Finding local maxima in large 3D Numpy arrays

Question

I'm processing some large volumetric image data that are present in three dimensional numpy arrays. I'll explain my task with two small 1D arrays. I have one image:

img = [5, 6, 70, 80, 3, 4, 80, 90]

and one segmented and labeled version of that image:

labels = [0, 0, 1, 1, 0, 0, 2, 2]

Each number in labels represents an object in img. Both arrays have the same dimensions. So in this example there's two objects in img:

[5, 6, 70, 80, 3, 4, 80, 90]

and what I'm trying to do now is finding the location of the maximum value of each object, which in this case would be the 3 and 7. Currently I loop over all labels, create a version of img which contains only the object corresponding to the current label, and look for the maximum value:

for label in range(1, num_labels + 1):
    imgcp = np.copy(img)
    imgcp[labels != label] = 0
    max_pos = np.argmax(imgcp)
    max_coords = np.unravel_index(pos, imgcp.shape)

One problem with this approach is that copying img in every step tends to create memory errors. I feel like memory management should prevent this, but is there a more memory efficient and possibly faster way to to this task?

Didn't understand the link between labels and img, nor what the 3 and 7 were. Try to look into np.argmax() might help you. — Mathieu
– Mathieu, Commented Mar 2, 2018 at 15:27
3 and 7 are the indices of the local maxima of his example array. — user3483203
– user3483203, Commented Mar 2, 2018 at 15:33
@Mathieu, basically each region of the same value in labels corresponds to a region of an object in img. I'm already using np.argmax() in my example; the difficulty is to apply it locally, i.e. once per label/object — smcs
– smcs, Commented Mar 2, 2018 at 15:38
Not sure, a 1D example will translate well into 3D. But maybe you can get some inspiration from here? Imho you have to be more specific, what you want to count as a local maximum in terms of area size and value difference. — Mr. T
– Mr. T, Commented Mar 2, 2018 at 16:34
@MrT, I am looking for the absolute maximum of each single object in img, as delimited by the labels array. Your link looks promising, I'll be able to try it next week! — smcs
– smcs, Commented Mar 2, 2018 at 16:42

Paul Panzer · Accepted Answer · 2018-03-02 17:10:03Z

2

Here is a method using argpartition.

# small 2d example
>>> data = np.array([[0,1,4,0,0,2,1,0],[0,4,1,3,0,0,0,0]])
>>> segments = np.array([[0,1,1,0,0,2,2,0],[0,1,1,1,0,0,0,0]])
>>> 
# discard zeros
>>> nz = np.where(segments)
>>> segc = segments[nz]
>>> dac = data[nz]

# count object sizes
>>> cnts = np.bincount(segc)
>>> bnds = np.cumsum(cnts)
# use counts to partition into objects
>>> idx = segc.argpartition(bnds[1:-1])
>>> dai = dac[idx]
# find maxima per object
>>> mx = np.maximum.reduceat(dai, bnds[:-1])
# find their positions
>>> am, = np.where(dai==mx.repeat(cnts[1:]))
# translate positions back to coordinate space
>>> im = idx[am]
>>> am = *(n[im] for n in nz),
>>> 
>>> 
# result
# coordinates, note that there are more points than objects because
# the maximum 4 occurs twice in object 1
>>> am
(array([1, 0, 0]), array([1, 2, 5]))
# maxima
>>> data[am]
array([4, 4, 2])
# labels
>>> segments[am]
array([1, 1, 2])

answered Mar 2, 2018 at 17:10

Paul Panzer

53.3k3 gold badges60 silver badges103 bronze badges

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

1 Comment

smcs Over a year ago

Thanks a bunch! Brilliant, this taught me a lot. It's about 15 times faster than my approach.

Tony · Accepted Answer · 2025-09-05 20:19:09Z

Did you ever check out maximum_position from scipy.ndimage?

The examples straight from their docs

>>> from scipy import ndimage
>>> import numpy as np
>>> a = np.array([[1, 2, 0, 0],
...               [5, 3, 0, 4],
...               [0, 0, 0, 7],
...               [9, 3, 0, 0]])
>>> ndimage.maximum_position(a)
(3, 0)

Features to process can be specified using `labels` and `index`:

>>> lbl = np.array([[0, 1, 2, 3],
...                 [0, 1, 2, 3],
...                 [0, 1, 2, 3],
...                 [0, 1, 2, 3]])
>>> ndimage.maximum_position(a, lbl, 1)
(1, 1)

If no index is given, non-zero `labels` are processed:

>>> ndimage.maximum_position(a, lbl)
(2, 3)

If there are no maxima, the position of the first element is returned:

>>> ndimage.maximum_position(a, lbl, 2)
(0, 2)

Collectives™ on Stack Overflow

Finding local maxima in large 3D Numpy arrays

2 Answers 2

1 Comment

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related