Ellipse Fitting

일반적인 conic section 피팅은 주어진 데이터 $\{ (x_i, y_i)\}$를 가장 잘 기술하는 이차식

$$F(x, y) = ax^2 + bxy +cy^2 + dx +ey +f=0 $$

의 계수 ${\bf u^T}= (a,b,c,d,e,f)$을 찾는 문제이다. 이 conic section이 타원이기 위해서는 2차항의 계수 사이에 다음과 같은 조건을 만족해야 한다.

$$\text{ellipse constraint:}~~ ac - b^2/4 >0$$

그리고 얼마나 잘 피팅되었난가에 척도가 필요한데 여기서는 주어진 데이터의 대수적 거리 $F(x,y)$을 이용하자. 주어진 점이 타원 위의 점이면 이 값은 정확히 0이 된다. 물론 주어진 점에서 타원까지의 거리를 사용할 수도 있으나 이는 훨씬 복잡한 문제가 된다.  따라서 해결해야 하는 문제는

\begin{gather}L = \sum _{i}  \left( ax_i^2 + bx_i y_i + cy_i^2 +dx_i + e y_i +f\right)^2 - \lambda( 4ac-b^2-1) \\= \left|\begin{pmatrix}x_0^2& x_0y_0 & y_0^2 & x_0 & y_0 & 1\\ x_1^2 & x_1 y_1& y_1^2 & x_1 & y_1 & 1 \\ x_2^2 & x_2y_2& y_2^2 & x_2& y_2 & 1\\ &&\vdots \\\end{pmatrix}\begin{pmatrix}a\\b\\c\\d\\e\\f \end{pmatrix}  \right|^2 -\lambda \left({\bf  u^T} \begin{pmatrix} 0& 0& 2&0&0&0\\ 0 &-1&0 &0 &0 &0\\ 2&0&0&0&0&0\\0&0&0&0&0&0 \\0&0&0&0&0&0\\0&0&0&0&0&0&  \end{pmatrix} \bf u -1\right) \\ =\bf u^T D^TD u -\lambda (u^T C u -1)\\ = \bf u^T S u -\lambda (u^T C u-1)\end{gather}

을 최소화시키는 계수 벡터 $\bf u$를 찾는 것이다. 여기서 제한조건으로 $4ac - b^2 =1= \bf u^T C u$로 설정했다. 

$\bf u^T$에 대해서 미분을 하면 

$$ \frac{\partial L}{\partial \bf u^T} =  \bf S u -\lambda C u=0$$

즉, 주어진 제한조건 $4ac - b^2=1$하에서 대수적 거리를 최소화시키는 타원방정식의 계수 $\bf u$를 구하는 문제는 scattering matrix $\bf S=D^T D$에 대한 일반화된 고유값 문제로 환원이 된다.

$$  \bf S u =\lambda C u \\ u^T C u =1$$

이 문제의 풀이는 직전의 포스팅에서 다른 바 있는데 $\bf S$의 제곱근 행렬 $\bf Q=S^{1/2}$를 이용하면 된다. 주어진 고유값 $\lambda$와 고유벡터 $\bf u$가 구해지면 대수적 거리는 $$\bf u^T S u = \lambda$$

이므로 이를 최소화시키기 위해서는 양의 값을 갖는 고유값 중에 최소에 해당하는 고유벡터를 고르면 된다. 그런데 고유값 $\lambda$의 부호별 개수는 $\bf C$의 고유값 부호별 개수와 동일함을 보일 수 있는데 (Sylverster's law of inertia),  $\bf C$의 고유값이 $\{-2,-1,2,0,0,0\}$이므로 $\lambda>0$인 고유값은 1개 뿐임을 알 수 있다. 따라서 $\bf S u = \lambda C u$를 풀어서 얻은  유일한 양의 고유값에 해당하는 고유벡터가 원하는 답이 된다.

https://kipl.tistory.com/370

Least Squares Fitting of Ellipses

https://kipl.tistory.com/565

Generalized eigenvalues problem

 Ref: https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/ellipse-pami.pdf

 

 

double FitEllipse(std::vector<CPoint>& points, double einfo[6] ) {     
    if ( points.size() < 6 ) return -1;
    double eigvals[6];
    std::vector<double> D(6 * points.size());
    double S[36];/*  S = ~D * D  */
    double C[36];
    double EIGV[36];/* R^T; transposed orthogonal matrix;*/

    double offx = 0, offy = 0;
    /* shift all points to zero */
    for(int i = points.size(); i--> 0; ) {	
        offx += points[i].x;
        offy += points[i].y;        	
    }
    offx /= points.size(); 
    offy /= points.size();

    /* for the sake of numerical stability, scale down to [-1:1];*/
    double scale = 0; 
    for (int i = points.size(); i-->0; ) {
        if (points[i].x > scale) scale = points[i].x;
        if (points[i].y > scale) scale = points[i].y;
    }
    double invscale = 1 / scale;
    /* ax^2 + bxy + cy^2 + dx + ey + f = 0*/
    /* fill D matrix rows as (x*x, x*y, y*y, x, y, 1 ) */
    for(int i = points.size(); i--> 0; ) {	
        double x = points[i].x - offx; x *= invscale; 
        double y = points[i].y - offy; y *= invscale;
        D[i*6 + 0] = x*x; D[i*6 + 1] = x*y;
        D[i*6 + 2] = y*y; D[i*6 + 3] = x;
        D[i*6 + 4] = y;   D[i*6 + 5] = 1;		
    }			

    /* scattering matrix: S = ~D * D (6x6)*/
    for (int i = 0; i < 6; i++) 
        for (int j = i; j < 6; j++) { /*upper triangle;*/
            double s = 0;
            for (int k = points.size(); k-- > 0; ) 
                s += D[k*6 + i] * D[k*6 + j];
            S[i*6 + j] = s;
        }
    for (int i = 1; i < 6; i++) /*lower triangle;*/
        for (int j = 0; j < i; j++) 	
            S[i*6 + j] = S[j*6 + i] ;
    
    /* fill constraint matrix C */
    for (int i = 0; i < 36 ; i++ ) C[i] = 0;
    C[12] =  2 ;//2x0 
    C[2 ] =  2 ;//0x2 
    C[7 ] = -1 ;//1x1

    /* find eigenvalues/vectors of scattering matrix; */
    double RT[36];	/* each row contains eigenvector; */
    JacobiEigens ( S, RT, eigvals, 6, 0 );
    /* create R and INVQ;*/
    double R[36];
    for (int i = 0; i < 6 ; i++) {
        eigvals[i] = sqrt(eigvals[i]);
        for ( int k = 0; k < 6; k++ ) {
            R[k*6 + i] = RT[i*6 + k];  /* R = orthogonal mat = transpose(RT);*/
            RT[i*6 + k] /= eigvals[i]; /* RT /= sqrt(eigenvalue) row-wise)*/
        }
    }
    /* create INVQ=R*(1/sqrt(eigenval))*RT;*/
    double INVQ[36];
    _MatrixMul(R, RT, 6, INVQ);

    /* create matrix INVQ*C*INVQ */
    double TMP1[36], TMP2[36];
    _MatrixMul(INVQ, C, 6, TMP1 );
    _MatrixMul(TMP1, INVQ, 6, TMP2 );
    
    /* find eigenvalues and vectors of INVQ*C*INVQ:*/
    JacobiEigens ( TMP2, EIGV, eigvals, 6, 0 );
    /* eigvals stores eigenvalues in descending order of abs(eigvals);*/
    /* search for a unique positive eigenvalue;*/
    int index = -1, count = 0;
    for (int i = 0 ; i < 3; i++ ) {
        if (eigvals[i] > 0) {
            index = i; // break;
            count++;
        }
    }
    /* only 3 eigenvalues must be non-zero 
    ** and only one of them must be positive;*/
    if ((count != 1) || (index == -1)) 
        return -1;
     
    /* eigenvector what we want: u = INVQ * v */
    double u[6]; 
    double *vec = &EIGV[index*6];
    for (int i = 0; i < 6 ; i++) {
        double s = 0;
        for (int k = 0; k < 6; k++) s += INVQ[i*6 + k] * vec[k];
        u[i] = s;
    }
    /* extract shape infos;*/
    PoseEllipse(u, einfo);
    /* recover original scale; center(0,1) and radii(2,3)*/
    for (int i = 0; i < 4; i++) einfo[i] *= scale;
    /* recover center */
    einfo[0] += offx; 
    einfo[1] += offy;
    return FitError(points, offx, offy, scale, u);
};