Блог

Скажем, меня интересует

 y = f(n, m, o)
  

Я создаю дискретизированные сетки для n , m , o :

 import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')
  

Для простоты давайте предположим f(n,m,o) = n m p :

 Y = N M O
  

Итак, теперь у меня есть аппроксимация по сетке f внутри Y . Из моего взгляда на документацию,

• Мне нужно взглянуть на многомерную интерполяцию
• scipy.interpolate.Griddata выглядит великолепно, но он предназначен для неструктурированных сеток. Это означает, что я должен перевести все мои результаты в длинные массивы — это не использует преимущества того, что у меня действительно есть структурированные сетки
• Другие методы, перечисленные в разделе «структурированный«, не поддерживают кубические и многомерные (более 2 измерений) входные данные одновременно.

Учитывая переменные, которые я уже создал, каков хороший способ интерполяции по этому вопросу?

Обновить

Это реализация первого ответа:

 import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')
Y = N M O

n,m,o = 5, 103, 1007

m_i = np.argmin(np.abs(np.floor(m)-mGrid))
m_f = m - m_i

n_i = np.argmin(np.abs(np.floor(n)-nGrid))
n_f = n - n_i

o_i = np.argmin(np.abs(np.floor(o)-oGrid))
o_f = o - o_i

A = Y[m_i-1:m_i 3, n_i-1:n_i 3, o_i-1:o_i 3]

# cubic smoothing polynome in 1D
# Catmull-Rom style
def v(p0, p1, p2, p3, f):
    return 0.5 * (p0 * (-f**3   2*f**2 - f)  
                  p1 * (3*f**3 - 5*f**2   2)  
                  p2 * (-3*f**3   4*f**2   f)  
                  p3 * (f**3 - f**2))



B = v(A[0], A[1], A[2], A[3], m_f)
C = v(B[0], B[1], B[2], B[3], n_f)
D = v(C[0], C[1], C[2], C[3], o_f)

# D is the final interpolated array

print m n o, D
  

К сожалению, f(n,m,o) = 1115 , пока D = 2215 . Поскольку нет ссылки на методологию, мне трудно понять, что именно происходит и почему приближение так далеко.

Если вы не нашли подходящую, создайте свою собственную!

Но имейте в виду, что трехкубическая интерполяция немного сложна. Ниже описывается простая трехкубическая интерполяция. В зависимости от ваших потребностей доступны гораздо более сложные и быстрые методы (на основе матрицы). Следует также знать о том факте, что существует несколько методов кубической экстраполяции, я просто использую тот, который, вероятно, наиболее распространен (Catmull-Rom).

При поиске точки (m, n, o) в 3D-массиве M ее окружение для трехкубической интерполяции состоит из сетки размером 4x4x4 точки (или 3x3x3 ячейки). Точки определяются на каждой оси с помощью: floor(a)-1 , floor(a) , floor(a) 1 , floor(a) 2 и т.д. для каждой оси. (Подумайте о кубике Рубика, где точка (m, n, o) находится в невидимом центральном кубе.)

Результатом трехкубической интерполяции является средневзвешенное значение этих 64 точек. Веса зависят от расстояния точки от точек сетки.

Давайте определим:

 m_i = np.floor(m).astype('int')    # integer part
m_f = m - m_i                      # fraction
  

И аналогичные определения для n и o. Теперь нам нужно поработать с массивом:

 A = M[m_i-1:m_i 3, n_i-1:n_i 3, o_i-1:o_i 3]
  

Теперь наша точка находится в (1 m_f, 1 n_f, 1 o_f). (Обратите внимание, что это предполагает, что матрица бесконечна. В противном случае необходимо учитывать краевые эффекты.)

64 коэффициента массива могут быть вычислены с помощью некоторой умной матричный математики, но здесь достаточно знать, что интерполяция является ассоциативной, то есть мы можем делать это по одной оси за раз:

 # cubic smoothing polynome in 1D
# Catmull-Rom style
def v(p0, p1, p2, p3, f):
    return 0.5 * (p0 * (-f**3   2*f**2 - f)  
                  p1 * (3*f**3 - 5*f**2   2)  
                  p2 * (-3*f**3   4*f**2   f)  
                  p3 * (f**3 - f**2))


# interpolate axis-by-axis
B = v(A[0], A[1], A[2], A[3], m_f)
C = v(B[0], B[1], B[2], B[3], n_f)
D = v(C[0], C[1], C[2], C[3], o_f)

# D is the final interpolated value
  

Таким образом, интерполяция вычисляется по одной оси за раз.

Давайте превратим это во что-то более практичное. Во-первых, нам нужно выполнить некоторое вычисление адреса из нашего собственного пространства координат в массив координат. Если определена ось m:

 m0: position of the first grid plane on the axis
m1: position of the last grid plane on the axis
n_m: number of samples
  

мы можем вычислить эквивалентную позицию массива из позиции «реального мира» m :

 m_arr = (n_m - 1) * (m -  m0) / (m1 - m0)
  

Значение -1 происходит из-за явления fencepost (если у нас есть 10 выборок, расстояние между первой и последней равно 9).

Теперь мы можем объединить все в функцию:

 # cubic smoothing polynome in 1D
# Catmull-Rom style
def interp_catmull_rom(p0, p1, p2, p3, f):
    return 0.5 * (p0 * (-f**3   2*f**2 - f)  
                  p1 * (3*f**3 - 5*f**2   2)  
                  p2 * (-3*f**3   4*f**2   f)  
                  p3 * (f**3 - f**2))

# linear interpolation
# only needs p1 and p2, the rest are there for compatibility
def interp_linear(p0, p1, p2, p3, f):
    return (1-f)*p1   f*p2


# don't interpolate, use the nearest point
# only needs p1 and p2, the rest are there for compatibility
def interp_nearest(p0, p1, p2, p3, f):
    if f > .5:
        return p2
    else:
        return p1

# 3D interpolation done axis-by-axis
def tri_interp(M, f, m, n, o, m0, m1, n0, n1, o0, o1):
    # M: 3D array
    # f: interpolation function to use
    # m,n,o: coordinates where to interpolate
    # m0: real world minimum for m
    # m1: real world maximum for m
    # n0,n1,o0,o0: as with m0 and m1

    # calculate the array coordinates
    m_arr = (M.shape[0] - 1) * (m - m0) / (m1 - m0)
    n_arr = (M.shape[1] - 1) * (n - n0) / (n1 - n0)
    o_arr = (M.shape[2] - 1) * (o - o0) / (o1 - o0)

    # if we are out of our data box, return a nan
    if m_arr < 0 or m_arr > M.shape[0] or 
       n_arr < 0 or n_arr > M.shape[1] or 
       o_arr < 0 or o_arr > M.shape[2]:
        return np.nan

    # calculate the integer parts and fractions
    m_i = np.floor(m_arr).astype('int')
    n_i = np.floor(n_arr).astype('int')
    o_i = np.floor(o_arr).astype('int')
    m_f = m_arr - m_i
    n_f = n_arr - n_i
    o_f = o_arr - o_i

    # edge effects may be nasty, we may need elements outside of the array
    # there may be more efficient ways to avoid it, but we'll create lists of
    # coordinates:

    n_coords, m_coords, o_coords = np.mgrid[m_i-1:m_i 3, n_i-1:n_i 3, o_i-1:o_i 3]

    # these coordinate arrays are clipped so that we are in the available data
    # for example, point (-1,3,7) will use the point (0,3,7) instead
    m_coords = m_coords.clip(0, M.shape[0]-1)
    n_coords = n_coords.clip(0, M.shape[1]-1)
    o_coords = o_coords.clip(0, M.shape[2]-1)

    # get the 4x4x4 cube:
    A = M[m_coords, n_coords, o_coords]

    # interpolate along the first axis (3D to 2D)
    B = f(A[0], A[1], A[2], A[3], m_f)

    # interpolate along the second axis (2D to 1D)
    C = f(B[0], B[1], B[2], B[3], n_f)

    # interpolate along the third axis (1D to scalar)
    D = f(C[0], C[1], C[2], C[3], o_f)

    return D
  

Теперь эта функция позволяет использовать любой метод интерполяции, который может выполняться по оси для четырех соседних точек. Теперь отображаются кубические, линейные и «ближайшие».

Но работает ли это? Давайте протестируем это, сгенерировав некоторые случайные данные и плоскость, пересекающую данные. (Последнее на удивление сложно.)

 # random data
M = np.random.random((10, 12, 16))
m0,m1 = -10.,10.
n0,n1 = -7.,5.
o0,o1 = -4.,5.

# create a grid (grid from -15..15,-15..15 in n-m plane)
gr = mgrid[-15:15.01:.1, -15:15.01:.1]

# create two perpendicular unit vectors (forming a plane)
# vn: normal vector of the plane
# vp: some vector, whose projection determines one direction
# v0: unit vector on the plane (perpendicular to vn and vp)
# v1: unit vector on the plane (perpendicular to vn and v0)
vn = np.array([-.2, .3, 1])
vp = np.array([0, -1, 0])
v1 = np.cross(vn, vp)
v2 = np.cross(vn, v1)
v1 /= numpy.linalg.norm(v1)
v2 /= numpy.linalg.norm(v2)

# the grid and the unit vectors define the 3d points on the plane
gr3d = gr[0][:,:,None] * v1   gr[1][:,:,None] * v2

# now we can fetch the points at grid points, just must flatten and reshape back
res = [ tri_interp(M, interp_catmull_rom, p[0], p[1], p[2], m0,m1,n0,n1,o0,o1) for p in gr3d.reshape(-1,3) ]
res = np.array(res).reshape((gr3d.shape[0], gr3d.shape[1]))
  

Итак, у нас есть плоскость, которая проходит через наш блок данных. Теперь мы действительно можем увидеть разницу между тремя различными методами интерполяции, выбрав interp_linear , interp_nearest или interp_catmull_rom в качестве функции одномерной интерполяции. Все приведенные ниже изображения выполнены в одинаковом масштабе как по размерам, так и по цвету.

1. Ближайший: найдите ближайшую известную точку в сетке (quick’n’dirty)

2. Линейный: возьмите куб 2x2x2 и используйте линейную интерполяцию для вычисления результата; результаты всегда находятся между существующими точками данных, никаких хитростей с краями не требуется. Однако первые различия не являются гладкими, что видно на изображении и иногда является проблемой при обработке сигналов.

3. Кубическая: возьмите куб 4x4x4 и используйте кубическую интерполяцию для вычисления результата. Результат имеет непрерывные первые дифференциалы, т. Е. результат «гладкий» во всех направлениях. Могут быть краевые эффекты из-за большого объема выборки. (Алгоритм интерполяции повторяет самые крайние вокселы, но они могут быть зеркальными, приниматься за некоторую константу и т.д.)

Различия в природе этих методов интерполяции могут быть еще лучше видны, если мы проведем отрезок линии через блок:

  # line segment from (-12,-12,-12) to (12,12,12)
 v = np.array([np.linspace(-12,12,1000)]*3).T

 res = [ tri_interp(M, interp_catmull_rom, p[0], p[1], p[2], m0,m1,n0,n1,o0,o1) for p in v ]
  

Для различных функций интерполяции мы имеем следующие значения вдоль отрезка линии:

Блочность «ближайшего» и перегибы в линейной интерполяции хорошо видны. Может показаться удивительным, что «линейная» интерполяция не дает линейных результатов. Это лучше всего понять, представив 2D-интерполяцию с квадратом со значениями углов 1,0,1,0 (т. Е. диагонально противоположные углы имеют одинаковые значения). Нет способа интерполировать это так, чтобы все разделы были линейными.

Приведенный здесь алгоритм является очень тривиальным алгоритмом, который можно ускорить с помощью нескольких трюков. Особенно, если в одной ячейке требуется несколько точек, подходы на основе матрицы намного быстрее. К сожалению, не существует очень быстрых методов, а методы на основе матрицы сложнее объяснить.

Один из способов, который я нашел, это ndimage.map_coordinates . Немного жаль, что он скрыт и не связан в модуле interpolate, но он работает удивительно хорошо и предлагает множество порядков приближения.

Два больших недостатка:

• это немного подозрительно на концах сетки, mode='nearest' заставит его выбрать значение границы, в то время как в противном случае оно принимало бы значения 0 для точек за пределами сетки.
• это не очень эффективно.

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

 import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')
Y = N M O
n,m,o = 0.01, 103, 1007
from scipy import ndimage
data = Y
# subtract lower end of grid
coords = np.array([[n-0,m-100,o-1000]])
ndimage.map_coordinates(data, coords.T, order=2, mode='nearest')