728x90

RANSAC 알고리즘을 써서 주어진 2차원 점집합에서 원을 추정한다. 원을 만들기 위해서는 최소한 3점이 필요하고, 또 일직선에 있지 않아야 한다. 이렇게 만들어진 원은 세 점을 꼭짓점으로 하는 삼각형의 외접원이다. 주어진 외접원에서 크게 벗어나지 않는 inliers를 찾고, 이 inliers에 대해 최소자승법으로 원의 중심과 반지름을 다시 추정한다. 무작위로 선택된 세 점에 대해 위의 과정을 반복 시행해서 구한 원 중에서 inliers의 벗어남 편차가 최소인 것을 결과로 선택한다.
// 참고: 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;(참조: kipl.tistory.com/207)
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 findInlier(double xp[], double yp[], int N, 
                double cx, double cy, double rad,
                double dist_th,
                double consensusx[], double consensusy[], double *var) {
    int ninlier = 0;
    double err = 0;
    *var = 0;           // variance of distance deviation;
    for (int k = 0; k < N; ++k) {
        double dist = sqrt(SQR(xp[k] - cx) + SQR(yp[k] - cy));
        double distdeviate = fabs(dist - rad) / rad ;
        if (distdeviate <= dist_th) {  //collect maybe_inliers;
            consensusx[ninlier] = xp[k];
            consensusy[ninlier] = yp[k] ;
            *var += SQR(dist - rad);
            ninlier++ ;
        }
    }
    return ninlier;
}
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));
    int iter = min(max_iter, trials99);
    int found = 0;   
    //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) { 
        int tri[3]; 
        double tx[3], ty[3];
        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
EM: Binarization  (0) 2008.07.01
Posted by helloktk

댓글을 달아 주세요