Skip to content

Commit b5e42ac

Browse files
IvanRenisonIvanRenison
authored
Add half planes intersection (#29)
* Use halfPlaneIntersection from el vasito * Fixes, improvements and better stress-tests * Add doubles stress test * Add test problems
1 parent db2a9c1 commit b5e42ac

File tree

12 files changed

+729
-3
lines changed

12 files changed

+729
-3
lines changed

‎content/geometry/ConvexHull.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ Points on the edge of the hull between two other points are not considered part
2020

2121
#include "Point.h"
2222

23-
typedef Point<ll> P;
23+
template<class P>
2424
vector<P> convexHull(vector<P> pts) {
2525
if (SZ(pts) <= 1) return pts;
2626
sort(ALL(pts));

‎content/geometry/HalfplaneIntersection.h‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
/**
2+
* Author: Iván Renison
3+
* Date: 2024-09-04
4+
* License: CC0
5+
* Source: notebook el vasito
6+
* Description: Computes the intersection of a set of half-planes.
7+
* Input is given as a set of planes, facing left.
8+
* The intersection must form a convex polygon or be empty.
9+
* Output is the convex polygon representing the intersection in CCW order.
10+
* The points may have duplicates and be collinear.
11+
* Time: O(n \log n)
12+
* Status: stress-tested ans problem tested
13+
*/
14+
#pragma once
15+
16+
#include "Point.h"
17+
#include "sideOf.h"
18+
#include "lineIntersection.h"
19+
20+
typedef Point<double> P;
21+
struct Line {
22+
P p, q;
23+
double a;
24+
Line() {}
25+
Line(P p, P q) : p(p), q(q), a((q - p).angle()) {}
26+
bool operator<(Line o) const { return a < o.a; }
27+
};
28+
#define L(a) a.p, a.q
29+
#define PQ(a) (a.q - a.p)
30+
31+
vector<P> halfPlaneIntersection(vector<Line> v) {
32+
sort(ALL(v));
33+
ll n = SZ(v), q = 1, h = 0;
34+
const double eps = 1e-9;
35+
vector<Line> c(n+2);
36+
#define I(j, k) lineInter(L(c[j]), L(c[k])).snd
37+
fore(i, 0, n) {
38+
while (q < h && sideOf(L(v[i]), I(h, h-1), eps) < 0) h--;
39+
while (q < h && sideOf(L(v[i]), I(q, q+1), eps) < 0) q++;
40+
c[++h] = v[i];
41+
if (q < h && abs(PQ(c[h]).cross(PQ(c[h-1]))) < eps) {
42+
if (PQ(c[h]).dot(PQ(c[h-1])) <= 0) return {};
43+
if (sideOf(L(v[i]), c[--h].p, eps) < 0) c[h] = v[i];
44+
}
45+
}
46+
while (q < h - 1 && sideOf(L(c[q]), I(h, h-1), eps) < 0) h--;
47+
while (q < h - 1 && sideOf(L(c[h]), I(q, q+1), eps) < 0) q++;
48+
if (h - q <= 1) return {};
49+
c[++h] = c[q];
50+
vector<P> s;
51+
fore(i, q, h) s.pb(I(i, i+1));
52+
return s;
53+
}

‎content/geometry/chapter.tex‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ \section{Geometric primitives}
88
\kactlimport{lineIntersection.h}
99
\kactlimport{sideOf.h}
1010
\kactlimport{OnSegment.h}
11+
\kactlimport{HalfplaneIntersection.h}
1112
\kactlimport{linearTransformation.h}
1213
\kactlimport{LineProjectionReflection.h}
1314
\kactlimport{Angle.h}

‎stress-tests/geometry/ConvexHull.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,8 @@
33
#include "../../content/geometry/ConvexHull.h"
44
#include "../utilities/bench.h"
55

6+
typedef Point<ll> P;
7+
68
namespace old {
79
pair<vi, vi> ulHull(const vector<P>& S) {
810
vi Q(SZ(S)), U, L;

‎stress-tests/geometry/FastDelaunay.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@
99
#include "../../content/geometry/circumcircle.h"
1010
#undef P
1111

12+
typedef Point<ll> P;
13+
1214
P2 top(P x) { return P2((double)x.x, (double)x.y); }
1315

1416
struct Bumpalloc {

‎stress-tests/geometry/HalfplaneIntersection.cpp‎

Lines changed: 186 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,186 @@
1+
#include "../utilities/template.h"
2+
#include "./utilities.h"
3+
4+
#include "../../content/geometry/HalfplaneIntersection.h"
5+
#include "../../content/geometry/ConvexHull.h"
6+
7+
#pragma GCC optimize ("trapv")
8+
9+
// Test against brote force with fractions
10+
namespace slow {
11+
12+
typedef __int128_t i128;
13+
14+
struct frac {
15+
i128 num, den;
16+
frac(i128 num_ = 0, i128 den_ = 1) : num(num_), den(den_) {
17+
i128 g = __gcd(num, den);
18+
num /= g;
19+
den /= g;
20+
if (den < 0) {
21+
num = -num;
22+
den = -den;
23+
}
24+
assert(den > 0);
25+
}
26+
bool operator<(frac f) const {
27+
return num * f.den < f.num * den;
28+
}
29+
bool operator==(frac f) const {
30+
return num == f.num && den == f.den;
31+
}
32+
bool operator<=(frac f) const {
33+
return num * f.den <= f.num * den;
34+
}
35+
bool operator>(frac f) const {
36+
return num * f.den > f.num * den;
37+
}
38+
frac operator+(frac o) const {
39+
return frac(num * o.den + o.num * den, den * o.den);
40+
}
41+
frac operator-(frac o) const {
42+
return frac(num * o.den - o.num * den, den * o.den);
43+
}
44+
frac operator*(frac o) const {
45+
return frac(num * o.num, den * o.den);
46+
}
47+
frac operator/(frac o) const {
48+
return frac(num * o.den, den * o.num);
49+
}
50+
};
51+
52+
typedef Point<frac> Pf;
53+
54+
vector<Pf> slow(const vector<pair<Pf,Pf>>& t) {
55+
ll n = SZ(t);
56+
vector<Pf> points;
57+
fore(i, 0, n) fore(j, 0, i) {
58+
auto [si, ei] = t[i];
59+
auto [sj, ej] = t[j];
60+
auto [x, p] = lineInter(si, ei, sj, ej);
61+
if (x == 1) {
62+
points.push_back(p);
63+
}
64+
}
65+
66+
vector<Pf> ans;
67+
for (Pf p : points) {
68+
bool valid = true;
69+
for (auto [s, e] : t) {
70+
ll side = sideOf(s, e, p);
71+
if (side == -1) {
72+
valid = false;
73+
break;
74+
}
75+
}
76+
if (valid) {
77+
ans.push_back(p);
78+
}
79+
}
80+
81+
ans = convexHull(ans);
82+
return ans;
83+
}
84+
}
85+
86+
typedef slow::Pf Pf;
87+
88+
P Pf_to_P(Pf p) {
89+
return P((double)p.x.num / p.x.den, (double)p.y.num / p.y.den);
90+
}
91+
92+
Pf randIntPt(ll lim) {
93+
return Pf{rand() % (2 * lim + 1) - lim, rand() % (2 * lim + 1) - lim};
94+
}
95+
96+
slow::frac randFrac(ll lim) {
97+
ll den = rand() % lim + 1;
98+
ll num = rand() % den;
99+
return slow::frac(num, den);
100+
}
101+
102+
Pf randDoublePt(ll lim) {
103+
Pf ans = randIntPt(lim);
104+
ans.x = ans.x + randFrac(lim / 2);
105+
ans.y = ans.y + randFrac(lim / 2);
106+
return ans;
107+
}
108+
109+
110+
const slow::frac INF = 500;
111+
void addInf(vector<pair<Pf,Pf>> &ans, slow::frac INF = INF) {
112+
slow::frac nINF = slow::frac() - INF;
113+
vector<Pf> infPts({Pf(INF, INF), Pf(nINF, INF), Pf(nINF, nINF), Pf(INF, nINF)});
114+
fore(i, 0, 4) {
115+
ans.push_back({infPts[i], infPts[(i + 1) % 4]});
116+
}
117+
}
118+
119+
void test(const vector<pair<Pf,Pf>>& t) {
120+
const double eps = 1e-11;
121+
vector<Line> t_(SZ(t));
122+
fore(i, 0, SZ(t)) {
123+
t_[i] = {Pf_to_P(t[i].first), Pf_to_P(t[i].second)};
124+
}
125+
vector<P> ans = halfPlaneIntersection(t_);
126+
assert(isConvexCCW(ans, eps));
127+
ans = convexHull(ans); // Remove colinear
128+
vector<Pf> ansf = slow::slow(t);
129+
vector<P> ans2(SZ(ansf));
130+
fore(i, 0, SZ(ansf)) {
131+
ans2[i] = Pf_to_P(ansf[i]);
132+
}
133+
assert(polygonEq(ans, ans2, eps));
134+
}
135+
136+
void testRandomInt() {
137+
ll n = rand() % 10 + 1;
138+
vector<pair<Pf,Pf>> t;
139+
fore(i, 0, n) {
140+
Pf p = randIntPt(10);
141+
Pf q = randIntPt(10);
142+
if (p == q) continue;
143+
t.push_back({p, q});
144+
}
145+
addInf(t);
146+
test(t);
147+
}
148+
149+
void testRandomDouble() {
150+
ll n = rand() % 10 + 1;
151+
vector<pair<Pf,Pf>> t;
152+
fore(i, 0, n) {
153+
Pf p = randDoublePt(10);
154+
Pf q = randDoublePt(10);
155+
if (p == q) continue;
156+
t.push_back({p, q});
157+
}
158+
addInf(t);
159+
test(t);
160+
}
161+
162+
int main() {
163+
164+
vector<vector<pair<Pf,Pf>>> handmade = {
165+
{{Pf(0, 0), Pf(5, 0)}, {Pf(5, -2), Pf(5, 2)}, {Pf(5, 2), Pf(2, 2)}, {Pf(0, 3), Pf(0, -3)}},
166+
{{Pf(0, 0), Pf(5, 0)}},
167+
{{Pf(0, 0), Pf(5, 0)}, {Pf(5, 0), Pf(0, 0)}}, // Line
168+
{{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
169+
{{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
170+
{{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
171+
{{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
172+
{{Pf(0, 0), Pf(1, 0)}, {Pf(0, 0), Pf(2, 0)}, {Pf(0, 0), Pf(3, 0)}}
173+
};
174+
175+
for (auto& t : handmade) {
176+
addInf(t);
177+
test(t);
178+
}
179+
180+
fore(_, 0, 1000) {
181+
testRandomInt();
182+
testRandomDouble();
183+
}
184+
185+
cout << "Tests passed!" << endl;
186+
}

‎stress-tests/geometry/HalfplaneIntersection2.cpp‎

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
#include "../utilities/template.h"
2+
#include "../utilities/randGeo.h"
3+
#include "./utilities.h"
4+
5+
#include "../../content/geometry/HalfplaneIntersection.h"
6+
#include "../../content/geometry/ConvexHull.h"
7+
8+
const double eps = 1e-11;
9+
10+
// Test against brote force with float128
11+
typedef Point<__float128> Pf;
12+
vector<Pf> slow(const vector<pair<Pf,Pf>>& t) {
13+
14+
ll n = SZ(t);
15+
vector<Pf> points;
16+
fore(i, 0, n) fore(j, 0, i) {
17+
auto [si, ei] = t[i];
18+
auto [sj, ej] = t[j];
19+
auto [x, p] = lineInter(si, ei, sj, ej);
20+
if (x == 1) {
21+
points.push_back(p);
22+
}
23+
}
24+
25+
vector<Pf> ans;
26+
for (Pf p : points) {
27+
bool valid = true;
28+
for (auto [s, e] : t) {
29+
ll side = sideOf(s, e, p, eps);
30+
if (side == -1) {
31+
valid = false;
32+
break;
33+
}
34+
}
35+
if (valid) {
36+
ans.push_back(p);
37+
}
38+
}
39+
40+
ans = convexHull(ans);
41+
return ans;
42+
}
43+
44+
const double INF = 1000;
45+
void addInf(vector<Line> &ans, double INF = INF) {
46+
vector<P> infPts({P(INF, INF), P(-INF, INF), P(-INF, -INF), P(INF, -INF)});
47+
fore(i, 0, 4) {
48+
ans.push_back({infPts[i], infPts[(i + 1) % 4]});
49+
}
50+
}
51+
52+
void test(const vector<Line>& t) {
53+
vector<P> ans = halfPlaneIntersection(t);
54+
assert(isConvexCCW(ans, eps));
55+
ans = convexHull(ans); // Remove colinear
56+
vector<pair<Pf,Pf>> t2(SZ(t));
57+
fore(i, 0, SZ(t)) {
58+
auto [p, q, _] = t[i];
59+
auto [px, py] = p;
60+
auto [qx, qy] = q;
61+
t2[i] = {Pf(px, py), Pf(qx, qy)};
62+
}
63+
vector<Pf> ans2_ = slow(t2);
64+
vector<P> ans2(SZ(ans2_));
65+
fore(i, 0, SZ(ans2_)) {
66+
auto [x, y] = ans2_[i];
67+
ans2[i] = P((double)x, (double)y);
68+
}
69+
assert(polygonEq(ans, ans2, eps));
70+
}
71+
72+
void testRandomInt() {
73+
ll n = rand() % 100 + 1;
74+
vector<Line> t;
75+
fore(i, 0, n) {
76+
P p = randIntPt(100);
77+
P q = randIntPt(100);
78+
if (p == q) continue;
79+
t.push_back({p, q});
80+
}
81+
addInf(t);
82+
test(t);
83+
}
84+
85+
void testRandomDouble() {
86+
ll n = rand() % 100 + 1;
87+
vector<Line> t;
88+
fore(i, 0, n) {
89+
P p = randDoublePt(100);
90+
P q = randDoublePt(100);
91+
if ((p - q).dist2() < eps) continue;
92+
t.push_back({p, q});
93+
}
94+
addInf(t);
95+
test(t);
96+
}
97+
98+
int main() {
99+
100+
vector<vector<Line>> handmade = {
101+
{{P(0, 0), P(5, 0)}, {P(5, -2), P(5, 2)}, {P(5, 2), P(2, 2)}, {P(0, 3), P(0, -3)}},
102+
{{P(0, 0), P(5, 0)}},
103+
{{P(0, 0), P(5, 0)}, {P(5, 0), P(0, 0)}}, // Line
104+
{{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
105+
{{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
106+
{{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
107+
{{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
108+
};
109+
110+
for (auto& t : handmade) {
111+
addInf(t);
112+
test(t);
113+
}
114+
115+
fore(_, 0, 1000) {
116+
testRandomInt();
117+
testRandomDouble();
118+
}
119+
120+
cout << "Tests passed!" << endl;
121+
}

0 commit comments

Comments
 (0)