Skip to content

Commit 0a4e8d3

Browse files
ChilleeChillee
committed
Added fuzz tests with SJTU code
1 parent e5b85b7 commit 0a4e8d3

File tree

1 file changed

+169
-119
lines changed

1 file changed

+169
-119
lines changed

‎fuzz-tests/geometry/HalfPlane.cpp‎

Lines changed: 169 additions & 119 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
#include <bits/stdc++.h>
22
using namespace std;
33

4-
#define rep(i, a, b) for(int i = a; i < int(b); ++i)
5-
#define trav(a, v) for(auto& a : v)
4+
#define rep(i, a, b) for (int i = a; i < int(b); ++i)
5+
#define trav(a, v) for (auto &a : v)
66
#define all(x) x.begin(), x.end()
77
#define sz(x) (int)(x).size()
88

@@ -16,179 +16,229 @@ typedef Point<double> P;
1616

1717
ostream &operator<<(ostream &os, P p) { return cout << "(" << p.x << "," << p.y << ")"; }
1818

19-
namespace MIT {
20-
#define eps 1e-8
21-
typedef Point<double> P;
19+
namespace sjtu {
20+
const double EPS = 1e-8;
21+
inline int sign(double a) { return a < -EPS ? -1 : a > EPS; }
22+
struct Point {
23+
double x, y;
24+
Point() {}
25+
Point(double _x, double _y) : x(_x), y(_y) {}
26+
Point operator+(const Point &p) const { return Point(x + p.x, y + p.y); }
27+
Point operator-(const Point &p) const { return Point(x - p.x, y - p.y); }
28+
Point operator*(double d) const { return Point(x * d, y * d); }
29+
Point operator/(double d) const { return Point(x / d, y / d); }
30+
bool operator<(const Point &p) const {
31+
int c = sign(x - p.x);
32+
if (c)
33+
return c == -1;
34+
return sign(y - p.y) == -1;
35+
}
36+
double dot(const Point &p) const { return x * p.x + y * p.y; }
37+
double det(const Point &p) const { return x * p.y - y * p.x; }
38+
double alpha() const { return atan2(y, x); }
39+
double distTo(const Point &p) const {
40+
double dx = x - p.x, dy = y - p.y;
41+
return hypot(dx, dy);
42+
}
43+
double alphaTo(const Point &p) const {
44+
double dx = x - p.x, dy = y - p.y;
45+
return atan2(dy, dx);
46+
}
47+
// clockwise
48+
Point rotAlpha(const double &alpha, const Point &o = Point(0, 0)) const {
49+
double nx = cos(alpha) * (x - o.x) + sin(alpha) * (y - o.y);
50+
double ny = -sin(alpha) * (x - o.x) + cos(alpha) * (y - o.y);
51+
return Point(nx, ny) + o;
52+
}
53+
Point rot90() const { return Point(-y, x); }
54+
Point unit() { return *this / abs(); }
55+
void read() { scanf("%lf%lf", &x, &y); }
56+
double abs() { return hypot(x, y); }
57+
double abs2() { return x * x + y * y; }
58+
void write() { cout << "(" << x << "," << y << ")" << endl; }
59+
};
60+
#define cross(p1, p2, p3) ((p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x) * (p2.y - p1.y))
61+
#define crossOp(p1, p2, p3) sign(cross(p1, p2, p3))
2262

23-
template<class P>
24-
pair<int, P> lineIntersection(P s1, P e1, P s2, P e2) {
25-
auto d = (e1-s1).cross(e2-s2);
26-
if (d == 0) //if parallel
27-
return {-((e1-s1).cross(s2-s1)==0 || s2==e2), P(0,0)};
28-
else
29-
return {1, s2-(e2-s2)*(e1-s1).cross(s2-s1)/d};
63+
Point isSS(Point p1, Point p2, Point q1, Point q2) {
64+
double a1 = cross(q1, q2, p1), a2 = -cross(q1, q2, p2);
65+
return (p1 * a2 + p2 * a1) / (a1 + a2);
3066
}
31-
32-
struct Line {
33-
P P1, P2;
34-
// Right hand side of the ray P1 -> P2
35-
explicit Line(P a = P(), P b = P()) : P1(a), P2(b) {};
36-
P intpo(Line y) {
37-
auto res = lineIntersection(P1, P2, y.P1, y.P2);
38-
assert(res.first == 1);
39-
return res.second;
40-
}
41-
P dir() {
42-
return P2 - P1;
43-
}
44-
bool contains(P x) {
45-
return (P2 - P1).cross(x - P1) < eps;
46-
}
47-
bool out(P x) {
48-
return !contains(x);
49-
}
67+
struct Border {
68+
Point p1, p2;
69+
double alpha;
70+
void setAlpha() { alpha = atan2(p2.y - p1.y, p2.x - p1.x); }
71+
void read() {
72+
p1.read();
73+
p2.read();
74+
setAlpha();
75+
}
5076
};
5177

52-
template<class T>
53-
bool mycmp(Point<T> a, Point<T> b) {
54-
// return atan2(a.y, a.x) < atan2(b.y, b.x);
55-
if (a.x * b.x < 0) return a.x < 0;
56-
if (abs(a.x) < eps) {
57-
if (abs(b.x) < eps) return a.y > 0 && b.y < 0;
58-
if (b.x < 0) return a.y > 0;
59-
if (b.x > 0) return true;
60-
}
61-
if (abs(b.x) < eps) {
62-
if (a.x < 0) return b.y < 0;
63-
if (a.x > 0) return false;
64-
}
65-
return a.cross(b) > 0;
78+
int n;
79+
const int MAX_N_BORDER = 20000 + 10;
80+
Border border[MAX_N_BORDER];
81+
82+
bool operator<(const Border &a, const Border &b) {
83+
int c = sign(a.alpha - b.alpha);
84+
if (c != 0)
85+
return c == 1;
86+
return crossOp(b.p1, b.p2, a.p1) >= 0;
87+
}
88+
89+
bool operator==(const Border &a, const Border &b) { return sign(a.alpha - b.alpha) == 0; }
90+
91+
void add(double x, double y, double nx, double ny) {
92+
border[n].p1 = Point(x, y);
93+
border[n].p2 = Point(nx, ny);
94+
border[n].setAlpha();
95+
n++;
6696
}
6797

68-
bool cmp(Line a, Line b) {
69-
return mycmp(a.dir(), b.dir());
98+
Point isBorder(const Border &a, const Border &b) { return isSS(a.p1, a.p2, b.p1, b.p2); }
99+
100+
Border que[MAX_N_BORDER];
101+
int qh, qt;
102+
103+
bool check(const Border &a, const Border &b, const Border &me) {
104+
Point is = isBorder(a, b);
105+
return crossOp(me.p1, me.p2, is) > 0;
70106
}
71107

72-
double Intersection_Area(vector <Line> b) {
73-
sort(b.begin(), b.end(), cmp);
74-
int n = b.size();
75-
int q = 1, h = 0, i;
76-
vector <Line> ca(b.size() + 10);
77-
for (i = 0; i < n; i++) {
78-
while (q < h && b[i].out(ca[h].intpo(ca[h - 1]))) h--;
79-
while (q < h && b[i].out(ca[q].intpo(ca[q + 1]))) q++;
80-
ca[++h] = b[i];
81-
if (q < h && abs(ca[h].dir().cross(ca[h - 1].dir())) < eps) {
82-
h--;
83-
if (b[i].out(ca[h].P1)) ca[h] = b[i];
84-
}
85-
}
86-
while (q < h - 1 && ca[q].out(ca[h].intpo(ca[h - 1]))) h--;
87-
while (q < h - 1 && ca[h].out(ca[q].intpo(ca[q + 1]))) q++;
88-
// Intersection is empty. This is sometimes different from the case when
89-
// the intersection area is 0.
90-
if (h - q <= 1) return 0;
91-
ca[h + 1] = ca[q];
92-
vector <P> s;
93-
for (i = q; i <= h; i++) s.push_back(ca[i].intpo(ca[i + 1]));
94-
s.push_back(s[0]);
95-
// for (auto i: s) cout<<i<<' ';
96-
// cout<<endl;
97-
double ans = 0;
98-
for (i = 0; i < (int) s.size() - 1; i++) ans += s[i].cross(s[i + 1]);
99-
return ans / 2;
108+
void convexIntersection() {
109+
qh = qt = 0;
110+
sort(border, border + n);
111+
n = unique(border, border + n) - border;
112+
for (int i = 0; i < n; ++i) {
113+
Border cur = border[i];
114+
while (qh + 1 < qt && !check(que[qt - 2], que[qt - 1], cur))
115+
--qt;
116+
while (qh + 1 < qt && !check(que[qh], que[qh + 1], cur))
117+
++qh;
118+
que[qt++] = cur;
119+
}
120+
while (qh + 1 < qt && !check(que[qt - 2], que[qt - 1], que[qh]))
121+
--qt;
122+
while (qh + 1 < qt && !check(que[qh], que[qh + 1], que[qt - 1]))
123+
++qh;
100124
}
125+
126+
double calcArea() {
127+
static Point ps[MAX_N_BORDER];
128+
int cnt = 0;
129+
130+
if (qt - qh <= 2) {
131+
return 0;
132+
}
133+
134+
for (int i = qh; i < qt; ++i) {
135+
int next = i + 1 == qt ? qh : i + 1;
136+
ps[cnt++] = isBorder(que[i], que[next]);
137+
}
138+
139+
double area = 0;
140+
for (int i = 0; i < cnt; ++i) {
141+
area += ps[i].det(ps[(i + 1) % cnt]);
142+
}
143+
area /= 2;
144+
area = fabsl(area);
145+
return area;
101146
}
102-
vector<MIT::Line> convert(vector<Line> x) {
103-
vector<MIT::Line> res;
104-
for (auto i: x)
105-
res.push_back(MIT::Line(i[1], i[0]));
106-
return res;
147+
148+
double halfPlaneIntersection(vector<Line> cur) {
149+
n = 0;
150+
for (auto i: cur)
151+
add(i[0].x, i[0].y, i[1].x, i[1].y);
152+
convexIntersection();
153+
return calcArea();
107154
}
155+
} // namespace sjtu
108156

109-
const double INF = 1e1;
157+
const double INF = 100;
110158
const double EPS = 1e-8;
111-
void addInf(vector<Line> &res, double INF=INF) {
159+
void addInf(vector<Line> &res, double INF = INF) {
112160
vector<P> infPts({P(INF, INF), P(-INF, INF), P(-INF, -INF), P(INF, -INF)});
113-
for (int i=0; i<4; i++)
114-
res.push_back({infPts[i], infPts[(i+1)%4]});
115-
}
116-
P randPt(int lim = INF) {
117-
return P(rand() % lim, rand() % lim);
161+
for (int i = 0; i < 4; i++)
162+
res.push_back({infPts[i], infPts[(i + 1) % 4]});
118163
}
164+
P randPt(int lim = INF) { return P(rand() % lim, rand() % lim); }
119165
void test1() {
120-
vector<Line> t({{P(0,0),P(5,0)}, {P(5,-2), P(5,2)}, {P(5,2),P(2,2)}, {P(0,3),P(0,-3)}});
166+
vector<Line> t({{P(0, 0), P(5, 0)}, {P(5, -2), P(5, 2)}, {P(5, 2), P(2, 2)}, {P(0, 3), P(0, -3)}});
121167
auto res = halfPlaneIntersection(t);
122168
assert(polygonArea2(res) == 20);
123169
addInf(t);
124170
res = halfPlaneIntersection(t);
125171
assert(polygonArea2(res) == 20);
126172
}
127173
void testInf() {
128-
vector<Line> t({{P(0,0), P(5,0)}});
174+
vector<Line> t({{P(0, 0), P(5, 0)}});
129175
addInf(t);
130176
auto res = halfPlaneIntersection(t);
131-
assert(polygonArea2(res)/(4*INF*INF) == 1);
132-
t = vector<Line>({{P(0,0), P(5,0)}, {P(0,0), P(0,5)}});
177+
assert(polygonArea2(res) / (4 * INF * INF) == 1);
178+
t = vector<Line>({{P(0, 0), P(5, 0)}, {P(0, 0), P(0, 5)}});
133179
addInf(t);
134180
res = halfPlaneIntersection(t);
135-
assert(polygonArea2(res)/(2*INF*INF) == 1);
181+
assert(polygonArea2(res) / (2 * INF * INF) == 1);
136182
}
137183
void testLine() {
138-
vector<Line> t({{P(0,0), P(5,0)}, {P(5,0), P(0,0)}});
184+
vector<Line> t({{P(0, 0), P(5, 0)}, {P(5, 0), P(0, 0)}});
139185
addInf(t);
140186
auto res = halfPlaneIntersection(t);
141187
assert(sz(res) >= 2);
142188
assert(polygonArea2(res) == 0);
143189
}
144190
void testPoint() {
145-
vector<Line> t({{P(0,0), P(5,0)}, {P(5,0), P(0,0)}, {P(0,0), P(0,5)}, {P(0,5), P(0,0)}});
191+
vector<Line> t({{P(0, 0), P(5, 0)}, {P(5, 0), P(0, 0)}, {P(0, 0), P(0, 5)}, {P(0, 5), P(0, 0)}});
146192
addInf(t);
147193
auto res = halfPlaneIntersection(t);
148194
assert(sz(res) >= 1);
149195
assert(polygonArea2(res) == 0);
150196
}
151197
void testEmpty() {
152-
vector<Line> t({{P(0,0), P(5,0)}, {P(5,0), P(0,0)}, {P(0,0), P(0,5)}, {P(0,5), P(0,0)},
153-
{P(0,2), P(5,2)}});
198+
vector<Line> t(
199+
{{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)}});
154200
addInf(t);
155201
auto res = halfPlaneIntersection(t);
156202
assert(sz(res) == 0);
157203
}
158204
void testRandom() {
159-
for (int i=0; i<1000; i++) {
205+
int lim = 10;
206+
for (int i = 0; i < 1000; i++) {
160207
vector<Line> t;
161-
for (int i=0; i<3; i++){
162-
Line cand{P(0,0), P(0,0)};
208+
for (int i = 0; i < 6; i++) {
209+
Line cand{P(0, 0), P(0, 0)};
163210
while (cand[0] == cand[1])
164-
cand = {randPt(), randPt()};
211+
cand = {randPt(lim), randPt(lim)};
165212
t.push_back(cand);
166213
}
167-
addInf(t);
214+
addInf(t, lim);
168215
auto res = halfPlaneIntersection(t);
169-
double area = MIT::Intersection_Area(convert(t));
170-
double diff = abs(2*area - polygonArea2(res));
216+
double area = sjtu::halfPlaneIntersection(t);
217+
double diff = abs(2 * area - polygonArea2(res));
171218
if (diff > EPS) {
172-
cout<<diff<<' '<<area<<' '<<endl;
173-
for (auto i: t) cout<<i[0]<<"->"<<i[1]<<' ';
174-
cout<<endl;
175-
for (auto i: res) cout<<i<<',';
176-
cout<<endl;
219+
cout << diff << ' ' << area << ' ' << endl;
220+
for (auto i : t)
221+
cout << i[0] << "->" << i[1] << ' ';
222+
cout << endl;
223+
for (auto i : res)
224+
cout << i << ',';
225+
cout << endl;
177226
}
178227
assert(diff < EPS);
179228
}
180229
}
181230

182231
int main() {
183-
// test1();
184-
// testInf();
185-
// testLine();
186-
// testPoint();
187-
// testEmpty();
188-
// testRandom();
189-
vector<Line> t({{P(9,8), P(9,1)}, {P(3,3),P(3,5)},{P(7,6),P(0,8)}});
190-
addInf(t);
191-
cout<<polygonArea2(halfPlaneIntersection(t))<<endl;
192-
cout<<MIT::Intersection_Area(convert(t))<<endl;
193-
cout<<"Tests passed!"<<endl;
232+
test1();
233+
testInf();
234+
testLine();
235+
testPoint();
236+
testEmpty();
237+
testRandom();
238+
// Failure case for MIT
239+
// vector<Line> t({{P(9, 8), P(9, 1)}, {P(3, 3), P(3, 5)}, {P(7, 6), P(0, 8)}});
240+
// addInf(t);
241+
// cout << polygonArea2(halfPlaneIntersection(t)) << endl;
242+
// cout << MIT::Intersection_Area(convert(t)) << endl;
243+
cout << "Tests passed!" << endl;
194244
}

0 commit comments

Comments
 (0)