function [moments] = Hu_moments(I)
% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-04
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
% Epage 
% 
%-----

% Make sure that the pixels of I have values that are indexed from 1,2,...N_g this way we can use MATLAB one based indexing
uElts = unique( I(:) );
N_g   = length(uElts);
Icp   = zeros(size(I));
for ii=1:N_g
  inds = find( I(:)==uElts(ii) );
  Icp(inds) = ii;
end
I = Icp;

[M,N] = size(I);

if( true )
  % the horizontal image axis is *x* and is mapped to [-1,+1] (these are the bin centers):
  dx     = 2/N;
  xs     = linspace(-1+0.5*dx,+1-0.5*dx,N);
  xArray = repmat( xs, [M,1] );
  
  % the vertical image axis is *y* and is mapped to [-1,+1] (these are the bin centers):
  dy     = 2/M;
  ys     = linspace(-1+0.5*dy,+1-0.5*dy,M); ys = ys';
  yArray = repmat( ys, [1,N] ); 
  yArray = flipud( yArray );
else
  xArray = repmat( 1:N, [M,1] );
  yArray = repmat( (1:M)', [1,N] ); yArray = flipud(yArray);
end
  
% compute m_pq for 0 <= p <= 3 and 0 <= q <= 3:
max_p = 3;
max_q = 3;
M_PQ = zeros(max_p+1,max_q+1);
for pp=0:max_p
  for qq=0:max_q
    xAToP = xArray.^pp;
    yAToQ = yArray.^qq;
    tmp   = I .* xAToP .* yAToQ;
    M_PQ(pp+1,qq+1) = sum(tmp(:));
  end
end

% compute the expressions for xbar and ybar:
xbar = M_PQ(2,1)/M_PQ(1,1);
ybar = M_PQ(1,2)/M_PQ(1,1);

% compute the expressions mu_pq for 0 <= p <= 3 and 0 <= q <= 3:
MU_PQ = zeros(max_p+1,max_q+1);
for pp=0:max_p
  for qq=0:max_q
    xAToP = (xArray - xbar).^pp;
    yAToQ = (yArray - ybar).^qq;
    tmp   = I .* xAToP .* yAToQ;
    MU_PQ(pp+1,qq+1) = sum(tmp(:));
  end
end

% compute the expressions for eta_pq for 0 <= p <= 3 and 0 <= q <= 3:
ETA_PQ = zeros(max_p+1,max_q+1);
for pp=0:max_p
  for qq=0:max_q
    gamma = ( pp + qq + 2 )/2;
    ETA_PQ(pp+1,qq+1) = MU_PQ(pp+1,qq+1)/(MU_PQ(1,1)^gamma);
  end
end

% compute the Hu moments:
%
phi1 =  ETA_PQ(3,1) + ETA_PQ(1,3);

phi2 = (ETA_PQ(3,1) - ETA_PQ(1,3))^2 + 4*ETA_PQ(2,2)^2;

phi3 = (ETA_PQ(4,1) - 3*ETA_PQ(2,3))^2 + (ETA_PQ(1,4) - 3*ETA_PQ(3,2))^2;

phi4 = (ETA_PQ(4,1) + ETA_PQ(2,3))^2 + (ETA_PQ(1,4) + ETA_PQ(3,2))^2;

phi5 = (ETA_PQ(4,1) - 3*ETA_PQ(2,3))*(ETA_PQ(4,1) + ETA_PQ(2,3))*( (ETA_PQ(4,1) + ETA_PQ(2,3))^2 - 3*(ETA_PQ(3,2) + ETA_PQ(1,4))^2 ) ...
     + (ETA_PQ(1,4) - 3*ETA_PQ(3,2))*(ETA_PQ(1,4) + ETA_PQ(3,2))*( (ETA_PQ(1,4) + ETA_PQ(3,2))^2 - 3*(ETA_PQ(2,3) + ETA_PQ(4,1))^2 );

phi6 = (ETA_PQ(3,1) - ETA_PQ(1,3))*( (ETA_PQ(4,1) + ETA_PQ(2,3))^2 - (ETA_PQ(3,2) + ETA_PQ(1,4))^2 ) ...
     + 4*ETA_PQ(2,2)*(ETA_PQ(4,1) + ETA_PQ(2,3))*(ETA_PQ(1,4) + ETA_PQ(3,2));

phi7 = (3*ETA_PQ(3,2) - ETA_PQ(1,4))*(ETA_PQ(4,1) + ETA_PQ(2,3))*( (ETA_PQ(4,1) + ETA_PQ(2,3))^2 - 3*(ETA_PQ(3,2) + ETA_PQ(1,4))^2 ) ...
     + (ETA_PQ(4,1) - 3*ETA_PQ(2,3))*(ETA_PQ(3,2) + ETA_PQ(1,4))*( (ETA_PQ(1,4) + ETA_PQ(3,2))^2 - 3*(ETA_PQ(4,1) + ETA_PQ(2,3))^2 );

moments = [ phi1, phi2, phi3, phi4, phi5, phi6, phi7 ];