Trilinear Interpolation

This section describes the trilinear interpolation algorithm using pseudo-Fortran code.

First, number the eight vertices of a hexahedron.

   x1 = x(i  ,j  ,k  )
   x2 = x(i+1,j  ,k  )
   x3 = x(i  ,j+1,k  )
   x4 = x(i+1,j+1,k  )
   x5 = x(i  ,j  ,k+1)
   x6 = x(i+1,j  ,k+1)
   x7 = x(i  ,j+1,k+1)
   x8 = x(i+1,j+1,k+1)

Defined as below, notice that xa = x1 when a = −1 and that xa = x2 when a = 1 and similar observations for the rest of the equations.

   xa = [(1-a)*x1 + (1+a)*x2] / 2
   xb = [(1-a)*x3 + (1+a)*x4] / 2
   xc = [(1-a)*x5 + (1+a)*x6] / 2
   xd = [(1-a)*x7 + (1+a)*x8] / 2
   xe = [(1-b)*xa + (1+b)*xb] / 2
   xf = [(1-b)*xc + (1+b)*xd] / 2
   x  = [(1-g)*xe + (1+g)*xf] / 2

Take three steps of substitutions.

   x = {(1-g) * [(1-b)*xa + (1+b)*xb] +
        (1+g) * [(1-b)*xc + (1+b)*xd]} / 4

   x = {(1-g) * [(1-b) * ((1-a)*x1 + (1+a)*x2)  +
                 (1+b) * ((1-a)*x3 + (1+a)*x4)] +
        (1+g) * [(1-b) * ((1-a)*x5 + (1+a)*x6)  +
                 (1+b) * ((1-a)*x7 + (1+a)*x8)]} / 8

   x = {(1-g)*(1-b)*(1-a)*x1 + (1-g)*(1-b)*(1+a)*x2 +
        (1-g)*(1+b)*(1-a)*x3 + (1-g)*(1+b)*(1+a)*x4 +
        (1+g)*(1-b)*(1-a)*x5 + (1+g)*(1-b)*(1+a)*x6 +
        (1+g)*(1+b)*(1-a)*x7 + (1+g)*(1+b)*(1+a)*x8} / 8

Regrouping

   x =       (x8 + x7 + x6 + x5 + x4 + x3 + x2 + x1) / 8 +
           a*(x8 - x7 + x6 - x5 + x4 - x3 + x2 - x1) / 8 +
           b*(x8 + x7 - x6 - x5 + x4 + x3 - x2 - x1) / 8 +
           g*(x8 + x7 + x6 + x5 - x4 - x3 - x2 - x1) / 8 +
         a*b*(x8 - x7 - x6 + x5 + x4 - x3 - x2 + x1) / 8 +
         a*g*(x8 - x7 + x6 - x5 - x4 + x3 - x2 + x1) / 8 +
         b*g*(x8 + x7 - x6 - x5 - x4 - x3 + x2 + x1) / 8 +
       a*b*g*(x8 - x7 - x6 + x5 - x4 + x3 + x2 - x1) / 8

Define

   f0 = (x8 + x7 + x6 + x5 + x4 + x3 + x2 + x1) / 8 - x
   f1 = (x8 - x7 + x6 - x5 + x4 - x3 + x2 - x1) / 8
   f2 = (x8 + x7 - x6 - x5 + x4 + x3 - x2 - x1) / 8
   f3 = (x8 + x7 + x6 + x5 - x4 - x3 - x2 - x1) / 8
   f4 = (x8 - x7 - x6 + x5 + x4 - x3 - x2 + x1) / 8
   f5 = (x8 - x7 + x6 - x5 - x4 + x3 - x2 + x1) / 8
   f6 = (x8 + x7 - x6 - x5 - x4 - x3 + x2 + x1) / 8
   f7 = (x8 - x7 - x6 + x5 - x4 + x3 + x2 - x1) / 8

Then define

   f = f0     + f1*a   + f2*b   + f3*g     +
       f4*a*b + f5*a*g + f6*b*g + f7*a*b*g

Exactly similar relations can be derived for y and z.

   g = g0     + g1*a   + g2*b   + g3*g     +
       g4*a*b + g5*a*g + g6*b*g + g7*a*b*g

   h = h0     + h1*a   + h2*b   + h3*g     +
       h4*a*b + h5*a*g + h6*b*g + h7*a*b*g

   f = 0
   g = 0
   h = 0

So the three equations can be solved for the unknowns a, b and g. A solution can be obtained from Newton's method by letting

  df = -f
  dg = -g
  dh = -h

   df = f1*da     + f2*db     + f3*dg     +
        f4*b*da   + f4*a*db
        f5*g*da   + f5*a*dg
        f6*g*db   + f6*b*dg
        f7*b*g*da + f7*a*g*db + f7*a*b*dg

Regrouping

   df = (f1 + f4*b + f5*g + f7*b*g)*da +
        (f2 + f4*a + f6*g + f7*a*g)*db +
        (f3 + f5*a + f6*b + f7*a*b)*dg
      = -f

Similarly

   dg = (g1 + g4*b + g5*g + g7*b*g)*da +
        (g2 + g4*a + g6*g + g7*a*g)*db +
        (g3 + g5*a + g6*b + g7*a*b)*dg
      = -g

   dh = (h1 + h4*b + h5*g + h7*b*g)*da +
        (h2 + h4*a + h6*g + h7*a*g)*db +
        (h3 + h5*a + h6*b + h7*a*b)*dg
      = -h

Using an initial guess of a = b = g = 0 the three equations can be solved for da, db and dg. Then one has

   anew = aold + da
   bnew = bold + db
   gnew = gold + dg

If the matrix of the coefficients are zero it is because the hexahedron is degenerate. The new a, b and g can be used to compute new coefficients and solve for new deltas and iteration should continue util the deltas are negligible. The interpolating program has a hard-wired limit of 20 iterations. If after 20 iterations or if the absolute values of any of the three parameters exceed 5.0, then the method fails and the starting stencil is returned as the "best" one.

After a, b and g have been computed the interplated value of any variable known at the vertices can be computed by

   v  = [(1-g)*(1-b)*(1-a)*v1 + (1-g)*(1-b)*(1+a)*v2 +
         (1-g)*(1+b)*(1-a)*v3 + (1-g)*(1+b)*(1+a)*v4 +
         (1+g)*(1-b)*(1-a)*v5 + (1+g)*(1-b)*(1+a)*v6 +
         (1+g)*(1+b)*(1-a)*v7 + (1+g)*(1+b)*(1+a)*v8] / 8

This equation was copied from the equation for X and can easily be verified by letting a, b and g take on their extreme values of ±1. If (x,y,z) is inside the hexahedron, then the parameters a, b and g will be between ±1. If one of them is outside those limits then the given point is outside the hexahedron and the above formula becomes an extrapolation, not an interpolation.