从点云到网格（二）VRIP介绍

VRIP（Volumetric Range Image Processing），顾名思义，是从深度图重建网格的一种方法。VRIP是Brian Curless和Marc Levoy在1996年提出来的方法，距今已经有20年的历史了，依然属于最好的方法之一。

VRIP的核心问题是

已知世界坐标系下，某物体表面ff在不同视角下的深度图f1^,...,fK^f1^,...,fK^，求ff。这里隐含深度图在世界坐标系下的位姿是已知的。

许多三维测量方法，比如激光、TOF、结构光、双目视觉等，都可以得到深度图。因此这是一个非常有意义的问题。如何将这些深度图融合成为一个平滑的单一网格呢？这就是VRIP要解决的问题。
下面给出一个结构光成像的例子。左上图是将不同视角拍摄到的深度图（已转化为三维网格）匹配在一起后的情形。右上图是左图中某一部分的切面，可以看到很多层网格重叠在一起，噪声、匹配误差、采样率等都反映在这个局部中。左下图是VRIP的结果，右下图是VRIP融合的结果和深度图放在一起（棕色是深度图，蓝色是融合后的网格，只能看到一点点）。

VRIP的基本假设

VRIP最重要的假设是测量误差沿着光线传播方向（投影方向）并服从高斯分布。假设fk^(u,v)fk^(u,v)是在第kk个视角下，光线从传感器(u,v)(u,v)位置沿着传播方向到达ff的距离的测量值，fk(u,v)fk(u,v)为其真实值。那么条件概率满足

P(fk(u,v)|fk^(u,v))=ck(u,v)exp[−12(fk^(u,v)−fk(u,v)σk(u,v))2]P(fk(u,v)|fk^(u,v))=ck(u,v)exp⁡[−12(fk^(u,v)−fk(u,v)σk(u,v))2]

因此，VRIP算法和视角非常相关。

VRIP的模型

VRIP尝试从概率的角度来描述核心问题。对任意曲面ff，P(f|f1^,...,fK^)P(f|f1^,...,fK^)为ff的条件概率。那么求解核心问题转化为一个最大似然问题

maxfP(f|f1^,...,fK^)maxfP(f|f1^,...,fK^)

经过一系列独立性假设，

P(f|f1^,...,fK^)=∏Kk=1∏Mi=1∏Nj=1P(fk(i,j)|fk^(i,j))P(f|f1^,...,fK^)=∏k=1K∏i=1M∏j=1NP(fk(i,j)|fk^(i,j))

取对数，并转化为求和

E(f)=−ΣΣΣlog[P(fk(i,j)|fk^(i,j))]E(f)=−ΣΣΣlog⁡[P(fk(i,j)|fk^(i,j))]

将离散情形转换为连续情形

E(f)=−Σ∫∫Aklog[P(fk(u,v)|fk^(u,v))]dudvE(f)=−Σ∫∫Aklog⁡[P(fk(u,v)|fk^(u,v))]dudv

定义d(u,v,fk)=fk^(u,v)−fk(u,v)d(u,v,fk)=fk^(u,v)−fk(u,v)，代入PP的表达式，则有

E(f)=12Σ∫∫Akwk(u,v)dk(u,v,fk)2dudvE(f)=12Σ∫∫Akwk(u,v)dk(u,v,fk)2dudv

上式还是建立在各个视角下的局部坐标系（传感器坐标系）。将上式转换到世界坐标系下，有

E(f)=12Σ∫∫Awk(x,y,z)dk(x,y,z)2[vk(x,y,z)⋅(∂z∂x),∂z∂y,−1)]dxdyE(f)=12Σ∫∫Awk(x,y,z)dk(x,y,z)2[vk(x,y,z)⋅(∂z∂x),∂z∂y,−1)]dxdy

其中z=f(x,y)z=f(x,y) ，vkvk是第kk个视角的投影方向，dk(x,y,z)=f^(x,y,z,vk)−f(x,y)dk(x,y,z)=f^(x,y,z,vk)−f(x,y) ，为三维点(x,y,z)(x,y,z)沿着第kk个视角的投影方向与深度图对应点的测量值f^f^的距离，即Signed Distance Function。积分中点乘那项是从uvuv坐标系转换到xyzxyz坐标系的JacobianJacobian矩阵。
下图是SDF的示意图。注意靠近传感器的方向距离为正，远离传感器的方向距离为负。

下图是同一个三维点在两个视角下的SDF示意图。注意SDF是沿着视角方向（光线传播反方向）的。图中d1<0,d2>0d1<0,d2>0。

于是，问题转换为

求，使得z=f(x,y)，使得，使得z=f(x,y)，使得E(f)E(f)最小。

求解这个问题涉及到很多偏导数、方向导数的知识，本人也并没有完全搞明白。但问题的解却是出乎意料的简单

若z=f(x,y)z=f(x,y)满足Σwk(x,y,z)dk(x,y,z)=0Σwk(x,y,z)dk(x,y,z)=0 ，则zz是最优解。

VRIP算法流程

定义D(x)=Σwk(x)dk(x)Σwk(x)D(x)=Σwk(x)dk(x)Σwk(x)，那么D(x)=0D(x)=0就是我们要重构的三维表面。实际应用中，D(x)D(x)可看作一个三维体数据（volume），D(x)=0D(x)=0通过提取DD的零等值面即可得到。因此，算法先增量构建DD，然后通过等值面提取方法得到三维网格。例如，Marching Cube就是一种高效提取等值面的算法，而且非常适合在GPU上实现。

需要指出的是，作者在论文中采用了TSDF（Truncated signed distance function），即在一条光线上，只考虑测量值z^z^附近一定范围z^±δzz^±δz内的体素。一方面是因为一条光线可能会穿过物体不止一次。另一方面这样也可以减小搜索范围，加速算法。

VRIP算法的框架如下：

/* 初始化 */ 
对每个三维体素，设其权重为0。
/* 深度图融合 */ 
对每个深度图 {
  /* 准备 */ 
  网格化深度图；
  计算每个点的权重；
  /* 更新体素 */
  对该视角下FOV中的体素 {
    沿投影方向找到深度网格中的对应点；
    计算其沿投影方向的SDF；
    插值得到其权重；
    更新这个体素的权重和SDF。
  }
}
/* 表面提取 */ 
提取零等值面。

着重强调一点，深度图的权重需要尽量准确，特别是（1）对噪声大的点能够赋予较小的权重，（2）对于法向量和视角方向角度比较大的三维点，可以降低其权重（曲率大，采样率不够，重建误差大）。

下图是两个视角下的TSDF的融合的示意图。

下图是多个视角下TSDF的融合（两颗真实的牙齿）

VRIP的优缺点

VRIP的优点主要有

• 它是一定意义下的最大似然解。这保证了解的精度。
• 它是一个增量方法，每一次得到新的深度图后，可简单快速地加入到TSDF中。
• 适合并行处理，可用GPU加速

VRIP的主要缺点有

• VRIP生成的网格会附加一定的平滑效果，在存在噪声和匹配误差的情况下有时不能重建出细微的结构。
• 如果深度图存在匹配误差，VRIP并不能消除这些误差。这些误差会反映在融合后的网格中（分层、噪声等）。

VRIP的加速

VRIP算法的时间复杂度比较高。作者在CPU端做了很多优化和加速工作。假设深度图平行于xyxy平面，且投影方向为正交投影，那么沿zz方向的所有体素在深度图上的投影相同，因此其TSDF可以在zz方向上简单计算得到，而且权重相等，不必要重复运算。作者的主要思路是，将不同视角下的深度图和坐标系，通过仿射变换和重采样，映射到相对标准的位置（深度图平行于xyxy平面），从而减少投影和权重的计算量。
如下图所示，（a）某个视角下的深度图及正交投影方向，投影方向和Voxel slices有一个夹角。（b）通过一个仿射变换，将投影方向变换为与Voxel slices相垂直。一般而言，仿射变换后深度图和Voxel slices仍然有夹角，因此还要将它在平行于slices的平面上重采样。（c）将计算得到的TSDF变换回原坐标系下的距离。

VRIP在SLAM中的应用

VRIP算法在RGBD SLAM中有着广泛的应用。Kinect Fusion的作者并没有采用CPU加速的方法，而是将VRIP算法移植到GPU上，利用Ray casing算法做到了网格的实时显示，效果非常好。后来的Kintinuous、Elastic Fusion、Dynamic Fusion也都采用类似的架构去生成网格。

转载请注明作者和出处(http://www.cnblogs.com/luyb)，未经允许请勿用于商业用途。 COPYRIGHT@CNBLOGS.COM/LUYB 联系方式：luyanbin7 at gmail.com