merged

Author: Chillee

content/geometry/Point.h

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

content/geometry/SegmentIntersection.h

@@ -6,8 +6,8 @@
  * Description:\\
 \begin{minipage}{75mm}
 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.
-The wrong position will be returned if P is Point<int> and the intersection point does not have integer coordinates.
+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.
 \end{minipage}
 \begin{minipage}{15mm}
@@ -16,31 +16,24 @@ Products of three coordinates are used in intermediate steps so watch out for ov
  * Status: Well tested with fuzz-test and with Kattis problem intersection.
  * Usage:
  * vector<P> inter = segInter(s1,e1,s2,e2);
- * if (inter.size()==1)
+ * if (sz(inter)==1)
  *   cout << "segments intersect at " << inter[0] << endl;
  */
 #pragma once
 
 #include "Point.h"
 #include "onSegment.h"
 
-template <class P>
-bool segInterProper(P a, P b, P c, P d, P& out) {
-	double oa = c.cross(d, a), ob = c.cross(d, b),
-	       oc = a.cross(b, c), od = a.cross(b, d);
-	if (sgn(oa) * sgn(ob) < 0 && sgn(oc) * sgn(od) < 0) {
-		out = (a * ob - b * oa) / (ob - oa);
-		return true;
-	}
-	return false;
-}
 template<class P> vector<P> segInter(P a, P b, P c, P d) {
-	P out;
-	if (segInterProper(a, b, c, d, out)) return {out};
+	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 vector<P>(all(s));
+	return {all(s)};
 }

content/geometry/onSegment.h

@@ -11,5 +11,5 @@
 #include "Point.h"
 
 template <class P> bool onSegment(P s, P e, P p) {
-	return s.cross(e, p) == 0 && (s - p).dot(e - p) <= 0;
+	return p.cross(s, e) == 0 && (s - p).dot(e - p) <= 0;
 }

fuzz-tests/geometry/SegmentIntersection.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;
+}