Sample Positions, Trilinear Interpolation, and the Half-Cell Offset

A grid of samples is a set of numbers with positions attached, and the positions are where most of the ambiguity lives. Two files can describe the same physical box with the same dimensions and produce images offset by half a cell, mirrored, or transposed, because the header recorded dimensions but not the convention that relates an index to a coordinate. Everything downstream, interpolation, contouring, and volume rendering, consumes that relation, so a wrong convention is not a small error but a systematic one that survives to the final image.

The other half of the problem is what the numbers mean between the samples. A scalar field is defined everywhere; the grid holds a finite sampling of it. Any visualization that needs a value at a position that is not a sample must reconstruct one, and the reconstruction is a model with its own error, its own smoothness, and its own artefacts.

Node-centred and cell-centred data differ by half a spacing

The two conventions are both common and neither is wrong. In node-centred data, sometimes called point data, values sit on the vertices of the grid. A grid with $n_x$ points along an axis has $n_x - 1$ cells, and the first and last points lie on the domain boundary. In cell-centred data, values sit at the centres of the cells: $n_x$ cells, $n_x$ values, and the first value is half a spacing inside the boundary.

The relationship is

$$x_i = x_0 + i,\Delta x \quad\text{(node-centred)}, \qquad x_i = x_0 + \left(i + \tfrac{1}{2}\right)\Delta x \quad\text{(cell-centred)}.$$

A file format that stores only dimensions and spacing cannot distinguish the two, and a reader that assumes the wrong one shifts the field by $\Delta x/2$ relative to geometry loaded from the same directory. The shift is invisible in a lone volume render and becomes obvious the moment an isosurface is compared against a CAD model, a mesh, or a slice plane. Conventions are usually declared in the format’s metadata rather than inferred, and where a format has no such field, the convention is part of the format’s documentation and not a matter of taste.

Finite-volume codes often produce cell-centred output by construction: the conserved quantities are averages over the cell volume, and a pressure at a vertex would require an extra interpolation that changes the meaning of the number. Finite-element codes produce node-centred output because the basis functions are attached to vertices. Converting between the two is not a resampling problem with a free choice of filter; it is a restriction or prolongation operator, and both directions lose or smear information. Averaging the eight cells around a node is a restriction that acts as a low-pass filter, and evaluating at cell centres from nodes is a prolongation that cannot recover detail finer than the node spacing.

Trilinear reconstruction and its cost

Inside a cell, the trilinear interpolant is the tensor product of three one-dimensional linear interpolants. With local coordinates $u, v, w \in [0,1]$ inside the cell, the value is a weighted sum of the eight corner values with weights $(1-u)(1-v)(1-w)$, and so on. The interpolant is continuous across cell faces, because the face value depends only on the four corners shared by both cells, and its gradient is discontinuous across faces. That C0-only continuity matters for the images: iso-surfaces extracted from trilinear data have visible creases along cell boundaries unless the normals are computed from a smoothed field rather than from the interpolant’s analytic gradient.

The evaluation cost is seven linear interpolations, which in the common formulation is 7 multiply-add pairs along the leading dimensions, plus eight corner loads. On a 64-byte cache line and with 32-bit floats, a line holds sixteen values, so a stride-1 traversal along the fastest axis reads each line once for sixteen cells; a traversal along the slowest axis touches a new line for every step. This asymmetry, not the arithmetic, is why a naive trilinear sampler on a large volume runs far below the rate the floating-point unit allows.

/* Trilinear sample on a cell-centred grid, with bounds clamped to the
   outermost cell centres. nx, ny, nz are the number of samples per axis. */
static inline float sample_trilinear(const float *v, int nx, int ny, int nz,
                                     float fx, float fy, float fz)
{
    float gx = fx - 0.5f, gy = fy - 0.5f, gz = fz - 0.5f;  /* index space */
    if (gx < 0) gx = 0; if (gy < 0) gy = 0; if (gz < 0) gz = 0;
    if (gx > nx - 1) gx = nx - 1;
    if (gy > ny - 1) gy = ny - 1;
    if (gz > nz - 1) gz = nz - 1;

    int i = (int)gx, j = (int)gy, k = (int)gz;
    if (i > nx - 2) i = nx - 2;
    if (j > ny - 2) j = ny - 2;
    if (k > nz - 2) k = nz - 2;
    float u = gx - i, v = gy - j, w = gz - k;

    const float *p = v + ((size_t)k * ny + j) * nx + i;
    size_t sx = 1, sy = nx, sz = (size_t)nx * ny;

    float c00 = p[0]        + u * (p[sx]        - p[0]);
    float c10 = p[sy]       + u * (p[sy + sx]   - p[sy]);
    float c01 = p[sz]       + u * (p[sz + sx]   - p[sz]);
    float c11 = p[sz + sy]  + u * (p[sz + sy + sx] - p[sz + sy]);
    float c0  = c00 + v * (c10 - c00);
    float c1  = c01 + v * (c11 - c01);
    return c0 + w * (c1 - c0);
}

The clamp at the boundary is a policy choice. Clamping holds the edge value constant beyond the last sample, which is right for a field that continues past the domain but wrong for one that ends. Returning a no-data marker, or refusing to sample, is the alternative; silently returning zero creates a dark shell around every volume, a failure that is easy to mistake for attenuation.

Error scales with the square of the spacing

For a smooth field, trilinear interpolation has an error term proportional to $h^2 \max|\partial^2 f/\partial x_i^2|$ over the cell, with the constant differing by direction because the three axes contribute independently. Halving the spacing reduces the error by roughly a factor of four in smooth regions. Near a discontinuity, such as a shock in a compressible flow solution or a material interface, the error does not follow that law: the interpolant creates a linear ramp across the cell that contains the jump, so an iso-surface at a level inside the jump lands somewhere in that cell and its position is determined by the sampling, not by the physics. This is the mechanism behind the stepped surfaces that appear when a sharp interface is contoured at low resolution.

Higher-order reconstructions change the trade. Tricubic interpolation uses 64 samples for a cell and is C1, so the extracted surfaces have continuous normals. It can overshoot, producing values outside the range of the local samples, which matters when the quantity is physically bounded, such as a concentration that cannot be negative. Monotone variants such as the monotone cubic Hermite construction constrain the slopes to avoid introducing new extrema, at the cost of flattening genuine peaks. The choice is between fidelity to smooth data and the guarantee that the reconstruction does not invent values the data does not contain.

Resampling to a second grid is filtering, not copying

Resampling arises whenever data and display grids disagree: a simulation on a stretched, curvilinear, or unstructured mesh needs to be shown on a uniform image grid, or a coarse field needs to be shown at display resolution. The operation is a reconstruction of the continuous field followed by a resampling of it, and the quality is decided by the filter used for the reconstruction.

Nearest-neighbour resampling is a box reconstruction and produces blocky output with the original cells visible. Linear resampling is the trilinear interpolant above and produces smooth output but suppresses detail at the new grid’s Nyquist limit. When the target grid is coarser than the source, the correct operation is a low-pass filter followed by decimation; skipping the low-pass step aliases high-frequency content into low-frequency patterns, which in a scientific image can look like structure in the data. A common approximation is to average over the footprint of the target cell in source space, which is exact for a box filter and cheap for integer ratios of spacing.

Two-pass separable resampling is the practical structure: filter and resample along $x$, then along $y$, then along $z$, each pass one-dimensional. A separable filter of width $k$ along each axis costs $O(k)$ per output sample per pass, versus $O(k^3)$ for the equivalent non-separable kernel. The result is identical for a separable kernel, which covers box, linear, and Gaussian filters, and differs for kernels designed to be non-separable.

Memory is the other constraint. A $512^3$ grid of 32-bit floats is 512 MiB before any derived array, and a resampling pass that needs both source and destination resident doubles that. Algorithms that touch the volume slice by slice, keeping only a few slices of each array, are the usual answer, and they interact with the traversal order chosen for interpolation.

Where interpolation stops being the issue

Interpolation of a scalar field is a solved problem in the sense that the error is understood and bounded; the unresolved part is usually the model, not the mathematics. A field with a no-data mask must be handled before interpolation, because a sentinel value such as $-9999$ or a NaN propagates through the weighted sum and produces a value that is neither data nor missing. Interpolating across a hole invents data, and the honest alternatives are to mark the output as missing, to reduce the weights to the valid corners and renormalize, or to fill the hole explicitly and record that it was filled.

The next question, what to do with the interpolated field, splits into two mechanisms. Extracting surfaces at a chosen level from the reconstructed field is covered in Marching Cubes and the Case Table, which depends directly on the interpolation being trilinear inside each cell. Integrating paths through a vector field, where interpolation error changes the trajectory rather than a single value, is the subject of Streamline Integration in Vector Fields. Both consume the same reconstruction, and both inherit its convention errors.