Skip to content

Instantly share code, notes, and snippets.

@lassoan
Last active September 3, 2019 16:08
Show Gist options
  • Save lassoan/2f5071c562108dac8efe277c78f2620f to your computer and use it in GitHub Desktop.
Save lassoan/2f5071c562108dac8efe277c78f2620f to your computer and use it in GitHub Desktop.
Histogram plot of volume masked by segments
# Generate input data
################################################
import SampleData
import numpy as np
# Load master volume
sampleDataLogic = SampleData.SampleDataLogic()
masterVolumeNode = sampleDataLogic.downloadMRBrainTumor1()
# Create segmentation
segmentationNode = slicer.vtkMRMLSegmentationNode()
slicer.mrmlScene.AddNode(segmentationNode)
segmentationNode.CreateDefaultDisplayNodes() # only needed for display
segmentationNode.SetReferenceImageGeometryParameterFromVolumeNode(masterVolumeNode)
# Create seed segment inside tumor
tumorSeed = vtk.vtkSphereSource()
tumorSeed.SetCenter(-6, 30, 28)
tumorSeed.SetRadius(10)
tumorSeed.Update()
segmentationNode.AddSegmentFromClosedSurfaceRepresentation(tumorSeed.GetOutput(), "Tumor", [1.0,0.0,0.0])
# Create seed segment inside tumor 2
referenceSeed = vtk.vtkSphereSource()
referenceSeed.SetCenter(-6, -50, -10)
referenceSeed.SetRadius(20)
referenceSeed.Update()
segmentationNode.AddSegmentFromClosedSurfaceRepresentation(referenceSeed.GetOutput(), "Reference", [0.0,0.0,1.0])
# Create seed segment outside tumor
backgroundSeedPositions = [[0,65,32], [1, -14, 30], [0, 28, -7], [0,30,64], [31, 33, 27], [-42, 30, 27]]
append = vtk.vtkAppendPolyData()
for backgroundSeedPosition in backgroundSeedPositions:
backgroundSeed = vtk.vtkSphereSource()
backgroundSeed.SetCenter(backgroundSeedPosition)
backgroundSeed.SetRadius(10)
backgroundSeed.Update()
append.AddInputData(backgroundSeed.GetOutput())
append.Update()
backgroundSegmentId = segmentationNode.AddSegmentFromClosedSurfaceRepresentation(append.GetOutput(), "Background", [0.0,1.0,0.0])
# Perform analysis
################################################
# Create segment editor to get access to effects
segmentEditorWidget = slicer.qMRMLSegmentEditorWidget()
# To show segment editor widget (useful for debugging): segmentEditorWidget.show()
segmentEditorWidget.setMRMLScene(slicer.mrmlScene)
segmentEditorNode = slicer.vtkMRMLSegmentEditorNode()
slicer.mrmlScene.AddNode(segmentEditorNode)
segmentEditorWidget.setMRMLSegmentEditorNode(segmentEditorNode)
segmentEditorWidget.setSegmentationNode(segmentationNode)
segmentEditorWidget.setMasterVolumeNode(masterVolumeNode)
# Set up masking parameters
segmentEditorWidget.setActiveEffectByName("Mask volume")
effect = segmentEditorWidget.activeEffect()
# set fill value to be outside the valid intensity range
intensityRange = masterVolumeNode.GetImageData().GetScalarRange()
effect.setParameter("FillValue", str(intensityRange[0]-1))
# Blank out voxels that are outside the segment
effect.setParameter("Operation", "FILL_OUTSIDE")
# Create a volume that will store temporary masked volumes
maskedVolume = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode", "Temporary masked volume")
effect.self().outputVolumeSelector.setCurrentNode(maskedVolume)
# Create chart
plotChartNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLPlotChartNode", "Histogram")
# Create histogram plot data series for each masked volume
for segmentIndex in range(segmentationNode.GetSegmentation().GetNumberOfSegments()):
# Set active segment
segmentID = segmentationNode.GetSegmentation().GetNthSegmentID(segmentIndex)
segmentEditorWidget.setCurrentSegmentID(segmentID)
# Apply mask
effect.self().onApply()
# Compute histogram values
histogram = np.histogram(arrayFromVolume(maskedVolume), bins=100, range=intensityRange)
# Save results to a new table node
segment = segmentationNode.GetSegmentation().GetNthSegment(segmentIndex)
tableNode=slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTableNode", segment.GetName() + " histogram table")
updateTableFromArray(tableNode, histogram)
tableNode.GetTable().GetColumn(0).SetName("Count")
tableNode.GetTable().GetColumn(1).SetName("Intensity")
# Create new plot data series node
plotSeriesNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLPlotSeriesNode", segment.GetName() + " histogram")
plotSeriesNode.SetAndObserveTableNodeID(tableNode.GetID())
plotSeriesNode.SetXColumnName("Intensity")
plotSeriesNode.SetYColumnName("Count")
plotSeriesNode.SetPlotType(slicer.vtkMRMLPlotSeriesNode.PlotTypeScatter)
plotSeriesNode.SetMarkerStyle(slicer.vtkMRMLPlotSeriesNode.MarkerStyleNone)
plotSeriesNode.SetUniqueColor()
# Add plot to chart
plotChartNode.AddAndObservePlotSeriesNodeID(plotSeriesNode.GetID())
# Show chart in layout
slicer.modules.plots.logic().ShowChartInLayout(plotChartNode)
# Delete temporary node
slicer.mrmlScene.RemoveNode(maskedVolume)
slicer.mrmlScene.RemoveNode(segmentEditorNode)
@lassoan
Copy link
Author

lassoan commented Sep 3, 2019

Line-71: plotSeriesNode is not defined yet.
Need be deleted.

Good catch, thank you, fixed in rev3.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment