深度濾波器詳細解讀
極品巧克力
前言
空間定位是VR\AR中的一項關鍵技術，計算機創建出來的虛擬圖像必須要依賴空間定位技術才能與現實（Reality）聯結在一起，所以它直接決定了用戶體驗的好壞。高精度的空間定位可以讓VR用戶體驗到完全的沉浸感，AR中的虛擬物體更加逼真；而如果空間定位的精度比較差的話，則VR用戶會頭暈想吐，AR中的虛擬物體會漂移。
空間定位技術又可以分為兩類，outside-in和inside-out。outside-in方案，需要在環境中架設定位基站，對用戶進行定位跟蹤，比如現在市面上的HTC Vive和Oculus Rift。
而inside-out方案，不需要在環境中架設基站，它只通過自身攜帶的傳感器來進行定位，比如微軟的HoloLens、WindowsMR頭盔系列、HTC Vive Focus、蘋果的ARKit、谷歌的ARCore和Project Tango，還有預計將在2018年發售的Oculus Santa Cruz、Magic Leap One。（當然，還有室內機器人、無人機、無人車）。
outside-in方案要架設定位基站，並且對用戶的活動范圍有限制，相當於"PC互聯網時代"，只適合PCVR。
而insiede-out方案更加輕便自然，就像人用眼睛觀察，就能知道自己在環境中的位置一樣。用戶只需要戴着頭盔，就能隨意走動，沒有運動范圍限制，相當於"移動互聯網時代"。從目前的發展趨勢來看，inside-out方案正在逐漸走向成熟，被應用在越來越多的VR一體機和AR眼鏡上。
而我接下來要探討的是，如何讓建立的地圖點更加准確。因為無論哪種算法，准確跟蹤定位的前提都是，要有准確的地圖點。但是，由於傳感器噪聲和環境噪聲的存在，總是會不可避免地創建出一些誤差大的甚至假的地圖點，這些錯誤的地圖點會對跟蹤定位的准確性造成很大的影響。
所以，我接下來就介紹一種可以讓地圖點更加准確的方法——深度濾波，它可以通過多次觀測，計算出一個地圖點是假的地圖點還是真實地圖點的概率，並且計算出其作為真實地圖點的最有可能的位置。它最初是用在三維重建中的(http://george-vogiatzis.org/)，后來被用在SVO（Semi-Direct Monocular Visual Odometry，https://github.com/uzh-rpg/rpg_svo）中，都取得了很好的效果。
本文主要是參考《Video-based, Real-Time Multi View Stereo》，按照我自己的理解，從最基礎的貝葉斯公式開始推導，一步步推導出最后的迭代公式，用最通俗易懂的方式，把深度濾波方法的來龍去脈講清楚。
1.測量值區間累加的方法
首先，是以前的基於投票的獲取深度的方法。要有一個特征點，然后在其它幀找匹配。
左圖是該特征點在一幅圖像上面的匹配，隨着Z的變化，對應在極線上的不同投影點位置。該特征點圖塊經過仿射變化之后，與不同的投影點位置的圖塊進行NCC匹配，獲取不同的匹配分數，然后記錄極大值點，落在哪個Z值區間。從Zmin到Zmax分了50段區間。右圖是把60幅圖像這樣的數組累加，把極大值點累加，記錄在每個Z值區間出現的極大值的頻數。
最后，會得到這樣的結果，最后會收斂。分別對應5、10、20、100張圖的時候的樣子。
這種方法比較准確，但是需要記錄大量的數據。比如上面（a），單是一幅圖上，從Zmin到Zmax就分了50段，於是要用50個元素的數組來記錄每段有無極大值。然后把60幅圖像的數組都合並起來，得到（b）。
而且，即便能知道累加后最高的那個直方圖對應的深度值，但不能知道，根據所有的測量結果所推算出來的深度值，也不知道，那個特征點對應的外點（假地圖點）的概率。
2.能推算深度值也能算外點概率的方法——窮舉法
所以，要想出其它的辦法，既可以綜合所有的測量結果推算出深度值，也可以推算出這個特征點對應的外點的概率。
那么，就用另外一種方法，不再記錄單幅圖像上的所有極大值，而是只記錄單幅圖像上的最大值。同樣地，也把所有的數組累加起來。積累了150幅圖像的測量值之后，可以得到下圖。
那么，在這樣的測量數據
的情況下，從這些測量數據中，可以反映出這個點的真正的深度值應該是多少呢？是外點的概率又是多少呢？
用
表示該特征點真正的深度值，用
表示內點的概率。
所以，可以用概率方法表示如下，也就是后驗概率。
要求出一組
和
，使得上面的
最大。但是，直接算，算不出來。用貝葉斯公式來轉換。
又因為
和
的區間分布的先驗值，就是每處都相同。
這里的
都是指概率密度，就像是把1千克的質量分散在這么大的體積里，每個地方的密度。要變成概率，還需要乘以一個微分段
，但因為在等式兩邊都乘以同樣的微分段，所以可以約掉了。貝葉斯公式也就可以同樣適用於概率密度的推算了。概率密度與實際概率的關系，參考《從貝葉斯到beta分布》。
所以，
又因為
之間互相獨立，所以，
也就是，找出一組
和
，使得
這個累乘結果最大。
而對其中的每一個元素，
可以認為是，內點的高斯分布概率，加上，外點的概率，表示成如下公式。其中，
是匹配點在極線上的一個像素的擾動，造成的對應點在深度上的變化。參考《視覺SLAM十四講》13.2.3的《高斯分布的深度濾波器》。這里的
都是指概率密度。
在
上，隨着
和
的調節，會變成如下不同的概率曲線，
所以，現在已知有這樣的確定的測量值。
每個測量值都對應一個
，所有的累乘。目標是，通過調節
和
，調節曲線的形狀，使得
這個累乘結果最大。其實，即便
和
相同，
所在的區間相同，
也是不同的，因為
里面的
是由兩幀相機的相對位置以及特征點的像素位置決定的，所以每個測量值的
都不相同。因為在該幀上面，假設在極線上測到的極大值點
，滿足方差為1個像素的高斯分布，即真實的極大值點與測量點
間滿足這樣的高斯分布。則計算出來的深度點
也會滿足高斯分布，方差用
表示，參考《視覺SLAM十四講》的《深度濾波器》。
但由於實際上
差別不大，為了思考上的方便，這里就先忽略掉
的差別。但實際計算上，還是會根據公式計算，考慮公式中
的差別。
將
展開，得到，
通過思考，會發現，測量值都是已知的固定的，就是上面的灰色的直方圖，為了使得累乘結果最大，需要調節曲線，使得直方圖高的地方概率高，就是曲線的形狀接近於直方圖的形狀的時候，所有測量值對應的概率的累乘結果會是最大的。
還是直接用窮舉的方法來找，怎么樣的
和
可以讓累乘結果最大。通過窮舉所有的
和
，得到每個組合對應的累乘結果
的值，畫成三維圖表示如下。分別對應5、10、20、100張圖的時候的樣子。
最后，可以收斂。根據窮舉的結果，可以找到讓累乘結果最大時的
和
，也就是最終想得到的
和
。
其實，上面的方法的本質是基於這樣的一個假設，如果那個特征點是內點的話，最后形成的直方圖會是中間非常凸起的形狀。如果那個特征點是外點的話，則最后形成的直方圖會是非常平坦的長方形的形狀。所以，上面的方法也可以近似理解為，看最后形成的直方圖的形狀，分布與凸起形狀和平坦的長方形形狀的相似度。
但是，對一張直方圖都要進行這樣的窮舉，計算量太大（
分50段，
分100段，要算5000種組合），而且還只是對應單個特征點。而且，每新來一張圖片，新來一個測量值
，調整直方圖后，針對這個新的測量值
，再窮舉出一幅三維圖
，然后之前算出來的三維圖
再乘以新的這幅三維圖
，也就是每種組合的計算結果對應相乘，得到新的累乘三維圖
，再在上面找峰值。公式如下。
問:上面的這種已知所有的測量值（直方圖），然后窮舉所有的
，求概率最大的地方，與卡爾曼濾波的用測量值更新得到真值的方法，有什么不同呢？首先，卡爾曼濾波只針對高斯分布成立，所以只對
有效。那么就把
都當成零，認為全是內點的情況下，兩者的區別呢？
答：本質上是一樣的，結果上也會是一樣的，只是使用場景不同。這種方法的特點是，知道每一個測量值各自對應的方差
，然后窮舉所有的
，使得累乘結果最大，就認為得到最大后驗概率。這種窮舉的方法，適合離線方法。而卡爾曼濾波，一步步迭代的，通過第一個測量值得到，預測值和對應的協方差。認為真值是不動的。然后用后面的測量值和各自對應的協方差，來更新。最后得到一個較准確的值。卡爾曼濾波的方法，適合在線的迭代的方法。在只有一個測量值的情況下，兩者是幾乎相同的。
如果在只有第一個測量值的情況下，上面的窮舉方法，窮舉的結果會發現，當
，
等於測量值的情況下，后驗概率會是最大的。也就是說，在只有一個測量值的情況下，就直接認為它是內點。所以，需要多一些測量值才比較好。
3.窮舉法的改進方法，已知是否內點，推導出
形式，用迭代加速計算
但是，每新來一個測量值都要針對新的測量值窮舉，然后再兩幅三維圖相乘，計算量還是蠻大的，還有沒有更方便計算的方法呢？
首先，對於這種窮舉找最大值的方法，現在一般都是用凸優化方法來代替了。如果能計算出優雅的雅克比解析式的話，就可以用梯度下降法，高斯牛頓法或LM方法來快速找到極大值點（可以轉換成極小值點來算，比如被1減）。但是，從展開的公式上看，也沒有方便地計算雅克比的方法。如果不能方便地計算雅克比，那么就用數值擾動的方法來求雅克比，這就接近於窮舉法了，但只是用數值擾動來算雅克比，沿着梯度變化找到極大值點，計算量應該還是比窮舉法要小的。這就是一般的窮舉法，凸優化法（有雅克比解析式，和，無雅克比解析式用數值擾動算雅克比）的關系。（http://www.cnblogs.com/ilekoaiq）但這些凸優化方法的前提都是，三維曲面上，只有一個極大值點。
通過對
展開，
發現並沒有什么方便計算的方法。主要是，其中的每個大括號里，都是由兩項加在一起的，這樣子，整個公式拆開來算的話，會非常龐大，而且找不出規律（雖然有點像快速傅里葉變換的蝶形乘法）。那么，有沒有辦法，把其中的加法消除掉呢？通過觀察發現，其中的加法，是因為不知道這個點是內點還是外點，於是就兩種情況都算。但是，如果通過其它方法，已經可以知道對應的這個點
是真實點高斯分布產生的，還是外點均勻分布產生的。那樣子，就可以把大括號里面的加法消除掉，只保留對應的一項。要么是真實點高斯分布產生的，要么是外點均勻分布產生的。為了方便理解，可以先想象一種情況，就是它們全都是真實點高斯分布產生的。在論文中，是在極線上搜索，在一定搜索范圍內能否找到極大值，如果能找到極大值，就把極大值測量點轉換成的
認為是真實點高斯分布產生的；如果不能找到極大值，就把最大值測量點轉換成的
認為是外點均勻分布產生的，其實這時候，也無所謂選取的最大值測量點的位置和轉換后的
的位置，因為產生它們的都是外點的均勻分布
。
再回到原來的公式。
設
表示這個測量點轉換成的
是否是真實點高斯分布產生的。
，表示是真實點高斯分布產生的。
，表示是外點的均勻分布產生的。則其中的每一項，可以表示如下，這也是最根本的推導。其中，
，來自於
。單個字母
表示某個取值的時候，
表示所有取值的概率累加，且
。要思考清楚各個變量之間的關系，參考《貝葉斯與維恩圖》。
（其中，
（或者
），也可以先看右邊的這項
（或者
），再看左邊的這項
（
），會更好理解。）
那么，現在，我們已經通過了之前的那種策略，知道了當前點是否內點，記為
，相比只知道
，又增加了信息。則在已知
、
的情況下，
的概率如以下公式所示。因為最根本的貝葉斯公式就是，
。
上式中，有些是內點，有些是外點，根據
來決定，各項的實際表達式會不同。為了數學形式上的簡潔與統一，上式中的各項可以用以下公式來表示，同時表示了內點和外點的情況。
所以，上上式可以簡潔地表示成，
我們要算的公式中，自變量是
，其它的量都是已知的。
是已知的，用
表示累加結果，
。
都是已知的常數，就算
會隨着
的變化而變化，那也沒關系，變量
的變化不經過它們影響因變量結果，所以它們的累乘結果可以認為是一個已知的常數。參考《從貝葉斯到Beta分布》，其實最后要算的是各個（
組合）概率之間的比例，也可以理解為其實是要歸一化的。所以，上式還可以寫成更簡潔的形式。（《Video-based, Real-Time Multi View Stereo》的附加材料的第一部分講的就是這樣的推導，但是它轉成了
，又轉回來，反而有些麻煩）。
是一個
分布。
是高斯分布的累乘，根據高斯分布的性質（以及2中的那個"隱含"高斯分布概率的推導），仍然還是高斯分布。所以，就用
的形式來表示。
然后，就可以很方便地迭代計算了。其中，
比較好算，
可以根據高斯分布乘法的性質，從
開始，一步步迭代計算出來。
當
時，上式中
展開相乘，結合高斯分布乘法，可以得到，
由高斯分布乘法可得，參考《高斯分布的加減乘除法》《直觀理解高斯函數相乘》《卡爾曼濾波基礎》。
如果
，則
，
。
所以，迭代一下后，得到的新的表達式為，
其中，
，
。這樣子，每新來一個測量值，只需要迭代一下就可以，不用再兩幅三維圖相乘。
（而另一方面，如果不知道
的話，是沒辦法明確推導出上面的這樣的Beta和高斯分離出來的公式的，只能用近似的方法。通過觀察最后的收斂的三維圖像，思考發現。當
一定時，累乘結果峰值的位置主要是由
來決定的。也就是
決定了曲線凸出來的地方，當曲線凸出來的地方的位置與直方圖的高值地方接近的時候，累乘結果就能取到較大的值。當
值一定時，
的分布也是一個中間凸的曲線。
的分布和
的分布是分別窮舉的，但不是獨立的，因為在
不同的地方，
的分布會發生變化，
的頂點位置會發生變化。但是，為了簡化計算，可以認為
和
是獨立的，所以分別用曲線去近似
的分布和
的分布。）
4.更完整形式的迭代計算
因為根據我們的策略，我們可以知道這個測量點對應的是內點還是外點。但是，我們之前的策略只能判定明顯的外點，如果不是明顯的外點的話，就認為是內點，但是，真實的情況是，這時候認為的內點，其實還是有可能是外點的。那么，就要考慮這種更真實的情況。
當用之前的策略判定當前測量點為外點產生的時候，就認為它是由外點產生的。直接使用第3部分中的迭代方法，更新這個外點。但是，當判定不是外點的時候，就仍然還是要用第2部分中的方法。
設之前的部分已經滿足
的形式，用公式表示如下。
則，如果新來一個測量點，沒有被判斷為外點。則公式更新如下。
根據高斯分布乘法的性質，處理上式中的
，參考《高斯分布的加減乘除法》。
用
表示
，用
表示
。（注意，這里與《Video-based, Real-Time Multi View Stereo》的附加材料里的公式不一樣，因為我認為《Video-based, Real-Time Multi View Stereo》的附加材料里的
的表達式寫錯了。）則上上式可以表示為。
這時候，結果是，兩個
形式的相加。相當於兩幅三維圖疊加（想象兩個高斯分布在概率軸方向上疊加在一起，形成了兩個山峰）。相加之后的結果，不會再是
的形式，或者其它的單一模型的形式。在這里，只能近似地認為，相加后的結果接近於
的形式，這要依賴於上式中的，（1）
和
很接近使得
接近於1；（2）
和
，
和
很接近，使得
和
很接近；（3）
和
兩個數都很大，使得
接近於
；（4）
很小，使得左邊項可以忽略。
所以，如果近似地認為，相加后的結果接近於
的形式的話，則用公式表示如下。
用一階矩和二階矩近似的方法，把新的參數估計出來。
就是用新的單一模型的概率密度三維圖去近似舊的多模型疊加的概率密度的三維圖。讓兩者的自變量軸（x、y軸）的一階矩和二階矩近似。用
表示舊的多模型疊加的概率密度的三維圖，用
表示新的單一模型的概率密度三維圖。其實這個三維圖的
軸，就是表示
的那一個軸，是無窮大的，只是最后顯示的時候，取了中間的
部分。
參考《高斯分布的各階矩》，
圖上的
的一階矩為
，
圖上的
的二階矩為
，
參考《常用分布函數及特征函數》，
圖上的
的一階矩為
，
圖上的
的二階矩為
，
圖上的
的一階矩為，
圖上的
的二階矩為，
圖上的
的一階矩為，
圖上的
的二階矩為，
所以，聯立方程，讓
圖的矩去等於
圖的矩，也就是
圖上的參數所要滿足的約束。
經過計算，最后得到迭代公式，
相當於是用一個Beta分布來近似
的分布，用一個高斯分布來近似
的分布。
於是，最終就推導出了這么一種簡潔的迭代公式。每新來一個測量值，迭代一下能快速計算出新的峰值點。相比第二部分，不用再針對新的測量值窮舉所有的
組合后再兩幅三維圖相乘；相比第三部分，更加接近真實的情況。
5.總結
深度濾波方法，理論推導非常優美，全部都由最簡潔的貝葉斯公式推導而來，它盡可能地把所有的觀測信息都利用上；還能進行快速迭代，實際運行的效果也很棒。應用在三維重建中，可以讓三維重建出來的模型更加真實;應用在VR\AR空間定位中，可以讓跟蹤定位更加准確。相信它還可以被應用在更多的跟蹤、建圖算法中。
6.求贊賞
您覺得，本文值多少？
7.參考文獻
-
George Vogiatzis, Carlos Hernandez. Video-based, real-time multi-view stereo [J]. Image & Vision Computing, 2011, 29(7):434-441.
-
Forster C, Pizzoli M, Scaramuzza D. SVO: Fast semi-direct monocular visual odometry[C]// IEEE International Conference on Robotics and Automation. IEEE, 2014:15-22.
-
Pizzoli M, Forster C, Scaramuzza D. REMODE: Probabilistic, monocular dense reconstruction in real time[C]// IEEE International Conference on Robotics and Automation. IEEE, 2014:2609-2616.
-
高翔.視覺SLAM十四講[M].北京:電子工業出版社,2017:325-327.