Surface Fitting

In basic terms, suppose you can measure some quantity which depends on several variables--say you can measure the hardness of concrete, which depends on the amount of various materials added to it. You test the hardness for a few material amounts to get some sample data. You want to somehow figure out the hardness at points in between your sample data. This is where surface fitting comes in--you want to fit a surface to your data (possibly in more than 3 dimensions, depending on your number of variables).

Lower dimensional interpolation (fitting) is pretty common. For instance, VBAhack's Cubic Spline Tutorial gives an example of single-variable interpolation. This post extends a single-variable interpolation to generate a surface fitting function. It does nothing to minimize error--that's another type of surface fitting you'd use if you had a theoretical model that predicted the relationships between your measured quantity and your variables.

This type of surface fitting would be useful if you wanted a "smooth" or "reasonable" fit when you only have a 1-dimensional interpolation handy. To keep the description reasonably short, I've used a fair amount of notation and jargon. Sorry; I included pseudocode.

n-dimensional Surface Fitting through Successive Coordinate Line Interpolation on a Finite Lattice

Suppose you have a real function F of finitely many real variables, x[i], i = 1, ..., n. For instance, this function could describe the temperature at each point in a fish tank, using three position coordinates. As another example, F could also describe the hardness of concrete as various amounts of different materials are added to the mixture.

Further suppose that the value of F is known at the points of a finite lattice grid. In 2D, this lattice could be the points on a 2x3 grid given by {(0, -1), (0.5, -1), (0, 0), (0.5, 0), (0, 1), (0.5, 1)}. This example gives x[1] a uniform spacing of 0.5 and x[2] a uniform spacing of 1. In general, the lattice could use any spacing in each direction, and even vary the spacing between points, so long as a general rectangular grid is formed. In each direction x[i], make a partition with points x[i,1] < x[i,2] < ... < x[i,m[i]], using a total of m_i points in the ith coordinate direction. In the example above, m[1] = 2, m[2] = 3, x[1,1] = 0, x[1,2] = 0.5, x[2,1] = -1, x[2,2] = 0, and x[2,3] = 1. The overall lattice is made of points which take values for their ith coordinate from the partition {x[i, c], c from 1 to m[i]}.

Finally, suppose you have a way to do real 1D interpolation. That is, given a sequence of points a[1] < a[2] < ... < a[m] which forms a partition of [a[1], a[m]] and associated values at those points v[1], ..., v[m], you can generate a function f passing through each point, so that (a[i], v[i]) = (a[i], f(a[i])), where you can evaluate f at any point x between a[1] and a[m]. A linear interpolation is a simple example of such a function. Between a[k] and a[k+1], f(x) = (a[k+1] - x) / (a[k+1] - a[k]) * v[k] + (x - a[k]) / (a[k+1] - a[k]) * v[k+1].


We wish to find the value of F at an arbitrary point P = (p[1], p[2], ..., p[n]) inside the lattice, which will usually not be on the lattice. Extending the temperature example above, we could equip a fish tank with a 3D grid of temperature sensors, and try to use this information to generate a reasonable graph of the temperature of the tank at other points. First, find the nearest lattice points to P, which generate an n-dimensional cell with P inside. That is, find indecies c[1] through c[n] such that lattice points x[i,c[i]] <= p[i] < x[i,c[i+1]] for each i from 1 to n. If you've followed this so far, finding such points should be simple.

Now, calculate F(P) = F(p[1], p[2], ..., p[n]) through successive interpolation. Suppose we have an interpolation function f[1] which approximates f[1](x) ~= F(x, p[2], ..., p[n]). If we have this function, we can calculate F(P) as simply f[1](p[1]). So, our problem becomes to generate such a function. Do so by taking a[i] = x[1, i] for i = 1 to m[1], corresponding to points (x[1, i], p[2], ..., p[n]). However, we don't have values v[i] at these points yet. We'll obtain them by recursively. That is, to find the value of F at (x[1, 1], p[2], ..., p[n]), we'll make a new interpolation function f[2] which approximates f[2](x) ~= F(x[1, 1], x, p[3], ..., p[n]). This will be generated by a[i] = x[2, i] for i = 1 to m[2], with recursively obtained values. Eventually, we'll get rid of all p[i] terms, to be replaced by lattice coordinates x[i, j]. This is where the recursion ends.

I've written this algorithm in Python-esque pseudocode

Code:

---

'P is an array from 1 to n representing the point you want to evaluate
'F is a lookup table mapping points on the lattice to values of the function to be interpolated
'xij is a 2D array. xij[i] gives the lattice partition in the x_i direction
'1DInterpolation is a function which finds an interpolation which passes through (a[i], v[i]), and returns the value of that interpolation at x. a[] is a partition, and v[] gives associated values at each point of a[].
'k is an accumulator, and an optimization. It tracks which indecies are in lattice partitions, so they don't have to be checked each recursive call

'No reasonable-length comment could explain this in any sane detail. Sorry.
'To use:
'1. Implement a 1DInterpolation function. Natural cubic spline interpolation is a decent choice.
'2. Decide on a rectangular lattice for your function. Use these partitions to generate xij.
'3. Evaluate your function at points on the lattice, and put those values in the lookup table F{}.
'4. Pick a point P you wish to interpolate.
'5. Interpolate by calling nDInterpolation(P, F, xij, 1DInterpolation).
function nDInterpolation(P[], F{}, xij[][], 1DInterpolation(a[], v[], x), Optional k=0)
    'Base case; the value of P is known since all of P's coordinates are on the lattice
    'For speed: ensure len() and F{} are constant-time calls.
    if k = len(P) then return F{P}

    'Vary the k+1'th coordinate over that coordinate's lattice partition.
    as = xij[k+1]
    dim vs[1 to len(as)]

    'Interpolate along the k+1'th coordinate recursively
    Pk1 = P[k+1] 'optimization; allows P to be passed ByRef throughout, avoiding lengthy copies
    for i = 1 to len(vs)
        P[k+1] = as[i]
        vs[i] = nDInterpolation(P, F, xij, 1DInterpolation, k+1)
    next i
    P[k+1] = Pk1

    return 1DInterpolation(as, vs, P[k+1])
end function

---

Conceptual Example Evaluation:
F(x, y, z) = x + y*z
Lattice = [0, 1, 2] X [1, 2] X [2, 2.5, 4] (X = Cartesian Product)
implies xij = [[0, 1, 2], [1, 2], [2, 2.5, 4]]

Take P = [0.5, 1.5, 2.6]. Each nDInterpolation call needs a P[]. These are listed below. For example, the first nDInterpolation takes the initial P, [0.5, 1.5, 2.6]. The second call, done inside the for loop of the first call, changes the 1st coordinate of P to 0, leaving the rest unchanged. That call recurses, forming a tree of calls. In the end, every point of the lattice will be needed to generate this interpolation, since the 1DInterpolation function may have global behavior.

Supposing 1DInterpolation is a standard linear interpolation, this will generate a multi-dimensional linear interpolation over the entire function. However, since this type of interpolation has purely local behavior, this general method will be inefficient.


k=number of spaces in indent, arranged in order of evaluation:

Code:

---

[0.5, 1.5, 2.6]
 [0, 1.5, 2.6]
  [0, 1, 2.6]
  [0, 1, 2]
  [0, 1, 2.5]
  [0, 1, 4]
  [0, 2, 2.6]
  [0, 2, 2]
  [0, 2, 2.5]
  [0, 2, 4]
 [1, 1.5, 2.6]
  [1, 1, 2.6]
  [1, 1, 2]
  [1, 1, 2.5]
  [1, 1, 4]
  [1, 2, 2.6]
  [1, 2, 2]
  [1, 2, 2.5]
  [1, 2, 4]
 [2, 1.5, 2.6]
  [2, 1, 2.6]
  [2, 1, 2]
  [2, 1, 2.5]
  [2, 1, 4]
  [2, 2, 2.6]
  [2, 2, 2]
  [2, 2, 2.5]
  [2, 2, 4]

---