Geometry & Recognition

주어진 데이터를 가장 잘 피팅하는 직선을 찾기 위해서는 일반적으로 각 데이터 점의 \(y\)값과 직선이 예측하는 \(y\)값의 차이, 즉 residual의 제곱합을 최소화하는 최소자승법을 사용한다. 이 방법은 계산이 간단하고 대부분의 경우 좋은 결과를 제공하지만, 오차를 \(y\)축 방향으로만 측정한다는 한계가 있다. 따라서 직선의 기울기가 매우 큰 경우에는 실제 점과 직선 사이의 최단거리는 작더라도 \(y\,\)방향 residual은 크게 나타날 수 있으며, 데이터에 포함된 outlier의 영향도 더욱 커져 올바른 직선 피팅을 방해할 수 있다.

이러한 문제를 완화하는 간단한 방법 중 하나는 \(y\,\)방향 residual 대신 데이터 점과 직선 사이의 수직거리(최단거리)를 오차로 사용하는 것이다. 이 경우 오차는 직선의 기울기에 의존하지 않는 기하학적 거리로 정의되므로, 기울기가 큰 직선에서도 보다 안정적인 피팅 결과를 얻을 수 있다. 이러한 방법을 Orthogonal Least Squares 또는 Total Least Squares(TLS)라 하며, 특히 \(x\)축과 \(y\)축 모두에 측정 오차가 존재하는 경우에 적합하다.

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

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

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

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

이다. 따라서 주어진 데이터에서 떨어진 거리 제곱의 합이 최소인 직선을 구하기 위해서는 다음을 최소화시키는 $(\theta, s)$을 구해야 한다. $$ L = \sum_i d_i^2 = \sum_i \big( x_i \sin \theta - y_i \cos \theta +s \big)^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_i x_i y_i + s \cos\theta \sum_i x_i + s \sin \theta \sum_i y_i = 0$$

$$ \frac{1}{2}\frac{\partial L}{\partial s}=\sin \theta \sum_i x_i -\cos \theta \sum_i y_i  + N s=0$$

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

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

두 번째 식은 원점에서 직선까지 거리 $s$가 0임을 의미하므로 직선이 질량중심(질량중심계에서 원점)을 통과함을 의미한다. 첫번째 식을 풀면

$$ \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$$이 정해진다. 질량중심계에서는 원점을 지나지만 원좌표계로 돌아오면 데이터의 질량중심을 통과하도록 평행이동시키면 된다.

$$  \boxed{\left(-A+ \sqrt{A^2+ (2B)^2} \right)  (x-\bar{x}) = 2B (y - \bar{y})   }  $$

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

$$ \bar{x} = \frac{1}{N} \sum_i x_i, \quad \bar{y} = \frac{1}{N} \sum_i y_i$$

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

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

단, $B=0$이면 위 식을 그대로 사용하지 말고 별도로 처리해야 한다. $A>0$이면 최적 직선은 $y=\bar y$인 수평선이고, $A<0$이면 $x=\bar x$인 수직선이다. 특히 $A=B=0$이면 중심화한 데이터의 $x$, $y$ 방향 분산이 같고 공분산이 0이다. 즉 데이터가 2차 모멘트 수준에서 모든 방향으로 동일하게 퍼져 있어 특별한 주축이 존재하지 않는다. 이 경우 질량중심을 지나는 모든 직선이 동일한 수직거리 제곱합을 가지므로, 유일한 최적 직선을 정할 수 없다. 원 위에 균등하게 놓인 점들, 정사각형 또는 정삼각형의 꼭짓점들이 대표적인 예다. 이 결과는 PCA를 적용한 결과와 동일하다. 중심화한 데이터의 공분산 행렬에서 가장 큰 고윳값에 대응하는 고유벡터가 최적 직선의 방향을 결정하며, 가장 작은 고윳값의 고유벡터는 그 직선의 법선 방향을 결정한다. 따라서 PCA 기반 피팅은 수직선을 포함한 모든 방향의 직선을 자연스럽게 처리한다.  https://kipl.tistory.com/211.

다만 TLS 역시 거리의 제곱합을 최소화하므로 outlier에 본질적으로 강인한 방법은 아니다. outlier의 영향을 줄이려면 Huber 또는 Tukey 가중치를 이용한 반복 가중 최소자승법, 혹은 RANSAC과 같은 robust fitting 방법을 사용해야 한다. 또한 통계에서 Deming regression은 일반적으로 $x$, $y$ 측정오차의 분산비를 반영하는 직교회귀를 가리킨다. 여기서 다룬 TLS는 두 축의 오차 분산이 같다고 보는 Deming regression의 특수한 경우로 해석할 수 있다.

일반적인 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$를 구하는 문제는 scatter matrix $\bf S=D^T D$에 대한 일반화된 고유값 문제로 환원이 된다.

$$  \bf S u =\lambda C u, \qquad 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

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;		
    }			

    /* scatter 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 scatter 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);
};

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

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

$$ \bf S = R \Lambda R^T, ~~~~\Lambda=\text{diag}(\sigma_1, ...,\sigma_n)$$

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

이제 $\bf S$의 제곱근 행렬을 $\bf Q, ~Q^2 =S$라면 

$$ \bf Q= R \Lambda^{1/2} R^T,~~~~\Lambda^{1/2} = \text{diag}( \sqrt{\sigma_1},...,\sqrt{\sigma_n})$$

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

$$ \bf QQ u = \lambda Cu ~\to ~  Qu = \lambda Q^{-1} C Q^{-1}Qu~~\to ~~ v = \lambda Q^{-1} C Q^{-1} v$$

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

$$\bf \Lambda^{-1/2} = \text{diag}( 1/\sqrt{\sigma_1},...,1/\sqrt{\sigma_n})$$