function [MV, Img] = MVestimate(F2, F1, blockLength,SearchArea)
%MVestimation computes the motion vectors of target frame F2 with
%respect to reference frame F1 
% [MV,  Img] = MVestimate(F2, F1, blockLength, SearhArea), uses logarithmic 
% search to find the motion vectors between frames F2 and F1
% using a block size of length blockLength 
% Img is the predicted image 
% Typical blocklengths are 16 or 8,
% SearchArea is the radius of the neighbourhood within which the best match
% is searched
% Example of use:
%   [MV, Img] = MVestimation('experJASP03\ET01049.JPG','experJASP03\ET01O45.JPG');

% For Color Images
[M,N,c]=size(F1);

if c>1
    ImageIn = rgb2ycbcr(F2);
    ImageIn = ImageIn(:,:,1);
    ImageOut = rgb2ycbcr(F1);
    ImageOut = ImageOut(:,:,1);
    % Color Output
    Img=zeros(M,N,3);
else
    ImageIn = F2;
    ImageOut = F1;
    Img=zeros(M,N);
end
    
p = SearchArea;
n = blockLength;

ImageIn = double(ImageIn);
ImageOut = double(ImageOut);

% Initialize min value
minD = 256*256;

MV = zeros(ceil(M/n), ceil(N/n));

for x = 1:n:M-n+1,
    for y = 1:n:N-n+1,
        b=ImageIn(x:x+n-1,y:y+n-1);
        mcx = 0; mcy = 0;
        offset = p;
        step = n;
        
        while offset>1,
            minDxy = minD;
            mx=0;  my=0;
            offset = ceil(offset/2);
            step = ceil(step/2);
            indLx = mcx - offset;
            indHx = mcx + offset;
            [x y];
            MV(mod(x,n),mod(y,n));
            for inx = indLx:step:indHx,
                lowLimX = x + inx;
                highLimX = lowLimX + n-1;
                
                if ((lowLimX>0)&&(highLimX<=M)),
                    indLy = mcy - offset;
                    indHy = mcy + offset;
                    
                    for iny = indLy:step:indHy,
                        lowLimY = y + iny;
                        highLimY = lowLimY + n-1;

                        if ((lowLimY>0)&&(highLimY<=N)),
                            Diff = abs(b-ImageOut(x+inx:x+n+inx-1,y+iny:y+iny+n-1));
                            MADv = sum(Diff(:));
                            
                            if ( MADv < minDxy),
                                minDxy = MADv;
                                mx=inx;
                                my=iny;
                            end                            
                        end                        
                    end                   
                end
            end
            mcx = mx;  mcy = my;
        end
        MV(fix(x/n)+mod(x,n),fix(y/n)+mod(y,n)) = mcx+sqrt(-1)*mcy;
        if c>1
            Img(x:x+n-1,y:y+n-1,:) = F1(x+mcx:x+n-1+mcx,y+mcy:y+n-1+mcy,:);
        else
            Img(x:x+n-1,y:y+n-1) = F1(x+mcx:x+n-1+mcx,y+mcy:y+n-1+mcy);
        end            
    end
end                                        
Img = uint8(Img);