双目立体匹配之SAD算法改正错误

话说我今天上午写的博客忘了保存，直接关了电脑，然后我就XX了，没办法，晚上补回来吧，还是老算法，还是熟悉的味道

SAD算法的基本流程：

1.构造一个小窗口，类似与卷积核。

2.用窗口覆盖左边的图像，选择出窗口覆盖区域内的所有像素点。

3.同样用窗口覆盖右边的图像并选择出覆盖区域的像素点。

4.左边覆盖区域减去右边覆盖区域，并求出所有像素点差的绝对值的和。

5.移动右边图像的窗口，重复3，4的动作。（这里有个搜索范围，超过这个范围跳出）

6.找到这个范围内SAD值最小的窗口，即找到了左边图像的最佳匹配的像素块。

下面我写一下对于这个算法的理解：

我们要构造一个窗口D分别用于覆盖左边图像和右边图像，窗口的size大小可以自己去定义，然后用D窗口圈出了AB两幅图像后，我们需要要左右窗口内的

选定的像素做差，然后求取其绝对值，并且对绝对值求和。在视差范围内移动窗口，并重复作差，求取绝对值，并求和。找到最小的匹配块，并且标记记录下来

废话不多说我们看代码，这份代码不是我写的，我写完后发现运行效率超级低，我以为我写错了，

赋值粘贴了别人的代码发现还是很慢，不知道什么原因。求大神指教。

for(i=1+win:1:wL-win)
    for(j=1+win:1:hL-win-dispMax)
        preSAD = 10000000;        temp=0.0;        OptimalDisp = dispMin;
        for(dispRange=dispMin:1:dispMax)
            curSAD=0.0;
            for(x=-win:1:win)
                for(y=-win:1:win)
                    if (j+y+dispRange <= hL)
                        temp=imR(i+x,j+y)-imL(i+x,j+y+dispRange);
                        if(temp<0.0)
                            temp=temp*-1.0;
                        end
                        curSAD=curSAD+temp;
                    end
                end
            end
            %Finding a best disaparty
            if (preSAD > curSAD)
                preSAD = curSAD;
                OptimalDisp = dispRange;
            end
        end
        %Final disparity
        dispMap_SAD(i,j) = OptimalDisp;
    end
end

（其实不难发现，算法中存在大量的冗余计算，主要体现在，你需要反复地计算出，已经计算过的差值，绝对值，求和，所以，能不能只计算一次插值呢，能不能先求取插值，比如生成一个插值图像（imgL-imgR），然后利用窗口D在插值图像上面滑动，计算和，然后在针对每个点进行匹配）= =这份代码也是在网上找到的，发现写代码还是一件费劲的事，自己还是很垃圾，慢慢努力吧

im1=imread('left2.png');
 im2=imread('right2.png');
 
 
if isrgb(im1)
    im1=rgb2gray(im1);
end
%imshow(im1);
im1=double(im1);
 
if isrgb(im2)
    im2=rgb2gray(im2);
end
%figure
%imshow(im2);
im2=double(im2);
 
D=20; %最大视差
N=9; %窗口大小的一半
[H,W]=size(im1);
 
%计算右图减去左图，相减产生D个矩阵放到imgDiff中 
imgDiff=zeros(H,W,D);
e=zeros(H,W);
tic
for i=1:D
    fprintf('%g\n',i)
    e(:,1:(W-i))=abs(im2(:,1:(W-i))- im1(:,(i+1):W));
   %e=conv2(e,e,'same');
    e2=zeros(H,W);%计算窗口内的和
    for y=(N+1):(H-N)
        for x=(N+1):(W-N)
            e2(y,x)=sum(sum(e((y-N):(y+N),(x-N):(x+N))));
        end
    end
    imgDiff(:,:,i)=e2;
end
%20
%找到最小的视差，到dispMap
 
dispMap=zeros(H,W);
 
for x=1:W
   for y=1:H
       %[val,id]=min(imgDiff(y,x,:));
       [val,id]=sort(imgDiff(y,x,:));
        if abs(val(1)-val(2))>10
            dispMap(y,x)=id(1);
        end
    end
end
toc
%显示
dispMap=dispMap*200/D;
dispMap=uint8(dispMap);
%toc
imshow(dispMap)
%imwrite(dispMap,'dispMap.jpg')