Как выполнить объемную визуализацию 3D-массива numpy с помощью VTK (MarchingCubes) в Python?

#python #image #vtk #pydicom #medical-imaging

Вопрос:

У меня есть массив 3D numpy, и я пытаюсь выполнить его объемную визуализацию с помощью VTK. Однако при визуализации я получаю совершенно другой объемный рендеринг. Я подозреваю, что это как-то связано с моим преобразованием массива numpy в формат изображения VTK, но я не могу понять, где я ошибаюсь. Я загрузил массив numpy сюда.

Может кто-нибудь помочь мне понять, где я ошибаюсь?

Это мой код:

 #!/usr/bin/env python

import os
import numpy as np

ArrayDicom = np.load('test3.npy')
data_matrix = ArrayDicom
w, d, h = ArrayDicom.shape

colors = vtkNamedColors()
iso_value = 200

reader = vtkImageImport()
data_string = data_matrix.tobytes()
reader.CopyImportVoidPointer(data_string, len(data_string))
reader.SetDataScalarTypeToUnsignedChar()
reader.SetNumberOfScalarComponents(1)
reader.SetDataExtent(0, w-1, 0, d-1, 0, h-1)
reader.SetWholeExtent(0, w-1, 0, d-1, 0, h-1)
reader.Update()

volume = vtkImageData()
volume.DeepCopy(reader.GetOutput())

surface = vtkMarchingCubes()
surface.SetInputData(volume)
surface.ComputeNormalsOn()
surface.SetValue(0, iso_value)

renderer = vtkRenderer()
renderer.SetBackground(colors.GetColor3d('DarkSlateGray'))

render_window = vtkRenderWindow()
render_window.AddRenderer(renderer)
render_window.SetWindowName('MarchingCubes')

interactor = vtkRenderWindowInteractor()
interactor.SetRenderWindow(render_window)

mapper = vtkPolyDataMapper()
mapper.SetInputConnection(surface.GetOutputPort())
mapper.ScalarVisibilityOff()

actor = vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetColor(colors.GetColor3d('MistyRose'))

renderer.AddActor(actor)

render_window.Render()
interactor.Start()
 

Это мой объемный рендеринг:

введите описание изображения здесь

Это мой ожидаемый объемный рендеринг:

введите описание изображения здесь

Ответ №1:

Numpy использует другой порядок массивов, чем VTK. Вы должны быть в состоянии изменить порядок w, h и d, чтобы получить нужную вещь.

Вот как ты этого хочешь:

 h, d, w = ArrayDicom.shape
 

Хорошо, вот сценарий преобразования, который я использовал для преобразования в файл VTK:

 import numpy as np
import SimpleITK as sitk

x = np.load("test3.npy")
y = sitk.GetImageFromArray(x)
sitk.WriteImage(y, "test3.vtk")
 

Это не так приятно, как правильно настроить импорт изображений VTK, но, ну, я парень из SimpleITK, и я знаю, что преобразование numpy работает в SimpleITK.

Комментарии:

1. Привет, Дэйв, спасибо за твой ответ. Я попытался изменить h, d и w, но, похоже, это не работает, я все равно получаю объем той же квадратной формы, что и на изображении. Может ли это быть связано с SetDataExtent и SetWholeExtent или с интервалом?

2. Привет, Дэйв, я тоже попробую, еще раз спасибо.