room.m

立体图像对编码的MATLAB实现程序。包含了匹配跟小波变换2种算法

function func_SPIHT_Main
I1=imread('roomL.BMP');
I2=imread('roomR.BMP');
f1=double(I2);
f2=double(I1);
[height,width]=size(f1)
fp=zeros(height,width);
fr=zeros(height,width);
N=16;R=28;dy=0;dv=0;mvx=0;mvy=0; 

for j=1:N:width-N+1
    for i=1:N:height-N+1,
MAD_min=256*N*N;dv=0;
if(i<=1+N)
    for l=0:R
        for k=-16:16  
   if(j==1)
       MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+l:i+l+N-1,j+abs(k):j+abs(k)+N-1))));
   elseif(j==width-N+1)
   MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+l:i+l+N-1,j-abs(k):j-abs(k)+N-1))));
   else
        MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+l:i+l+N-1,j+k:j+k+N-1))));
    end
   if MAD<MAD_min
    MAD_min=MAD;
    dx=l;dy=k;
   end
end;end

elseif(i==1+2*N) 
     for l=-R:R
         for k=-15:15
   if(j==1)
       MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+l:i+l+N-1,j+abs(k):j+abs(k)+N-1))));
   elseif(j==width-N+1)
       MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+l:i+l+N-1,j-abs(k):j-abs(k)+N-1))));
   else
        MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+l:i+l+N-1,j+k:j+k+N-1))));
    end
   if MAD<MAD_min
    MAD_min=MAD;
    dx=l;dy=k;
   end
end;end
dv=dx; %%%%%%%%%  将上一块的视差矢量储存
elseif(i<=width-2*N+1)
    t=12;
    for ll=-t:1:t
        for k=-15:1:15%for every search candidate
   if(j==1)
MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+ll+dv:i+ll+dv+N-1,j+abs(k):j+abs(k)+N-1)))); 
   elseif(j==width-N+1) MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+ll+dv:i+ll+dv+N-1,j-abs(k):j-abs(k)+N-1))));
   else
     MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+ll+dv:i+ll+dv+N-1,j+k:j+k+N-1))));
   end
  
if MAD<MAD_min
MAD_min=MAD;
dx=ll+dv;  
 dy=k;
end;end;end
dv=dx ;
else 
  for l=-R:0
     for k=-15:15
    if(j==1)
 MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+l:i+l+N-1,j+abs(k):j+abs(k)+N-1))));
    elseif(j==width-N+1) MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+l:i+l+N-1,j-abs(k):j-abs(k)+N-1))));
    else
         MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+l:i+l+N-1,j+k:j+k+N-1))));
    end
        
if MAD<MAD_min
   MAD_min=MAD;
    dx=l; dy=k;
   end
end ; end ;end
ii=round(i/2);jj=round(j/2);dxx=round(dx/2);dyy=round(dy/2);
if(j==1)
      fp(i:i+N-1,j:j+N-1)= f2(i+dx:i+dx+N-1,j+abs(dy):j+abs(dy)+N-1);  
 
elseif(j==width-N+1)
    fp(i:i+N-1,j:j+N-1)= f2(i+dx:i+dx+N-1,j-abs(dy):j-abs(dy)+N-1);
   
else    
fp(i:i+N-1,j:j+N-1)= f2(i+dx:i+dx+N-1,j+dy:j+dy+N-1);

end

%put the best matching block in the predicted image
iblk=floor((i-1)/N+1); jblk=floor((j-1)/N+1); %block index
mvx(iblk,jblk)=dx; mvy(iblk,jblk)=dy; %record the estimated MV
end;end;
subplot(2,2,1)
imshow(f1,[0 255])
title('右图像')
subplot(2,2,2)
imshow(f2,[0 255])
title('左图像')
subplot(2,2,3)
imshow(fp,[0 255])
title('预测右图像')
fr=imsubtract(fp,f1);
subplot(2,2,4)
imshow(fr,[0 255])
title('残差图像')
%-----------  输入 ----------------
Orig_I1 = I1;
Orig_I1=double(Orig_I1);
rate = 0.5;
%-----------  预处理 ----------------
OrigSize = size(Orig_I1, 1);
max_bits = floor(rate * OrigSize^2);
OutSize = OrigSize;
image_spiht1 = zeros(size(Orig_I1));
image_spiht2 = zeros(size(Orig_I1));
% "image " is the input of codec
[nRow, nColumn] = size(Orig_I1);
%-----------   Wavelet Decomposition   ----------------
n = size(Orig_I1,1);
n_log = log2(n); 
level = n_log;
% wavelet decomposition level can be defined by users manually.
type = 'bior4.4';
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(type);
[I_W1, S1] = func_DWT(Orig_I1, level, Lo_D, Hi_D);
[I_W2, S2] = func_DWT(fr, level, Lo_D, Hi_D);

%-----------   编码 ----------------
img_enc1 = func_SPIHT_Enc(I_W1, max_bits, nRow*nColumn, level);   
img_enc2 = func_SPIHT_Enc(I_W2, max_bits, nRow*nColumn, level);  

%-----------   解码  ----------------
img_dec1 = func_SPIHT_Dec(img_enc1);
img_dec2 = func_SPIHT_Dec(img_enc2);

%-----------   小波重建   ----------------
img_spiht1 = func_InvDWT(img_dec1, S1, Lo_R, Hi_R, level);
img_spiht2 = func_InvDWT(img_dec2, S2, Lo_R, Hi_R, level);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%右图像重建%%%%%
for j=1:N:width-N+1
    for i=1:N:height-N+1,
        iblk=floor((i-1)/N+1); jblk=floor((j-1)/N+1);
if(j==1)
      fpp(i:i+N-1,j:j+N-1)= f2(i+mvx(iblk,jblk):i+mvx(iblk,jblk)+N-1,j+abs(mvy(iblk,jblk)):j+abs(mvy(iblk,jblk))+N-1);  
  elseif(j==width-N+1)
    fpp(i:i+N-1,j:j+N-1)= f2(i+mvx(iblk,jblk):i+mvx(iblk,jblk)+N-1,j-abs(mvy(iblk,jblk)):j-abs(mvy(iblk,jblk))+N-1);
else    
fpp(i:i+N-1,j:j+N-1)= f2(i+mvx(iblk,jblk):i+mvx(iblk,jblk)+N-1,j+mvy(iblk,jblk):j+mvy(iblk,jblk)+N-1);
end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------   PSNR 分析  ----------------
Q = 255;
MSE1 = sum(sum((img_spiht1 - Orig_I1) .^ 2)) / nRow / nColumn;
MSE2 = sum(sum((fpp -f1) .^ 2)) / nRow / nColumn;
psnr1 = 10*log10(Q*Q/MSE1)
psnr2 = 10*log10(Q*Q/MSE2)
figure(2)
subplot(2,2,1)
imshow(Orig_I1,[0 255])
title('左图像')
subplot(2,2,2)
imshow(img_spiht1,[0 255])
title('解码后左图像')
subplot(2,2,3)
imshow(fr,[0 255])
title('残差图像')
subplot(2,2,4)
imshow(img_spiht2,[0 255])
title('解码后残差图像')
figure(3)
imshow(uint8(fpp))
title('重建右图像')
% subplot(1,2,2)
% imhist(fr,128)
% psnr=(20.3664,23.6279,26.2184,27.9317)
% rate1=0.1:0.1:0.4
% plot(rate1,psnr,'*')