Least Squares Fitting of Ellipses

일반적인 이차곡선은 다음의 이차식으로 표현이 된다:

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

6개의 계수는 모두 독립적이지 않고 어떤 종류의 이차곡선인가에 따라 제약조건이 들어온다. 주어진 점 데이터 $\{ (x_i, y_i) | i=1,2,...,n\}$를 피팅하는 이차곡선을 각각의 데이터 점에서 대수적 거리 ($=|F(x_i, y_i)|$)의 제곱을 최소로 하는 조건하에서 찾도록 하자:

$$ \epsilon = \sum_{i} |F(x_i , y_i )|^2 \rightarrow \text{min}$$

타원으로 피팅하는 경우 여러 가지 제약조건을 줄 수 있지만(e.g.: $a+c=1$, $\sqrt{a^2+b^2+...+f^2 }=1$ 등등) 여기서는 이차 곡선이 타원이기 위해서는 2차 항의 계수로 만드는 행렬 $\left(\begin{array}{cc} a & b/2\\b/2 & c\end{array}\right)$의 determinant (= $ac- b^2/4>0$)가 양수이어야 한다는 조건에서 (타원을 회전+평행 이동시키면 표준형 타원으로 만들 수 있는데, 이때 두 주축의 계수 곱이 determinant다. 따라서 타원이기 위해서는 값이 양수여야 된다) 

$$ 4ac - b^2 = 1$$

로 선택하도록 하자. 이 제약조건을 넣은 타원 피팅은 다음 식을 최소화하는 계수 벡터 $\mathbf{a}=[a,b,c,d,e,f]^T$을 찾는 문제가 된다:

\begin{gather}    \epsilon = \mathbf{a}^T. \mathbf{S}. \mathbf{a}   -\lambda ( \mathbf{a}^T.\mathbf{C}.\mathbf{a}-1)\end{gather}

여기서 $6\times 6$ scattering matrix $\mathbf{S}$는

$$ \mathbf{S}= \mathbf{D}^T. \mathbf{D}.$$

$n\times 6$ design matrix $\mathbf{D}$는 

$$ \mathbf{D}=\left[ \begin{array}{cccccc} 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\\ & & \vdots & &  \\ x_n^2& x_n y_n & y_n^2 & x_n & y_n & 1 \end{array}\right] .$$

그리고, $6\times 6$ 행렬 $\mathbf{C}$는

$$ \mathbf{C}= \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]  . $$

$\epsilon$을 $\mathbf{a}^T$에 대해 미분하면 최소 제곱 문제는 Lagrange multiplier를 고윳값으로 하는 고유 방정식의 해를 구하는 문제로 환원이 된다.

$$ \mathbf{S}.\mathbf{a} = \lambda \mathbf{C}.\mathbf{a}$$

그리고 주어진 고유값 $\lambda$에 해당하는 normalized 고유 벡터가 $\mathbf{a}$일 때 

$$\epsilon = \mathbf{a}^T. \mathbf{S}.\mathbf{a} = \lambda$$

이므로 최소의 양의 고윳값에 해당하는 고유 벡터가 해가 된다. Silverster의 law of inertia를 이용하면 위의 고유 방정식에서 양의 고윳값은 딱 1개만 존재함을 보일 수 있다. 고유값 $\lambda$에 해당하는 고유 벡터가 $\mathbf{a}$일 때 임의의 상수 $\mu$ 배를 한 $\mu\mathbf{a}$도 같은 고유값을 갖는 고유벡터다. normalized 고유벡터는 제약조건 $\mu^2 \mathbf{a}^T . \mathbf{C}. \mathbf{a} =1$을 만족하도록 크기를 조정하면 $ \mathbf{\tilde{a}} = \mathbf{a} / \sqrt{\mathbf{a}^T. \mathbf{C}. \mathbf{a}}$로 주어진다. 

 

이 일반화된 고유방정식의 풀이는 먼저 positive definite인 $\bf S$의 제곱근 행렬 $\bf Q=S^{1/2}$을 이용하면 쉽다. $\bf S$의 고유벡터를 구하면 eigendecomposition에 의해 $\bf  S = R\Lambda R^T$로 쓸 수 있으므로 제곱근 행렬은 $\bf  Q = R \Lambda ^{1/2} R^T$임을 쉽게 확인할 수 있다. 원래의 고유값 문제를 $\tt Q$을 이용해서 표현하면 

$$ \bf Q a = \lambda Q^{-1} C a = \lambda Q^{-1} C Q^{-1} Qa$$ 이므로 더 다루기 쉬운 $\bf Q^{-1} C Q^{-1}$의 고유값 문제로 환원이 됨을 알 수 있다.

 

추가로, 고유 방정식은 $\mathbf{C}$가 $3\times 3$ 크기의 block으로 나누어질 수 있으므로 $[a,b,c]^T$ 에 대한 고유 방정식으로 줄일 수 있어서 쉽게 해결할 수 있다. 물론 scattering matrix을 구성하는 moment의 개수가 많아서 matrix 연산 패키지를 사용하지 않는 경우 코드가 길어지게 된다. 

 

Note: 이 내용은 다음 논문을 정리한 것이다. Andrew Fitzgibbon, Maurizio Pilu, and Robert B. Fisher, Direct Least Square Fitting of Ellipses". IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 21, NO. 5, MAY 1999. 

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

 

구현은  https://kipl.tistory.com/566

Ellipse Fitting