ENH: Improve vtkMRMLVolumeArchetypeStorageNode for left handed volumes

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

Libs/MRML/Core/Testing/vtkMRMLVolumeArchetypeStorageNodeTest1.cxx

@@ -16,6 +16,8 @@
 #include "vtkMRMLVolumeArchetypeStorageNode.h"
 
 #include "vtkImageData.h"
+#include "vtkMatrix3x3.h"
+#include "vtkMatrix4x4.h"
 #include "vtkNew.h"
 #include "vtkPointData.h"
 #include <vtksys/SystemTools.hxx>
@@ -124,6 +126,67 @@ int TestVoxelVectorType(const std::string& tempDir, const std::string& fileExten
   return EXIT_SUCCESS;
 }
 
+int TestFlipsLeftHandedVolumes(const std::string& tempDir)
+{
+  std::cout << "TestFlipsLeftHandedVolumes" << std::endl;
+  // Create a flipped image data with arbitrary K values
+  vtkNew<vtkMRMLScene> scene;
+
+  auto volumeNode = vtkMRMLScalarVolumeNode::SafeDownCast(scene->AddNewNodeByClass("vtkMRMLScalarVolumeNode"));
+  CHECK_NOT_NULL(volumeNode);
+
+  int n_i{ 10 }, n_j{ 20 }, n_k{ 30 };
+  vtkNew<vtkImageData> imageData;
+  imageData->SetDimensions(n_i, n_j, n_k);
+  imageData->AllocateScalars(VTK_FLOAT, 1);
+  imageData->GetPointData()->GetScalars()->Fill(0);
+
+  for (int i = 0; i < n_i; i++)
+  {
+    for (int j = 0; j < n_j; j++)
+    {
+      for (int k = 0; k < n_k; k++)
+      {
+        imageData->SetScalarComponentFromDouble(i, j, k, 0, k);
+      }
+    }
+  }
+  volumeNode->SetIJKToRASDirections(-1., 0., 0., 0., -1., 0., 0., 0., -1.);
+
+  const auto fileName = tempFilename(tempDir, "flipped_volume", "mha", true);
+  auto storageNode =
+    vtkMRMLVolumeArchetypeStorageNode::SafeDownCast(scene->AddNewNodeByClass("vtkMRMLVolumeArchetypeStorageNode"));
+  CHECK_NOT_NULL(storageNode);
+  storageNode->SetSingleFile(true);
+  volumeNode->SetAndObserveImageData(imageData);
+  volumeNode->SetAndObserveStorageNodeID(storageNode->GetID());
+
+  // Check that when loading, the K axis is flipped
+  storageNode->SetFileName(fileName.c_str());
+  CHECK_BOOL(storageNode->WriteData(volumeNode), true);
+  CHECK_BOOL(storageNode->ReadData(volumeNode), true);
+
+  vtkNew<vtkMatrix4x4> matrix;
+  volumeNode->GetIJKToRASDirectionMatrix(matrix);
+
+  CHECK_DOUBLE(matrix->GetElement(0, 0), -1.);
+  CHECK_DOUBLE(matrix->GetElement(1, 1), -1.);
+  CHECK_DOUBLE(matrix->GetElement(2, 2), 1.);
+
+  for (int i = 0; i < n_i; i++)
+  {
+    for (int j = 0; j < n_j; j++)
+    {
+      for (int k = 0; k < n_k; k++)
+      {
+        CHECK_DOUBLE(volumeNode->GetImageData()->GetScalarComponentAsDouble(i, j, k, 0), n_k - k - 1.);
+      }
+    }
+  }
+
+  return EXIT_SUCCESS;
+}
+
 int vtkMRMLVolumeArchetypeStorageNodeTest1(int argc, char* argv[])
 {
   if (argc != 2)
@@ -143,6 +206,7 @@ int vtkMRMLVolumeArchetypeStorageNodeTest1(int argc, char* argv[])
   CHECK_EXIT_SUCCESS(TestVoxelVectorType(tempDir, "png",  false,     false,   true,  true));
   CHECK_EXIT_SUCCESS(TestVoxelVectorType(tempDir, "tif",  false,     false,   true,  false));
   CHECK_EXIT_SUCCESS(TestVoxelVectorType(tempDir, "jpg",  false,     false,   true,  false));
+  CHECK_EXIT_SUCCESS(TestFlipsLeftHandedVolumes(tempDir));
 
   std::cout << "Test passed." << std::endl;
   return EXIT_SUCCESS;

Libs/MRML/Core/vtkMRMLVolumeArchetypeStorageNode.cxx

@@ -39,6 +39,7 @@ Version:   $Revision: 1.6 $
 #include <vtkDataArray.h>
 #include <vtkErrorCode.h>
 #include <vtkImageChangeInformation.h>
+#include <vtkImageFlip.h>
 #include <vtkMatrix3x3.h>
 #include <vtkNew.h>
 #include <vtkPointData.h>
@@ -281,6 +282,32 @@ void ApplyImageSeriesReaderWorkaround(vtkMRMLVolumeArchetypeStorageNode * storag
     }
   }
 }
+
+//----------------------------------------------------------------------------
+bool IsIJKCoordinateSystemLeftHanded(vtkMatrix4x4* rasToIjkMatrix)
+{
+  // Check if the determinant of the orientation matrix is less than 0.
+  vtkNew<vtkMatrix3x3> orientation;
+  vtkAddonMathUtilities::GetOrientationMatrix(rasToIjkMatrix, orientation);
+  return orientation->Determinant() < 0.;
+}
+
+//----------------------------------------------------------------------------
+void FlipIJKCoordinateSystemHandedness(vtkImageData* imageData, vtkMatrix4x4* rasToIjkMatrix)
+{
+  // Flip K Axis
+  vtkNew<vtkImageFlip> flip;
+  flip->SetFilteredAxes(2);
+  flip->SetInputData(imageData);
+  flip->Update();
+  imageData->ShallowCopy(flip->GetOutput());
+
+  // Flip K Direction
+  for (int i = 0; i < 3; i++)
+  {
+    rasToIjkMatrix->SetElement(i, 2, -rasToIjkMatrix->GetElement(i, 2));
+  }
+}
 } // end of anonymous namespace
 
 //----------------------------------------------------------------------------
@@ -512,13 +539,21 @@ int vtkMRMLVolumeArchetypeStorageNode::ReadDataInternal(vtkMRMLNode *refNode)
     << ". Number of components: " << outputImage->GetNumberOfScalarComponents()
     << ". Pixel type: " << vtkImageScalarTypeNameMacro(outputImage->GetScalarType()) << ".");
 
-  vtkMatrix4x4* mat = reader->GetRasToIjkMatrix();
-  if ( mat == nullptr )
+  vtkMatrix4x4* rasToIjkMatrix = reader->GetRasToIjkMatrix();
+  if (rasToIjkMatrix == nullptr)
   {
     vtkErrorToMessageCollectionMacro(this->GetUserMessages(), "vtkMRMLVolumeArchetypeStorageNode::ReadDataInternal",
       "Reader returned nullptr RasToIjkMatrix");
   }
-  volNode->SetRASToIJKMatrix(mat);
+
+  // If volume is left-handed coordinates, modify it to right-handed coordinate
+  // to have support for every algorithms in 3D Slicer
+  if (IsIJKCoordinateSystemLeftHanded(rasToIjkMatrix))
+  {
+    FlipIJKCoordinateSystemHandedness(outputImage, rasToIjkMatrix);
+  }
+
+  volNode->SetRASToIJKMatrix(rasToIjkMatrix);
 
   if (volNode->IsA("vtkMRMLDiffusionTensorVolumeNode"))
   {