一、简介

        随机采样一致性（Random Sample Consensus，RANSAC）由斯坦福国际研究院的Fischler和Bolles于1981年首次提出[1]。RANSAC算法是一种随机参数估计迭代算法；从一组包含异常数据的样本数据集中，通过迭代的方式，估计已知数学模型的参数，并得到有效样本数据的算法；也可以将其理解为一种点集中离散值的检测方法。

         RANSAC存在2个假设：

①数据由“内点”（inliers）和“外点”（outliers）组成。“内点”是参与模型拟合的数据；“外点”是异常数据，不适用于或不符合模型的数据，噪声等。

②给定一组含有较少“内点”的数据，存在一个可以估计模型参数的过程或程序，得到的模型参数可以解释或适用于局内点。

        RANSAC是数理统计规律和几何规律结合的算法，它对粗差有较好的抵抗力；但是，RANSAC是一种不确定的算法，不一定得到最佳结果，为了得到更佳的结果，必须提高迭代次数。

二、方法流程

        RANSAC通过反复选择数据集，估计模型，直到迭代出符合要求的模型。

 

 2.1 设置先验参数

• 假设有N个点，记为P(P1,P2,…,PN)；拟合模型最少需要n个点。
• 假设拟合2D直线Y=kX+b，最少需要2个点，即n=2；
• 数据适用于模型的阈值，即点到直线的距离DisT；
• “内点”比例阈值为ProT；
• 最大迭代次数为Iter_Number。

2.2 RANSAC的计算流程

①随机选择用于模型估计的最小数据集Q：即随机选择M个点用于拟合2D直线（n≤M≤N），一般令M=n。

②使用数据集Q计算模型L，得到k和b。

③将所有数据P带入模型L，计算“内点”比例：即计算点到直线的距离D，若D≤DisT，则该点为“内点”；计算得到该模型的“内点”比例为Pro_current。

④之前最好模型的“内点”比例与当前模型的“内点”比例比较大小，记录较大者的模型参数和“内点”。

⑤重复1-4步，直到迭代结束或者当前模型的“内点”比例大于ProT。得到模型参数和“内点”等。

        以上是RANSAC的基本流程，为了增加结果精确度，可以和最小二乘等方法结合，例如，利用最后的“内点”数据，通过最小二乘法或其他方法，拟合模型，得到更佳的模型参数。

三、迭代次数的推导[2][3]

        我们需要根据特定的问题和数据集通过实验来确定参数；然而迭代次数可以从理论结果推断。假设“内点”在数据中的占比为Pro_current

        通常情况下，我们事先并不知道Pro_current的值，但是可以给出一些鲁棒的值。假设估计模型需要M个点，Pro_currentM是M个点均是内点的概率。1-Pro_currentM是M个点中至少有一个点为局外点的概率。此时表明我们从数据集中估计出了一个不好的模型。

        (1-Pro_currentM)iter_number表示算法永远都不会选择到n个点均为局内点的概率，就代表永远无法达到期望的内点比例，即

        上式取对数，得：

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

四、代码（python，拟合2D直线）

分享给有需要的人，代码质量勿喷

import numpy as np
import random
import math
import matplotlib.pyplot as pltfig = plt.figure()
ax = fig.add_subplot(1, 1, 1)#region data
N = 40
X = np.linspace(start = 0, stop = 20, num = N)
k = b = 2
Y = k * X + b + np.random.randn(N) * 3data_X = []
data_Y = []
for i in range(N):data_X.append(X[i])data_Y.append(Y[i])#添加粗差
for i in range(5):data_X.append(i+15)data_Y.append(1)N = len(data_X)
#endregion#region 直接用最小二乘拟合直线
# https://blog.csdn.net/xinjiang666/article/details/103782544
kLS = bLS = 0
sumX = sumY = 0
sumXX = sumXY = 0
for i in range(len(data_X)):xi = round(data_X[i])yi = round(data_Y[i])sumX += xisumY += yisumXX += xi * xisumXY += xi * yi
deltaX = sumXX * N - sumX * sumXif (abs(deltaX) > 1e-15):kLS = (sumXY * N - sumX * sumY) / deltaXbLS = (sumXX * sumY - sumX * sumXY) / deltaXYLS_pred = kLS * np.asarray(data_X) + bLS
#endregion#region RANSAC
#region 设置先验参数或期望
n = M = 2           #模型拟合所需的最小点数
Iter_Number = 1000  #最大迭代次数
DisT = 3            #点到直线的距离阈值
ProT = 0.95         #期望的最小内点比例
#endregionbest_k = 0
best_b = 0
best_inliers_number = N*0.1
best_X_inliers=[]
best_Y_inliers=[]
iter = 0
existed_k_b=[] #已经使用过的k、b#region 5 迭代
for i in range(Iter_Number):iter += 1#region 1 随机取M个点用于模型估计indx1, indx2 = random.sample(range(len(data_X)), M)point1_x = data_X[indx1]point1_y = data_Y[indx1]point2_x = data_X[indx2]point2_y = data_Y[indx2]#endregion#region 2 计算模型参数k = round((point2_y - point1_y) / (point2_x - point1_x), 3)b = round(point2_y - k * point2_x, 3)current_k_b = [k,b]if current_k_b in existed_k_b:#防止重复的两个点或参数用于计算continueelse:existed_k_b.append([k,b])#endregion#region 3 计算内点数量（或比例）X_inliers = []Y_inliers = []inliers_current = 0 #当前模型的内点数量for j in range(len(data_X)):x_current = data_X[j]y_current = data_Y[j]dis_current = abs(k * x_current - y_current + b) / math.sqrt(k * k + 1) # 点到拟合直线的距离dis_currentif (dis_current <= DisT):inliers_current += 1X_inliers.append(x_current)Y_inliers.append(y_current)print("第{}次, 内点数量={}, 最佳内点数量={}".format(iter, inliers_current, best_inliers_number))#endregion#region 4 如果当前模型的内点数量大于之前最好的模型的内点数量，跟新模型参数和迭代次数if (N>inliers_current >= best_inliers_number):Pro_current = inliers_current / N       #当前模型的内点比例Pro_currentbest_inliers_number = inliers_current   #更新最优内点的数量best_k = k  #更新模型参数best_b = b  #更新模型参数i = 0       #当前迭代置为0Iter_Number = math.log(1 - ProT) / math.log(1 - pow(Pro_current, M))    #更新迭代次数best_X_inliers = X_inliers  #更新内点best_Y_inliers = Y_inliers  #更新内点print("更新结果：k={}, b={}, 当前内点比例={}, 最佳内点比例={}, 新的迭代次数={}".format(best_k, best_b, Pro_current, best_inliers_number/N, Iter_Number))# endregion#region plotax.cla()ax.set_title("RANSAC and LS")ax.set_xlabel("x")ax.set_ylabel("y")plt.xlim((0, 20))plt.ylim((0, 50))plt.grid(True)#pointsax.scatter(data_X, data_Y, color='k')#LS fit lineax.plot(data_X, YLS_pred, color='b', linewidth=3)#RANSAC fit lineY_R_pred = best_k * np.asarray(data_X) + best_bax.plot(data_X, Y_R_pred, color='g', linewidth=5)plt.legend(['LS fit line','RANSAC fit line','point data',])plt.pause(0.001)#endregion#region XXX 如果最佳的内点比例大于期望比例，则跳出if ((best_inliers_number / N) > ProT):print("终止：内点比例=", (best_inliers_number / N), "大于 期望内点比例=", ProT)break#endregion#endregion
#endregion#region plot inliers
plt.plot(best_X_inliers, best_Y_inliers, "ro")
plt.legend(['LS','RANSAC','inlier','outlier',
])
plt.show()
#endregion

五、结果对比

参考

[1] Fischler M A, Bolles R C. Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography[J]. Communications of the ACM, 1981, 24(6): 381-395.

[2] RANSAC算法理解

[3] RANSAC算法详解(附Python拟合直线模型代码)