Geometry & Recognition

#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인  경우에..    
}

주어진 점집합에서 세 점을 뽑아서 만든 삼각형의 외접원에 다른 점이 하나도 포함하지 않으면 triangulation의 기본 삼각형 cell이 된다. 주어진 점으로 만들 수 있는 삼각형의 총개수가 ${}_nC_3$ 이므로, 기본 삼각형을 찾기 위해서는 이들 각각의 삼각형에 대서 나머지 점을 가지고 incircle 테스트를 수행해야 한다. 따라서 이 알고리즘은 ${\cal O}(n^4)$의 스텝이 필요하게 된다.

/*brute force attack*/
foreach point p1
    foreach point p2 
        foreach point p3
            foreach point p4 
                if(incircle(p1,p2,p3,p4))
                    iscell=false;
                    break ;
            endfor;
            if(iscell) 
                 add(triangle(p1,p2,p3));
        endfor;
    endfor;
endfor;

세 점에 의해서 형성이 되는 외접원은 대수적으로 쉽게 구할 수 있다. 여기서는 좀 더 기하학적인 접근을 쓰면, 평면의 점은 

$$(x, y)\rightarrow (x, y, z=x^2 + y^2)$$

의 mapping에 의해서 3차원 paraboloid 곡면의 점으로 옮길 수 있다. paraboloid 위의 세 점이 형성하는 3차원에서 평면이 paraboloid를 절단하는 곡선을 $x-y$ 평면으로 정사영하면 원이 된다는 것을 쉽게 알 수 있다.(incircle 포스팅 참조). 따라서 주어진 점이 세 점의 외접원에 포함되는가를 테스트하는 작업은 이 점을 paraboloid로 올렸을 때의 점과 (paraboloid로 올려진) 외접원의 3점이 형성하는 3차에서의 평면과 관계를 조사하는 것으로 바꿀 수 있다.

주어진 점이 외접원에 포함되면 paraboloid로 변환된 점은 평면의 아래에 놓이고, 외접원 밖의 점이면 평면 위에 놓이게 된다. 물론 외접원 위의 점은 평면에 놓인다. 따라서 평면의 법선 벡터를 구하고, 삼각형의 한 꼭짓점을 기준한 주어진 점의 변위 벡터와 내적을 구하면 내적의 값은 평면 위인지, 아래인지 또는 평면에 놓인 점인가에 따라서 부호가 달라진다. 평면의 수직 벡터를 고정하면(예제는 아래 방향: $n_z < 0$), 평면 위에 놓인 점과의 내적은 음수, 평면 아래에 놓인 점과의 내적은 양수가 되고, 평면의 점과의 내적은 0이다. 

주어진 세 점이 만드는 외접원 내부(and 경계)에 들어가는 점이 없으면 이 삼각형을 선택한다.

** 참고 : Computational Geometry in C(2nd Edition) by Joseph O'Rourke

std::vector<Triple> dt4(const std::vector<double>& x, const std::vector<double>& y) {
    const int n = x.size();
    if (n < 3) return std::vector<Triple> (); // null_vec;
    std::vector<double> z(n);
    for (int i = 0; i < n; i++) 
        z[i] = x[i] * x[i] + y[i] * y[i] ;

    std::vector<Triple> triples;
    /* For each triple (i,j,k) */
    for (int i = 0; i < n - 2; i++ )
        for (int j = i + 1; j < n; j++ )
            for (int k = i + 1; k < n; k++ )
                if ( j != k ) {
                    /* Compute normal to triangle (i,j,k)::  outter_product(j-i, k-i)*/
                    double nx = (y[j] - y[i]) * (z[k] - z[i]) - (y[k] - y[i]) * (z[j] - z[i]); 
                    double ny = (x[k] - x[i]) * (z[j] - z[i]) - (x[j] - x[i]) * (z[k] - z[i]);
                    double nz = (x[j] - x[i]) * (y[k] - y[i]) - (x[k] - x[i]) * (y[j] - y[i]);
                    
                    /* Only examine faces on bottom of paraboloid: nz < 0. */
                    int flag = (nz < 0);
                    if (flag) {
                        /* For each other point m */
                        for (int m = 0; m < n; m++) {
                            /* Check if m above (i,j,k)::i점을 기준으로 m 과 
                            ** normal 간의 내적으로 체크(내적이 양수이면 
                            ** m이 원의 내부에 완전히 들어 있는 경우가 된다. 
                            ** 0인 경우는 원주상에 놓인 경우이므로 배제한다
                            */
                            flag &= ((x[m]-x[i])*nx + (y[m]-y[i])*ny + (z[m]-z[i])*nz <= 0);
                        }
                    }
                    if (flag) {
                        // (i, j, k)의 외접원이 다른 점을 포함하지 않으므로 이 삼각형은 
                        // 삼각분할의 한 면을 형성하게 된다.
                        triples.push_back(Triple(i, j, k));
                    }
                }
    return triples;
}

이미지에서 Voronoi diagram으로 영역을 분할할 때 각 픽셀이 어느 Voronoi cell에 포함되는가를 알아야 하는 경우가 있다. 보통은 Voronoi 다이어그램으로 구한 cell을 폴리곤으로 표현하고, 해당 픽셀이 어느 폴리곤에 들어가는 가는 체크 해야 한다. 그러나, 이 과정은 복잡하고 계산이 많이 발생한다. 이미지에 만들어진 Voronoi diagram의 경우 cell mask를 이용하면 해당 픽셀이 어느 cell에 들어있는지를 바로 판단할 수 있다. 특히, cell의 개수가 적은 경우 mask를 gray 이미지로 처리할 수 있어서 메모리 사용도 줄일 수 있다.

Voronoi diagram의 이미지화 과정은 Voronoi 알고리즘을 이용할 필요는 없고 단지 각 cell을 형성하는 픽셀들은 그 cell의 중심까지 거리가 다른 cell보다 가깝다는 사실만 이용한다.

void rasterize_voronoi(std::vector<CPoint>& vorocenter, 
                       BYTE *image, int width, int height) {
    std::vector<BYTE> red(vorocenter.size()), 
                      green(vorocenter.size()), 
                      blue(vorocenter.size());
    for (int i = vorocenter.size(); i-->0;) {
    	red[i]   = rand() % 256;
        green[i] = rand() % 256;
        blue[i]  = rand() % 256;
    }
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            int min_id = 0; 
            int dist2_min = INT_MAX;
            for (int k = vorocenter.size(); k-->0;) {
                int dx = x - vorocenter[k].x;
                int dy = y - vorocenter[k].y;
                int dist2 = dx * dx + dy * dy;
                if (dist2 < dist2_min) {
                    dist2_min = dist2;
                    min_id = k;
                }
            }
            *image++ = blue[min_id];
            *image++ = green[min_id];
            *image++ = red[min_id];
        }
    }
    // draw cell center;
}