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()