Algorithm: how calculate INVERSE of bilinear interpolation? INVERSE of mapping on to an arbitrary quadrilateral?

Question

UPDATE: My terminology below is wrong. The "forward" algorithm I describe in "Lerp2D" (which I need inverse-of) takes 4 arbitrary corners. It is linear along each edge, but all 4 edges can independently stretch; it is not bilinear.

I've left bilinear in the title - if you come here looking for "inverse of bilinear", e.g. independent stretching in x and y, see Spektre's answer.

If you need a more general case (stretching defined by an arbitrary quadrilateral), then see the accepted answer.

Also see links that people have given, in comments on this question.

---

ORIGINAL QUESTION:

Bilinear interpolation is trivial to compute. But I need an algorithm that does the INVERSE operation. (algorithm will be useful to me in pseudo-code, or any widely-used computer language)

For example, here is a Visual Basic implementation of bilinear interpolation.

' xyWgt ranges (0..1) in x and y. (0,0) will return X0Y0,
(0,1) will return X0Y1, etc.
' For example, if xyWgt is relative location within an image,
' and the XnYn values are GPS coords at the 4 corners of the image,
' The result is GPS coord corresponding to xyWgt.
' E.g. given (0.5, 0.5), the result will be the GPS coord at center of image.
Public Function Lerp2D(xyWgt As Point2D, X0Y0 As Point2D, X1Y0 As Point2D, X0Y1 As Point2D, X1Y1 As Point2D) As Point2D
    Dim xY0 As Point2D = Lerp(X0Y0, X1Y0, xyWgt.X)
    Dim xY1 As Point2D = Lerp(X0Y1, X1Y1, xyWgt.X)

    Dim xy As Point2D = Lerp(xY0, xY1, xyWgt.Y)
    Return xy
End Function

where

' Weighted Average of two points.
Public Function Lerp(ByVal a As Point2D, ByVal b As Point2D, ByVal wgtB As Double) As Point2D
    Return New Point2D(Lerp(a.X, b.X, wgtB), Lerp(a.Y, b.Y, wgtB))
End Function

and

' Weighted Average of two numbers.
' When wgtB==0, returns a, when wgtB==1, returns b.
' Implicitly, wgtA = 1 - wgtB. That is, the weights are normalized.
Public Function Lerp(ByVal a As Double, ByVal b As Double, ByVal wgtB As Double) As Double
    Return a + (wgtB * (b - a))
End Function

---

In 1-D, I have determined the inverse function of Lerp:

' Calculate wgtB that would return result, if did Lerp(a, b, wgtB).
' That is, where result is, w.r.t. a and b.
' < 0 is before a, > 1 is after b.
Public Function WgtFromResult(ByVal a As Double, ByVal b As Double, ByVal result As Double) As Double

    Dim denominator As Double = b - a

    If Math.Abs(denominator) < 0.00000001 Then
        ' Avoid divide-by-zero (a & b are nearly equal).
        If Math.Abs(result - a) < 0.00000001 Then
            ' Result is close to a (but also to b): Give simplest answer: average them.
            Return 0.5
        End If
        ' Cannot compute.
        Return Double.NaN
    End If

    ' result = a + (wgt * (b - a))   =>
    ' wgt * (b - a) = (result - a)   =>
    Dim wgt As Double = (result - a) / denominator

    'Dim verify As Double = Lerp(a, b, wgt)
    'If Not NearlyEqual(result, verify) Then
    '    Dim test = 0    ' test
    'End If

    Return wgt
End Function

Now I need to do the same in 2-D:

' Returns xyWgt, which if given to Lerp2D, would return this "xy".
' So if xy = X0Y0, will return (0, 0). if xy = X1Y0, will return (1, 0), etc.
' For example, if 4 corners are GPS coordinates in corners of an image,
' and pass in a GPS coordinate,
' returns relative location within the image.
Public Function InverseLerp2D(xy As Point2D, X0Y0 As Point2D, X1Y0 As Point2D, X0Y1 As Point2D, X1Y1 As Point2D) As Point2D
    ' TODO ????
End Function

Answer 1

To simplify, let's begin by just considering a single intepolated value z.
Assume four values z₀₀, z₀₁, z₁₀, z₁₀, and two weights w₀ and w₁ applied to the first and second index, giving

z₀ = z₀₀ + w₀ × (z₁₀ - z₀₀)
z₁ = z₀₁ + w₀ × (z₁₁ - z₀₁)

and finally

z = z₀ + w₁ × (z₁ - z₀)
   = z₀₀ + w₀ × (z₁₀ - z₀₀) + w₁ × (z₀₁ - z₀₀) + w₁ × w₀ × (z₁₁ - z₁₀ - z₀₁ + z₀₀)

So, for your problem you will have to invert a two dimensional quadratic equation

x = x₀₀ + w₀ × (x₁₀ - x₀₀) + w₁ × (x₀₁ - x₀₀) + w₁ × w₀ × (x₁₁ - x₁₀ - x₀₁ + x₀₀)
y = y₀₀ + w₀ × (y₁₀ - y₀₀) + w₁ × (y₀₁ - y₀₀) + w₁ × w₀ × (y₁₁ - y₁₀ - y₀₁ + y₀₀)

Unfortunately, there isn't a simple formula to recover w₀ and w₁ from x and y. You can, however, treat it as a non-linear least squares problem and minimise

(x_w(w₀,w₁) - x)² + (y_w(w₀,w₁) - y)²

which you can do efficiently with the Levenberg–Marquardt algorithm.

Edit: Further Thoughts

It has occurred to me that you might be satisfied with an interpolation from (x, y) to (w₀, w₁) rather than the actual inverse. This will be less accurate in the sense that rev(fwd(w₀, w₁)) will likely be further from (w₀, w₁) than the actual inverse.

The fact that you're interpolating over an irregular mesh rather than a regular grid is going to make this a trickier proposition. Ideally you should join up your (x, y) points with non-overlapping triangles and use barycentric coordinates to linearly interpolate.
For numerical stability you should avoid shallow, pointy triangles. Fortunately, the Delaunay triangulation satifies this requirement and isn't too difficult to construct in two dimensions.

If you would like your reverse interpolation to take a similar form to your forward interpolation you can use the basis functions

1
x
y
x × y

and compute coefficients a_i, b_i, c_i and d_i (i equal to 0 or 1) such that

w₀ = a₀ + b₀ × x + c₀ × y + d₀ × x × y
w₁ = a₁ + b₁ × x + c₁ × y + d₁ × x × y

By substituting the relevant known values of x, y, w₀ and w₁ you'll get four simultaneous linear equations for each w that you can solve to get its coefficients.
Ideally you should use a numerically stable matrix inversion algorithm that can cope with near singular matrices (e.g. SVD), but you may be able to get away with Gaussian elimination if you're careful.

Sorry I can't give you any simpler options, but I'm afraid that this really is a rather tricky problem!

Answer 2

Define inverse of bilinear interpolation.

Your code is not very readable for me (not a VB coder) may be some additional text about the main idea behind will be better but from your code you are returning some weight I assume. From mine point of view it looks like this:

Bilinear interpolation is 2D image/matrix resize

• input is 2D image/matrix of resolution xs0,ys0 and new resolution xs1,ys1
• output is 2D image/matrix of resolution xs1,ys1

Inverse Bilinear interpolation is getting the original 2D image/matrix from resized image. This can be done only if (xs0<=xs1) AND (ys0<=ys1) else needed information is lost.

Inverse Bilinear interpolation algorithm

1. find the original raster in image

   first process lines only and group points with similar slopes together (green lines at the diagram raster image below). Compute intersection between neighboring lines and compute the most common smallest distance between intersections.

   Use histogram of intersection distances for that there should be more candidates which are a multiply of the original raster size so choose the smallest one. Chose only from these multiplies to avoid compression distortions problems. These are the points on grid of the original image (half bilinear interpolation)

2. interpolate raster grid lines

   group points by computed grid from bullet #1 and obtain the color of all 'blue' points.

3. obtain raster grid points

   apply bullet #1 and #2 on blue points (process columns) The result is original image. Do not forget to compute the grid size again because columns can have different one then lines.

   

   • graph x axis is processed line/row axis
   • graph y axis is color/component intensity value

[Edit1] was curious about it so I did a little testing

This approach is usable for (bi)linear scaled images with zoom >= 2.0 for smaller zooms there is no accuracy boost at least for current state of the code below (it could use some tweaking to be better).

Be careful to match de-interpolate to your interpolation (if not use mine below) because there can be some differences in coordinates mapping +/- 1 pixel position

if you have the source resolution computed then here is some code in C++ for the rest:

//---------------------------------------------------------------------------
const int dmax=5;   // max difference of derivation (threshold)
//---------------------------------------------------------------------------
float divide(float a,float b)
    {
    if (fabs(b)<1e-6) return 0.0;
    return a/b;
    }
//---------------------------------------------------------------------------
// (x,y) = intersection of line(xa0,ya0,xa1,ya1) and line(xb0,yb0,xb1,yb1)
// return true if lines intersect
// ta , tb are parameters for intersection point inside line a,b
int line_intersection(float &x ,float &y ,
                      float xa0,float ya0,float xa1,float ya1,float &ta,
                      float xb0,float yb0,float xb1,float yb1,float &tb,float _zeroacc=1e-6)
    {
    float   dxa,dya,dxb,dyb,q;
    dxa=xa1-xa0; dya=ya1-ya0; ta=0;
    dxb=xb1-xb0; dyb=yb1-yb0; tb=0;
    q=(dxa*dyb)-(dxb*dya);
    if (fabs(q)<=_zeroacc) return 0;            // no intersection
    ta=divide(dxb*(ya0-yb0)+dyb*(xb0-xa0),q);
    tb=divide(dxa*(ya0-yb0)+dya*(xb0-xa0),q);
    x=xa0+(dxa*ta);
    y=ya0+(dya*ta);
    return 1;                                   // x,y = intersection ta,tb = parameters
    }
//---------------------------------------------------------------------------
void lin_interpolate(int *dst,int dsz,int *src,int ssz)
    {
    int x,x0,x1,c0,c1;
    int cnt,d0=ssz,d1=dsz;
    for (cnt=0,x0=0,x1=1,x=0;x<dsz;x++)
        {
        c0=src[x0];
        c1=src[x1];
        dst[x]=c0+(((c1-c0)*cnt)/d1);
        cnt+=d0; while (cnt>=d1) { cnt-=d1; x0++; x1++; if (x1>=ssz) x1=ssz-1; }
        }
    }
//---------------------------------------------------------------------------
void lin_deinterpolate(int *dst,int dsz,int *src,int ssz)
    {
    float px,py,ta,tb;
    int x,x0,x1,x2,x3,x4,x5;
    int   d0,d1,cnt;
    int  *der=new int[ssz]; // derivation by 'x' (address in array) ... slopes
    int *dmap=new int[ssz]; // corresponding x-positions in dst[]
    // init derivation table
    for (x0=0,x1=1;x1<ssz;x0++,x1++) der[x1]=src[x1]-src[x0]; der[0]=der[1];
    // init coordinate conversion table
    for (d0=dsz,d1=ssz,cnt=0,x0=0,x=0;x<ssz;x++) { dmap[x]=x0; cnt+=d0; while (cnt>=d1) { cnt-=d1; x0++; } }
    // init original lines
    int lins=0;
    int *lin0=new int[ssz<<1];
    int *lin1=new int[ssz<<1];
    for (x0=0,x1=0,x=0;x<ssz;x++)
        {
        d0=der[x0]-der[x]; if (d0<0) d0=-d0;
        if (d0<=dmax) x1=x;
        if ((d0>dmax)||(x+1>=ssz))
            {
            if (x0!=x1)
                {
                lin0[lins]=x0;
                lin1[lins]=x1;
                lins++;
                x0=x1; x=x1;
                }
            else{
                x0=x; x1=x;
                }
            }
        }

    // first naive interpolation
    lin_interpolate(dst,dsz,src,ssz);

    // update inaccurate points
    for (d0=0,d1=1;d1<lins;d0++,d1++)
        {
        x=lin0[d1]-lin1[d0];            // distance betwen lines
        if ((x<1)||(x>2)) continue;     // use only inacurate points
        // if lines found and intersect in the right place (x1<=px<=x2) copy result py to dst
        x0=lin0[d0];
        x1=lin1[d0];
        x2=lin0[d1];
        x3=lin1[d1];
        if (line_intersection(px,py,x0,src[x0],x1,src[x1],ta,x2,src[x2],x3,src[x3],tb))
         if ((px>=x1)&&(px<=x2))
            {
            dst[dmap[int(ceil(px))]]=int(py);
//          der[int(px)]=-300;
            }
        }
    delete[]  der;
    delete[] lin0;
    delete[] lin1;
    delete[] dmap;
    }
//---------------------------------------------------------------------------

• lin_interpolate is 1D linear interpolation src[ssz] -> dst[dsz]
• lin_deinterpolate is 1D inverse linear interpolation src[ssz] -> dst[dsz]

Now to do it in 2D just do this:

Bilinear Interpolate:

first rescale all rows by interpolate 1D and then rescale columns by interpolate 1D. Either rewrite both functions to step through column not row or transpose image before and after the column rescale. You can also copy each column to row buffer rescale and then copy back so choose what you like the most.

Inverse or Reverse Bilinear Interpolate:

It is the same as Bilinear Interpolate but in reverse order !!! so first rescale columns by interpolate 1D and then rescale all rows by interpolate 1D. If you do not then accuracy may be a bit worse.

For integer 10-bit gray-scale images I tested is the accuracy for inverse bilinear about 3 times better then naive reverse bilinear.

OK here is some real gfx line for tested image:

• red arrows shows where is max accuracy boost
• white big dots are points from original image
• green lines/dots are points after linear interpolation
• Aqua lines are detected accurate points with the same slope (difference<=treshold)
• blue dots are output points after de-interpolate (intersection of near usable aqua lines)
• in the best case is blue dot in the middle of white dot but not always because of integer rounding errors

PS. provided code is not optimized

for 2D you do not need to resize/alloc all buffers per each row/column also init dmap[] can be done once per all rows and once per all columns