ZernikePolynomials.jl

A package to work with Zernike Polynomials in the Julia programming language.
Author rdoelman
Popularity
2 Stars
Updated Last
2 Years Ago
Started In
August 2019

ZernikePolynomials.jl

This package provides functionality to work with Zernike polynomials.

It features the following:

  • functions to convert between common serial indices of Zernike polynomials
  • functions to evaluate the polynomials on a grid
  • functions to estimate Zernike coefficients in a least squares sense from a user-provided input.

Conversion to and between (Noll & OSA/ANSI) sequential indices

This package provides conversion utility functions between three common ways of specifying a Zernike polynomial.

The first is a tuple (m,n) of integers to specify the polynomial Z_n^m(ρ,θ), where n ≥ m, and (ρ,θ) are polar coordinates (radius and angle).

The second (sequential) index is the OSA/ANSI standard index (0,1,2,3,...).

The third is Noll's sequential index (1,2,3,...).

An invalid combination (m,n) returns NaN in conversion.

>> using ZernikePolynomials
>> [mn2OSA(n,m) for n in 0:4, m in -5:5]
>> [mn2Noll(n,m) for n in 0:4, m in -5:5]
>> [OSA2mn(j) for j in 0:5]
>> [Noll2mn(j) for j in 0:5]
>> [Noll2OSA(OSA2Noll(j)) for j = 0:5]

Evaluating Zernike polynomials

The Zernike polynomials (ρ,θ) -> Z_n^m(ρ,θ) or (x,y) -> Z_n^m(x,y) can be obtained as a function as follows:

>> Zernike(1,1)
>> Zernike(1,1,coord=:polar)
>> Z = Zernike(5,index=:Noll,coord=:cartesian) # 5th polynomial by Noll's numbering
>> Z = Zernike(1,1,coord=:cartesian)
>> Z(0.5,0.2) # evaluate the function at cartesian coordinate (0.5,0.2)

Polar coordinates are used by default. The functions return 0 for ρ > 1.

Zernike polynomials and affine combination thereof can easily be evaluated on a grid of points using the function evaluateZernike().

For example, to evaluate 0.5Z_2 + 0.3Z_3 (by OSA numbering) on a 256x256 grid, use the following

>> evaluateZernike(256, [2, 3], [0.5, 0.3], index=:OSA)

To evaluate it on a (square) grid with predefined coordinates, use for example

>> x = LinRange(-2,2,256)
>> ϕ = evaluateZernike(x, [2, 3], [0.5, 0.3], index=:OSA)
>> using Plots
>> heatmap(ϕ)

Here the range x is a range that gives the x and y coordinates.

Zernike polynomial normalization

The Zernike polynomials are normalized according to Thibos et al. - "Standards for Reporting the Optical Aberrations of Eyes", i.e. N_n^m = sqrt(2 (n+1) / (1 + δ(m,0))) where δ(m,0) = 1 for m = 0 and 0 otherwise. These normalization constants can be obtained by normalization(m,n):

>> [normalization(OSA2mn(i)...) for i in 0:5]

Note that the definition on the Wikipedia page on Zernike polynomials is different. Here the polynomials are normalized between [-1,1]. When using or reporting (estimated) Zernike coefficients, it is important to be aware of which normalization has been used.

Estimating Zernike coefficients in a least squares sense

A common use for Zernike polynomials is to approximate a given 2D input. Since Zernike polynomials are often used in optics to approximate a 2D phase, this is the term used in the function documentation.

The function Zernikecoefficients() estimates (in a least-squares sense) the optimal coefficients for a sum of Zernike polynomials. A vector of (sequential) indices should be provided to specify which coefficients should be estimated

>> ϕ = evaluateZernike(256, [2, 3], [0.5, 0.3], index=:OSA)
>> Zernikecoefficients(ϕ, [2, 3]) # ≈ [0.5, 0.3]

It is also possible to specify on which coordinates the phase is defined:

>> x = LinRange(-2,2,256)
>> ϕ = evaluateZernike(x, [2, 3], [0.5, 0.3], index=:Noll)
>> Zernikecoefficients(x, ϕ, [1, 3], index=:Noll) # ≈ [0.0, 0.3]

Note that in this last example the coefficient for the Zernike polynomial with Noll index 1 is estimated, but this is not present in ϕ. Therefore the estimated coefficient is 0.

Plotting example

using Plots
using ZernikePolynomials

x = LinRange(-1,1,256)
ϕ = evaluateZernike(x, [2, 3], [0.5, 0.3], index=:Noll)
f(x,y) = (x^2 + y^2) <= 1 ? 1. : NaN # NaNs are not plotted
heatmap( ϕ .* [f(X,Y) for X in x, Y in x])