一. 计算视差原理

1.1 立体图像

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

假设两幅图像经过了矫正，那么对应点的寻找限制在图像的同一行上。一旦找到 对应点，由于深度是和偏移成正比的，那么深度（Z 坐标）可以直接由水平偏移来 计算

其中，f 是经过矫正图像的焦距，b 是两个照相机中心之间的距离，xl 和 xr 是左右两 幅图像中对应点的 x 坐标。分开照相机中心的距离称为基线。矫正后的立体照相机：

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

1.2 计算视差图

在该立体重建算法中，我们将对于每个像素尝试不同的偏移，并按照局部图像周围 归一化的互相关值，选择具有最好分数的偏移，然后记录下该最佳偏移。因为每个偏移在某种程度上对应于一个平面，所以该过程有时称为扫平面法。

1.3 NCC视差匹配方法

当密集地应用在图像中时，归一化的互相关值可以很快地计算出来。这和我们在第 2 章中应用于稀疏点对应的不同。我们使用每个像素周围的图像块（根本上说，是 局部周边图像）来计算归一化的互相关。对于这里的情形，我们可以在像素周围重 新写出公式（2.3）中的== NCC==

1.4 双目立体匹配的步骤

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

2. 极线校正：校正的目的是使两帧图像极线处于水平方向，或者说是使两帧图像的光心处于同一水平线上。
   2.1 由标定得到的内参中畸变信息中可以对图像去除畸变
   2.2 通过校正函数校正以后得到相机的矫正变换R和新的投影矩阵P，接下来是要对左右视图进行去畸变，并得到重映射矩阵。
   2.3 根据上述得到的重映射参数map1,map2，我们需要进一步对原始图像进行重映射到新的平面中才能去除图像畸变。
   2.4 通过上述两步操作，我们成功地对图像去除了畸变，并且校正了图像极线。

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

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

二. 代码实现

2.1 图片集

左相机：

右相机：

2.2 代码实现

# -*- 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('11.png').convert('L'), 'f')
im_r = array(Image.open('12.png').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()

2.3 运行结果及分析

窗口值wid = 15 时

窗口值wid = 9 时

窗口值wid = 5 时

窗口值wid = 3 时

图片集2：

窗口值wid = 15 时：

窗口值wid = 9 时：

窗口值wid = 3时：

三. 实验结果及小结

1. 从实验结果可以看出，窗口值的大小改变对结果影响较大，wid值越小，匹配结果的区分度越低，并且成散粉分裂状；随着窗口值增大，匹配区分度逐渐清晰，匹配精度变高。
2. 当窗口值过大时，效果也并没有更好，暗区会越来越黑，导致匹配效果也不太好。
3. 实验遇到的错误： numpy版本问题，依旧需要卸载重装方可解决；使用numpy里的函数需要在函数名前面加np.前缀方可。