Does anyone know how to do an "inverse" trilinear interpolation?

Question

Trilinear interpolation approximates the value of a point (x, y, z) inside a cube using the values at the cube vertices. I´m trying to do an "inverse" trilinear interpolation. Knowing the values at the cube vertices and the value attached to a point how can I find (x, y, z)? Any help would be highly appreciated. Thank you!

Answer 1

I'm not sure you can for all cases. For example using tri-linear filtering for colours where each colour (C) at each point is identical means that wherever you interpolate to you will still get the colour C returned. In this situation ANY x,y,z could be valid. As such it would be impossible to say for definite what the initial interpolation values were.

I'm sure for some cases you can reverse the maths but, i imagine, there are far too many cases where this is impossible to do without knowing more of the input information.

Good luck, I hope someone will prove me wrong :)

Answer 2

The problem as you're describing it somewhat ill-defined.
What you're asking for basically translates to this: I have a 3D function and I know its values in 8 known points. I'd like to know what is the point in which the function received value V.
The trouble is that in most likelihood there is an infinite number of such points which make a set of surfaces, lines or points, depending on the data.
One way to find this set is to use an iso-surfacing algorithm like Marching cubes.

Answer 3

The wikipedia page for trilinear interpolation has link to a NASA page which allegedly describes the inversing process - have you had a look at that?

Answer 4

You are solving for 3 unknowns given 1 piece of data, and as you are using a linear interpolation your answer will typically be a plane (2 free variables). Depending on the cube there may be no solutions or a 3D solution space.

I would do the following. Let v be the initial value. For each "edge" of the 12 edges (pair of adjacent vertices) of the cube look to see if 1 vertex is >=v and the other <=v - call this an edge that crosses v.

If no edges cross v, then there are no possible solutions.

Otherwise, for each edge that crosses v, if both vertices for the edge equal v, then the whole edge is a solution. Otherwise, linearly interpolate on the edge to find the point that has a value of v. So suppose the edge is (x1, y1, z1)->v1 <= v <= (x2, y2, z2)->v2.

s = (v-v1)/(v2-v1)
(x,y,z) = (s*(x2-x1)+x1, (s*(y2-y1)+y1, s*(z2-z1)+z1)

This will give you all edge points that are equal to v. This is a solution, but possibly you want an internal solution - be aware that if there is an internal solution there will always be an edge solution.

If you want an internal solution then just take any point linearly between the edge solutions - as you are linearly interpolating then the result will also be v.

Answer 5

Let's start with 2d: think of a bilinear hill over a square km, with heights say 0 10 20 30 at the 4 corners and a horizontal plane cutting the hill at height z. Draw a line from the 0 corner to the 30 corner (whether adjacent or diagonal). The plane must cut this line, for any z, so all points x,y,z fall on this one line, right ? Hmm.

OK, there are many solutions -- any z plane cuts the hill in a contour curve. Say we want solutions to be spread out over the whole hill, i.e. minimize two things at once:

• vertical distance z - bilin(x,y),
• distance from x,y to some point in the square.

Scipy.optimize.leastsq is one way of doing this, sample code below; trilinear is similar.

(Optimizing any two things at once requires an arbitrary tradeoff or weighting: food vs. money, work vs. play ... Cf. Bounded rationality )

""" find x,y so bilin(x,y) ~ z  and  x,y near the middle """
from __future__ import division
import numpy as np
from scipy.optimize import leastsq

zmax = 30
corners = [ 0, 10, 20, zmax ]
midweight = 10

def bilin( x, y ):
    """ bilinear interpolate
        in: corners at 0 0  0 1  1 0  1 1 in that order (binary)
        see wikipedia Bilinear_interpolation ff.
    """
    z00,z01,z10,z11 = corners  # 0 .. 1
    return (z00 * (1-x) * (1-y)
        + z01 * (1-x) * y
        + z10 * x * (1-y)
        + z11 * x * y)

vecs = np.array([ (x, y) for x in (.25, .5, .75) for y in (.25, .5, .75) ])
def nearvec( x, vecs ):
    """ -> (min, nearest vec) """
    t = (np.inf,)
    for v in vecs:
        n = np.linalg.norm( x - v )
        if n < t[0]:  t = (n, v)
    return t

def lsqmin( xy ):  # z, corners
    x,y = xy
    near = nearvec( np.array(xy), vecs )[0] * midweight
    return (z - bilin( x, y ), near )
        # i.e. find x,y so both bilin(x,y) ~ z  and  x,y near a point in vecs

#...............................................................................
if __name__ == "__main__":
    import sys
    ftol = .1
    maxfev = 10
    exec "\n".join( sys.argv[1:] )  # ftol= ...
    x0 = np.array(( .5, .5 ))
    sumdiff = 0
    for z in range(zmax+1):
        xetc = leastsq( lsqmin, x0, ftol=ftol, maxfev=maxfev, full_output=1 )
            # (x, {cov_x, infodict, mesg}, ier)
        x,y = xetc[0]  # may be < 0 or > 1
        diff = bilin( x, y ) - z
        sumdiff += abs(diff)
        print  "%.2g  %8.2g  %5.2g %5.2g" % (z, diff, x, y)

print "ftol %.2g maxfev %d midweight %.2g  =>  av diff %.2g" % (
    ftol, maxfev, midweight, sumdiff/zmax)