Как ускорить расширение 3D-области в массиве boolean numpy? - PullRequest
Новая
       36

Как ускорить расширение 3D-области в массиве boolean numpy?

0 голосов
/ 16 июня 2020

У меня есть логическая маска массива 3D numpy, которая была сегментирована из тома мозга МРТ.
Воксели мозга = True. Все остальное = Ложь.

Я бы хотел увеличить эту маску так, чтобы она охватывала окружающие ткани в объеме МРТ, а не только сегментированный орган, возможно, 10-миллиметровую корку, не относящуюся к мозгу. вокруг мозга.

Я пробовал использовать двухмерную дилатацию, используя skimage.morphology.dilation с алмазным фильтром. Хотя это хорошо и быстро для одного изображения, мне нужно повторить это на нескольких срезах объема и как минимум в двух плоскостях, чтобы приблизиться к равномерному расширению 3D-маски.

Я в основном взял свой код отсюда: https://scipy-lectures.org/packages/scikit-image/index.html

типичная форма объема = 512, 512, 270 

# 1st pass in axial plane
(x, y, z) = np.shape(3dMask)
for slice_number in range(z):
    image_slice = 3dMask[:, :, slice_number]
    3dMask[:, :, slice_number] = morphology.binary_dilation(image_slice, morphology.diamond(30))
# repeat in coronal plane...

Это очень хорошо работает с желаемым эффектом в каждом фрагменте, но очень медленно для 3D. Я могу ускорить процесс, расширив только те фрагменты, которые содержат хотя бы один «Истинный», но это неизбежно оставляет 100+ фрагментов в каждой плоскости. Все еще медленный.

В надежде, что побочный цикл python все замедляет, я искал трехмерную эквивалентную одиночную функцию в numpy и skimage, но не нашел ничего, что я мог бы признать полезным.

Я поигрался с идеей найти геометрический c центр и просто увеличить объем на 5%, но в маске обязательно будут дыры (пространство между двумя половинами мозга) который больше не будет соответствовать объему МРТ и поэтому бесполезен ...

Полагаю, это означает, что я делаю это неправильно, поскольку я новичок как в numpy, так и в skimage.

Есть ли быстрый способ сделать это? Возможно, это трехмерная альтернатива расширению 2D-изображения?

1 Ответ

1 голос
/ 17 июня 2020

Этот вопрос на самом деле имеет небольшую тонкость, которую я постараюсь раскрыть.

Первое, что следует отметить, это то, что большинство функций scikit-image действительно отлично работают в 3D, включая binary_dilation! Итак, вы должны в идеальном мире уметь: 

dilated = morphology.binary_dilation(
    mask3d, morphology.ball(radius=30)
)

Я говорю в идеальном мире, потому что это дает сбой на моей машине, вероятно, потому что эта давняя ошибка SciPy предотвращает работу фильтров SciPy (которые scikit-image использует под капотом) с большими размерами окрестностей. с радиусом 30 на самом деле то же самое, что и 30-кратное расширение ромба с радиусом 1! Вы можете сделать это вручную в for-l oop или можете использовать scipy.ndimage.binary_dilation, используя аргумент ключевого слова iterations. (См. этот выпуск для обсуждения этого вопроса.) 

from scipy import ndimage as ndi

# make a little 3D diamond:
diamond = ndi.generate_binary_structure(rank=3, connectivity=1)
# dilate 30x with it
dilated = ndi.binary_dilation(mask3d, diamond, iterations=30)

На самом деле вы можете довольно далеко продвинуться с этой стратегией. Например, если ваш набор данных не имеет одинакового разрешения по осям x, y и z, возможно, вы захотите расширить, скажем, вдвое больше по x и y. Вы можете сделать это в два этапа: 

dilated1 = ndi.binary_dilation(mask3d, diamond, iterations=15)
flat = np.copy(diamond)
flat[:, :, 0] = 0
flat[:, :, -1] = 0
dilated2 = ndi.binary_dilation(mask3d, flat, iterations=15)

Наконец, обратите внимание, что двоичное расширение эквивалентно (недвоичной) свертке с последующим пороговым значением выше 0. Итак, я обнаружил, что это также работает: 

from scipy import signal

b = morphology.ball(radius=30)
dilated = signal.fftconvolve(mask3d, b, mode='same') > 0

Однако для этого размера изображения и на моей машине это было медленнее, чем повторное расширение. Но об этом стоит помнить, потому что производительность будет разной для разных наборов данных.

В качестве примечания я рекомендую публиковать полный рабочий код в ваших вопросах StackOverflow, как описано здесь . В вашем случае np.shape(3dMask) является синтаксической ошибкой, поскольку 3dMask не является допустимым идентификатором Python! =)

Надеюсь, это поможет!

...