2D Faltung in Python mit fehlenden Daten

Question

Ich weiß, es gibt scipy.signal.convolve2d Funktion, um 2-dimensionale Convolution für 2d numpy Array zu behandeln, und es gibt numpy.ma-Modul, um fehlende Daten zu behandeln, aber diese beiden Methoden scheinen nicht um miteinander kompatibel zu sein (was bedeutet, dass selbst wenn Sie ein 2d-Array in numpy maskieren, der Prozess in convolve2d nicht beeinflusst wird). Gibt es eine Möglichkeit, fehlende Werte in der Faltung mit nur numpy und scipy Paketen zu behandeln?

Zum Beispiel:

  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 

Wunschergebnis für Faltung (Array, kernel boundary = 'wrap'):

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

---

Danke für die Anregung von Aguy, das ist ein wirklich guter Weg, um die Berechnung des Ergebnisses nach der Faltung zu helfen. Lassen Sie uns jetzt sagen, dass wir die Maske des Array von Array.mask bekommen können, die uns ein Ergebnis von

    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 

Wie geben würde, kann ich diese Maske verwenden das Ergebnis nach Faltung in eine maskierte Array zu konvertieren?

Answer 1

Ich glaube nicht, das Ersetzen mit 0s ist die richtige Art, dies zu tun, Sie stoßen die coolved Werte auf 0 an. Diese Missings sollten buchstäblich als "Missings" behandelt werden. Weil sie die fehlenden Informationen repräsentieren und es nicht gerechtfertigt ist, anzunehmen, dass sie 0 sein könnten, und sie sollten überhaupt nicht in irgendeine Berechnung einbezogen werden.

Ich hat versucht, fehlende Werte numpy.nan Einstellung und dann falten, das Ergebnis deutet darauf hin, dass Überschneidungen zwischen dem Kernel und jedem fehlenden einem nan im Ergebnis ergeben, auch wenn die Überlappung mit einem 0 im Kernel ist, so dass Sie eine bekommen vergrößertes Loch von Missings im Ergebnis. Abhängig von Ihrer Anwendung könnte dies das gewünschte Ergebnis sein.

Aber in einigen Fällen möchten Sie nicht so viele Informationen nur für 1 einzelne zu verwerfen (vielleicht < = 50% der fehlenden ist immer noch tolerierbar). In solchen Fällen habe ich ein anderes Modul astropy gefunden, das eine bessere Implementierung hat: numpy.nan s werden ignoriert (oder durch interpolierte Werte ersetzt?).

So astropy verwenden, würden Sie wie folgt vor:

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) 

Aber noch, Sie haben keine Kontrolle darüber, wie viel von tolerierbaren fehlt. Um dies zu erreichen, habe ich eine Funktion erstellt, die die scipy.ndimage.convolve() für die anfängliche Faltung verwendet, sondern manuell neu berechnen Werte, wenn Missings (numpy.nan) beteiligt sind:

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 

Im Folgenden eine Figur ist, um die Ausgabe zu demonstrieren.Auf der linken Seite mit 3 numpy.nan Löcher mit Größen von einem 30x30 zufällige Karte ist:

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

Auf der rechten Seite ist die convolved Ausgang, durch einen 5x5 Kernel (mit allen 1s) und eine Toleranz von 50% (max_missing=0.5).

Also die ersten 2 kleineren Löcher werden mit nahe gelegenen Werten gefüllt, und in der letzten, weil die Anzahl der Missings>0.5x5x5 = 12.5, numpy.nan s platziert sind, um fehlende Informationen darzustellen.