BUG: Force loading images from DICOM using right-handed coordinate sy…

…stem

Automatically flip left handed volumes when loading into 3D Slicer from DICOM to make them compatible with all algorithms.
Fixes flipped normals in segmentation node for left handed volumes.

Renamed FlipIJKCoordinateSystemHandedness method to ReverseSliceOrder to make it explicit that flipping of handedness is achieved
by reversing slice order. The method was introduced recently (no Slicer Stable Release has been created since then), therefore
the old name was not deprecated but removed.

See related commit that applies this regularization when the image is read using volume storage node:
Slicer@3802d81

See discussion at https://discourse.slicer.org/t/dentalsegmentator-result-mesh-is-inverted/37408/9

Libs/MRML/Core/vtkMRMLVolumeNode.cxx

@@ -1362,11 +1362,11 @@ bool vtkMRMLVolumeNode::IsIJKCoordinateSystemRightHanded(vtkMatrix4x4* ijkToRasM
 }
 
 //----------------------------------------------------------------------------
-void vtkMRMLVolumeNode::FlipIJKCoordinateSystemHandedness(vtkImageData* imageData, vtkMatrix4x4* ijkToRasMatrix)
+void vtkMRMLVolumeNode::ReverseSliceOrder(vtkImageData* imageData, vtkMatrix4x4* ijkToRasMatrix)
 {
   if (ijkToRasMatrix == nullptr)
   {
-    vtkGenericWarningMacro("vtkMRMLVolumeNode::FlipIJKCoordinateSystemHandedness failed: ijkToRasMatrix is invalid");
+    vtkGenericWarningMacro("vtkMRMLVolumeNode::ReverseSliceOrder failed: ijkToRasMatrix is invalid");
     return;
   }
   int imageDimensions[3] = { 0, 0, 0 };
@@ -1399,7 +1399,7 @@ void vtkMRMLVolumeNode::FlipIJKCoordinateSystemHandedness(vtkImageData* imageDat
       else
       {
         // There has been other data arrays, log a warning because we only flip the first one
-        vtkGenericWarningMacro("vtkMRMLVolumeNode::FlipIJKCoordinateSystemHandedness: Multiple types of point data arrays were found,"
+        vtkGenericWarningMacro("vtkMRMLVolumeNode::ReverseSliceOrder: Multiple types of point data arrays were found,"
           " only the " << vtkDataSetAttributes::GetAttributeTypeAsString(pointDataType) << " array will be flipped");
       }
     }
@@ -1453,6 +1453,6 @@ void vtkMRMLVolumeNode::SetIJKCoordinateSystemToRightHanded()
     // already right-handed
     return;
   }
-  vtkMRMLVolumeNode::FlipIJKCoordinateSystemHandedness(this->GetImageData(), ijkToRAS);
+  vtkMRMLVolumeNode::ReverseSliceOrder(this->GetImageData(), ijkToRAS);
   this->SetIJKToRASMatrix(ijkToRAS);
 }

Libs/MRML/Core/vtkMRMLVolumeNode.h

@@ -262,7 +262,7 @@ class VTK_MRML_EXPORT vtkMRMLVolumeNode : public vtkMRMLDisplayableNode
 
   /// Switch the IJK coordinate system handedness between left-handed and right-handed
   /// by inverting the K axis direction. The physical location of voxels do not change.
-  static void FlipIJKCoordinateSystemHandedness(vtkImageData* imageData, vtkMatrix4x4* ijkToRasMatrix);
+  static void ReverseSliceOrder(vtkImageData* imageData, vtkMatrix4x4* ijkToRasMatrix);
 
 protected:
   vtkMRMLVolumeNode();

Modules/Scripted/DICOMPlugins/DICOMScalarVolumePlugin.py

@@ -483,12 +483,26 @@ def setVolumeNodeProperties(self, volumeNode, loadable):
             #
             self.addSeriesInSubjectHierarchy(loadable, volumeNode)
 
+            # File for each slice. Will be reversed if reading a left-handed volume.
+            fileList = loadable.files
+
+            ijkToRAS = vtk.vtkMatrix4x4()
+            volumeNode.GetIJKToRASMatrix(ijkToRAS)
+            if not slicer.vtkMRMLVolumeNode.IsIJKCoordinateSystemRightHanded(ijkToRAS):
+                # This volume is in left-handed coordinate system. Reorder the slices to make it right-handed
+                # to ensure all computational algorithms work well (e.g., avoid turning generated surfaces inside out).
+                logging.info(f"Reverse slice order to force volume '{volumeNode.GetName()}' ({volumeNode.GetID()}) to use right-handed IJK coordinate system")
+                slicer.vtkMRMLVolumeNode.ReverseSliceOrder(volumeNode.GetImageData(), ijkToRAS)
+                volumeNode.SetIJKToRASMatrix(ijkToRAS)
+                # Reorder the file list so that i-th slice corresponds to i-th DICOM.instanceUIDs attribute value
+                fileList = list(reversed(fileList))
+
             #
             # add list of DICOM instance UIDs to the volume node
             # corresponding to the loaded files
             #
             instanceUIDs = ""
-            for file in loadable.files:
+            for file in fileList:
                 uid = slicer.dicomDatabase.fileValue(file, self.tags["instanceUID"])
                 if uid == "":
                     uid = "Unknown"
@@ -499,7 +513,7 @@ def setVolumeNodeProperties(self, volumeNode, loadable):
             # Choose a file in the middle of the series as representative frame,
             # because that is more likely to contain the object of interest than the first or last frame.
             # This is important for example for getting a relevant window/center value for the series.
-            file = loadable.files[int(len(loadable.files) / 2)]
+            file = fileList[int(len(fileList) / 2)]
 
             #
             # automatically select the volume to display