Inverse Bilinear Interpolation?

Question

I have four 2d points, p0 = (x0,y0), p1 = (x1,y1), etc. that form a quadrilateral. In my case, the quad is not rectangular, but it should at least be convex.

  p2 --- p3
  |      |
t |  p   |
  |      |
  p0 --- p1
     s

I'm using bilinear interpolation. S and T are within [0..1] and the interpolated point is given by:

bilerp(s,t) = t*(s*p3+(1-s)*p2) + (1-t)*(s*p1+(1-s)*p0)

Here's the problem.. I have a 2d point p that I know is inside the quad. I want to find the s,t that will give me that point when using bilinear interpolation.

Is there a simple formula to reverse the bilinear interpolation?

---

Thanks for the solutions. I posted my implementation of Naaff's solution as a wiki.

Answer 1

Well, if p is a 2D point, yes, you can get it easily. In that case, S is the fractional component of the total width of the quadrilateral at T, T is likewise the fractional component of the total height of the quadrilateral at S.

If, however, p is a scalar, it's not necessarily possible, because the bilinear interpolation function is not necessarily monolithic.

Answer 2

If all you have is a single value of p, such that p is between the minimum and maximum values at the four corners of the square, then no, it is not possible in general to find a SINGLE solution (s,t) such that the bilinear interpolant will give you that value.

In general, there will be an infinite number of solutions (s,t) inside the square. They will lie along a curved (hyperbolic) path through the square.

If your function is a vector valued one, so you have two known values at some unknown point in the square? Given known values of two parameters at each corner of the square, then a solution MAY exist, but there is no assurance of that. Remember that we can view this as two separate, independent problems. The solution to each of them will lie along a hyperbolic contour line through the square. If the pair of contours cross inside the square, then a solution exists. If they do not cross, then no solution exists.

You also ask if a simple formula exists to solve the problem. Sorry, but not really that I see. As I said, the curves are hyperbolic.

One solution is to switch to a different method of interpolation. So instead of bilinear, break the square into a pair of triangles. Within each triangle, we can now use truly linear interpolation. So now we can solve the linear system of 2 equations in 2 unknowns within each triangle. There may be one solution in each triangle, except for a rare degenerate case where the corresponding piecewise linear contour lines happen to be co-incident.

Answer 3

Some of the responses have slightly misinterpreted your question. ie. they are assuming you are given the value of an unknown interpolated function, rather than an interpolated position p(x,y) inside the quad that you want to find the (s,t) coordinates of. This is a simpler problem and there is guaranteed to be a solution that is the intersection of two straight lines through the quad.

One of the lines will cut through the segments p0p1 and p2p3, the other will cut through p0p2 and p1p3, similar to the axis-aligned case. These lines are uniquely defined by the position of p(x,y), and will obviously intersect at this point.

Considering just the line that cuts through p0p1 and p2p3, we can imagine a family of such lines, for each different s-value we choose, each cutting the quad at a different width. If we fix an s-value, we can find the two endpoints by setting t=0 and t=1.

So first assume the line has the form: y = a0*x + b0

Then we know two endpoints of this line, if we fix a given s value. They are:

(1-s)p0 + (s)p1

(1-s)p2 + (s)p3

Given these two end points, we can determine the family of lines by plugging them into the equation for the line, and solving for a0 and b0 as functions of s. Setting an s-value gives the formula for a specific line. All we need now is to figure out which line in this family hits our point p(x,y). Simply plugging the coordinates of p(x,y) into our line formula, we can then solve for the target value of s.

The corresponding approach can be done to find t as well.

Answer 4

Since you're working in 2D, your bilerp function is really 2 equations, 1 for x and 1 for y. They can be rewritten in the form:

x = t * s * A.x + t * B.x + s * C.x + D.x
y = t * s * A.y + t * B.y + s * C.y + D.y

Where:

A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0

Rewrite the first equation to get t in terms of s, substitute into the second, and solve for s.

Answer 5

I think it's easiest to think of your problem as an intersection problem: what is the parameter location (s,t) where the point p intersects the arbitrary 2D bilinear surface defined by p0, p1, p2 and p3.

The approach I'll take to solving this problem is similar to tspauld's suggestion.

Start with two equations in terms of x and y:

x = (1-s)*( (1-t)*x0 + t*x2 ) + s*( (1-t)*x1 + t*x3 )
y = (1-s)*( (1-t)*y0 + t*y2 ) + s*( (1-t)*y1 + t*y3 )

Solving for t:

t = ( (1-s)*(x0-x) + s*(x1-x) ) / ( (1-s)*(x0-x2) + s*(x1-x3) )
t = ( (1-s)*(y0-y) + s*(y1-y) ) / ( (1-s)*(y0-y2) + s*(y1-y3) )

We can now set these two equations equal to each other to eliminate t. Moving everything to the left-hand side and simplifying we get an equation of the form:

A*(1-s)^2 + B*2s(1-s) + C*s^2 = 0

Where:

A = (p0-p) X (p0-p2)
B = ( (p0-p) X (p1-p3) + (p1-p) X (p0-p2) ) / 2
C = (p1-p) X (p1-p3)

Note that I've used the operator X to denote the 2D cross product (e.g., p0 X p1 = x0*y1 - y0*x1). I've formatted this equation as a quadratic Bernstein polynomial as this makes things more elegant and is more numerically stable. The solutions to s are the roots of this equation. We can find the roots using the quadratic formula for Bernstein polynomials:

s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )

The quadratic formula gives two answers due to the +-. If you're only interested in solutions where p lies within the bilinear surface then you can discard any answer where s is not between 0 and 1. To find t, simply substitute s back into one of the two equations above where we solved for t in terms of s.

I should point out one important special case. If the denominator A - 2*B + C = 0 then your quadratic polynomial is actually linear. In this case, you must use a much simpler equation to find s:

s = A / (A-C)

This will give you exactly one solution, unless A-C = 0. If A = C then you have two cases: A=C=0 means all values for s contain p, otherwise no values for s contain p.

Answer 6

Here's my implementation of Naaff's solution, as a community wiki. Thanks again.

This is a C implementation, but should work on c++. It includes a fuzz testing function.

---

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

int equals( double a, double b, double tolerance )
{
    return ( a == b ) ||
      ( ( a <= ( b + tolerance ) ) &&
        ( a >= ( b - tolerance ) ) );
}

double cross2( double x0, double y0, double x1, double y1 )
{
    return x0*y1 - y0*x1;
}

int in_range( double val, double range_min, double range_max, double tol )
{
    return ((val+tol) >= range_min) && ((val-tol) <= range_max);
}

/* Returns number of solutions found.  If there is one valid solution, it will be put in s and t */
int inverseBilerp( double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double x, double y, double* sout, double* tout, double* s2out, double* t2out )
{
    int t_valid, t2_valid;

    double a  = cross2( x0-x, y0-y, x0-x2, y0-y2 );
    double b1 = cross2( x0-x, y0-y, x1-x3, y1-y3 );
    double b2 = cross2( x1-x, y1-y, x0-x2, y0-y2 );
    double c  = cross2( x1-x, y1-y, x1-x3, y1-y3 );
    double b  = 0.5 * (b1 + b2);

    double s, s2, t, t2;

    double am2bpc = a-2*b+c;
    /* this is how many valid s values we have */
    int num_valid_s = 0;

    if ( equals( am2bpc, 0, 1e-10 ) )
    {
     if ( equals( a-c, 0, 1e-10 ) )
     {
      /* Looks like the input is a line */
      /* You could set s=0.5 and solve for t if you wanted to */
      return 0;
     }
     s = a / (a-c);
     if ( in_range( s, 0, 1, 1e-10 ) )
      num_valid_s = 1;
    }
    else
    {
     double sqrtbsqmac = sqrt( b*b - a*c );
     s  = ((a-b) - sqrtbsqmac) / am2bpc;
     s2 = ((a-b) + sqrtbsqmac) / am2bpc;
     num_valid_s = 0;
     if ( in_range( s, 0, 1, 1e-10 ) )
     {
      num_valid_s++;
      if ( in_range( s2, 0, 1, 1e-10 ) )
       num_valid_s++;
     }
     else
     {
      if ( in_range( s2, 0, 1, 1e-10 ) )
      {
       num_valid_s++;
       s = s2;
      }
     }
    }

    if ( num_valid_s == 0 )
     return 0;

    t_valid = 0;
    if ( num_valid_s >= 1 )
    {
     double tdenom_x = (1-s)*(x0-x2) + s*(x1-x3);
     double tdenom_y = (1-s)*(y0-y2) + s*(y1-y3);
     t_valid = 1;
     if ( equals( tdenom_x, 0, 1e-10 ) && equals( tdenom_y, 0, 1e-10 ) )
     {
      t_valid = 0;
     }
     else
     {
      /* Choose the more robust denominator */
      if ( fabs( tdenom_x ) > fabs( tdenom_y ) )
      {
       t = ( (1-s)*(x0-x) + s*(x1-x) ) / ( tdenom_x );
      }
      else
      {
       t = ( (1-s)*(y0-y) + s*(y1-y) ) / ( tdenom_y );
      }
      if ( !in_range( t, 0, 1, 1e-10 ) )
       t_valid = 0;
     }
    }

    /* Same thing for s2 and t2 */
    t2_valid = 0;
    if ( num_valid_s == 2 )
    {
     double tdenom_x = (1-s2)*(x0-x2) + s2*(x1-x3);
     double tdenom_y = (1-s2)*(y0-y2) + s2*(y1-y3);
     t2_valid = 1;
     if ( equals( tdenom_x, 0, 1e-10 ) && equals( tdenom_y, 0, 1e-10 ) )
     {
      t2_valid = 0;
     }
     else
     {
      /* Choose the more robust denominator */
      if ( fabs( tdenom_x ) > fabs( tdenom_y ) )
      {
       t2 = ( (1-s2)*(x0-x) + s2*(x1-x) ) / ( tdenom_x );
      }
      else
      {
       t2 = ( (1-s2)*(y0-y) + s2*(y1-y) ) / ( tdenom_y );
      }
      if ( !in_range( t2, 0, 1, 1e-10 ) )
       t2_valid = 0;
     }
    }

    /* Final cleanup */
    if ( t2_valid && !t_valid )
    {
     s = s2;
     t = t2;
     t_valid = t2_valid;
     t2_valid = 0;
    }

    /* Output */
    if ( t_valid )
    {
     *sout = s;
     *tout = t;
    }

    if ( t2_valid )
    {
     *s2out = s2;
     *t2out = t2;
    }

    return t_valid + t2_valid;
}

void bilerp( double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double s, double t, double* x, double* y )
{
    *x = t*(s*x3+(1-s)*x2) + (1-t)*(s*x1+(1-s)*x0);
    *y = t*(s*y3+(1-s)*y2) + (1-t)*(s*y1+(1-s)*y0);
}

double randrange( double range_min, double range_max )
{
    double range_width = range_max - range_min;
    double rand01 = (rand() / (double)RAND_MAX);
    return (rand01 * range_width) + range_min;
}

/* Returns number of failed trials */
int fuzzTestInvBilerp( int num_trials )
{
    int num_failed = 0;

    double x0, y0, x1, y1, x2, y2, x3, y3, x, y, s, t, s2, t2, orig_s, orig_t;
    int num_st;
    int itrial;
    for ( itrial = 0; itrial < num_trials; itrial++ )
    {
     int failed = 0;
     /* Get random positions for the corners of the quad */
     x0 = randrange( -10, 10 );
     y0 = randrange( -10, 10 );
     x1 = randrange( -10, 10 );
     y1 = randrange( -10, 10 );
     x2 = randrange( -10, 10 );
     y2 = randrange( -10, 10 );
     x3 = randrange( -10, 10 );
     y3 = randrange( -10, 10 );
     /*x0 = 0, y0 = 0, x1 = 1, y1 = 0, x2 = 0, y2 = 1, x3 = 1, y3 = 1;*/
     /* Get random s and t */
     s = randrange( 0, 1 );
     t = randrange( 0, 1 );
     orig_s = s;
     orig_t = t;
     /* bilerp to get x and y */
     bilerp( x0, y0, x1, y1, x2, y2, x3, y3, s, t, &x, &y );
     /* invert */
     num_st = inverseBilerp( x0, y0, x1, y1, x2, y2, x3, y3, x, y, &s, &t, &s2, &t2 );
     if ( num_st == 0 )
     {
      failed = 1;
     }
     else if ( num_st == 1 )
     {
      if ( !(equals( orig_s, s, 1e-5 ) && equals( orig_t, t, 1e-5 )) )
       failed = 1;
     }
     else if ( num_st == 2 )
     {
      if ( !((equals( orig_s, s , 1e-5 ) && equals( orig_t, t , 1e-5 )) ||
          (equals( orig_s, s2, 1e-5 ) && equals( orig_t, t2, 1e-5 )) ) )
         failed = 1;
     }

     if ( failed )
     {
      num_failed++;
      printf("Failed trial %d\n", itrial);
     }
    }

    return num_failed;
}

int main( int argc, char** argv )
{
    int num_failed;
    srand( 0 );

    num_failed = fuzzTestInvBilerp( 100000000 );

    printf("%d of the tests failed\n", num_failed);
    getc(stdin);

    return 0;
}