6

スカラー値の大きな 3D numpy 配列があります (必要に応じて「ボリューム」と呼んでください)。私は、一連の不規則な、前もって知られていないすべての非整数 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 つの追加ポイントでボリュームをサンプリングし (これは少なくとも を 1 回呼び出すだけで実行できますmap_coordinates)、たとえば から x 導関数を計算し(i(x+h,y,z)-i(x-h,y,z))/(2*h)ます。(はい、追加のタップ数を 3 に減らして「片側」の違いを実行できることはわかっていますが、非対称性は私を悩ませます。)

私の本能は、勾配を取得するためのより直接的な方法があるはずですが、それを理解したり、Scipy 実装の内部で何が起こっているかを理解するのに十分なスプライン数学を (まだ) 知りません: scipy/scipy/ndimage/src/ni_interpolation.c.

中央差分よりも「より直接的に」勾配を取得するより良い方法はありますか? できれば、Scipy の内部をハッキングするのではなく、既存の機能を使用してそれらを取得できるようにするものです。

4

1 に答える 1

1

Aha: numpy コードで引用されているスプラインに関する古典的な論文によると、次数 n のスプラインとその導関数は次のように関連付けられています。

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

したがって、SciPy のスプライン補間を使用して、低次の事前フィルター処理されたボリュームも維持し、導関数ごとに数回クエリを実行することで導関数を取得できました。これは、かなりの量のメモリを追加することを意味します (キャッシュの「メイン」ボリュームとの競合かもしれません) が、おそらく低次スプラインの評価の方が高速であるため、中央の差分よりも全体的に高速かどうかは明らかではありません。私が現在行っている小さなオフセットを使用しています。まだ試していません。

于 2011-08-16T13:50:55.013 に答える