In this article
Imagine three maps sampled on the same rings. Their colours differ because they contain different values. Their source coordinates are identical, and we want to display all three on the same output grid.
A general interpolation call can rebuild a triangulation, locate every output point, compute its interpolation weights, and combine the source values. Repeating that whole sequence for each map is valid, but several steps depend on inputs that have not changed.
This is a pattern I work with in CorneaForge’s polar-to-Cartesian conversion. The interesting part is the dependency analysis: which work belongs to the geometry, and which work belongs to the changing measurements? Once that boundary is clear, the code and its correctness conditions become easier to explain.
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. This does not establish a particular speedup, but it explains exactly which repeated calculations have been removed.
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.
The CorneaForge path reviewed for this article uses a fixed angular grid. Within that constrained domain, the radial shape and target resolution determine the coordinate grids. The cache for maps with gaps additionally identifies the exact mask of non-finite source values. Two maps can share that 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.
The implementation uses a bounded least-recently-used cache for these patterns. 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.
Those numbers describe the inspected array layout, not a general requirement of interpolation. Other representations can make different tradeoffs. They are enough to show why a cache limit copied from a large server deserves reconsideration on a smaller machine.
Software reuse and CPU caches answer different questions
There are two meanings of “cache” here. The software cache retains a computed operator so the program can skip work. A processor’s L1 cache retains recently accessed memory close to a CPU core.
The gather operation may repeatedly access a compact source array, which makes locality worth investigating. But the indices, weights, output, and temporary arrays occupy their own memory. A small input alone does not prove that the complete working set fits in L1.
Nor does it prove that L1 behavior explains an observed speedup. The geometric construction may simply have dominated the original runtime. Attributing a gain to hardware locality needs appropriate measurements, including the relevant memory-access behavior, alongside the elapsed times. The reuse argument is already useful without an unsupported explanation involving a particular CPU cache or instruction set.
Check the calculation before timing it
I would first compare the cached result with the intended reference using identical coordinates, source ordering, missing-value rules, domain masks, and numerical precision. A planar function provides a useful check: linear interpolation should reproduce it inside admitted triangles, within floating-point tolerance.
A curved synthetic surface tests a different property. The interpolant need not reproduce the generating surface exactly; the comparison is between the cached implementation and the same interpolation method without reuse. Otherwise, an algorithm change can be mistaken for a faster implementation of the old calculation.
Edges and gaps deserve their own cases. Outside points must not accidentally index the last triangle through a -1 lookup. Degenerate or insufficient point sets need a defined response. The cached tables should also be rejected when a source coordinate, target location, ordering, or validity pattern changes.
Only then would I report three timings: table construction, warm-cache application, and end-to-end processing of a representative sequence of maps. The report should include resolutions, pattern reuse, repetitions, hardware, numerical error, memory use, and a distribution of timings. The historical timing comments in this code are not a substitute for reproducing that benchmark.
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.