基于python3计算机视觉编程(七)NCC实现视差图匹配

一. 立体视觉的研究背景

立体视觉是计算机视觉领域的一个重要课题，它的目的在于重构场景的三维几何信息。立体视觉的研究具有重要的应用价值，其应用包括移动机器人的自主导航系统，航空及遥感测量，工业自动化系统等。

立体视觉的研究方法之一，利用多幅图象来恢复三维信息的方法，它是被动方式的。根据图象获取方式的区别又可以划分成普通立体视觉和通常所称的光流(optical
flow)两大类。普通立体视觉研究的是由两摄像机同时拍摄下的两幅图象，而光流法中研究的是单个摄像机沿任一轨道运动时顺序拍下的两幅或更多幅图象。前者可以看作后者的一个特例，它们具有相同的几何构形，研究方法具有共同点。双目立体视觉是它的一个特例。

其中立体匹配是立体视觉系统的核心，是建立图像间的对应从而计算视差的过程，是极为重要的。

二. 双目立体匹配NCC

2.1 NCC算法简介

归一化相关性，normalization cross-correlation，因此简称$NCC$。顾名思义，就是用于归一化待匹配目标之间的相关程度，注意这里比较的是原始像素。通过在待匹配像素位置$p$（$p_x$，$p_y$）构建$3\ast 3$邻域匹配窗口，与目标像素位置$p^{\prime }$($px+d$,$py$）同样构建邻域匹配窗口的方式建立目标函数来对匹配窗口进行度量相关性，注意这里构建相关窗口的前提是两帧图像之间已经校正到水平位置，即光心处于同一水平线上，此时极线是水平的，否则匹配过程只能在倾斜的极线方向上完成，这将消耗更多的计算资源。

2.1 NCC基本原理

NCC是一种基于统计学计算两组样本数据相关性的算法，其取值范围为[-1, 1]之间，而对图像来说，每个像素点都可以看出是RGB数值，这样整幅图像就可以看成是一个样本数据的集合，如果它有一个子集与另外一个样本数据相互匹配则它的ncc值为1，表示相关性很高，如果是-1则表示完全不相关，基于这个原理，实现图像基于模板匹配识别算法，其中第一步就是要归一化数据。$NCC$计算公式如下图所示：

• 其中$NCC(p,d)$$NCC(p,d)NCC(p,d)$得到的值得范围将在$[-1,1]$之间$[-1,1]$之间$[-1,1]$之间。
• $W_p$为之前提到的匹配窗口。
• $I_1(x,y)$为原始图像的像素值。
• $I_1(px,py)$为原始窗口内像素的均值。
• $I_2(x+d,y)$为原始图像在目标图像上对应点位置在xxx方向上偏移ddd后的像素值。
• $I_2(p_x+d,p_y)$为目标图像匹配窗口像素均值。
• 若$NCC=-1$，则表示两个匹配窗口完全不相关，相反，若$NCC=1$时，表示两个匹配窗口相关程度非常高。

2.2 匹配流程

1. 采集图像：通过标定好的双目相机采集图像，当然也可以用两个单目相机来组合成双目相机。
2. 极线校正：校正的目的是使两帧图像极线处于水平方向，或者说是使两帧图像的光心处于同一水平线上。通过校正极线可以方便后续的$NCC$操作。
   

• 由标定得到的内参中畸变信息中可以对图像去除畸变。
• 通过校正函数校正以后得到相机的矫正变换$R$和新的投影矩阵$P$，接下来是要对左右视图进行去畸变，并得到重映射矩阵。

3. 特征匹配：这里便是我们利用$NCC$做匹配的步骤啦，匹配方法如上所述，右视图中与左视图待测像素同一水平线上相关性最高的即为最优匹配。完成匹配后，我们需要记录其视差$d$，即待测像素水平方向$x_l$与匹配像素水平方向$x_r$之间的差值$d=x_r-x_l$ ，最终我们可以得到一个与原始图像尺寸相同的视差图$D$。
4. 深度恢复：通过上述匹配结果得到的视差图$D$，我们可以很简单的利用相似三角形反推出以左视图为参考系的深度图。计算原理如下图所示：
   

如图，$Tx$为双目相机基线，$f$为相机焦距，这些可以通过相机标定步骤得到。而$x_r-x_l$就是视差$d$。通过公式$z=\frac{f\ast T_x}{d}$可以很简单地得到以左视图为参考系的深度图了。

至此，我们便完成了双目立体匹配。倘若只是用于图像识别，那么到步骤3时已经可以结束了。
参考博文https://blog.csdn.net/qq_38898129/article/details/93074387?depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromBaidu-5&utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromBaidu-5

三. 视差算法立体匹配

3.1 视差算法简介

立体匹配算法也称为视差算法,目的是求解两幅图片之间的视差值。主要是通过找出每对图像间的对应关系，根据三角测量原理，得到视差图；在获得了视差信息后，根据投影模型很容易地可以得到原始图像的深度信息和三维信息。

影响视差算法的因素：

• 光学失真和噪声（亮度、色调、饱和度等失衡）
• 平滑表面的镜面反射
• 投影缩减（Foreshortening）
• 透视失真（Perspective distortions）
• 低纹理（Low texture）
• 重复纹理（Repetitive/ambiguous patterns）
• 重叠和非连续
• 透明物体

3.2 算法步骤

1. 匹配代价计算： `

匹配代价计算的目的是衡量待匹配像素与候选像素之间的相关性。两个像素无论是否为同名点，都可以通过匹配代价函数计算匹配代价，代价越小则说明相关性越大，是同名点的概率也越大。匹配代价的算法由很多种，本文采用$NCC$算法匹配。

`
2. 代价聚合：

代价聚合的根本目的是让代价值能够准确的反映像素之间的相关性。上一步匹配代价的计算往往只会考虑局部信息，通过两个像素邻域内一定大小的窗口内的像素信息来计算代价值，这很容易受到影像噪声的影响，而且当影像处于弱纹理或重复纹理区域，这个代价值极有可能无法准确的反映像素之间的相关性，直接表现就是真实同名点的代价值非最小。

3.视差获取：

对于区域算法来说，在完成匹配代价的叠加以后，视差的获取就很容易了，只需在一定范围内选取叠加匹配代价最优的点（$SAD$和$SSD$取最小值，$NCC$取最大值）作为对应匹配点，如胜者为王算法$WTA$（Winner-take-all）。而全局算法则直接对原始匹配代价进行处理，一般会先给出一个能量评价函数，然后通过不同的优化算法来求得能量的最小值，同时每个点的视差值也就计算出来了。

4. 视差细化（亚像素级）：

大多数立体匹配算法计算出来的视差都是一些离散的特定整数值，可满足一般应用的精度要求。但在一些精度要求比较高的场合，如精确的三维重构中，就需要在初始视差获取后采用一些措施对视差进行细化，如匹配代价的曲线拟合、图像滤波、图像分割等。

四. NCC视差图匹配代码实现

代码：

import stereo
import scipy.misc
from PIL import Image
from pylab import *
from scipy.ndimage import *


def plane_sweep_ncc(im_l, im_r, start, steps, wid):
    # 使用归一化的互相关计算视差图像 """
    m, n = im_l.shape
    # 保存不同求和值的数组
    mean_l = zeros((m, n))
    mean_r = zeros((m, n))
    s = zeros((m, n))
    s_l = zeros((m, n))
    s_r = zeros((m, n))
    # 保存深度平面的数组
    dmaps = zeros((m, n, steps))
    # 计算图像块的平均值
    filters.uniform_filter(im_l, wid, mean_l)
    filters.uniform_filter(im_r, wid, mean_r)
    # 归一化图像
    norm_l = im_l - mean_l
    norm_r = im_r - mean_r
    # 尝试不同的视差
    for displ in range(steps):
        # 将左边图像移动到右边，计算加和
        filters.uniform_filter(roll(norm_l, -displ - start) * norm_r, wid, s)  # 和归一化
        filters.uniform_filter(roll(norm_l, -displ - start) * roll(norm_l, -displ - start), wid, s_l)
        filters.uniform_filter(norm_r * norm_r, wid, s_r)  # 和反归一化
        # 保存 ncc 的分数
        dmaps[:, :, displ] = s / sqrt(s_l * s_r)
    # 为每个像素选取最佳深度
    return argmax(dmaps, axis=2)


def plane_sweep_gauss(im_l, im_r, start, steps, wid):
    # 使用带有高斯加权周边的归一化互相关计算视差图像 """
    m, n = im_l.shape
    # 保存不同加和的数组
    mean_l = zeros((m, n))
    mean_r = zeros((m, n))
    s = zeros((m, n))
    s_l = zeros((m, n))
    s_r = zeros((m, n))
    # 保存深度平面的数组
    dmaps = zeros((m, n, steps))
    # 计算平均值
    filters.gaussian_filter(im_l, wid, 0, mean_l)
    filters.gaussian_filter(im_r, wid, 0, mean_r)
    # 归一化图像
    norm_l = im_l - mean_l
    norm_r = im_r - mean_r
    # 尝试不同的视差
    for displ in range(steps):
        # 将左边图像移动到右边，计算加和
        filters.gaussian_filter(roll(norm_l, -displ - start) * norm_r, wid, 0, s)  # 和归一化

        filters.gaussian_filter(roll(norm_l, -displ - start) * roll(norm_l, -displ - start), wid, 0, s_l)
        filters.gaussian_filter(norm_r * norm_r, wid, 0, s_r)  # 和反归一化
    # 保存 ncc 的分数
    dmaps[:, :, displ] = s / sqrt(s_l * s_r)
    # 为每个像素选取最佳深度
    return argmax(dmaps, axis=2)


im_l = array(Image.open('G:/QQdata/tsukuba/scene1.row3.col1.ppm').convert('L'), 'f')
im_r = array(Image.open('G:/QQdata/tsukuba/scene1.row3.col5.ppm').convert('L'), 'f')
# 开始偏移，并设置步长
steps = 50
start = 4

# ncc 的宽度
wid = 12

res = stereo.plane_sweep_ncc(im_l, im_r, start, steps, wid)

imsave('depth.png', res)

4.1 探究不同窗口值对结果的影响

测试数据集：

测试数据可以从这里下载http://vision.middlebury.edu/stereo/data/scenes2003/本文使用的测试数据是cones

窗口大小wid = 3

窗口大小wid = 5

这两组的对比可以看出当窗口值过低时低纹理区域容易出现误匹配，匹配精度较低，并且噪声也很多。随着窗口值的增大，其中错误的匹配的区域也随之慢慢被矫正。不过这只是猜想，接下来继续增大窗口值以用来验证正确性。

窗口大小wid = 8

窗口大小wid = 10

对比wid=8与wid=10

随着窗口值增大可以印证上面的判断纹理低的区域错误的匹配基本上得到矫正，但是在信息的细节上越来越模糊可以从上图勾出的红圈与篮圈中看出。

窗口大小wid = 20

从wid=20可以明显看出，图片中的噪声几乎被矫正完，但随之而然的图片细节鲁棒也越来越差。由此可以得出结论。wid值越小图片对噪声越敏感，在低纹理区会出现很多错误匹配，但细节鲁棒好，值越大噪声鲁棒好，错误匹配随之被纠正但细节随之会被忽略掉，导致图片的精度小。

4.2 结论

• 不同的窗口值会影响匹配的结果，值越小，图像的噪声点会很多，但细节信息保存的很好。值越大，通过上图的对比可以看出图像的边缘部分过于平滑丢失掉了细节信息，但抗噪声的能力比较强。
• 通过上图实现可以发现，面具表面的细节纹理都没有被算法所体现出来，这是上文介绍到的光学失真(光线亮度)所引起的，所以在选取图片时，图片是否光线充足，是否透明，是否重叠都会影响算法的匹配精确度。