Быстрая интерполяция данных сетки
У меня есть большой 3d np.ndarray данных, который представляет физическую переменную, выбранную по объему в регулярной сетке (как в значении в массиве [0,0,0] представляет значение в физических координатах (0,0,0)).
Я хотел бы перейти к более мелкому интервалу сетки путем интерполяции данных в грубой сетке. В настоящее время я использую линейную интерполяцию scipy griddata, но она довольно медленная (~90 сек для массива 20x20x20). Это немного переоценено для моих целей, позволяя случайную выборку объемных данных. Есть ли что-нибудь, что может использовать мои регулярно расположенные данные и тот факт, что есть только ограниченный набор конкретных точек, к которым я хочу интерполировать?
2 ответа
Конечно! Есть два варианта, которые делают разные вещи, но оба используют регулярный характер исходных данных.
Первый scipy.ndimage.zoom, Если вы просто хотите создать более плотную регулярную сетку, основанную на интерполяции исходных данных, это путь.
Второй scipy.ndimage.map_coordinates, Если вы хотите интерполировать несколько (или много) произвольных точек в ваших данных, но при этом использовать регулярную сетку исходных данных (например, не требуется квадродерево), это путь.
"Масштабирование" массива (scipy.ndimage.zoom)
В качестве быстрого примера (Это будет использовать кубическую интерполяцию. Использование order=1 для билинейного, order=0 для ближайших и т. д.):
import numpy as np
import scipy.ndimage as ndimage
data = np.arange(9).reshape(3,3)
print 'Original:\n', data
print 'Zoomed by 2x:\n', ndimage.zoom(data, 2)
Это дает:
Original:
[[0 1 2]
[3 4 5]
[6 7 8]]
Zoomed by 2x:
[[0 0 1 1 2 2]
[1 1 1 2 2 3]
[2 2 3 3 4 4]
[4 4 5 5 6 6]
[5 6 6 7 7 7]
[6 6 7 7 8 8]]
Это также работает для 3D (и nD) массивов. Однако имейте в виду, что если вы увеличите, например, в 2 раза, вы увеличите масштаб по всем осям.
data = np.arange(27).reshape(3,3,3)
print 'Original:\n', data
print 'Zoomed by 2x gives an array of shape:', ndimage.zoom(data, 2).shape
Это дает:
Original:
[[[ 0 1 2]
[ 3 4 5]
[ 6 7 8]]
[[ 9 10 11]
[12 13 14]
[15 16 17]]
[[18 19 20]
[21 22 23]
[24 25 26]]]
Zoomed by 2x gives an array of shape: (6, 6, 6)
Если у вас есть что-то вроде 3-полосного RGB-изображения, которое вы хотите увеличить, вы можете сделать это, указав последовательность кортежей в качестве коэффициента увеличения:
print 'Zoomed by 2x along the last two axes:'
print ndimage.zoom(data, (1, 2, 2))
Это дает:
Zoomed by 2x along the last two axes:
[[[ 0 0 1 1 2 2]
[ 1 1 1 2 2 3]
[ 2 2 3 3 4 4]
[ 4 4 5 5 6 6]
[ 5 6 6 7 7 7]
[ 6 6 7 7 8 8]]
[[ 9 9 10 10 11 11]
[10 10 10 11 11 12]
[11 11 12 12 13 13]
[13 13 14 14 15 15]
[14 15 15 16 16 16]
[15 15 16 16 17 17]]
[[18 18 19 19 20 20]
[19 19 19 20 20 21]
[20 20 21 21 22 22]
[22 22 23 23 24 24]
[23 24 24 25 25 25]
[24 24 25 25 26 26]]]
Произвольная интерполяция данных с регулярной сеткой с использованием map_coordinates
Первое, что нужно понять о map_coordinates в том, что он работает в пиксельных координатах (например, так же, как вы бы индексировали массив, но значения могут быть плавающими). Из вашего описания, это именно то, что вы хотите, но если часто смущает людей. Например, если у вас есть x, y, z "реальные" координаты, вам нужно преобразовать их в основанные на индексе "пиксельные" координаты.
Во всяком случае, скажем, мы хотели интерполировать значение в исходном массиве в позиции 1.2, 0.3, 1.4.
Если вы думаете об этом с точки зрения более раннего случая изображения RGB, первая координата соответствует "полосе", вторая - "строке", а последняя - "столбцу". Какой порядок соответствует тому, что полностью зависит от того, как вы решите структурировать свои данные, но я собираюсь использовать их как координаты "z, y, x", поскольку это облегчает визуализацию сравнения с печатным массивом.
import numpy as np
import scipy.ndimage as ndimage
data = np.arange(27).reshape(3,3,3)
print 'Original:\n', data
print 'Sampled at 1.2, 0.3, 1.4:'
print ndimage.map_coordinates(data, [[1.2], [0.3], [1.4]])
Это дает:
Original:
[[[ 0 1 2]
[ 3 4 5]
[ 6 7 8]]
[[ 9 10 11]
[12 13 14]
[15 16 17]]
[[18 19 20]
[21 22 23]
[24 25 26]]]
Sampled at 1.2, 0.3, 1.4:
[14]
Еще раз, это кубическая интерполяция по умолчанию. Использовать order kwarg для управления типом интерполяции.
Здесь стоит отметить, что все scipy.ndimage Операции сохраняют dtype исходного массива. Если вы хотите получить результаты с плавающей запятой, вам нужно привести исходный массив в число с плавающей точкой:
In [74]: ndimage.map_coordinates(data.astype(float), [[1.2], [0.3], [1.4]])
Out[74]: array([ 13.5965])
Еще одна вещь, которую вы можете заметить, это то, что формат интерполированных координат довольно громоздок для одной точки (например, он ожидает массив 3xN вместо массива Nx3). Тем не менее, это возможно лучше, когда у вас есть последовательности координат. Например, рассмотрим случай выборки по линии, проходящей через "куб" данных:
xi = np.linspace(0, 2, 10)
yi = 0.8 * xi
zi = 1.2 * xi
print ndimage.map_coordinates(data, [zi, yi, xi])
Это дает:
[ 0 1 4 8 12 17 21 24 0 0]
Это также хорошее место, чтобы упомянуть, как обрабатываются граничные условия. По умолчанию все, что находится за пределами массива, установлено в 0. Таким образом, последние два значения в последовательности 0, (т.е. zi >> 2 для последних двух элементов).
Если бы мы хотели, чтобы точки вне массива были, скажем, -999 (Мы не можем использовать nan так как это целочисленный массив. Если ты хочешь nan Вам нужно будет бросить в поплавки.)
In [75]: ndimage.map_coordinates(data, [zi, yi, xi], cval=-999)
Out[75]: array([ 0, 1, 4, 8, 12, 17, 21, 24, -999, -999])
Если бы мы хотели, чтобы он возвращал ближайшее значение для точек вне массива, мы бы сделали:
In [76]: ndimage.map_coordinates(data, [zi, yi, xi], mode='nearest')
Out[76]: array([ 0, 1, 4, 8, 12, 17, 21, 24, 25, 25])
Вы также можете использовать "reflect" а также "wrap" в качестве граничных мод, в дополнение к "nearest" и по умолчанию "constant", Это довольно очевидно, но попробуйте немного поэкспериментировать, если вы запутались.
Например, давайте интерполируем линию вдоль первой строки первой полосы в массиве, которая простирается на двойное расстояние массива:
xi = np.linspace(0, 5, 10)
yi, zi = np.zeros_like(xi), np.zeros_like(xi)
По умолчанию дают:
In [77]: ndimage.map_coordinates(data, [zi, yi, xi])
Out[77]: array([0, 0, 1, 2, 0, 0, 0, 0, 0, 0])
Сравните это с:
In [78]: ndimage.map_coordinates(data, [zi, yi, xi], mode='reflect')
Out[78]: array([0, 0, 1, 2, 2, 1, 2, 1, 0, 0])
In [78]: ndimage.map_coordinates(data, [zi, yi, xi], mode='wrap')
Out[78]: array([0, 0, 1, 2, 0, 1, 1, 2, 0, 1])
Надеюсь, это немного прояснит ситуацию!
Отличный ответ от Джо. По его предложению я создал пакет регулярной сетки ( https://pypi.python.org/pypi/regulargrid/, источник по адресу https://github.com/JohannesBuchner/regulargrid).
Он обеспечивает поддержку n-мерных декартовых сеток (при необходимости здесь) с помощью очень быстрого scipy.ndimage.map_coordinates для произвольных координатных масштабов.