728x90

일직선에 있지 않은 세 점 $A$, $B$, $C$가 만드는 외접원 (circumcircle)에 주어진 점 $D$가 포함되는가를 판별하는 것은 기하 알고리즘에서 매우 중요하다. 직접 외접원의 중심과 반지름을 구해서 판단을 할 수 있지만, 이 경우 수치 계산의 오류가 들어갈 여지가 매우 많다 (외접원의 중심과 반지름을 구하는 식을 만들어 보면 쉽게 알 수 있다). 여기서는 이런 문제에 부딪히지 않고 보다 robust 한 판단을 주는 식을 구해보자.

이는 공간 상의 네 점 $A(a_x, a_y, a_z)$, $B(b_x, b_y, b_z)$, $C(c_x, c_y, c_z)$, $D(d_x, d_y, d_z)$만드는 사면체의 (부호를 가지는) 부피를 구하는 공식이 다음 행렬식으로 주어진다는 사실에서 출발한다.

사용자 삽입 이미지

이 식을 이용하면 평면에서 3점의 외접원에 주어진 점이 포함이 되는가를 판별하는 식 역시 3차원 공간에서의 체적으로 표현이 가능하다.

$z= x^2 + y^2$로 표현이 되는 포물면(paraboloid)의 한 점 $(a, b, a^2 + b^2)$에서 접평면의 방정식을 구하면,

$$z=2ax + 2bx - (a^2 + b^2)$$

로 주어진다. 이 평면을 $z$-방향으로 $r^2$만큼 이동시키면

$$z=2ax + 2by -(a^2 +b^2) + r^2$$

로 주어지고, 이 평면이 포물면을 절단하는 곡선을 바닥에 정사영하면 중심이 $(a,b)$이고 반지름이 $r^2$인 원이 된다:

$$(x-a)^2 + (y-b)^2 = r^2$$

3점을 포물면 위로 올리면 3차원의 점이 되고, 이 3점이 만드는 평면과 포물면의 절단 곡선을 바닥으로 정사영하면 원래의 외접원이 된다. 따라서 외접원의 내부점은 평면 아래에, 외부점은 평면 위에 놓인다. 원주점은 평면에 놓인다. 동일한 (평면에 놓인) 바닥 삼각형에 대해서 꼭지점(테스트하려는 점)의 상대적 위치가 정반대에 놓이므로 행렬식으로 표현된 이들 4점이 만드는 사면체의 부호가 달라진다. 즉, 외접원 내부인가 외부인가를 포물면 위로 올려진 외접원 3점과 테스트하려는 점이 구성하는 사면체 체적의 부호로 판별이 가능하다

 

다시 정리하면, 세 점 $A(a_x, a_y)$, $B(b_x, b_y)$, $C(c_x, c_y)$, $D(d_x, d_y)$가 주어진 경우에 $A, B, C$(반시계 방향)에 의해서 형성이 되는 원에 점 $D$가 포함이 되는지의 여부는 공간상의 네 점

$$A = (a_x, a_y, a_x^2+a_y^2),$$

$$B = (b_x, b_y, b_x^2+b_y^2),$$

$$C = (c_x, c_y, c_x^2+c_y^2),$$

$$D = (d_x, d_y, d_x^2+d_y^2)$$

가 형성하는 체적의 부호에 의해서 결정이 된다.

사용자 삽입 이미지


위의 행렬식은 좀 더 간단히 정리할 수 있다. 1, 2, 3행에서 4행을 행단위로 빼도 값이 변하지 않으므로 오른쪽 식을 얻고(4열에 대해서 전개 후), 1열에 $2d_x$를 곱한 것과 2열에 $2d_y$를 곱한 것을 각각 3열에서 빼면 아래줄의 식을 얻는다.

사용자 삽입 이미지

이 3x3 행렬식은 쉽게 계산이 된다. $ccw(A, B, C)$를 계산하여 전체의 부호를 고정시킨다.

// A, B, C가 만드는 외접원에 D가 들어가는가?;
// det = 0인 경우는 세 점이 일직선이나, D가 원주에 놓인 경우다;
// 마지막 판단을 수정하여서 원주에 놓인 경우도 해결할 수 있다.
int incircle(POINT A, POINT B, POINT C, POINT D) {
    double ccw = CCW(A, B, C);
    double adx = A.x - D.x, ady = A.y - D.y,
           bdx = B.x - D.x, bdy = B.y - D.y,
           cdx = C.x - D.x, cdy = C.y - D.y,
           bdxcdy = bdx * cdy, cdxbdy = cdx * bdy,
           cdxady = cdx * ady, adxcdy = adx * cdy,
           adxbdy = adx * bdy, bdxady = bdx * ady,
           alift = adx * adx + ady * ady,
           blift = bdx * bdx + bdy * bdy,
           clift = cdx * cdx + cdy * cdy;
    double det = alift * (bdxcdy - cdxbdy)
               + blift * (cdxady - adxcdy)
               + clift * (adxbdy - bdxady);

    if (ccw > 0) return det > 0; // CCW인 경우에..
    else         return det < 0; // CW인  경우에..     
};

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

단순 다각형의 무게중심  (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
Chain Hull  (2) 2012.09.16
Posted by helloktk

댓글을 달아 주세요

728x90

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

사용자 삽입 이미지

 

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

Circle Drawing Algorithm  (1) 2008.06.03
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
Posted by helloktk

댓글을 달아 주세요

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

    계속 보게 되네요 ㅋㅋ