UVa 12310 - Point Location

contents

  1. 1. Problem
  2. 2. Sample Input
  3. 3. Sample Output
  4. 4. Solution
    1. 4.1. Previously
    2. 4.2. Detail
    3. 4.3. KD tree
    4. 4.4. 線段樹

Problem

給定一平面上 $n$ 個點,接著拉 $m$ 個線段拼湊成不相交的區域,並將標記數個點位置所在的區域。接著有數個詢問「點在哪一個先前標記的同一區域內」。

這一問題常在幾何處理中出現,詳細可查閱 Wikipedia - Point location,在此不額外說明。

Sample Input

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
14 16 5 5
0 0
30 0
40 0
60 0
60 50
0 50
20 40
10 10
30 10
30 20
40 20
50 30
50 40
30 40
1 2
2 9
9 8
8 7
7 9
9 10
10 11
11 3
3 4
4 5
5 6
6 1
12 13
13 14
14 12
2 3
20 20
10 20
35 10
45 39
1 60
28 11
29 14
34 7
40 38
70 1
0 0 0 0

Sample Output

1
2
3
4
5
1
2
3
4
5

Solution

Previously

在大學修計算幾何的時候,知道這一類的題目要透過梯形剖分的方式,搭配的資料結構類似四分樹,分隔的基準是一個線段,將線段兩個端點分別拉出無限延長的垂直線,但這方法對於線段本身是垂直時,處理情況比較複雜,那時構思了好幾天沒個底,只好放棄。現在,在 Cadence 公司上班快滿一年,解決幾何操作不在少數,成長的我應該能解決這一題了吧。

Detail

題目已經保證給予的線段除了端點外,彼此之間不會相交任何一點。便可以對所有線段建立空間分割相關的資料結構,好讓我們對於解決,從詢問點出發往 $y$ 無窮遠的射線,最後會打到哪一個線段。

處理的流程如下:

Step 1. 輸入線段

Step 2. 套上邊界,需要能包含所有詢問點

Step 3. 對每一個相連的線段集合,找出 $y$ 軸最高的點,放出輔助射線

Step 4. 對每一個頂點進行極角排序,進行繞行分組

首先,對於幾何剖分的問題,為處理方便,都會套上一個無限大的外框,把所有頂點包在裡面。你可以透過 GeoGebra 這套數學軟體,將處理過程中的資料丟進去,方便你除錯。如果要批次輸入一大筆,請透過上方的建立按鈕 (Button),將 command 的語法以行為單位放入,語法類似 javascript。

接著,我們需要對每一個線段集合的最高點放出輔助射線 (這部分已經縮減了不少,原則上對每一個頂點放出射線也可以),這些輔助射線用來解決內部的孤立形,要把整個圖串成一個連通圖。否則,當詢問屬於內部的孤立形狀的外部時,我們會缺少足夠的資訊連到包含這個孤立形的外框。

最後,對每一個頂點相連的邊進行極角排序,接著才能決定相鄰的邊,而任兩個相鄰的邊所夾的區域屬於同一個集合,因此我們需要對每一個有方向的邊進行編號,將邊與邊合併到同一個集合。對於每一個詢問,只需要對詢問點放出射線找到接觸的線段,便可以知道所在的集合。

  • 若使用 KD tree,即使交替選擇 $x$-$y$ 軸分割,也沒辦法保證樹高。因為處理的是線段,而不是點,每一次挑選的分隔軸,將會產生三個子樹,中間的子樹為與分隔軸相交的線段。在大部分的情況下,整體時間複雜度落在 $O(n \log n)$,空間複雜度 $O(n)$
  • 若使用線段樹,將在 $x$ 軸上劃分,每一個線段至多被拆成 $O(\log n)$ 個,詢問射線第一個碰觸的線段時,單一詢問的時間複雜度 $O(\log n)$,常數相較於 KD tree 多。整體時間複雜度落在 $O(n \log n)$,空間複雜度 $O(n \log n)$

還有許多 spatial data structure 可以考慮,這裡只挑選兩個出來實驗,上述皆為在線算法。也可以將所有操作離線,這時候將套用掃描線算法,但對於垂直線段如何在算法中維護二元樹,目前遇到一些實作上的問題,這將在未來才能給各位參考。

KD tree

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
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
#include <bits/stdc++.h>
using namespace std;
const float eps = 1e-6;
const float BOX_MX = 2000000;
struct Pt {
float x, y;
Pt() {}
Pt(float a, float b): x(a), y(b) {}
Pt operator-(const Pt &a) const { return Pt(x - a.x, y - a.y); }
Pt operator+(const Pt &a) const { return Pt(x + a.x, y + a.y); }
Pt operator*(const double a) const { return Pt(x * a, y * a); }
bool operator<(const Pt &a) const {
if (fabs(x - a.x) > eps) return x < a.x;
if (fabs(y - a.y) > eps) return y < a.y;
return false;
}
bool operator==(const Pt &a) const {
return fabs(x - a.x) < eps && fabs(y - a.y) < eps;
}
void println() const {
printf("(%lf, %lf)\n", x, y);
}
};
struct Seg {
Pt s, e;
int i;
Seg() {}
Seg(Pt a, Pt b, int i):s(a), e(b), i(i) {}
void println() {
printf("Segment((%lf, %lf), (%lf, %lf))\n", s.x, s.y, e.x, e.y);
}
};
double dot(Pt a, Pt b) {
return a.x * b.x + a.y * b.y;
}
double cross(Pt o, Pt a, Pt b) {
return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
double cross2(Pt a, Pt b) {
return a.x * b.y - a.y * b.x;
}
int between(Pt a, Pt b, Pt c) {
return dot(c - a, b - a) >= -eps && dot(c - b, a - b) >= -eps;
}
int onSeg(Pt a, Pt b, Pt c) {
return between(a, b, c) && fabs(cross(a, b, c)) < eps;
}
int cmpZero(float v) {
if (fabs(v) > eps) return v > 0 ? 1 : -1;
return 0;
}
Pt getIntersect(Seg a, Seg b) {
Pt u = a.s - b.s;
double t = cross2(b.e - b.s, u)/cross2(a.e - a.s, b.e - b.s);
return a.s + (a.e - a.s) * t;
}
struct AngleCmp {
Pt o;
AngleCmp(Pt o = Pt()):o(o) {}
bool operator() (const pair<Pt, int>& ppa, const pair<Pt, int>& ppb) {
Pt pa = ppa.first, pb = ppb.first;
Pt p1 = pa - o, p2 = pb - o;
if (p1.y == 0 && p2.y == 0 && p1.x * p2.x <= 0) return p1.x > p2.x;
if (p1.y == 0 && p1.x >= 0 && p2.y != 0) return true;
if (p2.y == 0 && p2.x >= 0 && p1.y != 0) return false;
if (p1.y * p2.y < 0) return p1.y > p2.y;
double c = cross2(p1, p2);
return c > 0 || (c == 0 && fabs(p1.x) < fabs(p2.x));
}
};
static Pt pts[10005];
static Seg segs[30005];
struct BSP {
static const int MAXN = 60005;
static const int MAXNODE = 60005;
struct Node {
float lx, ly, rx, ry;
Node *ls, *ms, *rs;
void extend(Node *u) {
if (u == NULL)
return ;
lx = min(lx, u->lx), ly = min(ly, u->ly);
rx = max(rx, u->rx), ry = max(ry, u->ry);
}
} nodes[MAXNODE];
int sn[MAXNODE];
Seg *seg[MAXNODE];
float axis[MAXN];
Node *root;
int size;
Node* newNode() {
Node *p = &nodes[size++];
assert(size < MAXNODE);
p->ls = p->ms = p->rs = NULL;
sn[p-nodes] = 0;
return p;
}
Node* _build(int k, Seg segs[], int n) {
if (n == 0)
return NULL;
if (k == 2)
k = 0;
int Lsize = 0, Msize = 0, Rsize = 0;
Seg *L = NULL, *M = NULL, *R = NULL;
if (k == 0) {
for (int i = 0; i < n; i++)
axis[i<<1] = segs[i].s.x, axis[i<<1|1] = segs[i].e.x;
nth_element(axis, axis+n, axis+2*n);
const float mval = axis[n];
L = segs;
R = std::partition(segs, segs+n, [mval](const Seg &s) {
return max(s.s.x, s.e.x) <= mval;
});
M = std::partition(R, segs+n, [mval](const Seg &s) {
return min(s.s.x, s.e.x) <= mval;
});
Msize = segs+n - M;
Rsize = M - R;
Lsize = R - segs;
} else {
for (int i = 0; i < n; i++)
axis[i<<1] = segs[i].s.y, axis[i<<1|1] = segs[i].e.y;
nth_element(axis, axis+n, axis+2*n);
const float mval = axis[n];
L = segs;
R = std::partition(segs, segs+n, [mval](const Seg &s) {
return max(s.s.y, s.e.y) <= mval;
});
M = std::partition(R, segs+n, [mval](const Seg &s) {
return min(s.s.y, s.e.y) <= mval;
});
Msize = segs+n - M;
Rsize = M - R;
Lsize = R - segs;
}
Node *u = newNode();
u->lx = BOX_MX, u->ly = BOX_MX;
u->rx = -BOX_MX, u->ry = -BOX_MX;
if (Lsize == n || Rsize == n || Msize == n) {
sn[u - nodes] = n, seg[u - nodes] = segs;
for (int i = 0; i < n; i++) {
u->lx = min(u->lx, min(segs[i].s.x, segs[i].e.x));
u->rx = max(u->rx, max(segs[i].s.x, segs[i].e.x));
u->ly = min(u->ly, min(segs[i].s.y, segs[i].e.y));
u->ry = max(u->ry, max(segs[i].s.y, segs[i].e.y));
}
} else {
u->ls = _build(k+1, L, Lsize), u->extend(u->ls);
u->ms = _build(k+1, M, Msize), u->extend(u->ms);
u->rs = _build(k+1, R, Rsize), u->extend(u->rs);
}
return u;
}
void build_tree(Seg s[], int m) {
size = 0;
root = _build(0, s, m);
}
Pt q_st, q_ed;
int q_si;
void rayhit(Seg &seg) {
if (seg.s.x == seg.e.x) {
if (cmpZero(seg.s.x - q_st.x) == 0) {
double low = min(seg.s.y, seg.e.y);
if (low > q_st.y && low < q_ed.y) {
q_ed.y = low;
q_si = seg.i;
}
}
return ;
}
if (max(seg.s.x, seg.e.x) < q_st.x || min(seg.s.x, seg.e.x) > q_st.x)
return ;
float y = seg.s.y + (float) (seg.e.y - seg.s.y) * (q_st.x - seg.s.x) / (seg.e.x - seg.s.x);
if (y > q_st.y && y < q_ed.y) {
q_ed.y = y;
q_si = seg.i;
}
}
void search(Node *u) {
if (u == NULL)
return ;
if (u->lx > q_st.x || u->rx < q_st.x || u->ry <= q_st.y || u->ly >= q_ed.y)
return ;
for (int i = 0; i < sn[u - nodes]; i++)
rayhit(seg[u - nodes][i]);
search(u->ls);
search(u->ms);
search(u->rs);
}
pair<int, Pt> raycast(Pt st) {
q_st = st;
q_ed = Pt(st.x, BOX_MX+1);
q_si = -1;
search(root);
return {q_si, q_ed};
}
} tree;
struct Disjoint {
static const int MAXN = 65536;
uint16_t parent[MAXN], weight[MAXN];
void init(int n) {
if (n >= MAXN)
exit(0);
for (int i = 0; i <= n; i++)
parent[i] = i, weight[i] = 1;
}
int findp(int x) {
return parent[x] == x ? x : parent[x] = findp(parent[x]);
}
int joint(int x, int y) {
x = findp(x), y = findp(y);
if (weight[x] >= weight[y])
parent[y] = x, weight[x] += weight[y];
else
parent[x] = y, weight[y] += weight[x];
}
} egroup, sgroup;
int main() {
int n, m, p, q;
while (scanf("%d %d %d %d", &n, &m, &p, &q) == 4 && n) {
for (int i = 0; i < n; i++) {
int x, y;
scanf("%d %d", &x, &y);
pts[i] = Pt(x, y);
}
sgroup.init(n);
for (int i = 0; i < m; i++) {
int st_i, ed_i;
scanf("%d %d", &st_i, &ed_i);
st_i--, ed_i--;
segs[i] = Seg(pts[st_i], pts[ed_i], i);
sgroup.joint(st_i, ed_i);
}
segs[m] = Seg(Pt(BOX_MX, BOX_MX), Pt(BOX_MX, -BOX_MX), m), m++;
segs[m] = Seg(Pt(BOX_MX, -BOX_MX), Pt(-BOX_MX, -BOX_MX), m), m++;
segs[m] = Seg(Pt(-BOX_MX, -BOX_MX), Pt(-BOX_MX, BOX_MX), m), m++;
segs[m] = Seg(Pt(-BOX_MX, BOX_MX), Pt(BOX_MX, BOX_MX), m), m++;
static map<Pt, vector<pair<Pt, uint8_t>>> g; g.clear();
static vector<vector<Pt>> on_seg; on_seg.clear(); on_seg.resize(m);
for (int i = 0; i < m; i++) {
on_seg[segs[i].i].reserve(4);
on_seg[segs[i].i].push_back(segs[i].s);
on_seg[segs[i].i].push_back(segs[i].e);
}
// for (int i = 0; i < n; i++)
// pts[i].println();
// for (int i = 0; i < m; i++)
// segs[i].println();
tree.build_tree(segs, m);
{
Pt top[10005];
for (int i = 0; i < n; i++)
top[i] = pts[i];
for (int i = 0; i < n; i++) {
int gid = sgroup.findp(i);
if (top[gid].y < pts[i].y)
top[gid] = pts[i];
}
for (int i = 0; i < n; i++) {
if (sgroup.findp(i) != i)
continue;
auto p = top[i];
auto hit = tree.raycast(p);
if (hit.first >= 0) {
on_seg[hit.first].emplace_back(hit.second);
g[p].emplace_back(hit.second, 1);
g[hit.second].emplace_back(p, 1);
}
}
}
for (int i = 0; i < m; i++) {
vector<Pt> &a = on_seg[i];
sort(a.begin(), a.end());
a.resize(unique(a.begin(), a.end()) - a.begin());
auto *prev = &g[a[0]];
for (int j = 1; j < a.size(); j++) {
prev->emplace_back(a[j], 0);
prev = &g[a[j]];
prev->emplace_back(a[j-1], 0);
}
}
for (auto &e : g)
sort(e.second.begin(), e.second.end(), AngleCmp(e.first));
static map<Pt, map<Pt, int>> R; R.clear();
int Rsize = 0;
for (auto &e : g) {
int sz = e.second.size();
map<Pt, int> &Rg = R[e.first];
for (auto &f : e.second) {
int &eid = Rg[f.first];
if (eid == 0)
eid = ++Rsize;
}
}
egroup.init(Rsize);
for (auto &e : g) {
int sz = e.second.size();
map<Pt, int> &Rg = R[e.first];
for (int i = sz-1, j = 0; j < sz; i = j++) {
int l = R[e.second[i].first][e.first];
int r = Rg[e.second[j].first];
egroup.joint(l, r);
if (e.second[i].second != 0) {
r = Rg[e.second[i].first];
assert(l > 0 && r > 0);
egroup.joint(l, r);
}
}
}
for (auto &e : g) {
int sz = e.second.size();
int n = 0;
for (int i = 0; i < sz; i++) {
if (e.second[i].second == 0)
e.second[n++] = e.second[i];
}
e.second.resize(n);
}
static int region[65536]; memset(region, 0, sizeof(0)*Rsize);
for (int i = 0; i < p; i++) {
float x, y;
scanf("%f %f", &x, &y);
pair<int, Pt> hit = tree.raycast(Pt(x, y));
if (hit.first < 0)
continue;
if (g.find(hit.second) != g.end()) {
auto &adj = g[hit.second];
AngleCmp cmp(hit.second);
int pos = 0;
pair<Pt, int> q(Pt(x, y), i+1);
pos = lower_bound(adj.begin(), adj.end(), q, cmp) - adj.begin() - 1;
assert(pos >= 0);
int l = R[adj[pos].first][hit.second];
assert(l > 0);
l = egroup.findp(l);
region[l] = i+1;
} else {
auto &adj = on_seg[hit.first];
Pt lpt = adj[0], rpt = adj[1];
int l;
if (cmpZero(cross(lpt, rpt, Pt(x, y))) <= 0)
l = R[lpt][rpt];
else
l = R[rpt][lpt];
assert(l > 0);
l = egroup.findp(l);
region[l] = i+1;
}
}
for (int i = 0; i < q; i++) {
float x, y;
scanf("%f %f", &x, &y);
pair<int, Pt> hit = tree.raycast(Pt(x, y));
if (hit.first < 0) {
printf("0\n");
continue;
}
if (g.find(hit.second) != g.end()) {
auto &adj = g[hit.second];
AngleCmp cmp(hit.second);
int pos = 0;
pair<Pt, int> q(Pt(x, y), i+1);
pos = lower_bound(adj.begin(), adj.end(), q, cmp) - adj.begin() - 1;
assert(pos >= 0);
int l = R[adj[pos].first][hit.second];
assert(l > 0);
l = egroup.findp(l);
printf("%d\n", region[l]);
} else {
auto &adj = on_seg[hit.first];
Pt lpt = adj[0], rpt = adj[1];
int l;
if (cmpZero(cross(lpt, rpt, Pt(x, y))) <= 0)
l = R[lpt][rpt];
else
l = R[rpt][lpt];
assert(l > 0);
l = egroup.findp(l);
printf("%d\n", region[l]);
}
}
}
return 0;
}

線段樹

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
401
402
403
404
405
406
407
408
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
const double BOX_MX = (3000000);
const int AUXID = 200000;
struct Pt {
double x, y;
Pt(double a = 0, double b = 0): x(a), y(b) {}
Pt operator-(const Pt &a) const { return Pt(x - a.x, y - a.y); }
Pt operator+(const Pt &a) const { return Pt(x + a.x, y + a.y); }
Pt operator*(const double a) const { return Pt(x * a, y * a); }
bool operator<(const Pt &a) const {
if (fabs(x - a.x) > eps) return x < a.x;
if (fabs(y - a.y) > eps) return y < a.y;
return false;
}
bool operator==(const Pt &a) const {
return fabs(x - a.x) < eps && fabs(y - a.y) < eps;
}
void println() const {
printf("(%lf, %lf)\n", x, y);
}
};
struct Seg {
Pt s, e;
int i;
Seg() {}
Seg(Pt a, Pt b, int i):s(a), e(b), i(i) {}
void println() const {
printf("Segment((%lf, %lf), (%lf, %lf))\n", s.x, s.y, e.x, e.y);
}
};
double dot(Pt a, Pt b) {
return a.x * b.x + a.y * b.y;
}
double cross(Pt o, Pt a, Pt b) {
return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
double cross2(Pt a, Pt b) {
return a.x * b.y - a.y * b.x;
}
int between(Pt a, Pt b, Pt c) {
return dot(c - a, b - a) >= -eps && dot(c - b, a - b) >= -eps;
}
int onSeg(Pt a, Pt b, Pt c) {
return between(a, b, c) && fabs(cross(a, b, c)) < eps;
}
int cmpZero(double v) {
if (fabs(v) > eps) return v > 0 ? 1 : -1;
return 0;
}
Pt getIntersect(Seg a, Seg b) {
Pt u = a.s - b.s;
double t = cross2(b.e - b.s, u)/cross2(a.e - a.s, b.e - b.s);
return a.s + (a.e - a.s) * t;
}
struct AngleCmp {
Pt o;
AngleCmp(Pt o = Pt()):o(o) {}
bool operator() (const pair<Pt, int>& ppa, const pair<Pt, int>& ppb) {
Pt pa = ppa.first, pb = ppb.first;
Pt p1 = pa - o, p2 = pb - o;
if (p1.y == 0 && p2.y == 0 && p1.x * p2.x <= 0) return p1.x > p2.x;
if (p1.y == 0 && p1.x >= 0 && p2.y != 0) return true;
if (p2.y == 0 && p2.x >= 0 && p1.y != 0) return false;
if (p1.y * p2.y < 0) return p1.y > p2.y;
double c = cross2(p1, p2);
return c > 0 || (c == 0 && fabs(p1.x) < fabs(p2.x));
}
};
static Pt pts[10005];
static Seg segs[30005];
static double interpolate(const Pt& p1, const Pt& p2, double& x) {
if (cmpZero(p1.x - p2.x) == 0) return min(p1.y, p2.y);
return p1.y + (double)(p2.y - p1.y) / (p2.x - p1.x) * (x - p1.x);
}
struct CMP {
static double x;
bool operator() (const Seg &i, const Seg &j) {
double l = interpolate(i.s, i.e, x);
double r = interpolate(j.s, j.e, x);
return l < r;
}
};
double CMP::x;
struct BSP {
struct Node {
double ly, ry;
void init() {
ly = BOX_MX, ry = -BOX_MX;
}
void extend(Node *u) {
if (u == NULL)
return ;
ly = min(ly, u->ly);
ry = max(ry, u->ry);
}
} bbox[262144];
vector<double> x;
set<Seg, CMP> tree[524288];
Seg segs[262144];
void clear(int k, int l, int r) {
tree[k].clear();
bbox[k].init();
if (l+1 >= r)
return ;
int m = (l+r)/2;
clear(k<<1, l, m);
clear(k<<1|1, m, r);
}
void insert(Seg &s, int k, int l, int r) {
if (s.s.x <= x[l] && x[r] <= s.e.x) {
CMP::x = (x[l] + x[r])/2;
tree[k].insert(s);
double ly = interpolate(s.s, s.e, x[l]);
double ry = interpolate(s.s, s.e, x[r]);
bbox[k].ly = min(bbox[k].ly, min(ly, ry));
bbox[k].ry = max(bbox[k].ry, max(ly, ry));
return ;
}
if (l+1 >= r)
return ;
int m = (l+r)/2;
if (s.s.x <= x[m]) {
insert(s, k<<1, l, m);
bbox[k].extend(&bbox[k<<1]);
}
if (s.e.x > x[m]) {
insert(s, k<<1|1, m, r);
bbox[k].extend(&bbox[k<<1|1]);
}
}
void build_tree(Seg s[], int m) {
memcpy(segs, s, sizeof(segs[0])*m);
x.clear();
for (int i = 0; i < m; i++)
x.push_back(segs[i].s.x), x.push_back(segs[i].e.x);
sort(x.begin(), x.end());
x.resize(unique(x.begin(), x.end()) - x.begin());
clear(1, 0, x.size()-1);
for (int i = 0; i < m; i++) {
if (segs[i].s.x > segs[i].e.x)
swap(segs[i].s, segs[i].e);
if (segs[i].s.x < segs[i].e.x)
insert(segs[i], 1, 0, x.size()-1);
}
}
Pt q_st, q_ed;
int q_si;
void search(int k, int l, int r) {
if (bbox[k].ly >= q_ed.y || bbox[k].ry <= q_st.y)
return ;
if (x[l] <= q_st.x && q_st.x <= x[r]) {
CMP::x = q_st.x;
auto it = tree[k].upper_bound(Seg(q_st, q_st, 0));
while (it != tree[k].end()) {
double y = interpolate(it->s, it->e, CMP::x);
if (y > q_st.y) {
if (y < q_ed.y) {
q_ed.y = y;
q_si = it->i;
}
break;
}
it++;
}
}
if (l+1 >= r)
return ;
int m = (l+r)/2;
if (q_st.x <= x[m])
search(k<<1, l, m);
if (q_st.x >= x[m])
search(k<<1|1, m, r);
}
pair<int, Pt> raycast(Pt st) {
q_st = st;
q_ed = Pt(st.x, BOX_MX+1);
q_si = -1;
search(1, 0, x.size()-1);
return {q_si, q_ed};
}
} tree;
struct Disjoint {
static const int MAXN = 65536;
uint16_t parent[MAXN], weight[MAXN];
void init(int n) {
if (n >= MAXN)
exit(0);
for (int i = 0; i <= n; i++)
parent[i] = i, weight[i] = 1;
}
int findp(int x) {
return parent[x] == x ? x : parent[x] = findp(parent[x]);
}
int joint(int x, int y) {
x = findp(x), y = findp(y);
if (weight[x] >= weight[y])
parent[y] = x, weight[x] += weight[y];
else
parent[x] = y, weight[y] += weight[x];
}
} egroup, sgroup;
int main() {
int n, m, p, q;
while (scanf("%d %d %d %d", &n, &m, &p, &q) == 4 && n) {
for (int i = 0; i < n; i++) {
int x, y;
scanf("%d %d", &x, &y);
pts[i] = Pt(x, y);
}
sgroup.init(n);
for (int i = 0; i < m; i++) {
int st_i, ed_i;
scanf("%d %d", &st_i, &ed_i);
st_i--, ed_i--;
segs[i] = Seg(pts[st_i], pts[ed_i], i);
sgroup.joint(st_i, ed_i);
}
segs[m] = Seg(Pt(BOX_MX, BOX_MX), Pt(BOX_MX, -BOX_MX), m), m++;
segs[m] = Seg(Pt(BOX_MX, -BOX_MX), Pt(-BOX_MX, -BOX_MX), m), m++;
segs[m] = Seg(Pt(-BOX_MX, -BOX_MX), Pt(-BOX_MX, BOX_MX), m), m++;
segs[m] = Seg(Pt(-BOX_MX, BOX_MX), Pt(BOX_MX, BOX_MX), m), m++;
static map<Pt, vector<pair<Pt, uint8_t>>> g; g.clear();
static vector<vector<Pt>> on_seg; on_seg.clear(); on_seg.resize(m);
for (int i = 0; i < m; i++) {
on_seg[segs[i].i].reserve(4);
on_seg[segs[i].i].push_back(segs[i].s);
on_seg[segs[i].i].push_back(segs[i].e);
}
// for (int i = 0; i < n; i++)
// pts[i].println();
// for (int i = 0; i < m; i++)
// segs[i].println();
tree.build_tree(segs, m);
{
Pt top[10005];
for (int i = 0; i < n; i++)
top[i] = pts[i];
for (int i = 0; i < n; i++) {
int gid = sgroup.findp(i);
if (top[gid].y < pts[i].y)
top[gid] = pts[i];
}
for (int i = 0; i < n; i++) {
if (sgroup.findp(i) != i)
continue;
auto p = top[i];
auto hit = tree.raycast(p);
if (hit.first >= 0) {
on_seg[hit.first].emplace_back(hit.second);
g[p].emplace_back(hit.second, 1);
g[hit.second].emplace_back(p, 1);
}
}
}
for (int i = 0; i < m; i++) {
vector<Pt> &a = on_seg[i];
sort(a.begin(), a.end());
a.resize(unique(a.begin(), a.end()) - a.begin());
auto *prev = &g[a[0]];
for (int j = 1; j < a.size(); j++) {
prev->emplace_back(a[j], 0);
prev = &g[a[j]];
prev->emplace_back(a[j-1], 0);
}
}
for (auto &e : g)
sort(e.second.begin(), e.second.end(), AngleCmp(e.first));
static map<Pt, map<Pt, int>> R; R.clear();
int Rsize = 0;
for (auto &e : g) {
int sz = e.second.size();
map<Pt, int> &Rg = R[e.first];
for (auto &f : e.second) {
int &eid = Rg[f.first];
if (eid == 0)
eid = ++Rsize;
}
}
egroup.init(Rsize);
for (auto &e : g) {
int sz = e.second.size();
map<Pt, int> &Rg = R[e.first];
for (int i = sz-1, j = 0; j < sz; i = j++) {
int l = R[e.second[i].first][e.first];
int r = Rg[e.second[j].first];
egroup.joint(l, r);
if (e.second[i].second != 0) {
r = Rg[e.second[i].first];
assert(l > 0 && r > 0);
egroup.joint(l, r);
}
}
}
for (auto &e : g) {
int sz = e.second.size();
int n = 0;
for (int i = 0; i < sz; i++) {
if (e.second[i].second == 0)
e.second[n++] = e.second[i];
}
e.second.resize(n);
}
static int region[65536]; memset(region, 0, sizeof(0)*Rsize);
for (int i = 0; i < p; i++) {
float x, y;
scanf("%f %f", &x, &y);
pair<int, Pt> hit = tree.raycast(Pt(x, y));
if (hit.first < 0)
continue;
if (g.find(hit.second) != g.end()) {
auto &adj = g[hit.second];
AngleCmp cmp(hit.second);
int pos = 0;
pair<Pt, int> q(Pt(x, y), i+1);
pos = lower_bound(adj.begin(), adj.end(), q, cmp) - adj.begin() - 1;
assert(pos >= 0);
int l = R[adj[pos].first][hit.second];
assert(l > 0);
l = egroup.findp(l);
region[l] = i+1;
} else {
auto &adj = on_seg[hit.first];
Pt lpt = adj[0], rpt = adj[1];
int l;
if (cmpZero(cross(lpt, rpt, Pt(x, y))) <= 0)
l = R[lpt][rpt];
else
l = R[rpt][lpt];
assert(l > 0);
l = egroup.findp(l);
region[l] = i+1;
}
}
for (int i = 0; i < q; i++) {
float x, y;
scanf("%f %f", &x, &y);
pair<int, Pt> hit = tree.raycast(Pt(x, y));
if (hit.first < 0) {
printf("0\n");
continue;
}
if (g.find(hit.second) != g.end()) {
auto &adj = g[hit.second];
AngleCmp cmp(hit.second);
int pos = 0;
pair<Pt, int> q(Pt(x, y), i+1);
pos = lower_bound(adj.begin(), adj.end(), q, cmp) - adj.begin() - 1;
assert(pos >= 0);
int l = R[adj[pos].first][hit.second];
assert(l > 0);
l = egroup.findp(l);
printf("%d\n", region[l]);
} else {
auto &adj = on_seg[hit.first];
Pt lpt = adj[0], rpt = adj[1];
int l;
if (cmpZero(cross(lpt, rpt, Pt(x, y))) <= 0)
l = R[lpt][rpt];
else
l = R[rpt][lpt];
assert(l > 0);
l = egroup.findp(l);
printf("%d\n", region[l]);
}
}
}
return 0;
}