2016-07-12 3 views
1

2d numpy 배열에 2 차원 컨볼 루션을 처리하는 scipy.signal.convolve2d 함수가 있으며, 누락 된 데이터를 처리하기위한 numpy.ma 모듈이 있지만이 두 가지 메소드가없는 것 같습니다. 서로 호환 될 수 있습니다 (즉, numpy로 2 차원 배열을 마스크해도 convolve2d의 프로세스는 영향을받지 않습니다). numpy 및 scipy 패키지 만 사용하여 회선 누락 값을 처리 할 수있는 방법이 있습니까? 예를 들어파이썬에서 누락 된 데이터가있는 2 차원 컨볼 루션

:

  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 - 

덕분에, 즉 인 convolution 후에 결과 계산을 돕는 정말 좋은 방법. 이제 어떤이 어떻게 마스크 배열로 회선 후 결과를 변환하려면이 마스크를 사용하여 우리에게

    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에서 배열의 마스크를 얻을 수 가정 해 봅시다?

+0

누락 값을 처리하기 위해 회선은 어떻게되어 있습니까? '결과 [0,1]'이'-'보다는'0 '이되어야한다고 생각합니다 ... – Julien

+1

마스크 된 값을 0으로 대체 한 다음, 다시 컨벌루션하고 마스크를 결과에 다시 적용하려고합니다. – Aguy

답변

1

나는 0으로 대체하는 것이 올바른 방법이라고 생각하지 않는다. 당신은 0을 향한 값을 누그러 뜨리고있다. 이러한 누락은 문자 그대로 "누락"으로 처리되어야한다. 그들은 누락 된 정보를 나타 내기 때문에 0 일 수 있다고 가정하는 것은 정당화되지 않으며 전혀 계산에 관여해서는 안됩니다.

누락 된 값을 numpy.nan으로 설정 한 다음 convolve를 시도한 결과 커널과 누락 사이의 겹침은 결과가 nan이고 커널에서 0이 겹친 경우에도 마찬가지입니다. 그 결과에 놓친 구멍을 확대. 응용 프로그램에 따라 원하는 결과가 될 수 있습니다.

그러나 경우에 따라 1 개의 누락 된 부분에 대해 너무 많은 정보를 삭제하고 싶지 않을 수도 있습니다 (어쩌면 < = 누락 된 부분의 50 %는 여전히 허용됩니다). 이 경우 더 좋은 구현을 가진 astropy 모듈을 발견했습니다 : numpy.nan은 무시됩니다 (또는 보간 된 값으로 대체 되었습니까?).

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) 

을하지만 여전히, 당신은 허용하는 방법을 누락 많이 제어 할 수 없습니다 : 그래서 astropy를 사용

, 다음을 수행 할 것입니다. 내가 초기 회선에 대한 scipy.ndimage.convolve()를 사용하는 함수를 만든 것을 달성하지만, 수동에서 missing (numpy.nan)를 포함 할 때마다 값 - 계산을 다시하려면

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 구멍이 30x30 랜덤 맵 :

  1. 의 1x1
  2. × 3
  3. 5 ×

오른쪽 enter image description here

컨 볼빙 된 (convolved) 출력은이 5x5 커널 (모두 1), 허용 수준은 50 % (max_missing=0.5)입니다.

누락 된 숫자가>0.5x5x5 = 12.5, numpy.nan 인 누락 된 정보를 나타 내기 때문에 처음 2 개의 작은 구멍은 가까운 값을 사용하여 채워집니다.

+0

포트란에서 컴파일 된 파이썬 모듈을 만든이 저장소 (https://github.com/Xunius/Random-Fortran-scripts)도 확인해보십시오. – Jason

+0

그리고 여기 내 자신의 솔루션을 업데이트했습니다. https://github.com/Xunius/python_convolution – Jason