점집합을 일반적인 2차 곡선으로 피팅하는 경우에 방정식은
$$ a x^2 + by^2 + cxy +d x + ey +f = 0$$
의 계수를 주어진 데이터를 이용하여서 구해야 한다. 실제 문제에서는 타원, 포물선 쌍곡 선등의 타입에 따라 몇 가지 제약 조건을 넣어 피팅을 한다. 원은 타원의 특별한 경우로 일반적으로 $a = b$, $c = 0$의 제약 조건이 필요하다. 그러나 보다 엄밀하게 제약을 하게 되면 $a = b = 1$의 추가 조건을 줄 수 있다. 이 경우는 점들이 모두 일직선에 있는 경우를 ($a = b = 0$) 취급할 수 없게 된다. 이 예외적인 경우를 제외하고는 최소자승법을 사용하면 계수를 매우 쉽게 구할 수 있기 때문에 많이 이용된다.
문제: 주어진 데이터를 fitting 하는 이차곡선(원)
$$x^2 + y^2 + A x + B y + C = 0$$
의 계수 $A, B, C$를 최소자승법을 사용해서 구하라.
주어진 점집합이 원 위의 점이면 우변이 0이 되어야 하나, 실제 데이터를 얻는 과정에서 여러 노이즈에 노출되므로 일반적으로 0이 되지 않는다. 최소자승법은 주어진 점들이 원에서 벗어나는 정도의 제곱 합이 최소가 되도록 하는 계수 $A, B, C$를 결정한다. 원과 점의 편차의 제곱합
$$ L=\sum_ i \left |x_i^2 + y_i^2 + A x_i + B y_i + C \right|^2 , $$
의 극값을 찾기 위해서 $A, B,$ 그리고 $C$에 대해 미분을 하면
$$\frac{\partial L}{\partial A} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) x_i = 0, $$
$$\frac{\partial L}{\partial B} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) y_i = 0, $$
$$\frac{\partial L}{\partial C} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) = 0. $$
이 연립방정식을 풀면 $A, B, C$를 구할 수 있다. 우선 세 번째 식에서
$$ CN = -S_{x^2} - S_{y^2} - AS_x - BS_y ,$$
을 얻고, 이를 첫번째와 두 번째 식에 각각 대입하면
$$A ( NS_{x^2} - S_x^2) + B ( N S_{xy} - S_x S_y ) =-N S_{x^3} - N S_{xy^2} + S_{x^2} S_x + S_{y^2} S_x, $$
$$A ( NS_{xy} - S_x S_y ) + B ( N S_{y^2} - S_y^2) = -N S_{y^3} - N S_{x^2 y} +S_{x^2} S_y +S_{y^2} S_y, $$
을 얻을 수 있다. 다시 정리하면 두 개의 연립방정식
$$-a_1 A - a_2 B = 2c_1,$$
$$-b_1 A - b_2 B = 2c_2,$$
을 얻는다. $a_1, a_2 = b_1, b_2, c_1, c_2$는 코드에서 정의되어 있다. 그리고
따라서, 추정된 원의 중심 $(c_x, c_y)$는
$$ c_x = - \frac{A}{2} = \frac{c_1 b_2 - c_2 a_2}{ a_1 b_2 - a_2 b_1},$$
$$ c_y = - \frac{B}{2} = \frac{-c_1 b_1 + c_2 a_1}{a_1 b_2 -a_2 b_1},$$
로 주어지고, 반지름은
$$r^2 = c_x^2 +c_y^2 - C = c_x^2 + c_y^2 + \frac{1}{N}( S_{x^2}+S_{y^2}- 2c_x S_x - 2 c_y S_y)$$
로 주어진다.
Note: 질량중심 좌표계로 옮긴 후 moment를 계산하면 $S_x = S_y =0$이어서 식이 더 간결해지고, 수치적인 안정성도 좋아질 수 있다.
Ref: I. Kasa, A curve fitting procedure and its error analysis. IEEE Trans. Inst. Meas., 25:8-14, 1976
/* 구현 코드:*/
double circleFit_LS(std::vector<CPoint>& Q, double& cx, double& cy, double& radius) {
if (Q.size() < 3) return -1;
double sx = 0.0, sy = 0.0;
double sx2 = 0.0, sy2 = 0.0, sxy = 0.0;
double sx3 = 0.0, sy3 = 0.0, sx2y = 0.0, sxy2 = 0.0;
double mx = 0, my = 0; /* center of mass;*/
for (int k = Q.size(); k-->0;)
mx += Q[k].x, my += Q[k].y;
mx /= Q.size(); my /= Q.size();
/* compute moments; */
for (int k = Q.size(); k-->0;) { /* offset (mx, my)*/
double x = Q[k].x - mx, xx = x * x;
double y = Q[k].y - my, yy = y * y;
sx += x; sy += y;
sx2 += xx; sy2 += yy; sxy += x * y;
sx3 += x * xx; sy3 += y * yy;
sx2y += xx * y; sxy2 += yy * x;
}
/* compute a's,b's,c's; */
const int N = Q.size();
double a1 = 2.0 * (sx * sx - sx2 * N);
double a2 = 2.0 * (sx * sy - sxy * N);
double b1 = a2;
double b2 = 2.0 * (sy * sy - sy2 * N);
double c1 = (sx2 + sy2) * sx - (sx3 * N + sxy2) * N;
double c2 = (sx2 + sy2) * sy - (sy3 * N + sx2y) * N;
double det = a1 * b2 - a2 * b1;
if (fabs(det) < 1.e-10) return -1; /*collinear한 경우임;*/
/* center; */
cx = (c1 * b2 - c2 * b1) / det;
cy = (a1 * c2 - a2 * c1) / det;
/* radius squared */
double radsq = cx * cx + cy * cy + (sx2 + sy2 - 2 * (sx * cx + sy * cy)) / N;
radius = sqrt(radsq);
cx += mx; cy += my; /* recover offset; */
return FitError(Q, cx, cy, radius);
}
'Image Recognition > Fundamental' 카테고리의 다른 글
PCA Line Fitting (0) | 2020.11.12 |
---|---|
Histogram Equalization (0) | 2020.11.12 |
Integer Sqrt (0) | 2020.11.11 |
Parabolic Interpolation in Peak Finding (3) | 2020.11.10 |
Histogram Matching (0) | 2012.11.03 |