ENH: Add DICOM grouping frames by orientation and harden transform op…

…tion.

Co-authored-by: Jean-Christophe Fillion-Robin <jchris.fillionr@kitware.com>

Modules/Scripted/DICOMPlugins/DICOMScalarVolumePlugin.py

@@ -22,6 +22,42 @@
 # plugin architecture.
 #
 
+class TagValueList(list):
+    def __init__(self, tag, epsilon):
+        self.tag = tag
+        self.epsilon = epsilon
+        # Duplicate values for these tags will compared numerically,
+        # instead of checking the values of the dictionary subseriesValues stored as strings.
+        self.subseriesTagsToVectorValues = [
+            "imageOrientationPatient",
+            "diffusionGradientOrientation",
+        ]
+
+    def tagValueToVector(self, value):
+        if value == "":
+            return numpy.array(0)
+        return numpy.array([float(element) for element in value.split("\\")])
+
+    def __contains__(self, value):
+        # The value is a vector, use custom numeric comparison
+        if self.tag in self.subseriesTagsToVectorValues:
+            valueVector = self.tagValueToVector(value)
+            for subseriesValue in self:
+                subseriesVector = self.tagValueToVector(subseriesValue)
+                # The 'eps' value default is set to 1.e-6, which reproduces the same behavior as
+                # the ITK numeric comparison logic (i.e. absolute difference).
+                # Reference:
+                #   Class: ITK/Modules/Numerics/Optimizersv4/include/itkObjectToObjectMetric.hxx
+                #   Method: VerifyDisplacementFieldSizeAndPhysicalSpace
+                #   URL: https://github.com/hinerm/ITK/blob/master/Modules/Numerics/Optimizersv4/include/itkObjectToObjectMetric.hxx#L471-L475.
+                if numpy.allclose(valueVector, subseriesVector, rtol=0.0, atol=self.epsilon):
+                    return True, subseriesValue
+            return False, value
+
+        # The value is a string, use default list __contains__ method
+        value = value.replace(",", "_")  # remove commas so it can be used as an index
+        return super().__contains__(value), value
+
 class DICOMScalarVolumePluginClass(DICOMPlugin):
     """ ScalarVolume specific interpretation code
     """
@@ -90,7 +126,7 @@ def settingsPanelEntry(panel, parent):
         importFormatsComboBox.addItem(_("default (apply regularization transform)"), "default")
         importFormatsComboBox.addItem(_("none"), "none")
         importFormatsComboBox.addItem(_("apply regularization transform"), "transform")
-        # In the future additional option, such as "resample" (harden the applied transform) may be added.
+        importFormatsComboBox.addItem(_("harden regularization transform"), "hardenTransform")
 
         importFormatsComboBox.currentIndex = 0
         formLayout.addRow(_("Acquisition geometry regularization:"), importFormatsComboBox)
@@ -134,6 +170,10 @@ def compareVolumeNodes(volumeNode1, volumeNode2):
             comparison += _("Pixel data mismatch") + "\n"
         return comparison
 
+    def hardenAcquisitionGeometryRegularization(self):
+        settings = qt.QSettings()
+        return (settings.value("DICOM/ScalarVolume/AcquisitionGeometryRegularization", "default") == "hardenTransform")
+
     def acquisitionGeometryRegularizationEnabled(self):
         settings = qt.QSettings()
         return (settings.value("DICOM/ScalarVolume/AcquisitionGeometryRegularization", "default") != "none")
@@ -171,10 +211,12 @@ def cleanNodeName(self, value):
         cleanValue = cleanValue.replace("\\", "-")
         return cleanValue
 
-    def examineFiles(self, files):
+    def examineFiles(self, files, eps=1.e-6):
         """ Returns a list of DICOMLoadable instances
         corresponding to ways of interpreting the
         files parameter.
+
+        eps value default to 1.e-6 as same behaviour of ITK numeric comparison logic.
         """
 
         seriesUID = slicer.dicomDatabase.fileValue(files[0], self.tags['seriesUID'])
@@ -223,10 +265,10 @@ def examineFiles(self, files):
             # check for subseries values
             for tag in subseriesTags:
                 value = slicer.dicomDatabase.fileValue(file, self.tags[tag])
-                value = value.replace(",", "_")  # remove commas so it can be used as an index
                 if tag not in subseriesValues:
-                    subseriesValues[tag] = []
-                if not subseriesValues[tag].__contains__(value):
+                    subseriesValues[tag] = TagValueList(tag=tag, epsilon=eps)
+                contained, value = subseriesValues[tag].__contains__(value)
+                if not contained:
                     subseriesValues[tag].append(value)
                 if (tag, value) not in subseriesFiles:
                     subseriesFiles[tag, value] = []
@@ -539,7 +581,8 @@ def load(self, loadable, readerApproach=None):
         if volumeNode:
             self.acquisitionModeling = self.AcquisitionModeling()
             self.acquisitionModeling.createAcquisitionTransform(volumeNode,
-                                                                addAcquisitionTransformIfNeeded=self.acquisitionGeometryRegularizationEnabled())
+                                                                addAcquisitionTransformIfNeeded=self.acquisitionGeometryRegularizationEnabled(),
+                                                                hardenAcquisitionTransform=self.hardenAcquisitionGeometryRegularization())
 
         return volumeNode
 
@@ -827,7 +870,7 @@ def cornersToWorld(self, volumeNode, corners):
                         volumeNode.TransformPointToWorld(corners[slice, row, column], worldCorners[slice, row, column])
             return worldCorners
 
-        def createAcquisitionTransform(self, volumeNode, addAcquisitionTransformIfNeeded=True):
+        def createAcquisitionTransform(self, volumeNode, addAcquisitionTransformIfNeeded=True, hardenAcquisitionTransform=False):
             """Creates the actual transform if needed.
             Slice corners are cached for inspection by tests
             """
@@ -846,6 +889,8 @@ def createAcquisitionTransform(self, volumeNode, addAcquisitionTransformIfNeeded
                     self.fixedCorners = self.cornersToWorld(volumeNode, self.originalCorners)
                     if not numpy.allclose(self.fixedCorners, self.targetCorners):
                         raise Exception("Acquisition transform didn't fix slice corners!")
+                    if hardenAcquisitionTransform:
+                        volumeNode.HardenTransform()
                 else:
                     logging.warning(warningText + "  Regularization transform is not added, as the option is disabled.")
             elif maxError > 0 and maxError > self.zeroEpsilon: