% M = MOMENTS (IM, P, Q, CENTRAL, SCALED)
%
% Calculates moments of order (P+Q) (can be arrays of indentical length)
% on image IM. If CENTRAL is set to 1 (default: 0), returns translation-
% invariant moments; if SCALED is set to 1 (default: 0), returns scale-
% invariant moments.
%
% After: M. Sonka et al., Image processing, analysis and machine vision.

function m = moments (im,p,q,central,scaled)

	if (nargin < 5), scaled = 0; 	end;
	if (nargin < 4), central = 0; end;
  if (nargin < 3)
  	error ('Insufficient number of parameters.');
  end;
   
	if (length(p) ~= length(q))
  	error ('Arrays P and Q should have equal length.');
  end;
   
  if (scaled & ~central)
  	error ('Scale-invariant moments should always be central.');
  end;

  % xx, yy are grids with co-ordinates
  [xs,ys] = size(im);
  [xx,yy] = meshgrid(-(ys-1)/2:1:(ys-1)/2,-(xs-1)/2:1:(xs-1)/2);
   
	if (central)
      
  	% Calculate zeroth and first order moments
	  m00 = sum(sum(im));
	  m10 = sum(sum(im.*xx));
	  m01 = sum(sum(im.*yy));
      
    % This gives the center of gravity
    xc  = m10/m00;
    yc  = m01/m00;
      
    % Subtract this from the grids to center the object
    xx  = xx - xc;
    yy  = yy - yc;
      
  end;
   
  % Calculate moment(s) (p,q).
  for i = 1:length(p)
		m(i) = sum(sum((xx.^p(i)).*(yy.^q(i)).*im));
  end;
   
  if (scaled)
      
  	c = 1 + (p+q)/2;
      
    % m00 should be known, as scaled moments are always central
    m = m ./ (m00.^c);
      
	end;
	      
return;