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


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

leda_mutex segment_rep::mutex_id_counter;
unsigned long segment_rep::id_counter = 0;

segment_rep::segment_rep()  {id  = id_counter++; }

segment_rep::segment_rep(const point& p, const point& q) : start(p), end(q)
{ dx  = q.xcoord() - p.xcoord(); 
  dy  = q.ycoord() - p.ycoord(); 
  mutex_id_counter.lock();
  id  = id_counter++; 
  mutex_id_counter.unlock();
}


segment::segment() { PTR = new segment_rep; }

segment::segment(const point& x, const point& y) 
{ PTR = new segment_rep(x,y); }

segment::segment(const point& x, const vector& v) 
{ PTR = new segment_rep(x,x+v); }

segment::segment(double x1, double y1, double x2, double y2) 
{ PTR = new segment_rep(point(x1,y1), point(x2,y2)); }

segment::segment(const point& p, double alpha, double length)
{ point q = p.translate_by_angle(alpha,length);
  PTR  = new segment_rep(p,q); 
 }
  

segment segment::translate_by_angle(double alpha, double d) const
{ point p = ptr()->start.translate_by_angle(alpha,d);
  point q = ptr()->end.translate_by_angle(alpha,d);
  return segment(p,q);
 }


segment segment::translate(double dx, double dy) const
{ point p = ptr()->start.translate(dx,dy);
  point q = ptr()->end.translate(dx,dy);
  return segment(p,q);
 }


segment segment::translate(const vector& v) const
{ point p = ptr()->start.translate(v);
  point q = ptr()->end.translate(v);
  return segment(p,q);
 }


segment segment::rotate(double alpha) const
{ point p = start();
  point q = end();
  return segment(p,q.rotate(p,alpha));
}

segment segment::rotate(const point& origin, double alpha) const
{  point p = start().rotate(origin,alpha);
   point q = end().rotate(origin,alpha);
   return segment(p,q);
}

segment segment::rotate90() const
{ return segment(start(),point(xcoord1()-dy(),ycoord1()+dx())); }

segment segment::rotate90(const point& origin) const
{  return segment(start().rotate90(origin),end().rotate90(origin)); }


segment segment::reflect(const point& p, const point& q) const
{ return segment(start().reflect(p,q),end().reflect(p,q)); }

segment segment::reflect(const point& p) const
{ return segment(start().reflect(p),end().reflect(p)); }



double segment::slope() const 
{ double dx = ptr()->dx;
  double dy = ptr()->dy;
  return (dx == 0) ? MAXDOUBLE : dy/dx; 
}



double segment::sqr_length() const
{ double dx = ptr()->dx;
  double dy = ptr()->dy;
  return (dx*dx+dy*dy);
 }

double segment::length() const { return sqrt(sqr_length()); }


double segment::angle() const
{ double dx = ptr()->dx;
  double dy = ptr()->dy;
  if (dx != 0 || dy != 0)  
      return atan2(dy,dx); 
  else
      return 0;
}


double segment::angle(const segment& s) const
{ double dx1 =   ptr()->dx;
  double dy1 =   ptr()->dy;
  double dx2 = s.ptr()->dx;
  double dy2 = s.ptr()->dy;
  
  double norm  = (dx1*dx1+dy1*dy1)*(dx2*dx2+dy2*dy2);

  if (norm == 0) return 0;

  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;
}


double  segment::y_proj(double x)  const
{ return  ycoord1() - slope() * (xcoord1() - x); }

double  segment::x_proj(double y)  const
{ if (is_vertical())  
       return  xcoord1();
  else
       return  xcoord1() - (ycoord1() - y)/slope(); 
}

double segment::y_abs() const 
{ if (ptr()->dx == 0) 
     return -MAXDOUBLE;
  else
     return ycoord1() - slope()*xcoord1();
 }


bool segment::contains(const point& p) const
{ point a = source();
  point b = target();
  if (a == b) return a == p;
  return orientation(a,b,p) == 0 && (b-a)*(p-a) >=0 && (a-b)*(p-b) >= 0;
 }
 
 


bool segment::intersection(const segment& s) const
{ // decides whether this and |s| intersect 
  return (orientation(*this,s.start()) != orientation(*this,s.end()))
      && (orientation(s,start())       != orientation(s,end()));
}


bool segment::intersection(const segment& s, point& inter) const
{ 
  if (!intersection(s)) return false;

  point a = source();
  point b = target();
  point c = s.source();
  point d = s.target();

  if (is_trivial())
  { bool C = s.contains(a);
    if (C) inter = a;
    return C;
   }

  if ( s.is_trivial() )
  { bool C = contains(c);
    if (C) inter = c;
    return C;
   }


  if (start() == s.start() || start() == s.end())
  { inter = start();
    return true;
   }

  if (end() == s.start() || end() == s.end())
  { inter = end();
    return true;
   }


  double cx,cy;

  if (is_vertical())
     cx = xcoord1();
  else
     if (s.is_vertical())
        cx = s.xcoord1();
     else
        cx = (s.y_abs()-y_abs())/(slope()-s.slope());

  if (is_vertical())
     cy = s.slope() * cx + s.y_abs();
  else
     cy = slope() * cx + y_abs();

  inter = point(cx,cy);

  return true;
}


bool segment::intersection_of_lines(const segment& s, point& inter) const
{ double cx,cy;

  if (slope() == s.slope()) return false;

  if (start() == s.start() || start() == s.end())
  { inter = start();
    return true;
   }

  if (end() == s.start() || end() == s.end())
  { inter = end();
    return true;
   }


  if (is_vertical())
     cx = xcoord1();
  else
     if (s.is_vertical())
        cx = s.xcoord1();
     else
        cx = (s.y_abs()-y_abs())/(slope()-s.slope());

  if (is_vertical())
     cy = s.slope() * cx + s.y_abs();
  else
     cy = slope() * cx + y_abs();

  inter = point(cx,cy);

  return true;
}


segment segment::perpendicular(const point& q) const
{ point p = source();
  point r;
  intersection_of_lines(translate(q-p).rotate90(),r);
  return segment(q,r);
 }


double segment::distance(const point& p) const
{ point a = source();
  point b = target();
  if (orientation(a,b.rotate90(a),p) > 0) return a.distance(p);
  if (orientation(b,a.rotate90(b),p) > 0) return b.distance(p);
  return perpendicular(p).length();
}



ostream& operator<<(ostream& out, const segment& s) 
{ out << "[" << s.start() << "===" << s.end() << "]"; 
  return out;
 } 


istream& operator>>(istream& in, segment& s) 
{ // syntax: {[} p {===} q {]}

  point p,q; 
  char c;

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

  in >> p;

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

  in >> q; 

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

  s = segment(p,q); 
  return in; 

 } 


 
int cmp_segments_at_xcoord(const segment& s1, const segment& s2, const point& p)
{
  double s1yd = s2.dx()*(s1.dy()*(p.xcoord() - s1.xcoord1()) + 
			 s1.dx()*s1.ycoord1());

  double s2yd = s1.dx()*(s2.dy()*(p.xcoord() - s2.xcoord1()) + 
			 s2.dx()*s2.ycoord1());

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

