$O(log n)$ 시간에 점 위치 찾기¶
다음 문제를 생각해 봅시다. 차수가 1이거나 0인 정점이 없는 평면 분할과 많은 수의 쿼리가 주어집니다. 각 쿼리는 하나의 점이며, 그 점이 속하는 분할의 면을 구해야 합니다. 각 쿼리에 오프라인으로 $O(\log n)$에 답하겠습니다.
이 문제는 Voronoi 다이어그램이나 어떤 단순 다각형에서 몇몇 점의 위치를 찾아야 할 때 발생할 수 있습니다.
알고리즘¶
먼저, 각 쿼리 점 $p\ (x_0, y_0)$에 대해 다음과 같은 간선을 찾고자 합니다. 점이 어떤 간선에 속한다면 그 점은 우리가 찾은 간선 위에 놓여야 하고, 그렇지 않다면 이 간선은 직선 $x = x_0$과 유일한 점 $(x_0, y)$에서 교차해야 하며, 이때 $y < y_0$이고 이 $y$는 그러한 모든 간선 중에서 최대여야 합니다. 다음 그림은 두 경우를 모두 보여 줍니다.
이 문제는 스위프 라인 알고리즘을 사용하여 오프라인으로 풀 것입니다. 쿼리 점과 간선 끝점의 x좌표를 오름차순으로 순회하면서 간선의 집합 $s$를 유지합시다. 각 x좌표에 대해 미리 몇 가지 이벤트를 추가해 둘 것입니다.
이벤트는 add, remove, vertical, get의 네 가지 유형입니다. 각 수직 간선(두 끝점의 x좌표가 같은 간선)에 대해 해당 x좌표에 vertical 이벤트를 하나 추가합니다. 그 밖의 모든 간선에 대해서는 끝점 x좌표의 최솟값에 add 이벤트를 하나, 끝점 x좌표의 최댓값에 remove 이벤트를 하나 추가합니다. 마지막으로, 각 쿼리 점에 대해 그 x좌표에 get 이벤트를 하나 추가합니다.
각 x좌표에 대해 이벤트를 유형에 따라 (vertical, get, remove, add) 순서로 정렬합니다. 다음 그림은 각 x좌표에 대해 정렬된 순서의 모든 이벤트를 보여 줍니다.
스위프 라인 과정 동안 두 개의 집합을 유지합니다. 하나는 수직이 아닌 모든 간선을 위한 집합 $t$이고, 다른 하나는 수직 간선만을 위한 집합 $vert$입니다. 각 x좌표의 처리를 시작할 때마다 집합 $vert$를 비웁니다.
이제 고정된 x좌표에 대한 이벤트를 처리해 봅시다.
- vertical 이벤트가 들어오면, 해당 간선 끝점의 y좌표 중 최솟값을 $vert$에 삽입하기만 하면 됩니다.
- remove 또는 add 이벤트가 들어오면, 해당 간선을 $t$에서 제거하거나 $t$에 추가합니다.
- 마지막으로, 각 get 이벤트에 대해 $vert$에서 이진 탐색을 수행하여 점이 어떤 수직 간선 위에 있는지 확인해야 합니다. 점이 어떤 수직 간선 위에도 놓이지 않는다면, 이 쿼리의 답을 $t$에서 찾아야 합니다. 이를 위해 다시 이진 탐색을 수행합니다. 몇몇 퇴화 경우(예를 들어 삼각형 $(0,~0)$, $(0,~2)$, $(1, 1)$에서 점 $(0,~0)$을 쿼리하는 경우)를 처리하기 위해, 이 x좌표의 모든 이벤트를 처리한 후 모든 get 이벤트에 다시 답하고 두 답 중 더 나은 것을 선택해야 합니다.
이제 집합 $t$를 위한 비교자를 정해 봅시다. 이 비교자는 두 간선이 공통으로 덮는 모든 x좌표에서 한 간선이 다른 간선보다 위에 놓이지 않는지를 확인해야 합니다. 두 간선 $(a, b)$와 $(c, d)$가 있다고 가정합시다. 그러면 비교자는 다음과 같습니다(의사 코드):
$val = sgn((b - a)\times(c - a)) + sgn((b - a)\times(d - a))$
if $val \neq 0$
then return $val > 0$
$val = sgn((d - c)\times(a - c)) + sgn((d - c)\times(b - c))$
return $val < 0$
이제 모든 쿼리에 대해 대응하는 간선을 구했습니다. 면은 어떻게 찾을 수 있습니까? 간선을 찾지 못했다면 점이 외부 면에 있다는 뜻입니다. 점이 찾은 간선 위에 있다면 면이 유일하게 정해지지 않습니다. 그렇지 않다면 후보는 두 개이며, 이 간선을 경계로 하는 두 면입니다. 어느 쪽이 답인지 어떻게 확인할 수 있습니까? 이 간선은 수직이 아니라는 점에 주목하십시오. 그러면 답은 이 간선의 위쪽에 있는 면입니다. 수직이 아닌 각 간선에 대해 이러한 면을 찾아봅시다. 각 면을 반시계 방향으로 순회한다고 생각해 봅시다. 이 순회 중에 간선을 지나면서 x좌표가 증가했다면, 이 면이 바로 이 간선에 대해 찾아야 하는 면입니다.
참고¶
사실 퍼시스턴트 트리(persistent tree)를 사용하면 이 방법으로 쿼리에 온라인으로 답할 수 있습니다.
구현¶
다음 코드는 정수에 대해 구현되어 있지만, 비교 메서드와 점 타입을 바꾸면 double에서도 동작하도록 쉽게 수정할 수 있습니다. 이 구현은 분할이 DCEL 안에 올바르게 저장되어 있고 외부 면의 번호가 $-1$이라고 가정합니다.
각 쿼리에 대해 점이 $i$번 면의 내부에 엄격하게 놓여 있으면 쌍 $(1, i)$가 반환되고, 점이 $i$번 간선 위에 놓여 있으면 쌍 $(0, i)$가 반환됩니다.
typedef long long ll;
bool ge(const ll& a, const ll& b) { return a >= b; }
bool le(const ll& a, const ll& b) { return a <= b; }
bool eq(const ll& a, const ll& b) { return a == b; }
bool gt(const ll& a, const ll& b) { return a > b; }
bool lt(const ll& a, const ll& b) { return a < b; }
int sgn(const ll& x) { return le(x, 0) ? eq(x, 0) ? 0 : -1 : 1; }
struct pt {
ll x, y;
pt() {}
pt(ll _x, ll _y) : x(_x), y(_y) {}
pt operator-(const pt& a) const { return pt(x - a.x, y - a.y); }
ll dot(const pt& a) const { return x * a.x + y * a.y; }
ll dot(const pt& a, const pt& b) const { return (a - *this).dot(b - *this); }
ll cross(const pt& a) const { return x * a.y - y * a.x; }
ll cross(const pt& a, const pt& b) const { return (a - *this).cross(b - *this); }
bool operator==(const pt& a) const { return a.x == x && a.y == y; }
};
struct Edge {
pt l, r;
};
bool edge_cmp(Edge* edge1, Edge* edge2)
{
const pt a = edge1->l, b = edge1->r;
const pt c = edge2->l, d = edge2->r;
int val = sgn(a.cross(b, c)) + sgn(a.cross(b, d));
if (val != 0)
return val > 0;
val = sgn(c.cross(d, a)) + sgn(c.cross(d, b));
return val < 0;
}
enum EventType { DEL = 2, ADD = 3, GET = 1, VERT = 0 };
struct Event {
EventType type;
int pos;
bool operator<(const Event& event) const { return type < event.type; }
};
vector<Edge*> sweepline(vector<Edge*> planar, vector<pt> queries)
{
using pt_type = decltype(pt::x);
// collect all x-coordinates
auto s =
set<pt_type, std::function<bool(const pt_type&, const pt_type&)>>(lt);
for (pt p : queries)
s.insert(p.x);
for (Edge* e : planar) {
s.insert(e->l.x);
s.insert(e->r.x);
}
// map all x-coordinates to ids
int cid = 0;
auto id =
map<pt_type, int, std::function<bool(const pt_type&, const pt_type&)>>(
lt);
for (auto x : s)
id[x] = cid++;
// create events
auto t = set<Edge*, decltype(*edge_cmp)>(edge_cmp);
auto vert_cmp = [](const pair<pt_type, int>& l,
const pair<pt_type, int>& r) {
if (!eq(l.first, r.first))
return lt(l.first, r.first);
return l.second < r.second;
};
auto vert = set<pair<pt_type, int>, decltype(vert_cmp)>(vert_cmp);
vector<vector<Event>> events(cid);
for (int i = 0; i < (int)queries.size(); i++) {
int x = id[queries[i].x];
events[x].push_back(Event{GET, i});
}
for (int i = 0; i < (int)planar.size(); i++) {
int lx = id[planar[i]->l.x], rx = id[planar[i]->r.x];
if (lx > rx) {
swap(lx, rx);
swap(planar[i]->l, planar[i]->r);
}
if (lx == rx) {
events[lx].push_back(Event{VERT, i});
} else {
events[lx].push_back(Event{ADD, i});
events[rx].push_back(Event{DEL, i});
}
}
// perform sweep line algorithm
vector<Edge*> ans(queries.size(), nullptr);
for (int x = 0; x < cid; x++) {
sort(events[x].begin(), events[x].end());
vert.clear();
for (Event event : events[x]) {
if (event.type == DEL) {
t.erase(planar[event.pos]);
}
if (event.type == VERT) {
vert.insert(make_pair(
min(planar[event.pos]->l.y, planar[event.pos]->r.y),
event.pos));
}
if (event.type == ADD) {
t.insert(planar[event.pos]);
}
if (event.type == GET) {
auto jt = vert.upper_bound(
make_pair(queries[event.pos].y, planar.size()));
if (jt != vert.begin()) {
--jt;
int i = jt->second;
if (ge(max(planar[i]->l.y, planar[i]->r.y),
queries[event.pos].y)) {
ans[event.pos] = planar[i];
continue;
}
}
Edge* e = new Edge;
e->l = e->r = queries[event.pos];
auto it = t.upper_bound(e);
if (it != t.begin())
ans[event.pos] = *(--it);
delete e;
}
}
for (Event event : events[x]) {
if (event.type != GET)
continue;
if (ans[event.pos] != nullptr &&
eq(ans[event.pos]->l.x, ans[event.pos]->r.x))
continue;
Edge* e = new Edge;
e->l = e->r = queries[event.pos];
auto it = t.upper_bound(e);
delete e;
if (it == t.begin())
e = nullptr;
else
e = *(--it);
if (ans[event.pos] == nullptr) {
ans[event.pos] = e;
continue;
}
if (e == nullptr)
continue;
if (e == ans[event.pos])
continue;
if (id[ans[event.pos]->r.x] == x) {
if (id[e->l.x] == x) {
if (gt(e->l.y, ans[event.pos]->r.y))
ans[event.pos] = e;
}
} else {
ans[event.pos] = e;
}
}
}
return ans;
}
struct DCEL {
struct Edge {
pt origin;
Edge* nxt = nullptr;
Edge* twin = nullptr;
int face;
};
vector<Edge*> body;
};
vector<pair<int, int>> point_location(DCEL planar, vector<pt> queries)
{
vector<pair<int, int>> ans(queries.size());
vector<Edge*> planar2;
map<intptr_t, int> pos;
map<intptr_t, int> added_on;
int n = planar.body.size();
for (int i = 0; i < n; i++) {
if (planar.body[i]->face > planar.body[i]->twin->face)
continue;
Edge* e = new Edge;
e->l = planar.body[i]->origin;
e->r = planar.body[i]->twin->origin;
added_on[(intptr_t)e] = i;
pos[(intptr_t)e] =
lt(planar.body[i]->origin.x, planar.body[i]->twin->origin.x)
? planar.body[i]->face
: planar.body[i]->twin->face;
planar2.push_back(e);
}
auto res = sweepline(planar2, queries);
for (int i = 0; i < (int)queries.size(); i++) {
if (res[i] == nullptr) {
ans[i] = make_pair(1, -1);
continue;
}
pt p = queries[i];
pt l = res[i]->l, r = res[i]->r;
if (eq(p.cross(l, r), 0) && le(p.dot(l, r), 0)) {
ans[i] = make_pair(0, added_on[(intptr_t)res[i]]);
continue;
}
ans[i] = make_pair(1, pos[(intptr_t)res[i]]);
}
for (auto e : planar2)
delete e;
return ans;
}