%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% grad.m subroutine for Affine Parameter Calculation
% --------------------------------------------------
%  AUTHOR: Maher Khoury
%    DATE: March 1, 1999
% PURPOSE:
%          Caluclates the image gradients in 3 dimensions given
%          two input frames
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ix, Iy, It] = grad(im1, im2);

%%%% Defime the size of the gradients
sizey = size(im1,1);
sizex = size(im1,2);

%%%% Get neighbouring pixels %%%%
im1_right = im1(1:sizey-1, 2:sizex);
im1_down = im1(2:sizey,1:sizex-1);
im1_down_right = im1(2:sizey,2:sizex);
im2_right = im2(1:sizey-1, 2:sizex);
im2_down = im2(2:sizey,1:sizex-1);
im2_down_right = im2(2:sizey,2:sizex);

%%%% Calculate Ix %%%%
Ix = zeros(sizey,sizex);
Ix(1:sizey-1,1:1:sizex-1) = 0.25.*(im1_right + im2_right + im1_down_right + im2_down_right - ...
   (im1(1:sizey-1, 1:sizex-1) + im2(1:sizey-1,1:sizex-1) + im1_down + im2_down));
Ix(sizey,:) = 0;
Ix(:,sizex) = 0;

%%%% Calculate Iy %%%%
Iy = zeros(sizey,sizex);
Iy(1:sizey-1,1:1:sizex-1) = 0.25.*(im1_down + im2_down + im1_down_right + im2_down_right - ...
    (im1(1:sizey-1, 1:sizex-1) + im2(1:sizey-1, 1:sizex-1) + im1_right + im2_right));
Iy(sizey,:) = 0;
Iy(:,sizex) = 0;

%%%% Calculate It %%%%
It = zeros(sizey,sizex);
It(1:sizey-1,1:1:sizex-1) = 0.25.*(im2(1:sizey-1,1:sizex-1) + im2_right + im2_down + im2_down_right - ...
   (im1(1:sizey-1,1:sizex-1) + im1_right + im1_down + im1_down_right));
It(sizey,:) = 0;
It(:,sizex) = 0;