Interpolation is a very important subject in computing and comes up regularly in mathematics. The shortest distance between 2 points is a straight line. Interpolation gives the value of points along this line. I use the term point lightly, it could be a point on the number line or a two dimensional coordinate or any number system defining addition, subtraction and multiplication by a scaler (ie .5 for 50%)

      Source code is presented below. It implements Nth dimensional interpolation on arrays essentially giving them floating point indexing. And is presented in C Sharp and Python.

      Given points A and B where A is the start and B is the end of the line and using a variable ‘q’ to represent the position on the line from 0 (A) through 1 (B). The fundamental formula for linear interpolation is ‘A + (B – A) * q’.

      The dimensionality I speak of is an extension of bi-linear interpolation (where q is of two dimensions). Usually these dimensions are ‘x’ and ‘y’. And there are 4 points in the interpolation (top-left,top-right,bottom-left, and bottom right). Lets start with y.

      first, interpolate from top to bottom for each side.

      lets call top-left C, top-right D, bottom-left E and bottom-right F .

            Left = C + (E – C) * Qy
            Right = D + (F – D) * Qy

            Middle = Left + (Right – Left) * Qx

            The middle point is the final result.
      This technique of working through the dimensions continues for any number of dimensions. The algorithm below accomplishes this through folding.

      Using 2 dimensions, there are 4 points. So listing the values in the order of higher dimensions you get a list containing C,D,E,F.
To interpolate across multiple dimensions, the list is recursively halved and linear interpolation is performed between the halves until a final result is obtained.

            G = C + (E-C) *Q[0]             H = D + (F-D) *Q[0]

      The list is then G,H. Continuing, once again combining halves, it comes to G + (H-G) * Q[1] (each step, the pertinent dimension increases by one while the list is halved to the final, single, element.

      Using 3 dimensions there are 8 points (or 2 to the 3’rd power) the list would be

            V[0,0,0], V[0,0,1], V[0, 1, 0], V[0, 1, 1], V[1,0,0], V[1,0,1], V[1, 1, 0], V[1, 1, 1] .

      Higher dimensions have a lower index, so V[0,*,*]->V[1,*,*] becomes V[0,0,*],V[0,1,*],V[1,0,*],V[1,1,*] etc and these are the initial values of the ‘flat’ array in the code.

      Note that with sizes of a dimension larger than two, the algorithm works by interpolating off the coordinate’s whole number portion towards the next.

      The most common use of bi linear interpolation is computer graphics where an image looks far less digital if interpolated (taking rough edges off). In this case C,D,E and F would be color values. Which are then interpolated across color components (channels). All the red components interpolated to the results red component, etc.

      My favorite use of Nth dimensional interpolation is color space transforms in rgb. The rgb components make for 3 dimensions and the transform adds a 4th dimension. I then usually initialize with the actual color values so that [0,0,0,0] is initialized to black and [0,255,255,255] is initialized to white so the lower dimensions of the index are the color value. Then I initialize [1,*,*,*] to smoothed colors from random. Then when interpolated [0,r,g,b] is not transformed [1,r,g,b] is completely transformed, and [.5,r,g,b] is half way in between. I have achieved some awesome metallic transforms with this scheme. Note, this can use quite a bit of memory and might not be viable for 32bit platforms.

public class Nth
{
    public delegate T InterpolateDelegate(T a, T b, double q);

    public double InterpolateDouble(double a, double b, double q)
    {
        return a + (b - a) * q;
    }

    public static T Interpolate(System.Array inputArray, double[] coords, InterpolateDelegate interpol)
    {
        int dimension;
        int numDimensions = coords.Length;

        if (inputArray.Rank != numDimensions)
            throw new System.ArgumentException("inputArray and coords must have the same number of dimensions");

        int stackHeight = 1 << numDimensions;

        T[] flat = new T[stackHeight];

        int[] baseCoords = new int[numDimensions];
        int[] interpCoords = new int[numDimensions];

        double[] _q = new double[numDimensions];
        for (dimension = 0; dimension < numDimensions; dimension++)
        {
            baseCoords[dimension] = (int) coords[dimension];
            _q[dimension] = coords[dimension] - baseCoords[dimension];
        }

        int dim = 0;

        for (int i = 0; i < stackHeight; i++)
        {
            int ii = i;

            dim = 0;
            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) inputArray.GetValue(interpCoords);
        }

        int foldedStackHeight = stackHeight;
        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];
    }
}

# ChatGPT helped me convert from C Sharp!
# see https://johncornell.net/nth-dimensional-interpolation-on-array/

from numpy import array
from typing import List


class Nth:

    def general_interpolate(a, b, q):
        return a + (b - a) * q

    def interpolate(input_array: array, coords: List[float], interpol: callable):
        num_dimensions = len(coords)

        if input_array.ndim != num_dimensions:
            raise ValueError("input_array and coords must have the same number of dimensions")

        stack_height = 1 << num_dimensions flat = [None] * stack_height base_coords = [0] * num_dimensions interp_coords = [0] * num_dimensions _q = [0.0] * num_dimensions for dimension in range(num_dimensions): base_coords[dimension] = int(coords[dimension]) _q[dimension] = coords[dimension] - base_coords[dimension] for i in range(stack_height): ii = i for dimension in range(num_dimensions): p = base_coords[dimension] + (ii % 2) if p >= input_array.shape[dimension]: p = input_array.shape[dimension] - 1

                interp_coords[dimension] = p

                ii >>= 1

            flat[i] = input_array[tuple(interp_coords)]

        folded_stack_height = stack_height

        dim = num_dimensions - 1;

        while folded_stack_height != 1:
            folded_stack_height >>= 1
            for position in range(folded_stack_height):
                flat[position] = interpol(flat[position], flat[position + folded_stack_height], _q[dim])
                flat[position + folded_stack_height] = None

            dim -= 1

        return flat[0]