Batch processing
Iterate through dicom series
This examples shows how to perform an operation on each series in the dicom database.
db = slicer.dicomDatabase
patients = db.patients()
patientCount = 0
for patient in patients:
patientCount += 1
print(f"Patient {patient} ({patientCount} of {len(patients)})")
for study in db.studiesForPatient(patient):
print(f"Study {study}")
for series in db.seriesForStudy(study):
print(f"Series {series}")
temporaryDir = qt.QTemporaryDir()
for instanceUID in db.instancesForSeries(series):
qt.QFile.copy(db.fileForInstance(instanceUID), f"{temporaryDir.path()}/{instanceUID}.dcm")
patientID = slicer.dicomDatabase.instanceValue(instanceUID, '0010,0020')
outputPath = os.path.join(convertedPath, patientID, study, series, "BatchResult")
if not os.path.exists(outputPath):
os.makedirs(outputPath)
# do an operation here that processes the series into the outputPath
Note
It can be helpful for debugging to include a comment with python commands that can be pasted into the console to run the script. With this approach any global variables, such as vtk class instances, defined in the script will exist after the script runs and you can easily inspect them or call methods on them.
"""
filePath = "/data/myscript.py"
exec(open(filePath).read())
"""
Extracting volume patches around segments
For machine learning apps, such as MONAI Label it can be helpful to reduce the size of the problem by extracting out subsets of data. This example shows how to iterate through a directory of segmentations, compute their bounding boxes, and save out new volumes and segmentations centered around the segmentation.
Here the ROI is aligned with the volume. See Segmentations for examples using oriented bounding boxes and other options.
import glob
import os
def segROI(segmentationNode):
# Compute bounding boxes
import SegmentStatistics
segStatLogic = SegmentStatistics.SegmentStatisticsLogic()
segStatLogic.getParameterNode().SetParameter("Segmentation", segmentationNode.GetID())
segStatLogic.getParameterNode().SetParameter("LabelmapSegmentStatisticsPlugin.obb_origin_ras.enabled",str(True))
segStatLogic.getParameterNode().SetParameter("LabelmapSegmentStatisticsPlugin.obb_diameter_mm.enabled",str(True))
segStatLogic.getParameterNode().SetParameter("LabelmapSegmentStatisticsPlugin.obb_direction_ras_x.enabled",str(True))
segStatLogic.getParameterNode().SetParameter("LabelmapSegmentStatisticsPlugin.obb_direction_ras_y.enabled",str(True))
segStatLogic.getParameterNode().SetParameter("LabelmapSegmentStatisticsPlugin.obb_direction_ras_z.enabled",str(True))
segStatLogic.computeStatistics()
stats = segStatLogic.getStatistics()
# Draw ROI for each oriented bounding box
import numpy as np
for segmentId in stats["SegmentIDs"]:
# Get bounding box
obb_origin_ras = np.array(stats[segmentId,"LabelmapSegmentStatisticsPlugin.obb_origin_ras"])
obb_diameter_mm = np.array(stats[segmentId,"LabelmapSegmentStatisticsPlugin.obb_diameter_mm"])
obb_direction_ras_x = np.array(stats[segmentId,"LabelmapSegmentStatisticsPlugin.obb_direction_ras_x"])
obb_direction_ras_y = np.array(stats[segmentId,"LabelmapSegmentStatisticsPlugin.obb_direction_ras_y"])
obb_direction_ras_z = np.array(stats[segmentId,"LabelmapSegmentStatisticsPlugin.obb_direction_ras_z"])
# Create ROI
segment = segmentationNode.GetSegmentation().GetSegment(segmentId)
roi=slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsROINode")
roi.SetName(segment.GetName() + " OBB")
roi.GetDisplayNode().SetHandlesInteractive(False) # do not let the user resize the box
roi.SetSize(obb_diameter_mm * 2) # make the ROI twice the size of the segmentation
# Position and orient ROI using a transform
obb_center_ras = obb_origin_ras+0.5*(obb_diameter_mm[0] * obb_direction_ras_x + obb_diameter_mm[1] * obb_direction_ras_y + obb_diameter_mm[2] * obb_direction_ras_z)
boundingBoxToRasTransform = np.row_stack((np.column_stack(((1,0,0), (0,1,0), (0,0,1), obb_center_ras)), (0, 0, 0, 1)))
boundingBoxToRasTransformMatrix = slicer.util.vtkMatrixFromArray(boundingBoxToRasTransform)
roi.SetAndObserveObjectToNodeMatrix(boundingBoxToRasTransformMatrix)
return roi
labelFiles = glob.glob("/data/imagesTr/labels/final/*.nii.gz")
for labelFile in labelFiles:
slicer.mrmlScene.Clear()
print(labelFile)
baseName = os.path.basename(labelFile)
ctFile = os.path.join("/data/imagesTr", baseName)
print(ctFile)
ct = slicer.util.loadVolume(ctFile)
seg = slicer.util.loadSegmentation(labelFile)
roi = segROI(seg)
cropVolumeParameters = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLCropVolumeParametersNode")
cropVolumeParameters.SetInputVolumeNodeID(ct.GetID())
cropVolumeParameters.SetROINodeID(roi.GetID())
slicer.modules.cropvolume.logic().Apply(cropVolumeParameters)
croppedCT = cropVolumeParameters.GetOutputVolumeNode()
seg.SetReferenceImageGeometryParameterFromVolumeNode(croppedCT)
segLogic = slicer.modules.segmentations.logic()
labelmap = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode")
segLogic.ExportAllSegmentsToLabelmapNode(seg, labelmap, slicer.vtkSegmentation.EXTENT_REFERENCE_GEOMETRY)
slicer.util.saveNode(croppedCT, f"/data/crops/{baseName}")
slicer.util.saveNode(labelmap, f"/data/crops/labels/final/{baseName}")
slicer.app.processEvents() # to watch progress