2d畳み込みPythonでデータが欠落している

Question

2次元のnumpy配列の2次元畳み込みを処理するscipy.signal.convolve2d関数があり、欠落したデータを処理するnumpy.maモジュールがありますが、これらの2つの方法は見えません（つまり、numpyで2次元配列をマスクしても、convolve2dのプロセスは影響を受けません）。 numpyとscipyパッケージのみを使用してコンボリューションで欠損値を処理する方法はありますか？例えば

：畳み込み用

  1 - 3 4 5 
      1 2 - 4 5 
    Array = 1 2 3 - 5 
      - 2 3 4 5 
      1 2 3 4 - 

    Kernel = 1 0 
      0 -1 

所望の結果（アレイ、カーネル、境界= 'ラップ'）：Aguyから提案ため

   -1 - -1 -1 4 
       -1 -1 - -1 4 
    Result = -1 -1 -1 - 5 
       - -1 -1 4 4 
       1 -1 -1 -1 - 

---

おかげで、すなわち畳み込み後の結果の計算を助ける本当に良い方法です。今度は、私たちに

    False True False False False      
        False False True False False 
    Array.mask == False False False True False 
        True False False False False 
        False False False False True 

どのようにマスクされた配列にコンボリューション後に結果を変換するために、このマスクを使用することができますか？の結果を与えることになる、我々はArray.maskからアレイのマスクを得ることができるとしましょうか

Answer 1

これを行う正しい方法は0と置き換えてはいけないと思います.0に向かって積み重ねられた値を微調整しています。これらのミスは、文字通り「ミッシング」として扱われるはずです。彼らは欠けている情報を表しており、0であると仮定するだけの正当性はないので、計算には全く関与してはいけません。

私はnumpy.nanに欠損値を設定しようとしたし、その後、畳み込み、結果はあなたが得るので、カーネルと欠落している間のいずれかの重複が、重複がカーネルで0である場合でも、結果にnanを与えることを示唆しています結果のミスの拡大した穴。アプリケーションによっては、これが望ましい結果になる可能性があります。

しかし、場合によっては、1つの不足分の情報だけを捨てることは望ましくありません（おそらく、< =行方不明の50％は許容されます）。このような場合、astropyの実装がより良い別のモジュールが見つかりました。numpy.nanは無視されます（または補間された値で置き換えられますか？）。

のでastropyを使用するには、次の操作を行います。

from astropy.convolution import convolve 
inarray=numpy.where(inarray.mask,numpy.nan,inarray) # masking still doesn't work, has to set to numpy.nan 
result=convolve(inarray,kernel) 

しかし、それでもまだ、あなたは許容されるどのくらい不足しているのを制御することはできません。 missings（numpy.nan）が関与している時はいつでも、私は最初のコンボリューションのためscipy.ndimage.convolve()を使用する関数を作成しますが、手動で再計算値ました、ことを達成するために：以下

def convolve2d(slab,kernel,max_missing=0.5,verbose=True): 
    '''2D convolution with missings ignored 

    <slab>: 2d array. Input array to convolve. Can have numpy.nan or masked values. 
    <kernel>: 2d array, convolution kernel, must have sizes as odd numbers. 
    <max_missing>: float in (0,1), max percentage of missing in each convolution 
        window is tolerated before a missing is placed in the result. 

    Return <result>: 2d array, convolution result. Missings are represented as 
        numpy.nans if they are in <slab>, or masked if they are masked 
        in <slab>. 

    ''' 

    from scipy.ndimage import convolve as sciconvolve 

    assert numpy.ndim(slab)==2, "<slab> needs to be 2D." 
    assert numpy.ndim(kernel)==2, "<kernel> needs to be 2D." 
    assert kernel.shape[0]%2==1 and kernel.shape[1]%2==1, "<kernel> shape needs to be an odd number." 
    assert max_missing > 0 and max_missing < 1, "<max_missing> needs to be a float in (0,1)." 

    #--------------Get mask for missings-------------- 
    if not hasattr(slab,'mask') and numpy.any(numpy.isnan(slab))==False: 
     has_missing=False 
     slab2=slab.copy() 

    elif not hasattr(slab,'mask') and numpy.any(numpy.isnan(slab)): 
     has_missing=True 
     slabmask=numpy.where(numpy.isnan(slab),1,0) 
     slab2=slab.copy() 
     missing_as='nan' 

    elif (slab.mask.size==1 and slab.mask==False) or numpy.any(slab.mask)==False: 
     has_missing=False 
     slab2=slab.copy() 

    elif not (slab.mask.size==1 and slab.mask==False) and numpy.any(slab.mask): 
     has_missing=True 
     slabmask=numpy.where(slab.mask,1,0) 
     slab2=numpy.where(slabmask==1,numpy.nan,slab) 
     missing_as='mask' 

    else: 
     has_missing=False 
     slab2=slab.copy() 

    #--------------------No missing-------------------- 
    if not has_missing: 
     result=sciconvolve(slab2,kernel,mode='constant',cval=0.) 
    else: 
     H,W=slab.shape 
     hh=int((kernel.shape[0]-1)/2) # half height 
     hw=int((kernel.shape[1]-1)/2) # half width 
     min_valid=(1-max_missing)*kernel.shape[0]*kernel.shape[1] 

     # dont forget to flip the kernel 
     kernel_flip=kernel[::-1,::-1] 

     result=sciconvolve(slab2,kernel,mode='constant',cval=0.) 
     slab2=numpy.where(slabmask==1,0,slab2) 

     #------------------Get nan holes------------------ 
     miss_idx=zip(*numpy.where(slabmask==1)) 

     if missing_as=='mask': 
      mask=numpy.zeros([H,W]) 

     for yii,xii in miss_idx: 

      #-------Recompute at each new nan in result------- 
      hole_ys=range(max(0,yii-hh),min(H,yii+hh+1)) 
      hole_xs=range(max(0,xii-hw),min(W,xii+hw+1)) 

      for hi in hole_ys: 
       for hj in hole_xs: 
        hi1=max(0,hi-hh) 
        hi2=min(H,hi+hh+1) 
        hj1=max(0,hj-hw) 
        hj2=min(W,hj+hw+1) 

        slab_window=slab2[hi1:hi2,hj1:hj2] 
        mask_window=slabmask[hi1:hi2,hj1:hj2] 
        kernel_ij=kernel_flip[max(0,hh-hi):min(hh*2+1,hh+H-hi), 
            max(0,hw-hj):min(hw*2+1,hw+W-hj)] 
        kernel_ij=numpy.where(mask_window==1,0,kernel_ij) 

        #----Fill with missing if not enough valid data---- 
        ksum=numpy.sum(kernel_ij) 
        if ksum<min_valid: 
         if missing_as=='nan': 
          result[hi,hj]=numpy.nan 
         elif missing_as=='mask': 
          result[hi,hj]=0. 
          mask[hi,hj]=True 
        else: 
         result[hi,hj]=numpy.sum(slab_window*kernel_ij) 

     if missing_as=='mask': 
      result=numpy.ma.array(result) 
      result.mask=mask 

    return result 

すると、出力を実証する図です。左側上のサイズの3つのnumpy.nan穴を有する30×30ランダムマップされる：右側に

1. 1x1の
2. の3x3
3. 5x5の

は、畳み込み出力します5x5カーネル（すべて1）、および許容レベル50％（max_missing=0.5）です。

したがって、最初の2つの小さな穴は、近くの値を使用して塗りつぶされます。最後の1つは、ミッシングの数>0.5x5x5 = 12.5,numpy.nanが不足している情報を表すために配置されます。