Incremental Delaunay Triangulation

#define REMOVE_FLAG    (-1)
struct TRIANGLE {
   int a, b, c;
   TRIANGLE() {}
   TRIANGLE(int va, int va, int vc) : a(va), b(vb), c(vc) {}
};
struct EDGE {
    int a, b;
    EDGE() {}
    EDGE(int va, int vb) : a(va), b(vb) {}
    void Unflag() {a = REMOVE_FLAG; b = REMOVE_FLAG;}
    bool Unflagged() {a == REMOVE_FLAG || b == REMOVE_FLAG;}
    bool IsSame(EDGE e) {
    	return ((a == e.a && b == e.b) || (a == e.b && b == e.a)) ;    
    }
};
// 수퍼 삼각형은 입력 점집합을 모두 포함하도록 충분히 큰 삼각형을 선택한다;
// 삼각화한 후 최외각이 convex가 아니면 수퍼 삼각형을 더 크게 다시 잡아야 함.
// 마지막 수정: 2021.4.7, 코드 간결화; 2024.4.5, 단순화;
// max num of triangles = 2 * (N - 2) - 1;
// max num of edegs     = 3 * (N - 2);
int inc_delaunay_triangulation(std::vector<CPoint>& pts, std::vector<TRIANGLE>& V) {
    const int N = pts.size()
    if (N < 3) return 0;
    int tmax = 2 * ((N + 3) - 2);
    int emax = 3 * ((N + 3) - 2);
    std::vector<EDGE> E(emax);
    V.resize(tmax);
    std::vector<CPoint> P(pts.size());
    for (int i = points.size(); i-->0; ) P[i] = pts[i];
    set_supertriangle(P); // P.size()은 +3 만큼 증가함;
    // First element of triangle index is supertriangle vertex;
    int nt = 0;
    V[nt++] = TRIANGLE(N, N + 1, N + 2);
    for (int i = 0; i < N; i++) {
        int ne = 0;
        /* 추가점(P[i])를 외접원에 포함시키는 삼각형을 모두 찾아 제거한다.
        ** 이들 삼각형의 공유되는 내부 edge를 제외한 나머지 외부 edge와
        ** 추가점으로 새로운 삼각형을 만들어 리스트에 추가한다. */ 
        for (int j = 0; j < nt; j++) {
            TRIANGLE &t = V[j];
            if (inCircle(P[t.a], P[t.b], P[t.c], P[i])) {
                // 제거되는 삼각형의 edge는 백업한다; 내부 edge는 중복된다.
                E[ne++] = EDGE(t.a, t.b);
                E[ne++] = EDGE(t.b, t.c);
                E[ne++] = EDGE(t.c, t.a);
                // P[i]를 포함하는 외접원 삼각형은 제거;
                // swap(j, nt-1) to remove triangle at j
                V[j--] = V[--nt];  
            }
        }
        // 중복되어 있는 내부 edge를 제거;
        remove_double_edges(E, ne);
        
        // 외부 edge와 P[i]가 만드는 새 삼각형을 리스트에 추가;
        for (int j = 0; j < ne; j++) {
            if (E[j].Unflagged()) continue;
            V[nt++] = TRIANGLE(E[j].a, E[j].b, i);
        }
    }
    // done with super_triangle;
    V.resize(nt);
    return remove_supertriangle(N, V) ;
}
// double edge 제거;
void remove_double_edges(std::vector<EDGE>& E, int ne) {
    for (int j = 0; j < ne - 1; j++)
    	for (int k = j + 1; k < ne; k++)
            if (E[j].IsSame(E[k])) {
            	E[j].Unflag(); E[k].Unflag();
            }
}
// 삼각형 중에서 수퍼 삼각형의 꼭지점(N,N+1,N+2)을 가지고 있는 것을 제거;
int remove_supertriangle(int N, std::vector<TRIANGLE>& V) {
    int nt = V.size();
    for (int i = 0; i < nt; i++)
        if (V[i].a >= N || V[i].b >= N || V[i].c >= N)
            // swap(i, nt - 1) to remove triangle at i;
            V[i--] = V[--nt];
    V.resize(nt);
    return nt;
}
// A,B,C가 만드는 외접원에 D가 들어가는가?;
int inCircle(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인  경우에..    
}