Zernike Functions

class galsim.zernike.Zernike(coef, R_outer=1.0, R_inner=0.0)[source]

A class to represent a Zernike polynomial series (http://en.wikipedia.org/wiki/Zernike_polynomials#Zernike_polynomials).

Zernike polynomials form an orthonormal basis over the unit circle. The convention used here is for the normality constant to equal the area of integration, which is pi for the unit circle. I.e.,

\[\int_\mathrm{unit circle} Z_i Z_j dA = \pi \delta_{i, j}.\]

Two generalizations of the unit circle Zernike polynomials are also available in this class: annular Zernike polynomials, and polynomials defined over non-unit-radius circles.

Annular Zernikes are orthonormal over an annulus instead of a circle (see Mahajan, J. Opt. Soc. Am. 71, 1, (1981)). Similarly, the non-unit-radius polynomials are orthonormal over a region with outer radius not equal to 1. Taken together, these generalizations yield the orthonormality condition:

\[\int_\mathrm{annulus} Z_i Z_j dA = \pi \left(R_\mathrm{outer}^2 - R_\mathrm{inner}^2\right) \delta_{i, j}\]

where \(0 <= R_\mathrm{inner} < R_\mathrm{outer}\) indicate the inner and outer radii of the annulus over which the polynomials are orthonormal.

The indexing convention for i and j above is that from Noll, J. Opt. Soc. Am. 66, 207-211(1976). Note that the Noll indices begin at 1; there is no Z_0. Because of this, the series coefficients argument coef effectively begins with coef[1] (coef[0] is ignored). This convention is used consistently throughout GalSim, e.g., OpticalPSF, OpticalScreen, zernikeRotMatrix, and zernikeBasis.

As an example, the first few Zernike polynomials in terms of Cartesian coordinates x and y are

Noll index	Polynomial
1	1
2	2 x
3	2 y
4	sqrt(3) (2 (x^2 + y^2) - 1)

A few mathematical convenience operations are additionally available. Zernikes can be added, subtracted, or multiplied together, or multiplied by scalars. Note, however, that two Zernikes can only be combined this way if they have matching R_outer and R_inner attributes. Zernike gradients, Laplacians and the determinant of the Hessian matrix are also available as properties that return new Zernike objects.

Parameters:

coef – Zernike series coefficients. Note that coef[i] corresponds to Z_i under the Noll index convention, and coef[0] is ignored. (I.e., coef[1] is ‘piston’, coef[4] is ‘defocus’, …)
R_outer – Outer radius. [default: 1.0]
R_inner – Inner radius. [default: 0.0]

__add__(rhs)[source]

Add two Zernikes.

Requires that each operand’s R_outer and R_inner attributes are the same.

__sub__(rhs)[source]

Subtract two Zernikes.

Requires that each operand’s R_outer and R_inner attributes are the same.

__neg__()[source]: Negate a Zernike.

__mul__(rhs)[source]

Multiply two Zernikes, or multiply a Zernike by a scalar.

If both operands are Zernikes, then the R_outer and R_inner attributes of each must be the same.

__rmul__(rhs)[source]: Equivalent to obj * rhs. See __mul__ for details.

__call__(x, y)[source]

Evaluate this Zernike polynomial series at Cartesian coordinates x and y. Synonym for evalCartesian.

Parameters:

x – x-coordinate of evaluation points. Can be list-like.
y – y-coordinate of evaluation points. Can be list-like.

Returns:

Series evaluations as numpy array.

evalCartesian(x, y)[source]

Evaluate this Zernike polynomial series at Cartesian coordinates x and y.

Parameters:

x – x-coordinate of evaluation points. Can be list-like.
y – y-coordinate of evaluation points. Can be list-like.

Returns:

Series evaluations as numpy array.

evalCartesianGrad(x, y)[source]

Evaluate the gradient of this Zernike polynomial series at cartesian coordinates x and y.

Parameters:

x – x-coordinate of evaluation points. Can be list-like.
y – y-coordinate of evaluation points. Can be list-like.

Returns:

Tuple of arrays for x-gradient and y-gradient.

evalPolar(rho, theta)[source]

Evaluate this Zernike polynomial series at polar coordinates rho and theta.

Parameters:

rho – radial coordinate of evaluation points. Can be list-like.
theta – azimuthal coordinate in radians (or as Angle object) of evaluation points. Can be list-like.

Returns:

Series evaluations as numpy array.

rotate(theta)[source]

Return new Zernike polynomial series rotated by angle theta.

For example:

>>> Z = Zernike(coefs)
>>> Zrot = Z.rotate(theta)
>>> Z.evalPolar(r, th) == Zrot.evalPolar(r, th + theta)

Parameters:: theta – angle in radians.
Returns:: A new Zernike object.

class galsim.zernike.DoubleZernike(coef, *, uv_outer=1.0, uv_inner=0.0, xy_outer=1.0, xy_inner=0.0)[source]

The Cartesian product of two (annular) Zernike polynomial series. Each double Zernike term is the product of two single Zernike polynomials over separate annuli:

\[DZ_{k,j}(u, v, x, y) = Z_k(u, v) Z_j(x, y)\]

The double Zernike’s therefore satisfy the orthonormality condition:

\[\int_{annuli} DZ_{k,j} DZ_{k',j'} = A_1 A_2 \delta_{k, k'} \delta_{j, j'}\]

where \(A_1\) and \(A_2\) are the areas of the two annuli.

Double Zernikes series are useful for representing the field and pupil dependent wavefronts of a telescope. We adopt the typical convention that the first index (the k index) corresponds to the field-dependence, while the second index (the j index) corresponds to the pupil-dependence.

Parameters:

coef – [kmax, jmax] array of double Zernike coefficients. Like the single Zernike class, the 0th index of either axis is ignored.
uv_outer – Outer radius of the first annulus. [default: 1.0]
uv_inner – Inner radius of the first annulus. [default: 0.0]
xy_outer – Outer radius of the second annulus. [default: 1.0]
xy_inner – Inner radius of the second annulus. [default: 0.0]

galsim.zernike.noll_to_zern(j)[source]

Convert linear Noll index to tuple of Zernike indices. j is the linear Noll coordinate, n is the radial Zernike index and m is the azimuthal Zernike index.

c.f. https://oeis.org/A176988

Parameters:: j – Zernike mode Noll index
Returns:: (n, m) tuple of Zernike indices

galsim.zernike.zernikeRotMatrix(jmax, theta)[source]

Construct Zernike basis rotation matrix. This matrix can be used to convert a set of Zernike polynomial series coefficients expressed in one coordinate system to an equivalent set of coefficients expressed in a rotated coordinate system.

For example:

>>> Z = Zernike(coefs)
>>> R = zernikeRotMatrix(jmax, theta)
>>> rotCoefs = R.dot(coefs)
>>> Zrot = Zernike(rotCoefs)
>>> Z.evalPolar(r, th) == Zrot.evalPolar(r, th + theta)

Note that not all values of jmax are allowed. For example, jmax=2 raises an Exception, since a non-zero Z_2 coefficient will in general rotate into a combination of Z_2 and Z_3 coefficients, and therefore the needed rotation matrix requires jmax=3. If you run into this kind of Exception, raising jmax by 1 will eliminate it.

Also note that the returned matrix is intended to multiply a vector of Zernike coefficients coef indexed following the Noll (1976) convention, which starts at 1. Since python sequences start at 0, we adopt the convention that coef[0] is unused, and coef[i] corresponds to the coefficient multiplying Z_i. As a result, the size of the returned matrix is [jmax+1, jmax+1].

Parameters:

jmax – Maximum Zernike index (in the Noll convention) over which to construct matrix.
theta – angle of rotation in radians.

Returns:

Matrix of size [jmax+1, jmax+1].

galsim.zernike.zernikeBasis(jmax, x, y, R_outer=1.0, R_inner=0.0)[source]

Construct basis of Zernike polynomial series up to Noll index jmax, evaluated at a specific set of points x and y.

Useful for fitting Zernike polynomials to data, e.g.:

>>> x, y, z = myDataToFit()
>>> basis = zernikeBasis(11, x, y)
>>> coefs, _, _, _ = np.linalg.lstsq(basis.T, z)
>>> resids = Zernike(coefs).evalCartesian(x, y) - z

or equivalently:

>>> resids = basis.T.dot(coefs).T - z

Note that since we follow the Noll indexing scheme for Zernike polynomials, which begins at 1, but python sequences are indexed from 0, the length of the leading dimension in the result is jmax+1 instead of jmax. We somewhat arbitrarily fill the 0th slice along the first dimension with 0s (result[0, …] == 0) so that it doesn’t impact the results of np.linalg.lstsq as in the example above.

Parameters:

jmax – Maximum Noll index to use.
x – x-coordinates (can be list-like, congruent to y)
y – y-coordinates (can be list-like, congruent to x)
R_outer – Outer radius. [default: 1.0]
R_inner – Inner radius. [default: 0.0]

Returns:

[jmax+1, x.shape] array. Slicing over first index gives basis vectors corresponding to individual Zernike polynomials.

galsim.zernike.zernikeGradBases(jmax, x, y, R_outer=1.0, R_inner=0.0)[source]

Construct bases of Zernike polynomial series gradients up to Noll index jmax, evaluated at a specific set of points x and y.

Note that since we follow the Noll indexing scheme for Zernike polynomials, which begins at 1, but python sequences are indexed from 0, the length of the leading dimension in the result is jmax+1 instead of jmax. We somewhat arbitrarily fill the 0th slice along the first dimension with 0s (result[0, …] == 0).

Parameters:

jmax – Maximum Noll index to use.
x – x-coordinates (can be list-like, congruent to y)
y – y-coordinates (can be list-like, congruent to x)
R_outer – Outer radius. [default: 1.0]
R_inner – Inner radius. [default: 0.0]

Returns:

[2, jmax+1, x.shape] array. The first index selects the gradient in the x/y direction, slicing over the second index gives basis vectors corresponding to individual Zernike polynomials.

galsim.zernike.doubleZernikeBasis(kmax, jmax, u, v, x, y, *, uv_outer=1.0, uv_inner=0.0, xy_outer=1.0, xy_inner=0.0)[source]

Construct basis of DoubleZernike polynomial series up to Noll indices (kmax, jmax), evaluated at (u, v, x, y).

Parameters:

kmax – Maximum Noll index for first annulus.
jmax – Maximum Noll index for second annulus.
u – Coordinates in the first annulus.
v – Coordinates in the first annulus.
x – Coordinates in the second annulus.
y – Coordinates in the second annulus.
uv_outer – Outer radius of the first annulus.
uv_inner – Inner radius of the first annulus.
xy_outer – Outer radius of the second annulus.
xy_inner – Inner radius of the second annulus.

Returns:

[kmax+1, jmax+1, u.shape] array. Slicing over the first two dimensions gives the basis vectors corresponding to individual DoubleZernike terms.

galsim.zernike.describe_zernike(j)[source]

Create a human-readable string describing the jth (unit circle) Zernike mode as a bivariate polynomial in x and y.

Parameters:: j – Zernike mode Noll index
Returns:: string description of jth mode.