双目立体视觉 - David-ustc's blog

双目立体视觉———《人工智能导论之computer vision》实验报告

实验题目

在middlebury数据集上实现双目立体视觉图像重建

实验目的

了解计算机视觉基本知识。
学习计算机数字图像处理基本概念。
双目视觉图像重建算法的学习与实践。

实验原理

双目立体视觉图像
大体步骤有：双目标定 $\Rightarrow$ 立体校正（含消除畸变）$\Rightarrow$ 立体匹配 $\Rightarrow$ 视差计算 $\Rightarrow$ 深度/3D坐标计算
第三、四步包括：
(1)必须从一幅图像中选出位于场景中一个表面上的某一特定位置；
(2)必须在另一幅图像中鉴别出同一个位置；
(3)测出这两个对应像点之间的视差。在一组重复摄影的两张照片上(立体模式)，同一地物的影象，沿着摄影基线(摄影地点和下一个摄影地点之间的飞行方向线)方向位置变换，这个变化量叫立体视差。
双目标定
双目标定的目标是获得左右两个相机的内参、外参和畸变系数，其中内参包括左右相机的fx，fy，cx，cy，外参包括左相机相对于右相机的旋转矩阵和平移向量，畸变系数包括径向畸变系数（k1， k2，k3）和切向畸变系数（p1，p2）。实验采用的middlebury数据集，包含校正效果已经很好的图像和参数表：
Here is a sample calib.txt file for one of the full-size training image pairs:
cam0=[3997.684 0 1176.728; 0 3997.684 1011.728; 0 0 1]
cam1=[3997.684 0 1307.839; 0 3997.684 1011.728; 0 0 1]
doffs=131.111
baseline=193.001
width=2964
height=1988
ndisp=280
isint=0
vmin=31
vmax=257
dyavg=0.918
dymax=1.516
极线校正
立体校正的目的是将拍摄于同一场景的左右两个视图进行数学上的投影变换，使得两个成像平面平行于基线，且同一个点在左右两幅图中位于同一行，简称共面行对准。只有达到共面行对准以后才可以应用三角原理计算距离。
如下图所示：
立体匹配与视差图计算
假设两幅图像经过了校正，那么对应点的寻找限制在图像的同一行上。一旦找到对应点，由于深度是和偏移成正比的，那么深度（Z坐标）可以直接由水平偏移来计算，
其中，$f$是经过校正图像的焦距，$b$(baseline)是两个照相机中心之间的距离，$x_l,x_r$是左右两幅图像中对应点的横坐标，$x_l - x_r = disparity$为视差。
根据上图可以得到相似三角形：$\Delta_{p_l, p, p_r} \sim \Delta_{O_l, p, O_r}$
以及对应的性质：
$$\frac{b+ x_l - x_r}{b} = \frac{Z-f}{Z}$$
计算公式在校正后必要时加上修正，即$d = x_l - x_r +doffs=x_l-x_r+cx_0-cx_1$。

双目视觉的重点在于匹配算法，接下来介绍几种基础的普通计算方法，并选取部分实践评估。常用的立体匹配方法基本上可以分为两类：局部方法例如BM、SGM等，非局部的，即全局方法例如Dynamic Programming、Graph Cut等。局部方法计算量小，匹配质量相对较低，全局方法省略了代价聚合而采用了优化能量函数的方法，匹配质量较高，但是计算量也比较大。目前OpenCV中已经实现的方法有BM、binaryBM、SGBM这几种方法。
此外在立体匹配生成视差图之后，还可以对视差图进行滤波后处理，例如Guided Filter、Bilatera Filter等。视差图滤波能够将稀疏视差转变为稠密视差，并在一定程度上降低视差图噪声，改善视差图的视觉效果，但是比较依赖初始视差图的质量。

匹配像素算法———归一化互相关(NCC:Normalized cross-correlation)
在立体重建算法中，我们将对于每个像素尝试不同的偏移，并按照局部图像周围归一化的互相关值，选择具有最佳评估结果的视差。因为每个偏移在某种程度上对应于一个平面，所以该过程有时称为扫平面法。
如有校正过的两帧图像，NCC算法对图像$I_1$一个待匹配像素构建$n*n$匹配窗口，在图像$I_2$极线上对每一个像素构建匹配窗口与待匹配像素匹配窗口计算相关性，相关性最高的视为最优匹配。由于NCC匹配流程是通过在同一行查找最优匹配，因此可以并行处理，对于运行效率有一定的提升。
我们使用每个像素周围的图像块（根本上说，是局部周边图像）来计算归一化的互相关。针对图像周围像素NCC公式：
$$ncc(I_1, I_2) =\frac{\sum_{x}(I_1(x) -\mu_1)(I_2(x)- \mu_2)}{\sqrt{\sum_{x}(I_1(x) -\mu_1)^2 \sum_{x}(I_2(x)- \mu_2)^2}}$$
匹配像素算法————SGBM算法
(1).预处理
SGBM采用水平Sobel算子，把图像做处理，公式为：
$$Sobel (x,y)=2[P(x+1,y)-P(x-1,y)]+ P(x+1,y-1)-P(x-1,y-1)+ P(x+1,y+1)-P(x-1,y+1)$$
用一个函数将经过水平Sobel算子处理后的图像上每个像素点（P表示其像素值）映射成一个新的图像：P_new表示新图像上的像素值，prefiltercap作为传递参数可定义。
$$P_{new} = \begin{cases}
```
\qquad\qquad 0 \qquad\qquad p <\text{-prefiltercap};\\
p+prefiltercap \qquad \text{-prefiltercap}\leq p \leq \text{-prefiltercap};\\
 2*prefiltercap,\quad\quad p >\text{-prefiltercap}.
\end{cases}
```
$$
预处理实际上是得到图像的梯度信息,经预处理的图像保存起来，将会用于计算代价。

(2).代价计算
大部分匹配算法中就是在校正图像上的每一行，寻找一条使得全局能量函数最小的最优匹配路径。在论文中，作者提出了一个具体的能量函数：
$$L(D) = \sum_p [C(p,D_p) +\sum_{q\in N_p}p_1\cdot I(|D_p -D_q|=1)+\sum_{q\in N_p}p_2 I(|D_p -D_q|>1)$$
opencv中两个参数P1，P2是这样设定的：
P1 =8*cn*sgbm.SADWindowSize*sgbm.SADWindowSize;
P2 = 32*cn*sgbm.SADWindowSize*sgbm.SADWindowSize。
Np 指像素p的相邻像素点；C(p, Dp)指当前像素点disparity为Dp时的cost函数值，可以为BT代价或者MI代价等，SGBM中能量函数一般定义为C(original)和C(smooth)之和，其中 C(data)为图像数据项用于判断匹配像素点之间的相似性，C(smooth)为sobel预处理之后像素点与邻域像素的平滑约束项，用于保证视差的连续性；
P1 是一个惩罚系数，它适用于像素p相邻像素中disparity值与p的disparity值相差1的那些像素；P2 是一个惩罚系数，它适用于像素p相邻像素中disparity值与p的disparity值相差大于1的那些像素；$I[]$示性函数如果函数中的参数为真返回1，否则返回0。
按照上述过程，只要依次求解每个阶段的最优值$min L_k()$，该行像素的视差也就求出来了。从这个过程中可以发现，该计算流程实际上就是接下来所要讲解的SGM在水平方向上的聚合。
(3).代价聚合和动态规划
动态规划即在水平扫描线的视差空间切面上寻找最优路径，如下图所示，以像素点的行方向为横坐标，视差值d 为纵坐标，依次将整个过程分为 1、2、3、…、k 共k个阶段。将上述各个阶段所处的匹配阶段用不同的状态表示。如下图所示，每一个阶段的状态就是对应像素点的匹配情况，共有三种状态：相互匹配记为 M、左可见右遮挡为 L、右可见左遮挡为 R，状态的选择满足无后效性。
针对该思想，Hirschmuller提出了semi-global matching方法，通过在多个方向1维路径上平等地进行代价聚合，然后使用WTA求解视差，来作为一个近似求解2维能量最小化的过程。此时能量（损失）函数化简成：
$$ L(r,d) = C(\bold{p},d) +min(L_r(\bold{p-r},d),L_r(\bold{p-r},d-1)+P_1,L_r(\bold{p-r},d+1)+P_1,min_iL_r(\bold{p-r},i)+P_2)-min_kL_r(\bold{p-r},k)$$
代价聚合的路径可以为8或16个，这样可以全面覆盖整张图像，最终的聚合代价是所有路径上聚合代价的和：$S(\bold{p},d) = \sum_r L_r(\bold{r},d)$
(4).视差计算
经过代价聚合后，根据winner takes all的思想，取最大值$d^* = argmin S(\bold{p},d)$得到视差图，即可进一步得到深度图。
(5).后处理proprocessing，这里只选取两步骤简要介绍。
A.唯一性检测
采用比率测试匹配筛选方法，即最低代价是次低代价的(1+ uniquenessRatio/100)倍的时候时，最低代价匹配的视差值才认为是正确视差值，否则这个像素点的视差值设定为负值。其中uniquenessRatio是一个常数参数，一般设置为0.8左右。
B.左右一致性检测
设左图中像素p的视差为d，右图中对应像素的视差为d，那么一般认为时，$|d-d^|<threshold$才最有可能为正确的视差，否则将其置为一个无效值。threshold一般取为1或者2。函数cv2.stereoSGBM_create()代码中已经介绍过该参数。

算法结构

ncc算法有两种处理方法，分别是采用平均滤波器和高斯滤波器计算归一化互相关系数ncc。两者除了filter不同，在算法上无差异，流程如下：

def plane_sweep_ncc(im_l,im_r,start,steps,wid):
    mean_l = filter(im_l)
    mean_r = filter(im_r)
    norm_l = im_l - mean_l
    norm_r = im_r - mean_r
    for displ in range(steps):
        s = filter(norm_l.roll(-displ - start) * norm_r)  
        s_l = filter(norm_l.roll(-displ - start) * norm_l.roll(-displ - start))
        s_r = filter(norm_r * norm_r) 
        dmaps[:, :, displ] = s / sqrt((s_l * s_r)+1e-5)
    #选取最佳深度
    return disparity = dmaps.argmax(, axis=2)

Python语言程序包opencv（cv2）有与C++的接口，内部定义了SGBM对象，可直接操作。

import cv2
stereo = cv2.StereoSGBM_create(
    #最小视差，默认为0。此参数决定左图中的像素点在右图匹配搜索的起点。
    minDisparity=0,
    numDisparities=256,  # max_disp has to be dividable by 16 f. E. HH 192, 256
    blockSize=3,
    #控制视差变化平滑性的参数。P1越大，视差越平滑，P2越大，边缘越差。windowsize默认3; 5; 7 ; 本实验设置为3
    P1=8 * 3 * window_size ** 2,
    P2=32 * 3 * window_size ** 2,
    #一致性
    disp12MaxDiff=1,
    #唯一性
    uniquenessRatio=15,
    #视差连通区域像素点个数的大小。对于每一个视差点，当其连通区域的像素点个数小于size时，认为该视差值无效，是噪点。
    speckleWindowSize=0,
    speckleRange=2,
    #水平sobel预处理后，映射滤波器大小。默认为15
    preFilterCap=63,
    mode=cv2.STEREO_SGBM_MODE_SGBM_3WAY
)
    disparity = stereo.compute(img_left, img_right).astype(np.float32)

结果分析

middlebury数据集提供了真实视差图图pfm文件，通过文件IO操作将字节序稍作变换，即可得到png图片。视差图作为数组读入后，运用一步公式和calib.txt标定文件参数，计算可得深度图。本小节全部以Umbrella-perfect为分析对象，视差图和深度图真值如下：
ncc算法提出时期较早，效果不如新算法好，主要通过调节wid参数实现最佳匹配。wid越大，图像越平滑，但细节信息更少，wid越小，噪声越多，缺乏稳健性。下图分别是平均滤波器（wid=9）和高斯滤波器(wid=3)的视差图，可以明显看出物体轮廓不清晰。
OpenCV的SGBM算法出彩之处在于后处理，受噪音、图像轮廓等因素影响小，从而成为主流的立体视觉匹配算法。实验成图和数据集真值图对比如下，虽然仍有差距，但是视觉上的匹配效果差强人意。

真值视差	SGBM视差
真值深度	SGBM深度

参考文献：
[1]vision.middlebury.edu/stereo/data/2014
[2]Hirschm. Stereo Processing by Semiglobal Matching and Mutual Information[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2007, 30(2):328-341.
[3]blog.csdn.net/dulingwen/article/details/104142149
[4]blog.csdn.net/weixin_45617915/article/details/105763749

class report

SSE可搜索加密上一篇

eigenface方法下一篇