ENH: Add CLI module for DID calculations

AutoscoperM/AutoscoperMLib/RadiographGeneration.py

@@ -318,12 +318,34 @@ def optimizeCameras(
         def progressCallback(_x):
             return None
 
+    # Parallel calls to cliModule
+    cliModule = slicer.modules.calculatedataintensitydensity
+    cliNodes = []
     for i in range(len(cameras)):
         camera = cameras[i]
         vrgFName = glob.glob(os.path.join(cameraDir, f"cam{camera.id}", "*.tif"))[0]
-        _calculateDataIntensityDensity(camera, vrgFName)
+        # _calculateDataIntensityDensity(camera, vrgFName)
+        cliNode = slicer.cli.run(
+            cliModule,
+            None,
+            {"whiteRadiographFName": vrgFName},
+            wait_for_completion=False,
+        )
+        cliNodes.append(cliNode)
+
+    for i in range(len(cameras)):
+        camera = cameras[i]
+        while cliNodes[i].IsBusy():
+            slicer.app.processEvents()
+        if cliNode.GetStatus() & cliNode.ErrorsMask:
+            # error
+            errorText = cliNode.GetErrorText()
+            slicer.mrmlScene.RemoveNode(cliNode)
+            raise ValueError("CLI execution failed: " + errorText)
+        camera.DID = float(cliNodes[i].GetOutputText())
         progress = ((i + 1) / len(cameras)) * 50 + 40
         progressCallback(progress)
+        slicer.mrmlScene.RemoveNode(cliNodes[i])
 
     cameras.sort(key=lambda x: x.DID)
 

CMakeLists.txt

@@ -37,6 +37,7 @@ endif()
 # Extension modules
 add_subdirectory(AutoscoperM)
 add_subdirectory(TrackingEvaluation)
+add_subdirectory(calculateDataIntensityDensity)
 ## NEXT_MODULE
 
 #-----------------------------------------------------------------------------

calculateDataIntensityDensity/CMakeLists.txt

@@ -0,0 +1,6 @@
+#-----------------------------------------------------------------------------
+set(MODULE_NAME calculateDataIntensityDensity)
+
+SlicerMacroBuildScriptedCLI(
+  NAME ${MODULE_NAME}
+  )

calculateDataIntensityDensity/calculateDataIntensityDensity.py

@@ -0,0 +1,69 @@
+import sys
+
+
+def calculateDataIntensityDensity(whiteRadiographFName: str) -> float:
+    """
+    Calculates the data intensity density of the given camera on its corresponding white radiograph.
+    Internal function used by :func:`optimizeCameras`.
+
+    :param whiteRadiographFName: White radiograph file name
+    :type whiteRadiographFName: str
+
+    return Data intensity density
+    :rtype: float
+    """
+
+    import numpy as np
+    import SimpleITK as sitk
+
+    MEAN_COMPARISON = 170  # 255 / 3 * 2
+
+    # Read in the white radiograph
+    whiteRadiograph = sitk.ReadImage(whiteRadiographFName)
+
+    # Superpixel Segmentation
+    slicImageFilter = sitk.SLICImageFilter()
+    slicImageFilter.SetSuperGridSize([15, 15, 15])  # smaller grid size = finer grid overall default is [50,50,50]
+    labelImage = slicImageFilter.Execute(whiteRadiograph)
+
+    # Get the mean pixel value for each label
+    labelStatsFilter = sitk.LabelStatisticsImageFilter()
+    labelStatsFilter.Execute(whiteRadiograph, labelImage)
+    N = labelStatsFilter.GetNumberOfLabels()
+    meanColor = np.zeros((N, 1))
+    m, n = labelImage.GetSize()
+    labels = list(labelStatsFilter.GetLabels())
+    labels.sort()
+    for i, label in enumerate(labels):
+        meanColor[i, 0] = labelStatsFilter.GetMean(label)
+
+    # Create a binary label from the labelImage where all '1' are labels whose meanColor are < 255/3
+    labelShapeFilter = sitk.LabelShapeStatisticsImageFilter()
+    labelShapeFilter.Execute(labelImage)
+    binaryLabels = np.zeros((m, n))
+    for i, label in enumerate(labels):
+        if label == 0:
+            continue
+        if meanColor[i, 0] < MEAN_COMPARISON:
+            pixels = list(labelShapeFilter.GetIndexes(label))
+            for j in range(0, len(pixels), 2):
+                y = pixels[j]
+                x = pixels[j + 1]
+                binaryLabels[x, y] = 1
+
+    # Calculate the Data Intensity Density
+    # Largest Region based off of https://discourse.itk.org/t/simpleitk-extract-largest-connected-component-from-binary-image/4958/2
+    binaryImage = sitk.Cast(sitk.GetImageFromArray(binaryLabels), sitk.sitkUInt8)
+    componentImage = sitk.ConnectedComponent(binaryImage)
+    sortedComponentImage = sitk.RelabelComponent(componentImage, sortByObjectSize=True)
+    largest = sortedComponentImage == 1
+
+    return np.sum(sitk.GetArrayFromImage(largest))
+
+
+if __name__ == "__main__":
+    EXPECTED_ARGS = 2
+    if len(sys.argv) < EXPECTED_ARGS:
+        print("Usage: calculateDataIntensityDensity <input FName>")
+        sys.exit(1)
+    print(calculateDataIntensityDensity(sys.argv[1]))

calculateDataIntensityDensity/calculateDataIntensityDensity.xml

@@ -0,0 +1,23 @@
+<?xml version="1.0" encoding="utf-8"?>
+<executable>
+  <category>Examples</category>
+  <index>0</index>
+  <title>calculateDataIntensityDensity</title>
+  <description><![CDATA[Apply a Gaussian blur to an image]]></description>
+  <version>0.1.0.</version>
+  <documentation-url>https://github.com/username/project</documentation-url>
+  <license/>
+  <contributor>Andras Lasso (PerkLab)</contributor>
+  <acknowledgements><![CDATA[This work is part of the National Alliance for Medical Image Computing (NAMIC), funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 EB005149.]]></acknowledgements>
+  <parameters>
+    <label>IO</label>
+    <description><![CDATA[Input/output parameters]]></description>
+    <string>
+      <name>whiteRadiographFName</name>
+      <label>White Radiograph File Name</label>
+      <index>1</index>
+      <description><![CDATA[White radiograph file name]]></description>
+      <channel>input</channel>
+    </string>
+  </parameters>
+</executable>