ENH: Issue warning when NIFTI IO coerces matrix to orthonormal

NIFTI IO silently coerced non-orhonormal matrices to be
orthonormal. We now issue a warning when this is done.

addresses #3476

Author: zivy

Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

@@ -2151,6 +2151,24 @@ NiftiImageIO::SetNIfTIOrientationFromImageIO(unsigned short origdims, unsigned s
   mat44 matrix =
     nifti_make_orthog_mat44(dirx[0], dirx[1], dirx[2], diry[0], diry[1], diry[2], dirz[0], dirz[1], dirz[2]);
   matrix = mat44_transpose(matrix);
+  // Check if matrix is orthogonal and issue a warning if it was
+  // coerced to orthogonal. Use same epsilon as used in SetImageIOOrientationFromNIfTI
+  const unsigned int matDim(this->GetDirection(0).size());
+  vnl_matrix<float>  imageMat(matDim, matDim);
+  for (unsigned c = 0; c < matDim; ++c)
+  {
+    auto col = this->GetDirection(c);
+    for (unsigned r = 0; r < matDim; ++r)
+    {
+      imageMat[r][c] = col[r];
+    }
+  }
+  auto candidateIdentity = imageMat * imageMat.transpose();
+  if (!candidateIdentity.is_identity(1.0e-4))
+  {
+    itkWarningMacro("Non-orthogonal direction matrix coerced to 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/test/CMakeLists.txt

@@ -16,6 +16,7 @@ itkNiftiImageIOTest13.cxx
 itkNiftiReadAnalyzeTest.cxx
 itkNiftiReadWriteDirectionTest.cxx
 itkExtractSlice.cxx
+itkNiftiWriteCoerceOrthogonalDirectionTest.cxx
 )
 
 # For itkNiftiImageIOTest.h.
@@ -101,3 +102,8 @@ itk_add_test(NAME itkNiftiReadWriteDirectionSmallVoxelTest
 itk_add_test(NAME itkNiftiSmallVoxelsAffinePrecisionTest
 COMMAND ITKIONIFTITestDriver itkNiftiImageIOTest13
       DATA{Input/SmallVoxels_AffinePrecision.nii.gz})
+
+itk_add_test(NAME itkNiftiWriteCoerceOrthogonalDirectionTest
+        COMMAND ITKIONIFTITestDriver itkNiftiWriteCoerceOrthogonalDirectionTest
+        ${ITK_TEST_OUTPUT_DIR}
+        )

Modules/IO/NIFTI/test/itkNiftiWriteCoerceOrthogonalDirectionTest.cxx

@@ -0,0 +1,66 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         https://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#include <iostream>
+#include "itkIOTestHelper.h"
+#include "itkTestingMacros.h"
+#include "itkNiftiImageIO.h"
+
+int
+itkNiftiWriteCoerceOrthogonalDirectionTest(int argc, char * argv[])
+{
+  if (argc < 2)
+  {
+    std::cerr << "Missing Parameters." << std::endl;
+    std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv) << "testOutputDir" << std::endl;
+    std::cerr << "1 argument required, received " << argc << std::endl;
+    for (int i = 0; i < argc; ++i)
+    {
+      std::cerr << "\t" << i << " : " << argv[i] << std::endl;
+    }
+    return EXIT_FAILURE;
+  }
+
+  const unsigned int dim = 2;
+  using ImageType = itk::Image<unsigned char, dim>;
+
+  ImageType::IndexType  startIndex = { { 0, 0 } };
+  ImageType::SizeType   imageSize = { { 2, 2 } };
+  ImageType::RegionType region;
+  region.SetSize(imageSize);
+  region.SetIndex(startIndex);
+  auto image1 = ImageType::New();
+  image1->SetRegions(region);
+  image1->Allocate();
+
+  ImageType::DirectionType mat1;
+  mat1.SetIdentity();
+  mat1[0][0] = 0.5; // make matrix non-orthonormal
+  image1->SetDirection(mat1);
+
+  const std::string outputDir = argv[1];
+  const std::string outputFilename = outputDir + "/coercedDirection.nii.gz";
+  itk::IOTestHelper::WriteImage<ImageType, itk::NiftiImageIO>(image1, outputFilename);
+  auto image2 = itk::IOTestHelper::ReadImage<ImageType>(outputFilename);
+  // Nifti image writing coerces the direction matrix to be
+  // orthonormal, so the matrices are expected to be different
+  if (image1->GetDirection() == image2->GetDirection())
+  {
+    return EXIT_FAILURE;
+  }
+  return EXIT_SUCCESS;
+}