BUG: Fix isAffine insuffcient precision in caluclations

Fixes #3333

Co-authored-by: Vladimir S. FONOV <vladimir.fonov@gmail.com>

Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

@@ -1753,10 +1753,15 @@ Normalize(std::vector<double> & x)
 static bool
 IsAffine(const mat44 & nifti_mat)
 {
-  const vnl_matrix_fixed<float, 4, 4> mat(&(nifti_mat.m[0][0]));
+  vnl_matrix_fixed<double, 4, 4> mat;
+
+  for (int i = 0; i < 4; ++i)
+    for (int j = 0; j < 4; ++j)
+      mat[i][j] = double{ nifti_mat.m[i][j] };
+
   // First make sure the bottom row meets the condition that it is (0, 0, 0, 1)
   {
-    float bottom_row_error = itk::Math::abs(mat[3][3] - 1.0f);
+    double bottom_row_error = itk::Math::abs(mat[3][3] - 1.0);
     for (int i = 0; i < 3; ++i)
     {
       bottom_row_error += itk::Math::abs(mat[3][i]);
@@ -1767,28 +1772,29 @@ IsAffine(const mat44 & nifti_mat)
     }
   }
 
-  const float condition = vnl_matrix_inverse<float>(mat.as_matrix()).well_condition();
+  const double condition = vnl_matrix_inverse<double>(mat.as_matrix()).well_condition();
   // Check matrix is invertible by testing condition number of inverse
-  if (!(condition > std::numeric_limits<float>::epsilon()))
+  if (!(condition > std::numeric_limits<double>::epsilon()))
   {
     return false;
   }
 
   // Calculate the inverse and separate the inverse translation component
   // and the top 3x3 part of the inverse matrix
-  const vnl_matrix_fixed<float, 4, 4> inv4x4Matrix = vnl_matrix_inverse<float>(mat.as_matrix()).as_matrix();
-  const vnl_vector_fixed<float, 3>    inv4x4Translation(inv4x4Matrix[0][3], inv4x4Matrix[1][3], inv4x4Matrix[2][3]);
-  const vnl_matrix_fixed<float, 3, 3> inv4x4Top3x3 = inv4x4Matrix.extract(3, 3, 0, 0);
+  const vnl_matrix_fixed<double, 4, 4> inv4x4Matrix = vnl_matrix_inverse<double>(mat.as_matrix()).as_matrix();
+  const vnl_vector_fixed<double, 3>    inv4x4Translation(inv4x4Matrix[0][3], inv4x4Matrix[1][3], inv4x4Matrix[2][3]);
+  const vnl_matrix_fixed<double, 3, 3> inv4x4Top3x3 = inv4x4Matrix.extract(3, 3, 0, 0);
 
   // Grab just the top 3x3 matrix
-  const vnl_matrix_fixed<float, 3, 3> top3x3Matrix = mat.extract(3, 3, 0, 0);
-  const vnl_matrix_fixed<float, 3, 3> invTop3x3Matrix = vnl_matrix_inverse<float>(top3x3Matrix.as_matrix()).as_matrix();
-  const vnl_vector_fixed<float, 3>    inv3x3Translation = -(invTop3x3Matrix * mat.get_column(3).extract(3));
+  const vnl_matrix_fixed<double, 3, 3> top3x3Matrix = mat.extract(3, 3, 0, 0);
+  const vnl_matrix_fixed<double, 3, 3> invTop3x3Matrix =
+    vnl_matrix_inverse<double>(top3x3Matrix.as_matrix()).as_matrix();
+  const vnl_vector_fixed<double, 3> inv3x3Translation = -(invTop3x3Matrix * mat.get_column(3).extract(3));
 
   // Make sure we adhere to the conditions of a 4x4 invertible affine transform matrix
-  const float     diff_matrix_array_one_norm = (inv4x4Top3x3 - invTop3x3Matrix).array_one_norm();
-  const float     diff_vector_translation_one_norm = (inv4x4Translation - inv3x3Translation).one_norm();
-  constexpr float normed_tolerance_matrix_close = 1e-2;
+  const double     diff_matrix_array_one_norm = (inv4x4Top3x3 - invTop3x3Matrix).array_one_norm();
+  const double     diff_vector_translation_one_norm = (inv4x4Translation - inv3x3Translation).one_norm();
+  constexpr double normed_tolerance_matrix_close = 1e-2;
   return !((diff_matrix_array_one_norm > normed_tolerance_matrix_close) ||
            (diff_vector_translation_one_norm > normed_tolerance_matrix_close));
 }

Modules/IO/NIFTI/test/CMakeLists.txt

@@ -12,6 +12,7 @@ itkNiftiImageIOTest9.cxx
 itkNiftiImageIOTest10.cxx
 itkNiftiImageIOTest11.cxx
 itkNiftiImageIOTest12.cxx
+itkNiftiImageIOTest13.cxx
 itkNiftiReadAnalyzeTest.cxx
 itkNiftiReadWriteDirectionTest.cxx
 itkExtractSlice.cxx
@@ -96,3 +97,7 @@ itk_add_test(NAME itkNiftiReadWriteDirectionSmallVoxelTest
         DATA{Input/SmallVoxels_nosform.nii.gz}
         ${ITK_TEST_OUTPUT_DIR}
         )
+
+itk_add_test(NAME itkNiftiSmallVoxelsAffinePrecisionTest
+COMMAND ITKIONIFTITestDriver itkNiftiImageIOTest13
+      DATA{Input/SmallVoxels_AffinePrecision.nii.gz})

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

@@ -0,0 +1 @@
+5afad4b7b6fe872a896ed7eb555a572dc893e8074555973aca1bff3d2c12f6549a6d79b88d8c869c4e0b514de8b049a92b74fc09765fa82cbb21494f4e6608d6

Modules/IO/NIFTI/test/itkNiftiImageIOTest13.cxx

@@ -0,0 +1,52 @@
+/*=========================================================================
+ *
+ *  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
+ *
+ *         http://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 "itkNiftiImageIOTest.h"
+
+// Test to read small voxel NIFTI file which was triggering numerical instability
+// in IsAffine loading code
+int
+itkNiftiImageIOTest13(int argc, char * argv[])
+{
+  if (argc != 2)
+  {
+    return EXIT_FAILURE;
+  }
+
+  constexpr unsigned int Dimension = 3;
+
+  using PixelType = float;
+  using ImageType = itk::Image<PixelType, Dimension>;
+
+  try
+  {
+    using ReaderType = itk::ImageFileReader<ImageType>;
+    ReaderType::Pointer reader = ReaderType::New();
+    reader->SetFileName(argv[1]);
+    reader->Update();
+    ImageType::Pointer image = reader->GetOutput();
+  }
+  catch (const itk::ExceptionObject & err)
+  {
+    std::cerr << "ExceptionObject caught !" << std::endl;
+    std::cerr << err << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  return EXIT_SUCCESS;
+}