1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
#include "../util.h"
using pt = complex<lll>;
// Kreuzprodukt, 0, falls kollinear.
auto cross(pt a, pt b) {return imag(conj(a) * b);}
auto cross(pt p, pt a, pt b) {return cross(a - p, b - p);}
#pragma GCC diagnostic ignored "-Wunused-variable"
#include <geometry/delaunay.cpp>
vector<pt> convexHull(vector<pt> pts){
sort(all(pts), [](const pt& a, const pt& b){
return real(a) == real(b) ? imag(a) < imag(b)
: real(a) < real(b);
});
pts.erase(unique(all(pts)), pts.end());
int k = 0;
vector<pt> h(2 * sz(pts));
auto half = [&](auto begin, auto end, int t) {
for (auto it = begin; it != end; it++) {
while (k > t && cross(h[k-2], h[k-1], *it) < 0) k--; //allow collinear points!
h[k++] = *it;
}};
half(all(pts), 1); // Untere Hülle.
half(next(pts.rbegin()), pts.rend(), k); // Obere Hülle.
h.resize(k);
return h;
}
lll area(const vector<pt>& poly) { //poly[0] == poly.back()
lll res = 0;
for (int i = 0; i + 1 < sz(poly); i++)
res += cross(poly[i], poly[i + 1]);
return res;
}
// Liegt p auf der Strecke a-b?
bool pointInLineSegment(pt a, pt b, pt p) {
if (cross(a, b, p) != 0) return false;
auto dist = norm(a - b);
return norm(a - p) < dist && norm(b - p) < dist;
}
// Test auf Streckenschnitt zwischen a-b und c-d.
// (nur intern)
bool lineSegmentIntersection(pt a, pt b, pt c, pt d) {
if (cross(a, b, c) == 0 && cross(a, b, d) == 0) {
return pointInLineSegment(a,b,c) ||
pointInLineSegment(a,b,d) ||
pointInLineSegment(c,d,a) ||
pointInLineSegment(c,d,b);
}
return cross(a, b, c) * cross(a, b, d) < 0 &&
cross(c, d, a) * cross(c, d, b) < 0;
}
// 1 => c links von a->b
// 0 => a, b und c kolliniear
// -1 => c rechts von a->b
int ccw(pt a, pt b, pt c) {
auto orien = cross(b - a, c - a);
return (orien > 0) - (orien < 0);
}
bool inOutCirc(pt a, pt b, pt c, pt p) {
lll p2 = norm(p);
lll A = norm(a)-p2;
lll B = norm(b)-p2;
lll C = norm(c)-p2;
return ccw(a, b, c) * (cross(p, a, b)*C + cross(p, b, c)*A + cross(p, c, a)*B) > 0;
}
void stress_test(ll LIM, ll range) {
ll queries = 0;
for (int tries = 0; tries < LIM; tries++) {
int n = Random::integer<int>(3, 30);
auto ps = Random::points<lll>(n, -range, range);
bool skip = true;
for (int i = 2; i < n; i++) skip &= cross(ps[i-2], ps[i-1], ps[i]) == 0;
if (skip) continue;
for (int i = 0; i < n; i++) {
for (int j = 0; j < i; j++) {
skip |= ps[i] == ps[j];
}
}
if (skip) continue;
auto hull = convexHull(ps);
lll expectedArea = area(hull);
hull.pop_back();
auto got = delaunay(ps);
if (sz(got) % 3 != 0) cerr << "error: not triangles" << FAIL;
if (sz(got) / 3 + sz(hull) - 3 + 1 != 2 * sz(ps) - 4) cerr << "error: wrong number" << FAIL;
//all triangles should be oriented ccw
lll gotArea = 0;
for (int i = 0; i < sz(got); i += 3) gotArea += cross(got[i], got[i+1], got[i+2]);
if (gotArea != expectedArea) cerr << "error: wrong area" << FAIL;
for (int i = 0; i < sz(got); i++) {
int ii = i + 1;
if (i / 3 != ii / 3) ii -= 3;
for (int j = 0; j < i; j++) {
int jj = j + 1;
if (j / 3 != jj / 3) jj -= 3;
if (got[i] == got[j] && got[ii] == got[jj]) cerr << "error: dublicate" << FAIL;
if (lineSegmentIntersection(got[i], got[ii], got[j], got[jj])) cerr << "error: intersection" << FAIL;
}
bool seen = false;
for (pt p : ps) seen |= p == got[i];
if (!seen) cerr << "error: invalid point" << FAIL;
}
for (int i = 0; i < sz(got); i += 3) {
for (pt p : ps) {
if (p == got[i]) continue;
if (p == got[i+1]) continue;
if (p == got[i+2]) continue;
if (inOutCirc(got[i], got[i+1], got[i+2], p)) cerr << "error: not delaunay" << FAIL;
}
}
queries += n;
}
cerr << "tested random queries: " << queries << endl;
}
constexpr int N = 100'000;
void performance_test() {
timer t;
auto ps = Random::points<lll>(N, -1'000'000'000, 1'000'000'000);
t.start();
auto got = delaunay(ps);
t.stop();
hash_t hash = sz(got);
if (t.time > 500) cerr << "too slow: " << t.time << FAIL;
cerr << "tested performance: " << t.time << "ms (hash: " << hash << ")" << endl;
}
int main() {
if (!sanitize) {
stress_test(100'000, 10);
stress_test(100'000, 10'000);
stress_test(100'000, 1'000'000'000);
performance_test();
} else {
stress_test(10'000, 10);
stress_test(10'000, 10'000);
stress_test(10'000, 1'000'000'000);
}
}
|
