双目Stereo重建算法SGM（1） - 互信息（Mutual Information）

1 概述

对于机器视觉，我目前还是刚起步，可能步还没起来。

近日开始深入学习双目Stereo相关算法，并主要参考了Heiko Hirschmüller于2008年发表的文献(Hirschmüller 2008)，其中提出了相关算法的主要思想与方法。该文献的主要贡献在于利用互信息（Mutual Information）作为计算matching cost的基础进行立体重建工作，在计算机视觉领域产生了深远的影响。值得注意的是，他提出的方法不仅具有创新性，在后续研究中也得到了广泛应用与改进。具体而言，在当前主流的OpenCV框架中（当前我使用的是版本3.4.1），其核心思想已经被集成到SGBM实现中作为基础模块运用。这些日子致力于复现实验中的关键算法，并希望以此为契机进行进一步的研究探索与技术积累。然而由于个人基础较为薄弱，在理解和复现过程中仍存在诸多疑问与不足之处；因此选择将此作为一种研究过程中的记录方式，并便于后续回顾与同行交流讨论

为了复现算法的过程，在Ubuntu系统上搭建开发环境并进行代码调试研究中发现了一些问题

本研究将详细阐述算法实现的具体步骤

2 基本假设

为了便于分析，在此假设两台构成双目相机的摄像头是完全相同的，并且已经经过立体校准处理。经rectification处理后的灰度图像在OpenCV中具有数据类型CV_8UC1。

3 Warp

为了计算matching cost, 我们需要利用互信息. 根据(Hirschmüller 2008)在文中的描述, 计算双目成像中两个图片的互信息时, 需要针对当前的disparity map, 对其中一个图片进行变形处理. 这里的"变形"也是首次 encounter. 我们通常将双目系统中左相机获取到的图象定义为基准图象（保持不变）, 右相机获取到的图象作为匹配图象, 而disparity map则表示右图象向左图象进行匹配时, 右图象中每个像素与左图象中最佳匹配像素在水平(x)坐标上的差值. 在计算互信息的过程中, 我们需要对右图象进行变形处理. 具体来说, 就是对左图象中的每一个像素坐标(x₁i, y₁i), 根据当前的disparity map确定在右图中对应的像素坐标(x₂i, y₁i), 然后从右图中提取该像素并将其放置到新生成的变形后图象（warped image）相应的位置(x₁i, y₁i). 这样就得到了一个变形后的warped image. 在后续计算互信息时, 我们将基于左图片和这个变形后的warped image来进行比较和分析, 而不再直接使用原始右图片.

确实存在一个潜在的问题。原文没有明确说明warp的具体实现方法,尤其是在左右图像匹配过程中无法获取到对应点的像素区域的处理方式,该情况同样没有提供相应的解决方案。由于左右相机在拍摄同一场景时必定存在一个重叠区域,但该重叠区域一般不可能覆盖整个Field of View（FOV）。因此,对于处于非重叠区域中的像素点,它们也无法获得对应的disparity值,而优化算法在这种情况下可能会给出不可靠的结果并认为这些像素不具备有效的disparity信息。另外,对于无法获取到disparity值的像素区域,同样没有提供有效的处理办法

我目前采用了Naïve的方法：对于右图中拥有相同的视差值的像素点，则选择不进行 warping 处理并保留其在原始位置上。这种方法未必完全准确，并且现阶段还不能确保所有像素点都能实现良好的覆盖效果。后续仍需进一步优化和完善

4 Joint Entropy and Mutual Information

该研究指出，在计算mutual information的过程中，其主要步骤是对joint entropy进行估计。这一过程可被表示为引用公式（4）、（5）及（6）。

[

该算法通过基于相似性计算的方法来实现目标图像的位置估计。

[

h_{I_1,I_2}(i,k)= -\frac{1}{n} \cdot \log\left(P_{I_1,I_2}(i,k)\otimes g(i,k)\right)\otimes g(i,k)\%0

[

该算法通过计算每个位置的匹配度...

其中， H I1 ,I2 即为joint entroy。这里我花了蛮久时间理解它们。

在本研究中采用了清晰的符号语言。其中_I₁_ 和_I₂_分别代表基于强度的左右两张图像。而_I₁^p_ 和_I₂^p_则代表左右两张图像中对应像素点处的强度。在这里,I₂
代表经过变形校正后的图像。由此,H(I₁,I₂)形成一个二维表格。这个表格具有固定大小为
256×256。无论_I₁
和_I₂
自身的大小如何变化,_H(I₁,I₂)_始终维持这一固定大小。

在原始公式（5）中,对该变量P_{I1,I2}执行了两次卷积操作.所有卷积核均为高斯内核.鉴于本人在概率统计理论基础方面的薄弱,耗费大量时间试图弄懂公式（5）的核心意图.原始描述指出,第一次卷 conv相当于执行了一次 Parzen 估计.其中 Parzen 估计实际上等同于内核密度估计（KDE）.作为概率统计中的一个基本概念,KDE是一种无需预先设定特定的概率密度函数形式的非参数化概率密度估计方法.具体而言,"无参数化"意味着在对随机变量的概率密度进行估计时,完全依赖所获得的实际数据而非预先假设任何特定的概率密度函数形式.因此,在公式（5）中,第一次卷 conv相当于执行了一次内核密度估计.

Kernel density estimation: https://en.wikipedia.org/wiki/Kernel_density_estimation

核函数（统计学）：参考文献链接

Would you mind explaining Parzen window (kernel) density estimation in plain language? https://stats.stackexchange.com/questions/244012/can-you-explain-parzen-window-kernel-density-estimation-in-laymans-terms

Interpretation of the bandwidth value within the framework of kernel density estimation.

KDE与卷积运算的关系可参考

http://ieor.berkeley.edu/~ieor165/ 页面上的第十节内容为《概率分布表征》

原文式(5)描述的第二个高斯卷积操作仍然未能让我完全掌握其具体目的。根据原文的解释，该过程将KDE计算出的结果进一步转化为lookup table形式，然而我仍然未能完全理解其深层含义。

在这个式(5)中存在两个不确定性问题：其一是_n_的具体取值范围是什么？我认为_n_指的是图像中的像素总数（number of all correspondences）。其二是式(5)中出现了对数计算（logarithmic computation），这一操作不允许自变量等于零（zero）。然而，在实际应用中可能存在自变量等于零的情况（包括必然出现的情况），因此建议将这些情况下的零替换为一个极小正数（small positive number）来避免违反对数函数定义域的要求。那么具体采用多小的一个极小正数值？后面采用1e-30作为替代值进行处理。这个替代值的选择是基于个人研究习惯的结果（arbitrary choice）。在互信息计算过程中使用单精度浮点型变量来进行数据运算（single precision floating-point variable）[citation:Single-precision floating-point format]。需要注意的是，在旧版本中的OpenCV库中对于Mat对象的操作与新版本有所不同。（In older versions of OpenCV, the logarithmic operation on Mat objects differs from newer versions.）

v 2.4版本操作内容包括对输入数组源（InputArray src）执行对数运算（log），并将其结果存储到输出数组 destination（OutputArray dst）。

v3.4:

该代码库中的核心数组实现细节可参考此标签：https://docs.opencv.org/3.4/d2/de8/group__core__array.html#ga937ecdce46...

总体而言，在3.4版本中对于自变量出现负值以及零的情况都更加符合对log()函数预期。该版本中log()函数对自变量为负值以及零的同时还包含非数字（NaN）以及无限大（Inf）的处理均被设定为"未定义"状态, 因此必须迫使用户对这些情况进行特殊处理, 这一改进确实增强了程序的整体鲁棒性和调节灵活性

在探讨完前面两个问题后，接下来是如何实现卷积运算。实际上，卷积运算可被视为一种图像模糊处理方法（smoothing）。当采用Gaussian kernel作为卷积核时，则是对图像执行高斯模糊处理。该功能的设计可见于OpenCV源码中，并具体实现位于opencv2/core/base.hpp文件中的第268行。

我认为将采用BORDER_ISOLATED策略，并以避免引入超出图像范围的数据。但ROI的具体定义是什么呢？然而，在OpenCV关于GaussianBlur函数的官方文档中并未对此进行明确说明。目前看来我已经采用了这一策略

而关键的问题是原文并未说明进行卷积时，图像的边界如何处理。

进行Gaussian smoothing时，在OpenCV中通常会使用一个7x7大小的核来完成平滑处理；这个过程仅仅是一个参数设置的问题。

卷积运算部分讨论结束后，则是原文式(6)的具体计算过程。实际上最初我想到的是, 式(6)实际上是基于两个单通道灰度图完成了一个二维直方图的计算。不过这个直方图具有特殊性, 即其取值范围均为0至255, 并将每个维度划分为256个区间。这一描述在以下网页中也有相应的表述

Can someone explain how to determine the joint entropy measure for two images using OpenCV? The link here provides a detailed explanation: [https://stackoverflow.com/questions/14898349/how-can-i-calculate-the-joint-entropy-of-two-images-in-opencv]

于是式(6)的运算就变为在OpenCV中进行一次2D直方图运算。OpenCV提供了calcHist( )函数。这个函数的channels参数的定义非常迷。经过几次调试之后，可以确定当calcHist( )的输入图像为两个单通道的灰度图像时，channels定义为{0, 1}，这表示使用第一幅图像的第一个通道，和第二个图像的第一个通道。然后剩下一个问题是与OpenCV自身的设计有关的，即所生成直方图的数据type。OpenCV在这里的处理也比较隐晦，经过测试，使用两个CV_8UC1类型的灰度图像，生成的直方图的数据type为CV_32FC1，也就是单精度浮点数（float）。这是正确的方法，若生成直方图仍使用某些整数类型，可能导致类型的数值溢出。

以上便是对原文式(4)、(5)和(6)的讨论。

5 Entropy of a Single Image

为了最终获得mutual information值,仍需计算两个图像各自的熵.这些熵已被建议通过联合熵来推导.具体而言,即可通过对二维矩阵_P(I₁,I₂)_沿行或列求和,可得两幅图的概率分布.进而根据原文中的公式（7a）与（7b）计算各幅图的熵.我的观点是继续采用整个图像像素的数量作为计数依据.

OpenCV为此功能提供了一个直接执行的函数，在调用reduce()函数时，其中所使用的数据类型仍为CV_32FC1。

6 Mutual Information

基于之前的讨论以及对原文中式（8a）、（8b）、（9a）和（9b）的详细分析，在当前阶段我们已经能够进行初步的互信息计算了。然而，在当前阶段仍需假设存在一个disparity map。

基于现有的一个自定义立体摄像头系统，在一个静止的墙面上进行成像拍摄，并通过捕获到的图像作为样本来计算互信息值。

包含如图所示的左右摄像头所获取的画面

图表 1 中的左边图像（经过缩放调整后）

图 2

通过计算获得的结果显示为互信息图像（如图3所示），此处采用了恒定的disparity map，并使用了随机选择的数值值进行计算。为了方便地以灰度形式展示互信息值分布情况，在对每个位置的数据值进行了归一化处理之后得到了浮点型Mat对象。该对象经过归一化处理后能够被OpenCV imshow()函数正确地以灰度形式进行展示。如果希望将该数据保存为一张灰度图，则需要将该浮点型Mat对象与255相乘并进行适当缩放才能实现正确保存

第3幅图展示了 Mutual information

通过图形观察可知，在(https://paperpile.com/c/ShN2Fm/DvzD)文献中给出了一个类似的实例。值得注意的是，在本文中图形显示与默认图片坐标系存在区别：本文图形中的y轴方向与常规设置不同。此外，在处理不理想的disparity map时（如图3所示），mutual information的分布呈现出沿主对角线的扩散状态。而在理想情况下（如图4所示），mutual information在主对角线上应显示为明亮区域（白色），远离主对角线则较为暗淡（灰色）。为了进一步验证该方法的有效性，在实验中采用仅使用左图替代右图的方法（即两幅图像完全相同），所得mutual information灰度图像结果如图4所示。此时disparity map全为零值，在图像上确实显示出较为明显的主对角线亮带。

图 4

后期会针对原文的算法进行进一步研究和复现，现在先做到这里。

参考文献

Hirschmüller, Heiko, 在2008年发表于《IEEE Transactions on Pattern Analysis and Machine Intelligence》期刊上的一篇文章中提出了“立体处理通过半全局匹配与互信息方法”的研究