Как выполнить билинейную интерполяцию в Python
Я хотел бы выполнить блинейную интерполяцию с использованием Python.
Пример точки GPS, для которой я хочу интерполировать высоту:
B = 54.4786674627
L = 17.0470721369
используя четыре смежные точки с известными координатами и значениями высоты:
n = [(54.5, 17.041667, 31.993), (54.5, 17.083333, 31.911), (54.458333, 17.041667, 31.945), (54.458333, 17.083333, 31.866)]
z01 z11
z
z00 z10
и вот моя примитивная попытка:
import math
z00 = n[0][2]
z01 = n[1][2]
z10 = n[2][2]
z11 = n[3][2]
c = 0.016667 #grid spacing
x0 = 56 #latitude of origin of grid
y0 = 13 #longitude of origin of grid
i = math.floor((L-y0)/c)
j = math.floor((B-x0)/c)
t = (B - x0)/c - j
z0 = (1-t)*z00 + t*z10
z1 = (1-t)*z01 + t*z11
s = (L-y0)/c - i
z = (1-s)*z0 + s*z1
где z0 и z1
z01 z0 z11
z
z00 z1 z10
Я получаю 31,964, но из другого программного обеспечения я получаю 31,961.
Мой сценарий правильный?
Можете ли вы предоставить другой подход?
9 ответов
Вот функция многократного использования, которую вы можете использовать. Включает в себя документы и проверки данных:
def bilinear_interpolation(x, y, points):
'''Interpolate (x,y) from values associated with four points.
The four points are a list of four triplets: (x, y, value).
The four points can be in any order. They should form a rectangle.
>>> bilinear_interpolation(12, 5.5,
... [(10, 4, 100),
... (20, 4, 200),
... (10, 6, 150),
... (20, 6, 300)])
165.0
'''
# See formula at: http://en.wikipedia.org/wiki/Bilinear_interpolation
points = sorted(points) # order points by x, then by y
(x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points
if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2:
raise ValueError('points do not form a rectangle')
if not x1 <= x <= x2 or not y1 <= y <= y2:
raise ValueError('(x, y) not within the rectangle')
return (q11 * (x2 - x) * (y2 - y) +
q21 * (x - x1) * (y2 - y) +
q12 * (x2 - x) * (y - y1) +
q22 * (x - x1) * (y - y1)
) / ((x2 - x1) * (y2 - y1) + 0.0)
Вы можете запустить тестовый код, добавив:
if __name__ == '__main__':
import doctest
doctest.testmod()
Запустив интерполяцию для вашего набора данных, вы получите:
>>> n = [(54.5, 17.041667, 31.993),
(54.5, 17.083333, 31.911),
(54.458333, 17.041667, 31.945),
(54.458333, 17.083333, 31.866),
]
>>> bilinear_interpolation(54.4786674627, 17.0470721369, n)
31.95798688313631
Не уверен, что это сильно поможет, но я получаю другое значение при выполнении линейной интерполяции с использованием scipy:
>>> import numpy as np
>>> from scipy.interpolate import griddata
>>> n = np.array([(54.5, 17.041667, 31.993),
(54.5, 17.083333, 31.911),
(54.458333, 17.041667, 31.945),
(54.458333, 17.083333, 31.866)])
>>> griddata(n[:,0:2], n[:,2], [(54.4786674627, 17.0470721369)], method='linear')
array([ 31.95817681])
Вдохновленный отсюда, я придумал следующий фрагмент. API оптимизирован для многократного использования одной и той же таблицы:
from bisect import bisect_left
class BilinearInterpolation(object):
""" Bilinear interpolation. """
def __init__(self, x_index, y_index, values):
self.x_index = x_index
self.y_index = y_index
self.values = values
def __call__(self, x, y):
# local lookups
x_index, y_index, values = self.x_index, self.y_index, self.values
i = bisect_left(x_index, x) - 1
j = bisect_left(y_index, y) - 1
x1, x2 = x_index[i:i + 2]
y1, y2 = y_index[j:j + 2]
z11, z12 = values[j][i:i + 2]
z21, z22 = values[j + 1][i:i + 2]
return (z11 * (x2 - x) * (y2 - y) +
z21 * (x - x1) * (y2 - y) +
z12 * (x2 - x) * (y - y1) +
z22 * (x - x1) * (y - y1)) / ((x2 - x1) * (y2 - y1))
Вы можете использовать это так:
table = BilinearInterpolation(
x_index=(54.458333, 54.5),
y_index=(17.041667, 17.083333),
values=((31.945, 31.866), (31.993, 31.911))
)
print(table(54.4786674627, 17.0470721369))
# 31.957986883136307
В этой версии нет проверки ошибок, и вы столкнетесь с проблемами, если попытаетесь использовать ее на границах индексов (или за их пределами). Полную версию кода, включая проверку ошибок и дополнительную экстраполяцию, смотрите здесь.
Простая реализация на основе этой формулы:
def bilinear_interpolation(x,y,x_,y_,val):
a = 1 /((x_[1] - x_[0]) * (y_[1] - y_[0]))
xx = np.array([[x_[1]-x],[x-x_[0]]],dtype='float32')
f = np.array(val).reshape(2,2)
yy = np.array([[y_[1]-y],[y-y_[0]]],dtype='float32')
b = np.matmul(f,yy)
return a * np.matmul(xx.T, b)
Ввод: Здесьx_ список [x0,x1] а также y_ список [y0,y1]
bilinear_interpolation(x=54.4786674627,
y=17.0470721369,
x_=[54.458333,54.5],
y_=[17.041667,17.083333],
val=[31.993,31.911,31.945,31.866])
Выход:
array([[31.95912739]])
Это то же решение, которое определено здесь, но применяется к некоторой функции и сравнивается с
interp2dдоступно в Scipy. Мы используем библиотеку numba, чтобы сделать функцию интерполяции даже быстрее, чем реализация Scipy.
import numpy as np
from scipy.interpolate import interp2d
import matplotlib.pyplot as plt
from numba import jit, prange
@jit(nopython=True, fastmath=True, nogil=True, cache=True, parallel=True)
def bilinear_interpolation(x_in, y_in, f_in, x_out, y_out):
f_out = np.zeros((y_out.size, x_out.size))
for i in prange(f_out.shape[1]):
idx = np.searchsorted(x_in, x_out[i])
x1 = x_in[idx-1]
x2 = x_in[idx]
x = x_out[i]
for j in prange(f_out.shape[0]):
idy = np.searchsorted(y_in, y_out[j])
y1 = y_in[idy-1]
y2 = y_in[idy]
y = y_out[j]
f11 = f_in[idy-1, idx-1]
f21 = f_in[idy-1, idx]
f12 = f_in[idy, idx-1]
f22 = f_in[idy, idx]
f_out[j, i] = ((f11 * (x2 - x) * (y2 - y) +
f21 * (x - x1) * (y2 - y) +
f12 * (x2 - x) * (y - y1) +
f22 * (x - x1) * (y - y1)) /
((x2 - x1) * (y2 - y1)))
return f_out
Мы делаем довольно большой массив интерполяции для оценки производительности каждого метода.
Пример функции:
x = np.linspace(0, 4, 13)
y = np.array([0, 2, 3, 3.5, 3.75, 3.875, 3.9375, 4])
X, Y = np.meshgrid(x, y)
Z = np.sin(np.pi*X/2) * np.exp(Y/2)
x2 = np.linspace(0, 4, 1000)
y2 = np.linspace(0, 4, 1000)
Z2 = bilinear_interpolation(x, y, Z, x2, y2)
fun = interp2d(x, y, Z, kind='linear')
Z3 = fun(x2, y2)
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(10, 6))
ax[0].pcolormesh(X, Y, Z, shading='auto')
ax[0].set_title("Original function")
X2, Y2 = np.meshgrid(x2, y2)
ax[1].pcolormesh(X2, Y2, Z2, shading='auto')
ax[1].set_title("bilinear interpolation")
ax[2].pcolormesh(X2, Y2, Z3, shading='auto')
ax[2].set_title("Scipy bilinear function")
plt.show()
Тест производительности
Python без библиотеки numba
bilinear_interpolation функция в данном случае такая же, как
numba версия за исключением того, что мы меняем
prange с нормальным питоном
range в цикле for и удалите декоратор функций
jit
%timeit bilinear_interpolation(x, y, Z, x2, y2)
Дает 7, 15 с ± 107 мс на цикл (среднее ± стандартное отклонение из 7 прогонов, по 1 циклу)
Python с numba numba
%timeit bilinear_interpolation(x, y, Z, x2, y2)
Дает 2,65 мс ± 70,5 мкс на цикл (среднее ± стандартное отклонение из 7 прогонов, по 100 циклов в каждом)
Scipy реализация
%%timeit
f = interp2d(x, y, Z, kind='linear')
Z2 = f(x2, y2)
Дает 6,63 мс ± 145 мкс на цикл (среднее ± стандартное отклонение из 7 прогонов, по 100 циклов в каждом)
Тесты производительности выполняются на процессоре Intel(R) Core(TM) i7-8700K @ 3,70 ГГц.
Я думаю, что смысл делать floor Функция заключается в том, что обычно вы хотите интерполировать значение, координата которого лежит между двумя дискретными координатами. Однако у вас, кажется, уже есть реальные значения координат ближайших точек, что упрощает математические вычисления.
z00 = n[0][2]
z01 = n[1][2]
z10 = n[2][2]
z11 = n[3][2]
# Let's assume L is your x-coordinate and B is the Y-coordinate
dx = n[2][0] - n[0][0] # The x-gap between your sample points
dy = n[1][1] - n[0][1] # The Y-gap between your sample points
dx1 = (L - n[0][0]) / dx # How close is your point to the left?
dx2 = 1 - dx1 # How close is your point to the right?
dy1 = (B - n[0][1]) / dy # How close is your point to the bottom?
dy2 = 1 - dy1 # How close is your point to the top?
left = (z00 * dy1) + (z01 * dy2) # First interpolate along the y-axis
right = (z10 * dy1) + (z11 * dy2)
z = (left * dx1) + (right * dx2) # Then along the x-axis
В переводе с вашего примера может быть немного ошибочной логики, но суть в том, что вы можете взвесить каждую точку на основе того, насколько она ближе к целевой точке интерполяции, чем ее другие соседи.
Предлагаю следующее решение:
def bilinear_interpolation(x, y, z01, z11, z00, z10):
def linear_interpolation(x, z0, z1):
return z0 * x + z1 * (1 - x)
return linear_interpolation(y, linear_interpolation(x, z01, z11),
linear_interpolation(x, z00, z10))
Решение здесь показывает билинейную интерполяцию, я применил его метод здесь . Но моя адаптация этого метода, просто вычисляющая близость ко всем четырем углам на основе евклидова расстояния и использующая это как простое средневзвешенное значение, работает намного лучше (моя адаптация находится по той же ссылке).
import numpy as np
def func(x, y):
s1 = 0.2; x1 = 36.5; y1 = 32.5
s2 = 0.4; x2 = 36.1; y2 = 32.8
g1 = np.exp( -4 *np.log(2) * ((x-x1)**2+(y-y1)**2) / s1**2)
g2 = np.exp( -2 *np.log(2) * ((x-x2)**2+(y-y2)**2) / s2**2)
return g1 + g2
D = 20
x = np.linspace(36,37,D)
y = np.linspace(32,33,D)
xx,yy = np.meshgrid(x,y)
zz = func(xx,yy)
def find_corners(xi,yi):
idx1 = np.searchsorted(x, xi, side="left")
idx2 = np.searchsorted(y, yi, side="left")
cs = [(idx2,idx1),(idx2-1,idx1),(idx2,idx1-1),(idx2-1,idx1-1)]
return cs
def cdist(p1,p2):
distances = np.linalg.norm(p1 - p2, axis=1)
return distances
def cell_interp(x, y, points):
a = np.array([x,y]).reshape(-1,2)
b = np.array(points)[:,:2]
ds = cdist(a,b)
ds = ds / np.sum(ds)
ds = 1. - ds
c = np.array(points)[:,2]
iz = np.sum(c * ds) / np.sum(ds)
return iz
def grid_interp(intx,inty):
cs = find_corners(intx,inty)
i,j = cs[0][0],cs[0][1]
i,j = cs[1][0],cs[1][1]
i,j = cs[2][0],cs[2][1]
i,j = cs[3][0],cs[3][1]
i0,j0 = cs[0][0],cs[0][1]
i1,j1 = cs[1][0],cs[1][1]
i2,j2 = cs[2][0],cs[2][1]
i3,j3 = cs[3][0],cs[3][1]
introw = [(xx[i0,j0],yy[i0,j0],zz[i0,j0]),
(xx[i1,j1],yy[i1,j1],zz[i1,j1]),
(xx[i2,j2],yy[i2,j2],zz[i2,j2]),
(xx[i3,j3],yy[i3,j3],zz[i3,j3])]
return cell_interp(intx,inty,introw)
x2 = np.linspace(36.0001,36.9999,D*2)
y2 = np.linspace(32.0001,32.9999,D*2)
xx2,yy2 = np.meshgrid(x2,y2)
zz2 = func(xx2,yy2)
grid_interp_vec = np.vectorize(grid_interp,otypes=[np.float])
zz2_grid = grid_interp_vec(xx2,yy2)
print (np.mean(np.square(zz2-zz2_grid)))