intermotion.m

these m files contain the solution of some problems in the book "communication system" by Haykin

%This function first calculate block motion vectors (with integer pel accuracy), using EBMA(), 
%and then interpolate the block motion vectors to obtain a dense motion field
%An example of main function calling this function is "Prob6_14.m", which can be entered on the command window.
%!!!The function also use the function "EBMA" for motion estimation of every macroblock

function [Target_Img,Anchor_Img,Predict_Img,ox,oy,px,py,PSNR]=InterMotion(TargetName,AnchorName,Img_Height,Img_Width,BlockSize,rangs,range)



%	TargetName,AnchorName: 
%		File Names of Target Frame and Anchor Frame
%	Img_Height,Img_Width:  
%		Image Height and Width of a Frame
%  BlockSize:             
%		The size of Macro Block in Frame is BlockSize(1) by BlockSize(2)
%	rangs,range:      
%		The Search Field in Frame A is from (rangs(1),rangs(2)) to (range(1),range(2))
%  Target_Img,Anchor_Img,Predict_Img:
%		Image Matrix for Target Frame, Anchor Frame, Predicted Frame
%	ox,oy,px,py:
%		The location of Motion vector is (ox,oy), (px,py) for the direction 
%	PSNR
%		The peak signal and noise ratio between original image and predicted image
%	Author: Xiaofeng Xu, Polytechnic University  4/21/2002



%Read images from files
fid = fopen(TargetName,'r');
Target_Img= fread(fid,[Img_Height,Img_Width]);
fclose(fid);

fid = fopen(AnchorName,'r');
Anchor_Img= fread(fid,[Img_Height,Img_Width]);
fclose(fid);

figure;
imshow(uint8(Target_Img));

figure;
imshow(uint8(Anchor_Img));


m=1;

imv=1;
jmv=1;

Predict_Img=Target_Img;
t0=clock;

%First perform Integer pixel BMA
for i=1:BlockSize(1):Img_Height-BlockSize(1)+1
   RangeStart(1)=i+rangs(1);
   RangeEnd(1)=i+BlockSize(1)-1+range(1);
   if RangeStart(1)<1
      RangeStart(1)=1;
   end   
   if RangeEnd(1)>Img_Height
      RangeEnd(1)=Img_Height;
   end
   jmv=1;
   for j=1:BlockSize(2):Img_Width-BlockSize(2)+1
      RangeStart(2)=j+rangs(2);
      RangeEnd(2)=j+BlockSize(2)-1+range(2);
	   if RangeStart(2)<1
   	   RangeStart(2)=1;
	   end   
   	if RangeEnd(2)>Img_Width
      	RangeEnd(2)=Img_Width;
	   end
      [px(m), py(m),Predict_Img(i:i+BlockSize(1)-1,j:j+BlockSize(1)-1)]=EBMA(Target_Img,Anchor_Img,BlockSize,[i,j],RangeStart,RangeEnd); 
      ox(m)=j;
      oy(m)=i;
      %Predict_Img(i:i+BlockSize(1)-1,j:j+BlockSize(1)-1)=Target_Img(py(m)+i:py(m)+i+BlockSize(1)-1,px(m)+j:px(m)+j+BlockSize(1)-1);
      mvfx(imv,jmv)=px(m);
      mvfy(imv,jmv)=py(m);
      jmv=jmv+1;
      m=m+1;
   end
   imv=imv+1;
end

hold on
quiver(ox,oy,px,py);

hold off

%show prediction results using block MV
figure;
imshow(uint8(Predict_Img));


%Caculate the error image
Error_Img=Anchor_Img-Predict_Img;

%Caculate PSNR
PSNR=10*log10(255*255/mean(mean((Error_Img.^2))));
PSNR
etime(clock,t0)


t0=clock;

%Interpolation for motion field
%accomplish factor of 16 interpolation by repeating factor of 2 interpolation 4 times.

for i=1:4
[m_Height,m_Width]=size(mvfx);
Up_mvfx=zeros(m_Height*2,m_Width*2);
Up_mvfx(1:2:m_Height*2,1:2:m_Width*2)=mvfx;
Up_mvfx(1:2:m_Height*2-1,2:2:m_Width*2-1)=(mvfx(:,1:m_Width-1)+mvfx(:,2:m_Width))/2;
Up_mvfx(2:2:m_Height*2-1,1:2:m_Width*2-1)=(mvfx(1:m_Height-1,:)+mvfx(2:m_Height,:))/2;
Up_mvfx(2:2:m_Height*2-1,2:2:m_Width*2-1)=(mvfx(1:m_Height-1,1:m_Width-1)+mvfx(1:m_Height-1,2:m_Width)+mvfx(2:m_Height,1:m_Width-1)+mvfx(2:m_Height,2:m_Width))/4;
mvfx=Up_mvfx;

[m_Height,m_Width]=size(mvfy);
Up_mvfy=zeros(m_Height*2,m_Width*2);
Up_mvfy(1:2:m_Height*2,1:2:m_Width*2)=mvfy;
Up_mvfy(1:2:m_Height*2-1,2:2:m_Width*2-1)=(mvfy(:,1:m_Width-1)+mvfy(:,2:m_Width))/2;
Up_mvfy(2:2:m_Height*2-1,1:2:m_Width*2-1)=(mvfy(1:m_Height-1,:)+mvfy(2:m_Height,:))/2;
Up_mvfy(2:2:m_Height*2-1,2:2:m_Width*2-1)=(mvfy(1:m_Height-1,1:m_Width-1)+mvfy(1:m_Height-1,2:m_Width)+mvfy(2:m_Height,1:m_Width-1)+mvfy(2:m_Height,2:m_Width))/4;
mvfy=Up_mvfy;

end

%obtained predicted image based on interpolated dense motion field. Shift the interpolated motion field by 8,8 

m=1;
Predict_Img=Target_Img;
for i=8:Img_Height-8-1
   for j=8:Img_Width-8-1
      Predict_Img(i,j)=Target_Img(round(mvfy(i-7,j-7)+i),round(mvfx(i-7,j-7)+j));
      %round the generally fractional motion vectors in the dense motion field into integers
      m=m+1;
   end
end

etime(clock,t0)

Error_Img=Anchor_Img-Predict_Img;


%show prediction results using interpolated motion field

PSNR=10*log10(255*255/mean(mean((Error_Img.^2))));
PSNR
%hold on
%quiver(ox,oy,px,py);

%hold off
%axis image

figure;
imshow(uint8(Anchor_Img));
figure;
imshow(uint8(Predict_Img));

figure;
imshow(uint8(Error_Img));