Ransac 算法的探索和应用

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 + 5fig = 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.99for 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，ba = (y_2 - y_1) / (x_2 - x_1)b = y_1 - a * x_1# 算出内点数目total_inlier = 0for index in range(SIZE + SIZE_N): # SIZE * 2 is because add 2 times noise of SIZEy_estimate = a * RANDOM_X[index] + bif abs(y_estimate - RANDOM_Y[index]) < sigma:total_inlier = total_inlier + 1# record these points that between +-sigmarandom_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_inlierbest_a = abest_b = b# update the latest better pointsrandom_x2b = np.array(pretotal) # 散点图的横轴。random_y2b = np.array(pretotal) # 散点图的纵轴。random_x2b = random_x2random_y2b = random_y2SIZE2 = 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，ba2 = -1 / ab2 = y_12 - (a2 * x_12)# 算出内点数目total_inlier2 = 0for index in range(SIZE2):    # SIZE * 2 is because add 2 times noise of SIZEy_estimate2 = a2 * RANDOM_XX[index] + b2if abs(y_estimate2 - RANDOM_YY[index]) < sigma2:total_inlier2 = total_inlier2 + 1# record these points that between +-sigmarandom_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_inlier2best_a2 = a2best_b2 = b2# update the latest better pointsrandom_x22b = np.array(pretotal2) # 散点图的横轴。random_y22b = np.array(pretotal2) # 散点图的纵轴。random_x22b = random_x22random_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()