Numpy自定义移动窗口计数范围内的单元格

Question

我有一个 2d numpy 数组，我想从中知道落在两个堆叠内核中指定的两个值之间的所有像素的计数。边界应从计数中排除。

例如：

input = np.array([[1, 0, 0, 2, 1],
   [0, 2, 1, 1, 7],
   [2, 0, 6, 4, 1],
   [1, 2, 3, 0, 5],
   [6, 5, 4, 3, 2]])

kernel_min=np.array([[0, 1, 0],
   [0, 1, 2],
   [0, 0, 1]])

kernel_max=np.array([[2, 3, 2],
   [2, 2, 4],
   [2, 2, 3]])

min_max = np.dstack((kernel_min,kernel_max))

outcome= [[0, 0, 0, 0, 0],
   [0, 7, 7, 9, 0],
   [0, 8, 8, 8, 0],
   [0, 8, 8, 8, 0],
   [0, 0, 0, 0, 0]]

为此，我创建了一个循环遍历输入数组中所有元素的脚本。在每个元素周围提取内核大小的区域，并对内核范围内的单元格进行计数。

def create_heatmap(input,min_max):
  input_selection=np.zeros(shape(min_max ),dtype='f')
  heatmap = np.zeros(shape(input),dtype='uint16')
  length,width=shape(input)
  center=shape(min_max )[0]/2
  #exclude edge pixels
  for i in range(0+center,length-center):
    for j in range(0+center,width-center):
      # Mask area the size of the kernel around the selected cell
      input_selection= input[i-center:i+center+1,j-center:j+center+1]
      # Count the number of cells within kernel range:
      heatmap[i,j]= shape(np.where((input_selection>=min_max [:,:,0]) & (input_selection<=min_max [:,:,1])))[1]
  return heatmap`

然而，这种遍历所有元素的循环非常耗时（我有一个巨大的数组）。有没有办法加快两个内核范围内的像素计数？例如，计算所有值的移动窗口函数。

Answer 1

您可以通过使用数组步幅技巧在没有 Python 循环的情况下做到这一点。简单的方法是使用现成的功能，例如view_as_windows 来自 Scikit Image 或 extract_patches 来自 Scikit Learn。

from skimage.util import view_as_windows

def with_strides(img, kernel_min, kernel_max):
    win_w, win_h = kernel_min.shape
    heatmap = np.zeros(img.shape, dtype='uint16')
    target = heatmap[win_h//2:-win_h//2+1, win_w//2:-win_w//2+1]

    windowed = view_as_windows(img, (win_h, win_w))
    mask = ((kernel_min <= windowed) & (windowed <= kernel_max))
    np.sum(mask, axis=(2,3), out=target)

    return heatmap

但是，这种方法对内存要求很高。

另一种方法是仍然使用 Python 循环，但循环遍历内核而不是映像。这个想法是通过在更大的数组上调用 Numpy 函数来减少开销。

def kernel_loops(img, kernel_min, kernel_max):
    img_h, img_w = img.shape
    win_h, win_w = kernel_min.shape
    heatmap = np.zeros(img.shape, dtype='uint16')
    target = heatmap[win_h//2:-win_h//2+1, win_w//2:-win_w//2+1]

    for i in range(win_h):
        for j in range(win_w):
            # Negative index fails when slicing to the end
            source = img[i:i+img_h-win_h+1, j:j+img_w-win_w+1]
            target += (kernel_min[i,j] <= source) & (source <= kernel_max[i,j])

    return heatmap

目前我无法提供有代表性的时间安排，但加速应该是显着的。看看窗口大小对性能的影响会很有趣。