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