Geometry & Recognition

영상에 담고 있는 object을 간단히 근사를 할 때 타원으로 많이 기술한다(e.g: head tracking). 타원으로 기술하면 장축의 방향으로 object의 기울어진 방향을, 장축과 단축의 길이로 object의 크기를 가늠할 수 있다. object의 픽셀 분포에서 형상에 대한 정보는 2차 moment를 계산해서 얻을 수 있다. 이는 타원이 2차 곡선이기 때문에 가능하다. 그리고 질량중심을 원점으로 잡으면 2차 central moment를 계산해야 한다. 그런데 통계적인 의미를 부여하기 위해서는 central moment를 object의 픽셀로 나눈 normalized central moment로 구성한 covariance matrix를 사용하면 된다.

$$\mathrm{\Sigma }=\left(\begin{array}{cc}{\tilde{\mu }}_{20}& {\tilde{\mu }}_{11}\\ {\tilde{\mu }}_{11}& {\tilde{\mu }}_{02}\end{array}\right),{\tilde{\mu }}_{pq}\equiv \frac{{\mu }_{pq}}{{\mu }_{00}},(p+q=2)$$ 

(Note: ${\mu }_{00}$로 정규화를 하지 않더라도 문제는 없다).

$\mathrm{\Sigma }$는 영상에서 object pixel의 $x$ 뱡향 분산(${\mu }_{20}$: $x$축에 대한 회전관성), $y$ 방향 분산(${\mu }_{02}$: $y$축에 대한 회전관성), $x$-$y$의 correlation을 나타낸다. $\mathrm{\Sigma }$가 대칭행렬이므로 두 개의 음이 아닌 고윳값을 가진다. 

$$\begin{array}{c}{\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{array}$$

큰 고유값에 해당하는 고유벡터의 방향이 타원의 장축 방향에 해당하고 (픽셀 변동이 심하므로) , 작은 고윳값의 고유벡터 방향은 단축 방향이다. 그리고 고윳값은 각각 장축과 단축의 반지름의 제곱에 비례한다(${\tilde{\mu }}_{pq}$는 단위가 거리 제곱이다). object의 orientation인 타원의 장축 방향은 

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

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

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

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

로 쓸 수 있다. 그리고 이 두 축에 대한 object의 회전관성은 정의에 의해서 다음 식으로 구할 수 있다: $(\overline{x},\overline{y})=\text{center of mass}$

$$\begin{array}{rl}\text{major axis:}{I}_{\text{min}}& =\sum _{(x,y)\in \text{object}}|-(x-\overline{x})\sin \theta +(y-\overline{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{array}$$

$$\begin{array}{rl}\text{minor axis:}{I}_{\text{max}}& =\sum _{(x,y)\in \text{object}}|(x-\overline{x})\cos \theta +(y-\overline{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{array}$$

(note: object의 orientation 각은 $I_{\text{min}}$을 최소화시키는 값이다)

표준 타원의 장축 반지름이 $a$고 단축 반지름이 $b$일 때 ($x^2/a^2+y^2/b^2=1$) 2차 central moment(회전관성)는 간단한 계산에 의해서

$${\mu }_{20}^{\text{(e)}}=\frac{\pi }{4}{a}^{3}b,{\mu }_{02}^{\text{(e)}}=\frac{\pi }{4}a{b}^3,{\mu }_{11}^{\text{(e)}}=0$$

으로 구해짐을 알 수 있다.

주어진 타원이 object을 잘 피팅하려면 타원의 두 축에 대한 회전관성이 각각 object의 주축에 대한  회전관성과 같은 값을 가져야 할 것이다:

$$I_{\text{min}}={\mu }_{20}^{\text{(e)}},I_{\text{max}}={\mu }_{02}^{\text{(e)}}$$

이 두 식을 풀면 타원의 장축과 단축의 반지름을 구할 수 있다. 

$$a=(\frac{4}{\pi }{)}^{1/4}(\frac{I_{\text{max}}^3}{I_{\text{min}}}{)}^{1/8},b=(\frac{4}{\pi }{)}^{1/4}(\frac{I_{\text{min}}^3}{I_{\text{max}}}{)}^{1/8}$$

보통 object을 타원 피팅할 때 윤곽선 정보를 이용하는데, 이 방법은 윤곽선을 추출할 필요가 없어서 편리하다. 단 내부에 빈 곳이 있는 object의 경우 회전관성을 감소시키므로 좋은 결과를 기대할 수 없다. 영상이 다수의 object를 담고 있을 때는 connected component labeling을 한 후 각각의 component에 대해서 fitting을 수행하면 된다.

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

원뿔을 평면으로 잘랐을 때 나타나는 곡선인 conic section은 직교 좌표계에서 $(x,y)$에 대한 2차 형식으로 쓰인다.

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

이 conic section이 타원을 기술할 때 parameter {$a,b,c,d,e,f$}를 이용해서 타원의 중심, 장축과 단축의 길이, 그리고 회전각에 대한 공식을 구해보자. 2차 항을 행렬을 써서 표현하면

$$(x,y)\left(\begin{array}{cc}a& b/2\\ b/2& c\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\cdots =0$$

따라서, 적절한 회전 변환을 하면 두 좌표의 곱으로 주어지는 $xy$-항을 없앨 수 있다. 회전 변환을 시행하면 행렬은 eigenvalue를 성분으로 하는 대각 행렬이 된다. 회전 변환이 determinant를 보존하므로 determinant는 행렬의 두 eigenvalue의 곱으로 주어짐을 알 수 있다.

$$\begin{array}{c}\left(\begin{array}{cc}a& b/2\\ b/2& c\end{array}\right)⟹\left(\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right)\text{det}=ac-b^2/4\end{array}$$

회전 후의 방정식이 타원의 방정식(원점이 이동된)을 기술하기 위해서는 $\text{det}>0$ 이어야 한다 (회전시킨 후 식에서 $x^2$과 $y^2$의 계수는 두 eigenvalue로 주어지므로 같은 부호를 가져야 한다.)

$$F=\text{ellipse}\iff b^2-4ac<0$$

conic section $F$가 타원을 기술한다면, 다음과 같이 평행이동을 시켜서 타원의 중심 $(x_0,y_0)$이 원점에 놓이게 하고, 회전변환을 시켜서 $xy$ 항을 없애도록 하자.

$$x\to x_0+x\cos \theta -y\sin \theta ,y\to y_0+x\sin \theta +y\cos \theta$$

여기서 회전각 $\theta$는 $x$을 기준으로 측정된다. $F$에 적용해서 $xy$-항이 없어지는 조건과 1차 항이 사라지도록 하는 조건을 찾으면 타원의 중심 $(x_0,y_0)$와 $\theta$는 

$$\begin{array}{c}\tan 2\theta =\frac{b}{a-c}\\ \\ \text{ellipse center: }(x_0,y_0)=\left(\begin{array}{c}\frac{2cd-be}{b^2-4ac},\frac{2ae-bd}{b^2-4ac}\end{array}\right)\end{array}$$

로 주어짐을 확인할 수 있다. Eigenvalue의 제곱근의 역수가 두 축의 반지름을 결정하므로 음수가 되어서는 안된다 (둘 중 하나가 0이면 직선이고, 둘 모두 0이면 한 점에 해당). 이는 대칭행렬의 고유값이 항상 0보다 작지 않다는 사실에서 기인한다. 위에서 구한 회전각 $\theta$를 이용해서 두 고유값을 표현하면 ,

$${\lambda }_1=a{\cos }^2\theta +b\cos \theta \sin \theta +c{\sin }^2\theta$$

$${\lambda }_2=a{\sin }^2\theta -b\cos \theta \sin \theta +c{\cos }^2\theta$$

위의 회전변환식을 대입했을 떄 나머지 상수항은 (첫 번째 등호는 계산해서 확인할 수 있음)

$$a{x}_0^2+b{x}_0{y}_0^2+c{y}_0^2+d{x}_0+e{y}_0+f=f-(a{x}_0^2+b{x}_0{y}_0+c{y}_0^2)\equiv -{\text{scale}}^{-1}$$

로 주어짐을 알 수 있다. 따라서 회전시킨 타원은 표준형 꼴

$${\lambda }_1{x}^2+{\lambda }_2{y}^2={\text{scale}}^{-1}$$

로 표현된다. 이 표준형 타원의 두 축의 반지름은 각각

$$r_x=\sqrt{\frac{{\text{scale}}^{-1}}{{\lambda }_1}},r_y=\sqrt{\frac{{\text{scale}}^{-1}}{{\lambda }_2}}$$ 로 주어진다.

// conic_params(a, b, c, d, e, f): ax^2 + bxy + cy^2 + dx + ey + f = 0;
// ellipse_params(radius_x, radius_y, center_x, center_y, tilt_angle w.r.t x-axis);
bool conic_to_ellipse(double conic_params[6], double ellipse_params[5]) {
    const double a = conic_params[0];
    const double b = conic_params[1];
    const double c = conic_params[2];
    const double d = conic_params[3];
    const double e = conic_params[4];
    const double f = conic_params[5];
    // get ellipse orientation w.r.t x-axis;
    const double theta = 0.5 * atan2(b, a - c);
    // get scaled x/y radius;
    const double ct = cos(theta);
    const double st = sin(theta);
    const double ap = a * ct * ct + b * ct * st + c * st * st;
    const double cp = a * st * st - b * ct * st + c * ct * ct;
    // get center of ellipse;
    const double cx = (2 * c * d - b * e) / (b * b - 4 * a * c);
    const double cy = (2 * a * e - b * d) / (b * b - 4 * a * c);
    // get scale factor
    const double val = a * cx * cx + b * cx * cy + c * cy * cy;
    const double scale_inv = val - f;
    if (scale_inv / ap <= 0 || scale_inv / cp <= 0) {
        TRACE("Error! ellipse parameters are imaginary\n");
        return 0;
    }
    ellipse_params[0] = sqrt(scale_inv / ap);
    ellipse_params[1] = sqrt(scale_inv / cp);
    ellipse_params[2] = cx;
    ellipse_params[3] = cy;
    ellipse_params[4] = theta;
    return 1;
};