RANSAC原理及直线拟合(python动态图解)
一、簡介
????????隨機采樣一致性(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擬合直線模型代碼)
總結
以上是生活随笔為你收集整理的RANSAC原理及直线拟合(python动态图解)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: Awt学习总结
- 下一篇: Android学习笔记 56. TabL