【原】手写梯度下降《二》之

From WIKI：

Random sample consensus (RANSAC) is an iterative method to estimate parameters of a mathematical model from a set of observed data that contains outliers, when outliers are to be accorded no influence on the values of the estimates. Therefore, it also can be interpreted as an outlier detection method.^[1] It is a non-deterministic algorithm in the sense that it produces a reasonable result only with a certain probability, with this probability increasing as more iterations are allowed. The algorithm was first published by Fischler and Bolles at SRI International in 1981. They used RANSAC to solve the Location Determination Problem (LDP), where the goal is to determine the points in the space that project onto an image into a set of landmarks with known locations.

A basic assumption is that the data consists of "inliers", i.e., data whose distribution can be explained by some set of model parameters, though may be subject to noise, and "outliers" which are data that do not fit the model. The outliers can come, for example, from extreme values of the noise or from erroneous measurements or incorrect hypotheses about the interpretation of data. RANSAC also assumes that, given a (usually small) set of inliers, there exists a procedure which can estimate the parameters of a model that optimally explains or fits this data.

 伪代码：

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框架下计算相机位姿变换、或者数据拟合等等。其中，最大迭代次数，这个不是固定的，他与当前迭代内点比例等等参数有关，请参考:https://en.wikipedia.org/wiki/Random_sample_consensus   或者 https://www.cnblogs.com/cfantaisie/archive/2011/06/09/2076864.html

下面例子，我们在进行直线拟合，手动添加一些“较为夸张”的噪点，基于RANSAC框架，进行直线拟合，需要说明的是：循环迭代结束的时候，我们可以进一步采取最小二乘等手段进行曲线拟合，这样效果更好。同时，有关迭代次数问题，它涉及到程序的执行时效，所以，实际应用RANSAC框架来处理数据的时候，最大迭代次数，最好每次都更新。

#define _CRT_SECURE_NO_WARNINGS
#include<iostream>
#include<vector>
#include<algorithm>
#include<ctime>
using namespace std;

#include<opencv2/opencv.hpp>
using namespace cv;


// 功能：画点
void drawPoints(vector<Point2d> points, Mat& image, bool is_green)
{
    if (is_green)
    {
        for (int i = 0; i < points.size(); i++)
        {
            circle(image, Point2i(int(points[i].x), int(points[i].y)), 2, cv::Scalar(0, 255, 0), 2, CV_AA);
        }
    }
    else
    {
        for (int i = 0; i < points.size(); i++)
        {
            circle(image, Point2i(int(points[i].x), int(points[i].y)), 2, cv::Scalar(0, 0, 255), 2, CV_AA);
        }
    }

}

// 功能：画线
void drawLine(Point2d begin, Point2d end, Mat& image)
{
    line(image, Point2i(int(begin.x), int(begin.y)),
        Point2i(int(end.x), int(end.y)), Scalar(255, 0, 0), 1, CV_AA);
}

// 功能：将一张图经行纵向镜像
void upDownMirror(Mat& image)
{
    Mat image_cpy = image.clone();
    for (int j = 0; j < image.rows; j++)
    {
        for (int i = 0; i < image.cols; i++)
        {
            image.ptr<cv::Vec3b>(j)[i] = image_cpy.ptr<cv::Vec3b>(image.rows - 1 - j)[i];
        }
    }
}

// 功能：最简单的直线拟合 Ransac 框架

class Ransac
{
public:
    Ransac(vector<double> x, vector<double> y);
    ~Ransac();
    static int getRand(int min, int max);//[min, min + max - 1]
public:
    vector<double> x_, y_;
};

Ransac::Ransac(vector<double> x, vector<double> y) :x_(x), y_(y)
{

}

Ransac::~Ransac()
{

}

int Ransac::getRand(int min, int max)//[min, min + max - 1]
{
    return rand() % max + min;
}

int main()
{
    // time seed 
    srand((unsigned int)time(nullptr));

    // <1>产生随机数
    vector<Point2d> points;
    // 精确解
    const int k = 1;
    const int b = 10;

    // x [0, 390]; y [10, 400] total_numbers = 40
    for (int x = 0; x < 400; x += 10)
    {
        points.push_back(Point2d(x, k*x + b + Ransac::getRand(0, 60)));
    }
    //添加离群点
    
    points[35].y = 100;
    points[36].y = 200;
    points[37].y = 200;
    points[38].y = 120;
    points[39].y = 110;


    const int itr_nums = 300;
    double k_estimate = 0;
    double b_estimate = 0;
    Point2d best_parameters; // [k, b]

                             // 统计纵向距离 y
    int count = 0; // 内点计数器
    int y_threshold = 50; // 判定内点的阈值
    int total = 0; // 内点总数

                   //样本点
    Point2d sample_first;
    Point2d sample_second;
    // 最佳样本点
    Point2d best_sample_first;
    Point2d best_sample_second;

    for (int i = 0; i < itr_nums; i++)
    {
        // <1>、随机抽取两个样本
        int index_first = Ransac::getRand(0, points.size() - 2);
        int index_second = Ransac::getRand(0, points.size() - 2);
        sample_first = points[index_first];
        sample_second = points[index_second];

        if (sample_first == sample_second)
        {
            continue;
        }
        // 计算斜率 k = (y2 - y1)/(x2 - x1)
        k_estimate = (sample_second.y - sample_first.y) / (sample_second.x - sample_first.x);
        // 计算截距 b = y1 - k * x1
        b_estimate = sample_first.y - k_estimate * sample_first.x;
        // <2>、根据距离，来找出所有样本点中的内点，并统计数量
        for (int i = 0; i < points.size(); i++)
        {
            // delta = k * x  + b - y
            double y_error = abs(k_estimate * points[i].x + b_estimate - points[i].y) / sqrt((pow(k_estimate, 2) + 1));
            //cout << "y_error = " << y_error << endl;
            if (y_error < y_threshold)
            {
                count++;
            }
        }
        //  <3>、找出内点数量最多的那一组
        if (count > total)
        {
            total = count;

            best_sample_first = sample_first;
            best_sample_second = sample_second;

            best_parameters.x = k_estimate;
            best_parameters.y = b_estimate;
        }
        count = 0;
    }
    cout << "内点数 = " << total << endl;
    cout << "斜率 k = " << best_parameters.x << "截距 b = " << best_parameters.y << endl;

    // 统计内点
    vector<Point2d> inliners_points;
    for (int i = 0; i < points.size(); i++)
    {
        // delta = k * x  + b - y
        double y_error = abs(best_parameters.x * points[i].x + best_parameters.y - points[i].y) / sqrt((pow(k_estimate, 2) + 1));
        //cout << "y_error = " << y_error << endl;
        if (y_error < y_threshold)
        {
            inliners_points.push_back(points[i]);
        }
    }

    Mat image = Mat::zeros(500, 500, CV_8UC3);

    drawLine(best_sample_first, best_sample_second, image);
    drawPoints(points, image, true);
    drawPoints(inliners_points, image, false); //画内点
    upDownMirror(image);
    imshow("image", image);
    waitKey(0);
    return 1;
}

总的样本数是40，内点数35，比例已经算很高了。绿色是外点，红色是内点。当然迭代次数越多，得到的解越精确。