Find maxima of segmented regions

user_10352072 0
Viewed 68

I have 2 3D grayscale images of equal dimensions.

The first image (img) is the raw data of a 3D volume (obtained using a microscope) and contains various cells. Example slice of raw 3d image

The second image (imglab) is a masked version of the raw image where each identified cell is filled with a single, unique value (i.e. cell 1 = all 1s, cell 2 = all 2s). All non cell regions are zeros. Example slice of masked 3d image

I am now trying to find the coordinates of the maximum value of each cell from the raw data that corresponds to the labeled mask array.

Currently, I have a loop which is extremely inefficient. I suspect there is a way to setup multiple conditions using a single np.where call, but I can't figure out how to do this. The current for loop method is below:

coordinates = []
for i in range(1, int(imglab.max())): # Number of cells = max value of label image
    max_val = np.max(img[imglab == i])
    max_coord = np.where((img == max_val) & (imglab == i))
    coordinates.append(max_coord)
image-processing numpy python
1 Answers

When it is difficult to find a way of coding efficient and compatible with numpy, but when the code with for loops is trivial, you can use njitfrom numba.

It works best to process flat arrays, so first let's write a function in numba that does what you ask but in 1d:

from numba import njit, int64

@njit
def fast_max_flat(img_flat, imglab_flat):
    n_cells =int(imglab_flat.max())  # number of cells
    max_values = np.full(n_cells, - np.inf)  # stores the n_cells max values seen so far
    max_coords = np.zeros(n_cells, dtype=int64)  # stores the corresponding coordinate
    n_pixels = len(img)
    for i in range(n_pixels):
        label = imglab_flat[i]
        value = img_flat[i]
        if max_values[label] < value:
            max_values[label] = value
            max_coords[label] = i
    return max_coords

And then write a python wrapper that ravels the array, applies the previous function, and retrieves the coordinates as a list:

def wrapper(img, imglab):
    dim = img.shape
    coords = fast_max_flat(img.ravel(), imglab.ravel())
    return [np.unravel_index(coord, dim) for coord in coords]

On my machine, with a 100 x 100 x 100 image of 3 cells, this is ~50 times faster than your method.

user_72634050
Related