Last active
September 3, 2019 16:08
-
-
Save lassoan/2f5071c562108dac8efe277c78f2620f to your computer and use it in GitHub Desktop.
Histogram plot of volume masked by segments
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Good catch, thank you, fixed in rev3.