【问题标题】:2d convolution in python with missing datapython中缺少数据的二维卷积
【发布时间】:2021-04-06 17:06:57
【问题描述】:

我知道有 scipy.signal.convolve2d 函数来处理 2d numpy 数组的二维卷积,并且有 numpy.ma 模块来处理丢失的数据,但这两种方法似乎彼此不兼容(意味着即使您在 numpy 中屏蔽 2d 数组,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

卷积的期望结果(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 -

感谢Aguy的建议,这是帮助卷积后计算结果的好方法。现在假设我们可以从 Array.mask 中获取 Array 的掩码,这会给我们一个结果

                   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

如何使用此掩码将卷积后的结果转换为掩码数组?

【问题讨论】:

  • 卷积应该如何处理缺失值?我觉得Result[0,1] 应该是0 而不是-...
  • 您想将掩码值替换为 0,然后进行卷积,然后将掩码重新应用于结果。

标签: python numpy scipy convolution


【解决方案1】:

我不认为用 0 替换是这样做的正确方法,您正在将 covolved 值推向 0。这些缺失应该被视为“缺失”。因为它们代表了缺失的信息片段,没有理由仅仅假设它们可能是 0,它们根本不应该参与任何计算。

我尝试将缺失值设置为numpy.nan,然后进行卷积,结果表明内核之间的任何重叠和任何缺失都会在结果中给出nan,即使重叠是内核中的0,所以你会在结果中得到一个扩大的缺失洞。根据您的应用,这可能是您想要的结果。

但在某些情况下,您不希望仅仅因为 1 次缺失就丢弃这么多信息(也许 astropy 具有更好的实现:numpy.nans 被忽略(或替换为插值?)。

因此,使用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)

但是,您仍然无法控制可以容忍的缺失程度。为了实现这一点,我创建了一个函数,它使用scipy.ndimage.convolve() 进行初始卷积,但只要涉及缺失(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

下图展示了输出。左边是一个 30x30 的随机地图,有 3 个numpy.nan 洞,大小为:

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

右侧是卷积输出,由 5x5 内核(全为 1)和 50% 的容差水平 (max_missing=0.5)。

所以前 2 个较小的孔使用附近的值填充,而在最后一个中,因为缺失的数量 > 0.5x5x5 = 12.5numpy.nans 被放置来表示缺失的信息。

【讨论】:

  • 请查看这个 repo (github.com/Xunius/Random-Fortran-scripts),我在其中创建了一个从 Fortran 编译的 python 模块,这将大大加快速度。
  • 我更新了自己的解决方案,在这里github.com/Xunius/python_convolution
  • 我对@Jason 的解决方案有两个问题:您能否发布您如何创建/使用您的 fxn 的代码?你为什么不使用 convolve2d,我怎么能使用它的属性,比如不同的模式,例如“有效”,即不保持图像的相同形状?
【解决方案2】:

我发现了一个黑客。而不是 nan 使用虚数(将其更改为 1i)运行卷积并设置虚值高于阈值的任何地方都是 nan。每当它在下面时,就取真正的价值。这是一个代码sn-p:

frames_complex = np.zeros_like(frames_, dtype=np.complex64)
frames_complex[np.isnan(frames_)] = np.array((1j))
frames_complex[np.bitwise_not(np.isnan(frames_))] =                         
frames_[np.bitwise_not(np.isnan(frames_))]
convolution = signal.convolve(frames_complex, gaussian_window, 'valid')
convolution[np.imag(convolution)>0.2] = np.nan
convolution = convolution.astype(np.float32)

【讨论】:

  • 制作滚动窗口功能的绝妙技巧。
  • 这不是把 nan 的实际值归零了吗?
  • 我想直接卷积方法不使用虚数,所以幸运的是没问题。但是如果将 FFT 用于卷积,这会混淆吗?
猜你喜欢
  • 2013-10-16
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2017-02-24
  • 2018-04-09
  • 2015-04-11
  • 2011-01-27
  • 2010-11-09
相关资源
最近更新 更多