Merge pull request InsightSoftwareConsortium#4847 from cookpa/sform_p…

…ermissive_svd_fix

BUG: permissive cosine computation from non-orthogonal sform had reve…

Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

@@ -2097,7 +2097,7 @@ NiftiImageIO::SetImageIOOrientationFromNIfTI(unsigned short dims, double spacing
 
         if (svd.singularities() == 0)
         {
-          mat = svd.V() * svd.U().conjugate_transpose();
+          mat = svd.U() * svd.V().conjugate_transpose();
           mat44 _mat;
 
           for (int i = 0; i < 3; ++i)

Modules/IO/NIFTI/test/CMakeLists.txt

@@ -224,7 +224,10 @@ itk_add_test(
   DATA{${ITK_DATA_ROOT}/Input/LPSLabels_noqform.nii.gz}
   DATA{${ITK_DATA_ROOT}/Input/LPSLabels_nosform.nii.gz}
   DATA{${ITK_DATA_ROOT}/Input/NonOrthoSform.nii.gz}
-  ${ITK_TEST_OUTPUT_DIR})
+  ${ITK_TEST_OUTPUT_DIR}
+  8.0 # Tolerance in degrees for direction recovered from non-orthogonal sform
+  )
+
 
 itk_add_test(
   NAME
@@ -236,7 +239,9 @@ itk_add_test(
   DATA{Input/SmallVoxels_noqform.nii.gz}
   DATA{Input/SmallVoxels_nosform.nii.gz}
   DATA{Input/SmallVoxelsNonOrthoSform.nii.gz}
-  ${ITK_TEST_OUTPUT_DIR})
+  ${ITK_TEST_OUTPUT_DIR}
+  4.0 # Tolerance in degrees for direction recovered from non-orthogonal sform
+  )
 
 itk_add_test(
   NAME

Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx

@@ -17,6 +17,7 @@
  *=========================================================================*/
 #include <iostream>
 #include <fstream>
+#include "itkMath.h"
 #include "itkNiftiImageIO.h"
 #include "itkImageFileReader.h"
 #include "itkImageFileWriter.h"
@@ -79,8 +80,8 @@ itkNiftiReadWriteDirectionTest(int argc, char * argv[])
   {
     std::cerr << "Missing Parameters." << std::endl;
     std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv)
-              << "imageWithBothQAndSForms imageWithNoQform imageWithNoSform  imageWithNonOrthoSform testOutputDir"
-              << std::endl;
+              << "imageWithBothQAndSForms imageWithNoQform imageWithNoSform imageWithNonOrthoSform testOutputDir "
+              << " [permissive_sform_threshold_deg=8.0]" << std::endl;
     std::cerr << "6 arguments required, received " << argc << std::endl;
     for (int i = 0; i < argc; ++i)
     {
@@ -89,13 +90,21 @@ itkNiftiReadWriteDirectionTest(int argc, char * argv[])
     return EXIT_FAILURE;
   }
 
+  // Threshold for sform test - how much can the sform matrix recovered from the non-ortho image differ from the
+  // original
+  double sformPermissiveAngleThreshold = 8.0;
+  if (argc > 6)
+  {
+    sformPermissiveAngleThreshold = std::stod(argv[6]);
+  }
+
   using TestImageType = itk::Image<float, 3>;
   TestImageType::Pointer inputImage = itk::ReadImage<TestImageType>(argv[1]);
   TestImageType::Pointer inputImageNoQform = itk::ReadImage<TestImageType>(argv[2]);
   TestImageType::Pointer inputImageNoSform = itk::ReadImage<TestImageType>(argv[3]);
 
 
-  // check if rotation matrix is orthogonal
+  // Check if rotation matrix is orthogonal
   if (!CheckRotation<TestImageType>(inputImage))
   {
     std::cerr << "Rotation matrix of " << argv[1] << " is not orthgonal" << std::endl;
@@ -188,34 +197,65 @@ itkNiftiReadWriteDirectionTest(int argc, char * argv[])
     return EXIT_FAILURE;
   }
 
-  // should not read the image with non-orthogonal sform
+  // Should not read the image with non-orthogonal sform
   ITK_TRY_EXPECT_EXCEPTION(ReadImage<TestImageType>(argv[4], false));
 
-  // check if environment flag is used properly
-  // check without permissive option
+  // Check if environment flag is used properly
+  // Check without permissive option
   itksys::SystemTools::PutEnv("ITK_NIFTI_SFORM_PERMISSIVE=NO");
   ITK_TRY_EXPECT_EXCEPTION(itk::ReadImage<TestImageType>(argv[4]));
-  // check with permissive option
+  // Check with permissive option
   itksys::SystemTools::PutEnv("ITK_NIFTI_SFORM_PERMISSIVE=YES");
   ITK_TRY_EXPECT_NO_EXCEPTION(itk::ReadImage<TestImageType>(argv[4]));
 
-  // this should work
+  // This should work
   TestImageType::Pointer inputImageNonOrthoSform = ReadImage<TestImageType>(argv[4], true);
   dictionary = inputImageNonOrthoSform->GetMetaDataDictionary();
   if (!itk::ExposeMetaData<std::string>(dictionary, "ITK_sform_corrected", temp) || temp != "YES")
   {
     std::cerr << "ITK_sform_corrected metadata flag was not properly set" << std::endl;
-    std::cerr << " expected YES, recieved:" << temp.c_str() << std::endl;
+    std::cerr << " expected YES, received:" << temp.c_str() << std::endl;
     return EXIT_FAILURE;
   }
 
-  // check the resulting rotation matrix is orthogonal
+  itksys::SystemTools::PutEnv("ITK_NIFTI_SFORM_PERMISSIVE=NO");
+
+  // Check the resulting rotation matrix is orthogonal
   if (!CheckRotation<TestImageType>(inputImageNonOrthoSform))
   {
     std::cerr << "Rotation matrix after correcting is not orthgonal" << std::endl;
     return EXIT_FAILURE;
   }
 
+  // Check similarity between direction without shear (argv[2]) and direction with shear (argv[4])
+  // inputImageNoQformDirection should be somewhat similar to inputImageNonOrthoSformDirection.
+  // The shear in the test data is pretty large, so they are off by a bit, but check angle to catch
+  // rotation flips introduced by previous bugs
+  auto inputImageNonOrthoSformDirection = inputImageNonOrthoSform->GetDirection();
+
+  // Compare the direction matrices
+  auto nonOrthoDirectionInv = inputImageNonOrthoSformDirection.GetTranspose();
+
+  // Compute R_relative = R1 * R2^-1
+  auto rRelative = inputImageNoQformDirection * nonOrthoDirectionInv;
+
+  // Calculate the trace of the relative rotation matrix
+  double trace = rRelative[0][0] + rRelative[1][1] + rRelative[2][2];
+
+  // Calculate the angle of rotation between the two matrices
+  double angle = std::acos((trace - 1.0) / 2.0) * 180.0 / vnl_math::pi;
+
+  // The angle between the two matrices will depend on the amount of shear in the sform. Some test images
+  // have relatively large shear to make sure they trigger the sform correction. In practice, permissive mode
+  // should really only be used to account for small numerical errors in the sform matrix
+  if (angle > sformPermissiveAngleThreshold)
+  {
+    std::cerr << "Error: direction matrices are not similar" << std::endl;
+    std::cerr << "Rotation angle: " << angle << " degrees, exceeds threshold: " << sformPermissiveAngleThreshold
+              << std::endl;
+    return EXIT_FAILURE;
+  }
+
   std::cout << "Test finished." << std::endl;
   return EXIT_SUCCESS;
 }