立体视觉

单目照相机将三维的事物转化为二维的图像,丢失了距离信息,也就是视深。和我们双眼获知物体距离的原理类似,(对极几何)可以参见下图:

记cc'距离为B,xx'距离为视差d,空间点到基线的距离为Z,相机的焦距为f,则根据相似三角形定理得:d=Bf/Z

也就是距离信息可以转化为根据不同方向的图片求视差d的问题。

局部的方法(计算匹配代价):

  • SAD(Sum of Absolute Differences):绝对差值的和
  • SSD(Sum of Squared Differences):差值平方的和
  • NCC(Normalized Cross Correlation):图像的相关性

原理简述:其实视差的求取关键在于找到对应的x与x',则问题转换为如何找到左图中某点在右图中的对应点?以图像点x为例,在左图选取一个以x点为中心的特定大小的矩阵块,然后在右图中遍历,选取同样大小的矩阵块使得匹配代价最小,这时矩阵的中心点也就是右图中对应的x'点所在的位置。匹配代价的计算是针对矩阵中对应元素分别操作的,每个元素间的误差有两种定义方式,一种是用绝对差值,一种是用差值的平方,然后要将矩阵中所有元素误差加和,这样寻找最佳匹配点的方法分别叫做SAD和SSD。因为是对矩阵中对应元素操作,既可以按行列遍历操作求解(即下面每种方法给出的第一种代码形式),也可以直接借助矩阵操作函数像absdiff、sum、multiply等进行整体运算。(即每种方法的第二种代码形式)。

NCC计算图像匹配区域的互相关性,与上面两种的计算方法不同,公式如下:

其中Wp为以像素p=(px,py)为中心的匹配窗口。I1为左图匹配窗口内的像素,I2为右图匹配窗口内的像素。注意这里是要找相关性最大的匹配区域。

全局的方法:动态规划、图割法等。

OPENCV3.4.4中提供的方法:

  • BM(Block Matching)
  • SGBM(Semi-Global Block Matching)

原理简述:BM块匹配,可以看作封装好的SAD方法,原理同上。

SGBM因为是半全局算法,涉及到动态规划,效果比BM要好。它的原理较为复杂,具体可参考博客:

https://blog.csdn.net/A_L_A_N/article/details/81490043

这里对全局方法未作展开,后面的文章主要第一借助opencv完成相机的标定工作,第二结合相机参数和视差图,完成视深图和稀疏图的三维重建。

以下代码在VS2017+OPENCV3.4.4的环境下运行成功,参考博客:

https://blog.csdn.net/liyingjiang22/article/details/53156331

https://www.cnblogs.com/riddick/p/8318997.html 

https://blog.csdn.net/cxgincsu/article/details/74451940

SGBM效果比较好,所以后期将采用该种方法继续研究。

下面为可执行代码:

SAD

Mat SAD(Mat left, Mat right, bool turn, int size) {
	int i, j, d, count, x, y, size_s, k, index;
	double temp, temp1, temp2, min = 1000000000;
	int row = left.rows;
	int col = left.cols;
	Mat output(row, col, CV_8UC1);
	vector<pair<double, int>> v;
	size_s = size / 2;
	for (i = 0; i < row; i++) {
		for (j = 0; j < col; j++) {
			v.clear();
			for (d = 0; d <= 79; d++) {//在-79~79范围内搜索
				if (j - d < 0 && turn)break;
				else if (j + d >= col && (!turn))break;
				temp = 0;
				count = 0;//自适应patch,计算窗口像素个数
				if (!turn)d = -d;//用于复用循环
				for (x = i - size_s; x <= i + size_s; x++) {
					for (y = j - d - size_s; y <= j - d + size_s; y++) {
						if (x < 0 || y < 0 || y + d >= col || y + d < 0 || x >= row || y >= col)continue;
						if (turn) {
							temp1 = static_cast<double>(right.ptr<uchar>(x)[y]);
							temp2 = static_cast<double>(left.ptr<uchar>(x)[y + d]);
						}
						else {
							temp1 = static_cast<double>(right.ptr<uchar>(x)[y + d]);
							temp2 = static_cast<double>(right.ptr<uchar>(x)[y]);
						}
						count++;
						temp += temp1 > temp2 ? (temp1 - temp2) : (temp2 - temp1);
					}
				}
				if (!turn)d = -d;
				v.push_back(make_pair(temp / count, d));
			}
			index = 0;
			for (k = 0; k < v.size(); k++) {
				if (v[k].first < min) {
					min = v[k].first;
					index = v[k].second;
				}
			}
			output.ptr<uchar>(i)[j] = static_cast<uchar>(index * 3);
		}
	}
	return output;
}

利用矩阵并行计算,并将方法用类封装。 

#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++)  //左图从DSR开始遍历
	{
		for (int j = 0; j < Height - winSize; j++)
		{
			Kernel_L = L(Rect(i, j, winSize, winSize));
			Mat MM(1, DSR, CV_32F, Scalar(0)); //
			for (int k = 0; k < DSR; k++)
			{
				int x = 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 = static_cast<float>(ADD[0]);
					MM.at<float>(k) = a;
				}
			}
			Point minLoc;
			minMaxLoc(MM, NULL, NULL, &minLoc, NULL);
			int loc = minLoc.x;
			//int loc=DSR-loc;
			Disparity.at<char>(j, i) = loc * 16;
		}
		double rate = double(i) / (Width);
		cout << "SAD已完成" << setprecision(2) << rate * 100 << "%" << endl; //处理进度
	}
	return Disparity;
}

SSD

Mat SSD(Mat left, Mat right, bool turn, int size) {
	int i, j, d, count, x, y, size_s, k, index;
	double temp, temp1, temp2, min=1000000000;
	int row = left.rows;
	int col = left.cols;
	Mat output(row, col, CV_8UC1);
	vector<pair<double, int>> v;
	size_s = size / 2;
	for (i = 0; i < row; i++) {
		for (j = 0; j < col; j++) {
			v.clear();
			for (d = 0; d <= 79; d++) {//在-79~79范围内搜索
				if (j - d < 0 && turn)break;
				else if (j + d >= col && (!turn))break;
				temp = 0;
				count = 0;//自适应patch,计算窗口像素个数
				if (!turn)d = -d;//用于复用循环
				for (x = i - size_s; x <= i + size_s; x++) {
					for (y = j - d - size_s; y <= j - d + size_s; y++) {
						if (x < 0 || y < 0 || y + d >= col || y + d < 0 || x >= row || y >= col)continue;
						if (turn) {
							temp1 = static_cast<double>(right.ptr<uchar>(x)[y]);
							temp2 = static_cast<double>(left.ptr<uchar>(x)[y + d]);
						}
						else {
							temp1 = static_cast<double>(right.ptr<uchar>(x)[y+d]);
							temp2 = static_cast<double>(right.ptr<uchar>(x)[y]);
						}
						count++;
						temp += pow(temp1 - temp2, 2);
					}
				}
				if (!turn)d = -d;
				v.push_back(make_pair(temp / count, d));
			}
			index = 0;
			for (k = 0; k < v.size(); k++) {
				if (v[k].first < min) {
					min = v[k].first;
					index = v[k].second;
				}
			}
			output.ptr<uchar>(i)[j] = static_cast<uchar>(index * 3);
		}
	}
	return output;
}

利用矩阵并行计算,并将方法用类封装。

#include"iostream"
#include"opencv2/opencv.hpp"
#include"iomanip"
using namespace std;
using namespace cv;

class SSD
{
public:
	SSD() :winSize(7), DSR(30) {}
	SSD(int _winSize, int _DSR) :winSize(_winSize), DSR(_DSR) {}
	Mat computerSSD(Mat &L, Mat &R); //计算SAD
private:
	int winSize; //卷积核的尺寸
	int DSR;     //视差搜索范围
};

Mat SSD::computerSSD(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++)  //左图从DSR开始遍历
	{
		for (int j = 0; j < Height - winSize; j++)
		{
			Kernel_L = L(Rect(i, j, winSize, winSize));
			Mat MM(1, DSR, CV_32F, Scalar(0)); //
			for (int k = 0; k < DSR; k++)
			{
				int x = i - k;
				if (x >= 0)
				{
					Kernel_R = R(Rect(x, j, winSize, winSize));
					Mat Dif;
					absdiff(Kernel_L, Kernel_R, Dif);//
					Dif = Dif.mul(Dif);
					Scalar ADD = sum(Dif);
					float a = static_cast<float>(ADD[0]);
					MM.at<float>(k) = a;
				}
			}
			Point minLoc;
			minMaxLoc(MM, NULL, NULL, &minLoc, NULL);
			int loc = minLoc.x;
			//int loc=DSR-loc;
			Disparity.at<char>(j, i) = loc * 16;
		}
		double rate = double(i) / (Width);
		cout << "SSD已完成" << setprecision(2) << rate * 100 << "%" << endl; //处理进度
	}
	return Disparity;
}

Tip: A.mul(B)是矩阵A与矩阵B对应位相乘得到的新矩阵,A*B则对应着矩阵乘法,二者不是同一概念。

  NCC

Mat NCC(Mat left, Mat right, bool turn, int size) {
	int i, j, d, x, y, size_s, k, index;
	double tempUp, tempDown, temp1, temp2, result, max = 0;
	int row = left.rows;
	int col = left.cols;
	Mat output(row, col, CV_8UC1);
	vector<pair<double, int>> v;
	size_s = size / 2;
	for (i = 0; i < row; i++) {
		for (j = 0; j < col; j++) {
			v.clear();
			for (d = 0; d <= 79; d++) {//在-79~79范围内搜索
				if (j - d < 0 && turn)break;
				else if (j + d >= col && (!turn))break;
				tempUp = 0;
				tempDown = 0;
				if (!turn)d = -d;//用于复用循环
				for (x = i - size_s; x <= i + size_s; x++) {
					for (y = j - d - size_s; y <= j - d + size_s; y++) {
						if (x < 0 || y < 0 || y + d >= col || y + d < 0 || x >= row || y >= col)continue;
						if (turn) {
							temp1 = static_cast<double>(right.ptr<uchar>(x)[y]);
							temp2 = static_cast<double>(left.ptr<uchar>(x)[y + d]);
						}
						else {
							temp1 = static_cast<double>(right.ptr<uchar>(x)[y + d]);
							temp2 = static_cast<double>(right.ptr<uchar>(x)[y]);
						}
						tempUp += temp1 * temp2;
					}
				}
				tempUp = 0;
				tempDown = 0;
				for (x = i - size_s; x <= i + size_s; x++) {
					for (y = j - d - size_s; y <= j - d + size_s; y++) {
						if (x < 0 || y < 0 || y + d >= col || y + d < 0 || x >= row || y >= col)continue;
						if (turn) {
							temp1 = pow(static_cast<double>(right.ptr<uchar>(x)[y]), 2);
							temp2 = pow(static_cast<double>(left.ptr<uchar>(x)[y + d]), 2);
						}
						else {
							temp1 = pow(static_cast<double>(right.ptr<uchar>(x)[y + d]), 2);
							temp2 = pow(static_cast<double>(right.ptr<uchar>(x)[y]), 2);
						}
					}
				}
				tempUp += sqrt(temp1 * temp2);
				result = tempUp / tempDown;
				if (!turn)d = -d;
				v.push_back(make_pair(result, d));
			}
			index = 0;
			for (k = 0; k < v.size(); k++) {
				if (v[k].first >max) {
					max = v[k].first;
					index = v[k].second;
				}
			}
			output.ptr<uchar>(i)[j] = static_cast<uchar>(index * 3);
		}
	}
	return output;
}

利用矩阵并行计算,并将方法用类封装。

#include"iostream"
#include"opencv2/opencv.hpp"
#include"iomanip"
using namespace std;
using namespace cv;

class NCC
{
public:
	NCC() :winSize(7), DSR(30) {}
	NCC(int _winSize, int _DSR) :winSize(_winSize), DSR(_DSR) {}
	Mat computerNCC(Mat &L, Mat &R); //计算SAD
private:
	int winSize; //卷积核的尺寸
	int DSR;     //视差搜索范围
};

Mat NCC::computerNCC(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++)  //左图从DSR开始遍历
	{
		for (int j = 0; j < Height - winSize; j++)
		{
			Kernel_L = L(Rect(i, j, winSize, winSize));
			Mat MM(1, DSR, CV_32F, Scalar(0)); //
			for (int k = 0; k < DSR; k++)
			{
				int x = i - k;
				if (x >= 0)
				{
					Kernel_R = R(Rect(x, j, winSize, winSize));
					Mat Dif = Kernel_L.mul(Kernel_R);
					Scalar ADD = sum(Dif);
					float a = static_cast<float>(ADD[0]);
					MM.at<float>(k) = a;
				}
			}
			for (int k = 0; k < DSR; k++)
			{
				int x = i - k;
				if (x >= 0)
				{
					Kernel_R = R(Rect(x, j, winSize, winSize));
					Scalar ADD2 = sum(Kernel_L)*sum(Kernel_R);
					float b = static_cast<float>(ADD2[0]);
					MM.at<float>(k) /=sqrt(b);
				}
			}
			Point minLoc;
			minMaxLoc(MM, NULL, NULL, &minLoc, NULL);
			int loc = minLoc.x;
			//int loc=DSR-loc;
			Disparity.at<char>(j, i) = loc * 16;
		}
		double rate = double(i) / (Width);
		cout << "NCC已完成" << setprecision(2) << rate * 100 << "%" << endl; //处理进度
	}
	return Disparity;
}

OPENCV3.4.4中,StereoBM不能直接声明实例,需要放到智能指针里cv::Ptr<cv::StereoBM>声明才行,且必须通过setter和

getter方法来设置和获取参数。StereoSGBM同理。(这里的参数要根据实际情况调整,注意图像格式问题)

使用样例如下: 

BM

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

int main() {
	//读入灰度图left和right
	Mat left = imread("D:/project_test/stereo_matching/images/im1.png",0);
	Mat right = imread("D:/project_test/stereo_matching/images/im2.png",0);
	if (left.empty()||right.empty()) {
		printf("could not found image!");
		return -1;
	}
	imshow("left", left);
	imshow("right", right);
	//存储视差
	Mat disp(7, 30, CV_8U, Scalar(0)); //视差图
	//bm
	int mindisparity = 0;//最小视差,默认为0。此参数决定左图中的像素点在右图匹配搜索的起点
	int numdisparities = 64;//视差搜索范围长度,其值为16的整数倍。最大视差= minDisparity + numDisparities -1;
	int blocksize =15;//SAD代价计算窗口大小,默认为5。窗口大小为奇数,一般在3*3 到21*21之间
	//智能指针,利用setter和getter对bm操作
	cv::Ptr<cv::StereoBM> bm = cv::StereoBM::create(numdisparities, blocksize);
	//setter设置参数
	bm->setBlockSize(blocksize);
	bm->setMinDisparity(mindisparity);
	bm->setNumDisparities(numdisparities);

	bm->setPreFilterSize(9);
	bm->setPreFilterCap(31);

	bm->setTextureThreshold(10); 
	bm->setUniquenessRatio(15); 
	bm->setSpeckleWindowSize(100);
	bm->setSpeckleRange(32);
	bm->setDisp12MaxDiff(1); 
	
	copyMakeBorder(left, left, 0, 0, 80, 0, IPL_BORDER_REPLICATE);  //防止黑边
	copyMakeBorder(right, right, 0, 0, 80, 0, IPL_BORDER_REPLICATE);
	bm->compute(left, right, disp);
	/*U<S<F*/
	//bm计算出的视差都是CV_16S格式的,使用32位float格式可以得到真实的视差值,所以我们需要除以16
	disp.convertTo(disp, CV_32F, 1.0 / 16); //除以16得到真实视差值
	disp = disp.colRange(80, disp.cols);

	//归一化时,normType应该为cv::NORM_MINMAX,alpha为归一化后的最大值,beta为归一化后的最小值
	Mat disp8U = Mat(disp.rows, disp.cols, CV_8UC1);
	normalize(disp, disp8U, 0, 255, NORM_MINMAX, CV_8UC1);

	imshow("disp", disp8U);
	waitKey(0);
	return 0;
}

minDisparity:最小视差,默认为0。
numDisparities:视差搜索范围长度,其值必须为16的整数倍。
blockSize:SAD代价计算窗口大小,默认为5。窗口大小为奇数,一般在3*3 到21*21之间;
P1、P2:能量函数参数,P1是相邻像素点视差增/减 1 时的惩罚系数;P2是相邻像素点视差变化值大于1时的惩罚系数。

SGBM

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

int main()
{
	Mat left = imread("D:/project_test/stereo_matching/images/im1.png", 0);
	Mat right = imread("D:/project_test/stereo_matching/images/im2.png", 0);
	if (!left.data || !right.data) {
		cout << "cannot found images!\n";
	}
	Mat disp(7, 30, CV_8U, Scalar(0));
	int mindisparity = 0;
	int ndisparities = 64;
	int SADWindowSize = 15;
	//SGBM
	Ptr<StereoSGBM> sgbm = StereoSGBM::create(mindisparity, ndisparities, SADWindowSize);
	int P1 = 8 * left.channels() * SADWindowSize* SADWindowSize;
	int P2 = 32 * left.channels() * SADWindowSize* SADWindowSize;
	sgbm->setP1(P1);
	sgbm->setP2(P2);
	sgbm->setPreFilterCap(31);
	sgbm->setUniquenessRatio(15);
	sgbm->setSpeckleRange(32);
	sgbm->setSpeckleWindowSize(100);
	sgbm->setDisp12MaxDiff(1);
	sgbm->compute(left, right, disp);
	copyMakeBorder(left, left, 0, 0, 80, 0, IPL_BORDER_REPLICATE);  //防止黑边
	copyMakeBorder(right, right, 0, 0, 80, 0, IPL_BORDER_REPLICATE);
	disp.convertTo(disp, CV_32F, 1.0 / 16);                //除以16得到真实视差值
	disp = disp.colRange(80, disp.cols);
	Mat disp8U = Mat(disp.rows, disp.cols, CV_8UC1);       //显示
	normalize(disp, disp8U, 0, 255, NORM_MINMAX, CV_8UC1);
	imshow("SBGM_disp",disp8U);
	waitKey(0);
	return 0;
}

这里给出矩阵方法SAD、SSD、NCC和BM、SGBM的运行结果。

得到的视差图如下:

SAD
SSD
NCC
BM
SGBM

上周总结:

  1. 由对极几何基本原理,了解常见的局部立体匹配算法SAD、SSD、NCC,并加以实现,调试运行。 ​​
  2. Opencv3.4提供了BM和半全局的SGBM方法。只用调用指针调节参数就可得到视差图,参数设置参考相关博客。
  3. 全局有动态规划等方法暂未作展开。
  4. 利用opencv3.4附带的相机标定样例函数,传入自己的棋盘图像,得出相机参数矩阵。
  5. 继续看有关计算机视觉相关视频。

本周目标:

  1. 利用SBGM(较高效)得到的视差图和相机参数矩阵获得稀疏图的三维重建。
  2. 了解orb-slam

 后期补充:部分视差图效果并不理想且用时较长。

附原始左右图:

左图
右图

 

 

 

 

 

 

  • 3
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值