单目照相机将三维的事物转化为二维的图像,丢失了距离信息,也就是视深。和我们双眼获知物体距离的原理类似,(对极几何)可以参见下图:
记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,并加以实现,调试运行。
- Opencv3.4提供了BM和半全局的SGBM方法。只用调用指针调节参数就可得到视差图,参数设置参考相关博客。
- 全局有动态规划等方法暂未作展开。
- 利用opencv3.4附带的相机标定样例函数,传入自己的棋盘图像,得出相机参数矩阵。
- 继续看有关计算机视觉相关视频。
本周目标:
- 利用SBGM(较高效)得到的视差图和相机参数矩阵获得稀疏图的三维重建。
- 了解orb-slam
后期补充:部分视差图效果并不理想且用时较长。
附原始左右图: