《CV中的多視圖幾何》-RANSAC

08-10

《CV中的多視圖幾何》-RANSAC

4 人贊了文章

最近一直在啃《多視圖幾何》，厚厚的一本，不知何年何月能啃完，然而大神們都已經啃了好幾次了，慚愧慚愧。。。。。。

前幾天讀完了RANSAC那一章節，對這個魯棒演算法還是很感興趣的，第一次見到這個演算法還是在高博的《SLAM十四講》中的代碼看到的。。。那時候也不知道這是個什麼玩意，後來在給導師翻譯文章的時候遇到了這個演算法的介紹，算是了解了一點點，直到看完了多視圖幾何中所講的RANSAC，才算是理出了一些頭緒。在知乎上也看了很多大佬的文章，加深了理解，自己也想寫點文章，作為學習筆記，同時也希望能對初學者有一些幫助吧。

RANSAC，隨機採樣一致性，是一種魯棒演算法。我們以直線擬合作為例子，在擬合一條直線時，我們最常用，最先想到的，肯定是最小二乘法了，簡單又實用！簡單地說一下最小二乘法，以擬合直線為例子。

記 $Y=HX$ ，其中， $Y=left[ y_{1}, y_{2},y_{3},...,y_{n} ight]_{1 imes n}$ ， $H=left[ h_{1},h_{2} ight]_{1 imes2}$ ， $X=left[ x_{1} , x_{2} ,..., x_{n} ight]_{2 imes n}$ ，其中， $x_{i}=left[ x^{(i)},1 ight]_{2 imes n}^{T}$ 。

常用的代價函數就是 $L_{2}$ 代價函數：

$J=frac{1}{2}||Y-HX||_{2}^{2}=frac{1}{2}left( Y-HX ight)^{T}( Y-HX )$

很容易算出來，在我的這種矩陣形式的下（也就說，矩陣 $X，Y$ 的定義形式不同，比如把 $X$ 的行數定義為點的總數，列數為特徵的數目，解的形式也會有形式上的差異），解析解為：

$H=YX^{T}(XX^{T})^{-1}$

所以，只要給定數據集，我們套用上面的公式很容易就算出來最小二乘解，但是，當數據中有一些很大的雜訊，導致個別點偏離正確模型非常遠時（將這種點稱為「外點（outliers）」）,會導致我們的最小二乘解會偏向這些外點，從而失去了精度，比如，下面這個圖的情況：

其中，綠線是用最小二乘估計出來的，紅線為真解。

因此，在存在外點的情況下，最小二乘法往往很容易失去精度。這個時候，我們就要採用魯棒代價函數，來儘可能地減小外點的干擾。如Huber函數，還有本文要講的RANSAC。下面，我們就來看什麼是RANSAC。

RANSAC的基本思路：

（1）.從數據中採集最小子集（最小子集的大小由所研究的問題決定，如擬合直線，則最小子集數目為2，因為兩點確定一條直線），然後用這個最小子集擬合一條直線 $l_{i}$ （得到 $ilde{H}$ ）；

（2）.計算每個點到直線 $l_{1}$ 的距離 $d_{ot}$ （在擬合直線的情況下，我採用的是點 $x_{i}$ 通過 $ilde{H}$ 得到的 $ilde{y_{i}}= ilde{H}x_{i}$ 與真值的差的絕對值 $|y_{i}- ilde{y_{i}}|$ ），如果 $d_{ot}<t$ ， $t$ 是距離閾值，則記為「內點（inlier）」。統計本次得到的內點的數目 $n_{i}$ 。如果 $n_{i}>T$ ， $T$ 是內點數目閾值，則說明當前所估計的 $ilde{H}$ 足夠好，利用這些內點再重新估計 $ilde{H}$ ，然後返回 $ilde{H}$ ，退出。否則返回（1）繼續執行。

不停地重複上面兩步，直到估計的 $ilde{H}$ 可以得到最大的內點數，然後利用這些內點再重新估計 $ilde{H}$ 即可。但有一個問題還沒有解決，那就是我們究竟需要重複上面步驟多少次？即採集多少次子集？一種最直觀的想法就是遍曆數據中的所有點，這固然是一個容易想到的方法，但是，如果數據集非常大，則這種方法就是非常費時間的了，效率會變低。為了解決這個問題，給出了這樣的一個公式，用來決定需要採樣次數 $N$ ：

$N=log(1-p)/log(1-(1-varepsilon)^{S})$

其中， $S$ 子集的大小， $p$ 表示由 $S$ 個點組成的隨機樣本中至少有一次沒有外點的概率，通常取 $p=0.99$ ， $varepsilon$ 表示數據集中出現外點的概率（我們可以用數據集中外點所佔數據集的比例來近似）。有了這個公式我們就可以估計出我們需要採樣多少次了。

然而，這就又引出了另一個問題—— $varepsilon$ 我們不知道啊！通常我們並不知道我們的數據集中有多少個外點。所以，這就引出了自適應決定採樣次數版本的RANSAC：

$N=infty$ （或者取100000，總之很大就行），sample_count = 0
當sample_count< $N$

選取一個樣本，估計 $ilde{H}$ ，並計算內點數。若內點數大於閾值，則利用這些內點重新估計 $ilde{H}$ ，然後返回 $ilde{H}$ ，終止。
令 $varepsilon$ = $1-$ （內點數）/（總點數）
取 $p=0.99$ ，由 $N=log(1-p)/log(1-(1-varepsilon)^{S})$ ，求出 $N$
sample_count++

用得到的最多的內點，重新估計 $ilde{H}$ ，然後返回，終止。

以上就是RANSAC的演算法的基本介紹了。總而言之，RANSAC的核心思想在於從帶有外點的數據中儘可能地找出所有的內點，然後利用內點來擬合出最佳的模型。

下面，給出我自己用python寫的RANSAC演算法的實現，以及小例子（本人代碼功底堪憂，沒啥風格，望大佬輕噴/(ㄒoㄒ)/~~）

import numpy as npimport randomimport matplotlib.pyplot as pltimport win_unicode_consolewin_unicode_console.enable()def curve_fitting(X,Y,t): """輸入X:數據集，其中X.shape[0] = 特徵的數目，X.shape[1] = 數據的數目；輸入Y:真值集，其中Y.shape[0] = 特徵的數目，Y.shape[1] = 數據的數目；輸入t:模型估計的值與真值的距離輸出：學習到的參數矩陣 """ T = int(X.shape[1]*0.9) #內點數目的閾值 t_ = t #點到生成的模型的距離閾值 e = 1.0 #外點的比例的初始估計 p = 0.99 #選擇的子集中至少有一次沒有外點的概率 N = 1000000 #選擇次數的初始估計，這個數給的非常大也無所謂，後面很快就降下來了 sample_count = 0 #採樣的次數 num = range(X.shape[1]) #樣本的總數 min_points = X.shape[0] #每次選取的點的數目 X_batch = np.zeros([X.shape[0],min_points]) Y_batch = np.zeros([Y.shape[0],min_points]) inliers = [] #用來存放每個內點子集 index = 0. #inliers中最大內點子集的下標 max_num = 0. #inliers中最大內點子集的數目 ###開始運行ransac while(sample_count <= N): randnum = random.sample(num,min_points) #隨機選出min_points個點 for i in range(min_points): X_batch[:,i] = X[:,randnum[i]] Y_batch[:,i] = Y[:,randnum[i]] A = np.dot(Y_batch,X_batch.T) B = np.dot(X_batch,X_batch.T) C = np.linalg.inv(B) H = np.dot(A,C) #得到的H估計 ###開始尋找內點 inliers_index = [] #用來存放內點的位置 for i in range(X.shape[1]): x = X[:,i].reshape(X.shape[0],-1) dist = abs(np.dot(H,x)-Y[0,i]) #計算模型估計出的值與真值的距離 if dist < t: inliers_index.append(i) print(length: ,len(inliers_index)) if len(inliers_index) > max_num: max_num = len(inliers_index) #更新內點的最大數目 index = sample_count #更新最大內點子集的位置 if len(inliers_index) >= T: #當前採集到的內點子集的數目大於給定閾值，則說明本次估計的H足夠好了，返回H break e = 1.0 - len(inliers_index)/X.shape[1] inliers.append(inliers_index) #將內點子集放入inliers中 N = np.log(1-p+1e-14)/np.log(1-(1-e)**min_points+1e-14) #更新選擇次數N sample_count += 1 print(the amount of inlier_set: , len(inliers)) print(the max number of inlier_set: , len(inliers[index])) ###用得到的含有最多內點的子集重新估計模型 sub_point = inliers[index] X_batch = np.zeros([X.shape[0],len(inliers[index])]) Y_batch = np.zeros([Y.shape[0],len(inliers[index])]) for i in range(len(inliers[index])): X_batch[:,i] = X[:,sub_point[i]] Y_batch[:,i] = Y[:,sub_point[i]] A = np.dot(Y_batch,X_batch.T) B = np.dot(X_batch,X_batch.T) C = np.linalg.inv(B) H = np.dot(A,C) return H if __name__ == __main__: ###構造數據集 X = np.arange(5,20,1).reshape(1,-1) X_ = np.concatenate([X,np.ones([X.shape[0],X.shape[1]])]) #這一步的目的是將形如Y = HX+B整合為形如Y=HX的等價形式 H = np.array([[5,4]]) ###真實的模型 Y = np.dot(H,X_) + np.random.randn(1,X.shape[1]) #Y = HX_ ###構造外點 Y[0,2] += 100 H_estimation = curve_fitting(X_,Y,1.0) print(Estimation: , H_estimation) Y_estimation = np.dot(H_estimation,X_) #用最小二乘直接估計 A = np.dot(Y,X_.T) B = np.dot(X_,X_.T) C = np.linalg.inv(B) H_ = np.dot(A,C) Y_ = np.dot(H_,X_) ###繪製圖像 plt.scatter(X.T,Y.T) #真實值 plt.plot(X.T,Y_estimation.T,r) #RANSAC估計得到的H plt.plot(X.T,Y_.T,g) #最小二乘法得到的H plt.show()

運行結果如下：

補充： $N$ 初始取得非常大，這個不需要在意，只是一個估計，後面非常快地就會變成兩位數。

第一次寫知乎文章，寫的不好，還望見諒。。。。。。希望能對大家有所幫助吧！共同進步O(∩_∩)O