Semi-Global Matching(SGM)演算法原文理解
參考:@迷霧forest http://blog.csdn.net/wsj998689aa/article/details/49464017,原博主對SGM演算法的精髓理解的很透,我是在參考他文章的基礎上,才能看懂SGM演算法幾處關鍵的地方。本文的不同在於加入了一些我自己的理解,並且調整了一下整個演算法闡述的思路,當是自己的一個閱讀筆記。後邊打算再做一下SGM原始演算法與OpenCV的SGBM演算法實現的對比,並且爭取自己實現一下SGM演算法。因為我的導師不搞這個方向,很多東西都是我自己的理解,恐怕會有不對的地方,還請大家多指正。下面直接從第二章開始。
2.半全域性立體匹配
半全域性立體匹配演算法基於一種逐畫素匹配的方法,該方法使用互資訊來評價匹配代價,並通過組合很多一維的約束來近似一個全域性的二維平滑約束。本方法分成下面幾個步驟,其中有一些並不是必須的。
2.1 逐畫素匹配代價的計算
輸入的左右影象必須已經知道對極幾何模型,但是不一定是校正過的,因為有些影象是難以校正的,比如推掃式影象。要計算參考圖(base image)某點P的匹配代價,需要用到其灰度為Ibp,及在待匹配圖(match image)的疑似匹配點q,其灰度為Imq,p和q之間有極線方程q=ebm(p,d),d是極線的引數。如果影象已經被校正了,那麼待匹配圖在參考圖的右邊,且d代表視差。
一個重要的問題是用於匹配的區域的大小和形狀,區域越大,魯棒性越好,但是大了以後會導致物體的邊界模糊。這裡不使用P鄰域內的視差是連續的這一假設,也就是說,只有Ibp和Imq這兩個灰度值被用來計算點p匹配代價,其他的點不考慮。
互資訊M的定義式,互資訊是基於熵來計算的,下面的三項分別是圖1的熵,圖2的熵,以及圖12的聯合熵,緊接著給出熵和聯合熵的定義:
什麼是熵?參考@迷霧forest:先說說熵,熵是用來表徵隨機變數的不確定性(可以理解為變數的資訊量),不確定性越強那麼熵的值越大(最大為1),那麼影象的熵其實就代表影象的資訊量。互資訊度量的是兩個隨機變數之間的相關性,相關性越大,那麼互資訊就越大。可以想想看,兩幅影象如果匹配程度非常高,說明這兩幅影象相關性大還是小?顯然是大,知道一幅影象,另外一幅影象馬上就知道了,相關性已經不能再大了!!!反之,如果兩幅影象配準程度很低,那麼兩幅影象的互資訊就會非常小。所以,立體匹配的目的當然就是互資訊最大化。這就是為什麼使用互資訊的原因。
此時,為了計算互資訊,就要先計算熵,而根據上的定義式,為了計算HI和HI1I2,首先需要知道PI和PI1I2是啥,其中PI代表某個點i的概率分佈,也就是灰度直方圖中灰度為i的點出現的概率。PI1I2是二維概率分佈,二維概率分佈是什麼意思呢,比如現在左圖是100*100,則一共有10000個畫素點,對於左圖中的每一個點p,它在左圖的灰度值為i,它在右圖的對應點的灰度值為j,它們構成一個點對(i,j),此時點對(i,j)對應的計數加1;然後遍歷其他所有點,如果某個點對出現過,就把該點對對應的計數值加1,最後把所有點對的計數值除以10000,就得到某個點對(i,j)的概率。不難理解,因為i和j的範圍都是0-255,所以這個二維概率分佈圖PI1I2的尺寸為256*256。
現在知道了PI和PI1I2,就要返回去計算HI和HI1I2。對於HI1I2,作者根據Kim等人的成果,直接用泰勒展開式,把原來的聯合熵公式轉換為下面這個:
也就是各畫素點對應的h值之和的形式,那麼hI1I2又是什麼呢?hI1I2是根據PI1I2來計算的,其公式是:
也就是說,得到二維概率分佈圖(其尺寸永遠是256*256),然後與高斯核進行卷積,然後求log,然後再卷積,就得到這個hI1I2,其實它也是個256*256的圖,這樣就可以建一個256*256的表,然後HI1I2的解釋就是,對於逐個點p,其左影象中灰度值為I1P,右影象灰度值為I2p,則其對應的hI1I2可以直接查表得出,然後整個圖的HI1I2就是把所有點p的hI1I2加起來。
這裡列出由兩個圖得到hI1I2的過程,這裡有一個warp,它的意思是,在求二維概率分佈之前,要根據一個初始的視差圖,對待匹配影象進行修正,修正後,左圖中某一點大概率能有一個正確的對應點,則此時兩個點的灰度值基本相同,這樣就使得(5,5)(10,10)這樣的點對佔大多數,所以P的影象基本是一個45度的直線。但是這個作為參考的視差圖是粗糙的,它是怎麼得到的作者沒有提到,只知道遮擋點的視差都被設為0,只有非遮擋點才有視差,用這樣粗糙的視差圖作為參考去對待匹配影象做修正,就可能使對應點的視差不一致,因此就會出現不在45°線上的情況,這就需要高斯模糊來去躁了。至此,就討論完了HI1I2的計算。還有HI1和HI2。
對於HI,其計算方法本來跟HI1I2是不同的。但是這裡有個問題,如果考慮遮擋的問題,那麼有些點根本沒有可以匹配的點,那麼這些點不應該被考慮,為此,我們又一次想到聯合熵,聯合熵是基於聯合概率分佈的,其基於的點都可以保證是匹配點,因此,這裡考慮對於HI的計算,也使用HI1I2即聯合熵的求法。公式如下:
這裡hI的計算就跟上面很相似了,需要注意的是計算P的時候不要把遮擋點統計進來。
最後,互資訊MI的最終定義就是:
相應的,基於互資訊的某點p的匹配代價CMI就是:
這就是說,計算某一點的匹配代價,就是先根據極線方程求出在右影象中的灰度,根據這兩個點的灰度構成的灰度值對直接在互資訊圖M中查詢,然後取負。
這裡還有一個問題,求互資訊之前要對待匹配影象進行校正,根據Kim等,採用一個迭代的過程,一開始採用一個隨機的視差圖來進行第一次的互資訊計算,然後根據這個互資訊就可以得到新的視差圖,再根據這個視差圖再求互資訊,如此迭代,一般也不要太多次數,比如3次可以。為什麼可以這麼“隨意”,因為畫素點很多,所以即使隨機生成的視差圖有錯誤,體現在概率分佈上也沒有特別大的影響。
最後就是分層互資訊。既然一個粗糙的初始視差圖也可以得到較精準的概率分佈,那麼在前面的迭代中可以使用一個基於相關的快速方法,並且只在最後一次迭代時使用更精準和耗時的方法。下面就介紹這種方法,也就是分層互資訊法(HMI)。
用HMI來計算匹配代價,該方法遞迴的使用降取樣過的參考圖和待匹配圖,取樣率是1/2,然後得到的視差圖的長寬也是上一幅視差圖的一半(注意,因為我們只對影象降取樣,而影象上的點仍然是可能佔滿0-255空間的,因此概率分佈圖仍然是256*256)。那麼如果迭代四次,第一次使用的是1/16的圖,下一次是1/8,再1/4和1/2,最後一次是1,這就得到較高精度的匹配代價。單次迭代的時間複雜度是O(WHD),所以總的運算時間是 ,可見覆雜度並不是很高,之比BT演算法多了14%。注意低一層的圖只作為下一層的輸入,而不作為最後結果的輸入。上面的公式中為啥要乘3,根據@迷霧forest,乘以3的原因是隨機生成的視差圖十分不靠譜,需要反覆迭代3次才能得到同樣解析度下的靠譜視差圖,然後再參與後續高解析度的計算。
2.2 代價聚合
由2.1已經得到基於互資訊的匹配代價,但是這個逐畫素的匹配代價容易受到誤匹配和噪聲點的影響,因此這裡還要考慮使用某點的鄰域視差資料來構造懲罰函式以增加平滑性約束。也就是說,把代價聚合裡面加入平滑性約束的考慮,也就是設計一個全域性的能量函式E(D),用它來綜合匹配代價和平滑性約束。能量函式的定義式如下:
(注意這裡懲罰值的意義,懲罰就代表著不想選這種情況,如果有一個視差圖它的懲罰值比另一幅大,那麼就不選這個,也就是通過懲罰篩掉了我們不想要的情況,這個彎要轉過來)
式中第一項是基於互資訊的代價計算項;後兩項是當前畫素p和其鄰域Np內所有畫素q之間的約束,如果p和q的視差只差1,那麼懲罰P1,否則懲罰係數P2,P2要大於P1。如果當前點和鄰域點的差值比較小,就選擇一個小的懲罰值,這麼做保證可以對斜面和曲面有一定的適應性(也就是說不是完全要篩掉它們,而是還有容忍)。
現在問題就變成找到一個視差圖,它能使這個全域性的能量函式最小。關於求這個視差圖的方法,@迷霧forest文章說的講比較清楚:這個時候問題來了,這個E對p是不可導的,這意味著我們常看到的梯度下降,牛頓高斯等等演算法在這裡都不適用,作者於是採用了動態規劃來解決這一問題。簡單地說,p的代價想要最小,那麼前提必須是鄰域內的點q的代價最小,q想要代價最小,那麼必須保證q的領域點m的代價最小,如此傳遞下去。怎麼利用動態規劃來求解E,其實這個求解問題是NP完全問題,想在2D影象上直接利用動態規劃求解是不可能的,只有沿著每一行或者每一列求解才能夠滿足多項式時間,但是這裡問題來了,如果我們只沿著每一行求解,那麼行間的約束完全考慮不到,q是p的領域的點其實這個時候被弱化到了q是p的左側點或者右側點,這樣的求優效果肯定很差。於是,大招來了!!我們索性不要只沿著橫或者縱來進行優化,而是沿著一圈8個或者16個方向進行優化。
什麼叫“沿著一圈8個或者16個方向進行優化”呢,結合下面的執行步驟詳細解釋。現在可以採用一種新的匹配代價聚合方法,它使用某點p上鄰域內所有方向上的代價構成一維的DP問題。S(p,d)代表某一點(畫素點名p,視差值d)的代價聚合值,它的求法是,所以Lr是什麼?以p為中心,可以找到16個(或8個,下同)鄰域方向,然後在每個鄰域方向上有一個指向p的向量,而r表示16個向量中的某一個。因此Lr就是代表當前向量路徑r上的代價聚合值。所以,S(p,d)就是把16個代價聚合值相加,而這16個又都是各自路徑上的最小代價。所以求出來的S(p,d)就是p點的最小代價。注意r是一條路,它只能是這個方向,但是長度可以延長。也就是說,16個方向分別算到底。那麼下面就是考慮,Lr(p,d)到底如何計算呢。且看下式:
P點上某個鄰域方向的代價聚合值分為三項,第一是當前點的匹配代價,就是前面的CMI;第二項是min(當前鄰域方向上p-r這個畫素點的當前視差代價聚合值,p-r點的視差差值為1的代價聚合值 + P1,p-r點的視差插值大於1的最小代價聚合值 + P2);最後一項是p-r點的視差插值大於1的最小代價聚合值,最後一項是為了防止計算結果太大,做的統一的調整。
關於這個p-r,它的意義不是很好理解,我的理解的是,因為每次要沿著一個向量方向算到底,那麼p-r就是某個向量上距離p為一個單位、兩個單位等等點,直到算完(addsthe lowest cost of the previous pixel p-r ofthe path)。
把各個方向上的Lr都求和得到S(p,d),作者建議方向數最好用16,至少為8,這些方向最好是水平、垂直或位於某對角線上,如果不是,那就多做一步,調整為這三種。
最後考慮一下計算量和實現方法。作者推薦的是先計算所有點的C(p,d),然後把值規範化到0-211的取值範圍內。所有的匹配代價都存在C[]數組裡面,該陣列大小W*H*D,然後用一個與C同樣大小的陣列S[]來儲存代價聚合值。然後對每個點b,計算所有鄰域向量方向上的代價聚合值,計算過程中體現動態規劃中記憶化搜尋的思想。
2.5 視差求精
後處理這一塊我感覺http://blog.csdn.net/zhubaohua_bupt/article/details/51866567所做的介紹比本文要更詳細,因為視差求精本來就是精益求精嘛,所以越全面參考價值越大。
2.6 大尺寸圖片的處理
SGM演算法需要儲存臨時變數包括每個點的匹配代價、聚合代價、融合前的視差圖等等,因此中間步驟需要的儲存空間非常大,很可能導致空間佔滿的情況。為此,提出一種方法,把基準圖分成片,以片的形式參與到2.1-2.3步的運算,然後再組合起來進行2.4的運算。劃分片的時候,它們需要是重疊的,也就是每個片取大一點。這種方法對於大尺寸的照片非常適用。
2.7 視差圖融合
注意這裡不是左右視差圖的融合,而是針對不同視角獲取到的視差圖進行融合。我的應用只用了一個視角,所以這裡沒怎麼看。