Bilinear interpolation with non-aligned input points

Here's my own technique, along with code for deriving the resulting value. It requires three lerps of the output values (and three percentage calculations to determine the lerp percentages):

Note that this is not bilinear interpolation. It does not remap the quad of input points to the quad of output values, as some input points can result in output values outside the output quad.

Here I'm showing the non-aligned input values on a Cartesian plane (using the sample input values from the question above, multiplied by 10 for simplicity).

To calculate the 'north' point (top green dot), we calculate the percentage across the X axis as

  (inputX - northwestX) / (northeastX - northwestX)
= (-4.2 - -19) / (10 - -19)
= 0.51034

We use this percentage to calculate the intercept at the Y axis by lerping between the top Y values:

  (targetValue - startValue) * percent + startValue
= (northeastY  - northwestY) * percent + northwestY
= (-8 - -7) * 0.51034 + -7
= -7.51034

We do the same on the 'south' edge:

  (inputX - southwestX) / (southeastX - southwestX)
= (-4.2 - -11) / (9 - -11)
= 0.34

  (southeastY - southwestY) * percent + southwestY
= (7 - 4) * 0.34 + 4
= 5.02

Finally, we use these two values to calculate the final percentage between the north and south edges:

  (inputY - southY) / (northY - southY)
= (1 - 5.02) / (-7.51034 - 5.02)
= 0.3208

With these three percentages in hand we can calculate our final output values by lerping between the points:

nw = Vector(-150,-100)
ne = Vector( 150,-100)
sw = Vector(-150, 100)
se = Vector( 150, 100)

north  = lerp( nw, ne, 0.51034)       --> (  3.10, -100.00)
south  = lerp( sw, se, 0.34)          --> (-48.00,  100.00)
result = lerp( south, north, 0.3208)  --> (-31.61,   35.84)

Finally, here is some (Lua) code performing the above. It uses a mutable Vector object that supports the ability to copy values from another vector and lerp its values towards another vector.

-- Creates a bilinear interpolator
-- Corners should be an object with nw/ne/sw/se keys,
-- each of which holds a pair of mutable Vectors
-- { nw={inp=vector1, out=vector2}, … }
function tetragonalBilinearInterpolator(corners)
  local sides = {
    n={ pt=Vector(), pts={corners.nw, corners.ne} },
    s={ pt=Vector(), pts={corners.sw, corners.se} }
  }

  for _,side in pairs(sides) do
    side.minX = side.pts[1].inp.x
    side.diff = side.pts[2].inp.x - side.minX
  end

  -- Mutates the input vector to hold the result
  return function(inpVector)
    for _,side in pairs(sides) do
      local pctX = (inpVector.x - side.minX) / side.diff
      side.pt:copyFrom(side.pts[1].inp):lerp(side.pts[2].inp,pctX)
      side.inpY = side.pt.y
      side.pt:copyFrom(side.pts[1].out):lerp(side.pts[2].out,pctX)
    end
    local pctY = (inpVector.y-sides.s.inpY)/(sides.n.y-sides.s.inpY)
    return inpVector:copyFrom(sides.s.pt):lerp(sides.n.pt,pctY)
  end
end

local interp = tetragonalBilinearInterpolator{
  nw={ inp=Vector(-19,-7),  out=Vector(-150,-100) },
  ne={ inp=Vector( 10,-8),  out=Vector( 150,-100) },
  sw={ inp=Vector(-11, 4),  out=Vector(-150, 100) },
  se={ inp=Vector(  9, 7),  out=Vector( 150, 100) }
}
print(interp(Vector(-4.2, 1))) --> <-31.60 35.84>

One way is to convert your 4 input points to 4 new points that are in a rectangle.

You can use the min and max x and y values to create a rectangle that is outside the 4 original points. These are your new input points.

inputs:

import numpy as np
input = np.array([[-1.9,-0.7],[1.0,-0.8],[0.9,0.7],[-1.1,0.4]])
output = np.array([[-150,-100],[150,-100],[150,100],[-150,100]])

New points:

import matplotlib.pyplot as plt

x = input[:,0]
y = input[:,1]
plt.scatter(x, y,label="input")

x_sorted = np.sort(x)
y_sorted = np.sort(y)

x_grid = [x_sorted[0],x_sorted[0],x_sorted[-1],x_sorted[-1]]
y_grid = [y_sorted[0],y_sorted[-1],y_sorted[-1],y_sorted[0]]
plt.scatter(x_grid, y_grid,label="new_input",marker='+')

plt.legend(loc=(1.04, 0))
plt.show() 

result:

Use inverse bilinear interpolation/extrapolation to find predicted output values for your new input points.

The first step is to find the alpha (relative x) values and beta(relative y value) values for the original points:

alpha = (x-x_sorted[0])/(x_sorted[-1]-x_sorted[0])
beta = (y-y_sorted[0])/(y_sorted[-1]-y_sorted[0])

The next step is to use symbolic math or matrix inversion to solve the inverse bilinear interpolation/extrapolation.

matrix = np.zeros((4,4))

for i in range(4):
    matrix[i][0] =  alpha[i] * beta[i] - alpha[i] - beta[i] + 1
    matrix[i][1] = -alpha[i] * beta[i] + beta[i]
    matrix[i][2] =  alpha[i] * beta[i]
    matrix[i][3] = -alpha[i] * beta[i] + alpha[i]

output_grid = np.array([np.linalg.solve(matrix, output[:,0]), np.linalg.solve(matrix, output[:,1])])

When using bilinear interpolation using these new input points and output values you get the same output values on the original points.

from scipy.interpolate import RegularGridInterpolator

z1 = np.array([
    [output_grid[0][0], output_grid[0][3]],
    [output_grid[0][1], output_grid[0][2]]
]) 

z2 = np.array([
    [output_grid[1][0], output_grid[1][3]],
    [output_grid[1][1], output_grid[1][2]]
]) 

interp_func1 = RegularGridInterpolator(([y_sorted[0],y_sorted[-1]], [x_sorted[0],x_sorted[-1]]), z1)
interp_func2 = RegularGridInterpolator(([y_sorted[0],y_sorted[-1]], [x_sorted[0],x_sorted[-1]]), z2)

value = interp_func1([[input[0][1], input[0][0]]])[0]
print(value) # print x output of first original point
value = interp_func2([[input[0][1], input[0][0]]])[0]
print(value) # print y output of first original point