文章目录
Overview
1. 图像采集
2. 双目标定
3. 双目矫正
stereoRectify
CameraInfo DKRP
initUndistortRectifyMap
Remap
4. 立体匹配
视差计算
5. 三维重建
深度图 图像类型
总结

Overview
双目立体视觉的整体流程包括：

图像采集
双目标定
双目矫正
立体匹配
三维重建


  

1. 图像采集
双目相机采集 左右目图像
2. 双目标定
通过 双目标定工具 对双目相机进行标定，得到如下结果参数：




内参
外参




相机矩阵 $K_1, K_2$
旋转矩阵 $R$


畸变系数 $D_1, D_2$
平移向量 $t$


《Learning OpenCV》中对于 Translation 和 Rotation 的图示是这样的:

  

示例代码：
cv::Matx33d K1, K2, R;
cv::Vec3d T;
cv::Vec4d D1, D2;

int flag = 0;
flag |= cv::fisheye::CALIB_RECOMPUTE_EXTRINSIC;
flag |= cv::fisheye::CALIB_CHECK_COND;
flag |= cv::fisheye::CALIB_FIX_SKEW;

cv::fisheye::stereoCalibrate(
        obj_points_, img_points_l_, img_points_r_,
        K1, D1, K2, D2, img_size_, R, T,
        flag, cv::TermCriteria(3, 12, 0));

3. 双目矫正
双目矫正 主要包括两方面：畸变矫正 和 立体矫正 。
利用 OpenCV的函数，主要分为

stereoRectify
initUndistortRectifyMap
remap

stereoRectify
根据双目标定的结果 $K_1, K_2, D_1, D_2, R, t$，利用 OpenCV函数 stereoRectify，计算得到如下参数

左目 矫正矩阵(旋转矩阵) $R_1$ (3x3)
右目 矫正矩阵(旋转矩阵) $R_2$ (3x3)
左目 投影矩阵 $P_1$ (3x4)
右目 投影矩阵 $P_2$ (3x4)
disparity-to-depth 映射矩阵 $Q$ (4x4)

其中，
左右目投影矩阵（horizontal stereo, ${c_x}_1'={c_x}_2'$ if CV_CALIB_ZERO_DISPARITY is set）
$P_1 =
	\begin{bmatrix}
	f' & 0 & {c_x}_1' & 0 \\
	0 & f' & c_y' & 0 \\
	0 & 0 & 1   & 0
	\end{bmatrix}$
$P_2 =
	\begin{bmatrix}
	f' & 0 & {c_x}_2' & t_x' \cdot f' \\
	0 & f' & c_y' & 0 \\
	0 & 0 & 1   & 0
	\end{bmatrix}$
where
$t_x' = -B$
disparity-to-depth 映射矩阵
$Q =
	\begin{bmatrix}
	1 & 0 & 0 & -{c_x}_1' \\
	0 & 1 & 0 & -c_y'     \\
	0 & 0 & 0 & f'        \\
  0 & 0 & -\frac{1}{t_x'} & \frac{ {c_x}_1'-{c_x}_2'}{t_x'}
	\end{bmatrix}$
通过 $P_2$ 可计算出 基线 长度:
$\begin{aligned}
    \text{baseline} = B = - t_x' = - \frac{ {P_2}_{(03)} }{f'}
\end{aligned}$
示例代码：
cv::Mat R1, R2, P1, P2, Q;
cv::fisheye::stereoRectify(
        K1, D1, K2, D2, img_size_, R, T,
        R1, R2, P1, P2, Q,
        CV_CALIB_ZERO_DISPARITY, img_size_, 0.0, 1.1);

CameraInfo DKRP
参考：sensor_msgs/CameraInfo Message


D: distortion parameters.

For “plumb_bob”, the 5 parameters are: (k1, k2, t1, t2, k3)



K: Intrinsic camera matrix for the raw (distorted) images.

Projects 3D points in the camera coordinate frame to 2D pixel coordinates using the focal lengths (fx, fy) and principal point (cx, cy).

$\mathbf{K} =
\begin{bmatrix}
f_x & 0 & c_x \\
0 & f_y & c_y \\
0 & 0 & 1
\end{bmatrix}$



R: Rectification matrix (stereo cameras only).

A rotation matrix aligning the camera coordinate system to the ideal stereo image plane so that epipolar lines in both stereo images are parallel.
For monocular cameras $\mathbf{R} = \mathbf{I}$



P: Projection/camera matrix.

For monocular cameras

$\mathbf{P} =
  \begin{bmatrix}
  f_x & 0 & c_x & 0 \\
  0 & f_y & c_y & 0 \\
  0 & 0 & 1 & 0
  \end{bmatrix}$
For a stereo pair, the fourth column [Tx Ty 0]’ is related to the position of the optical center of the second camera in the first camera’s frame. We assume Tz = 0 so both cameras are in the same stereo image plane.

The first camera

$\mathbf{P} =
  \begin{bmatrix}
  f_x' & 0 & c_x' & 0 \\
  0 & f_y' & c_y' & 0 \\
  0 & 0 & 1 & 0
  \end{bmatrix}$
The second camera

$\mathbf{P} =
  \begin{bmatrix}
  f_x' & 0 & c_x' & -f_x' \cdot B \\
  0 & f_y' & c_y' & 0 \\
  0 & 0 & 1 & 0
  \end{bmatrix}$


Given a 3D point $[X Y Z]'$, the projection $(x, y)$ of the point onto the rectified image is given by:

$\begin{bmatrix} u \\ v \\ w \end{bmatrix} =
\mathbf{P} \cdot \begin{bmatrix} X \\ Y \\ Z \\ 1 \end{bmatrix}
, \quad
\begin{cases} x = \frac{u}{w} \\ y = \frac{v}{w} \end{cases}$


initUndistortRectifyMap
左右目 分别利用 OpenCV函数 initUndistortRectifyMap 计算 the undistortion and rectification transformation map，得到

左目map: $map^l_1, map^l_2$
右目map: $map^r_1, map^r_2$

示例代码：
cv::fisheye::initUndistortRectifyMap(K1, D1, R1, P1, img_size, CV_16SC2, rect_map_[0][0], rect_map_[0][1]);
cv::fisheye::initUndistortRectifyMap(K2, D2, R2, P2, img_size, CV_16SC2, rect_map_[1][0], rect_map_[1][1]);

Remap
左右目 分别利用 OpenCV函数 remap 并根据 左右目map 对左右目图像进行 去畸变 和 立体矫正，得到 左右目矫正图像
示例代码：
cv::remap(img_l, img_rect_l, rect_map_[0][0], rect_map_[0][1], cv::INTER_LINEAR);
cv::remap(img_r, img_rect_r, rect_map_[1][0], rect_map_[1][1], cv::INTER_LINEAR);

4. 立体匹配
根据双目矫正图像，通过 BM或SGM等立体匹配算法 对其进行立体匹配，计算 视差图

  

视差计算
通过 OpenCV函数 stereoBM (block matching algorithm)，生成 视差图(Disparity Map) (CV_16S or CV_32F)

disparity map from stereoBM of OpenCV :

It has the same size as the input images. When disptype == CV_16S, the map is a 16-bit signed single-channel image, containing disparity values scaled by 16. To get the true disparity values from such fixed-point representation, you will need to divide each disp element by 16. If disptype == CV_32F, the disparity map will already contain the real disparity values on output.

So if you’ve chosen disptype = CV_16S during computation, you can access a pixel at pixel-position (X,Y) by: short pixVal = disparity.at<short>(Y,X);, while the disparity value is float disparity = pixVal / 16.0f;; if you’ve chosen disptype = CV_32F during computation, you can access the disparity directly: float disparity = disparity.at<float>(Y,X);

Disparity Map
Disparity map post-filtering

5. 三维重建
（1）算法1：根据视差图，利用 $f'$ 和 $B$ 通过几何关系计算 深度值，并利用相机内参计算 三维坐标

  

根据上图相似三角形关系，得
$\frac{Z}{B} = \frac{Z-f}{B-d_w}
\quad \Longrightarrow \quad
Z = \frac{Bf}{d_w}$
其中，$f$ 和 $d_w$ 分别为 成像平面的焦距和视差，单位均为 物理单位，将其转换为 像素单位，上式写为
$Z = \frac{B f'}{d_p}$
其中，
$d_p = (O_r - u_r) + (u_l - O_l) = (u_l - u_r) + (O_r - O_l)$
最终，深度计算公式如下，通过遍历图像可生成 深度图
$Z = \text{depth} = \frac{B \cdot f'}{d_p}
\quad \text{with} \quad
d_p = \text{disp}(u,v) + ({c_x}_2' - {c_x}_1')$
根据 小孔成像模型，已知 $Z$ 和 相机内参 可计算出 三维点坐标，从而可生成 三维点云
$\begin{aligned}
	\begin{cases}
	Z = \text{depth} = \frac{f' \cdot B}{d_p} \\
	X = \frac{u-{c_x}_1'}{f'} \cdot Z \\
	Y = \frac{v-{c_y}'}{f'} \cdot Z
	\end{cases}
\end{aligned}
\text{或}
\begin{aligned}
	\begin{cases}
  \text{bd} = \frac{B}{d_p}\\
	Z = \text{depth} = f' \cdot \text{bd} \\
	X = (u-{c_x}_1') \cdot \text{bd} \\
	Y = (u-{c_y}') \cdot \text{bd}
	\end{cases}
\end{aligned}$
其中，$\text{disp}(u,v)$ 代表 视差图 坐标值
（2）算法2：根据视差图，利用 $Q$ 矩阵 计算 三维点坐标（reprojectImageTo3D）
$\begin{bmatrix} X' \\ Y' \\ Z' \\ W \end{bmatrix} =
Q \cdot
\begin{bmatrix} u \\ v \\ \text{disp}(u,v) \\ 1 \end{bmatrix}$
最终，三维点坐标为
$\begin{bmatrix} X \\ Y \\ Z \end{bmatrix} =
\begin{bmatrix}
\frac{X'}{W} \\[2ex]
\frac{Y'}{W} \\[2ex]
\frac{Z'}{W}
\end{bmatrix}$
深度图 图像类型

单位meter --> 32FC1
单位millimeter --> 16UC1

总结

  

ZYNQ上基于FPGA加速的 双目立体视觉 三维重建 生成 实时深度图彩虹图 （zedboard板）
1.在ZYNQ上完成双目采集成像，PS+PL协同处理。

2.在ZYNQ上完成相机标定，镜头校准。

3.在ZYNQ上完成立体匹配算法，目前已实现局部立体匹配BM以及半全局立体匹配SGBM

4.在ZYNQ上完成了立体匹配后的深度图转换彩虹图，即彩色深度图

5.在ZYNQ上也完成了测距算法，但是目前没有识别算法，只能取图像一点计算距离。精度在10mm内。有效距离目前测试在5m内。

6.优化了所有的算法，目前接近720p的图像质量，并且帧率在30fps左右。
详细的实现过程会出一个专栏，每周更一次，从零开始的项目搭建ZYNQ上的实时三维重建。
部分实验过程记录：
标定过程：



得到相机的内参矩阵和外参矩阵


校准与立体匹配：

伪彩色化：（深度彩虹图）

运行测试图：

转载 双目视觉三维重建 作者tiemaxiaosu
http://blog.csdn.net/tiemaxiaosu/article/details/51734667
一、三维重建概述
         三维重建主要是研究如何从得到的匹配点中计算出相机的投影矩阵(如果是外部标定的话，就是求出相机的外部参数)以及如何计算出匹配点的三维坐标。
        目前研究的进展与之还相差较远。研究人员为了能够表达三维空间信息，目前较多地采用三维矢量图形来替代三维位图。主要的重构方法有如下几种：
        (1)、空间点的重建
        这是三维重构中的最基本的方法。空间任一点在两个摄像机中分别成像，得到该点在两个图像中的对应坐标，在知道两摄像机的参数矩阵的条件下，通过建立以该点的世界坐标为未知数的 4 个线性方程，可以用最小二乘法求解得该点的世界坐标。有一种简化计算的办法，选取两个内部参数完全相同的摄像机，将他们平行配置，使他们的光轴相互平行。这种方法计算相当简单，不需要进行外部参数标定。但一般情况下，两个摄像机的内部参数完全相同(包括两摄像机的拍摄状态也相同)是很难选到。摄像机安装时无法看到光轴和成像平面，故难以实施。
       (2)、空间直线、空间二次曲线的重建
       空间直线和空间曲线也是组成空间图形的主要基元。根据空间直线射影变换前后仍是直线的性质，以此来分析空间直线与图像中的直线之间的关系。空间直线在CCD上的成像可以这样认为，空间直线与摄像机光心构成的平面与成像平面的交线，两台摄像机有两个这样的平面，这两个平面的交线就是空间直线。空间二次曲线是空间二次曲面与平面的交线，因此二次曲线的平面曲线，在CCD上的成像可以认为是由空间二次曲线与光心组成的锥面与成像平面的交线。两台摄像机有两个这样的锥面，求这两个空间锥面的交线，就是空间曲线。
       (3)、全像素的三维重建
       对图像中的每个像素都进行三维重建，这当然是最理想的情况，目前要实施这种重构，需要相当严格的测量条件，且仅限于对某一具体对象，但效果并不理想。对此人们正在进行大量的研究，还有待于进一步的努力。
二、空间点的三维重建
       1、基本模型
       空间点三维重建的基本模型如图(4.1)所示，对于空间物体表面任意一点P，如果用C1摄像机观察，看到它在C1摄像机的图像点位于p1，但他们无法由p1得知P的三维位置，事实上，在O1P(O1为C1摄像机的光心)连线上任意一点P’的图像点都是p1，因此，由p1点的位置，我们只知道空间点位于O1P1与O2P2两条直线的交点，即它们的三维位置是唯一确定的。
                    
       假如我们能够得到物体表面上所有点的三维坐标，则三维物体的形状与位置就是唯一确定的，在某些简单场合，例如三维物体时一个多面体，我们只需要知道它的各个顶点的三维坐标与相邻关系，则该多面体的形状与位置是唯一确定的。因此，用立体视觉的方法获取三维点的坐标是十分基础的。
       2、最小二乘法求解三维坐标值
       下面我们来具体介绍求解三维点坐标的算法。在进行讨论之前，我们要明确一下几点基本假设：空间任意点P在两个摄像机C1与C2上的图像点平p1与p2已经从两个图像中分别检测出来。即已知p1与p2为空间同一点P的对应点。至于对应点是如何找到的，这属于图像配准问题。还有一个假设就是，摄像机C1与C2已标定，他们的投影矩阵分别为M1与M2，于是有：
                   
        其中，(u1,v1,1)与(u2,v2,1)分别为P1与P2点在各自图像中的图像齐次坐标；(X,y,Z,1)为P点在世界坐标系下的齐次坐标；mijk为Mk的第i行第j列元素。式(4.1)与式(4.2)消去Zc1与Zc2得到关于X,Y,Z的四个线性方程：
              
       由于空间点P是O1p1与O2p2的交点，它必然同时满足式(4.3)，式(4.4)，我们可以联立这四个式子，从而求出P点的坐标(X,Y,Z)。四个方程，三个未知数，也就是说这个方程组只有三个独立的方程，这是因为我们已经假设P1与P2点是空间同一点P的对应点，因此已经假设了直线O1p1与O2p2一定相交，或者说，这四个方程必定有解，而且解是唯一的。而实际应用中，由于数据总是有噪声的，我们可以用最小二乘法求出P的三维点坐标(X,Y,Z)。
       3、视差测距法
       以上我们介绍了空间点三维重建的一般方法，为了更清楚的了解上述过程的物理意义，下面我们讨论在一种简单的摄像机配置下的空间点重建方法，视差测距法。此方法要求选用两个内部参数完全相同的摄像机，将它们平行放置，使它们的光轴相互平行，另有一对坐标轴共线，两个成像平面共面，两摄像机的光心有一个固定的距离d，这样求解图像点的世界坐标时只涉及摄像机内参数。原理如图(4.2)所示：
                
       在上述特殊摄像机配置下，假设C1坐标系为O1x1y1z1，C2坐标系为O2x2y2z2，焦距为f，摄像机间距为d，任何空间点P的坐标在C1下为(x1,y1,z1)，在C2坐标系下为(x2,y2,z2)，在左摄像机中的图像点左边为(u1,v1)，在右摄像机中的图像点坐标为(u2,v2)。
       根据式(4.5)，式(4.6)的摄影比例关系，对于左右摄像机有如下关系：
                       
       世界坐标系与摄像机坐标系的关系可表述如下：
                        
       联立上两式可得：
                        
       进而可得：
                       
       由式(4.5)和式(4.8)可计算的空间点的三维坐标：
                 
      式(4.9)即为双目立体视觉重构三维空间点的方法，可见，通过图像对视差和像点坐标就可以恢复出物体的三维坐标。但是由于对设备配置的要求比较苛刻，只是被应用在一些简单的视觉系统。本文采用的还是针对一般摄像机参数的三维重建方法，这样大大扩展了带云台的摄像机的视野范围，从而增强了移动机器人视觉系统的灵活性。


。



 

双目 MATLAB标定 ，查阅博主的【计算机视觉】摄像机标定 matlab toolbox_calib工具箱（单目标定和双目标定）

1  基本原理



 

得到了立体标定参数之后，就可以把参数放入xml文件，然后用cvLoad读入OpenCV了。具体的方法可以参照Learning OpenCV第11章的例子，上面就是用cvSave保存标定结果，然后再用cvLoad把之前的标定结果读入矩阵的



2.14. xml文件示例

 

 

这里需要注意的是Matlab标定结果中的om向量，这个向量是旋转矩阵通过Rodrigues变换之后得出的结果，如果要在cvStereoRectify中使用的话，需要首先将这个向量用cvRodrigues转换成旋转矩阵。关于Rodrigues变换，Learning OpenCV的第11章也有说明。

 



2.15. 旋转矩阵的Rodrigues形式表示

 

 

 

 

 

 

3. 立体校准和匹配

有了标定参数，校准的过程就很简单了。

我使用的是OpenCV中的cvStereoRectify，得出校准参数之后用cvRemap来校准输入的左右图像。这部分的代码参考的是Learning OpenCV 十二章的例子。

 

校准之后，就可以立体匹配了。立体匹配OpenCV里面有两种方法，一种是Block Matching，一种是Graph Cut。Block Matching用的是SAD方法，速度比较快，但效果一般。Graph Cut可以参考Kolmogrov03的那篇博士论文，效果不错，但是运行速度实在是慢到不能忍。所以还是选择BM。

 

以下是我用BM进行立体匹配的参数设置

 

 

[cpp] view plain copy

BMState = cvCreateStereoBMState(CV_STEREO_BM_BASIC,0);  
	
assert(BMState != 0);  
	
BMState->preFilterSize=13;  
	
BMState->preFilterCap=13;  
	
BMState->SADWindowSize=19;  
	
BMState->minDisparity=0;  
	
BMState->numberOfDisparities=unitDisparity*16;  
	
BMState->textureThreshold=10;  
	
BMState->uniquenessRatio=20;  
	
BMState->speckleWindowSize=13;  
 

 

其中minDisparity这个参数我设置为0是由于我的两个摄像头是前向平行放置，相同的物体在左图中一定比在右图中偏右，如下图3.1所示。所以没有必要设置回搜的参数。

如果为了追求更大的双目重合区域而将两个摄像头向内偏转的话，这个参数是需要考虑的。

 



3.1. 校正后的左右视图

 

 

另外需要提的参数是uniquenessRatio，实验下来，我感觉这个参数对于最后的匹配结果是有很大的影响。uniquenessRatio主要可以防止误匹配，其主要作用从下面三幅图的disparity效果比对就可以看出。在立体匹配中，我们宁愿区域无法匹配，也不要误匹配。如果有误匹配的话，碰到障碍检测这种应用，就会很麻烦。

 



3.2. UniquenessRatio为0时的匹配图，可以看到大片的误匹配区域

 



3.3. UniquenessRatio为10时的disparity map, 可以看到误匹配被大量减少了, 但还是有噪点

 



3.4. UniquenessRatio为20时的disparity map, 可以看到误匹配基本被去除了, 点云干净了很多

 

 

 

关于cvFindStereoCorrespondenceBM这个函数的源代码，曾经做过比较详细的研究，过一段时间也会把之前写的代码注释整理一下，发篇博文。

 

 

 

 

4. 实际距离的测量

在用cvFindStereoCorrespondenceBM得出disparity map之后，还需要通过cvReprojectImageTo3D这个函数将单通道Disparity Map转换成三通道的实际坐标矩阵。

具体的数学原理可以参考下面这个公式(from chenyusiyuan http://blog.csdn.net/chenyusiyuan/archive/2009/12/25/5072597.aspx,实际深度的一些问题这篇博文中也有提到)



4.1 距离转换公式

 

 

但是在实际操作过程中，用cvReprojectImageTo3D得到的数据并未如实际所想，生成深度矩阵所定义的世界坐标系我就一直没弄清楚。这在下面的例子中会详细说明，希望这方面的专家能帮忙解答一下：

 

 

图4.2是测量时的实际场景图，场景中主要测量的三个物体就是最前面的利乐包装盒、中间的纸杯、和最远的塑料瓶。

 

 



4.2. 实际场景中三个待测物体的位置

 

 

图4.3是校准后的左右图和匹配出来的disparity map，disparity窗口中是实际的点云，object窗口是给disparity map加了个阈值之后得到的二值图，主要是为了分割前景和背景。可以看到要测的三个物体基本被正确地分割出来了

 



4.3. 双目摄像头得到的disparity map

 

图4.4是在disparity窗口中选取一个点后然后在实际坐标矩阵中得到的对应三维信息，在这里，我在三个物体的点云上各选一个点来代表一个物体实际的坐标信息。(这里通过鼠标获取一点坐标信息的方法参考的是opencv sample里的watershed.cpp)

 

 



4.4. 对应点的三维坐标

 

 

在这里可以看到，(265, 156)也就是利乐包装盒的坐标是(13, 12, -157)，(137, 142)纸杯的坐标是(77, 30, -312)，(95, 115)塑料瓶的坐标是(144, 63, -482)。

补充一下：为了方便显示，所以视差图出来之后进行了一个0-255的normalize，所以value值的前一个是normalize之后点的灰度值，后一个是normalize之前点的实际视差图。

由cvFindStereoCorrespondenceBM算法的源代码：

 

dptr[y*dstep] = (short)(((ndisp - mind - 1 + mindisp)*256 + (d != 0 ? (p-n)*128/d : 0) + 15) >> 4);

其中

ndisp是ndisp = state->numberOfDisparities;

mindisp是mindisp = state->minDisparity;

mind就是sad得出的视差

实际视差大约是(64-mind-1)*256=1163, 基本是对的, 后面一项修正值在通常情况下可以忽略

 

 

 

目前我还是不是很清楚立体坐标系原点(左相机坐标原点)和尺度，但是从这三个点的z坐标可以大致看出这三个物体的距离差大概是1:2:3，基本与实际场景中物体的位置一致。因此，可以通过这种方法确定出物体的大致距离信息。

 

但是，如果就从摄像头参数本身来测量距离的话，就不是很明白了，还求这方面的大牛解答。

 

5.  一些问题

5.1 关于StereoCalibrate

OpenCV自带的cvStereoCalibrate感觉不怎么好用，用这个函数求出的内参外参和旋转平移矩阵进行校准，往往无法达到行对准，有时甚至会出现比较可怕的畸变。在看了piao的http://www.opencv.org.cn/forum/viewtopic.php?f=1&t=4603帖子之后，也曾经尝试过现用cvCalibrateCamera2单独标定(左右各20幅图)，得出的结果基本和Matlab单独标定的相同，然后再在cvStereoCalibrate中将参数设成CV_CALIB_USE_INTRINSIC_GUESS，用来细化内参数和畸变参数，结果得出的标定结果就又走样了。

不知道有谁在这方面有过成功经验的，可以出来分享一下。毕竟用Matlab工具箱还是麻烦了些。

 

 

 

5.2 Translation向量以及立体匹配得出的世界坐标系

Learning OpenCV中对于Translation和Rotation的图示是这样的



5.1. Learning OpenCV中的图示

 

可是在实验过程中发现，如果将Translation向量按尺度缩放，对于StereoRectify之后的左右视图不会有变化，比如将T = [ -226.73817   -0.62302  8.93984 ] ，变成T = [ -22.673817   -0.062302  0.893984 ]，在OpenCV中显示的结果不会有任何变化。而且我如果修改其中的一个参量的话，左右视图发生的变化也不是图5.1中所示的那种变化(比如把x缩小，那么视图发生的变化不是往x轴方向的平移)。

 

因此又回到了老问题，这里这些坐标的尺度究竟是什么？通过ReprojectTo3D那个函数得到的三维坐标又是以哪个点为原点，那三个方向为x,y,z轴的？ 

 

补充: 对这个问题的解答来自于和maxwellsdemon的讨论

他的解释如下:rotation是两者的旋转角度的关系，但是你要把它矫正平行，也是需要translation matrix的。你可以设想，两个看似已经平行了的摄像头，但是深度上放置的有差距，那么在矫正的时候会议translation matrix所对应的角度或者直线为基准，二者旋转一个小角度，使得完全平行。

 

 

 

 

 

目前立体匹配算法是计算机视觉中的一个难点和热点，算法很多，但是一般的步骤是：

 

 

A、匹配代价计算

匹配代价计算是整个立体匹配算法的基础，实际是对不同视差下进行灰度相似性测量。常见的方法有灰度差的平方SD（squared intensity differences），灰度差的绝对值AD（absolute intensity differences）等。另外，在求原始匹配代价时可以设定一个上限值，来减弱叠加过程中的误匹配的影响。以AD法求匹配代价为例，可用下式进行计算，其中T为设定的阈值。



 

图18

 

 

B、 匹配代价叠加

一般来说，全局算法基于原始匹配代价进行后续算法计算。而区域算法则需要通过窗口叠加来增强匹配代价的可靠性，根据原始匹配代价不同，可分为：



 

图19

 

 

C、 视差获取

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

 

 

D、视差细化（亚像素级）

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

 

有关立体匹配的介绍和常见匹配算法的比较，推荐大家看看Stefano Mattoccia 的讲义 Stereo Vision: algorithms and applications，190页的ppt，讲解得非常形象详尽。

 

 

 

 

 

1． opencv2.1和opencv2.0在做stereo vision方面有什么区别了？

2.1版增强了Stereo Vision方面的功能：

(1) 新增了 SGBM 立体匹配算法（源自Heiko Hirschmuller的《Stereo Processing by Semi-global Matching and Mutual Information》），可以获得比 BM 算法物体轮廓更清晰的视差图（但低纹理区域容易出现横/斜纹路，在 GCstate->fullDP 选项使能时可消减这种异常纹路，但对应区域视差变为0，且运行速度会有所下降），速度比 BM 稍慢， 352*288的帧处理速度大约是 5 帧/秒；

(2) 视差效果：BM < SGBM < GC；处理速度：BM > SGBM > GC ；

(3) BM 算法比2.0版性能有所提升，其状态参数新增了对左右视图感兴趣区域 ROI 的支持（roi1 和 roi2，由stereoRectify函数产生）；

(4) BM 算法和 GC 算法的核心代码改动不大，主要是面向多线程运算方面的（由 OpenMP 转向 Intel TBB）；

(5) cvFindStereoCorrespondenceBM 函数的disparity参数的数据格式新增了 CV_32F 的支持，这种格式的数据给出实际视差，而 2.0 版只支持 CV_16S，需要除以 16.0 才能得到实际的视差数值。

 

 

2． 用于立体匹配的图像可以是彩色的吗？

在OpenCV2.1中，BM和GC算法只能对8位灰度图像计算视差，SGBM算法则可以处理24位（8bits*3）彩色图像。所以在读入图像时，应该根据采用的算法来处理图像：

 

 

int color_mode = alg == STEREO_SGBM ? 1 : 0;

//

// 载入图像

cvGrabFrame( lfCam );

cvGrabFrame( riCam );

frame1 = cvRetrieveFrame( lfCam );

frame2 = cvRetrieveFrame( riCam );

if(frame1.empty()) break;

resize(frame1, img1, img_size, 0, 0);

resize(frame2, img2, img_size, 0, 0);

// 选择彩色或灰度格式作为双目匹配的处理图像

if (!color_mode && cn>1)

{

cvtColor(img1, img1gray, CV_BGR2GRAY);

cvtColor(img2, img2gray, CV_BGR2GRAY);

img1p = img1gray;

img2p = img2gray;

}

else

{

img1p = img1;

img2p = img2;

}

 

 

3． 怎样获取与原图像有效像素区域相同的视差图？

在OpenCV2.0及以前的版本中，所获取的视差图总是在左侧和右侧有明显的黑色区域，这些区域没有有效的视差数据。视差图有效像素区域与视差窗口（ndisp，一般取正值且能被16整除）和最小视差值（mindisp，一般取0或负值）相关，视差窗口越大，视差图左侧的黑色区域越大，最小视差值越小，视差图右侧的黑色区域越大。其原因是为了保证参考图像（一般是左视图）的像素点能在目标图像（右视图）中按照设定的视差匹配窗口匹配对应点，OpenCV 只从参考图像的第 (ndisp - 1 + mindisp) 列开始向右计算视差，第 0 列到第 (ndisp - 1 + mindisp) 列的区域视差统一设置为 (mindisp - 1) *16；视差计算到第 width + mindisp 列时停止，余下的右侧区域视差值也统一设置为 (mindisp - 1) *16。  

00177 static const int DISPARITY_SHIFT = 4;
…
00411     int ndisp = state->numberOfDisparities;
00412     int mindisp = state->minDisparity;
00413     int lofs = MAX(ndisp - 1 + mindisp, 0);
00414     int rofs = -MIN(ndisp - 1 + mindisp, 0);
00415     int width = left->cols, height = left->rows;
00416     int width1 = width - rofs - ndisp + 1;
…
00420     short FILTERED = (short)((mindisp - 1) << DISPARITY_SHIFT);
…
00466     // initialize the left and right borders of the disparity map
00467     for( y = 0; y < height; y++ )
00468     {
00469         for( x = 0; x < lofs; x++ )
00470             dptr[y*dstep + x] = FILTERED;
00471         for( x = lofs + width1; x < width; x++ )
00472             dptr[y*dstep + x] = FILTERED;
00473     }
00474     dptr += lofs;
00475
00476     for( x = 0; x < width1; x++, dptr++ )

…


 

这样的设置很明显是不符合实际应用的需求的，它相当于把摄像头的视场范围缩窄了。因此，OpenCV2.1 做了明显的改进，不再要求左右视图和视差图的大小（size）一致，允许对视差图进行左右边界延拓，这样，虽然计算视差时还是按上面的代码思路来处理左右边界，但是视差图的边界得到延拓后，有效视差的范围就能够与对应视图完全对应。具体的实现代码范例如下：

 

//
// 对左右视图的左边进行边界延拓，以获取与原始视图相同大小的有效视差区域
copyMakeBorder(img1r, img1b, 0, 0, m_nMaxDisp, 0, IPL_BORDER_REPLICATE);
copyMakeBorder(img2r, img2b, 0, 0, m_nMaxDisp, 0, IPL_BORDER_REPLICATE);

//
// 计算视差
if( alg == STEREO_BM )
{
	bm(img1b, img2b, dispb);
	// 截取与原始画面对应的视差区域（舍去加宽的部分）
	displf = dispb.colRange(m_nMaxDisp, img1b.cols);	
}
else if(alg == STEREO_SGBM)
{
	sgbm(img1b, img2b, dispb);
	displf = dispb.colRange(m_nMaxDisp, img1b.cols);
}


 

 

 

 

4． cvFindStereoCorrespondenceBM的输出结果好像不是以像素点为单位的视差？

“@scyscyao：在OpenCV2.0中，BM函数得出的结果是以16位符号数的形式的存储的，出于精度需要，所有的视差在输出时都扩大了16倍(2^4)。其具体代码表示如下：

dptr[y*dstep] = (short)(((ndisp - mind - 1 + mindisp)*256 + (d != 0 ? (p-n)*128/d : 0) + 15) >> 4);

可以看到，原始视差在左移8位(256)并且加上一个修正值之后又右移了4位，最终的结果就是左移4位。

因此，在实际求距离时，cvReprojectTo3D出来的X/W,Y/W,Z/W都要乘以16 (也就是W除以16)，才能得到正确的三维坐标信息。”

 

在OpenCV2.1中，BM算法可以用 CV_16S 或者 CV_32F 的方式输出视差数据，使用32位float格式可以得到真实的视差值，而CV_16S 格式得到的视差矩阵则需要 除以16 才能得到正确的视差。另外，OpenCV2.1另外两种立体匹配算法 SGBM 和 GC 只支持 CV_16S 格式的 disparity 矩阵。

 

 

 

5． 如何设置BM、SGBM和GC算法的状态参数？

 

（1）StereoBMState

// 预处理滤波参数

preFilterType：预处理滤波器的类型，主要是用于降低亮度失真（photometric distortions）、消除噪声和增强纹理等, 有两种可选类型：CV_STEREO_BM_NORMALIZED_RESPONSE（归一化响应） 或者 CV_STEREO_BM_XSOBEL（水平方向Sobel算子，默认类型）, 该参数为 int 型；
	
preFilterSize：预处理滤波器窗口大小，容许范围是[5,255]，一般应该在 5x5..21x21 之间，参数必须为奇数值, int 型
	
preFilterCap：预处理滤波器的截断值，预处理的输出值仅保留[-preFilterCap, preFilterCap]范围内的值，参数范围：1 - 31（文档中是31，但代码中是 63）, int
// SAD 参数

SADWindowSize：SAD窗口大小，容许范围是[5,255]，一般应该在 5x5 至 21x21 之间，参数必须是奇数，int 型
	
minDisparity：最小视差，默认值为 0, 可以是负值，int 型
	
numberOfDisparities：视差窗口，即最大视差值与最小视差值之差, 窗口大小必须是 16 的整数倍，int 型
// 后处理参数

textureThreshold：低纹理区域的判断阈值。如果当前SAD窗口内所有邻居像素点的x导数绝对值之和小于指定阈值，则该窗口对应的像素点的视差值为 0（That is, if the sum of absolute values of x-derivatives computed over SADWindowSize by SADWindowSize pixel neighborhood is smaller than the parameter, no disparity is computed at the pixel），该参数不能为负值，int 型
	
uniquenessRatio：视差唯一性百分比， 视差窗口范围内最低代价是次低代价的(1 + uniquenessRatio/100)倍时，最低代价对应的视差值才是该像素点的视差，否则该像素点的视差为 0 （the minimum margin in percents between the best (minimum) cost function value and the second best value to accept the computed disparity, that is, accept the computed disparity d^ only if SAD(d) >= SAD(d^) x (1 + uniquenessRatio/100.) for any d != d*+/-1 within the search range ），该参数不能为负值，一般5-15左右的值比较合适，int 型
	
speckleWindowSize：检查视差连通区域变化度的窗口大小, 值为 0 时取消 speckle 检查，int 型
	
speckleRange：视差变化阈值，当窗口内视差变化大于阈值时，该窗口内的视差清零，int 型
// OpenCV2.1 新增的状态参数

roi1, roi2：左右视图的有效像素区域，一般由双目校正阶段的 cvStereoRectify 函数传递，也可以自行设定。一旦在状态参数中设定了 roi1 和 roi2，OpenCV 会通过cvGetValidDisparityROI 函数计算出视差图的有效区域，在有效区域外的视差值将被清零。
	
disp12MaxDiff：左视差图（直接计算得出）和右视差图（通过cvValidateDisparity计算得出）之间的最大容许差异。超过该阈值的视差值将被清零。该参数默认为 -1，即不执行左右视差检查。int 型。注意在程序调试阶段最好保持该值为 -1，以便查看不同视差窗口生成的视差效果。具体请参见《使用OpenGL动态显示双目视觉三维重构效果示例》一文中的讨论。
在上述参数中，对视差生成效果影响较大的主要参数是 SADWindowSize、numberOfDisparities 和 uniquenessRatio 三个，一般只需对这三个参数进行调整，其余参数按默认设置即可。

在OpenCV2.1中，BM算法有C和C++ 两种实现模块。

 

（2）StereoSGBMState

SGBM算法的状态参数大部分与BM算法的一致，下面只解释不同的部分：

SADWindowSize：SAD窗口大小，容许范围是[1,11]，一般应该在 3x3 至 11x11 之间，参数必须是奇数，int 型
	
P1, P2：控制视差变化平滑性的参数。P1、P2的值越大，视差越平滑。P1是相邻像素点视差增/减 1 时的惩罚系数；P2是相邻像素点视差变化值大于1时的惩罚系数。P2必须大于P1。OpenCV2.1提供的例程 stereo_match.cpp 给出了 P1 和 P2 比较合适的数值。
	
fullDP：布尔值，当设置为 TRUE 时，运行双通道动态编程算法（full-scale 2-pass dynamic programming algorithm），会占用O(W*H*numDisparities)个字节，对于高分辨率图像将占用较大的内存空间。一般设置为 FALSE。
注意OpenCV2.1的SGBM算法是用C++ 语言编写的，没有C实现模块。与H. Hirschmuller提出的原算法相比，主要有如下变化：

算法默认运行单通道DP算法，只用了5个方向，而fullDP使能时则使用8个方向（可能需要占用大量内存）。
	
算法在计算匹配代价函数时，采用块匹配方法而非像素匹配（不过SADWindowSize=1时就等于像素匹配了）。
	
匹配代价的计算采用BT算法（"Depth Discontinuities by Pixel-to-Pixel Stereo" by S. Birchfield and C. Tomasi），并没有实现基于互熵信息的匹配代价计算。
	
增加了一些BM算法中的预处理和后处理程序。
 

（3）StereoGCState

GC算法的状态参数只有两个：numberOfDisparities 和 maxIters ，并且只能通过 cvCreateStereoGCState 在创建算法状态结构体时一次性确定，不能在循环中更新状态信息。GC算法并不是一种实时算法，但可以得到物体轮廓清晰准确的视差图，适用于静态环境物体的深度重构。

注意GC算法只能在C语言模式下运行，并且不能对视差图进行预先的边界延拓，左右视图和左右视差矩阵的大小必须一致。

 

 

6． 如何实现视差图的伪彩色显示？

首先要将16位符号整形的视差矩阵转换为8位无符号整形矩阵，然后按照一定的变换关系进行伪彩色处理。我的实现代码如下： 

 

// 转换为 CV_8U 格式，彩色显示
dispLfcv = displf, dispRicv = dispri, disp8cv = disp8;
if (alg == STEREO_GC)
{
	cvNormalize( &dispLfcv, &disp8cv, 0, 256, CV_MINMAX );
} 
else
{
	displf.convertTo(disp8, CV_8U, 255/(m_nMaxDisp*16.));
}
F_Gray2Color(&disp8cv, vdispRGB);


 

 

灰度图转伪彩色图的代码，主要功能是使灰度图中 亮度越高的像素点，在伪彩色图中对应的点越趋向于 红色；亮度越低，则对应的伪彩色越趋向于 蓝色；总体上按照灰度值高低，由红渐变至蓝，中间色为绿色。其对应关系如下图所示：



 

图20

 

 

void F_Gray2Color(CvMat* gray_mat, CvMat* color_mat)
{
	if(color_mat)
		cvZero(color_mat);
		
	int stype = CV_MAT_TYPE(gray_mat->type), dtype = CV_MAT_TYPE(color_mat->type);
	int rows = gray_mat->rows, cols = gray_mat->cols;

	// 判断输入的灰度图和输出的伪彩色图是否大小相同、格式是否符合要求
	if (CV_ARE_SIZES_EQ(gray_mat, color_mat) && stype == CV_8UC1 && dtype == CV_8UC3)
	{
		CvMat* red = cvCreateMat(gray_mat->rows, gray_mat->cols, CV_8U);
		CvMat* green = cvCreateMat(gray_mat->rows, gray_mat->cols, CV_8U);
		CvMat* blue = cvCreateMat(gray_mat->rows, gray_mat->cols, CV_8U);
		CvMat* mask = cvCreateMat(gray_mat->rows, gray_mat->cols, CV_8U);

		// 计算各彩色通道的像素值
		cvSubRS(gray_mat, cvScalar(255), blue);	// blue(I) = 255 - gray(I)
		cvCopy(gray_mat, red);			// red(I) = gray(I)
		cvCopy(gray_mat, green);			// green(I) = gray(I),if gray(I) < 128
		cvCmpS(green, 128, mask, CV_CMP_GE );	// green(I) = 255 - gray(I), if gray(I) >= 128
		cvSubRS(green, cvScalar(255), green, mask);
		cvConvertScale(green, green, 2.0, 0.0);

		// 合成伪彩色图
		cvMerge(blue, green, red, NULL, color_mat);

		cvReleaseMat( &red );
		cvReleaseMat( &green );
		cvReleaseMat( &blue );
		cvReleaseMat( &mask );
	}
}


 

 

 

7． 如何将视差数据保存为 txt 数据文件以便在 Matlab 中读取分析？

由于OpenCV本身只支持 xml、yml 的数据文件读写功能，并且其xml文件与构建网页数据所用的xml文件格式不一致，在Matlab中无法读取。我们可以通过以下方式将视差数据保存为txt文件，再导入到Matlab中。 

 

void saveDisp(const char* filename, const Mat& mat)		
{
	FILE* fp = fopen(filename, "wt");
	fprintf(fp, "%02d/n", mat.rows);
	fprintf(fp, "%02d/n", mat.cols);
	for(int y = 0; y < mat.rows; y++)
	{
		for(int x = 0; x < mat.cols; x++)
		{
			short disp = mat.at<short>(y, x); // 这里视差矩阵是CV_16S 格式的，故用 short 类型读取
			fprintf(fp, "%d/n", disp); // 若视差矩阵是 CV_32F 格式，则用 float 类型读取
		}
	}
	fclose(fp);
}


 

相应的Matlab代码为：

  

function img = txt2img(filename)
data = importdata(filename);
r = data(1);    % 行数
c = data(2);    % 列数
disp = data(3:end); % 视差
vmin = min(disp);
vmax = max(disp);
disp = reshape(disp, [c,r])'; % 将列向量形式的 disp 重构为 矩阵形式
%  OpenCV 是行扫描存储图像，Matlab 是列扫描存储图像
%  故对 disp 的重新排列是首先变成 c 行 r 列的矩阵，然后再转置回 r 行 c 列
img = uint8( 255 * ( disp - vmin ) / ( vmax - vmin ) );
mesh(disp);
set(gca,'YDir','reverse');  % 通过 mesh 方式绘图时，需倒置 Y 轴方向
axis tight; % 使坐标轴显示范围与数据范围相贴合，去除空白显示区


 

显示效果如下：



 

图21