RANSAC

Posted on Edited on In Waline:

RANSAC（RAndom SAmple Consensus，随机采样一致）算法是从一组含有“外点”（outliers）的数据中正确估计数学模型参数的迭代算法。本文是关于 Sinkhorn 相关网络资料的整理。

RANSAC

1 简介

定义：

  • RANSAC（RAndom SAmple Consensus，随机采样一致）算法是从一组含有“外点”（outliers）的数据中正确估计数学模型参数的迭代算法。
  • “外点”一般指数据中的噪声，比如匹配中的误匹配和估计曲线中的离群点。所以，RANSAC 也是一种“外点”检测算法。
  • RANSAC 算法是一种不确定算法，它只能在一种概率下产生结果，并且这个概率会随着迭代次数的增加而加大。

基本假设：

  • 数据是由“内点”和“外点”组成的。“内点”就是组成模型参数的数据，“外点”就是不适合模型的数据。同时，RANSAC 假设：给定一组“内点”数据集合（通常较小），存在一个程序可以估计出符合“内点”的模型。

2 算法基本思想

基本思想：通过反复选择数据集去估计模型，一直迭代到估计出认为比较好的模型。

算法步骤-1：

  1. 随机地从数据集 SS 中选择 ss 个数据点，一般取 ss 能刚好组成一个最小子集（满足模型估计的最小数据集，如直线拟合为 2 个点），以 ss 个数据点估计待求模型；
  2. 根据前一步的模型，遍历 SS 中的每一个数据点，判断是否满足这个模型，即是否在距离阈值 tt 内：如果在，则为 inlier，将这个点加入 inlier 集合（一致集）SiS_i 中；
  3. 如果 SiS_i 的大小大于阈值 TT，则用 SiS_i 中的所有点重新估计模型，更新迭代次数 kk，并结束算法循环；
  4. 如果 SiS_i 的大小小于阈值 TT，选择一个新的最小子集 ss，重复上述过程；
  5. 经过 kk 次实验选择一个最大一致集 SiS_i，并用 SiS_i 的所有点重新估计模型。

算法步骤-2：

  1. 选择出可以估计出模型的最小数据集（对于直线拟合来说就是两个点，对于计算单应矩阵就是 4 个点）；
  2. 使用这个数据集来计算出数据模型；
  3. 将所有数据带入这个模型，计算出 inlier 的数目（在一定误差范围内的适合当前迭代推出的模型的数据）；
  4. 比较当前模型和之前推出的最好的模型的 inlier 的数量，记录最大 inlier 数的模型参数和 inlier 数；
  5. 重复 1-4 步，直到迭代结束或者当前模型足够好（ inlier 数目大于一定数量）。

3 迭代次数推导

假设 inlier 在数据中的占比为 tt：

t=ninliersninliers+noutlierst = \frac{n_{\text{inliers}}}{n_{\text{inliers}} + n_{\text{outliers}}}

那么每次计算模型使用 NN 个点的情况下，选取的点至少有一个 outlier 的概率为：

1−tN1 - t^N

即在迭代 kk 次的情况下，(1−tN)k(1 - t^N)^k 就是 kk 次迭代计算模型都至少采样到一个 outlier 去计算模型的概率。那么能采样到正确的 NN 个点去计算出正确模型的概率就是：

P=1−(1−tN)kP = 1 - (1 - t^N)^k

从而求得迭代次数为：

k=log⁡(1−P)log⁡(1−tN)k = \frac{\log (1 - P)}{\log (1 - t^N)}

  • tt 通常是一个先验值。若事先不知道 tt 的取值，则使用自适应迭代次数的方法：一开始设定一个无穷大的迭代次数，然后每次更新模型参数时，用当前的 inlier 比值作为 tt 来估算迭代次数。

4 算法伪码

输入：

  • data：观测数据
  • model：适应于数据的模型
  • n：适用于模型的最小数据个数
  • k：算法的迭代次数
  • t：用于决定数据是否适应于模型的阈值（判断是否为 inlier 的阈值）
  • d：半段模型是否适用于数据集的 inlier 数目

输出：

  • best_model：与数据最匹配的模型参数，没找到则返回 NULL
  • best_consensus_set：估计出模型的数据点
  • best_error：与数据相关的模型误差

算法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
iterations = 0;
best_model = NULL;
best_consensus_set = NULL;
best_error = DBL_MAX;  // 无穷大

while (iterations < k) {
    maybe_inliers = 从数据集中随机选择 n 个点;
    maybe_model = 求得适合于 maybe_inliers 的模型参数;
    consensus_set = maybe_inliers;
    
    for (data 中不属于 maybe_inliers 的点 p) {
        if (数据点 p 适合 maybe_model，且误差小于 t) {
            将点 p 添加到 consensus_set;
        }
    }
    
    if (|consensus_set| > d) {
        // 已经找到可能较好的模型，检测模型误差
        better_model = 适合于 consensus_set 中所有点的模型参数;
        this_error = bette_model 产生的误差;
        
        if (this_error < best_error) {
            // 发现了更好的模型
            best_model = better_model;
            best_consensus_set = consensus_set;
            best_error = this_error;
        }
    }
    
    更新迭代总次数 k;
    
    iterations++;
}

return best_model, best_consensus_set, best_error;

5 用 Python 拟合直线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
import numpy as np
import matplotlib.pyplot as plt
import random
import math

# 数据量。
SIZE = 50
# 产生数据。np.linspace 返回一个一维数组，SIZE指定数组长度。
# 数组最小值是 0，最大值是 10。所有元素间隔相等。
X = np.linspace(0, 10, SIZE)
Y = 3 * X + 10

fig = plt.figure()
# 画图区域分成1行1列。选择第一块区域。
ax1 = fig.add_subplot(1,1, 1)
# 标题
ax1.set_title("RANSAC")

# 让散点图的数据更加随机并且添加一些噪声。
random_x = []
random_y = []
# 添加直线随机噪声
for i in range(SIZE):
    random_x.append(X[i] + random.uniform(-0.5, 0.5)) 
    random_y.append(Y[i] + random.uniform(-0.5, 0.5)) 
# 添加随机噪声 outliers
for i in range(SIZE):
    random_x.append(random.uniform(0,10))
    random_y.append(random.uniform(10,40))
RANDOM_X = np.array(random_x) # 散点图的横轴。
RANDOM_Y = np.array(random_y) # 散点图的纵轴。

# 画散点图。
ax1.scatter(RANDOM_X, RANDOM_Y)
# 横轴名称。
ax1.set_xlabel("x")
# 纵轴名称。
ax1.set_ylabel("y")

# 使用RANSAC算法估算模型
# 迭代最大次数，每次得到更好的估计会优化iters的数值
iters = 100000
# 数据和模型之间可接受的差值
sigma = 0.25
# 最好模型的参数估计和内点数目
best_a = 0
best_b = 0
pretotal = 0
# 希望的得到正确模型的概率
P = 0.99

i = 0
while i < iters:
    i += 1
    
    # 随机在数据中红选出两个点去求解模型
    sample_index = random.sample(range(SIZE * 2), 2)
    x_1 = RANDOM_X[sample_index[0]]
    x_2 = RANDOM_X[sample_index[1]]
    y_1 = RANDOM_Y[sample_index[0]]
    y_2 = RANDOM_Y[sample_index[1]]

    # y = ax + b 求解出a，b
    a = (y_2 - y_1) / (x_2 - x_1)
    b = y_1 - a * x_1

    # 算出内点数目
    total_inlier = 0
    for index in range(SIZE * 2):
        y_estimate = a * RANDOM_X[index] + b
        if abs(y_estimate - RANDOM_Y[index]) < sigma:
            total_inlier = total_inlier + 1

    # 判断当前的模型是否比之前估算的模型好
    if total_inlier > pretotal:
        iters = int(math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE * 2), 2)))
        pretotal = total_inlier
        best_a = a
        best_b = b

    # 判断是否当前模型已经符合超过一半的点
    # 注意有 2 * SIZE 个点，一半为 outliers
    if total_inlier > SIZE:
        break

# 用我们得到的最佳估计画图
Y = best_a * RANDOM_X + best_b

# 直线图
ax1.plot(RANDOM_X, Y)
text = "best_a = " + str(best_a) + "\nbest_b = " + str(best_b)
plt.text(5,10, text,
         fontdict={'size': 8, 'color': 'r'})
plt.show()

6 ORB 中的 RANSAC 调用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include <opencv2/opencv.hpp>
#include <iostream>
#include <opencv2/xfeatures2d.hpp>
#include <opencv2/features2d/features2d.hpp>

void extracte_orb(
    cv::Mat                   input, 
    std::vector<cv::KeyPoint> &keypoint, 
    cv::Mat                   &descriptor) {
    cv::Ptr<cv::ORB> f2d = cv::ORB::create(500);
    f2d->detect(input, keypoint);
    cv::Mat image_with_kp;
    f2d->compute(input, keypoint, descriptor);
    cv::drawKeypoints(input, keypoint, image_with_kp, cv::Scalar::all(-1), cv::DrawMatchesFlags::DRAW_RICH_KEYPOINTS);
	cv::imwrite("orb" + std::to_string(random()) + ".png", image_with_kp);
}

void match_two_image(
    cv::Mat                   image1, 
    cv::Mat                   image2, 
    std::vector<cv::KeyPoint> keypoint1,
    std::vector<cv::KeyPoint> keypoint2, 
    cv::Mat                   descriptor1, 
    cv::Mat                   descriptor2) {
    cv::BFMatcher matcher(cv::NORM_HAMMING);
    std::vector<cv::DMatch> matches;
    matcher.match(descriptor1,descriptor2, matches);
    cv::Mat good_matches_image;
    cv::drawMatches(image1, keypoint1, image2, keypoint2,
                    matches, good_matches_image, cv::Scalar::all(-1), cv::Scalar::all(-1),
                    std::vector<char>(), cv::DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS);
    cv::imwrite("good_matches_image.png", good_matches_image);
    {
        std::vector<cv::KeyPoint> RAN_KP1, RAN_KP2;
        std::vector<cv::Point2f> keypoints1, keypoints2;
        for (int i = 0; i < matches.size(); i++) {
            keypoints1.push_back(keypoint1[matches[i].queryIdx].pt);
            keypoints2.push_back(keypoint2[matches[i].trainIdx].pt);
            RAN_KP1.push_back(keypoint1[matches[i].queryIdx]);
            RAN_KP2.push_back(keypoint2[matches[i].trainIdx]);
        }

        std::vector<uchar> RansacStatus;
        cv::findFundamentalMat(keypoints1, keypoints2, RansacStatus, cv::FM_RANSAC);
        std::vector<cv::KeyPoint> ransac_keypoints1, ransac_keypoints2;
        std::vector<cv::DMatch> ransac_matches;
        int index = 0;
        for (size_t i = 0; i < matches.size(); i++) {
            if (RansacStatus[i] != 0) {
                ransac_keypoints1.push_back(RAN_KP1[i]);
                ransac_keypoints2.push_back(RAN_KP2[i]);
                matches[i].queryIdx = index;
                matches[i].trainIdx = index;
                ransac_matches.push_back(matches[i]);
                index++;
            }
        }
        cv::Mat after_ransac_sift_match;
        cv::drawMatches(image1, ransac_keypoints1, image2, ransac_keypoints2,
                        ransac_matches, after_ransac_sift_match, cv::Scalar::all(-1), cv::Scalar::all(-1),
                        std::vector<char>(), cv::DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS);
        cv::imwrite("after_ransac_orb_match.png", after_ransac_sift_match);
    }
}

int main(int argc, char *argv[]) {
    cv::Mat image1 = cv::imread(argv[1]);
    cv::Mat image2 = cv::imread(argv[2]);
    std::vector<cv::KeyPoint> keypoint1,keypoint2;
    cv::Mat descriptor1, descriptor2;
    extracte_orb(image1, keypoint1, descriptor1);
    extracte_orb(image2, keypoint2, descriptor2);
    match_two_image(image1, image2, keypoint1, keypoint2, descriptor1, descriptor2);
    return 0;
}

7 OpenCV 对 RANSAC 的应用

7.1 findHomography()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
CV_EXPORTS_W 
Mat findHomography( 
    InputArray   srcPoints, 
    InputArray   dstPoints,
    int          method = 0, 
    double       ransacReprojThreshold = 3,
    OutputArray  mask = noArray(), 
    const int    maxIters = 2000,
    const double confidence = 0.995
);

/** @overload */
CV_EXPORTS 
Mat findHomography( 
    InputArray srcPoints, 
    InputArray dstPoints,
    OutputArray mask, 
    int method = 0, 
    double ransacReprojThreshold = 3 
);

参数说明：

  • srcPoints：源平面的坐标矩阵，可以是 CV_32FC2 类型，也可以是 vector<Point2f> 类型。

  • dstPoints：目标平面的坐标矩阵，可以是 CV_32FC2 类型，也可以是 vector<Point2f> 类型。

  • method：计算单应矩阵所使用的方法。不同的方法对应不同的参数，具体如下：

    • 0：利用所有点的常规方法；

      RANSAC：基于 RANSAC 的鲁棒算法；

      LMEDS：最小中值鲁棒算法；

      RHO：基于 PROSAC 的鲁棒算法。

  • ransacReprojThreshold：将点对视为内点的最大允许重投影错误阈值（仅用于 RANSAC 和 RHO 方法）。如果

    ∥dstPointsi−convertPointsHomogeneous⁡(H∗srcPointsi)∥>ransacReprojThreshold\| dstPoints_i - \operatorname{convertPointsHomogeneous}(H * srcPoints_i) \| \\ > ransacReprojThreshold

    则点被认为是个 outlier（即错误匹配点对）。若 srcPoints 和 dstPoints 是以像素为单位的，则该参数通常设置在 1 到 10 的范围内。

  • mask：可选输出掩码矩阵，通常由鲁棒算法（RANSAC 或 LMEDS）设置。注意，输入掩码矩阵是不需要设置的。

  • maxIters：RANSAC 算法的最大迭代次数，默认为 2000。

  • confidence：可信度，取值为 [0,1][0, 1]。

7.2 findFundamental()

1
2
3
4
5
6
7
8
9
10
CV_EXPORTS_W 
Mat findFundamentalMat(
    InputArray  points1, 
    InputArray  points2,
    // FM_7POINT, FM_8POINT, FM_LMEDS, FM_RANSAC
    int         method = FM_RANSAC,  
    double      param1 = 3., 
    double      param2 = 0.99,                              
    OutputArray mask = noArray()
);

参数说明：

  • points1：第一张图像的点；
  • points2：第二章图像的点；
  • param1：该参数用于 RANSAC 算法，它是从点到对极线的最大距离（以像素为单位），超过该距离，该点被视为异常值，并且不用于计算最终的基础矩阵。 可以将其设置为 1-3 ，具体取决于点定位的精度，图像分辨率和图像噪声。
  • param2：该参数仅在 RANSAC 算法以及 LMedS 算法中， 它指定了估计矩阵正确的期望置信度（概率）。

7.3 solvePnPRansac()solvePnPRansac()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
bool cv::solvePnPRansac(
    InputArray  objectPoints,
    InputArray  imagePoints,
    InputArray  cameraMatrix,
    InputArray  distCoeffs,
    OutputArray rvec,
    OutputArray tvec,
    bool        useExtrinsicGauss = false,
    int         iterationsCount = 100,
    float       reprojectionError = 8.0,
    double      confidence = 0.99,
    OutputArray inliers = noArray(),
    int         flags = SOLVEPNP_INTERATIVE
);

参数说明：

  • objectPoints：目标坐标系下的目标点，3×N/N×3,1 channel3 \times N / N \times 3, 1\ \text{channel} 或者 1×N/N×1,3 channel1 \times N / N \times 1, 3\ \text{channel}，NN 为点的数量，也可以是 vector<Point3f>。
  • imagePoints：目标对应的图像点，2×N/N×2,1 channel2 \times N / N \times 2, 1\ \text{channel} 或者 1×N/N×1,2 channel1 \times N / N \times 1, 2\ \text{channel}，也可以是 vector<Point2f>。
  • cameraMatrix：相机内参矩阵。
  • distCoeffs：畸变参数向量，可以是 NULL 或 4、5、8、12 或 14 个变量组成的向量。
  • rvec：输出的旋转向量，结合 tvec 将点从模型坐标系转换到相机坐标系。
  • tvec：输出的位移向量。
  • useExtrinsicGauss：使用在 SOLVEPNP_INTERATIVE 的参数。若为 true，则函数使用提供的 rvec 和 tvec 作为旋转和平移的初始值。
  • iterationsCount：迭代次数。
  • reprojectionError：RANSAC 算法的 inlier 阈值，表示观测和估计点投影的最大距离阈值。
  • confidence：算法结果的置信度。
  • inliers：包含在 objectPoints 和 imagePoints 的 inlier 的索引。
  • flags：求解 PnP 的方法（见 solvePnP()）。

参考