传统立体匹配相似性测度笔记（更新中）

概述

主要参考文章《Evaluation of Stereo Matching Costs on Images with Radiometric Differences》进行传统立体匹配领域匹配相似性测度的梳理。论文链接。

匹配测度

文章比较了15种测度以及三种立体匹配方法的优劣。
将匹配测度分为三种类型：
1. non-parametric costs
2. parametric costs
3. mutual information

parametric costs

1. absolute difference(AD) && SAD

隐含假设是，认为同名点的亮度一致。
对于左图的一个像素$\mathbf{p}$而言，右图对应的像素为$\mathbf{p}-\mathbf{d}$。
$$C_{AD}(\mathbf{p},\mathbf{d})=|I_L(\mathbf{p})-I_R(\mathbf{p}-\mathbf{d})|$$
通常，对于局部的匹配方法来说，AD通常表现为SAD的方式，即：
$$C_{SAD}(\mathbf{p},\mathbf{d})=\sum _{\mathbf{q}\in N_{\mathbf{p}}}|I_L(\mathbf{p})-I_R(\mathbf{p}-\mathbf{d})|$$

2. BT测度

此外，还测试了 BT(Birchfield and Tomasi) 测度，原文为：

… the sampling-insensitive absolute difference of Birchfield and Tomasi…

该测度计算了the extrema of linear interpolations of the corresponding pixels of intrest with their neighbors. 这个测度通常用于pixel-wise的全局方法。 实际上，BT代价也是像素灰度值差值的绝对值，和AD有相似性，只不过BT利用了亚像素的灰度信息。

$$\begin{array}{rl}{C}_{BT}(\mathbf{p},\mathbf{d})& =\min (A,B)\\ A& =\max \left(0,I_L(\mathbf{p})-I_R^{\max }(\mathbf{p}-\mathbf{d}),I_R^{\min }(\mathbf{p}-\mathbf{d})-I_L(\mathbf{p})\right)\\ B& =\max \left(0,I_R(\mathbf{p}-\mathbf{d})-I_L^{\max }(\mathbf{p}),I_L^{\min }(\mathbf{p})-I_R(\mathbf{p}-\mathbf{d})\right)\\ I^{\min }(\mathbf{p})& =\min \left(I^{-}(\mathbf{p}),I(\mathbf{p}),I^{+}(\mathbf{p})\right)\\ I^{\max }(\mathbf{p})& =\max \left(I^{-}(\mathbf{p}),I(\mathbf{p}),I^{+}(\mathbf{p})\right)\\ I^{-}(\mathbf{p})& =\left(I\left(\mathbf{p}-{\left[\begin{array}{ll}1& 0\end{array}\right]}^T\right)+I(\mathbf{p})\right)/2\\ I^{+}(\mathbf{p})& =\left(I\left(\mathbf{p}+{\left[\begin{array}{ll}1& 0\end{array}\right]}^T\right)+I(\mathbf{p})\right)/2\end{array}$$

在opencv中的源码为：

for( int d = minD; d < maxD; d++ )
 {
 		// 右图逆序取出 
      int v = prow2[width-x-1 + d];
      int v0 = buffer[width-x-1 + d];// 最小值
      int v1 = buffer[width-x-1 + d + width2];// 最大值
      int c0 = std::max(0, u - v1); c0 = std::max(c0, v0 - u);
      int c1 = std::max(0, v - u1); c1 = std::max(c1, u0 - v);
	// 不同于论文，opencv综合了sobel滤波结果的BT以及原始灰度值的BT
      cost[x*D + d] = (CostType)(cost[x*D+d] + (std::min(c0, c1) >> diff_scale));
  }

参考链接

3. 滤波型的测度(mean, log,bilateral)

上图为滤波型的测度作用之前的原始图像，以下在介绍测度公式后所给出的结果，均是对该图进行的操作结果。

首先是均值型的，直接减去邻居像素的均值然后再加个常数，这里是128，加常数的目的是使其不为负值。
$$I_{\mathrm{mean}}(\mathbf{p})=I(\mathbf{p})-\frac{1}{\left|N_{\mathbf{p}}\right|}\sum _{\mathbf{q}\in N_{\mathbf{p}}}I(\mathbf{q})+128$$
处理结果为：

其次是LoG，目的是smoothing for removing noise and removes an offeset in intensities. 通常用于local real-time methods.这里采用$5\times 5$大小的LoG核。
$$I_{LoG}=I\otimes K_{LoG},K_{LoG}(x,y)=-\frac{1}{\pi {\sigma }^4}\left(1-\frac{x^2+y^2}{2{\sigma }^2}\right)e^{-\frac{x^2+y^2}{2{\sigma }^2}}$$
处理结果为：

最后，考虑了background subtraction by bilateral filtering(BilSub)。每一个像素都减去了双边滤波后的值，通过这个方式来做background subtraction.使用的LoG的核是$15\times 15$大小的，空间域参数${\sigma }_s=3$，以及灰度域参数${\sigma }_r=20$。
$$\begin{array}{l}{I}_{\text{BilSub }}(\mathbf{p})=I(\mathbf{p})-\frac{{\sum }_{\mathbf{q}\in N_{\mathrm{p}}}I(\mathbf{p})e^s{e}^r}{{\sum }_{\mathbf{q}\in N_{\mathrm{p}}}{e}^s{e}^r}\\ s=-\frac{(\mathbf{q}-\mathbf{p}{)}^2}{2{\sigma }_s^2}r=-\frac{(I(\mathbf{q})-I(\mathbf{p}){)}^2}{2{\sigma }_r^2}\end{array}$$
处理结果为：

4. ZSAD

ZASD全称呼是zero-mean sum of absolute differences，相对于SAD来说，改变在于：在计算SAD之前，像素值先减去窗口内的均值，其实本质上还是SAD。注意，被减去的窗口内的均值其实对于窗口内的每一个像素来说都是一样的，而不是像median filter那样，每一个像素要减去的值都不一样。
对于BilSub以及ZSAD这些测度，文中的表述是"…there are further common costs for removing an offset in intensities",我对其中"an offset in intensities"的理解是相当于去除类似于背景值。
$$\begin{array}{rl}{C}_{ZSAD}(\mathbf{p},\mathbf{d})& =\sum _{\mathbf{q}\in N_{\mathbf{p}}}\left|I_L(\mathbf{q})-{\bar{I}}_L(\mathbf{p})-I_R(\mathbf{q}-\mathbf{d})+{\bar{I}}_R(\mathbf{p}-\mathbf{d})\right|\\ \bar{I}(\mathbf{p})& =\frac{1}{N_{\mathbf{p}}}\sum _{\mathbf{q}\in N_{\mathbf{p}}}I(\mathbf{q})\end{array}$$

5. 大名鼎鼎的NCC

据我所知，在许多实际应用中，传统的匹配算法里用的测度不是NCC就是Census。NCC的全名是Normallized cross-correlation，文章中对其的评价是：“NCC compensates gain changes and is statistically the optimal method for dealing with Gaussian noise.”。在统计上来说，是最适应于高斯噪声的。公式为：
$$C_{NCC}(\mathbf{p},\mathbf{d})=\frac{{\sum }_{\mathbf{q}\in N_{\mathrm{p}}}{I}_L(\mathbf{q})I_R(\mathbf{q}-\mathbf{d})}{\sqrt{{\sum }_{\mathbf{q}\in N_{\mathrm{p}}}{I}_L(\mathbf{q}{)}^2{\sum }_{\mathbf{q}\in N_{\mathrm{p}}}{I}_R(\mathbf{q}-\mathbf{d}{)}^2}}$$

此外，在这篇博客中，对NCC也有相对详实的介绍，可供参考。

5.1 常见变种1—MNCC

这是一个基于Moravec算子的NCC变种。很接近于NCC，但是计算的速度更快，但是文中表示，MNCC似乎会对实验结果产生一定的影响。

5.2 常见变种2—ZNCC

类似于ZSAD之于SAD，ZNCC也是相当于NCC加了个zero-mean。ZNCC是唯一一个在gain和offset上都能进行补偿的parametric cost。公式为：
$$\begin{array}{rl}& C_{\text{ZNCC }}(\mathbf{p},\mathbf{d})=\\ & \frac{{\sum }_{\mathbf{q}\in N_{\mathrm{p}}}\left(I_L(\mathbf{q})-{\bar{I}}_L(\mathbf{p})\right)\left(I_R(\mathbf{q}-\mathbf{d})-{\bar{I}}_R(\mathbf{p}-\mathbf{d})\right)}{\sqrt{{\sum }_{\mathbf{q}\in N_{\mathbf{p}}}{\left(I_L(\mathbf{q})-{\bar{I}}_L(\mathbf{p})\right)}^2{\sum }_{\mathbf{q}\in N_{\mathbf{p}}}{\left(I_R(\mathbf{q}-\mathbf{d})-{\bar{I}}_R(\mathbf{p}-\mathbf{d})\right)}^2}}\end{array}$$

Non-Parametric Matching Costs

所谓的Non-parametric匹配测度依赖的都是影像灰度的相对顺序，即"… based on the local order of intensities".

1. Rank filter

Rank filter最先被提出时，是为了提高窗口匹配的鲁棒性，换言之，是为了起到对窗口内的粗差不敏感的作用。然而，再仔细想想，如果对窗口内的粗差不敏感了，那么，那些边缘处，本就是看似“粗差”的真实情况，是不是也就不能够被有效识别出来呢？确实是这样的，对于这种所谓的鲁棒的方法来说，他们也会进一步地导致物体边缘的blur情况。在论文中，使用了的Rank fillter的窗口大小是$15\times 15$大小的，具体公式为：

$$I_{Rank}(\mathbf{p})=\sum _{\mathbf{q}\in N_{\mathrm{p}}}T[I(\mathbf{q})<I(\mathbf{p})]$$

其中，$T[]$表示，如果是真则返回1，如果是假则返回0.

Rank filter在没有纹理的区域容易受到噪声影响，比如上图中小熊的右边墙壁。

1.1 改进 — Soft Rank filter

刚才提到，在rank filter中存在受弱纹理影响的问题，针对该问题，Zitnick等人提出了Soft rank filter。Soft rank filter的思想是，通过定义一个线性的，soft的过渡区（在0-1之间）来避免弱纹理的影响，而实验结果证明，这个想法是行之有效的。究其原因，弱纹理的影响本质上是因为像素的数值太过于相近的时候仍然表现为一大一小，而这是不符合实际情况的。具体的公式为：

$$I_{\text{SoftRank }}(\mathbf{p})=\sum _{\mathbf{q}\in N_{\mathrm{p}}}\min \left(1,\max \left(0,\frac{I(\mathbf{p})-I(\mathbf{q})}{2t}+\frac{1}{2}\right)\right)$$
在论文实验中，阈值$t=8$
结果如下图所示，可以看到，红框框出来的地方，明显保真多了。

2. 大名鼎鼎的census filter

census变换就是窗口中心的元素和邻居们一个一个的比较，然后表示为一个位串，两幅影像都经过这种变换后，对于每个像素来说，都有一个位串，两个影像之间就通过hamming距离来计算相似度。是不是很像rank filter？但是rank filter是变换了之后直接AD了，而census是计算汉明距离，而且更重要的是，rank不管空间邻居的结构相对位置信息是怎么样的，census管。这也就导致了census的效果，要比rank好。当然，速度也要比rank慢一些。论文中的实验使用的窗口大小是$9\times 7$的。

3. ordinal measure

这个方法原文是这样讲的：“…based on the distance of rank permutations of corresponding matching windows.”
这个方法不能够用filter来实现，而且需要基于窗口的匹配。其相对于Rank还有Census的优势可能就是他不依赖于兴趣点的数值。

Mutual Information

互信息。这个测度其实也是作者自己提出来的，对了，这个作者就是提出SGM的那位大佬。
类似于$9\times 9$,$11\times 11$这样的窗口，体现不出meaningful probability distribution.这就需要窗口大一些，然而如果窗口过大的话，就会导致边缘的blur。因此，采取了一种能够利用逐像素的匹配方法的MI的计算方式（其实就是SGM的方法）。首先用一个随机的视差图，通过这个视差图，我们可以知道两张影像同名点的关系，进而，我们就可以计算需要的probability distributions. 因为这个计算是对于整张影像来说的，所以这个得到的probability distributions是可靠的。通过MI的泰勒扩展，使得代价矩阵是可微的。初始的视差图可以随便给，之后迭代精化就可以了，文中表示，一般来说迭代个三次就够了。为了效率考虑，文章使用HMI计算，即Hierarchical MI，意思就是先从低分辨率开始算，然后逐层传递回高分辨率去。其中有一个细节值得注意，就是"It should be noted that the
disparity image of the lower-resolution level is used only for
calculating the matching costs of the higher-resolution level,
but not for restricting the disparity range, as this could easily lead to missing small objects"。这种分层的计算方式和迭代的效果几乎一样好，但是用的时间却更少了，奈斯。

其他

AD-Census

前面介绍了AD，介绍了Census，那么AD-Census顾名思义就是AD和Census两种测度的结合。显然，AD和Census的尺度并不一样，当我们想要结合两个还不错的东西的时候，或许首先要想到的就是这两个东西在结合的时候是否是同权重的？而这，就自然而然的会涉及到尺度，试想，如果一个测度的范围在10-100000之间，另一个只在0-1之间，那么如果简单的相加的话，自然会削弱那个0-1范围的测度所起到的作用。因此，测度尺度的范围一定要解决，在Ad-Census中，是通过引入一个自然指数函数来解决的：

$$\rho (c,\lambda )=1-\exp \left(-\frac{c}{\lambda }\right)$$
其中，$c$是代价值，而$\lambda$是一个控制参数（大于0），显然，通过这个自然指数参数，可以将任何的代价$c$归一化到0-1之间。
将这个自然指数参数应用至AD代价以及Census代价，并将作用后的代价相加，即为Ad-Census的计算公式了：

$$\begin{array}{rl}C(\mathbf{p},d)=& \rho \left(C_{\text{census }}(\mathbf{p},d),{\lambda }_{\text{census }}\right)+\\ & \rho \left(C_{AD}(\mathbf{p},d),{\lambda }_{AD}\right)\end{array}$$

这个范围在0-2之间。
参考李博博客中的代码为：

	// ad代价
	const float32 cost_ad = abs(img_left_[y*width+x] - img_right_[y*width+x_r]);

	// census代价
	const auto& census_val_r = census_right_[y * width_ + xr];
	const float32 cost_census = static_cast<float32>(Hamming64(census_val_l, census_val_r));

	// ad-census代价
	cost = 1 - exp(-cost_ad / lambda_ad) + 1 - exp(-cost_census / lambda_census);

AD-Census即综合了AD在重复纹理区域的优势，Census在弱纹理区域的优势。
理论参考链接

AD-Census出自文章《On building an accurate stereo matching system on graphics hardware》

AD + gradient

实际上就是sad损失与梯度损失的加权配比。