N’th-Dimensional Interpolation Revisited: When a Point Becomes a Sample

When I first wrote about N’th-dimensional interpolation on an array, the goal was straightforward:
given a coordinate with fractional parts, find the neighboring integer coordinates and interpolate
between them. In two dimensions this becomes bilinear interpolation. In three dimensions it becomes
trilinear interpolation. In N dimensions, the same idea generalizes naturally.

The core pattern is simple:

coordinate
    → integer base coordinate
    → fractional offset per dimension
    → collect 2^n neighboring corner values
    → fold those values through interpolation
    → final interpolated value

That original version worked by treating each corner as a point. For each generated corner coordinate,
the array was sampled directly:

value = inputArray[cornerCoordinate]

In C# terms, the key line was essentially:

flat[i] = inputArray.GetValue(interpCoords);

That line was correct, but it also hid something important. It made a quiet assumption:
a coordinate sample means one array cell.

The new realization is that this does not have to be true.

The Sampling Seam

The important change is replacing direct array access with a sample delegate:

flat[i] = sample(inputArray, interpCoords);

This small change opens the algorithm up. Each corner no longer has to mean “read one value from
the array.” Each corner can now mean “sample this location according to some rule.”

The original behavior still exists as the default sample:

SampleDefault(array, coords):
    return array[coords]

But once sampling is abstracted, the interpolation algorithm becomes more than a point interpolator.
It becomes a framework where the meaning of a sample can be changed.

PointSample     → read one cell
BoxSample       → average a local rectangle
CubeSample      → average a local cube
HyperBoxSample  → average a local N-dimensional region
KernelSample    → use weighted neighborhood logic

This is the conceptual upgrade:

Interpolation does not have to interpolate points.
It can interpolate samples.

The Half-Pixel Thought

The discovery that led me back to this code was the idea of a half-pixel.
In two-dimensional image terms, a half-pixel location can be thought of as the space between neighboring
pixels. At exactly halfway between four pixels, ordinary bilinear interpolation naturally averages the
surrounding 2×2 rectangle.

That gives a helpful way to think about the coordinate:

(x, y)       → sample at the pixel/cell position
(x+.5, y+.5) → sample halfway into the neighboring rectangle

In N dimensions, the same idea generalizes:

coords[d] + 0.5

This shifts the sampling coordinate by half a cell in each dimension. Then the existing interpolation
logic does what it already knows how to do: it finds the surrounding 2^n corners and folds them down
into a final value.

In 2D, this means the half-shift samples across a rectangle.
In 3D, it samples across a cube.
In N dimensions, it samples across a hyper-rectangle.

Point Sample vs. Region Sample

There are now two related but distinct ideas:

1. Shift the coordinate
   coords → coords + 0.5

2. Change the sample meaning
   point sample → region/kernel sample

The coordinate shift changes where interpolation happens.
The sample delegate changes what each corner means.

Together, they create a very flexible structure:

coordinate transform
    → corner generation
    → sample delegate
    → interpolation fold

That separation matters. The interpolator does not need to know whether a sample is a point,
a rectangle, a cube, or a weighted neighborhood. It only needs a value for each corner.
The sample function owns the meaning of that value.

Pseudo-Code: The Sample Delegate

The delegate idea can be expressed like this:

SampleDelegate(array, coords):
    return some value of type T from the array at or around coords

The original point sample:

PointSample(array, coords):
    return array[coords]

A simple 2D box sample might look like:

BoxSample2D(array, coords):
    sum = 0
    count = 0

    for dy in 0..1:
        for dx in 0..1:
            p = clamp(coords + (dx, dy))
            sum += array[p]
            count += 1

    return sum / count

A 3D cube sample follows the same pattern:

CubeSample3D(array, coords):
    sum = 0
    count = 0

    for dz in 0..1:
        for dy in 0..1:
            for dx in 0..1:
                p = clamp(coords + (dx, dy, dz))
                sum += array[p]
                count += 1

    return sum / count

And the N-dimensional version is the natural continuation:

HyperBoxSampleND(array, coords):
    sum = 0
    count = 0

    for each offset in all binary offsets for N dimensions:
        p = clamp(coords + offset)
        sum += array[p]
        count += 1

    return sum / count

For N dimensions, the number of offsets in a 2-wide hyper-box is:

2^n

That mirrors the interpolation itself, which also gathers 2^n neighboring corner values.
This symmetry is part of what makes the idea feel natural.

The Updated Interpolation Shape

With the sampling delegate in place, the high-level interpolation algorithm becomes:

Interpolate(array, coords, interpolator, sample):
    split coords into base coordinates and fractional q values

    for each corner among 2^n corners:
        cornerCoord = baseCoord + cornerOffset
        flat[i] = sample(array, cornerCoord)

    while more than one value remains:
        fold values together using q for the current dimension

    return final folded value

The old version is still available by passing the default point sample.
The new version allows richer sampling without rewriting the interpolation fold.

Place for Updated Code

Below is the updated C# implementation.

    public class Nth
    {
        public delegate T SampleDelegate<T>(System.Array inputArray, int[] coords);
        public static T SampleDefault<T>(Array inputArray, int[] coords)
        {
            return (T)inputArray.GetValue(coords);
        }
        
        public delegate T InterpolateDelegate<T>(T a, T b, double q);
        public static double InterpolateDouble(double a, double b, double q)
        {
            return a + (b - a) * q;
        }

        //if you like param arrays here is a nice convenience wrapper
        public static T Interpolate<T>(System.Array inputArray, int[] coords, InterpolateDelegate<T> interpol,
            bool half = false, SampleDelegate<T>? sample = null)
        {
            double[] newCoords = new double[coords.Length];
            for (int i=0; i<coords.Length; i++)
            {
                newCoords[i] = coords[i];
            }
            
            return Interpolate(inputArray, newCoords, interpol, half, sample);
        }

        public static T Interpolate_HalfUnit<T>(System.Array inputArray, int[] coords, InterpolateDelegate<T> interpol,
            SampleDelegate<T>? sample = null)
        {
            return Interpolate(inputArray, coords, interpol, true, sample);
        }
        
        public static T Interpolate<T>(System.Array inputArray, double[] coords, InterpolateDelegate<T> interpol, bool half = false, SampleDelegate<T>? sample = null)
        {
            int dimension;
            int numDimensions = coords.Length;

            if (inputArray.Rank != numDimensions)
                throw new System.ArgumentException("inputArray and coords must have the same number of dimensions");
            
            if (sample == null) 
                sample = SampleDefault<T>;
            
            int stackHeight = 1 << numDimensions;
            
            T[] flat = new T[stackHeight];

            int[] baseCoords = new int[numDimensions];
            int[] interpCoords = new int[numDimensions];
            
            double[] _q = new double[numDimensions];
            if (!half)
            {
                for (dimension = 0; dimension < numDimensions; dimension++)
                {
                    baseCoords[dimension] = (int)Math.Floor(coords[dimension]);
                    _q[dimension] = coords[dimension] - baseCoords[dimension];
                }
            }
            else
            {
                for (dimension = 0; dimension < numDimensions; dimension++)
                {
                    double shifted = coords[dimension] + 0.5;
                    if (shifted >= inputArray.GetLength(dimension))
                        shifted = inputArray.GetLength(dimension) - 1;

                    baseCoords[dimension] = (int) Math.Floor(shifted);
                    _q[dimension] = shifted - baseCoords[dimension];
                }
            }
            
            for (int i = 0; i < stackHeight; i++)
            {
                int ii = i;
                
                for (dimension = 0; dimension < numDimensions; dimension++)
                {
                    int p =  baseCoords[dimension] + (ii % 2);
                    if (p >= inputArray.GetLength(dimension)) p = inputArray.GetLength(dimension) - 1;

                    interpCoords[dimension] = p;

                    ii >>= 1;
                }

                flat[i] = (T) sample(inputArray, interpCoords);
            }

            int foldedStackHeight = stackHeight;
            int dim = numDimensions-1;

            while (foldedStackHeight != 1)
            {
                foldedStackHeight >>= 1;
                for (int position = 0; position < foldedStackHeight; position++)
                {
                    flat[position] = interpol(flat[position], flat[position + foldedStackHeight], _q[dim]);
                    flat[position + foldedStackHeight] = default(T);
                }

                dim--;
            }

            return flat[0];
        }
    }

Why This Matters

The original algorithm answered the question:

How do I interpolate between neighboring points in an N-dimensional array?

The revised version asks a broader question:

What should a sample mean before interpolation happens?

That is a much more powerful question.

For image data, a sample might mean a pixel, a half-pixel blend, or a small rectangle.
For volume data, it might mean a voxel or a cube of voxels.
For procedural fields, it might mean a local kernel.
For generalized numerical arrays, it might mean a neighborhood summary.

The algorithm did not need to become complicated to support this.
It only needed one seam:

array.GetValue(coords)
    → sample(array, coords)

That is the moment where a point becomes a sample.

Closing Thought

This update is exciting to me because it shows the original N’th-dimensional interpolation routine
becoming more general without losing its original simplicity.

The interpolation fold still does the same elegant work:
it reduces 2^n corner values down to one final value.

But now those corner values can carry more meaning.
They can be raw points, half-shifted blends, local regions, or eventually weighted kernels.

In short:

Point → Sample → Region → Kernel

That is a small change in code, but a large change in what the algorithm can express.