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

…tion.

Modules/Scripted/DICOMPlugins/DICOMScalarVolumePlugin.py

@@ -90,7 +90,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 +134,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 +175,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'])
@@ -212,6 +218,13 @@ def examineFiles(self, files):
             "diffusionGradientOrientation",
         ]
 
+        # Duplicate values for these tags will compared numerically,
+        # instead of checking the values of the dictionary subseriesValues stored as strings.
+        vectorTags = [
+            "imageOrientationPatient",
+            "diffusionGradientOrientation",
+        ]
+
         #
         # first, look for subseries within this series
         # - build a list of files for each unique value
@@ -224,10 +237,34 @@ def examineFiles(self, files):
             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):
+
+                if tag in vectorTags:
+                    if value != "":
+                        vector = numpy.array([float(element) for element in value.split("\\")])
+
+                        found = False
+                        for subseriesValue in subseriesValues[tag]:
+                            subseriesVector = numpy.array([float(element) for element in subseriesValue.split("\\")])
+
+                            # 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(vector, subseriesVector, rtol = 0., atol = eps):
+                                found = True
+                                value = subseriesValue
+                                break
+
+                        if not found:
+                            subseriesValues[tag].append(value)
+                elif not subseriesValues[tag].__contains__(value):
                     subseriesValues[tag].append(value)
+
                 if (tag, value) not in subseriesFiles:
                     subseriesFiles[tag, value] = []
                 subseriesFiles[tag, value].append(file)
@@ -539,7 +576,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 +865,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 +884,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: