RANSAC: Circle Fit

RANSAC 알고리즘을 써서 주어진 2차원 점집합에서 원을 추정한다. 원을 만들기 위해서는 최소한 3점이 필요하고, 또 일직선에 있지 않아야 한다. 이렇게 만들어진 원은 세 점을 꼭짓점으로 하는 삼각형의 외접원(circumcircle)이다. 주어진 외접원에서 크게 벗어나지 않는 inliers를 찾는다(추가로 이 inliers에 대해 최소자승법으로 원의 중심과 반지름을 다시 구해서 보다 정밀하게 추정하는 과정을 넣을 수도 있다). 무작위로 선택된 세 점에 대해 위의 과정을 반복 시행해서 구한 원 중에서 가장 많은 inliers를 포함하는 원을 결과로 사용한다.

// 참고: http://en.wikipedia.org/wiki/RANSAC

red = inliers

// 2024.5.31 재작성;
double dist_deviate(const CfPt& pts, double cparam[3]) {
    double dx = pts.x - cparam[0];
    double dy = pts.y - cparam[1];
    return fabs(hypot(dx, dy) - cparam[2]);
}
int circumcircle(CfPt pts[3], double cparam[3]) {
    double x1 = pts[0].x, x2 = pts[1].x, x3 = pts[2].x;
    double y1 = pts[0].y, y2 = pts[1].y, y3 = pts[2].y;
    double bax = x2 - x1, bay = y2 - y1;
    double cax = x3 - x1, cay = y3 - y1;
    double E = bax * (x1 + x2) + bay * (y1 + y2);
    double F = cax * (x1 + x3) + cay * (y1 + y3);
    double G = 2. * (bax * (y3 - y2) - bay * (x3 - x2));
    if (G == 0.) return 0;    //error;
    //assert(fabs(G)>small_epsilon); //to prevent collinear or degenerate case;
    cparam[0] = (cay * E - bay * F) / G; //cx;
    cparam[1] = (bax * F - cax * E) / G; //cy;
    cparam[2] = hypot(cparam[0]-x1, cparam[1]-y1); //rad;
    return 1;
};
int num_sampling3(double prob_fail, double inlier_ratio) {
    return int(log(prob_fail)/log(1-pow(inlier_ratio, 3))); 
}
std::vector<int> Ransac_CircleFit(std::vector<CfPt>& points, double circle_param[3]) {
    if (points.size() < 3)
        return std::vector<int> (); //return null_vector;

    CfPt center; double inv_scale;
    // normalize input points for the sake of numerical stability;
    std::vector<CfPt> nor_pts = normalize(points, inv_scale, center);
    // distance threshold;
    double distance_thresh = sqrt(double(points.size())) * inv_scale;
    //ransac
    int sample_num = 1000;	//number of sample
    int ransac_count = 0;
    const double prob_fail = 0.01;
    double best_cparam[3] = {0};
    std::vector<int> best_inliers;
    while (sample_num > ransac_count) {
        // pick random 3 indices:[0,points.size()-1];
        int triple[3];
        random_triple(points.size()-1, triple);
        CfPt selected[3];
        for (int i = 0; i < 3; i++) 
            selected[i] = nor_pts[triple[i]];
        // circumcircle of 3 points;
        if (circumcircle(selected, circle_param)) {
            // find inliers;
            std::vector<int> inliers;
            inliers.reserve(points.size());
            for (int i = nor_pts.size(); i-->0;) {
                // error measure = algebric distance;
                double deviation = dist_deviate(nor_pts[i], circle_param);
                if (fabs(deviation) < distance_thresh)
                    inliers.push_back(i);
            }
            if (inliers.size() > best_inliers.size()) {			
                // update sampling_num;
                sample_num = num_sampling3(prob_fail, double(inliers.size())/points.size());
                // update best_inliers;
                best_inliers.swap(inliers);
                // update best circle param;
                for (int i = 0; i < 3; i++) 
                    best_cparam[i] = circle_param[i];
            }
        }
        if (++ransac_count > 1500) {
            TRACE("error! ransac_count exceed!\n");
            break;
        }
    }
    // recover original coordinate and scale;
    denormalize(best_cparam, best_cparam, inv_scale, center);
    if (best_cparam[0] > 0 && best_cparam[1] > 0) {
        for (int i = 0; i < 3; i++)
            circle_param[i] = best_cparam[i];
        TRACE("circle_found(%d, %d)\n", sample_num, ransac_count);
        // more accurate estimation needed at this stage;
    } else 
        best_inliers.clear();
    return best_inliers;
}

https://kipl.tistory.com/207

Least Squares Fitting of Circles

https://kipl.tistory.com/357

Circle Fitting: Pratt