How to find the local minima of a smooth multidimensional array in NumPy efficiently?

Question

Say I have an array in NumPy containing evaluations of a continuous differentiable function, and I want to find the local minima. There is no noise, so every point whose value is lower than the values of all its neighbors meets my criterion for a local minimum.

I have the following list comprehension which works for a two-dimensional array, ignoring potential minima on the boundaries:

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

However, this is quite slow. I would also like to get this to work for any number of dimensions. For example, is there an easy way to get all the neighbors of a point in an array of any dimensions? Or am I approaching this problem the wrong way altogether? Should I be using numpy.gradient() instead?

Finding the global maxima: http://stackoverflow.com/questions/3584243/python-get-the-position-of-the-biggest-item-in-an-numpy-array/3584260#3584260 — endolith, Oct 06 '11 at 03:45

unutbu · Accepted Answer · 2019-11-06T16:08:03.780

The location of the local minima can be found for an array of arbitrary dimension using Ivan's detect_peaks function, with minor modifications:

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

def detect_local_minima(arr):
    # https://stackoverflow.com/questions/3684484/peak-detection-in-a-2d-array/3689710#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)

which you can use like this:

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

This says the minima occur at indices [0,0,3], [0,4,0], [1,1,1] and [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]

Nice! It runs about 65 times as fast as my original, and works for any number of dimensions. — ptomato, Oct 21 '10 at 14:22
For newer version of numpy. It will have this error: `numpy boolean subtract, the - operator, is deprecated, use the bitwise_xor, the ^ operator, or the logical_xor function instead.`. Maybe change `detected_minima = local_min - eroded_background` to `bitwise_xor`, `^` or `logical_xor ` — Echan, Nov 06 '19 at 14:49

eumiro · Answer 2 · 2010-10-21T14:18:06.207

5

Try this for 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)))

This will return you an array2d-like array with True/False where local minima (four neighbors) are located.

edited Oct 21 '10 at 14:18

answered Oct 21 '10 at 10:29

eumiro

207,213
34
299
261

Well, this actually finds local maxima, and needs `&`'s instead of `&&`'s, and requires parentheses around the comparisons, but it runs thirty times faster than my original. – ptomato Oct 21 '10 at 14:17
what is happening in here? – john k Mar 18 '18 at 03:52
@johnktejik its comparing an element in the array to the elements above, below, to the left, and to the right. if its smaller than all, it replaces the element with "True" to let you know there is a local minimum there. – Joseph Farah Jul 09 '18 at 17:19
It depends on your definition of local minima. If you are in a "valley", but the "valley" has a plateau at the minima, the entire plateau will be ignored. No local minima for that valley will be recorded – Thomas Hubregtsen May 03 '20 at 10:00

How to find the local minima of a smooth multidimensional array in NumPy efficiently?

2 Answers2

Linked