728x90

평면 상에 주어진 점들의 convex hull을 구하는 알고리즘 중의 하나인 Graham scan에서는 먼저 주어진 점들을 한 점을 기준으로 각도로 정렬하는 과정이 필요했다. 그러면 점들이 순차적으로 연결된 단순 다각형에서는 sorting과정이 없이 Graham scan 알고리즘을 바로 적용이 가능한 것일까? 이에 대한 counter example은 쉽게  찾을 수 있다. 단순 폴리곤에 대해서도 항상 각도에 대한 정렬 과정이 필요한 것인가? 답은 아니다. 정렬 과정이 없이도 단순 폴리곤에 대해서는 쉽게 convex hull을 구할 수 있는 알고리즘이 존재한다. 정렬 과정이 없이 단순 폴리곤의 convex hull을 찾는 알고리즘에 대해서 알아보자.

Melkman Algorithm;

우선 폴리곤의 이웃하는 세 꼭짓점을 잡아서, 반시계 방향의 삼각형을 구성한다. 이 삼각형을 deque에 넣는다 (bottom = top). 폴리곤을 순환하면서 새 점이 들어올 때 이미 만들어진 convex hull의 내부점이 아니면 이 점이 포함되도록 convex hull을 업데이트한다: Graham scan 알고리즘의 scanning 과정을 bottom을 기준으로 반시계 방향으로 convexity를 만족시킬 때까지 bottom을 제거하고, top을 기준으로는 시계방향으로 convexity를 만족시킬 때까지 top을 제거한다. 이 과정이 끝나면 새 점을 deque에 추가한다. 이 과정을 나머지 모든 점들에 대해서도 진행한다.

int LEFTSIDE(CPoint A, CPoint B, CPoint C){ // cross(AB, AC) > 0: C가 AB 대해서 왼쪽
	return ((B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x)) > 0;
}
int convexhull_spolygon(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
    int N = pts.size();
    if (N < 3) return 0;
    std::vector<CPoint> Q(2 * N + 1);  // deque;
    int bot = N - 2, top = bot + 3;
    // |0|1|......|N-2|N-1|N+0|N+1|N+2|..............|2N|;
    //              bot    top    
    Q[bot] = Q[top] = pts[2];
    if (LEFTSIDE(pts[0], pts[1], pts[2])) { //2->0->1->2; Q에 ccw로 입력;
        Q[bot + 1] = pts[0]; Q[bot + 2] = pts[1];
    } else {
        Q[bot + 1] = pts[1]; Q[bot + 2] = pts[0];
    }
    int i = 2; // last touched index;
    while (++i < N) {
        // 기존의 convex_hull에 들어있으면 제외.
        if (LEFTSIDE(Q[bot + 0], Q[bot + 1], pts[i]) && 
            LEFTSIDE(Q[top - 1], Q[top + 0], pts[i]))
            continue;
        // bottom에 대해서 ccw 방향으로 체크(bot 증가 방향)
        // pts[i]가 (bot)->(bot+1)라인의 오른쪽에 있으면, 기존의 bottom을 제외;
        while (!LEFTSIDE(Q[bot], Q[bot + 1], pts[i]) && (bot < top)) bot++;
        Q[--bot] = pts[i];
        // 추가점에서 top->top-1의 ccw 방향으로 convexity 체크하여 만족하지 않은 경우
        // top을 감소
        while (!LEFTSIDE(Q[top - 1], Q[top], pts[i]) && (top > bot)) top-- ;
        Q[++top] = pts[i];
    }
    hull.resize(top - bot);
    for (int i = hull.size() - 1; i >= 0; --i) hull[i] = Q[i + bot];
    return hull.size();
};

**네이버 블로그 이전;

 

 

'Computational Geometry' 카테고리의 다른 글

Minimum Enclosing Circle  (0) 2021.03.01
단순 다각형 만들기  (0) 2021.01.25
단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Posted by helloktk

댓글을 달아 주세요

728x90

단순 다각형의 무게중심(center of gravity or center of mass)은 다각형을 균일한 밀도의 판으로 생각했을 때 판의 무게중심과 같다. 가장 단순한 다각형인 삼각형의 무게중심은 세 꼭짓점의 산술평균으로 표현된다.

$$ \text{CoG} = \frac{1}{3} ({{\bf P} + {\bf Q} + {\bf R}}).$$

다각형은 삼각형으로 분할되므로 이 분할된 삼각형의 무게중심을 이용하면 쉽게 계산할 수 있다. 분할된  삼각형의 무게중심을 면적으로 가중치를 준 평균값이 다각형의 무게중심이 된다. 

 

실제 계산에서는 다각형을 삼각 분할하지 않고도 간단한 방법에 의해서 무게중심을 구할 수 있다. 원점과 다각형의 각 변의 꼭짓점을 이용해서 삼각형들을 구성하면 원래의 다각형을 겹치게(원점이 내부에 있으면 겹침이 없다) 분할할 수 있다. 분할된 삼각형으로 무게중심을 구할 때 겹치는 영역의 기여를 제거해야 한다. 그런데 다각형 밖의 영역을 분할하는 삼각형은 다각형 내부를 분할하는 삼각형과는 다른 orientation을 가지게 된다. 삼각형의 면적은 한 꼭짓점을 공유하는 두 변의 외적에 비례하므로, 반대의 orientation을 갖는 삼각형은 자동으로 반대 부호의 면적을 가지게 된다. 따라서 분할된 삼각형 면적 가중치를 외적으로 주어서 무게중심을 구하면 겹치는 영역이 자동으로 상쇄되는 효과를 얻을 수 있다.

노란색 영역에 있는 분할 삼각형은 음의 면적

$$\begin{align} \text{CoG} &= \frac{1}{\text{다각형 면적}} \sum (\text{삼각형 CoG})  (\text{면적})  \\ &= \frac{1}{\text{다각형 면적}} \sum \frac{1}{3} \left( {\bf V}_i + {\bf V}_{i+1} + {\bf O}\right ) \frac{ {\bf V}_{i} \times {\bf V}_{i+1} }{2} \\ &= \frac{1}{3}\frac{1}{   \text{다각형 면적} }\sum ( {\bf V}_{i} + {\bf V}_{i+1}) \frac{{\bf V}_{i} \times {\bf V}_{i+1}}{2} \end{align}$$ 

다각형의 면적($=\sum \frac{1}{2}({\bf V}_i \times {\bf V}_{i+1})$)을 구할 때 삼각형과 동일하게 orientation에 따라 부호를 포함하도록 설정하면 다각형의 면적 부호가 삼각형의 면적 부호로 상쇄되므로 다각형의 orientation에 무관하게 성립하는 공식이 된다.

POINT polygon_centroid(POINT V[], int N) {
    double cx = 0, cy = 0, area2 = 0;
    for(int i = 0, j = N - 1; i < N; j = i++) {
        double tri_area2 = V[i].x * V[j].y - V[i].y * V[j].x;
        cx += (V[i].x + V[j].x) * tri_area2;
        cy += (V[i].x + V[j].x) * tri_area2;
        area2 += tri_area2;
    }
    cx /= 3 * area2;
    cy /= 3 * area2;
    return POINT(cx, cy)
};

** 네이버 블로그에서 수정 이전;

'Computational Geometry' 카테고리의 다른 글

단순 다각형 만들기  (0) 2021.01.25
단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Point in Polygon  (2) 2020.12.14
Posted by helloktk

댓글을 달아 주세요

728x90

평면상의 다각형(모서리의 교차가 없는 단순 다각형)의 면적을 구하는 것은 단순하지 않을 것처럼 보이지만 계산식은 무척이나 간단하게 주어진다. 기본적인 아이디어는 다각형에 임의의 점을 찍으면 이 점과 이웃한 두 개의 꼭짓점으로 형성이 되는 삼각형의 합으로 다각형을 분할할 수 있다. 분할된 삼각형의 면적을 구하여 합산하면 다각형의 면적을 구할 수 있다.

 

세 점 ${\bf P, Q, R}$(이 순서대로 반시계방향으로 배열)이 만드는 삼각형의 면적은

$$\text {삼각형의 면적}=  \frac{1}{2} ({\bf R} - {\bf  Q})  \times ( {\bf P} - {\bf Q}); \quad  (\Rightarrow \text {시계 방향이면 면적이 음수})$$

로 주어지므로, 꼭짓점이 ${\bf P}_0(x_0, y_0), {\bf P}_1(x_1, y_1),....$(반시계 방향)으로 주어지는 $N$각형의 면적은 아래와 같이 주어진다. 

$$\begin{align} \text{다각형 면적} &= \sum \text{각 삼각형의 면적} \\ &= \frac{1}{2}\sum ({\bf P}_{i+1}-{\bf Q})\times ({\bf P}_{i}-{\bf Q})\quad \quad ({\bf Q}\text{는  임의의 점})\\ &= \frac{1}{2} \sum\left(  {\bf P}_{i+1} \times {\bf P}_{i} - {\bf P}_{i+1}\times {\bf Q} + {\bf P}_{i} \times {\bf Q}\right) \\ &=\frac{1}{2} \sum {\bf P}_{i+1} \times {\bf P}_{i} \\ &= \frac{1}{2}\sum \left( x_{i+1} y_{i} -x_{i}y_{i+1} \right) \end{align}$$

이 결과는 $\bf Q$에 무관하다. 다각형의 꼭짓점이 시계 방향으로 정렬이 된 경우는 면적이 음수로 나온다(윈도우 DC는 위-아래가 역전되어 있으므로 orientation이 반대로 보인다). 그리고 이 공식은 단순 다각형에만 적용이 되고 모서리의 교차가 있는 경우에는 적용이 되지 않는다.

double simple_polygon_area2D(POINT point[], int N) {
    double area = 0;
    for (int i = 0, j = N - 1; i < N; j = i++) 
        area += (point[i].x * point[j].y) - (point[i].y * point[j].x);
    area /= 2;
    // return area;  // signed area;
    return area < 0 ? -area: area;
}

**네이버 블로그 수정 이전;

'Computational Geometry' 카테고리의 다른 글

단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Point in Polygon  (2) 2020.12.14
Incremental Delaunay Triangulation  (1) 2020.12.01
Posted by helloktk

댓글을 달아 주세요

728x90

세 점 $A$, $B$, $C$가 만드는 삼각형의 외접원의 반지름을 구해보자.  세 변을 나타내는 벡터는

$$\vec{a}= \vec{OC}-\vec{OB}, \quad \vec{b}=\vec{OA}-\vec{OC},\quad\vec{c}=\vec{OB}-\vec{OA}.$$

세 점이 일직선 위에 놓이지 않기 위해서는 

$$\vec{a}\times \vec{b} \ne  0$$

를 만족해야 한다. 

그림의 삼각형에서 $\angle(BOC)$가 $\angle(BAC)$의 두 배이므로,

$$\sin \theta =\frac{a/2}{R}.$$

이제 삼각형 면적은 세 변의 길이 ($a, b, c$)와 외접원의 반지름($R$)으로 표현된다:

$$\text{area}=\frac{1}{2}bc \sin \theta ~~\longrightarrow~~ \frac{abc}{4R}.$$

또한, 삼각형의 면적을 세 점이 일직선 위에 있는가를 테스트하는 식으로 표현하면 

$$\text{area}=\frac{1}{2} | \vec{a}\times \vec{b}|$$

처럼 쓸 수 있다. 따라서 삼각형의 외접원의 반지름은 아래처럼 주어진다.

$$R = \frac{abc}{4\cdot \text{area}} = \frac{| \vec{a}||\vec{b} ||\vec{c} |}{2|\vec{a}\times \vec{b}| }.$$

#define SQR(x) ((x) * (x))
double circumradius(CPoint A, CPoint B, CPoint C) {
    double ax = C.x - B.x, ay = C.y - B.y;
    double bx = A.x - C.x, by = A.y - C.y;
    double crossab = ax * by - ay * bx;
    if (crossab != 0.) { 
        double a = sqrt(SQR(ax) + SQR(ay));
        double b = sqrt(SQR(bx) + SQR(by)); 
        double cx = B.x - A.x, cy = B.y - A.y;       
        double c = sqrt(SQR(cx) + SQR(cy));
        return (0.5 * a * b * c / fabs(crossab));
    } else 
        return ( -1.0 ); 
};
Posted by helloktk

댓글을 달아 주세요

Chain Hull

Computational Geometry 2012. 9. 16. 19:28
728x90

참고: http://cm.bell-labs.com/who/clarkson/2dch.c에 구현된 알고리즘을 수정함:

  1. 점집합을 x-좌표의 값(같으면, y좌표의 값) 순으로 정렬한다. 시작점은 hull의 꼭짓점이 된다.
  2. x-좌표 순서대로 현재 점이 hull의 꼭짓점인가를 검사한다: hull의 꼭짓점이 되기 위해서는 직전까지 구성된 hull의 에지 왼편에 있어야 한다. 오른편에 있으면 직전까지 구성된 hull에서 현재 점이 왼편에 놓일 때까지 꼭짓점을 제거한다.
  3. 남는 점들을 1번 과정의 역순으로 정렬한다. 최대 x-좌표, 최소 y-좌표점이 다시 윗 hull의 시작점이 된다.
  4. x-좌표가 작아지는 순서대로 2번 과정을 수행한다.

/* A->B->C가 반시계방향으로 정렬되거나(< 0), 일직선상에 있으면(=) 참이다*/
static int ccw(POINT A, POINT B, POINT C);

더보기
/* A->B->C가 반시계방향으로 정렬되거나(< 0), 일직선상에 있으면(=) 참이다*/
static int ccw(POINT A, POINT B, POINT C) {
    double a = A.x - B.x, b = A.y - B.y, c = C.x - B.x, d = C.y - B.y;
    return (a * d - b * c) <= 0; /* true if points A,B,C counterclockwise */
}
static int cmpl(const void *a, const void *b) {
    POINT *A = (POINT* )a, *B = (POINT* )b;
    int v = (A)->x - (B)->x ;       //lower-hull은 x 증가: y감소 순으로 정렬;
    if (v > 0) return 1; 
    if (v < 0) return -1;
    v = (B)->y - (A)->y ;
    if (v > 0) return 1; 
    if (v < 0) return -1;
    return 0;
}
static int cmph(const void *a, const void *b) {
    return cmpl(b, a);      // upper-hull은 x감소, y-증가 순으로 정렬;
}
static int makeChain(POINT *V, int n, int (*cmp)(const void*, const void*)) {
    qsort(V, n, sizeof(POINT), cmp);
    // 꼭짓점(i(>j))가 현재까지 만들어진 hull(s까지)의 에지보다 더 오른쪽에 
    // 있으면(ccw(i,j,j-1)), 왼쪽에 위치할 때까지 기존 hull을 줄인다. 
    // 이 점을 hull에 추가.
    int s = 1;
    for (int i = 2; i < n; i++) {
        int j = s;
        while (j >= 1 && ccw(V[i], V[j], V[j - 1])) j--; // ccw의 =때문에 일직선상의 
                                                         // 마지막 점만 들어감;
        s = j + 1; // j = last index of hull;
        // (j + 1)에 새 hull-point을 추가;
        POINT t = V[i]; V[i] = V[s]; V[s] = t;         // swap: V[i] <-> V[j + 1]
    }
    return s;
}
int chainHull(std::vector<POINT>& P, std::vector<POINT>& chull) {
    int n = P.size();
    if (n < 3) return 0;
    // copy to working array;
    chull.resize(n + 1);
    for (int i = 0; i < n; i++) chull[i] = P[i];
    /* make lower hull */
    int u = makeChain(&chull[0], n, cmpl);      
    /* make upper hull */   
    chull[n] = chull[0];
    u += makeChain(&chull[u], n - u + 1, cmph); 
    chull.resize(u);
    return u;
};

 

 

 

 

'Computational Geometry' 카테고리의 다른 글

Point in Polygon  (2) 2020.12.14
Incremental Delaunay Triangulation  (1) 2020.12.01
Chain Hull  (2) 2012.09.16
Quick Hull  (2) 2012.09.16
Monotone Cubic Interpolation  (0) 2012.09.09
Brute Force Convex Hull Algorithm  (0) 2012.09.06
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2020.12.30 15:10 신고  댓글주소  수정/삭제  댓글쓰기

    C++ STL에는 계산 기하학에서 유용하게 활용할 수 있는 <complex> 클래스가 있어요!
    어려워질 수록 struct로 정의하는 것이 나아보이기는 하지만,
    대부분의 경우에 간단하게 구현할 수 있어 추천드립니다.

  2. hgmhc 2020.12.30 15:11 신고  댓글주소  수정/삭제  댓글쓰기

    chain hull을 해주셨으니...
    rotating calipers를 막연히 기다려보겠습니다...

Quick Hull

Computational Geometry 2012. 9. 16. 02:07
728x90

파란선: 직전 convex hull을 표시함

Input = a set S of n points 
    Assume that there are at least 2 points in the input set S of points
QuickHull (S) { 
    // Find convex hull from the set S of n points
    Convex Hull := {} 
    가장 왼쪽(A)과 오른쪽(B)에 있는 두 점을 찾아 convec hull에 추가;
    선분 AB을 기준으로 나머지 (n-2)개 점들을 두 그룹 S1과 S2로 나눔;
        S1 = 선분 AB의 오른쪽에 놓인 점; 
        S2 = 선분 BA의 오른쪽에 놓인 점;
    FindHull (S1, A, B) 
    FindHull (S2, B, A) 
};
FindHull (Sk, P, Q) { 
   If Sk has no point, 
        then  return. 
   선분 PQ에서 가장 멀리 떨어진 점(C)을 찾아 convex hull에서 P와 Q의 사이에 추가;
   세점 P, C, Q는 나머지 Sk의 나머지 점들을 세 그룹S0, S1, and S2 
        S0 = 삼각형 PCQ에 포함되는 점;
        S1 = 선분 PC의 오른쪽에 놓인 점;
        S2 = 선분 CQ의 오른쪽에 놓인 점; 
    FindHull(S1, P, C) 
    FindHull(S2, C, Q) 
}
Output = convex hull

참고: http://www.cse.yorku.ca/~aaw/Hang/quick_hull/Algorithm.html

 

// P가 line(A, B)의 완전히 왼쪽에 있는가?

static bool leftSide(CPoint P, CPoint A, CPoint B) ;

더보기
static bool leftSide(CPoint P, CPoint A, CPoint B) {
   int ax = P.x - A.x,  ay = P.y - A.y;
   int bx = B.x - A.x,  by = B.y - A.y;
   return (ax * by - ay * bx) >= 0;
};

double dist2ToLine(CPoint P, CPoint A, CPoint B);

더보기
double dist2ToLine(CPoint P, CPoint A, CPoint B) {
    double dx = (B.x - A.x), dy = (B.y - A.y);
    double den = (dx * dx + dy * dy);
    if (den == 0.) return 0; //degenerate;
    double du = (P.x - A.x), dv = (P.y - A.y);
    double dist = (du * dy - dv * dx);
    return dist * dist / den; //distance-square!
};

#define SWAP_POINT(A, B)

static void findHull(CPoint *S, int n, CPoint A, CPoint B, std::vector <CPoint>& hull) ;

더보기
static void findHull(CPoint *S, int n, CPoint A, CPoint B, std::vector<CPoint>& hull) {
    if (n < 1) return ;
    //Find furtherst point 
    double maxdist = 0;
    for (int i = 0; i < n; i++) {
        double dist = dist2ToLine(S[i], A, B);
        if (dist > maxdist) {
            SWAP_POINT(S[0], S[i]);
            maxdist = dist;
        }
    }
    if (maxdist == 0.) return ; //collinear case;
    hull.push_back(S[0]);
    int j = 1;
    for (int i = 1; i < n; i++) {
        if (leftSide(S[i], A, S[0])) {
            SWAP_POINT(S[i], S[j]);
            j++;
        };
    };
    int k = j;
    for (int i = j; i < n; i++) {
        if (leftSide(S[i], S[0], B)) {
            SWAP_POINT(S[i], S[k]);
            k++;
        }
    };
    //1,...,j-1;
    findHull(&S[1], j-1, A, S[0], hull);
    //j,...,k-1;
    findHull(&S[j], k-j, S[0], B, hull);
};

void findMaxMin(std::vector<CPoint>& pts) ;

더보기
void findMaxMin(std::vector<CPoint>& pts) {
    int minx = pts[0].x, maxx = minx;
    int minid = 0, maxid = 0;
    for (int i = 0; i < pts.size(); i++) {
        if (pts[i].x > maxx)
        if (pts[i].x < minx)
    }
    SWAP_POINT(pts[0], pts[minid]); if (maxid == 0) maxid = minid;
    SWAP_POINT(pts[1], pts[maxid]);
};

void hullToPolyline(std::vector<CPoint>& hull);

void quickHull(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
    if (pts.size() < 3) return;
    hull.clear();
    //Find left and right most points, say A & B, and add A & B to convex hull ;
    findMaxMin(pts);
    hull.push_back(pts[0]);
    hull.push_back(pts[1]);
    int j = 2;
    for (int i = 2; i < pts.size(); i++) {
        if (leftSide(pts[i], pts[0], pts[1])) {
            SWAP_POINT(pts[i], pts[j]);
            j++;
        }
    }
    //2,3,...,j-1;
    findHull(&pts[2], j-2, pts[0], pts[1], hull);
    //j,j+1,...,pts.size()-1;
    if (j < pts.size()) // in order to avoid dereferencing &pts[pts.size()];
        findHull(&pts[j], pts.size()-j, pts[1], pts[0], hull);
    hullToPolyline(hull);
};

출력 hull은 단순히 convex hull을 구성하는 정렬이 안된 점들만 주므로, hull의 에지를 구하고 싶으면 추가적인 수정이 필요함.

더보기
static int cmph(const void *a, const void *b) {
    CPoint *A = (CPoint *)a , *B = (CPoint *)b ;
    int v = (A->x - B->x) ;
    if (v > 0 ) return 1;
    if (v < 0 ) return -1;
    v = B->y - A->y ;
    if (v > 0 ) return 1;
    if (v < 0 ) return -1;
    return 0;
};
static int cmpl(const void * a, const void *b) {
    return cmph(b, a);
};
void hullToPolyline(std::vector<CPoint>& hull) {
    CPoint A = hull[0];
    CPoint B = hull[1];
    // 가장 멀리 떨어진 두 점(hull[0], hull[1])을 연결하는 직선을 기준으로 프로젝션을 구하여서 순차적으로 
    int k = 2;
    for (int i = 2; i < hull.size(); i++) {
        if (leftSide(hull[i], A, B)) {
            SWAP_POINT(hull[i], hull[k]); k++;
        };
    };
    // k-1; last index of hull left side of line(A,B);
    // upper part reordering: 
    qsort(&hull[0], k, sizeof(CPoint), cmph);
    //lower part reordering;
    if (k < hull.size())
        qsort(&hull[k], hull.size()-k, sizeof(CPoint), cmpl);
    }
};

또한 입력점들의 순서를 그대로 유지하고 싶으면, double pointer를 이용하거나 복사복은 이용하여야 한다.

 

'Computational Geometry' 카테고리의 다른 글

Incremental Delaunay Triangulation  (1) 2020.12.01
Chain Hull  (2) 2012.09.16
Quick Hull  (2) 2012.09.16
Monotone Cubic Interpolation  (0) 2012.09.09
Brute Force Convex Hull Algorithm  (0) 2012.09.06
Curvature of Curve  (0) 2012.08.07
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2020.12.30 15:12 신고  댓글주소  수정/삭제  댓글쓰기

    quick hull이 평균 O(n log n)이였던가요?

728x90

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

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

void rasterize_voronoi(POINT vorocenter [], int N, BYTE *image, int width, int height) {
    // RGB 컬러이미지이므로 image는 적어도 (3 * width * height) 크기의 버퍼를 가져야 한다.
    // make color lookup table;
    DWORD *color = new DWORD [N]; //RGBQUAD;
    // 랜덤 컬러 LUT;
    for (int i = 0; i < N; i++) 
        color[i] = (rand() % 256) << 16 | (rand() % 256) << 8 | (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; //caution overflow;
            for (int k = 0; k < N; k++) {
                int dx = x - vorocenter[k].x;
                int dy = y - vorocenter[k].y;
                int ds2 = dx * dx + dy * dy;
                if (ds2 < dist2_min) {
                    dist2_min = ds2;
                    min_id = k;
                }
            }
            // set pixel value(here,we use RGB-color) wih verocenter[i]'s id
            *image++ = (color[min_id] >> 16) & 0xFF;   //B;
            *image++ = (color[min_id] >> 8)  & 0xFF;   //G;
            *image++ = (color[min_id])       & 0xFF;   //R;
        }
    };
    // draw cell center;
    delete[] color ;
};

사용자 삽입 이미지

 

'Image Recognition' 카테고리의 다른 글

EM Algorithm: Line Fitting 예  (0) 2008.06.29
Gaussian Mixture Model  (2) 2008.06.07
Rasterizing Voronoi Diagram  (0) 2008.05.26
RANSAC Algorithm  (0) 2008.05.24
Contour Tracing  (0) 2008.05.22
Gausssian Scale Space  (0) 2008.05.22
Posted by helloktk

댓글을 달아 주세요

728x90

단순 다각형(simple polygon)의 삼각화 알고리즘의 다른 구현이다. 이전 구현과 같은 Ear-Clipping 알고리즘을 써서 구현한 것이다. 폴리곤의 각 꼭짓점을 돌면서 현재 꼭짓점의 직전과 직후의 꼭짓점들로 구성한 삼각형이 Ear인가를 체크해서 맞으면 이 삼각형을 추가하고, 현 꼭짓점은 주어진 폴리곤에서 제거를 한 후 나머지에 대해 테스트를 계속한다. $N$ 개의 꼭짓점을 가지는 단순 폴리곤의 경우에 $N-2$ 개의 삼각형으로 분해가 된다.

주어진 폴리곤이 시계방향으로 정렬된 경우는 반시계 방향으로 바꾸어서 정렬시킨 후 작업을 한다. 구현된 알고리즘의 복잡도는 ${\cal O}(N^2)$이다. 

last update: 2012.10.03;

static bool leftSide(CPoint *P, CPoint *A, CPoint *B) ;
더보기
static bool leftSide(CPoint *P, CPoint *A, CPoint *B) {
    // P가 선분(A,B)의 왼쪽 또는 위에 있는가?
    return ((B->x - A->x) * (P->y - A->y) - ( B->y - A->y) * (P->x - A->x)) >=0;
}
static bool insideTriangle (CPoint *P, CPoint *A, CPoint *B, CPoint *C) ;
더보기
static bool insideTriangle (CPoint *P, CPoint *A, CPoint *B, CPoint *C) {    
    // 반시계방향으로 정렬된 삼각형(A,B,C)의 안에 점 P가 있는가?(경계포함);
    if (!leftSide(P, A, B)) return false;
    if (!leftSide(P, B, C)) return false;
    if (!leftSide(P, C, A)) return false;
    return true;
};​
static bool earTest(int a, int b, int c, int nv, int *V, CPoint *points) ;
더보기
static bool earTest(int a, int b, int c, int nv, int *V, CPoint *points) {
    // Check the triangle {V[a], V[b], V[c]} is an Ear; 
    CPoint *A  = &points[V[a]]; 
    CPoint *B  = &points[V[b]]; 
    CPoint *C  = &points[V[c]]; 
    // colinear(A->B->C) or concave인가를 체크; 반시계 방향으로 정렬된 입력점이므로 왼편에
    // 있으면 concave이므로 ear가 안된다;
    if (leftSide(B, A, C)) return false; 
    // 이미 만들어진 삼각형의 꼭지점이 포함되는지 체크한다;
    for (int k = 0; k < nv; ++k) {
        if ((k == a) || (k == b) || (k == c)) continue;
        if (insideTriangle(&points[V[k]], A, B, C)) return false;
    }
    return true;
};
static double polygonArea2(std::vector<CPoint>& pts) ;
더보기
static double polygonArea2(std::vector<CPoint>& pts) {
    double area2 = 0.;
    for (int p = pts.size() - 1, q = 0; q < pts.size(); p = q++)
        area2 += pts[p].x * pts[q].y - pts[p].y * pts[q].x;
    return (area2);
}
/* Polygon should be simple(ccw or cw);*/
int polygonTriangulation(std::vector<CPoint>& points, std::vector<Triangle>& triplets) {
    int nv = points.size();
    if (nv < 3) return 0;
    triplets.clear(); triplets.reserve(nv - 2);
    std::vector<int> V(nv);  // contains vertex indices;
    // 주어진 단순폴리곤을 반시계방향으로 정렬;
    if (polygonArea2(points) > 0)
        for (int v = 0; v < nv; ++v) V[v] = v;
    else 
        for (int v = 0; v < nv; ++v) V[v] = nv - 1 - v;
    // (nv-2)개 꼭지점을 차례로 제거한다. 한개의 꼭지점이 제거될때마다 
    // 삼각형이 하나씩 생긴다;
    int err_detect = 2 * nv;   /* error detection */
    for (int v = nv - 1; nv > 2;  ) {
        // 단순폴리곤이 아닌 경우에 무한루프를 도는 것을 방지;
        if (0 >= (err_detect--)) {
            TRACE("triangulate::ERROR\n");
            return 0;
        }
        // <u,v,w>는 연속하는 세 꼭지점의 인덱스임;
        int u = v % nv; 
        v = (u + 1) % nv;
        int w = (v + 1) % nv;
        if (earTest(u, v, w, nv, &V[0], &points[0])) { 
            triplets.push_back(Triangle(points[V[u]], points[V[v]], points[V[w]]));  
            // 꼭지점 V[v]를 제거함;
            for (int k = v; k < nv - 1; ++k) V[k] = V[k + 1];
            --nv;       
            /* resest error detection counter */
            err_detect = 2 * nv;
        }
    }
    return triplets.size(); // number of triangle = (nv-2);
};

최적화 후:

'Computational Geometry' 카테고리의 다른 글

Wu's Line Algorithm  (1) 2008.06.02
Brute Force Triangulation  (1) 2008.05.28
Polygon Triangulation (II)  (0) 2008.05.26
Polygon Triangulation  (4) 2008.05.25
Polygon Fill  (0) 2008.05.22
Fortune's Sweep Algorithm  (0) 2008.05.22
Posted by helloktk

댓글을 달아 주세요

728x90

코드 제거: 2012년 9월 4일;  코드 구현: http://kipl.tistory.com/13

Polygon의 내부를 겹치지 않게 분할하는 것을 polygon의 삼각화라고 한다. N개의 꼭짓점이 있는 polygon의 경우에 N-2개의 서로 겹치지 않은 삼각형을 내부에 가지게 되며, polygon의 경계와 겹치지 않는  N-3개의 내부 경계 라인을(diagonal)을 가지게 된다. 

사용자 삽입 이미지


분할된 삼각형의 외접원 반지름에 편차가 심한 경우 삼각형의 면적이 균일하게 분할되지 않는다. 이런 경우 인접하는 두 개의 삼각형으로 형성된 사변형에서 현재의 대각 변에 교차되는 대각 변으로 재분할을 시도한다. 이렇게 새로이 만든 삼각형을 검사하여 이전보다 면적이 균일하면 이 분할을 사용한다. 물론 삼각형의 외접원의 반경이 작을 때는 거의 균일한 분할을 얻을 수 있다. Polygon Optimizing을 참고하기 바란다.

*simple 폴리곤이 아닌 경우는 아래의 링크를 참조하기 바란다.

1). O(n log(n)) 복잡도를 갖는 polygon triangulation algorithm;
==> 'Fast Polygon Triangulation based on Seidel's Algorithm'.
==> http://www.cs.unc.edu/~dm/CODE/GEM/chapter.html

2). Poly2Tri 홈페이지:
==> Fast and Robust Simple Polygon Triangulation With/Without Holes by Sweep Line Algorithm
==> http://www.mema.ucl.ac.be/~wu/Poly2Tri/poly2tri.html

 

 

triangulation notes: (몇 개의 사이트는 없어졌다).
Nice lecture notes:
http://arachne.ics.uci.edu/~eppstein/junkyard/godfried.toussaint.html
▶ Narkhede & Manocha's description/code of Seidel's alg:
http://www.cs.unc.edu/~dm/CODE/GEM/chapter.html
▶ Some school project notes w/ triangulation overview & diagrams:
http://www.mema.ucl.ac.be/~wu/FSA2716-2002/project.html
▶ Toussaint paper about sleeve-following, including interesting
description & opinion on various other algorithms:
http://citeseer.ist.psu.edu/toussaint91efficient.html
▶ Toussaint outline & links:
http://cgm.cs.mcgill.ca/~godfried/teaching/cg-web.html
http://geometryalgorithms.com/algorithms.htm
▶ History Of Triangulation Algorithms
http://cgm.cs.mcgill.ca/~godfried/teaching/cg-projects/97/Thierry/thierry507webprj/complexity.html
▶ Ear Cutting For Simple Polygons
http://cgm.cs.mcgill.ca/~godfried/teaching/cg-projects/97/Ian//cutting_ears.html
▶ Intersections for a set of 2D segments
http://geometryalgorithms.com/Archive/algorithm_0108/algorithm_0108.htm
▶ Simple Polygon Triangulation
http://cgafaq.info/wiki/Simple_Polygon_Triangulation
▶ KKT O(n log log n) algo
http://portal.acm.org/citation.cfm?id=150693&dl=ACM&coll=portal&CFID=11111111&CFTOKEN=2222222
▶ Poly2Tri implemenation, good notes and looks like good code, sadly the
license is non-commercial only:(최근에 사이트가 존재하지 않음)
http://www.mema.ucl.ac.be/~wu/Poly2Tri/poly2tri.html
▶ FIST
http://www.cosy.sbg.ac.at/~held/projects/triang/triang.html
▶ Nice slides on monotone subdivision & triangulation:
http://www.cs.ucsb.edu/~suri/cs235/Triangulation.pdf

▶ Interesting forum post re monotone subdivision in Amanith:
http://www.amanith.org/forum/viewtopic.php?pid=43
 

'Computational Geometry' 카테고리의 다른 글

Brute Force Triangulation  (1) 2008.05.28
Polygon Triangulation (II)  (0) 2008.05.26
Polygon Triangulation  (4) 2008.05.25
Polygon Fill  (0) 2008.05.22
Fortune's Sweep Algorithm  (0) 2008.05.22
Triangulating Monotone Polygon  (0) 2008.05.22
Posted by helloktk

댓글을 달아 주세요

  1. 이상민 2008.05.26 12:55  댓글주소  수정/삭제  댓글쓰기

    잘 보았습니다. 이 분야로 졸업논문을 쓰고 있는데 많은 도움이 되었습니다.

  2. 이동기 2009.02.05 16:27  댓글주소  수정/삭제  댓글쓰기

    좋은글 공유해주셔서 감사합니다 잘 보고있습니다^^

    그런데 궁금한게 있습니다.

    이 알고리즘이 시계방향으로 그릴때와 반시계방향으로 그릴때 결과가 다르게 나오나요??

    저만 그런건지;;; 아니면 원래 그런건지 궁금해서요;;

    아. 그리고 마지막 점은 첫번째 점하고 같아야 하나요?

  3. helloktk 2009.02.07 00:09  댓글주소  수정/삭제  댓글쓰기

    코드에서 comment했듯이 위의 코드는 반시계방향의 입력에 대해서만 제대로 작동합니다(CCW 함수를 그렇게 만들었습니다.--> 물론 윈도우 그래픽은 y-방향이 위-아래가 우리가 보통 쓰는 좌표와 반대로 되어 있으므로 윈도우 그래픽상으로는 시계방향으로 나타났니다). 코드를 좀 수정하면 시계방향의 입력 폴리곤은 항상 다시 반시계방향으로 순서를 뒤집어서 만들수 있습니다.


    처음점과 마지막 점은 같지 않습니다. 같은 점(중복점)이 있으면 문제가 발생할 수 있습니다-->이것도 사전에 제거 되어야 합니다.

    제대로 할려면

    1. 입력 폴리곤의 중복점 제거
    2. 입력 폴리곤이 시계방향인지 반시계방향인지를 먼저 체크.
    3. 시계방향이면 점의 배열의 순서를 역으로 함.
    4. ear-clipping 알고리즘 적용

    이런 프로시져를 만들어야 합니다.

  4. YJ 2009.09.29 14:23  댓글주소  수정/삭제  댓글쓰기

    좋은 글 감사합니다. 제 작업에 많은 도움이 되었습니다.

728x90

 

 

Fortune2.zip
다운로드
사용자 삽입 이미지

 

 

사용자 삽입 이미지

 

평면에서 보로노이 다이어그램을 계산하는 알고리즘 중에 1980년대에 Fortune에 의해서 발견이 된 Sweep line 알고리즘에 대해서 살펴보겠다. 이 알고리즘은 보로노이 다이어그램 계산문제의 optimal 한 답 중에 하나이다(~O(n * log(n)). Fortune이 자신의 알고리즘을 직접 C로 구현한 것이 있는데 이것은 웹상에서 쉽게 구할 수 있다. 그러나 알고리즘의 얼개는 간단하나 구현에 들어간 자료구조 등을 간단히 않아서 코드를 이해하는 것이 쉽지 않았다(이에 비하면 incremental delaunay triangulation은 매우 단순하다). 몇 번을 간단한 자료구조를 이용해서 재 작성을 해보려고 했는데 잘 되지 않았다. 인터넷에서 구할 수 있는 다양한 구현 소스를 분석해 보았으나 쉬운 것을 찾기가 힘들었다(triangle.c에서도 구현되어 있고, 자바로도 구현이 된 것을 찾을 수 있다.) 원래의 목적은 Fortune의 코드가 메모리 처리를 잘못하여 메모리 leak이 발생하여서 이 문제를 해결하려고 시도한 것인데, 메모리 문제는 자바의 경우에 처럼 가비지 콜렉터를 만들어서 처리할 수는 있다. 이 알고리즘의 구현에 들어가는 기본적인 자료구조는 priority queue와 balanced binary tree이다. 이들의 기본적인 구현은 이미 잘 알려져 있으므로 이것들을 직접적으로 알고리즘에 적용하는 문제만 남는다. balanced binary tree로 만들어진 자료구조를 쓰면 자료를 찾는 시간이 O(log(n))의 시간 복잡도를 주지만, 이것이 알고리즘을 직관적으로 보는 것을 방해하므로 여기서는 이중 연결 리스트를 이용하도록 한다. 전체적으로 알고리즘은 O(n^2)의 복잡도를 가진다.

1. priority queue 구성: 주어진 입력점들을 가지고 구성한다. 알고리즘 중간에 세 점으로 원이 구성이 되는 circle event가 생성이 되는데 이것은 따로 event 큐를 구성한다. event의 우선순위는  x 좌표의 순서에 의해서 구성하고 동일한 x 좌표값에 대해서는 y값이 작은 순서대로 구성한다(compare 구조체에서 point와 event에 대한 비교 연산자를 정의한다).

std::priority_queue <point, std::vector <point>, compare> points ; // 입력점(point-구조체)들로
                                                                   // 만들어지는 priority_queue;
std::priority_queue <event*, std::vector <event*>, compare> events;//circle events(event-구조체);

2. site event 처리: sweep line(x=const)은 전체 평면을 반으로 나누는 역할을 한다. 알고리즘은 sweep 라인이 쓸고 지나간 영역에서만 관심 영역이다. sweep라인과 주어진 점은 하나의 포물선을 정의할 수 있다(주어짐 점이 포물선의 초점을 구성함). 따라서 sweep 라인이 진행하면 sweep line의 왼편의 입력점들은 각각 하나의 포물선을 형성하고, 이 포물선들의 구간 중에서 sweep라인에 가장 가까운 부분은 포물선 arc로 연결인 되는 하나의 beach line을 형성한다.  sweep라인이 새로운 입력점을 지나면 여기서 새로운 포물선 arc가 생기는데, 이것은 이미 만들어진 beach라인을 구성하는 포물선 arc 중 하나를 둘로 가르고, 그 자신이 새로운 beach 라인의 일부를 구성하게 된다; sweep라인이 전진하면서 비치라인을 구성하는 arc가 다른 arc에 파 묻혀서 없어지는 경우가 생긴다. 이 경우가 circle event가 만들어지는 시점이다.

.......................

*           BEFORE                                  AFTER
*
*           new point                               new arc
*               |                                       |
*     __        |       _____                   __      |       _____
*       \      \|/        |                       \    \|/       arc a
*        \      `         |                        \    `       __|__
*         \     X         |                         p------X    __|__<-- arc c
*      a   |             arc a                  a    |            |
*         /               |                         /            arc a
*        /                |                        /              |
*      /__              __|__                     /__           __|__
*   __/    \              |                    __/   \            |
*           \             |                           \           |
*            \            |                            \          |
*       b     |          arc b                  b       |        arc b
*            /            |                            /          |
*           /             |                           /           |
*        __/            __|__                      __/          __|__

 

3. circle event 처리; sweep라인이 각각의 site에 도달하면 새로운 포물선의 arc가 비치라인에 추가가 되는데 sweep 라인으로부터 멀어진 arc들은 어느 시점에서 없어져야 한다. 이것은 인접하는 세 개의 arc들의 교점이 현재의 sweepline 위치에서 하나로 만들어져서 중앙의 arc가 없어지는 경우로 이것을 circle event라고 한다. 이때 만들어지는 교점은 보로노이 에지의 꼭짓점을 구성하게 된다. circle event에서는 나머지 두 개의 남은 arc의 교점이 trace 하는 새로운 에지가 생기게 된다.

*    __
*       \
*   a    \
*         \
*   __     |
*     `   /
*      \ /
*   b   X           <-- Arc b is overtaken at point X (this is a circle center)
*     / \
*   __,   \
*          \
*           \
*   c        |
*           /
*          /
*         /

*
*       BEFORE                                              AFTER
*
*      \        arc a                                   \
*        \                                               \              a.s0
*         \   <-- a.s0                                    \               |
*          \                                               \             \|/
*           \                                               \             `
*            \                                               \
*    arc b   X     <-- termination point                      .--------------------<-- new segment
*           /                                                /          
*          /                                                /            .
*         /   <-- c.s1                                     /            /|\ 
*        /                                                /              |
*       /    arc c                                       /              c.s1
*      /      

                                 

//포물선 arc를 정의하기 위한 클래스;

struct arc {
    point p;                //focus of parabola;
    arc *prev, *next;       //double linked-list;
    event *e;
    ////////////////////////////////////////////////////
    seg *s0 ;               //edge of voronoi starting from break point;
    seg *s1;                //edge of voronoi starting from break point;
    arc(point _p, arc *_a=0, arc *_b=0)
        : p(_p), prev(_a), next(_b), e(0), s0(0), s1(0) {}
};

// voronoi edge를 정의하는 segment class;
struct seg {
    point start, end;                           //defines segment;
    bool done;
    int ref ;                                   //referece count in order to distingush the double edges;
    seg(point _p, int _ref=0)
        : start(_p), end(0,0), done(false), ref(_ref)
    { output.push_back(this); }                 //garbage collector;
    void finish(point p) { if (done) return; end = p; done = true; }
};

//메인 wrapper;

int fortune_main(std::vector<POINT>& pts,          //voronoi points;
                       std::vector<EDGE>& edges)     //voronoi edges;

{
      //point priority_queue구성;
      std::vector<POINT>::iterator iter=pts.begin();
      for(; iter !=pts.end(); iter++) {
          point P(iter->x, iter->y) ;
          points.push_back(P);
      }
      //body of algorithm ;
      while(!points.empty()) {
           if(!events.empty() && events.top()->x <= points.top().x) 
                 process_event() ; //circle event ;
           else 
                 process_site() ;   //site event;
      }
      //for remaining circle events if any;
      while(!events.empty()) 
           process_event() ;
                                                                      
      //2. prepare output edges and clean memory;
}

arc * root=NULL;                              //global variable;
void process_site() {
    point P= points.top(); points.pop(); //because popping return void in STL;
    if(!root) { // root of double linked list for arcs(global);
        root = new arc(P) ; return  ;
    }
    //
    for(arc* a = root; a; a=a->next) {
         point Q;
         if(intersect(P, a, &Q)) { //arc a와 점P에서의 포물선(degenerate 된)이 만나는 점 Q;
              duplicate_arc(P, a) ; //교차하는 arc를 둘로 만든다 
                                    //(다음 arc가 있고, 만약에 이것과도 교차하면 circle event 임으로 제외)
              insert_new_arc(P, a, a->next) ;//새 arc를 중간에 삽입한다;
              //새 arc의 에지 세팅 ::교차점을 출발점으로 하는 두개의 반직선을 만든다: 
              //도착점은 circle event에 대부분 결정되고, 나머지는 후처리 과정에서 결정됨)
              a = a->next ;
              a->prev->s1 =a->s0 = new seg(Q, ref_count) ;
              a->next->s0 =a->s1 = new seg(Q, ref_count++) ; //쌍으로 생성되는 반직선을 identify하기 
                                                             //위해서 동일번호를 부여함.
              //새 site 추가로 비치라인을 구성하는 arc의 초점들과 circle event를 만들 수 있는가 체크함;
              check_circle_event(a, P.x) ;
              check_circle_event(a->prev, P.x) ;
              check_circle_event(a->next, P.x) ;
         }
    }//for-;
};
//
void process_event(){
    event *e = events.top(); events.pop();
    if (e->valid) {
        // start a new edge.(single edges)
        seg *s = new seg(e->p, ref_count++);
        // remove the associated arc from the front. and attach a new segment;
        arc *a = e->a;
        remove_arc(a, s);
        // finish the edges before and after a==>new voronoi vertex;
        if (a->s0) a->s0->finish(e->p);
        if (a->s1) a->s1->finish(e->p);        
        // recheck circle events on either side of p:
        if (a->prev) check_circle_event(a->prev, e->x);
        if (a->next) check_circle_event(a->next, e->x);
        delete a;
    }
    delete e; 
};
// Look for a new circle event for arc a.
void check_circle_event(arc *a, double x0/*=current sweep-line*/) { 
    // Invalidate any old event.
    if (a->e && a->e->x != x0)  a->e->valid = false;
    a->e = NULL;
    if (!a->prev || !a->next) return;
    double maxx;
    point O;
    point& A = a->prev->p ;
    point& B = a->p ;
    point& C = a->next->p ;
    //collinear가 아니고, 최대값이 현재의 site event위치보다도 큰 경우에;
    if (circle(A, B, C, &maxx, &O) && maxx > x0) { //점 A,B,C에 의해서 형성이 되는 
                                                   //원 중심(O)과 최우측x(maxx) 좌표를 얻음;
        // create new event.
        a->e = new event(maxx, O, a);
        events.push(a->e);
    }
};

나머지 함수들은 모두 단순한 구현이므로 생략한다;

보로노이 에지는 seg collector에서 각각의 segment를 끄집어내어서 그리면 된다. 그러나 각각의 segment는 보로노이 에지를 전체를 커버하는 것이 아니다. site event인 경우에는 항상 듀얼로 반직선을 생성하는데. 이 두 개의 반직선이 하나의 보로노이 에지를 정의한다(따라서 ref를 참조하면 온전한 하나의 에지를 찾을 수 있다). circle event의 경우에는 에지가 듀얼로 생성하지 않았으므로 이 경우에는 하나의 segment가 그 에지를 표현한다. 에지 액세스를 쉽게 하기 위해서는 모든 에지를 듀얼 구조로 만들어서 사용할 수 있다.

n 개의 점들의 보로노이 다이어그램은 얼마나 많은 꼭짓점과 에지를 가질까?

이것은 오일러 공식 V-E+F=2를 사용하면 된다. 보로노이 다이어그램은 경우 바깥쪽의 에지들은 무한대로 연결이 되어 있다. 이것을 무한대에서 가상의 꼭짓점을 가정하고 그것에 연결이 되어 있다고 생각하면 된다. 따라서 V --> V+1로 계산을 해야 한다. 오일러 공식 : V+1-E + F= 2; 여기서 F=n임을 알 수 있다. (n개의 점들이 하나의 face상에 놓여 있음)그리고 하나의 에지가 두 개의 꼭짓점을 연결하므로 각각의 꼭짓점에서 나간 에지의 합(=deg(v))을 계산하면 2*E 개 임을 알 수 있다. sum deg(v) = 2 * E;그런데 각각의 꼭짓점에서는 적어도 3개 이상의 에지와 연결이 되어 있으므로 위식의 좌변은 (V+1) * 3 <= 2 * E;를 준다. 따라서 오일러 공식과 이 부등식을 연관시키면 V <=  2*n - 5; E <=  3*n - 6;임을 알 수 있다. 즉 필요한 메모리의 양은 입력점의 수의 선형적으로 비례한다.

따라서 전체 이벤트의 개수는 site 이벤트와 꼭짓점에 해당하는 circle event 만큼이 있으므로 O(3*n) 정도이다. 각각의 event에 대해서 arc노드를 검색해야 하므로 O(n) 번 탐색을 해야 한다(balanced binary tree의 경우에는 log(n)). 따라서 알고리즘의 복잡도는 O(n^2 ( n*log(n))이다./**** http://blog.naver.com/helloktk/80041603288*/

** 첨부된 실행파일로 알고리즘을 테스트해 볼 수 있다(중복된 입력점은 제거하였고, 세 점이 일직선상에 놓인 것을 방지하기 위해서 작은 랜덤 값을 첨가하였다)

 

'Computational Geometry' 카테고리의 다른 글

Polygon Triangulation  (4) 2008.05.25
Polygon Fill  (0) 2008.05.22
Fortune's Sweep Algorithm  (0) 2008.05.22
Triangulating Monotone Polygon  (0) 2008.05.22
Trapezoidalization  (0) 2008.05.22
Optimizing Polygon Triangulation  (0) 2008.05.22
Posted by helloktk

댓글을 달아 주세요

728x90

(1) closest pair in the planar point set
     * brute force : O(n^2) --> O(n log n), see also divide and conquer algorithm.
     * ref : http://www.cs.mcgill.ca/~cs251/ClosestPair/ClosestPairPS.html 

사용자 삽입 이미지

(2) simplicity of polygon :

Shamos-Hoey Algorithm

 

 

 

사용자 삽입 이미지

(3) voronoi diagram :Fortune algorithm


 

 

사용자 삽입 이미지

(4) segment intersection : Bentley-Ottmann Algorithm

      * 파란색 점 : sweep-line과 교차하는 세그먼트들 간의(붉은색으로 표시)의 교차점으로(before sweep-line)로 point 이벤트 큐에 들어간다.
      * 청록색 점 : sweep-line이 지나간 후의 교차점.

사용자 삽입 이미지

sweep-line알고리즘을 구현하는데서 스텝별로 쪼개서 구현하려고 하니 조금 귀찮은 작업들이 많아진다 (DC에 그려지는 것을 자동으로 저장하는 부분을 만들어서 매 스텝마다 캡처하는 작업은 줄였다). 그러나 근본적인 문제는 입력되는 데이터에 degeneracy가 있는 경우다. 이것 때문에 가장 기본적인 세그먼트 intersection 부분부터 다시 살펴보아야 했다.

/**
** http://blog.naver.com/helloktk/80051980882 에서 옮긴 글.
*/

 

'Computational Geometry' 카테고리의 다른 글

Polygon Fill  (0) 2008.05.22
Fortune's Sweep Algorithm  (0) 2008.05.22
Triangulating Monotone Polygon  (0) 2008.05.22
Trapezoidalization  (0) 2008.05.22
Optimizing Polygon Triangulation  (0) 2008.05.22
Sweep Line Algorithm  (0) 2008.05.22
Posted by helloktk

댓글을 달아 주세요