Geometry & Recognition

마우스 클릭으로 생긴 이미지 상의 입력 점을 가지고 triangulation을 한 후에 각각의 삼각형 내의 영역을 그 영역의 평균 컬러 값으로 채우는 효과이다. voronoi 다이어그램을 이용해서 tiling 하는 효과와 유사하다.

그림: Lena사진(일부)을 이용한 예:

사용 algorithm
  - incremental delaunay triangulation;
  - polygon fill;

Voronoi Tiling:

주어진 점집합에서 세 점을 뽑아서 만든 삼각형의 외접원에 다른 점이 하나도 포함하지 않으면 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

// input  x[0], x[1],....,x[n-1],
// input  y[0], y[1],....,y[n-1];
// calculate z[0]=x[0]^2+y[0]^2, z[1]=x[1]^2+y[1]^2,...;
int dt4(int n, double x[], double y[]) {
    double *z = new double [n] ;
    for(int i = 0; i < n; i++) 
        z[i] = x[i] * x[i] + y[i] * y[i] ;

    int ntri = 0;
    /* 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) {
                        ntri++;
                        // (i, j, k)의 외접원이 다른 점을 포함하지 않으므로 이 삼각형은 
                        // 삼각분할의 한 면을 형성하게 된다.
                        // addTriangle(tri(i, j, k));
                    }
                }
    delete[] z ;
    return ntri;
}

monotone 폴리곤은 한 축 방향으로 꼭짓점의 좌표값이 단조 증가하는 형태를 갖는 단순 폴리곤이다. 단순 폴리곤을 monotone 폴리곤으로 분할하는 방법은 잘 알려져 있고 O(n log n)의 복잡도를 갖는다. monotone 폴리곤은 한 방향에 대해서 이미 정렬이 되어 있기 때문이 그 방향으로 sweep-line 기법을 쓰면 선형 시간 내에 triangulation이 가능하다. 이러한 알고리즘은 이미지 처리에 적용하면 매우 효율적으로 이용이 가능하다. raster 상의 데이터는 scanline 방향으로 항상 정렬이 되어 있기 때문이다.

참고 : http://www.personal.kent.edu/~rmuhamma/Compgeometry/MyCG/PolyPart/polyPartition.htm

Algorithm Monotone_Triangulation
Input: Monotone polygon P (ccw-ness needed!)
Output: Set of triangles
(1) 주어진 폴리곤의 꼭지점을 x-좌표를 기준으로 정렬하여 V에 저장 (O(n) due to monotonicity).
(2) sweep-line 리스트(L)을 초기화 :: V의 처음 두 포인트(이벤트)를 넣음.
(3) while(V is not empty) do
       get next vertex (p) from V ;
       if (p is opposite to points in L) then
            while(size(L)>1) do
                out triangle { first(L), second(L), p }
                remove(first(L))
            end while
            add p to L;
       else 
            while ( size(L)>1 && ccw(last(L), previous(last(L)), p)) do
                 out triangle { last(L), previous(last(L)), p }
                 remove(last(L))
            end while
            add p to L
       end if
     end while/*입력이 ccw이면 ccw-test도 ccw로, cw면 test도 cw로 해야한다.*/ 

/* 구현 코드(C++)의 일부이다.*/
struct endpoint {
    int     vid;                /* polygon vertex id */
    double* x;                  /* X coordinate of the endpoint */
    double* y;                  /* Y coordinate of the endpoint */
} ;
//x-크기순서대로--> y-크기순 서대로; (ascending-order);
static 
int point_compare(double* x1, double* y1, double* x2, double* y2) {
    if (*x1 > *x2) return  1;
    if (*x1 < *x2) return -1;
    if (*y1 > *y2) return  1;
    if (*y1 < *y2) return -1;
    return 0;
}
//qsort에서 사용;
static 
int endpoint_compare(const void* p1, const void* p2) {
    endpoint* e1 = *(endpoint**) p1;
    endpoint* e2 = *(endpoint**) p2;
    return point_compare(e1->x, e1->y, e2->x, e2->y);
}
//
class endpointqueue {
    int n;                      /* total number of endpoints */
    int next;                   /* index of next endpoint on queue */
    endpoint* data;             /* array of all endpoints */
    endpoint** sorted;          /* sorted list of endpoint pointers */
    bool create(int npoly, double polyx[], double polyy[]);
public:
    ~endpointqueue() {
        if(data) delete[] data ;
        if(sorted) delete[] sorted;
    }
    endpointqueue() : data(NULL), sorted(NULL), n(0), next(0) {}
    endpointqueue(int npoly, double polyx[], double polyy[]) {
        create(npoly, polyx, polyy);
    } ;
    endpoint * getnext() { return (next >= n) ? NULL : sorted[next++];}
} ;
bool endpointqueue::create(int n, double polyx[], double polyy[]) {
    data = new endpoint [n] ;
    this->n = n ;
    this->next = 0;
    for (int i = 0; i < n; i++) {
        data[i].x = &polyx[i] ;
        data[i].y = &polyy[i] ;
        data[i].vid = i;
    };
    sorted = new endpoint* [n] ;
    for (i = 0; i < n; i++) {
        sorted[i] = &data[i];
    };
    //non-optimal sorting; use monotonicity of input;
    qsort(sorted, n, sizeof(endpoint*), endpoint_compare);
    return true ;
}
//deque의 행동을 재정의할 필요가 있음.

struct sweepline {
    int     n;                      /* number of vertices in polygon */
    endpoint* deque;
    int bot, top;
    sweepline() : n(0), deque(NULL) {}
    sweepline(int n_) {
        n = n_; bot = -1; top = -1; 
        deque = new endpoint [n] ;
    }
    ~sweepline(){ if(deque) delete [] deque;}
    endpoint* add(endpoint* e) ;
    void removefirst() { if (bot < top) bot++;} //safe::size()>1
    void removelast() { if (top > bot) top--;}  //safe::size()>1;
    int size() { if (top >= 0 && bot >= 0 && top >= bot) return top - bot + 1; else return 0;}
    endpoint* getlast()     { ASSERT(size() > 0); return &deque[top];}
    endpoint* getprelast()  { ASSERT(size() > 1); return &deque[top - 1];}
    endpoint* getfirst()    { ASSERT(size() > 0); return &deque[bot  ];}
    endpoint* getsecond()   { ASSERT(size() > 1); return &deque[bot + 1];}
    bool opposite(endpoint* e) ;
    void draw(CDC *pDC, DWORD color = RGB(0xFF, 0, 0)) ;
};
endpoint* sweepline::add(endpoint* e) {
    if (bot == -1 && top == -1) //first adding;
        bot = 0;
    top++;
    deque[top].x = e->x ;
    deque[top].y = e->y ;
    deque[top].vid = e->vid ;
    return &deque[top];
}
bool sweepline::opposite(endpoint* e) {
    if (e->vid == ((deque[top].vid + 1) % this->n))
        return false;
    else
        return true ;
}
void monotone_triangulate(CPoint P[], int n, 
                          void (*print_triangle)(int a, int b, int c)) {
    double *polyx = new double [n] ;
    double *polyy = new double [n] ;
    for (int i = 0; i < n; i++) {
        polyx[i] = P[i].x;
        polyy[i] = P[i].y;
    };
    //
    endpointqueue eq(n, polyx, polyy);
    sweepline sl(n);
    endpoint *e ;
    sl.add(eq.getnext());
    sl.add(eq.getnext());
    while ((e = eq.getnext())) {
        TRACE("event=%d(x=%5.0f, y=%5.0f)\n", e->vid, *e->x, *e->y);
        if (sl.opposite(e)) {
            while (sl.size() > 1) {
                //(first, second, p);
                print_triangle(sl.getfirst()->vid, sl.getsecond()->vid, e->vid);
                sl.removefirst();
            }
            sl.add(e);
        } else {
            while (sl.size() > 1 && isccw(sl.getlast(), sl.getprelast(), e)) {
                //triangle(last, previous(last), p);
                print_triangle(sl.getlast()->vid, sl.getprelast()->vid, e->vid);
                sl.removelast();
            }
            sl.add(e);
        }
    }
    //
    delete [] polyx;
    delete [] polyy;
};

/*
**  http://blog.naver.com/helloktk/80051120809 에서 옮김.
*/