RANSAC 알고리즘을 이용하여서 원을 추정한다. 원을 만들기 위해서는 최소한 3개의 일직선상에 있지 않은 점을 잡아야 하면, 이때 원은 세점의 외접원으로 주어진다. 이 원을 기본으로 해서 다시 원의 중심과 반지름을 추정한다. 추정된 원들 중에서 inlier 들의 벗어남의 편차가 가장 작은 것을 결과로 이용한다.
// 참고: http://en.wikipedia.org/wiki/RANSAC 
//

// [0,max_num) 사이에서 3개의 숫자를 무작위로 선택함;
void GetRandomTriplet(int max_num, int triplet[3]) {

// 세 점의 외접원을 구함;
int CircumCircle(double x1, double x2, double x3,
                 double y1, double y2, double y3,
                 double *cx, double *cy, double *rad) {

// 3x3 선형방정식을 푼다: A * x = b;
bool SolveLinearEQ3x3(double A[9], double bb[3], double x[3]) {

// x^2 + y^2 + A*x + B*y + C = 0;
// Least square fit for A, B, C;
int CircleFit_LS(int N, double xp[], double yp[],
                   double *cx, double *cy, double *rad) {  

// 데이터 중에서 주어진 원에서 거리 임계값 이내의 데이터만 골라냄;
int findInlier(double xp[], double yp[], int N,
                double cx, double cy, double rad,
                double dist_th,
                double consensusx[], double consensusy[], double *var) {

int RansacCircleFit(int N, double xp[], double yp[],
                    int sample_th,      //# of inliers; >= 50% of data(N), 66%;
                    double dist_th,     // 25% of radius;   |dist-rad|/rad< dist_th
                    int max_iter,
                    double *centerx,
                    double *centery,
                    double *radius)
{    
    double pr = double(sample_th) / double(N);
    double trials99 = log(1. - 0.99)/log(1. - pow(pr, 3));
    double tx[3], ty[3];
    int tri[3];
    int found = 0;  
    int iter = min(max_iter, trials99);
    //inlier buffer;
    std::vector<double> consensusx(N), consensusy(N) ;
    double min_dev = 1.e+10, var, sdev;
    if (sample_th < 3) sample_th = 3;
    while (iter) {
        GetRandomTriplet(N, tri);
        for (int i = 0; i < 3; i++) {
            tx[i] = xp[tri[i]];
            ty[i] = yp[tri[i]];
        }
        double cx, cy, rad;
        if (!CircumCircle(tx[0], ty[0], tx[1], ty[1], tx[2], ty[2], &cx, &cy, &rad)) {
            // if tree points are degenerate or on the sample line, discard them!
            continue ;
        }
        int ninlier = findInlier(xp, yp, N, cx, cy, rad, dist_th, &consensusx[0], &consensusy[0], &var);
        if (ninlier >= sample_th) {          
            // estimate model params using maybe-inliers;
            if (!CircleFit_LS(ninlier, &consensusx[0], &consensusy[0], &cx, &cy, &rad))
                continue; //least-square fitting fails;
            // collect real-inliers;
            ninlier = findInlier(xp, yp, N, cx, cy, rad, dist_th / 2, &consensusx[0], &consensusy[0], &var);
            if (ninlier < sample_th)
                continue;
            sdev = sqrt(var / ninlier);
            // estimate model params using inliers again;
            if (!CircleFit_LS(ninlier, &consensusx[0], &consensusy[0], &cx, &cy, &rad))
                continue;            
            TRACE("cx = %f, cy = %f, rad=%f\n", cx, cy, rad);
            if (sdev < min_dev) {
                *centerx = cx;
                *centery = cy;
                *radius  = rad;
                min_dev = sdev;
                found = 1;
                // we need more elaborate termination conditions;
#if _DEBUG
#endif
            }
        }
        --iter;
    }
    return found;
};

  • sample_th = 2 * N / 3;
  • dist_deviate_th = 0.25;
  • 파란색: 최소자승법을 이용한 피팅 (outlier의 영향을 많이 받음);
  • 붉은색: RANSAC 피팅 결과 (2017.01.04 수정)


see also http://blog.naver.com/helloktk/80029898273


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

Chamfer Match  (0) 2008.08.01
Retinex Algorithm  (2) 2008.07.26
RANSAC : Circle Fit  (0) 2008.07.21
KMeans Algorithm  (0) 2008.07.19
Robust Line Fitting  (0) 2008.07.08
Bayesian Spam Filtering  (0) 2008.07.03
Posted by helloktk

댓글을 달아 주세요