Geometry & Recognition

RGB 컬러 이미지의 gray level의 cdf를 이용해서 histogram equalization을 수행한다. 컬러 이미지의 gray level은 

$$gray = \frac{r + g+ b}{3}$$

으로 정의한다.

void test_color_equalize_new(CRaster& raster, CRaster& out) {
    int hist[256] = {0};
    int chist[256], lut[256], partition[256 + 1];
    CSize sz = raster.GetSize();
    out = raster; // clone; 
    // pdf of gray level: g = (r + g + b) / 3
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++){
            int a = *p++; a += *p++; a += *p++;
            hist[a / 3]++;
        }
    };
    // cdf;
    chist[0] = hist[0];
    for (int i = 1; i < 256; i++)
        chist[i] = hist[i] + chist[i - 1];
    
    /* assign equal number of pixels in each gray levels;*/
    int num_pixels = sz.cx * sz.cy;
    double pixels_per_level = double(num_pixels) / 256;

    /* first and last in partition */
    partition[0]   = 0;
    partition[256] = 256;
    /* intermediate; */
    for (int j = 0, i = 1; i < 256; i++) {
        double desired = i * pixels_per_level;			
        while (chist[j + 1] <= desired) j++;
        /* nearest interpolation */
        if ((desired - chist[j]) < (chist[j + 1] - desired)) partition[i] = j;
        else partition[i] = j + 1;
    } 
    /* fill equalization LUT */
    for (int j = 0; j < 256; j++) {
        int i = 0;
        while (partition[i + 1] <= j) i++;
        lut[j] = i;
    } 
    // needs hue preserving processes;
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *q = (BYTE *)out.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            *q++ = lut[*p++]; *q++ = lut[*p++]; *q++ = lut[*p++];
        }
    }
}

void test_color_equalize_HSV(CRaster& raster, CRaster& out) {
    int hist[256] = {0};
    int chist[256] = {0}, lut[256] = {0}, partition[256 + 1];
    CSize sz = raster.GetSize();
    out = raster; // cloning;
    std::vector<double> fH(sz.cx * sz.cy);
    std::vector<double> fV(sz.cx * sz.cy);
    std::vector<double> fS(sz.cx * sz.cy);

    int n = sz.cx * sz.cy;
    double r, g, b;
    for (int y = 0, k = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++, k++){
            b = *p++, g = *p++, r = *p++;
            RGBtoHSV(r / 255., g / 255., b / 255., fH[k], fS[k], fV[k]);
        }
    };	
    // make histogram of V-color;
    for (int k = 0; k < n; k++)
        hist[int(fV[k] * 255)]++;

    // cdf;
    chist[0] = hist[0];
    for (int i = 1; i < 256; i++) {
        chist[i] = hist[i] + chist[i - 1];
    }
    /* assign equal number of pixels in each V-color;*/
    int num_pixels = sz.cx * sz.cy;
    double pixels_per_level = double(num_pixels) / 256;

    /* first and last in partition */
    partition[0]   = 0;
    partition[256] = 256;
    /* intermediate; */
    for (int j = 0, i = 1; i < 256; i++) {
        double desired = i * pixels_per_level;			
        while (chist[j + 1] <= desired) j++;
        /* nearest interpolation */
        if ((desired - chist[j]) < (chist[j + 1] - desired)) partition[i] = j;
        else partition[i] = j + 1;
    } 
    /* fill equalization LUT */
    for (int j = 0; j < 256; j++) {
        int i = 0;
        while (partition[i + 1] <= j) i++;
        lut[j] = i;
    } 
    for (int k = 0; k < n; k++)
        fV[k]= lut[int(fV[k] * 255)] / 255.;

    for (int y = 0, k = 0; y < sz.cy; y++) {
        BYTE *q = (BYTE *)out.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++, k++) {
            HSVtoRGB(fH[k], fS[k], fV[k], r, g, b);
            *q++ = int(b * 255); 
            *q++ = int(g * 255);
            *q++ = int(r * 255);
        }
    }
}

/* fR Red component, range: [0, 1]
** fG Green component, range: [0, 1]
** fB Blue component, range: [0, 1]
** fH Hue component, range: [0, 360]
** fS Hue component, range: [0, 1]
** fV Hue component, range: [0, 1] */
void RGBtoHSV(double fR, double fG, double fB, double& fH, double& fS, double& fV);
void HSVtoRGB(double fH, double fS, double fV, double& fR, double& fG, double& fB);

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

$$ 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

주어진 점집합을 원으로 피팅하기 위해 이차식

$$A(x^2 + y^2) + Bx + Cy +D =0$$

을 이용하자. 원의 경우는 $A=0$인 경우는 직선을 나타내고, $A\ne0$인 경우가 원을 표현한다. 물론 $A=1$로 설정을 할 수 있으나 이 경우는 주어진 데이터가 원의 일부 부분만 나타내는 경우 문제가 생긴다. 원은 중심과 반지름을 알면 결정되므로 3개의 변수가 필요한데 위의 이차 다항식은 4개의 변수를 가진다. 원을 제대로 표현하기 위해서는 제한 조건이 들어오는데 

$$ \Big(x + \frac {B}{2A} \Big)^2 + \Big(y + \frac {C}{2A} \Big)^2 = \frac {B^2 +C^2 -4A  D}{4A^2}$$

으로 쓸 수 있으므로 좌변이 양수이므로 

$$ B^2+C^2 - 4A  D>0$$이어야 한다. 따라서 $$ B^2 + C^2 - 4A  D=1$$인 제한 조건을 둘 수 있다.

주어진 각 점들이 정확히 추정하는 원 상에 있지 않으므로 이차식의 좌변의 값이 반드시 0이 되지는 않는다. 잘 피팅하기 위해서는 제한조건을 만족하면 각 점들에서 차이를 최소로 만드는 파라미터를 찾으면 된다:

$$ \epsilon(A,B, C, D)= \sum | A(x_i ^2 + y_i^2) + Bx_i + Cy_i +D |^2 - \lambda (B^2 +C^2 - 4A D-1) \longrightarrow \text {min}$$

여기서 $\lambda$는 Lagrange multiplier이다. 식을 행렬로 표현하기 위해

$$z_i= x_i^2 + y_i^2 , ~M_{zz}= \frac {1}{n} \sum z_i^2 , ~M_{xx} = \frac{1}{n} \sum x_i^2, ~M_{yy}=\frac{1}{n}\sum y_i^2,  $$

$$M_{zx}= \frac {1}{n} \sum z_i x_i  , ~M_{zy} = \frac{1}{n} \sum z_i y_i , ~M_{xy}=\frac{1}{n}\sum x_i y_i, $$

$$M_{z}= M_{xx}+M_{yy} , ~M_{x} = \frac{1}{n} \sum x_i, ~M_{y}=\frac {1}{n}\sum y_i,  $$

$${\tt x}= \left(\begin {array}{c} A \\ B \\ C \\D\end {array}\right) , ~ \tt M = \left( \begin{array}{cccc} M_{zz}& M_{zx}& M_{zy} & M_z \\ M_{zx} & M_{xx}& M_{xy} & M_{x}\\ M_{zy} & M_{xy} & M_{yy} &M_y \\M_z &M_x &M_y& 1 \end{array} \right), ~\tt  N=\left( \begin{array}{cccc} 0& 0& 0& -2 \\ 0& 1& 0 & 0\\ 0& 0& 1& 0\\ -2&0&0& 0 \end{array} \right)$$

로 놓으면 ($\tt M=\text {scattering matrix}$) $$ n^{-1} \epsilon= {\tt x}^T. {\tt M}. {\tt x} - \lambda ( {\tt x}^T. {\tt N}. {\tt x} - 1)$$

로 표현된다. 질량중심 좌표계에서 계산을 하면 $M_x= M_y =0$이므로 행렬이 더 간단해진다.

$\epsilon$을 $\tt x^T$에 대해서 미분하면

$$ \tt M.x = \lambda \tt N.x$$

인 lagrange multiplier가 고윳값이 되는 일반화된 고윳값 방정식을 얻는다. 이 고윳값 식이 해를 가지려면

$$\text {det}( \tt M -\lambda \tt N)=0$$

이어야 한다. 이 특성 방정식은 4차 다항식이므로 closed form 형태의 해를 쓸 수 있으나 실질적으로 수치해석적으로 구하는 것이 더 편하다(Newton-Raphson). $\text {min}(\epsilon)$은 제한조건 때문에

$$ n^{-1}\epsilon = \lambda \tt x^T. N. x = \lambda > 0$$

이므로 $0$보다 큰 최소 고윳값이 원하는 해이다. ($\lambda = 0$은 $\text {det}(\tt M)=0$, 즉, $\tt M$이 singular 한 경우)

그런데 주어진 고윳값 $\lambda$ 해당하는 고유 벡터가 $\tt x$일 때, 임의의 0이 아닌 상수 $\mu$에 대해 $\mu \tt x$도 역시 고유벡터이다. 고유 벡터의 크기를 고정하기 위해서 제한조건을 적용하면 $1=\mu^2 \tt x^T .N .x = \mu^2 \tt x^T. M. x/\lambda$이므로 $\mu = \sqrt{\lambda/ \tt x^T. M. x}$로 선택하면 된다. 그렇지만 해가 원을 나타내는 경우 원의 중심 $(-B/2A, -C/2A)$, 반지름 $\sqrt{ (B/2A)^2 + (C/2A)^2 - D/A}$은 계수의 비에만 의존하므로 굳이 normalized 된 고유벡터를 구할 필요가 없다. 따라서, 해 중에 직선인 경우는 무시하고($A=0$) 원만을 선택할 때는 $A=1$로 놓고 계산을 해도 영향을 받지 않는다. 또한 이 경우 고유벡터 성분을 구체적으로 적을 수 있어 수치해석에 의존할 필요도 없어진다.

최소의 고윳값에 대해서 고유 방정식을 풀면 고유벡터의 성분은

$$ C_{xy} \equiv  M_{xx} M_{yy}- M_{xy}^2 , ~~\Delta \equiv \lambda^2 - M_z   \lambda +C_{xy}$$

$$A=1$$

$$B = \frac {M_{zx}(M_{yy} - \lambda)-M_{zy} M_{xy}}{ \Delta }$$

$$C = \frac {M_{zy}(M_{xx} - \lambda)-M_{zx} M_{xy}}{ \Delta }$$

$$D=-M_z -2\lambda$$로 주어진다. 따라서 원의 중심은

$$ c_x = -\frac {B}{2A}= \frac{    M_{zx}(M_{yy} -\lambda) - M_{zy}M_{xy}}{2  \Delta}, \quad c_y= -\frac{C}{2A}= \frac{ M_{zy}(M_{xx} -\lambda )- M_{zx} M_{xy}}{2  \Delta} $$

원의 반지름은 

$$ r =\sqrt { (B/2A)^2+(C/2A)^2 -D/A} = \sqrt{ c_x^2 + c_y^2 +M_{z} +2\lambda}$$ 

로 주어진다.

** 추가: Sylvester's law of inertia에 의하면 위에서 구하는 고유값 $\lambda$의 부호 개수는  $\tt N$의 고유값 부호별 개수와 같아야 하는데, $\tt N$의 고유값이 $\{-2,1,1,2\}$이므로 양수인 고유값 $\lambda$는 3개가 존재한다.

Ref: V. Pratt, Direct least-squares fitting of algebraic surfaces, Computer Graphics 21, 1987, 145–152.

double circleFit(std::vector<CPoint>& data, double& centerx, double& centery, double& radius) {
    const int maxIter = 99;
    double mx = 0, my = 0;
    for (int i = data.size(); i-- > 0;)
        mx += data[i].x, my += data[i].y;

    // center of mass;
    mx /= data.size(); my /= data.size();

    // moment calculation;
    double Mxy, Mxx, Myy, Mzx, Mzy, Mzz;
    Mxx = Myy = Mxy = Mzx = Mzy = Mzz = 0.;
    for (int i = data.size(); i-- > 0;) {
        double xi = data[i].x - mx;   //  center of mass coordinate
        double yi = data[i].y - my;   //  center of mass coordinate
        double zi = xi * xi + yi * yi;
        Mxy += xi * yi; Mxx += xi * xi;
        Myy += yi * yi; Mzx += xi * zi;
        Mzy += yi * zi; Mzz += zi * zi;
    }
    Mxx /= data.size(); Myy /= data.size();
    Mxy /= data.size(); Mzx /= data.size();
    Mzy /= data.size(); Mzz /= data.size();
    double Mz = Mxx + Myy;
    double Cxy = Mxx * Myy - Mxy * Mxy;
    double Vz = Mzz - Mz * Mz;
    // coefficients of characteristic polynomial;
    double C2 = 4 * Cxy - 3 * Mz * Mz - Mzz; 
    double C1 = Vz * Mz + 4 * Cxy * Mz - Mzx * Mzx - Mzy * Mzy;
    double C0 = Mzx * (Mzx * Myy - Mzy * Mxy) + Mzy * (Mzy * Mxx - Mzx * Mxy) - Vz * Cxy;
    // Newton's method starting at lambda = 0
    double lambda = 0;
    double y = C0;  // det(lambda = 0)
    for (int iter = 0; iter < maxIter; iter++) {
        double Dy = C1 + lambda * (2. * C2 + 16.*lambda * lambda);
        double lambdaNew = lambda - y / Dy;
        if ((lambdaNew == lambda)) break;
        double ynew = C0 + lambdaNew * (C1 + lambdaNew * (C2 + 4 * lambdaNew * lambdaNew));
        if (fabs(ynew) >= fabs(y))  break;
        lambda = lambdaNew;
        y = ynew;
    }
    double DEL = lambda * lambda - lambda * Mz + Cxy;
    double cx = (Mzx * (Myy - lambda) - Mzy * Mxy) / DEL / 2;
    double cy = (Mzy * (Mxx - lambda) - Mzx * Mxy) / DEL / 2;
    centerx = cx + mx;   // recover origianl coordinate;
    centery = cy + my;
    radius = sqrt(cx * cx + cy * cy + Mz + 2 * lambda);
    return fitError(data, centerx, centery, radius);
}

표준형 타원

$$ \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$

의 두 주축에 대한 central moment (moment of inertia)는 각각

$$ \mu_{20} = \frac{\pi}{4} a^3 b, \quad \mu_{02} = \frac{\pi}{4} ab^3$$ 

으로 주어짐은 쉽게 계산할 수 있다. 그런데 타원의 면적인 0차 central moment가 $\mu_{00}= \pi ab$이므로 normalized central moment (위치의 분산을 의미한다)는 각각

$$\tilde\mu_{20} = \frac{\mu_{20}}{\mu_{00}} = \frac{1}{4}a^2, \quad  \tilde\mu_{02} = \frac{\mu_{02}}{\mu_{00}} = \frac{1}{4}b^2$$

따라서 장축과 단축의 반지름은 주축에 대한 2차 normalized central moment을 구하면 얻을 수 있다.

기울어진 타원의 경우는 $\tilde\mu_{pq}$가 픽셀의 분포를 알려주므로 이를 이용한 covariance 행렬의 고유값을 구하면 장축과 단축의 반지름을 알 수 있고, 기울어진 정도도 알 수 있다. 공변행렬

$$ \Sigma = \left( \begin{array}{cc} \tilde\mu_{20} & \tilde\mu_{11} \\ \tilde\mu_{11} & \tilde\mu_{02} \end{array} \right)$$

의 두 고윳값은

$$ \lambda_1 = \frac{1}{2}\left( \tilde\mu_{20}+ \tilde\mu_{02} + \sqrt{ (\tilde\mu_{20} - \tilde\mu_{02})^2 + 4 \tilde\mu_{11}^2} \right),$$

$$ \lambda_2 = \frac{1}{2}\left( \tilde\mu_{20}+ \tilde\mu_{02} - \sqrt{ (\tilde\mu_{20} - \tilde\mu_{02})^2 + 4 \tilde\mu_{11}^2} \right).$$

따라서 기울어진 타원의 장축과 단축의 반지름, 그리고 회전각은

$$ a = 2\sqrt{ \lambda_1}, \quad b = 2\sqrt{ \lambda_2}, \quad \tan(2\theta)=\frac{2 \tilde\mu_{11}}{ \tilde\mu_{20}-\tilde\mu_{02}   }$$

로 주어지므로 회전된 좌표계 $(u, v)$에서 타원의 방정식은

$$ \frac{ u^2 }{4 \lambda_1} + \frac{v^2 }{4\lambda_2} = 1$$

로 쓰인다. 

다시 원 좌표계 $(x, y)$로 돌아가기 위해서 (물론 질량중심이 원점인 좌표계이다)

$$ u =  \cos \theta x+  \sin \theta y, \quad v = -\sin \theta x + \cos \theta y$$

를 대입하면 타원을 2차식 형태

$$ A x^2 + 2B xy + C y^2 = 1$$

로 쓸 수 있는데, 그 계수는

$$ A= \frac{1}{8}\left( \frac{1}{\lambda_1} + \frac{1}{\lambda_2} +\Big( \frac{1}{\lambda_1} - \frac{1}{\lambda_2} \Big) \cos (2\theta) \right)=\frac{\tilde\mu_{02}}{ 4(\tilde\mu_{20}\tilde\mu_{02}-\tilde\mu_{11}^2 )} $$

$$ C= \frac{1}{8}\left( \frac{1}{\lambda_1} + \frac{1}{\lambda_2} -\Big( \frac{1}{\lambda_1} - \frac{1}{\lambda_2} \Big) \cos (2\theta) \right)=\frac{\tilde\mu_{20}}{ 4(\tilde\mu_{20}\tilde\mu_{02}-\tilde\mu_{11}^2 )} $$

$$ B= \frac{1}{8}\left( \frac{1}{\lambda_1} - \frac{1}{\lambda_2} \right)\sin (2\theta)=\frac{-\tilde\mu_{11}}{ 4(\tilde\mu_{20}\tilde\mu_{02}-\tilde\mu_{11}^2 )} $$

로 주어진다. 주어진 object의 normalized central moment를 구하면 타원 fitting에 대한 정보를 구할 수 있음을 보였다.

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

$\Sigma$는 영상에서 object pixel의 $x$ 뱡향 분산($\mu_{20}$: $x$축에 대한 회전관성), $y$ 방향 분산($\mu_{02}$: $y$축에 대한 회전관성), $x$-$y$의 correlation을 나타낸다. $\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}

(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, \quad \mu_{02}^\text{(e)} = \frac{\pi}{4} ab^3,\quad \mu_{11}^\text{(e)}= 0$$

으로 구해짐을 알 수 있다.

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

$$ I_\text{min} = \mu_{20}^\text{(e)}, \quad I_\text{max} =\mu_{02}^\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을 타원 피팅할 때 윤곽선 정보를 이용하는데, 이 방법은 윤곽선을 추출할 필요가 없어서 편리하다. 단 내부에 빈 곳이 있는 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);
}