% M = ZERNIKE_MOMENTS (IM, ORDER)
%
% Calculates Zernike moments up to and including ORDER (<= 12) on image IM.
% Default: ORDER = 12.

function m = zernike_moments (im, order)

  if (nargin < 2),             order = 12;                      end;
  if (order < 1 | order > 12), error ('order should be 1..12'); 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);

  % Calculate center of mass and distance of any pixel to it

  m  = moments (im,[0 1 0],[0 0 1],0,0);
  xc = m(2)/m(1); yc = m(3)/m(1);
  xx = xx - xc; yy = yy - yc;

  len     = sqrt(xx.^2+yy.^2);
  max_len = max(max(len));

  % Map pixels to unit circle; prevent divide by zero.

  rho        = len/max_len;
  rho_tmp    = rho; rho_tmp(find(rho==0)) = 1;
  theta      = acos((xx/max_len)./rho_tmp);

  % Flip angle for pixels above center of mass

  yneg            = length(find(yy(:,1)<0));
  theta(:,1:yneg) = 2*pi - theta(:,1:yneg);

  % Calculate coefficients

  c = zeros(order,order);
  s = zeros(order,order);

  i = 1;
  for n = 2:order
    for l = n:-2:0
      r    = polynomial (n,l,rho);
      c    = sum(sum(r.*cos(l*theta)))*((n+1)/(pi*max_len^2));
      s    = sum(sum(r.*sin(l*theta)))*((n+1)/(pi*max_len^2));
      m(i) = sqrt(c^2+s^2);
      i    = i + 1;
    end;
  end;

return

function p = polynomial (n,l,rho)

  switch (n)
    case 2, switch (l)
        case 0, p = 2*(rho.^2)-1;
        case 2, p =   (rho.^2);
      end;
    case 3, switch (l)
        case 1, p = 3*(rho.^3)-2*rho;
        case 3, p =   (rho.^3);
      end;
    case 4, switch (l)
        case 0, p = 6*(rho.^4)-6*(rho.^2)+1;
        case 2, p = 4*(rho.^4)-3*(rho.^2);
        case 4, p =   (rho.^4);
      end;
    case 5, switch (l)
        case 1, p = 10*(rho.^5)-12*(rho.^3)+3*rho;
        case 3, p =  5*(rho.^5)- 4*(rho.^3);
        case 5, p =    (rho.^5);
      end;
    case 6, switch (l)
        case 0, p = 20*(rho.^6)-30*(rho.^4)+12*(rho.^2)-1;
        case 2, p = 15*(rho.^6)-20*(rho.^4)+ 6*(rho.^2);
        case 4, p =  6*(rho.^6)- 5*(rho.^4);
        case 6, p =    (rho.^6);
      end;
    case 7, switch (l)
        case 1, p = 35*(rho.^7)-60*(rho.^5)+30*(rho.^3)-4*rho;
        case 3, p = 21*(rho.^7)-30*(rho.^5)+10*(rho.^3);
        case 5, p =  7*(rho.^7)- 6*(rho.^5);
        case 7, p =    (rho.^7);
      end;
    case 8, switch (l)
        case 0, p = 70*(rho.^8)-140*(rho.^6)+90*(rho.^4)-20*(rho.^2)+1;
        case 2, p = 56*(rho.^8)-105*(rho.^6)+60*(rho.^4)-10*(rho.^2);
        case 4, p = 28*(rho.^8)- 42*(rho.^6)+15*(rho.^4);
        case 6, p =  8*(rho.^8)-  7*(rho.^6);
        case 8, p =    (rho.^8);
      end;
    case 9, switch (l)
        case 1, p = 126*(rho.^9)-280*(rho.^7)+210*(rho.^5)-60*(rho.^3)+5*rho;
        case 3, p =  84*(rho.^9)-168*(rho.^7)+105*(rho.^5)-20*(rho.^3);
        case 5, p =  36*(rho.^9)- 56*(rho.^7)+ 21*(rho.^5);
        case 7, p =   9*(rho.^9)-  8*(rho.^7);
        case 9, p =     (rho.^9);
      end;
    case 10, switch (l)
        case  0, p = 252*(rho.^10)-630*(rho.^8)+560*(rho.^6)-210*(rho.^4)+30*(rho.^2)-1;
        case  2, p = 210*(rho.^10)-504*(rho.^8)+420*(rho.^6)-140*(rho.^4)+15*(rho.^2);
        case  4, p = 129*(rho.^10)-252*(rho.^8)+168*(rho.^6)- 35*(rho.^4);
        case  6, p =  45*(rho.^10)- 72*(rho.^8)+ 28*(rho.^6);
        case  8, p =  10*(rho.^10)-  9*(rho.^8);
        case 10, p =     (rho.^10);
      end;
    case 11, switch (l)
        case  1, p = 462*(rho.^11)-1260*(rho.^9)+1260*(rho.^7)-560*(rho.^5)+105*(rho.^3)-6*rho;
        case  3, p = 330*(rho.^11)- 840*(rho.^9)+ 756*(rho.^7)-280*(rho.^5)+ 35*(rho.^3);
        case  5, p = 165*(rho.^11)- 360*(rho.^9)+ 252*(rho.^7)- 56*(rho.^5);
        case  7, p =  55*(rho.^11)-  90*(rho.^9)+  36*(rho.^7);
        case  9, p =  11*(rho.^11)-  10*(rho.^9);
        case 11, p =     (rho.^11);
      end;
    case 12, switch (l)
        case  0, p = 924*(rho.^12)-2772*(rho.^10)+3150*(rho.^8)-1680*(rho.^6)+420*(rho.^4)-42*(rho.^2)+1;
        case  2, p = 792*(rho.^12)-2310*(rho.^10)+2520*(rho.^8)-1260*(rho.^6)+280*(rho.^4)-21*(rho.^2);
        case  4, p = 495*(rho.^12)-1320*(rho.^10)+1260*(rho.^8)- 504*(rho.^6)+ 70*(rho.^4);
        case  6, p = 220*(rho.^12)- 495*(rho.^10)+ 360*(rho.^8)-  84*(rho.^6);
        case  8, p =  66*(rho.^12)- 110*(rho.^10)+  45*(rho.^8);
        case 10, p =  12*(rho.^12)-  11*(rho.^10);
        case 12, p =     (rho.^12);
      end;
  end;

return