fig6_21_hbma.m

these are the Mtlab functions to plot the figures in the book Communication system (by Haykin)

%for performing HBMA between two images
%main file

%function[]=hbma(FileName1, FileName2,dx,dy, MVname, PImgName);
%Img1=fread(fopen(FileName1),[dx,dy]);
%Img2=fread(fopen(FileName2),[dx,dy]);

dx=352;dy=288;
Img1=fread(fopen('foreman69.Y'),[dx,dy]);
Img2=fread(fopen('foreman72.Y'),[dx,dy]);

img13=zeros(dy,dx);img23=zeros(dy,dx);
img13=Img1'; img23=Img2';
figure(1);
image(img13), colormap(gray(256)); axis image;
figure(2);
image(img23), colormap(gray(256)); axis image;

%level 2 image, use averaging filter
dy2=dy/2; dx2=dx/2;
img12=zeros(dy2,dx2);img22=zeros(dy2,dx2);
img12(1:dy2,1:dx2)=(img13(1:2:dy,1:2:dx)+img13(1:2:dy,2:2:dx)+img13(2:2:dy,1:2:dx)+img13(2:2:dy,2:2:dx))/4;
img22(1:dy2,1:dx2)=(img23(1:2:dy,1:2:dx)+img23(1:2:dy,2:2:dx)+img23(2:2:dy,1:2:dx)+img23(2:2:dy,2:2:dx))/4;


figure(3);
image(img12), colormap(gray(256)); axis image;
figure(4);
image(img22), colormap(gray(256)); axis image;

%level 1 image (top level)  use averaging filter from level 2
dy1=dy2/2; dx1=dx2/2;
img11=zeros(dy1,dx1);img21=zeros(dy1,dx1);
img11(1:dy1,1:dx1)=(img12(1:2:dy2,1:2:dx2)+img12(1:2:dy2,2:2:dx2)+img12(2:2:dy2,1:2:dx2)+img12(2:2:dy2,2:2:dx2))/4;
img21(1:dy1,1:dx1)=(img22(1:2:dy2,1:2:dx2)+img22(1:2:dy2,2:2:dx2)+img22(2:2:dy2,1:2:dx2)+img22(2:2:dy2,2:2:dx2))/4;
figure(5);
image(img11), colormap(gray(256)); axis image;
figure(6);
image(img21), colormap(gray(256)); axis image;


% level 1 matching 
bsize=16; srange=4;
pimg1=zeros(dy1,dx1);
pimg1=img11;
bdy1=ceil(dy1/bsize); bdx1=ceil(dx1/bsize);
mvy1=zeros(bdy1,bdx1); mvx1=zeros(bdy1,bdx1);
[mvy1,mvx1,pimg1]=matching(img11, img21, dy1, dx1, bsize, srange, srange, mvy1, mvx1, pimg1);


figure(7);
image(pimg1), colormap(gray(256)); axis image;axis off;
print -deps2 foreman_HBMA1_pimg.eps
figure(8);
imagesc(abs(img21-pimg1)), colormap(gray(256)); axis image;axis off;
figure(9);
plot_MV_function(dy,dx,bsize*4,mvy1*4,mvx1*4);set(gca,'XTick',[],'YTick',[]);
print -deps2 foreman_HBMA1_MV.eps

% level 2 matching
% interpolate MV from level 1 by replication at block level
bsize=16; srange=4;
pimg2=zeros(dy2,dx2);
pimg2=img12;
bdy2=ceil(dy2/bsize); bdx2=ceil(dx2/bsize);
mvy2=zeros(bdy2,bdx2); mvx2=zeros(bdy2,bdx2);

mvy2(1:2:bdy1*2,1:2:bdx1*2)=2*mvy1(1:bdy1,1:bdx1);
mvy2(1:2:bdy1*2,2:2:bdx1*2)=2*mvy1(1:bdy1,1:bdx1);
mvy2(2:2:bdy1*2,1:2:bdx1*2)=2*mvy1(1:bdy1,1:bdx1);
mvy2(2:2:bdy1*2,2:2:bdx1*2)=2*mvy1(1:bdy1,1:bdx1);

mvx2(1:2:bdy1*2,1:2:bdx1*2)=2*mvx1(1:bdy1,1:bdx1);
mvx2(1:2:bdy1*2,2:2:bdx1*2)=2*mvx1(1:bdy1,1:bdx1);
mvx2(2:2:bdy1*2,1:2:bdx1*2)=2*mvx1(1:bdy1,1:bdx1);
mvx2(2:2:bdy1*2,2:2:bdx1*2)=2*mvx1(1:bdy1,1:bdx1);

[mvy2,mvx2,pimg2]=matching(img12, img22, dy2, dx2, bsize, srange, srange, mvy2, mvx2, pimg2);

figure(10);
image(pimg2), colormap(gray(256)); axis image;axis off;
print -deps2 foreman_HBMA2_pimg.eps
figure(11);
imagesc(abs(img22-pimg2)), colormap(gray(256)); axis image;axis off;
figure(12);
plot_MV_function(dy,dx,bsize*2,mvy2*2,mvx2*2);set(gca,'XTick',[],'YTick',[]);
print -deps2 foreman_HBMA2_MV.eps

% level 3 matching
% interpolate MV from level 1 by replication at block level
bsize=16; srange=4;
pimg3=zeros(dy,dx);
pimg3=img13;
bdy3=ceil(dy/bsize); bdx3=ceil(dx/bsize);
mvy3=zeros(bdy3,bdx3); mvx3=zeros(bdy3,bdx3);

mvy3(1:2:bdy2*2,1:2:bdx2*2)=2*mvy2(1:bdy2,1:bdx2);
mvy3(1:2:bdy2*2,2:2:bdx2*2)=2*mvy2(1:bdy2,1:bdx2);
mvy3(2:2:bdy2*2,1:2:bdx2*2)=2*mvy2(1:bdy2,1:bdx2);
mvy3(2:2:bdy2*2,2:2:bdx2*2)=2*mvy2(1:bdy2,1:bdx2);

mvx3(1:2:bdy2*2,1:2:bdx2*2)=2*mvx2(1:bdy2,1:bdx2);
mvx3(1:2:bdy2*2,2:2:bdx2*2)=2*mvx2(1:bdy2,1:bdx2);
mvx3(2:2:bdy2*2,1:2:bdx2*2)=2*mvx2(1:bdy2,1:bdx2);
mvx3(2:2:bdy2*2,2:2:bdx2*2)=2*mvx2(1:bdy2,1:bdx2);

%integer level
[mvy3,mvx3,pimg3]=matching(img13, img23, dy, dx, bsize, srange, srange, mvy3, mvx3, pimg3);
figure(13);
image(pimg3), colormap(gray(256)); axis image; axis off;
figure(14);
imagesc(abs(img23-pimg3)), colormap(gray(256)); axis image;axis off;
figure(15);
plot_MV_function(dy,dx,bsize,mvy3,mvx3);axis off;


%refine using half-pel at the last level
srange=1;
[mvy3,mvx3,pimg3]=matching_half(img13, img23, dy, dx, bsize, srange, srange, mvy3, mvx3, pimg3);

figure(16);
image(pimg3), colormap(gray(256)); axis image;axis off;
print -deps2 foreman_HBMA3_pimg.eps
figure(17);
imagesc(abs(img23-pimg3)), colormap(gray(256)); axis image;axis off;
figure(18);
plot_MV_function(dy,dx,bsize,mvy3,mvx3);set(gca,'XTick',[],'YTick',[]);
print -deps2 foreman_HBMA3_MV.eps

mse=sum(sum((img23-pimg3).*(img23-pimg3)))/dx/dy;
psnr=10*log10(255*255/mse)