C++ 解法, 执行用时: 5ms, 内存消耗: 376KB, 提交时间: 2020-11-10

#include<bits/stdc++.h>
using namespace std;

#define real float
struct Node{
    real x, y;
};

vector<real> ransac(vector<Node>& arr,
                   real outlier_prob=0.1,
                   real accept_prob=1e-3,
                   real threshold=10.0){
    default_random_engine generator;
    int n_sample = arr.size();
    int n = 2;   // 拟合模型需要的最小数据量
    // 计算理论最大迭代次数
    real inlier_prob = 1 - outlier_prob;
    real sample_fail = 1 - inlier_prob*inlier_prob;
    int k = log(accept_prob) / log(sample_fail);
    
    real a_res, b_res, c_res;
    real min_error = numeric_limits<real>::max();
    while(k--){
        // 随机采样 n 个样本
        uniform_int_distribution<int> sampler(0, n_sample-1);
        int idx1=0, idx2=0;
        while(idx1==idx2){
            idx1 = sampler(generator);
            idx2 = sampler(generator);
        }
        Node p1=arr[idx1], p2=arr[idx2];
        // 拟合模型：a*x1+b*y1+c=0 和 a*x2+b*y2+c=0
        // 两式相减：a*(x1-x2) = b*(y2-y1)
        //    解得：a=z*(y2-y1), b=z*(x1-x2)
        // 归一化时 z 会被约去，令 z=1，得 a=y2-y1, b=x1-x2
        // 把上述a,b代入 a*x1+b*y1+c=0 解得 c=x2*y1-y2*x1
        real a = p2.y - p1.y;
        real b = p1.x - p2.x;
        real c = (p2.x * p1.y) - (p2.y * p1.x);
        // 归一化到 a^2 + b^2 = 1
        real coef = sqrt(a*a + b*b);
        a /= coef;
        b /= coef;
        c /= coef;
        // 测试数据，计算可能的局内点
        real error = 0.0;
        int n_inlier = 0;
        for(int i=0; i<n_sample; ++i){
            real err_i = fabs(a*arr[i].x + b*arr[i].y + c);
            if(err_i < threshold){ // 若低于阈值，则为局内点
                ++n_inlier;
                error += err_i;
            }
        }
        if(n_inlier) error /= n_inlier; 
        // 若有足够多的点被归类为局内点
        if(static_cast<real>(n_inlier)/static_cast<real>(n_sample) > 0.7){
            if(error < min_error){ // 若新模型更好
                min_error = error;
                a_res = a;
                b_res = b;
                c_res = c;
            }
        }
    }
    return {a_res, b_res, c_res};
}

int main(){
    int N=0;
    cin>> N;
    vector<Node> arr(N, {0,0});
    for(int i=0; i<N; ++i){
        cin>> arr[i].x >> arr[i].y;
    }
    auto res = ransac(arr, 0.1, 1e-2, 10);
    cout<< res[0] << " " << res[1] << " " << res[2] << endl;
}

C++ 解法, 执行用时: 5ms, 内存消耗: 376KB, 提交时间: 2020-08-19

#include<bits/stdc++.h>
using namespace std;

#define real float
struct Node{
    real x, y;
};

vector<real> ransac(vector<Node>& arr,
                   real outlier_prob=0.1,
                   real accept_prob=1e-3,
                   real threshold=10.0){
    default_random_engine generator;
    int n_sample = arr.size();
    int n = 2;   // 拟合模型需要的最小数据量
    // 计算理论最大迭代次数
    real inlier_prob = 1 - outlier_prob;
    real sample_fail = 1 - inlier_prob*inlier_prob;
    int k = log(accept_prob) / log(sample_fail);
    
    real a_res, b_res, c_res;
    real min_error = numeric_limits<real>::max();
    while(k--){
        // 随机采样 n 个样本
        uniform_int_distribution<int> sampler(0, n_sample-1);
        int idx1=0, idx2=0;
        while(idx1==idx2){
            idx1 = sampler(generator);
            idx2 = sampler(generator);
        }
        Node p1=arr[idx1], p2=arr[idx2];
        // 拟合模型：a*x1+b*y1+c=0 和 a*x2+b*y2+c=0
        // 两式相减：a*(x1-x2) = b*(y2-y1)
        //    解得：a=z*(y2-y1), b=z*(x1-x2)
        // 归一化时 z 会被约去，令 z=1，得 a=y2-y1, b=x1-x2
        // 把上述a,b代入 a*x1+b*y1+c=0 解得 c=x2*y1-y2*x1
        real a = p2.y - p1.y;
        real b = p1.x - p2.x;
        real c = (p2.x * p1.y) - (p2.y * p1.x);
        // 归一化到 a^2 + b^2 = 1
        real coef = sqrt(a*a + b*b);
        a /= coef;
        b /= coef;
        c /= coef;
        // 测试数据，计算可能的局内点
        real error = 0.0;
        int n_inlier = 0;
        for(int i=0; i<n_sample; ++i){
            real err_i = fabs(a*arr[i].x + b*arr[i].y + c);
            if(err_i < threshold){ // 若低于阈值，则为局内点
                ++n_inlier;
                error += err_i;
            }
        }
        // 若有足够多的点被归类为局内点
        if(static_cast<real>(n_inlier)/static_cast<real>(n_sample) > 0.7){
            if(error < min_error){ // 若新模型更好
                min_error = error;
                a_res = a;
                b_res = b;
                c_res = c;
            }
        }
    }
    return {a_res, b_res, c_res};
}

int main(){
    int N=0;
    cin>> N;
    vector<Node> arr(N, {0,0});
    for(int i=0; i<N; ++i){
        cin>> arr[i].x >> arr[i].y;
    }
    auto res = ransac(arr, 0.1, 1e-4, 10);
    cout<< res[0] << " " << res[1] << " " << res[2] << endl;
}