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인  경우에..     
};

삼각형의 외접원을 기하학적인 방법이 아닌 대수적인 방법을 이용해서 구해보자. 중심이 $(x_c, y_c)$이고 반지름이 $R$인 원의 방정식은

$$(x-x_c)^2 + (y- y_c)^2 = R^2$$

으로 주어진다. 이 식에서 세 개의 변수 $x_c$, $y_c$, $R$를 결정해야는데, 각 변수에 대해서 2차 방정식의 형태로 주어져 있으므로 좀 더 쉬운 1차식 형태가 되도록 변형을 해보자. 방정식을 전개하면

$$(2x) x_c + (2y) y_c + R^2 - x_c^2 - y_c^2 = x^2 + y^2  \\ \longrightarrow ~(2x) x_c + (2y) y_c +d = x^2 +y^2,\quad\quad d=R^2 -x_c^2 -y_c^2.$$

이제 원의 방정식은 $x_c$, $y_c$, $d$에 대한 1차식으로 정리되었다. 삼각형의 외접원의 방정식은 세 꼭짓점을 통과하므로 세 꼭지 짐 $A(x_1, y_1)$, $B(x_2, y_2)$, $C(x_3, y_3)$을 위의 식에 넣으면 3개의 방정식을 얻는다. 이 방정식을 행렬로 표현하면

따라서, Cramer의 규칙을 적용하면 다음과 같이 외접원의 중심 (또한 $d$을 구해서 외접원의 반지름을 구할 수 있다)을 얻을 수 있다.

determinant를 정리하면,

을 얻을 수 있다. 이 식의 분모는 세 꼭짓점이 일직선 위에 있지 않으면 0이 아니다. 한 꼭짓점과 다른 꼭짓점 간 좌표의 차이 4개 $(x_2-x_1, y_2-y_1, x_3-x_1, y_3-y_1)$와 좌표의 합 4개 $(x_1+x_2, y_1+y_2, x_1+x_3, y_1+y_3)$를 계산하여 중심점의 좌표를 얻을 수 있다.

구현된 코드는 http://kipl.tistory.com/113; 또는

int circumcenter2(CPoint P, CPoint Q, CPoint R, double *xc, double *yc) {
    double A = Q.x - P.x, B = Q.y - P.y;
    double C = R.x - P.x, D = R.y - P.y;
    double E = A * (P.x + Q.x) + B * (P.y + Q.y);
    double F = C * (P.x + R.x) + D * (P.y + R.y);
    double G = 2. * ( A * D - B * C);
    if (G != 0.) { // 세 점이 일직선에 놓이지 않는 경우; 이 경우만 외접원이 정의된다;
        *xc = (D * E - B * F) / G;
        *yc = (A * F - C * E) / G;
        return 1;
    } else return 0;
}

분모 $G$의 계산을 $G = 2 ( A  (R.y - Q.y) - B  (R.x - Q.x))$ 로 쓰는 경우가 많다. 위의 행렬식을 보면 같은 결과지만 식을 좀 더 대칭적으로 보이게 한다: $G = 2.(A  (D - B) - B  (C - A)) = 2 (A D - B C)$.

평면 위 세 점 $A$, $B$, $C$가 삼각형을 만들 때 외접원 중심을 구해보자. 점 $C$을 원점으로 하면, $A$, $B$는 각각

의 두 벡터로 표현이 된다 (원점을 $C$점으로 평행 이동했다고 생각하면 된다). 이 두 벡터에 각각 수직이고 삼각형과 같은 평면에 놓인 ($\vec a$, $\vec b$를 각각 시계 방향으로 회전시켜 만든) 두 벡터 $\vec m$, $\vec n$를 다음과 같이 정의하자.

그러면, 변 $CA$의 중점을 지나면서 수직인 직선의 방정식은

또, $CB$의 중점을 지나면서 수직인 직선의 방정식은 

이 두 직선의 교점이 외접원의 중심이 된다. 매개변수 $s$, $t$를 구하기 위해서 두 식을 빼면

 매개변수 $s$는 $\vec b$와 $\vec n$이 수직인 조건을 이용하면

따라서, 외접원의 중심은  

여기서, 3개 벡터의 외적의 성질을 이용해서 (평면에서 외적은 숫자이다)

가 항등식으로 성립하고, 또한 $\vec a$, $\vec m$이 서로 수직이면서 길이가 같고, $\vec m$ 도 $\vec b$에 수직이면 길이가 같으므로

따라서, 

여기서, 윗 식의 분모를 보면

이므로, 이 값이 0이 아니려면 세 점이 일직선에 있지 않으면 된다.

그런데 이 벡터는 점 $C$를 원점으로 하여서 계산을 한 것이므로, 원래의 좌표계에 대한 식으로 바꾸려면 $(C_x, C_y)$를 더해 주어야 한다.

또는, 성분으로 표현하면

int circumcenter ( CPoint A, CPoint B, CPoint C, double *xc, double *yc) {
    double ax = A.x - C.x, ay = A.y - C.y ;
    double bx = B.x - C.x, by = B.y - C.y ;
    double asq = ax * ax + ay * ay;
    double bsq = bx * bx + by * by;
    double ccw = ax * by - ay * bx;
    if ( ccw != 0. ) { //세 점임 일직선 위에 있지 않는 경우; 이 경우만 외접원이 정의됨;
        *xc = C.x + ( by * asq - ay * bsq ) / ( 2 * ccw ) ;
        *yc = C.y + ( -bx * asq + ax * bsq ) / ( 2 * ccw ) ;
        return 1;
    } else return 0;
}

세 점 $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}| }.$$

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 = hypot(ax, ay);
        double b = hypot(bx, by); 
        double cx = B.x - A.x, cy = B.y - A.y;       
        double c = hypot(cx, cy);
        return (0.5 * a * b * c / fabs(crossab));
    } else 
        return 0;   //collinear
};