home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / ledar34.zip / leda-r-3_4_tar / LEDA-3.4 / src / plane / _circle.c next >
C/C++ Source or Header  |  1996-09-03  |  8KB  |  319 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA 3.4
  4. +
  5. +  _circle.c
  6. +
  7. +  This file is part of the LEDA research version (LEDA-R) that can be 
  8. +  used free of charge in academic research and teaching. Any commercial
  9. +  use of this software requires a license which is distributed by the
  10. +  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
  11. +  (fax +49 681 31104).
  12. +
  13. +  Copyright (c) 1991-1996  by  Max-Planck-Institut fuer Informatik
  14. +  Im Stadtwald, 66123 Saarbruecken, Germany     
  15. +  All rights reserved.
  16. + 
  17. *******************************************************************************/
  18.  
  19. #include <LEDA/circle.h>
  20. #include <math.h>
  21.  
  22.  
  23. //------------------------------------------------------------------------------
  24. // circles
  25. //
  26. // S. N"aher (1996)
  27. //------------------------------------------------------------------------------
  28.  
  29.  
  30. circle_rep::circle_rep(const point& p1, const point& p2, const point& p3) : 
  31. a(p1), b(p2), c(p3), cp(0), rp(0), first_side_of(true) 
  32. { orient = orientation(p1,p2,p3); }
  33.  
  34. circle_rep::~circle_rep() 
  35. { if (cp) delete cp;
  36.   if (rp) delete rp;
  37.  }
  38.  
  39.  
  40.  
  41. circle::circle(const point& a, const point& b, const point& c)
  42. { PTR = new circle_rep(a,b,c); 
  43.   if (a == b && a == c)
  44.   { ptr()->cp = new point(a);
  45.     ptr()->rp = new double(0);
  46.    }
  47. }
  48.  
  49.  
  50. circle::circle() 
  51. { PTR = new circle_rep(point(1,0),point(0,1),point(-1,0));
  52.   ptr()->cp = new point(0,0);
  53.   ptr()->rp = new double(1);
  54.  }
  55.  
  56. circle::circle(const point& m, double r)       
  57. { point a(m.xcoord(),  m.ycoord()+r); 
  58.   point b(m.xcoord(),  m.ycoord()-r); 
  59.   point c(m.xcoord()+r,m.ycoord()); 
  60.   PTR = new circle_rep(a,b,c); 
  61.   ptr()->cp = new point(m);
  62.   ptr()->rp = new double(r);
  63. }
  64.  
  65. circle::circle(double x, double y, double r)  
  66. { PTR = new circle_rep(point(x,y+r),point(x+r,y),point(x,y-r)); 
  67.   ptr()->cp = new point(x,y);
  68.   ptr()->rp = new double(r);
  69. }
  70.  
  71.  
  72.  
  73. point circle::center() const
  74. { if (ptr()->cp == 0)
  75.   { if (collinear(ptr()->a,ptr()->b,ptr()->c)) 
  76.        error_handler(1,"circle::center(): points are collinear.");
  77.     line l1 = p_bisector(ptr()->a,ptr()->b);
  78.     line l2 = p_bisector(ptr()->b,ptr()->c);
  79.     point m;
  80.     l1.intersection(l2,m);
  81.     ptr()->cp = new point(m);
  82.    }
  83.   return *(ptr()->cp);
  84.  }
  85.  
  86.  
  87. double circle::radius() const
  88. { if (ptr()->rp == 0)
  89.   { double r = center().distance(ptr()->a);
  90.     ptr()->rp = new double(r);
  91.    }
  92.  return *(ptr()->rp);
  93. }
  94.  
  95.   
  96.  
  97. bool circle::operator==(const circle& c)  const
  98. { if (!contains(c.ptr()->a)) return false;
  99.   if (!contains(c.ptr()->b)) return false;
  100.   if (!contains(c.ptr()->c)) return false;
  101.   return true;
  102. }
  103.  
  104.  
  105. double circle::distance(const point& p) const
  106. { double d = p.distance(center());
  107.   return (d - radius());
  108.  }
  109.  
  110.  
  111. double circle::distance(const line& l) const
  112. { double d = l.distance(center());
  113.   return (d - radius());
  114.  }
  115.  
  116. double circle::distance(const circle& c) const
  117. { double d = center().distance(c.center());
  118.   return (d - radius() - c.radius());
  119.  }
  120.  
  121. int circle::side_of(const point& p) const
  122. { 
  123.   double ax = ptr()->a.xcoord();
  124.   double ay = ptr()->a.ycoord();
  125.  
  126.   // compute D1 D2 D3 (if not done before)
  127.  
  128.   if (ptr()->first_side_of) 
  129.   { double bx = ptr()->b.xcoord() - ax;
  130.     double by = ptr()->b.ycoord() - ay;
  131.     double bw = bx*bx + by*by;
  132.  
  133.     double cx = ptr()->c.xcoord() - ax;
  134.     double cy = ptr()->c.ycoord() - ay;
  135.     double cw = cx*cx + cy*cy;
  136.  
  137.     ptr()->D1 = cy*bw - by*cw;
  138.     ptr()->D2 = bx*cw - cx*bw;
  139.     ptr()->D3 = by*cx - bx*cy;
  140.  
  141.     ptr()->first_side_of = false;
  142.    }
  143.  
  144.    double px = p.xcoord() - ax;
  145.    double py = p.ycoord() - ay;
  146.    double pw = px*px + py*py;
  147.  
  148.    double D  = px*ptr()->D1  + py*ptr()->D2 + pw*ptr()->D3;
  149.  
  150.    if (D != 0)
  151.       return (D > 0) ? 1 : -1;
  152.    else
  153.       return 0;
  154. }
  155.  
  156.  
  157.  
  158. bool circle::outside(const point& p) const 
  159. { return (ptr()->orient * side_of(p)) < 0; }
  160.  
  161.  
  162. bool circle::inside(const point& p) const 
  163. { return (ptr()->orient * side_of(p)) > 0; }
  164.  
  165.  
  166. bool circle::contains(const point& p) const 
  167. { return side_of(p) == 0; }
  168.  
  169.  
  170. circle circle::translate(double alpha, double d) const
  171. { point a1 = ptr()->a.translate(alpha,d);
  172.   point b1 = ptr()->b.translate(alpha,d);
  173.   point c1 = ptr()->c.translate(alpha,d);
  174.   return circle(a1,b1,c1);
  175.  }
  176.  
  177. circle circle::translate(const vector& v) const
  178. { point a1 = ptr()->a.translate(v);
  179.   point b1 = ptr()->b.translate(v);
  180.   point c1 = ptr()->c.translate(v);
  181.   return circle(a1,b1,c1);
  182.  }
  183.  
  184. circle  circle::rotate(const point& o, double alpha) const
  185. { point a1 = ptr()->a.rotate(o,alpha);
  186.   point b1 = ptr()->b.rotate(o,alpha);
  187.   point c1 = ptr()->c.rotate(o,alpha);
  188.   return circle(a1,b1,c1);
  189.  }
  190.  
  191. circle  circle::rotate(double alpha)  const
  192. { point o(0,0);
  193.   point a1 = ptr()->a.rotate(o,alpha);
  194.   point b1 = ptr()->b.rotate(o,alpha);
  195.   point c1 = ptr()->c.rotate(o,alpha);
  196.   return circle(a1,b1,c1);
  197. }
  198.  
  199.  
  200. circle  circle::rotate90(const point& o)  const
  201. { point a1 = ptr()->a.rotate90(o);
  202.   point b1 = ptr()->b.rotate90(o);
  203.   point c1 = ptr()->c.rotate90(o);
  204.   return circle(a1,b1,c1);
  205. }
  206.  
  207.  
  208. circle  circle::reflect(const point& p, const point& q)  const
  209. { point a1 = ptr()->a.reflect(p,q);
  210.   point b1 = ptr()->b.reflect(p,q);
  211.   point c1 = ptr()->c.reflect(p,q);
  212.   return circle(a1,b1,c1);
  213. }
  214.  
  215. circle  circle::reflect(const point& p)  const
  216. { point a1 = ptr()->a.reflect(p);
  217.   point b1 = ptr()->b.reflect(p);
  218.   point c1 = ptr()->c.reflect(p);
  219.   return circle(a1,b1,c1);
  220. }
  221.  
  222.  
  223.  
  224. list<point> circle::intersection(const line& l) const
  225. { list<point> result;
  226.  
  227.   segment s = l.perpendicular(center());
  228.   double  d = s.length();
  229.   double  r = radius();
  230.  
  231.   if (d==r) result.append(s.end());
  232.   
  233.   if (d < r)
  234.   { segment sl = l.ptr()->seg;
  235.     double t = sqrt(r*r - d*d)/sl.length();
  236.     double dx = t*sl.dx();
  237.     double dy = t*sl.dy();
  238.     result.append(s.end().translate(vector(dx,dy)));
  239.     result.append(s.end().translate(vector(-dx,-dy)));
  240.   }
  241.  
  242.   return result;
  243. }
  244.  
  245.  
  246. list<point> circle::intersection(const segment& s) const
  247. { list<point> result,L;
  248.  
  249.   L = intersection(line(s));
  250.  
  251.   point p;
  252.   double  d  = s.length();
  253.  
  254.   forall(p,L)
  255.   { double d1 = s.ptr()->start.distance(p);
  256.     double d2 = s.ptr()->end.distance(p);
  257.     if (d1 <= d && d2 <= d) result.append(p);
  258.    }
  259.  
  260.   return result;
  261. }
  262.  
  263. list<point> circle::intersection(const circle& c) const
  264. { list<point> result;
  265.   segment s(center(), c.center());
  266.   double d  = s.length();
  267.   double r1 = radius();
  268.   double r2 = c.radius();
  269.  
  270.   if (d > (r1+r2) || (d+r2) < r1 || (d+r1) < r2) return result;
  271.  
  272.  
  273.   double x = (d*d + r1*r1 - r2*r2)/(2*d);
  274.   double alpha = acos(x/r1);
  275.   double beta  = s.angle() + alpha;
  276.   double gamma = s.angle() - alpha;
  277.  
  278.   result.append(center().translate(beta,r1));
  279.   if (alpha!=0) result.append(center().translate(gamma,r1));
  280.  
  281.   return result;
  282.  }
  283.  
  284.  
  285. segment circle::left_tangent(const point& p) const
  286. { if (inside(p)) error_handler(1,"left_tangent:: point inside circle");
  287.   segment s(p,center());
  288.   double d = s.length();
  289.   double alpha = asin(radius()/d) + s.angle();
  290.   point touch = p.translate(alpha,sqrt(d*d - radius()*radius()));
  291.   return segment(p,touch);
  292. }
  293.  
  294. segment circle::right_tangent(const point& p) const
  295. { if (inside(p)) error_handler(1,"right_tangent:: point inside circle");
  296.   segment s(p,center());
  297.   double d = s.length();
  298.   double alpha = s.angle() - asin(radius()/d);
  299.   point touch = p.translate(alpha,sqrt(d*d - radius()*radius()));
  300.   return segment(p,touch);
  301. }
  302.  
  303.  
  304. ostream& operator<<(ostream& out, const circle& c) 
  305. { out << c.center() << " "<< c.radius(); 
  306.   return out;
  307.  } 
  308.  
  309. istream& operator>>(istream& in,  circle& c) 
  310. { point cent;
  311.   double rad;
  312.   if (in) in >> cent;
  313.   if (in) in >> rad;
  314.   c = circle(cent,rad);
  315.   return in;
  316. }
  317.  
  318.  
  319.