/*******************************************************************************
+
+  LEDA 3.5
+
+  _d3_plane.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
#include <LEDA/d3_plane.h>
#include <math.h>
#include <ctype.h>

//------------------------------------------------------------------------------
// d3_plane 
//
// by S. Naeher (1997)
//------------------------------------------------------------------------------


leda_mutex d3_plane_rep::mutex_id_counter;

unsigned long d3_plane_rep::id_counter = 0;


d3_plane_rep::d3_plane_rep(const d3_point& p, const double& x, const double& y,
                                                               const double& z)
: a(p),nx(x),ny(y),nz(z)

{ if (nx == 0 && ny == 0 && nz == 0)
     error_handler(1,"d3_plane::d3_plane(): cannot construct plane.");
  mutex_id_counter.lock();
  id  = id_counter++; 
  mutex_id_counter.unlock();
}


d3_plane::d3_plane(const d3_point& a, const d3_point& b, const d3_point& c)
{ double X1 = b.X() - a.X();
  double Y1 = b.Y() - a.Y();
  double Z1 = b.Z() - a.Z();
  double X2 = c.X() - a.X();
  double Y2 = c.Y() - a.Y();
  double Z2 = c.Z() - a.Z();
  double nx = Z1*Y2 - Y1*Z2;
  double ny = X1*Z2 - Z1*X2;
  double nz = Y1*X2 - X1*Y2;
  PTR = new d3_plane_rep(a,nx,ny,nz); 
}


d3_plane::d3_plane(const d3_point& p, const vector& v) 
{ PTR = new d3_plane_rep(p,v.hcoord(0),v.hcoord(1),v.hcoord(2)); }


d3_plane d3_plane::translate(double dx, double dy, double dz) const
{ d3_point p = point1().translate(dx,dy,dz);
  return d3_plane(p,normal());
}


d3_plane d3_plane::translate(const vector& v) const 
{ d3_point p = point1().translate(v);
  return d3_plane(p,normal());
}


// Distances

d3_plane d3_plane::reflect(const d3_point& q) const
{ // reflect plane across point q
  d3_point a = point1();
  d3_point b = a.translate(normal());
  d3_point a1 = a.reflect(q);
  d3_point b1 = b.reflect(q);
  return d3_plane(a1,b1-a1);
 }


vector d3_plane::normal_project(const d3_point& p) const
{ vector v = p - point1();
  double nx = ptr()->nx;
  double ny = ptr()->ny;
  double nz = ptr()->nz;

  double W = nx*nx+ny*ny+nz*nz;
  double A = -(nx*v.hcoord(0)+ny*v.hcoord(1)+nz*v.hcoord(2))/W;
  return vector(A*nx,A*ny,A*nz);
}

d3_point d3_plane::reflect_point(const d3_point& p) const
{ return p.translate(2*normal_project(p)); } 


d3_plane  d3_plane::reflect(const d3_plane& Q) const
{ // reflect plane across plane q
  d3_point a = point1();
  d3_point b = a.translate(normal());
  d3_point a1 = Q.reflect_point(a);
  d3_point b1 = Q.reflect_point(b);
  return d3_plane(a1,b1-a1);
 }


double d3_plane::sqr_dist(const d3_point& p)  const
{ vector vec = normal_project(p);
  return vec.sqr_length();
 }

double d3_plane::distance(const d3_point& p) const 
{ return sqrt(sqr_dist(p)); }


bool d3_plane::parallel(const d3_plane& Q) const
{ double nx1 = ptr()->nx;
  double ny1 = ptr()->ny;
  double nz1 = ptr()->nz;
  double nx2 = Q.ptr()->nx;
  double ny2 = Q.ptr()->ny;
  double nz2 = Q.ptr()->nz;
  return nz1*ny2 == ny1*nz2 && nx1*nz2 == nz1*nx2 && ny1*nx2 == nx1*ny2;
}


int d3_plane::side_of(const d3_point& q) const
{ 
  d3_point a = point1();

  double qx = q.X()*a.W() - a.X()*q.W();
  double qy = q.Y()*a.W() - a.Y()*q.W();
  double qz = q.Z()*a.W() - a.Z()*q.W();

  double nx = ptr()->nx;
  double ny = ptr()->ny;
  double nz = ptr()->nz;

  double D = nx*qx + ny*qy + nz*qz;

  if (D > 0) return +1;
  if (D < 0) return -1;
  return 0;
}


d3_plane d3_plane::to_d3_plane() const
{ return d3_plane(point1(),normal()); }


int d3_plane::operator==(const d3_plane& P) const 
{ return parallel(P) && contains(P.point1()); }

   
ostream& operator<<(ostream& out, const d3_plane& P)
{ cout << P.point1() << " " << P.normal(); return out; } 

istream& operator>>(istream& in, d3_plane& P) 
{ d3_point p;
  vector n(3);
  cin >> p >> n;
  P = d3_plane(p,n);
  return in; 
} 
