I made a 3D ROI with Segment Editor module in an MRI image… How can I get the area of each sagittal slice?
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.
You can use the Segment Cross Section Area module in Sandbox extension (listed under the Examples category in the Extension manager).
1 Like