matlab opencv 计算视差disparity

0.前言 
视差与深度的关系戳这里。在opencv里应该有三种求视差的方法和一些问题集锦，可参考：http://blog.csdn.net/chenyusiyuan/article/details/5967291/ 
1.代码

void saveDisp(const char* filename, Mat mat);
int main()
{
    char* filename = "../disparity_value.txt";

    Mat left = imread("../image/left_199.png",CV_LOAD_IMAGE_GRAYSCALE);
    Mat right = imread("../image/right_199.png",CV_LOAD_IMAGE_GRAYSCALE);
    double t = (double)getTickCount();
    Mat disp;
    StereoBM bm;  
    CvStereoBMState *BMState = cvCreateStereoBMState();
    bm.state->SADWindowSize = 15;  
    bm.state->minDisparity = 0;  
    bm.state->numberOfDisparities = 80;  
    bm.state->textureThreshold = 10;  
    bm.state->uniquenessRatio = 8;  
    bm.state->speckleWindowSize = 10;  
    bm.state->speckleRange = 32;  
    bm.state->disp12MaxDiff = -1;  
    /*防止左侧黑边*/
    copyMakeBorder(left, left, 0, 0, 80, 0, IPL_BORDER_REPLICATE);  
    copyMakeBorder(right, right, 0, 0, 80, 0, IPL_BORDER_REPLICATE);  

    bm.operator()(left,right,disp,CV_32F);
    disp = disp.colRange(80, left.cols);  

    t = ((double)getTickCount() - t)/getTickFrequency();
    cout<<"time:"<<t<<endl;
    imshow("disp",disp);
    saveDisp(filename,disp);
    cvWaitKey(0);
    return 0;
}

/**
将视差保存到txt,用于MATLAB读取
**/
void saveDisp(const char* filename,  Mat mat)       
{
    int ih = mat.type();
    //float disp = mat.at<float>(0,0); 
    FILE* fp = fopen(filename, "wt");
    fprintf(fp, "%02d\n", mat.rows);
    fprintf(fp, "%02d\n", mat.cols);
    for(int y = 0; y < mat.rows; y++)
    {
        for(int x = 0; x < mat.cols; x++)
        {
            float disp = mat.at<float>(y, x); // 这里视差矩阵是CV_16S 格式的，故用 short 类型读取
            fprintf(fp, "%f\n", disp); // 若视差矩阵是 CV_32F 格式，则用 float 类型读取
        }
    }
    fclose(fp);
}
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

相应的Matlab代码为：
function img = txt2img(filename)
data = importdata(filename);
r = data(1);    % 行数
c = data(2);    % 列数
disp = data(3:end); % 视差
vmin = min(disp);
vmax = max(disp);
disp = reshape(disp, [c,r])'; % 将列向量形式的 disp 重构为 矩阵形式
%  OpenCV 是行扫描存储图像，Matlab 是列扫描存储图像
%  故对 disp 的重新排列是首先变成 c 行 r 列的矩阵，然后再转置回 r 行 c 列
img = uint8( 255 * ( disp - vmin ) / ( vmax - vmin ) );
mesh(disp);
set(gca,'YDir','reverse');  % 通过 mesh 方式绘图时，需倒置 Y 轴方向
axis tight; % 使坐标轴显示范围与数据范围相贴合，去除空白显示区

http://blog.csdn.net/u010025211/article/details/53197337

https://cn.mathworks.com/help/vision/ref/reconstructscene.html

http://www.visionbib.com/bibliography/contents.html

http://docs.opencv.org/3.0-beta/doc/py_tutorials/py_calib3d/py_depthmap/py_depthmap.html

http://docs.opencv.org/trunk/d3/d14/tutorial_ximgproc_disparity_filtering.html

https://dzone.com/articles/computing-disparity-map-opencv

看求disparity的算法挺简单的，本来想自己实现下，

先用matlab写个简单的，就是http://vision.deis.unibo.it/~smatt/Seminars/StereoVision.pdf里面的1领域的window search

clear
clc
img1 = rgb2gray(imread('left1.jpg'));
img2 = rgb2gray(imread('right1.jpg'));
[n,m] = size(img1);


sigma = 1;
gausFilter = fspecial('gaussian', [5,5], sigma);
img1= imfilter(img1, gausFilter, 'replicate');
img2= imfilter(img2, gausFilter, 'replicate');


LARGE = 10000000;
max = 60;
d = zeros(n,m);


for i=1:n
    fprintf('processing %d \n', i);
    for j=1:m
        minDiff = LARGE;
        t = -1;
        for k=0:max
            if j+k<=m
                diff = abs(img1(i,j)-img2(i,j+k));
                if diff < minDiff
                    minDiff = diff ;
                    t = k;
                end
            else
                break;
            end
        end
        d(i, j) = t;
    end
end

但是得不到讲义里面的答案，一个基本的形状也得不到（难道是显示数据的问题？？转成图片这一步没做好？？）

Anyway，用Matlab和opencv自带的函数都OK，自己搞不定就想看下源码，Matlab的源码看着很蛋疼，一下又链到一个不知所云的函数，那就看opencv源码把，找是找到了，可以没什么注释，也咩有说具体算法参考的什么资料，

说不得，还是先看片基础的论文吧，因为网上的博客时候都是讲的差不多，具体讲的也比较少，可找论文又出现问题，根本不知道什么搜索关键字，因为想实现个基本的算法，所以想找最开始提出这种算法的论文，但是年代太久有找不到，真是GG

说不得，看书籍把，一些Computer Vision的书籍有讲这些算法，learning opencv也有讲，但是经历上面的本末倒置后，竟然看不下去，浑浑噩噩过了几天，最好还是找到opencv的sample代码，调用opencv现成的函数

后面也想写全局求disparity矩阵的模拟退火算法，先copy了一个网上可以求极值的模拟退火代码，但是应用不上求disparity矩阵的实际问题，而且速度很慢，估计是哪里想错了吧

教训：明白自己要做什么：是自己写算法还是调用现成的

自己写就要先弄清楚算法步骤，先找简单的博客，先实现简单原始的

调用现成的就可以参考opencv自带的sample例子