/*******************************************************************************
+
+  LEDA 3.5
+
+  _d3_point.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_point.h>
#include <math.h>
#include <ctype.h>

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

leda_mutex d3_point_rep::mutex_id_counter;

unsigned long d3_point_rep::id_counter = 0;

d3_point_rep::d3_point_rep(double a, double b, double c) 
{ x = a; 
  y = b; 
  z = c; 
  mutex_id_counter.lock();
  id  = id_counter++; 
  mutex_id_counter.unlock();
}


d3_point::d3_point() { PTR = new d3_point_rep; }

d3_point::d3_point(double x, double y, double z) 
{ PTR = new d3_point_rep(x,y,z); }

d3_point::d3_point(double x, double y, double z,double w)
{ PTR = new d3_point_rep(x/w,y/w,z/w); }

d3_point::d3_point(vector v) 
{ PTR = new d3_point_rep(v[0], v[1], v[2]); }

// Projections

point d3_point::project_xy() const { return point(xcoord(),ycoord()); }
point d3_point::project_yz() const { return point(ycoord(),zcoord()); }
point d3_point::project_xz() const { return point(xcoord(),zcoord()); }


// Translations

d3_point d3_point::translate(double dx, double dy, double dz) const
{ return d3_point(xcoord()+dx,ycoord()+dy,zcoord()+dz); }


d3_point d3_point::translate(const vector& v) const 
{ return translate(v[0],v[1],v[2]); }


// Distances

double d3_point::sqr_dist(const d3_point& p)  const
{ double dx = p.ptr()->x - ptr()->x; 
  double dy = p.ptr()->y - ptr()->y;
  double dz = p.ptr()->z - ptr()->z;
  return dx*dx + dy*dy + dz*dz;
 }

double d3_point::xdist(const d3_point& q) const 
{ return fabs(xcoord() - q.xcoord()); }

double d3_point::ydist(const d3_point& q) const 
{ return fabs(ycoord() - q.ycoord()); }

double d3_point::zdist(const d3_point& q) const 
{ return fabs(zcoord() - q.zcoord()); }


d3_point center(const d3_point& a, const d3_point& b)
{ double x = (a.xcoord()+b.xcoord())/2;
  double y = (a.ycoord()+b.ycoord())/2;
  double z = (a.zcoord()+b.zcoord())/2;
  return d3_point(x,y,z);
 } 


d3_point d3_point::reflect(const d3_point& q) const
{ // reflect point across point q
  return d3_point(2*q.xcoord()-xcoord(), 2*q.ycoord()-ycoord(),
                                         2*q.zcoord()-zcoord());
 }


d3_point  d3_point::reflect(const d3_point& a, const d3_point& b,
                                               const d3_point& c) const
{  
  // reflect point across plane through a, b, and c

  double x1 = b.xcoord() - a.xcoord();
  double y1 = b.ycoord() - a.ycoord();
  double z1 = b.zcoord() - a.zcoord();

  double x2 = c.xcoord() - a.xcoord();
  double y2 = c.ycoord() - a.ycoord();
  double z2 = c.zcoord() - a.zcoord();

  double x3 = xcoord() - a.xcoord();
  double y3 = ycoord() - a.ycoord();
  double z3 = zcoord() - a.zcoord();

  double x = (z1*y2-y1*z2);
  double y = (x1*z2-z1*x2);
  double z = (y1*x2-x1*y2);

  if (x == 0 && y == 0 && z == 0)
      error_handler(1,"d3_point::reflect(a,b,c): a,b,c are coplanar");


  double f = -2*(x*x3+y*y3+z*z3)/(x*x+y*y+z*z);

  return translate(f*x,f*y,f*z);
}

  


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


int d3_point::operator==(const d3_point& p) const 
{ return ptr()->x == p.ptr()->x && 
         ptr()->y == p.ptr()->y && 
         ptr()->z == p.ptr()->z; 
}

   

int side_of_sphere(const d3_point& a, const d3_point& b, const d3_point& c, 
                                                         const d3_point& d,
                                                         const d3_point& e)

{  double ax = a.xcoord();
   double ay = a.ycoord();
   double az = a.zcoord();

   double bx = b.xcoord() - ax;
   double by = b.ycoord() - ay;
   double bz = b.zcoord() - az;
   double bw = bx*bx + by*by + bz*bz;

   double cx = c.xcoord() - ax;
   double cy = c.ycoord() - ay;
   double cz = c.zcoord() - az;
   double cw = cx*cx + cy*cy + cz*cz;

   double dx = d.xcoord() - ax;
   double dy = d.ycoord() - ay;
   double dz = d.zcoord() - az;
   double dw = dx*dx + dy*dy + dz*dz;

   double ex = e.xcoord() - ax;
   double ey = e.ycoord() - ay;
   double ez = e.zcoord() - az;
   double ew = ex*ex + ey*ey + ez*ez;

   return -orientation(d3_point(bx/bw,by/bw,bz/bw),
                       d3_point(cx/cw,cy/cw,cz/cw),
                       d3_point(dx/dw,dy/dw,dz/dw),
                       d3_point(ex/ew,ey/ew,ez/ew));
 }





ostream& operator<<(ostream& out, const d3_point& p)
{ out << "(" << p.xcoord() << "," << p.ycoord() << "," << p.zcoord() << ")";
  return out;
 } 

istream& operator>>(istream& in, d3_point& p) 
{ // syntax: {(} x {,} y {,} z{)}

  double x,y,z; 
  char c;

  do in.get(c); while (in && isspace(c));

  if (!in) return in;

  if (c != '(') in.putback(c);

  in >> x;

  do in.get(c); while (isspace(c));
  if (c != ',') in.putback(c);

  in >> y; 

  do in.get(c); while (isspace(c));
  if (c != ',') in.putback(c);

  in >> z; 

  do in.get(c); while (c == ' ');
  if (c != ')') in.putback(c);

  p = d3_point(x,y,z); 
  return in; 

 } 


int d3_point::cmp_xyz(const d3_point& a, const d3_point& b)
{ int d = compare(a.xcoord(),b.xcoord());
  if (d) return d;
  d = compare(a.ycoord(),b.ycoord());
  if (d) return d;
  d = compare(a.zcoord(),b.zcoord());
  return d;
 }



double volume(const d3_point& a, const d3_point& b, const d3_point& c,
                                                    const d3_point& d)
{ double AX = a.xcoord();
  double AY = a.ycoord();
  double AZ = a.zcoord();
  double BX = b.xcoord();
  double BY = b.ycoord();
  double BZ = b.zcoord();
  double CX = c.xcoord();
  double CY = c.ycoord();
  double CZ = c.zcoord();
  double DX = d.xcoord();
  double DY = d.ycoord();
  double DZ = d.zcoord();
   
  double D = - AZ*BY*CX + AY*BZ*CX + AZ*BX*CY - AX*BZ*CY 
             - AY*BX*CZ + AX*BY*CZ + AZ*BY*DX - AY*BZ*DX 
             - AZ*CY*DX + BZ*CY*DX + AY*CZ*DX - BY*CZ*DX 
             - AZ*BX*DY + AX*BZ*DY + AZ*CX*DY - BZ*CX*DY 
             - AX*CZ*DY + BX*CZ*DY + AY*BX*DZ - AX*BY*DZ 
             - AY*CX*DZ + BY*CX*DZ + AX*CY*DZ - BX*CY*DZ;

  return D/6;
}



int orientation(const d3_point& a, const d3_point& b, const d3_point& c,
                                                      const d3_point& d)
{ //sign of (d-a) * cross_prod(b-a,c-a)  

  double x1 = b.xcoord() - a.xcoord();
  double y1 = b.ycoord() - a.ycoord();
  double z1 = b.zcoord() - a.zcoord();

  double x2 = c.xcoord() - a.xcoord();
  double y2 = c.ycoord() - a.ycoord();
  double z2 = c.zcoord() - a.zcoord();

  double x3 = d.xcoord() - a.xcoord();
  double y3 = d.ycoord() - a.ycoord();
  double z3 = d.zcoord() - a.zcoord();
    
  double D = (z1*y2-y1*z2)*x3 + (x1*z2-z1*x2)*y3 + (y1*x2-x1*y2)*z3;

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



bool collinear(const d3_point& a, const d3_point& b, const d3_point& c)
{ // cross product(b-a, c-a) == 0
  double X1 = b.xcoord() - a.xcoord();
  double Y1 = b.ycoord() - a.ycoord();
  double Z1 = b.zcoord() - a.zcoord();
  double X2 = c.xcoord() - a.xcoord();
  double Y2 = c.ycoord() - a.ycoord();
  double Z2 = c.zcoord() - a.zcoord();

  return  Z1*Y2 == Y1*Z2 &&
          X1*Z2 == Z1*X2 &&
          Y1*X2 == X1*Y2;
 }




inline int d2_orientation(double ax,double ay, double bx,double by,
                                               double cx,double cy)
{ double E = (ax-bx) * (ay-cy) - (ay-by) * (ax-cx);
  if (E > 0) return  1;
  else if (E < 0) return -1;
  else return  0;
}
 

int orientation_xy(const d3_point& a, const d3_point& b, const d3_point& c)
{ return d2_orientation(a.xcoord(),a.ycoord(),b.xcoord(),b.ycoord(), 
                                              c.xcoord(),c.ycoord()); }

int orientation_yz(const d3_point& a, const d3_point& b, const d3_point& c)
{ return d2_orientation(a.ycoord(),a.zcoord(),b.ycoord(),b.zcoord(), 
                                              c.ycoord(),c.zcoord()); }

int orientation_xz(const d3_point& a, const d3_point& b, const d3_point& c)
{ return d2_orientation(a.xcoord(),a.zcoord(),b.xcoord(),b.zcoord(), 
                                              c.xcoord(),c.zcoord()); }



d3_point point_on_positive_side(const d3_point& a, const d3_point& b, 
                                                   const d3_point& c)
{ 
  double x1 = b.xcoord() - a.xcoord();
  double y1 = b.ycoord() - a.ycoord();
  double z1 = b.zcoord() - a.zcoord();

  double x2 = c.xcoord() - a.xcoord();
  double y2 = c.ycoord() - a.ycoord();
  double z2 = c.zcoord() - a.zcoord();

  double cx = z1*y2-y1*z2;
  double cy = x1*z2-z1*x2; 
  double cz = y1*x2-x1*y2;

  return a.translate(cx,cy,cz);
}

