基于Stereo Processing by Semiglobal Matching and Mutual Information对SGM算法的理解

看大佬的分析，醍醐灌顶！大家对算法的一些困惑也贴上来。
(转载http://blog.csdn.net/wsj998689aa/article/details/49464017, 作者：迷雾forest)
SGM算法源于《Stereo Processing by Semi-Global Matching and Mutual Information》一文，我认为这篇文章是立体匹配算法中最给力的，放眼KITTI，可以发现目前排名前五十的算法几乎一半都是对SGM的改进，具有最强的实用价值。SGM中文名称“半全局匹配”，顾名思义，其介于局部算法和全局算法之间，所谓半全局指的是算法既没有只考虑像素的局部区域，也没有考虑所有的像素点。例如，BM计算某一点视差的时候，往往根据目标像素周围的矩形区域进行代价聚合计算；DoubleBP在计算目的像素视差的时候，考虑的则是图像所有的像素点。抛开具体的方法不说，SGM中考虑到的只有非遮挡点，这正是定义为半全局的原因。
有的童鞋这时可能会想到Non-Local，为啥它要叫“非局部”，难道也是一种半全局算法？我认为不是，NLCA的目的是构建代价聚合项，它没有能量函数构建优化过程，而构建代价聚合基本上是局部算法才要做的事情，然而基于的却是所有的像素点，所以被称为非局部，并且NLCA和SGM的套路也完全不同。

我总是喜欢在博客开头扯一点之前写过的算法，回归正题，SGM的文献比较晦涩，作者Heiko Hirschmuller进行了大篇幅的论述，其中涉及到了很多的细节，最开始读这篇文献的时候我完全不理解算法的细节，但是SGM又经典的不行，没办法只能在反复理解论文的同时，结合代码进行理解。这篇文章主要涉及三个部分：分层互信息的代价计算+基于动态规划的代价聚合+其他，我打算分别基于这三个部分写三篇博客，与童鞋们分享我对SGM的理解。
和所有的立体匹配算法一样，SGM的第一步也是代价计算，它采用了基于互信息的计算方法，互信息是一种熵。那我们就先说说熵，熵是用来表征随机变量的不确定性（可以理解为变量的信息量），不确定性越强那么熵的值越大（最大为1），**那么图像的熵其实就代表图像的信息量。互信息度量的是两个随机变量之间的相关性，相关性越大，那么互信息就越大。可以想想看，两幅图像如果匹配程度非常高，说明这两幅图像相关性大还是小？显然是大，知道一幅图像，另外一幅图像马上就知道了，相关性已经不能再大了！！！反之，如果两幅图像配准程度很低，那么两幅图像的互信息就会非常小。所以，立体匹配的目的当然就是互信息最大化。**这就是为什么使用互信息的原因。熵和互信息的定义分别如下所示。

其中，互信息的前两项是图像的熵，最后一项叫做联合熵，联合熵是熵的变种，其依赖的自然是联合概率分布（熵依赖的是概率分布），联合熵的公式如下所示：

这时，一个最自然的问题是，图像的概率分布P是什么意思？答案一句话，图像的灰度直方图。图像的灰度值是0_{255，每个灰度值对应的像素个数除以图像像素个数就是该灰度值对应的概率，单幅图像的概率密度是一维的，那么自然地，两幅图像的联合概率密度就是二维的，它的定义域取值就是（0，0）}（255，255），公式如下所示：

其中，（i，k）指的是像素灰度值对，I1和I2分别代表了左图和矫正之后的右图，p代表像素点，n代表着像素的总数，公式的含义是统计不同灰度值对的个数并归一化。因此，两幅图像对应的概率分布可以用一幅图表示，这幅图的大小一定是256*256，弄清楚这点很重要，我之前一直误以为上述公式中的（i，k）指的是像素坐标，其实上面的公式只想告诉我们一件事情：联合概率分布就是归一化之后的统计直方图。

其中，作者称h为数据项，注意看它的自变量是像素灰度值，所以可以事先建立一个查表，将每个灰度值对的h值都保存下来（一共有256*256），那么联合熵的计算速度就会很快，h的计算公式如下所示：

（概率P是对所有的corresponding points（假设有n个）构成灰度值对的统计，可以用联合直方图计算。如果两个图像匹配得越好，计算得H就会非常小，加负号后就是非常大，也就是互信息会越大。

等号右边前两个H分别是对I1和I2而言，可以由之前联合分布直方图快速求到。至于为什么加上他们的原因，结合论文中的阐述，我认为忽略遮挡情况（可以在选点阶段避免）某幅图像的像素点P的信息熵，代表了这个点的灰度值在图中信息量的多少，HI越大，说明越信息量越多，越有价值。）

其中，g(i,k)指的是高斯平滑，正如前面所说，P是一幅表征概率值的图像（分辨率恒为256*256），正是针对这幅图像进行平滑处理，但是作者却在这里称之为Parzen估计，我想来想去觉得是一个意思，Parzen估计是一种非参数估计方法，不需要事先假设数学模型的参数形式而直接估计单点的概率值，最常用的Parzen窗函数就是高斯核函数。作者提供了一个图示也可以说明这个问题：

其中，第一个坐标系就是联合概率分布图，由于事先要对匹配图像（一般是右图）基于视差图进行修正（warp），调整右图像素的位置，使得左右两图尽可能相同之后计算同一像素位置的联合概率，由此导致了诸如（10,10）、（25,25）、（125,125）这样的灰度值对的概率值更大一些。再对这样的图像进行77窗口大小的高斯平滑处理，高斯平滑的目的是消除噪音，因为这里基于的视差图是很粗糙的（只有非遮挡点才有视差值，遮挡点的视差都为0），根据这样的视差图对右图进行处理，难免一些位置所对应的两图灰度值不相等，体现在P图上就是噪音，这个时候当然要平滑。例如可能有（10,100）这样的灰度值对出现，这说明右图该点经过warp之后很可能是错误的，但是这种情况肯定是少数，对应的概率值也就是很小，体现在P图上就是一个噪音点，用高斯平滑处理，一下子就没了。其实这就是检测outlier，去除outlier的另外一种手段，只不过SGM是从概率分布的角度处理罢了，而其他的算法一般都是从空间分布来处理outlier。进一步，取负对数便得到了h值的图示，最后得到了图像的联合熵图。
其实，平滑归平滑，也只是基于256*256的图像进行平滑，*只需根据像素位置计算出权重即可，同样可以根据查表的方式实现，这和具体的双目图像大小仍旧无关，一个不小的优点哦！
作者为什么要用互信息来计算代价匹配值？根据文献给出的信息，作者是考虑到了互信息对光照具有鲁棒性，我们注意看下面的公式，这个计算方式很有意思，我们已经得到了互信息图，剩下要做的事情只是根据左图和右图挑选出来的像素点的灰度值对，在互信息图中直接查找就行（又是一个速度优势），注意：要取个负号，这点直觉上很好理解，互信息越大->相关性越大->两个点的匹配程度越高->代价计算值理应越小，童鞋们千万不要忘记WTA是怎么运作的啊！！！！
通过利用MI计算的某个像素点P的matching cost可以如下表示

其中q是满足对极约束的在Im图像中的点，d是通过初始视差图D的数据。想在2D的图像找到使得cost aggregation 最小的视差图是个NP-hard的问题，那么单独在某方向r上（比如从左到右，从右到左等）选择匹配点的路径使得cost aggregation 最小（即最终L最小），该问题就变成多项式时间内可以解决的。
大家有没有觉得作者在此处的思维很是特别？先根据一个视差图作为先验来确定互信息图，然后基于这样的互信息图去计算各个视差的代价，这在理论上显然会导致第一步得到的视差图非常近似于初始视差图！！！难怪作者花了大篇幅介绍HMI，但是却没有给出中间结果。好比是法官审判一个杀人罪犯，他的判罚依据竟然是另一个抢劫罪犯的判罚结果，而不是根据其他杀人罪犯的判罚结果。此处我觉得有点怪，拿出来说说，大家都想想看哈。
本来熵和联合熵的计算方法不同，但是由于遮挡点的问题，将二者都设置为相同，都采用泰勒近似的方式。具体来说，图像的熵是基于图像的概率分布来计算的，但由于遮挡点的存在，有些像素根本就不会有匹配点，如果也将这样的像素考虑在内恐怕不合适，别忘了我们的**目的是令能够匹配的点尽量匹配，不能匹配的点我们根本就没必要让它们匹配，**于是将目光放在了联合熵的定义上，联合熵是基于联合概率分布的，其基于的点可以保证都是匹配点，这样的概率公式如下所示：

进一步，基于这样的概率分布得到的图像熵公式如下所示：

最后，作者采用分层互信息（HMI）进行代价匹配计算，由于概率分布图和图像大小无关（一直都是256*256），所以可以采用分层计算的方式来加速（反正不同层都是256*256），每一层的计算对应一次迭代，别忘记首次迭代需要基于视差图对右图进行warp，没有视差图怎么办？文章里面说是随机生成就好，并且由于像素个数很多，随机生成的个别错误可以被容忍，warp之后的右图和左图的匹配程度还是会比较高，迭代次数也相对较低，这也是SGM的一大优点。
  作者还对HMI的时间复杂度进行了计算，由于单次迭代的时间复杂度是O(WHD)，每次下采样都会令三个参数同时减半，所以上次迭代将会是当前迭代速度的1/8，假设我们有4次迭代，那么相比较于BT算法，只比它慢了14%，注意，**算法首次迭代使用的是最小的视差图，并且乘以3的原因是随机生成的视差图十分不靠谱，需要反复迭代3次才能得到同样分辨率下的靠谱视差图，然后再参与后续高分辨率的计算。**

总结：HMI是SGM的重要一步骤，通过HMI计算得到的代价匹配值对光照鲁棒，一旦这步做好了，后续的代价聚合，迭代求精等等操作弄起来就轻松多了。

SGM的代价聚合，其实仔细看看，这并不是严格意义上的代价聚合，因为SGM是为了优化一个能量函数，这和一般的全局算法一样，如何利用优化算法求解复杂的能量函数才是重中之重，其能量函数如下所示：（匹配代价是有噪音的，有时候错误的深度对应的匹配代价反而会比正确深度对应的匹配代价小。因此需要进行优化，也就是代价聚合。）

其中，C(p, Dp)代表的就是基于互信息的代价计算项，后面两项指的是当前像素p和其邻域内所有像素q之间的约束，如果q和p的视差只差了1，那么惩罚P1，如果大于1，那么惩罚P2，这么做基本上是机器学习中的常用方法，即所谓的正则化约束。这里需要注意的是，P2要大于P1，（第一项是匹配代价，第二项表示如果该像素的视差与周围像素相差1，则要加一个值为P1的惩罚，注意这里是累加符号，表明周围有几个视差相差1的像素，则要加几个P1，第三项表示如果该像素的视差与周围像素相差大于1，则要加一个值为P2的惩罚。也是周围有多少个像素满足就要加几个P2.。为什么要有第二项和第三项呢？这两项往往也被称为平滑项，可以这样理解，一个像素点的视差值取决于他自身的cost以及周围像素的视差情况，如果完全相信自身，那就是第一项，但是自身有误差，所以要考虑周围的像素的视差，那么就是后两项，但是要想相信周围的视差必须要加一个惩罚项，表明来自周围的信息的可靠程度没有自身cost的可靠程度高，P1的惩罚主要是为了在弯曲或者倾斜的表面进行调整（深度变化小），P2的惩罚则主要是为了处理深度不连续的情况（深度变化大）。）。

1. 假如不考虑像素之间的制约关系，不假设领域内像素应该具有相同的视差值，那么最小化E(D)就是最小化每一个C，这样的视差图效果很差，因为图像总会收到光照，噪声等因素的影响，最小的代价对应的视差往往是“假的”，并且这样做全然不考虑相邻之间的像素关系，例如，一个桌面的视差明显应该相同，但是可能由于倾斜光照的影响，每个像素的最小代价往往会不同，所以看起来就会乱七八糟，东一块西一块。这就是加上约束的目的。

2. 添加两个正则化项一方面是为了保证视差图平滑，另一方面为了保持边缘（我理解文献这一段中的 changes 是2种不同的“差异”，一个是相邻像素点视差值的差异，另一个是相邻像素点灰度值的差异。independent of their size，指的是P2与相邻像素点的视差差异的大小无关；而文献中这句话下面3行给出的P2公式，是说P2的值与相邻像素点的灰度差异的大小成反比。 也就是说，对于视差差异较大情况的惩罚，要视情况分别对待：如果此处相邻像素点灰度变化不大，则此处是连续物体的可能性更大，那么视差差异较大就是不合理的，因此要给予较大的P2来惩罚；但如果此处相邻像素点灰度的变化也很大，则此处是物体边缘的可能性就较大，那么视差差异较大就是合理的，因此要用较小的P2予以保留。 不知道这么理解是否合理呢~另一分析保持边缘保持的是object boundaries，要在smooth的同时保持object boundaries，这样叫做discontinuity preserving）。惩罚的越大，说明越不想看到这种情况发生，具体来说，如果q和p之间的视差有所差异但又不大，那么就要付出代价，你不是想最小化能量函数么？那么二者都要小，如果没有第二项，那么求出来的视差图将会有明显的锯齿现象，如果只有第三项，那么求出来的视差图边缘部分将会得到保持，但由于没有对相差为1的相邻像素进行惩罚，物体内部很可能出现一个“斜面”。

3. 这事情还没完，本文中有对这两项的解释，原文内容如下所示：

这句话的隐含意思是，如果我们让P1<P2，那么会允许出现小的斜面，也会保持边缘，前面一句我理解，惩罚的力度不大，就会导致这种事情还会发生，这也正是作者想看到的，水至清而无鱼嘛，不过，后一句中的P2并不是常数项，是根据相邻像素的差距来决定的，括号里面的“与大小无关”看起来就更加矛盾了，不知道哪位可以给好好解释一下这句话？
（ 这句引用了Graph Cut 论文, 里面有这么一句话: For example,in regularization-based vision [6], E smooth makes f smooth everywhere. This leads to poor results at object boundaries. Energy functions that do not have this problem are called discontinuity-preserving. 好像是用一个constant的penalty, 会使得所有地方都smooth, 这使边缘的结果不好, 所以P2是一个适应灰度梯度的值, 使得梯度很大的边缘, smooth不那么强, 来保持边缘. 另外, 指向II-B 会不会就是需要3个箭头? 算当前位置的cost的时候, 需要知道当前位置及对应点的灰度值, 然后把2个灰度值带入MI矩阵里面?）

有了能量函数，下面要做的就是求解它了，这个时候问题来了，这个E对p是不可导的，这意味着我们常看到的梯度下降，牛顿高斯等等算法在这里都不适用，作者于是采用了动态规划来解决这一问题，动态规划相信大家都知道了，但是其真正的精髓却是深藏不露，我早在大三期间就接触到了动态规划算法，这么多年过去了，虽然时而会用到这个算法，但到现在仍旧不敢说自己彻底懂它。。。。

简单地说，p的代价想要最小，那么前提必须是邻域内的点q的代价最小，q想要代价最小，那么必须保证q的领域点m的代价最小，如此传递下去。

本文只说说作者是怎么利用动态规划来求解E，其实这个求解问题是NP完全问题，想在2D图像上直接利用动态规划求解是不可能的，只有沿着每一行或者每一列求解才能够满足多项式时间（又叫做扫描线优化），但是这里问题来了，如果我们只沿着每一行求解，那么行间的约束完全考虑不到，q是p的领域的点其实这个时候被弱化到了q是p的左侧点或者右侧点，这样的求优效果肯定很差。于是，大招来了！！我们索性不要只沿着横或者纵来进行优化，而是沿着一圈8个或者16个方向进行优化。

这是一幅神奇的图示，我一直没有弄明白它到底是什么意思，笨死了，直到有一天我终于领悟它的真谛（仰天长啸）。我们先来看看优化求解过程：

（P1是对物体不平整表面的适应，P2是对灰度不连续（可能是多个物体）的适应，这样可以起到平滑视差图的效果。最后一项减去r方向在前一个像素点的L，添加这一项并不会改变最终视差图结果，只是为了限制L（不断增长至极大值）的规模。p1,p2的选取可以根据具体的应用进行调整，一般p1是一个定值，p2可以随着图像梯度变化而变化，往往要保证p2>=p1）
每一个点的代价聚合值是“当前代价+min（路径相邻点的当前视差代价聚合值，路径相邻点的视差差值为1的代价聚合值 + P1，路径相邻点的视差插值大于1的最小代价聚合值 + P2）- 路径相邻点的视差插值大于1的最小代价聚合值 ”，听起来够绕口的，其实就好比最小代价的蔓延，当前代价聚合值由当前代价和路径上一点的加了惩罚的最小代价聚合值所决定（最后那一项纯粹是为了防止数字过大，这是常用手段）。

其实为什么分解为8个方向想想看也很正常，能量函数E中每个p的能量是“自身代价本身+周围像素q带来的惩罚”，周围像素足足有8个，想求它们和的最小化十分难，最朴素的想法就是“分而求之”，我们就规定一个方向r，这个方向上p的邻居q只有一个，那么沿着这一方向的p的代价聚合值就成为了上面公式的样子。进一步，将8个方向的代价聚合值都加起来，就形成了p的最终代价聚合值。然后用WTA搞一下得到的视差图可以得到一个较小的能量E，目的就达到了。

我们来想想SGM的优化过程和DoubleBP有什么区别。

1. 先看能量函数，DoubleBP是每个像素自身代价加上周围像素的一个二元势函数值。SGM呢？是自身的代价加上周围像素带来的惩罚。其实二者是一个意思。
2. 再看优化过程，DoubleBP靠的是置信度传播算法，最后WTA的目标是一个置信度向量，这个置信度向量其实和向量没关系，每个分量都是去当前视差d的代价+周围像素的消息，这一点和SGM简直是太像了。
3. 再说说二者的区别，消息的每个分量可以理解为q对p取每个视差的支持力度，而SGM索性直接求取最小的惩罚，这点比DoubleBP要直接许多，所以SGM很快，DoubleBP很慢。

我认为这块内容非常值得单独拉出来说说，以后有时间好好的写写。

最后，我们可以看看SGM的整体流程图，这么长的流程图！！这个没啥好进一步解释的，唯一想说的就是我认为II-A那里多画了两个箭头，指向II-B的箭头应该只有MI一个。

Ib是参考图像(base image)，Im是待匹配的图像（matching image），第一次迭代的初始视差图Dinit可以随机产生。我们首先选择在Ib中的n个点，使他们尽量避免在两幅图的遮挡情况，然后找到Im中的对应位置通过Dinit变换后的n个点，按照立体视觉视差的计算步骤——matching cost，cost aggregation，disparity computation完成一次迭代，再把本次计算得到的视差图作为下一次的初始视差图进入下一轮迭代。

惊讶有两点还需要解释：在C阶段，作者交换“主”“客”位置，可以得到两个视差图Db和Dm，然后用一致性检测保证一对一的映射，否则该点作为黑点（255）。

其次，为了缩短迭代的时间，作者采用了层次计算matching cost（MI）的方法。Dinit和Ib，Im可以通过三次分辨率为原图的1/16（scale factor）来迭代得到，之后依次改为1/8，1/4、1/2的分辨率则迭代结束，Run time 比用BT（openCV中使用此方法）计算matching cost只多出14%。流程图中Down scale 和 Up scale 保证了图像在输入输出的大小不变。前面的1.14.

1. 对互信息基于泰勒展开进行近似计算的具体过程是什么？

2. 联合熵中的两个卷积是根据泰勒展开得到的，还是故意加上的？

因为此部分内容涉及到公式推理，在《Visual Correspondence Using Energy Minimization and Mutual Information》一文中有比较详细的推理过程，感兴趣的朋友可以去仔细看看，本文中就不详细说明了，希望和大家继续探讨。
大家的问题集锦：
问：始视差图该怎么获取，如果说用别的方法获取了，这不就相当于用别的立体匹配做了一遍又用SGM再做一遍？ 假设，可以用其他任何方法粗略获取视差图D，那么这个D如果不精确，就会对联合概率分布存在很大影响啊？因为两边的遮挡点并没有完全去除。 以上是我的理解，哪里出问题了啊？
迷雾forest回复:初始视差图的获取不用纠结，随机生成就好，并且因为这个视差图是随机的，所以在1/16要反复迭代3次，至于遮挡点，遮挡点是根据左右一致检测来确定的，和初始视差图关系不大。

参考
https://blog.csdn.net/tangli_sea/article/details/88218634
https://blog.csdn.net/lcc921528642/article/details/42550819
https://blog.csdn.net/wsj998689aa/article/details/50488249