source: birds/moments.m @ 150

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