如何在NumPy中找到平滑多维数组的局部最小值

Question

假设我有一个包含连续可微函数评估结果的NumPy数组，我想要找到局部最小值。由于没有噪音，所以任何值低于其所有邻居的点都符合我对局部最小值的标准。

我有以下的列表推导式，适用于二维数组，忽略边界上的潜在最小值：

import numpy as N

def local_minima(array2d):
    local_minima = [ index 
                     for index in N.ndindex(array2d.shape)
                     if index[0] > 0
                     if index[1] > 0
                     if index[0] < array2d.shape[0] - 1
                     if index[1] < array2d.shape[1] - 1
                     if array2d[index] < array2d[index[0] - 1, index[1] - 1]
                     if array2d[index] < array2d[index[0] - 1, index[1]]
                     if array2d[index] < array2d[index[0] - 1, index[1] + 1]
                     if array2d[index] < array2d[index[0], index[1] - 1]
                     if array2d[index] < array2d[index[0], index[1] + 1]
                     if array2d[index] < array2d[index[0] + 1, index[1] - 1]
                     if array2d[index] < array2d[index[0] + 1, index[1]]
                     if array2d[index] < array2d[index[0] + 1, index[1] + 1]
                   ]
    return local_minima

然而，这个速度相当慢。我还希望能够使其适用于任意维度的情况。例如，是否有一种简单的方法可以获取数组中任意维度的点的所有邻居？或者我完全错误地解决了这个问题？我应该使用`numpy.gradient()`吗？

- ptomato

---

寻找全局最大值：https://dev59.com/4nA65IYBdhLWcg3w4CwL#3584260 - endolith

Answer 1

可以使用Ivan的detect_peaks函数（稍作修改）来查找任意维度数组的局部最小值的位置：

import numpy as np
import scipy.ndimage.filters as filters
import scipy.ndimage.morphology as morphology

def detect_local_minima(arr):
    # https://dev59.com/XnA65IYBdhLWcg3wsgyA#3689710
    """
    Takes an array and detects the troughs using the local maximum filter.
    Returns a boolean mask of the troughs (i.e. 1 when
    the pixel's value is the neighborhood maximum, 0 otherwise)
    """
    # define an connected neighborhood
    # http://www.scipy.org/doc/api_docs/SciPy.ndimage.morphology.html#generate_binary_structure
    neighborhood = morphology.generate_binary_structure(len(arr.shape),2)
    # apply the local minimum filter; all locations of minimum value 
    # in their neighborhood are set to 1
    # http://www.scipy.org/doc/api_docs/SciPy.ndimage.filters.html#minimum_filter
    local_min = (filters.minimum_filter(arr, footprint=neighborhood)==arr)
    # local_min is a mask that contains the peaks we are 
    # looking for, but also the background.
    # In order to isolate the peaks we must remove the background from the mask.
    # 
    # we create the mask of the background
    background = (arr==0)
    # 
    # a little technicality: we must erode the background in order to 
    # successfully subtract it from local_min, otherwise a line will 
    # appear along the background border (artifact of the local minimum filter)
    # http://www.scipy.org/doc/api_docs/SciPy.ndimage.morphology.html#binary_erosion
    eroded_background = morphology.binary_erosion(
        background, structure=neighborhood, border_value=1)
    # 
    # we obtain the final mask, containing only peaks, 
    # by removing the background from the local_min mask
    detected_minima = local_min ^ eroded_background
    return np.where(detected_minima)       

你可以像这样使用：

arr=np.array([[[0,0,0,-1],[0,0,0,0],[0,0,0,0],[0,0,0,0],[-1,0,0,0]],
              [[0,0,0,0],[0,-1,0,0],[0,0,0,0],[0,0,0,-1],[0,0,0,0]]])
local_minima_locations = detect_local_minima(arr)
print(arr)
# [[[ 0  0  0 -1]
#   [ 0  0  0  0]
#   [ 0  0  0  0]
#   [ 0  0  0  0]
#   [-1  0  0  0]]

#  [[ 0  0  0  0]
#   [ 0 -1  0  0]
#   [ 0  0  0  0]
#   [ 0  0  0 -1]
#   [ 0  0  0  0]]]

这意味着最小值出现在索引[0,0,3]，[0,4,0]，[1,1,1]和[1,3,3]：

print(local_minima_locations)
# (array([0, 0, 1, 1]), array([0, 4, 1, 3]), array([3, 0, 1, 3]))
print(arr[local_minima_locations])
# [-1 -1 -1 -1]

- unutbu

---

不错！它的运行速度大约比我的原始版本快65倍，并且适用于任何维数。 - ptomato

对于较新版本的numpy，它会出现以下错误：numpy boolean subtract, the - operator, is deprecated, use the bitwise_xor, the ^ operator, or the logical_xor function instead.。也许可以将detected_minima = local_min - eroded_background更改为bitwise_xor、^或logical_xor。 - Echan

Answer 2

尝试以下方法实现2D效果：

import numpy as N

def local_minima(array2d):
    return ((array2d <= N.roll(array2d,  1, 0)) &
            (array2d <= N.roll(array2d, -1, 0)) &
            (array2d <= N.roll(array2d,  1, 1)) &
            (array2d <= N.roll(array2d, -1, 1)))

这将返回一个类似于2D数组的数组，其中包含局部最小值（四个相邻值）的True/False值。

- eumiro

---

这实际上是找到本地最大值，并需要使用&而不是&&，并需要在比较周围加上括号，但它的运行速度比我的原始代码快30倍。 - ptomato

这里发生了什么？ - john k

@johnktejik 它在将数组中的元素与其上方、下方、左侧和右侧的元素进行比较。如果该元素小于所有元素，则用“True”替换该元素，以告知您此处存在局部最小值。 - Joseph Farah

这取决于您对局部最小值的定义。如果您处于一个“山谷”中，但是该“山谷”在最小值处有一个高原，则整个高原将被忽略。该山谷的任何局部最小值都不会被记录。 - Thomas Hubregtsen