Hello everybody.
I am trying to visualize the slices of a DICOM file together with a VTU within the same scene. Said VTU is a file resulting from performing FEM simulations using an STL geometry that I obtained from manually segmenting the bone in the corresponding DICOM. I am taking as a base code the example MedicalDemo3 (https://examples.vtk.org/site/Python/Medical/MedicalDemo3/) removing, to simplify the code, everything related to the extraction of the ‘bone’ and the ‘skin’ (leaving only what has to do with the reader and the three anatomical planes). The idea is that the VTU is visualized in the same way as the ‘skin’: a 3D geometry overlapped with the DICOM slices.
Above the example code, I modified the reader so that it can read existing DICOMs in a directory and added the VTU reader, along with a couple more variables for display:
# VTK
# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.util.numpy_support import vtk_to_numpy
from vtkmodules.vtkIOImage import vtkDICOMImageReader
from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera
from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridReader
from vtkmodules.vtkCommonCore import (
vtkLookupTable,
)
from vtkmodules.vtkFiltersModeling import vtkOutlineFilter
from vtkmodules.vtkImagingCore import vtkImageMapToColors
from vtkmodules.vtkRenderingCore import (
vtkActor,
vtkCamera,
vtkImageActor,
vtkPolyDataMapper,
vtkDataSetMapper,
vtkRenderWindow,
vtkRenderWindowInteractor,
vtkRenderer
)
vtu_path="/path/to/file.vtu"
dicom_dir="/path/to/dicom/directory"
pixel_spacing=(0.25, 0.25)
slice_thickness=0.25
# -----------------------------------------------------------------------------
# VTK pipeline
# -----------------------------------------------------------------------------
# Set colors (in RGBA)
colors = vtkNamedColors()
colors.SetColor("BkgColor", [82, 87, 110, 255]) # paraview background color
colors.SetColor("SkinColor", [240, 184, 160, 255])
# Create the renderer, the render window, and the interactor
renderer = vtkRenderer()
render_window = vtkRenderWindow()
render_window.AddRenderer(renderer)
render_window_interactor = vtkRenderWindowInteractor()
interactor_style = vtkInteractorStyleTrackballCamera()
render_window_interactor.SetInteractorStyle(interactor_style)
render_window_interactor.SetRenderWindow(render_window)
# Set a background color for the renderer and set the size of the
# render window (expressed in pixels).
renderer.SetBackground(colors.GetColor3d("BkgColor"))
render_window.SetSize(640, 480)
# ***** Read VTU *****
vtu_reader = vtkXMLUnstructuredGridReader()
vtu_reader.SetFileName(vtu_path)
vtu_reader.Update()
# Dataset mapper
dataset_mapper = vtkDataSetMapper()
dataset_mapper.SetInputConnection(vtu_reader.GetOutputPort())
dataset_mapper.ScalarVisibilityOn()
scalar_range = vtu_reader.GetOutput().GetScalarRange()
dataset_mapper.SetScalarRange(scalar_range)
# Color map
dataset_lut = dataset_mapper.GetLookupTable()
dataset_lut.SetHueRange(0.666, 0.0) # H - hue
dataset_lut.SetSaturationRange(1.0, 1.0) # S - saturation
dataset_lut.SetValueRange(1.0, 1.0) # V - value
dataset_lut.Build()
dataset_mapper.SetLookupTable(dataset_lut)
# Mesh
mesh = vtkActor()
mesh.SetMapper(dataset_mapper)
# Mesh: Setup default representation to surface
mesh.GetProperty().SetRepresentationToSurface()
mesh.GetProperty().SetPointSize(1)
mesh.GetProperty().EdgeVisibilityOff()
print(f"Mesh bounds: {mesh.GetBounds()}")
# ***** Read DICOM slices *****
dicom_reader = vtkDICOMImageReader()
dicom_reader.SetDirectoryName(dicom_dir)
dicom_reader.Update()
# Compute bounds
outline_data = vtkOutlineFilter()
outline_data.SetInputConnection(dicom_reader.GetOutputPort())
outline_data.Update()
print(f"DICOM bounds: {outline_data.GetOutput().GetBounds()}")
map_outline = vtkPolyDataMapper()
map_outline.SetInputConnection(outline_data.GetOutputPort())
outline = vtkActor()
outline.SetMapper(map_outline)
outline.GetProperty().SetColor(colors.GetColor3d("Black"))
print(f"Slice thickness [mm]: {slice_thickness}")
print(f"Pixel spacing [mm]: {pixel_spacing}")
# Check images array
image = dicom_reader.GetOutput()
rows, cols, slices = image.GetDimensions()
image_data = vtk_to_numpy(image.GetPointData().GetScalars()).reshape(rows, cols, -1)
print(f"DICOM image data: [shape={image_data.shape}, dtype={image_data.dtype}, "
f"min={image_data.min()}, max={image_data.max()}]")
# Get the number of slices for each anatomical plane
n_sagittal_slices = image_data.shape[0]
n_coronal_slices = image_data.shape[1]
n_axial_slices = image_data.shape[2]
# Get the bounds to use in SetDisplayExtent()
dicom_sagittal_min = 0
dicom_sagittal_max = int(n_sagittal_slices - 1)
dicom_sagittal = int((dicom_sagittal_max + dicom_sagittal_min) / 2) # displays the middle slice in sagittal plane
dicom_axial_min = 0
dicom_axial_max = int(n_axial_slices - 1)
dicom_axial = int((dicom_axial_max + dicom_axial_min) / 2) # displays the middle slice in sagittal plane
dicom_coronal_min = 0
dicom_coronal_max = int(n_coronal_slices - 1)
dicom_coronal = int((dicom_coronal_max + dicom_coronal_min) / 2) # displays the middle slice in sagittal plane
# Start by creating a black/white lookup table.
bw_lut = vtkLookupTable()
bw_lut.SetTableRange(0, 2000)
bw_lut.SetSaturationRange(0, 0) # no color saturation
bw_lut.SetHueRange(0, 0)
bw_lut.SetValueRange(0, 1) # from black to white
bw_lut.Build() # effective built
# Create SAGITTAL slices
sagittal_colors = vtkImageMapToColors()
sagittal_colors.SetInputConnection(dicom_reader.GetOutputPort())
sagittal_colors.SetLookupTable(bw_lut)
sagittal_colors.Update()
sagittal = vtkImageActor()
sagittal.GetMapper().SetInputConnection(sagittal_colors.GetOutputPort())
sagittal.SetDisplayExtent(dicom_sagittal, dicom_sagittal,
dicom_coronal_min, dicom_coronal_max,
dicom_axial_min, dicom_axial_max)
sagittal.ForceOpaqueOn()
# Create AXIAL slices
axial_colors = vtkImageMapToColors()
axial_colors.SetInputConnection(dicom_reader.GetOutputPort())
axial_colors.SetLookupTable(bw_lut)
axial_colors.Update()
axial = vtkImageActor()
axial.GetMapper().SetInputConnection(axial_colors.GetOutputPort())
axial.SetDisplayExtent(dicom_sagittal_min, dicom_sagittal_max, dicom_coronal_min, dicom_coronal_max, dicom_axial,
dicom_axial)
axial.ForceOpaqueOn()
# Create CORONAL slices
coronal_colors = vtkImageMapToColors()
coronal_colors.SetInputConnection(dicom_reader.GetOutputPort())
coronal_colors.SetLookupTable(bw_lut)
coronal_colors.Update()
coronal = vtkImageActor()
coronal.GetMapper().SetInputConnection(coronal_colors.GetOutputPort())
coronal.SetDisplayExtent(dicom_sagittal_min, dicom_sagittal_max, dicom_coronal, dicom_coronal, dicom_axial_min,
dicom_axial_max)
coronal.ForceOpaqueOn()
# It is convenient to create an initial view of the data. The
# FocalPoint and Position form a vector direction. Later on
# (ResetCamera() method) this vector is used to position the camera
# to look at the data in this direction.
a_camera = vtkCamera()
a_camera.SetViewUp(0, 0, -1)
a_camera.SetPosition(0, -1, 0)
a_camera.SetFocalPoint(0, 0, 0)
a_camera.ComputeViewPlaneNormal()
a_camera.Azimuth(30.0)
a_camera.Elevation(30.0)
# Actors are added to the renderer.
renderer.AddActor(mesh)
renderer.AddActor(outline)
renderer.AddActor(sagittal)
renderer.AddActor(axial)
renderer.AddActor(coronal)
# An initial camera view is created. The Dolly() method moves
# the camera towards the FocalPoint, thereby enlarging the image.
renderer.SetActiveCamera(a_camera)
# Calling Render() directly on a vtkRenderer is strictly forbidden.
# Only calling Render() on the vtkRenderWindow is a valid call.
render_window.SetWindowName("MedicalDemo3")
render_window.Render()
renderer.ResetCamera()
a_camera.Dolly(1.5)
# Note that when camera movement occurs (as it does in the Dolly()
# method), the clipping planes often need adjusting. Clipping planes
# consist of two planes: near and far along the view direction. The
# near plane clips out objects in front of the plane; the far plane
# clips out objects behind the plane. This way only what is drawn
# between the planes is actually rendered.
renderer.ResetCameraClippingRange()
# Interact with the data.
render_window.Render()
render_window_interactor.Initialize()
render_window_interactor.Start()
where the outputs of the prints are:
Mesh bounds: (-0.04964113235, 0.04310790077, -0.06312587112, 0.01191986725, 0.06700585037, 0.1173627079)
DICOM bounds: (0.0, 149.75, 0.0, 149.75, 0.0, 77.5)
Slice thickness [mm]: 0.25
Pixel spacing [mm]: (0.25, 0.25)
DICOM image data: [shape=(600, 600, 311), dtype=int16, min=-1000, max=3095]
I manage to show both things separately:
(In the image on the left, I commented the ‘outline’, ‘sagittal’, ‘coronal’ and ‘axial’ actors. In the image on the right, only the ‘mesh’ actor is commented)
The problem comes when I try to show everything at once, since only the DICOM slices can be seen. I understand it has to do with the scale of each thing: DICOM slices are in pixels and the associated VTU is in m. It seems reasonable then, to use the pixel spacing and the slice thickness to rescale the DICOM images to the size of the VTU scale. However, it is not possible since the SetDisplayExtent() method requires data of type integer, not float.
Another way I can think of, given the values of the ‘DICOM bounds’ output, is to rescale “something” (I’m not quite sure what) so that those limits are (almost) the same as the ‘Mesh bounds’. Does anyone know how this can be done? Should I approach it differently? I don’t want to rescale the VTU itself, as it will have realistic measurements in its final version.
Thanks in advance.
P.D.: I am using Python v3.8 and vtk==9.3.0