之前隻是簡單了解RANSAC模型，知道它是幹什麼的。然後今天有個課程設計的報告，上去講了一下RANSAC，感覺這個東西也沒那麼複雜，是以今天就總結一些RASAC并用Python實作一下直線拟合。
RANSAC簡介
RANSAC(RAndom SAmple Consensus,随機采樣一緻)算法是從一組含有“外點”(outliers)的資料中正确估計數學模型參數的疊代算法。“外點”一般指的的資料中的噪聲，比如說比對中的誤比對和估計曲線中的離群點。是以，RANSAC也是一種“外點”檢測算法。RANSAC算法是一種不确定算法，它隻能在一種機率下産生結果，并且這個機率會随着疊代次數的增加而加大（之後會解釋為什麼這個算法是這樣的）。RANSAC算最早是由Fischler和Bolles在SRI上提出用來解決LDP(Location Determination Proble)問題的。
對于RANSAC算法來說一個基本的假設就是資料是由“内點”和“外點”組成的。“内點”就是組成模型參數的資料，“外點”就是不适合模型的資料。同時RANSAC假設：在給定一組含有少部分“内點”的資料，存在一個程式可以估計出符合“内點”的模型。
以上文字翻譯于wiki，并加有一些部落客自己的了解。
算法基本思想和流程
RANSAC是通過反複選擇資料集去估計出模型，一直疊代到估計出認為比較好的模型。
具體的實作步驟可以分為以下幾步：
1.選擇出可以估計出模型的最小資料集；(對于直線拟合來說就是兩個點，對于計算Homography矩陣就是4個點)
2.使用這個資料集來計算出資料模型；
3.将所有資料帶入這個模型，計算出“内點”的數目；(累加在一定誤差範圍内的适合目前疊代推出模型的資料)
4.比較目前模型和之前推出的最好的模型的“内點“的數量，記錄最大“内點”數的模型參數和“内點”數；
5.重複1-4步，直到疊代結束或者目前模型已經足夠好了(“内點數目大于一定數量”)。
疊代次數推導
這裡有一點就是疊代的次數我們應該選擇多大呢？這個值是否可以事先知道應該設為多少呢？還是隻能憑經驗決定呢？
這個值其實是可以估算出來的。下面我們就來推算一下。
假設“内點”在資料中的占比為t
那麼我們每次計算模型使用N個點的情況下，選取的點至少有一個外點的情況就是
也就是說，在疊代k次的情況下，(1−tn)k就是k次疊代計算模型都至少采樣到一個“外點”去計算模型的機率。那麼能采樣到正确的N個點去計算出正确模型的機率就是
通過上式，可以求得
“内點”的機率t通常是一個先驗值。然後P是我們希望RANSAC得到正确模型的機率。如果事先不知道t的值，可以使用自适應疊代次數的方法。也就是一開始設定一個無窮大的疊代次數，然後每次更新模型參數估計的時候，用目前的“内點”比值來當成是t估算出疊代次數。
用Python實作直線拟合
import numpy as np
import matplotlib.pyplot as plt
import random
import math
# 資料量。
SIZE = 50
# 産生資料。np.linspace 傳回一個一維數組，SIZE指定數組長度。
# 數組最小值是0，最大值是10。所有元素間隔相等。
X = np.linspace(0, 10, SIZE)
Y = 3 * X + 10
fig = plt.figure()
# 畫圖區域分成1行1列。選擇第一塊區域。
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(0,10))
random_y.append(random.uniform(10,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的數值
iters = 100000
# 資料和模型之間可接受的內插補點
sigma = 0.25
# 最好模型的參數估計和内點數目
best_a = 0
best_b = 0
pretotal = 0
# 希望的得到正确模型的機率
P = 0.99
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 > pretotal:
iters = math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE * 2), 2))
pretotal = 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(5,10, text,
fontdict={'size': 8, 'color': 'r'})
plt.show()
最後得到結果如下圖：