In this article
I implemented CorneaForge’s polar-to-Cartesian interpolation so that maps sharing a sampling grid can reuse the same geometric calculation. The code constructs a triangulation, finds the triangle containing each output pixel, and stores three source indices and three interpolation weights. Applying the next map becomes a gather, multiplication, and sum.
This supports a production pipeline that generates corneal maps from device exports. Different maps often share the same rings and angles while carrying different measurements. Rebuilding their geometry each time repeats work that the changing values do not affect.
The engineering problem is deciding exactly what can be reused. Missing measurements change the point set; output resolution changes the tables; and retaining many tables consumes memory. I made those dependencies explicit in the cache design.
One triangle, two maps
Take three source positions, A = (0, 0), B = (1, 0), and C = (0, 1). We want a value at q = (0.25, 0.25).
The point has barycentric weights 0.5, 0.25, and 0.25 relative to those vertices. If the source values are 10, 20, and 40, interpolation gives:
0.5 × 10 + 0.25 × 20 + 0.25 × 40 = 20Now change the values to 20, 30, and 50:
0.5 × 20 + 0.25 × 30 + 0.25 × 50 = 30The weights did not need to change. We changed the numbers attached to the vertices, not the vertices or the query point. We could have determined those weights before either map arrived.
Try changing the values below. Then remove vertex C. The second operation changes which observations are available, so the previously valid three-point calculation can no longer be used as it stands.
0.50 × 10 + 0.25 × 20 + 0.25 × 40 = 20.00
The weights stay fixed while measurements change. Removing a valid vertex changes the interpolation support, so this cached result can no longer be used.
Store the reusable operator
For a complete image, repeat the geometric construction for every output location. Each valid output point needs three source indices and three weights. Together, they define an operation:
y = W xHere x contains the source values and y the interpolated outputs. W describes the geometry. For linear interpolation within a two-dimensional triangulation, each valid row refers to at most three source vertices.
There is no reason to allocate a large dense matrix full of zeros. Store the three indices and weights per output point instead:
# Built once for a compatible geometry:
# indices: (n_outputs, 3), integer source indices
# weights: (n_outputs, 3), interpolation weights
# valid: (n_outputs,), points admitted by the geometry
# Applied to each new map:
# values: (n_sources,)
result = np.full(indices.shape[0], np.nan, dtype=np.float64)
selected = values[indices[valid]]
result[valid] = (selected * weights[valid]).sum(axis=1)This snippet is only the application stage. It assumes the tables are already built and compatible with the source ordering. It does not implement triangulation, periodic boundaries, missing-source policies, or handling for a degenerate point set.
SciPy already permits passing a precomputed triangulation to LinearNDInterpolator. With fixed output locations, we can reuse more: the triangle selected for each output and its barycentric weights. SciPy exposes the required triangle lookup and coordinate transforms in its Delaunay API.
The resulting hot path is a gather, multiplication, and reduction: triangle construction, lookup, and weight calculation have moved out of repeated map processing.
Reuse changes the cost model
Let building the geometric tables cost C, and applying them to one map cost A. For K maps, rebuilding every time costs:
K(C + A)Building once and reusing the tables costs approximately:
C + KAThe approximation leaves out lookup and cache-management overhead. It also assumes the same geometry really can be reused.
For a teaching example, choose C = 10, A = 1, and K = 10 in arbitrary time units. Rebuilding costs 110; reuse costs 20. These are constructed values, not measurements from CorneaForge.
For one map, both expressions give 11. The initial construction has not disappeared. A benchmark that times only the second call measures a warm-cache application, which may be appropriate for a repeated workload but says little about startup or a stream of unseen geometries.
Missingness is an input to the geometry
A cache is correct only while the facts behind its result remain true.
If a source value becomes invalid and the interpolation policy excludes invalid points, the point set changes. The old table might refer to a removed point. Replacing its value with zero does not fix the table; it invents a source measurement.
In our three-point example, removing C leaves no triangle spanning the query. In a larger map, removing a point may produce a different triangulation from the remaining points. That new triangulation may bridge a gap, so an additional support-quality policy can still be necessary. Rebuilding geometry and deciding that the resulting estimate is scientifically acceptable are separate checks.
In CorneaForge, I use a fixed angular grid. The radial shape and target resolution determine the coordinate grids. For maps with gaps, the cache key also includes the exact mask of non-finite source values. Two maps reuse an entry when their valid points occupy the same locations in the same ordering.
That key would be inadequate for a general library accepting arbitrary coordinates. Two arrays can have the same shape and validity mask while representing different positions. A general cache must identify the source and target coordinates, their ordering, and all settings that affect the operator. Coordinate rescaling or a changed boundary rule can matter even when array dimensions stay unchanged.
“Same shape” is a useful implementation condition. “Same mathematical problem” is the condition we actually need.
A cache has a memory budget
Repeated validity masks make reuse possible, but the number of masks is not universally small. A workload with mostly unique patterns can spend memory retaining tables that will never be requested again.
I bounded the missing-value cache to 128 entries with a least-recently-used policy. When it reaches its entry limit, it removes the entry that has gone unused for the longest time. This controls entry count; the memory cost still depends strongly on output resolution.
We can calculate that cost without timing anything. Three 32-bit indices and three 64-bit weights require 12 + 24 = 36 bytes per output location. Add two one-byte masks and the table storage is about 38 bytes per location. At 512 × 512, that is roughly 9.50 MiB per entry. Keeping 128 such entries retains about 1.19 GiB for those arrays alone, before object overhead and temporary working arrays.
A saved benchmark from 31 August 2026 confirms that allocation: 128 generated missing-value patterns at 512 × 512 retained exactly 1,275,068,416 bytes, or 1.19 GiB, in cached arrays. This was a controlled software fixture, separate from patient data. It gives the cache limit a concrete memory cost; a smaller machine may need fewer entries.
Verify the operator and measure the workload
The software cache skips repeated calculations. That is distinct from a processor’s L1 cache, which keeps recently accessed memory close to a CPU core. The optimization here follows from reusing the geometric operator; explaining a gain through a particular hardware cache would require separate profiling.
For numerical verification, compare cached and uncached interpolation using identical coordinates, ordering, missing-value rules, domain masks and precision. A planar function should be reproduced inside admitted triangles. A curved surface tests agreement between the two implementations, rather than exact recovery of the generating surface. Points outside the triangulation and changes to the validity mask deserve explicit cases.
For performance, distinguish table construction, warm-cache application, and end-to-end processing of a representative map sequence. Resolution, mask reuse and memory use determine whether the cache helps. The retained benchmark above measures cold construction and allocation under distinct masks. A matched sequence containing cache hits and misses would quantify the application speedup.
Put the optimization back into the pipeline
Suppose interpolation accounts for 20% of a pipeline and becomes ten times faster. In a simplified model, the new normalized runtime is:
0.80 + 0.20 / 10 = 0.82The overall speedup is 1 / 0.82, about 1.22×. The calculation is illustrative. It applies the familiar limit associated with Amdahl’s work: improving one component leaves the remaining work to be done.
I look for this separation of fixed structure and changing values before reaching for a different language or accelerator. It appears in repeated image warps, projection operators, and sensor geometry. The opportunity is worth testing when the operator is expensive, reusable, and affordable to retain. If its inputs change constantly, rebuilding may be the simpler and better choice.
Check your cache rule
The values change; coordinates, ordering, validity, and interpolation rules stay fixed. Rebuild? No. The same operator can be applied to the new values.
The array shape stays fixed, but one source position moves. Rebuild? Yes, unless you establish that the operator is unchanged. Matching dimensions does not establish that.
A warm-cache microbenchmark is ten times faster. Is the application ten times faster? Measure the application. Construction, misses, other computations, and I/O still count.