% block_match.m
% Motion Estimation using block matching 
% Match criterion: SSE (sum of square errors)
% Author: Dr. David Pan, for EE 604, Digital Image Processing
% Dept. of ECE, Univ. of Alabama in Huntsville (UAH)

clear all;
close all;
BS = 16;  % Block size is 16 x 16
SW = 8;   % Search window size is +/- 8 pixels
I_prev = vcif2('tennis_cif.yuv',1); % previous frame
I_curr = vcif2('tennis_cif.yuv',2); % current frame

DIM = size(I_curr);
y_dim = DIM(1);
x_dim = DIM(2);

curr_blk = zeros(BS, BS);
candi_blk = zeros(BS, BS);

Num_blk_y = y_dim/BS;
Num_blk_x = x_dim/BS;

mv_y = zeros(Num_blk_y, Num_blk_x);  % Motion vector matrix
mv_x = zeros(Num_blk_y, Num_blk_x);  % Motion vector matrix

SSE_tot = 0;  % SSE of the entire frame

for y = 1: Num_blk_y
    for x = 1: Num_blk_x
        % Get the current block (y, x)
        for i = 1: BS
            for j = 1: BS
                curr_blk (i, j) = I_curr((y-1)*BS + i, (x-1)*BS + j);
            end
        end
        
        % For each current block, search for a match block in the previous frame
        SSE_min = BS*BS*255*255;   % Largest possible sum square error value
        dy_min = 0;
        dx_min = 0;
        
        for dy = -SW:SW
            for dx = -SW:SW
                % Boundary check of the upper left (UL) corner pixel locations
                UL_y = 1 + (y-1)*BS + dy; 
                UL_x = 1 + (x-1)*BS + dx;
                if ((UL_y < 1 || UL_x < 1)) 
                    continue
                end
                % Boundary check of the lower right (LR) corner pixel locations
                LR_y = y*BS + dy; 
                LR_x = x*BS + dx;
                if ((LR_y > y_dim || LR_x > x_dim)) 
                    continue
                end

                % Get the candidate block (y, x)
                for i = 1: BS
                    for j = 1: BS
                        prev_blk (i, j) = I_prev((y-1)*BS + dy + i, (x-1)*BS + dx + j);
                    end
                end
                
                % MAD calculation
                SSE = 0;
                for i = 1: BS
                    for j = 1: BS
                        SSE = SSE + (curr_blk(i,j) - prev_blk(i,j))^2;
                    end
                end
                
                if (SSE < SSE_min)
                    SSE_min = SSE;
                    dy_min = dy;
                    dx_min = dx;
                end
            end
        end
        
        mv_y(y, x) = dy_min;
        mv_x(y, x) = dx_min;
        
        SSE_tot = SSE_tot + SSE_min;
    end
end

% Mean square error without motion compensation 
SSE = 0;
for y = 1: y_dim
    for x = 1: x_dim
        SSE = SSE + (I_curr(y,x) - I_prev(y,x))^2;
    end
end    
MSE = SSE / (y_dim*x_dim)

% Mean square error after motion compensation
MSE_mc = SSE_tot / (y_dim*x_dim)



                        
                        

