RANSAC算法应用

问题描述

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

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

RANSAC算法原理

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

1. 随机抽取数据集中的一些样本 $S_i$ ，估计一个模型 $M_i$ 。
2. 通过一定标准检验能较好的符合模型$M_i$的数据集$S_i'$，称为内点（inliner）。
3. 若$size(S_i)>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;
}