Index: birds/hu_moments.m
===================================================================
--- birds/hu_moments.m	(revision 74)
+++ birds/hu_moments.m	(revision 74)
@@ -0,0 +1,31 @@
+% M = HU_MOMENTS (IM)
+%
+% Calculates 7 moments of Hu on image IM, invariant to translation, 
+% rotation and scale.
+%
+% After: M. Sonka et al., Image processing, analysis and machine vision.
+
+function m = hu_moments (im)
+
+	p = [ 1 0 2 1 2 0 3 ];
+	q = [ 1 2 0 2 1 3 0 ];
+
+  n = moments(im,p,q,1,1);
+   
+  m(1) = n(2) + n(3);
+  m(2) = (n(3) - n(2))^2   + 4*n(1)^2;
+  m(3) = (n(7) - 3*n(4))^2 + (3*n(5) - n(6))^2;
+  m(4) = (n(7) +   n(4))^2 + (  n(5) + n(6))^2;
+  m(5) = (  n(7) - 3*n(4)) * (n(7) + n(4)) * ...
+           (  (n(7) + n(4))^2 - 3*(n(5) + n(6))^2) + ...
+         (3*n(5) -   n(6)) * (n(5) + n(6)) * ...
+           (3*(n(7) + n(4))^2 -   (n(5) + n(6))^2);
+  m(6) = (n(3) - n(2)) * ((n(7) + n(4))^2 - (n(5) + n(6))^2) + ...
+          4*n(1) * (n(7)+n(4)) * (n(5)+n(6));      
+  m(7) = (3*n(5) -   n(6)) * (n(7) + n(4)) * ...
+           (  (n(7) + n(4))^2 - 3*(n(5) + n(6))^2) - ...
+         (  n(7) - 3*n(4)) * (n(5) + n(6)) * ...
+           (3*(n(7) + n(4))^2 -   (n(5) + n(6))^2);
+           
+return;
+
Index: birds/moments.m
===================================================================
--- birds/moments.m	(revision 74)
+++ birds/moments.m	(revision 74)
@@ -0,0 +1,61 @@
+% 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;
Index: birds/zernike_moments.m
===================================================================
--- birds/zernike_moments.m	(revision 74)
+++ birds/zernike_moments.m	(revision 74)
@@ -0,0 +1,129 @@
+% 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
+
