Run Python code in each subfolder automatically

Support
1

I have this file structure:

  1. Main Folder: W001-W100.

    • Subfolder: W001

      • Around 500 DICOM files (I am working with CT data).

    • Subfolder: W002

      • Around 500 DICOM files (I am working with CT data).

    • Subfolder: ...

    • Subfolder: W100

      • Around 500 DICOM files (I am working with CT data).

enter image description here

I am using Slicer software and Python to automate the process of getting surface meshes automatically.

This is the Python script I wrote that needs to be run in each subfolder (please, don’t take into account the code itself)

As you can see, there are two lines with the location of the subfolder:

  • Second line, specifying the location of the subfolder
  • Almost the last line, indicating where to save the .ply file.

As I have around 1000 subfolders, I don’t want to run the script changing those lines. In other words, I wan to automate the process.

I am much more familiar with R than with Python. So:

  • could please indicate me how can I automate the process in detail by using Python?
  • could you please edit the code to adapt to my needs?
  • could you please show me a step-by-step process?
# Load DICOM files
dicomDataDir = "C:/Users/mario.modesto/Desktop/DICOM/W001"  # input folder with DICOM files
import os 
baboon_skull  = os.path.basename(dicomDataDir)
loadedNodeIDs = []  # this list will contain the list of all loaded node IDs

from DICOMLib import DICOMUtils
with DICOMUtils.TemporaryDICOMDatabase() as db:
  DICOMUtils.importDicom(dicomDataDir, db)
  patientUIDs = db.patients()
  for patientUID in patientUIDs:
    loadedNodeIDs.extend(DICOMUtils.loadPatientByUID(patientUID))

# Display volume rendering
logic = slicer.modules.volumerendering.logic()
volumeNode = slicer.mrmlScene.GetNodeByID('vtkMRMLScalarVolumeNode1')
displayNode = logic.CreateVolumeRenderingDisplayNode()
displayNode.UnRegister(logic)
slicer.mrmlScene.AddNode(displayNode)
volumeNode.AddAndObserveDisplayNodeID(displayNode.GetID())
logic.UpdateDisplayNodeFromVolumeNode(displayNode, volumeNode)

# find the files NodeID
volumeNode = getNode('2: Facial Bones  0.75  H70h')

#create a blank Markup ROI
roiNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsROINode")

#set the new markup ROI to the dimensions of the volume
cropVolumeParameters = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLCropVolumeParametersNode")
cropVolumeParameters.SetInputVolumeNodeID(volumeNode.GetID())
cropVolumeParameters.SetROINodeID(roiNode.GetID())
slicer.modules.cropvolume.logic().SnapROIToVoxelGrid(cropVolumeParameters)  # optional (rotates the ROI to match the volume axis directions)
slicer.modules.cropvolume.logic().FitROIToInputVolume(cropVolumeParameters)
slicer.mrmlScene.RemoveNode(cropVolumeParameters)

#set the cropping parameters
cropVolumeLogic = slicer.modules.cropvolume.logic()
cropVolumeParameterNode = slicer.vtkMRMLCropVolumeParametersNode()
cropVolumeParameterNode.SetIsotropicResampling(True)

#set the output resolution to 2 millimeters. units in slicer is always in mm. 
cropVolumeParameterNode.SetSpacingScalingConst(2)
cropVolumeParameterNode.SetROINodeID(roiNode.GetID())
cropVolumeParameterNode.SetInputVolumeNodeID(volumeNode.GetID())

#do the cropping
cropVolumeLogic.Apply(cropVolumeParameterNode)

#obtain the nodeID of the cropped volume 
croppedVolume = slicer.mrmlScene.GetNodeByID(cropVolumeParameterNode.GetOutputVolumeNodeID())

# Segmentation
segmentationNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLSegmentationNode")
segmentationNode.CreateDefaultDisplayNodes() # only needed for display
segmentationNode.SetReferenceImageGeometryParameterFromVolumeNode(croppedVolume)
addedSegmentID = segmentationNode.GetSegmentation().AddEmptySegment("skull")

# Create segment editor to get access to effects
segmentEditorWidget = slicer.qMRMLSegmentEditorWidget()
segmentEditorWidget.setMRMLScene(slicer.mrmlScene)
segmentEditorNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLSegmentEditorNode")
segmentEditorWidget.setMRMLSegmentEditorNode(segmentEditorNode)
segmentEditorWidget.setSegmentationNode(segmentationNode)
segmentEditorWidget.setMasterVolumeNode(croppedVolume)

# Thresholding
segmentEditorWidget.setActiveEffectByName("Threshold")
effect = segmentEditorWidget.activeEffect()
effect.setParameter("MinimumThreshold","115")
effect.setParameter("MaximumThreshold","3071")
effect.self().onApply()

# Clean up
segmentEditorWidget = None
slicer.mrmlScene.RemoveNode(segmentEditorNode)

# Make segmentation results visible in 3D
segmentationNode.CreateClosedSurfaceRepresentation()

# Creatint surface mesh and saving

surfaceMesh = segmentationNode.GetClosedSurfaceInternalRepresentation(addedSegmentID)
writer = vtk.vtkPLYWriter()
writer.SetInputData(surfaceMesh)
writer.SetFileName("C:/Users/mario.modesto/Desktop/DICOM/"+baboon_skull+"_surfaceMesh.ply")
writer.Update()

I would really appreciate your help.

UPDATE

I tried this code:

from DICOMLib import DICOMUtils
import os
from pathlib import Path
import shutil

for root, dirnames, _ in os.walk("C:/Users/mario.modesto/Desktop/DICOM/"):
    root_path = Path(root)
    for dirname in dirnames:
        dir_path = root_path / dirname
        dicomDataDir = dirname  # input folder with DICOM files
        baboon_skull  = os.path.basename(dicomDataDir)
        loadedNodeIDs = []  # this list will contain the list of all loaded node IDs
        
        with DICOMUtils.TemporaryDICOMDatabase() as db:
          DICOMUtils.importDicom(dicomDataDir, db)
          patientUIDs = db.patients()
          for patientUID in patientUIDs:
            loadedNodeIDs.extend(DICOMUtils.loadPatientByUID(patientUID))
        
        # Display volume rendering
        # https://slicer.readthedocs.io/en/latest/developer_guide/script_repository.html#display-volume-using-volume-rendering
        logic = slicer.modules.volumerendering.logic()
        volumeNode = slicer.mrmlScene.GetNodeByID('vtkMRMLScalarVolumeNode1')
        displayNode = logic.CreateVolumeRenderingDisplayNode()
        displayNode.UnRegister(logic)
        slicer.mrmlScene.AddNode(displayNode)
        volumeNode.AddAndObserveDisplayNodeID(displayNode.GetID())
        logic.UpdateDisplayNodeFromVolumeNode(displayNode, volumeNode)

I inserted some code to load DICOM files and to show volume rendering in the for loop. However, now I see this error:

Traceback (most recent call last):
  File "<string>", line 27, in <module>
AttributeError: 'NoneType' object has no attribute 'AddAndObserveDisplayNodeID'

However, when I tried the original code with only one folder, it works. Any idea?
2

I solved this by using this code (answered in SO):

import os
from DICOMLib import DICOMUtils

yourpath = r"<enter your path>/W001-W100"

#walk through DICOM directory
for dir in os.scandir(yourpath):
    # Load DICOM files
    dicomDataDir = dir.path  # path to input folder with DICOM files
    baboon_skull  = dir.name
    loadedNodeIDs = []  # this list will contain the list of all loaded node IDs

    with DICOMUtils.TemporaryDICOMDatabase() as db:
        DICOMUtils.importDicom(dicomDataDir, db)
        patientUIDs = db.patients()
        for patientUID in patientUIDs:
            loadedNodeIDs.extend(DICOMUtils.loadPatientByUID(patientUID))


    < rest of your code >

    writer.SetFileName(fr"{dicomDataDir}_surfaceMesh.ply")
    writer.Update()
3

Looking great.

By the way, if you do not need the metada contents of DICOM files, you can convert the DICOMs to NRRD outside of Slicer using a tool like DCM2NIIX (there is even a slicer extension). Often clinical research needs those DICOM fields for their project, but I don’t think that’s probably the case for those baboon scans. That may simplify your code and workflow.

But regardless you are almost there I think.