1

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?

smcs
  • 1,772
  • 3
  • 18
  • 48
  • 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 Mar 02 '18 at 15:27
  • 2
    3 and 7 are the indices of the local maxima of his example array. – user3483203 Mar 02 '18 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 Mar 02 '18 at 15:38
  • 1
    Not sure, a 1D example will translate well into 3D. But maybe you can get some inspiration [from here](https://stackoverflow.com/questions/9111711/get-coordinates-of-local-maxima-in-2d-array-above-certain-value)? 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 Mar 02 '18 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 Mar 02 '18 at 16:42
  • @speedymcs If you don't set a threshold for your local maxima, you will catch a lot of noise. If you do, you might cut out some of the signal. Threshold is everything for signal processing. – Mr. T Mar 02 '18 at 16:44

1 Answers1

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])
Paul Panzer
  • 51,835
  • 3
  • 54
  • 99
  • Thanks a bunch! Brilliant, this taught me a lot. It's about 15 times faster than my approach. – smcs Mar 07 '18 at 16:41