Updated onSegment.h and both segment intersections (#78)

Author: Chillee

content/geometry/Point.h

@@ -9,6 +9,7 @@
  */
 #pragma once
 
+template <class T> int sgn(T x) { return (x > 0) - (x < 0); }
 template<class T>
 struct Point {
 	typedef Point P;

content/geometry/SegmentIntersection.h

@@ -1,53 +1,39 @@
 /**
- * Author: Ulf Lundstrom
- * Date: 2009-03-21
+ * Author: Victor Lecomte, chilli
+ * Date: 2019-04-27
  * License: CC0
- * Source:
+ * Source: https://vlecomte.github.io/cp-geo.pdf
  * Description:\\
 \begin{minipage}{75mm}
-If a unique intersetion point between the line segments going from s1 to e1 and from s2 to e2 exists r1 is set to this point and 1 is returned.
-If no intersection point exists 0 is returned and if infinitely many exists 2 is returned and r1 and r2 are set to the two ends of the common line.
-The wrong position will be returned if P is Point<int> and the intersection point does not have integer coordinates.
+If a unique intersection point between the line segments going from s1 to e1 and from s2 to e2 exists then it is returned.
+If no intersection point exists an empty vector is returned. If infinitely many exist a vector with 2 elements is returned, containing the endpoints of the common line segment.
+The wrong position will be returned if P is Point<ll> and the intersection point does not have integer coordinates.
 Products of three coordinates are used in intermediate steps so watch out for overflow if using int or long long.
-Use segmentIntersectionQ to get just a true/false answer.
 \end{minipage}
 \begin{minipage}{15mm}
 \includegraphics[width=\textwidth]{content/geometry/SegmentIntersection}
 \end{minipage}
- * Status: Well tested with unitTest and with Kattis problem intersection.
+ * Status: Well tested with fuzz-test and with Kattis problem intersection.
  * Usage:
- * Point<double> intersection, dummy;
- * if (segmentIntersection(s1,e1,s2,e2,intersection,dummy)==1)
- *   cout << "segments intersect at " << intersection << endl;
+ * vector<P> inter = segInter(s1,e1,s2,e2);
+ * if (sz(inter)==1)
+ *   cout << "segments intersect at " << inter[0] << endl;
  */
 #pragma once
 
 #include "Point.h"
+#include "onSegment.h"
 
-template<class P>
-int segmentIntersection(const P& s1, const P& e1,
-		const P& s2, const P& e2, P& r1, P& r2) {
-	if (e1==s1) {
-		if (e2==s2) {
-			if (e1==e2) { r1 = e1; return 1; } //all equal
-			else return 0; //different point segments
-		} else return segmentIntersection(s2,e2,s1,e1,r1,r2);//swap
-	}
-	//segment directions and separation
-	P v1 = e1-s1, v2 = e2-s2, d = s2-s1;
-	auto a = v1.cross(v2), a1 = v1.cross(d), a2 = v2.cross(d);
-	if (a == 0) { //if parallel
-		auto b1=s1.dot(v1), c1=e1.dot(v1),
-		     b2=s2.dot(v1), c2=e2.dot(v1);
-		if (a1 || a2 || max(b1,min(b2,c2))>min(c1,max(b2,c2)))
-			return 0;
-		r1 = min(b2,c2)<b1 ? s1 : (b2<c2 ? s2 : e2);
-		r2 = max(b2,c2)>c1 ? e1 : (b2>c2 ? s2 : e2);
-		return 2-(r1==r2);
-	}
-	if (a < 0) { a = -a; a1 = -a1; a2 = -a2; }
-	if (0<a1 || a<-a1 || 0<a2 || a<-a2)
-		return 0;
-	r1 = s1-v1*a2/a;
-	return 1;
+template<class P> vector<P> segInter(P a, P b, P c, P d) {
+	auto oa = c.cross(d, a), ob = c.cross(d, b),
+	     oc = a.cross(b, c), od = a.cross(b, d);
+	// Checks if intersection is single non-endpoint point.
+	if (sgn(oa) * sgn(ob) < 0 && sgn(oc) * sgn(od) < 0)
+		return {(a * ob - b * oa) / (ob - oa)};
+	set<P> s;
+	if (onSegment(c, d, a)) s.insert(a);
+	if (onSegment(c, d, b)) s.insert(b);
+	if (onSegment(a, b, c)) s.insert(c);
+	if (onSegment(a, b, d)) s.insert(d);
+	return {all(s)};
 }

content/geometry/SegmentIntersectionQ.h

@@ -6,15 +6,13 @@ \section{Geometric primitives}
 	\kactlimport{SegmentDistance.h}
 	\columnbreak
 	\kactlimport{SegmentIntersection.h}
-	\kactlimport{SegmentIntersectionQ.h}
-	\columnbreak
 	\kactlimport{lineIntersection.h}
 	\kactlimport{sideOf.h}
 	\kactlimport{onSegment.h}
 	\kactlimport{linearTransformation.h}
 	\kactlimport{Angle.h}
 
-\section{Circles}	
+\section{Circles}
 	\kactlimport{CircleIntersection.h}
 	\kactlimport{circleTangents.h}
 	\kactlimport{circumcircle.h}

content/geometry/chapter.tex

@@ -1,17 +1,15 @@
 /**
- * Author: Ulf Lundstrom
- * Date: 2009-04-09
+ * Author: Victor Lecomte, chilli
+ * Date: 2019-04-26
  * License: CC0
- * Source: Basic geometry
- * Description: Returns true iff p lies on the line segment from s to e. Intended for use with e.g. Point<long long> where overflow is an issue. Use (segDist(s,e,p)<=epsilon) instead when using Point<double>.
+ * Source: https://vlecomte.github.io/cp-geo.pdf
+ * Description: Returns true iff p lies on the line segment from s to e. Intended for use with e.g. Point<long long> where overflow is an issue.  Use (segDist(s,e,p)<=epsilon) instead when using Point<double>.
  * Status:
  */
 #pragma once
 
 #include "Point.h"
 
-template<class P>
-bool onSegment(const P& s, const P& e, const P& p) {
-	P ds = p-s, de = p-e;
-	return ds.cross(de) == 0 && ds.dot(de) <= 0;
+template <class P> bool onSegment(P s, P e, P p) {
+	return p.cross(s, e) == 0 && (s - p).dot(e - p) <= 0;
 }

content/geometry/onSegment.h

@@ -27,9 +27,9 @@ int main() {
 			P a = ps[i], b = ps[(i+1)%N];
 			P c = ps[j], d = ps[(j+1)%N];
 			P r1, r2;
-			int r = segmentIntersection(a, b, c, d, r1, r2);
-			if (r == 2) goto fail;
-			if (r == 1) {
+			auto r = segInter(a, b, c, d);
+			if (sz(r) == 2) goto fail;
+			if (sz(r) == 1) {
 				if (i+1 == j || (j+1) % N == i) ;
 				else goto fail;
 			}

fuzz-tests/geometry/PolygonCut.cpp

@@ -0,0 +1,71 @@
+#include <bits/stdc++.h>
+using namespace std;
+
+#define rep(i, a, b) for(int i = a; i < int(b); ++i)
+#define trav(a, v) for(auto& a : v)
+#define all(x) x.begin(), x.end()
+#define sz(x) (int)(x).size()
+
+typedef long long ll;
+typedef pair<int, int> pii;
+typedef vector<int> vi;
+
+
+#include "../content/geometry/SegmentIntersection.h"
+namespace oldImpl {
+template<class P>
+int segmentIntersection(const P& s1, const P& e1,
+        const P& s2, const P& e2, P& r1, P& r2) {
+    if (e1==s1) {
+        if (e2==s2) {
+            if (e1==e2) { r1 = e1; return 1; } //all equal
+            else return 0; //different point segments
+        } else return segmentIntersection(s2,e2,s1,e1,r1,r2);//swap
+    }
+    //segment directions and separation
+    P v1 = e1-s1, v2 = e2-s2, d = s2-s1;
+    auto a = v1.cross(v2), a1 = v1.cross(d), a2 = v2.cross(d);
+    if (a == 0) { //if parallel
+        auto b1=s1.dot(v1), c1=e1.dot(v1),
+             b2=s2.dot(v1), c2=e2.dot(v1);
+        if (a1 || a2 || max(b1,min(b2,c2))>min(c1,max(b2,c2)))
+            return 0;
+        r1 = min(b2,c2)<b1 ? s1 : (b2<c2 ? s2 : e2);
+        r2 = max(b2,c2)>c1 ? e1 : (b2>c2 ? s2 : e2);
+        return 2-(r1==r2);
+    }
+    if (a < 0) { a = -a; a1 = -a1; a2 = -a2; }
+    if (0<a1 || a<-a1 || 0<a2 || a<-a2)
+        return 0;
+    r1 = s1-v1*a2/a;
+    return 1;
+}
+}
+typedef Point<double> P;
+ostream &operator<<(ostream &os, P p) { return cout << "(" << p.x << "," << p.y << ")"; }
+bool eq(P a, P b) {
+    return (a-b).dist()<1e-8;
+}
+int main() {
+    rep(t,0,1000000) {
+        const int GRID=6;
+        P a(rand()%GRID, rand()%GRID), b(rand()%GRID, rand()%GRID), c(rand()%GRID, rand()%GRID), d(rand()%GRID, rand()%GRID);
+        P tmp1, tmp2;
+        auto res = oldImpl::segmentIntersection(a,b,c,d, tmp1, tmp2);
+        auto res2 = segInter(a,b,c,d);
+        if (res != sz(res2)) {
+            cout<<a<<' '<<b<<' '<<c<<' '<<d<<endl;
+            cout<<"old: "<<res<<" new: "<<sz(res2)<<endl;
+        }
+        assert(res==sz(res2));
+        if (res==1) {
+            assert(eq(*res2.begin(), tmp1));
+        } else if (res==2) {
+            vector<P> a(res2.begin(), res2.end());
+            vector<P> b({tmp1, tmp2});
+            sort(all(b));
+            assert(eq(a[0], b[0]) && eq(a[1],b[1]));
+        }
+    }
+    cout<<"Tests passed!"<<endl;
+}