Zernike Polynomials for Circular and Annular Domains

Written by Brett Patterson, UK Astronomy Technology Centre, The Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ. October 2001. Email: brett.patterson@physics.org

Updated: March 2002.

Introduction

This Mathematica Notebook sets up functions that describe the circular Zernike polynomials as described by Noll, and the corresponding annular Zernike polynomials described by Mahajan (see below for references). Each set of polynomials is plotted and tested for othogonality and normalisation, over a circular or annular domain as appropriate. The non-orthogonality and non-normalisation of the circular Zernikes over an annular domain is demonstrated. Finally, an aberration is defined on an annular domain and it is shown that the circular Zernikes, while fitting the data correctly, do not describe the aberration correctly.

For illustration, the Zernike functions are defined in this Notebook, but the functions are more carefully defined in the accompanying Package file, Zernike.m, with consideration for numerical performance and accuracy, and memory useage. For normal usage, the command Get["Zernike.m"] (or <<"Zernike.m") should be used in order to make the Zernike polynomial functions available to a Mathematica session.

Circular Zernikes

Noll, R. J., Zernike polynomials and atmospheric turbulence, J. Opt. Soc. Am., 66, 207, 1976.

Definition

Firstly, define the radial function, Equation (2) in Noll:

R[m_,  
        n_,  ρ_]  :=  (∑_(s  =  0)^((n - Abs[m])/2))(((((-1))^s) (n - s)!)/(s! (((n + m)/2) - 
                      s)! (((n - m)/2) - 
                      s)!)) (ρ^(n - 2  s))

Circular Zernikes defined using Noll's ordering from (Subscript[Z, n]^m) ->  (Subscript[Z, j]), and normalisation over the unit disc (0, 1).
The function selects the jth term by starting with the first term (n = m = 0) and proceeding by incrementing m and n until it reaches the desired jth term.
Note the test "i < 100" sets an arbitrary limit to the order of Zernike that can be calculated.

Zernike[j_,  ρ_,  θ_]  := 
      Module[{m,  n,  i}, 
{m,  n}  =  
          Catch[
For[ 
              n = m = i = 0, 
i < 
                100, 
n += 1, 
Do[ 
                If[EvenQ[n - m], i += 1; If[i == j,  Throw[{m,  n}]]; 
                  If[m ≠ 0,  i += 1;  
                    If[i == j,  
                      Throw[{m,  n}]]]], 
{m, 0, 
                  n}]
]
]; 
(Sqrt[n + 1]) R[m,  n,  ρ] 
          Which[m == 0,  1,  EvenQ[j], (Sqrt[2]) Cos[m θ],  
            OddQ[j], (Sqrt[2]) Sin[m θ]]
]

For convenience add the condition that the Zernikes are zero outside the unit disc (0, 1):

Zernike[
        j_,  ρ_  /;  ρ > 1,  θ_]  :=  0

Test Definition

Print a table of Zernike polynomials:

\(TraditionalForm\`\(aberrations\  = \ {"piston", \ "tip", \ "tilt", \ "defocus", \ "astigmatism", \ "astigmatism", \ "coma", \ "coma", \ "trefoil", \ "trefoil", \ "spherical", \ "5th-order astigmatism", \ "5th-order astigmatism", \ "quadrafoil", \ "quadrafoil", \ "5th-order coma", \ "5th-order coma", \ "5th-order trefoil", \ "5th-order trefoil", \ "pentafoil", \ "pentafoil", \ "5th-order spherical"};\)\)

Transpose[Append[Transpose[Prepend[Table[{j,  
                        Zernike[j,  
                          r,  θ]
,  (cl  =  
                            CoefficientList[
                              az  =  
                                Chop[Expand[
                                    N[Zernike[j,  r,  θ]]]],  r];  
                          lcl  =  Length[cl];  
                          m  =  
                            Which[Length[cl[[lcl]]] == 0,  
                              cl[[lcl]],  Length[cl[[lcl]]] == 1,  
                              cl[[lcl]],  
                              True, cl[[lcl]][[1]]];   m),  
                        Collect[
                            Collect[Expand[az/m],  trig_[__ θ]],  
                            trig_[θ]] //. { 
                            1.  x_  ->  
                              x, -1.  x_ -> -x,  
                            n_ θ  ->  
                              IntegerPart[n] θ}},  {j,  1,  22}], {"Order, \!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\)", "Circular Zernike, (Subscript[Z, j])", "Prefactor to (r^n)", "Polynomial"}]], Prepend[
                aberrations[[Range[22]]], "Aberration"]]]//TableForm

Plot the Zernikes over the unit disc:

Show[
        GraphicsArray[
          Partition[
            Table[DensityPlot[
                Zernike[i,  (Sqrt[(x^2) + (y^2)]),  
                  ArcTan[y,  x]],  {x,  -1,  1},  {y,  -1,  
                  1}, 
Mesh -> False,  
                PlotPoints -> 48,  PlotRange -> All,  
                PlotLabel -> aberrations[[i]], 
                DisplayFunction -> Identity],  {i,  2,  11}],  
            5]]];

Check orthonormalisation, following Noll Equation (3), using the weight 1/π:

Table[
        Integrate[
            Zernike[i,  ρ,  θ] Zernike[
                j,  ρ,  θ] ρ,  {ρ,  0,  
              1},  {θ,  -π,  π}]/π,  {i,  1,  
          11},  {j,  1,  11}] // MatrixForm

Annular Zernikes

Mahajan, V. N., Zernike annular polynomials for imaging systems with annular pupils, J. Opt. Soc. Am., 71, 75, 1981.

Definition

Note that for consistency with Noll's notation, the azimuthal frequency parameter is designated by "m" (Mahajan uses "l").

First define two iterative functions Q and h.

Mahajan, Equation (A4) (This is a specific case for m = 0, and must be defined first)

Q[0,  j_,  u_,  ϵ_]  :=  
      R[0,  2  j,  (Sqrt[u]),  ϵ]

Mahajan, Equation (A13) (This is the general case)

Q[m_,  j_ /; j ≥ 0,  
        u_,  ϵ_]  :=  ((2 (2  j + 2  m - 1) h[m - 1, 
                j,  ϵ])/((j + m) (1 - (ϵ^2)) 
              Q[m - 1, j,  
                0,  ϵ])) (∑_(i  =  0)^j)(( Q[m - 1, i, 0,  ϵ] Q[m - 1, i, u,  ϵ])/h)[m - 1, i,  ϵ]

Mahajan, Equation (A5) (This is a specific case for m = 0, and must be defined first)

h[0, 
        j_,  ϵ_]  :=  ((1 - (ϵ^2))/(2 (2  j + 1)))

Mahajan, Equation (A14) (This is the general case)

h[m_,  
        j_ /; j ≥ 
            0,  ϵ_]  :=  -((2 (2  j + 2  m - 1) 
                Q[m - 1, j + 1, 
                  0,  ϵ])/((j + m) (1 - (ϵ^2)) 
                Q[m - 1, j, 0,  ϵ])) 
        h[m - 1, j,  ϵ]

Now define the specific cases for the radial function:

Mahajan, Equation (13b) (Case for ϵ = 0, reduces to the circular radial function)

R[m_,  n_,  ρ_,  0]  :=  
      R[m,  n,  ρ]

Mahajan, Equation (18) (Case for m = 0 and even n, reduces to a radial function)

R[0,  
        n_,  ρ_,  ϵ_]  :=  R[0, 
          n,  (Sqrt[(((ρ^2) - (ϵ^2))/(1 - (ϵ^2)))])] /; EvenQ[n]

Mahajan, Equation (25)

R[n_,  
        n_,  ρ_,  ϵ_]  :=  (ρ^n)/(Sqrt[(∑_(i  =  0)^n) (ϵ^(2  i))])

Mahajan, Equation (A1) (This is the general case, in terms of the iterative functions Q and h defined above.)

R[m_,  n_,  ρ_,  ϵ_]  :=  
      Module[{j},  
        j = (n - m)/
            2;  (Sqrt[((1 - (ϵ^2))/(2 (2  j + m + 1) 
                    h[m,  j,  ϵ]))]) (ρ^m) 
          Q[m,  j,  (ρ^2),  ϵ]]

Finally, define the full expression for the annular Zernikes, including Noll's ordering from (Subscript[Z, n]^m) ->  (Subscript[Z, j]), and the normalisation over the unit annulus (ϵ, 1).
The function selects the jth term by starting with the first term (n=m=0) and proceeding by incrementing m and n until it reaches the desired jth term.
Note the test "i < 100" sets an arbitrary limit to the order of Zernike that can be calculated.

Mahajan, Equation (8)

Zernike[j_,  ρ_,  θ_,  ϵ_]  := 
      Module[{m,  n,  i}, 
{m,  n}  =  
          Catch[
For[ 
              n = m = i = 0, 
i < 
                100, 
n += 1, 
Do[ 
                If[EvenQ[n - m], i += 1; If[i == j,  Throw[{m,  n}]]; 
                  If[m ≠ 0,  i += 1;  
                    If[i == j,  
                      Throw[{m,  n}]]]], 
{m, 0, 
                  n}]
]
]; 
(Sqrt[n + 1]) 
          Simplify[R[m,  n,  ρ,  ϵ],  
            0 ≤ ϵ ≤ 1] 
          Which[m == 0,  1,  EvenQ[j], (Sqrt[2]) Cos[m θ],  
            OddQ[j], (Sqrt[2]) Sin[m θ]]
]

For convenience, define the annular Zernikes to be zero outside the annulus (ϵ, 1):

Zernike[
        j_,  ρ_,  θ_,  ϵ_]  :=  
      0  /; (ϵ > ρ  ∨  ρ > 1)

Test Definition

Print a table of annular Zernike polynomials,
using the specific example of the obscuration ratio for the Primary Mirror of the Visible & Infrared Survey Telescope (VISTA), ϵ = 1.2 / 4.1 = 0.29.

ϵ =. ;Transpose[Append[Transpose[Prepend[Table[{j,  
                          Zernike[j,  r,  θ,  ϵ],  
                          Expand[
                            Zernike[j,  r,  θ,  1.2/4.1]]},  {j,  
                          1,  11}], {"Order, \!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\)", "Annular Zernike, (Subscript[Z, j])", "VISTA M1"}]], Prepend[
                  aberrations[[Range[11]]], "Aberration"]]]//TableForm

Extend the previous table, excluding the general form (which becomes bulky):

ϵ =. ;
Transpose[Append[Transpose[Prepend[Table[{j,  
                        Collect[
                          Expand[Zernike[j,  r,  θ,  1.2/4.1]],  
                          trig_[θ]]
,  (cl  =  
                            CoefficientList[
                              az  =  
                                Chop[Expand[
                                    Zernike[j,  r,  θ,  
                                      1.2/4.1]]],  r];  
                          lcl  =  Length[cl];  
                          m  =  
                            Which[Length[cl[[lcl]]] == 0,  
                              cl[[lcl]],  Length[cl[[lcl]]] == 1,  
                              cl[[lcl]],  
                              True, cl[[lcl]][[1]]];   m),  
                        Collect[
                            Collect[Expand[az/m],  trig_[__ θ]],  
                            trig_[θ]] /. {1.  x_ -> x}},  {j,  
                        1,  22}], {"Order, \!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\)", "Annular Zernike, (Subscript[Z, j]) (VISTA M1: ϵ = 1.2/4.1)", "Prefactor to (r^n)", "Polynomial"}]], Prepend[
                aberrations[[Range[22]]], "Aberration"]]]//TableForm

This next table is the same as above, but with the general annular form included:

ϵ =. ;
Transpose[Append[Transpose[Prepend[Table[{j,  
                        Zernike[j,  r,  θ,  ϵ],  
                        Collect[
                          Expand[Zernike[j,  r,  θ,  1.2/4.1]],  
                          trig_[θ]]
,  (cl  =  
                            CoefficientList[
                              az  =  
                                Chop[Expand[
                                    Zernike[j,  r,  θ,  
                                      1.2/4.1]]],  r];  
                          lcl  =  Length[cl];  
                          m  =  
                            Which[Length[cl[[lcl]]] == 0,  
                              cl[[lcl]],  Length[cl[[lcl]]] == 1,  
                              cl[[lcl]],  
                              True, cl[[lcl]][[1]]];   m),  
                        Collect[
                            Collect[Expand[az/m],  trig_[__ θ]],  
                            trig_[θ]] /. {1.  x_ -> x}},  {j,  
                        1,  22}], {"Order, \!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\)", "Annular Zernike, (Subscript[Z, j])", "VISTA M1: ϵ = 1.2/4.1", "Prefactor to (r^n)", "Polynomial"}]], Prepend[
                aberrations[[Range[22]]], "Aberration"]]]//TableForm

Plot the annular Zernikes over an annulus with obscuration 1.2 / 4.1 = 0.29 (VISTA Primary Mirror), for example.

Show[
        GraphicsArray[
          Partition[
            Table[DensityPlot[
                Zernike[i,  (Sqrt[(x^2) + (y^2)]),  ArcTan[y,  x],  
                  1.2/4.1],  {x,  -1,  1},  {y,  -1,  
                  1}, 
Mesh -> False,  
                PlotPoints -> 48,  PlotRange -> All,  
                PlotLabel -> aberrations[[i]], 
                DisplayFunction -> Identity],  {i,  2,  11}],  
            5]]];

Check orthonormalisation, following Mahajan Equation (9)

Table[
          Integrate[
              Zernike[i,  ρ,  θ,  1.2/4.1] Zernike[
                  j,  ρ,  θ,  1.2/4.1] ρ,  {ρ,  
                1.2/4.1,  1},  {θ,  -π,  π}]/
            NIntegrate[ρ,  {ρ,  1.2/4.1,  
                1},  {θ,  -π,  π}],  {i,  1,  
            11},  {j,  1,  11}] // Chop // MatrixForm

Compare with the orthonormalisation of circular Zernikes over the same annulus (using standard circular normalisation weight of 1/π):

Table[
        Integrate[
            Zernike[i,  ρ,  θ] Zernike[
                j,  ρ,  θ] ρ,  {ρ,  1.2/4.1,  
              1},  {θ,  -π,  π}]/π,  {i,  1,  
          11},  {j,  1,  11}] // MatrixForm

This result demonstrates the non-orthogonality, and non-normalization of the circular Zernikes over an annulus.

Check convergence of annular to circular for ϵ -> 0

\(TraditionalForm\`Prepend[
        Table[{j, \ 
            Limit[Zernike[
                j, \ ρ, \ θ, \ ϵ], \ ϵ -> 
                0], \ Zernike[j, \ ρ, \ θ]}, \ {j, \ 1, \ 
            22}], \ {"order", \ "annular Zernike with e->0", \ "circular Zernike"}]
 // \ TableForm\)

Try fitting to a test wavefront

Define the wavefront function using annular Zernikes:

wavefront[x_, y_] := 
      Zernike[4, (Sqrt[(x^2) + (y^2)]), ArcTan[y, x], 1.2/4.1] + 
        2 Zernike[11, (Sqrt[(x^2) + (y^2)]), ArcTan[y, x], 1.2/4.1]

Make a data array of wavefront value, making sure to avoid the origin where θ is undefined, and plot the data:

data  =  
        Table[{x,  y,  wavefront[x,  y]},  {x,  -1.01,  1,  
            0.1},  {y,  -1.01,  1,  
            0.1}];
Dimensions[data]

ListDensityPlot[
        Map[Map[##[[3]] &,  #] &, data],  Mesh -> False,  
        PlotRange -> All];

Select the annular subset of the data:

dataDisc  =  
        Select[Flatten[data, 
            1],  (1.2/4.1 < #[[1]]^2 + #[[2]]^2 < 
                1) &];
Dimensions[dataDisc]

Use the SimpleFit function to determine the least-squares coefficients.

Load the SimpleFit function, which is an modified version of the Mathematica Fit function.

<< "SimpleFit.m"

The SimpleFit function determines the least-squares coefficients to the supplied set of basis functions.
If the built-in Fit function was used, the premultipliers in the Zernike polynomials would be included in the coefficients.

SimpleFit[dataDisc,  
        Table[Zernike[i, (Sqrt[(x^2) + (y^2)]),  ArcTan[y,  x],  
            1.2/4.1],  {i,  11}],  {x,  y}] // Chop

As expected, the fit reproduces the true wavefront parameters.

Now do the same using circular Zernikes:

SimpleFit[dataDisc,  
        Table[Zernike[i, (Sqrt[(x^2) + (y^2)]),  ArcTan[y,  x]],  {i,  
            11}],  {x,  y}] // Chop

The result fits the wavefront, but does not describe it using the correct combination of aberrations.