Using Python Module to Render an Image in 3D

[vertebrae]

Hello,

I am looking to use a 3D slicer python module to render thresholded spinal data in 3D and to potentially save as an STL file. Does anyone know how to do this?

[vertebra]

Hey! I would take a look at this repository. It has the code: https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository

[lassoan (Andras Lasso)]

More specifically, these examples show how to automatically segment structures, display it in 3D and save as STL file: https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#How_to_run_segment_editor_effects_from_a_script

[vertebra]

One more question Andras!

Whenever I attempt to run my code with any slicer.util or other commands, it says slicer is not defined. Even though I have imported slicer.

What’s going on?

[vertebra]

One more question Andras. Also, do you have any sources for coding the keepLargestIsland function?

[lassoan (Andras Lasso)]

@vertebra

Whenever I attempt to run my code with any slicer.util or other commands, it says slicer is not defined. Even though I have imported slicer.

How do you run your script?

[vertebra]

I have a module in 3D Slicer but the code is just written in IDLE. When running it from IDLE to print something, it says you have to import from slicer.

[vertebrae]

I also have the same question regarding the islands function.

[lassoan (Andras Lasso)]

@vertebrae

I also have the same question regarding the islands function.

Entire source code of Slicer is open. Implementaitno of all built-in Segment editor effects are available here: Slicer/Modules/Loadable/Segmentations/EditorEffects at main · Slicer/Slicer · GitHub

@vertebra

When running it from IDLE to print something, it says you have to import from slicer.

You need to run Python scripts that use Slicer classes from Slicer’s Python environment.

[vertebra]

I’m running my python module through slicer. Even when I import slicer before the class, an error saying there is no module named slicer appears. What should I do about this?

Also, kind of a long shot… Do you have or could you create a completed python module for vertebra segmentation using flood fill, threshold, smoothing, and fiducial nodes selected by the user? That would be a great help especially if you have comments explaining each function.

[vertebra]

Thanks Andras once again for all the help!

[lassoan (Andras Lasso)]

@vertebra

I’m running my python module through slicer. Even when I import slicer before the class, an error saying there is no module named slicer appears. What should I do about this?

I still don’t have a clear idea of what scripts you are trying to run where. Maybe you could record a screen capture video and share that, but probably it would be even better to follow steps described in Slicer programming tutorials, such as this.

@vertebra

Do you have or could you create a completed python module for vertebra segmentation using flood fill, threshold, smoothing, and fiducial nodes selected by the user? That would be a great help especially if you have comments explaining each function.

Slicer community members has been working a lot on vertebra segmentation, so we know quite a lot about potential approaches and capabilities and limitations of each.

There is no such thing as a perfect segmentation. The right solution depends on what is the clinical goal, expertise and time constraints of users, input image modality, quality, resolution, computational resources, amount of data sets, etc. If you can tell more about your project then we can give advice on potential approaches.

[vertebra]

So the general approach we have decided upon is:

1. The user would mark a fiducial point on the lumbar vertebra
2. The code creates a segment naming it lumbar vertebra
3. The threshold segment editor is used to set it between (150, 400) to isolate the bone
4. The islands segment editor is used to keep the largest island
5. Around the fiducial point, the flood fill segment editor is used to encompass the singular lumbar vertebra
6. The smoothing/median segment editor is used to smooth the image
7. Render the image in 3D and save to an STL file

We are open to any suggestions on how to improve the code or alter the algorithm

[lassoan (Andras Lasso)]

@vertebra

The threshold segment editor is used to set it between (150, 400) to isolate the bone

Isolation of vertebrae is the most challenging problem. If you are so to have so high-quality, high-resolution images that you can separate a vertebra by simple thresholding then the whole segmentation process can be very easy (potentially easier than the steps you described). Do you segment human clinical CTs? What is the image resolution? Can you post a few screenshots?

The two main improvements that I would consider:

• use local thresholding effect (provided by SegmentEditorExtraEffect extension): it implements step 3-5 and it allows using GrowFromSeed or Watershed transforms for more robust connectivity computation
• solidify the bone (fill in internal holes caused by lower density of cancellous bone) using Wrap Solidify effect (provided by SurfaceWrapSolidify extension), as it fill all holes without risk of leaking or missing some holes

[vertebra]

I mean we have been experimenting on CT Chest and CT Cardio from the sample data. If we only use one segment, how would grow from seeds or watershed be implemented? We will definitely try the wrap solidify effect. And we will start experimenting with localized thresholding.

[vertebrae]

Hello Andras,

I cannot find the SurfaceWrapSolidify extension in the extension manager. How can I access this?

[lassoan (Andras Lasso)]

@vertebra

If we only use one segment, how would grow from seeds or watershed be implemented?

For Grow from seeds and Watershed effects you always need to specify at least two segments. In case of vertebra segmentation, it would be a “bone” and “other” segment. Local thresholding effect automatically creates a temporary “other” segment from regions that have different intensity or they are disconnected.

@vertebra

we have been experimenting on CT Chest and CT Cardio from the sample data

On these images, you cannot separate a single vertebra by simple thresholding (but of course, as I mentioned, there is no such things “perfect segmentation”, so maybe what you can get from these images is good enough).

Here is an example of what you can get with a single-click segmentation using Local thresholding effect on CTACardio example (threshold range of 265-1009, features size = 6mm, segmentation algorithm = growcut):

[vertebra]

Thank you so much! This helps astoundingly. The localized threshold will be implemented instead of steps 3-5. You’re the greatest!

[vertebra]

Hello! We are still beginners at the coding process in python. Of course, this doesn’t work yet. We compiled some of the codes from different resources on GitHub that we identified. Could you please give us some tips or edits? We don’t have the rendering in yet.

import os
import vtk, qt, ctk, slicer
import logging
from SegmentEditorEffects import *
import vtkITK
import SimpleITK as sitk
import sitkUtils
import math
import vtkSegmentationCorePython as vtkSegmentationCore
import vtkSlicerSegmentationsModuleLogicPython as vtkSlicerSegmentationsModuleLogic
import SampleData

class VertebraSegmentation()

Load master volume

sampleDataLogic = SampleData.SampleDataLogic()
masterVolumeNode = sampleDataLogic.downloadCTACardio()

##gets the node coordinates to run the grow cut from later
hierarchy = slicer.vtkMMRLSubjectHierarchyNode(slicer.mmrlScene)
sceneItemID = hierarchy.GetSceneItemID()
subjectItemID = hierarchy.GetItemChildWithName(sceneItemID,‘Fiducial Nodes’)
fidList = slicer.util.getNode(‘FiducialNodes’)

Create segmentation

segmentationNode = slicer.vtkMRMLSegmentationNode()
slicer.mrmlScene.AddNode(segmentationNode)
segmentationNode.CreateDefaultDisplayNodes() # only needed for display
segmentationNode.name = ‘Lumbar’
segmentationNode.SetReferenceImageGeometryParameterFromVolumeNode(masterVolumeNode)

Create seed segment inside lumbar and name

lumbarSeed = fidList[0]
lumbarSeed.SetRadius(3)
lumbarSeed.Update()
segmentationNode.AddSegmentFromClosedSurfaceRepresentation(lumbarSeed.GetOutput(), “Lumbar Vertebra”, [1.0,0.0,0.0])

def init(self, scriptedEffect):
SegmentEditorThresholdEffect.init(self, scriptedEffect)
scriptedEffect.name = ‘Local Threshold’
self.previewSteps = 4

def updateMRMLFromGUI(self): ## sets the parameters for local threshold - grow cut at 6mm for feature size and 265.00 to 1009.00 for threshold range
SegmentEditorThresholdEffect.updateMRMLFromGUI(self)
featureSizeMm = 6.000 ##self.minimumMinimumFeatureSize.value
self.scriptedEffect.setParameter(MINIMUM_FEATURE_MM_PARAMETER_NAME, featureSizeMm)

segmentationAlgorithm = "GrowCut" ##self.segmentationAlgorithmSelector.currentText
self.scriptedEffect.setParameter(SEGMENTATION_ALGORITHM_PARAMETER_NAME, segmentationAlgorithm)

def runGrowCut(self, masterImageData, seedLabelmap, outputLabelmap): ## runs the grow cut - local threshold on the segment

self.clippedMaskImageData = slicer.vtkOrientedImageData()
intensityBasedMasking = self.scriptedEffect.parameterSetNode().GetMasterVolumeIntensityMask()
segmentationNode = self.scriptedEffect.parameterSetNode().GetSegmentationNode()
success = segmentationNode.GenerateEditMask(self.clippedMaskImageData,
  self.scriptedEffect.parameterSetNode().GetMaskMode(),
  masterImageData, # reference geometry
  "", # edited segment ID
  self.scriptedEffect.parameterSetNode().GetMaskSegmentID() if self.scriptedEffect.parameterSetNode().GetMaskSegmentID() else "",
  masterImageData if intensityBasedMasking else None,
  self.scriptedEffect.parameterSetNode().GetMasterVolumeIntensityMaskRange() if intensityBasedMasking else None)

import vtkSlicerSegmentationsModuleLogicPython as vtkSlicerSegmentationsModuleLogic
self.growCutFilter = vtkSlicerSegmentationsModuleLogic.vtkImageGrowCutSegment()
self.growCutFilter.SetIntensityVolume(masterImageData)
self.growCutFilter.SetSeedLabelVolume(seedLabelmap)
self.growCutFilter.SetMaskVolume(self.clippedMaskImageData)
self.growCutFilter.Update()
outputLabelmap.ShallowCopy(self.growCutFilter.GetOutput())

def apply(self, ijkPoints):
kernelSizePixel = self.getKernelSizePixel()
if kernelSizePixel[0]<=0 and kernelSizePixel[1]<=0 and kernelSizePixel[2]<=0:
return

qt.QApplication.setOverrideCursor(qt.Qt.WaitCursor)

# Get parameter set node
parameterSetNode = self.scriptedEffect.parameterSetNode()

# Get parameters
minimumThreshold = self.scriptedEffect.doubleParameter("MinimumThreshold")
maximumThreshold = self.scriptedEffect.doubleParameter("MaximumThreshold")

# Get modifier labelmap
modifierLabelmap = self.scriptedEffect.defaultModifierLabelmap()

# Get master volume image data
masterImageData = self.scriptedEffect.masterVolumeImageData()

# Set intensity range
oldMasterVolumeIntensityMask = parameterSetNode.GetMasterVolumeIntensityMask()
parameterSetNode.MasterVolumeIntensityMaskOn()
oldIntensityMaskRange = parameterSetNode.GetMasterVolumeIntensityMaskRange()
intensityRange = [265.00, 1009.00]
if oldMasterVolumeIntensityMask:
  intensityRange = [max(oldIntensityMaskRange[0], minimumThreshold), min(oldIntensityMaskRange[1], maximumThreshold)]
parameterSetNode.SetMasterVolumeIntensityMaskRange(intensityRange)

roiNode = lumbarSeed ##self.roiSelector.currentNode()
clippedMasterImageData = masterImageData
if roiNode is not None:
  worldToImageMatrix = vtk.vtkMatrix4x4()
  masterImageData.GetWorldToImageMatrix(worldToImageMatrix)

  bounds = [0,0,0,0,0,0]
  roiNode.GetRASBounds(bounds)
  corner1RAS = [bounds[0], bounds[2], bounds[4], 1]
  corner1IJK = [0, 0, 0, 0]
  worldToImageMatrix.MultiplyPoint(corner1RAS, corner1IJK)

  corner2RAS = [bounds[1], bounds[3], bounds[5], 1]
  corner2IJK = [0, 0, 0, 0]
  worldToImageMatrix.MultiplyPoint(corner2RAS, corner2IJK)

  extent = [0, -1, 0, -1, 0, -1]
  for i in range(3):
      lowerPoint = min(corner1IJK[i], corner2IJK[i])
      upperPoint = max(corner1IJK[i], corner2IJK[i])
      extent[2*i] = int(math.floor(lowerPoint))
      extent[2*i+1] = int(math.ceil(upperPoint))

  imageToWorldMatrix = vtk.vtkMatrix4x4()
  masterImageData.GetImageToWorldMatrix(imageToWorldMatrix)
  clippedMasterImageData = slicer.vtkOrientedImageData()
  self.padder = vtk.vtkImageConstantPad()
  self.padder.SetInputData(masterImageData)
  self.padder.SetOutputWholeExtent(extent)
  self.padder.Update()
  clippedMasterImageData.ShallowCopy(self.padder.GetOutput())
  clippedMasterImageData.SetImageToWorldMatrix(imageToWorldMatrix)

# Pipeline
self.thresh = vtk.vtkImageThreshold()
self.thresh.SetInValue(LABEL_VALUE)
self.thresh.SetOutValue(BACKGROUND_VALUE)
self.thresh.SetInputData(clippedMasterImageData)
self.thresh.ThresholdBetween(minimumThreshold, maximumThreshold)
self.thresh.SetOutputScalarTypeToUnsignedChar()
self.thresh.Update()

self.erode = vtk.vtkImageDilateErode3D()
self.erode.SetInputConnection(self.thresh.GetOutputPort())
self.erode.SetDilateValue(BACKGROUND_VALUE)
self.erode.SetErodeValue(LABEL_VALUE)
self.erode.SetKernelSize(
  kernelSizePixel[0],
  kernelSizePixel[1],
  kernelSizePixel[2])

self.erodeCast = vtk.vtkImageCast()
self.erodeCast.SetInputConnection(self.erode.GetOutputPort())
self.erodeCast.SetOutputScalarTypeToUnsignedInt()
self.erodeCast.Update()

# Remove small islands
self.islandMath = vtkITK.vtkITKIslandMath()
self.islandMath.SetInputConnection(self.erodeCast.GetOutputPort())
self.islandMath.SetFullyConnected(False)
self.islandMath.SetMinimumSize(125)  # remove regions smaller than 5x5x5 voxels

self.islandThreshold = vtk.vtkImageThreshold()
self.islandThreshold.SetInputConnection(self.islandMath.GetOutputPort())
self.islandThreshold.ThresholdByLower(BACKGROUND_VALUE)
self.islandThreshold.SetInValue(BACKGROUND_VALUE)
self.islandThreshold.SetOutValue(LABEL_VALUE)
self.islandThreshold.SetOutputScalarTypeToUnsignedChar()
self.islandThreshold.Update()

# Points may be outside the region after it is eroded.
# Snap the points to LABEL_VALUE voxels,
snappedIJKPoints = self.snapIJKPointsToLabel(ijkPoints, self.islandThreshold.GetOutput())
if snappedIJKPoints.GetNumberOfPoints() == 0:
  qt.QApplication.restoreOverrideCursor()
  return

# Convert points to real data coordinates. Required for vtkImageThresholdConnectivity.
seedPoints = vtk.vtkPoints()
origin = masterImageData.GetOrigin()
spacing = masterImageData.GetSpacing()
for i in range(snappedIJKPoints.GetNumberOfPoints()):
  ijkPoint = snappedIJKPoints.GetPoint(i)
  seedPoints.InsertNextPoint(
    origin[0]+ijkPoint[0]*spacing[0],
    origin[1]+ijkPoint[1]*spacing[1],
    origin[2]+ijkPoint[2]*spacing[2])

segmentationAlgorithm = self.scriptedEffect.parameter(SEGMENTATION_ALGORITHM_PARAMETER_NAME)
if segmentationAlgorithm == SEGMENTATION_ALGORITHM_MASKING:
  self.runMasking(seedPoints, self.islandThreshold.GetOutput(), modifierLabelmap)

else:
  self.floodFillingFilterIsland = vtk.vtkImageThresholdConnectivity()
  self.floodFillingFilterIsland.SetInputConnection(self.islandThreshold.GetOutputPort())
  self.floodFillingFilterIsland.SetInValue(SELECTED_ISLAND_VALUE)
  self.floodFillingFilterIsland.ReplaceInOn()
  self.floodFillingFilterIsland.ReplaceOutOff()
  self.floodFillingFilterIsland.ThresholdBetween(LABEL_VALUE, LABEL_VALUE)
  self.floodFillingFilterIsland.SetSeedPoints(seedPoints)
  self.floodFillingFilterIsland.Update()

  self.maskCast = vtk.vtkImageCast()
  self.maskCast.SetInputData(self.thresh.GetOutput())
  self.maskCast.SetOutputScalarTypeToUnsignedChar()
  self.maskCast.Update()

  self.imageMask = vtk.vtkImageMask()
  self.imageMask.SetInputConnection(self.floodFillingFilterIsland.GetOutputPort())
  self.imageMask.SetMaskedOutputValue(OUTSIDE_THRESHOLD_VALUE)
  self.imageMask.SetMaskInputData(self.maskCast.GetOutput())
  self.imageMask.Update()

  imageMaskOutput = slicer.vtkOrientedImageData()
  imageMaskOutput.ShallowCopy(self.imageMask.GetOutput())
  imageMaskOutput.CopyDirections(clippedMasterImageData)

  imageToWorldMatrix = vtk.vtkMatrix4x4()
  imageMaskOutput.GetImageToWorldMatrix(imageToWorldMatrix)

  segmentOutputLabelmap = slicer.vtkOrientedImageData()
  if segmentationAlgorithm == SEGMENTATION_ALGORITHM_GROWCUT:
    self.runGrowCut(clippedMasterImageData, imageMaskOutput, segmentOutputLabelmap)
  elif segmentationAlgorithm == SEGMENTATION_ALGORITHM_WATERSHED:
    self.runWatershed(clippedMasterImageData, imageMaskOutput, segmentOutputLabelmap)
  else:
    logging.error("Unknown segmentation algorithm: \"" + segmentationAlgorithm + "\"")

  segmentOutputLabelmap.SetImageToWorldMatrix(imageToWorldMatrix)

  self.selectedSegmentThreshold = vtk.vtkImageThreshold()
  self.selectedSegmentThreshold.SetInputData(segmentOutputLabelmap)
  self.selectedSegmentThreshold.ThresholdBetween(SELECTED_ISLAND_VALUE, SELECTED_ISLAND_VALUE)
  self.selectedSegmentThreshold.SetInValue(LABEL_VALUE)
  self.selectedSegmentThreshold.SetOutValue(BACKGROUND_VALUE)
  self.selectedSegmentThreshold.SetOutputScalarType(modifierLabelmap.GetScalarType())
  self.selectedSegmentThreshold.Update()
  modifierLabelmap.ShallowCopy(self.selectedSegmentThreshold.GetOutput())

self.scriptedEffect.saveStateForUndo()
self.scriptedEffect.modifySelectedSegmentByLabelmap(modifierLabelmap, slicer.qSlicerSegmentEditorAbstractEffect.ModificationModeAdd)

parameterSetNode.SetMasterVolumeIntensityMask(oldMasterVolumeIntensityMask)
parameterSetNode.SetMasterVolumeIntensityMaskRange(oldIntensityMaskRange)

qt.QApplication.restoreOverrideCursor()

[lassoan (Andras Lasso)]

This code is way too long to be discussed here. If you need feedback, then create an extension using Extension Wizard module, upload it to github and just post links here. Give a high-level overview of what you are trying to achieve and ask specific questions related to specific parts of the code.