Converting between different Zernike polynomial index systems

The Zernike polynomials are a very popular basis to describe aberrations in the wavefront of optical beams. They are used in ophthalmology  [1], microscopy, and laser metrology, among others  [2]. Because they are orthogonal on the unit disk, they are most suitable to describe deviations from a flat phase front in a single transverse plane of the beam.

They are given by

where we can verify that they are indeed orthogonal on the unit disk:

For a given value of \(n\in\{\{0\}\cup\mathbb{N}\}=\mathbb{N}_0\), the possible \(m\) values are in the range \(m\in\{-n,n\}\) and increase in increments of 2. To simplify the enumeration of Zernike polynomials, various linear index schemes have been devised. In short, each possible pair of \((n,m)\) is assigned to a monotonically increasing index that we will denote \(j\). In this post, we’ll go over some index schemes and write Python scripts that will allow us to map form the linear index \(j\) to the quantum pair \((n,m)\). First, however, it will be useful to investigate the structure of the pairs.

The Triangular Structure of the Quantum Pairs

By simply enumerating the quantum pairs and arranging them on a grid where rows indicate the value of \(n\) and columns the value of \(m\), their triangular structure becomes readily apparent.

Since each row contains one element more than the previous row, the number of elements contained in the first \(k\) rows is given by the \(k\)th triangular number \(T_k\), defined by

For index systems that are monotonic in \(n\), i.e. most of them, we can directly relate the quantum number \(n\) to the index of the smallest triangular number larger than \(j\). Fortunately, the index can easily be determined via the formula  [3]

The remainder \(r=j-T_k+n\) can then be used to determine what column corresponds to that particular \(j\) value. The various linear schemes use different rules to assign a column to the remainder.

Noll Indices

The rules to assign the Noll index is to have it monotonic in \(n\) and then assign in order of increasing \(|m|\), where odd indices correspond to negative values of \(m\) and even indices to positive values of \(m\). The correspondence between the Noll indices and the quantum pairs is the subject of OEIS sequence A176988  [4].

Since the Noll indices start at 1, the first quantum number can be written as one less than the triangular number index, \(n=k-1\).

The \(m\) value can be computed from the remainder \(r=j-T_k+n\). The latter can serve as an index into the list m = [k for k in range(-n,n+2,2)], as it guaranteed to lie in the range \([0,n-1]\). Because of the rules governing the assignment of signs, using the remainder directly as an index into this list is rather tedious. A simpler way is to create a list that generates all the possible values of \(|m|\) and sort them in ascending order. For \(n=3\), that list would be [1,1,3,3]. \(r\) can be used as an index in this list to determine the value of \(|m|\). The sign is then decided by the parity of \(j\) directly, as stated by the Noll rules.

Here’s the Python code that implements this (a C++ implementation is available as part of my Zernike library)

import numpy as np
def NollToQuantum(j):
    triangular_numbers_idx = np.array(np.ceil((1+np.sqrt(1+8*j))/2),dtype=int)-1
    triangular_numbers     = np.array(triangular_numbers_idx*(triangular_numbers_idx+1)/2).astype(int)
    n = triangular_numbers_idx - 1

    r     = j - triangular_numbers
    rpn   = r+n

    # -- Have to work value-per-value here.
    for i in range(n.size):
        m_vec = np.array([j for j in range(-n[i],n[i]+2,2)])
        m_vec = np.sort(np.abs(m_vec))
        m[i]  = (-1)**(j[i]%2)*m_vec[rpn[i]]

    return n, m

\(m\) could be computed directly via the one-liner

m = (-1)**j * ((n % 2) + 2 * np.array((r + ((n+1)%2))/2).astype(int))

but in practice I found the code above to be faster, even though it generated a larger bytecode than the one-liner version.

Let’s work through an example value. For instance, suppose that we have \(j=8\), then the formula above \(k=4\) and thus \(n=3\). The remainder is \(r=8-10+3=1\). From our list above, we have that \(|m|=1\). Since \(j\) is even, \(m\) is positive, and we have \(m=1\). And indeed, \((j=8)\mapsto(3,1)\).

Phasics Indices

The Phasics index system is monotonic in \(n\), as Noll’s is, and also monotonic in \(|m|\). However, this system elects to always assign positive values of \(m\) to lower linear indices, regardless of parity, as stated in the SID4 6.4 User Handbook, Chapter 6.

From these characteristics, we see that Noll and Phasics only differ on rows where \(T_n\) is even, as We can thus use the same basic function above, but add the single-line

m[i] *= (-1)**(triangular_numbers[i]%2+1)

at the end of the loop.

Sources

1. L. Thibos, R. A. Applegate, J. T. Schwiegerling, and R. Webb, "Standards for reporting the optical aberrations of eyes," in Vision Science and Its Applications (Optical Society of America, 2000), p. SuC1.
2. J. C. Wyant and K. Creath, "Basic Wavefront Aberration Theory for Optical Metrology," in Applied Optics and Optical Engineering, Volume XI, R. R. Shannon and J. C. Wyant, eds. (1992), Vol. XI.
3. J. D’Aurizio, "Largest Triangular Number less than a Given Natural Number," https://math.stackexchange.com/q/1417583", accessed Sep. 26th, 2018.
4. R. J. Mathar, "A176988: Triangle read by rows which contains Noll’s indices of Zernike polynomials in row n sorted along increasing index of the azimuthal quantum number," https://oeis.org/A176988, accessed Sep. 26th, 2018.

Appendix: Table of index values generated by the algorithm