本文详细记录了EMVS这篇论文的一些细节。主要用于个人记录，转载请注明出处。
本人近期准备深入研究这个算法，并尝试做相应改进，所以总结了此文，梳理自己的思路。本文并没有把论文分析的很完整，不适合完全没有了解过EMVS的人阅读。如果你准备深入研究这个方法，有自己的想法，欢迎和我交流：jfjxxdy2007@126.com
同时指出，本文主要围绕ETH于2018年在IJCV发表的"EMVS: Event-Based Multi-View Stereo—3D Reconstruction with an Event Camera in Real-Time"[1] 一文进行展开，而不是在此之前发表在BMVC发表的会议论文"EMVS: Event-based Multi-View Stereo"[2] 。虽然两者本质上是一样的，但[1]更加详细。
算法的思路非常清晰明了，对应的代码也很简洁。但论文的安排个人感觉有些杂乱，我在此进行一下梳理。 首先我们要了解EMVS的基本思想：我们需要提前知道相机真实的运动轨迹，在运动过程中，从每个触发事件的像素进行反向投影投射出一条射线到所在的空间中，并对这个空间进行了网格划分（这个空间在论文中用DSI进行表示，虽然DSI的本意并不是这个意思，请看这篇博客：DSI到底是指什么？）。当将所有事件反向投影为射线以后，在DSI空间中由同一个事件源产生的事件应该是多条射线的相交，我们在DSI网格中选取局部极大值，从而确定了事件源的深度。
可以看出，这种方法基于的一个重要假设是：相机运动轨迹完全已知，即为实现了相机自定位情况下的建图。
完整的算法非常清晰明了，只有3步，对应的代码也是主要有3个函数：1. 对每个事件反向投影寻找射线在DSI中穿过的网格；2. 选取网格的局部最大值；3. 点云的后期处理。
从论文安排上看，1-4章介绍了一些意义不大的背景知识，5章介绍了完整的方法，包括：在DSI中寻找反向投影射线穿过网格的Space-Sweep方法、Confident Map的建立以及局部最大值的寻找、地图融合与参考点的选择、点云后期滤波等，第6章分析了为什么采用projective的方式进行网格划分而不是传统的定尺寸（我觉得这部分放在这里不太合适，虽然作者也提到可以跳过，我感觉是为了期刊发表而在[2]上新增的内容），第7章给出了第5章计算反投影射线的具体方法，往后是实验总结。
EMVS的核心工作我觉得就可以算是这个Space-Sweep算法。原始的算法请参考论文[3]。基本思想是下图：首先在轨迹上一点选择一个位置放置一个虚拟的相机，假设待测量的点在这个虚拟相机的一定深度范围内 [ Z m i n , Z m a x ] [Z_{min}, Z_{max}] [Zmin,Zmax]，把这个区间均匀切片，每片距离相机的距离记为 Z i Z_i Zi 。我们要做的是，找到某一个位置上相机平面的一个事件在 Z i Z_i Zi 平面上的坐标。显然平面到平面的映射是一个单应性矩阵，用 H Z i H_{Z_i} HZi 表示“从真实相机相平面，到虚拟相机 Z i Z_i Zi平面”的映射。笨方法就是，依次计算相机相平面到所有的 Z i Z_i Zi 的单应性，但这样效率太低，而是采用计算出到 Z 0 Z_0 Z0（第一个切片）的单应性矩阵后，通过 H Z i H Z 0 − 1 H_{Z_i}H_{Z_0}^{-1} HZiHZ0−1 计算到剩下所有的位置。为什么这样计算？后面我们进行推导可以发现， H Z i H Z 0 − 1 H_{Z_i}H_{Z_0}^{-1} HZiHZ0−1 是一个很方便计算的东西。
现在我们进行一下推导，这个计算过程到底是怎么计算的。在此之前，解释一下为什么是叫“改进的Space-Sweep算法”。原始的算法，在每一层切片中按照恒定大小划分 X Y XY XY网格，可以理解成从 Z 0 Z_0 Z0 到任意 Z i Z_i Zi 上的坐标都是统一的物理量（例如真实世界的物理单位），是uniformed的，而这里改进为了projective的，即每一个网格切片都有恒定数量的 X Y XY XY网格，而每个网格的对应的真实物理大小是不同的。为什么这样选取将在下面进行介绍，这也是论文第6章的主要内容。
现在带入第 Z i Z_i Zi 个切片的距离，以及法向量 e 3 e_3 e3 ，得到论文式(8)： H Z i − 1 ∼ R 1 Z i t e 3 T H_{Z_i}^{-1}\sim R\frac{1}{Z_i}te_3^T HZi−1∼RZi1te3T （本文最后数学部分解释了这个单应性矩阵的由来），之所以是“逆”是因为不带逆是从真实相机相平面到切片，而这个公式计算的是统一的虚拟相机坐标系下的 R , t R,t R,t 对应的单应性变换。当 i i i 取0和任意值后，便得到了从 Z 0 Z_0 Z0 切片到 Z i Z_i Zi 切片的更新： H Z i H Z 0 − 1 = ( R + 1 Z i t e 3 T ) − 1 ( R + 1 Z 0 t e 3 T ) H_{Z_i}H_{Z_0}^{-1}=(R+\frac{1}{Z_i}te_3^T)^{-1}(R+\frac{1}{Z_0}te_3^T) HZiHZ0−1=(R+Zi1te3T)−1(R+Z01te3T)，这就是论文的公式(13)，利用矩阵求逆引理（补充在了博客最后）化简第一部分，得到论文公式(14)： H Z i H Z 0 − 1 ∼ I + Z 0 − Z i Z 0 ( Z i − C z ) C e 3 T H_{Z_i}H_{Z_0}^{-1}\sim I+\frac{Z_0-Z_i}{Z_0(Z_i-C_z)}Ce_3^T HZiHZ0−1∼I+Z0(Zi−Cz)Z0−ZiCe3T，注意这里用到了 C = − R T t C=-R^Tt C=−RTt 即光心与变换参数的关系。我为了方便书写，不带角标的字母则表示完整的向量（例如 C = [ C x , C y , C z ] T C=[C_x, C_y, C_z]^T C=[Cx,Cy,Cz]T ）。那么得到 H Z i H Z 0 − 1 H_{Z_i}H_{Z_0}^{-1} HZiHZ0−1 后，右乘投影射线在第0个切片 Z 0 Z_0 Z0 的坐标 [ x ( Z 0 ) , y ( Z 0 ) , Z 0 ] T [x(Z_0), y(Z_0), Z_0]^T [x(Z0),y(Z0),Z0]T 后，进行归一化，得到在第 i i i个切片的坐标，即公式(15)： x ( Z i ) = Z 0 Z i δ x ( Z 0 ) + 1 Z i ( 1 − δ ) C x x(Z_i)=\frac{Z_0}{Z_i}\delta x(Z_0)+\frac{1}{Z_i}(1-\delta)C_x x(Zi)=ZiZ0δx(Z0)+Zi1(1−δ)Cx和对应的 y ( Z i ) y(Z_i) y(Zi)，其中 δ = ( Z i − Z 0 ) / ( Z 0 − C z ) \delta=(Z_i-Z_0)/(Z_0-C_z) δ=(Zi−Z0)/(Z0−Cz)。这里进行耐心的计算和合并同类项即可，注意到这个与原始的Space-Sweep算法不同，多了两个系数因子 Z 0 / Z i , 1 / Z i Z_0/Z_i, 1/Z_i Z0/Zi,1/Zi，就是因为切片的坐标是与 Z i Z_i Zi 成比例的，如果都是统一的物理坐标，完全可以通过相似三角形变换得到（直接用相似就可以得到参考文献[3]中的公式(6)，而不需要原文的公式(5)）。
最后完整的捋一下思路：对于每个实际的相机位姿，通过单应性计算相机平面到第0个切片的单应性，然后通过公式(15)快速计算剩下所有切片中的坐标。
在完成计算反投影射线穿过的DSI的栅格之后，我们需要从这个栅格中寻找局部最大值。论文采用的方法是，首先建立一个Confident Map，然后寻找局部最大值。
所谓Confident Map，就是计算在虚拟相机每个像素的反投影射线穿过的每个栅格投票的最大值，得到Confident Map。这个名称的含义，可以理解成，这个视角下最可能是事件源的点。
之后在Confident Map这个二维平面中，通过自适应阈值的方法，选取区域最大值。论文采用5*5的区域，常数 C = 10 C=10 C=10。保留最大值，获得在虚拟相机像平面的坐标，再索引得到对应的深度。
论文中在5.2.4节提到，每当运动超过一定距离（例如当前参考的平均点深度的15%~40%）后，插入新的参考帧，即新的虚拟相机；而在代码中，选择了一段儿轨迹中间位置作为虚拟相机；
5.2.5小节提到，通过每个点附近一片区域内是否有其他的点，判断是否为离群点，从而删除；
7.2小节提到，每个事件计算一次真实相机的位置再反投影的计算代价是极大的，采用累积一批事件后计算一个位置，进行批量反投影计算。代码中采用1024个事件为一组，相机位置采用时间戳进行线性插值得到；
最简单的方式是，对最近的一个坐标网格进行投票；而论文和代码中，采用这个点附近4个整数点，进行双线性插值的方式对4个顶点进行投票。
论文的第6章详细介绍了projective的网格，而不是uniformed，两者之间的区别，以及projective的优势。说实话我搞不懂6.1这么一大段，尤其是正中间极其显眼的示意图，除了便于理解外有什么意义。既然作者提到了，我也介绍一下。
首先左上角为虚拟相机采用projective的方式进行划分切片的网格，保证每个切片的网格数量是一样的，但实际的物理坐标不同。(b)图表示uniformed方式，在DSI空间反投影的结果，注意到坐标单位是m，而且投射的结果是直线，很好理解；右下图 ( c)表示如果用projective的方式，投射的反投影在Y轴的切片上，投票路径是曲线（左上我绘制了绿色的网格，对应a图的绿色射线，绿色射线依次穿过1-10，而在projective的网格下实际的路径是条曲线），右下角(d)图表示，如果用逆深度绘制这个DSI空间，那么反投影射线基本是直线。然后用了6.1整个小节进行了证明，为什么( c)是曲线，(d)是直线，式(4)对 X , Y X,Y X,Y 分别除以了 Z Z Z 即意味着实际距离与深度有关，而式(5)(6)还用到了消逝点、直线空间表示等几何知识，真的是醉了。
更有趣的是，6.2小节理论证明了如果用uniformed方式，像素平面的一个点反投影后在远处的切片会对应更大的面积（这么一个很显然的事情），用到了附录的大量证明，以及矩阵行列式引理等，就是为了说明一个结论：如果采用uniformed方式，一条射线在更远的地方对应了更多的网格，而用projective的方式，对应的数量基本就是1个网格。
那么就要分析一下，为什么projective的方式比uniformed要好？作者指出有以下原因：1. 物理意义上，一个事件源投射到相机平面，再投射回去，用uniformed会对应更多的网格，这个意义有点儿奇怪；2. 自适应阈值选取最大值时，不需要考虑尺度问题；3. 节省内存；4. 精度足够VO使用。而我觉得，最重要的一点，是作者在进行 Confident Map 时，采用了投影的方式将DSI投影到平面，如果是uniformed的DSI，那么极其不方便投影，将会是一个下采样，也不方便索引。
作者开源了代码：https://github.com/uzh-rpg/rpg_emvs 有开源代码才是我认真研究这个方法的主要原因。
这个代码主要由3个函数，即对应：反投影与DSI计数、寻找局部最大值、点云生成与处理。其中包含了一些技巧，例如类的派生继承与动态绑定，以及Eigen库对4x4矩阵的加速的利用，L1缓存的利用等。如果是我写，至少效率低1个数量级吧。
单应性矩阵的基础知识：假设有两个坐标系，分别记为1系和2系，1到2系的变化参数是 R , t R,t R,t即三维坐标点 X X X 有如下计算式： X 2 = R X 1 + t X_2=RX_1+t X2=RX1+t假设1系中有一个平面，其法向量是 n n n，距离坐标原点的距离是 d d d，那么这个平面上任意一点 X 1 X_1 X1 的方程为： n ⋅ X 1 = d n \cdot X_1=d n⋅X1=d，即 1 d n T X 1 = 1 \frac{1}{d}n^TX_1=1 d1nTX1=1。这时可以将 上式写成 X 2 = R X 1 + t = R X 1 + t ( 1 d n T X 1 ) = ( R + t n T d ) X 1 = H X 1 X_2=RX_1+t=RX_1+t(\frac{1}{d}n^TX_1)=(R+\frac{tn^T}{d})X_1=HX_1 X2=RX1+t=RX1+t(d1nTX1)=(R+dtnT)X1=HX1从而得到单应性矩阵 H = R + t n T d X 1 H=R+\frac{tn^T}{d}X_1 H=R+dtnTX1。
( A + x y H ) − 1 = A − 1 − A − 1 x y H A − 1 1 + y H A − 1 x (A+xy^H)^{-1}=A^{-1}-\frac{A^{-1}xy^HA^{-1}}{1+y^HA^{-1}x} (A+xyH)−1=A−1−1+yHA−1xA−1xyHA−1 式(23)带入 x = 1 d C , y = n , A = I x=\frac{1}{d}C, y=n, A=I x=d1C,y=n,A=I，得到式(24)
https://blog.csdn.net/tfb760/article/details/108361375 从而由式(23)得到式(25)，来支撑正文第6章中的公式(9)。
[1].Rebecq, H., Gallego, G., Mueggler, E. et al. EMVS: Event-Based Multi-View Stereo—3D Reconstruction with an Event Camera in Real-Time. Int J Comput Vis 126, 1394–1414 (2018). https://doi.org/10.1007/s11263-017-1050-6
[2]. Rebecq, Henri; Gallego, Guillermo; Scaramuzza, Davide (2016). EMVS: Event-based Multi-View Stereo. In: British Machine Vision Conference (BMVC), York, UK, 19 September 2016 - 22 September 2016, 1-111
[3]. R. T. Collins, “A space-sweep approach to true multi-image matching,” Proceedings CVPR IEEE Computer Society Conference on Computer Vision and Pattern Recognition, San Francisco, CA, USA, 1996, pp. 358-363, doi: 10.1109/CVPR.1996.517097.
