function [Ux,Uy] = imFluxoOpticoDiferencial(im0, im1, Sigmafiltro, s)
% im0, im1: imagens correspondente ao frame 0 e 1 respectivamente. 
% Sigmafiltro: parametro que define o tamanho do filtro Gaussiano.
% s: parametro que estabelece cada cuantos pixels sera calculado
% o fluxo optico.

%Determinando as mascara gaussiana e sua derivada
[gFilt,gxFilt] = MaskGaussiana(Sigmafiltro);

%Suavizando os frames para eliminar os possiveis defeitos
im0blur = conv2(conv2(im0,gFilt','same'),gFilt,'same');
im1blur = conv2(conv2(im1,gFilt','same'),gFilt,'same');

%Determinando as derivadas espaciais e temporais
fx	= conv2(conv2(im0blur,gFilt','same'),gxFilt,'same');
fy	= conv2(conv2(im0blur,gFilt,'same'),gxFilt','same');
ft	= im1blur - im0blur;

%Determinando os componente da matriz G
G(:,:,1) = conv2(conv2(fx.*fx,gFilt','same'),gFilt,'same');
G(:,:,2) = conv2(conv2(fx.*fy,gFilt','same'),gFilt,'same');
G(:,:,3) = G(:,:,2); 
G(:,:,4) = conv2(conv2(fy.*fy,gFilt','same'),gFilt,'same');

%Determinando os componente da matriz b
b(:,:,1) = conv2(conv2(fx.*ft,gFilt','same'),gFilt,'same');
b(:,:,2) = conv2(conv2(fy.*ft,gFilt','same'),gFilt,'same');

[dimy,dimx]=size(im0);
Ux = zeros(dimy/s, ceil(dimx/s));
Uy = zeros(dimy/s, ceil(dimx/s));
cx = 1; 
for x = 1 : s : dimx
        cy = 1;
        for y = 1 : s : dimy
                GG = [G(y,x,1) G(y,x,2);G(y,x,3) G(y,x,4)];
                bb = [b(y,x,1);b(y,x,2)];
                if( rank(GG) < 2 )
                        Ux(cy,cx) = 0;
                        Uy(cy,cx) = 0;
                else
                        u = inv(GG) * bb; 
                        Ux(cy,cx) = u(1);
                        Uy(cy,cx) = u(2);
                end
                cy = cy + 1;
        end
        cx = cx + 1;
end

if(1)
    fg = figure(2);set(fg,'color','w');
    subplot(121);stem(gFilt);axis('square');
    title('Filtro Gaussiano');
    subplot(122);stem(gxFilt);axis('square');
    title('Filtro Gaussiano Derivado'); 
    
    fg = figure(3);set(fg,'color','w');
    subplot(121);
    imshow(im0,[min(im0(:)) max(im0(:))]); 
    title('Frame 0');
    subplot(122);
    imshow(im1,[min(im1(:)) max(im1(:))]); 
    title('Frame 1');
    
    fg = figure(4);set(fg,'color','w');
    subplot(121);
    imshow(im0blur,[min(im0blur(:)) max(im0blur(:))]); 
    title('Filtro pasa-bajo Frame 0');
    subplot(122);
    imshow(im1blur,[min(im1blur(:)) max(im1blur(:))]); 
    title('Filtro pasa-bajo Frame 1');    
    
    fg = figure(5);set(fg,'color','w');    
    subplot(121);
    imshow(fx,[min(fx(:)) max(fx(:))]); 
    title('derivada horizontal');
    subplot(122);
    imshow(fy,[min(fy(:)) max(fy(:))]); 
    title('derivada vertical');

    [xramp,yramp] = meshgrid(1:s:dimx,1:s:dimy);
    fg = figure(6); set(fg,'color','w');
    imagesc(im0(:,:,1)); colormap gray; 
    hold on;quiver( xramp, yramp, Ux, Uy, 5 );axis equal;hold on;
end


function [gFilt,gxFilt] = MaskGaussiana(sigmaBlur)
%Longitude do filtro em fun玢o do sigma 
gBlurSize = 2 * round(2.5 * sigmaBlur) + 1; 
x = [1:gBlurSize] - round((gBlurSize+1)/2);  
%Filtro de suavizamento
gFilt = exp(- x .* x / (2.0*sigmaBlur*sigmaBlur));
gFilt = gFilt / sum(gFilt(:));
%Filtro da derivada
gxFilt = (-x/sigmaBlur^2) .* gFilt;