% Function : Half Pixel EBMA
%f1: anchor frame; f2: target frame, fp: predicted image;
%mvx,mvy: store the MV image
%widthxheight: image size; N: block size, R: search range
f1 = fopen('foreman72.Y');
f2 = fopen('foreman66.Y');
f1 = fread(f1,[352,288]);
f2 = fread(f2,[352,288]);
f1 = double(f1)';
f2 = double(f2)';
[height,width] = size(f1);
N = 16;
R = input('\nPlease indicate the search range R\n');
% Target Frame Interpolation
f3 = zeros(2*height,2*width);
tic
for k = 1:height-1
for kk = 1:width-1
f3(2*k-1,2*kk-1) = f2(k,kk);
f3(2*k-1,2*kk) = (f2(k,kk) + f2(k+1,kk))/2;
f3(2*k,2*kk-1) = (f2(k,kk) + f2(k,kk+1))/2;
f3(2*k-1,2*kk-1) = (f2(k,kk) + f2(k,kk+1) + f2(k+1,kk) + f2(k+1,kk+1))/4;
end
end
for k = 1:height-1
f3(2*k-1,2*width-1) = f2(k,width);
f3(2*k,2*width-1) = (f2(k,width)+f2(k+1,width))/2;
end
for k = 1:width-1
f3(2*height-1,2*k-1) = f2(height,k);
f3(2*height-1,2*k) = (f2(height,k)+f2(height,k+1))/2;
end
for k = 1:2*height-1
f3(k,2*width) = f3(k,2*width-1);
end
for k = 1:2*width-1
f3(2*height,k) = f3(2*height-1,k);
end
f3(2*height,2*width) = f3(2*height-1,2*width-1);
mvx = zeros(height/N,width/N);
mvy = zeros(height/N,width/N);
for i = 1 : N : height - N +1 % for every block in the anchor frame : height
for j = 1 : N : width - N +1 % for every block in the anchor frame : width
MAD_min = 256 * N * N;
for k = -R : 0.5 : R
for kk = -R : 0.5 : R % for every search candidate
if (2*(i+k)-1>0)&&(2*(j+kk)-1>0)&&(2*(i+k+N-1)-1<2*height+1)&&(2*(j+kk+N-1)-1<2*width+1)
% evaluate before MAD calculation
flag1_1 = 2*(i+k)-1;flag1_2 = 2*(i+k+N-1)-1;
flag2_1 = 2*(j+kk)-1;flag2_2 = 2*(j+kk+N-1)-1;
flag = f3(flag1_1:2:flag1_2,flag2_1:2:flag2_2);
MAD = sum(sum(abs(f1(i:i+N-1,j:j+N-1)-flag)));
% calculate MAD for this candidate
if MAD < MAD_min
MAD_min = MAD;
dy = k;
dx = kk;
end
end
end
end
fp(i:i+N-1,j:j+N-1) = f3(2*(i+dy)-1:2:2*(i+dy+N-1)-1,2*(j+dx)-1:2:2*(j+dx+N-1)-1);
% put the best matching block in the predicted image
iblk = floor((i-1)/N)+1; % block index
jblk = floor((j-1)/N)+1; % block index
mvx(iblk,jblk) = dx; % record the estimated MV
mvy(iblk,jblk) = dy; % record the estimated MV
end
end
toc
figure;imshow(f1,[]);title('Anchor Frame - Half Pixel EBMA')
figure;imshow(f2,[]);title('Target Frame - Half Pixel EBMA')
figure;quiver(mvx,mvy);title('Move Vector - Half Pixel EBMA')
figure;imshow(fp,[]);title('Predicted Image - Half Pixel EBMA')
figure;imshow(f1-fp,[]);title('Predicted Error Image - Half Pixel EBMA')
% PSNR
f_mse = sum(sum((f1-fp).^2))/(height*width);
maximum=double(max(f1(:)));
fsnr = (maximum^2)/f_mse;
fpsnr = 10*log10(fsnr);
fprintf('\nPSNR(Half Pixel EBMA) is %f.\n',fpsnr);
评论0