How to find local maxima of 3D array in python?

Question

Hi I'm trying to find local maxima in a 3D numpy array, but I can't seem to find a easy way to do that using numpy, scipy, or anything else.

For now I implemented it using scipy.signal.argrelexrema. But it's very long to process large arrays, and only works on separated axis.

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

I know this solution is far from optimal and only works with order=1. Is there any better way to find local maxima in a 3D array in python ?

EDIT :

I use the following method now, which is actually much faster and also works when 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

code optimalisation belongs elsewhere in the SO sister site universe. This is also about statistics. See: code review here or stats here — ZF007
– ZF007, Commented Apr 1, 2019 at 10:44
Look at this code here. stackoverflow.com/questions/49072148/… — Farhood ET
– Farhood ET, Commented Apr 1, 2019 at 10:53
@ZF007 Thank you for the pointer, I asked the question on code review — theobdt
– theobdt, Commented Apr 1, 2019 at 13:59
@FarhoodET Thank you but it's actually different from what I'm looking for — theobdt
– theobdt, Commented Apr 1, 2019 at 14:00
@theobdt, suppose you have two adjacent elements of "data" with the same large value, call it v0, at locations x0,y0,z0 and x0,y0+1,z0. Then wouldn't filtered contain those same values at the two locations, and therefore you would not see "data > filtered" and thus the "fat" maximum would be missed? — Mark Lavin
– Mark Lavin, Commented Sep 1, 2019 at 22:40

NanoBennett · Accepted Answer · 2019-04-01 18:13:40Z

1

Assuming some statistical representation of your data you should be able to perform a 3D local max like this. Hopefully this answers your question.

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))

edited Apr 1, 2019 at 18:13

answered Apr 1, 2019 at 14:00

NanoBennett

1,9001 gold badge13 silver badges13 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

theobdt Over a year ago

Thank you very much for your answer but I think your method is more for detecting global maxima than local maxima, since you are using a global threshold.

Collectives™ on Stack Overflow

How to find local maxima of 3D array in python?

1 Answer 1

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related