SGBM算法获取视差图
立体校正后的左右两幅图像得到后，匹配点是在同一行上的，可以使用OpenCV中的BM算法或者SGBM算法计算视差图。由于SGBM算法的表现要远远优于BM算法，因此采用SGBM算法获取视差图。SGBM中的参数设置如下：
int numberOfDisparities = ((imgSize.width / 8) + 15) & -16;
cv::Ptr sgbm =cv::StereoSGBM::create(0, 16, 3);
sgbm->setPreFilterCap(32);
int SADWindowSize = 9;
int sgbmWinSize = SADWindowSize > 0 ? SADWindowSize : 3;
sgbm->setBlockSize(sgbmWinSize);
intcn = imgL.channels();
sgbm->setP1(8* cn*sgbmWinSize*sgbmWinSize);
sgbm->setP2(32* cn*sgbmWinSize*sgbmWinSize);
sgbm->setMinDisparity(0);
sgbm->setNumDisparities(numberOfDisparities);
sgbm->setUniquenessRatio(10);
sgbm->setSpeckleWindowSize(100);
sgbm->setSpeckleRange(32);
sgbm->setDisp12MaxDiff(1);
intalg = STEREO_SGBM;
if(alg == STEREO_HH)
sgbm->setMode(cv::StereoSGBM::MODE_HH);
else if(alg == STEREO_SGBM)
sgbm->setMode(cv::StereoSGBM::MODE_SGBM);
else if(alg == STEREO_3WAY)
sgbm->setMode(cv::StereoSGBM::MODE_SGBM_3WAY);
sgbm->compute(imgL, imgR, disp);
默认计算出的是左视差图，如果需要计算右视差图，则将上面加粗的三条语句替换为下面前三条语句。由于视差值计算出来为负值，disp类型为16SC1，因此需要取绝对值，然后保存：
sgbm->setMinDisparity(-numberOfDisparities);
sgbm->setNumDisparities(numberOfDisparities);
sgbm->compute(imgR, imgL, disp);
disp = abs(disp);
SGBM算法得到的左、右视差图如下，左视差图的数据类型为CV_16UC1，右视差图的数据类型为CV_16SC1 （SGBM中视差图中不可靠的视差值设置为最小视差（mindisp-1）*16。因此在此例中，左视差图中不可靠视差值设置为-16，截断值为0；右视差图中不可靠视差值设置为（-numberOfDisparities-1）*16，取绝对值后为（numberOfDisparities+1）*16，所以两幅图会有较大差别）：
左视差图（不可靠视差值为0）
右视差图（不可靠视差值为 （numberOfDisparities+1）*16 ）
如果将右视差图不可靠视差值也设置为0，则如下
至此，左视差图和右视差图遥相呼应。
2. 视差图空洞填充
视差图中视差值不可靠的视差大多数是由于遮挡引起，或者光照不均匀引起。既然牛逼如SGBM也觉得不可靠，那与其留着做个空洞，倒不如用附近可靠的视差值填充一下。
空洞填充也有很多方法，在这里我检测出空洞区域，然后用附近可靠视差值的均值进行填充。填充后的视差图
填充后左视差图
填充后右视差图
3. 视差图转换为深度图
视差的单位是像素（pixel），深度的单位往往是毫米（mm）表示。而根据平行双目视觉的几何关系（此处不再画图推导，很简单），可以得到下面的视差与深度的转换公式：
depth = ( f * baseline) / disp
上式中，depth表示深度图；f表示归一化的焦距，也就是内参中的fx； baseline是两个相机光心之间的距离，称作基线距离；disp是视差值。等式后面的均已知，深度值即可算出。
在上面我们用SGBM算法获取了视差图，接下来转换为深度图，函数代码如下：
/*
函数作用：视差图转深度图输入：dispMap ----视差图，8位单通道，CV_8UC1
K ----内参矩阵，float类型输出：depthMap ----深度图，16位无符号单通道，CV_16UC1
*/
voiddisp2Depth(cv::Mat dispMap, cv::Mat&depthMap, cv::Mat K)
{
inttype = dispMap.type();
float fx = K.at<float>(0, 0);
float fy = K.at<float>(1, 1);
float cx = K.at<float>(0, 2);
float cy = K.at<float>(1, 2);
float baseline = 65; //基线距离65mm
if(type == CV_8U)
{
const float PI = 3.14159265358;
intheight = dispMap.rows;
intwidth = dispMap.cols;
uchar* dispData = (uchar*)dispMap.data;
ushort* depthData = (ushort*)depthMap.data;
for (int i = 0; i < height; i++)
{
for (int j = 0; j < width; j++)
{
intid = i*width + j;
if(!dispData[id]) continue; //防止0除
depthData[id] =ushort( (float)fx *baseline / ((float)dispData[id]) );
}
}
}
else
{
cout <<"please confirm dispImg's type!"<< endl;
cv::waitKey(0);
}
}
注：png的图像格式可以保存16位无符号精度，即保存范围为0-65535，如果是mm为单位，则最大能表示约65米的深度，足够了。
上面代码中我设置深度图的精度为CV_16UC1，也就是ushort类型，将baseline设置为65mm，转换后保存为png格式即可。如果保存为jpg或者bmp等图像格式，会将数据截断为0-255。所以保存深度图，png格式是理想的选择。（如果不是为了获取精确的深度图，可以将baseline设置为1，这样获取的是相对深度图，深度值也是相对的深度值）
转换后的深度图如下;
左深度图
右深度图
空洞填充后的深度图，如下：
左深度图（空洞填充后）
右深度图（空洞填充后）
视差图到深度图完成。
注：视差图和深度图中均有计算不正确的点，此文意在介绍整个流程，不特别注重算法的优化，如有大神望不吝赐教。
附：视差图和深度图的空洞填充
步骤如下：
① 以视差图dispImg为例。计算图像的积分图integral，并保存对应积分图中每个积分值处所有累加的像素点个数n（空洞处的像素点不计入n中，因为空洞处像素值为0，对积分值没有任何作用，反而会平滑图像）。
② 采用多层次均值滤波。首先以一个较大的初始窗口去做均值滤波（积分图实现均值滤波就不多做介绍了，可以参考我之前的一篇博客），将大区域的空洞赋值。然后下次滤波时，将窗口尺寸缩小为原来的一半，利用原来的积分图再次滤波，给较小的空洞赋值（覆盖原来的值）；依次类推，直至窗口大小变为3x3，此时停止滤波，得到最终结果。
③ 多层次滤波考虑的是对于初始较大的空洞区域，需要参考更多的邻域值，如果采用较小的滤波窗口，不能够完全填充，而如果全部采用较大的窗口，则图像会被严重平滑。因此根据空洞的大小，不断调整滤波窗口。先用大窗口给所有空洞赋值，然后利用逐渐变成小窗口滤波覆盖原来的值，这样既能保证空洞能被填充上，也能保证图像不会被过度平滑。
空洞填充的函数代码如下，仅供参考：
voidinsertDepth32f(cv::Mat& depth)
{
const intwidth = depth.cols;
const intheight = depth.rows;
float* data = (float*)depth.data;
cv::Mat integralMap = cv::Mat::zeros(height, width, CV_64F);
cv::Mat ptsMap = cv::Mat::zeros(height, width, CV_32S);
double* integral = (double*)integralMap.data;
int* ptsIntegral = (int*)ptsMap.data;
memset(integral,0, sizeof(double) * width * height);
memset(ptsIntegral,0, sizeof(int) * width * height);
for (int i = 0; i < height; ++i)
{
intid1 = i * width;
for (int j = 0; j < width; ++j)
{
intid2 = id1 + j;
if (data[id2] > 1e-3)
{
integral[id2] = data[id2];
ptsIntegral[id2] =1;
}
}
}
// 积分区间
for (int i = 0; i < height; ++i)
{
intid1 = i * width;
for (int j = 1; j < width; ++j)
{
intid2 = id1 + j;
integral[id2] += integral[id2 -1];
ptsIntegral[id2] +=ptsIntegral[id2 -1];
}
}
for (int i = 1; i < height; ++i)
{
intid1 = i * width;
for (int j = 0; j < width; ++j)
{
intid2 = id1 + j;
integral[id2] += integral[id2 -width];
ptsIntegral[id2] +=ptsIntegral[id2 - width];
}
}
intwnd;
double dWnd = 2;
while (dWnd > 1)
{
wnd =int(dWnd);
dWnd /=2;
for (int i = 0; i < height; ++i)
{
intid1 = i * width;
for (int j = 0; j < width; ++j)
{
intid2 = id1 + j;
int left = j - wnd - 1;
intright = j + wnd;
int top = i - wnd - 1;
intbot = i + wnd;
left = max(0, left);
right = min(right,width -1);
top = max(0, top);
bot = min(bot, height -1);
intdx = right - left;
intdy = (bot - top) * width;
intidLeftTop = top * width + left;
intidRightTop = idLeftTop + dx;
intidLeftBot = idLeftTop + dy;
intidRightBot = idLeftBot + dx;
intptsCnt = ptsIntegral[idRightBot] + ptsIntegral[idLeftTop] -(ptsIntegral[idLeftBot] + ptsIntegral[idRightTop]);
doublesumGray = integral[idRightBot] +integral[idLeftTop] - (integral[idLeftBot] + integral[idRightTop]);
if (ptsCnt <= 0)
{
continue;
}
data[id2] =float(sumGray / ptsCnt);
}
}
int s = wnd / 2 * 2 + 1;
if (s > 201)
{
s =201;
}
cv::GaussianBlur(depth, depth,cv::Size(s, s), s, s);
}
}
