Ransac 算法python 应用和实现

Ransac 算法是一种常用的图像匹配算法，在参数估计领域也经常被使用到。针对估计各种曲线的鲁棒模型参数，效果显著。这里对ransac算法进行某些探索。

python program:

import numpy as np
import matplotlib.pyplot as plt
import random
import math

# 数据量。
SIZE = 60
SIZE_N = 10 # the numbe of noise
# 产生数据。np.linspace 返回一个一维数组，SIZE指定数组长度。
# 数组最小值是0，最大值是10。所有元素间隔相等。
X = np.linspace(0, 10, SIZE)
Y = -2 * X + 5

fig = plt.figure()
# 画图区域分成1行1列。选择第一块区域。
ax1 = fig.add_subplot(111)
# 标题
ax1.set_title("title ")


# 让散点图的数据更加随机并且添加一些噪声。
random_x = []
random_y = []

random_x2 = []
random_y2 = []

random_x2b = []
random_y2b = []

random_x22 = []
random_y22 = []

random_x22b = []
random_y22b = []
# 添加直线随机噪声
for i in range(SIZE):
    random_x.append(X[i] + random.uniform(-1, 1)) 
    random_y.append(Y[i] + random.uniform(-1, 1)) 
# 添加随机噪声
for i in range(SIZE_N):
    random_x.append(random.uniform(-SIZE,SIZE))
    random_y.append(random.uniform(-SIZE,SIZE))
RANDOM_X = np.array(random_x) # 散点图的横轴。
RANDOM_Y = np.array(random_y) # 散点图的纵轴。



# 使用RANSAC算法估算模型
# 迭代最大次数，每次得到更好的估计会优化iters的数值
iters = 1000
iters2 = int(iters/2)
# 数据和模型之间可接受的差值
sigma = 3
sigma2 = 10
# 最好模型的参数估计和内点数目
best_a = 0
best_b = 0
best_a2 = 0
best_b2 = 0
pretotal = 0
pretotal2 = 0
# 希望的得到正确模型的概率
P = 0.99

for i in range(iters):
    # update the record position for seconde RANSAC 
    random_x2 = []
    random_y2 = []
    # 随机在数据中红选出两个点去求解模型
    sample_index = random.sample(range(SIZE + SIZE_N),2)
    x_1 = RANDOM_X[sample_index[0]]
    x_2 = RANDOM_X[sample_index[1]]
    y_1 = RANDOM_Y[sample_index[0]]
    y_2 = RANDOM_Y[sample_index[1]]

    # y = ax + b 求解出a，b
    a = (y_2 - y_1) / (x_2 - x_1)
    b = y_1 - a * x_1

    # 算出内点数目
    total_inlier = 0
    for index in range(SIZE + SIZE_N): # SIZE * 2 is because add 2 times noise of SIZE
        y_estimate = a * RANDOM_X[index] + b
        if abs(y_estimate - RANDOM_Y[index]) < sigma:
            total_inlier = total_inlier + 1
            # record these points that between +-sigma
            random_x2.append(RANDOM_X[index])
            random_y2.append(RANDOM_Y[index])

    # 判断当前的模型是否比之前估算的模型好
    if total_inlier > pretotal:
        iters = math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE + SIZE_N), 2))
        pretotal = total_inlier
        best_a = a
        best_b = b
        # update the latest better points
        random_x2b = np.array(pretotal) # 散点图的横轴。
        random_y2b = np.array(pretotal) # 散点图的纵轴。
        random_x2b = random_x2
        random_y2b = random_y2
        SIZE2 = pretotal
 
    # 判断是否当前模型已经超过八成的点
    if total_inlier > 0.8 * SIZE:
        break

# 用我们得到的最佳估计画图
# 横轴名称。
ax1.set_xlabel("top view x-axis")
# 纵轴名称。
ax1.set_ylabel("top view y-axis")

Y = best_a * RANDOM_X + best_b

# show the ransac2 points:
ax1.scatter(random_x2b, random_y2b, c='b', marker='v')

# 直线图
ax1.scatter(RANDOM_X, RANDOM_Y, c='r', marker='^')

ax1.plot(RANDOM_X, Y, c='b',)
text = "best_a = " + str(best_a) + "\nbest_b = " + str(best_b)
plt.text(5,50, text,
         fontdict={'size': 12, 'color': 'b'})


# the seconde ransac call the point that cover the largest area
RANDOM_XX = np.array(random_x2b) # 散点图的横轴。
RANDOM_YY = np.array(random_y2b) # 散点图的纵轴。

for i in range(iters2):
    random_x22 = []
    random_y22 = []
    # 随机在数据中红选出一个点去求解模型
    sample_index2 = random.sample(range(SIZE2),1)
    x_12 = RANDOM_XX[sample_index2[0]]
    y_12 = RANDOM_YY[sample_index2[0]]


    # y = ax + b 求解出a，b
    a2 = -1 / a
    b2 = y_12 - (a2 * x_12)

    # 算出内点数目
    total_inlier2 = 0
    for index in range(SIZE2):    # SIZE * 2 is because add 2 times noise of SIZE
        y_estimate2 = a2 * RANDOM_XX[index] + b2
        if abs(y_estimate2 - RANDOM_YY[index]) < sigma2:
            total_inlier2 = total_inlier2 + 1
            # record these points that between +-sigma
            random_x22.append(RANDOM_XX[index])
            random_y22.append(RANDOM_YY[index])
            

    # 判断当前的模型是否比之前估算的模型好
    if total_inlier2 > pretotal2:
        print("total_inlier2:", total_inlier2)
        print("SIZE2:", SIZE2)
        iters = math.log(1 - P) / math.log(1 - pow(total_inlier2 / SIZE2, 2))
        pretotal2 = total_inlier2
        best_a2 = a2
        best_b2 = b2
        
        # update the latest better points
        random_x22b = np.array(pretotal2) # 散点图的横轴。
        random_y22b = np.array(pretotal2) # 散点图的纵轴。
        random_x22b = random_x22
        random_y22b = random_y22
 
    # 判断是否当前模型已经超过八成的点
    if total_inlier2 > 0.8 * SIZE2:
        break
    
# 用我们得到的最佳估计画图
YY = best_a2 * RANDOM_XX + best_b2

# show the ransac2 points:
ax1.scatter(random_x22b, random_y22b, c='g', marker='o')

ax1.set_aspect('equal', adjustable='box')
# 直线图
ax1.plot(RANDOM_XX, YY, c='g' )
text = "best_a2 = " + str(best_a2) + "\nbest_b2 = " + str(best_b2)
plt.text(1,30, text,
         fontdict={'size': 12, 'color': 'g'})
plt.show()

ptyhon results:

References：

ransac实现参考：
scatter()使用方法
Matplotlib 绘制等轴正方形图
random.uniform( ) 函数教程与实例