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

//------------------------------------------------------------------------------
// points 
//
// by S. Naeher (1995)
//------------------------------------------------------------------------------

leda_mutex point_rep::mutex_id_counter;

unsigned long point_rep::id_counter = 0;

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


point::point()                            { PTR = new point_rep; }
point::point(double x, double y)          { PTR = new point_rep(x,y); }
point::point(double x, double y, double w){ PTR = new point_rep(x/w,y/w); }
point::point(vector v)                    { PTR = new point_rep(v[0], v[1]); }

double point::angle(const point& q, const point& r) const
{
  double dx1 = q.ptr()->x - ptr()->x; 
  double dy1 = q.ptr()->y - ptr()->y; 
  double dx2 = r.ptr()->x - ptr()->x; 
  double dy2 = r.ptr()->y - ptr()->y; 

  if ((dx1 == 0 && dy1 == 0) || (dx2 == 0 && dy2 == 0)) 
     error_handler(1,"point::angle:  zero vector input.");
  
  double norm  = (dx1*dx1+dy1*dy1)*(dx2*dx2+dy2*dy2);

  double cosfi = (dx1*dx2+dy1*dy2) / sqrt(norm);

  if (cosfi >=  1.0 ) return 0;
  if (cosfi <= -1.0 ) return LEDA_PI;
  
  double fi = acos(cosfi);

  return (dx1*dy2 > dy1*dx2) ? fi : -fi;

}
  


// Rotations 

point point::rotate90(const point& p) const
{ double px = p.xcoord();
  double py = p.ycoord();
  double dx = xcoord() - px;
  double dy = ycoord() - py;
  return point(px-dy,py+dx);
}

point point::rotate90() const
{ return point(-ycoord(),xcoord()); }


point point::rotate(const point& origin, double fi) const
{ double cx = origin.xcoord();
  double cy = origin.ycoord();
  double sinfi = sin(fi);
  double cosfi = cos(fi);
  double dx = xcoord() - cx;
  double dy = ycoord() - cy;
  return point(cx+dx*cosfi-dy*sinfi,cy+dx*sinfi+dy*cosfi);
}


point point::rotate(double fi) const
{ double sinfi = sin(fi);
  double cosfi = cos(fi);
  double x = xcoord();
  double y = ycoord();
  return point(x*cosfi-y*sinfi,x*sinfi+y*cosfi);
}


point point::reflect(const point& p, const point& q) const
{ // reflect point across line through p and q

  double px = p.xcoord();
  double py = p.ycoord();

  double x1 = xcoord()   - px;
  double y1 = ycoord()   - py;
  double x2 = q.xcoord() - px;
  double y2 = q.ycoord() - py;

  double L = (x1*x1 + y1*y1) * (x2*x2 + y2*y2);

  double cosfi = (x1*x2 + y1*y2);
  double sinfi = (x1*y2 - x2*y1); 
  double cos2 = (cosfi*cosfi - sinfi*sinfi)/L;
  double sin2 = 2*cosfi*sinfi/L;

  return point(px + x1*cos2-y1*sin2, py + x1*sin2+y1*cos2);
}


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


// Translations

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


point point::translate_by_angle(double phi, double d) const
{ double dx = cos(phi) * d;
  double dy = sin(phi) * d;
  if (fabs(dx) < 1e-12) dx = 0; 
  if (fabs(dy) < 1e-12) dy = 0; 
  return point(xcoord()+dx,ycoord()+dy);
 }

point point::translate(const vector& v) const 
{ return point(xcoord()+v[0],ycoord()+v[1]); }


// Distances

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

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

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

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


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

   

int side_of_circle(const point& a, const point& b, const point& c, 
                                                   const point& d)
{  double ax = a.xcoord();
   double ay = a.ycoord();

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

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

   double D1 = cy*bw - by*cw;
   double D2 = bx*cw - cx*bw;
   double D3 = by*cx - bx*cy;


   double dx = d.xcoord() - ax;
   double dy = d.ycoord() - ay;

   double D  = D1*dx  + D2*dy + D3*(dx*dx + dy*dy);


   if (D != 0) 
      return (D > 0) ? 1 : -1;
   else
      return 0;
 }




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

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

  double x,y; 
  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 (c == ' ');
  if (c != ')') in.putback(c);

  p = point(x,y); 
  return in; 

 } 


int point::cmp_xy(const point& a, const point& b)
{ int r = compare(a.xcoord(),b.xcoord());
  return (r!=0) ? r : compare(a.ycoord(),b.ycoord());
 }

int point::cmp_yx(const point& a, const point& b)
{ int r = compare(a.ycoord(),b.ycoord());
  return (r!=0) ? r : compare(a.xcoord(),b.xcoord());
 }
