Как найти локальные максимумы 3D-массива в python?
Привет, я пытаюсь найти локальные максимумы в трехмерном массиве 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
Я знаю, что это решение далеко от оптимального и работает только с порядком = 1. Есть ли лучший способ найти локальные максимумы в трехмерном массиве в python?
РЕДАКТИРОВАТЬ :
Сейчас я использую следующий метод, который на самом деле намного быстрее, а также работает, когда порядок > 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
Посмотрите на этот код здесь. stackoverflow.com/questions/49072148/…
@ ZF007 ZF007 Спасибо за указание, я задал вопрос о проверке кода.
@FarhoodET Спасибо, но на самом деле это отличается от того, что я ищу
@theobdt, предположим, у вас есть два соседних элемента «данных» с одинаковым большим значением, назовите его v0, в местах x0, y0, z0 и x0, y0 + 1, z0. Тогда не будет ли filtered содержать одни и те же значения в двух местах, и, следовательно, вы не увидите «данные> отфильтровано», и, таким образом, будет пропущен «жирный» максимум?
Ответы 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))
Большое спасибо за ваш ответ, но я думаю, что ваш метод больше предназначен для обнаружения глобальных максимумов, чем локальных максимумов, поскольку вы используете глобальный порог.
оптимизация кода относится к другим местам во вселенной родственных сайтов SO. Это тоже про статистику. См.: обзор кода здесь или статистика здесь.