Skip to content

Commit e61acd5

Browse files
ChilleeChillee
committed
Doing some precision analysis
1 parent 1aa4b39 commit e61acd5

File tree

2 files changed

+137
-19
lines changed

2 files changed

+137
-19
lines changed

‎content/geometry/HalfPlane.h

+7-1
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
#include "Point.h"
1515
#include "sideOf.h"
1616
#include "lineIntersection.h"
17+
#include "lineDistance.h"
1718

1819
typedef Point<double> P;
1920
typedef array<P, 2> Line;
@@ -25,6 +26,7 @@ bool cmp(Line a, Line b) {
2526
auto s = angDiff(a, b);
2627
return s == 0 ? sideOf(sp(b), a[0]) >= 0 : s < 0;
2728
}
29+
double mxErr1 = 0, mxErr2 = 0;
2830
vector<P> halfPlaneIntersection(vector<Line> vs) {
2931
sort(all(vs), cmp);
3032
vector<Line> deq(sz(vs) + 5);
@@ -38,7 +40,11 @@ vector<P> halfPlaneIntersection(vector<Line> vs) {
3840
while (i!=n && ah<at && sideOf(sp(vs[i]),ans[ah])<0) ah++;
3941
auto res = lineInter(sp(vs[i]), sp(deq[at]));
4042
if (res.first != 1) continue;
41-
ans[at++] = res.second;
43+
mxErr1 = max(mxErr1, (lineDist(sp(vs[i]), res.second)));
44+
mxErr2 = max(mxErr2, (lineDist(sp(vs[i]), lineInter2(sp(vs[i]), sp(deq[at])))));
45+
assert((res.second - lineInter2(sp(vs[i]), sp(deq[at]))).dist() < 1e-8);
46+
// ans[at++] = res.second;
47+
ans[at++] = lineInter2(sp(vs[i]), sp(deq[at]));
4248
deq[at] = vs[i];
4349
}
4450
if (at - ah <= 2) return {};

‎fuzz-tests/geometry/HalfPlane.cpp

+130-18
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,6 @@ typedef Point<double> P;
1515
ostream &operator<<(ostream &os, P p) { return cout << "(" << p.x << "," << p.y << ")"; }
1616
#include "../../content/geometry/HalfPlane.h"
1717

18-
1918
namespace sjtu {
2019
typedef double T;
2120
const T EPS = 1e-8;
@@ -148,12 +147,112 @@ T calcArea() {
148147

149148
T halfPlaneIntersection(vector<Line> cur) {
150149
n = 0;
151-
for (auto i: cur)
150+
for (auto i : cur)
152151
add(i[0].x, i[0].y, i[1].x, i[1].y);
153152
convexIntersection();
154153
return calcArea();
155154
}
155+
#undef cross
156156
} // namespace sjtu
157+
namespace mit {
158+
#define eps 1e-8
159+
typedef Point<double> P;
160+
161+
template <class P> pair<int, P> lineIntersection(P s1, P e1, P s2, P e2) {
162+
auto d = (e1 - s1).cross(e2 - s2);
163+
if (d == 0) // if parallel
164+
return {-((e1 - s1).cross(s2 - s1) == 0 || s2 == e2), P(0, 0)};
165+
else
166+
return {1, s2 - (e2 - s2) * (e1 - s1).cross(s2 - s1) / d};
167+
}
168+
169+
struct Line {
170+
P P1, P2;
171+
// Right hand side of the ray P1 -> P2
172+
explicit Line(P a = P(), P b = P()) : P1(a), P2(b){};
173+
P intpo(Line y) {
174+
auto res = lineIntersection(P1, P2, y.P1, y.P2);
175+
assert(res.first == 1);
176+
return res.second;
177+
}
178+
P dir() { return P2 - P1; }
179+
bool contains(P x) { return (P2 - P1).cross(x - P1) < eps; }
180+
bool out(P x) { return !contains(x); }
181+
};
182+
183+
template <class T> bool mycmp(Point<T> a, Point<T> b) {
184+
// return atan2(a.y, a.x) < atan2(b.y, b.x);
185+
if (a.x * b.x < 0)
186+
return a.x < 0;
187+
if (abs(a.x) < eps) {
188+
if (abs(b.x) < eps)
189+
return a.y > 0 && b.y < 0;
190+
if (b.x < 0)
191+
return a.y > 0;
192+
if (b.x > 0)
193+
return true;
194+
}
195+
if (abs(b.x) < eps) {
196+
if (a.x < 0)
197+
return b.y < 0;
198+
if (a.x > 0)
199+
return false;
200+
}
201+
return a.cross(b) > 0;
202+
}
203+
204+
bool cmp(Line a, Line b) { return mycmp(a.dir(), b.dir()); }
205+
206+
double Intersection_Area(vector<Line> b) {
207+
sort(b.begin(), b.end(), cmp);
208+
int n = b.size();
209+
int q = 1, h = 0, i;
210+
vector<Line> ca(b.size() + 10);
211+
for (i = 0; i < n; i++) {
212+
while (q < h && b[i].out(ca[h].intpo(ca[h - 1])))
213+
h--;
214+
while (q < h && b[i].out(ca[q].intpo(ca[q + 1])))
215+
q++;
216+
ca[++h] = b[i];
217+
if (q < h && abs(ca[h].dir().cross(ca[h - 1].dir())) < eps) {
218+
if (ca[h].dir().dot(ca[h - 1].dir()) > 0) {
219+
h--;
220+
if (b[i].out(ca[h].P1))
221+
ca[h] = b[i];
222+
} else {
223+
// The area is either 0 or infinite.
224+
// If you have a bounding box, then the area is definitely 0.
225+
return 0;
226+
}
227+
}
228+
}
229+
while (q < h - 1 && ca[q].out(ca[h].intpo(ca[h - 1])))
230+
h--;
231+
while (q < h - 1 && ca[h].out(ca[q].intpo(ca[q + 1])))
232+
q++;
233+
// Intersection is empty. This is sometimes different from the case when
234+
// the intersection area is 0.
235+
if (h - q <= 1)
236+
return 0;
237+
ca[h + 1] = ca[q];
238+
vector<P> s;
239+
for (i = q; i <= h; i++)
240+
s.push_back(ca[i].intpo(ca[i + 1]));
241+
s.push_back(s[0]);
242+
// for (auto i: s) cout<<i<<' ';
243+
// cout<<endl;
244+
double ans = 0;
245+
for (i = 0; i < (int)s.size() - 1; i++)
246+
ans += s[i].cross(s[i + 1]);
247+
return ans / 2;
248+
}
249+
} // namespace mit
250+
vector<mit::Line> convert(vector<Line> x) {
251+
vector<mit::Line> res;
252+
for (auto i : x)
253+
res.push_back(mit::Line(i[1], i[0]));
254+
return res;
255+
}
157256

158257
const double INF = 100;
159258
const double EPS = 1e-8;
@@ -203,8 +302,8 @@ void testEmpty() {
203302
assert(sz(res) == 0);
204303
}
205304
void testRandom() {
206-
int lim = 1e1;
207-
for (int i = 0; i < 10000000; i++) {
305+
int lim = 3;
306+
for (int i = 0; i < 1000000; i++) {
208307
vector<Line> t;
209308
for (int i = 0; i < 6; i++) {
210309
Line cand{P(0, 0), P(0, 0)};
@@ -215,16 +314,17 @@ void testRandom() {
215314
addInf(t, lim);
216315
auto res = halfPlaneIntersection(t);
217316
double area = sjtu::halfPlaneIntersection(t);
218-
double resArea = polygonArea2(res)/2;
219-
if (isnan(resArea)) resArea = 0;
317+
// double resArea = polygonArea2(res)/2;
318+
double resArea = mit::Intersection_Area(convert(t));
220319
double diff = abs(area - resArea);
221320
if (diff > EPS) {
222-
cout << diff << ' ' << area << ' ' <<resArea<< endl;
321+
cout << diff << ' ' << area << ' ' << resArea << endl;
223322
for (auto i : t)
224323
cout << i[0] << "->" << i[1] << ' ';
225324
cout << endl;
226325
for (auto i : t)
227-
cout << "{P"<<i[0] << "," << "P"<<i[1] << "},";
326+
cout << "{P" << i[0] << ","
327+
<< "P" << i[1] << "},";
228328
cout << endl;
229329
for (auto i : res)
230330
cout << i << ',';
@@ -235,18 +335,30 @@ void testRandom() {
235335
}
236336

237337
int main() {
238-
test1();
239-
testInf();
240-
testLine();
241-
testPoint();
242-
testEmpty();
243-
// testRandom();
338+
// test1();
339+
// testInf();
340+
// testLine();
341+
// testPoint();
342+
// testEmpty();
343+
srand(time(0));
344+
testRandom();
244345
// Case that messes with precision
245-
vector<Line> t({{P(8,9),P(8,2)},{P(3,9),P(5,2)},{P(8,2),P(8,3)},{P(7,2),P(1,7)},{P(1,0),P(7,1)},{P(9,2),P(5,6)},{P(10,10),P(-10,10)},{P(-10,10),P(-10,-10)},{P(-10,-10),P(10,-10)},{P(10,-10),P(10,10)}});
346+
vector<Line> t({{P(8, 9), P(8, 2)},
347+
{P(3, 9), P(5, 2)},
348+
{P(8, 2), P(8, 3)},
349+
{P(7, 2), P(1, 7)},
350+
{P(1, 0), P(7, 1)},
351+
{P(9, 2), P(5, 6)},
352+
{P(10, 10), P(-10, 10)},
353+
{P(-10, 10), P(-10, -10)},
354+
{P(-10, -10), P(10, -10)},
355+
{P(10, -10), P(10, 10)}});
246356
auto res = halfPlaneIntersection(t);
247-
cout<<polygonArea2(res)<<endl;
248-
cout << sjtu::halfPlaneIntersection(t)<<endl;
249-
// Failure case for MIT's half plane cod
357+
cout<<fixed<<setprecision(30)<<mxErr1<<' '<<mxErr2<<endl;
358+
cout << polygonArea2(res) << endl;
359+
cout << sjtu::halfPlaneIntersection(t) << endl;
360+
cout << mit::Intersection_Area(convert(t))<<endl;
361+
// Failure case for mit's half plane cod
250362
// vector<Line> t({{P(9, 8), P(9, 1)}, {P(3, 3), P(3, 5)}, {P(7, 6), P(0, 8)}});
251363
// Failure case for old code.
252364
// vector<Line> t({{P(3,0),P(3,3)},{P(5,3), P(5,0)}, {P(4,-2), P(0,1)}, {P(3,-1), P(0,-1)}});

0 commit comments

Comments
 (0)