Colored volume rendering

coloredvolumerendering.py

# This script can be used with this example scene:
# https://github.com/lassoan/PublicTestingData/releases/download/data/ColoredVolumeRenderingScene.mrb

volumeNode = getNode('Panoramix-cropped')
segmentationNode = getNode('Panoramix-cropped segmentation')
boostSegmentOpacity = False  # Increase voxel value inside segments to allow making them them more opaque

volumesLogic = slicer.modules.volumes.logic()

segmentIds = vtk.vtkStringArray()
labelmapVolumeNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode", "__temp__")
if not slicer.modules.segmentations.logic().ExportSegmentsToLabelmapNode(segmentationNode, segmentIds, labelmapVolumeNode, volumeNode):
    raise RuntimeError("Export of segment failed.")
colorTableNode = labelmapVolumeNode.GetDisplayNode().GetColorNode()

mapToRGB = vtk.vtkImageMapToColors()
mapToRGB.SetOutputFormatToRGB()
mapToRGB.SetInputData(labelmapVolumeNode.GetImageData())
mapToRGB.SetLookupTable(colorTableNode.GetLookupTable())
mapToRGB.Update()

smooth = vtk.vtkImageGaussianSmooth()
smooth.SetInputData(mapToRGB.GetOutput())
smooth.SetRadiusFactor(3.0)
smooth.Update()

shiftScale = vtk.vtkImageShiftScale()
shiftScale.SetOutputScalarType(vtk.VTK_UNSIGNED_CHAR)
shiftScale.SetScale(0.08)
shiftScale.SetShift(1000)
shiftScale.ClampOverflowOn()

if boostSegmentOpacity:
    # Create volume where segment region voxel intensities are increased
    boostedVolumeNode = volumesLogic.CloneVolume(slicer.mrmlScene, volumeNode, f"{volumeNode.GetName()} boosted")
    voxels = slicer.util.arrayFromVolume(volumeNode)
    labelVoxels = slicer.util.arrayFromVolume(labelmapVolumeNode)
    boostedVoxels = slicer.util.arrayFromVolume(boostedVolumeNode)
    boostedVoxels[labelVoxels>0] = voxels[labelVoxels>0] + 500
    slicer.util.arrayFromVolumeModified(boostedVolumeNode)
    # Smooth the boosted volume to avoid sharp edges at the segment boundary
    smooth2 = vtk.vtkImageGaussianSmooth()
    smooth2.SetInputData(boostedVolumeNode.GetImageData())
    smooth2.SetRadiusFactor(3.0)
    smooth2.Update()
    boostedVolumeNode.SetAndObserveImageData(smooth2.GetOutput())
    # Use boosted volume as alpha
    shiftScale.SetInputData(boostedVolumeNode.GetImageData())
else:
    # Use original volume as alpha
    shiftScale.SetInputData(volumeNode.GetImageData())

shiftScale.Update()

appendComponents = vtk.vtkImageAppendComponents()
appendComponents.AddInputData(smooth.GetOutput())
appendComponents.AddInputData(shiftScale.GetOutput())
appendComponents.Update()

coloredVolumeNode = volumesLogic.CloneVolume(slicer.mrmlScene, volumeNode, f"{volumeNode.GetName()} colored")
coloredVolumeNode.SetAndObserveImageData(appendComponents.GetOutput())

#slicer.mrmlScene.RemoveNode(colorTableNode)
#slicer.mrmlScene.RemoveNode(labelmapVolumeNode)
#slicer.mrmlScene.RemoveNode(boostedVolumeNode)