Geometry & Recognition

주어진 데이터를 잘 피팅하는 직선을 찾기 위해서는 데이터를 이루는 각 점의 $y$ 값과 같은 위치에서 구하려는 직선과의 차이 (residual)의 제곱을 최소화시키는 조건을 사용한다. 그러나 직선의 기울기가 매우 커지는 경우에는  데이터에 들어있는 outlier에 대해서는 그 차이가 매우 커지는 경우가 발생할 수도 있어 올바른 피팅을 방해할 수 있다. 이를 수정하는 간단한 방법 중 하나는 $y$값의 차이를 이용하지 않고 데이터와 직선 간의 최단거리를 이용하는 것이다. 

수평에서 기울어진 각도가 $\theta$이고 원점에서 거리가 $s$인 직선의 방정식은 

$$x\sin \theta -y\cos \theta +s=0$$

이고, 한 점 $(x_i,y_i)$에서 이 직선까지 거리는

$$d_i=|\sin \theta y_i-\cos \theta x_i+s|$$

이다. 따라서 주어진 데이터에서 떨어진 거리 제곱의 합이 최소인 직선을 구하기 위해서는 다음을 최소화시키는 $\theta ,s$을 구해야 한다. $$L=\sum _i(\sin \theta x_i-\cos \theta y_i+s{)}^2$$$$(\theta ,s)=\text{argmin}(L)$$

$\theta$와 $s$에 대한 극값 조건에서 

$$\frac{1}{2}\frac{\partial L}{\partial \theta }=\frac{1}{2}\sin 2\theta \sum _i(x_i^2-y_i^2)-\cos 2\theta \sum x_i{y}_i+s\sin \theta \sum _i{x}_i+s\cos \theta \sum _i{y}_i=0$$

$$\frac{1}{2}\frac{\partial L}{\partial s}=\cos \theta \sum y_i-\sin \theta \sum _i{x}_i-Ns=0$$

주어진 데이터의 질량중심계에서 계산을 수행하면 $\sum _i{x}_i=\sum _i{y}_i=0$ 이므로 데이터의 2차 모멘트를 $$A=\sum _i(x_i^2-y_i^2),B=\sum _i{x}_i{y}_i$$로 놓으면 직선의 파라미터를 결정하는 식은

$$\begin{gathered} \frac{1}{2}A\sin 2\theta -B\cos 2\theta =0\to \tan 2\theta =\frac{2B}{A} \\ s=0 \end{gathered}$$

두 번째 식은 직선이 질량중심(질량중심계에서 원점)을 통과함을 의미한다. 첫번째 식을 풀면

$$\tan \theta =\frac{-A\pm \sqrt{A^2+(2B{)}^2}}{2B}$$

두 해 중에서 극소값 조건을 만족시키는 해가 직선을 결정한다. 그런데

$$\frac{1}{2}\frac{{\partial }^{2}L}{\partial {\theta }^2}=A\cos 2\theta +2B\sin 2\theta =\pm \sqrt{A^2+(2B{)}^2}>0$$

이므로 위쪽 부호로 직선($x\sin \theta =y\cos \theta$)이 정해진다. 질량중심계에서는 원점을 지나지만 원좌표계로 돌아오면 데이터의 질량중심을 통과하도록 평행이동시키면 된다.

$$(-A+\sqrt{A^2+(2B{)}^2})(x-\overline{x})=2B(y-\overline{y})$$

여기서 주어진 데이터의 질량중심은 원좌표계에서

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

이다. 또한 원좌표계에서 $A$와 $B$의 계산은 

$$A=\sum _i[(x_i-\overline{x}{)}^2-(y_i-\overline{y}{)}^2],B=\sum (x_i-\overline{x})(y_i-\overline{y})$$

이 결과는 데이터 분포에 PCA를 적용해서 얻은 결과와 동일하다. PCA에서는 공분산 행렬의 고유값이 큰 쪽에 해당하는 고유벡터가 직선의 방향을 결정했다. https://kipl.tistory.com/211.  또한 통계에서는 Deming regression이라고 불린다.

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

$$F(x,y)=a{x}^2+bxy+c{y}^2+dx+ey+f=0$$

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

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

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

$$\begin{array}{c}L=\sum _i{(a{x}_i^2+b{x}_i{y}_i+c{y}_i^2+d{x}_i+e{y}_i+f)}^2-\lambda (4ac-b^2-1)\\ ={\left|\left(\begin{array}{cccccc}{x}_0^2& x_0{y}_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_2{y}_2& y_2^2& x_2& y_2& 1\\ & & ⋮\end{array}\right)\left(\begin{array}{c}a\\ b\\ c\\ d\\ e\\ f\end{array}\right)\right|}^2-\lambda ({\mathbf{u}}^{\mathbf{T}}\left(\begin{array}{cccccc}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{array}\right)\mathbf{u}\mathbf{-}\mathbf{1})\\ ={\mathbf{u}}^{\mathbf{T}}{\mathbf{D}}^{\mathbf{T}}\mathbf{D}\mathbf{u}\mathbf{-}\lambda \mathbf{(}{\mathbf{u}}^{\mathbf{T}}\mathbf{C}\mathbf{u}\mathbf{-}\mathbf{1}\mathbf{)}\\ ={\mathbf{u}}^{\mathbf{T}}\mathbf{S}\mathbf{u}\mathbf{-}\lambda \mathbf{(}{\mathbf{u}}^{\mathbf{T}}\mathbf{C}\mathbf{u}\mathbf{-}\mathbf{1}\mathbf{)}\end{array}$$

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

${\mathbf{u}}^{\mathbf{T}}$에 대해서 미분을 하면 

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

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

$$\begin{gathered} \mathbf{S}\mathbf{u}\mathbf{=}\lambda \mathbf{C}\mathbf{u} \\ {\mathbf{u}}^{\mathbf{T}}\mathbf{C}\mathbf{u}\mathbf{=}\mathbf{1} \end{gathered}$$

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

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

https://kipl.tistory.com/370

https://kipl.tistory.com/565

 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 smax = points[0].x, smin = points[0].y;
    for (int i = points.size(); i-->1; ) {
        smax = max(smax, max(points[i].x, points[i].y));
        smin = min(smin, min(points[i].x, points[i].y));
    }
    double scale = smax - smin; 
    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);
};

$\mathbf{S}$가 positive definite 행렬이고, $\mathbf{C}$는 대칭행렬일 때 아래의 일반화된 eigenvalue 문제를 푸는 방법을 알아보자. 
$$\mathbf{S}\mathbf{u}\mathbf{=}\lambda \mathbf{C}\mathbf{u}$$

타원을 피팅하는 문제에서 이런 형식의 고유값 문제에 부딛히게 된다. 이 경우 $\mathbf{S}$는 scattering matrix이고, $\mathbf{C}$는 타원피팅에 걸리는 제한조건때문에 나온다. $\mathbf{S}$가 positive definite이므로  $\mathbf{S}$의 eigenvalues $\{{\sigma }_i>0\}$는 모두 0보다 크고, eigenvector를 이용하면 

$$\mathbf{S}\mathbf{=}\mathbf{R}\mathbf{\Lambda }{\mathbf{R}}^{\mathbf{T}}\mathbf{,}\mathbf{\Lambda }\mathbf{=}\mathbf{\text{diag}}\mathbf{(}{\sigma }_{\mathbf{1}}\mathbf{,}\mathbf{.}\mathbf{.}\mathbf{.}\mathbf{,}{\sigma }_{\mathbf{n}}\mathbf{)}$$

처럼 분해할 수 있다. $\mathbf{R}$은 $\mathbf{S}$의 eigenvector를 열로 가지는 행렬로 orthogonal 행렬이다: ${\mathbf{R}}^{\mathbf{-}\mathbf{1}}\mathbf{=}{\mathbf{R}}^{\mathbf{T}}$.

이제 $\mathbf{S}$의 제곱근 행렬을 $\mathbf{Q}\mathbf{,}{\mathbf{Q}}^{\mathbf{2}}\mathbf{=}\mathbf{S}$라면 

$$\mathbf{Q}\mathbf{=}\mathbf{R}{\mathbf{\Lambda }}^{\mathbf{1}\mathbf{/}\mathbf{2}}{\mathbf{R}}^{\mathbf{T}}\mathbf{,}{\mathbf{\Lambda }}^{\mathbf{1}\mathbf{/}\mathbf{2}}\mathbf{=}\mathbf{\text{diag}}\mathbf{(}\sqrt{{\sigma }_{\mathbf{1}}}\mathbf{,}\mathbf{.}\mathbf{.}\mathbf{.}\mathbf{,}\sqrt{{\sigma }_{\mathbf{n}}}\mathbf{)}$$

임을 쉽게 확인할 수 있다. $\mathbf{Q}$을 이용하면 구하려는 고유값 문제는 

$$\mathbf{Q}\mathbf{Q}\mathbf{u}\mathbf{=}\lambda \mathbf{C}\mathbf{u}\to \mathbf{Q}\mathbf{u}\mathbf{=}\lambda {\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{C}{\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{Q}\mathbf{u}\to \mathbf{v}\mathbf{=}\lambda {\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{C}{\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{v}$$

이므로 ${\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{C}{\mathbf{Q}}^{\mathbf{-}\mathbf{1}}$의 고유값 문제 $(1/\lambda ,\mathbf{Q}\mathbf{u}\mathbf{)}$로 단순화됨을 알 수 있다. $\mathbf{Q}$의 역행렬이 ${\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{=}\mathbf{R}{\mathbf{\Lambda }}^{\mathbf{-}\mathbf{1}\mathbf{/}\mathbf{2}}{\mathbf{R}}^{\mathbf{T}}$임을 쉽게 체크할 수 있으므로 직접적으로 역행렬을 계산할 필요가 없어진다.

$${\mathbf{\Lambda }}^{\mathbf{-}\mathbf{1}\mathbf{/}\mathbf{2}}\mathbf{=}\mathbf{\text{diag}}\mathbf{(}\mathbf{1}\mathbf{/}\sqrt{{\sigma }_{\mathbf{1}}}\mathbf{,}\mathbf{.}\mathbf{.}\mathbf{.}\mathbf{,}\mathbf{1}\mathbf{/}\sqrt{{\sigma }_{\mathbf{n}}}\mathbf{)}$$