RANSAC算法:随机抽样一致算法(random sample consensus,RANSAC)
一个简单的例子是从一组观测数据中找出合适的二维直线。假设观测数据中包含局内点和局外点,其中局内点近似的被直线所通过,而局外点远离于直线。简单的最小二乘法不能找到适应于局内点的直线,原因是最小二乘法尽量去适应包括局外点在内的所有点。相反,RANSAC能得出一个仅仅用局内点计算出模型,并且概率还足够高。但是,RANSAC并不能保证结果一定正确,为了保证算法有足够高的合理概率,我们必须小心的选择算法的参数。
RANSAC通过反复选择数据中的一组随机子集来达成目标。被选取的子集被假设为局内点,并用下述方法进行验证:
1.有一个模型适应于假设的局内点,即所有的未知参数都能从假设的局内点计算得出。
2.用1中得到的模型去测试所有的其它数据,如果某个点适用于估计的模型,认为它也是局内点。
3.如果有足够多的点被归类为假设的局内点,那么估计的模型就足够合理。
4.然后,用所有假设的局内点去重新估计模型,因为它仅仅被初始的假设局内点估计过。
5.最后,通过估计局内点与模型的错误率来评估模型。
这个过程被重复执行固定的次数,每次产生的模型要么因为局内点太少而被舍弃,要么因为比现有的模型更好而被选用。
下面用一张图来模拟一下数据拟合的过程。
与最小二乘法的区别:
RANSA通过多次迭代,计算一个模型包含的局内点,模型好坏的判别依据是:模型包含越多的局内点,则这个模型越好。
而最小二乘法则是计算出一条直线,让所有的样本距离这条直线的距离最短。所以判别模型的标准就是样本距离直线的距离和。
优缺点:
RANSAC 算法的优点是能鲁棒的估计模型参数。例如,他能从包含大量局外点的数据集中估计出高精度的参数。
缺点是它计算参数的迭代次数没有上限,如果设置迭代次数的上限,得到的结果可能不是最优的结果,甚至可能得到错误的结果。
RANSAC只有一定的概率得到的可信的模型,概率与迭代次数成正比。另一个缺点是它要求设置跟问题相关的阈值,
RANSAC职能从特定的数据集中估计出一个模型,如果存在两个(或多个)模型,RANSAC不能找到别的模型。
输入:
data —— 一组观测数据
model —— 适应于数据的模型
n —— 适用于模型的最少数据个数
k —— 算法的迭代次数
t —— 用于决定数据是否适应于模型的阀值
d —— 判定模型是否适用于数据集的数据数目
输出:
best_model —— 跟数据最匹配的模型参数(如果没有找到好的模型,返回null)
best_consensus_set —— 估计出模型的数据点
best_error —— 跟数据相关的估计出的模型错误
附上一个可以实现的代码:
# -*- coding: utf-8 -*-
import numpy
import scipy # use numpy if scipy unavailable
import scipy.linalg # use numpy if scipy unavailable
import pylab
def ransac(data, model, n, k, t, d, debug=False, return_all=False):
'''fit model parameters to data using the RANSAC algorithm
This implementation written from pseudocode found at
Given:
data - a set of observed data points # 可观测数据点集
model - a model that can be fitted to data points #
n - the minimum number of data values required to fit the model # 拟合模型所需的最小数据点数目
k - the maximum number of iterations allowed in the algorithm # 最大允许迭代次数
t - a threshold value for determining when a data point fits a model #确认某一数据点是否符合模型的阈值
d - the number of close data values required to assert that a model fits well to data
Return:
bestfit - model parameters which best fit the data (or nil if no good model is found)
'''
iterations = 0
bestfit = None
besterr = numpy.inf
best_inlier_idxs = None
while iterations < k:
maybe_idxs, test_idxs = random_partition(n, data.shape[0])
maybeinliers = data[maybe_idxs, :]
test_points = data[test_idxs]
maybemodel = model.fit(maybeinliers)
test_err = model.get_error(test_points, maybemodel)
also_idxs = test_idxs[test_err < t] # select indices of rows with accepted points
alsoinliers = data[also_idxs, :]
if debug:
print
'test_err.min()', test_err.min()
print
'test_err.max()', test_err.max()
print
'numpy.mean(test_err)', numpy.mean(test_err)
print
'iteration %d:len(alsoinliers) = %d' % (
iterations, len(alsoinliers))
if len(alsoinliers) > d:
betterdata = numpy.concatenate((maybeinliers, alsoinliers))
bettermodel = model.fit(betterdata)
better_errs = model.get_error(betterdata, bettermodel)
thiserr = numpy.mean(better_errs)
if thiserr < besterr:
bestfit = bettermodel
besterr = thiserr
best_inlier_idxs = numpy.concatenate((maybe_idxs, also_idxs))
iterations += 1
if bestfit is None:
raise ValueError("did not meet fit acceptance criteria")
if return_all:
return bestfit, {'inliers': best_inlier_idxs}
else:
return bestfit
def random_partition(n, n_data):
#返回n行数据
all_idxs = numpy.arange(n_data)
numpy.random.shuffle(all_idxs)
idxs1 = all_idxs[:n]
idxs2 = all_idxs[n:]
return idxs1, idxs2
#线性系统可以使用最小二乘法求解,这个类是定义最小二乘法求解的一个示例
class LinearLeastSquaresModel:
def __init__(self, input_columns, output_columns, debug=False):
self.input_columns = input_columns
self.output_columns = output_columns
self.debug = debug
def fit(self, data):
A = numpy.vstack([data[:, i] for i in self.input_columns]).T
B = numpy.vstack([data[:, i] for i in self.output_columns]).T
x, resids, rank, s = scipy.linalg.lstsq(A, B)
return x
def get_error(self, data, model):
A = numpy.vstack([data[:, i] for i in self.input_columns]).T
B = numpy.vstack([data[:, i] for i in self.output_columns]).T
B_fit = scipy.dot(A, model)
err_per_point = numpy.sum((B - B_fit) ** 2, axis=1) # sum squared error per row
return err_per_point
def test():
#创建一些数据(局内点)使其在一条直线附近
n_samples = 500
n_inputs = 1
n_outputs = 1
A_exact = 20 * numpy.random.random((n_samples, n_inputs)) # x坐标
perfect_fit = 60 * numpy.random.normal(size=(n_inputs, n_outputs)) # the model(斜率)
B_exact = scipy.dot(A_exact, perfect_fit) # y坐标
#验证y坐标数组的大小,在数组计算中尺寸大小很重要,如果尺寸大小不匹配,则导致计算错误,或者无法计算
assert B_exact.shape == (n_samples, n_outputs)
#pylab.plot( A_exact, B_exact, 'b.', label='data' ) #这两行是图片显示创建出来的数据,不含噪声的数据
#pylab.show()
#添加一点点高斯噪声(仅用最小二乘的方法即可很好的去解决小噪声的线性回归问题)
A_noisy = A_exact + numpy.random.normal(size=A_exact.shape) # x坐标添加高斯噪声
B_noisy = B_exact + numpy.random.normal(size=B_exact.shape) # y坐标....
#pylab.plot( A_noisy, B_noisy, 'b.', label='data') #这两句是图片显示带有噪声的数据。
#pylab.show()
#添加局外点,对模型进行干扰
n_outliers = 100 # 500个数据点有100个是putliers
all_idxs = numpy.arange(A_noisy.shape[0])
numpy.random.shuffle(all_idxs) # 索引随机排列
outlier_idxs = all_idxs[:n_outliers] # 选取all_idxs前100个做outlier_idxs
non_outlier_idxs = all_idxs[n_outliers:] # 后面的不是outlier_idxs
A_noisy[outlier_idxs] = 20 * numpy.random.random((n_outliers, n_inputs)) # 外点的横坐标
B_noisy[outlier_idxs] = 50 * numpy.random.normal(size=(n_outliers, n_outputs)) # 外点的纵坐标
pylab.plot( A_noisy, B_noisy, 'b.', label='data' )
pylab.show()
#建立模型
all_data = numpy.hstack((A_noisy, B_noisy)) # 组成坐标对
input_columns = range(n_inputs) # the first columns of the array
output_columns = [n_inputs + i for i in range(n_outputs)] # the last columns of the array
debug = False
model = LinearLeastSquaresModel(input_columns, output_columns, debug=debug)
linear_fit, resids, rank, s = scipy.linalg.lstsq(all_data[:, input_columns],all_data[:, output_columns])
#使用RANSAC算法
ransac_fit, ransac_data = ransac(all_data, model,50, 1000, 7000, 300, # misc. parameters
debug=debug, return_all=True)
sort_idxs = numpy.argsort(A_exact[:, 0]) # 对A_exact排序, sort_idxs为排序索引
A_col0_sorted = A_exact[sort_idxs] # maintain as rank-2 array
pylab.plot(A_noisy[:, 0], B_noisy[:, 0], 'k.', label='data')
pylab.plot(A_noisy[ransac_data['inliers'], 0], B_noisy[ransac_data['inliers'], 0], 'bx',label='RANSAC data')
pylab.plot(A_col0_sorted[:, 0],numpy.dot(A_col0_sorted, ransac_fit)[:, 0] ,label='RANSAC fit')
pylab.plot(A_col0_sorted[:, 0],numpy.dot(A_col0_sorted, perfect_fit)[:, 0],label='exact system')
pylab.plot(A_col0_sorted[:, 0],numpy.dot(A_col0_sorted, linear_fit)[:, 0] ,label='linear fit')
pylab.legend()
pylab.show()
#运行主函数
if __name__ == '__main__':
test()
这是我自己创建的数据,数据点多的那一大块是我创建的局内点,下面那一小排是我创建的局外点。
实验结果如下图所示:
应用:
RANSAC可以实现特征点检测,通过特征点匹配去辨识人脸。达到一个人脸识别的效果。
同样,通过特征点匹配也可以对多张图进行一个图像融合,比如全景照片的拍摄,用的就是图像融合技术。
参考博客:https://blog.csdn.net/laobai1015/article/details/51682596
https://blog.csdn.net/robinhjwy/article/details/79174914
https://blog.csdn.net/fandq1223/article/details/53175964