「CodeVs 2819」无尽的毁灭-Voronoi 图

给出 $n$ 个关键点,求出满足以下条件的点的个数及坐标:

  • 距离它最近的关键点有至少 $k$ 个

链接

CodeVs 2819

Voronoi 图

Voronoi diagram

Voronoi 图就是根据 $n$ 个在平面上不重合种子点,把平面分成 $n$ 个区域,使得每个区域内的点到它所在区域的种子点的距离比到其它区域种子点的距离近。

题解

由 Voronoi 图的定义可知,Voronoi 边上的点(Voronoi 顶点)到相邻区域内的种子点的距离相等,所以答案一定在这些点中,它们只要满足相邻至少 $k$ 个区域就可以了。

具体做法,我们先读入点,建出 Delaunay 三角网,由于 Voronoi 图是 Delaunay 三角网的对偶图,然后用最左转线算法求出其对偶图,并记录每个 Voronoi 顶点相邻哪些 Delaunay 边,由于 Voronoi 顶点是 Delaunay 边构成的封闭区域的中心。

所以我们在每个点相邻的 Delaunay 边中枚举相邻边,求出这两条边的中垂线,然后求出其交点就是 Voronoi 顶点的坐标,由于答案要按水平序输出且方便去重,我们可以用 set 来维护。

由于数据范围较大,计算几何的常数也较大,这里 Delaunay 三角剖分采用分治法实现,当然也可以用扫描线等 $O(n \log n)$ 的算法实现。

时间复杂度 $O(n \log n)$。

一些坑点

Delaunay 后建边的时候要注意判断平行的边,对于当前点的双向链表中,若出现平行边,选择离当前点最近的点,否则最左转线后可能会出现重边和一些跳跃链接的边,导致最后求出的交点在区域外,被这个坑 WA 了一页

代码

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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
/**
* Copyright (c) 2017, xehoth
* All rights reserved.
* 「CodeVS 2819」无尽的毁灭 18-11-2017
* Voronoi 图
* @author xehoth
*/
#include <bits/stdc++.h>
namespace {
const double EPS = 1e-6;
const int MAXN = 40010;
const int MAXM = MAXN * 2;
const double PI2 = M_PI * 2;
int n, k;
template <typename T>
inline T square(const T &x) {
return x * x;
}
struct Point {
double x, y;
int id;
Point(double x = 0, double y = 0, int id = -1) : x(x), y(y), id(id) {}
inline bool operator<(const Point &a) const {
return x < a.x || (fabs(x - a.x) < EPS && y < a.y);
}
inline bool operator==(const Point &a) const {
return fabs(x - a.x) < EPS && fabs(y - a.y) < EPS;
}
inline double dist2(const Point &a) {
return (x - a.x) * (x - a.x) + (y - a.y) * (y - a.y);
}
inline Point operator-(const Point &p) const {
return Point(x - p.x, y - p.y);
}
inline double dis(const Point &p) const {
return sqrt(square(x - p.x) + square(y - p.y));
}
inline Point operator+(const Point &p) const {
return Point(x + p.x, y + p.y);
}
inline double operator*(const Point &p) const { return x * p.y - y * p.x; }
inline Point operator*(const double i) const { return Point(x * i, y * i); }
inline double operator^(const Point &p) const { return x * p.x + y * p.y; }
inline Point operator/(const double i) const { return Point(x / i, y / i); }
};
namespace Delaunay {
struct Point3D {
double x, y, z;
Point3D(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
Point3D(const Point &p) { x = p.x, y = p.y, z = p.x * p.x + p.y * p.y; }
inline Point3D operator-(const Point3D &a) const {
return Point3D(x - a.x, y - a.y, z - a.z);
}
inline double dot(const Point3D &a) { return x * a.x + y * a.y + z * a.z; }
};
struct Edge {
int id;
std::list<Edge>::iterator c;
Edge(int id = 0) { this->id = id; }
};
inline int cmp(double v) { return fabs(v) > EPS ? (v > 0 ? 1 : -1) : 0; }
inline double cross(const Point &o, const Point &a, const Point &b) {
return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}
inline Point3D cross(const Point3D &a, const Point3D &b) {
return Point3D(a.y * b.z - a.z * b.y, -a.x * b.z + a.z * b.x,
a.x * b.y - a.y * b.x);
}
inline int inCircle(const Point &a, Point b, Point c, const Point &p) {
if (cross(a, b, c) < 0) std::swap(b, c);
Point3D a3(a), b3(b), c3(c), p3(p);
b3 = b3 - a3, c3 = c3 - a3, p3 = p3 - a3;
Point3D f = cross(b3, c3);
return cmp(p3.dot(f));
}
inline int intersection(const Point &a, const Point &b, const Point &c,
const Point &d) {
return cmp(cross(a, c, b)) * cmp(cross(a, b, d)) > 0 &&
cmp(cross(c, a, d)) * cmp(cross(c, d, b)) > 0;
}
class Delaunay {
public:
std::list<Edge> head[MAXN];
Point p[MAXN];
int n, rename[MAXN];
inline void init(int n, Point *p) {
memcpy(this->p, p, sizeof(Point) * n);
std::sort(this->p, this->p + n);
for (register int i = 0; i < n; i++) rename[p[i].id] = i;
this->n = n;
divide(0, n - 1);
}
inline void addEdge(int u, int v) {
head[u].push_front(Edge(v));
head[v].push_front(Edge(u));
head[u].begin()->c = head[v].begin();
head[v].begin()->c = head[u].begin();
}
void divide(int l, int r) {
if (r - l <= 2) {
for (register int i = l; i <= r; i++)
for (register int j = i + 1; j <= r; j++) addEdge(i, j);
return;
}
register int mid = (l + r) / 2;
divide(l, mid), divide(mid + 1, r);
std::list<Edge>::iterator it;
register int nowl = l, nowr = r;
for (register int update = 1; update;) {
update = 0;
Point ptL = p[nowl], ptR = p[nowr];
for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
Point t = p[it->id];
register double v = cross(ptR, ptL, t);
if (cmp(v) > 0 ||
(cmp(v) == 0 && ptR.dist2(t) < ptR.dist2(ptL))) {
nowl = it->id, update = 1;
break;
}
}
if (update) continue;
for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
Point t = p[it->id];
register double v = cross(ptL, ptR, t);
if (cmp(v) < 0 ||
(cmp(v) == 0 && ptL.dist2(t) < ptL.dist2(ptR))) {
nowr = it->id, update = 1;
break;
}
}
}
addEdge(nowl, nowr);
for (;;) {
Point ptL = p[nowl], ptR = p[nowr];
register int ch = -1, side = 0;
for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
if (cmp(cross(ptL, ptR, p[it->id])) > 0 &&
(ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0))
ch = it->id, side = -1;
}
for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
if (cmp(cross(ptR, p[it->id], ptL)) > 0 &&
(ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0))
ch = it->id, side = 1;
}
if (ch == -1) break;
if (side == -1) {
for (it = head[nowl].begin(); it != head[nowl].end();) {
if (intersection(ptL, p[it->id], ptR, p[ch])) {
head[it->id].erase(it->c);
head[nowl].erase(it++);
} else
it++;
}
nowl = ch, addEdge(nowl, nowr);
} else {
for (it = head[nowr].begin(); it != head[nowr].end();) {
if (intersection(ptR, p[it->id], ptL, p[ch])) {
head[it->id].erase(it->c), head[nowr].erase(it++);
} else {
it++;
}
}
nowr = ch, addEdge(nowl, nowr);
}
}
}
inline bool parallel(const Point &a, const Point &b, const Point &c) {
return fabs((b - a) * (c - a)) < EPS;
}
inline void getEdge(std::vector<std::pair<int, int> > &ret) {
ret.reserve(n);
std::list<Edge>::iterator it, itp;
std::set<std::pair<int, int> > vis;
for (register int i = 0; i < n; i++) {
for (it = head[i].begin(); it != head[i].end(); it++) {
if (it->id < i) continue;
Point now = p[it->id];
for (itp = head[i].begin(); itp != head[i].end(); ++itp) {
if (itp->id < i) continue;
if (parallel(p[i], p[it->id], p[itp->id]) &&
p[itp->id].dist2(p[i]) < now.dist2(p[i])) {
now = p[itp->id];
}
}
if (vis.find(std::make_pair(p[i].id, now.id)) == vis.end()) {
vis.insert(std::make_pair(p[i].id, now.id));
ret.push_back(std::make_pair(p[i].id, now.id));
}
}
}
}
};
} // namespace Delaunay
namespace PlanarGraph {
Point p[MAXN];
struct Edge {
int u, v;
double angle;
Edge(int u, int v) : u(u), v(v) {
angle = atan2(p[v - 1].y - p[u - 1].y, p[v - 1].x - p[u - 1].x);
if (angle < 0) angle += PI2;
}
Edge() {}
} edge[MAXM];
bool vis[MAXM + 1];
int regionCnt, infArea, rank[MAXM + 1], near[MAXM + 1];
std::vector<int> et[MAXN + 1];
std::vector<Edge> vc[MAXN];
inline void findRegion(int x, int id) {
if (vis[id]) return;
double area = 0;
while (!vis[id]) {
area += p[x - 1] * p[edge[id].v - 1];
vis[id] = true, near[id] = regionCnt, x = edge[id].v;
vc[regionCnt].push_back(edge[id]);
if (!rank[id ^ 1])
id = et[x].back();
else
id = et[x][rank[id ^ 1] - 1];
}
if (area < 0) infArea = regionCnt;
regionCnt++;
}
inline void findDualGraph(const int n, const int m) {
static std::pair<double, int> tmp[MAXM + 1];
for (register int i = 0; i != m << 1; i++)
tmp[i] = std::make_pair(edge[i].angle, i);
std::sort(tmp, tmp + (m << 1));
for (int i = 0, id; i != m << 1; i++) {
id = tmp[i].second;
const Edge &e = edge[id];
rank[id] = et[e.u].size(), et[e.u].push_back(id);
}
for (register int i = 1; i <= n; i++)
for (register int j = 0; j != (int)et[i].size(); j++)
findRegion(i, et[i][j]);
}
} // namespace PlanarGraph
namespace Voronoi {
Delaunay::Delaunay delaunay;
PlanarGraph::Edge vorE[MAXM];
struct Line {
Point s, t;
Line(const Point &s, const Point &t) : s(s), t(t) {}
inline Point intersect(const Line &l) const {
return s +
(t - s) * (((s - l.s) * (l.t - l.s)) / ((l.t - l.s) * (t - s)));
}
inline Line getPerpendicularBisector() {
return Line(
Point((s.x + t.x) / 2, (s.y + t.y) / 2),
Point((s.x + t.x) / 2 + s.y - t.y, (s.y + t.y) / 2 + t.x - s.x));
}
};
struct Cmp {
inline bool operator()(const std::pair<double, double> &a,
const std::pair<double, double> &b) const {
if (fabs(a.first - b.first) < EPS) {
if (fabs(a.second - b.second) < EPS) return false;
return a.second < b.second;
}
return a.first < b.first;
}
};
const double INF = 1e7;
inline void insert(PlanarGraph::Edge *p1, PlanarGraph::Edge *p2,
std::set<std::pair<double, double>, Cmp> &ret) {
// parallel (p1, p2)
if (!(fabs((PlanarGraph::p[p1->v - 1] - PlanarGraph::p[p1->u - 1]) *
(PlanarGraph::p[p2->v - 1] - PlanarGraph::p[p2->u - 1])) <
EPS)) {
// circle center
Point o = Line(PlanarGraph::p[p1->u - 1], PlanarGraph::p[p1->v - 1])
.getPerpendicularBisector()
.intersect(Line(PlanarGraph::p[p2->u - 1],
PlanarGraph::p[p2->v - 1])
.getPerpendicularBisector());
// valid
if (fabs(o.x) <= INF && fabs(o.y) <= INF && !isinf(o.x) &&
!isnan(o.x) && !isinf(o.y) && !isnan(o.y)) {
ret.insert(std::make_pair(o.x, o.y));
}
}
}
inline void getCenter(std::set<std::pair<double, double>, Cmp> &ret, int i) {
PlanarGraph::Edge *p1 = &PlanarGraph::vc[i][0],
*p2 = &PlanarGraph::vc[i].back();
insert(p1, p2, ret);
if (i == PlanarGraph::infArea) return;
for (register int j = 0; j < (int)PlanarGraph::vc[i].size() - 1; j++) {
p1 = &PlanarGraph::vc[i][j], p2 = &PlanarGraph::vc[i][j + 1];
insert(p1, p2, ret);
}
}
inline void buildVoronoi() {
delaunay.init(n, PlanarGraph::p);
std::vector<std::pair<int, int> > edge;
delaunay.getEdge(edge);
register int m = edge.size();
for (register int i = 0; i < m; i++) {
if (edge[i].first > edge[i].second)
std::swap(edge[i].first, edge[i].second);
PlanarGraph::edge[i << 1] =
PlanarGraph::Edge(edge[i].first + 1, edge[i].second + 1);
PlanarGraph::edge[i << 1 | 1] =
PlanarGraph::Edge(edge[i].second + 1, edge[i].first + 1);
}
PlanarGraph::findDualGraph(n, m);
for (register int i = 0, a; i != m; i++) {
a = i << 1;
vorE[i].u = PlanarGraph::near[a];
vorE[i].v = PlanarGraph::near[a ^ 1];
}
static int deg[MAXN + 1];
for (register int i = 0; i < m; i++) deg[vorE[i].u]++, deg[vorE[i].v]++;
std::set<std::pair<double, double>, Cmp> ans;
for (register int i = 0; i < PlanarGraph::regionCnt; i++) {
if (deg[i] >= k) getCenter(ans, i);
}
std::cout << ans.size() << '\n';
for (std::set<std::pair<double, double>, Cmp>::iterator it = ans.begin();
it != ans.end(); ++it) {
std::cout << std::fixed << std::setprecision(4) << it->first << ' '
<< it->second << '\n';
}
}
} // namespace Voronoi
inline void solve() {
std::cin >> n >> k;
for (register int i = 0; i < n; i++) {
std::cin >> PlanarGraph::p[i].x >> PlanarGraph::p[i].y;
PlanarGraph::p[i].id = i;
}
Voronoi::buildVoronoi();
}
} // namespace
int main() {
std::ios::sync_with_stdio(false), std::cin.tie(NULL), std::cout.tie(NULL);
solve();
return 0;
}
分享到