RANSAC算法

RANSAC算法

RANSAC是一种用于估计模型参数的迭代方法。当根据观测到的数据进行模型拟合时，如果数据存在outliers时，会对最终模型的拟合结果造成影响。

因此，RANSAC可以认为是一种outliers检测方法。

比如，在对观察到的点集进行直线拟合时，采用最小二乘对所有观测点进行拟合，outliers的存在会对误差占据主导地位，导致拟合失败。这时，我们就需要在拟合之前，利用RANSAC对outliers进行剔除，然后再利用最小二乘得到鲁棒的拟合结果。

整体概述

参考源
RANSAC算法在迭代过程中，主要包括两个步骤，即：

• 从数据集中随机选择一组能够拟合模型的数据，且数据量刚好能够拟合模型。
• 给定一个误差阈值，检查所有数据是否符合上述模型，符合就是内点，否则就是外点。

RANSAC不断迭代上述两个步骤，直到获得足够多的内点。

主要流程

根据RANSAC在二维运动估计中的运用案例，写一下其主要流程：

1. 随机从所有的点集中选择3个不重复的点作为一个子集（之所以选择3个点，因为要想求解仿射变换至少需要3个点才能提供足够的约束，来求解6个未知参数）。
2. 利用上述子集进行模型拟合。
3. 利用上述模型，检查所有点，记录内点的个数。
4. 更新最优模型及对应的内点个数。
5. 重复上述1-4步，直到满足最大迭代次数。
6. 利用上述最优模型，对所有点进行筛选，得到所有内点。
7. 利用上述所有内点，再次进行模型拟合，获得鲁棒结果。（之所以再次拟合，是因为在迭代过程中获得模型结果，都是用最少的点进行拟合的结果）

代码

作为一种参考，下面给出opencv中一种实现，结合上述主要流程食用更佳。
在运动估计前，首先利用RANSAC算法，对光流追踪的点集进行筛选后，再利用estimateGlobalMotionLeastSquares()进行运动估计。

Mat estimateGlobalMotionRansac(
        InputArray points0, InputArray points1, int model, const RansacParams &params,
        float *rmse, int *ninliers)
{
    CV_INSTRUMENT_REGION();

    CV_Assert(model <= MM_AFFINE);
    CV_Assert(points0.type() == points1.type());
    const int npoints = points0.getMat().checkVector(2);
    CV_Assert(points1.getMat().checkVector(2) == npoints);

    if (npoints < params.size)
        return Mat::eye(3, 3, CV_32F);

    const Point2f *points0_ = points0.getMat().ptr<Point2f>();
    const Point2f *points1_ = points1.getMat().ptr<Point2f>();
    const int niters = params.niters(); // 估算推导出的迭代数

    // current hypothesis
    std::vector<int> indices(params.size);
    std::vector<Point2f> subset0(params.size);
    std::vector<Point2f> subset1(params.size);

    // best hypothesis
    std::vector<int> bestIndices(params.size);

    Mat_<float> bestM;
    int ninliersMax = -1;

    RNG rng(0);
    Point2f p0, p1;
    float x, y;

    for (int iter = 0; iter < niters; ++iter)
    {
        // params.size--Minimum number of data points required to estimate model parameters.
        for (int i = 0; i < params.size; ++i)
        {
            bool ok = false;
            while (!ok)
            {
                ok = true;
                indices[i] = static_cast<unsigned>(rng) % npoints; // randomly selected values from data
                // The randomly selected data cannot be the same as the previous data.Otherwise, random selection fails
                for (int j = 0; j < i; ++j) 
                    if (indices[i] == indices[j])
                        { ok = false; break; } // if fails, continue to randomly select
            }
        }
        // based the above index, extract the data used to fit the models
        for (int i = 0; i < params.size; ++i)
        {
            subset0[i] = points0_[indices[i]];
            subset1[i] = points1_[indices[i]];
        }

        // cal models
        Mat_<float> M = estimateGlobalMotionLeastSquares(subset0, subset1, model, 0);

        // put all the data into the model above, then cal the nums of maybeInliers
        int numinliers = 0;
        for (int i = 0; i < npoints; ++i)
        {
            p0 = points0_[i];
            p1 = points1_[i];
            x = M(0,0)*p0.x + M(0,1)*p0.y + M(0,2);
            y = M(1,0)*p0.x + M(1,1)*p0.y + M(1,2);
            if (sqr(x - p1.x) + sqr(y - p1.y) < params.thresh * params.thresh)
                numinliers++;
        }
        // record the model corresponding to the max nums of maybeInliers
        if (numinliers >= ninliersMax)
        {
            bestM = M;
            ninliersMax = numinliers;
            bestIndices.swap(indices);
        }
    }

    // after iters, the nums of maybeInliers less then the nums of points required,
    // use Least squares
    if (ninliersMax < params.size)
    {
        // compute RMSE
        for (int i = 0; i < params.size; ++i)
        {
            subset0[i] = points0_[bestIndices[i]];
            subset1[i] = points1_[bestIndices[i]];
        }
        bestM = estimateGlobalMotionLeastSquares(subset0, subset1, model, rmse);
    }
    else
    {
        subset0.resize(ninliersMax);
        subset1.resize(ninliersMax);
        // select maybeInliers
        for (int i = 0, j = 0; i < npoints && j < ninliersMax ; ++i)
        {
            p0 = points0_[i];
            p1 = points1_[i];
            x = bestM(0,0)*p0.x + bestM(0,1)*p0.y + bestM(0,2);
            y = bestM(1,0)*p0.x + bestM(1,1)*p0.y + bestM(1,2);
            if (sqr(x - p1.x) + sqr(y - p1.y) < params.thresh * params.thresh)
            {
                subset0[j] = p0;
                subset1[j] = p1;
                j++;
            }
        }
        bestM = estimateGlobalMotionLeastSquares(subset0, subset1, model, rmse);
    }

    if (ninliers)
        *ninliers = ninliersMax;

    return std::move(bestM);
}