RANSAC是“Random Sample Consensus（随机抽样一致）”的缩写。它可以从一组包含“局外点”的观测数据集中，通过迭代方式估计数学模型的参数。它是一种不确定的算法——它有一定的概率得出一个合理的结果；为了提高概率必须提高迭代次数。该算法最早由Fischler和Bolles于1981年提出。

RANSAC的基本假设是：

  1. 数据由“局内点”组成，例如：数据的分布可以用一些模型参数来解释；
  2. “局外点”是不能适应该模型的数据；
  3. 除此之外的数据属于噪声。

局外点产生的原因有：

  1. 噪声的极值；
  2. 错误的测量方法；
  3. 对数据的错误假设.

RANSAC也做了以下假设：给定一组（通常很小的）局内点，存在一个可以估计模型参数的过程；而该模型能够解释或者适用于局内点。

示例

一个简单的例子是从一组观测数据中找出合适的2维直线。假设观测数据中包含局内点和局外点，其中局内点近似的被直线所通过，而局外点远离于直线。简单的最小二乘法不能找到适应于局内点的直线，原因是最小二乘法尽量去适应包括局外点在内的所有点。相反，RANSAC能得出一个仅仅用局内点计算出模型，并且概率还足够高。但是，RANSAC并不能保证结果一定正确，为了保证算法有足够高的合理概率，我们必须小心的选择算法的参数。

算法分析

设，N为样本点个数，n为求解模型需要的最少点个数, w表示内点的概率。

\(w^n\)n个点都是内点的概率

\(1 - w^n\)n个至少有一个外点(采样失败)的概率

\((1 - w^n)^K\)K次采样全部失败的概率

\(p = 1 - (1 - w^n)^K\)K次采样至少有一次成功的概率（置信度）

\[ K = \frac{log(1 - p)}{log(1 - w^n)} \]

值得注意的是，这个结果假设n个点都是独立选择的；也就是说，某个点被选定之后，它可能会被后续的迭代过程重复选定到。这种方法通常都不合理，由此推导出的k值被看作是选取不重复点的上限。例如，要从上图中的数据集寻找适合的直线，RANSAC算法通常在每次迭代时选取2个点，计算通过这两点的直线，要求这两点必须唯一。 ​ 为了得到更可信的参数，标准偏差或它的乘积可以被加到k上。k的标准偏差定义为：

\[ SD(K) = \frac{\sqrt[2]{1 - w^n}}{w^n} \]

算法步骤

  1. .随机采样 n 个点
  2. 对该 n个点拟合模型
  3. 计算其它点到拟合模型的距离，小于一定阈值，当作内点，统计内点个数
  4. .重复K次，选择内点数最多的模型
  5. 利用所有的内点重新估计模型(可选)
Given:
    data – a set of observations
    model – a model to explain observed data points
    n – minimum number of data points required to estimate model parameters
    k – maximum number of iterations allowed in the algorithm
    t – threshold value to determine data points that are fit well by model 
    d – number of close data points required to assert that a model fits well to data

Return:
    bestFit – model parameters which best fit the data (or nul if no good model is found)

iterations = 0
bestFit = nul
bestErr = something really large
while iterations < k {
    maybeInliers = n randomly selected values from data
    maybeModel = model parameters fitted to maybeInliers
    alsoInliers = empty set
    for every point in data not in maybeInliers {
        if point fits maybeModel with an error smaller than t
             add point to alsoInliers
    }
    if the number of elements in alsoInliers is > d {
        % this implies that we may have found a good model
        % now test how good it is
        betterModel = model parameters fitted to all points in maybeInliers and alsoInliers
        thisErr = a measure of how well betterModel fits these points
        if thisErr < bestErr {
            bestFit = betterModel
            bestErr = thisErr
        }
    }
    increment iterations
}
return bestFit

RANSAC的优点是它能鲁棒的估计模型参数。例如，它能从包含大量局外点的数据集中估计出高精度的参数。RANSAC的缺点是它计算参数的迭代次数没有上限；如果设置迭代次数的上限，得到的结果可能不是最优的结果，甚至可能得到错误的结果。RANSAC只有一定的概率得到可信的模型，概率与迭代次数成正比。RANSAC的另一个缺点是它要求设置跟问题相关的阀值。 ​ RANSAC只能从特定的数据集中估计出一个模型，如果存在两个（或多个）模型，RANSAC不能找到别的模型。

循环终止条件

\[ K = \frac{log(1 - p)}{log(1 - w^n)} \]

这里置信度p一般取[0.95, 0.99]，而一般情况下\(\omega^n\)是未知的，因此，可以取最坏条件下类内点的比例，或者在初始状态下设置为最坏条件下的比例，然后随着迭代次数，不断更新为当前最大的类内点比例。

Code

#include <iostream>
#include <cmath>
#include <random>
#include <vector>
#include <array>
#include <set>
#include <memory>
#include <algorithm>
#include <Eigen/Eigen>
#include <opencv2/opencv.hpp>

template <class ModelT, class SampleT, int min_sample>
class RandomSampleConsensus
{
protected:
    typedef ModelT ModelT;
    typedef SampleT SampleT;
    enum{ min_sample_ = min_sample};
public:
    class IModel {
    protected:
        enum{min_sample_ = min_sample};
        typedef SampleT SampleT;
    public:
        IModel() {}
        virtual ~IModel(){}
    public:
        virtual bool generate(const std::array<std::shared_ptr<SampleT>, min_sample_> &samples) = 0;
        virtual bool evaluate(const std::set<std::shared_ptr<SampleT>> &samples, float threshold) = 0;
    public:
        const std::set<std::shared_ptr<SampleT>> &getInliers() const { return inliers_; }
        int getWeight() { return inliers_.size(); }
    protected:
        std::set<std::shared_ptr<SampleT>> inliers_;
    };

public:
    RandomSampleConsensus(float threshold, float confidence = 0.995, float probability_of_inliers = 0.5, int iteration_total = 0)
    {
        reset(threshold, confidence, probability_of_inliers, iteration_total);
    }
    ~RandomSampleConsensus() {}
public:
    void reset(float threshold, float confidence = 0.995, float probability_of_inliers = 0.5, int iteration_total = 0)
    {
        threshold_ = threshold;
        iteration_total_ = iteration_total;
        if (iteration_total_ <= 0) {
            iteration_total_ = std::log(1 - confidence) / std::log(1 - std::pow(probability_of_inliers, min_sample_));
        }
        models_.clear();
        for (int pos = 0; pos < iteration_total_; pos++) {
            models_.push_back(std::make_shared<ModelT>());
        }
        std::srand((unsigned int)(std::time(nullptr)));
    }

    bool evaluate(const std::set<std::shared_ptr<SampleT>> &samples)
    {
        if (!samples.size() || samples.size() < min_sample) {
            return false;
        }
        for (int pos = 0; pos < iteration_total_; pos++) {
            auto &model = models_.at(pos);
            std::set<std::shared_ptr<SampleT>> remaining = samples;
            std::array<std::shared_ptr<SampleT>, min_sample_> model_sample;
            randomGenerateMinSample(remaining, model_sample);
            if (!model->generate(model_sample)) {
                return false;
            }
            if (!model->evaluate(remaining, threshold_)) {
                return false;
            }
        }
        // sort models by inliers count.
        std::sort(models_.begin(), models_.end(), [](const auto &left, const auto &right)
        {
            if (left->getWeight() > right->getWeight()) {
                return true;
            }
            return false;
        });
        return true;
    }
    const std::vector<std::shared_ptr<ModelT>> getModels(void) {return models_;}
protected:
    bool randomGenerateMinSample(std::set<std::shared_ptr<SampleT>> samples, std::array<std::shared_ptr<SampleT>, min_sample_> &model_sample)
    {
        int index = std::rand() % samples.size();
        for(auto iter = model_sample.begin(); iter != model_sample.end(); iter ++){
            int index = std::rand() % samples.size();
            auto sample = samples.begin();
            std::advance(sample, index);
            (*iter) = (*sample);
            samples.erase(sample);
        }
        return true;
    }
private:
    int iteration_total_;
    float threshold_;
    std::vector<std::shared_ptr<ModelT>> models_;
};


class LineModel : public RandomSampleConsensus<LineModel, Eigen::Vector2d, 2>::IModel
{
public:
    LineModel() {}
    virtual ~LineModel() {}
public:
    virtual bool generate(const std::array<std::shared_ptr<SampleT>, min_sample_> &samples)
    {
        if (!samples[0] || !samples[1]) {
            return false;
        }
        // y = kx + b;
        k_ = (samples[0]->y() - samples[1]->y()) / (samples[0]->x() - samples[1]->x());
        b_ = samples[0]->y() - k_ * samples[0]->x();
        model_sample_ = samples;
        return true;
    }
    virtual bool evaluate(const std::set<std::shared_ptr<SampleT>> &samples, float threshold)
    {
        for (auto iter : samples) {
            if (!iter) {
                return false;
            }
            if (threshold > distance(iter)) {
                inliers_.insert(iter);
            }
        }
        return true;
    }
    float distance(std::shared_ptr<SampleT> sample)
    {
        // http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
        // kx - y + b = 0
        return std::fabs(k_ * sample->x() + (-1.) * sample->y() + b_) / std::sqrt(k_ * k_ + (-1. * -1.));
    }
    const std::array<std::shared_ptr<SampleT>, min_sample_> getModelSample() const { return model_sample_; }
private:
    float k_, b_;
    std::array<std::shared_ptr<SampleT>, min_sample_> model_sample_;
};

double slope(int x0, int y0, int x1, int y1)
{
    return (double)(y1 - y0) / (x1 - x0);
}
void draw_line(cv::Mat& output, cv::Point a, cv::Point b, cv::Scalar color, int LineWidth)
{
    double s = slope(a.x, a.y, b.x, b.y);

    cv::Point p(0, 0), q(output.cols, output.rows);

    p.y = -(a.x - p.x) * s + a.y;
    q.y = -(b.x - q.x) * s + b.y;

    cv::line(output, p, q, color, LineWidth, cv::LINE_AA, 0);
}

int main(int argc, char * argv[])
{
    while (true) {
        int canvas_size = 800;
        int points_size = 500;
        int outlier_size = 50;
        cv::Mat canvas(canvas_size, canvas_size, CV_8UC3, cv::Scalar::all(255));

        std::random_device seed;
        std::mt19937 random = std::mt19937(seed());
        std::uniform_int_distribution<int> range(0, canvas_size - 1);

        std::set<std::shared_ptr<Eigen::Vector2d>> observations;
        {
            std::normal_distribution<double> noise(0, 25);
            for (int pos = 0; pos < points_size; pos++)
            {
                int value = range(random);
                cv::Point point(floor(value + noise(random)), floor(value + noise(random)));
                cv::circle(canvas, point, floor(canvas_size / 100) + 3, cv::Scalar(0, 0, 0), 2, cv::LINE_AA);
                observations.insert(std::make_shared<Eigen::Vector2d>(point.x, point.y));
            }
        }
        {
            std::normal_distribution<double> noise(0, 200);
            for (int pos = 0; pos < outlier_size; pos++) {
                int value = range(random);
                cv::Point point(floor(value + noise(random)), floor(value + noise(random))); 
                cv::circle(canvas, point, floor(canvas_size / 100) + 3, cv::Scalar(0, 0, 0), 2, cv::LINE_AA);
                observations.insert(std::make_shared<Eigen::Vector2d>(point.x, point.y));
            }
        }
        RandomSampleConsensus<LineModel, Eigen::Vector2d, 2> line_ransac(20);
        line_ransac.evaluate(observations);

        auto model = line_ransac.getModels().front();
        auto inliers = model->getInliers();
        if (inliers.size())
        {
            for (auto& inlier : inliers)
            {
                cv::Point point(floor(inlier->x()), floor(inlier->y()));
                cv::circle(canvas, point, floor(canvas_size / 100), cv::Scalar(0, 255, 0), -1, cv::LINE_AA);
            }
        }
        auto line = model->getModelSample();
        {
            auto begin = line[0];
            auto end = line[1];
            if (line[0] && line[1])
            {
                cv::Point begin(line[0]->x(), line[0]->y());
                cv::Point end(line[1]->x(), line[1]->y());
                draw_line(canvas, begin, end, cv::Scalar(0, 0, 255), 2);
            }
        }

        cv::imshow("RANSAC", canvas);
        char key = cv::waitKey(1);
        if (key == 27) {
            return 0;
        }
    }
    return 0;
}

demo中使用算法检测直线。下一步可以通过最小二乘拟合出精确度更高的直线。

Reference

https://en.wikipedia.org/wiki/Random_sample_consensus https://en.wikipedia.org/wiki/Hough_transform

Raguram R, Chum O, Pollefeys M, et al. USAC: A Universal Framework for Random Sample Consensus[J]. Pattern Analysis & Machine Intelligence IEEE Transactions on, 2013, 35(8):2022-2038.

O. Chum, J. Matas, and J. Kittler, “Locally Optimized RANSAC,”Proc. DAGM-Symp. Pattern Recognition, pp. 236-243, 2003.

Martin A. Fischler and Robert C. Bolles (June 1981). "Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography". Comm. of the ACM 24: 381–395. doi:10.1145/358669.358692.

David A. Forsyth and Jean Ponce (2003). Computer Vision, a modern approach. Prentice Hall. ISBN 0-13-085198-1.

Richard Hartley and Andrew Zisserman (2003). Multiple View Geometry in Computer Vision (2nd ed.). Cambridge University Press.

P.H.S. Torr and D.W. Murray (1997). "The Development and Comparison of Robust Methods for Estimating the Fundamental Matrix". International Journal of Computer Vision 24: 271–300. doi:10.1023/A:1007927408552.

Ondrej Chum (2005). "Two-View Geometry Estimation by Random Sample and Consensus". PhD Thesis. http://cmp.felk.cvut.cz/~chum/Teze/Chum-PhD.pdf

Sunglok Choi, Taemin Kim, and Wonpil Yu (2009). "Performance Evaluation of RANSAC Family". In Proceedings of the British Machine Vision Conference (BMVC). http://www.bmva.org/bmvc/2009/Papers/Paper355/Paper355.pdf.