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

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

Теперь Scipy поддерживает это просто превосходно: я фильтрую объем с помощью

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 (xh, y, z)) / (2 * h) . (Да, я знаю, что могу уменьшить количество дополнительных нажатий до 3 и сделать «односторонние» различия, но асимметрия меня бы раздражала.)

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

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