numpy 矢量化函数以通过严格比较在2D阵列中找到局部最小值和最大值

np8igboo  于 2023-01-09  发布在  其他
关注(0)|答案(1)|浏览(160)

我正在尝试提高一个函数的性能,该函数返回输入2D NumPy数组的局部最小值和最大值。该函数按预期工作,但对于我的用例来说太慢了。我想知道是否可以创建该函数的矢量化版本来提高其性能。
以下是定义元素是否为局部最小值(最大值)的正式定义:

其中

A=[a_m,n]是2D矩阵,mn分别是行和列,w_hw_w分别是滑动窗口的高度和宽度。
我试过使用skimage.morphology.local_minimumskimage.morphology.local_maxima,但是如果一个元素的值小于或等于(大于或等于)它的所有邻居,它们就认为这个元素是最小值(最大值)。
在我的例子中,如果一个元素严格小于(大于)它的所有邻居,我需要函数将它视为最小(最大)。
当前的实现使用numpy.lib.stride_tricks.sliding_window_view的滑动窗口方法,但是函数不一定要使用这种方法。
下面是我当前的实现:

import numpy as np

def get_local_extrema(array, window_size=(3, 3)):
    # Check if the window size is valid
    if not all(size % 2 == 1 and size >= 3 for size in window_size):
        raise ValueError("Window size must be odd and >= 3 in both dimensions.")

    # Create a map to store the local minima and maxima
    minima_map = np.zeros_like(array)
    maxima_map = np.zeros_like(array)

    # Save the shape and dtype of the original array for later
    original_size = array.shape
    original_dtype = array.dtype
    # Get the halved window size
    half_window_size = tuple(size // 2 for size in window_size)

    # Pad the array with NaN values to handle the edge cases
    padded_array = np.pad(array.astype(float),
                         tuple((size, size) for size in half_window_size),
                         mode='constant', constant_values=np.nan)

    # Generate all the sliding windows
    windows = np.lib.stride_tricks.sliding_window_view(padded_array, window_size).reshape(
        original_size[0] * original_size[1], *window_size)

    # Create a mask to ignore the central element of the window
    mask = np.ones(window_size, dtype=bool)
    mask[half_window_size] = False

    # Iterate through all the windows
    for i in range(windows.shape[0]):
        window = windows[i]
        # Get the value of the central element
        center_val = window[half_window_size]
        # Apply the mask to ignore the central element
        masked_window = window[mask]

        # Get the row and column indices of the central element
        row = i // original_size[1]
        col = i % original_size[1]

        # Check if the central element is a local minimum or maximum
        if center_val > np.nanmax(masked_window):
            maxima_map[row, col] = center_val
        elif center_val < np.nanmin(masked_window):
            minima_map[row, col] = center_val

    return minima_map.astype(original_dtype), maxima_map.astype(original_dtype)

a = np.array([[8, 8, 4, 1, 5, 2, 6, 3],
              [6, 3, 2, 3, 7, 3, 9, 3],
              [7, 8, 3, 2, 1, 4, 3, 7],
              [4, 1, 2, 4, 3, 5, 7, 8],
              [6, 4, 2, 1, 2, 5, 3, 4],
              [1, 3, 7, 9, 9, 8, 7, 8],
              [9, 2, 6, 7, 6, 8, 7, 7],
              [8, 2, 1, 9, 7, 9, 1, 1]])

(minima, maxima) = get_local_extrema(a)

print(minima)
# [[0 0 0 1 0 2 0 0]
#  [0 0 0 0 0 0 0 0]
#  [0 0 0 0 1 0 0 0]
#  [0 1 0 0 0 0 0 0]
#  [0 0 0 1 0 0 3 0]
#  [1 0 0 0 0 0 0 0]
#  [0 0 0 0 6 0 0 0]
#  [0 0 1 0 0 0 0 0]]

print(maxima)
# [[0 0 0 0 0 0 0 0]
#  [0 0 0 0 7 0 9 0]
#  [0 8 0 0 0 0 0 0]
#  [0 0 0 4 0 0 0 8]
#  [6 0 0 0 0 0 0 0]
#  [0 0 0 0 0 0 0 8]
#  [9 0 0 0 0 0 0 0]
#  [0 0 0 9 0 9 0 0]]

expected_minima = np.array([[0, 0, 0, 1, 0, 2, 0, 0],
                            [0, 0, 0, 0, 0, 0, 0, 0],
                            [0, 0, 0, 0, 1, 0, 0, 0],
                            [0, 1, 0, 0, 0, 0, 0, 0],
                            [0, 0, 0, 1, 0, 0, 3, 0],
                            [1, 0, 0, 0, 0, 0, 0, 0],
                            [0, 0, 0, 0, 6, 0, 0, 0],
                            [0, 0, 1, 0, 0, 0, 0, 0]])

expected_maxima = np.array([[0, 0, 0, 0, 0, 0, 0, 0],
                            [0, 0, 0, 0, 7, 0, 9, 0],
                            [0, 8, 0, 0, 0, 0, 0, 0],
                            [0, 0, 0, 4, 0, 0, 0, 8],
                            [6, 0, 0, 0, 0, 0, 0, 0],
                            [0, 0, 0, 0, 0, 0, 0, 8],
                            [9, 0, 0, 0, 0, 0, 0, 0],
                            [0, 0, 0, 9, 0, 9, 0, 0]])

np.testing.assert_array_equal(minima, expected_minima)
np.testing.assert_array_equal(maxima, expected_maxima)

print('All tests passed')

关于如何向量化这个函数的任何建议或想法将非常感谢。
先谢了!

    • 编辑#1**

在使用NumPy玩了一会儿之后,我设法让下面的代码几乎可以工作了,如果我理解正确的话,这是一种完全矢量化的方式:

def get_local_extrema_2(img):
  minima_map = np.zeros_like(img)
  maxima_map = np.zeros_like(img)

  minima_map[1:-1, 1:-1] = np.where(
    (a[1:-1, 1:-1] < a[:-2, 1:-1]) &
    (a[1:-1, 1:-1] < a[2:, 1:-1]) &
    (a[1:-1, 1:-1] < a[1:-1, :-2]) &
    (a[1:-1, 1:-1] < a[1:-1, 2:]) &
    (a[1:-1, 1:-1] < a[2:, 2:]) &
    (a[1:-1, 1:-1] < a[:-2, :-2]) &
    (a[1:-1, 1:-1] < a[2:, :-2]) &
    (a[1:-1, 1:-1] < a[:-2, 2:]),
    a[1:-1, 1:-1],
    0)
  
  maxima_map[1:-1, 1:-1] = np.where(
    (a[1:-1, 1:-1] > a[:-2, 1:-1]) &
    (a[1:-1, 1:-1] > a[2:, 1:-1]) &
    (a[1:-1, 1:-1] > a[1:-1, :-2]) &
    (a[1:-1, 1:-1] > a[1:-1, 2:]) &
    (a[1:-1, 1:-1] > a[2:, 2:]) &
    (a[1:-1, 1:-1] > a[:-2, :-2]) &
    (a[1:-1, 1:-1] > a[2:, :-2]) &
    (a[1:-1, 1:-1] > a[:-2, 2:]),
    a[1:-1, 1:-1],
    0)

  return minima_map, maxima_map

get_local_extrema_2的输出为:
最小值Map:

[[0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0]
 [0 0 0 0 1 0 0 0]
 [0 1 0 0 0 0 0 0]
 [0 0 0 1 0 0 3 0]
 [0 0 0 0 0 0 0 0]
 [0 0 0 0 6 0 0 0]
 [0 0 0 0 0 0 0 0]]

最大值图:

[[0 0 0 0 0 0 0 0]
 [0 0 0 0 7 0 9 0]
 [0 8 0 0 0 0 0 0]
 [0 0 0 4 0 0 0 0]
 [0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0]]

上述方法的问题在于没有检测到边界上的最小值或最大值像素。

    • 编辑#2**

即使在输出阵列中存在1而不是局部最小值(最大值)的值,即0和1(或假和真)的2d阵列,也是好的。

wqlqzqxt

wqlqzqxt1#

您对局部最大值的定义有缺陷。例如,在一维数组[1,2,3,4,4,3,2,1]中,存在一个局部最大值,但您的定义忽略了它。skimage.morphology.local_maxima将正确标识此局部最大值。
如果你真的需要实现你的定义,我会使用一个窗口大小的正方形结构元素的膨胀(腐 eclipse ),但不包括中心像素。原始图像中大于(小于)过滤图像中的任何像素都将满足你的局部最大值(最小值)的定义。
我使用scikit-image实现了这个功能,但发现它在图像边缘会做一件奇怪的事情,所以它不会检测边缘附近的局部极大值或极小值:

se = np.ones((3, 3))
se[1, 1] = 0
minima = a * (skimage.morphology.erosion(a, footprint=se) > a)
maxima = a * (skimage.morphology.dilation(a, footprint=se) < a)

使用DIPlib(公开:我是一个作者)这将正确工作,也在图像边缘:

import diplib as dip

se = np.ones((3, 3), dtype=np.bool_)
se[1, 1] = False
minima = a * (dip.Erosion(a, se) > a)
maxima = a * (dip.Dilation(a, se) < a)

查看skimage.morphology.dilation的源代码,它使用默认的边界扩展'reflect'调用scipy.ndimage.grey_dilation,这意味着图像边缘的每个局部最大值都有一个具有相同值的邻居,因此在此定义中不会被检测为局部最大值,而是应该使用'constant'扩展。cval设置为数据类型的最小可能值。例如,对于uint8输入数组,它应该执行ndi.grey_dilation(image, footprint=footprint, output=out, mode='constant', cval=0). GitHub issue

相关问题