Delaunay 삼각분할과 Voronoi 다이어그램¶
평면 위의 점들의 집합 $\{p_i\}$를 생각해 봅시다. $\{p_i\}$의 Voronoi 다이어그램 $V(\{p_i\})$는 평면을 $n$개의 영역 $V_i$로 분할한 것이며, 여기서 $V_i = \{p\in\mathbb{R}^2;\ \rho(p, p_i) = \min\ \rho(p, p_k)\}$입니다. Voronoi 다이어그램의 셀은 (무한할 수도 있는) 다각형입니다. $\{p_i\}$의 Delaunay 삼각분할 $D(\{p_i\})$는 모든 점 $p_i$가 각 삼각형 $T \in D(\{p_i\})$의 외접원 바깥 또는 경계 위에 있는 삼각분할입니다.
Voronoi 다이어그램이 연결되어 있지 않고 Delaunay 삼각분할이 존재하지 않는 까다로운 퇴화 경우가 있습니다. 이는 모든 점이 한 직선 위에 있는 경우입니다.
성질¶
Delaunay 삼각분할은 가능한 모든 삼각분할 중에서 최소 각을 최대화합니다.
점 집합의 최소 유클리드 신장 트리는 그 점 집합의 Delaunay 삼각분할의 간선들의 부분집합입니다.
쌍대성¶
$\{p_i\}$가 한 직선 위에 있지 않고, $\{p_i\}$ 중 어떤 네 점도 한 원 위에 있지 않다고 가정합시다. 그러면 $V(\{p_i\})$와 $D(\{p_i\})$는 쌍대이므로, 둘 중 하나를 구하면 나머지 하나를 $O(n)$에 구할 수 있습니다. 이 조건을 만족하지 않으면 어떻게 해야 합니까? 모든 점이 한 직선 위에 있는 경우는 쉽게 처리할 수 있습니다. 그렇지 않은 경우에는 $V$와 $D'$가 쌍대이며, 여기서 $D'$는 $D$에서 간선을 사이에 둔 두 삼각형이 같은 외접원을 공유하는 간선을 모두 제거하여 얻은 것입니다.
Delaunay 삼각분할과 Voronoi 다이어그램 구하기¶
쌍대성 덕분에 $V$와 $D$ 중 하나만 계산하는 빠른 알고리즘만 있으면 됩니다. $D(\{p_i\})$를 $O(n\log n)$에 구하는 방법을 설명하겠습니다. 삼각분할은 Guibas와 Stolfi의 분할 정복 알고리즘으로 구합니다.
쿼드 간선(quad-edge) 자료 구조¶
알고리즘을 수행하는 동안 $D$는 쿼드 간선 자료 구조에 저장됩니다. 이 자료 구조는 그림에 설명되어 있습니다:
알고리즘에서는 간선에 대해 다음 함수들을 사용합니다:
make_edge(a, b)
이 함수는 점a에서 점b로 가는 고립된 간선을 그 역방향 간선 및 두 쌍대 간선과 함께 생성합니다.splice(a, b)
이 함수는 알고리즘의 핵심 함수입니다. 이 함수는a->Onext와b->Onext를 서로 바꾸고,a->Onext->Rot->Onext와b->Onext->Rot->Onext를 서로 바꿉니다.delete_edge(e)
이 함수는 삼각분할에서 e를 삭제합니다.e를 삭제하려면splice(e, e->Oprev)와splice(e->Rev, e->Rev->Oprev)를 호출하기만 하면 됩니다.connect(a, b)
이 함수는a,b,e가 모두 같은 왼쪽 면을 갖도록a->Dest에서b->Org로 가는 새 간선e를 생성합니다. 이를 위해e = make_edge(a->Dest, b->Org),splice(e, a->Lnext),splice(e->Rev, b)를 호출합니다.
알고리즘¶
이 알고리즘은 삼각분할을 계산하여 두 개의 quad-edge를 반환합니다. 하나는 가장 왼쪽 정점에서 나가는 반시계 방향 볼록 껍질 간선이고, 다른 하나는 가장 오른쪽 정점에서 나가는 시계 방향 볼록 껍질 간선입니다.
모든 점을 x 기준으로 정렬하고, $x_1 = x_2$이면 y 기준으로 정렬합시다. 어떤 구간 $(l, r)$에 대해 문제를 풀어 봅시다(처음에는 $(l, r) = (0, n - 1)$입니다). $r - l + 1 = 2$이면 간선 $(p[l], p[r])$를 추가하고 반환합니다. $r - l + 1 = 3$이면 먼저 간선 $(p[l], p[l + 1])$과 $(p[l + 1], p[r])$를 추가합니다. 또한 splice(a->Rev, b)를 사용하여 이 간선들을 연결해야 합니다. 이제 삼각형을 닫아야 합니다. 다음 동작은 $p[l], p[l + 1], p[r]$의 방향에 따라 달라집니다. 세 점이 한 직선 위에 있으면 삼각형을 만들 수 없으므로 단순히 (a, b->Rev)를 반환합니다. 그렇지 않으면 connect(b, a)를 호출하여 새 간선 c를 만듭니다. 점들이 반시계 방향으로 배치되어 있으면 (a, b->Rev)를 반환합니다. 그렇지 않으면 (c->Rev, c)를 반환합니다.
이제 $r - l + 1 \ge 4$라고 가정합시다. 먼저 $L = (l, \frac{l + r}{2})$와 $R = (\frac{l + r}{2} + 1, r)$를 재귀적으로 풉시다. 이제 이 삼각분할들을 하나의 삼각분할로 병합해야 합니다. 점들이 정렬되어 있으므로, 병합하는 동안 L에서 R로 가는 간선(이른바 교차 간선)을 추가하고 L에서 L로, R에서 R로 가는 간선 일부를 제거한다는 점에 유의하십시오. 교차 간선은 어떤 구조를 가집니까? 이 간선들은 모두 분할 기준 x 값에 놓인, y축에 평행한 직선과 교차해야 합니다. 이로써 교차 간선들 사이에 선형 순서가 정해지므로, 연속한 교차 간선, 가장 아래에 있는 교차 간선 등을 논할 수 있습니다. 알고리즘은 교차 간선을 오름차순으로 추가합니다. 인접한 두 교차 간선은 항상 공통 끝점을 가지며, 이들이 정의하는 삼각형의 세 번째 변은 L에서 L로, 또는 R에서 R로 이어진다는 점에 유의하십시오. 현재 교차 간선을 base라고 부릅시다. base의 다음 간선은 base의 왼쪽 끝점에서 오른쪽 끝점의 R-이웃 중 하나로 가거나, 그 반대로 갑니다. base와 이전 교차 간선의 외접원을 생각해 봅시다. 이 원이 base를 현으로 가지면서 Oy 방향으로 더 멀리 놓인 다른 원들로 변형된다고 가정합시다. 원은 한동안 위로 올라가지만, base가 L과 R의 위쪽 공통접선이 아닌 한 L 또는 R에 속한 점을 만나게 되며, 이 점은 외접원 안에 아무 점도 없는 새로운 삼각형을 만들어 냅니다. 이 삼각형의 새로운 L-R 간선이 다음으로 추가되는 교차 간선입니다. 이를 효율적으로 수행하기 위해 두 간선 lcand와 rcand를 계산하는데, lcand는 이 과정에서 처음 만나는 L의 점을, rcand는 처음 만나는 R의 점을 가리키도록 합니다. 그런 다음 먼저 만나게 될 쪽을 선택합니다. 처음에 base는 L과 R의 아래쪽 공통접선을 가리킵니다.
구현¶
in_circle 함수의 구현은 GCC 전용이라는 점에 유의하십시오.
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& a) { return a >= 0 ? a ? 1 : 0 : -1; }
struct pt {
ll x, y;
pt() { }
pt(ll _x, ll _y) : x(_x), y(_y) { }
pt operator-(const pt& p) const {
return pt(x - p.x, y - p.y);
}
ll cross(const pt& p) const {
return x * p.y - y * p.x;
}
ll cross(const pt& a, const pt& b) const {
return (a - *this).cross(b - *this);
}
ll dot(const pt& p) const {
return x * p.x + y * p.y;
}
ll dot(const pt& a, const pt& b) const {
return (a - *this).dot(b - *this);
}
ll sqrLength() const {
return this->dot(*this);
}
bool operator==(const pt& p) const {
return eq(x, p.x) && eq(y, p.y);
}
};
const pt inf_pt = pt(1e18, 1e18);
struct QuadEdge {
pt origin;
QuadEdge* rot = nullptr;
QuadEdge* onext = nullptr;
bool used = false;
QuadEdge* rev() const {
return rot->rot;
}
QuadEdge* lnext() const {
return rot->rev()->onext->rot;
}
QuadEdge* oprev() const {
return rot->onext->rot;
}
pt dest() const {
return rev()->origin;
}
};
QuadEdge* make_edge(pt from, pt to) {
QuadEdge* e1 = new QuadEdge;
QuadEdge* e2 = new QuadEdge;
QuadEdge* e3 = new QuadEdge;
QuadEdge* e4 = new QuadEdge;
e1->origin = from;
e2->origin = to;
e3->origin = e4->origin = inf_pt;
e1->rot = e3;
e2->rot = e4;
e3->rot = e2;
e4->rot = e1;
e1->onext = e1;
e2->onext = e2;
e3->onext = e4;
e4->onext = e3;
return e1;
}
void splice(QuadEdge* a, QuadEdge* b) {
swap(a->onext->rot->onext, b->onext->rot->onext);
swap(a->onext, b->onext);
}
void delete_edge(QuadEdge* e) {
splice(e, e->oprev());
splice(e->rev(), e->rev()->oprev());
delete e->rev()->rot;
delete e->rev();
delete e->rot;
delete e;
}
QuadEdge* connect(QuadEdge* a, QuadEdge* b) {
QuadEdge* e = make_edge(a->dest(), b->origin);
splice(e, a->lnext());
splice(e->rev(), b);
return e;
}
bool left_of(pt p, QuadEdge* e) {
return gt(p.cross(e->origin, e->dest()), 0);
}
bool right_of(pt p, QuadEdge* e) {
return lt(p.cross(e->origin, e->dest()), 0);
}
template <class T>
T det3(T a1, T a2, T a3, T b1, T b2, T b3, T c1, T c2, T c3) {
return a1 * (b2 * c3 - c2 * b3) - a2 * (b1 * c3 - c1 * b3) +
a3 * (b1 * c2 - c1 * b2);
}
bool in_circle(pt a, pt b, pt c, pt d) {
// If there is __int128, calculate directly.
// Otherwise, calculate angles.
#if defined(__LP64__) || defined(_WIN64)
__int128 det = -det3<__int128>(b.x, b.y, b.sqrLength(), c.x, c.y,
c.sqrLength(), d.x, d.y, d.sqrLength());
det += det3<__int128>(a.x, a.y, a.sqrLength(), c.x, c.y, c.sqrLength(), d.x,
d.y, d.sqrLength());
det -= det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), d.x,
d.y, d.sqrLength());
det += det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), c.x,
c.y, c.sqrLength());
return det > 0;
#else
auto ang = [](pt l, pt mid, pt r) {
ll x = mid.dot(l, r);
ll y = mid.cross(l, r);
long double res = atan2((long double)x, (long double)y);
return res;
};
long double kek = ang(a, b, c) + ang(c, d, a) - ang(b, c, d) - ang(d, a, b);
if (kek > 1e-8)
return true;
else
return false;
#endif
}
pair<QuadEdge*, QuadEdge*> build_tr(int l, int r, vector<pt>& p) {
if (r - l + 1 == 2) {
QuadEdge* res = make_edge(p[l], p[r]);
return make_pair(res, res->rev());
}
if (r - l + 1 == 3) {
QuadEdge *a = make_edge(p[l], p[l + 1]), *b = make_edge(p[l + 1], p[r]);
splice(a->rev(), b);
int sg = sgn(p[l].cross(p[l + 1], p[r]));
if (sg == 0)
return make_pair(a, b->rev());
QuadEdge* c = connect(b, a);
if (sg == 1)
return make_pair(a, b->rev());
else
return make_pair(c->rev(), c);
}
int mid = (l + r) / 2;
QuadEdge *ldo, *ldi, *rdo, *rdi;
tie(ldo, ldi) = build_tr(l, mid, p);
tie(rdi, rdo) = build_tr(mid + 1, r, p);
while (true) {
if (left_of(rdi->origin, ldi)) {
ldi = ldi->lnext();
continue;
}
if (right_of(ldi->origin, rdi)) {
rdi = rdi->rev()->onext;
continue;
}
break;
}
QuadEdge* basel = connect(rdi->rev(), ldi);
auto valid = [&basel](QuadEdge* e) { return right_of(e->dest(), basel); };
if (ldi->origin == ldo->origin)
ldo = basel->rev();
if (rdi->origin == rdo->origin)
rdo = basel;
while (true) {
QuadEdge* lcand = basel->rev()->onext;
if (valid(lcand)) {
while (in_circle(basel->dest(), basel->origin, lcand->dest(),
lcand->onext->dest())) {
QuadEdge* t = lcand->onext;
delete_edge(lcand);
lcand = t;
}
}
QuadEdge* rcand = basel->oprev();
if (valid(rcand)) {
while (in_circle(basel->dest(), basel->origin, rcand->dest(),
rcand->oprev()->dest())) {
QuadEdge* t = rcand->oprev();
delete_edge(rcand);
rcand = t;
}
}
if (!valid(lcand) && !valid(rcand))
break;
if (!valid(lcand) ||
(valid(rcand) && in_circle(lcand->dest(), lcand->origin,
rcand->origin, rcand->dest())))
basel = connect(rcand, basel->rev());
else
basel = connect(basel->rev(), lcand->rev());
}
return make_pair(ldo, rdo);
}
vector<tuple<pt, pt, pt>> delaunay(vector<pt> p) {
sort(p.begin(), p.end(), [](const pt& a, const pt& b) {
return lt(a.x, b.x) || (eq(a.x, b.x) && lt(a.y, b.y));
});
auto res = build_tr(0, (int)p.size() - 1, p);
QuadEdge* e = res.first;
vector<QuadEdge*> edges = {e};
while (lt(e->onext->dest().cross(e->dest(), e->origin), 0))
e = e->onext;
auto add = [&p, &e, &edges]() {
QuadEdge* curr = e;
do {
curr->used = true;
p.push_back(curr->origin);
edges.push_back(curr->rev());
curr = curr->lnext();
} while (curr != e);
};
add();
p.clear();
int kek = 0;
while (kek < (int)edges.size()) {
if (!(e = edges[kek++])->used)
add();
}
vector<tuple<pt, pt, pt>> ans;
for (int i = 0; i < (int)p.size(); i += 3) {
ans.push_back(make_tuple(p[i], p[i + 1], p[i + 2]));
}
return ans;
}