Algorithm: how to calculate the INVERSION of bilinear interpolation?

Bilinear interpolation is trivial to calculate. But I need an algorithm that performs an INVERSE operation. (the algorithm will be useful to me in pseudocode or any widely used computer language)

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

' 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 defined the inverse Lerp function:

' 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
+4
source share
2 answers

, intepolated z.
, z 00, z 01, z 10, z 10 w 0 w 1 ,

z 0= z 00 + w 0 & times; (z 10 - z 00)
 z 1= z 01 + w 0 & times; (z 11 - z 01)

z = z 0 + w 1 & times; (z 1 - z 0)
   = z 00 + w 0 & times; (z 10 - z 00) + w 1 & times; (z 01 - z 00) + w 1 & times; w 0 & times; (z 11 - z 10 - z 01 + z 00)

,

x = x 00 + w 0 & times; (x 10 - x 00) + w 1 & times; (x 01 - x 00) + w 1 & times; w 0 & times; (x 11 - x 10 - x 01 + x 00)
 y = y 00 + w 0 & times; (y 10 - y 00) + w 1 & times; (y 01 - y 00) + w 1 & times; w 0 & times; (y 11 - y 10 - y 01 + y 00)

, w 0 w 1 x y.

(x w (w 0, w 1)) x) 2 + (y w (w 0, w 1) - y) 2

-.

:

, (x, y) (w 0, w 1), . , rev (fwd (w 0, w 1)), , (w 0, w 1), .

, , , . (x, y) .
, . , Delaunay .

, ,

1


x & times;

a i, b i, c i d i (i, 0 1) ,

w 0= a 0 + b 0 & times; x + c 0 & times; y + d 0 & times; x & times;
w 1= a 1 + b 1 & times; x + c 1 & times; y + d 1 & times; x & times;

x, y, w 0 w 1, w, , . , (, SVD), , .

, , , !

+1

.

( VB), , - , , . :

- 2D- /

  • 2D-/ xs0, ys0 xs1, ys1
  • 2D-/ xs1, ys1

2D-/ . , (xs0<=xs1) AND (ys0<=ys1) .

  • ( ). .

    , , , . , . ( )

  • # 1 "" .

  • bullet # 1 # 2 ( ). . , .

    Inverse bilinear interpolation

    • x /
    • y - /

[Edit1] ,

() zoom >= 2.0 , , ( , ).

, ( ), +/- 1

, , - ++:

//---------------------------------------------------------------------------
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 - 1D src[ssz] -> dst[dsz]
  • lin_deinterpolate - 1D src[ssz] -> dst[dsz]

, 2D, :

:

1D, 1D. , . , , .

:

, Bilinear Interpolate, !!! 1D, 1D. , .

10- , , 3 , .

OK - gfx :

real data

  • ,
  • -
  • / -
  • (difference<=treshold)
  • ( -).
  • - , - .

PS.,

2D resize/alloc / init dmap[]

+4

Source: https://habr.com/ru/post/1536712/


All Articles