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

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

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

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

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

\begin{gather} \begin{pmatrix} a & b/2 \\ b/2 & c \end{pmatrix} \Longrightarrow \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix}  \quad \text{det} = ac - b^2/4\end{gather}

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

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

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

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

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

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

로 주어짐을 확인할 수 있다. 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$$

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

$$ ax_0^2 + b x_0 y_0^2 + cy_0^2 + dx_0 +e y_0 +f = f - (ax_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}} , \quad 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;
};
728x90
,

타원은 원뿔을 평면으로 잘랐을 때 생기는 곡선 중의 하나로 다음 이차 형식(quadratic form)으로 표현된다.

$$\text{conic eq: }   ax^2 + bxy + cy^2 +d x+ ey + f=0$$

이차 형식의 계수를 구하기 위해서는 타원 위의 서로 다른 5개의 점이 필요하다 (타원이기 위해서는 $b^2 - 4a c < 0$). 주어진 5개의 평면상의 점 $\{(x_i, y_i)|i=1,2,..,5\}$가 결정하는 타원은

\begin{gather}  ax_1^2 + bx_1 y_1 + c y_1^2 + d x_1 +e y_1 + f =0 \\  ax_2^2 + bx_2 y_2 + c y_2^2 + d x_2 +e y_2 + f =0 \\  ax_3^2 + bx_3 y_3 + c y_3^2 + d x_3 +e y_3 + f =0 \\ ax_4^2 + bx_4 y_4 + c y_4^2 + d x_4 +e y_4 + f =0 \\ ax_5^2 + bx_5 y_5 + c y_5^2 + d x_5 +e y_5 + f =0 \end{gather}

이를 만족시키는 해는 다음 행렬식으로 쓰인 방정식의 해와 같다:

$$ \begin{pmatrix}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_3^ 2& x_3 y_3 & y_3 ^2 & x_3 & y_3 & 1\\ x_4^ 2& x_4 y_4 & y_4 ^2 & x_4 & y_4 & 1\\ x_5^ 2& x_5 y_5 & y_5 ^2 & x_5 & y_5 & 1   \end{pmatrix} \begin{pmatrix} a\\b\\c\\d\\e\\f  \end{pmatrix}=0 $$

또는 간단히

$$ \mathbf{A}.\mathbf{x} =\mathbf{0}$$

의 형태로 쓰여진다. 6개의 conic-parameters를 결정해야 하므로, $\mathbf{A}$는 마지막 행을 0으로 채워서 $6\times6$ 행렬로 만들어서 사용한다.  이 경우 $\mathbf{A}_{6\times 6}$는 determinant가 0인 singular 행렬이 되고 적어도 null 공간의 차원은 1 이상이다(1인 경우만 의미있는 해를 준다). $\mathbf{A}$의 singular value decomposition(SVD)을 이용하면위 방정식의 nontrivial 한 해를 얻을 수 있다. SVD에 의해 $\mathbf{A}$는 다음처럼 쓸 수 있다( singular value 중 1개가 0인 경우).

$$\mathbf{A}= \mathbf{U}. \text{diag}(w_0,w_1, w_2, w_3, w_4, 0).\mathbf{V}^T= \mathbf{U}.\begin{pmatrix}  w_0 \mathbf{v}_0^T \\ w_1 \mathbf{v}_1^T\\  w_2 \mathbf{v}_2^T\\ w_3 \mathbf{v}_3^T\\ w_4 \mathbf{v}_4^T\\ 0\end{pmatrix}$$

행렬 $\mathbf{V}$의 각 열벡터 $\{\mathbf{v}_i\}$는 6차원 벡터 공간의 기저를 형성한다. 그리고 $\mathbf{U}$는 orthogonal 행렬이므로 오른쪽에서 곱해지는 벡터의 크기를 변화시키지 않고 단지 회전만 시킨다. singular value가 $0$에 해당하는 $\mathbf{V}$의 열벡터 $\mathbf{v}_5$를 대입하면

$$ \mathbf{A} \cdot \mathbf{v}_5    =   \mathbf{U}. (0,0,0,0,0)^T  = 0$$

이어서 $\mathbf{A}.\mathbf{x}=0$의 nontrivial solution이 됨을 알 수 있다. 물론 singular value가 0인 경우가 여럿이 생기면 해당하는 열벡터의 선형결합도 해가 되므로 원하는 답이 아니다. 이 경우는 주어진 5개의 점들의 일부가 일직선상에 있는 경우에 해당한다(또는 겹치는 경우). 이때는 다시 5개의 점을 선택해야 한다.

red = inliers

// 2024.05.30에 새롭게 작성;
double alg_distance(const CfPt& pts, double conic_param[6]) {
    return conic_param[0] * pts.x * pts.x + 
           conic_param[1] * pts.x * pts.y +
           conic_param[2] * pts.y * pts.y + 
           conic_param[3] * pts.x + 
           conic_param[4] * pts.y + 
           conic_param[5];
}
int num_sampling5(double prob_fail, double inlier_ratio) {
    return int(log(prob_fail)/log(1-pow(inlier_ratio, 5))); 
}
std::vector<int> Ransac_EllipseFit(std::vector<CfPt>& points, double ellipse_param[5]) {
    if (points.size() < 5) return std::vector<int> (); //return null_vector;

    CfPt center; double inv_scale;
    // normalize input points for the sake of numerical stability;
    std::vector<CfPt> nor_pts = normalize(points, inv_scale, center);
    // RANSAC;
    const double prob_fail = 0.01;
    int sample_num = 1000;	//number of sample
    int ransac_count = 0;
    const double distance_thresh = sqrt(3.84) * inv_scale;
    double best_eparam[5] = {0};
    std::vector<int> best_inliers;
    while (sample_num > ransac_count) {
        // pick random 5 indices:[0,points.size()-1];
        int quintet[5];
        random_quintet(points.size()-1, quintet);
        CfPt selected[5];
        for (int i = 0; i < 5; i++) 
            selected[i] = nor_pts[quintet[i]];
        // ellipse parameters with 5 points;
        double conic_param[6];
        solve_conic(selected, conic_param);
        // find inliers;
        std::vector<int> inliers;
        inliers.reserve(points.size());
        for (int i = nor_pts.size(); i-->0;) {
            // error measure = algebric distance;
            double deviation = alg_distance(nor_pts[i], conic_param);
            if (fabs(deviation) < distance_thresh)
                inliers.push_back(i);
        }
        if ((inliers.size() > best_inliers.size()) && 
             solve_ellipse(conic_param, ellipse_param)) {
            // update sampling_num;
            sample_num = num_sampling5(prob_fail, double(inliers.size())/points.size());
            // update best_inliers;
            best_inliers.swap(inliers);
            // update best ellipse param;
            for (int i = 0; i < 5; i++) 
                best_eparam[i] = ellipse_param[i];    
        }
        if (++ransac_count > 1500) {
            TRACE("error! ransac_count exceed!\n");
            break;
        }
    }
    // recover original scale and coordinate;
    denormalize(best_eparam, best_eparam, inv_scale, center);
    if (best_eparam[0] > 0 && best_eparam[1] > 0) {
        for (int i = 0; i < 5; i++)
            ellipse_param[i] = best_eparam[i];
        TRACE("ellipse found(%d, %d)\n", sample_num, ransac_count);
    } else 
        best_inliers.clear();
    return best_inliers;
}
728x90
,

1. Derivative Based Algorithms:

  • Thresholded absolute gradient: $$ f  =\sum_x \sum _y |I (x+1, y) - I(x, y)| $$  where   $|I(x+1, y) - I(x, y)| > \text{threshold};$
  • Squared gradient: $$f=\sum_x \sum_y ( I(x+1, y)- I(x, y)) ^2 $$ where $(I(x+1, y) - I(x, y))^2 > \text{threshold};$
  • Brenner gradient: $$f=\sum_x \sum_y ( I(x+2, y)- I(x, y)) ^2 $$ where $(I(x+2, y) - I(x, y))^2 > \text{threshold};$
  • Tenenbaum gradient: $$f = \sum_x \sum_y (\text{SobelX}^2 (x, y)+ \text{SobelY}^2(x,y));$$
  • Energy Laplace: $$f = \sum_x \sum_y ( (L * I) (x, y))^2$$  where $$L =\left[\begin{array}{ccc} -1& -4 & -1 \\  -4 & 20 & -4 \\ -1 & -4 & -1 \end{array}  \right]$$ 
  • Sum of modified Laplace: $$f = \sum_x \sum_y (\text{LaplaceX}(x, y)| + |\text{LaplaceY}(x, y)|) $$
  • Sum of squared Gaussian derivatives $$f = \frac{1}{\text{total pixels}} \sum_x \sum_y ((\text{Gaussian derivative X}(x,y)) ^2 + (\text{Gaussian derivative Y}(x,y))^2 )$$ where $\sigma  = d/(2\sqrt{3})$, $d$ = dimension of the smallest feature;

 

2. Statistical Algorithms:

  •  Variance Measure ; $$f = \frac{1}{\text{total pixels}}\sum _x \sum _y ( I  (x, y) - \overline{I})^2;$$ where $\overline{I}$ is the mean of image.
  • Normalized Variance Measure ; $$f = \frac{1}{\text{total pixels}\times \overline{I}} \sum_x \sum_y ( I(x, y)- \overline{I})^2;$$
  • Auto-correlation Measure: $$f = \sum_x\sum_y I(x, y) I(x+1, y) - \sum _x \sum_y I(x, y) I(x+2, y);$$
  • Standard Deviation-based Correlation: $$f = \sum_x \sum_y I(x, y) I(x+1, y) - \overline{I}^2(x, y)\times \text{total pixels};$$

3. Histogram-based Algorithms:

  • Range Algorithm: $$f =\text{max}\{i| \text{histogram}(i)>0\} - \text{min}\{ i| \text{histogram}(i)>0\};$$
  • Entropy Algorithm: $$f=-\sum_{i=0}^{255} p(i) \log_{2} p(i)$$ where $p(i)= \text{histogram}(i)/\text{total pixels};$

4. Other Algorithms:

  • Threshold Contents: $$f=\sum_x \sum_y I(x, y)$$ where $I(x,y) \ge \text{threshold};$
  • Thresholded Pixel Count: $$f=\sum_x\sum_y \theta(\text{threshold}- I(x, y));$$ where $\theta(x)$ is the step-function.
  • Image Power: $$f= \sum_{x}\sum_{y} I^2(x, y)$$ where $I(x, y) \ge\text{threshold};$

Ref: Dynamic evaluation of autofocusing for automated microscopic analysis of blood smear and pap smear, J. Microscopy,227, 15(2007).                           

 

728x90

'Image Recognition' 카테고리의 다른 글

Union-Find Connected Component Labeling  (0) 2012.11.01
RANSAC: Ellipse Fitting  (1) 2012.10.07
Statistical Region Merging  (2) 2012.03.25
Local Histogram Equalization  (0) 2012.03.10
2차원 Savitzky-Golay Filters 응용  (0) 2012.02.28
,