full search.m

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

function func_SPIHT_Main
% 立体图像对。

I1=imread('roomL.BMP');
I2=imread('roomr.BMP');

f1=im2double(I2);
f2=im2double(I1);
[height,width]=size(f1)
fp=zeros(height,width);
fr=zeros(height,width);
N=16;R=16;
for i=1:N:height-N+1,
for j=1:N:width-N+1
MAD_min=256*N*N;mvx=0;mvy=0;
for k=-R:1:R,
for l=-R:1:R %for every search candidate
 if(or((i==1&j<width-N+1),(j==1&i<height-N+1))) 
     MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+abs(k):i+abs(k)+N-1,j+abs(l):j+abs(l)+N-1))));

 elseif(i>1&i<height-N+1&j>1&j<width-N+1)
MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+k:i+k+N-1,j+l:j+l+N-1))));

 elseif(or((i==height-N+1&j>1),(j==width-N+1&i>1)))      
     MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i-abs(k):i-abs(k)+N-1,j-abs(l):j-abs(l)+N-1))));

   elseif(i==1&j==width-N+1)      
       MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i+abs(k):i+abs(k)+N-1,j-abs(l):j-abs(l)+N-1))));

   elseif(i==height-N+1&j==1)      
       MAD=sum(sum(abs(f1(i:i+N-1,j:j+N-1)-f2(i-abs(k):i-abs(k)+N-1,j+abs(l):j+abs(l)+N-1))));

   else;
   end

% calculate MAD for this candidate
if MAD<MAD_min
MAD_min=MAD;
dy=k;dx=l;
end

end;end;

ii=round(i/2);jj=round(j/2);dxx=round(dx/2);dyy=round(dy/2);

%%%%%%%%%%%%%%%%  进行图像重建
if(or((i==1&j<width-N+1),(j==1&i<height-N+1)))
  fp(i:i+N-1,j:j+N-1)= f2(i+abs(dy):i+abs(dy)+N-1,j+abs(dx):j+abs(dx)+N-1);  
 

elseif(i>1&i<height-N+1&j>1&j<width-N+1)
     fp(i:i+N-1,j:j+N-1)= f2(i+dy:i+dy+N-1,j+dx:j+dx+N-1);
   
elseif(or((i==height-N+1&j>1),(j==width-N+1&i>1)))   
   fp(i:i+N-1,j:j+N-1)= f2(i-abs(dy):i-abs(dy)+N-1,j-abs(dx):j-abs(dx)+N-1); 
 
elseif(i==1&j==width-N+1)
  fp(i:i+N-1,j:j+N-1)= f2(i+abs(dy):i+abs(dy)+N-1,j-abs(dx):j-abs(dx)+N-1); 
 
elseif(i==height-N+1&j==1)
    fp(i:i+N-1,j:j+N-1)= f2(i-abs(dy):i-abs(dy)+N-1,j+abs(dx):j+abs(dx)+N-1);
  
else;
end
subplot(2,2,1)
imshow(f1)
subplot(2,2,2)
imshow(f2)
subplot(2,2,3)
imshow(fp)
%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;
fr=imsubtract(f1,fp);
subplot(2,2,4)
imshow(fr)
%-----------  输入 ----------------
Orig_I1 = imread('roomL.BMP');
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);
%-----------右图像重建－－－－－－－－－－
fe=imadd(fp,img_spiht2)
%-----------   PSNR 分析  ----------------
Q = 255;
MSE1 = sum(sum((img_spiht1 - Orig_I1) .^ 2)) / nRow / nColumn;
MSE2 = sum(sum((fe -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])
subplot(2,2,2)
imshow(img_spiht1,[0 255])
subplot(2,2,3)
imshow(fr)
subplot(2,2,4)
imshow(img_spiht2)
figure(3)
imshow(fe)

% figure,plot(rate,psnr,'*')