728x90

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

.

이 이차 형식의 계수를 구하기 위해서는 타원 위의 5개의 서로 다른 점이 필요하다 $(a c > 0)$. 그러나 주어진 점 정보를 얻는 과정에서 노이즈가 포함되어 conic eq에 대입했을 때 우변이 0이 안될 수 있다. 이런 경우를 고려해서, 5개 점 $\{(x_i,y_i)|i=1,...,5\}$에 대해서 conic eq의 좌변 (대수적인 거리^2) 최소가 되도록 하는 conic parameter 구하도록 하자:

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

또는 간단히

의 형태로 쓰여진다. 6개의 conic-parameters를 결정해야 하므로, A는 마지막 행을 0으로 채워서 6x6 행렬로 만들어서 사용한다.  이 경우 $\bf A$(->6x6)는 determinatn가 0인 singular 행렬이 된다. 근사해를 구하기 위해 $\bf A$를 singular value decomposition(SVD)해서 아래처럼 쓰면

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

 

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

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

더보기
static double fitting_error(double conic_params[6], double x, double 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) {
    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

댓글을 달아 주세요

  1. 오현경 2019.05.05 22:59  댓글주소  수정/삭제  댓글쓰기

    죄송하지만 전체 소스를 볼 수 있을까요??

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

댓글을 달아 주세요

RANSAC Algorithm

Image Recognition 2008. 5. 24. 09:02
728x90

어떤 현상을 설명하는 이론 모델을 만들려면 현상에서 관측 데이터를 모아야 한다. 그런데 관측 데이터에는 모델에 대한 잘못된 가정이나 측정의 오차 등에서 생기는 여러 가지 형태의 노이즈가 들어 있게 마련이다. 이러한 노이즈는 모델을 정확히 예측하는 것을 방해한다. 모델 파라미터의 예측을 방해하는 데이터(outliners)가 들어 있는 관측 데이터로부터 어떻게 하면 모델을 구축할 수 있는가라는 질문에 대한 한 해결책을 RANSAC("RANdom SAmple Consensus") 알고리즘이 제공한다. 이 알고리즘은 1981년 Fischler과 Bolles에 의해서 제안되었다.

RANSAC 알고리즘에 들어가는 입력은 관측된 데이터, 관측 결과를 피팅하거나 설명할 수 있는 모델 파라미터, 그리고 신뢰성을 나타내는 파라미터를 필요로 한다.

RANSAC 알고리즘은 주어진 원본 데이터에서 일부를 임의로 선택하는 과정을 반복하여서 최적의 파라미터를 예측하는 프러시저의 형태를 가진다. 전체의 관측 데이터가 M개 있고, 모델 파라미터를 예측하는데 N개의 데이터가 필요한 경우, 알고리즘의 동작은

     1. 임의로 관측 데이터에서 N개의 부분 데이터를 선택한다.
     2. 선택 데이터를 (가상의) inlier로 생각하고 모델을 예측한다.
     3. 원본 데이터(M) 중에서 예측된 모델에 잘 맞는가를 체크한다. 잘 맞는 데이터 수를 K라고 한다.
     4. 만약에 K가 충분하면, 이 모델을 확정하고 알고리즘을 종료한다.
     5. 모델이 좋지 않으면 14 과정을  L번 반복한다.
     6. 충분한 반복 후에도 좋은 모델을 찾지 못하면 예측에 실패한 것으로 간주한다.

1. 임계값 K는 주어진 전체 데이터 중에서 어느 정도의 비율로 찾는 모델에 잘 맞는가로 판단하는 기준인데, 사용자가 결정해야 한다. 대략적으로 주어진 샘플에서 inlier의 비율을 Pg정도라고 생각되면 다음 정도로 잡으면 된다:

 
   K = M * (1- P
g)

2. 얼마나 많은 반복을 해야하는가? 주어진 관측 데이터에서 inlier일 확률이 Pg인 경우에 L번의 모델 예측 시도가 실패할 확률을 계산하여서 이것이 주어진 설정값, pfail 보다도 작은 경우에 모델 예측의 실패로 처리하므로,

  pfail   = L번의 모델예측 실패 확률
          = (한 번 시도가 실패할 확률)L
          = (1 - 한 번 시도가 성공할 확률))L
          = (1 - (랜덤 데이터가 모델을 예측할 확률)N))L
          = (1- (Pg) N)L

이 사실로부터 최대 반복 횟수는 
    
   L = log(pfail)/log(1-(Pg)N)

로 주어진다.
inlier가 반이 섞인 경우 Pg (=주어진 데이터중에서 inlier일 확률) = 0.5,
pfail = 0.01 인 경우, N = 3 (세 점만 있으면 모델 구성이 가능한 원의 피팅이 한 예)이면 최대 반복 횟수는 윗 식에 적용하면,
     
                   L = 35 회

RANSAC 알고리즘은 주어진 관측 데이터에 많은 outlier가 존재하더라도 매우 정확히 모델 파라미터를 예측할 수 있는 능력이 있는 robust 한 알고리즘이지만, 얼마나 많은 반복을 해야 하는가에 대한 상한 값이 제공되지 않는다(사용자가 pfail 를 설정한다). 따라서 사용자 설정에 의한 상한 값은 최적의 모델이 아닐 수 있다.

"The RANSAC procedure is opposite to that of conventional smoothing techniques: Rather than using as much of the data as possible to obtain an initial solution and then attempting to eliminate the invalid data points, RANSAC uses as small an initial data set as feasible and enlarges this set with consistent data when possible".

Martin A. Fischler and Robert C. Bolles (June 1981). "Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography". Comm. of the ACM 24: 381–395.

예제 코드:
ransac을 이용한 라인 피팅: http://blog.naver.com/helloktk/80029006029
ransac을 이용한 원 피팅: http://kipl.tistory.com/32

 

RANSAC: Circle Fit

RANSAC 알고리즘을 써서 주어진 2차원 점집합에서 원을 추정한다. 원을 만들기 위해서는 최소한 3점이 필요하고, 또 일직선에 있지 않아야 한다. 이렇게 만들어진 원은 세 점을 꼭짓점으로 하는

kipl.tistory.com

ransac을 이용한 타원 피팅: kipl.tistory.com/110

 

RANSAC Ellipse Fitting

타원은 원뿔을 평면으로 잘랐을 때 생기는 곡선 중의 하나로 다음 이차 형식(quadratic form)으로 표현된다. . 이 이차 형식의 계수를 구하기 위해서는 타원 위의 5개의 서로 다른 점이 필요하다 $(a

kipl.tistory.com

 

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

Gaussian Mixture Model  (2) 2008.06.07
Rasterizing Voronoi Diagram  (0) 2008.05.26
RANSAC Algorithm  (0) 2008.05.24
Contour Tracing  (0) 2008.05.22
Gausssian Scale Space  (0) 2008.05.22
Watershed Algorithm 적용의 예  (2) 2008.05.21
Posted by helloktk

댓글을 달아 주세요