| 1 | % M = ZERNIKE_MOMENTS (IM, ORDER)
|
|---|
| 2 | %
|
|---|
| 3 | % Calculates Zernike moments up to and including ORDER (<= 12) on image IM.
|
|---|
| 4 | % Default: ORDER = 12.
|
|---|
| 5 |
|
|---|
| 6 | function m = zernike_moments (im, order)
|
|---|
| 7 |
|
|---|
| 8 | if (nargin < 2), order = 12; end;
|
|---|
| 9 | if (order < 1 | order > 12), error ('order should be 1..12'); end;
|
|---|
| 10 |
|
|---|
| 11 | % xx, yy are grids with co-ordinates
|
|---|
| 12 |
|
|---|
| 13 | [xs,ys] = size(im);
|
|---|
| 14 | [xx,yy] = meshgrid(-(ys-1)/2:1:(ys-1)/2,-(xs-1)/2:1:(xs-1)/2);
|
|---|
| 15 |
|
|---|
| 16 | % Calculate center of mass and distance of any pixel to it
|
|---|
| 17 |
|
|---|
| 18 | m = moments (im,[0 1 0],[0 0 1],0,0);
|
|---|
| 19 | xc = m(2)/m(1); yc = m(3)/m(1);
|
|---|
| 20 | xx = xx - xc; yy = yy - yc;
|
|---|
| 21 |
|
|---|
| 22 | len = sqrt(xx.^2+yy.^2);
|
|---|
| 23 | max_len = max(max(len));
|
|---|
| 24 |
|
|---|
| 25 | % Map pixels to unit circle; prevent divide by zero.
|
|---|
| 26 |
|
|---|
| 27 | rho = len/max_len;
|
|---|
| 28 | rho_tmp = rho; rho_tmp(find(rho==0)) = 1;
|
|---|
| 29 | theta = acos((xx/max_len)./rho_tmp);
|
|---|
| 30 |
|
|---|
| 31 | % Flip angle for pixels above center of mass
|
|---|
| 32 |
|
|---|
| 33 | yneg = length(find(yy(:,1)<0));
|
|---|
| 34 | theta(:,1:yneg) = 2*pi - theta(:,1:yneg);
|
|---|
| 35 |
|
|---|
| 36 | % Calculate coefficients
|
|---|
| 37 |
|
|---|
| 38 | c = zeros(order,order);
|
|---|
| 39 | s = zeros(order,order);
|
|---|
| 40 |
|
|---|
| 41 | i = 1;
|
|---|
| 42 | for n = 2:order
|
|---|
| 43 | for l = n:-2:0
|
|---|
| 44 | r = polynomial (n,l,rho);
|
|---|
| 45 | c = sum(sum(r.*cos(l*theta)))*((n+1)/(pi*max_len^2));
|
|---|
| 46 | s = sum(sum(r.*sin(l*theta)))*((n+1)/(pi*max_len^2));
|
|---|
| 47 | m(i) = sqrt(c^2+s^2);
|
|---|
| 48 | i = i + 1;
|
|---|
| 49 | end;
|
|---|
| 50 | end;
|
|---|
| 51 |
|
|---|
| 52 | return
|
|---|
| 53 |
|
|---|
| 54 | function p = polynomial (n,l,rho)
|
|---|
| 55 |
|
|---|
| 56 | switch (n)
|
|---|
| 57 | case 2, switch (l)
|
|---|
| 58 | case 0, p = 2*(rho.^2)-1;
|
|---|
| 59 | case 2, p = (rho.^2);
|
|---|
| 60 | end;
|
|---|
| 61 | case 3, switch (l)
|
|---|
| 62 | case 1, p = 3*(rho.^3)-2*rho;
|
|---|
| 63 | case 3, p = (rho.^3);
|
|---|
| 64 | end;
|
|---|
| 65 | case 4, switch (l)
|
|---|
| 66 | case 0, p = 6*(rho.^4)-6*(rho.^2)+1;
|
|---|
| 67 | case 2, p = 4*(rho.^4)-3*(rho.^2);
|
|---|
| 68 | case 4, p = (rho.^4);
|
|---|
| 69 | end;
|
|---|
| 70 | case 5, switch (l)
|
|---|
| 71 | case 1, p = 10*(rho.^5)-12*(rho.^3)+3*rho;
|
|---|
| 72 | case 3, p = 5*(rho.^5)- 4*(rho.^3);
|
|---|
| 73 | case 5, p = (rho.^5);
|
|---|
| 74 | end;
|
|---|
| 75 | case 6, switch (l)
|
|---|
| 76 | case 0, p = 20*(rho.^6)-30*(rho.^4)+12*(rho.^2)-1;
|
|---|
| 77 | case 2, p = 15*(rho.^6)-20*(rho.^4)+ 6*(rho.^2);
|
|---|
| 78 | case 4, p = 6*(rho.^6)- 5*(rho.^4);
|
|---|
| 79 | case 6, p = (rho.^6);
|
|---|
| 80 | end;
|
|---|
| 81 | case 7, switch (l)
|
|---|
| 82 | case 1, p = 35*(rho.^7)-60*(rho.^5)+30*(rho.^3)-4*rho;
|
|---|
| 83 | case 3, p = 21*(rho.^7)-30*(rho.^5)+10*(rho.^3);
|
|---|
| 84 | case 5, p = 7*(rho.^7)- 6*(rho.^5);
|
|---|
| 85 | case 7, p = (rho.^7);
|
|---|
| 86 | end;
|
|---|
| 87 | case 8, switch (l)
|
|---|
| 88 | case 0, p = 70*(rho.^8)-140*(rho.^6)+90*(rho.^4)-20*(rho.^2)+1;
|
|---|
| 89 | case 2, p = 56*(rho.^8)-105*(rho.^6)+60*(rho.^4)-10*(rho.^2);
|
|---|
| 90 | case 4, p = 28*(rho.^8)- 42*(rho.^6)+15*(rho.^4);
|
|---|
| 91 | case 6, p = 8*(rho.^8)- 7*(rho.^6);
|
|---|
| 92 | case 8, p = (rho.^8);
|
|---|
| 93 | end;
|
|---|
| 94 | case 9, switch (l)
|
|---|
| 95 | case 1, p = 126*(rho.^9)-280*(rho.^7)+210*(rho.^5)-60*(rho.^3)+5*rho;
|
|---|
| 96 | case 3, p = 84*(rho.^9)-168*(rho.^7)+105*(rho.^5)-20*(rho.^3);
|
|---|
| 97 | case 5, p = 36*(rho.^9)- 56*(rho.^7)+ 21*(rho.^5);
|
|---|
| 98 | case 7, p = 9*(rho.^9)- 8*(rho.^7);
|
|---|
| 99 | case 9, p = (rho.^9);
|
|---|
| 100 | end;
|
|---|
| 101 | case 10, switch (l)
|
|---|
| 102 | case 0, p = 252*(rho.^10)-630*(rho.^8)+560*(rho.^6)-210*(rho.^4)+30*(rho.^2)-1;
|
|---|
| 103 | case 2, p = 210*(rho.^10)-504*(rho.^8)+420*(rho.^6)-140*(rho.^4)+15*(rho.^2);
|
|---|
| 104 | case 4, p = 129*(rho.^10)-252*(rho.^8)+168*(rho.^6)- 35*(rho.^4);
|
|---|
| 105 | case 6, p = 45*(rho.^10)- 72*(rho.^8)+ 28*(rho.^6);
|
|---|
| 106 | case 8, p = 10*(rho.^10)- 9*(rho.^8);
|
|---|
| 107 | case 10, p = (rho.^10);
|
|---|
| 108 | end;
|
|---|
| 109 | case 11, switch (l)
|
|---|
| 110 | case 1, p = 462*(rho.^11)-1260*(rho.^9)+1260*(rho.^7)-560*(rho.^5)+105*(rho.^3)-6*rho;
|
|---|
| 111 | case 3, p = 330*(rho.^11)- 840*(rho.^9)+ 756*(rho.^7)-280*(rho.^5)+ 35*(rho.^3);
|
|---|
| 112 | case 5, p = 165*(rho.^11)- 360*(rho.^9)+ 252*(rho.^7)- 56*(rho.^5);
|
|---|
| 113 | case 7, p = 55*(rho.^11)- 90*(rho.^9)+ 36*(rho.^7);
|
|---|
| 114 | case 9, p = 11*(rho.^11)- 10*(rho.^9);
|
|---|
| 115 | case 11, p = (rho.^11);
|
|---|
| 116 | end;
|
|---|
| 117 | case 12, switch (l)
|
|---|
| 118 | case 0, p = 924*(rho.^12)-2772*(rho.^10)+3150*(rho.^8)-1680*(rho.^6)+420*(rho.^4)-42*(rho.^2)+1;
|
|---|
| 119 | case 2, p = 792*(rho.^12)-2310*(rho.^10)+2520*(rho.^8)-1260*(rho.^6)+280*(rho.^4)-21*(rho.^2);
|
|---|
| 120 | case 4, p = 495*(rho.^12)-1320*(rho.^10)+1260*(rho.^8)- 504*(rho.^6)+ 70*(rho.^4);
|
|---|
| 121 | case 6, p = 220*(rho.^12)- 495*(rho.^10)+ 360*(rho.^8)- 84*(rho.^6);
|
|---|
| 122 | case 8, p = 66*(rho.^12)- 110*(rho.^10)+ 45*(rho.^8);
|
|---|
| 123 | case 10, p = 12*(rho.^12)- 11*(rho.^10);
|
|---|
| 124 | case 12, p = (rho.^12);
|
|---|
| 125 | end;
|
|---|
| 126 | end;
|
|---|
| 127 |
|
|---|
| 128 | return
|
|---|
| 129 |
|
|---|
