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