Day 109-Q1 | Randomized Incremental Delaunay Triangulation

問題

2次元平面上に $N$ 個の点が与えられる。座標は相異なり、一般位置（同一直線上に3点が並ばず、同一円周上に4点が乗らない）が保証される。

これらの点集合の Delaunay三角形分割 を構築し、その三角形分割に含まれる辺の本数と、全ての辺の長さの総和を出力せよ。

Delaunay三角形分割とは、「どの三角形の外接円の内部にも、他のどの点も含まれない」という性質（外接円空性）を満たす三角形分割である。

入力形式

N
x_1 y_1
x_2 y_2
...
x_N y_N

制約

$3 \le N \le 500$

$0 \le x_i, y_i \le 10^4$（実数、小数点以下2桁まで）

点は一般位置（同一直線上3点なし、同一円周上4点なし）

入出力例

入力例1

出力例1

5
5.519201

4点は凸四角形$A(0,0),B(1,0),C(1,1),D(0,1.1)$をなす。対角線$AC$を選ぶ分割$\{ABC,ACD\}$が外接円空性を満たす（三角形$ABC$の外接円半径$\approx0.707$に対し$D$までの距離$\approx0.781$で外側、三角形$ACD$の外接円半径$\approx0.711$に対し$B$までの距離$\approx0.778$で外側）。辺は$AB,BC,CD,DA,AC$の5本、長さの総和は$1+1+\sqrt{1.01}+1.1+\sqrt2\approx5.519201$。

概念図: 点挿入とedge flip

ヒント（段階的開示）

ヒント1: 方向性

$N$個の点から作れる三角形分割の候補は複数あるが、外接円空性（どの三角形の外接円の内部にも他の点がない）を満たすものを見つけたい。全ての三角形の組み合わせについて外接円空性を愚直にチェックする方法は$O(N^4)$かかり、$N$が大きくなると現実的ではない。三角形分割を1つずつ改善しながら外接円空性へ近づけていく、という発想が必要になる。

ヒント2: アプローチ

点をランダムな順序で1つずつ挿入していく乱択増分アルゴリズムを使う。

1. 全ての点を内部に含む十分大きな「super triangle（番兵三角形）」を用意し、そこから三角形分割を始める。
2. 点をランダム順に1つずつ挿入する。挿入する点$p$を含む三角形を見つけ、その三角形を$p$を頂点とする3つの三角形に分割する。
3. 分割によって新しくできた辺の「向かい側」にある隣接三角形との間で、外接円空性が崩れていないかを調べる。崩れていれば、その2つの三角形が共有する辺（対角線）をflip（架け替え）して局所的に外接円空性を回復する。
4. 全点を挿入し終えたら、super triangleの頂点を含む三角形を取り除く。

ランダムな挿入順序を使うことで、期待計算量が悪化しにくいことが理論的に保証される。

ヒント3: 誘導（コード骨格）

def in_circle(a, b, c, d):
    # a,b,c が反時計回りのとき、dがa,b,cの外接円の内部にあるか
    ax, ay = a[0]-d[0], a[1]-d[1]
    bx, by = b[0]-d[0], b[1]-d[1]
    cx, cy = c[0]-d[0], c[1]-d[1]
    det = ((ax*ax+ay*ay)*(bx*cy-by*cx)
           - (bx*bx+by*by)*(ax*cy-ay*cx)
           + (cx*cx+cy*cy)*(ax*by-ay*bx))
    return det > 0  # 正なら内部

# 隣接する2三角形 (u,v,w) と (u,v,x) が辺(u,v)を共有しているとき、
# xがtriangle(u,v,w)の外接円内部にあるなら、辺(u,v)をflipして
# (w,x)を新しい対角線にすると外接円空性が改善する。

模範解答 (Python)

import sys
import random

def solve():
    data = sys.stdin.read().split()
    idx = 0
    n = int(data[idx]); idx += 1
    pts = []
    for _ in range(n):
        x = float(data[idx]); idx += 1
        y = float(data[idx]); idx += 1
        pts.append((x, y))

    minx = min(p[0] for p in pts); maxx = max(p[0] for p in pts)
    miny = min(p[1] for p in pts); maxy = max(p[1] for p in pts)
    d = max(maxx - minx, maxy - miny, 1.0) * 20
    midx, midy = (minx + maxx) / 2, (miny + maxy) / 2
    super_pts = [
        (midx - d, midy - d),
        (midx + d, midy - d),
        (midx, midy + d),
    ]
    P = pts + super_pts
    N = len(pts)
    SA, SB, SC = N, N + 1, N + 2

    def cross(o, a, b):
        return (a[0]-o[0])*(b[1]-o[1]) - (a[1]-o[1])*(b[0]-o[0])

    def ccw(tri):
        a, b, c = tri
        return tri if cross(P[a], P[b], P[c]) > 0 else (a, c, b)

    def in_circle(tri, d_idx):
        a, b, c = tri
        ax, ay = P[a][0]-P[d_idx][0], P[a][1]-P[d_idx][1]
        bx, by = P[b][0]-P[d_idx][0], P[b][1]-P[d_idx][1]
        cx, cy = P[c][0]-P[d_idx][0], P[c][1]-P[d_idx][1]
        det = ((ax*ax+ay*ay)*(bx*cy-by*cx)
               - (bx*bx+by*by)*(ax*cy-ay*cx)
               + (cx*cx+cy*cy)*(ax*by-ay*bx))
        return det > 1e-9

    def in_triangle(p, tri):
        a, b, c = tri
        d1 = cross(P[a], P[b], p)
        d2 = cross(P[b], P[c], p)
        d3 = cross(P[c], P[a], p)
        neg = d1 < -1e-9 or d2 < -1e-9 or d3 < -1e-9
        pos = d1 > 1e-9 or d2 > 1e-9 or d3 > 1e-9
        return not (neg and pos)

    triangles = [ccw((SA, SB, SC))]
    order = list(range(N))
    random.shuffle(order)

    for pi in order:
        p = P[pi]
        ti = next(i for i, t in enumerate(triangles) if in_triangle(p, t))
        a, b, c = triangles.pop(ti)
        triangles.extend([ccw((a, b, pi)), ccw((b, c, pi)), ccw((c, a, pi))])

        changed = True
        while changed:
            changed = False
            for i in range(len(triangles)):
                for j in range(i + 1, len(triangles)):
                    t1, t2 = triangles[i], triangles[j]
                    shared = set(t1) & set(t2)
                    if len(shared) != 2:
                        continue
                    u, v = tuple(shared)
                    w = [x for x in t1 if x not in shared][0]
                    x = [x for x in t2 if x not in shared][0]
                    if in_circle(t1, x):
                        triangles[i] = ccw((u, w, x))
                        triangles[j] = ccw((v, w, x))
                        changed = True

    result = [t for t in triangles if SA not in t and SB not in t and SC not in t]

    edges = set()
    for a, b, c in result:
        for u, v in [(a, b), (b, c), (c, a)]:
            edges.add((min(u, v), max(u, v)))

    total_len = sum(((P[u][0]-P[v][0])**2 + (P[u][1]-P[v][1])**2) ** 0.5 for u, v in edges)

    print(len(edges))
    print(f"{total_len:.6f}")

solve()

計算量: 理論的な乱択増分Delaunay構築は期待$O(N \log N)$（点位置探索をDAGで管理する場合）。本解答は簡略実装（最悪$O(N^3)$程度）だが$N\le500$では十分高速。30回のランダム点配置で$O(N^4)$の愚直な外接円空性チェックと突き合わせるストレステストで一致を確認済み。

Step-by-Step 解説

1super triangleで境界処理を回避する
全ての入力点を内部に含む十分大きな仮想三角形（super triangle）を最初に用意し、実在しない3頂点SA, SB, SCを持つ1つの三角形からスタートする。

2ランダム順に点を挿入し、含む三角形を3分割する
挿入する点$p$を含む三角形を線形探索で見つけ（in_triangle）、その三角形を$p$を頂点とする3つの小さな三角形に置き換える。

3in_circle判定でDelaunay性を局所的に回復する（flip）
隣接三角形の「向かいの頂点」が自分の外接円の内部にあれば、共有辺をflip（対角線を架け替える）ことで解消する。これを全ての辺について崩れがなくなるまで繰り返す（Lawsonのflipアルゴリズム）。

4super triangleの頂点を含む三角形を除去する
全ての実点を挿入し終えたら、番兵に接している三角形を結果から取り除く。

5辺集合を復元し、本数と総延長を計算する
残った各三角形から3本の辺を取り出し、重複を除いた辺の集合を作り、本数と長さの総和を出力する。

よくあるミス

ミス	原因	正しい書き方
`in_circle`の外接円行列式の符号を逆にしてしまう	三角形の頂点順序が反時計回り(CCW)であることを前提とした式なのに、頂点順序を正規化せずに使う	三角形を必ず`ccw()`で正規化してから`in_circle`を呼ぶ
点が三角形の辺上にちょうど乗るケースの扱いを誤る	`in_triangle`の3つのcross積が1つだけ0になる境界ケースの判定漏れ	`neg`と`pos`の両方が真になったときだけ「外」と判定するロジックにする
浮動小数点の座標をそのまま厳密比較し、flipが無限ループする	数値誤差でわずかに外接円境界上にある点の判定が振動する	`in_circle`にイプシロン(`1e-9`程度)の許容誤差を設ける
super triangleを小さく作ってしまい、一部の入力点を内部に含められない	バウンディングボックスの何倍にするかの係数が小さすぎる	幅・高さの最大値の10〜20倍程度余裕を持たせる

次のステップ

発展: 点位置探索を「履歴DAG」に置き換え、真の期待$O(N \log N)$を実現する
発展: Delaunay三角形分割の双対グラフを取るとVoronoi図が得られる。最近傍点探索クエリへの応用を考える
発展: 3次元に拡張し、4面体分割（Delaunay tetrahedralization）として点群補間に応用する

自己評価

完全理解おおむね理解要復習

Day 109-Q1 — Randomized Incremental Delaunay Triangulation

問題