Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Appearance settings

Commit 61036aa

Browse filesBrowse files
bqi343bqi343
committed
change default halfplane isect
1 parent 26d1428 commit 61036aa
Copy full SHA for 61036aa

File tree

Expand file treeCollapse file tree

5 files changed

+156
-148
lines changed
Filter options
Expand file treeCollapse file tree

5 files changed

+156
-148
lines changed

‎Implementations/content/contest/TemplateMinimal.cpp

Copy file name to clipboard
+7-2Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,11 @@
1-
#include<bits/stdc++.h>
1+
#include <bits/stdc++.h>
22
using namespace std;
33

4+
template<class T> using V = vector<T>;
5+
#define all(x) begin(x), end(x)
6+
47
int main() {
5-
8+
ios::sync_with_stdio(false);
9+
cin.tie(nullptr);
10+
611
}

‎Implementations/content/data-structures/1D Range Queries (9.2)/LazySeg (15.2).h

Copy file name to clipboardExpand all lines: Implementations/content/data-structures/1D Range Queries (9.2)/LazySeg (15.2).h
+1-1Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66

77
tcT, int SZ> struct LazySeg {
88
static_assert(pct(SZ) == 1); // SZ must be power of 2
9-
const T ID = 0; T cmb(T a, T b) { return a+b; }
9+
const T ID{}; T cmb(T a, T b) { return a+b; }
1010
T seg[2*SZ], lazy[2*SZ];
1111
LazySeg() { F0R(i,2*SZ) seg[i] = lazy[i] = ID; }
1212
void push(int ind, int L, int R) { /// modify values for current node

‎Implementations/content/geometry (13)/Polygons/HalfPlaneIsect.h

Copy file name to clipboard
+54-83Lines changed: 54 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -1,94 +1,65 @@
11
/**
2-
* Description: half-plane intersection area
2+
* Description: Returns vertices of half-plane intersection.
3+
* A half-plane is the area to the left of a ray, which is defined
4+
* by a point \texttt{p} and a direction \texttt{dp}.
5+
* Area of intersection should be sufficiently precise when all inputs
6+
* are integers with magnitude $\le 10^5$. Probably works with floating
7+
* point too (but not well tested). Assumes intersection is bounded
8+
* (easiest way to ensure this is to uncomment the code below).
39
* Time: O(N\log N)
410
* Source: Own
5-
* http://www.cs.umd.edu/class/spring2020/cmsc754/Lects/lect06-duality.pdf might be of interest
6-
* I would be interested if someone has a simpler implementation ...
7-
* Verification:
11+
* Verification:
812
* https://open.kattis.com/problems/bigbrother
913
* https://icpc.kattis.com/problems/domes
10-
* (planes through two points with integer coordinates <= 10^7)
14+
* https://codeforces.com/gym/100962 B
1115
*/
1216

13-
#include "Point.h"
17+
#include "AngleCmp.h"
1418

15-
using Half = AR<T,3>; // half-plane
16-
using vH = V<Half>;
17-
P hp_point(const Half& h) { return {h[0],h[1]}; } // direction of half-plane
18-
P isect(const Half& h0, const Half& h1) { // Cramer's rule to intersect half-planes
19-
AR<T,3> vals;
20-
FOR(i,-1,2) {
21-
int x = (i == 0 ? 2 : 0), y = (i == 1 ? 2 : 1);
22-
vals[1+i] = h0[x]*h1[y]-h0[y]*h1[x];
23-
}
24-
assert(vals[0] != 0); return {vals[1]/vals[0],vals[2]/vals[0]};
25-
}
26-
T eval(const Half& h, T x) { assert(h[1] != 0); // evaluate half-plane at x-coordinate
27-
return (h[2]-h[0]*x)/h[1]; }
28-
T x_isect(const Half& h0, const Half& h1) { return isect(h0,h1).f; } // x-coordinate of intersection
29-
30-
vH construct_lower(P x, vH planes) { // similar to convex hull (by duality)
31-
sort(all(planes),[](const Half& a, const Half& b) {
32-
return cross(hp_point(a),hp_point(b)) > 0; });
33-
vH res{{1,0,x.f}}; // >= x.f
34-
planes.pb({-1,0,-x.s}); // <= x.s
35-
auto lst_x = [&](Half a, Half b) {
36-
if (cross(hp_point(a),hp_point(b)) == 0) // parallel half-planes, remove lower one
37-
return a[2]/a[1] <= b[2]/b[1] ? x.f : x.s;
38-
return x_isect(a,b);
39-
};
40-
each(t,planes) {
41-
while (sz(res) > 1 && lst_x(res.bk,t) <= lst_x(end(res)[-2],res.bk))
42-
res.pop_back();
43-
res.pb(t);
44-
}
45-
return res;
46-
}
19+
struct Ray {
20+
P p, dp; // origin, direction
21+
P isect(const Ray& L) const {
22+
return p+dp*(cross(L.dp,L.p-p)/cross(L.dp,dp)); }
23+
bool operator<(const Ray& L) const {
24+
return angleCmp(dp,L.dp); }
25+
};
4726

48-
T isect_area(vH planes) {
49-
const T BOUND = 1e9; P x{-BOUND,BOUND};
50-
planes.pb({0,1,-BOUND}); // y >= -BOUND
51-
planes.pb({0,-1,-BOUND}); // -y >= -BOUND
52-
vH upper, lower;
53-
each(t,planes) {
54-
if (t[1] == 0) { // vertical line
55-
T quo = t[2]/t[0];
56-
if (t[0] > 0) ckmax(x.f,quo);
57-
else ckmin(x.s,quo); // -x >=
58-
} else if (t[1] > 0) lower.pb(t);
59-
else upper.pb(t);
27+
vP halfPlaneIsect(V<Ray> rays_orig) {
28+
// int DX = 1e9, DY = 1e9; // bound input by rectangle [0,DX] x [0,DY]
29+
// rays_orig.pb({P{0,0},P{1,0}});
30+
// rays_orig.pb({P{DX,0},P{0,1}});
31+
// rays_orig.pb({P{DX,DY},P{-1,0}});
32+
// rays_orig.pb({P{0,DY},P{0,-1}});
33+
sor(rays_orig); // sort rays by angle
34+
V<Ray> rays; // without parallel rays
35+
each(t,rays_orig) {
36+
if (!sz(rays) || rays.bk < t) { rays.pb(t); continue; }
37+
// WARNING: use cross(rays.bk,t) > EPS instead of rays.bk < t if working with floating point dp
38+
if (cross(t.dp,t.p-rays.bk.p) > 0) rays.bk = t; // last two rays are parallel, remove one
6039
}
61-
if (x.f >= x.s) return 0;
62-
lower = construct_lower(x,lower);
63-
64-
each(t,upper) t[0] *= -1, t[1] *= -1;
65-
upper = construct_lower({-x.s,-x.f},upper);
66-
each(t,upper) t[0] *= -1, t[1] *= -1;
67-
reverse(all(upper));
68-
int iu = 1, il = 1;
69-
T lst = x.f, lst_dif = eval(upper[1],lst)-eval(lower[1],lst);
70-
T ans = 0;
71-
while (iu < sz(upper)-1 && il < sz(lower)-1) { // sweep vertical line through lower and upper hulls
72-
T nex_upper = x_isect(upper[iu],upper[iu+1]);
73-
T nex_lower = x_isect(lower[il],lower[il+1]);
74-
T nex = min(nex_upper,nex_lower);
75-
T nex_dif = eval(upper[iu],nex)-eval(lower[il],nex);
76-
auto avg_val = [](T a, T b) -> T {
77-
if (a > b) swap(a,b);
78-
if (b <= 0) return 0;
79-
if (a >= 0) return (a+b)/2;
80-
return b/(b-a)*b/2;
81-
};
82-
ans += (nex-lst)*avg_val(lst_dif,nex_dif);
83-
assert(x.f <= nex && nex <= x.s);
84-
lst = nex, lst_dif = nex_dif;
85-
iu += lst == nex_upper;
86-
il += lst == nex_lower;
40+
auto bad = [&](const Ray& a, const Ray& b, const Ray& c) {
41+
P p1 = a.isect(b), p2 = b.isect(c);
42+
if (dot(p2-p1,b.dp) <= EPS) { // this EPS is required ...
43+
if (cross(a.dp,c.dp) <= 0) return 2; // isect(a,b,c) = empty
44+
return 1; // isect(a,c) == isect(a,b,c), remove b
45+
}
46+
return 0; // all the rays matter
47+
};
48+
#define reduce(t) \
49+
while (sz(poly) > 1) { \
50+
int b = bad(poly.at(sz(poly)-2),poly.bk,t); \
51+
if (b == 2) return {}; \
52+
if (b == 1) poly.pop_back(); \
53+
else break; \
54+
}
55+
deque<Ray> poly;
56+
each(t,rays) { reduce(t); poly.pb(t); }
57+
for(;;poly.pop_front()) {
58+
reduce(poly[0]);
59+
if (!bad(poly.bk,poly[0],poly[1])) break;
8760
}
88-
return ans;
89-
}
90-
91-
Half plane_right(P a, P b) { // half-plane to right of a -> b
92-
return {b.s-a.s,a.f-b.f,(b.s-a.s)*a.f+(a.f-b.f)*a.s}; }
93-
Half plane_through(P p, P dir) { // half-plane through p in direction dir
94-
return {dir.f,dir.s,dot(p,dir)}; }
61+
assert(sz(poly) >= 3); // if you reach this point, area should be nonzero
62+
vP poly_points; F0R(i,sz(poly))
63+
poly_points.pb(poly[i].isect(poly[(i+1)%sz(poly)]));
64+
return poly_points;
65+
}

‎Implementations/content/geometry (13)/Polygons/HalfPlaneIsect2.h

Copy file name to clipboardExpand all lines: Implementations/content/geometry (13)/Polygons/HalfPlaneIsect2.h
-62Lines changed: 0 additions & 62 deletions
This file was deleted.

‎Implementations/content/geometry (13)/Polygons/_old_HalfPlaneIsect.h

Copy file name to clipboard
+94Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
/**
2+
* Description: half-plane intersection area
3+
* Time: O(N\log N)
4+
* Source: Own
5+
* http://www.cs.umd.edu/class/spring2020/cmsc754/Lects/lect06-duality.pdf might be of interest
6+
* I would be interested if someone has a simpler implementation ...
7+
* Verification:
8+
* https://open.kattis.com/problems/bigbrother
9+
* https://icpc.kattis.com/problems/domes
10+
* (planes through two points with integer coordinates <= 10^7)
11+
*/
12+
13+
#include "Point.h"
14+
15+
using Half = AR<T,3>; // half-plane
16+
using vH = V<Half>;
17+
P hp_point(const Half& h) { return {h[0],h[1]}; } // direction of half-plane
18+
P isect(const Half& h0, const Half& h1) { // Cramer's rule to intersect half-planes
19+
AR<T,3> vals;
20+
FOR(i,-1,2) {
21+
int x = (i == 0 ? 2 : 0), y = (i == 1 ? 2 : 1);
22+
vals[1+i] = h0[x]*h1[y]-h0[y]*h1[x];
23+
}
24+
assert(vals[0] != 0); return {vals[1]/vals[0],vals[2]/vals[0]};
25+
}
26+
T eval(const Half& h, T x) { assert(h[1] != 0); // evaluate half-plane at x-coordinate
27+
return (h[2]-h[0]*x)/h[1]; }
28+
T x_isect(const Half& h0, const Half& h1) { return isect(h0,h1).f; } // x-coordinate of intersection
29+
30+
vH construct_lower(P x, vH planes) { // similar to convex hull (by duality)
31+
sort(all(planes),[](const Half& a, const Half& b) {
32+
return cross(hp_point(a),hp_point(b)) > 0; });
33+
vH res{{1,0,x.f}}; // >= x.f
34+
planes.pb({-1,0,-x.s}); // <= x.s
35+
auto lst_x = [&](Half a, Half b) {
36+
if (cross(hp_point(a),hp_point(b)) == 0) // parallel half-planes, remove lower one
37+
return a[2]/a[1] <= b[2]/b[1] ? x.f : x.s;
38+
return x_isect(a,b);
39+
};
40+
each(t,planes) {
41+
while (sz(res) > 1 && lst_x(res.bk,t) <= lst_x(end(res)[-2],res.bk))
42+
res.pop_back();
43+
res.pb(t);
44+
}
45+
return res;
46+
}
47+
48+
T isect_area(vH planes) {
49+
const T BOUND = 1e9; P x{-BOUND,BOUND};
50+
planes.pb({0,1,-BOUND}); // y >= -BOUND
51+
planes.pb({0,-1,-BOUND}); // -y >= -BOUND
52+
vH upper, lower;
53+
each(t,planes) {
54+
if (t[1] == 0) { // vertical line
55+
T quo = t[2]/t[0];
56+
if (t[0] > 0) ckmax(x.f,quo);
57+
else ckmin(x.s,quo); // -x >=
58+
} else if (t[1] > 0) lower.pb(t);
59+
else upper.pb(t);
60+
}
61+
if (x.f >= x.s) return 0;
62+
lower = construct_lower(x,lower);
63+
64+
each(t,upper) t[0] *= -1, t[1] *= -1;
65+
upper = construct_lower({-x.s,-x.f},upper);
66+
each(t,upper) t[0] *= -1, t[1] *= -1;
67+
reverse(all(upper));
68+
int iu = 1, il = 1;
69+
T lst = x.f, lst_dif = eval(upper[1],lst)-eval(lower[1],lst);
70+
T ans = 0;
71+
while (iu < sz(upper)-1 && il < sz(lower)-1) { // sweep vertical line through lower and upper hulls
72+
T nex_upper = x_isect(upper[iu],upper[iu+1]);
73+
T nex_lower = x_isect(lower[il],lower[il+1]);
74+
T nex = min(nex_upper,nex_lower);
75+
T nex_dif = eval(upper[iu],nex)-eval(lower[il],nex);
76+
auto avg_val = [](T a, T b) -> T {
77+
if (a > b) swap(a,b);
78+
if (b <= 0) return 0;
79+
if (a >= 0) return (a+b)/2;
80+
return b/(b-a)*b/2;
81+
};
82+
ans += (nex-lst)*avg_val(lst_dif,nex_dif);
83+
assert(x.f <= nex && nex <= x.s);
84+
lst = nex, lst_dif = nex_dif;
85+
iu += lst == nex_upper;
86+
il += lst == nex_lower;
87+
}
88+
return ans;
89+
}
90+
91+
Half plane_right(P a, P b) { // half-plane to right of a -> b
92+
return {b.s-a.s,a.f-b.f,(b.s-a.s)*a.f+(a.f-b.f)*a.s}; }
93+
Half plane_through(P p, P dir) { // half-plane through p in direction dir
94+
return {dir.f,dir.s,dot(p,dir)}; }

0 commit comments

Comments
0 (0)
Morty Proxy This is a proxified and sanitized view of the page, visit original site.