Circle Fitting: Pratt

주어진 점집합을 원으로 피팅하기 위해 이차식

$$A(x^2 + y^2) + Bx + Cy +D =0$$

을 이용하자. 원의 경우는 $A=0$인 경우는 직선을 나타내고, $A\ne0$인 경우가 원을 표현한다. 물론 $A=1$로 설정을 할 수 있으나 이 경우는 주어진 데이터가 원의 일부 부분만 나타내는 경우 문제가 생긴다. 원은 중심과 반지름을 알면 결정되므로 3개의 변수가 필요한데 위의 이차 다항식은 4개의 변수를 가진다. 원을 제대로 표현하기 위해서는 제한 조건이 들어오는데 

$$ \Big(x + \frac {B}{2A} \Big)^2 + \Big(y + \frac {C}{2A} \Big)^2 = \frac {B^2 +C^2 -4A  D}{4A^2}$$

으로 쓸 수 있으므로 좌변이 양수이므로 

$$ B^2+C^2 - 4A  D>0$$이어야 한다. 따라서 $$ B^2 + C^2 - 4A  D=1$$인 제한 조건을 둘 수 있다.

주어진 각 점들이 정확히 추정하는 원 상에 있지 않으므로 이차식의 좌변의 값이 반드시 0이 되지는 않는다. 잘 피팅하기 위해서는 제한조건을 만족하면 각 점들에서 차이를 최소로 만드는 파라미터를 찾으면 된다:

$$ \epsilon(A,B, C, D)= \sum | A(x_i ^2 + y_i^2) + Bx_i + Cy_i +D |^2 - \lambda (B^2 +C^2 - 4A D-1) \longrightarrow \text {min}$$

여기서 $\lambda$는 Lagrange multiplier이다. 식을 행렬로 표현하기 위해

$$z_i= x_i^2 + y_i^2 , ~M_{zz}= \frac {1}{n} \sum z_i^2 , ~M_{xx} = \frac{1}{n} \sum x_i^2, ~M_{yy}=\frac{1}{n}\sum y_i^2,  $$

$$M_{zx}= \frac {1}{n} \sum z_i x_i  , ~M_{zy} = \frac{1}{n} \sum z_i y_i , ~M_{xy}=\frac{1}{n}\sum x_i y_i, $$

$$M_{z}= M_{xx}+M_{yy} , ~M_{x} = \frac{1}{n} \sum x_i, ~M_{y}=\frac {1}{n}\sum y_i,  $$

$${\tt x}= \left(\begin {array}{c} A \\ B \\ C \\D\end {array}\right) , ~ \tt M = \left( \begin{array}{cccc} M_{zz}& M_{zx}& M_{zy} & M_z \\ M_{zx} & M_{xx}& M_{xy} & M_{x}\\ M_{zy} & M_{xy} & M_{yy} &M_y \\M_z &M_x &M_y& 1 \end{array} \right), ~\tt  N=\left( \begin{array}{cccc} 0& 0& 0& -2 \\ 0& 1& 0 & 0\\ 0& 0& 1& 0\\ -2&0&0& 0 \end{array} \right)$$

로 놓으면 ($\tt M=\text {scattering matrix}$) $$ n^{-1} \epsilon= {\tt x}^T. {\tt M}. {\tt x} - \lambda ( {\tt x}^T. {\tt N}. {\tt x} - 1)$$

로 표현된다. 질량중심 좌표계에서 계산을 하면 $M_x= M_y =0$이므로 행렬이 더 간단해진다.

$\epsilon$을 $\tt x^T$에 대해서 미분하면

$$ \tt M.x = \lambda \tt N.x$$

인 lagrange multiplier가 고윳값이 되는 일반화된 고윳값 방정식을 얻는다. 이 고윳값 식이 해를 가지려면

$$\text {det}( \tt M -\lambda \tt N)=0$$

이어야 한다. 이 특성 방정식은 4차 다항식이므로 closed form 형태의 해를 쓸 수 있으나 실질적으로 수치해석적으로 구하는 것이 더 편하다(Newton-Raphson). $\text {min}(\epsilon)$은 제한조건 때문에

$$ n^{-1}\epsilon = \lambda \tt x^T. N. x = \lambda > 0$$

이므로 $0$보다 큰 최소 고윳값이 원하는 해이다. ($\lambda = 0$은 $\text {det}(\tt M)=0$, 즉, $\tt M$이 singular 한 경우)

그런데 주어진 고윳값 $\lambda$ 해당하는 고유 벡터가 $\tt x$일 때, 임의의 0이 아닌 상수 $\mu$에 대해 $\mu \tt x$도 역시 고유벡터이다. 고유 벡터의 크기를 고정하기 위해서 제한조건을 적용하면 $1=\mu^2 \tt x^T .N .x = \mu^2 \tt x^T. M. x/\lambda$이므로 $\mu = \sqrt{\lambda/ \tt x^T. M. x}$로 선택하면 된다. 그렇지만 해가 원을 나타내는 경우 원의 중심 $(-B/2A, -C/2A)$, 반지름 $\sqrt{ (B/2A)^2 + (C/2A)^2 - D/A}$은 계수의 비에만 의존하므로 굳이 normalized 된 고유벡터를 구할 필요가 없다. 따라서, 해 중에 직선인 경우는 무시하고($A=0$) 원만을 선택할 때는 $A=1$로 놓고 계산을 해도 영향을 받지 않는다. 또한 이 경우 고유벡터 성분을 구체적으로 적을 수 있어 수치해석에 의존할 필요도 없어진다.

최소의 고윳값에 대해서 고유 방정식을 풀면 고유벡터의 성분은

$$ C_{xy} \equiv  M_{xx} M_{yy}- M_{xy}^2 , ~~\Delta \equiv \lambda^2 - M_z   \lambda +C_{xy}$$

$$A=1$$

$$B = \frac {M_{zx}(M_{yy} - \lambda)-M_{zy} M_{xy}}{ \Delta }$$

$$C = \frac {M_{zy}(M_{xx} - \lambda)-M_{zx} M_{xy}}{ \Delta }$$

$$D=-M_z -2\lambda$$로 주어진다. 따라서 원의 중심은

$$ c_x = -\frac {B}{2A}= \frac{    M_{zx}(M_{yy} -\lambda) - M_{zy}M_{xy}}{2  \Delta}, \quad c_y= -\frac{C}{2A}= \frac{ M_{zy}(M_{xx} -\lambda )- M_{zx} M_{xy}}{2  \Delta} $$

원의 반지름은 

$$ r =\sqrt { (B/2A)^2+(C/2A)^2 -D/A} = \sqrt{ c_x^2 + c_y^2 +M_{z} +2\lambda}$$ 

로 주어진다.

 

** 추가: Sylvester's law of inertia에 의하면 위에서 구하는 고유값 $\lambda$의 부호 개수는  $\tt N$의 고유값 부호별 개수와 같아야 하는데, $\tt N$의 고유값이 $\{-2,1,1,2\}$이므로 양수인 고유값 $\lambda$는 3개가 존재한다.

 

Ref: V. Pratt, Direct least-squares fitting of algebraic surfaces, Computer Graphics 21, 1987, 145–152.

double circleFit(std::vector<CPoint>& data, double& centerx, double& centery, double& radius) {
    const int maxIter = 99;
    double mx = 0, my = 0;
    for (int i = data.size(); i-- > 0;)
        mx += data[i].x, my += data[i].y;

    // center of mass;
    mx /= data.size(); my /= data.size();

    // moment calculation;
    double Mxy, Mxx, Myy, Mzx, Mzy, Mzz;
    Mxx = Myy = Mxy = Mzx = Mzy = Mzz = 0.;
    for (int i = data.size(); i-- > 0;) {
        double xi = data[i].x - mx;   //  center of mass coordinate
        double yi = data[i].y - my;   //  center of mass coordinate
        double zi = xi * xi + yi * yi;
        Mxy += xi * yi; Mxx += xi * xi;
        Myy += yi * yi; Mzx += xi * zi;
        Mzy += yi * zi; Mzz += zi * zi;
    }
    Mxx /= data.size(); Myy /= data.size();
    Mxy /= data.size(); Mzx /= data.size();
    Mzy /= data.size(); Mzz /= data.size();
    double Mz = Mxx + Myy;
    double Cxy = Mxx * Myy - Mxy * Mxy;
    double Vz = Mzz - Mz * Mz;
    // coefficients of characteristic polynomial;
    double C2 = 4 * Cxy - 3 * Mz * Mz - Mzz; 
    double C1 = Vz * Mz + 4 * Cxy * Mz - Mzx * Mzx - Mzy * Mzy;
    double C0 = Mzx * (Mzx * Myy - Mzy * Mxy) + Mzy * (Mzy * Mxx - Mzx * Mxy) - Vz * Cxy;
    // Newton's method starting at lambda = 0
    double lambda = 0;
    double y = C0;  // det(lambda = 0)
    for (int iter = 0; iter < maxIter; iter++) {
        double Dy = C1 + lambda * (2. * C2 + 16.*lambda * lambda);
        double lambdaNew = lambda - y / Dy;
        if ((lambdaNew == lambda)) break;
        double ynew = C0 + lambdaNew * (C1 + lambdaNew * (C2 + 4 * lambdaNew * lambdaNew));
        if (fabs(ynew) >= fabs(y))  break;
        lambda = lambdaNew;
        y = ynew;
    }
    double DEL = lambda * lambda - lambda * Mz + Cxy;
    double cx = (Mzx * (Myy - lambda) - Mzy * Mxy) / DEL / 2;
    double cy = (Mzy * (Mxx - lambda) - Mzx * Mxy) / DEL / 2;
    centerx = cx + mx;   // recover origianl coordinate;
    centery = cy + my;
    radius = sqrt(cx * cx + cy * cy + Mz + 2 * lambda);
    return fitError(data, centerx, centery, radius);
}

https://kipl.tistory.com/32

RANSAC: Circle Fit

https://kipl.tistory.com/207

Least Squares Fitting of Circles