Geometry & Recognition

일직선에 있지 않은 세 점 $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 행렬식은 쉽게 계산이 된다. $\text{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인  경우에..     
};