728x90

점집합을 일반적인 2차 곡선으로 피팅하는 경우에 방정식은

$$ a x^2 + by^2 + cxy +d x + ey +f = 0$$

의 계수를 주어진 데이터를 이용하여서 구해야 한다. 실제 문제에서는 타원, 포물선 쌍곡 선등의 타입에 따라 몇 가지 제약 조건을 넣어 피팅을 한다. 원은 타원의 특별한 경우로 일반적으로 $a = b, c = 0$의 제약 조건이 필요하다. 그러나 보다 엄밀하게 제약을 하게 되면 $a = b = 1$의 추가 조건을 줄 수 있다. 이 경우는 점들이 모두 일직선에 있는 경우를 ($a = b = 0$) 취급할 수 없게 된다. 이 예외적인 경우를 제외하고는 최소자승법을 사용하면 계수를 매우 쉽게 구할 수 있기 때문에 많이 이용된다.

 

문제: 

$$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)$$

로 주어진다.

/* 구현 코드 */
BOOL circleFit_LS(POINT Q[], int N, POINT *center, double *radius) {
    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;
    /* compute summations */
    for (int k = 0; k < N; k++) {
        double x = Q[k].x, xx = x * x;
        double y = Q[k].y, 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 */
    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)                /*collinear한 경우임;*/
        return FALSE;
    /* floating point center */
    double cx = (c1 * b2 - c2 * b1) / det;
    double cy = (a1 * c2 - a2 * c1) / det;
    /* compute radius squared */
    double radsq = cx * cx + cy * cy + (sx2 + sy2 - 2 * (sx * cx  + sy * cy)) / N;
    *radius = sqrt(radsq);
    /* integer center */
    center->x = int(cx + .5) ;
    center->y = int(cy + .5) ;
    return TRUE;
}

 

'Image Recognition > Fundamental' 카테고리의 다른 글

PCA Line Fitting  (0) 2020.11.12
Histogram Equalization  (0) 2020.11.12
Least Square Fitting of Circle  (0) 2020.11.11
Integer Sqrt  (0) 2020.11.11
Parabolic Interpolation in Peak Finding  (3) 2020.11.10
Histogram Matching  (0) 2012.11.03
Posted by helloktk

댓글을 달아 주세요