立体匹配算法

转载请注明出处：http://www.cnblogs.com/adong7639/p/4267326.html

立体匹配算法最新动态：

http://vision.middlebury.edu/stereo/eval/

http://www.cvlibs.net/datasets/kitti/eval_stereo_flow.php?benchmark=stereo

相关文献：http://blog.csdn.net/xuyuhua1985/article/details/26283389

介绍立体匹配的基本原理: http://vision.deis.unibo.it/~smatt/Seminars/StereoVision.pdf（比较清晰）

立体匹配综述性文章 ： http://wenku.baidu.com/view/5b359d7d5acfa1c7aa00cc7b.html

立体匹配算法的基本目标：找出图像的每个像素点在另一个视角的图像上对应的像素点，算出视差图像，估算出景深图像。

最简单的SAD块匹配算法

//Stereo Match By SAD
#include <opencv2/opencv.hpp>
#include <vector>
#include <algorithm>
#include <iostream>  
#include <windows.h>  
#include <string>  

using namespace std;
using namespace cv;

DWORD t1;  
DWORD t2;  

void timebegin()  
{  
    t1 = GetTickCount();  
}  

void timeend(string str)  
{  
    t2 = GetTickCount();  
    cout << str << " is "<< (t2 - t1)/1000 << "s" << endl;  
}  


float sadvalue(const Mat &src1, const Mat &src2)
{
    Mat  matdiff = cv::abs(src1 -src2);
    int  saddiff = cv::sum(matdiff)[0];
    return saddiff;
}

float GetMinSadIndex(std::vector<float> &sad)
{
    float minsad = sad[0];
    int index = 0;
    int len = sad.size();
    for (int i = 1; i < len; ++i)
    {
        if (sad[i] < minsad)
        {
            minsad = sad[i];
            index = i;
        }
    }
    return index;
}

void MatDataNormal(const Mat &src, Mat &dst)
{
    normalize(src, dst, 255, 0, NORM_MINMAX );
    dst.convertTo(dst, CV_8UC1);
}


void GetPointDepthRight(Mat &disparity, const Mat &leftimg, const Mat  &rightimg, 
    const int MaxDisparity, const  int winsize)
{
    int row = leftimg.rows;
    int col = leftimg.cols;
    if (leftimg.channels() == 3 && rightimg.channels() == 3)
    {
        cvtColor(leftimg, leftimg, CV_BGR2GRAY);
        cvtColor(rightimg, rightimg, CV_BGR2GRAY);
    }

    //Mat disparity = Mat ::zeros(row,col, CV_32S);
    int w = winsize;
    int rowrange = row - w;
    int colrange = col - w - MaxDisparity;

    for (int i = w; i < rowrange; ++i)
    {
        int *ptr = disparity.ptr<int>(i);
        for (int j = w; j < colrange; ++j)
        {
            //Rect rightrect;
            Mat rightwin = rightimg(Range(i - w,i + w + 1),Range(j - w,j + w + 1)); 
            std::vector<float> sad(MaxDisparity);
            for (int d = j; d < j + MaxDisparity; ++d)
            {
                //Rect leftrect;
                Mat leftwin = leftimg(Range(i - w,i + w + 1),Range(d - w,d + w + 1));
                sad[d - j] = sadvalue(leftwin, rightwin);
            }
            *(ptr + j) = GetMinSadIndex(sad);
        }
    }
}

void GetPointDepthLeft(Mat &disparity, const  Mat &leftimg, const Mat  &rightimg, 
    const int MaxDisparity, const  int winsize)
{
    int row = leftimg.rows;
    int col = leftimg.cols;
    if (leftimg.channels() == 3 && rightimg.channels() == 3)
    {
        cvtColor(leftimg, leftimg, CV_BGR2GRAY);
        cvtColor(rightimg, rightimg, CV_BGR2GRAY);
    }

    //Mat disparity = Mat ::zeros(row,col, CV_32S);
    int w = winsize;
    int rowrange = row - w;
    int colrange = col - w;

    for (int i = w; i < rowrange; ++i)
    {
        int *ptr = disparity.ptr<int>(i);
        for (int j = MaxDisparity + w; j < colrange; ++j)
        {
            //Rect leftrect;
            Mat leftwin = leftimg(Range(i - w,i + w + 1),Range(j - w,j + w + 1)); 
            std::vector<float> sad(MaxDisparity);
            for (int d = j; d >  j -  MaxDisparity; --d)
            {
                //Rect rightrect;
                Mat rightwin = rightimg(Range(i - w,i + w + 1),Range(d - w,d + w + 1));
                sad[j - d] = sadvalue(leftwin, rightwin);
            }
            *(ptr + j) = GetMinSadIndex(sad);
        }
    }
}

//(Left-Right Consistency (LRC)
void CrossCheckDiaparity(const Mat &leftdisp, const Mat &rightdisp, Mat &lastdisp, 
    const int MaxDisparity, const int winsize)
{
    int row = leftdisp.rows;
    int col = rightdisp.cols;
    int w = winsize;
    int rowrange = row - w;
    int colrange = col - MaxDisparity - w;
    int diffthreshold = 2;
    for (int i = w; i < row -w; ++i)
    {
        const int *ptrleft = leftdisp.ptr<int>(i);
        const int *ptrright = rightdisp.ptr<int>(i);
        int *ptrdisp = lastdisp.ptr<int>(i);
        for (int j = MaxDisparity + w; j < col - MaxDisparity - w; ++j)
        {
            int leftvalue = *(ptrleft + j);
            int rightvalue = *(ptrright + j - leftvalue );
            int diff = abs(leftvalue - rightvalue);
            if (diff > diffthreshold)
            {
                *(ptrdisp + j) = 0;
            }else
            {
                *(ptrdisp + j) = leftvalue;
            }
        }
    }

}


int main()
{
    Mat leftimg = imread("left1.png",0);   
    Mat rightimg = imread("right1.png",0); 

    if (leftimg.channels() == 3 && rightimg.channels() == 3)
    {
        cvtColor(leftimg, leftimg, CV_BGR2GRAY);
        cvtColor(rightimg, rightimg, CV_BGR2GRAY);
    }

    float scale = 1;
    int row = leftimg.rows * scale;
    int col = leftimg.cols * scale;
    resize(leftimg, leftimg, Size( col, row));
    resize(rightimg,rightimg, Size(col, row));
    Mat depthleft = Mat ::zeros(row,col, CV_32S);
    Mat depthright = Mat ::zeros(row,col, CV_32S);
    Mat lastdisp = Mat ::zeros(row,col, CV_32S);
    int MaxDisparity = 60 * scale;
    int winsize = 31*scale;

    timebegin();
    GetPointDepthLeft(depthleft, leftimg, rightimg, MaxDisparity,  winsize);
    GetPointDepthRight(depthright, leftimg, rightimg, MaxDisparity,  winsize);
    CrossCheckDiaparity(depthleft,depthright, lastdisp, MaxDisparity, winsize);
    timeend("time ");

    MatDataNormal(depthleft,depthleft);
    MatDataNormal(depthright, depthright);
    MatDataNormal(lastdisp, lastdisp);
    namedWindow("left", 0);
    namedWindow("right", 0);
    namedWindow("depthleft", 0);
    namedWindow("depthright", 0);
    namedWindow("lastdisp",0);
    imshow("left", leftimg);
    imshow("right", rightimg);
    imshow("depthleft", depthleft);
    imshow("depthright", depthright);
    imshow("lastdisp",lastdisp);

    string strsave = "result_";
    imwrite(strsave +"depthleft.jpg", depthleft);
    imwrite(strsave +"depthright.jpg", depthright);
    imwrite(strsave +"lastdisp.jpg",lastdisp);
    waitKey(0);
    return 0;
}

  

　　left.png 　　　　　　　　　　　　right.png      

   

　　　　left1_depthleft.jpg                        right1_depthleft.jpg                              lastdisp

OpenCv中实现了三种立体匹配算法：

BM算法

SGBM算法 Stereo Processing by Semiglobal Matching and Mutual Information

GC算法 算法文献：Realistic CG Stereo Image Dataset with Ground Truth Disparity Maps

参考：http://blog.csdn.net/wqvbjhc/article/details/6260844

BM算法：速度很快，效果一般

void BM()
{
　　IplImage * img1 = cvLoadImage("left.png",0);
    IplImage * img2 = cvLoadImage("right.png",0);
    CvStereoBMState* BMState=cvCreateStereoBMState();
    assert(BMState);
    BMState->preFilterSize=9;
    BMState->preFilterCap=31;
    BMState->SADWindowSize=15;
    BMState->minDisparity=0;
    BMState->numberOfDisparities=64;
    BMState->textureThreshold=10;
    BMState->uniquenessRatio=15;
    BMState->speckleWindowSize=100;
    BMState->speckleRange=32;
    BMState->disp12MaxDiff=1;

    CvMat* disp=cvCreateMat(img1->height,img1->width,CV_16S);
    CvMat* vdisp=cvCreateMat(img1->height,img1->width,CV_8U);
    int64 t=getTickCount();
    cvFindStereoCorrespondenceBM(img1,img2,disp,BMState);
    t=getTickCount()-t;
    cout<<"Time elapsed:"<<t*1000/getTickFrequency()<<endl;
    cvSave("disp.xml",disp);
    cvNormalize(disp,vdisp,0,255,CV_MINMAX);
    cvNamedWindow("BM_disparity",0);
    cvShowImage("BM_disparity",vdisp);
    cvWaitKey(0);
    //cvSaveImage("cones\\BM_disparity.png",vdisp);
    cvReleaseMat(&disp);
    cvReleaseMat(&vdisp);
    cvDestroyWindow("BM_disparity");
}

 

　　left.png 　　　　　　　　　　　　right.png                                    disparity.jpg

SGBM算法，作为一种全局匹配算法，立体匹配的效果明显好于局部匹配算法，但是同时复杂度上也要远远大于局部匹配算法。算法主要是参考Stereo Processing by Semiglobal Matching and Mutual Information。

opencv中实现的SGBM算法计算匹配代价没有按照原始论文的互信息作为代价，而是按照块匹配的代价。

参考：http://www.opencv.org.cn/forum.php?mod=viewthread&tid=23854

#include <highgui.h>
#include <cv.h>
#include <cxcore.h>
#include <iostream>
using namespace std;
using namespace cv;
int main()
{

    IplImage * img1 = cvLoadImage("left.png",0);
    IplImage * img2 = cvLoadImage("right.png",0);
    cv::StereoSGBM sgbm;
    int SADWindowSize = 9;
    sgbm.preFilterCap = 63;
    sgbm.SADWindowSize = SADWindowSize > 0 ? SADWindowSize : 3;
    int cn = img1->nChannels;
    int numberOfDisparities=64;
    sgbm.P1 = 8*cn*sgbm.SADWindowSize*sgbm.SADWindowSize;
    sgbm.P2 = 32*cn*sgbm.SADWindowSize*sgbm.SADWindowSize;
    sgbm.minDisparity = 0;
    sgbm.numberOfDisparities = numberOfDisparities;
    sgbm.uniquenessRatio = 10;
    sgbm.speckleWindowSize = 100;
    sgbm.speckleRange = 32;
    sgbm.disp12MaxDiff = 1;
    Mat disp, disp8;
    int64 t = getTickCount();
    sgbm((Mat)img1, (Mat)img2, disp);
    t = getTickCount() - t;
    cout<<"Time elapsed:"<<t*1000/getTickFrequency()<<endl;
    disp.convertTo(disp8, CV_8U, 255/(numberOfDisparities*16.));

    namedWindow("left", 1);
    cvShowImage("left", img1);
    namedWindow("right", 1);
    cvShowImage("right", img2);
    namedWindow("disparity", 1);
    imshow("disparity", disp8);
    waitKey();
    imwrite("sgbm_disparity.png", disp8);   
    cvDestroyAllWindows();
    return 0;
}

　　　　　　　　left.png 　　　　　　right.png                                    disparity.jpg　　　　　　　　　

GC算法 效果最好，速度最慢

void GC()
{
    IplImage * img1 = cvLoadImage("left.png",0);
    IplImage * img2 = cvLoadImage("right.png",0);
    CvStereoGCState* GCState=cvCreateStereoGCState(64,3);
    assert(GCState);
    cout<<"start matching using GC"<<endl;
    CvMat* gcdispleft=cvCreateMat(img1->height,img1->width,CV_16S);
    CvMat* gcdispright=cvCreateMat(img2->height,img2->width,CV_16S);
    CvMat* gcvdisp=cvCreateMat(img1->height,img1->width,CV_8U);
    int64 t=getTickCount();
    cvFindStereoCorrespondenceGC(img1,img2,gcdispleft,gcdispright,GCState);
    t=getTickCount()-t;
    cout<<"Time elapsed:"<<t*1000/getTickFrequency()<<endl;
    //cvNormalize(gcdispleft,gcvdisp,0,255,CV_MINMAX);
    //cvSaveImage("GC_left_disparity.png",gcvdisp);
    cvNormalize(gcdispright,gcvdisp,0,255,CV_MINMAX);
    cvSaveImage("GC_right_disparity.png",gcvdisp);


    cvNamedWindow("GC_disparity",0);
    cvShowImage("GC_disparity",gcvdisp);
    cvWaitKey(0);
    cvReleaseMat(&gcdispleft);
    cvReleaseMat(&gcdispright);
    cvReleaseMat(&gcvdisp);
}

 

           

　　　　　left.png 　　　　　　right.png                disparity.jpg　       

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

参看：http://blog.csdn.net/chenyusiyuan/article/details/5967291