Number of voxels dependent on HU value

Support
, ,
1

Hi all,

I created a segment of volume and got the histogram of segment. What I need to do next is to have statistical data about segment such as; How many voxels above/under 2000 HU value? And their percentage?
I could again create 2 segments with 2000 being the boundary, but I just need the number. Is there a way that I can assing the number of voxels under 2000 HU to a variable in possibly Python Interactor.

Best regards.

Code im using:

##just histogram
import numpy as np
import SampleData 

masterVolumeNode=getNode('inputVolumeNode')
segmentationNode=getNode('Segmentation')



# Create segment editor to get access to effects
segmentEditorWidget = slicer.qMRMLSegmentEditorWidget()   #segment editor modulüne erişmek için
# 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)    #segmentation
segmentEditorWidget.setMasterVolumeNode(masterVolumeNode)    #segmentvolume

# Set up masking parameters
segmentEditorWidget.setActiveEffectByName("Mask volume")   #mask volume kısmı histogram için segmentasyonu masklemelisin
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()):#for segmentIndex in segmentsayisi
  # 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=400, 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)

image
image1920×1030 225 KB

2

You can convert a volume masked by a segmentation to a numpy array with slicer.util.arrayFromVolume and compute the number from the array.

2 Likes
3

You have already solved the problem. You can get the histogram with any resolution and then sum values of the histogram to get total number of voxels in a certain intensity range.

Alternatively, you can find a complete implementation of this feature as a module with a nice GUI in LungCTAnalyzer extension: it computes volume of each region that is specified by voxel intensity range. It shows nice real-time preview, generates, PDF report, etc. To get started, you can use the module as is, just ignoring that takes “left lung” and “right lung” segments as inputs. You can put any other structure in those segments and you can specify any intensity range. If it is close enough to what you want to achieve, you can easily customize it to be more convenient and intuitive for your data sets - it is just a Python scripted module that you can modify with a text editor.

1 Like
5

import numpy as np
import SampleData
masterVolumeNode=getNode(‘inputVolumeNode_masked’)
segmentationNode=getNode(‘Segmentation’)

Create segment editor to get access to effects

segmentEditorWidget = slicer.qMRMLSegmentEditorWidget() segment editor modulüne erişmek için

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) segmentation
segmentEditorWidget.setMasterVolumeNode(masterVolumeNode) #segmentvolume

Set up masking parameters

segmentEditorWidget.setActiveEffectByName(“Mask volume”) #mask volume kısmı histogram için segmentasyonu masklemelisin
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()):#for segmentIndex in segmentsayisi

Set active segment

segmentID = segmentationNode.GetSegmentation().GetNthSegmentID(segmentIndex)
segmentEditorWidget.setCurrentSegmentID(segmentID)

Apply mask

effect.self().onApply()

Compute histogram values

histogram_upper = np.histogram(arrayFromVolume(maskedVolume), bins=1, range=(2000,np.amax(intensityRange)))

Save results to a new table node

segment = segmentationNode.GetSegmentation().GetNthSegment(segmentIndex)
tableNode=slicer.mrmlScene.AddNewNodeByClass(“vtkMRMLTableNode”, segment.GetName() + " histogram_upper table")
updateTableFromArray(tableNode, histogram_upper)
tableNode.GetTable().GetColumn(0).SetName(“Count”)
tableNode.GetTable().GetColumn(1).SetName(“Intensity”)

Create new plot data series node

#histogram_lower
histogram_lower= np.histogram(arrayFromVolume(maskedVolume), bins=1, range=(np.amin(intensityRange),1999))
tableNode=slicer.mrmlScene.AddNewNodeByClass(“vtkMRMLTableNode”, segment.GetName() + " histogram_lower table")
updateTableFromArray(tableNode, histogram_lower)
tableNode.GetTable().GetColumn(0).SetName(“Count”)
tableNode.GetTable().GetColumn(1).SetName(“Intensity”)

Create new plot data series node

Delete temporary node

slicer.mrmlScene.RemoveNode(maskedVolume)
slicer.mrmlScene.RemoveNode(segmentEditorNode)

Fig.1: Number of voxels under 2000HU

image
image1105×753 109 KB

Fig.2: Number of voxels above 2000HU
image
image1052×768 104 KB

Thank you

6

Since I have just started using this program, I am not very familiar yet. I don’t know how to use this.
Thank you