Area of each slice

Support
1

I made a 3D ROI with Segment Editor module in an MRI image… How can I get the area of each sagittal slice?

2

To do this I would convert the segmentation to a labelmap, count the labeled voxels across the sagittal plane, and multiply by the voxel cross sectional area. Here is code which would do this for you:

segmentationName= 'MySegmentation' # change this to whatever your segmentation is called
segmentationNode = getNode(segmentationName)
labelMapVolumeNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLLabelMapVolumeNode', segmentationNode.GetName()+'-label')
slicer.modules.segmentations.logic().ExportAllSegmentsToLabelmapNode(segmentationNode, labelMapVolumeNode , slicer.vtkSegmentation.EXTENT_REFERENCE_GEOMETRY)
labelMapNumpyArray = slicer.util.arrayFromVolume(labelMapVolumeNode )
voxelSizeMm = labelMapVolumeNode.GetSpacing()

Next, you can find the area by summing all the non-zero labeled voxels in each sagittal slice and multiplying by the voxel area in that slice plane. For this, you need to know which IJK dimension corresponds to your sagittal plane because there can be an arbitrary relationship between the I, J, and K dimensions and the anatomical directions. The relationship between the IJK directions and the RAS directions can be seen by examining the volume’s IJK to RAS directions matrix in the “Volumes” module, in the “Volume Information” section. If you were to determine that the sagittal slice varies in the K direction (i.e. by seeing that the K direction vector points in the anatomical “Right” direction), then the following code would find the areas you want.

mm2PerVoxel = voxelSizeMm[0] * voxelSizeMm[1] # the cross sectional area of the voxel in the IJ plane
# The numpy indexing is KJI rather than IJK, so to sum across the IJ plane, we need to sum across axis 1 and 2 rather than 0 and 1
areaInVoxels = np.sum(labelMapNumpyArray > 0, axis=(1,2))
areaInMm2 = areaInVoxels * mm2PerVoxel

If you saw that your sagittal slice varied in the J dimension instead, then you would want to do this:

mm2PerVoxel = voxelSizeMm[0] * voxelSizeMm[2] # the cross sectional area of the voxel in the IK plane
# The numpy indexing is KJI rather than IJK, so to sum across the IK plane, we need to sum across axis 2 and 0 rather than 0 and 2
areaInVoxels = np.sum(labelMapNumpyArray > 0, axis=(0,2))
areaInMm2 = areaInVoxels * mm2PerVoxel

If you have multiple segments, you could do np.sum(labelMapNumpyArray==segmentNumber) for each segment number instead of the >0, which will include all segments.

At the end of any of these, areaInMm2 would be a numpy vector with as many values as there are sagittal slices.

3

You can use the Segment Cross Section Area module in Sandbox extension (listed under the Examples category in the Extension manager).

1 Like
4

I like this but it is not clear which axis (axial, coronal or sagittal) is chosen or can be chosen. The video does not indicate it as well.
After reviewing it seems “column” is for coronal, row is for sagittal, and slice is axial, correct ?
Many thanks

5

You can choose IJK axis in Volume axis selector. It is not anatomical axis, but row, column, or slice. If you want to get cross-sections along an arbitrary axis direction (including any anatomical axis) then there are several other options. To summarize, if you need cross-section areas:

  • in axial slices: you can use the Segment Cross-Section Area module in Sandbox extension.
  • along a line in arbitrary direction: you can use SegmentGeometry module in SlicerBiomech extension.
  • along a curve: you can segment the region and use Cross-section analysis or Stenosis measurement modules in VMTK extension.
  • in a single slice (either axial oblique): you can use a closed curve markup (enable area measurement in Measurements section in Markups module).
2 Likes