3D Dicom Visualisation in Python

Question

I am new to 3D image processing . I would like to know how to view the dicom series in python. I tried using matplotlib and VTK. In matplot I am not able to view the volume like I view in matlab using volViewer. Regarding VTK I am not able to import VTKRAyCASt for viewing 3D. The version I am using is 8.2.0.

I am doing the processing using scipy.ndimages

Kindly suggest me some resources to my volume dicom files

Respirator

Answer 1

You can try ipyvolume https://github.com/maartenbreddels/ipyvolume for interactive plotting, I found it quite useful. Also, you can plot them with matplotlib by using marching cubes to obtain the surface mesh but it is quite slow though:

from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import numpy as np
from skimage import measure

def plot_3d(image, threshold=-300): 
    p = image.transpose(2,1,0)
    verts, faces, normals, values = measure.marching_cubes_lewiner(p, threshold)
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    mesh = Poly3DCollection(verts[faces], alpha=0.1)
    face_color = [0.5, 0.5, 1]
    mesh.set_facecolor(face_color)
    ax.add_collection3d(mesh)
    ax.set_xlim(0, p.shape[0])
    ax.set_ylim(0, p.shape[1])
    ax.set_zlim(0, p.shape[2])

    plt.show()

The threshold of -300 HU is fine for visualizing chest CT scans but change it if you going to use MRI (check your intensity values distribution) or binary volumes (threshold =0).

There is some example of visualization:

Litharge

Answer 2

With vtkplotter you should be able to do this easily:

from vtkplotter import *

volume = load(mydicomdir) #returns a vtkVolume object
show(volume, bg='white')

To install: pip install vtkplotter

Photodynamics

vtkplotter looked promising, but worked on it for a half an hour and couldn't get it to not crash on .load with Process finished with exit code 139 (interrupted by signal 11: SIGSEGV) – Huberman

uhmm does it work with some test dataset? E.g. from vtkplotter import datadir, load then load(datadir+'embryo.tif').show() – Photodynamics

On two of my machines it was crashing when trying to load DICOM data, and also- there was no test dataset with my install (via pip). Also- this is a bad forum for bug reporting- happy to go to a mailing list or git issues. – Huberman

Feel free to submit to github.com/marcomusy/vtkplotter/issues – Photodynamics

Answer 3

This code can help you to create a 3D model from DICOM data using the VTK library!

import sys
import os
import itk
import vtk

# Define input and output directories directly as variables
dicom_directory = "Patient_Data"
output_image = "knee_3D.nrrd"
series_name = None  # Series name is set to None by default

print("Start: Reading and writing images from the DICOM directory")

# Default to current directory
dirName = dicom_directory if dicom_directory else "."
print(f"Using DICOM directory: {dirName}")

PixelType = itk.ctype("signed short")
Dimension = 3

ImageType = itk.Image[PixelType, Dimension]

# Use GDCMSeriesFileNames to get the DICOM filenames
namesGenerator = itk.GDCMSeriesFileNames.New()
namesGenerator.SetUseSeriesDetails(True)
namesGenerator.AddSeriesRestriction("0008|0021")
namesGenerator.SetGlobalWarningDisplay(False)
namesGenerator.SetDirectory(dirName)

print("Searching for DICOM files...")
seriesUID = namesGenerator.GetSeriesUIDs()

if len(seriesUID) < 1:
    print("No DICOMs in: " + dirName)
    sys.exit(1)

print("The directory: " + dirName)
print("Contains the following DICOM Series: ")
for uid in seriesUID:
    print(uid)

seriesFound = False
for uid in seriesUID:
    seriesIdentifier = uid
    if series_name:
        seriesIdentifier = series_name
        seriesFound = True
    print("Reading: " + seriesIdentifier)
    fileNames = namesGenerator.GetFileNames(seriesIdentifier)

    # Read the DICOM series
    print(f"Reading series: {seriesIdentifier}...")
    reader = itk.ImageSeriesReader[ImageType].New()
    dicomIO = itk.GDCMImageIO.New()
    reader.SetImageIO(dicomIO)
    reader.SetFileNames(fileNames)
    reader.ForceOrthogonalDirectionOff()

    # Write the 3D image
    print("Writing the 3D image...")
    writer = itk.ImageFileWriter[ImageType].New()
    outFileName = os.path.join(dirName, seriesIdentifier + ".nrrd")
    if output_image:
        outFileName = output_image
    writer.SetFileName(outFileName)
    writer.UseCompressionOn()
    writer.SetInput(reader.GetOutput())
    print("Writing: " + outFileName)
    writer.Update()
    print("3D image successfully written.")

    if seriesFound:
        break

Elisavetpol