算法:如何计算双线性插值的INVERSE?映射到任意四边形的逆函数?

Too*_*eve 4 algorithm interpolation

更新:我下面的术语是错误的。我在“ Lerp2D”中描述的“正向”算法(需要逆运算)需要四个任意角。它沿着每个边都是线性的,但是所有4个都可以独立拉伸。它不是双线性的

我已经在标题中保留了双线性-如果您来这里是寻找“双线性的逆”,例如in x和中的独立拉伸y,请参见Spektre的答案

如果您需要更一般的情况(由任意四边形定义拉伸),请参见已接受的答案。

在此问题的评论中,另请参阅人们给出的链接。


原始问题:

双线性插值的计算很简单。但是我需要一种执行逆运算的算法。(算法对我来说可以是伪代码或任何广泛使用的计算机语言)

例如,这是双线性插值的Visual Basic实现。

' 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
Run Code Online (Sandbox Code Playgroud)

哪里

' 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
Run Code Online (Sandbox Code Playgroud)

' 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
Run Code Online (Sandbox Code Playgroud)

在1-D中,我确定了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
Run Code Online (Sandbox Code Playgroud)

现在我需要在二维中执行相同的操作:

' 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
Run Code Online (Sandbox Code Playgroud)

Spe*_*tre 5

定义双线性插值的逆。

您的代码对我来说不是很易读(不是VB编码器),可能有一些有关背后主要思想的附加文字会更好,但是我认为从您的代码中可以得到一些帮助。从我的角度来看,它看起来像这样:

双线性插值是2D图像/矩阵调整大小

  • 输入的是 分辨率为xs0,ys0和新分辨率为xs1,ys1的2D图像 /矩阵
  • 输出是 分辨率为xs1,ys1的2D图像 /矩阵

逆双线性插值从调整大小后的图像中获取原始2D图像/矩阵。仅当(xs0<=xs1) AND (ys0<=ys1)其他所需信息丢失时,才可以这样做。

逆双线性插值算法

  1. 在图像中找到原始栅格

    仅第一个处理线,并将具有相似斜率的点组合在一起(下图的光栅图上的绿线)。计算相邻线之间的交点并计算交点之间最常见的最小距离。

    使用相交距离的直方图表示应该有更多的候选对象,这些候选对象是原始栅格大小的倍数,因此请选择最小的一个。仅从这些乘法中进行选择,以避免出现压缩失真问题。这些是原始图像的网格上的点(半双线性插值)

  2. 内插栅格网格线

    根据项目符号#1的计算网格对点进行分组,并获取所有“蓝色”点的颜色。

  3. 获取栅格网格点

    在蓝色点(过程列)上应用项目符号#1#2,结果为原始图像。不要忘记再次计算网格大小,因为列可以具有与行不同的列。

    逆双线性插值

    • 图x轴为加工线/行轴
    • 图y轴是颜色/分量强度值

[Edit1]对此很好奇,所以我做了一些测试

这种方法可用于缩放zoom >= 2.0较小的(双)线性缩放图像,至少对于以下代码的当前状态没有精度提升(可以使用一些调整以使其更好)。

注意将反插值与您的插值匹配(如果不使用下面的方法),因为坐标映射+/- 1像素位置可能会有一些差异

如果您已计算出源分辨率,那么下面是C ++中的一些代码:

//---------------------------------------------------------------------------
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;
    }
//---------------------------------------------------------------------------
Run Code Online (Sandbox Code Playgroud)
  • lin_interpolate 是一维线性插值 src[ssz] -> dst[dsz]
  • lin_deinterpolate 是一维逆线性插值 src[ssz] -> dst[dsz]

现在要在2D模式下执行以下操作:

双线性插值:

首先通过插值1D重新缩放所有行,然后通过插值1D重新缩放列。重写这两个函数可以单步执行列而不是行,或者在列重新缩放之前和之后对图像进行转置。您还可以将每列复制到行缓冲区重新缩放,然后再复制回来,以便选择最喜欢的内容。

逆或逆双线性插值:

它与双线性插值相同,但顺序相反!因此,首先通过插值1D重新缩放列,然后通过插值1D重新缩放所有行。如果您不这样做,那么准确性可能会差一些。

对于我测试的10位整数灰度图像,逆双线性的准确度是纯逆双线性的3倍左右。

好的,这里是一些用于测试图像的真实gfx行:

真实数据

  • 红色箭头显示最大精度提升在哪里
  • 白色大点是原始图像中的点
  • 绿线/点是线性插值后的点
  • 以相同的坡度检测到水线的准确点 (difference<=treshold)
  • 蓝点是解插值后的输出点(接近可用的浅绿色线的交点)
  • 最好的情况是蓝点位于白点中间,但并非总是因为整数舍入误差

PS。提供的代码未优化

对于2D,不需要resize/alloc每个行/列的所有缓冲区,也init dmap[]可以对所有行一次和对所有列一次