Jak uzyskać gradient interpolantu scipy bezpośrednio?

Question

Mam dużą trójwymiarową tablicę wartości skalarnych (OK, nazwij to "tom", jeśli musisz). Chcę interpolować gładkie pole skalarne nad tym ciągiem nieregularnych, nie wszystkich znanych nie-integralnych współrzędnych xyz z góry.

wsparcie

& Pobierz scipy za to po prostu doskonałe: filtrować głośność

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

i powołać

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

dla (x, y, z) odsetek do uzyskania widocznie dobrze wychowane (gładkie itp.) wartości interpolowane.

Jak dotąd tak dobrze. Jednak moja aplikacja potrzebuje również lokalnych pochodnych pola interpolowanego. Obecnie uzyskuję je przez centralne różnicowanie: próbuję również objętości w 6 dodatkowych punktach (można to zrobić za pomocą tylko jednego połączenia z map_coordinates) i obliczyć np. Pochodną x z (i(x+h,y,z)-i(x-h,y,z))/(2*h). (Tak, wiem, że mogłem zmniejszyć liczbę dodatkowych odczepów do 3 i zrobić "jednostronne" różnice, ale irytowała mnie asymetria.)

Instynktem jest to, że powinien istnieć bardziej bezpośredni sposób uzyskiwania gradientu ale nie znam wystarczającej matematyki splajnowej (jeszcze), aby to zrozumieć, lub zrozumieć, co jest dzieje się w odwzorowaniach implementacji Scipy: scipy/scipy/ndimage/src/ni_interpolation.c.

Czy istnieje lepszy sposób uzyskania moich gradientów "bardziej bezpośrednio" niż centralne różnicowanie? Najlepiej taki, który pozwala na ich uzyskanie przy użyciu istniejącej funkcjonalności, a nie na hackowaniu na instrumencie Scipy.

Answer 1

Aha: Zgodnie z classic paper on splines cytowanego w kodzie numpy, wypusty porządku N i ich pochodne są związane przez

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

Więc za pomocą interpolacji spline scipy za I mógłby dostać moje pochodne również poprzez utrzymanie niższej - wolumin z prefiltrem i zapytania, które kilka razy na pochodną. Oznacza to dodanie sporej ilości pamięci (może konkurencji z "głównym" woluminem dla pamięci podręcznej), ale prawdopodobnie ocena splajnów niższego rzędu jest szybsza, więc nie jest dla mnie jasne, czy byłaby szybsza, czy nie ogólna niż centralna różnica. używając małych przesunięć, które robię obecnie. Jeszcze tego nie próbowałem.