Python script to get axis-aligned bounding box

Development
1

Could you tell me how to get an axis-aligned bounding box that fits a segment using python interactor? The script repository describes how to create an oriented bounding box, but not an axis-aligned bounding box. I would like to know a script to make a roi like attached image.
Thank you!

axis-aligned bounding box
axis-aligned bounding box699×987 79.9 KB

2

Probably the easiest is to get the segment as a numpy array as shown here and then use basic numpy function to get index range of non-zero elements.

3

Could you give a sample script because I am new to python.
Thank you in advance.

4

You can use the code snippet that I linked to get the numpy array, then search on Google/StackOverflow for code snippets to get range of non-zero elements in a numpy array.

5

I was able to make a script but it is very long. Can you make it simpler?

#Get volume and segmentation
volumeNode = slicer.mrmlScene.GetFirstNodeByClass("vtkMRMLScalarVolumeNode")
segmentationNode = slicer.mrmlScene. GetFirstNodeByClass ("vtkMRMLSegmentationNode")

# Export "Segment_1" to LabelmapNode
labelmapVolumeNode=slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode")
segmentIds = vtk.vtkStringArray()
segmentId = segmentationNode.GetSegmentation().GetSegmentIdBySegmentName("Segment_1")
segmentIds.InsertNextValue(segmentId)
slicer.modules.segmentations.logic().ExportSegmentsToLabelmapNode(segmentationNode, segmentIds, labelmapVolumeNode, volumeNode)

#Get 3 axis length and center coordinate (kJI coordinate) of a bounding box
import numpy as np
label = array("LabelMapVolume")
points  = np.where( label == 1 )
length_z=np.max(points[0])- np.min(points[0])+1
length_y=np.max(points[1])- np.min(points[1])+1
length_x=np.max(points[2])- np.min(points[2])+1
center=[(np.max(points[0])+ np.min(points[0]))/2, (np.max(points[1])+ np.min(points[1]))/2, (np.max(points[2])+ np.min(points[2]))/2]

#Transform KJI to RAS coordinate
ijkToRasMatrix = vtk.vtkMatrix4x4()
volumeNode.GetIJKToRASMatrix(ijkToRasMatrix)
ijkPoint = np.array([center[2], center[1], center[0]])  # This is IJK, opposite order compared to KJI.
rasPoint = ijkToRasMatrix.MultiplyFloatPoint(np.append(ijkPoint,1))

#Make an axis-alinged bounding box of "Segment_1". 
roi=slicer.mrmlScene.AddNewNodeByClass("vtkMRMLAnnotationROINode")
roi.SetXYZ(rasPoint[0], rasPoint[1], rasPoint[2])
spacing = volumeNode.GetSpacing()
roi.SetRadiusXYZ(length_x*spacing[0]/2, length_y*spacing[1]/2, length_z*spacing[2]/2)
6

This looks simple enough. Note that the Annotations module is deprecated and will be removed in Slicer-5.4, so I would recommend to use Markups module instead.

7

I changed to use markups module and leave the script below for other users.
Thank you for your help.

#Get volume and segmentation. 
volumeNode = slicer.mrmlScene.GetFirstNodeByClass("vtkMRMLScalarVolumeNode")
segmentationNode = slicer.mrmlScene. GetFirstNodeByClass ("vtkMRMLSegmentationNode")

# Export “Segment_1” to LabelmapNode
labelmapVolumeNode=slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode")
segmentIds = vtk.vtkStringArray()
segmentId = segmentationNode.GetSegmentation().GetSegmentIdBySegmentName("Segment_1")
segmentIds.InsertNextValue(segmentId)
slicer.modules.segmentations.logic().ExportSegmentsToLabelmapNode(segmentationNode, segmentIds, labelmapVolumeNode, volumeNode)

#Get three axis length and center coordinate (kJI coordinate) of a bounding box
import numpy as np
label = array("LabelMapVolume")
points  = np.where( label == 1 )
length_z=np.max(points[0])- np.min(points[0])+1
length_y=np.max(points[1])- np.min(points[1])+1
length_x=np.max(points[2])- np.min(points[2])+1
center=[(np.max(points[0])+ np.min(points[0]))/2, (np.max(points[1])+ np.min(points[1]))/2, (np.max(points[2])+ np.min(points[2]))/2]

#Transform KJI to RAS coordinate
ijkToRasMatrix = vtk.vtkMatrix4x4()
volumeNode.GetIJKToRASMatrix(ijkToRasMatrix)
ijkPoint = np.array([center[2], center[1], center[0]])  # This is IJK, opposite order compared to KJI.
rasPoint = ijkToRasMatrix.MultiplyFloatPoint(np.append(ijkPoint,1))
print(rasPoint)

#Make an axis-alinged bounding box of "Segment_1". 
roi=slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsROINode")
roi.SetXYZ(rasPoint[0], rasPoint[1], rasPoint[2])
spacing = volumeNode.GetSpacing()
roi.SetRadiusXYZ(length_x*spacing[0]/2, length_y*spacing[1]/2, length_z*spacing[2]/2)
1 Like