From measured points to Zernike coefficients

Put a dozen coloured points on a blank disk. Each point has a position and a measured value. Between the points, there is empty space.

Now fill the disk with a smooth colour map. It looks like a much richer observation, but the instrument has not acquired anything new. An algorithm has supplied values between the measurements. The image combines observations with assumptions about what happens between them.

This distinction matters in the corneal processing work I do in CorneaForge. A map can be useful for viewing a surface; a mathematical fit can be useful for describing its shape; a model can consume either. Those uses share source data, but they do not produce interchangeable evidence. We can see why with a small synthetic surface, without needing an ophthalmic dataset.

First locate the observations

Suppose an instrument samples a surface on concentric rings. Each observation contains a radius r, an angle θ, and a value z.

To place the observation on a Cartesian diagram, use:

x = r cos(θ)
y = r sin(θ)

For r = 2 and θ = π/2, the location is approximately (0, 2). We have changed the coordinate description of the point. We have not changed its measured value, created neighbouring measurements, or determined how the surface behaves nearby.

The next step depends on what we want. To draw an image, we might evaluate values on a regular pixel grid. To estimate shape coefficients, we can build a fit directly from the original point coordinates. There is no mathematical requirement to rasterize a surface before fitting Zernike functions to it.

That choice is useful in practice. A raster is convenient for visualization or a convolutional network. Direct fitting avoids treating the many interpolated pixels as though they were independent instrument measurements.

Interpolation begins with an assumption

Consider two measurements on a line: 10 at the left endpoint and 20 at the right. A linear interpolant gives 15 halfway between them. The arithmetic is straightforward because we have chosen a rule: between these observations, the interpolated value varies linearly.

A curved physical signal need not pass through that midpoint value. Interpolation produces an estimate consistent with a chosen construction, not a new observation guaranteed by the instrument.

In two dimensions, a triangle provides an equally simple construction. Let its vertices be:

Vertex	Position	Value
A	(0, 0)	10
B	(1, 0)	20
C	(0, 1)	40

At q = (0.25, 0.25), the barycentric weights are 0.5, 0.25, and 0.25. These numbers express the position itself:

q = 0.5 A + 0.25 B + 0.25 C

Use the same weights for the values:

z(q) = 0.5 × 10 + 0.25 × 20 + 0.25 × 40 = 20

Inside the triangle, the weights are nonnegative and sum to one. Along an edge, one weight becomes zero. At a vertex, its weight is one and the others are zero. These are useful checks on an implementation before examining a complicated map.

With many source points, a triangulation defines which triples provide the local interpolation. A Delaunay triangulation uses the source positions; the measured values are supplied afterwards. Its empty-circumcircle condition says that a triangle’s circumcircle contains no other source point strictly inside it. Cocircular configurations can admit more than one valid triangulation, so uniqueness should not be assumed for every grid.

SciPy’s LinearNDInterpolator implements this piecewise linear construction and accepts a precomputed Delaunay triangulation. Its handling of queries outside the convex hull is explicit: the default fill value is NaN.

The convex hull needs a little caution. A point can lie inside the outer boundary while sitting in a large gap between measurements. A triangulation can bridge that gap. Being inside the hull therefore does not establish that the estimate has adequate local support. A domain-specific quality rule may need to reject it.

Describe a surface with a few shapes

An interpolated grid answers “what value should I display here?” A basis expansion asks another question: “how much of each chosen shape explains these observations?”

Zernike polynomials form a family of functions on a disk, widely used in optics. For a deliberately simple demonstration, normalize the radius as ρ = r/R, where R is the analysis radius, and choose three real, unnormalized modes:

Z₀ = 1
Z₁ = 2ρ² − 1
Z₂ = ρ² cos(2θ)

The first adds a constant height. The second changes the surface radially. The third adds an oriented variation. We can construct a surface by choosing their coefficients:

z(ρ, θ) = 0.2 Z₀ + 0.4 Z₁ − 0.3 Z₂

At the center, ρ = 0, so the value is 0.2 − 0.4 = −0.2. At the right edge, ρ = 1, θ = 0, it becomes 0.2 + 0.4 − 0.3 = 0.3. At the upper edge, θ = π/2, the cosine changes sign and the value becomes 0.9.

These are arbitrary units and intentionally unnormalized modes. The coefficients teach the mechanism; they are not clinical measurements and should not be compared directly with a device’s Zernike output.

Surface values are also given in the text below.

Constant, a:

0.20

Radial, b:

0.40

Oriented, c:

−0.30

z = a + b(2ρ² − 1) + cρ² cos(2θ)

Centre: −0.20. Right edge: 0.30. Top edge: 0.90.

Non-normalized modes, arbitrary units, fixed colour scale from −3 to +3. This is a constructed surface, not an optical wavefront measurement or patient map.

Recover the mixture from samples

If we know the sample locations but not the coefficients, we can evaluate each basis function at each location. Arrange those evaluations into a matrix A: one row per sample and one column per mode.

Then A c gives the predicted sample values for a coefficient vector c. Fitting asks for coefficients that make those predictions close to the observed vector z:

choose c to minimize ||A c − z||²

Here is the complete toy calculation:

import numpy as np

rho = np.array([0.0, 1.0, 1.0, 0.5])
theta = np.array([0.0, 0.0, np.pi / 2, np.pi / 4])

A = np.column_stack([
    np.ones_like(rho),
    2 * rho**2 - 1,
    rho**2 * np.cos(2 * theta),
])
coefficients = np.array([0.2, 0.4, -0.3])
observations = A @ coefficients

estimated, _, rank, singular_values = np.linalg.lstsq(
    A, observations, rcond=None
)
residual = observations - A @ estimated

assert rank == 3
assert np.allclose(estimated, coefficients)
assert np.allclose(residual, 0.0)

A has shape (4, 3): four points, three modes. Both coefficient vectors have length three; the observation and residual vectors have length four. The points provide enough independent information to identify the three coefficients, and there is no noise, so the recovery works to numerical precision.

Real data need more scrutiny. A small residual shows that the chosen functions approximate the observed values under this objective. It does not establish accurate extrapolation, stability under noise, or diagnostic usefulness. Adding modes gives the fit more flexibility, which can also make coefficients harder to estimate reliably. NumPy’s least-squares documentation describes the returned rank and singular values, which help inspect the numerical problem.

The coordinate convention is part of the result

Even with identical input values, changing the analysis disk changes the fit. A physical point at radius 2 mm has ρ = 0.5 when R = 4 mm and ρ = 1 when R = 2 mm. The basis functions evaluate differently at that point.

Moving the center changes the coordinates too. Subtracting a constant height, often called removing piston, does not move the coordinate origin. Rotation can change the balance of oriented modes. Comparing two coefficient vectors therefore requires agreement on center, orientation, radius, mode ordering, and normalization.

The familiar orthogonality of Zernike functions also needs its domain. Orthogonality on a continuous disk under a specified integration measure does not make columns of an arbitrarily sampled design matrix orthogonal. Ring sampling, missing sectors, and repeated locations alter the discrete fitting problem. Equal weight per recorded point is not automatically equal weight per unit area.

In the fitting component I inspected, least squares is applied to the available points without an added area-weighting scheme. Describing the implementation that way matters: a mathematically attractive weighting method should not appear in an article as though the code already uses it.

A coefficient describes a specified quantity

One fitting path in CorneaForge works on elevation residuals after subtracting a fitted reference sphere and defining centered coordinates. Those choices belong to the result. A residual is a difference from that reference, not an absolute surface and not a disease label.

A geometric height, a residual elevation, and an optical path difference are also different quantities. They can all be expanded in a Zernike basis. Sharing the basis does not make them interchangeable. Optical calculations may involve multiple surfaces, refractive indices, and ray geometry; a constant multiplier is not a universal conversion between every height map and every wavefront.

The need to state optical conventions is longstanding. The VSIA taskforce’s standards for reporting optical aberrations address reference axes and describing functions so that results can be compared meaningfully.

For an ML pipeline, the coefficients are a compact representation chosen by the engineer. They retain some structure and discard other detail. Their usefulness has to be evaluated for the actual task, just as the usefulness of an image representation does. The scientific work includes making the transformations explicit enough that a later result can be interpreted.

Two quick checks

I double the image width and height. Have I collected four times as many measurements?

No. I have evaluated the interpolant at four times as many pixel locations. The original observations and their uncertainty are unchanged.

Two tools report a coefficient called “coma.” Can I compare the numbers?

Only after checking the fitted quantity, units, disk, coordinate convention, normalization, and mode definition. The shared label is the beginning of that comparison, not its conclusion.