算法：如何计算双线性插值的INVERSE？映射到任意四边形的逆函数？ [英] Algorithm: how calculate INVERSE of bilinear interpolation? INVERSE of mapping on to an arbitrary quadrilateral?

问题描述

更新：我下面的术语是错误的。我在 Lerp2D中描述的正向算法（需要逆运算）需要四个任意角。它沿着每个边都是线性的，但是所有 4 个边可以独立伸展；

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.

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

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.

原始问题：

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

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)

例如，这是Visual Basic的bilinear实现

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

其中

' 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的反函数：

---

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

现在我需要在2-D中做同样的事情：

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

推荐答案

为简单起见，让我们开始吧只需考虑单个内插值 z 。

假定四个值 z ₀₀ ， z ₀₁ ， z ₁₀ ， z ₁₀ 和两个权重 w ₀ 和 w ₁ 应用于第一和第二个索引，给出

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 ₀₁ ）



和最后

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 ₀₀ ）

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 ₀₀ ）

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₀₀)

不幸的是，没有简单的公式可以恢复 w ₀ 和 w ₁ 来自 x 和 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 ） ²

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

您可以使用 Levenberg-Marquardt算法。

编辑：进一步的思考

在我看来，您可能对（ x y ）到（ w ₀ ， w ₁ ）而不是实际的倒数。从rev（fwd（ w ₀ ， w ₁ ））可能的意义上讲，这将不太准确。比实际的逆数更远离（ w ₀ ， w ₁ ）。

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.

事实上，您是在不规则网格而不是规则网格上进行插值，这将使这一建议变得更加棘手。理想情况下，您应该将（ x ， y ）点与不重叠的三角形连接起来，并使用重心坐标进行线性插值。

为了保持数值稳定性，应避免使用尖的三角形。幸运的是， Delaunay三角剖分满足了这一要求，并且也不是那么困难构造为二维。

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

1
x
y
x × y

并计算系数 a _i ， b _i ， c _i 和 d _i （ i 等于0或1），这样

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

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

通过代入 x ， y ， w ₀ 和 w ₁ ，对于每个，您将获得四个同时出现的 linear 方程> w 您可以求解以获得其系数。

理想情况下，您应该使用数值稳定的矩阵求逆算法，该算法可以处理近似奇异矩阵（例如 SVD ），但是您也许可以摆脱高斯消除，如果您注意的话。

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!

这篇关于算法：如何计算双线性插值的INVERSE？映射到任意四边形的逆函数？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！