점 구름을 삼각형으로 덮는 두 철학 — Bowyer

대각선은 외접원 부호 한 번으로 뒤집히고, 프론트의 다음 삼각형은 형상과 크기의 곱이 고른다.

1934년, 보리스 델로네는 스승 보로노이의 다이어그램에서 쌍대(dual) 하나를 꺼냈다. 점을 삼각형으로 잇되, 어느 외접원 안에도 다른 점이 없게. 90년이 지났다. 격자 생성기는 아직도 이 규칙으로 돈다.

시드 3으로 뽑은 점 12개는 삼각형 17개가 된다. 그 출력을 읽기 전에, 네 점이 사각형일 때 대각선이 어느 쪽으로 가는지부터 손보자. 위쪽 정점을 위아래로 끌면, 공유 에지를 넘는 순간 대각선이 뒤집힌다.

Drag vertex T through the shared edge. When T enters the circumcircle of the lower triangle, the diagonal snaps from LR to BT — a single Lawson flip restoring the empty-circumcircle rule.

정점 T가 아래 삼각형의 외접원 안으로 들어가는 위치를 보자. 대각선이 LR에서 BT로 바뀌는 그 한 번이 로슨 플립(Lawson flip)이다. 빈 외접원을 국소적으로 되돌리는 최소 연산이다.

점 구름을 삼각형으로 덮는 두 철학#

같은 점 집합의 삼각화는 끝이 없다. CFD가 피하는 것은 가늘고 뾰족한 삼각형이다. 납작하면 수치 미분 오차가 커진다. 선형계 조건수도 같이 나빠진다.

Delaunay 삼각화는 가능한 것 가운데 최소 내각을 최대화한다. 가장 뾰족한 각을 가능한 한 둔하게 만든다. 점 위치만 주면 답이 하나라, 사후에 같은 점을 다시 이어도 연결이 흔들리지 않는다.

전진 프론트는 순서가 반대다. 점을 먼저 뿌리지 않는다. 경계선을 출발점으로 두고 안쪽으로 삼각형을 하나씩 붙인다. 새 점은 그때 만든다. 매 순간 제일 나은 요소를 집는 지역 탐욕이다.

빈 외접원 — Delaunay를 정의하는 성질#

규정 문장은 하나다. 어떤 삼각형의 외접원도 내부에 다른 정점을 품지 않는다. 빈 외접원 조건(empty circumcircle)이다. 점 $d$ 가 삼각형 $abc$ 의 외접원 내부에 있는지는 행렬식 하나다. $abc$ 를 반시계 방향으로 두고 다음을 계산한다.

\begin{vmatrix} a_x-d_x & a_y-d_y & (a_x-d_x)^2+(a_y-d_y)^2 \\ b_x-d_x & b_y-d_y & (b_x-d_x)^2+(b_y-d_y)^2 \\ c_x-d_x & c_y-d_y & (c_x-d_x)^2+(c_y-d_y)^2 \end{vmatrix} > 0

여기서 $a,b,c$ 는 삼각형 세 정점, $d$ 는 검사할 점이다. 행렬식이 양수면 $d$ 는 외접원 안, 음수면 밖, 0이면 원 위에 정확히 놓인다. 나눗셈도 제곱근도 없다. 부호만 보면 되니 부동소수점 오차에도 비교적 안전하다. 이 부호가 삽입 한 스텝에서 어떻게 깨지는지는 빈 외접원과 Bowyer–Watson 점 삽입 쪽에 적어 두었다.

Bowyer–Watson: 공동을 파고 다시 채운다#

로슨 플립을 반복해도 Delaunay 격자는 만들어진다. 더 곧은 길은 1981년, 에이드리언 바우어와 데이비드 왓슨이 같은 학술지에 나란히 낸 증분 삽입이다.

이미 Delaunay인 격자에 새 점 $p$ 를 넣으면, $p$ 를 외접원 안에 품는 삼각형은 전부 무효가 된다. 지우면 다각형 구멍, 즉 **공동(cavity)**이 생긴다. 공동의 경계 에지 각각과 $p$ 를 이어 새 삼각형을 채우면 끝이다.

절차는 세 단계다.

새 점 $p$ 를 외접원에 포함하는 삼각형(나쁜 삼각형)을 모두 찾는다.
나쁜 삼각형들이 이루는 공동의 경계 에지를 뽑는다. 경계 에지는 나쁜 삼각형 하나에만 속한 에지다.
나쁜 삼각형을 지우고, 경계 에지마다 $p$ 를 이어 삼각형을 새로 만든다.

시작할 때는 모든 입력 점을 감싸는 거대한 슈퍼 삼각형 하나에서 출발한다. 점을 다 넣은 뒤 슈퍼 삼각형 정점에 걸린 삼각형을 걷어내면 최종 Delaunay 격자만 남는다.

캔버스를 클릭해 점을 넣으면, 매 삽입마다 나쁜 삼각형이 지워지고 공동이 다시 채워진다.

Show circumcircles

Click on the canvas to insert a point. Each insertion rebuilds the Delaunay triangulation via Bowyer–Watson. Toggle the circumcircles: every one stays empty of other vertices — that is the Delaunay property.

외접원 표시를 켜 보자. 어느 원도 다른 정점을 담지 않아야 한다. 점을 촘촘히 찍어도 이 성질은 깨지지 않는다. 그게 Delaunay다.

Python — 증분 Delaunay를 처음부터 짠다#

외접원 판정부터 삽입 커널까지, 라이브러리 없이 순수 Python으로 옮긴 코드다.

import numpy as np
 
def circumcircle(a, b, c):
    """세 점의 외접원 중심과 반지름 제곱을 반환. 일직선이면 None."""
    ax, ay = a; bx, by = b; cx, cy = c
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if abs(d) < 1e-12:
        return None
    a2, b2, c2 = ax * ax + ay * ay, bx * bx + by * by, cx * cx + cy * cy
    ux = (a2 * (by - cy) + b2 * (cy - ay) + c2 * (ay - by)) / d
    uy = (a2 * (cx - bx) + b2 * (ax - cx) + c2 * (bx - ax)) / d
    return (ux, uy), (ax - ux) ** 2 + (ay - uy) ** 2
 
def in_circumcircle(p, tri, pts):
    """점 p가 삼각형 tri의 외접원 내부면 True."""
    cc = circumcircle(pts[tri[0]], pts[tri[1]], pts[tri[2]])
    if cc is None:
        return False
    (ux, uy), r2 = cc
    return (p[0] - ux) ** 2 + (p[1] - uy) ** 2 < r2 - 1e-12
 
def bowyer_watson(points):
    """증분 삽입으로 Delaunay 삼각화를 만든다."""
    pts = list(points)
    xs = [p[0] for p in pts]; ys = [p[1] for p in pts]
    dmax = 20 * max(max(xs) - min(xs), max(ys) - min(ys))
    cx, cy = (min(xs) + max(xs)) / 2, (min(ys) + max(ys)) / 2
    base = len(pts)                                    # 슈퍼 삼각형 정점 시작 인덱스
    pts += [(cx - dmax, cy - dmax), (cx + dmax, cy - dmax), (cx, cy + dmax)]
    tris = [(base, base + 1, base + 2)]
 
    for pi in range(base):                             # 점을 하나씩 삽입
        p = pts[pi]
        bad = [t for t in tris if in_circumcircle(p, t, pts)]
        edge_count = {}                                # 공동 경계 = 한 번만 등장한 에지
        for t in bad:
            for e in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]:
                key = tuple(sorted(e))
                edge_count[key] = edge_count.get(key, 0) + 1
        boundary = [e for e, n in edge_count.items() if n == 1]
        tris = [t for t in tris if t not in bad]
        tris += [(a, b, pi) for a, b in boundary]      # p를 경계 에지에 잇는다
 
    return [t for t in tris if all(i < base for i in t)]  # 슈퍼 삼각형 제거
 
if __name__ == "__main__":
    rng = np.random.default_rng(3)
    P = [tuple(xy) for xy in rng.random((12, 2))]
    T = bowyer_watson(P)
    print(f"점 {len(P)}개 -> 삼각형 {len(T)}개")
    print("첫 삼각형 정점 인덱스:", T[0])

실행하면 점 12개 -> 삼각형 17개가 나온다. 볼록 껍질 안의 삼각형 개수는 대략 $2n - 2 - h$ 개다( $n$ 은 점 수, $h$ 는 볼록 껍질 위 점 수). 순수 삽입 커널이 60줄 안에 들어온다.

열두 점을 그 식에 넣으면 $2\times 12-2=22$ 다. 출력이 17이니, 식이 딱 맞을 때 껍질 위의 점은 $22-17=5$ 다. 슈퍼 삼각형 스케일은 점 구름 폭의 20배다. 샘플은 단위 정사각형에서 뽑히므로 그 폭은 1을 넘지 않고, 세 꼭짓점은 상자 중심에서 그 스케일만큼 떨어져 입력 점 12개를 안에 둔다. 20은 포함을 위한 여유지, 삼각형이 뾰족한지와는 무관하다.

코드가 두는 슈퍼 삼각형의 밑변은 스케일의 두 배이고, 높이도 스케일의 두 배다. 넓이는 스케일 제곱의 두 배가 된다. 점 구름 폭을 1로 두면 스케일은 20, 넓이는 800이다. 점이 들어 있는 정사각형 넓이 1과 견주면 슈퍼 삼각형은 그 800배다. 마지막에 슈퍼 정점에 걸린 삼각형을 걷어내므로 이 800은 결과에 남지 않는다. 남는 것은 점 12개로 만든 삼각형 17개뿐이다.

본문 행렬식은 부호만 본다. 코드는 반지름 제곱에서 1e-12를 뺀 값보다 거리 제곱이 작을 때만 내부로 친다. 원 위에 거의 앉은 점은 밖으로 밀린다. 세 점이 거의 일직선이면 같은 1e-12로 중심 계산을 포기한다. 공동 쪽은 에지가 몇 번 나왔는지만 센다. 두 번 나온 에지는 구멍 안쪽이다. 한 번 나온 에지만 새 삼각형의 밑변이 되고, 그 개수만큼 삼각형이 추가된다.

경계는 어떻게 지키나 — 제약과 파이프#

여기까지는 점만 있으면 되는 문제였다. 실제 CFD 형상은 다르다. 익형의 표면, 실린더의 벽처럼 반드시 격자 에지로 존재해야 하는 경계선이 있다. 그런데 순수 Delaunay는 이 경계 에지를 마음대로 가로지를 수 있다.

해법이 제약 Delaunay 삼각화(constrained DT)다. 잃어버린 경계선 $AB$ 를 복구하려면, 먼저 $AB$ 가 가로지르는 삼각형들의 띠, 이른바 **파이프(pipe)**를 찾는다. 점 $A$ 에 붙은 삼각형에서 출발해, $AB$ 가 뚫고 지나가는 이웃 삼각형을 따라가다 $B$ 에 닿으면 파이프가 완성된다.

파이프가 정해지면 그 안을 다시 삼각화해 $AB$ 를 에지로 복원한다. 대각선 교환(swapping of diagonals)이나 분할 정복으로 처리한다. 경계가 살아난 대신, 그 근처 삼각형은 더 이상 완벽한 Delaunay가 아닐 수 있다. 형상 충실도와 격자 품질을 맞바꾸는 셈이다.

전진 프론트 — 탐욕으로 요소를 쌓는다#

전진 프론트 기법(AFT)은 이 경계 문제를 아예 뒤집어 푼다. 경계가 곧 출발선이다.

경계 세그먼트 목록에서 기준 세그먼트 $AB$ 를 하나 고른다. 이 $AB$ 와 가장 좋은 삼각형을 이루는 점 $C$ 를 찾는다. 후보는 두 종류다. 기존 프론트 위의 점을 쓰거나, 새 내부 점 $I$ 를 만들거나. 둘 중 품질이 높은 쪽을 택해 삼각형을 만들고, 프론트를 갱신한다. 새로 생긴 에지는 프론트에 추가하고, 이미 프론트에 있던 에지는 닫힌 것으로 보고 제거한다.

품질은 두 요소의 곱으로 잰다.

\lambda = \alpha \cdot \delta

여기서 $\alpha$ 는 삼각형의 형상 계수(정삼각형이면 1, 납작할수록 0), $\delta$ 는 요구된 요소 크기와 실제 에지 길이의 일치도(둘이 같으면 1)다. $\lambda$ 가 클수록 좋은 요소다.

형상 계수는 이렇게 짤 수 있다.

import math
 
def shape_quality(a, b, c):
    """정삼각형이면 1, 납작할수록 0에 가까운 형상 계수 alpha."""
    l1 = math.dist(a, b); l2 = math.dist(b, c); l3 = math.dist(c, a)
    area = abs((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) / 2
    denom = l1 * l1 + l2 * l2 + l3 * l3
    return 4 * math.sqrt(3) * area / denom if denom > 0 else 0.0

프론트가 비면(닫힌 경계가 안쪽으로 다 모이면) 격자가 완성된다. 매 단계에서 국소 최적을 고르므로 AFT의 초기 격자 품질은 대개 Delaunay보다 높다. 특히 경계층처럼 방향성이 강한 이방성 격자를 만들 때 강점이 크다.

한 변의 길이가 같은 정삼각형에서 넓이는 변 제곱에 $\sqrt{3}/4$ 를 곱한 값이고, 세 변 제곱의 합은 변 제곱의 세 배다. 코드의 $4\sqrt{3}$ 곱하기 넓이를 그 합으로 나누면 $4\sqrt{3}\cdot(\sqrt{3}/4)/(3)=1$ 이다. 한 변이 길어지고 높이가 줄면 분자는 줄고 분모는 커져 $\alpha$ 가 0 쪽으로 간다. 크기 일치 $\delta$ 가 1이어도 곱 $\lambda$ 는 같이 작아진다. 납작한 요소가 프론트에 남지 않는 이유다.

대가는 속도다. 새 요소마다 프론트의 다른 에지와 교차하는지 검사해야 한다. 소박하게 짜면 요소 하나에 $O(n)$ , 전체가 $O(n^2)$ 로 부풀 수 있다. 배경 격자나 사분트리(quadtree)로 이 검사를 국소화하지 않으면, 프론트가 길어질수록 교차 검사가 생성을 잡아먹는다.

언제 무엇을 고르나#

두 방법은 경쟁이 아니라 분업에 가깝다.

Delaunay는 점 분포가 주어졌을 때 유일하고 빠른 답을 준다. 볼록 껍질 채우기, 산란 데이터 보간, 사후 격자 개선에 강하다. 반면 경계를 스스로 지키지 못해 제약 처리가 따로 필요하다.

AFT는 경계 충실도와 초기 요소 품질이 뛰어나고 이방성 확장이 쉽다. 대신 구현이 까다롭고 자료구조 설계에 성능이 크게 좌우된다.

그래서 현대 격자 생성기는 둘을 섞는다. AFT로 요소 배치 순서를 정하되 Delaunay 연결 기준으로 점을 잇는 프론트–Delaunay 방식이 대표적이다.

익형 에지가 잘리면 삽입을 거기서 끝내지 않는다#

점 12개가 삼각형 17개로 끝나는 것은, 점이 이미 영역 안에 있을 때의 이야기다. 익형 표면이 대각선에 잘리면 그 17개 중 일부는 형상을 대표하지 못한다. 파이프로 그 띠를 다시 삼각화하거나, 잘린 에지를 프론트의 기준 세그먼트로 남겨 다음 점 $C$ 를 고르게 해야 한다.

$\lambda = \alpha\delta$ 가 큰 삼각형만 쌓는다고 잘린 경계가 돌아오지는 않는다. 제약 에지를 되돌릴 때 점을 어디에 더 찍는지는 Steiner 점으로 잘린 경계를 복구하는 절차에 적혀 있다. 삽입 커널 60줄을 먼저 맞춘 다음, 경계가 잘리는 형상에서만 그 절차로 넘어가면 된다.

점 구름을 삼각형으로 덮는 두 철학 — Bowyer–Watson과 전진 프론트