'algorithm'에 해당되는 글 7건

  1. 2010.03.24 Savitzky-Golay smoothing filter (2)
  2. 2008.05.28 Brute Force Triangulation
  3. 2008.05.26 Polygon Triangulation (II)
  4. 2008.05.25 Polygon Triangulation (4)
  5. 2008.05.24 RANSAC Algorithm
  6. 2008.05.22 Polygon Fill
  7. 2008.05.22 Fortune's Sweep Algorithm

The Savitzky–Golay method essentially performs a local polynomial regression (of degree k) on a series of values (of at least k+1 points which are treated as being equally spaced in the series) to determine the smoothed value for each point.
""

Savitzky–Golay 필터는 일정한 간격으로 주어진 데이터들이 있을 때(이들 데이터는 원래의 정보와 노이즈를 같이 포함한다), 각각의 점에서 주변의 점들을 가장 잘 피팅하는 k-차의 다항식을 최소자승법을 이용해서 찾아서 그 지점에서의 데이터 값을 결정하는 필터이다. 이 필터는 주어진 데이터에서의 극대나 극소, 또는 봉우리의 폭을 상대적으로 잘 보존한다.(주변 점들에 동등한 가중치를 주는 Moving Average Filter와 비교해 볼 수 있다).

간단한 예로, 2차의 다항식과 5개의 데이터 점

{(-2, d0), (-1, d1), (0, d2), (1, d3), (2, d4)}

을 이용해서 중앙에서의 값을 결정하는 방법을 살펴보자. 사용하려고 하는 다항식은

p[x] = a0 + a1*x + a2*x2

이다. 다항식의 계수들은 다항식의 값과 실제 데이터의 값과의 차이를 최소화시키는 것들로 잡아야 한다.. 즉, 최소자승의 원리를 적용하여서 구하면 된다. 계산된 다항식의 값과 실제 데이터 값 사이의 차의 제곱을 구하면:

L = |a0-2*a1+4*a2-d0|2+|a0 -a1+ a2-d1|2+|a0-d2|2+|a0+a1+a2-d3|2
                          
+|a0+2*a1+4*a2-d4|2

이 식의 a0, a1, a2에 대한 극값을 가질 조건은

5*a0 + 10*a2 = d0 + d1 + d2 + d3 + d4

10*a1 = -2d0– d1 + d3+ 2d4

10*a0 + 34*a2 = 4d0 + d1+ d3 + 4d4        

이 식을 만족시키는 a
0를 구하면, 필터의 출력값 (원도우 중앙에서 값)이 결정된다.

필터출력 = a0 = (-3*d0 + 12*d1 + 17*d2 + 12*d- 3*d4) / 35;

위에서 계수 a0, a1, a2를 결정하는 방정식은 행렬로 정리하면 아래의 식과 같이 표현할 수 있다.

""

좌변의 5행3열 행렬을 A, a = [a0, a1, a2]T, d = [d0, d1, d2, d3, d4]T로 놓으면, 이 행렬방정식은 A.a = d 형태로 쓸 수 있고, 양변에 AT를 곱해서 왼쪽을 대칭행렬로 바뀌면

(AT.A).a = AT.d

따라서 해는,

a = (AT.A)-1.(AT.d)

이 식은 임의의 k-차 다항식을 이용한 경우에도 사용할 수 있다. 이 경우에 행렬 (ATA) 는 (k+1)x(k+1)의 대칭행렬이 된다. 행렬 A는 다항식의 찻수와 피팅에 사용이 될 데이터의 구간의 크기가 주어지면 정해지므로,  윗 식에서 (AT.A)-1.AT의 첫행 (a0을 d로 표현하는 식의 계수들)을 구해서 프로그램상에서는 lookup table를 만들어서 사용할 수 있다. 아래표는 mathematica 를 이용해서 윈도우 크기가 7 (7개 점)인 경우 2차 다항식을 사용할 때 계수를 구하는 과정이다.

 

""

2차 다항식일 때, 같은 방식으로 다양한 윈도우 크기에 따른 계수를 구할 수 있다.
*크기(n)에 따른 필터값 결정계수 (중앙에 대해 좌우대칭이다);

n=5;  W[] = {-3, 12, 17, 12, -3};
n=7;  W[] = {-2, 3,  6, 7, 6, 3. -2};
n=9;  W[] = {-21, 14. 39, 54, 59, 54, 39, 14, -21};
n=11; W[] = {-36, 9, 44, 69, 84, 89, 84, 69, 44, 9, -33};

필터출력 =   ∑ W[i]d[i] / ∑ W[i]

저작자 표시 비영리 변경 금지
신고
크리에이티브 커먼즈 라이선스
Creative Commons License

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

Adaboost  (0) 2010.12.28
Blur Detection  (0) 2010.05.25
Savitzky-Golay smoothing filter  (2) 2010.03.24
Watershed Algorithm Code  (0) 2010.03.19
Retinex 알고리즘  (10) 2010.02.03
Gaussian Mixture Model & KMeans  (4) 2010.01.30
Posted by helloktk

주어진 점집합에서 세점을 뽑아서 이들이 이루는 삼각형의 외접원에 다른 점이 하나도 포함하지 않으면  triangulation의 기본삼각형 cell이 된다. nC3개의 주어진 점의 combination에 의해서 생성이 되는 삼각형을 나머지 점들에 대해서 incircle 테스트를 수행하기 떄문에 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)-->(x, y, x^2 + y^2);

의 매핑에 의해서 3차원의 파라볼로이드상의 점으로 옮길 수 있다. 파라볼로이드 상의 세점에 의해서 만들어지는 평면과 파라볼로이드 면이 만나는 곡선을 x-y  평면상으로 프로젝션하면 원이 된다는 것을 쉽게 알 수 있다.(incircle 포스팅 참조). 따라서 평면상에서 세점에 의해서 만들어진 원에 주어진 점이 포함된는가를 테스트하는것은 파라볼로이드 위의 점으로 변환했을 때((x, y)-> (x, y, x^2 + y^2)), 변환된 파라볼로이드 위의 세점이 만드는 평면에 주어진 점(파라볼로이드로 변환됨)이 포함된는가만 보면 된다. 

원의 내부점이면 파라볼로이드로 변환된 점은 평면의 아래에 놓이고, 원의 외부점이면 평면의 위에 놓이게 된다. 물론 원 상의 점은 평면에 놓인다. 따라서 평면의 수직벡터 구하고, 삼각형의 한 꼭지점을 기준한 주어진 점의 변위벡터와 내적을 구하면 내적의 값은 평면 위인지 아래인지 또는 평면상에 놓인 점인가에 따라서 부호가 달라진다. 평면의 수직벡터를 고정하면(예제는 아래방향: nz < 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 ii = 0; ii < n; ii++)
        z[ii] = x[ii] * x[ii] + y[ii] * y[ii] ;

    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;
}

사용자 삽입 이미지

신고
크리에이티브 커먼즈 라이선스
Creative Commons License

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

Wu's Line Algorithm  (1) 2008.06.02
삼각형 외접원의 Inclusion Test  (0) 2008.05.29
Brute Force Triangulation  (0) 2008.05.28
Polygon Triangulation (II)  (0) 2008.05.26
Polygon Triangulation  (4) 2008.05.25
Polygon Fill  (0) 2008.05.22
Posted by helloktk

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

알고리즘은 주어진 폴리곤이 시계방향인 경우에는 반시계방향으로 수정하여서 작업을 한다. 구현된 알고리즘의 복잡도는 O(N^2)이다. 

last update: 2012.10.03;

static bool leftSide(CPoint *P, CPoint *A, CPoint *B) {

더보기

static bool insideTriangle (CPoint *P, CPoint *A, CPoint *B, CPoint *C) {

더보기

static bool earTest(int a, int b, int c, int nv, int *V, CPoint *points) {

더보기

static double polygonArea2(std::vector<CPoint>& pts) {

더보기

/* 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) ;
    // 주어진 단순폴리곤을 반시계방향으로 정렬한다;
    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");
            triplets.clear();
            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);
}


최적화 후:





신고
크리에이티브 커먼즈 라이선스
Creative Commons License

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

삼각형 외접원의 Inclusion Test  (0) 2008.05.29
Brute Force Triangulation  (0) 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

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

폴리곤의 내부를 겹치지 않게 분할하는 것을 폴리곤의 삼각화라고 한다. N개의 꼭지점이 있는 폴리곤의 경우에 N-2개의 서로 겹치지 않은 삼각형을 내부에 가지게 되며, 폴리곤의 경계와 겹치지 않는  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
 

신고
크리에이티브 커먼즈 라이선스
Creative Commons License

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

Brute Force Triangulation  (0) 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

어떤 현상을 설명하는 이론 모델을 만들고자 하면 현상에서 관측데이터를 모아야한다. 그런데 관측 데이터에는 모델에 대한 잘못된 가정이나 측정장비의 오차등에서 생기는 여러가지 형태의 노이즈가 들어있게 마련이다. 이러한 노이즈는 모델을 정확히 예측하는 것을 방해한다. 이처럼 모델파라미터의 예측을 방해하는 데이터(outliner)가 들어 있는 관측데이터로부터 어떻게 하면 모델을 구축할 수 있는가라는 질문에 대한 한 대답을 RANSAC ("RANdom SAmple Consensus")  알고리즘이 제공한다. 이 알고리즘은 1981년 Fischler 과Bolles에 의해서 제안되었다.

RANSAC 알고리즘에 들어가는 입력은 관측된 데이터값과, 관측결과를 피팅하거나 설명할 수 있는 모델파라미터, 그리고 신뢰성을 나타내는 파라미터를 필요로 한다.

RANSAC알고리즘은 주어진 원본데이터에서 반복적으로 일부를 임의로 선택하는 과정을 반복하여서 최적의 파라미터를 예측하는 프로시져의 형태를 가진다. 전체의 관측데이터가 M개 있고, 모델파라미터를 예측하는데 N개의 데이터가 필요한 경우에, 알고리즘은 :

     1. 임의로 관측데이터에서 N개의 서브데이터를 선택한다.
     2. 이 데이터를 (가상의) inlier로 생각하고 모델을 예측한다.
     3. 원본데이터중에서 예측된 모델에 잘 맞는가를 체크한다. 잘 맞는 데이터 수를 K라고 한다.
     4. 만약에 K가 충분하면, 이 모델을 확정하고 알고리즘을 종료한다.
     5. 모델이 좋지 않으면 1-->4과정을  L 번  반복한다.
     6. 충분히 반복후에도 좋은 모델을 찾지 못하면 모델예측에 실패한 것으로 간주한다.

1. 임계값 K 는 모델이 얼마나 잘 맞는가에 대한 사용자입력으로부터 결정이 된다. 그러나 대략적으로 주어진 샘플에서 inlier의 비율 Pg 이라고 생각되는 정도를 잡으면 된다 :
 
   K = M * (1- P
g)

2. 얼마나 많은 반복을 해야하는가? 주어진 관측데이터에서 inlier일 확률이 Pg인 경우에 L번의 모델예측시도가 실패할 확률을 계산하여서 이것이 주어진 설정값, pfail 보다도 작은 경우에 모델예측의 실패로 처리하므로,

  pfail   = L번의 모델예측 실패확률
          = (한번 시도가 실패할 확률)L
          = (1 - 한번 시도가 성공할 확률))L
          = (1 - (랜덤데이터가 모델을 예측할 확률)N))L
          = (1- (Pg) N)L

이 사실로 부터 최대 반복횟수는 
    
   L = log(pfail)/log(1-(Pg)N)

로 주어진다.
inlier가 반이 섞인 경우 Pg (=주어진 데이터중에서 inlier일 확률)
= 0.5,
p
fail = 0.01 인 경우, 만약에, N = 3 (세 점만 있으면 모델구성이 가능한 원의 피팅이 한 예임) 일 때, 최대 반복횟수는 윗 식에 적용하면,
     
                   L = 35 회


RANSAC알고리즘은 주어진 관측데이터에 많은 outlier가 존재하더라도 매우 정확히 모델파라미터를 예측할 수 있는 능력이 있는 robust한 알고리즘이지만, 얼마나 많은 반복을 해야하는가에 대한 상한값이 제공되지 않는다(사용자가 pfail 를 설정한다). 따라서 사용자 설정에 의한 상한값은 최적의 모델이 아닐 수 있다.

"The RANSAC procedure is opposite to that of conventional smoothing techniques: Rather than using as much of the data as possible to obtain an initial solution and then attempting to eliminate the invalid data points, RANSAC uses as small an initial data set as feasible and enlarges this set with consistent data when possible".

Martin A. Fischler and Robert C. Bolles (June 1981). "Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography". Comm. of the ACM 24: 381–395.

예제코드:
ransac을 이용한 라인 피팅: http://blog.naver.com/helloktk/80029006029
ransac을 이용한 원 피팅: http://kipl.tistory.com/32

신고
크리에이티브 커먼즈 라이선스
Creative Commons License

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

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
Watershed Algorithm 적용의 예  (2) 2008.05.21
Posted by helloktk

MFC를 써서 폴리곤의 내부를 칠하기 위해서는 윈도우  GDI region을 다루는  CRgn  클래스의 CRgn::CreatePolygonRgn() 메쏘드를 써서 해당 폴리곤으로 bound된 영역을 만든 이후에 CDC::FillRgn() 함수를 이용하여서 원하는 색으로 칠할 수 있다. 직접  raster에 칠하려고 할 때는 이러한 방법이 너무 불편하다 (물론 라스터(DIB)를 윈도우 비트맵으로 변환한 후에 위의 방법대로 칠하고, 다시 윈도우 비트맵을 라스터로 바꾸는 작업을 하면 되는데 번거롭다).

쉽게 생각할 수 있는 방법은  point_in_polygon()  방법을 써서 각각의 픽셀이 주어진 폴리곤 내부점인지를 확인하면서 작업을 할 수 있다. 그러나 이 방법은 한 픽셀을 검사하기 위해서 주어진 폴리곤의 모든 에지를 다 건들어야 하는 낭비를 하게 된다. 보다 효율적인 방법은 point_in_polygon 알고리즘을 구현할 때 사용하는 원리를 이용하면 적어도 한 스캔라인의 모든 픽셀에 대해서는 딱 한번만 전체 폴리곤의 에지를 검사하게 만들 수 있다. inclusion 테스트는 해당점에서 출발하는 수평반직선이 폴리곤과 몇번 교차하는지를 세서 내부점인지 외부점인지를 판별한다. 이것을 생각하면 수평선이 폴리곤과 교차하는 점들을 모두 구하여서(한번 폴리곤의 전체 에지검사가 필요) 크기의 순서대로 정렬할 때, 짝수번째 구간(0->1, 2->3,.....)에 들어가는 스캔라인의 픽셀들은 모두 폴리곤의 내부점이 된다는 것을 알 수 있다.

사용자 삽입 이미지

알고리즘 순서:
(1) 각각의 스캔라인에 대해서 교차점들을 구한다:
       꼭지점-i와 꼭지점-j(=i-1)를 잇는 직선의 방정식: x = x[i] + (x[j] - x[i])*(y - y[i]) / (y[j] - y[i]) 이므로
       현재의 scanline 위치에서 y값이 두 꼭지점의 y값 (y[i], y[j]) 사이에 있으면, 교차점의 x-좌표는 윗식의
       좌변값이 된다.
(2) 교차점들을 크기의 순서대로 정렬한다.
(3) 짝수번째 구간에 들어가는 스캔라인상의 픽셀들은 해당 색깔로 칠한다.

/*코드 샘플
** 폴리곤을 플로팅타입으로 변환이 필요?
*/
// qsort용 비교함수;
static int comp(const void *a, const void *b) {
    return *(int*)a - *(int*)b;
};
void FillPolygon(POINT poly[], int npoly,
                 DWORD color,
                 BYTE *image, int width, int height, int bpp)
{
    int *nodeX = new int [npoly] ; //never exceeds npoly;
    for (int y=0; y<height; y++)
    {
        // find intersection nodes;
        int nodes=0;
        for (int i=0, j=npoly-1; i<npoly; j = i++)
        {
            if (poly[i].y<y && poly[j].y>=y || poly[j].y< y && poly[i].y>= y)
            {//수평인 에지는 그리지 않는다(부등식을 자세히 보라)
                nodeX[nodes++] = (int)(poly[i].x +
                                double(y-poly[i].y)*double(poly[j].x-poly[i].x)/double(poly[j].y-poly[i].y) +.5);
                                //round to integer!!.
            }
        }
        // sort nodes (ascending order);
        qsort(nodeX, nodes, sizeof(int), comp) ;

        //  fill the pixels between node pairs.
        for (i=0; i<nodes; i+=2)
        {
            if   (nodeX[i  ]>=width) break;
            if   (nodeX[i+1]> 0 )
            {
                if (nodeX[i  ]< 0 ) nodeX[i  ]=0 ;
                if (nodeX[i+1]> width) nodeX[i+1]=width;
                for (int x=nodeX[i]; x<nodeX[i+1]; x++)
                    SetPixel(x, y, color, image, width, height, bpp); //다른 프로토타입을 가질 수 있다.
            }
        }
    }
    delete[] nodeX;
}

사용자 삽입 이미지

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

신고
크리에이티브 커먼즈 라이선스
Creative Commons License

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

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
Trapezoidalization  (0) 2008.05.22
Posted by helloktk


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
*/

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


신고
크리에이티브 커먼즈 라이선스
Creative Commons License

'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


티스토리 툴바