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


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



이제 원의 방정식은 xc , yc, d에 대한 1차 식으로 정리되었다. 삼각형의 외접원의 방정식은 세 꼭지점을 통과하므로 세 꼭지짐 A(x1, y1), B(x2, y2), C(x3, y3)을 위의 식에 넣으면 3개의 방정식을 얻는다. 이 방정식을 행렬로 표현하면


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


determinant를 정리하면,


을 얻을 수 있다. 이 식들의 분모는 세 꼭지점에 동일직선 위에 있지 않으면 0이 아니다. 한 꼭지점과 다른 꼭지점간의 좌표의 차이 4개 (x2-x1, y2-y1, x3-x1, y3-y1)와 좌표의 합 4개(x1+x2, y1+y2, x1+x3, y1+y3)를 계산하여서 중심점의 좌표를 얻을 수 있다. 

//

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

int circumcenter2(CPoint P, CPoint Q, CPoint R, double *xc, double *yc) {

double A = Q.x - P.x;

double B = Q.y - P.y;

double C = R.x - P.x;

double 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).  단지,식이 보다 대칭적으로 보일 뿐이다.

저작자 표시 비영리 변경 금지
신고
Posted by helloktk

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


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


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


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

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


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


따라서, 외접원의 중심은  



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


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


따라서, 

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


이므로, 이 값이 0이 아니기 위해서는 세 꼭지점이 일직선상에 놓여있지 않으면 된다.


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


또는, 성분으로 표현하면



이 설명은 공간상에 놓인 삼각형에 대해서도, 삼각형이 놓인 평면을 xy평면으로 생각하면 적용된다. 단지 조심해야 할 것은 m, n벡터를 구하는 것이다.


int circumcenter(CPoint A, CPoint B, CPoint C, double *xc, double *yc){

double ax = A.x - C.x ;

double ay = A.y - C.y ;

double bx = B.x - C.x ;

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

}



저작자 표시 비영리 변경 금지
신고
Posted by helloktk

세 개의 점 A, B, C가 만드는 삼각형의 외접원의 반지름을 구해보자. 


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


이어야 한다. 

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


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


또한, 삼각형의 면적이 세 점이 직선상에 있는가를 테스트 하는 값으로 아래 식처럼 쓰여지므로,


외접원의 반지름은 다음과 같이 구할 수 있다:


#define SQR(x) ((x)*(x))

double circumradius(CPoint A, CPoint B, CPoint C) {

double ax = C.x - B.x; 

double ay = C.y - B.y;

double bx = A.x - C.x; 

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

double 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

동일 직선상에 놓이지 않은 세점 A, B, C에 의해서 만들어지는 삼각형의 외접원 (circumcircle)에 한점 D가 포함되는가를 안되는가를 판단하는 것은 평면상의 점들을 삼각화화는 알고리즘에서 매우 중요하다. 물론 직접 외접원의 중심과 반지름을 구하고 이것을 가지고 중심과 점 D와의 거리를 계산하여서 판단을 할 수 있지만 이 경우에 수치계산의 오류가 들어갈 여지가 매우 많다(직접 중심과 반지름을 구하는 식을 만들어 보면 쉽게 알 수 있다). 여기서는 이것을 회피하는 방법을 고안하여서 보다 robust한 판단을 주는 식을 유도하려고 한다.

이것은 공간상의 네점 A(ax, ay, az), B(bx, by, bz), C(cx, cy, cz), D(dx, dy, dz)로 만들어지는 사변형의 체적을 구하는 공식이 아래의 행렬식으로 주어진다는 사실에서 출발한다.

사용자 삽입 이미지

이 사실을 이용하면 평면상의 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;

으로 주어지고, 이 평면과 포물체가 만나는 경계선이 바닥에 정사영하면 아래의 방정식을 얻는다.

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

즉, 단면의 프로젝션은 중심이 (a,b)이고 반지름이 r^2인 원이 된다.

따라서, 원주상의 점들은 포물체상의 곡선상에 존재하고, 원의 내부점은 평면아래의 포물체면에 존재하게되고, 원주 밖의 점은 평면위의 포물체면에 존재하게 된다. 세점이 그리는 원인 경우에 원 내부와 외부점에 의한 사면체의 꼭지점이 평면을 기준으로 반대로 향하게 된다. 따라서 이 사면체의 체적도 부호를 달리하게 된다. 즉, 체적의 부호를 판별함으로써, 외접원내에 있는지 밖에 있는지를 판별할 수 있는 것이다.



다시 정리하면, 세점 a(ax, ay), b(bx, by) c(cx, cy) d(dx, dy)가 주어진 경우에 a, b, c(반시계방향)에 의해서 형성이 되는 원에 점 d가 포함이 되는지의 여부는 공간상의 네점

        A = (ax, ay, ax^2+ay^2),
        B = (bx, by, bx^2+by^2),
        C = (cx, cy, cx^2+cy^2),
        D = (dx, dy, dx^2+dy^2)

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

사용자 삽입 이미지


위의 행렬식은 좀 더 단순화 시킬 수 있다. 1, 2, 3행에서 4행을 행단위로 빼도 값이 변화지 않으므로 오른쪽 식을 얻고(4열에 대해서 전개 후), 1열에 2*dx를 곱한것과 2열에 2*dy를 곱한것을 각각 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' 카테고리의 다른 글

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