Многомерная сплайн-интерполяция в python / scipy? - PullRequest
27 голосов
/ 04 июня 2011

Существует ли библиотечный модуль или другой простой способ реализации многомерной сплайн-интерполяции в python?

В частности, у меня есть набор скалярных данных на равномерно распределенной трехмерной сетке, которые мне нужно интерполировать в небольшом количестве точек, разбросанных по всей области. Для двух измерений я использовал scipy.interpolate.RectBivariateSpline , и я, по сути, ищу расширение этого для трехмерных данных.

N-мерные процедуры интерполяции, которые я нашел, не достаточно хороши: я бы предпочел сплайны более LinearNDInterpolator для гладкости, и у меня слишком много точек данных (часто более миллиона) для, например , радиальная базисная функция для работы.

Если кто-нибудь знает библиотеку python, которая может это сделать, или, возможно, одну на другом языке, который я мог бы назвать или портировать, я был бы очень признателен.

Ответы [ 2 ]

43 голосов
/ 04 июня 2011

Если я правильно понимаю ваш вопрос, ваши входные данные "наблюдения" регулярно отображаются в сетке?

Если это так, scipy.ndimage.map_coordinates делает именно то, что вы хотите.

Это немного сложно понять при первом проходе, но, по сути, вы просто передаете ему последовательность координат, в которую вы хотите интерполировать значения сетки в координатах пикселя / вокселя / n-мерного индекса.

В качестве примера 2D:

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt

# Note that the output interpolated coords will be the same dtype as your input
# data.  If we have an array of ints, and we want floating point precision in
# the output interpolated points, we need to cast the array as floats
data = np.arange(40).reshape((8,5)).astype(np.float)

# I'm writing these as row, column pairs for clarity...
coords = np.array([[1.2, 3.5], [6.7, 2.5], [7.9, 3.5], [3.5, 3.5]])
# However, map_coordinates expects the transpose of this
coords = coords.T

# The "mode" kwarg here just controls how the boundaries are treated
# mode='nearest' is _not_ nearest neighbor interpolation, it just uses the
# value of the nearest cell if the point lies outside the grid.  The default is
# to treat the values outside the grid as zero, which can cause some edge
# effects if you're interpolating points near the edge
# The "order" kwarg controls the order of the splines used. The default is 
# cubic splines, order=3
zi = ndimage.map_coordinates(data, coords, order=3, mode='nearest')

row, column = coords
nrows, ncols = data.shape
im = plt.imshow(data, interpolation='nearest', extent=[0, ncols, nrows, 0])
plt.colorbar(im)
plt.scatter(column, row, c=zi, vmin=data.min(), vmax=data.max())
for r, c, z in zip(row, column, zi):
    plt.annotate('%0.3f' % z, (c,r), xytext=(-10,10), textcoords='offset points',
            arrowprops=dict(arrowstyle='->'), ha='right')
plt.show()

enter image description here

Чтобы сделать это в n-измерениях, нам просто нужно передать массивы соответствующего размера:

import numpy as np
from scipy import ndimage

data = np.arange(3*5*9).reshape((3,5,9)).astype(np.float)
coords = np.array([[1.2, 3.5, 7.8], [0.5, 0.5, 6.8]])
zi = ndimage.map_coordinates(data, coords.T)

Что касается масштабирования и использования памяти, map_coordinates создаст отфильтрованную копию массива, если вы используете порядок> 1 (т. Е. Не линейную интерполяцию).Если вы просто хотите интерполировать в очень небольшом количестве точек, это довольно большие издержки.Однако оно не увеличивается с количеством точек, в которые вы хотите интерполировать.Если у вас достаточно ОЗУ для одной временной копии массива входных данных, все будет в порядке.

Если вы не можете сохранить копию своих данных в памяти, вы можете либо: а) указать prefilter=False и order=1 и использовать линейную интерполяцию, либо б) заменить исходные данные отфильтрованной версией, используяndimage.spline_filter, а затем вызовите map_coordinates с prefilter=False.

Даже если у вас достаточно оперативной памяти, сохранение отфильтрованного набора данных может быть большим ускорением, если вам нужно несколько раз вызывать map_coordinates (например, интерактивное использование и т. Д.).

2 голосов
/ 04 июня 2011

Гладкую сплайн-интерполяцию в dim> 2 трудно реализовать, и поэтому не так много свободно доступных библиотек, способных сделать это (на самом деле, я не знаю ни одной).

Вы можете попробовать интерполяцию с обратным взвешенным расстоянием, см .: Интерполяция с обратным взвешенным расстоянием (IDW) с Python . Это должно привести к достаточно плавным результатам и лучше масштабироваться, чем RBF, к более крупным наборам данных.

...