深度学习 立体匹配 关系体计算 图像立体匹配

一 概念与步骤

1.概念： SAD(Sum of absolute differences)是一种图像匹配算法。基本思想：差的绝对值之和。此算法常用于图像块匹配，将每个像素对应数值之差的绝对值求和，据此评估两个图像块的相似度。该算法快速、但并不精确，通常用于多级处理的初步筛选。

2.步骤：（1）构造一个小窗口,类似于卷积核；
（2）用窗口覆盖左边的图像,选择出窗口覆盖区域内的所有像素点；
（3）同样用窗口覆盖右边的图像并选择出覆盖区域的像素点；
（4）左边覆盖区域减去右边覆盖区域，并求出所有像素点灰度差的绝对值之和；
（5）移动右边图像的窗口，重复（3）-（4）的处理（这里有个搜索范围,超过这个范围跳出）；
（6）找到这个范围内SAD值最小的窗口,即找到了左图锚点的最佳匹配的像素块。

3.说明：用SAD算法可以得出左右图像的视差，若再一步处理就可以得到深度图，深度与视差存在反比的关系，我们可以简单体会一下：将手指头放在离眼睛不同距离的位置，并轮换睁、闭左右眼，可以发现手指在不同距离的位置，视觉差也不同，且距离越近，视差越大，其中距离的远近就是深度了。而且可以注意到，当你用左眼看手指时，手指在你眼中的靠右位置，而用右眼看时，手指在你眼中靠左的位置。两只眼分别看到的图像是一样大的，如果用（x,y）表示左眼视图中某个位置的坐标，那么相应的该位置右眼视图的坐标应该为（x-d,y），其中d为视差。这时(x,y)和(x-d,y)就是最佳匹配点。但是实际情况我们并不知道d是多少。SAD算法就给出了如何求d.
SAD算法的思想:只需我们按视差搜索范围从0开始搜索，找到左右图像最匹配的点，对应的视差值就确定了。那么怎么确定最佳匹配点呢？试想一下，如果视差为0，也就是左右图像一样，那么这个点上下左右区域对应的点都应该相同，所以像素相减后都为0，由于视差的存在（简单理解为从不同的角度看物体，由于光照的影响像素值也会发生改变），该点上下左右区域的像素值不会完全相等，但是我们依然可以利用这个思想，设定一个小窗口，在左右两幅图中计算其像素值差的绝对值之和。假如视差搜索范围为0-50，那么就会得到51个结果。若在某个视差值d下该绝对值之和最小，那么d就为该中心点对应的视差。再由视差与深度的关系就可以得到深度图。

二.代码实现

关于SAD算法的opencv实现代码有很多，这里给出两个：

#include"iostream"
#include"opencv2/opencv.hpp"
#include"iomanip"
using namespace std;
using namespace cv;
 
class SAD
{
public:
	SAD() :winSize(7), DSR(30) {}
	SAD(int _winSize, int _DSR) :winSize(_winSize), DSR(_DSR) {}
	Mat computerSAD(Mat &L, Mat &R); //计算SAD
private:
	int winSize; //卷积核的尺寸
	int DSR;     //视差搜索范围
 
};
 
Mat SAD::computerSAD(Mat &L, Mat &R)
{
	int Height = L.rows;
	int Width = L.cols;
	Mat Kernel_L(Size(winSize, winSize), CV_8U, Scalar::all(0));
	Mat Kernel_R(Size(winSize, winSize), CV_8U, Scalar::all(0));
	Mat Disparity(Height, Width, CV_8U, Scalar(0)); //视差图
 
	for (int i = 0; i<Width - winSize; i++)  
	{
		for (int j = 0; j<Height - winSize; j++)
		{
			Kernel_L = L(Rect(i, j, winSize, winSize));
			Mat MM(1, DSR, CV_32F, Scalar(0)); //MM是一个1行DSR列的图像（矩阵）
 
			for (int k = 0; k<DSR; k++)
			{
				int x = i - k; //为什么是i-k参见我上面的叙述
				if (x >= 0)
				{
					Kernel_R = R(Rect(x, j, winSize, winSize));
					Mat Dif;
					absdiff(Kernel_L, Kernel_R, Dif);//
					Scalar ADD = sum(Dif);
					float a = ADD[0];//a为视差为k是相应窗口的像素差值的绝对值之和
					MM.at<float>(k) = a;//将a赋给MM的第k列，因为从0开始搜索，遍历结束后MM每一列为视差为列序号时对应的SAD值，我们取其最小即可
				}
 
			}
			Point minLoc; //point数据类型为二维点对象，有横纵xy两个坐标
			minMaxLoc(MM, NULL, NULL, &minLoc, NULL);//返回MM最小值的坐标
 
			int loc = minLoc.x;//取最小值坐标的横坐标x值，即为对应的列序号，也就是相应的视差值
			//int loc=DSR-loc;
			Disparity.at<char>(i, j) = loc * 16;//*16只是为了方便显示
 
		}
		double rate = double(i) / (Width);
		//cout << "已完成" << setprecision(2) << rate * 100 << "%" << endl; //处理进度
	}
	return Disparity;
}

我在其中做了一些注释
这个代码比较容易理解，更能体现SAD算法的过程和思想。最后的loc*16只是简单的把视差乘以一个倍数方便显示，所得的图是视差图而不是深度图，若要得到深度图还需进一步处理。
这段代码应该是没有问题的，但是我运行不出来，一直出现abort() has been called的错误，也许是opencv版本或者指针的什么问题，但是无伤大雅，理解更为重要。

接下来这个可以运行，就是代码比较粗暴，不是很容易理解。

#include <opencv2/opencv.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <iostream>
#include <ctime>
using namespace std;
using namespace cv;

const int n = 7; //窗口大小：2*n+1
const int range = 150;//视差范围

void SAD(uchar* Limg,
	uchar* Rimg,
	uchar* Oimg,int w,int h)
{
	for (int y = 0; y< h; y++)
	{
		for (int x = 0; x< w; x++){
			unsigned int bestCost = 999999;
			unsigned int bestDisparity = 0;
			for (int d = 0; d <= range; d++)
			{
				unsigned int cost = 0;
				for (int i = -n; i <= n; i++)
				{
					for (int j = -n; j <= n; j++)
					{
						int yy, xx, xxd;
						yy = y + i;
						if (yy < 0) yy = 0;
						if (yy >= h) yy = h-1;

						xx = x + j;
						if (xx < 0) xx = 0;
						if (xx >= w) xx = w-1;

						xxd = xx - d;
						if (xxd < 0) xxd = 0;
						if (xxd >= w) xxd = w-1;
						cost += abs((int)(Limg[yy*w + xx] - Rimg[yy*w + xxd]));
					}
				}
				if (cost < bestCost)
				{
					bestCost = cost;
					bestDisparity = d;
				}
				Oimg[y*w + x] = bestDisparity*4;
			}
		}
	}
}

int main()
{
	clock_t starttime = clock();
	Mat imL, imR, imO;
	imL = imread("im2.png", 0);
	if (imL.empty())
	{
		return -1;
	}
	imR = imread("im6.png", 0);
	if (imR.empty())
	{
		return -1;
	}
	imO.create(imL.rows, imL.cols, CV_8UC1);

	SAD(imL.data, imR.data, imO.data, imL.cols, imL.rows);
	namedWindow("left", WINDOW_AUTOSIZE);
	namedWindow("right", WINDOW_AUTOSIZE);
	namedWindow("Output", WINDOW_AUTOSIZE);
	imwrite("11.png", imO);
	imshow("Output", imO);
	imshow("left", imL);
	imshow("right", imR);
	clock_t endtime = clock();
	printf("%d\n", (endtime - starttime));
	cout << (endtime - starttime) << endl;
	waitKey(0);
	return 0;
}