立体匹配-----NCC视差匹配

目录

一、立体匹配算法

　　　1.立体匹配算法分类

二、NCC 视差匹配方法

           1.原理

 　　　2.NCC计算公式

　　　3.算法流程

　　　4.代码实现

　　    5.不同场景运行

三、结论

四、遇到的问题及解决方法

一、立体匹配算法

1.立体匹配算法分类

 在立体匹配中，匹配问题可以看成是寻找两组数据相关程度的过程。根据采用图像表示的基元不同，立体匹配算法有多种分类。

  ①根据算法运行时约束的作用范围：分为局部匹配算法和全局匹配算法。

  ②基于生成的视差图：可分为稠密匹配和稀疏匹配。稠密匹配：是基于生成的视差图，对于所有像素都能生成确定视差值，称为稠密匹配。稀疏匹配：只选择关键像素点[通常为角点或者边缘点]计算视差值的方法称为稀疏匹配，该算法计算速度较快，但后续还需要通过插值算法计算缺失像素点的视差值，因此应用场景上有很大限制。

1.1全局匹配算法

  全局（半全局）立体匹配算法主要是采用了全局的优化理论方法估计视差，建立一个全局能量函数，其包含一个数据项和平滑项，通过最小化全局能量函数得到最优的视差值。其中，图割、置信传播、动态规划,粒子群算法、遗传算法等优化算法都是常用的求解能量最小化的方法。（以后作补充，也可参考我之前的粒子群算法、遗传算法文章）

1.2局部匹配算法

（1）算法介绍：局部立体匹配算法又称基于窗口的方法或基于支持区域的方法。算法对参考图像中的每个像素计算一个合适大小、形状和权重的窗口，然后对这个窗口内的视差值进行加权平均。理想的支持窗口应该完全覆盖弱纹理区域，并在窗口内深度连续。与全局立体匹配算法相似，通过优化一个代价函数的方法计算最佳视差。但是，在局部立体匹配算法的能量函数中，只有基于局部区域的约束数据项，没有平滑项。局部匹配算法仅利用某一点邻域的灰度、颜色、梯度等信息进行计算匹配代价，计算复杂度较低，大多实时的立体匹配算法都属于局部立体匹配的范畴，但局部立体匹配算法对低纹理区域、重复纹理区域、视差不连续和遮挡区域匹配效果不理想。

（2）基本原理：在参考图像中选择一个点，选择该点邻域内一个支持窗口，然后依据一定的相似性判断准则，在待匹配图像中寻找与支持窗口最相似的子窗口，该子窗口所对应的像素点即为对应的匹配点。

（3）不同场景运用：固定窗口代价聚合使用固定大小和形状的窗口作为代价聚合的基元，通常是一个矩形，并假设支持窗口内的其它像素点与待匹配点具有相同的视差。固定窗口法精度不高，但易实现、耗时短，在一些对实时性要求极高的场合得到了应用。

  基于双边滤波的代价聚合算法仍然使用固定大小和形状的窗口，但窗口内的元素权重不同，权重由目标图像在该窗口内像素与窗口中心的灰度差和距离计算。基于双边滤波的代价聚合算法精度高，但计算复杂，实时性差，算法性能随窗口尺寸指数增加。

  基于分割的代价聚合算法的主要思想是：预先将作为参考图像的左图进行分割，对于支持窗口内与窗口中心处于同一分割的像素，对应的权值取1，否则为一个远小于1的正数。但是图像分割是一个非常耗时的操作，同样无法在实时性要求较高的场合使用。

  基于十字的代价聚合算法（Cross-based Cost Aggregation,CBCA）的支持窗口形状并不确定，会根据匹配点邻域的灰度值而改变，具体实现后面的更新中会介绍。该方法可以使用GPU并行计算，具有较好的实时性，现广泛应用于各种算法的代价聚合步骤。

（4）优点：算法成熟、计算简单、速度快，能进行图像实时处理，匹配精度较高，计算量比全局匹配算法更小。

（5）缺点：对噪声敏感，弱纹理区易造成误匹配

二、NCC 视差匹配方法

1.原理：对于原始的图像内任意一个像素 点(px, py )构建一个n x n的邻域作为匹配窗口。然后对于目标相素位置(px + d, py )同样构建一个n x n大小的匹配窗口，对两个窗口进行相似度度量,注意这里的有一个取值范围。对于两幅图像来说，在进行NCC计算之前要对图像处理，也就是将两帧图像校正到水平位置，即光心处于同一 水平线上，此时极线是水平的,否则匹配过程只能在倾斜的极线方向上完成,这将消耗更多的计算资源。

2.NCC计算公式：

 

 

 

求出的NCC(p,d)值在[−1,1]之间

3.算法流程：

 ①采集图像：通过标定好的双目相机采集图像，当然也可以用两个单目相机来组合成双目相机。 

 ②极线校正：校正的目的是使两帧图像极线处于水平方向，或者说是使两帧图像的光心处于同一水平线上。通过校正极线可以方便后续的NCC操作。 

 ③特征匹配：这里便是我们利用NCC做匹配的步骤，匹配方法如上所述，右视图中与左视图待测像素同一水平线上相关性最高的即为最优匹配。完成匹配后，需要记录其视差d，即待测像素水平方向xl与匹配像素水平方向xr之间的差值d=xr−xl，最终可以得到一个与原始图像尺寸相同的视差图D。

 ④深度恢复：通过上述匹配结果得到的视差图D，我们可以很简单的利用相似三角形反推出以左视图为参考系的深度图。计算原理如下图所示：

 

 

 

如图，Tx为双目相机基线，f为相机焦距，这些可以通过相机标定步骤得到。而xr - xl就是视差D。

通过公式 z = f * Tx / d可以很简单地得到以左视图为参考系的深度图了。

4.代码实现：

（1）源代码

# -*- coding: utf-8 -*-
from PIL import Image
from pylab import *
import cv2
from numpy import *
from numpy.ma import array
from scipy.ndimage import filters
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(np.roll(norm_l, -displ - start) * norm_r, wid, s) # 和归一化
        filters.uniform_filter(np.roll(norm_l, -displ - start) * np.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 np.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(np.roll(norm_l, -displ - start) * norm_r, wid, 0, s) # 和归一化
     filters.gaussian_filter(np.roll(norm_l, -displ - start) * np.roll(norm_l, -displ - start), wid, 0, s_l)
     filters.gaussian_filter(norm_r*norm_r,wid,0,s_r) # 和反归一化
     # 保存 ncc 的分数
     dmaps[:,:,displ] = s / np.sqrt(s_l * s_r)
 # 为每个像素选取最佳深度
 return np.argmax(dmaps, axis=2)

im_l = array(Image.open(r'C:/Users/LE/PycharmProjects/untitled/NCC/3.jpg').convert('L'), 'f')
im_r = array(Image.open(r'C:/Users/LE/PycharmProjects/untitled/NCC/4.jpg').convert('L'),'f')
# 开始偏移，并设置步长
steps = 12
start = 4
# ncc 的宽度
wid = 9
res = plane_sweep_ncc(im_l,im_r,start,steps,wid)
import scipy.misc
scipy.misc.imsave('depth.png',res)
show()