% Estimation du mouvement par technique de BlockMatching

%seq = aviread('foreman_color_qcif.avi');
% Rq0 : aviread remplace par VideoReader dans Matlab R20015a
% Rq2 : Impossible de lire la sequence 'foreman' avec VideoReader (??)

%vidObj = VideoReader('xylophone.mp4');
%vidObj = VideoReader('Tom2.mp4');
vidObj = VideoReader('Vectra21Frames.mpg');

% chargement de la sequence

vidHeight  = vidObj.Height;
vidWidth   = vidObj.Width;
seq = struct('cdata', zeros(vidHeight, vidWidth, 3, 'uint8'), ...
    'colormap', []); 

k = 1;
while hasFrame(vidObj)
    seq(k).cdata = readFrame(vidObj);
    k=k+1;
end

% Affichage de la sequence
figure (1); 
movie(seq);

%   Prelevement des 2 images successives
im1 = seq(15).cdata;
im2 = seq(16).cdata;

% Affichage des 2 images
figure;
imshow(im1);
figure;
imshow(im2);

% Conversion : RGB -> YUV et conservation du Y
YCBCR1 = rgb2ycbcr(im1);
YCBCR2 = rgb2ycbcr(im2);
Y1 = double( YCBCR1(:, :, 1) ); % la luminance de la r?f?rence : Y(t-1)
Y2 = double( YCBCR2(:, :, 1) ); % la luminance de la cible : Y(t)

% Prendre qu'une portion des images => moins de calcul
% Rq : Prendre des dimensions modulo taille du bloc (c.a.d 8)
Y1 = Y1(1:120, 81:160);
Y2 = Y2(1:120, 81:160);

[nblignes, nbcolonnes] = size( Y1 );

% Definition de la fonction MSE
MSE = @(x1, x2) mean2((x1 - x2).^2);

% Boucle sur les blocs "b" de la cible Y(t) soit ici Y2
%   => tenir compte de l'effet de bord 
%       (t.q fenetre de recherche toujours incluse dans la reference 
%       Y(t-1) soit ici Y1)           
%       ("enlever" 8 lignes et 8 colonnes autour de Y2)

% Pour le resultat : affichage de Y2
 
figure(2);
imshow( uint8( Y2 ) ); 

    
    for i = 9 : 8 : nblignes-9
        for j = 9 : 8 : nbcolonnes-9
            
            % prelevement de b
            b = Y2(i:i+7, j:j+7);
            
            % prelevement de la fenetre de recherche "B" dans la reference
            % Y(t-1) soit ici Y1 
            B = Y1(i-8:i+7, j-8:j+7);
            
            % calcul des EQM
            % Rq : nlfilter procede sur chacun de blocs (8x8) de B et en
            % les recouvrant
            EQM = nlfilter(B, [8 8], MSE, b);
            
            % retenir l'EQM minimale
            EQM_min = min( min( EQM ));
            
            % position des elements presentants le min dans EQM
            % Rq : s'il y en a plusieurs, on prendra la 1ere occurence
            [I, J] = find( EQM == EQM_min );
            
            % affichage du vecteur mouvement en sur-impression dans Y2
            hold on;
            X = [ i+4 i+I(1)-4 ];
            Y = [ j+4 j+J(1)-4 ];
            line(X, Y);                  
                       
        end
    end