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

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