OpenMVS中的PatchMatch算法笔记

本文主要记录多视立体几何中传统算法使用得较广的一个开源框架OpenMVS的一些实现流程、算法细节等，重点针对稠密重建部分，使用的是PatchMatch的方法，SGM的部分暂时不涉及。

基本流程

pipeline2.png

主要的数据结构

OpenMVS在整个稠密重建这一块，主要有几个比较大的数据结构，分别是Scene，DenseDepthMapData及DepthEstimator，PatchMatchCUDA。前两个为数据类，后两个为深度估计的操作类，分别对应CPU的实现和CUDA版本的实现。Scene主要是存储输入（相机内外参，图像，稀疏点云等）输出（稠密点云）的，DenseDepthMapData主要是存储中间过程的各个视角的深度图。

scene.png

DenseDepthMapData.png

具体模块介绍

Neighbor的选择：SelectViews

多视几何最核心的思想是确定同一个三维点在不同视角下的对应像素，再利用三角法可以估计出深度信息。因此，对于每个reference图像找到可以用来做重建的共视图像就是首要的一步。在OpenMVS中，有多个数据结构中有neighbor的出现，Scene中的neighbor包含所有共视的图像，而在DepthData中的neighbor则是从中做了一定挑选以后用来做重建的共视图像，可以选多个，也可以只选一个。

选择多个neighbor

neighbor的确认：在输入的稀疏点云中，每个点都存储了一个track，代表着哪些view可以看到它，于是对于当前视角而言，能共同看到三个以上三维点的视角即为neighbor。

neighbor的排序：针对每个neighbor，会有一个评分，这个评分的维度主要是有三个，一个是三角化的角度，一个是相对于reference视角的scale，另一个是两个视角的overlap。

对于稀疏重建中，当前视角与neighbor视角都能看到的三维点，根据对应点的关系，我们可以构造如图的三角化，对应的角度即为三角化的角度。

angle.png

我们对于这个角度值进行一个打分，我们设置最理想的角度为 $\theta_{opt} = 15$ 度，具体的公式如下：
$w_{\theta}= \left\{ \begin{matrix} \left(\frac{\theta}{\theta_{opt}}\right)^{1.5}, & \theta < \theta_{opt} \\ \exp\left(-\frac{(\theta - \theta_{opt})^2}{2(\theta_{opt}*1.3)^2}\right), & \theta \geq \theta_{opt} \end{matrix} \right.$

weight_angle.png

由 $w_{\theta}$ 的示意图像可以知道，倾向于选择接近理想角度的neighbor，同时对于小于理想角度的情形惩罚力度会比大于理想角度情形惩罚力度要大。这主要是说三角化只有达到一定的角度，基线才能够足够保证三维重建的可靠性及准确性。

而对于每个公共三维点，scale的定义为 $s = \frac{f_1}{f_2}\frac{d_2}{d_1}$ ，其中 $f$ 为相机焦距，而 $d$ 为该三维点在对应视角的深度。我们同样对scale进行打分，设置了一个期望区间 $[1,1.6]$ ，对处于该区间之外的scale进行惩罚，具体的公式为：
$w_{s}= \left\{ \begin{matrix} \left(\frac{1.6}{s}\right)^{2}, & s > 1.6 \\ 1.0, & 1.0 \leq s \leq 1.6 \\ s^2, & s < 1.0 \end{matrix} \right.$

weight_scale.png

同样，可以看到 $w_s$ 主要惩罚那些scale大于1.6及小于1.0的情形。

两个视角的overlap值 $w_{a}$ ，计算方式为找到两个视角下overlap的三维点，然后投影到reference视角，此时将refrence视角的图像化成16*16的格子，按投影点算像素的占有率，这个比例即当做neighbor视角与reference视角的overlap值。

综合上述三个维度，score的最终计算公式为
$score = w_a\cdot\sum_{common~nodes} w_{\theta}\cdot w_{s}$

根据score从大到小排序，我们可以截断选择成绩最好的若干共视图像用来重建。

选择最优neighbor

基于上面得到的共视图像集合，我们也可以选择最优neighbor来进行重建。选择最优neighbor的策略比较复杂，本质上是在求解一个马尔科夫随机场的优化问题。我们定义view的编号为 $x_i \in [0, N-1]$ ，view $x_i$ 对应的最优neighbor为 $y(x_i)\in[0, N]$ ，这里需要注意的是当 $y$ 取 $N$ 时，代表这个view没有找到适合的neighbor，它在接下来的重建中也将会被抛弃。求解最优neighbor其实是在做一个全局优化，优化的代价函数为
$\arg\min\limits_y ~E_{data}(x_i, y(x_i)) + E_{pairwise}(x_i, x_j)$
其中， $E_{data}$ 描述的是当view $x_i$ 选中某一个 $y$ 的代价，而 $E_{pairwise}$ 描述的则是，当 $x_j$ 是 $x_i$ 的其中一个neighbor时，两者对应的最优neighbor相关性的一个代价。总体的原则是，view $x_i$ 选中 $y$ 的自身代价要小，同时还要尽量避免它和它的neighbor互为最优选择。这两个代价的具体形式可以描述为
$E_{data}(x_i, y(x_i))= \left\{ \begin{matrix} avgScore / score(x_i, y(x_i)), & y(x_i)\neq N \\ C_{none} * avgScore*r_{none}, & y(x_i)=N~且~x_i没有neighbor\\ C_{none} * E_{data}(x_i, y(neighbor(x_i))), & y(x_i)=N~且~x_i有neighbor \end{matrix} \right.$

$E_{pairwise}(x_i, x_j) = \left\{ \begin{matrix} C_{penalty}, & 当~y(x_i)=x_j~且~y(x_j)=x_i \\ C_{empty} + r*\frac{w_a(x_i, x_j)}{w_a(x_i, y(x_i))}, & 当y(x_i)\neq N ~且~y(x_j)=N\\ C_{empty} + r*\frac{w_a(x_i, x_j)}{w_a(x_j, y(x_j))}, & 当y(x_i)= N ~且~y(x_j)\neq N\\ r*\left(\frac{w_a(x_i, x_j)}{w_a(x_i, y(x_i))} + \frac{w_a(x_i, x_j)}{w_a(x_j, y(x_j))}\right), & 其他情形 \end{matrix} \right.$

其中，avgScore对应了 $x_i$ 与所有neighbor的score值的平均，而 $C,r$ 等都为预设的常数，默认参数中 $C_{none} = 6.0, r_{none} = 0.01, C_{penalty}=7.2, C_{empty} = 2.4, r = 0.3$ 。 $E_{data}$ 项比较好理解， $E_{pairwise}$ 项其实是惩罚两个view互为最优neighbor的情形，也同时惩罚找不到最优neighbor的情形。

深度初始化：InitViews

初始化的时候主要干了如下几件事情：

按照之前计算的scale值，将neighbor对应的image重采样至对应大小，同时更新相应的相机参数；
去除一些score值过小的neighbor，这个阈值是动态设定的，默认的是排名第一的neighbor对应的score值的3%；
初始化当前视角的深度及法向信息

这里重点介绍第3件事情。深度初始化的目的是要给reference图像的每个像素找到一个初始的深度值及法向值，这里利用的在二维区域上的Delaunay分割来实现的。我们利用稀疏重建的点云，将其往图像的二维区域上投影，通过Delaunay分割可以将图像区域做一个三角剖分，同时保证上述二维投影点为三角形节点。如果加上图像区域的4个角，则可以实现整个图像区域的三角剖分。实际上，这些投影点都对应有深度信息，我们通过把投影点按深度进行提升，就可以把二维区域的三角剖分转化成一个三维的表面网格。在这里，由于图像区域的4个角不一定对应稀疏点云中的点，我们可以用平均深度来当做其深度。有了上述表面网格以后，我们可以估计出节点位置的法向方向。同时对于二维区域的任何一个像素，可以往该表面网格竖直投影找到其三维位置，即可得到初始的深度值，其实也就是找到该像素所在的三角形，利用线性插值可以得到深度值。而法向值则对应为该表面网格的三角片的法向方向。

init.png

深度估计CPU版本

深度估计的CPU版本主要是实现了如下参考文献的算法

Shen, Shuhan. "Accurate multiple view 3d reconstruction using patch-based stereo for large-scale scenes." IEEE transactions on image processing 22.5 (2013): 1901-1914.

像素处理的顺序

每个像素点都需要估计深度，CPU版本中，像素的处理时串行的，预先订制好了一个处理的顺序。在有些操作中，深度估计依赖于邻域的一些信息，所以，这个顺序对后续的处理还是有影响的。

zigzag.png

处理时，主要按照上图中所示的顺序来进行，当该像素点位于mask之外时，则会被忽略。

匹配代价的计算：ScorePixel

给定一个reference图像上的像素，patch对应的图像邻域大小，其对应的深度，对应的法向等，即可找到它在所有neighbor图像上的对应像素及对应patch，然后可以计算patch与patch之间的匹配代价，最后综合所有的neighbor的匹配代价可以得到一个最终的聚合匹配代价。需要注意的是，reference图像上的patch我们可以取成 $n\times n$ 的方形邻域，但经过变换以后，方形邻域每个像素都将在neighbor图像上找到对应点，这些对应点围成的patch就不再是方形的了，这是与基于极线校正方法计算匹配代价的比较大的不同。

按照上述说法，整个计算将分成几部分

给定reference图像的像素，及其深度，法向信息，如何确定变换矩阵（单应矩阵） $H$ ，找到匹配像素？
单个neighbor图像的匹配代价值是如何计算的？
聚合匹配代价是如何计算的？

单应矩阵 $H$ 的确定

给定两个视角，其投影矩阵（将世界坐标系的点投影到相机归一化平面坐标系的矩阵）分别为 $P_0 = [I|0]$ 和 $P_1 = [A|a]$ ，对于空间中的一个平面 $\pi^T X = 0$ ，其中 $\pi = (v^T, 1)^T$ ，则该平面上任意一点在两个相机上的投影之间的关系可以由一个单应矩阵来确定 $x'=Hx$ ，这个单应矩阵在齐次坐标下可以表达为^[1]
$H=A-av^T$
上式的推导如下。取平面上任一三维点 $X(X',Y',Z')$ ，则其在第0个相机下的投影点 $x_0(x_0',y_0')$ 满足

$\left( \begin{matrix} x_0'\\ y_0'\\ 1 \end{matrix} \right) = \rho_0 P_0 \left( \begin{matrix} X'\\ Y'\\ Z'\\ 1 \end{matrix} \right) = \rho_0 \left( \begin{matrix} X'\\ Y'\\ Z' \end{matrix} \right)$
点 $X$ 在平面上，所以有 $v^TX + 1=0$ ，用 $x_0$ 代入，则有 $\rho_0 = - v^T(x_0,1)$ .

该点在第一个相机下的投影点 $x_1(x_1', y_1')$ 满足
$\left( \begin{matrix} x_1\\ 1 \end{matrix} \right) = \left( \begin{matrix} x_1'\\ y_1'\\ 1 \end{matrix} \right) = \rho_1 P_1 \left( \begin{matrix} X'\\ Y'\\ Z'\\ 1 \end{matrix} \right) = \rho_1 \left( AX+a \right) = \frac{\rho_1}{\rho_0} \left(A\left(\begin{matrix} x_0\\ 1 \end{matrix} \right)+ a \rho_0\right) = s (A-av^T)\left( \begin{matrix} x_0\\ 1 \end{matrix} \right)$
这里的 $s$ 代表的是该三维点在两个相机坐标系下深度的比例。

在当前的场景中，我们取reference相机为0号相机，而neighbor相机为1号相机，不考虑内参矩阵 $K$ 的情况下，我们有
$P_0= \left( \begin{matrix} R_0 & - R_0C_0\\ 0 & 1 \end{matrix} \right) ,~~ P_1= \left( \begin{matrix} R_1 & - R_1C_1\\ 0 & 1 \end{matrix} \right)$
如果我们通过变换使得 $P_0 = [I|0]$ ，则 $P_1$ 可以写成
$\begin{eqnarray*} P_1&=& \left( \begin{matrix} R_1 & - R_1C_1\\ 0 & 1 \end{matrix} \right) \left( \begin{matrix} R_0 & - R_0C_0\\ 0 & 1 \end{matrix} \right)^{-1} \\ &=& \left( \begin{matrix} R_1 & - R_1C_1\\ 0 & 1 \end{matrix} \right) \left( \begin{matrix} R_0^{-1} & C_0\\ 0 & 1 \end{matrix} \right)\\ &=& \left( \begin{matrix} R_1R_0^{-1} & R_1(C_0-C_1)\\ 0 & 1 \end{matrix} \right)\\ &=& [R_1R_0^{-1}|R_1(C_0-C_1)] \end{eqnarray*}$
经过点 $X$ 的patch平面可以定义为
$n^TX+d=0$
其中， $n$ 为点 $X$ 的法线，而 $d$ 则代表了世界坐标系原点到该平面的距离，与上面平面的定义对比可知， $v=-\frac{n}{n^TX}$ .

根据上面的结果，可以得到 $P_0$ 与 $P_1$ 之间的单应矩阵为
$H=R_1R_0^{-1}-R_1(C_0-C_1)\cdot \frac{n^T}{-n^TX} =R_1R_0^{-1}+ \frac{R_1(C_0-C_1)\cdot n^T}{n^TX}$

再引入内参矩阵变可得到最终的单应矩阵
$H = K_1\left(R_1R_0^{-1}+ \frac{R_1(C_0-C_1)\cdot n^T}{n^TX}\right)K_0^{-1}$

这里的 $H$ 是没有考虑深度比例 $s$ 的，因为 $s$ 的值与具体的点是有关的。在OpenMVS的实现中， $s$ 相关的操作十分隐蔽，对于给定像素 $x=(x_0, y_0)^T$ ，其对应点应该满足
$H\left( \begin{matrix} x_0\\ y_0\\ 1 \end{matrix} \right) = \left( \begin{matrix} x_1\\ y_1\\ z_1 \end{matrix} \right) \Rightarrow \left( \begin{matrix} x_1/z_1\\ y_1/z_1\\ 1 \end{matrix} \right)$
上述箭头对应的操作，OpenMVS中是通过Point2f的构造函数实现的

explicit inline TPoint2(const cv::Point3_<TYPE>& pt) : Base(pt.x/pt.z,pt.y/pt.z) {}

匹配代价的计算

这里的匹配代价是通过ZNCC（zero-normalized cross correlation）来计算的。按照ZNCC的定义有
$ncc = \frac{\sum\limits_{q\in B} (I_q-\overline{I_q}) (I_{H(q)} - \overline{I_{H(q)}})}{\sqrt{\sum\limits_{q\in B} (I_q-\overline{I_q})^2 * \sum\limits_{q\in B} (I_{H(q)}-\overline{I_{H(q)}})^2 }}$
这里的 $B$ 代表在reference图像空间的邻域， $I$ 代表该像素的灰度值， $\bar{I}$ 代表邻域范围内的平均值。在具体实现时，这里的平均值其实是一个加权平均，权重为一个以中心像素的灰度及位置为中心的一个高斯分布
$w_q= \exp(w_c + w_s) = \exp\left(-\frac{(I_{ij}-I_c)^2}{2*\sigma_c^2} - \frac{i^2 +j^2}{2*\sigma_s^2}\right)$

$\overline{I_q} = \frac{\sum\limits_{q\in B} I_q * w_q}{\sum\limits_{q\in B} w_q}$

需要注意的是，计算 $\overline{I_{H(q)}}$ 时使用的权重仍然是 $w_q$ 。但从代码的实现上，匹配代价中ZNCC部分的计算并不完全遵循上述公式，具体可写成
$ncc=\frac{\sum\limits_{q\in B} (I_q-\overline{I_q})*I_{H(q)}*w_q}{\sqrt{\sum\limits_{q\in B}w_q * (I_q - \overline{I_q})^2}\sqrt{\sum\limits_{q \in B} w_q*(I_{H(q)}^2-\overline{I_{H(q)}}^2)}}$
ZNCC的意义是说两个patch越相近，值越大，而我们定义的匹配代价希望是越小，代表匹配效果越好，因此从ZNCC转换成匹配代价的过程中我们先定义 $score = 1-ncc$ 。同时，这个score值还可以进行多个维度的加权。第一个是光滑度的考量，主要是看当前像素与其四邻域的像素深度值及法向值是否差异不大，具体的公式如下
$w_d = \exp\left(-\frac{(d_p/d)^2}{2\sigma_d^2}\right),~~w_n=\exp\left(-\frac{\delta\theta^2}{2\sigma_n^2}\right)，~~w_{smooth}=\prod_{i=1}^4(1-c_d*w_d)(1-c_n*w_n)$
其中， $c_d,c_n,\sigma_d,\sigma_n$ 都为设定好的常数。另一个维度便是一致性的考量，即reference图像与neighbor图像对应点的深度一致性，这个一致性是通过将neighbor图像对应点的深度还原成世界坐标系的三维点，再投影回reference图像，计算投影值与原像素值的差异来描述的。具体公式为
$dist = \sqrt{\delta x^2 + \delta y^2}\\ w_{con} = c_{con}\min(\sqrt{dist*(dist+2.0)}, 4.0)$
综合起来，最终的匹配代价表达式为
$score=(1-ncc)*w_{smooth} + w_{con}$
ncc的取值范围是 $[-1,1]$ ，因此理论上，score的值最小为0，最大为 $2+\max(w_{con})$ ，而实际上我们会给score设一个阈值，默认设置是0.55，高于这个阈值的都会被抛弃，所以最后得到的score将在区间 $[0, 0.55]$ 。

聚合匹配代价的计算

该像素的聚合匹配代价，可以由每个neighbor的匹配代价平均而来，也可以直接使用最小的匹配代价。在具体实现中，当neighbor数小于2时，选取的是最小的匹配代价，而大于2时，在去掉一些过大的匹配代价值后再平均。

传播与随机优化：ProcessPixel

在传播优化时，像素的处理顺序一种是按照前节描述的那样从左上到右下，另一种是反过来从右下到左上。上述操作我们迭代运行，一次左上到右下，一次右下到左上，总共3次。对于每个待考察的像素，我们查看其四邻域的深度和法向值，以邻域像素的值来重新估计待考察像素的匹配代价，如果匹配代价变小，则用邻域像素的值来取代。实际操作中，并不是直接使用了邻域的深度，而是使用了邻域像素的patch平面，取当前像素的三维点与相机原点的连线与该patch平面的交点作为当前像素新的三维点位置。

interpolation.png

如上图，假设当前像素为 $P'$ ，传播之前的相机坐标系位置为 $X'$ ，对应法向 $n'$ 和深度 $d'$ ，传播的像素为 $P$ ，相机坐标系位置为 $X$ ，对应法向 $n$ 和深度 $d$ ，则传播后 $P'$ 对应的位置为 $CX'$ 与patch平面的交点 $X''$ ，法向为 $n$ ，深度为 $d''$ 。新的深度的计算公式可推导如下：
$\begin{eqnarray*} n^T(X''-X) &=& 0\\ n^T(d''K^{-1}P'-dK^{-1}P) &=& 0\\ n^TK^{-1}(d''P'-dP)&=&0 \\ (n_x, n_y, n_z)* \left( \begin{matrix} \frac{1}{f} & & -\frac{c_x}{f}\\ & \frac{1}{f} & -\frac{c_y}{f}\\ & & 1 \end{matrix} \right) * \left( \begin{matrix} d''x' - dx \\ d''y' - dy \\ d'' - d \end{matrix} \right) &=& 0 \end{eqnarray*}$
求解这个线性方程组可以得到
$d'' = d * \frac{\frac{n_x}{f}(x-c_x) + \frac{n_y}{f}(y-c_y) + n_z}{\frac{n_x}{f}(x'-c_x) + \frac{n_y}{f}(y'-c_y) + n_z}$
在进行完邻域传播以后，还有一步是随机优化操作。基本上按照score的值按大小分成三挡，score越大的随机范围越大，在最大的那一档中，基本上就是完全随机找一个深度和法向与当前去比，而在前两档中，则是给定不同的随机范围，并随着迭代次数增加而收紧这个范围。

深度估计CUDA版本

OpenMVS里CUDA版本的深度估计算法主要是基于如下参考文献

Xu, Qingshan, and Wenbing Tao. "Multi-view stereo with asymmetric checkerboard propagation and multi-hypothesis joint view selection." arXiv preprint arXiv:1805.07920 (2018).

算法的基本流程与CPU版本类似，都是先计算匹配代价，然后再传播优化，随机调整，但是具体步骤和公式上有很多不同之处。

匹配代价的计算

对于给定的reference像素，给定neighbor图像，给定深度和法向，匹配代价可以表示成
$score = 1-\frac{\sum\limits_{q\in B}w_qI_qI_{H(q)}*\sum\limits_{q\in B}w_q - \left(\sum\limits_{q\in B}w_qI_q\right)\left(\sum\limits_{q\in B}w_qI_{H(q)}\right)} {\sqrt{\left(\sum\limits_{q\in B}w_qI_q^2\right)*\left(\sum\limits_{q\in B}w_q\right)-\left(\sum\limits_{q\in B}w_qI_q\right)^2}*\sqrt{\left(\sum\limits_{q\in B}w_qI_{H(q)}^2\right)*\left(\sum\limits_{q\in B}w_q\right)-\left(\sum\limits_{q\in B}w_qI_{H(q)}\right)^2}} +w_{con}$
其中所有的记号都与CPU版本中的公式一样。每个像素匹配代价的计算都可以并行完成。而在计算聚合匹配代价时，与CPU版本的简单平均或是取最小值不同，CUDA版本有一套很复杂的加权方式，并且这个加权方式与传播融合在一起，将在下节介绍。

传播优化

在CPU版本中，传播优化时像素是顺序进行的，即后一个像素的处理是依赖于上一个像素的结果。为了能够并行地进行传播，这里采用了棋盘格红黑像素分隔的方法，具体的分法可以查看下图，黑色像素的深度法向更新只依赖于邻域内的红色像素，而红色像素的深度法向更新只依赖于邻域内的黑色像素。在每一轮的传播中，先是固定红色像素的取值，并行更新所有的黑色像素，然后再固定黑色像素的取值，并行更新所有的红色像素。

checkboard.png

对于顺序传播而言，每个像素的更新有可能从较远的地方传播而来，而对于并行传播而言，每个像素只能接受邻域像素的传播，因此为了提高传播的效率，并行传播时的邻域会选取得更大，而不仅限于上下左右的四邻域。具体的实现中，邻域的范围如下图所示

cuda_neighbor.png

传播过程中，邻域像素被分成了8条路径，每条路径对应上图中的一种颜色，需要注意的是，在实现中，淡色的直线路径不只4个像素，而是11个像素。对于每条路径，我们会挑选出confidence值最大的那个像素当做这条路径的候选者（论文中称为hypothesis），我们的目标就是确定这8个候选像素的深度法向值能否传播给当前像素。对于每个hypothesis而言，我们可以通过插值得到当前像素的新的深度和法向，新的值对应每个neighbor图像都会有一个匹配代价，我们需要将其加权平均得到最终这个hypothesis的聚合匹配代价，传播是否成功就取决于聚合匹配代价与当前像素的取值计算得到的聚合匹配代价相比是否变小。

下面我们来介绍权重的计算。如果一个neighbor图像对应的权重大，我们认为它起决定性作用的概率更大。如果我们能得到这样一种概率分布，也就大概能得到neighbor图像的权重了。假设有 $N$ 个neighbor图像，它们被选中的概率为 $P_0, P_1, ..., P_{N-1}$ ，满足 $\sum\limits_{i=0}^{N-1}P_i=1$ ，实现中并不是直接以这个概率来作为权重的，而是通过均匀分布采样的方式来确定。如下图，不同的概率布满了区间 $[0,1]$ ，线段的长度代表着概率的大小，我们通过 $[0,1]$ 区间的均匀分布产生了32个点，然后看这32个点落在哪个概率区间，权重即是落在该区间点的个数除以32。从大数定理的角度来看，只要采样点的数目足够多，得到的权重其实将趋于上述概率值，但通过采样的手段，可以增加一点随机性，比如，没有被sample到的图像在下一轮传播中会被剔除neighbor集合。

prob.png

问题现在转化成如何确定每个neighbor图像选中的概率。这个概率与两个因素有关，第一个是一致性，即该neighbor图像是否也是当前像素的四邻域像素的neighbor图像，第二个是与8个hypothesis在该neighbor图像下的匹配代价有关。这两个因素分别会产生两个系数 $w_1,w_2$ ，而每个neighbor图像对应的概率为 $p_j = w_{1j}*w_{2j}$ 归一化后的值。

因素1： $w_{1j}=0.9n+0.1(4-n)$ ，这里的 $n$ 代表四邻域像素中共享当前neighbor图像的个数，最大为4。

因素2： $w_{2j} = \sum\limits_{i=0}^K \exp\left(-\frac{c_{ij}^2}{2\sigma^2}\right)/K$ ，其中 $c_{ij}$ 代表第 $i$ 个hypothesis在第 $j$ 个neighbor图像上的匹配代价。需要注意的是并不是8个hypothesis都要参与上述平均，即 $K$ 不一定为8。实现中，对 $c_{ij}$ 设置了阈值，并且这个阈值是与迭代次数相关的，具体来说，只有 $c_{ij} < 0.8*\exp(-\frac{iter^2}{2*4^2})$ ，这个hypothesis才会参与平均，总体原则是说随着迭代次数增多，阈值收紧，更倾向于选择score较小的几个去确定概率。

在每一轮黑色或者红色像素做完空间传播以后，还是需要做几次随机的调整，以期获得更好的深度和法向信息。

深度优化

这一步主要是两个操作，一个去除孤立的小团体（RemoveSmallSegments），一个是填补小洞（GapInterpolation）。

RemoveSmallSegments:从当前像素出发，从四邻域生长，只要深度差异在一定范围内，即合并入小团体，直至无法生长。检测小团体的点数，当点数小于阈值（预设100）时，则该小团体中的所有点将被删除。

GapInterpolation:逐行扫描，如果有相隔小于（默认7）个像素没有深度值，则启动填补小洞操作，深度上取小洞边界两个像素的深度值线性插值，法向按角度线性插值，信任值则去两个端点最小的。逐列扫描，进行同样的操作。

深度滤波

这里的深度滤波近似于在多个视角上去做平滑，对于reference图像上每一个二维像素点，自身有一个重建的深度，而其在neighbor图像上的对应点，也有一个对应的深度，这个深度是针对neighbor图像相机坐标系的，通过变换，我们可以得到它相当于reference图像相机坐标系的深度，我们可以将reference图像的深度与neighbor图像变换得来的深度进行加权平均，作为滤波后的深度。实现时，先将neighbor图像的深度图转换到reference图像相机坐标系的深度图，对于每一个像素，查看深度深度差异，只有大多数neighbor的深度差异在一定合理范围内时，才会利用这些深度参与加权平均，而权重即是深度对应的confidence值。此时，新的confidence值则等于参与平均的那些深度对应的confidence值的和减去没参与平均的那些深度对应的confidence值的和。当大多数neighbor的深度与当前深度不一致的时候，则会将该点删除，将深度置0。

深度融合

fuse.png

整个融合的过程大体可以由上图来解释，首先每个视角对应的深度图经由该视角的相机内参及位姿，可以得到世界坐标系的三维坐标。然后对于neighbor之间的对应点，我们考察他们在reference相机坐标系下的深度，如图所示， $N_1, N_2$ 的深度比reference像素的深度要大，被认为是处于遮挡状态，会被删除掉，做法就是将对应像素在其相机坐标系的深度图中置0； $N_4$ 的深度与reference深度相差不大，最终它会被用来平均reference的深度； $N_3$ 上述两个条件都不满足，会被暂时保留。而对于reference像素本身，如果和他一致的neighbor数（像 $N_4$ 这样的点）很少，那么它也将被删除掉，深度置0。

Multiple View Geometry in Computer Vision, 2nd ed, p326, Result 13.1 ↩

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 212,294评论 6赞 493
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 90,493评论 3赞 385
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 157,790评论 0赞 348
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 56,595评论 1赞 284
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 65,718评论 6赞 386
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 49,906评论 1赞 290
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,053评论 3赞 410
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 37,797评论 0赞 268
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,250评论 1赞 303
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 36,570评论 2赞 327
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 38,711评论 1赞 341
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 34,388评论 4赞 332
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,018评论 3赞 316
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,796评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,023评论 1赞 266
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 46,461评论 2赞 360
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 43,595评论 2赞 350

OpenMVS中的PatchMatch算法笔记

基本流程

主要的数据结构

具体模块介绍

Neighbor的选择：SelectViews

选择多个neighbor

选择最优neighbor

深度初始化：InitViews

深度估计CPU版本

像素处理的顺序

匹配代价的计算：ScorePixel

单应矩阵的确定

匹配代价的计算

聚合匹配代价的计算

传播与随机优化：ProcessPixel

深度估计CUDA版本

匹配代价的计算

传播优化

深度优化

深度滤波

深度融合

推荐阅读更多精彩内容

单应矩阵 $H$ 的确定