RANSAC算法应用

问题描述

在对数据进行建模时，我们往往会遇到数据中存在异常数据（outlier，一般称为外点）的情况，如图，用一条直线拟合图中的点，其中左下和右上的点明显是外点。

此时使用最小二乘法，则会由于外点的存在偏离正确模型较多，如下图，因此我们需要找出并剔除这些外点，从而得到一个较好的模型。

RANSAC算法原理

RANSAC（Random Sample Consensus）基本思想是通过多次随机抽取部分数据进行模型估计，然后使用全部数据评估这些模型的正确性并进行比较，得出最好的模型。具体步骤如下。

1. 随机抽取数据集中的一些样本 Si ，估计一个模型 Mi 。
2. 通过一定标准检验能较好的符合模型Mi的数据集S′i，称为内点（inliner）。
3. 若size(Si)>N，N为人为设定的阈值，则算法结束，得到模型M
4. 否则重新进行步骤1～3.若一定重复次数后算法依然没结束，则选取内点数目最多的模型。
5. 使用得到模型的内点再次计算更优的模型M′

使用RANSAC算法得到的拟合直线如下

具体实现

RANSAC进行直线拟合的代码如下：

---

ransac.hpp

#include <opencv2\opencv.hpp>#include <vector>class RansacModelEstimator{public:    RansacModelEstimator() {};    cv::Mat run(std::vector<cv::Point3d> point_set);    int iterate(std::vector<int>& mask);    void estimateParm(std::vector<cv::Point3d>& initial_point);    void estimateInliners(std::vector<int>& mask);    cv::Mat estimateParm1(std::vector<cv::Point3d>& initial_point);    double estimateDeparture(cv::Point3d& p);private:    std::vector<cv::Point3d> point_set;    cv::Point3d theta;    std::vector<int> rand_index;};

---

ransac.cpp

#include "ransac.hpp"#include <random>#include <vector>using namespace std;using namespace cv;double RansacModelEstimator::estimateDeparture(Point3d& p){    return fabs(p.dot(theta)/ norm(Point2d(theta.x,theta.y)));}void RansacModelEstimator::estimateInliners(vector<int>& mask){    mask.clear();    for (size_t i = 0; i < point_set.size(); i++)    {        if (estimateDeparture(point_set[i]) < 0.3)            mask.push_back(i);    }}void RansacModelEstimator::estimateParm(vector<cv::Point3d>& initial_point){    theta = initial_point[0].cross(initial_point[1]);}Mat RansacModelEstimator::estimateParm1(vector<cv::Point3d>& inliner_point){    Mat constrains(inliner_point.size(),3,CV_64F);    Mat b(inliner_point.size(), 1, CV_64F);    b.setTo(0);    for (int i = 0;i<inliner_point.size();i++)    {        constrains.at<double>(i, 0) = inliner_point[i].x;        constrains.at<double>(i, 1) = inliner_point[i].y;        constrains.at<double>(i, 2) = inliner_point[i].z;    }    Mat tconstrains;    Mat result;    cv::SVD::solveZ(constrains, result);    return result;}int RansacModelEstimator::iterate(vector<int>& mask){    static std::random_device rd;    static std::default_random_engine e(rd());    static std::uniform_int_distribution<> u(0, point_set.size()-1);    vector<cv::Point3d> initial_point;    int last=-1;    for (size_t i = 0; i < 2; i++)    {        int index = u(e);        while (index == last)            index = u(e);        initial_point.push_back(point_set[index]);        last = index;    }    estimateParm(initial_point);    estimateInliners(mask);    return mask.size();}Mat RansacModelEstimator::run(vector<cv::Point3d> point_set_){    point_set = point_set_;    vector<int> optimal_mask;    vector<int> current_mask;    int optimal_inliner_num = 0;    int inliner_num;    for (size_t i = 0; i < 300; i++)    {        inliner_num = iterate(current_mask);        if (optimal_inliner_num < inliner_num)        {            optimal_mask = current_mask;            optimal_inliner_num = inliner_num;        }    }    vector<Point3d> inliners;    for (size_t i = 0; i < optimal_inliner_num; i++)    {        inliners.push_back(point_set[optimal_mask[i]]);    }    return estimateParm1(inliners);}int main(){    RansacModelEstimator est;    Mat result = est.run({ {1,1,1},{2,2,1},{3,4,1},{5,4,1},{7,6.9,1},{6,6.1,1},{1,111,1},{23,2,1},{11.8,12.1,1} });    std::cout << result << endl;    return 0;}