Как получить градиент скиппи интерполяции напрямую?

Question

У меня большой массив трехмерных скалярных значений (ОК, если нужно, назовите его «объем»). Я хочу интерполировать гладкое скалярное поле над этим в последовательности нерегулярных, не все известные заранее нецелые координаты XYZ.

Теперь поддержка Сципи просто превосходна: я фильтрую громкость с помощью

filtered_volume = scipy.ndimage.interpolation.spline_filter(volume)

и вызвать

scipy.ndimage.interpolation.map_coordinates(
    filtered_volume,
    [[z],[y],[x]],
    prefilter=False)

для (x, y, z) интереса для получения, по-видимому, хорошо поведенческих (гладких и т. Д.) Интерполированных значений.

Пока все хорошо. Однако моему приложению также нужны локальные производные интерполированного поля. В настоящее время я получаю их с помощью центральной разности: я также пробую объем в 6 дополнительных точках (это можно сделать по крайней мере одним вызовом map_coordinates) и вычисляю, например, производную x от (i(x+h,y,z)-i(x-h,y,z))/(2*h). (Да, я знаю, что мог бы уменьшить количество дополнительных отводов до 3 и сделать «односторонние» различия, но асимметрия меня раздражала.)

Мой инстинкт заключается в том, что должен быть более прямой способ получения градиента но я не знаю достаточно математики сплайнов (пока), чтобы понять это или понять, что происходит в сущности реализации Scipy: scipy/scipy/ndimage/src/ni_interpolation.c.

Есть ли лучший способ получения моих градиентов "более напрямую", чем центральное дифференцирование? Предпочтительно тот, который позволяет получить их, используя существующую функциональность, а не взламывая внутренности Сципи.

Answer 1

Ага: согласно классической статье о сплайнах , процитированной в числовом коде, сплайны порядка n и их производные связаны

   n          n-1           n-1
 dB (x)/dx = B   (x+1/2) - B   (x-1/2)

Таким образом, используя сплайн-интерполяцию SciPy, я мог бы получить мои производные, также поддерживая предварительно фильтрованный объем более низкого порядка и запрашивая его пару раз на производную. Это означает добавление достаточного объема памяти (возможно, конкуренция с «основным» томом для кэша), но, по-видимому, оценка сплайнов более низкого порядка выполняется быстрее, поэтому для меня не очевидно, будет ли он быстрее или не в целом, чем центральный дифференциал используя небольшие смещения, которые я делаю в настоящее время. Еще не пробовал.