Least Squares Fitting of Circles

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

$$a{x}^2+b{y}^2+cxy+dx+ey+f=0$$

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

 

문제: 주어진 데이터를 fitting 하는 이차곡선(원)

$$x^2+y^2+Ax+By+C=0$$

의 계수 $A,B,C$를 최소자승법을 사용해서 구하라. 

 

주어진 점집합이 원 위의 점이면 우변이 0이 되어야 하나, 실제 데이터를 얻는 과정에서 여러 노이즈에 노출되므로 일반적으로 0이 되지 않는다. 최소자승법은 주어진 점들이 원에서 벗어나는 정도의 제곱 합이 최소가 되도록 하는 계수 $A,B,C$를 결정한다.  원과 점의 편차의 제곱합
$$L=\sum _i{|x_i^2+y_i^2+A{x}_i+B{y}_i+C|}^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$를 구할 수 있다. 계산을 단순하게 만들고 수치적 안정성을 높이기 위해 입력점들의 질량중심 

$$m_x=\frac{1}{N}\sum _i{x}_i,m_y=\frac{1}{N}\sum _i{y}_i$$

계에서 계산을 하자. 이를 위해 입력점의 $x$, $y$ 성분에서 각각 $m_x$, $m_y$만큼을 빼준 값을 좌표값으로 대입하면 된다: 

$$x_i\to x_i-m_x,y_i\to y_i-m_y$$

그러면 질량중심 좌표계에서는 $S_x=\sum _i{x}_i=0$, $S_y=\sum _i{y}_i=0$이 된다.

우선 세 번째 식에서 

$$C=-\frac{S_{x^2}+S_{y^2}}{N}$$

을 얻을 수 있고, 첫번째와 두 번째 식에서는 각각

$$S_{x^2}A+S_{xy}B=-(S_{x^3}+S_{x{y}^2})$$

$$S_{xy}A+S_{y^2}B=-(S_{y^3}+S_{x^{2}y})$$

을 얻을 수 있다. 이를 풀면

$$A=\frac{-S_{y^2}(S_{x^3}+S_{x{y}^2})+S_{xy}(S_{y^3}+S_{x^{2}y})}{S_{x^2}{S}_{y^2}-S_{xy}^2}$$

$$B=\frac{-S_{x^2}(S_{y^3}+S_{x^{2}y})+S_{xy}(S_{x^3}+S_{x{y}^2})}{S_{x^2}{S}_{y^2}-S_{xy}^2}$$

여기서 주어진 데이터의 각 차수에 해당하는 moment는 아래처럼 계산된다:

추정된 원의 중심 $(c_x,c_y)$는 

$$c_x=-\frac{A}{2},c_y=-\frac{B}{2}$$

로 주어지고, 반지름은 

$$r^2=c_x^2+c_y^2-C=c_x^2+c_y^2+\frac{1}{N}(S_{x^2}+S_{y^2})$$

로 주어진다.

Ref: I. Kasa, A curve fitting procedure and its error analysis. IEEE Trans. Inst. Meas., 25:8-14, 1976

/* 구현 코드: 2024.04.01, typing error 수정 & 질량중심계 계산으로 수정;*/
double circleFit_LS(std::vector<CPoint> &Q, double &cx, double &cy, double &radius) {
    if (Q.size() < 3) return -1;
    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;
        sx2  += xx;      sy2  += yy;      sxy  += x * y;
        sx3  += x * xx;  sy3  += y * yy;
        sx2y += xx * y;  sxy2 += yy * x;
    }
    double det = sx2 * sy2 - sxy * sxy;
    if (fabs(det) < 1.e-10) return -1;    /*collinear한 경우임;*/
    /* center in cm frame; */
    double a = sx3 + sxy2;
    double b = sy3 + sx2y;
    cx = (sy2 * a - sxy * b) / det / 2;
    cy = (sx2 * b - sxy * a) / det / 2;
    /* radius squared */
    double radsq = cx * cx + cy * cy + (sx2 + sy2) / Q.size();
    radius = sqrt(radsq);
    cx += mx; cy += my; /* recover offset; */
    return fitError(Q, cx, cy, radius);
}

https://kipl.tistory.com/357

Circle Fitting: Pratt

https://kipl.tistory.com/32

RANSAC: Circle Fit