Code covered by the BSD License

# Fast computation of Zernike Radial Polynomials

### Mohammed Sadeq Al-Rawi (view profile)

18 Jan 2012 (Updated )

Q-recursive based iteration function

```function R = zernike_radial_polynomials(n,r)
%  An iterative QuRecursive method (using q-recursive) to generate Zernike
%  radial polynomials in matlab. It accepts as input the moment order n and
%  a vector of r values  (has been written using a vectorized
%  implementation: multiple r values)
%
%
%
%   Written by Mohammed S. Al-Rawi, 2007-2008
%   Last updated 2011.
%   rawi707@yahoo.com
%
% The fast method presented in this code:
%  The method implemented in this work is q-recurseive see Ref.
%  Way, C.-C., Paramesran, R., and Mukandan, R., A comparative analysis of algorithms for fast computation of Zernike moments, Pattern Recognition 36 (2003) 731-742.
%  It uses radial polynomials of fixed order p with a varying index q to
%  compute Zernike moments
%
%   Faster than
%   http://www.mathworks.com/matlabcentral/fileexchange/7687-zernike-polynomials/content/zernpol.m
%
%  Related publications
%  M. S. Al-Rawi, "Fast Computation of Pseudo Zernike Moments," Journal of Real Time Image Processing, vol. 5, pp. 3-10, 2010.
%  M. S. Al-Rawi, "Fast Zernike Moments," Journal of Real Time Image Processing, vol. 3, pp. 89-96, 2009.
%
%
%
%
%
%
%  R(n,n) is obtained from the result with R(n +1)
%  R(0,0) is written as R(0 +1) since the matrix always starts at
%  1 in matlab.
%  e.g. n=4, r=0.4, gives
%
%
%     0.1936         0   -0.3776         0    0.0256

% above, R(4, 4) is R(4 +1) = 0.0256, R(4,2) is R(2 +1) = -0.3776
% R(4, 0) is R(0 +1) = 0.1936
% results form direct computation Paul Fricker 11/13/2006  zernpol(4,0,0.4)
% zernpol(4,0,0.4)
%
% ans =
%
%     0.1936
%
%  performance comparison generating R's up to the 60th order, 0,2,4,..,60
% xx = 0:2:60;
% tic; R = zernpol(max(xx)*ones(length(xx),1) , xx, 0.4 ) ; toc
% Elapsed time is 0.006597 seconds.
% tic; R = zernike_radial_polynomials(60, .4);toc
% Elapsed time is 0.000355 seconds.
%
% so, it's 30 times faster....
% %
% %
% Example: using a vector of r values
%
%
%   tic; R = zernike_radial_polynomials(60, [0.1:0.01: 0.5 ]);toc
%   Elapsed time is 0.000643 seconds.
%   for 41 r values
%

%
% What are Zernike polynomials?
%   Zernike functions, which are an orthogonal basis on the unit
%   circle.  The series representation of the radial Zernike
%   polynomials is
%
%          (n-m)/2
%            __
%    m      \       s                                          n-2s
%   Z(r) =  /__ (-1)  [(n-s)!/(s!((n-m)/2-s)!((n+m)/2-s)!)] * r
%    n      s=0
%
%   The following table shows the first 12 polynomials.
%
%       n    m    Zernike polynomial    Normalization
%       ---------------------------------------------
%       0    0    1                        sqrt(2)
%       1    1    r                           2
%       2    0    2*r^2 - 1                sqrt(6)
%       2    2    r^2                      sqrt(6)
%       3    1    3*r^3 - 2*r              sqrt(8)
%       3    3    r^3                      sqrt(8)
%       4    0    6*r^4 - 6*r^2 + 1        sqrt(10)
%       4    2    4*r^4 - 3*r^2            sqrt(10)
%       4    4    r^4                      sqrt(10)
%       5    1    10*r^5 - 12*r^3 + 3*r    sqrt(12)
%       5    3    5*r^5 - 4*r^3            sqrt(12)
%       5    5    r^5                      sqrt(12)
%       ---------------------------------------------
%
%
%  Direct computation is very slow :(
%

if any( r>1 | r<0 | n<0)
error(':zernike_radial_polynomials either r is less than or greater thatn 1,   r must be between 0 and 1 or n is less than 0.')
end

r=r(:); % columninzing it   :)
R = zeros(n+1, length(r) ); % a mat to store results

%idx0 = logical(r==n);
idx0 = ~logical(r.^n) ;  % if any low r exist, and high n,then treat as 0

idx_nz = ~idx0; % any idx with r non zero

if any(idx0) && ~mod(n,2)   % this is an extreme case. for r=0 and m=0..so, ..
R(0+1, idx0) = (-1)^(n/2);
% exit out of this function, r is 0, and thats it
end

if any(idx_nz) % for any r above zero
r(idx0)=[]; % removing zeros if any
switch n

case 0
R(idx_nz) = 1;
case 1
R(1 +1, idx_nz)  = r;
case 2
R(0 +1, idx_nz) = 2*r.*r- 1;
R(2+1, idx_nz)  = r.*r;
case 3
R(1 +1, idx_nz)  = (3*r.*r-2).*r;
R(3 +1, idx_nz) = r.*r.*r;
otherwise

R(n +1, idx_nz)  = r.^n;
R(n-2 +1, idx_nz) = n*r.^n - (n-1)*r.^(n-2);
r22=(r.*r)';
for m = n-4: -2: mod(n,1)
p=n; q=m;
H3 =  -(4*(q+2)*(q+1))/((p+q +2.0)*(p-q));
H2 =  (H3*(p+q+4)*(p-q-2))/(4*(q+3.0)) + (q+2);
H1 =  ((q+4)*(q+3))/2.0 - H2*(q+4)+ (H3*(p+q+6)*(p-q-4))/8.0;

R(m  +1,idx_nz) =   H1*R(m+4  +1, idx_nz) +  R(m+2 +1, idx_nz) .*(H2+H3./r22) ;

end
end

end

```