이제 원의 방정식은 $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)$를 계산하여 중심점의 좌표를 얻을 수 있다.
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;
}
따라서, 적절한 회전 변환을 하면 두 좌표의 곱으로 주어지는 $xy$-항을 없앨 수 있다. 회전 변환을 시행하면 행렬은 eigenvalue를 성분으로 하는 대각 행렬이 된다. 회전 변환이 determinant를 보존하므로 determinant는 행렬의 두 eigenvalue의 곱으로 주어짐을 알 수 있다.
\begin{gather} \begin{pmatrix} a & b/2 \\ b/2 & c \end{pmatrix} \Longrightarrow \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} \quad \text{det} = ac - b^2/4\end{gather}
회전 후의 방정식이 타원의 방정식(원점이 이동된)을 기술하기 위해서는 $\text{det}>0$ 이어야 한다 (회전시킨 후 식에서 $x^2$과 $y^2$의 계수는 두 eigenvalue로 주어지므로 같은 부호를 가져야 한다.)
$$ F=\text{ellipse} \Leftrightarrow b^2 -4ac <0$$
conic section $F$가 타원을 기술한다면, 다음과 같이 평행이동을 시켜서 타원의 중심 $(x_0, y_0)$이 원점에 놓이게 하고, 회전변환을 시켜서 $xy$ 항을 없애도록 하자.
$$ x\to x_0 + x \cos \theta - y \sin \theta , \quad y\to y_0 + x \sin \theta + y \cos \theta$$
여기서 회전각 $\theta$는 $x$을 기준으로 측정된다. $F$에 적용해서 $xy$-항이 없어지는 조건과 1차 항이 사라지도록 하는 조건을 찾으면 타원의 중심 $(x_0, y_0)$와 $\theta$는
로 주어짐을 확인할 수 있다. Eigenvalue의 제곱근의 역수가 두 축의 반지름을 결정하므로 음수가 되어서는 안된다 (둘 중 하나가 0이면 직선이고, 둘 모두 0이면 한 점에 해당). 이는 대칭행렬의 고유값이 항상 0보다 작지 않다는 사실에서 기인한다. 위에서 구한 회전각 $\theta$를 이용해서 두 고유값을 표현하면 ,
$$\lambda_1 = a\cos^2 \theta + b \cos \theta\sin \theta + c \sin^2 \theta$$
$$\lambda_2 = a\sin^2 \theta - b \cos \theta\sin \theta + c \cos^2 \theta$$
위의 회전변환식을 대입했을 떄 나머지 상수항은 (첫 번째 등호는 계산해서 확인할 수 있음)
$$ ax_0^2 + b x_0 y_0^2 + cy_0^2 + dx_0 +e y_0 +f = f - (ax_0^2 + b x_0 y_0 + c y_0^2)\equiv -\text{scale}^{-1}$$
의 형태로 쓰여진다. 6개의 conic-parameters를 결정해야 하므로, $\mathbf{A}$는 마지막 행을 0으로 채워서 $6\times6$ 행렬로 만들어서 사용한다. 이 경우 $\mathbf{A}_{6\times 6}$는 determinant가 0인 singular 행렬이 되고 적어도 null 공간의 차원은 1 이상이다(1인 경우만 의미있는 해를 준다). $\mathbf{A}$의 singular value decomposition(SVD)을 이용하면위 방정식의 nontrivial한 해를 얻을 수 있다. SVD에 의해 $\mathbf{A}$는 다음처럼 쓸 수 있다( singular value 중 1개가 0인 경우).
행렬 $\mathbf{V}$의 각 열벡터 $\{\mathbf{v}_i\}$는 6차원 벡터 공간의 기저를 형성한다. 그리고 $\mathbf{U}$는 orthogonal 행렬이므로 오른쪽에서 곱해지는 벡터의 크기를 변화시키지 않고 단지 회전만 시킨다. singular value가 $0$에 해당하는 $\mathbf{V}$의 열벡터 $\mathbf{v}_5$를 대입하면
이어서 $\mathbf{A}.\mathbf{x}=0$의 nontrivial solution이 됨을 알 수 있다. 물론 singular value가 0인 경우가 여럿이 생기면 해당하는 열벡터의 선형결합도 해가 되므로 원하는 답이 아니다. 이 경우는 주어진 5개의 점들의 일부가 일직선상에 있는 경우에 해당한다(또는 겹치는 경우). 이때는 다시 5개의 점을 선택해야 한다.
static double fitting_error(CfPt q, double conic_params[6]) {
double x = q.x, y = q.y;
return conic_params[0] * x * x +
conic_params[1] * x * y +
conic_params[2] * y * y +
conic_params[3] * x +
conic_params[4] * y +
conic_params[5];
}
//
static int ransac_sample_number(int ndata, int ngood, double prob, int ntry);
