Hello,
please is there any way, how to get the histogram of intensities of only a bordered area by segmentation? I am working with MRI scans, one is 3D and serves for the segmantation of an organ and the second, which is only 3 slices, is a relaxation map. I would like to obtain intensities histograms of chosen organs and compare pathological and physiological values
many thanks
You can use Mask Volume effect (in SegmentEditorExtraEffects extension) to blank out areas of a volume that is outside of a segment. You can then use VTK or numpy to compute histogram (see example of computing and displaying histogram of a volume in Python) .
Additionally, if you want to compare healthy and pathological using statistics calculated on the histograms, you can use the SlicerRadiomics extension, which is the Slicer module for PyRadiomics. Features in class âFirstOrderâ are the histogram features and are calculated on the histogram of the Region of Interest, which you can define using segments or binary label maps.
Thank you very much for your reply!
I would like to ask one more thing: can I really blank out areas out of my volume of interest? They become black and that causes a big column in the histogram and I would like them to be of no value
You can set the âblankingâ intensity to any value and then ignore that value in the histogram.
Hi Diana,
I did something similar for my data. Assuming you start with a Volume node named âinputVolumeNodeâ and draw a segmentation ROI on it (default name âSegmentationâ) the following code should do it for you. Caveat: You need the latest nightly build of slicer to draw the plot :
import SimpleITK as sitk
import sitkUtils
import numpy as np
newScalarNode = slicer.mrmlScene.AddNewNodeByClass(âvtkMRMLScalarVolumeNodeâ) #to be used later
seg = getNode(âSegmentationâ)
#create a temporary labelmap
labelMapNode = slicer.mrmlScene.AddNewNodeByClass(âvtkMRMLLabelMapVolumeNodeâ)
#using reference volume will match the geometry
ids = vtk.vtkStringArray()
seg.GetSegmentation().GetSegmentIDs(ids)
slicer.modules.segmentations.logic().ExportSegmentsToLabelmapNode(seg,ids,labelMapNode,inputVolumeNode)
mif = sitk.MaskImageFilter()
mif.SetOutsideValue(-1234) #custom dummy value thatâs recognizable
roiSITK = sitkUtils.PullVolumeFromSlicer(inputVolumeNode)
maskSITK = sitkUtils.PullVolumeFromSlicer(labelMapNode)
final = mif.Execute(roiSITK,maskSITK)
sitkUtils.PushVolumeToSlicer(final,newScalarNode)
#creating histogram <-------------- from this onwards the code is from slicer scripts webpage
np_tmp = arrayFromVolume(newScalarNode)
np_tmp=np_tmp[np_tmp!=-1234] # flattens the array and removes selected outside label values
histogram = np.histogram(np_tmp, bins=50)
#Save results to a new table node
tableNode=slicer.mrmlScene.AddNewNodeByClass(âvtkMRMLTableNodeâ)
updateTableFromArray(tableNode, histogram)
tableNode.GetTable().GetColumn(0).SetName(âCountâ)
tableNode.GetTable().GetColumn(1).SetName(âIntensityâ)
#Create plot ------ DOES NOT EXIST ON PRE VERSION 4.9 !
plotSeriesNode = slicer.mrmlScene.AddNewNodeByClass(âvtkMRMLPlotSeriesNodeâ, newScalarNode.GetName() + â histogramâ)
plotSeriesNode.SetAndObserveTableNodeID(tableNode.GetID())
plotSeriesNode.SetXColumnName(âIntensityâ)
plotSeriesNode.SetYColumnName(âCountâ)
plotSeriesNode.SetPlotType(plotSeriesNode.PlotTypeScatterBar)
plotSeriesNode.SetColor(0, 0.6, 1.0)
#Create chart and add plot
plotChartNode = slicer.mrmlScene.AddNewNodeByClass(âvtkMRMLPlotChartNodeâ)
plotChartNode.AddAndObservePlotSeriesNodeID(plotSeriesNode.GetID())
plotChartNode.YAxisRangeAutoOff()
plotChartNode.SetYAxisRange(0, 500000)
#Show plot in layout
slicer.modules.plots.logic().ShowChartInLayout(plotChartNode)
Iâve also implemented this exact feature - see here: Histogram of volume in regions defined by segments
Probably it is very similar to what @drusmanbashir described above. Main difference is that I used SegmentEditorExtraEffects extensionâs Mask Volume effect to do the masking, which is somewhat more complicated, but works even if the segmentationâs master volume is not the same as the volume where you need compute the histograms on. There is also a small difference in the blanking value computation: I retrieve scalar range from the volume and select a blanking value that lies just outside this range; then specify the valid range in np.histogram to make sure the histogram bin distribution is not affected by the blanking value.
Thanks a lot Andras,
Your code appears to be a skeleton of a very useful module - one that would enable histogram visualisation of active segmentations.
Could you add a snipet that would update the histogram to reflect single axial view ROI (instead of entire segmentation map as current) stats as the user scrolls through the axial data? As a user I will also be interested to know the histogram of the current 2D circle I have just drawn over an object on the MRI image for example
Thanks
USMAN.
This is an example for slice-by-slice processing: https://gist.github.com/lassoan/d85be45b7284a3b4e580688fccdb1d02. It shows that you can extract a single slice of a volume using numpy indexing. You can apply the same technique for histogram computation: you donât pass the entire 3D array to np.histogram but you extract a slice first. So, with 2-3 extra lines of code you can get per-slice histograms.
Hello,
thank you for your replies! I am still working on my data, and you have already helped a lot, but there are still some things I would like to improve. When I want to compare I need some statistics of course, so I used the SlicerRadiomics extension. However, it works only if I re-segment the T1 map (it does not allow me to use the segmantation from the 3D volume on the T1 map), which I do in some cases anyway, because T1 maps and 3D volumes are a little bit shifted. Anyway, there is one more thing I would like to do: to unify the intensities range within the segments from 0 to 4000. That is the most frequent range I get, however sometimes it is different. I can set the number of bins, but that means the bins have different thickness. Can you help me work this out?
You can export segmentation to a labelmap in Segmentations module, Import/Export section.
To solve this, you can set bin range and width in numpy.histogram by bins and range parameters for all histograms, instead of relying on automatic bin generation.
Dear Diana,
What type of configuration are you using? On which version of slicer?
In the the latest version of SlicerRadiomics, the mask correction is set to True by default (allowing you to re-use your segmentation).
If that is not the case, but you are able to provide a parameter file customization, you can also use it that way (by specifying a configuration file which has the setting âcorrectMaskâ set to true). See here for more information on building a parameter file
Moreover, If you feed the image and segmentation to PyRadiomics, for histogram statistics, no discretization is applied (i.e. a bin width fixed to 1). Exception is uniformity and entropy, for that youâd need to set the bin width (default 25). Therefore, you would not need to unify your intensities.
example for a configuration file you can use (paste this in a text file and store with extension â.ymlâ):
settings:
binWidth: 1
correctMask: true
imageType:
Original: {}
featureClass:
firstorder:
Images may need smaller bin size (e.g., your entire scalar range may be between 0.0-1.0 if scalar type is float or double) or larger bin size (image is rescaled when loaded and there can be many intensity values that no voxel data is assigned to; so you could have large gaps, very noisy histogram if you use bin size = 1.0). So, in general, a fixed bin size of 1 is probably not an optimal choice.
@lassoan, you are correct. In that sense, mentioning bin width fixed to 1 is incorrect. In the actual code no discretization is applied, and it is therefore equivalent to calculating on a histogram with the number of bins equal to the number of unique gray values in the ROI. In fact, it just uses the values directly in calculating the features (it does not first make a histogram).
Hello, sorry for the belated reply. I am using the version 4.9. from the date 03.04. The reason why I think it does not work is that the relaxation map consists only of 3 slices, whereas the segmentation is from a 3d volume and it spatially exceeds the relaxation map. What I would like to do is to extract from the segmentation only the three slices, is it possible? So far, I resegmented the maps, but that is not ideal, because only the middle slice is segmented properly, the outer two are not (is it because there are no points above and below, so the algorithm âcanât growâ in these directions?) - I did not want the segmentation to exceed my ROI
thank you for your replies
When you export your segmentation to a labelmap (in Segmentations module, export section) then in the advanced section select your relaxation map as âReference volumeâ. The exported labelmap will have exactly the same size, spacing, and direction as the relaxation map.