Add half planes intersection (#29)

* Use halfPlaneIntersection from el vasito

* Fixes, improvements and better stress-tests

* Add doubles stress test

* Add test problems

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,53 @@
+/**
+ * Author: Iván Renison
+ * Date: 2024-09-04
+ * License: CC0
+ * Source: notebook el vasito
+ * 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 in CCW order.
+ * The points may have duplicates and be collinear.
+ * Time: O(n \log n)
+ * Status: stress-tested ans problem tested
+ */
+#pragma once
+
+#include "Point.h"
+#include "sideOf.h"
+#include "lineIntersection.h"
+
+typedef Point<double> P;
+struct Line {
+	P p, q;
+	double a;
+	Line() {}
+	Line(P p, P q) : p(p), q(q), a((q - p).angle()) {}
+	bool operator<(Line o) const { return a < o.a; }
+};
+#define L(a) a.p, a.q
+#define PQ(a) (a.q - a.p)
+
+vector<P> halfPlaneIntersection(vector<Line> v) {
+	sort(ALL(v));
+	ll n = SZ(v), q = 1, h = 0;
+	const double eps = 1e-9;
+	vector<Line> c(n+2);
+	#define I(j, k) lineInter(L(c[j]), L(c[k])).snd
+	fore(i, 0, n) {
+		while (q < h && sideOf(L(v[i]), I(h, h-1), eps) < 0) h--;
+		while (q < h && sideOf(L(v[i]), I(q, q+1), eps) < 0) q++;
+		c[++h] = v[i];
+		if (q < h && abs(PQ(c[h]).cross(PQ(c[h-1]))) < eps) {
+			if (PQ(c[h]).dot(PQ(c[h-1])) <= 0) return {};
+			if (sideOf(L(v[i]), c[--h].p, eps) < 0) c[h] = v[i];
+		}
+	}
+	while (q < h - 1 && sideOf(L(c[q]), I(h, h-1), eps) < 0) h--;
+	while (q < h - 1 && sideOf(L(c[h]), I(q, q+1), eps) < 0) q++;
+	if (h - q <= 1) return {};
+	c[++h] = c[q];
+	vector<P> s;
+	fore(i, q, h) s.pb(I(i, i+1));
+	return s;
+}

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,186 @@
+#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(5, 0)}, {Pf(5, 0), Pf(5, 5)}, {Pf(5, 5), Pf(0, 5)}, {Pf(0, 5), Pf(0, 0)}, {Pf(1, 5), Pf(1, 0)}}, // Parallel lines
+		{{Pf(0, 0), Pf(5, 0)}, {Pf(5, 0), Pf(5, 5)}, {Pf(5, 5), Pf(0, 5)}, {Pf(1, 5), Pf(1, 0)}, {Pf(0, 5), Pf(0, 0)}}, // Parallel lines
+		{{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,121 @@
+#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
+		{{P(0, 0), P(5, 0)}, {P(5, 0), P(5, 5)}, {P(5, 5), P(0, 5)}, {P(0, 5), P(0, 0)}, {P(1, 5), P(1, 0)}}, // Parallel lines
+		{{P(0, 0), P(5, 0)}, {P(5, 0), P(5, 5)}, {P(5, 5), P(0, 5)}, {P(1, 5), P(1, 0)}, {P(0, 5), P(0, 0)}} // Parallel lines
+	};
+
+	for (auto& t : handmade) {
+		addInf(t);
+		test(t);
+	}
+
+	fore(_, 0, 1000) {
+		testRandomInt();
+		testRandomDouble();
+	}
+
+	cout << "Tests passed!" << endl;
+}