1. 立体匹配

1.1 概述

立体匹配是立体视觉研究中的关键部分。其目标是在两个或多个视点中匹配相应像素点，计算视差。通过建立一个能量代价函数，对其最小化来估计像素点的视差，求得深度。

点P和Q，映射到左相机OR像面上的同一点p≡q，只要找到p和q在右相机OT像面上的对应点就可以通过三角计算估计深度。找到对应点的过程，即立体匹配。

为了找到对应点，需要增加约束，最常用的是极线约束。

P和Q映射到左相机QR像面上的同一点p≡q，直线pq上的点对应点一定位于右像面的直线p’q’上，p’q’即为直线pq的极线，这就是极线约束。

接下来就可以根据视差估计深度，然后通过Graph cuts算法给每一个像素点分配视差从而得到深度图，不再详细说明。

1.2 主要立体匹配算法分类

1）根据采用图像表示的基元不同，立体匹配算法分为：
A、区域立体匹配算法：可获取稠密视差图。缺点：受图像的仿射畸变和辐射畸变影响较大；像素点约束窗口的大小与形状选择比较困难，选择过大，在深度不连续处，视差图中会出现过度平滑现象；选择过小，对像素点的约束比较少，图像信息没有得到充分利用，容易产生误匹配。
B、基于特征的立体匹配算法：可获得稀疏的视差图，经差值估计可获得稠密视差图。可提取点、线、面等局部特征，也可提取多边形和图像结构等全局特征。缺点：特征提取易受遮挡、光线、重复纹理等影响较大；差值估计计算量大。
C、基于相位立体匹配算法：假定在图像对应点中，其频率范围内，其局部相位是相等的，在频率范围内进行视差估计。
2）依据采用最优化理论方法的不同，立体匹配算法可以分为：
A、局部的立体匹配算法
B、全局的立体匹配算法

还有立体匹配算法介绍的更详细内容：参考博客

1.3 立体匹配的基本步骤

立体匹配过程：

1. 匹配代价计算： 一般是通过计算左右两图对应像素3个通道的灰度值差来决定匹配代价的，常用的就是基于像素点匹配代价计算，一般有AD，SD，TAD什么的，基于区域的匹配代价计算一般有SAD，SSD， STAD之类的。匹配代价计算会生成一个disparity space image，也就是DSI。这个DSI是一个三维的空间，也就是每一个视差，得到一张代价图。假如视差范围是0～16，则会得到17幅代价图。视差搜索范围就是MiddleBurry网站上的stereopair值，也就是说在视差范围内（比如0-16）内搜索匹配代价，得到17张匹配代价图，然后找到匹配代价最小的对应的视差值就是此像素对应的视差。
2. 代价聚合：其实就是一个滤波的过程，对每一幅代价图进行聚合，最简单的就是采用boxfilter。第一步代价计算只是得到了图像上所有孤立像素的视差值，但是这些时差值都是孤立的，引入了过多噪声，比如一片区域的视差值都是10，可是引入噪声后就会导致这一片的视差值都不一样，那么就需要一个滤波的过程，也就是我们所说的局部立体匹配方法，即采用窗口卷积达到局部滤波的目的。
3. 计算视差：常用的方法就是WTA算法（局部），对于图像中的同一个点，选出17幅代价图中匹配代价最小的那张图，该幅图对应的视差值就选取为最终的视差。或者在全局立体匹配中采用能量函数的方法，分为数据项和平滑项，数据项其实就是代价计算，平滑项就是代价聚合，只不过窗口大小是整幅图像，也可以试试如果把平滑项前面的系数lamda设为0，那么得到的结果和单纯代价计算的局部立体匹配是一样的。
4. 视差精化：也就是对得到的视差进行优化的过程，如：左右一致性检测、区域投票等；这步其实是很多立体匹配的遮羞布，比如用遮挡处理，中值滤波，左右一致性检测等，都能使最后的是视差图提升1%左右，它是很多论文的遮羞布。但是不可否认的是，立体匹配最关键的步骤仍然是代价计算和代价聚合步骤。
   在立体匹配方法中，基于全局和局部的算法有些区别。不过基本步骤都差不多。有些时候，基于局部的算法，第一步和第二步是合并在一起进行的，基于全局的算法，会跳过第二步。

2. 归一化互相关（NCC）视差匹配法

2.1 原理

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

$NCC$计算公式如下：
$NCC(p,d)=\frac{{\sum }_{(x,y)\in W_p}(I_1(x,y)-\overline{I_1}(P_x,P_y))\cdot (I_1(x+d,y)-\overline{I_2}(P_x+d,P_y))}{{\sum }_{(x,y)\in W_p}(I_1(x,y)-\overline{I_1}(P_x,P_y){)}^2\cdot {\sum }_{(x,y)\in W_p}(I_2(x+d,y)-\overline{I_2}(P_x+d,P_y){)}^2}$

其中$NCC(p,d)$得到的值的范围将在[−1,1]之间。
$Wp$为之前提到的匹配窗口。
$I_1(x,y)$为原始图像的像素值。
$\overline{I_1}(p_x,p_y)$为原始窗口内像素的均值。
$I_2(x+d,y)$为原始图像在目标图像上对应点位置在$x$方向上偏移$d$后的像素值。
$\overline{I_2}(p_x+d,p_y)$为目标图像匹配窗口像素均值。

若$NCC=-1$，则表示两个匹配窗口完全不相关，相反，若$NCC=1$时，表示两个匹配窗口相关程度非常高。

2.2 匹配流程

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

2. 极线校正：校正的目的是使两帧图像极线处于水平方向，或者说是使两帧图像的光心处于同一水平线上。通过校正极线可以方便后续的$NCC$操作。
   1）由标定得到的内参中畸变信息中可以对图像去除畸变。
   2）通过校正函数校正以后得到相机的矫正变换R和新的投影矩阵P，接下来是要对左右视图进行去畸变，并得到重映射矩阵。

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

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

如图，$Tx$为双目相机基线，$f$为相机焦距，这些可以通过相机标定步骤得到。而$xr-xl$就是视差$d$。
通过公式$z=\frac{f\times Tx}{d}$可以很简单地得到以左视图为参考系的深度图了。
至此，我们便完成了双目立体匹配。倘若只是用于图像识别，那么到步骤3时已经可以结束了。

2.3 代码实现

NCCfaster.py

import numpy as np
import cv2
from PCV.localdescriptors import sift

im1 = 'C://Users//Garfield//PycharmProjects//untitled//NCC-master//im2.ppm'
im2 = 'C://Users//Garfield//PycharmProjects//untitled//NCC-master//im6.ppm'

img1 = cv2.imread(im1, cv2.CV_8UC1)
img2 = cv2.imread(im2, cv2.CV_8UC1)
rows, cols = img1.shape
print(img1.shape)

def translaton(image, shape):
    step = round((shape[0]-1)/2)
    print(step)
    shifted = []
    for i in range(0, step+1):
        for j in range(0, step+1):
            if i==0 and j==0:
                M1 = np.float32([[1, 0, i], [0, 1, j]])
                shifted.append(cv2.warpAffine(image, M1, (image.shape[1], image.shape[0])))
            elif i==0 and j!=0:
                M1 = np.float32([[1, 0, i], [0, 1, j]])
                M2 = np.float32([[1, 0, i], [0, 1, -j]])
                shifted.append(cv2.warpAffine(image, M1, (image.shape[1], image.shape[0])))
                shifted.append(cv2.warpAffine(image, M2, (image.shape[1], image.shape[0])))
            elif i!=0 and j==0:
                M1 = np.float32([[1, 0, i], [0, 1, j]])
                M2 = np.float32([[1, 0, -i], [0, 1, j]])
                shifted.append(cv2.warpAffine(image, M1, (image.shape[1], image.shape[0])))
                shifted.append(cv2.warpAffine(image, M2, (image.shape[1], image.shape[0])))
            else:
                M1 = np.float32([[1, 0, i], [0, 1, j]])
                M2 = np.float32([[1, 0, -i], [0, 1, -j]])
                M3 = np.float32([[1, 0, -i], [0, 1, j]])
                M4 = np.float32([[1, 0, i], [0, 1, -j]])
                shifted .append(cv2.warpAffine(image, M1, (image.shape[1], image.shape[0])))
                shifted.append(cv2.warpAffine(image, M2, (image.shape[1], image.shape[0])))
                shifted.append(cv2.warpAffine(image, M3, (image.shape[1], image.shape[0])))
                shifted.append(cv2.warpAffine(image, M4, (image.shape[1], image.shape[0])))

    print(len(shifted))
    return np.array(shifted)

#I(x,y)-avg(I(x,y))
def img_sub_avg(img_shifted, avg_img):
    len, height, width = img1_shifted.shape
    tmp_ncc1 = np.zeros([len, height, width])
    for i in range(len):
        tmp_ncc1[i] = img_shifted[i] - avg_img
    print(tmp_ncc1)
    return tmp_ncc1

def NCC(img1_sub_avg,img2_sub_avg, threshold, max_d):
    #设立阈值
    len, height, width = img1_sub_avg.shape
    thershould_shifted = np.zeros([len, height, width])
    ncc_max = np.zeros([height, width])
    ncc_d = np.zeros([height, width])
    for j in range(3, max_d):
        tmp_ncc1 = np.zeros([height, width])
        tmp_ncc2 = np.zeros([height, width])
        tmp_ncc3 = np.zeros([height, width])
        for k in range(len):
            M1 = np.float32([[1, 0, -j - 1], [0, 1, 0]])
            thershould_shifted[k] = cv2.warpAffine(img1_sub_avg[k], M1, (img1_sub_avg.shape[2], img1_sub_avg.shape[1]))
        for i in range(len):
            tmp_ncc1 += (img2_sub_avg[i])*(thershould_shifted[i])
            tmp_ncc2 += pow(img2_sub_avg[i], 2)
            tmp_ncc3 += pow(thershould_shifted[i], 2)

        tmp_ncc2 = tmp_ncc2*tmp_ncc3
        tmp_ncc2 = np.sqrt(tmp_ncc2)
        tmp_ncc4 = tmp_ncc1/tmp_ncc2
        for m in range(height):
            for n in range(width):
                if tmp_ncc4[m, n] > ncc_max[m ,n] and tmp_ncc4[m, n] > threshold:
                    ncc_max[m, n] = tmp_ncc4[m, n]
                    ncc_d[m , n] = j
    for i in ncc_d:
        print(i)
    return ncc_max, ncc_d

if __name__ == "__main__":
    disparity = np.zeros([rows, cols])
    NCC_value = np.zeros([rows, cols])
    deeps = np.zeros([rows, cols])
    # 用3*3卷积核做均值滤波
    avg_img1 = cv2.blur(img1, (7, 7))
    avg_img2 = cv2.blur(img2, (7, 7))
    fimg1 = img1.astype(np.float32)
    fimg2 = img2.astype(np.float32)
    avg_img1 = avg_img1.astype(np.float32)
    avg_img2  = avg_img2.astype(np.float32)
    img1_shifted = translaton(fimg1, [7, 7])
    img2_shifted = translaton(fimg2, [7, 7])
    img1_sub_avg = img_sub_avg(img1_shifted, avg_img1)
    img2_sub_avg = img_sub_avg(img2_shifted, avg_img2)
    ncc_max, ncc_d = NCC(img1_sub_avg,img2_sub_avg, threshold = 0.5, max_d = 64)
    print(img1_shifted.shape)
    disp = cv2.normalize(ncc_d, ncc_d, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX,
                          dtype=cv2.CV_8U)
    cv2.imshow("left", img1)
    cv2.imshow("right", img2)
    cv2.imshow("depth", disp)
    cv2.waitKey(0)  # 等待按键按下
    cv2.destroyAllWindows()#清除所有窗口

原始图像：

运行结果：


代码2：

# -*- coding: utf-8 -*-
import scipy.misc
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//Garfield//Desktop//towelmatch//jidian//1.jpg').convert('L'), 'f')
im_r = array(Image.open(r'C://Users//Garfield//Desktop//towelmatch//jidian//2.jpg').convert('L'),'f')
# 开始偏移，并设置步长
steps = 12
start = 4
# ncc 的宽度
wid = 2000
res = plane_sweep_ncc(im_l,im_r,start,steps,wid)

imsave('C://Users//Garfield//PycharmProjects//untitled//NCC-master//depth2000.png',res)
show()

此代码可以修改窗口值的大小，具体结果与分析放在下一小节

3. 不同窗口值对匹配结果的影响

原图：

以下分别是当窗口值设定为15，20，50，100，200，500，700，1000，2000的运行结果：

总的来看，窗口值设置的越大，区域化分得越为明显，而设置的越小噪声越大，区块划分的越不明显。同时，从本次实验结果图来看，wid设置为700时，结果最为理想。
另外可以看到

瓶身光照最为强烈的地方会与瓶身其他其他地方像素值有较大区分，以及窗口大小设置越大，则会更多的提取出有光照的像素点并连结成块，区分出亮部暗部。

4. 实验遇到的问题与解决

问题：

• ‘NoneType’ object has no attribute ‘shape’

求解cols,rows = img.shape时出现上述报错，是图片路径不对，程序无法获取到图片。

解决办法：
右键下图目录灰色处

—>copypath,再在原本的读取图片路径代码中粘贴
注意格式改成C://Users//…
参考博客：立体匹配_数据结构与算法
参考博客：立体匹配过程
参考博客:双目立体匹配算法–归一化互相关（NCC）详解和代码实现（python)
NCC的更详细相关
延伸学习：真实场景的双目立体匹配（Stereo Matching）获取深度图详解
原理详解：双目视觉（三）立体匹配算法

目录

1.立体视觉

2.视差

3.立体匹配

4.立体重建

5.归一化的互相关 (NCC)

6.实验

6.1 实验代码

6.2 实验分析

6.3 实验总结

1.立体视觉

一个多视图成像的特殊例子是立体视觉（或者立体成像），即使用两台只有水平（向一侧）偏移的照相机观测同一场景。当照相机的位置如上设置，两幅图像具有相同的图像平面，图像的行是垂直对齐的，那么称图像对是经过矫正的。该设置在机器人学中很常见，常被称为立体平台。

通过将图像扭曲到公共的平面上，使外极线位于图像行上，任何立体照相机设置都能得到矫正（通常构建立体平台来产生经过矫正的图像对）。

在已经确保两个摄像头的参数是完全一致的，并且两者的位置是平行之后，关注点就落到了如何计算物体的深度信息，这也是最重要最关键的地方。

标准立体视觉系统下的计算原理：这里所说的立体视觉指的是双目视觉，双目视觉是模拟人类视觉原理，使用计算机被动感知距离的方法。从两个或者多个点观察一个物体，获取在不同视角下的图像，根据图像之间像素的匹配关系，通过三角测量原理计算出像素之间的偏移来获取物体的三维信息。得到了物体的景深信息，就可以计算出物体与相机之间的实际距离，物体3维大小，两点之间实际距离。

实现一个立体视觉系统流程：

1. 获得图像
2. 相机标定
3. 图像校正
4. 立体匹配
5. 深度计算

2.视差

如上图，Z是深度，，也就是视差（disparity ）。

可以发现，深度 Z跟视差D成反比关系的，当视差D越小时，Z越大，物体离立体视觉系统也就越远， 当视差D越大，Z越小，物体离立体视觉系统也就越近。这一点和我们人眼系统是一样的，当我们观察离我们比较近的物体的时候，视差很大，可以获得的信息也就越多，当物体离我们很远的时候，视差很小，我们获得的信息也就很少了。

在图像处理中，我们通常用灰度值来表示视差信息，视差越大，其灰度值也就越大，在视差图像的视觉效果上表现出来就是图像越亮，物体离我们越远，其视差越小，灰度值也越小，视差图像也就越暗。

3.立体匹配

立体匹配主要是通过找出每对图像间的对应关系，根据三角测量原理，得到视差图；在获得了视差信息后，根据投影模型很容易地可以得到原始图像的深度信息和三维信息。立体匹配是建立立体视觉中最重要的一环，立体匹配的效果直接影响得到的三维信息。

对于图像的匹配是对两幅图片或者更多图片，找到它们相似的地方。而对于立体匹配而言，就是基于同一场景得到的多张二维图，通过找到相同点进一步还原场景的三维信息，一般采用的图像是双目图像，

基本的步骤：

1. 匹配代价计算
2. 代价聚合
3. 视差计算
4. 视差精化：对上一步得到的粗估计的视差图进行精确计算

立体匹配也分为全局匹配和局部匹配两种，这两种在匹配步骤上也有所差异。

全局匹配：

全局立体匹配算法会省略第二步，主要是采用了全局的优化理论方法估计视差，建立全局能量函数，通过最小化全局能量函数得到最优视差值。

全局匹配算法得到的结果比较准确，但是其运行时间比较长，不适合实时运行。主要的算法有图割、信念传播、动态规划等算法。

局部匹配：

给定在一幅图像上的某一点，选取该像素点邻域内的一个子窗口，在另一幅图像中的一个区域内，根据某种相似性判断依据，寻找与子窗口图像最为相似的子图，而其匹配的子图中对应的像素点就为该像素的匹配点。主要方法有SAD、SSD、NCC等。

4.立体重建

立体重建（又称致密深度重建）就是恢复深度图（或者相反，视差图），图像中每个像素的深度（或视差）都需要计算出来。

三维重建技术就是要在计算机中真实地重建出该物体表面的三维虚拟模型，构建一个物体完整的三维模型，大致可以分为三步：

1. 利用摄像机等图像采集设备对物体的点云数据从各个角度釆集，单个摄像机只能对物体的一个角度进行拍摄，要获得物体表面完整信息，需要从多个角度对物体拍摄；
2. 将第一步获得的各视角点云数据变换到同一个坐标系下，完成多视角点云数据的配准；
3. 根据配准好的点云数据构建出模型的网格表面。

流程如图：

5.归一化的互相关 (NCC)

归一化互相关匹配法NCC（Normalization cross correlation），是基于图像灰度信息的匹配方法，就是用于归一化待匹配目标之间的相关程度。归一化互相关匹配算法是图像匹配算法中较为经典的匹配算法。它是一种相似性度量或者匹配程度的表征，而不是一种图像匹配的完整方法，但是把互相关的思想作为度量测度，在许多匹配算法里都会用到。

在该立体重建算法中，将对于每个像素尝试不同的偏移，并按照局部图像周围归一化的互相关值，选择具有最好分数的偏移，然后记录下该最佳偏移。因为每个偏移在某种程度上对应于一个平面，所以该过程有时称为扫平面法。虽然该方法并不是立体重建中最好的方法，但是非常简单，通常会得出令人满意的结果。

当密集地应用在图像中时，归一化的互相关值可以很快地计算出来。这和在应用于稀疏点对应的不同，使用每个像素周围的图像块(根本上说，是局部周边图像）来计算归一化的互相关。

公式如下：

其中x为图像块的点，， 为图像块的均值。 

归一化的目的就是消除上述方法对于光照变化敏感的问题。  

优点：该方法步骤简单，较好的解决了对于光照变化敏感的问题。

缺点：矩形窗口的选用是一个缺憾。因此它只适合于具有平移和小角度旋转关系的图像配准。其次，它需要将模板图在搜索图上遍历所有可能的位置，当模板图很大的情况下，显然其计算量是相当大的。

6.实验

6.1 实验代码

计算视差图

# coding=utf-8
from PIL import Image
from pylab import *
import stereo
import cv2

im_l = array(Image.open('D:\\PycharmProjects\\ComputerVision\\ch06\\data\\im1.ppm').convert('L'), 'f')
im_r = array(Image.open('D:\\PycharmProjects\\ComputerVision\\ch06\\data\\im2.ppm').convert('L'), 'f')

# 开始偏移，并设置步长
steps = 12
start = 4

# ncc的宽度
wid = 3

# 标准版本
res = stereo.plane_sweep_ncc(im_l, im_r, start, steps, wid)
# 高斯版本
# res = stereo.plane_sweep_gauss(im_l, im_r, start, steps, wid)

import scipy.misc
scipy.misc.imsave('depth.png', res)

im = array(Image.open('D:\\PycharmProjects\\ComputerVision\\depth.png').convert('L'), 'f')
fig = plt.figure()
subplot(131)
imshow(im_l)
title('left')
axis('off')
subplot(132)
imshow(im_r)
title('right')
axis('off')
subplot(133)
imshow(im)
title('depth')
axis('off')
show()

均匀滤波器

# coding=utf-8
from pylab import *
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(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)

6.2 实验分析

6.2.1 使用均匀滤波器

（1）wid值不同

这两组图片都是使用均匀滤波器得到的，窗口值的取值都为1，3，5，7，9，11六个数，当窗口值为1时，图片中一片空白什么也没有，这是因为窗口太小，加上左右图片存在一定的偏移量，导致窗口匹配时，两张图片很难匹配到一起，所以得不到任何信息。

而随着窗口值增大为3时，视差图中的信息最多，是窗口匹配的最佳时候，窗口大小适宜，此时左右两张图的匹配信息最多。

接着窗口值再增大，视差图的信息就开始减少了，预测当窗口值达到一定数值时会得到和窗口值为1时的图像。所以窗口值的取值要适宜，不能过大或过小，所以当窗口值为3时，均匀滤波器得到最佳视差。

（2）wid = 3，steps取不同值

steps = 6

steps = 9

steps = 12

steps = 15

结论：在窗口值最佳的情况下，steps的取值大小对计算视差图并无明显影响 。

6.2.2 使用高斯滤波器

使用高斯滤波器得到的结果，当窗口值为1时就得到最佳结果，然后随着窗口值的增大，越来越少。

与均匀滤波器的最佳结果对比发现，高斯滤波器具有较少的噪声，但是缺少很多细节。

6.3 实验总结

• 窗口值的取值要适宜，在上面所选的样例图中，均匀滤波器的wid = 3最佳，高斯滤波器的wid = 1 最佳，均匀滤波器细节信息多，高斯滤波器噪声少，各有千秋。
• 深度跟偏移成正比，在视差图像的视觉效果上表现出来就是图像越亮，物体离我们越远，其视差越小，灰度值也越小，视差图像也就越暗。
• 深度跟视差成反比，距离较近时，视差很大，可以获得的信息也就越多，距离较近时，视差很小，我们获得的信息也就很少了。

一、NCC的基础概念

  NCC（normalized cross correlation）算法，归一化互相关匹配法，是基于图像灰度信息的匹配方法。

  图像匹配的方法主要有三种：基于灰度，基于特征，基于变换域。

二、公式介绍

  NCC算法的基础理论来可以讲，是将图像的相似性归结为2个向量的相似性。假如 a为向量1,b为向量2,根据点乘的定义得到

若 a与b相似,则它们的方向基本相同,其夹角近似为0,即cosθ≈1,这样,就可以根据cosθ的值判断2个向量的相似性。

  将其推广到多维的图像检测，则假设待搜索图像S的尺寸为M * M，模板r的尺寸为N * N（M>>N）。其中M，N代表图像象素。模板T在图像S上平移，搜索窗口所覆盖的子图记作S（i，j）．(i,j)为子图的左上角顶点在搜索图S中的坐标。通过相关函数计算子图与实时图的灰度相关值。对搜索图自上面下、自左面右遍历搜索，记录下每一个子图位置的互相关值。互相关值最大的子图位置即为匹配位置。

  在实际匹配应用中，搜索图和模板的相似性是通过度量函数来度量的，则归一化积相关匹配算法可定义为：

R(u,v)为点(u,v)处的NCC系数;M×N为匹配模板的大小;xi+u,j+v,yi,j分别为参加匹配的2幅图像中(i+u,j+v),(i,j)处的灰度值。R(u,v)的值越大,这2幅图像越相似。因此,可以根据R(u,v)值的大小判断2幅图像的相度。

{

每个尺寸点为一个单位向量，因此这个方程为每个点相对位置的两个点的的余弦值，这个值代表着两点的相似度。把这个余弦值累乘，就是两个图像的相似度了。

}

三、举例说明

  景象匹配实际上是比较 2 幅图像的相似性。也可以把图像展开为向量 , 这样 , 就可以归结为比较 2 个向量的相似性。根据向量点乘的定义a · b = | a | · | b | · cos θ . 若 2 个向量相似 , 则它们的方向相同 , 其夹角为 0 , 因此 , 可以根据 cos θ 的值来判断 2 个向量的相似性。把其推广到二维图像中 , 则下列式中

R (u, v) 为位置点 ( u, v) 的归一化互相关系数 ; N 1 ×N 2 为匹配模板的大小 ; x （i + u, j + v） , y（i, j） 分别为需匹配的 2 幅图像中( i+ u, j + v) , ( i, j) 处的灰度值。 R ( u, v) 的值越大 ,2 幅图像越相似。

四、NCC算法特点

  这种方法的优点是抗白噪声干扰能力强，且在灰度变化及几何畸变不大的情

况下精度很高，它的这种优点非常突出，但该方法受局部光照变化的影响，且匹

配速度较慢。

五、NCC闲谈

  在进行图片采集的时候，应该先进行聚焦。

  聚焦的步骤：首先，把图片设定为黑白图，然后通过粗聚焦，得到清晰地大概位置，然后再缩小范围（在粗聚焦得到的函数曲线的峰值附近）进行精聚焦，得到细聚焦函数曲线；

可用卷积积分来使函数更平滑。

  当遇到图像有旋转的时候如何呢？可以在原本目标中选一个参考区域，通过实时监控该区域内图像，让图像进行一定角度的实时转变，实现与实际监控图片的实时锁定。通过这种方法，可以得到整张图片的旋转角度，并且，极大的降低了运算量。

【作者】 孙卜郊； 周东华；

【Author】 SUN Bo-jiao,ZHOU Dong-hua(Department of Automation,Tsinghua University,Beijing 100084,China)

【机构】 清华大学自动化系；