Add half planes intersection

content/geometry/ConvexHull.h

@@ -20,7 +20,7 @@ Points on the edge of the hull between two other points are not considered part
 
 #include "Point.h"
 
-typedef Point<ll> P;
+template<class P>
 vector<P> convexHull(vector<P> pts) {
 	if (SZ(pts) <= 1) return pts;
 	sort(ALL(pts));

content/geometry/HalfplaneIntersection.h

@@ -0,0 +1,49 @@
+/**
+ * Author: chilli and Iván Renison
+ * Date: 2019-05-05
+ * License: CC0
+ * Source: https://github.com/zimpha/algorithmic-library/blob/master/computational-geometry/polygon.cc
+ * Description: Computes the intersection of a set of half-planes. Input is given as a set of planes, facing left.
+ * The intersection must form a convex polygon or be empty. Output is the convex polygon representing the intersection.
+ * The points may have duplicates and be collinear. Will not fail catastrophically if `eps > sqrt(2)(line intersection error)`.
+ * Likely to work for more ranges if 3 half planes are guaranteed to never intersect at the same point.
+ * Time: O(n \log n)
+ * Status: stress-tested
+ */
+#pragma once
+
+#include "Point.h"
+#include "sideOf.h"
+#include "lineIntersection.h"
+
+typedef Point<double> P;
+typedef pair<P, P> Line;
+#define L(a) a.fst, a.snd
+
+ll angDiff(Line a, Line b) {
+	return sgn((a.snd-a.fst).angle() - (b.snd-b.fst).angle());
+}
+vector<P> halfPlaneIntersection(vector<Line> v) {
+	const double eps = sqrt(2) * 1e-8;
+	sort(ALL(v), [&](Line a, Line b) {
+		ll s = angDiff(a, b);
+		return (s ? s : sideOf(L(a), b.fst)) < 0;
+	});
+	ll ah = 0, at = 0, n = SZ(v);
+	vector<Line> deq(n + 1);
+	vector<P> ans(n);
+	deq[0] = v[0];
+	fore(i, 1, n + 1) {
+		if (i == n) v.pb(deq[ah]);
+		if (angDiff(v[i], v[i - 1])) {
+			while (ah < at && sideOf(L(v[i]), ans[at-1], eps) < 0)
+				at--;
+			while (i < n && ah < at && sideOf(L(v[i]),ans[ah],eps)<0)
+				ah++;
+			auto [r, p] = lineInter(L(v[i]), L(deq[at]));
+			if (r == 1) ans[at++] = p, deq[at] = v[i];
+		}
+	}
+	if (at - ah < 3) return {};
+	return {ans.begin() + ah, ans.begin() + at};
+}

content/geometry/chapter.tex

@@ -8,6 +8,7 @@ \section{Geometric primitives}
 	\kactlimport{lineIntersection.h}
 	\kactlimport{sideOf.h}
 	\kactlimport{OnSegment.h}
+	\kactlimport{HalfplaneIntersection.h}
 	\kactlimport{linearTransformation.h}
 	\kactlimport{LineProjectionReflection.h}
 	\kactlimport{Angle.h}

stress-tests/geometry/ConvexHull.cpp

@@ -3,6 +3,8 @@
 #include "../../content/geometry/ConvexHull.h"
 #include "../utilities/bench.h"
 
+typedef Point<ll> P;
+
 namespace old {
 pair<vi, vi> ulHull(const vector<P>& S) {
 	vi Q(SZ(S)), U, L;

stress-tests/geometry/FastDelaunay.cpp

@@ -9,6 +9,8 @@
 #include "../../content/geometry/circumcircle.h"
 #undef P
 
+typedef Point<ll> P;
+
 P2 top(P x) { return P2((double)x.x, (double)x.y); }
 
 struct Bumpalloc {

stress-tests/geometry/HalfplaneIntersection.cpp

@@ -0,0 +1,184 @@
+#include "../utilities/template.h"
+#include "./utilities.h"
+
+#include "../../content/geometry/HalfplaneIntersection.h"
+#include "../../content/geometry/ConvexHull.h"
+
+#pragma GCC optimize ("trapv")
+
+// Test against brote force with fractions
+namespace slow {
+
+	typedef __int128_t i128;
+
+	struct frac {
+		i128 num, den;
+		frac(i128 num_ = 0, i128 den_ = 1) : num(num_), den(den_) {
+			i128 g = __gcd(num, den);
+			num /= g;
+			den /= g;
+			if (den < 0) {
+				num = -num;
+				den = -den;
+			}
+			assert(den > 0);
+		}
+		bool operator<(frac f) const {
+			return num * f.den < f.num * den;
+		}
+		bool operator==(frac f) const {
+			return num == f.num && den == f.den;
+		}
+		bool operator<=(frac f) const {
+			return num * f.den <= f.num * den;
+		}
+		bool operator>(frac f) const {
+			return num * f.den > f.num * den;
+		}
+		frac operator+(frac o) const {
+			return frac(num * o.den + o.num * den, den * o.den);
+		}
+		frac operator-(frac o) const {
+			return frac(num * o.den - o.num * den, den * o.den);
+		}
+		frac operator*(frac o) const {
+			return frac(num * o.num, den * o.den);
+		}
+		frac operator/(frac o) const {
+			return frac(num * o.den, den * o.num);
+		}
+	};
+
+	typedef Point<frac> Pf;
+
+	vector<Pf> slow(const vector<pair<Pf,Pf>>& t) {
+		ll n = SZ(t);
+		vector<Pf> points;
+		fore(i, 0, n) fore(j, 0, i) {
+			auto [si, ei] = t[i];
+			auto [sj, ej] = t[j];
+			auto [x, p] = lineInter(si, ei, sj, ej);
+			if (x == 1) {
+				points.push_back(p);
+			}
+		}
+
+		vector<Pf> ans;
+		for (Pf p : points) {
+			bool valid = true;
+			for (auto [s, e] : t) {
+				ll side = sideOf(s, e, p);
+				if (side == -1) {
+					valid = false;
+					break;
+				}
+			}
+			if (valid) {
+				ans.push_back(p);
+			}
+		}
+
+		ans = convexHull(ans);
+		return ans;
+	}
+}
+
+typedef slow::Pf Pf;
+
+P Pf_to_P(Pf p) {
+	return P((double)p.x.num / p.x.den, (double)p.y.num / p.y.den);
+}
+
+Pf randIntPt(ll lim) {
+	return Pf{rand() % (2 * lim + 1) - lim, rand() % (2 * lim + 1) - lim};
+}
+
+slow::frac randFrac(ll lim) {
+	ll den = rand() % lim + 1;
+	ll num = rand() % den;
+	return slow::frac(num, den);
+}
+
+Pf randDoublePt(ll lim) {
+	Pf ans = randIntPt(lim);
+	ans.x = ans.x + randFrac(lim / 2);
+	ans.y = ans.y + randFrac(lim / 2);
+	return ans;
+}
+
+
+const slow::frac INF = 500;
+void addInf(vector<pair<Pf,Pf>> &ans, slow::frac INF = INF) {
+	slow::frac nINF = slow::frac() - INF;
+	vector<Pf> infPts({Pf(INF, INF), Pf(nINF, INF), Pf(nINF, nINF), Pf(INF, nINF)});
+	fore(i, 0, 4) {
+		ans.push_back({infPts[i], infPts[(i + 1) % 4]});
+	}
+}
+
+void test(const vector<pair<Pf,Pf>>& t) {
+	const double eps = 1e-11;
+	vector<Line> t_(SZ(t));
+	fore(i, 0, SZ(t)) {
+		t_[i] = {Pf_to_P(t[i].first), Pf_to_P(t[i].second)};
+	}
+	vector<P> ans = halfPlaneIntersection(t_);
+	assert(isConvexCCW(ans, eps));
+	ans = convexHull(ans); // Remove colinear
+	vector<Pf> ansf = slow::slow(t);
+	vector<P> ans2(SZ(ansf));
+	fore(i, 0, SZ(ansf)) {
+		ans2[i] = Pf_to_P(ansf[i]);
+	}
+	assert(polygonEq(ans, ans2, eps));
+}
+
+void testRandomInt() {
+	ll n = rand() % 10 + 1;
+	vector<pair<Pf,Pf>> t;
+	fore(i, 0, n) {
+		Pf p = randIntPt(10);
+		Pf q = randIntPt(10);
+		if (p == q) continue;
+		t.push_back({p, q});
+	}
+	addInf(t);
+	test(t);
+}
+
+void testRandomDouble() {
+	ll n = rand() % 10 + 1;
+	vector<pair<Pf,Pf>> t;
+	fore(i, 0, n) {
+		Pf p = randDoublePt(10);
+		Pf q = randDoublePt(10);
+		if (p == q) continue;
+		t.push_back({p, q});
+	}
+	addInf(t);
+	test(t);
+}
+
+int main() {
+
+	vector<vector<pair<Pf,Pf>>> handmade = {
+		{{Pf(0, 0), Pf(5, 0)}, {Pf(5, -2), Pf(5, 2)}, {Pf(5, 2), Pf(2, 2)}, {Pf(0, 3), Pf(0, -3)}},
+		{{Pf(0, 0), Pf(5, 0)}},
+		{{Pf(0, 0), Pf(5, 0)}, {Pf(5, 0), Pf(0, 0)}}, // Line
+		{{Pf(0, 0), Pf(5, 0)}, {Pf(5, 0), Pf(0, 0)}, {Pf(0, 0), Pf(0, 5)}, {Pf(0, 5), Pf(0, 0)}}, // Point
+		{{Pf(0, 0), Pf(5, 0)}, {Pf(5, 0), Pf(0, 0)}, {Pf(0, 0), Pf(0, 5)}, {Pf(0, 5), Pf(0, 0)}, {Pf(0, 2), Pf(5, 2)}}, // Empty
+		{{Pf(0, 0), Pf(1, 0)}, {Pf(0, 0), Pf(2, 0)}, {Pf(0, 0), Pf(3, 0)}}
+	};
+
+	for (auto& t : handmade) {
+		addInf(t);
+		test(t);
+	}
+
+	fore(_, 0, 1000) {
+		testRandomInt();
+		testRandomDouble();
+	}
+
+	cout << "Tests passed!" << endl;
+}

stress-tests/geometry/HalfplaneIntersection2.cpp

@@ -0,0 +1,119 @@
+#include "../utilities/template.h"
+#include "../utilities/randGeo.h"
+#include "./utilities.h"
+
+#include "../../content/geometry/HalfplaneIntersection.h"
+#include "../../content/geometry/ConvexHull.h"
+
+const double eps = 1e-11;
+
+// Test against brote force with float128
+typedef Point<__float128> Pf;
+vector<Pf> slow(const vector<pair<Pf,Pf>>& t) {
+
+	ll n = SZ(t);
+	vector<Pf> points;
+	fore(i, 0, n) fore(j, 0, i) {
+		auto [si, ei] = t[i];
+		auto [sj, ej] = t[j];
+		auto [x, p] = lineInter(si, ei, sj, ej);
+		if (x == 1) {
+			points.push_back(p);
+		}
+	}
+
+	vector<Pf> ans;
+	for (Pf p : points) {
+		bool valid = true;
+		for (auto [s, e] : t) {
+			ll side = sideOf(s, e, p, eps);
+			if (side == -1) {
+				valid = false;
+				break;
+			}
+		}
+		if (valid) {
+			ans.push_back(p);
+		}
+	}
+
+	ans = convexHull(ans);
+	return ans;
+}
+
+const double INF = 1000;
+void addInf(vector<Line> &ans, double INF = INF) {
+	vector<P> infPts({P(INF, INF), P(-INF, INF), P(-INF, -INF), P(INF, -INF)});
+	fore(i, 0, 4) {
+		ans.push_back({infPts[i], infPts[(i + 1) % 4]});
+	}
+}
+
+void test(const vector<Line>& t) {
+	vector<P> ans = halfPlaneIntersection(t);
+	assert(isConvexCCW(ans, eps));
+	ans = convexHull(ans); // Remove colinear
+	vector<pair<Pf,Pf>> t2(SZ(t));
+	fore(i, 0, SZ(t)) {
+		auto [p, q] = t[i];
+		auto [px, py] = p;
+		auto [qx, qy] = q;
+		t2[i] = {Pf(px, py), Pf(qx, qy)};
+	}
+	vector<Pf> ans2_ = slow(t2);
+	vector<P> ans2(SZ(ans2_));
+	fore(i, 0, SZ(ans2_)) {
+		auto [x, y] = ans2_[i];
+		ans2[i] = P((double)x, (double)y);
+	}
+	assert(polygonEq(ans, ans2, eps));
+}
+
+void testRandomInt() {
+	ll n = rand() % 100 + 1;
+	vector<Line> t;
+	fore(i, 0, n) {
+		P p = randIntPt(100);
+		P q = randIntPt(100);
+		if (p == q) continue;
+		t.push_back({p, q});
+	}
+	addInf(t);
+	test(t);
+}
+
+void testRandomDouble() {
+	ll n = rand() % 100 + 1;
+	vector<Line> t;
+	fore(i, 0, n) {
+		P p = randDoublePt(100);
+		P q = randDoublePt(100);
+		if ((p - q).dist2() < eps) continue;
+		t.push_back({p, q});
+	}
+	addInf(t);
+	test(t);
+}
+
+int main() {
+
+	vector<vector<Line>> handmade = {
+		{{P(0, 0), P(5, 0)}, {P(5, -2), P(5, 2)}, {P(5, 2), P(2, 2)}, {P(0, 3), P(0, -3)}},
+		{{P(0, 0), P(5, 0)}},
+		{{P(0, 0), P(5, 0)}, {P(5, 0), P(0, 0)}}, // Line
+		{{P(0, 0), P(5, 0)}, {P(5, 0), P(0, 0)}, {P(0, 0), P(0, 5)}, {P(0, 5), P(0, 0)}}, // Point
+		{{P(0, 0), P(5, 0)}, {P(5, 0), P(0, 0)}, {P(0, 0), P(0, 5)}, {P(0, 5), P(0, 0)}, {P(0, 2), P(5, 2)}} // Empty
+	};
+
+	for (auto& t : handmade) {
+		addInf(t);
+		test(t);
+	}
+
+	fore(_, 0, 1000) {
+		testRandomInt();
+		testRandomDouble();
+	}
+
+	cout << "Tests passed!" << endl;
+}

stress-tests/geometry/LineHullIntersection.cpp

@@ -32,7 +32,9 @@ struct Point<double> {
 		return P(x*cos(a)-y*sin(a),x*sin(a)+y*cos(a)); }
 };
 
+
 #include "../../content/geometry/ConvexHull.h"
+typedef Point<ll> P;
 #include "../../content/geometry/LineHullIntersection.h"
 
 ll segmentIntersection(const P& s1, const P& e1,