如何在Python中找到3D数组的局部最大值？

Question

你好，我正在尝试在一个3D的numpy数组中找到局部最大值，但是使用numpy、scipy或其他方式似乎都没有简单的方法。

目前我使用了scipy.signal.argrelexrema进行了实现。但是处理大型数组非常耗时，而且只能在分离的轴上运行。

import numpy as np
from scipy.signal import argrelextrema


def local_maxima_3D(data, order=1):
    """Detects local maxima in a 3D array

    Parameters
    ---------
    data : 3d ndarray
    order : int
        How many points on each side to use for the comparison

    Returns
    -------
    coords : ndarray
        coordinates of the local maxima
    values : ndarray
        values of the local maxima
    """
    # Coordinates of local maxima along each axis
    peaks0 = np.array(argrelextrema(data, np.greater, axis=0, order=order))
    peaks1 = np.array(argrelextrema(data, np.greater, axis=1, order=order))
    peaks2 = np.array(argrelextrema(data, np.greater, axis=2, order=order))

    # Stack all coordinates 
    stacked = np.vstack((peaks0.transpose(), peaks1.transpose(),
                         peaks2.transpose()))

    # We keep coordinates that appear three times (once for each axis)
    elements, counts = np.unique(stacked, axis=0, return_counts=True)
    coords = elements[np.where(counts == 3)[0]]

    # Compute values at filtered coordinates
    values = data[coords[:, 0], coords[:, 1], coords[:, 2]]

    return coords, values

我知道这种解决方案远非最佳，并且仅在order=1时起作用。有没有更好的方法在Python中查找3D数组中的局部极大值？ 编辑: 我现在使用以下方法，它实际上比较快，而且当order> 1时也有效：

import numpy as np
from scipy import ndimage as ndi


def local_maxima_3D(data, order=1):
    """Detects local maxima in a 3D array

    Parameters
    ---------
    data : 3d ndarray
    order : int
        How many points on each side to use for the comparison

    Returns
    -------
    coords : ndarray
        coordinates of the local maxima
    values : ndarray
        values of the local maxima
    """
    size = 1 + 2 * order
    footprint = np.ones((size, size, size))
    footprint[order, order, order] = 0

    filtered = ndi.maximum_filter(data, footprint=footprint)
    mask_local_maxima = data > filtered
    coords = np.asarray(np.where(mask_local_maxima)).T
    values = data[mask_local_maxima]

    return coords, values

- theobdt

---

代码优化应该在SO姊妹站点宇宙的其他地方进行。这也涉及统计数据。请参见：此处的代码审查或此处的统计数据。 - ZF007

看看这里的代码。 https://stackoverflow.com/questions/49072148/finding-local-maxima-in-large-3d-numpy-arrays - Farhood ET

@ZF007 谢谢你的指引，我在代码审查上提出了问题。 - theobdt

@FarhoodET 谢谢你，但实际上与我所寻找的不同。 - theobdt

@theobdt，假设您有两个相邻的“数据”元素具有相同的大值v0，在位置x0，y0，z0和x0，y0 + 1，z0。那么过滤后的结果不会包含这两个位置上的相同值，因此您将看不到“data> filtered”，从而错过了“fat”最大值？ - Mark Lavin

Answer 1

假设您的数据有一些统计表示，您应该能够执行类似于这样的3D本地最大值。希望这回答了您的问题。

import numpy as np
import scipy.ndimage as ndimage

img = np.random.normal(size=(100, 256, 256))

# Get local maximum values of desired neighborhood
# I'll be looking in a 5x5x5 area
img2 = ndimage.maximum_filter(img, size=(5, 5, 5))

# Threshold the image to find locations of interest
# I'm assuming 6 standard deviations above the mean for the threshold
img_thresh = img2.mean() + img2.std() * 6

# Since we're looking for maxima find areas greater than img_thresh

labels, num_labels = ndimage.label(img2 > img_thresh)

# Get the positions of the maxima
coords = ndimage.measurements.center_of_mass(img, labels=labels, index=np.arange(1, num_labels + 1))

# Get the maximum value in the labels
values = ndimage.measurements.maximum(img, labels=labels, index=np.arange(1, num_labels + 1))

- NanoBennett

---

非常感谢您的回答，但我认为您的方法更适用于检测全局最大值而不是局部最大值，因为您使用了全局阈值。 - theobdt