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

Question

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

我在标题中留下了双线性 - 如果你来这里寻找“双线性的逆”，例如x 和 y 中的独立拉伸(stretch)，参见 Spektre's answer 。

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

另请参阅人们在此问题的评论中提供的链接。

原始问题:

双线性插值计算起来很简单。
但我需要一个执行 INVERSE 操作的算法。
(算法将在伪代码或任何广泛使用的计算机语言中对我有用)

例如，这里是双线性插值的 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

在哪里

' 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

在 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

现在我需要在二维中做同样的事情:

' 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

为简化起见，让我们首先考虑单个内插值 z。
假设四个值 z00、z01、z10、z10，以及应用于第一个和第二个索引的两个权重 w0 和 w1，给出

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)

不幸的是，没有一个简单的公式可以从 x 和 y 中恢复 w0 和 w1。但是，您可以将其视为非线性最小二乘问题并最小化

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

您可以使用 Levenberg–Marquardt algorithm 有效地做到这一点。

编辑 : 进一步的想法

我突然想到，您可能会对从 (x, y) 到 (w0, w1) 的插值感到满意，而不是实际的逆。这将不太准确，因为 rev(fwd(w0, w1)) 可能比实际的倒数更远离 (w0, w1)。

您在不规则网格而不是规则网格上进行插值的事实将使这成为一个更棘手的命题。理想情况下，您应该将 (x, y) 点与不重叠的三角形连接起来，并使用 barycentric coordinates 进行线性插值。
为了数值稳定性，您应该避免使用浅而尖的三角形。幸运的是，Delaunay triangulation 满足了这个要求，并且在二维中构建起来并不太难。

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

1
X
是
x × y

并计算系数 ai、bi、ci 和 di(i 等于 0 或 1)，使得

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

通过替换 x、y、w0 和 w1 的相关已知值，您将获得每个 w 的四个联立线性方程，您可以求解这些方程以获得其系数。
理想情况下，您应该使用可以处理近似奇异矩阵的数值稳定的矩阵求逆算法(例如 SVD )，但是如果您小心的话，您可能能够逃脱 Gaussian elimination 。

对不起，我不能给你任何更简单的选择，但恐怕这真的是一个相当棘手的问题!