다섯개의 점들 (x1,y1), (x2,y2), (x3,y3), (x4,y4), (x5,y5)에 대해서 conic eq의 왼쪽값(대수적인 거리^2)이 최소가 되도록 하는 conic parameter를 결정하는 식은 다음과 같다: 

이 식은 해는 다음의 방정식을 만족시키는 해를 구하는 것과 같다:

또는 행렬식으로 표현하면,

의 형태로 쓰여진다. 6개의 conic-parameters를 결정해야 하므로, A는 마지막 행을 0으로 채워서 6x6 행렬로 만들어서 사용한다.  따라서 A(->A6x6)는 singular 행렬이 된다. A를 다음과 같이 singular value decomposition(SVD)를 할 때

가장 작은 eigenvalue(w)에 해당하는 V의 열벡터가 에러를 최소로 하는 conic parameter를 준다 (행렬 V의 각 열들이 6차원 벡터공간의 기저를 형성한다. U는 orthogonal 행렬이어서 오른쪽에서 곱해지는 벡터의 크기를 변화시키지 않고 단지 회전만 시킨다. 따라서 가장 작은 고유값에 해당하는 V의 열벡터가 |A.x|를  최소화시킨다)



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

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

//

static int ransac_sample_number(int ndata, int ngood, double prob, int 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)

{

    const int npts = points.size(); 

    if (npts < 5) {

        return 0;

    }

    // 안정적인 conic fitting을 위해서 데이터를 (-1,1)구간으로 normalize함;

    std::vector<CfPt> nor_pts(npts);

    double ellipse_scale; //원점에서  타원점들까지의 평균거리;

    CfPt ellipse_center;  //타원의 중심 예측;

    normalize_points(points, nor_pts, &ellipse_scale, &ellipse_center);

    //

    std::vector<int> inliers_index(npts);

    best_index.resize(npts);

    // 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((npts - 1), rand_idx);

  //5 점을 이용해서 conic parameter estimate;

double conic_params[6];

get_conic_params(nor_pts, rand_idx, conic_params);


        //conic parameter를 이용해서 inlier후보들을 모음;

        int ninliers = 0;

        for (int i = 0; i < npts; i++) {

            double error = fitting_error(conic_params, nor_pts[i].x, nor_pts[i].y);

            if (fabs(error) < dist_threshold) {

                inliers_index[ninliers] = i;

                ninliers++;

            }

        }


        if (ninliers > 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 = ninliers;

                    //이 단계에서 inliers를 이용해서 다시 타원을 예측한다;

                    //if (ninliers >= 6) {

                    // std::vector<double> pts(2 * ninliers);

                    // for (int i = 0; i < ninliers; i++) {

                    // pts[2 * i + 0] = points[inliers_index[i]].x ;

                    // pts[2 * i + 1] = points[inliers_index[i]].y ;

                    // }

                    // EllipseFit(&pts[0], ninliers, conic_params);

                    // for (int i = 0; i < 5; i++) best_params[i] = conic_params[i];

                    //}

                    // if (ninliers == 5) 

                    for (int i = 0; i < 5; i++) best_params[i] = ellipse_params[i];

                    for (int i = 0; i < ninliers; i++) best_index[i] = inliers_index[i];

                    sample_num = ransac_sample_number(npts, ninliers, 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();

}

Posted by helloktk