function [u,v] = mvg_ex05(threshold,windowsize)

% Read input data:
I1 = double(imread('street1.pgm'));
I2 = double(imread('street2.pgm'));

[w,h] = size(I1);

% Compute gradients:
Ix  = I1(:, [2:end end]) - I1(:, [1 1:end-1]) / 2; % (I(x+1,y)-I(x-1,y)) / 2
Iy  = I1([2:end end], :) - I1([1 1:end-1], :) / 2; % (I(x,y+1)-I(x,y-1)) / 2
It  = I2 - I1;

% Elements of the structure tensor:
Ixx = Ix.*Ix;
Iyy = Iy.*Iy;
Ixy = Ix.*Iy;
Itx = It.*Ix;
Ity = It.*Iy;

for y=windowsize+1:h-windowsize
for x=windowsize+1:w-windowsize

    T = zeros(2,2);
    a = zeros(2,1);
    
    n = 1/(2*windowsize+1)^2;
    
    % Compute integral of structure tensors:
    for j=-windowsize:windowsize
    for i=-windowsize:windowsize
        xi = x+i;
        yj = y+j;
        
        T(1,1) = T(1,1) + Ixx(xi,yj) * n;
        T(1,2) = T(1,2) + Ixy(xi,yj) * n;
        T(2,2) = T(2,2) + Iyy(xi,yj) * n;

        a(1,1) = a(1,1) + Itx(xi,yj) * n;
        a(2,1) = a(2,1) + Ity(xi,yj) * n;
    end
    end
    T(2,1) = T(1,2);
    
    % Compute velocity:
    if det(T) > threshold
        b = -inv(T) * a;
        u(x,y) = b(1);
        v(x,y) = b(2);
    end
end
end

figure; colormap('gray'); imagesc(I1);
figure; colormap('gray'); imagesc(I2);
figure; colormap('gray'); imagesc(u);
figure; colormap('gray'); imagesc(v);