基于图像的三维重建——深度图计算方法1-SGM/tSGM(9)

立 体 算 法 流 程 介 绍

>全局立体匹配算法

  全局立体匹配算法主要是采用了全局的优化理论方法估计视差，建立全局能量函数，通过最小化全局能量函数得到最优视差值。通过二维相邻像素视差之间的约束(如平滑性约束)而得到更好的匹配效果，但是对内存的占用量大，速度慢不适合实时运行。主要的算法有图割(graph cuts)、置信传播(belief propagation)、动态规划等算法。

>局部立体匹配算法

  主要是采用局部优化方法进行视差值估计，局部立体匹配算法有SAD,SSD等算法，与全局立体匹配算法一样，也是通过能量最小化方法进行视差估计，但是在能量函数中，只有数据向，而没有平滑项。该算法由于每个像素计算互不干扰可以并行计算，所以可以实时。但由于所基于的局部窗口视差相同的假设在很多情况下并不成立导致匹配效果较差。

>半全局立体匹配算法SGM

  综合上述局部和全局算法的优缺点，半全局算法依旧采用全局框架，但是在计算能量函数最小化的步骤时使用高效率的一维路径聚合方法来代替全局算法中的二维最小化算法，使用一维最优来近似二维最优，得到的视差图在效果上和全局算法没有太大的差别，但是算法效率却有非常大的提升。

SGM 算 法 原 理 详 解

>深度图初始化

  深度图初始化主要有两种方法：1.Delaunay三角剖分；2.稀疏点投影到对应帧，得到稀疏深度图，然后三角化，将三角化内平面插值得到完整的深度图。
  1.Delaunay三角剖分
  首先，为什么用Delaunay三角剖分，因为对于给定点云的地形或其他对象建模，Delaunay三角剖分提供了一组很好的三角形，可以用作模型中的多边形。特别是，Delaunay三角剖分避免了狭窄的三角形(因为从面积上相比，它们有较大的外接圆。)下图是非Delaunay三角剖分与Delaunay三角剖分的对比：

  接下来介绍Delaunay三角化的三个重要特性：空圆特性、最大化最小角、对偶性。空圆特性是指点集P的Delaunay三角剖分满足P内任意一个点都不在P内任意一个三角面片的外接圆内。其中任意三个点的外接圆不包含第4个点即任意三个点的外接圆是空的；Delaunay三角剖分中的所有三角形都满足空圆特性；若四边形中的一个三角形不满足空圆特性，则将对角线互换之后构成的新的三角形一定满足空圆特性。空圆特性表现形式如下图所示：

  最大化最小角指的是Delaunay三角剖分最大化三角面片内三角形的最小角，Delaunay三角剖分中所有三角形的内角中最小角大于任何其他剖分的最小角，具体表现形式如下图所示：

  对偶性指的是Delaunay 三角剖分与Voronoi Diagram是互为对偶的关系，关系如下图所示：

  几种Delaunay生成算法时间复杂度的比较如下表所示：
  

  以Bowyer-Waston算法为例简单介绍一下Delaunay三角剖分的原理。Bowyer-Waston算法是一种增量算法，它的工作原理是向所需点子集的有效Delaunay三角剖分添加点，一次添加一个点。每次插入后，其外接圆包含新点的所有三角形都会被删除，留下一个星形多边形孔，然后使用新点重新对其进行三角剖分。通过利用三角剖分的连通性来有效定位要删除的三角形，该算法可以采用时间复杂度为O(N logN)的操作来对N个点进行三角剖分，尽管存在特殊的退化情况，在这种情况下，这会上升到O(N²)。下图是它的图形化解释：

  按我的理解就是说：在依次插入点时有两种情况，情况1就是新加入的点位于三角形内部，假设添加新的点p_r,连接p_rp_i,p_rp_j,p_rp_k，分别检查与三角形p_rp_ip_k,p_rp_jp_k,p_rp_ip_k相邻的三角形的空圆特性，如果不满足空圆特性，则进行翻转(即选择另一条对角线)，若p_jp_k需要翻转，则翻转之后需要对p_kp_m以及p_kp_j所在的四边形进行空圆特性检查。情况2就是新加入的点位于一条边上，假设添加新的点p_r,连接p_rp_l,p_rp_k,分别检查与三角形p_rp_jp_k,p_rp_jp_l,p_rp_lp_i,p_rp_kp_i相邻的三角形的空圆特性，如果不满足空圆特性，则进行翻转，若p_jp_k需要进行翻转，则翻转之后需要对p_kp_m以及p_kp_j所在的四边形进行空圆特性检查。

C++代码实现网址：https://github.com/Bl4ckb0ne/delaunay-triangulation，伪代码如下所示：

function BowyerWatson (pointList)
   // pointList is a set of coordinates defining the points to be triangulated
   triangulation := empty triangle mesh data structure
   add super-triangle to triangulation // must be large enough to completely contain all the points in pointList
   for each point in pointList do // add all the points one at a time to the triangulation
      badTriangles := empty set
      for each triangle in triangulation do // first find all the triangles that are no longer valid due to the insertion
         if point is inside circumcircle of triangle
            add triangle to badTriangles
      polygon := empty set
      for each triangle in badTriangles do // find the boundary of the polygonal hole
         for each edge in triangle do
            if edge is not shared by any other triangles in badTriangles
               add edge to polygon
      for each triangle in badTriangles do // remove them from the data structure
         remove triangle from triangulation
      for each edge in polygon do // re-triangulate the polygonal hole
         newTri := form a triangle from edge to point
         add newTri to triangulation
   for each triangle in triangulation // done inserting points, now clean up
      if triangle contains a vertex from original super-triangle
         remove triangle from triangulation
   return triangulation

  2.稀疏点投影到对应帧，得到稀疏深度图，然后三角化，将三角化内平面插值得到完整的深度图（具体形式如下图所示）

>匹配代价计算

  文献[1]中介绍的SGM的代价计算是基于互信息(Mutual Information,MI)的匹配测度计算算法来计算匹配代价，互信息是一种影像明暗变化不敏感的相关性测度。但由于原理复杂且计算需要迭代效率比较低，在实际应用中，更有效的方法如Census变换，故在此不再介绍MI。
  Census变换主要是使用像素邻域内的局部灰度差异将像素灰度转换为比特串即为census值，基于Census变换的匹配代价计算方法是计算左右影像对应的两个像素的Census变换值得汉明(Hamming)距离(两个比特串得对应位不相同的数量：先进性异或运算，再统计运算结果中1的个数)。

>代价聚合

  采用全局立体匹配算法，即找到每个像素的最优视差使得整体能量最小。能量方程如下：
  
  这是一个二维优化问题，为了提高优化效率，SGM将问题转化为使用一维路径聚合的方式来近似二维最优，提高效率的同时也保证了效果。SGM能量方程公式如下：
  
  式中，C为匹配代价，第二项和第三项是平滑项，我们期望的是视差是连续的。所以如果当前像素x_b与邻域像素x_N视差相差比较小(1个像素)我们会给一个比较小的惩罚P₁,如果大于一个像素则给一个大的惩罚P₂,但是实际场景中肯定会有一些视差不连续区域相差比较大(比如：前景和背景),如下图所示：

  为了处理这种情况，我们对P₂进行调整,公式如下所示，如果像素和它的邻域像素亮度差很大，那么该像素很可能是位于视差非连续区域，则一定程序上允许其和邻域像素的视差差值超过1个像素，对应超过1个像素的惩罚力度就适当减小一点。
  
  具体求解过程中，SGM路径代价聚合的思路就是将像素所有视差下的匹配代价进行像素周围所有路径上的一维聚合得到路径下的路径代价值，然后将所有路径代价值相加得到该像素聚合后的匹配代价值。一维聚合路径示意图如下所示(图中各个箭头方向就是聚合各个路径)：

  某一维路径r代价计算公式如下：
  
  式中L_ri(x_b,d) < C_max(x_b,d) + P₂ ;最终代价值(所有路径之和)为：
  

>视差计算

  代价聚合后，每个像素直接找聚合代价最小对应的视差值就是我们所要求的视差值。

>视差优化

  视差一致性检查
  将左右图互换，得到R-L视差图。与L-R视差图对比，根据左图的视差找到在右图中的对应视差，如果两者小于阈值则认为是准确的，反之是错的把该值剔除。如下图：

  亚像素计算
  上述计算的视差图是像素级别的，在实际使用中精度是无法满足需求的。SGM通过在上述计算的最小代价的视差层附近进行插值找到亚像素级的精度。如下图：

tSGM 原 理 介 绍

  与SGM的区别主要是在代价聚合的时候：1.使用金字塔图像计算视差(由粗糙到精即从低分辨率到高分辨率计算匹配代价)；2.每个像素的视差范围都不同，只在真值附近搜索大大减少了搜索范围和内存占用，如下图：

  每个像素的视差范围计算方法：如果当前像素值d(x)无效，搜索窗口win=31,反之win=7;以当前像素为中心，窗口大小win搜索所有有效视差都存储在disps中，disps的大小是numdisp,中值是disp，最大值是max，最小值min：


  3.代价聚合时，SGM由于每个像素视差范围都一样，所以各个搜索路径都有对应的视差，tSGM由于每个像素视差范围不一样有的可能都没有重叠范围，所以之前的代价聚合计算方法需要调整，把没有重合的情况考虑进去。

if d > dmax(xb - ri)
  Lri(xb,d) = Cri(xb,d) + P2;
if d < dmin(xb - ri)
  Lri(xb,d) = Cri(xb,d) + P2;
else
  Lri(xb,d) = SGM

参 考 文 献 与 代 码 整 体 框 架

  参考文献：[1]. Stereo Processing by Semi-global Matching and Mutual Information
       [2]. SURE: Photogrammetric Surface Reconstruction from Imagery
  代码框架：

文章参考：
1.https://en.wikipedia.org/wiki/Bowyer%E2%80%93Watson_algorithm
2.https://blog.csdn.net/electech6/article/details/116233056