Scipy Cookbook中RANSAC实现的不一致性

Question

最近我查看了Cookbook中的RANSAC实现：http://wiki.scipy.org/Cookbook/RANSAC，但发现它似乎和RANSAC算法本身不太一致。

看看那里的图表，为什么有些数据点离最佳模型很远（见底部的点），却被认为是“RANSAC数据”，而有些离模型更近的点却不是呢？

在我看来，这和RANSAC算法的主要思想相矛盾，因为在这个算法中，所有在预设阈值范围内的点都应该被视为有效点。

为什么在这个实现中不是这样？还有其他Python中的RANSAC实现吗？

谢谢你的帮助！

祝好，
Alexey

Answer 1

这并不与RANSAC的理念相矛盾。不过，图表有点误导人。

图中显示的蓝色十字是样本点（best_inlier_idxs = maybeinliers + alsoinliers），其中一些点（也就是alsoinliers中的点，称为共识集）是支持模型的（maybemodel），这个模型是通过随机选择的数据样本（maybeinliers）来拟合的。这意味着，alsoinliers中的所有点应该比那些不支持这个模型的最近点更接近maybemodel。

但是，图中并没有显示maybemodel的拟合结果。图中显示的是bettermodel（蓝线，“RANSAC拟合”），这是通过将模型参数拟合到best_inlier_idxs中的所有点（而不仅仅是maybeinliers中的点）得到的。

此外，best_inlier_idxs包含了alsoinliers和maybeinliers。在随机选择的样本maybeinliers中，可能有一些点实际上并不支持maybemodel的拟合（也就是说，它们与模型的距离超过了设定的阈值）。这些点也显示为蓝色十字，尽管它们比其他不在支持集中的点更远。

我稍微修改了一下绘图，以便在“RANSAC数据”中也指示出最佳建议模型（maybemodel）和随机样本（maybeinliers）。重要的是，某些十字周围的圆圈强调了随机样本包含在RANSAC数据中的事实。

以下是修改后的绘图代码：

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) )
    best_maybe_model = maybemodel
    best_random_set = maybe_idxs
    iterations+=1
if bestfit is None:
    raise ValueError("did not meet fit acceptance criteria")
if return_all:
    return bestfit, {'inliers':best_inlier_idxs, 'best_random_set':best_random_set,'best_maybe_model':best_maybe_model}
else:
    return bestfit

def test():
    # generate perfect input data

    n_samples = 500
    n_inputs = 1
    n_outputs = 1
    A_exact = 20*numpy.random.random((n_samples,n_inputs) )
    perfect_fit = 60*numpy.random.normal(size=(n_inputs,n_outputs) ) # the model
    B_exact = scipy.dot(A_exact,perfect_fit)
    assert B_exact.shape == (n_samples,n_outputs)

    # add a little gaussian noise (linear least squares alone should handle this well)
    A_noisy = A_exact + numpy.random.normal(size=A_exact.shape )
    B_noisy = B_exact + numpy.random.normal(size=B_exact.shape )

    if 1:
        # add some outliers
        n_outliers = 100
        all_idxs = numpy.arange( A_noisy.shape[0] )
        numpy.random.shuffle(all_idxs)
        outlier_idxs = all_idxs[:n_outliers]
        non_outlier_idxs = all_idxs[n_outliers:]
        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) )

    # setup model

    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 = True
    model = LinearLeastSquaresModel(input_columns,output_columns,debug=debug)

    linear_fit,resids,rank,s = scipy.linalg.lstsq(all_data[:,input_columns],
                                                  all_data[:,output_columns])

    # run RANSAC algorithm
    ransac_fit, ransac_data = ransac(all_data,model,
                                     50, 1000, 7e3, 300, # misc. parameters
                                     debug=debug,return_all=True)
    if 1:
        import pylab

        sort_idxs = numpy.argsort(A_exact[:,0])
        A_col0_sorted = A_exact[sort_idxs] # maintain as rank-2 array

        if 1:
            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_noisy[ransac_data['best_random_set'],0], B_noisy[ransac_data['best_random_set'],0], 'ro', mfc='none',label='best random set (maybeinliers)' )
        else:
            pylab.plot( A_noisy[non_outlier_idxs,0], B_noisy[non_outlier_idxs,0], 'k.', label='noisy data' )
            pylab.plot( A_noisy[outlier_idxs,0], B_noisy[outlier_idxs,0], 'r.', label='outlier 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.plot( A_col0_sorted[:,0],
                    numpy.dot(A_col0_sorted,ransac_data['best_maybe_model'])[:,0],
                    label='best proposed model (maybemodel)' )
        pylab.legend()
        pylab.show()

if __name__=='__main__':
    test()