Incremental Delaunay Triangulation

Incremental Delaunay Triangulation은 점을 하나씩 삽입하면서 영향을 받는 local region만 수정하여 Delaunay 조건을 유지하는 삼각분할 알고리즘이다. 전체 삼각분할을 다시 계산하지 않고도 새로운 점을 효율적으로 반영할 수 있기 때문에 동적 환경에서의 Delaunay Triangulation에 가장 널리 사용되는 방법 중 하나이다.

알고리즘:

• 모든 점을 포함하는 매우 큰 Super Triangle을 생성한다.
• 점들을 하나씩 삽입한다.
• 새 점이 영향을 주는 삼각형들을 찾는다.
• 해당 삼각형을 제거하여 빈 영역(cavity)을 만든다.
• 새 점과 cavity 경계를 연결하여 새로운 삼각형을 만든다.
• Delaunay 조건이 깨진 Edge는 Flip하여 복원한다.
• 모든 점을 삽입한 뒤 Super Triangle과 연결된 삼각형을 제거한다.

계산 복잡도: 점 위치 탐색 자료구조를 사용하는 경우

• 평균 시간복잡도: \( O(n\log n)\)
• 최악의 경우: \(O(n^2)\)

실제 무작위 점 분포에서는 대부분 \( O(n\log n)\)의 성능을 보인다. 

장점:

• 구현이 비교적 단순하다.
• 점을 동적으로 추가하기 쉽다.
• 평균적으로 매우 빠르다.
• 메모리 사용량이 적다.
• 병렬화 및 온라인 처리에 적합하다.

단점

• 점의 삽입 순서에 따라 수행 시간이 달라질 수 있다.
• 최악의 경우 \(O(n^2)\)까지 증가할 수 있다.
• 점이 거의 일직선으로 배열되거나 특수한 배치를 가지면 성능이 저하될 수 있다.
• 효율적인 점 위치 탐색(point location) 자료구조가 필요하다.

#define REMOVE_FLAG    (-1)
struct TRIANGLE {
   int a, b, c;
   TRIANGLE() {}
   TRIANGLE(int va, int vb, int vc) : a(va), b(vb), c(vc) {}
};
struct EDGE {
    int a, b;
    EDGE() {}
    EDGE(int va, int vb) : a(va), b(vb) {}
    void MarkRemoved() { a = REMOVE_FLAG; b = REMOVE_FLAG; }
    bool IsRemoved() const { return (a == REMOVE_FLAG || b == REMOVE_FLAG); }
    bool operator == (const EDGE& e) const {
    	return (a == e.a && b == e.b) || (a == e.b && b == e.a) ;    
    }
};
// 수퍼 삼각형은 입력 점집합을 모두 포함하도록 충분히 큰 삼각형을 선택한다;
// 삼각화한 후 최외각이 convex가 아니면 수퍼 삼각형을 더 크게 다시 잡아야 함.
// 마지막 수정: 2021.4.7, 코드 간결화; 2024.4.5, 단순화; 2026.8. 추가 단순화;
std::vector<TRIANGLE> inc_delaunay_triangulation(std::vector<CPoint>& pts) {
    const int N = pts.size();
    if (N < 3) return std::vector<TRIANGLE> ();
    std::vector<EDGE> E;
    E.reserve(3 * (N + 1));
    std::vector<TRIANGLE> V;
    V.reserve(2 * (N + 1));
    
    std::vector<CPoint> P = pts;  // 복사;
    set_supertriangle(P); // P.size()은 +3 만큼 증가함;
    // First element of triangle index is supertriangle vertex;
    V.push_back(TRIANGLE(N, N + 1, N + 2));
    for (int i = 0; i < N; i++) {
        E.clear(); // reset edge container;
        /* 추가점(P[i])를 외접원에 포함시키는 삼각형을 모두 찾아 제거한다.
        ** 이들 삼각형의 공유되는 내부 edge를 제외한 나머지 외부 edge와
        ** 추가점으로 새로운 삼각형을 만들어 리스트에 추가한다. */ 
        for (int j = 0; j < V.size(); ) {
            TRIANGLE &t = V[j];
            if (in_circle(P[t.a], P[t.b], P[t.c], P[i])) {
                // 제거되는 삼각형의 edge는 백업한다; 내부 edge는 중복된다.
                E.push_back(EDGE(t.a, t.b));
                E.push_back(EDGE(t.b, t.c));
                E.push_back(EDGE(t.c, t.a));
                
                // P[i]를 포함하는 외접원 삼각형은 제거;
                // swap(j, nt-1) to remove triangle at j
                V[j] = V.back();  
                V.pop_back();
            } else {
                ++j;
            }
        }
        // 중복되어 있는 내부 edge를 제거;
        remove_double_edges(E);
        
        // 외부 edge와 P[i]가 만드는 새 삼각형을 리스트에 추가;
        for (int j = 0; j < E.size(); j++) {
            if (E[j].IsRemoved()) continue;
            V.push_back(TRIANGLE(E[j].a, E[j].b, i));
        }
    }
    // remove triangles with supertriangle vertices;
    // having a vertex number greater than N;
    for (int i = 0; i < V.size(); ) {
        if (V[i].a >= N || V[i].b >= N || V[i].c >= N) {
            V[i] = V.back(); // i번째 원소를 마지막 원소와 swap해서 제거;
            V.pop_back();
        } else {
            i++;
        }
    }
    return V;
}
// double edge 제거 (삼중 이상으로 제거되는 경우를 방지);
void remove_double_edges(std::vector<EDGE>& E) {
    const int n = E.size();
    for (int j = 0; j < n - 1; j++) {
    	if (E[j].IsRemoved()) continue;  // 이미 제거된 에지
        for (int k = j + 1; k < n; k++) {
            if (E[k].IsRemoved()) continue;// 이미 제거된 에지
            if (E[j] == E[k]) {
            	E[j].MarkRemoved(); 
                E[k].MarkRemoved();
                break; //짝을 찾음: k-loop 벗어남;
            }
        }
    }
}
// A,B,C가 만드는 외접원에 D가 들어가는가?;
int in_circle(CPoint A, CPoint B, CPoint C, CPoint D) {
    CPoint AD = A - D, BD = B - D, CD = C - D;
    double ccw = CCW(A, B, C);
    double det = norm2(AD) * cross(BD, CD)
               + norm2(BD) * cross(CD, AD)
               + norm2(CD) * cross(AD, BD);
    if (ccw > 0) return det > 0; //CCW인 경우에..
    else         return det < 0; //CW인  경우에..    
}