Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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)
@AidenZhu186

This comment has been minimized.

Copy link

commented Sep 3, 2019

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

@lassoan

This comment has been minimized.

Copy link
Owner Author

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
You can’t perform that action at this time.