Best-fit Ellipse

영상에 담고 있는 object를 타원으로 근사할 때 타원으로 많이 기술한다(e.g: head tracking). 타원으로 기술하면 장축의 방향으로 object의 기울어진 방향(orientation)을, 장축과 단축의 길이로 object의 크기를 가늠할 수 있다. object의 형상 정보는 픽셀 분포의 2차 moment로부터 얻을 수 있다. 이는 타원이 2차 곡선으로 표현되므로 object의 2차 moment만으로 주요 형상 정보를 효과적으로 기술할 수 있기 때문이다. 그리고 질량중심(centroid)을 원점으로 잡은 후 2차 central moment를 계산해야 한다. 이를 object의 픽셀 개수(또는 면적)로 나누면 통계적 의미의 분산과 공분산을 얻을 수 있으며, 다음과 같은 covariance matrix를 구성할 수 있다. 

$$\Sigma = \left(\begin{array}{cc} \tilde\mu_{20} &  \tilde\mu_{11} \\ \tilde\mu_{11} & \tilde\mu_{02} \end{array} \right), \quad\quad\tilde\mu_{pq} \equiv \frac{\mu_{pq}}{\mu_{00}},  \quad (p+q=2)$$(Note: $\mu_{00}$로 정규화를 하지 않아도 이후의 orientation 계산에는 영향을 주지 않는다).

여기서 $\Sigma$는 영상에서 object 픽셀들의 $x$ 방향 분산(${\tilde{\mu}}_{20}$), $y$ 방향 분산 (${\tilde{\mu}}_{02}$) , $x$-$y$ covariance( (${\tilde{\mu}}_{11}$) )을 나타낸다. $\Sigma$가 대칭행렬이므로 두 개의 음이 아닌 고윳값을 가지며,

\begin{gather} \lambda_1 = \frac{\tilde\mu_{20} + \tilde\mu_{02}}{2} + \frac{ \sqrt{(\tilde\mu_{20} - \tilde\mu_{02})^2 + 4\tilde\mu_{11}^2 } }{2}, \\  \lambda_2 = \frac{\tilde\mu_{20} + \tilde\mu_{02}}{2} - \frac{ \sqrt{(\tilde\mu_{20} - \tilde\mu_{02})^2 + 4\tilde\mu_{11}^2 } }{2}\end{gather}로 주어진다.

큰 고유값에 해당하는 고유벡터는 픽셀 분포가 가장 넓게 퍼져 있는 방향, 즉 타원의 장축 방향이 되고, 작은 고윳값의 고유벡터 방향은 단축 방향이다. 그리고 균일한 밀도를 갖는 타원에서는 고윳값은 각각 장축과 단축의 반지름의 제곱에 비례한다($\tilde\mu_{pq}$는 단위가 거리 제곱이다). object의 orientation,  타원의 장축 방향은 

$$ \theta = \frac{1}{2} \tan^{-1} \Big( \frac{   2\tilde\mu_{11} }{\tilde\mu_{20} - \tilde\mu_{02}}\Big)= \frac{1}{2} \tan^{-1} \Big( \frac{   2\mu_{11} }{\mu_{20} - \mu_{02}}\Big) $$

로 계산된다. (https://kipl.tistory.com/58)

 

타원의 orientation 각도를 구했으므로 두 주축을 나타내는 단위벡터는

$$\text{major axis: }(\cos \theta, \sin \theta), \quad \text{minor axis: }(-\sin\theta, \cos \theta)$$

로 쓸 수 있다. 그리고 이 두 축에 대한 object의 회전관성은 회전축으로부터의 거리 제곱의 합으로 정의되므 다음 식으로 구할 수 있다: $(\bar{x}, \bar{y})=\text{center of mass}$

\begin{align}\text{major axis:} ~I_\text{min} &= \sum_{(x, y )\in \text{object}} | - (x-\bar{x})\sin \theta + (y - \bar{y})\cos \theta |^2 \\ &=\frac{\mu_{20}+\mu_{02}}{2} -\frac{\mu_{20} - \mu_{02}}{2} \cos (2\theta)-  \mu_{11} \sin(2\theta)\end{align}

\begin{align}\text{minor axis:} ~I_\text{max} &= \sum_{(x, y)\in \text{object}} | (x-\bar{x}) \cos \theta + (y - \bar{y}) \sin \theta|^2 \\&= \frac{\mu_{20}+\mu_{02}}{2}+\frac{\mu_{20}-\mu_{02}}{2}\cos (2\theta)+ \mu_{11}\sin(2\theta)\end{align} 장축 방향은 object가 가장 길게 퍼져 있는 방향이지만, 회전축으로 생각하면 대부분의 픽셀이 축 가까이에 위치하므로 회전관성은 가장 작다. 따라서 orientation 각도 $\theta$는 $I_\text{min}$을 최소로 만드는 방향이 된다.

 

장축 반지름이 $a$고 단축 반지름이 $b$인 표준 타원$$\frac{x^2}{a^2} + \frac{y^2}{b^2} =1$$의 2차 central moment는 간단한 계산에 의해서

$$\mu_{20}^\text{(e)} = \frac{\pi}{4} a^3 b, \quad \mu_{02}^\text{(e)} = \frac{\pi}{4} ab^3,\quad \mu_{11}^\text{(e)}= 0$$이다. 따라서 주어진 타원이 object 잘 피팅하려면 타원의 두 축에 대한 회전관성이 각각 object의 주축에 대한  회전관성과 각각 같아야 하므로,

$$ I_\text{min} = \mu_{02}^\text{(e)}, \quad I_\text{max} =\mu_{20}^\text{(e)}$$를 만족해야 한다.

이 두 식을 풀면 타원의 장축과 단축의 반지름은

$$a = \Big( \frac{4}{\pi}\Big)^{1/4} \Big( \frac{I_\text{max}^3}{I_\text{min}}\Big)^{1/8}, \quad b = \Big( \frac{4}{\pi}\Big)^{1/4} \Big( \frac{I_\text{min}^3}{I_\text{max}}\Big)^{1/8}$$으로 구할 수 있다.

 

일반적으로 object를 타원으로 근사할 때는 윤곽선(contour)을 추출한 후 타원을 피팅하는 방법이 많이 사용된다. 반면 위 방법은 object 내부의 픽셀 분포만으로 2차 moment를 계산하므로 윤곽선을 추출할 필요가 없다는 장점이 있다. 다만 내부에 큰 빈 영역(hole)이 존재하는 object에서는 회전관성이 감소하여 실제 형상을 충분히 반영하지 못할 수 있으므로 피팅 결과가 부정확해질 수 있다. 영상에 여러 개의 object가 포함되어 있는 경우에는 먼저 connected component labeling을 수행한 후 각 connected component에 대해 독립적으로 위 과정을 적용하면 된다.

void getEllipse(CRaster& raster) {
    const double four_pi = 1.0 / atan(1.0);
    CSize sz = raster.GetSize();
    double xsum = 0, ysum = 0;
    double x2sum = 0, y2sum = 0, xysum = 0;
    int count = 0;
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            if (*p++) {
                count++;
                xsum += x;      ysum += y;
                x2sum += x * x; y2sum += y * y; 
                xysum += x * y;
            }
        }
    }
    if (!count) return;
    double xm = xsum / count;
    double ym = ysum / count;
    double m20 = x2sum - xm * xm * count;
    double m02 = y2sum - ym * ym * count;
    double m11 = xysum - xm * ym * count;
    double theta2 = atan2(2 * m11, m20 - m02);
    double ct = cos(theta2), st = sin(theta2);
    double Imin = 0.5 * (m20 + m02) - 0.5 * (m20 - m02) * ct - m11 * st;
    double Imax = 0.5 * (m20 + m02) + 0.5 * (m20 - m02) * ct + m11 * st;
    double major = pow(four_pi, 0.25) * pow(Imax * Imax * Imax / Imin, 0.125);
    double minor = pow(four_pi, 0.25) * pow(Imin * Imin * Imin / Imax, 0.125);
    drawEllipse(raster, mx, my, major, minor, theta2 / 2);
}