타원은 원뿔을 평면으로 잘랐을 때 생기는 곡선 중의 하나로 다음 이차 형식(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개의 점을 선택해야 한다.

 

// 주어진 점에서 타원까지의 대수적 거리: conic-eq의 좌변의 값;

static double fitting_error(double conic_params[6], double x, double y);

더보기
static double fitting_error(CfPt q, double conic_params[6]) {
     double x = q.x, y = q.y;
     return conic_params[0] * x * x +
            conic_params[1] * x * y +
            conic_params[2] * y * y +
            conic_params[3] * x + 
            conic_params[4] * y + 
            conic_params[5];
}

//

static int ransac_sample_number(int ndata, int ngood, double prob, int ntry);

더보기
static int ransac_sample_number(int ndata, int ngood, double prob, int ntry) {     
      if (ndata == ngood) return 1;
      return (int)(log(1. - prob)/log(1. - pow((double)ngood/ndata, ntry)));
}
int ransac_ellipse_fitter(std::vector<CPoint>& points, int width, int height,
                          double ellipse_params[5],/*(a, b, cx, cy, theta)*/
                          std::vector<int> &best_index) {
    if (points.size() < 5) return 0;
    // 안정적인 conic fitting을 위해서 데이터를 (-1,1)구간으로 normalize함;
    std::vector<CfPt> nor_pts(points.size());
    double ellipse_scale; //원점에서  타원점들까지의 평균거리;
    CfPt ellipse_center;  //타원의 중심 예측;
    normalize_points(points, nor_pts, &ellipse_scale, &ellipse_center);
    //
    best_index.resize(nor_pts.size());
    // inverse_scale을 기준으로 만들어졌음;
    const double dist_threshold = sqrt(3.84) / ellipse_scale;
    double best_params[5];
#define MAX_ITERATION (1500)
    int sample_num = 1000; //number of sample
    int ransac_iter = 0;
    int max_inliers = 0;
    while (sample_num > ransac_iter) {
        int rand_idx[5];
        // 5개 점을 랜덤하게 선택;
        random_quintet((nor_pts.size() - 1), rand_idx);
        // 5개 점을 이용해서 conic parameter estimate;
        double conic_params[6];
        get_conic_params(nor_pts, rand_idx, conic_params);
        // conic parameter를 이용해서 inlier 후보들을 모음;
        std::vector<int> inliers; // inlier index을 저장;
        for (int i = nor_pts.size(); i-->0;) {
            double error = fitting_error(nor_pts[i], conic_params);
            if (fabs(error) < dist_threshold) inliers.push_back(i);
        }
        if (inliers.size() > max_inliers) {
            if (conic_to_ellipse(conic_params, ellipse_params)) {
                get_real_params(ellipse_params, ellipse_scale, ellipse_center, ellipse_params);
                // 주어진 영역안의 정상적인 타원만 선택함;
                if (regular_ellipse(ellipse_params, width, height, 0.5, 2.0)) {
                    max_inliers = inliers.size();
                    //이 단계에서 inliers를 이용해서 다시 타원을 예측한다;
                    //if (inliers.size() >= 6) { 
                    //   std::vector<CfPt> pts(inliers.size());
                    //   for (int i = inliers.size(); i-->0;) pts[i] = nor_pts[inliers[i]];
                    //   EllipseFit(pts, conic_params);
                    //   conic_to_ellipse(conic_params, ellipse_params);
                    //   get_real_params(ellipse_params, ellipse_scale, ellipse_center, ellipse_params);
                    //}
                    for (int i = 0; i < 5; i++) best_params[i] = ellipse_params[i];
                    for (int i = inliers.size(); i-->0;) best_index[i] = inliers[i];
                    sample_num = ransac_sample_number(nor_pts.size(), inliers.size(), 0.99, 5);
                }
            }
        }
        ransac_iter++;
        if (ransac_iter > MAX_ITERATION) {
            TRACE("Error! ransac_iter exceeds!\n");
            max_inliers = 0;
            break ;
        }
    }
    if (best_params[0] > 0 && best_params[1] > 0)
        for (int i = 0; i < 5; i++) ellipse_params[i] = best_params[i];
    else
        max_inliers = 0;
    best_index.resize(max_inliers); 
    return best_index.size();
}

 
 
 
 
 
728x90
Posted by helloktk
,