RANSAC算法做直线拟合
RANSAC算法之前了解過相關的原理,這兩天利用晚上閑暇的時間,看了一下RANSAC算法的Python代碼實現,這方面的資料很多了,這里就不在重復。在分析該RANSAC.py代碼之前,想用自己的對RANSAC的理解對其做下總結。
在實際應用中獲取到的數據,常常會包含有噪聲數據,這些噪聲數據會使對模型的構建造成干擾,我們稱這樣的噪聲數據點為outliers,那些對于模型構建起積極作用的我們稱它們為inliers,RANSAC做的一件事就是先隨機的選取一些點,用這些點去獲得一個模型(這個講得有點玄,如果是在做直線擬合的話,這個所謂的模型其實就是斜率),然后用此模型去測試剩余的點,如果測試的數據點在誤差允許的范圍內,則將該數據點判為inlier,否則判為outlier。inliers的數目如果達到了某個設定的閾值,則說明此次選取的這些數據點集達到了可以接受的程度,否則繼續前面的隨機選取點集后所有的步驟,不斷重復此過程,直到找到選取的這些數據點集達到了可以接受的程度為止,此時得到的模型便可認為是對數據點的最優模型構建。
在Cookbook/RANSAC中給出的是一個用RANSAC做直線擬合的例子。這個例子非常的直觀,而且代碼也很簡短易懂,為便于后面詳細解讀該代碼,這里把它貼出來:
1 # -*- coding: utf-8 -*- 2 import numpy 3 import scipy # use numpy if scipy unavailable 4 import scipy.linalg # use numpy if scipy unavailable 5 import pylab 6 7 ## Copyright (c) 2004-2007, Andrew D. Straw. All rights reserved. 8 9 def ransac(data,model,n,k,t,d,debug=False,return_all=False): 10 """fit model parameters to data using the RANSAC algorithm 11 12 This implementation written from pseudocode found at 13 http://en.wikipedia.org/w/index.php?title=RANSAC&oldid=116358182 14 15 Given: 16 data - a set of observed data points # 可觀測數據點集 17 model - a model that can be fitted to data points # 18 n - the minimum number of data values required to fit the model # 擬合模型所需的最小數據點數目 19 k - the maximum number of iterations allowed in the algorithm # 最大允許迭代次數 20 t - a threshold value for determining when a data point fits a model #確認某一數據點是否符合模型的閾值 21 d - the number of close data values required to assert that a model fits well to data 22 Return: 23 bestfit - model parameters which best fit the data (or nil if no good model is found) 24 """ 25 iterations = 0 26 bestfit = None 27 besterr = numpy.inf 28 best_inlier_idxs = None 29 while iterations < k: 30 maybe_idxs, test_idxs = random_partition(n,data.shape[0]) 31 maybeinliers = data[maybe_idxs,:] 32 test_points = data[test_idxs] 33 maybemodel = model.fit(maybeinliers) 34 test_err = model.get_error( test_points, maybemodel) 35 also_idxs = test_idxs[test_err < t] # select indices of rows with accepted points 36 alsoinliers = data[also_idxs,:] 37 if debug: 38 print 'test_err.min()',test_err.min() 39 print 'test_err.max()',test_err.max() 40 print 'numpy.mean(test_err)',numpy.mean(test_err) 41 print 'iteration %d:len(alsoinliers) = %d'%( 42 iterations,len(alsoinliers)) 43 if len(alsoinliers) > d: 44 betterdata = numpy.concatenate( (maybeinliers, alsoinliers) ) 45 bettermodel = model.fit(betterdata) 46 better_errs = model.get_error( betterdata, bettermodel) 47 thiserr = numpy.mean( better_errs ) 48 if thiserr < besterr: 49 bestfit = bettermodel 50 besterr = thiserr 51 best_inlier_idxs = numpy.concatenate( (maybe_idxs, also_idxs) ) 52 iterations+=1 53 if bestfit is None: 54 raise ValueError("did not meet fit acceptance criteria") 55 if return_all: 56 return bestfit, {'inliers':best_inlier_idxs} 57 else: 58 return bestfit 59 60 def random_partition(n,n_data): 61 """return n random rows of data (and also the other len(data)-n rows)""" 62 all_idxs = numpy.arange( n_data ) 63 numpy.random.shuffle(all_idxs) 64 idxs1 = all_idxs[:n] 65 idxs2 = all_idxs[n:] 66 return idxs1, idxs2 67 68 class LinearLeastSquaresModel: 69 """linear system solved using linear least squares 70 71 This class serves as an example that fulfills the model interface 72 needed by the ransac() function. 73 74 """ 75 def __init__(self,input_columns,output_columns,debug=False): 76 self.input_columns = input_columns 77 self.output_columns = output_columns 78 self.debug = debug 79 def fit(self, data): 80 A = numpy.vstack([data[:,i] for i in self.input_columns]).T 81 B = numpy.vstack([data[:,i] for i in self.output_columns]).T 82 x,resids,rank,s = scipy.linalg.lstsq(A,B) 83 return x 84 def get_error( self, data, model): 85 A = numpy.vstack([data[:,i] for i in self.input_columns]).T 86 B = numpy.vstack([data[:,i] for i in self.output_columns]).T 87 B_fit = scipy.dot(A,model) 88 err_per_point = numpy.sum((B-B_fit)**2,axis=1) # sum squared error per row 89 return err_per_point 90 91 def test(): 92 # generate perfect input data 93 n_samples = 500 94 n_inputs = 1 95 n_outputs = 1 96 A_exact = 20*numpy.random.random((n_samples,n_inputs) ) # x坐標 97 perfect_fit = 60*numpy.random.normal(size=(n_inputs,n_outputs) ) # the model(斜率) 98 B_exact = scipy.dot(A_exact,perfect_fit) # y坐標 99 assert B_exact.shape == (n_samples,n_outputs) #驗證y坐標數組的大小 100 #pylab.plot( A_exact, B_exact, 'b.', label='data' ) 101 #pylab.show() 102 103 # add a little gaussian noise (linear least squares alone should handle this well) 104 A_noisy = A_exact + numpy.random.normal(size=A_exact.shape ) # x坐標添加高斯噪聲 105 B_noisy = B_exact + numpy.random.normal(size=B_exact.shape ) # y坐標.... 106 #pylab.plot( A_noisy, B_noisy, 'b.', label='data' ) 107 108 if 1: 109 # add some outliers 110 n_outliers = 100 # 500個數據點有100個是putliers 111 all_idxs = numpy.arange( A_noisy.shape[0] ) 112 numpy.random.shuffle(all_idxs) # 索引隨機排列 113 outlier_idxs = all_idxs[:n_outliers] # 選取all_idxs前100個做outlier_idxs 114 non_outlier_idxs = all_idxs[n_outliers:] # 后面的不是outlier_idxs 115 A_noisy[outlier_idxs] = 20*numpy.random.random((n_outliers,n_inputs) ) # 外點的橫坐標 116 B_noisy[outlier_idxs] = 50*numpy.random.normal(size=(n_outliers,n_outputs) ) # 外點的縱坐標 117 #pylab.plot( A_noisy, B_noisy, 'b.', label='data' ) 118 #pylab.show() 119 120 121 # setup model 122 123 all_data = numpy.hstack( (A_noisy,B_noisy) ) # 組成坐標對 124 input_columns = range(n_inputs) # the first columns of the array 125 output_columns = [n_inputs+i for i in range(n_outputs)] # the last columns of the array 126 debug = False 127 model = LinearLeastSquaresModel(input_columns,output_columns,debug=debug) 128 129 linear_fit,resids,rank,s = scipy.linalg.lstsq(all_data[:,input_columns], 130 all_data[:,output_columns]) 131 132 # run RANSAC algorithm 133 ransac_fit, ransac_data = ransac(all_data,model, 134 50, 1000, 7e3, 300, # misc. parameters 135 debug=debug,return_all=True) 136 if 1: 137 import pylab 138 139 sort_idxs = numpy.argsort(A_exact[:,0]) # 對A_exact排序, sort_idxs為排序索引 140 A_col0_sorted = A_exact[sort_idxs] # maintain as rank-2 array 141 142 if 1: 143 pylab.plot( A_noisy[:,0], B_noisy[:,0], 'k.', label='data' ) 144 pylab.plot( A_noisy[ransac_data['inliers'],0], B_noisy[ransac_data['inliers'],0], 'bx', label='RANSAC data' ) 145 else: 146 pylab.plot( A_noisy[non_outlier_idxs,0], B_noisy[non_outlier_idxs,0], 'k.', label='noisy data' ) 147 pylab.plot( A_noisy[outlier_idxs,0], B_noisy[outlier_idxs,0], 'r.', label='outlier data' ) 148 pylab.plot( A_col0_sorted[:,0], 149 numpy.dot(A_col0_sorted,ransac_fit)[:,0], 150 label='RANSAC fit' ) 151 pylab.plot( A_col0_sorted[:,0], 152 numpy.dot(A_col0_sorted,perfect_fit)[:,0], 153 label='exact system' ) 154 pylab.plot( A_col0_sorted[:,0], 155 numpy.dot(A_col0_sorted,linear_fit)[:,0], 156 label='linear fit' ) 157 pylab.legend() 158 pylab.show() 159 160 if __name__=='__main__': 161 test()上面代碼跟原版的代碼相比,我刪除了一些冗余的東西。在test()中做的是直線擬合。在看test()部分之前,我們先來看看RANSAC部分的代碼,傳入RANSAC函數中的參數有8個,前面6個是比較重要的。data就是全部的數據點集,model注釋里給出的是擬合點集的模型,放到這個直線擬合的實例下,就是斜率,n就是擬合時所需要的最小數據點數目,放在這里直線擬合的例子中,就是用于選取的用于去做直線擬合的數據點數目,k就是最大允許的迭代次數,t是人為設定的用于判斷誤差接受許可的范圍。這幾個參數的含義知道了,剩下的就是理解while循環里面的內容了。在每一次循環中,選對所有的數據點做一個隨機的劃分,將數據點集分成兩堆,分別對應maybeinliers和test_points,maybeinliers這部分數據用于做直線擬合,這里直線擬合采用的是最小二乘法,得到擬合到的直線的斜率maybemodel,然后用該直線及測試數據的橫坐標去估計測試數據的縱坐標,也就是在該模型下測試數據的估計值,測試數據的估計值和測試數據的真實值做一個平方和便得到誤差,將得到的誤差分別和設定的可接受誤差進行判斷,在誤差范圍內的判定為inlier,否者判斷為outlier。當inliers的數目達到了設定的數目的要求是,再講inliers和maybeinliers放一下再做一下最小二乘擬合,便得到最終的最佳斜率了。
test()部分的內容很簡單,先生成在某條直線上的一些離散點,這里某條直線的斜率就是精確的模型:然后添加高斯平穩高斯噪聲:將其中的某些點變為outliers:最后用RANSAC擬合出來的結果如下:整個過程就醬紫,后面有時間繼續前面在BoW圖像檢索Python實戰用RANSAC做一個重排過程。
from:?http://yongyuan.name/blog/fitting-line-with-ransac.html
總結
以上是生活随笔為你收集整理的RANSAC算法做直线拟合的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: VLFeat SLIC超像素分割(Cpp
- 下一篇: GMM、fisher vector、SI