Быстрая интерполяция данных сетки

Ответ 1

Конечно! Есть два варианта, которые делают разные вещи, но оба используют регулярную структуру исходных данных.

Первый 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). Однако имейте в виду, что если вы увеличиваете масштаб на 2x, например, вы будете масштабировать по всем осям.

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 "real-world", вам нужно будет преобразовать их в координаты "пикселя" на основе индекса.

Во всяком случае, допустим, мы хотели интерполировать значение в исходном массиве в позиции 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 исходного массива. Если вы хотите получить результаты с плавающей запятой, вам нужно будет отличить исходный массив как float:

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, вам нужно будет сбрасывать на float. ):

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])

Надеюсь, это немного прояснит ситуацию!

Ответ 2

Отличный ответ Джо. Основываясь на его предположении, я создал пакет regulargrid (https://pypi.python.org/pypi/regulargrid/, источник в https://github.com/JohannesBuchner/regulargrid)

Он обеспечивает поддержку n-мерных декартовых сеток (по мере необходимости здесь) с помощью очень быстрого scipy.ndimage.map_coordinates для произвольных координатных масштабов.