Найти максимумы сегментированных областей

Question

У меня есть 2 3D-изображения в градациях серого одинакового размера.

Первое изображение (img) представляет собой необработанные данные трехмерного тома (полученные с помощью микроскопа) и содержит различные ячейки. Пример фрагмента необработанного трехмерного изображения

Второе изображение (imglab) является замаскированной версией необработанного изображения, где каждая идентифицированная ячейка заполнена одним уникальным значением (то есть ячейка 1 = все 1 с, ячейка 2 = все 2 с). Все не ячейки являются нулями. Пример фрагмента маскированного трехмерного изображения

Сейчас я пытаюсь найти координаты максимального значения каждой ячейки из необработанных данных, соответствующих помеченному массиву масок.

В настоящее время у меня есть цикл, который крайне неэффективен. Я подозреваю, что есть способ настроить несколько условий, используя один вызов np.where, но я не могу понять, как это сделать. Текущий метод петли for приведен ниже:

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)

Answer 1

Когда трудно найти способ кодирования, эффективный и совместимый с numpy, но когда код с петлями for тривиален, вы можете использовать njit из numba.

Лучше всего обрабатывать плоские массивы, поэтому сначала напишем функцию в numba, которая выполняет то, что вы просите, но в 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

А затем напишите оболочку Python, которая разбивает массив, применяет предыдущую функцию и получает координаты в виде списка:

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]

На моей машине с 100 x 100 x 100 изображением из 3 ячеек это примерно в 50 раз быстрее, чем ваш метод.