function [A_pq] = Zernike_moments(I, p,q)
% 
% 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

% Convert the (x,y) grid into a polar (rho,theta) grid:
%
rhoArray   = sqrt( xArray.^2 + yArray.^2 );
thetaArray = atan2( yArray, xArray ); % mapped to -\pi < theta < +\pi 
indsLTOne  = rhoArray(:) <= 1.0; % we only sum over these numbers

% Compute V_(rhoArray,thetaArray) (the Zernike polynomials evaluated at each point)
%
V = zeros(M,N);
for ii=1:M,
  for jj=1:N,
    V(ii,jj) = Zernike_polynomial(rhoArray(ii,jj),thetaArray(ii,jj), p,q);
  end
end
Vstar = conj(V);

ItimesV = I .* Vstar;

A_pq = (p+1)*sum( ItimesV(indsLTOne) )/pi;