本文詳細介紹一篇雙目事件相機里程計的論文，Event-based Stereo Visual Odometry，港科大沈邵劼團隊Yi Zhou和TU Berlin的Guillermo Gallego共同完成，並公佈了程式碼。我準備在接下來一段兒時間認真研究一下這篇論文與程式碼，歡迎想深入研究這個方法的朋友與我交流。轉載請註明出處

論文下載：arXiv, https://arxiv.org/abs/2007.15548
github程式碼：https://github.com/HKUST-Aerial-Robotics/ESVO
柯翰同學在公眾號的一篇簡單介紹：基於雙目事件相機的視覺里程計

1. 演算法流程

演算法總體流程如上圖所示，主要包括三個重要模組（虛線矩形框），分別是：資料處理、Tracking和Mapping。暫時不去管初始化部分，我們關注資料流的方向，可以發現，在Tracking時，傳入的是Time-Surface（本文往後縮寫為TS）和Local Map，得到運動變換位姿pose；之後，pose和TS等共同完成完整的Mapping過程。

2. 基礎知識

2.1 Time-Surface 時間表面

是一種事件相機資料表達方式，請看之前我部落格中對這個進行的總結：【事件相機整理】角點檢測與跟蹤總結 - Time surface / SAE

2.2 Student t 分佈

wiki介紹：https://en.wikipedia.org/wiki/Student%27s_t-distribution
作者通過對資料集測試，發現左右目相事件機中，匹配畫素的殘差$r$統計上符合“學生t分佈”（中文表述不清，原文為：$r$ denotes the temporal difference between two corresponding pixels $x_{i,1},x_{i,2}$ inside neighborhoods）由此實現了多幀影像對深度進行融合。

2.3 ZNCC 零均值歸一化相似性

是雙目畫素匹配時利用的一種度量方法，在極線約束上尋找匹配的畫素點。在演算法初始化部分使用。

2.4 IRLS 迭代重權重最小二乘

wiki介紹: https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares
一種非線性迭代優化演算法，在Mapping過程中利用IRLS得到最優逆深度

2.5 Compositional LK

LK光流跟蹤方法的一種變形，增量式而不是絕對式，避免了重複計算完整的雅克比矩陣。在Tracking過程中不斷更新計算Warpping Parameter，從而得到運動軌跡。注意每次只是計算更新的部分$W(x,\rho ;\theta )\leftarrow W(x,\rho ;\theta )\cdot W(x,\rho ;\Delta \theta )$，即式中的$\Delta \theta$

3. Tracking部分

3.1 TS negative

論文並沒有使用原始的TS進行Tracking，而是使用了 TS negative（後面簡寫為TSn)，$\overline{\tau }(x)=1-\tau (x)$。作者指出，由於TS的一側往往是斜坡狀，另一側是陡峭的，所以沿著斜坡可以搜尋到當前Edge在影像中的位置。由此，當使用$(1-\tau )$時，值越小意味著距離真實Edge越近。

3.2 目標函式

上面提到，當TSn越小，這意味著影像這個位置對應了真實世界的Edge。在Tracking時，我們需要根據已有的Mapping結果（Local Map），將已知逆深度的Edge點投影到Frame上，計算對應的TSn的數值，進行相加，越小這意味著投影的引數約正確。例如，當位姿完全正確時，真實世界的Edge投影到Frame時，對應的TSn應該全接近0。從而我們可以寫出目標函式為：$${\theta }^{\ast }=\underset{\theta }{\mathrm{argmin}}\sum _{x\in S}(\overline{\tau }(W(x,\rho ;\theta ){)}^2,$$ 即尋找一個最佳的Warpping Parameter $\theta$，使Frame中所有有逆深度的點($x\in S$)經過Warp後對應的TSn儘可能小。其中Warp就是定義為反投影后經過6DoF變換再投影，$$W(x,\rho ;\theta )={\pi }_{left}(T({\pi }_{ref}^{-1}(x,\rho ),G(\theta ))),$$ 即根據參考位置反投影，經過運動引數$G(\theta )$變換再投影到左目。論文指出，只利用左目進行計算就夠了，加上右目沒有明顯的提升。

3.3 優化過程

論文指出，利用Compositional LK演算法，不斷迭代更新，得到最佳的T。

4. Mapping

Mapping部分主要也有兩部分，首先對event的逆深度進行估計，之後再進行半稠密重建。

4.1 Event的深度估計

4.1.1 目標函式

深度估計時，需要知道Tracking的軌跡，進行運動補償（所以需要一個初始化部分，才能夠保證Tracking部分能夠首先進行），估計一個event的逆深度$\rho$。

當我們估計在$t$時刻一個event的深度$\rho$時，我們可以得到這個event在左右目的座標$x_1^t,x_2^t$，如果深度估計的正確，這兩者對應的TS應該是相同，且之前一段兒時間$\delta t$(100Hz)內對應的TS應該也是相同的。所以目標函式定義為： $${\rho }^{\ast }=\underset{\rho }{\mathrm{argmin}}C(x,\rho ,{\mathrm{T}}_{left}(,t),{\mathrm{T}}_{right}(,t),T_{t-\delta t,t})$$ 其中${\mathrm{T}}_{left/right}$為該畫素點在這段兒時間內的TS對應的值，而$C$具體為：$$C(...)=\sum _{x_{1,i}\in W_1,x_{2,i}\in W_2}|{\tau }_{left}^t(x_{1,i})-{\tau }_{right}^t(x_{2,i}){|}_2^2$$ 即對應的TS值的差。而$x_1,x_2$為：$$x_1=\pi (T_{c_{t-ϵ}}^{c_t}\cdot {\pi }^{-1}(x,{\rho }_k),$$ 即經過當前假設逆深度投影后，經過對應時間戳投影到指定時刻的Frame。

4.1.2 優化

優化採用任何一種非線性優化方法即可，例如GN方法

4.1.3 初始化

初始化時，由於不知道 $T$ 的軌跡，所以直接在極線上搜尋匹配，只用ZNCC即可。不知道$T$所以無法得到 $T_{c_{t-ϵ}}$，所以不是很精確，但作為初始化已經足夠了。

4.2 半稠密重建

由4.1，我們得到了一個稀疏的深度圖，下一步通過融合的方式，得到半稠密的深度圖。在這個過程中，並不是利用一次觀測，而是利用了多次觀測的結果，所以同時收斂了深度的取值。這部分內容類似深度濾波，但提出了新的融合方法。

如圖，每次觀測得到一個稀疏深度圖，經過一段兒時間後（M次觀測）融合得到更為精確、半稠密的深度圖。

作者通過真值，統計兩個資料集估計逆深度的殘差分佈，發現是個Student t分佈（下文縮寫為St）。當某個畫素有多次觀測時，多個逆深度進行融合，得到新的St分佈。

接下來用到了IRLS進行優化，但我還沒能理解這部分是在幹啥，就先跳過了。

在融合時，將每次觀測計算得到的有深度的畫素點投射到下一次觀測的像平面，由於不是整數畫素，所以在新的像平面中考慮周圍4畫素是否是同一個資料點。1. 如果是同一個資料點，則利用St分佈疊加方式將上一次的逆深度分佈作為先驗，疊加最新的觀測逆深度分佈獲得後驗逆深度分佈；2. 如果不是同一個資料點（由二者逆深度分佈的均值和方差判斷），則取方差較小的作為當前畫素點的逆深度；3. 如果最新的測量中沒有得到這個畫素點的逆深度，則直接定為上次的分佈。由此，實現了多次測量的逆深度的融合，也得到了更稠密的深度圖。

5. 小結

本文我認為最有創新之處就是多次測量的融合方法。

其他很多地方作者考慮的很細緻，也有很多技巧。參考文獻有五六十篇非常之多，只能說港科和Guillermo Gallego在EB和Stereo VO上有很深的積累。