RANSAC算法原理与代码实现

        以前一直把ransac算法误解成了类似于最小二乘的做法，直到今天专门搜了不少资料才搞懂其中的原理，ransac算法是从一组含有外点（也可认为是噪声点）的数据中通过迭代方式来估算出正确的数学模型，通过该模型来保留更多的内点。总的来说是一种基于概率的方式不断取随机抽取数据得到模型，并保留最好的模型，随着迭代次数的增加，得到的模型也最准确（但不一定是最准确的那个）。以直线拟合为例，已知给定一组数据（包含噪声点或叫离群点），具体的做法如下：

1.在一组数据中，随机选取两个点，以此来计算出一根直线参数，然后求剩余点到该直线的距离

2.记录距离小于epsilon的点的个数（小于epsilon的点为内点，epsilon是给定的一个误差值）

3.重新选取两个点，计算直线参数，并计算其余点到直线的距离epsilon的点数，若得到的点数大于历史最大点数，则保留本次模型（保留斜率a和截距b）。

4.重复1～3步，直到迭代次数大于设置的迭代次数或者内点数大于给定阈值，停止迭代。

所以，ransac不易受噪声点的影响，是一种稳健算法，而最小二乘是根据全局误差迭代得到模型，易受到噪声点影响。

使用python由一组含噪声点来拟合直线代码实现如下：

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


##################  生成数据  ################## 
# 数据量
SIZE = 50

# 生成数据 0～10所有元素间隔相等。
X = np.linspace(0, 10, SIZE)
Y = -5 * X + 10 #生成一根直线数据点

fig = plt.figure()
ax1 = fig.add_subplot(1, 1, 1)
# 标题
ax1.set_title("RANSAC")

random_x = []
random_y = []
# 添加直线随机噪声
for i in range(SIZE):
    random_x.append(X[i] + random.uniform(-0.5, 0.5)) 
    random_y.append(Y[i] + random.uniform(-0.5, 0.5)) 

# 添加随机噪声
for i in range(SIZE):
    random_x.append(random.uniform(-10,10))
    random_y.append(random.uniform(-40,40))

RANDOM_X = np.array(random_x) 
RANDOM_Y = np.array(random_y)

# 画散点图
ax1.scatter(RANDOM_X, RANDOM_Y)
ax1.set_xlabel("x")
ax1.set_ylabel("y")


##################  使用RANSAC算法估算模型  ################## 
# 迭代最大次数
iters = 1000
# 数据和模型之间可接受的差值
sigma = 0.25
# 最好模型的参数估计和内点数目
best_a = 0
best_b = 0
best_inner_num = 0
for i in range(iters):
    # 随机在数据中红选出两个点去求解模型
    sample_index = random.sample(range(SIZE * 2),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 * 2):
        y_estimate = a * RANDOM_X[index] + b
        if abs(y_estimate - RANDOM_Y[index]) < sigma:
            total_inlier = total_inlier + 1 #统计内点数

    # 判断内点数是否是历史最好
    if total_inlier > best_inner_num:
        best_inner_num = total_inlier
        best_a = a
        best_b = b

    # 判断是否当前模型已经符合超过一半的点
    if total_inlier > SIZE:
        break

# 用我们得到的最佳估计画图
Y = best_a * RANDOM_X + best_b

# 直线图
ax1.plot(RANDOM_X, Y)
text = "best_a = " + str(best_a) + "\nbest_b = " + str(best_b)
plt.text(4, -40, text,
         fontdict={'size': 8, 'color': 'r'})
plt.show()