zernike_radial_polynomials_kintner(n, m, r)

function R = zernike_radial_polynomials_kintner(n, m, r)

% Written by Mohammed S. Al-Rawi, 
% [e-mail]
% Last updated, June 2012.

% The method is an implementaiotn to the one shown in (known as the modified Kintner's method):
% A comparative analysis of algorithms for fast computation of Zernike moments
% Pattern Recognition
% Volume 36, Issue 3, March 2003, Pages 731–742
% 
% Chee-Way Chonga, ,  [Author Vitae], P. Raveendranb [Author Vitae], R. Mukundanc [Author Vitae]
% A comparative analysis of algorithms for fast computation of Zernike moments
% Chee-Way Chonga, ,  [Author Vitae], P. Raveendranb [Author Vitae], R. Mukundanc [Author Vitae]

% 
%  Example finding R_00, R_20, R_40, R_60
% R = zernike_radial_polynomials_kintner(6, 0, .1)
% 
% R =
% 
%     1.0000         0   -0.9800         0    0.9406         0   -0.8830
%     
%  Now, generating the same values using the q-recursive, we have to sratr from q=0 until p to get the same values   
%  q_recursive =
% 
%     1.0000         0         0         0         0         0         0
%          0    0.1000         0         0         0         0         0
%    -0.9800         0    0.0100         0         0         0         0
%          0   -0.1970         0    0.0010         0         0         0
%     0.9406         0   -0.0296         0    0.0001         0         0
%          0    0.2881         0   -0.0040         0    0.0000         0
%    -0.8830         0    0.0580         0   -0.0005         0    0.0000
%
% so, the first column of the resutls of the q-recursive is exactly the
% same as the one we generated using Kintner's (p-recursive) method

% to work with the vectorized version, many r values, here's an example
% showing this:
% R = zernike_radial_polynomials_kintner(6, 0, [.1 .2])
% 
% R =
% 
%     1.0000         0   -0.9800         0    0.9406         0   -0.8830
%     1.0000         0   -0.9200         0    0.7696         0   -0.5667
%     


if any( r>1 | r<0 | n<0 | m<0 | m>n | mod(n-m,2))
    error('zernike_radial_polynomials either r is less than or greater than 1,   r must be between 0 and 1 or n is less than 0. All N and M must differ by multiples of 2 (including 0)... ')
end

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

if(n==m)
    R(:, m +1) =  r.^m;
    return;
end


% finding the bases
q=m;
R(:, q +1) =  r.^q;

q2=q+2;
R(:, q2 +1) = q2*(r.^q2)-(q+1)*R(:, q +1);

% now,the rest
m4=m+4;
for p=m4:2:n
    pm2=p-2;
    pm1=p-1;
    k1= (p+q)*(p-q)*(pm2)/2;
    k2= 2*p*(pm1)*(pm2);
    k3= -q*q*(pm1)-p*(pm1)*(pm2);
    k4=(-p*(pm2+q)*(pm2-q))/2;
    R(:, p+1) = ( (k2*r2 +k3).*R(:, pm2 +1) + k4*R(:, p-4 +1) )/k1;
end