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

4

更新: 我下面使用的术语是错误的。我在“Lerp2D”中描述的“forward”算法(我需要它的反向操作)需要4个任意角落。它沿着每条边线性,但所有4条边缘都可以独立拉伸; 它不是双线性的。

我将双线性保留在标题中 - 如果您在寻找“双线性的反向操作”,例如在xy中进行独立拉伸,请参见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

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

并且

' 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

在一维情况下,我已经确定了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

现在我需要在二维中做同样的事情:
' 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

嗯,我认为解决两个未知数的线性方程组是解决问题的方法。正在调查中。最近没有做过这样的事情,所以对细节有些生疏。我没有可调用的数学包,因此需要为这个特定的情况编写代码细节。 - ToolmakerSteve
1
请参考以下链接:这个链接,以及这个更好的重复问题 - Phrogz
@Phrogz:谢谢 - 这正是我需要的链接,我在搜索中不知何故找不到它们(可能当我开始研究这个问题时没有正确的术语;我发现“反双线性插值”确实可以得到这两个链接)。 - ToolmakerSteve
1
请查看程序图形大师Iñigo Quilez的文章。https://iquilezles.org/www/articles/ibilinear/ibilinear.htm - Darkgaze
2个回答

5

定义双线性插值的反函数。

对我来说,您的代码不太易读(我不是VB程序员),可能更好的做法是加入一些关于主要思想的附加文本,但从您的代码中,我认为您正在返回一些权重。在我的观点中,它看起来像这样:

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

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

双线性插值的反函数是从调整大小的图像中获取原始2D图像/矩阵。只有当(xs0<=xs1) AND (ys0<=ys1)时才能完成此操作,否则需要的信息将丢失。

双线性插值的反函数算法

  1. 找到图像中的原始光栅

    首先仅处理线,并将具有相似斜率的点分组在一起(下图光栅图像中的绿色线)。计算相邻线之间的交点并计算最常见的最小交点距离。

    使用交点距离的直方图,应该有更多候选项是原始光栅大小的倍数,因此选择最小值。仅从这些倍数中选择以避免压缩失真问题。这些是原始图像网格上的点(半双线性插值)

  2. 插值光栅网格线

    根据bullet#1计算出的网格将点分组,并获取所有“蓝色”点的颜色。

  3. 获取光栅网格点

    在蓝色点上应用bullet #1#2(处理列)的结果是原始图像。不要忘记重新计算网格大小,因为列可以具有与行不同的网格大小。

    Inverse Bilinear interpolation

    • 图形x轴是处理的行/行轴
    • 图形y轴是颜色/组件强度值

[编辑1]对此很好奇,所以我进行了一些测试

对于缩放比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;
    }
//---------------------------------------------------------------------------
  • lin_interpolate是一维线性插值src[ssz] -> dst[dsz]
  • lin_deinterpolate是一维反线性插值src[ssz] -> dst[dsz]

现在要在二维中进行操作,只需执行以下操作:

双线性插值:

首先通过对一维插值进行缩放来重新调整所有行,然后通过插值 1D 对列进行重新调整。要么重写两个函数以使其遍历列而不是行,要么在列重新调整之前和之后转换图像。您还可以将每列复制到行缓冲区进行缩放,然后再复制回来,因此选择您最喜欢的方式。

反向双线性插值:

与双线性插值相同,但顺序相反!!! 因此,首先通过插值 1D 缩放列,然后通过插值 1D 调整所有行。如果您不这样做,精度可能会略微降低。

对于整数10位灰度图像,我测试的逆双线性精度约为朴素逆双线性的3倍。

好的,这里是一些经过测试的实际图像:

real data

  • 红色箭头显示了精度最大的提升位置
  • 白色大点是原始图像中的点
  • 绿线/点是线性插值后的点
  • Aqua行是检测到的具有相同斜率的准确点(difference<=treshold)
  • 蓝点是反插值后的输出点(靠近可用Aqua行的交点)
  • 在最佳情况下,蓝点在白点中间,但不总是因为整数舍入误差

PS。 提供的代码未经优化

对于二维,您无需每个行/列都resize/alloc所有缓冲区,init dmap[]可以在所有行和所有列上完成一次。


感谢您花时间处理此事。不幸的是,这并不适用于我的情况。我心中想要的反函数不是线性的,而你的是(我认为)。很抱歉没有以一种语言中立的方式澄清数学问题。 - ToolmakerSteve
@ToolmakerSteve 双线性插值就是2倍于线性插值,它运行良好。你只需先转换所有行,然后再转换所有列,或按相反顺序进行转换(顺序会影响编码顺序的精度)。在2D调整大小的图像上已经成功测试过了...请查看答案中的 Now to do it in 2D just do this - Spektre
1
我已经更新了我的问题,解释了我的术语错误,并引导人们访问这个答案,如果他们需要双线性的反函数。 - ToolmakerSteve

2
为了简化问题,我们首先只考虑单个插值值z
假设有四个值z00z01z10z11和两个权重w0w1应用于第一个和第二个索引,得到以下结果: z0 = z00 + w0 × (z10 - z00)
z1 = z01 + w0 × (z11 - z01)
最终得到 z = z0 + w1 × (z1 - z0)
   = z00 + w0 × (z10 - z00) + w1 × (z01 - z00) + w1 × w0 × (z11 - z10 - z01 + z00)
因此,对于您的问题,您将不得不反转一个二维二次方程。

x = x00 + w0 × (x10 - x00) + w1 × (x01 - x00) + w1 × w0 × (x11 - x10 - x01 + x00)
y = y00 + w0 × (y10 - y00) + w1 × (y01 - y00) + w1 × w0 × (y11 - y10 - y01 + y00)

很遗憾,没有一个简单的公式可以从xy中恢复出w0w1。但是,您可以将其视为非线性最小二乘问题,并最小化

(xw(w0,w1) - x)2 + (yw(w0,w1) - y)2

您可以使用Levenberg-Marquardt算法高效地完成此操作。

编辑:进一步思考

我觉得你可能会满意从(xy)到(w0w1)的插值而不是实际的反向插值。这种方法精度较低,因为rev(fwd(w0w1))距离(w0w1)可能比实际的反向插值更远。

您正在对不规则网格进行插值而不是正常网格,这将使问题更加棘手。理想情况下,您应该用非重叠三角形连接您的(xy)点,并使用重心坐标进行线性插值。barycentric coordinates 为了保证数值稳定性,应避免使用浅、尖的三角形。幸运的是,Delaunay triangulation 满足这个要求,在两个维度中构建起来并不是太困难。

如果您希望您的反向插值采用与前向插值类似的形式,您可以使用basis functions

1
x
y
x × y

并计算系数aibicidi(其中i等于0或1),使得

w0 = a0 + b0 × x + c0 × y + d0 × x × y
w1 = a1 + b1 × x + c1 × y + d1 × x × y

通过替换已知的xyw0w1的值,您将得到每个w的四个同时的线性方程,可以解出其系数。
理想情况下,您应该使用能够处理几乎奇异矩阵(例如SVD)的数值稳定矩阵求逆算法,但如果您小心一点,也可能使用高斯消元法

对不起,我无法为您提供更简单的选项,但恐怕这确实是一个相当棘手的问题!


如何从2个方程组中得到6个未知数(z00,z01,z10,z11,w0,w1)?您还需要使用相邻点……顺便说一句,如果已知原始分辨率和偏移量,则可以计算出w0和w1,但仍然会剩下2个无法计算的未知数…… - Spektre
1
@Spektre:z项纯粹是为了推导双线性插值的一般公式。将它们替换为x,以获得x的公式,并用y替换以获得y的公式。 - thus spake a.k.
谢谢,这正是我需要知道的。首先,知道我没有忽略任何明显的东西,而且这确实比它所反映的插值要复杂得多。在我的情况下,畸变与线性相差不远:我没有进行最小化或构建三角网格,而是通过猜测/迭代技术解决了问题,这对我来说更加熟悉。感谢您的数学知识和相关主题的指针,以防将来需要更强大的方法。 - ToolmakerSteve
对于任何偶然发现这个问题和答案的读者,Phrogz在我的问题评论中提供的两个链接很好,以及谷歌搜索“反双线性插值”。(但在谷歌搜索中,您必须忽略大多数链接--那些仅描述“双线性插值”,而不讨论如何获得反向。) - ToolmakerSteve

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接