ENH: Added permissive mode to NIFTI reader

CMakeLists.txt

@@ -338,6 +338,11 @@ if(NOT ("${ITK_NIFTI_IO_ANALYZE_FLAVOR}" STREQUAL "ITK4Warning" OR
 endif()
 mark_as_advanced(ITK_NIFTI_IO_ANALYZE_FLAVOR)
 
+# Allow non-orthogonal rotation matrix approximation in NIFTI files
+option(ITK_NIFTI_IO_SFORM_PERMISSIVE_DEFAULT "Allow non-orthogonal rotation matrix in NIFTI sform by default" OFF)
+mark_as_advanced(ITK_NIFTI_IO_SFORM_PERMISSIVE_DEFAULT)
+
+
 
 #-----------------------------------------------------------------------------
 # ITK build classes that are in the review process

Modules/IO/NIFTI/include/itkNiftiImageIO.h

@@ -219,6 +219,11 @@ class ITKIONIFTI_EXPORT NiftiImageIO : public ImageIOBase
   itkGetConstMacro(ConvertRASDisplacementVectors, bool);
   itkBooleanMacro(ConvertRASDisplacementVectors);
 
+  /** Allow to read nifti files with non-orthogonal sform*/
+  itkSetMacro(SFORM_Permissive, bool);
+  itkGetConstMacro(SFORM_Permissive, bool);
+  itkBooleanMacro(SFORM_Permissive);
+
 protected:
   NiftiImageIO();
   ~NiftiImageIO() override;
@@ -276,6 +281,9 @@ class ITKIONIFTI_EXPORT NiftiImageIO : public ImageIOBase
   IOComponentEnum m_OnDiskComponentType{ IOComponentEnum::UNKNOWNCOMPONENTTYPE };
 
   NiftiImageIOEnums::Analyze75Flavor m_LegacyAnalyze75Mode{};
+
+  bool m_SFORM_Permissive;
+  bool m_SFORM_Corrected{ false };
 };
 
 

Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

@@ -22,6 +22,8 @@
 #include <nifti1_io.h>
 #include "itkNiftiImageIOConfigurePrivate.h"
 #include "itkMakeUniqueForOverwrite.h"
+#include "itksys/SystemTools.hxx"
+#include "itksys/SystemInformation.hxx"
 
 namespace itk
 {
@@ -449,6 +451,7 @@ NiftiImageIO::NiftiImageIO()
   : m_NiftiImageHolder(new NiftiImageProxy(nullptr))
   , m_NiftiImage(*m_NiftiImageHolder.get())
   , m_LegacyAnalyze75Mode{ ITK_NIFTI_IO_ANALYZE_FLAVOR_DEFAULT }
+  , m_SFORM_Permissive{ ITK_NIFTI_IO_SFORM_PERMISSIVE_DEFAULT }
 {
   this->SetNumberOfDimensions(3);
   nifti_set_debug_level(0); // suppress error messages
@@ -460,6 +463,12 @@ NiftiImageIO::NiftiImageIO()
     this->AddSupportedWriteExtension(ext);
     this->AddSupportedReadExtension(ext);
   }
+  std::string envVar;
+  if (itksys::SystemTools::GetEnv("ITK_NIFTI_SFORM_PERMISSIVE", envVar))
+  {
+    envVar = itksys::SystemTools::UpperCase(envVar);
+    this->SetSFORM_Permissive(envVar != "NO" && envVar != "OFF" && envVar != "FALSE");
+  }
 }
 
 NiftiImageIO::~NiftiImageIO()
@@ -481,6 +490,7 @@ NiftiImageIO::PrintSelf(std::ostream & os, Indent indent) const
   os << indent << "ConvertRASDisplacementVectors: " << (m_ConvertRASDisplacementVectors ? "On" : "Off") << std::endl;
   os << indent << "OnDiskComponentType: " << m_OnDiskComponentType << std::endl;
   os << indent << "LegacyAnalyze75Mode: " << m_LegacyAnalyze75Mode << std::endl;
+  os << indent << "SFORM permissive: " << (m_SFORM_Permissive ? "On" : "Off") << std::endl;
 }
 
 bool
@@ -899,6 +909,8 @@ NiftiImageIO::SetImageIOMetadataFromNIfTI()
   // Necessary to clear dict if ImageIO object is re-used
   thisDic.Clear();
 
+  EncapsulateMetaData<std::string>(thisDic, "ITK_sform_corrected", m_SFORM_Corrected ? "YES" : "NO");
+
   std::ostringstream nifti_type;
   nifti_type << nim->nifti_type;
   EncapsulateMetaData<std::string>(thisDic, "nifti_type", nifti_type.str());
@@ -2031,6 +2043,7 @@ NiftiImageIO::SetImageIOOrientationFromNIfTI(unsigned short dims)
           // Only orthonormal matricies have transpose as inverse
           const vnl_matrix_fixed<float, 3, 3> candidate_identity = rotation * rotation.transpose();
           const bool                          is_orthonormal = candidate_identity.is_identity(1.0e-4);
+
           return is_orthonormal;
         }
       }();
@@ -2089,15 +2102,59 @@ NiftiImageIO::SetImageIOOrientationFromNIfTI(unsigned short dims)
           }();
           prefer_sform_over_qform = sform_and_qform_are_very_similar;
         }
-      } // sform is orthonormal
-    }   // sform not NIFTI_XFORM_UNKNOWN
+      }
+      else if (m_SFORM_Permissive) // sform is orthonormal
+      {
+        // Fix to deal with non-orthogonal matrixes
+        // this approach uses SVD decomposition
+        // maybe it is better to use nifti_make_orthog_mat44
+        // to be consistent with other software?
+        itkWarningMacro(<< this->GetFileName() << " has non-orthogonal sform");
 
+        vnl_matrix_fixed<double, 3, 3> mat;
+
+        for (int i = 0; i < 3; ++i)
+          for (int j = 0; j < 3; ++j)
+          {
+            mat[i][j] = double{ this->m_NiftiImage->sto_xyz.m[i][j] };
+          }
+
+        vnl_svd<double> svd(mat.as_ref(), 1e-8);
+
+        if (svd.singularities() == 0)
+        {
+          mat = svd.V() * svd.U().conjugate_transpose();
+          mat44 _mat;
+
+          for (int i = 0; i < 3; ++i)
+            for (int j = 0; j < 3; ++j)
+            {
+              _mat.m[i][j] = static_cast<float>(mat[i][j]);
+            }
+
+          // preserve origin
+          for (int i = 0; i < 3; ++i)
+          {
+            _mat.m[i][3] = this->m_NiftiImage->sto_xyz.m[i][3];
+          }
+          for (int i = 0; i < 4; ++i) // should be 0 0 0 1
+          {
+            _mat.m[3][i] = this->m_NiftiImage->sto_xyz.m[3][i];
+          }
+
+          this->m_SFORM_Corrected = true;
+          return _mat;
+        }
+      }
+    } // sform not NIFTI_XFORM_UNKNOWN
     if (prefer_sform_over_qform)
     {
+      this->m_SFORM_Corrected = false;
       return this->m_NiftiImage->sto_xyz;
     }
     else if (this->m_NiftiImage->qform_code != NIFTI_XFORM_UNKNOWN)
     {
+      this->m_SFORM_Corrected = false;
       return this->m_NiftiImage->qto_xyz;
     }
 
@@ -2279,6 +2336,15 @@ NiftiImageIO::SetNIfTIOrientationFromImageIO(unsigned short origdims, unsigned s
     itkWarningMacro("Non-orthogonal direction matrix coerced to orthogonal");
   }
 
+  // also issue a warning if original file had non-orthogonal matrix?
+  {
+    const MetaDataDictionary & thisDic = this->GetMetaDataDictionary();
+    std::string                temp;
+    if (itk::ExposeMetaData<std::string>(thisDic, "nifti_sform_corrected", temp) && temp == "YES")
+    {
+      itkWarningMacro("Non-orthogonal direction matrix in original nifti file was non-orthogonal");
+    }
+  }
   // Fill in origin.
   matrix.m[0][3] = static_cast<float>(-this->GetOrigin(0));
   matrix.m[1][3] = (origdims > 1) ? static_cast<float>(-this->GetOrigin(1)) : 0.0f;

Modules/IO/NIFTI/src/itkNiftiImageIOConfigurePrivate.h.in

@@ -23,5 +23,7 @@
 
 // Analyze default flavour
 #define ITK_NIFTI_IO_ANALYZE_FLAVOR_DEFAULT NiftiImageIOEnums::Analyze75Flavor::Analyze@ITK_NIFTI_IO_ANALYZE_FLAVOR@
+// Enable permissive treatment of SFORM
+#cmakedefine01 ITK_NIFTI_IO_SFORM_PERMISSIVE_DEFAULT
 
-#endif //itkNiftiImageIOConfigurePrivate_h
+#endif // itkNiftiImageIOConfigurePrivate_h

Modules/IO/NIFTI/test/CMakeLists.txt

@@ -93,6 +93,7 @@ itk_add_test(NAME itkNiftiReadWriteDirectionTest
         DATA{${ITK_DATA_ROOT}/Input/LPSLabels.nii.gz}
         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}
         )
 
@@ -101,6 +102,7 @@ itk_add_test(NAME itkNiftiReadWriteDirectionSmallVoxelTest
         DATA{Input/SmallVoxels.nii.gz}
         DATA{Input/SmallVoxels_noqform.nii.gz}
         DATA{Input/SmallVoxels_nosform.nii.gz}
+        DATA{Input/SmallVoxelsNonOrthoSform.nii.gz}
         ${ITK_TEST_OUTPUT_DIR}
         )
 

Modules/IO/NIFTI/test/Input/SmallVoxelsNonOrthoSform.nii.gz.sha512

@@ -0,0 +1 @@
+8ee6cb2fa2f8c95f33d8b34f1556779f2a5c793680bee872531beaaa7660ef0b7d3f50fc90fb905d79e96bb37dfc63158f9b39560241ca5cc55852d76921b745

Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx

@@ -21,19 +21,67 @@
 #include "itkImageFileReader.h"
 #include "itkImageFileWriter.h"
 #include "itkTestingMacros.h"
+#include "itksys/SystemTools.hxx"
+#include "itksys/SystemInformation.hxx"
 
 // debug
 #include <map>
 
+/** ReadFile
+ * read an image from disk
+ */
+template <typename TImage>
+typename TImage::Pointer
+ReadImage(const std::string & fileName, bool SFORM_Permissive)
+{
+  using ReaderType = itk::ImageFileReader<TImage>;
+
+  auto                                reader = ReaderType::New();
+  typename itk::NiftiImageIO::Pointer imageIO = itk::NiftiImageIO::New();
+  {
+    imageIO->SetSFORM_Permissive(SFORM_Permissive);
+    reader->SetImageIO(imageIO);
+    reader->SetFileName(fileName.c_str());
+    try
+    {
+      reader->Update();
+    }
+    catch (const itk::ExceptionObject & err)
+    {
+      std::cout << "Caught an exception: " << std::endl;
+      std::cout << err << ' ' << __FILE__ << ' ' << __LINE__ << std::endl;
+      throw;
+    }
+    catch (...)
+    {
+      std::cout << "Error while reading in image  " << fileName << std::endl;
+      throw;
+    }
+  }
+  typename TImage::Pointer image = reader->GetOutput();
+  return image;
+}
+
+
+template <typename TImage>
+bool
+CheckRotation(typename TImage::Pointer img)
+{
+  vnl_matrix_fixed<double, 3, 3>       rotation = img->GetDirection().GetVnlMatrix().extract(3, 3, 0, 0);
+  const vnl_matrix_fixed<double, 3, 3> candidate_identity = rotation * rotation.transpose();
+  return candidate_identity.is_identity(1.0e-4);
+}
+
 int
 itkNiftiReadWriteDirectionTest(int argc, char * argv[])
 {
-  if (argc < 5)
+  if (argc < 6)
   {
     std::cerr << "Missing Parameters." << std::endl;
     std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv)
-              << "imageWithBothQAndSForms imageWithNoQform imageWithNoSform testOutputDir" << std::endl;
-    std::cerr << "5 arguments required, received " << argc << std::endl;
+              << "imageWithBothQAndSForms imageWithNoQform imageWithNoSform  imageWithNonOrthoSform testOutputDir"
+              << std::endl;
+    std::cerr << "6 arguments required, received " << argc << std::endl;
     for (int i = 0; i < argc; ++i)
     {
       std::cerr << '\t' << i << " : " << argv[i] << std::endl;
@@ -46,6 +94,34 @@ itkNiftiReadWriteDirectionTest(int argc, char * argv[])
   TestImageType::Pointer inputImageNoQform = itk::ReadImage<TestImageType>(argv[2]);
   TestImageType::Pointer inputImageNoSform = itk::ReadImage<TestImageType>(argv[3]);
 
+
+  // check if rotation matrix is orthogonal
+  if (!CheckRotation<TestImageType>(inputImage))
+  {
+    std::cerr << "Rotation matrix of " << argv[1] << " is not orthgonal" << std::endl;
+    return EXIT_FAILURE;
+  }
+  if (!CheckRotation<TestImageType>(inputImageNoQform))
+  {
+    std::cerr << "Rotation matrix of " << argv[2] << " is not orthgonal" << std::endl;
+    return EXIT_FAILURE;
+  }
+  if (!CheckRotation<TestImageType>(inputImageNoSform))
+  {
+    std::cerr << "Rotation matrix of " << argv[3] << " is not orthgonal" << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  itk::MetaDataDictionary & dictionary = inputImage->GetMetaDataDictionary();
+  std::string               temp;
+
+  if (!itk::ExposeMetaData<std::string>(dictionary, "ITK_sform_corrected", temp) || temp != "NO")
+  {
+    std::cerr << "ITK_sform_corrected metadata flag was not properly set" << std::endl;
+    std::cerr << " expected NO, recieved:" << temp.c_str() << std::endl;
+    return EXIT_FAILURE;
+  }
+
   const auto inputImageDirection = inputImage->GetDirection();
   const auto inputImageNoQformDirection = inputImageNoQform->GetDirection();
   const auto inputImageNoSformDirection = inputImageNoSform->GetDirection();
@@ -74,7 +150,7 @@ itkNiftiReadWriteDirectionTest(int argc, char * argv[])
   }
 
   // Write image that originally had no sform direction representation into a file with both sform and qform
-  const std::string                            testOutputDir = argv[4];
+  const std::string                            testOutputDir = argv[5];
   const std::string                            testFilename = testOutputDir + "/test_filled_sform.nii.gz";
   itk::ImageFileWriter<TestImageType>::Pointer writer = itk::ImageFileWriter<TestImageType>::New();
   ITK_TRY_EXPECT_NO_EXCEPTION(itk::WriteImage(inputImageNoSform, testFilename));
@@ -112,6 +188,34 @@ itkNiftiReadWriteDirectionTest(int argc, char * argv[])
     return EXIT_FAILURE;
   }
 
+  // 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
+  itksys::SystemTools::PutEnv("ITK_NIFTI_SFORM_PERMISSIVE=NO");
+  ITK_TRY_EXPECT_EXCEPTION(itk::ReadImage<TestImageType>(argv[4]));
+  // 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
+  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;
+    return EXIT_FAILURE;
+  }
+
+  // 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;
+  }
+
   std::cout << "Test finished." << std::endl;
   return EXIT_SUCCESS;
 }

Testing/Data/Input/NonOrthoSform.nii.gz.sha512

@@ -0,0 +1 @@
+c12c74bdb7fa32dfbee00db9009110fa1734c7702ae5002f7dd5fb6740341b0fb94c4e5f1a53bde9f1bebd45ad8b64b7643f944a66fa25b44cff5dcc0054974b