Repository: InsightSoftwareConsortium/ITK (language: python)

Issue:
Support specification of Nifti sform tolerance via environmental variable
Per discussion with @effigies , it would be helpful to enable specification of a higher nifti sform tolerance via environmental variable.

There are cases where nifti headers, perhaps edited manually, that may have unintentional small shear, it would be helpful to specify an override via environmental variable a tolerance for orthogonality.

This would be used in docker images that use FSL, ANTs, etc.

Possibly related?

- https://github.com/InsightSoftwareConsortium/ITK/blob/fb61d2cadc33bf50d1f0ec5b842b70a1a7021d0e/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx#L1989-L2010

Todo:

- [ ] find a test case that demonstrates the issue
- [ ] create a patch (if needed) to enable specifying the related tolerance via environmental variable
- [ ] look into updated associated warnings to suggest setting the environmental variable if appropriate

CC @hjmjohnson 


Discussion / hints:
Just adding a note that I'm having trouble finding records of when we were hitting this. My memory was that there were small shears causing ITK to reject images, but it could also be differences between the sform column norms and `pixdim[1:4]`, which I believe is what this particular code is handling.

Will keep my eyes out for the next one. When we've encountered it, it's generally been due to a failure to validate and update headers before passing to ANTs, so it's possible we won't trigger the case again until the next refactor leaves a hole for a messy image to slip through.
I think this is already available via #4009 

https://github.com/InsightSoftwareConsortium/ITK/blob/fb61d2cadc33bf50d1f0ec5b842b70a1a7021d0e/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx#L2080-L2123

This works for me:

```
export ITK_NIFTI_SFORM_PERMISSIVE=1
PrintHeader bad_sform.nii
```
BTW, the expected values of `ITK_NIFTI_SFORM_PERMISSIVE` are true/false, on/off, yes/no. It cannot be toggled with 1/0. Is this ITK convention, or something to fix?
Sorry for not reading carefully before, I see the idea here is to add a numeric tolerance rather than the boolean we currently have - this seems like a good idea.

I had a go at creating a test image. However, I noticed that the permissive mode introduces a flip. 

It's very possible I'm doing something wrong here, so here's R code to create the dummy images, starting from a template

```R
library(RNifti)

# the nifti class contains pointers to memory, so read them
# independently to be safe
mni <- readNifti('tpl-MNI152NLin2009cAsym_res-02_T1w.nii.gz')
mni_rotated <- readNifti('tpl-MNI152NLin2009cAsym_res-02_T1w.nii.gz')
mni_sheared <- readNifti('tpl-MNI152NLin2009cAsym_res-02_T1w.nii.gz')

qoffset_x <- mni$qoffset_x
qoffset_y <- mni$qoffset_y
qoffset_z <- mni$qoffset_z

# Create the translation vector using qoffset values
translation_vector <- c(qoffset_x, qoffset_y, qoffset_z)

spacing <- diag(c(2, 2, 2))

# blank out qform
qform(mni) <- structure(diag(4), code=0L)
qform(mni_rotated) <- structure(diag(4), code=0L)
qform(mni_sheared) <- structure(diag(4), code=0L)

writeNifti(mni, 'mni_sform.nii.gz')

rotation <- matrix(c(
  0.99765903, -0.06763516,  0.01009717,
  0.06808255,  0.99622887, -0.05378478,
  -0.00642135,  0.05434632,  0.9985015
), nrow = 3, byrow = TRUE)

shear_matrix <- matrix(c(
  1, 0.0001, 0.0001,
  0.0, 1, 0.0,
  0.0, 0.0, 1
), nrow = 3, byrow = TRUE)

almost_rotation <- shear_matrix %*% rotation

sform_rotated <- rbind(cbind(rotation %*% spacing, translation_vector), c(0, 0, 0, 1))

sform(mni_rotated) <- structure(sform_rotated, code=1L)

writeNifti(mni_rotated, 'mni_rotated.nii.gz')

sform_sheared <- rbind(cbind(almost_rotation %*% spacing, translation_vector), c(0, 0, 0, 1))

sform(mni_sheared) <- structure(sform_sheared, code=1L)

writeNifti(mni_sheared, 'mni_sheared.nii.gz')
```

The shear is small but enough to make non-permissive mode fail. ITK-SNAP 4.2.0 refuses to read `mni_sheared.nii.gz`, and so does ANTs PrintHeader:

```bash
PrintHeader mni_sheared.nii.gz 4
 cant read mni_sheared.nii.gz
```

```bash
ITK_NIFTI_SFORM_PERMISSIVE=TRUE PrintHeader mni_sheared.nii.gz 4

WARNING: In /Users/pcook/tmp/NOT_BACKED_UP/antsMasterBuild/build/ITKv5/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx, line 2111
NiftiImageIO (0x7fccec328ce0): mni_sheared.nii.gz has non-orthogonal sform

WARNING: In /Users/pcook/tmp/NOT_BACKED_UP/antsMasterBuild/build/ITKv5/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx, line 2111
NiftiImageIO (0x7fccec328440): mni_sheared.nii.gz has non-orthogonal sform

-0.997662x-0.0680327x0.00647123x0.0675826x-0.996232x-0.0543497x0.0101444x-0.0537853x0.998501
```

To try to visualize what happens in permissive mode, I used antsApplyTransforms. This should read the image using the permissive mode code, and then resample it into the space of mni_rotated, which does not trigger the permissive mode sform formulation

```bash
ITK_NIFTI_SFORM_PERMISSIVE=TRUE antsApplyTransforms \
  -d 3 -i mni_sheared.nii.gz -r mni_rotated.nii.gz -o mni_shear_permissive.nii.gz  \
  -t Identity --verbose --float
```

But loading all these into ITK-SNAP, `mni_shear_permissive.nii.gz` appears rotated. I figured this was because it was going back to the qform, so I removed that

```R
library(RNifti)
permissive <- readNifti('mni_shear_permissive.nii.gz')
qform(permissive) <- structure(diag(4), code=0L)
c<- structure(diag(4), code=0L)
writeNifti(permissive, 'mni_shear_permissive_sform.nii.gz')
```

So I think the qform is not the issue, the sform created in permissive mode seems to be flipped somehow. Perhaps some sign indeterminancy in SVD? 

```
itksnap -g mni_sform.nii.gz -o mni_rotated.nii.gz mni_shear_permissive_sform.nii.gz 
```

cc: @vfonov 

Fix patch:
diff --git a/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx b/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx
index c018f88b90f..7fde5622799 100644
--- a/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx
+++ b/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)


Test patch:
diff --git a/Modules/IO/NIFTI/test/CMakeLists.txt b/Modules/IO/NIFTI/test/CMakeLists.txt
index de9efbf1f3e..07ce2583806 100644
--- a/Modules/IO/NIFTI/test/CMakeLists.txt
+++ b/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
diff --git a/Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx b/Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx
index 7c7d1f377de..081bb72407c 100644
--- a/Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx
+++ b/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;
 }