ZERNFUN.m and ZERNFUN2.m compute the Zernike functions Znm(r,theta). These functions, which form an orthogonal basis on the unit circle, are used in disciplines such as astronomy, optics, optometry, and ophthalmology to characterize functions and data on a circular domain. ZERNPOL.m computes the Zernike polynomials Znm(r), which are the radial portion of the Zernike functions.
A MATLAB Digest article describing the use of the Zernike functions for analyzing optics data (using a LASIK surgery data as an example) also is available, on the File Exchange as a PDF,
and in HTML at:
1.3.0.1  Updated license 

1.3.0.0  Removed an unneeded intermediate variable from the code. 

1.2.0.0  Updated Help documentation and screenshot. 

1.0.0.0  Need to point to a MATLAB Digest article which is based on this submission. 

Update to the algorithm for computing powers of R (vector of radius values).

Inspired by: showregions.m, quadcc
Inspired: Zernike decomposition, Zernike Polynomial Coefficients for a given Wavefront using Matrix Inversion in Matlab, zernike3, Fast computation of Zernike Radial Polynomials, Zernfun2a.m, showregions.m, quadcc, iahncajigas/nSTAT
Create scripts with code, output, and formatted text in a single executable document.
Paul Dumont (view profile)
Jiang You (view profile)
ming li (view profile)
Raghav Narasimhan (view profile)
Joseph Wang (view profile)
Hi, I am also wondering where can one find the eyedata used Joseph
Govi (view profile)
Hi Paul,
I am trying to do a reconstruction on the Corneal Topographer. Is the 'eyedata' that is used in 'Analyzing LASIK Optical Data Using Zernike Functions' the image from the topographer itself?
If not, what is it and where can I find it?
Thanks in Advance!
Mathieu Gauvin (view profile)
Dear Dr Fricker,
I really enjoyed reading your document "Analyzing LASIK Optical Data Using Zernike Functions". Amazing work!
I am trying to understand the format which your eyedata matrix should be and how to plot it. You probably use something fancier than imagesc here?
Is eye data simply a m x n image? I did save the topo from your example as a jpeg and cropped the axes to try to replicate your example, but I had trouble doing it. Is the raw data and code for your example available somewhere? Thanks!
Tim (view profile)
Where can I find the eyedata for "Analyzing LASIK Optical Data Using Zernike Functions" Everywhere online points to it being found here. HELP!
Eureka Ocean (view profile)
Nicolás Casaballe (view profile)
Hello. I think that the help of ZERNFUN has a typo. It states that "Each element k of M must be a positive integer, with possible values M(k) = N(k) to +N(k)", which is a contradiction. Even in the examples after that, negative values are used. After examining the code, I have arrived to the conclusion that is merely the help that has a typo. Am I right?
On the other hand, the help for the ZERNPOL function does indeed use strictly positive values for M(k). Perhaps this is the source of the problem.
Anael (view profile)
Well documented!
Mathieu (view profile)
Hi,
Thanks for providing these files to the comunity !
However, did you check the othogonality of the functions created?
It seems to me that something's wrong there.
I generated 2 polynoms like it is shown in the doc :
z(idx) = zernfun(3,1,r(idx),theta(idx));
z2(idx) = zernfun(4,0,r(idx),theta(idx));
When I check the orthogonality by simply doing :
sum(sum(z(idx).*z2(idx)))
the result is not 0.
Am I doing someting wrong or is there a pb here?
Thanks in adavnce
Konstantin (view profile)
12 Jun 2012 @JUAN CARDENAS
n = floor((1+sqrt(1+8*p))/2);
m=(n2*p+n.*(n+1));
Very good solution for using any p but
I think you were missing a ''sign for your definition of m. (can be seen when plotting the ZernikePolynomials)
JUAN CARDENAS (view profile)
Thank you Paul.
An update for enabling any P order, in zernfun2, instead of:
n = ceil((3+sqrt(9+8*p))/2);
m = 2*p  n.*(n+2);
change it by:
n = floor((1+sqrt(1+8*p))/2);
m=n2*p+n.*(n+1);
Shahab (view profile)
Hi,
Does anybody know about the discrete Zernike Transform (DZT)?
I am so eager to implement that.
thank you in advance
Vic (view profile)
Thanks for the code! I have a question about the unit of the zernike polynomials generated by this code, is it in [number of waves] or [micrometer]?
Thanks again!
Chauncey Graetzel (view profile)
I was running into problems when using my fitted terms in another program: the same coefficients did not give the same surface.
I found the cause: the cosine and sine terms are inverted in zernfun. M>0 should be a cosine, while M<0 should be a sine.
I've corrected this, starting at line 193
if any(idx_pos)
z(:,idx_pos) = y(:,idx_pos).*cos(theta*m(idx_pos)');
end
if any(idx_neg)
z(:,idx_neg) = y(:,idx_neg).*sin(theta*m(idx_neg)');
end
% note the required sign change "m" in the sine term
Bob (view profile)
Very great functions! But the normalization with the 'norm'  option seems to be wrong. I've changed in:
zernfun.m
% For the normalized polynomials the line 177:
y(:,j) = y(:,j)*sqrt((1+(m(j)~=0))*(n(j)+1)/pi);
% I have to replaced by:
if m(j)==0
y(:,j) = y(:,j)*sqrt((n(j)+1));
else
y(:,j) = y(:,j)*sqrt(2*(n(j)+1));
end
Elia (view profile)
I am still becoming accustomed to MATLAB. When I run the zernfun.m by itself without any modifications, I get this error
??? Input argument "n" is undefined.
Error in ==> zernfun at 88
if ( ~any(size(n)==1) )  ( ~any(size(m)==1) ).
However, when I initialized n, m, r, & theta, I get this error
??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N)
to change the limit. Be aware that exceeding your available stack space can
crash MATLAB and/or your computer.
Error in ==> meshgrid
Does anybody have any thoughts about making zernfun.m to work? Any help would be greatly appreciated.
Jeff (view profile)
Paul, Thanks very much for this code. It is very much appreciated.
I think, however, that you have a bug in how you are normalizing your Zernike polynomials. Noll's paper shows that for terms with m=0 the normalization term should be (n+1)**0.5 and for terms with m!=0 the normalization term should be (2*(n+1))**0.5. In the routine zernpol.m you are normalizing all terms (including those with m=0) by (2*(n+1))**0.5.
Perhaps you have already fixed this? I have just downloaded it on 12110 from this site. The creation date for the file I'm looking at is on Nov 13, 2006 12:56:48 PM.
Alex Kararg (view profile)
I am interested in extracting the 3D Zernike descriptors of 3D shapes? Do you happen to know if there is any Matlab code?
Michael Black (view profile)
Looks promising but I'm having a problem duplicating the LASIK example...I copied the lasik image and edited out the the lasik portion  but all I get from the following is a big red blob  not a decent reconstructed image...anybody know what's wrong?
image=imread('zk_fig2_w.jpg');
figure(1);
I=im2double(image);
imagesc(image);
% make grid coordinate matrices expressed in polar coordinates
L = size(I,1);
X = 1:2/(L1):1;
[x,y] = meshgrid(X);
x=x(:);
y=y(:);
[theta,r] = cart2pol(x,y);
% Compute the required degree and order values from n=07, inclusive
N = [];
M = [];
for n=0:7
N = [N n*ones(1,n+1)];
M = [M n:2:n];
end
is_in_circle = ( r <= 1);
Z = zernfun(N,M,r(is_in_circle),theta(is_in_circle));
a = Z\I(is_in_circle);
% Reconstruct image using Zernike coefficients
r=NaN(size(I));
r(is_in_circle) = Z*a;
% rescale to 0255 to display image
figure(2);
r = im2uint8(im2double(r));
imshow(r);
Valuable work, thank you.
very helpful!It will be perfect if include the zernike annular polynomials
Why this error? Thank you!
??? Error: File: E:\zernike\zernpol.m Line: 151 Column: 25
"identifier" expected, "(" found.
Very helpful, fast, and complete. Examples in the comments are very helpful too.
thank u
thank you very much
Thanks
this the function i was looking for past few weeks