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.
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.
ORIGINAL 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)
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
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 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
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
To simplify, let's begin by just considering a single intepolated value z.
Assume four values z00, z01, z10, z10, and two weights w0 and w1 applied to the first and second index, giving
z0 = z00 + w0 × (z10 - z00)
z1 = z01 + w0 × (z11 - z01)
and finally
z = z0 + w1 × (z1 - z0)
= z00 + w0 × (z10 - z00) + w1 × (z01 - z00) + w1 × w0 × (z11 - z10 - z01 + z00)
So, for your problem you will have to invert a two dimensional quadratic equation
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)
Unfortunately, there isn't a simple formula to recover w0 and w1 from x and y. You can, however, treat it as a non-linear least squares problem and minimise
(xw(w0,w1) - x)2 + (yw(w0,w1) - y)2
which you can do efficiently with the Levenberg–Marquardt algorithm.
Edit: Further Thoughts
It has occurred to me that you might be satisfied with an interpolation from (x, y) to (w0, w1) rather than the actual inverse. This will be less accurate in the sense that rev(fwd(w0, w1)) will likely be further from (w0, w1) than the actual inverse.
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
and compute coefficients ai, bi, ci and di (i equal to 0 or 1) such that
w0 = a0 + b0 × x + c0 × y + d0 × x × y
w1 = a1 + b1 × x + c1 × y + d1 × x × y
By substituting the relevant known values of x, y, w0 and w1 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!