TP FFT et filtrage linéaire

% Estimation du mouvement par technique de BlockMatching
% 
% %seq = aviread('foreman_color_qcif.avi');
% % Rq0 : aviread remplac?e par VideoReader dans Matlab R20015a
% % Rq2 : Impossible de lire la s?quence 'foreman' avec VideoReader (??)
% 
% vidObj = VideoReader('Vectra21Frames.mpg');

% chargement de la sequence

%vidHeight  = vidObj.Height; 
%vidWidth   = vidObj.Width;	
vidHeight  = 176;
vidWidth   = 320;

% 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);

%   Affichage des 2 images successives

im1 = imread('im1.ppm');
im2 = imread('im2.ppm');

% 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 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 des vecteurs en sur-impression sur 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 procède sur chacun de blocs (8x8) de B et en
            % les recouvrant
            EQM = nlfilter(B, [8 8], MSE, b);
            
            % retenir EQM minimale
            EQM_min = min( min( EQM ));
            
            % position des éléments présentants le min dans EQM
            % Rq : si impressionl 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);
            line(Y,X)
                       
        end
    end