compresse geometry, add test, and move HalfPlaneIntersection

code/Geometry.cc → code/Geometry.cpp

@@ -11,27 +11,23 @@ struct PT {
   PT operator * (double c)     const { return PT(x*c,   y*c  ); }
   PT operator / (double c)     const { return PT(x/c,   y/c  ); }
 };
-
 double dot(PT p, PT q)     { return p.x*q.x+p.y*q.y; }
 double dist2(PT p, PT q)   { return dot(p-q,p-q); }
 double cross(PT p, PT q)   { return p.x*q.y-p.y*q.x; }
 ostream &operator<<(ostream &os, const PT &p) {
   os << "(" << p.x << "," << p.y << ")";
 }
-
 // rotate a point CCW or CW around the origin
 PT RotateCCW90(PT p)   { return PT(-p.y,p.x); }
 PT RotateCW90(PT p)    { return PT(p.y,-p.x); }
 PT RotateCCW(PT p, double t) {
   return PT(p.x*cos(t)-p.y*sin(t), p.x*sin(t)+p.y*cos(t));
 }
-
 // project point c onto line through a and b
 // assuming a != b
 PT ProjectPointLine(PT a, PT b, PT c) {
   return a + (b-a)*dot(c-a, b-a)/dot(b-a, b-a);
 }
-
 // project point c onto line segment through a and b
 PT ProjectPointSegment(PT a, PT b, PT c) {
   double r = dot(b-a,b-a);
@@ -41,30 +37,25 @@ PT ProjectPointSegment(PT a, PT b, PT c) {
   if (r > 1) return b;
   return a + (b-a)*r;
 }
-
 // compute distance from c to segment between a and b
 double DistancePointSegment(PT a, PT b, PT c) {
   return sqrt(dist2(c, ProjectPointSegment(a, b, c)));
 }
-
 // compute distance between point (x,y,z) and plane ax+by+cz=d
 double DistancePointPlane(double x, double y, double z,
                           double a, double b, double c, double d)
 {
   return fabs(a*x+b*y+c*z-d)/sqrt(a*a+b*b+c*c);
 }
-
 // determine if lines from a to b and c to d are parallel or collinear
 bool LinesParallel(PT a, PT b, PT c, PT d) {
   return fabs(cross(b-a, c-d)) < EPS;
 }
-
 bool LinesCollinear(PT a, PT b, PT c, PT d) {
   return LinesParallel(a, b, c, d)
       && fabs(cross(a-b, a-c)) < EPS
       && fabs(cross(c-d, c-a)) < EPS;
 }
-
 // determine if line segment from a to b intersects with
 // line segment from c to d
 bool SegmentsIntersect(PT a, PT b, PT c, PT d) {
@@ -79,7 +70,6 @@ bool SegmentsIntersect(PT a, PT b, PT c, PT d) {
   if (cross(a-c, d-c) * cross(b-c, d-c) > 0) return false;
   return true;
 }
-
 // compute intersection of line passing through a and b
 // with line passing through c and d, assuming that unique
 // intersection exists; for segment intersection, check if
@@ -89,14 +79,12 @@ PT ComputeLineIntersection(PT a, PT b, PT c, PT d) {
   assert(dot(b, b) > EPS && dot(d, d) > EPS);
   return a + b*cross(c, d)/cross(b, d);
 }
-
 // compute center of circle given three points
 PT ComputeCircleCenter(PT a, PT b, PT c) {
   b=(a+b)/2;
   c=(a+c)/2;
   return ComputeLineIntersection(b, b+RotateCW90(a-b), c, c+RotateCW90(a-c));
 }
-
 // determine if point is in a possibly non-convex polygon (by William
 // Randolph Franklin); returns 1 for strictly interior points, 0 for
 // strictly exterior points, and 0 or 1 for the remaining points.
@@ -115,7 +103,6 @@ bool PointInPolygon(const vector<PT> &p, PT q) {
   }
   return c;
 }
-
 // determine if point is on the boundary of a polygon
 bool PointOnPolygon(const vector<PT> &p, PT q) {
   for (int i = 0; i < p.size(); i++)
@@ -123,7 +110,6 @@ bool PointOnPolygon(const vector<PT> &p, PT q) {
       return true;
     return false;
 }
-
 // compute intersection of line through points a and b with
 // circle centered at c with radius r > 0
 vector<PT> CircleLineIntersection(PT a, PT b, PT c, double r) {
@@ -140,7 +126,6 @@ vector<PT> CircleLineIntersection(PT a, PT b, PT c, double r) {
     ret.push_back(c+a+b*(-B-sqrt(D))/A);
   return ret;
 }
-
 // compute intersection of circle centered at a with radius r
 // with circle centered at b with radius R
 vector<PT> CircleCircleIntersection(PT a, PT b, double r, double R) {
@@ -155,7 +140,6 @@ vector<PT> CircleCircleIntersection(PT a, PT b, double r, double R) {
     ret.push_back(a+v*x - RotateCCW90(v)*y);
   return ret;
 }
-
 // This code computes the area or centroid of a (possibly nonconvex)
 // polygon, assuming that the coordinates are listed in a clockwise or
 // counterclockwise fashion.  Note that the centroid is often known as
@@ -168,11 +152,9 @@ double ComputeSignedArea(const vector<PT> &p) {
   }
   return area / 2.0;
 }
-
 double ComputeArea(const vector<PT> &p) {
   return fabs(ComputeSignedArea(p));
 }
-
 PT ComputeCentroid(const vector<PT> &p) {
   PT c(0,0);
   double scale = 6.0 * ComputeSignedArea(p);
@@ -182,7 +164,6 @@ PT ComputeCentroid(const vector<PT> &p) {
   }
   return c / scale;
 }
-
 // tests whether or not a given polygon (in CW or CCW order) is simple
 bool IsSimple(const vector<PT> &p) {
   for (int i = 0; i < p.size(); i++) {
@@ -196,71 +177,56 @@ bool IsSimple(const vector<PT> &p) {
   }
   return true;
 }
-
-int main() {
-
+// angle of vector y relative to x
+double angle(PT x, PT y) {
+  return atan2(y.y,y.x) - atan2(x.y,y.x);
+}
+int test() {
   // expected: (-5,2)
   cerr << RotateCCW90(PT(2,5)) << endl;
-
   // expected: (5,-2)
   cerr << RotateCW90(PT(2,5)) << endl;
-
   // expected: (-5,2)
   cerr << RotateCCW(PT(2,5),M_PI/2) << endl;
-
   // expected: (5,2)
   cerr << ProjectPointLine(PT(-5,-2), PT(10,4), PT(3,7)) << endl;
-
   // expected: (5,2) (7.5,3) (2.5,1)
   cerr << ProjectPointSegment(PT(-5,-2), PT(10,4), PT(3,7)) << " "
        << ProjectPointSegment(PT(7.5,3), PT(10,4), PT(3,7)) << " "
        << ProjectPointSegment(PT(-5,-2), PT(2.5,1), PT(3,7)) << endl;
-
   // expected: 6.78903
   cerr << DistancePointPlane(4,-4,3,2,-2,5,-8) << endl;
-
   // expected: 1 0 1
   cerr << LinesParallel(PT(1,1), PT(3,5), PT(2,1), PT(4,5)) << " "
        << LinesParallel(PT(1,1), PT(3,5), PT(2,0), PT(4,5)) << " "
        << LinesParallel(PT(1,1), PT(3,5), PT(5,9), PT(7,13)) << endl;
-
   // expected: 0 0 1
   cerr << LinesCollinear(PT(1,1), PT(3,5), PT(2,1), PT(4,5)) << " "
        << LinesCollinear(PT(1,1), PT(3,5), PT(2,0), PT(4,5)) << " "
        << LinesCollinear(PT(1,1), PT(3,5), PT(5,9), PT(7,13)) << endl;
-
   // expected: 1 1 1 0
   cerr << SegmentsIntersect(PT(0,0), PT(2,4), PT(3,1), PT(-1,3)) << " "
        << SegmentsIntersect(PT(0,0), PT(2,4), PT(4,3), PT(0,5)) << " "
        << SegmentsIntersect(PT(0,0), PT(2,4), PT(2,-1), PT(-2,1)) << " "
        << SegmentsIntersect(PT(0,0), PT(2,4), PT(5,5), PT(1,7)) << endl;
-
   // expected: (1,2)
   cerr << ComputeLineIntersection(PT(0,0), PT(2,4), PT(3,1), PT(-1,3)) << endl;
-
   // expected: (1,1)
   cerr << ComputeCircleCenter(PT(-3,4), PT(6,1), PT(4,5)) << endl;
-
   vector<PT> v;
-  v.push_back(PT(0,0));
-  v.push_back(PT(5,0));
-  v.push_back(PT(5,5));
-  v.push_back(PT(0,5));
-
+  v.push_back(PT(0,0)); v.push_back(PT(5,0)); v.push_back(PT(5,5)); v.push_back(PT(0,5));
   // expected: 1 1 1 0 0
   cerr << PointInPolygon(v, PT(2,2)) << " "
        << PointInPolygon(v, PT(2,0)) << " "
        << PointInPolygon(v, PT(0,2)) << " "
        << PointInPolygon(v, PT(5,2)) << " "
        << PointInPolygon(v, PT(2,5)) << endl;
-
   // expected: 0 1 1 1 1
   cerr << PointOnPolygon(v, PT(2,2)) << " "
        << PointOnPolygon(v, PT(2,0)) << " "
        << PointOnPolygon(v, PT(0,2)) << " "
        << PointOnPolygon(v, PT(5,2)) << " "
        << PointOnPolygon(v, PT(2,5)) << endl;
-
   // expected: (1,6)
   //           (5,4) (4,5)
   //           blank line
@@ -279,7 +245,6 @@ int main() {
   for (int i = 0; i < u.size(); i++) cerr << u[i] << " "; cerr << endl;
   u = CircleCircleIntersection(PT(1,1), PT(4.5,4.5), 5, sqrt(2.0)/2.0);
   for (int i = 0; i < u.size(); i++) cerr << u[i] << " "; cerr << endl;
-
   // area should be 5.0
   // centroid should be (1.1666666, 1.166666)
   PT pa[] = { PT(0,0), PT(5,0), PT(1,1), PT(0,5) };
@@ -287,57 +252,5 @@ int main() {
   PT c = ComputeCentroid(p);
   cerr << "Area: " << ComputeArea(p) << endl;
   cerr << "Centroid: " << c << endl;
-
   return 0;
 }
-
-// angle of vector y relative to x
-double angle(PT x, PT y) {
-  return atan2(y.y,y.x) - atan2(x.y,y.x);
-}
-
-// INTERSECTION DE DEMI-PLAN
-// demi-plan : à gauche de la ligne
-struct Line {
-  PT P;    // point
-  PT v;   // vecteur directeur
-  double ang;
-  Line() {}
-  Line(PT P, PT v):P(P),v(v){ ang = atan2(v.y, v.x); }
-  bool operator < (const Line& L) const {
-    return ang < L.ang;
-  }
-};
-bool OnLeft(const Line& L, const PT& p) {
-  return cross(L.v, p-L.P) > 0;
-}
-PT GetLineIntersection(const Line& a, const Line& b) {
-  return ComputeLineIntersection(a.P, a.P+a.v, b.P, b.P+b.v);
-}
-// INPUT : set of half-plane
-// OUPUT : polygon (if exists)
-// COMPLEXITY n log n ?
-vector<PT> HalfplaneIntersection(vector<Line> L) {
-  int n = L.size();
-  sort(L.begin(), L.end());
-  int first, last;
-  vector<PT> p(n);
-  vector<Line> q(n);
-  vector<PT> ans;
-  q[first=last=0] = L[0];
-  for(int i = 1; i < n; i++) {
-    while(first < last && !OnLeft(L[i], p[last-1])) last--;
-    while(first < last && !OnLeft(L[i], p[first])) first++;
-    q[++last] = L[i];
-    if(fabs(cross(q[last].v, q[last-1].v)) < EPS) {
-      last--;
-      if(OnLeft(q[last], L[i].P)) q[last] = L[i];
-    }
-    if(first < last) p[last-1] = GetLineIntersection(q[last-1], q[last]);
-  }
-  while(first < last && !OnLeft(q[first], p[last-1])) last--;
-  if(last - first <= 1) return ans;
-  p[last] = GetLineIntersection(q[last], q[first]);
-  for(int i = first; i <= last; i++) ans.push_back(p[i]);
-  return ans;
-}

code/Geometry_test.cpp

@@ -0,0 +1,69 @@
+#include <bits/stdc++.h>
+using namespace std;
+
+#include "Geometry.cpp"
+
+bool operator==(const PT &a, const PT &b) {
+  return a.x == b.x && a.y == b.y;
+}
+
+int main() {
+  assert(RotateCCW90(PT(2,5)) == PT(-5, 2));
+  assert(RotateCW90(PT(2,5)) == PT(5, -2));
+  assert(fabs(RotateCCW(PT(2,5),M_PI/2).x - (-5)) < 0.00001);
+  assert(fabs(RotateCCW(PT(2,5),M_PI/2).y - 2) < 0.00001);
+  assert(ProjectPointLine(PT(-5,-2), PT(10,4), PT(3,7)) == PT(5, 2));
+  assert(ProjectPointSegment(PT(-5,-2), PT(10,4), PT(3,7)) == PT(5, 2));
+  assert(ProjectPointSegment(PT(7.5,3), PT(10,4), PT(3,7)) == PT(7.5, 3));
+  assert(ProjectPointSegment(PT(-5,-2), PT(2.5,1), PT(3,7)) == PT(2.5, 1));
+  assert(fabs(DistancePointPlane(4,-4,3,2,-2,5,-8) - 6.78903) < 0.0001);
+  assert(LinesParallel(PT(1,1), PT(3,5), PT(2,1), PT(4,5)) == 1);
+  assert(LinesParallel(PT(1,1), PT(3,5), PT(2,0), PT(4,5)) == 0);
+  assert(LinesParallel(PT(1,1), PT(3,5), PT(5,9), PT(7,13)) == 1);
+  assert(LinesCollinear(PT(1,1), PT(3,5), PT(2,1), PT(4,5)) == 0);
+  assert(LinesCollinear(PT(1,1), PT(3,5), PT(2,0), PT(4,5)) == 0);
+  assert(LinesCollinear(PT(1,1), PT(3,5), PT(5,9), PT(7,13)) == 1);
+  assert(SegmentsIntersect(PT(0,0), PT(2,4), PT(3,1), PT(-1,3)) == 1);
+  assert(SegmentsIntersect(PT(0,0), PT(2,4), PT(4,3), PT(0,5)) == 1);
+  assert(SegmentsIntersect(PT(0,0), PT(2,4), PT(2,-1), PT(-2,1)) == 1);
+  assert(SegmentsIntersect(PT(0,0), PT(2,4), PT(5,5), PT(1,7)) == 0);
+  assert(ComputeLineIntersection(PT(0,0), PT(2,4), PT(3,1), PT(-1,3)) == PT(1, 2));
+  assert(ComputeCircleCenter(PT(-3,4), PT(6,1), PT(4,5)) == PT(1, 1));
+  vector<PT> v;
+  v.push_back(PT(0,0)); v.push_back(PT(5,0)); v.push_back(PT(5,5)); v.push_back(PT(0,5));
+  assert(PointInPolygon(v, PT(2,2)) == 1);
+  assert(PointInPolygon(v, PT(2,0)) == 1);
+  assert(PointInPolygon(v, PT(0,2)) == 1);
+  assert(PointInPolygon(v, PT(5,2)) == 0);
+  assert(PointInPolygon(v, PT(2,5)) == 0);
+  assert(PointOnPolygon(v, PT(2,2)) == 0);
+  assert(PointOnPolygon(v, PT(2,0)) == 1);
+  assert(PointOnPolygon(v, PT(0,2)) == 1);
+  assert(PointOnPolygon(v, PT(5,2)) == 1);
+  assert(PointOnPolygon(v, PT(2,5)) == 1);
+  // expected: (1,6)
+  //           (5,4) (4,5)
+  //           blank line
+  //           (4,5) (5,4)
+  //           blank line
+  //           (4,5) (5,4)
+  // vector<PT> u = CircleLineIntersection(PT(0,6), PT(2,6), PT(1,1), 5);
+  // for (int i = 0; i < u.size(); i++) cerr << u[i] << " "; cerr << endl;
+  // u = CircleLineIntersection(PT(0,9), PT(9,0), PT(1,1), 5);
+  // for (int i = 0; i < u.size(); i++) cerr << u[i] << " "; cerr << endl;
+  // u = CircleCircleIntersection(PT(1,1), PT(10,10), 5, 5);
+  // for (int i = 0; i < u.size(); i++) cerr << u[i] << " "; cerr << endl;
+  // u = CircleCircleIntersection(PT(1,1), PT(8,8), 5, 5);
+  // for (int i = 0; i < u.size(); i++) cerr << u[i] << " "; cerr << endl;
+  // u = CircleCircleIntersection(PT(1,1), PT(4.5,4.5), 10, sqrt(2.0)/2.0);
+  // for (int i = 0; i < u.size(); i++) cerr << u[i] << " "; cerr << endl;
+  // u = CircleCircleIntersection(PT(1,1), PT(4.5,4.5), 5, sqrt(2.0)/2.0);
+  // for (int i = 0; i < u.size(); i++) cerr << u[i] << " "; cerr << endl;
+  PT pa[] = { PT(0,0), PT(5,0), PT(1,1), PT(0,5) };
+  vector<PT> p(pa, pa+4);
+  PT c = ComputeCentroid(p);
+  assert(fabs(ComputeArea(p) - 5.) < 0.00001);
+  assert(fabs(c.x - 1.166666) < 0.00001);
+  assert(fabs(c.y - 1.166666) < 0.00001);
+  return 0;
+}

code/HalfPlaneIntersection.cpp

@@ -0,0 +1,45 @@
+// INTERSECTION DE DEMI-PLAN
+// demi-plan : à gauche de la ligne
+struct Line {
+  PT P;    // point
+  PT v;   // vecteur directeur
+  double ang;
+  Line() {}
+  Line(PT P, PT v):P(P),v(v){ ang = atan2(v.y, v.x); }
+  bool operator < (const Line& L) const {
+    return ang < L.ang;
+  }
+};
+bool OnLeft(const Line& L, const PT& p) {
+  return cross(L.v, p-L.P) > 0;
+}
+PT GetLineIntersection(const Line& a, const Line& b) {
+  return ComputeLineIntersection(a.P, a.P+a.v, b.P, b.P+b.v);
+}
+// INPUT : set of half-plane
+// OUPUT : polygon (if exists)
+// COMPLEXITY n log n ?
+vector<PT> HalfplaneIntersection(vector<Line> L) {
+  int n = L.size();
+  sort(L.begin(), L.end());
+  int first, last;
+  vector<PT> p(n);
+  vector<Line> q(n);
+  vector<PT> ans;
+  q[first=last=0] = L[0];
+  for(int i = 1; i < n; i++) {
+    while(first < last && !OnLeft(L[i], p[last-1])) last--;
+    while(first < last && !OnLeft(L[i], p[first])) first++;
+    q[++last] = L[i];
+    if(fabs(cross(q[last].v, q[last-1].v)) < EPS) {
+      last--;
+      if(OnLeft(q[last], L[i].P)) q[last] = L[i];
+    }
+    if(first < last) p[last-1] = GetLineIntersection(q[last-1], q[last]);
+  }
+  while(first < last && !OnLeft(q[first], p[last-1])) last--;
+  if(last - first <= 1) return ans;
+  p[last] = GetLineIntersection(q[last], q[first]);
+  for(int i = first; i <= last; i++) ans.push_back(p[i]);
+  return ans;
+}

main.tex

@@ -108,10 +108,13 @@ Temps de cuisson : $O(n)$
 {\scriptsize\lstinputlisting{code/ConvexHull.cc}} % OK
 
 \subsection{Miscellaneous geometry}
-{\scriptsize\lstinputlisting{code/Geometry.cc}}
+{\scriptsize\lstinputlisting{code/Geometry.cpp}}
 % OK RotateCCW
 % OK ComputeSignedArea
 
+\subsection{Half-plane intersection}
+{\scriptsize\lstinputlisting{code/HalfPlaneIntersection.cpp}}
+
 \subsection{Closest Pair}
 {\scriptsize\lstinputlisting{code/ClosestPair.cpp}}