diff --git a/CMakeLists.txt b/CMakeLists.txt index 54df5f0e9bc..8ad2077e95d 100644 --- a/CMakeLists.txt +++ b/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 "Allow non-orthogonal rotation matrix in NIFTI sform" ON) +mark_as_advanced(ITK_NIFTI_IO_SFORM_PERMISSIVE) + + #----------------------------------------------------------------------------- # ITK build classes that are in the review process diff --git a/Modules/IO/NIFTI/include/itkNiftiImageIO.h b/Modules/IO/NIFTI/include/itkNiftiImageIO.h index 4f2e4767250..32a2a8e96f0 100644 --- a/Modules/IO/NIFTI/include/itkNiftiImageIO.h +++ b/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 }; }; diff --git a/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx b/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx index c179e79b4bf..3e92a3edfa6 100644 --- a/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx +++ b/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx @@ -22,6 +22,7 @@ #include #include "itkNiftiImageIOConfigurePrivate.h" #include "itkMakeUniqueForOverwrite.h" +#include "vnl/algo/vnl_svd.h" namespace itk { @@ -449,6 +450,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 } { this->SetNumberOfDimensions(3); nifti_set_debug_level(0); // suppress error messages @@ -481,6 +483,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 +902,8 @@ NiftiImageIO::SetImageIOMetadataFromNIfTI() // Necessary to clear dict if ImageIO object is re-used thisDic.Clear(); + EncapsulateMetaData(thisDic, "ITK_sform_corrected", m_SFORM_Corrected? "YES":"NO"); + std::ostringstream nifti_type; nifti_type << nim->nifti_type; EncapsulateMetaData(thisDic, "nifti_type", nifti_type.str()); @@ -2031,6 +2036,7 @@ NiftiImageIO::SetImageIOOrientationFromNIfTI(unsigned short dims) // Only orthonormal matricies have transpose as inverse const vnl_matrix_fixed candidate_identity = rotation * rotation.transpose(); const bool is_orthonormal = candidate_identity.is_identity(1.0e-4); + return is_orthonormal; } }(); @@ -2089,15 +2095,48 @@ 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 mat; + + for (int i = 0; i < 4; ++i) + for (int j = 0; j < 4; ++j) + { + mat[i][j] = double{ this->m_NiftiImage->sto_xyz.m[i][j] }; + } + + vnl_svd svd(mat.as_matrix(), 1e-8); + if (svd.singularities()==0) + { + mat = svd.V() * svd.U().conjugate_transpose(); + mat44 _mat; + for (int i = 0; i < 4; ++i) + for (int j = 0; j < 4; ++j) + { + _mat.m[i][j] = float{ mat[i][j] }; + } + + 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 +2318,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(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(-this->GetOrigin(0)); matrix.m[1][3] = (origdims > 1) ? static_cast(-this->GetOrigin(1)) : 0.0f; diff --git a/Modules/IO/NIFTI/src/itkNiftiImageIOConfigurePrivate.h.in b/Modules/IO/NIFTI/src/itkNiftiImageIOConfigurePrivate.h.in index e860bc020c3..789405e33a0 100644 --- a/Modules/IO/NIFTI/src/itkNiftiImageIOConfigurePrivate.h.in +++ b/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 +#cmakedefine ITK_NIFTI_IO_SFORM_PERMISSIVE #endif //itkNiftiImageIOConfigurePrivate_h diff --git a/Modules/IO/NIFTI/test/CMakeLists.txt b/Modules/IO/NIFTI/test/CMakeLists.txt index c7df53a88e8..76a1566ffaa 100644 --- a/Modules/IO/NIFTI/test/CMakeLists.txt +++ b/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} ) diff --git a/Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx b/Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx index a9b91b27acc..ed083a5ccb4 100644 --- a/Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx +++ b/Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx @@ -25,15 +25,51 @@ // debug #include +/** ReadFile + * read an image from disk + */ +template +typename TImage::Pointer +ReadImage(const std::string & fileName, + bool SFORM_Permissive) +{ + using ReaderType = itk::ImageFileReader; + + 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; +} + 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 +82,16 @@ itkNiftiReadWriteDirectionTest(int argc, char * argv[]) TestImageType::Pointer inputImageNoQform = itk::ReadImage(argv[2]); TestImageType::Pointer inputImageNoSform = itk::ReadImage(argv[3]); + itk::MetaDataDictionary & dictionary = inputImage->GetMetaDataDictionary(); + std::string temp; + + if (!itk::ExposeMetaData(dictionary, "ITK_sform_corrected", temp) || temp!="NO") + { + std::cerr<<"ITK_sform_corrected metadata flag was not properly set"<GetDirection(); const auto inputImageNoQformDirection = inputImageNoQform->GetDirection(); const auto inputImageNoSformDirection = inputImageNoSform->GetDirection(); @@ -74,7 +120,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::Pointer writer = itk::ImageFileWriter::New(); ITK_TRY_EXPECT_NO_EXCEPTION(itk::WriteImage(inputImageNoSform, testFilename)); @@ -112,6 +158,19 @@ itkNiftiReadWriteDirectionTest(int argc, char * argv[]) return EXIT_FAILURE; } + // should not read the image with non-orthogonal sform + ITK_TRY_EXPECT_EXCEPTION(ReadImage(argv[4],false)); + + // this should work + TestImageType::Pointer inputImageNonOrthoSform = ReadImage(argv[4],true); + dictionary = inputImageNonOrthoSform->GetMetaDataDictionary(); + if (!itk::ExposeMetaData(dictionary, "ITK_sform_corrected", temp) || temp!="YES") + { + std::cerr<<"ITK_sform_corrected metadata flag was not properly set"<