基于图像的多视角立体视觉三维重建(四)——基于MVS算法的稠密点云重建_mvs三维重建

基础知识

  稠密点云重建的目的是在相机位姿已知的前提下，逐像素的计算图像中每一个像素点对应的三维点，得到场景物体表面密集的三维点云。
  稠密点云重建的前提相机的姿态已知,且待重建的三维点具有图像一致性,即三维点位于物体的表面上而非物体内部，那么这个三维点根据相机的内外参投影到不同视角的图像上，以投影点为中心取出两个方形的小patch，那么它们所包含的场景应该是比较类似的.

极线搜索

  两个视角的图像中，A视角中的一点在不同的深度的三维点投影到B视角中形成一条线，称为极线。当三维点的坐标未知，A图像中的一点，在B图像中可能匹配上对应的极线上的某点，在极线上找到B点的具体位置与相机的姿态之后，可以利用三角量量测的得到三位点的深度信息。

  如何判断极线上的某点是否为A点的同名点，主要利用以下约束：光度一致性约束、可视性约束。

---

光度一致性约束

  同一个空间的点在不同的视角的投影应当具有相同的光度，重建的核心在于恢复空间中具有光度一致性的点。（要求物体的表面为朗伯面，玻璃、镜子、水面不属于朗伯面）

光度一致性假设的度量方法
Sum of squared Differences(SSD):
$${\rho }_{SSD}(f,g)=||f-g|{|}^2$$
Sum of Absolute Differences(SAD):
$${\rho }_{SAD}(f,g)=||f-g|{|}_1$$
Normalized Cross Correlation(NCC):
$${\rho }_{NCC}(f,g)=\frac{(f-\overline{f})(g-\overline{g})}{{\delta }_f{\delta }_g}$$
其中f代表以A视图中某点的领域的像素组成的向量，g代表B视图中对应同名点邻域像素组成的向量

---

可视性约束

  1、图像上的点不能被遮挡。 2、重建的点前面不能出现点。 3、不能出现在物体的内部。

---

多视角立体技术

---

基于体素的方法

规则划分

  基于体素的方法等价于3D空间中的像素标定问题，在物体内部的点标记为1，在物体外部的点，标记为0，介于0-1之间的点就是物体的表面。
  标记完之后，使用光度一致性约束和可视性约束来对物体表面的点进行计算，评估标定质量。基于体素的方法的缺点是计算量大，在分辨率很大的时候，计算量会增长很多。

不规则划分

在点云较为稀疏的区域采用分辨率较小的体素，在点云较为稠密与细节比较精确的的区域采用分辨率较高的体素，可以保证物体重建的精度，另一方面可以自适应调整四面体的大小，减少计算量。

基于体素的方法的优化

典型的MRF离散化优化问题

参考图像上的每个像素都分配一个标签（内部或者外部）（无向图的优化）
$$E(f)=E_{data}(f)+E_{smoothness}(f)+E_{visibility}(f)$$
其中$E(f)$代表能量模型。
f表示无向图的所有顶点的向量，需要对每个顶点赋值
$E_{data}$是光度一致性假设，在不同视图上观察到该点的纹理或者光度尽量相似。
$E_{smoothness}$平滑项约束，相邻的体素（四面体）的标签趋于一致。
$E_{visibility}$不能被别的东西遮挡，不能遮挡别的体素。

基于体素的方法的优小结

优点：1、生成规则的点云
   2、便于提取物体的平面

缺点：1、精度受到空间划分分辨率的影响
   2、难以处理精度高、规模大的场景

---

基于空间patch扩散的方法(pmvs)

1、假设空间中的3d矩形patch
2、通过一定规则的扩张方法，使得pathch覆盖物体表面

过程：

1. 初始种子点生成采用sift、HOG等稀疏特征点，对这些点进行深度和法向量的重建
2. 扩张过程对已重建三维点的邻域进行匹配
3. 滤波过程采用光度一致性约束和可视性约束来去除噪点

3D patch定义：
位置、法向量（每个patch具有25个点， pathch的中心为点的位置，法向量为三维点邻域的法向量）

对三维空间上的一个pathch投影到不同视角中，计算NCC
1.在图像上均匀计算HOG/Harris特征
2.沿极线进行搜索找到匹配特征点
3.对匹配对进行三角化建立patch,
  中心: 三角化确定
  法向量: 指向参考图像
  可视图像$V(p)$: 法向量和视角的夹角足够小,通过光度一致性约束对可视图像进行筛选 V ∗ ( p ) = { I ∣ I ∈ V ( P ) , h ( p , I , R ( P ) ) ≤ a } V^{*}(p)=\{I|I\in V(P), h(p,I,R(P))\le a\} V∗(p)={I∣I∈V(P),h(p,I,R(P))≤a}
4.对patch位置和法向量进行优化$$\underset{c(p),n(p)}{\max }{g}^{\ast }(p)=\frac{1}{|V^{\ast }(p)\backslash R(p)|}\sum _{I\in V^{\ast }(p)\backslash R(p)}h(p,I,R(p))$$

---

深度图融合(重点)

基本步骤：

1. 为图像选择邻域图像构成立体图像对(通过图像pair进行深度图重建)
2. 计算每一幅参考图像的深度图，得到每幅参考图像的点云
3. 对所有深度图进行融合

全局视角的选择(建立候选视角池)

初步筛选出一些适合用于重建的图像,局部视角从全局视角中选择.
1、图像具有相同的内容、外观和尺度
2、图像具有较大的视差（宽基线）

其中，$F_V$表示V视角， $F_R$表示R视角，$w_N(f)$代表相机的三角测量夹角大小，$w_s(f)$代表两幅图像的尺度相似性。

图像分辨率的估计

用于衡量图像的分辨率，图像的一个像素的宽度对应的三维空间中的物体尺寸。

相关视角的选择(针对每个像素)

  全局视角的选择时固定的，局部视角的选择是针对每个像素的，每个像素在重建过程中不需要用到全部的全局视角，有些视角存在遮挡。局部视角选择是从全局视角中再一次进行筛选得到。一般从全局视角中选择四个当作局部视角进行优化。

局部视角的选择：
1、NCC值确定候选视角
2、实现要足够分散（不共面）
$$l_R(V)=g_R\ast \prod _{V^{{}^{\prime }}\in A}{w}_e(V,V^{{}^{\prime }})$$

深度图融合算法整体思路

区域生长法扩张（宽度优先搜索）：

1、根据重建的置信度（NCC）建立优先级队列（排序）
2、从初始的稀疏特征点开始深度估计
3、对每个种子点进行非线性深度优化（核心）
4、每次优化完成判断两种情况，满足的话将邻域像素添加到队列中当作种子点，并将当前像素patch的信息作为邻域像素的初始值：（1）邻域没有深度值（2）当前像素的置信度高于邻域像素一定范围（当前点因为重建优化完，置信度得到提升，和邻域拉开很大差距，说明这个区域很需要优化，有很大的优化空间）

  新建patch的可视图像和法向量的初始值等同于邻域patch
$$\begin{gathered} n(p^{{}^{\prime }})=n(p) \\ V(p^{{}^{\prime }})=V(P) \end{gathered}$$

深度值非线性优化

  假设参考视角为$R,k\in N$为与参考视角相关的局部邻域视角,$P_k$为三维点投影到第k个视角的投影矩阵. 每个像素提供了一个初始的深度和法向量和局部视角后，对深度值进行非线性优化。建立几何模型如下：

  其中，$o_R$表示参考视角的相机中心，(s, t)表示参考视角像素点,其初始深度为$h(s,t)$，已知相机姿态, 则对应的三维空间中的射线的单位向量为${\overrightarrow{r}}_R(s,t)$, 则该像素对应的三维空间的点为:
$$X_R(s,t)=O_R+h(s,t){\overrightarrow{r}}_R(s,t)$$

为了对像素点$(s,t)$的深度和三维中的法向量进行优化,以该像素为中心，建立nxn的patch, pathch对应三维中间中的一个很小的平面.同时引入两个变量来帮助表示patch中每个像素的三维坐标. 则像素$(s+i,t+j)$对应的深度为:
$$h(s+i,t+j)=h(s,t)+i\ast h_s(s,t)+j\ast h_t(s,t)$$

假设$(s+i,t+j)$处的射线方向近似为${\overrightarrow{r}}_R(s,t)$,则像素$(s+i,t+j)$对应的三维坐标为:
$$X_R(i,j)=O_R+{\overrightarrow{r}}_R(s,t)(h(s,t)+i\ast h_s+j\ast h_t)$$

光度模型：

  对朗伯反射进行简单假设,以提升颜色匹配的准确度. 为邻域内的每个视角K分配一个颜色尺度$c_k=[c_k^r,c_k^g,c_k^b{]}^T\in R^{3\times 1}$,如果深度$h(s,t)$估计准确，则有：
$$I_R(s+i,t+j)=c_k\cdot I_k(P_k(X_R(s+i,t+j)))$$

对所有的邻域视角和所有的patch中的patch中的三维点成立,其中$I_R(.)$表示取参考视角对应像素处的颜色(3通道,为向量形式), $I_k(.)$表示取第k个视角中图像上对应像素处的颜色,表示参考视角的一点和其他视角在对应同一点上的颜色一致.

结合几何模型和光度模型,可以得到深度值和法向量的优化数学模型如下：
$$E=\sum _{k\in N}\sum _{i=-\frac{n-1}{2}}^{\frac{n-1}{2}}\sum _{j=-\frac{n-1}{2}}^{\frac{n-1}{2}}[I_R(s+i,t+j)-c_k.I_k(P_k(X_R(s+i,t+j))){]}^2$$
  其中$X_R(s+i,t+j)$表示将参考图像中的$(s+i,t+k)$映射到三维空间中的坐标,乘以$P_k$表示将该点的三维坐标投影到第k个视角图像上的像素坐标,对这个像素坐标取$I_k$表示将第k个视角图像上这个点的像素颜色值取出来,再乘以前面的$c_k$,表示一个调节因子,增加鲁棒性

  为了书写方便$I_R(s+i,t+j)$简记为$I_R(i,j)$, 将$I_k(P_k(X(s+i,t+j)))$简记为$I_k(i,j)$,将$h_s(s,t),h_t(s,t)$分别简记为$h_s,h_t$于是,上式可以简记为:
$$E=\sum _{ijk}[I_R(i,j)-c_k.I_k(i,j){]}^2$$
$i,j\in [-\frac{n-1}{2},\frac{n-1}{2}]$是patch中的采样点，k表示视角的个数，需要优化的变量为三维坐标的深度$h(s,t),h_s,h_t$和颜色尺度因子 { c k } \{c_{k}\} {ck​}

颜色尺度的优化：

$$\begin{gathered} \frac{\partial E}{\partial c_k^r}=\sum _{ij}{r}_k^r(i,j)I_R^r(i,j)+c_k^r\sum _{ij}(I_k^r(i,j){)}^2 \\ \frac{\partial E}{\partial c_k^g}=\sum _{ij}{r}_k^g(i,j)I_R^g(i,j)+c_k^g\sum _{ij}(I_k^g(i,j){)}^2 \\ \frac{\partial E}{\partial c_k^b}=\sum _{ij}{r}_k^b(i,j)I_R^r(i,j)+c_k^b\sum _{ij}(I_k^b(i,j){)}^2 \\ {c}_k=[\frac{{\sum }_{ij}{I}_K^r(i,j)I_R^r(i,j)}{{\sum }_{ij}(I_K^r(i,j){)}^2},\frac{{\sum }_{ij}{I}_K^g(i,j)I_R^g(i,j)}{{\sum }_{ij}(I_K^g(i,j){)}^2},\frac{{\sum }_{ij}{I}_K^b(i,j)I_R^b(i,j)}{{\sum }_{ij}(I_K^b(i,j){)}^2}{]}^T \end{gathered}$$

颜色尺度的优化：

  分别对上式求三个偏导，并令结果为0，得：
$$\begin{gathered} c_k^r=\frac{{\sum }_{ij}{I}_K^r(i,j)I_R^r(i,j)}{{\sum }_{ij}(I_K^r(i,j){)}^2} \\ {c}_k^g=\frac{{\sum }_{ij}{I}_K^g(i,j)I_R^g(i,j)}{{\sum }_{ij}(I_K^g(i,j){)}^2} \\ {c}_k^b=\frac{{\sum }_{ij}{I}_K^b(i,j)I_R^b(i,j)}{{\sum }_{ij}(I_K^b(i,j){)}^2}{]}^T \end{gathered}$$

$h(s,t),h_s,h_t$的优化：

  待优化的模型为最小二乘的非线性优化问题，无法得到该函数的解析解，采用迭代的方式求取其近似解。首先引入$I_R(i,j)$关于$h(s,t),h_s,h_t$的线性表达式:
I R ( i , j ) ≈ c k ⋅ I k ( P k ( O R + r → R ( s , t ) ( h ( s , t ) + i ∗ h s + j ∗ h t ) ) ) + c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ( d h ( s , t ) + i ∗ d h s + j ∗ d h t ) = c k ⋅ I k ( i , j ) + c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ( d h ( s , t ) + i ∗ d h s + j ∗ d h t )

$$\begin{aligned} \boldsymbol{I}_{R}(i, j) \approx \boldsymbol{c}_{k} \cdot & \boldsymbol{I}_{k}\left(\boldsymbol{P}_{k}\left(\boldsymbol{O}_{R}+\overrightarrow{\boldsymbol{r}}_{R}(s, t)\left(h(s, t)+i * h_{s}+j * h_{t}\right)\right)\right)+\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\left(d h(s, t)+i * d h_{s}+j * d h_{t}\right) \\ =& \boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)+\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\left(d h(s, t)+i * d h_{s}+j * d h_{t}\right) \end{aligned}$$ IR​(i,j)≈ck​⋅=​Ik​(Pk​(OR​+r R​(s,t)(h(s,t)+i∗hs​+j∗ht​)))+ck​⋅∂h(s,t)∂Ik​(i,j)​(dh(s,t)+i∗dhs​+j∗dht​)ck​⋅Ik​(i,j)+ck​⋅∂h(s,t)∂Ik​(i,j)​(dh(s,t)+i∗dhs​+j∗dht​)​
由此,我们可以得到能量函数E的近似表达式为:
$$E=\sum _{ijk}{\left({\mathbf{I}}_R(i,j)-{\mathbf{c}}_k\cdot {\mathbf{I}}_k(i,j)+{\mathbf{c}}_k\cdot \frac{\partial {\mathbf{I}}_k(i,j)}{\partial h(s,t)}\left(dh(s,t)+i\ast d{h}_s+j\ast d{h}_t\right)\right)}^2$$

令:
b i j k = I R ( i , j ) − c k ⋅ I k ( i , j )    ∇ x = [ d h ( s , t ) , d h s , d h t ] T A i j k = ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) [ 1 i j ] T \boldsymbol{b}_{i j k}=\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)\\\;\\ \nabla \boldsymbol{x}=\left[d h(s, t), d h_{s}, d h_{t}\right]^{T}\\ \boldsymbol{A}_{i j k}=\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)\left[

$$\begin{array}{l} 1 \\ i \\ j \end{array}$$ \right]^{T} bijk​=IR​(i,j)−ck​⋅Ik​(i,j)∇x=[dh(s,t),dhs​,dht​]TAijk​=(ck​⋅∂h(s,t)∂Ik​(i,j)​) ​1ij​ ​T

则上式可以表达为:
$$E=\sum _{ijk}{\left(A_{ijk}\nabla x+b_{ijk}\right)}^2$$

求$E$关于$\nabla x$的梯度，可以得到：
$$\frac{\partial E}{\partial \nabla \mathbf{x}}=\sum _{ijk}{\mathbf{A}}_{ijk}^T\left({\mathbf{A}}_{ijk}\mathbf{x}+{\mathbf{b}}_{ijk}\right)$$

令$\frac{\partial E}{\partial \nabla x}=0$得：
$$\left(\sum _{ijk}{A}_{ijk}^T{A}_{ijk}\right)\nabla x+\sum _{ijk}{A}_{ijk}^T{b}_{ijk})=0$$

令$A={\sum }_{ijk}{A}_{ijk}^T{A}_{ijk},b={\sum }_{ijk}{A}_{ijk}{b}_{ijk}$，则有：
$$\nabla x=A^{-1}b$$

补充:
A i j k = ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) [ 1 i j ] T    A i j k T A i j k = [ 1 i j ] ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) T ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) [ 1 i j ] T = ∥ c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ∥ 2 [ 1 , i , j i , i 2 , i j j , i j , j 2 ]    A i j k T b i j k = [ 1 i j ] ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) T ( I R ( i , j ) − c k ⋅ I k ( i , j ) ) = s i j k [ 1 i j ]

$$\begin{array}{c} \boldsymbol{A}_{i j k}=\left(\boldsymbol{c}_{k} \cdot \frac{\partial I_{k}(i, j)}{\partial h(s, t)}\right)\left[\begin{array}{l} 1 \\ i \\ j \end{array}$$ \right]^{T}\\\;\\ \boldsymbol{A}_{i j k}^{T} \boldsymbol{A}_{i j k}=\left[ $$\begin{array}{l} 1 \\ i \\ j \end{array}$$ \right]\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)\left[ $$\begin{array}{l} 1 \\ i \\ j \end{array}$$ \right]^{T}=\left\|\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right\|^{2}\left[ $$\begin{array}{ccc} 1, & i, & j \\ i, & i^{2}, & i j \\ j, & i j, & j^{2} \end{array}$$ \right]\\\;\\ \boldsymbol{A}_{i j k}^{T} \boldsymbol{b}_{i j k}=\left[ $$\begin{array}{l} 1 \\ i \\ j \end{array}$$ \right]\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)\right)=s_{i j k}\left[ $$\begin{array}{l} 1 \\ i \\ j \end{array}$$ \right] \end{array} Aijk​=(ck​⋅∂h(s,t)∂Ik​(i,j)​) ​1ij​ ​TAijkT​Aijk​= ​1ij​ ​(ck​⋅∂h(s,t)∂Ik​(i,j)​)T(ck​⋅∂h(s,t)∂Ik​(i,j)​) ​1ij​ ​T= ​ck​⋅∂h(s,t)∂Ik​(i,j)​ ​2 ​1,i,j,​i,i2,ij,​jijj2​ ​AijkT​bijk​= ​1ij​ ​(ck​⋅∂h(s,t)∂Ik​(i,j)​)T(IR​(i,j)−ck​⋅Ik​(i,j))=sijk​ ​1ij​ ​​

其中$s_{ijk}={\left(c_k\cdot \frac{\partial I_k(i,j)}{\partial h(s,t)}\right)}^T\left(I_R(i,j)-c_k\cdot I_k(i,j)\right)$是标量，因此可以得到：
A = ∑ i j k A i j k T A i j k = ∑ i j k ∥ c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ∥ 2 [ 1 , i , j i , i 2 , i j j , i j , j 2 ] b = ∑ i j k A i j k T b i j k = ∑ i j k ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) T ( I R ( i , j ) − c k ⋅ I k ( i , j ) ) [ 1 i j ]

$$\begin{array}{l} \boldsymbol{A}=\sum\limits_{i j k} \boldsymbol{A}_{i j k}^{T} \boldsymbol{A}_{i j k}=\sum\limits_{i j k}\left\|\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right\|^{2}\left[\begin{array}{ccc} 1, & i, & j \\ i, & i^{2}, & i j \\ j, & i j, & j^{2} \end{array}$$ \right]\\ \boldsymbol{b}=\sum \limits_{i j k} \boldsymbol{A}_{i j k}^{T} \boldsymbol{b}_{i j k}=\sum\limits_{i j k}\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)\right)\left[ $$\begin{array}{l} 1 \\ i \\ j \end{array}$$ \right] \end{array} A=ijk∑​AijkT​Aijk​=ijk∑​ ​ck​⋅∂h(s,t)∂Ik​(i,j)​ ​2 ​1,i,j,​i,i2,ij,​jijj2​ ​b=ijk∑​AijkT​bijk​=ijk∑​(ck​⋅∂h(s,t)∂Ik​(i,j)​)T(IR​(i,j)−ck​⋅Ik​(i,j)) ​1ij​ ​​

算法流程：

1. 只进行$h(s,t)$的优化(迭代4次)
2. While(iteration <20)
3.   每间隔5次迭代，进行$h(s,t)，h_s,h_t$的优化，否则仅仅优化深度
4.   优化完成后跟踪判断每个视角的置信度
5.   如果视角的置信度太小，或者超过迭代14次尚未收敛，则从局部视角中移除该视角，并重新进行视角选择

---