26

2D 座標のリストと 3 番目の変数 (速度) から、サンプリングされた領域全体をカバーする 2D numpy 配列を作成しました。私の意図は、各ピクセルがその中にあるポイントの平均速度を含む画像を作成することです。その後、その画像をガウス フィルターでフィルター処理します。

問題は、領域が均一にサンプリングされていないことです。Nanしたがって、画像の中央に情報 ( ) のないピクセルがいくつかあります。配列をガウス フィルターでフィルター処理しようとすると、Nan伝播が画像全体を台無しにします。

この画像をフィルタリングする必要がありますが、情報のないすべてのピクセルを拒否します。つまり、ピクセルに情報が含まれていない場合は、フィルタリングで考慮されません。

平均化のための私のコードの例を次に示します。

Mean_V = np.zeros([len(x_bins), len(y_bins)])

for i, x_bin in enumerate(x_bins[:-1]):
    bin_x = (x > x_bins[i]) & (x <= x_bins[i+1])
    for j, y_bin in enumerate(y_bins[:-1]):
        bin_xy = (y[bin_x] > y_bins[j]) & (y[bin_x] <= y_bins[j+1])
        if (sum(x > 0 for x in bin_xy) > 0) :
            Mean_V[i,j]=np.mean(V[bin_x][bin_xy])
        else:
            Mean_V[i,j]=np.nan

編集：

2013 年に作成したこの質問に終わった Web サーフィン。この問題の解決策は、astropy ライブラリにあります。

http://docs.astropy.org/en/stable/convolution/

Astropy の畳み込みは、NaN ピクセルを隣接ピクセルからのカーネル加重補間に置き換えます。

ありがとうございます！

4

3 に答える 3

37

言葉で：

特定の配列Uの NaN を無視するガウス フィルターは、標準のガウス フィルターを 2 つの補助配列VとWに適用し、2 つの比率をとって結果Zを取得することで簡単に取得できます。

ここで、Vは元のUのコピーでNaN がゼロに置き換えられ、Wは元のUの NaN の位置を示すゼロを含む 1 の配列です。

NaN をゼロで置き換えると、フィルター処理された配列にエラーが発生しますが、同じガウス フィルターを別の補助配列に適用し、2 つを組み合わせることで補償できます。

Python で:

import numpy as np
import scipy as sp
import scipy.ndimage

sigma=2.0                  # standard deviation for Gaussian kernel
truncate=4.0               # truncate filter at this many sigmas

U=sp.randn(10,10)          # random array...
U[U>2]=np.nan              # ...with NaNs for testing

V=U.copy()
V[np.isnan(U)]=0
VV=sp.ndimage.gaussian_filter(V,sigma=sigma,truncate=truncate)

W=0*U.copy()+1
W[np.isnan(U)]=0
WW=sp.ndimage.gaussian_filter(W,sigma=sigma,truncate=truncate)

Z=VV/WW

数字で：

ここでは、デモンストレーションのためにガウス フィルターの係数を [0.25,0.50,0.25] に設定し、一般性を失うことなく合計すると 0.25+0.50+0.25=1 になります。

NaN をゼロで置き換え、ガウス フィルターを適用した後 (以下の VV を参照)、ゼロがエラーを導入することは明らかです。したがって、「真の」値を過小評価します。

ただし、これは、2 番目の補助配列 (以下の WW を参照) を使用することで補償できます。この配列は、同じガウス分布でフィルター処理した後、係数の合計のみを含みます。

したがって、フィルター処理された 2 つの補助配列を分割すると、合計が 1 になるように係数が再スケーリングされ、NaN の位置は無視されます。

array U         1   2   NaN 1   2    
auxiliary V     1   2   0   1   2    
auxiliary W     1   1   0   1   1
position        a   b   c   d   e

filtered VV_b   = 0.25*V_a  + 0.50*V_b  + 0.25*V_c
                = 0.25*1    + 0.50*2    + 0
                = 1.25

filtered WW_b   = 0.25*W_a  + 0.50*W_b  + 0.25*W_c
                = 0.25*1    + 0.50*1    + 0
                = 0.75

ratio Z         = VV_b / WW_b  
                = (0.25*1 + 0.50*2) / (0.25*1    + 0.50*1)
                = 0.333*1 + 0.666*2
                = 1.666

更新 - ゼロ除算:

以下は、以下のコメントからの @AndyL および @amain による有用な質問と回答を組み込んでいます。

NaN の大きな領域は、ガウス カーネルのサポート内に NaN エントリしかない場合、一部の位置でゼロ分母 (WW=0) になる可能性があります (理論的にはサポートは無限ですが、実際にはカーネルは通常切り捨てられます。上記のコード例の truncate' パラメータ)。その状況では、指定子もゼロ (VV=0) になるため、numpy は「RuntimeWarning: true_divide で検出された無効な値」をスローし、対応する位置で NaN を返します。

これはおそらく最も一貫性のある意味のある結果であり、numpy の警告に耐えられる場合は、それ以上の調整は必要ありません。

于 2016-03-30T11:21:06.540 に答える