用于平滑零区域的 Numpy 过滤器

Numpy filter to smooth out zero-regions

我有一个 0 和更大整数的 2D numpy 数组,其中的值代表区域标签。例如,

array([[9, 9, 9, 0, 0, 0, 0, 1, 1, 1],
       [9, 9, 9, 9, 0, 7, 1, 1, 1, 1],
       [9, 9, 9, 9, 0, 2, 2, 1, 1, 1],
       [9, 9, 9, 8, 0, 2, 2, 1, 1, 1],
       [9, 9, 9, 8, 0, 2, 2, 2, 1, 1],
       [4, 4, 4, 4, 0, 2, 2, 2, 1, 1],
       [4, 6, 6, 4, 0, 0, 0, 0, 0, 0],
       [4, 6, 6, 4, 0, 0, 0, 0, 0, 0],
       [4, 4, 4, 4, 5, 5, 5, 5, 5, 5],
       [4, 4, 4, 4, 5, 5, 5, 5, 5, 5]])

我希望等于 0 的索引(即零区域)取其附近最常见的值。该操作将实质上关闭零区域。我已经尝试了膨胀、腐蚀、grey-closing, and other morphology operations 的多种变体,但我无法完全消除零区域(不笨拙地混合其他区域)。一个体面的方法可能是定义一个仅对零进行卷积的内核,并使用过滤器区域中最常见的标签设置值。不过我不确定如何实现它。

基于卷积思路的可能解法

from scipy import stats
ar = #Your np array
blank = np.zeros(ar.shape)
#Size to search in for mode values
window_size = 3

for x,y in np.array(np.where(ar == 0)).T:
    window = ar[max(x-window_size,0):x+window_size,max(0,y-window_size):y+window_size]
    oneD = window.flatten()

    #fill blank array with modal value
    blank[x,y] = stats.mode(oneD[oneD != 0])[0]

#fill in the zeros
print ar + blank

我不确定这里是否可以避免循环

这是一个使用 Numba 的可行解决方案,我还没有介绍它,但应该非常快:

import numba
@numba.njit
def nn(arr):
    res = arr.copy()
    zeros = np.where(arr == 0)
    for n in range(len(zeros[0])):
        i = zeros[0][n]
        j = zeros[1][n]
        left = max(i-1, 0)
        right = min(i+2, arr.shape[1])
        top = max(j-1, 0)
        bottom = min(j+2, arr.shape[0])
        area = arr[left:right,top:bottom].ravel()
        counts = np.bincount(area[area != 0])
        res[i,j] = np.argmax(counts)
    return res

它产生:

array([[9, 9, 9, 9, 7, 1, 1, 1, 1, 1],
       [9, 9, 9, 9, 9, 7, 1, 1, 1, 1],
       [9, 9, 9, 9, 2, 2, 2, 1, 1, 1],
       [9, 9, 9, 8, 2, 2, 2, 1, 1, 1],
       [9, 9, 9, 8, 2, 2, 2, 2, 1, 1],
       [4, 4, 4, 4, 2, 2, 2, 2, 1, 1],
       [4, 6, 6, 4, 4, 2, 2, 2, 1, 1],
       [4, 6, 6, 4, 4, 5, 5, 5, 5, 5],
       [4, 4, 4, 4, 5, 5, 5, 5, 5, 5],
       [4, 4, 4, 4, 5, 5, 5, 5, 5, 5]])

此处的内核大小为 3x3,定义为 ij 减去 1 并加 2(加 2 是因为 Python 切片需要最后一个,例如 [0:3] 给你 3 个元素)。边界条件由 minmax.

处理

bincount 想法的功劳:

这里提出了一种矢量化方法。步骤是:

  1. 获取内核大小的 2D 滑动 windows,导致 4D 数组。我们可以用 skimage's view_as_windows 将它们作为视图,从而避免创建 任何额外的内存。

  2. Select windows 通过索引到 4D 数组以零为中心。这会强制复制。但是假设零的数量比输入数组中的元素总数相对较小,这应该没问题。

  3. 对于每一个被选中的windows,将每个window偏移一个合适的偏移量,以使用np.bincount来执行计数。因此,使用 bincount 并获得不包括零的最大计数。最大计数的argmax应该是我们的家伙!

这是涵盖这些步骤的实现 -

from skimage.util import view_as_windows as viewW

def fill_zero_regions(a, kernel_size=3):
    hk = kernel_size//2 # half_kernel_size    

    a4D = viewW(a, (kernel_size,kernel_size))
    sliced_a = a[hk:-hk,hk:-hk]
    zeros_mask = sliced_a==0
    zero_neighs = a4D[zeros_mask].reshape(-1,kernel_size**2)
    n = len(zero_neighs) # num_zeros

    scale = zero_neighs.max()+1
    zno = zero_neighs + scale*np.arange(n)[:,None] # zero_neighs_offsetted

    count = np.bincount(zno.ravel(), minlength=n*scale).reshape(n,-1)
    modevals = count[:,1:].argmax(1)+1
    sliced_a[zeros_mask] = modevals
    return a

样本运行-

In [23]: a
Out[23]: 
array([[9, 9, 9, 0, 0, 0, 0, 1, 1, 1],
       [9, 9, 9, 9, 0, 7, 1, 1, 1, 1],
       [9, 9, 9, 9, 0, 2, 2, 1, 1, 1],
       [9, 9, 9, 8, 0, 2, 2, 1, 1, 1],
       [9, 9, 9, 8, 0, 2, 2, 2, 1, 1],
       [4, 4, 4, 4, 0, 2, 2, 2, 1, 1],
       [4, 6, 6, 4, 0, 0, 0, 0, 0, 0],
       [4, 6, 6, 4, 0, 0, 0, 0, 0, 0],
       [4, 4, 4, 4, 5, 5, 5, 5, 5, 5],
       [4, 4, 4, 4, 5, 5, 5, 5, 5, 5]])

In [24]: fill_zero_regions(a)
Out[24]: 
array([[9, 9, 9, 0, 0, 0, 0, 1, 1, 1],
       [9, 9, 9, 9, 9, 7, 1, 1, 1, 1],
       [9, 9, 9, 9, 2, 2, 2, 1, 1, 1],
       [9, 9, 9, 8, 2, 2, 2, 1, 1, 1],
       [9, 9, 9, 8, 2, 2, 2, 2, 1, 1],
       [4, 4, 4, 4, 2, 2, 2, 2, 1, 1],
       [4, 6, 6, 4, 4, 2, 2, 2, 1, 0],
       [4, 6, 6, 4, 4, 5, 5, 5, 5, 0],
       [4, 4, 4, 4, 5, 5, 5, 5, 5, 5],
       [4, 4, 4, 4, 5, 5, 5, 5, 5, 5]])

正如所见,我们没有解决边界情况。如果需要这样做,请使用零填充数组作为输入数组,如下所示:np.pad(a, (k//2,k//2), 'constant'),其中 k 作为内核大小(样本为 =3)。