Fastest way to get bounding boxes around segments in a label map

algorithm, numpy, performance, python

Solution

If I understood your problem correctly, you have groups of voxels, and you would like to have the extremes of a group in each axis.

Let'd define:

`arr`: 3D array of integer labels

`labels`: list of labels (integers 0..labmax)

The code:

import numpy as np

# number of highest label:
labmax = np.max(labels)

# maximum and minimum positions along each axis (initialized to very low and high values)
b_first = np.iinfo('int32').min * np.ones((3, labmax + 1), dtype='int32')
b_last = np.iinfo('int32').max * np.ones((3, labmax + 1), dtype='int32')

# run through all of the dimensions making 2D slices and marking all existing labels to b
for dim in range(3):
    # create a generic slice object to make the slices
    sl = [slice(None), slice(None), slice(None)]

    bf = b_first[dim]
    bl = b_last[dim]

    # go through all slices in this dimension
    for k in range(arr.shape[dim]):
        # create the slice object
        sl[dim] = k
        # update the last "seen" vector
        bl[arr[sl].flatten()] = k

        # if we have smaller values in "last" than in "first", update
        bf[:] = np.clip(bf, None, bl)

After this operation we have six vectors giving the smallest and largest indices for each axis. For example, the bounding values along second axis of label `13` are `b_first[1][13]` and `b_last[1][13]`. If some label is missing, all corresponding `b_first` and `b_last` will be the maximum `int32` value.

I tried this with my computer, and for a (300,300,200) array it takes approximately 1 sec to find the values.

Problem

A 3D label map is matrix in which every pixel (voxel) has an integer label. These values are expected to be contiguous, meaning that a segment with label `k` will not be fragmented. Given such label map (segmentation), what is the fastest way to obtain the coordinates of a minimum bounding box around each segment, in Python? I have tried the following: - Iterate through the matrix using multiindex iterator (from `numpy.nditer`) and construct a reverse index dictionary. This means that for every label you get the 3 coordinates of every voxel where the label is present. - For every label get the max and min of each coordinate. The good thing is that you get all the location information in one O(N) pass. The bad thing is that I dont need this detailed information. I just need the extremities, so there might be a faster way to do this, using some numpy functions which are faster than so many list appends. Any suggestions? The one pass through the matrix takes about 8 seconds on my machine, so it would be great to get rid of it. To give an idea of the data, there are a few hundred labels in a label map. Sizes of the label map can be 700x300x30 or 300x300x200 or something similar. Edit: Now storing only updated max and min per coordinate for every label. This removes the need to maintain and store all these large lists (append).

Original source