[转载]三线性插值
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 = &S722;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
where the g's and h's are defined exactly similarly to the f's. Note that in the three equations above everything is known except for a, b and g. From the expression for x, and the similar equations for y and z, it is seen that
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
where
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.
