In this article
In CorneaForge, I implemented two ways of working with the same native corneal samples: interpolation into maps for viewing, and direct fitting of Zernike coefficients to describe surface shape. Both start from measured points on a polar grid. They serve different purposes, so I keep their transformations explicit.
The elevation-fitting path extracts valid coordinates, forms residuals against a conic reference surface, and solves a least-squares problem on those points. It returns coefficients in micrometers and a fit residual. The rendering path evaluates an interpolant on a pixel grid. Increasing that grid’s resolution changes the display, while the original measurements remain the same.
Here is how those calculations work, using a small constructed surface that makes each step visible.
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 CUse the same weights for the values:
z(q) = 0.5 × 10 + 0.25 × 20 + 0.25 × 40 = 20Inside 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.
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.
On real data, the residual measures how closely the chosen modes fit the available observations. More modes increase flexibility but can make the coefficients sensitive to noise and incomplete sampling. NumPy’s least-squares documentation describes the returned rank and singular values, which help inspect that numerical stability.
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.
My implementation applies ordinary least squares to the available points. Each retained sample has equal weight; I have not added area weighting. This makes the sampling pattern part of the fitted objective, alongside the chosen modes and analysis radius.
A coefficient describes a specified quantity
The elevation path in CorneaForge fits 36 normalized Zernike terms through radial order 7, within a default analysis radius of 4 mm. It uses a conic reference surface, with separate anterior and posterior shape parameters, and forms the residual as reference sag minus measured sag relative to the apex. Preserving that sign convention matters: reversing it reverses the fitted coefficients.
The solver evaluates the basis at the original valid point coordinates and converts its fitted coefficients from millimeters to micrometers. It also returns the root mean square of the remaining residual, in millimeters. These are explicit outputs of a specified transformation, rather than numbers whose meaning can be recovered from a label alone.
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, these coefficients offer a compact representation of surface shape. I can compare that representation with image-based features while retaining the reference surface, radius, normalization and units needed to interpret the experiment.
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.