ENH: Grouping frames by orientation in DICOM module and add harden tr…

…ansform option in the DICOM settings.

Libs/vtkITK/vtkITKArchetypeImageSeriesReader.cxx

@@ -1290,7 +1290,6 @@ void vtkITKArchetypeImageSeriesReader::AnalyzeDicomHeaders()
   this->ImageOrientationPatient.resize( 0 );
   this->ImagePositionPatient.resize( 0 );
 
-
   itk::GDCMImageIO::Pointer gdcmIO = itk::GDCMImageIO::New();
   if ( !gdcmIO->CanReadFile(this->Archetype) )
     {

Modules/Scripted/DICOMPlugins/DICOMScalarVolumePlugin.py

@@ -6,6 +6,7 @@
 import qt
 import vtk
 import vtkITK
+import math
 from slicer.i18n import tr as _
 
 import slicer
@@ -90,6 +91,7 @@ def settingsPanelEntry(panel, parent):
         importFormatsComboBox.addItem(_("default (apply regularization transform)"), "default")
         importFormatsComboBox.addItem(_("none"), "none")
         importFormatsComboBox.addItem(_("apply regularization transform"), "transform")
+        importFormatsComboBox.addItem(_("apply harden regularization transform"), "hardenTransform")
         # In the future additional option, such as "resample" (harden the applied transform) may be added.
 
         importFormatsComboBox.currentIndex = 0
@@ -134,6 +136,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,7 +177,7 @@ def cleanNodeName(self, value):
         cleanValue = cleanValue.replace("\\", "-")
         return cleanValue
 
-    def examineFiles(self, files):
+    def examineFiles(self, files, eps = 1.e-4):
         """ Returns a list of DICOMLoadable instances
         corresponding to ways of interpreting the
         files parameter.
@@ -224,10 +230,40 @@ 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 == "imageOrientationPatient" or tag == "diffusionGradientOrientation":
+                    if value != "":
+                        directionCosine = [float(i) for i in value.split("\\")]
+                        a = math.sqrt(directionCosine[0]*directionCosine[0] + directionCosine[1]*directionCosine[1] + directionCosine[2]*directionCosine[2])
+                        for k in range(3):
+                            directionCosine[k] /= a
+                        a = math.sqrt(directionCosine[3]*directionCosine[3] + directionCosine[4]*directionCosine[4] + directionCosine[5]*directionCosine[5])
+                        for k in range(3, 6):
+                            directionCosine[k] /= a
+
+                        found = False
+                        for subseriesValue in subseriesValues[tag]:
+                            subserDirCos = [float(i) for i in subseriesValue.split("\\")]
+                            a = math.sqrt(subserDirCos[0]*subserDirCos[0] + subserDirCos[1]*subserDirCos[1] + subserDirCos[2]*subserDirCos[2])
+                            b = (directionCosine[0]*subserDirCos[0] + directionCosine[1]*subserDirCos[1] + directionCosine[2]*subserDirCos[2]) / a
+                            if b < 1 - eps :
+                                continue
+
+                            a = math.sqrt(subserDirCos[3]*subserDirCos[3] + subserDirCos[4]*subserDirCos[4] + subserDirCos[5]*subserDirCos[5])
+                            b = (directionCosine[3]*subserDirCos[3] + directionCosine[4]*subserDirCos[4] + directionCosine[5]*subserDirCos[5]) / a
+                            if  b > 1 - 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 +575,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 +864,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 +883,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: