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

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

Solution

Previously

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

Detail

• 相關問題 演算法筆記 - Polygon Partition / Polygon Trapezoidalization

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

處理的流程如下：

首先，對於幾何剖分的問題，為處理方便，都會套上一個無限大的外框，把所有頂點包在裡面。你可以透過 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

#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;
}

線段樹

#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;
}