Geometry & Recognition

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

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

으로 정의한다.

std::vector<BYTE> color_equalize_new(std::vector<BYTE> &rgb) {
    int hist[256] = {0};
    int chist[256], lut[256], partition[256 + 1];
    std::vector<BYTE> out = raster; // clone; 
    // pdf of gray level: g = (r + g + b) / 3
    for (int k = 0; k < rgb.size(); k += 3) {
        int a = rgb[k+0]; a += rgb[k+1]; a += rgb[k+2];
        ++hist[a/3];
    };
    // cdf;
    for (int i = 0, s = 0; i < 256; i++) {
        s += hist[i]; chist[i] = s;
    }
    
    /* assign equal number of pixels in each gray levels;*/
    int num_pixels = rgb.size()/3;
    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++;
        /* chist[j] <= desired < chist[j+1];*/
        /* nearest interpolation */
        if ((desired - chist[j]) < (chist[j+1] - desired)) partition[i] = j;
        else partition[i] = j + 1;
    } 
    /* Find i s.t. partion[i] <= j < partition[i+1];*/
    for (int j = 0; j < 256; j++) {
        int i = 0;
        while (partition[i+1] <= j) i++;
        lut[j] = i;
    } 
    // needs hue preserving processes;
    for (k = rgb.size(); k-->0;)
        res[k] = lut[src[k]];
    return out;
}

std::vector<BYTE> color_equalize_HSV(std::vector<BYTE> &rgb) {
    int hist[256] = {0};
    int chist[256] = {0}, lut[256] = {0}, partition[256 + 1];
    std::vector<BYTE> out = raster; // cloning;
    std::vector<double> fH(rgb.size()/3);
    std::vector<double> fV(rgb.size()/3);
    std::vector<double> fS(rgb.size()/3);

    const int n = rgb.size()/3;
    double r, g, b;
    for (int k = 0; k < rgb.size(); k += 3) {
        b = rgb[k+0], g = rgb[k+1], r = rgb[k+2];
        RGBtoHSV(r / 255, g / 255, b / 255, fH[k], fS[k], fV[k]);
    };	
    // make histogram of V-color;
    for (int k = n; k-->0;)
        ++hist[int(fV[k] * 255)];

    // cdf;
    for (int i = 0, s = 0; i < 256; i++) {
        s += hist[i]; chist[i] = s;
    }
    /* assign equal number of pixels in each V-color;*/
    int num_pixels = rgb.size()/3;
    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 = n; k-->0;)
        fV[k]= lut[int(fV[k] * 255)] / 255.;

    for (int k = 0; k < rgb.size(); k += 3) {
        HSVtoRGB(fH[k], fS[k], fV[k], r, g, b);
        out[k+0] = int(b * 255); 
        out[k+1] = int(g * 255);
        out[k+2] = int(r * 255);
    }
    return out;
}

/* 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)=a{x}^2+bxy+c{y}^2+dx+ey+f=0$$

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

$$ϵ=\sum _i|F(x_i,y_i){|}^2\to \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{array}{c}ϵ={\mathbf{a}}^T.\mathbf{S}.\mathbf{a}-\lambda ({\mathbf{a}}^T.\mathbf{C}.\mathbf{a}-1)\end{array}$$

여기서 $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\\ & & ⋮& & \\ 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].$$

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

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

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

$$ϵ={\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$을 만족하도록 크기를 조정하면 $\stackrel{\mathbf{~}}{\mathbf{a}}=\mathbf{a}/\sqrt{{\mathbf{a}}^T.\mathbf{C}.\mathbf{a}}$로 주어진다. 

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

$$\mathbf{Q}\mathbf{a}\mathbf{=}\lambda {\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{C}\mathbf{a}\mathbf{=}\lambda {\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{C}{\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{Q}\mathbf{a}$$ 이므로 더 다루기 쉬운 ${\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{C}{\mathbf{Q}}^{\mathbf{-}\mathbf{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\ne 0$인 경우가 원을 표현한다. 물론 $A=1$로 설정을 할 수 있으나 이 경우는 주어진 데이터가 원의 일부 부분만 나타내는 경우 문제가 생긴다. 원은 중심과 반지름을 알면 결정되므로 3개의 변수가 필요한데 위의 이차 다항식은 4개의 변수를 가진다. 원을 제대로 표현하기 위해서는 제한 조건이 들어오는데 

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

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

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

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

$$ϵ(A,B,C,D)=\sum |A(x_i^2+y_i^2)+B{x}_i+C{y}_i+D{|}^2-\lambda (B^2+C^2-4AD-1)⟶\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,$$

$$\mathtt{x}=\left(\begin{array}{c}A\\ B\\ C\\ D\end{array}\right),\mathtt{M}\mathtt{=}\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)\mathtt{,}\mathtt{N}\mathtt{=}\left(\begin{array}{cccc}0& 0& 0& -2\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ -2& 0& 0& 0\end{array}\right)$$

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

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

$ϵ$을 ${\mathtt{x}}^{\mathtt{T}}$에 대해서 미분하면

$$\mathtt{M}\mathtt{.}\mathtt{x}\mathtt{=}\lambda \mathtt{N}\mathtt{.}\mathtt{x}$$

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

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

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

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

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

그런데 주어진 고윳값 $\lambda$ 해당하는 고유 벡터가 $\mathtt{x}$일 때, 임의의 0이 아닌 상수 $\mu$에 대해 $\mu \mathtt{x}$도 역시 고유벡터이다. 고유 벡터의 크기를 고정하기 위해서 제한조건을 적용하면 $1={\mu }^2{\mathtt{x}}^{\mathtt{T}}\mathtt{.}\mathtt{N}\mathtt{.}\mathtt{x}\mathtt{=}{\mu }^{\mathtt{2}}{\mathtt{x}}^{\mathtt{T}}\mathtt{.}\mathtt{M}\mathtt{.}\mathtt{x}\mathtt{/}\lambda$이므로 $\mu =\sqrt{\lambda /{\mathtt{x}}^{\mathtt{T}}\mathtt{.}\mathtt{M}\mathtt{.}\mathtt{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,\mathrm{\Delta }\equiv {\lambda }^2-M_z\lambda +C_{xy}$$

$$A=1$$

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

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

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

$$c_x=-\frac{B}{2A}=\frac{M_{zx}(M_{yy}-\lambda )-M_{zy}{M}_{xy}}{2\mathrm{\Delta }}$$

$$c_y=-\frac{C}{2A}=\frac{M_{zy}(M_{xx}-\lambda )-M_{zx}{M}_{xy}}{2\mathrm{\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$의 부호 개수는  $\mathtt{N}$의 고유값 부호별 개수와 같아야 하는데, $\mathtt{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);
}

https://kipl.tistory.com/32

https://kipl.tistory.com/207