ENH: Deprecate vnl_determinant/vnl_det for Eigen-backed itk::Math::Determinant

Modules/Core/Common/include/itkExtractImageFilter.hxx

@@ -19,6 +19,7 @@
 #define itkExtractImageFilter_hxx
 
 #include "itkImageAlgorithm.h"
+#include "itkMathDeterminant.h"
 #include "itkObjectFactory.h"
 #include "itkProgressReporter.h"
 
@@ -193,15 +194,15 @@ ExtractImageFilter<TInputImage, TOutputImage>::GenerateOutputInformation()
       break;
       case DirectionCollapseStrategyEnum::DIRECTIONCOLLAPSETOSUBMATRIX:
       {
-        if (vnl_determinant(outputDirection.GetVnlMatrix()) == 0.0)
+        if (Math::Determinant(outputDirection.GetVnlMatrix()) == 0.0)
         {
           itkExceptionStringMacro("Invalid submatrix extracted for collapsed direction.");
         }
       }
       break;
       case DirectionCollapseStrategyEnum::DIRECTIONCOLLAPSETOGUESS:
       {
-        if (vnl_determinant(outputDirection.GetVnlMatrix()) == 0.0)
+        if (Math::Determinant(outputDirection.GetVnlMatrix()) == 0.0)
         {
           outputDirection.SetIdentity();
         }

Modules/Core/Common/include/itkHexahedronCell.hxx

@@ -19,7 +19,7 @@
 #define itkHexahedronCell_hxx
 #include "itkMath.h"
 #include "vnl/vnl_matrix_fixed.h"
-#include "vnl/algo/vnl_determinant.h"
+#include "itkMathDeterminant.h"
 
 #include <algorithm> // For copy_n.
 
@@ -398,7 +398,7 @@ HexahedronCell<TCellInterface>::EvaluatePosition(CoordinateType *          x,
     }
 
     // ONLY 3x3 determinants are supported.
-    const double d = vnl_determinant(mat);
+    const double d = Math::Determinant(mat);
     // spell-check-disable
     // d=vtkMath::Determinant3x3(rcol,scol,tcol);
     // spell-check-enable
@@ -431,9 +431,9 @@ HexahedronCell<TCellInterface>::EvaluatePosition(CoordinateType *          x,
       mat3.put(2, i, fcol[i]);
     }
     double params[Self::CellDimension3D]{ 0.5, 0.5, 0.5 };
-    pcoords[0] = params[0] - vnl_determinant(mat1) / d;
-    pcoords[1] = params[1] - vnl_determinant(mat2) / d;
-    pcoords[2] = params[2] - vnl_determinant(mat3) / d;
+    pcoords[0] = params[0] - Math::Determinant(mat1) / d;
+    pcoords[1] = params[1] - Math::Determinant(mat2) / d;
+    pcoords[2] = params[2] - Math::Determinant(mat3) / d;
 
     if (pcoord)
     {

Modules/Core/Common/include/itkImageBase.hxx

@@ -34,6 +34,7 @@
 #include "itkSpatialOrientation.h"
 #include <cstring>
 #include "itkMath.h"
+#include "itkMathDeterminant.h"
 
 namespace itk
 {
@@ -133,7 +134,7 @@ ImageBase<VImageDimension>::SetDirection(const DirectionType & direction)
 {
   bool modified = false;
 
-  if (vnl_determinant(direction.GetVnlMatrix()) == 0.0)
+  if (Math::Determinant(direction.GetVnlMatrix()) == 0.0)
   {
     itkExceptionMacro("Bad direction, determinant is 0. Refusing to change direction from " << this->m_Direction
                                                                                             << " to " << direction);

Modules/Core/Common/include/itkMathDeterminant.h

@@ -0,0 +1,126 @@
+/*=========================================================================
+ *
+ *  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.
+ *
+ *=========================================================================*/
+#ifndef itkMathDeterminant_h
+#define itkMathDeterminant_h
+
+#include "itkMacro.h"
+#include "vnl/vnl_matrix.h"
+#include "vnl/vnl_matrix_fixed.h"
+
+#include "itk_eigen.h"
+#include ITK_EIGEN(Dense)
+
+namespace itk
+{
+// Forward declaration lets the Matrix overload parse under the circular include
+// with itkMatrix.h; its body instantiates only where itk::Matrix is complete.
+template <typename T, unsigned int VRows, unsigned int VColumns>
+class Matrix;
+
+namespace Math
+{
+namespace detail
+{
+// Eigen chooses direct cofactor formulas for VDim <= 4 and PartialPivLU beyond;
+// the determinant is transpose-invariant, so the row-major map matches ITK
+// storage without affecting the value.
+template <unsigned int VDim, typename TReal>
+TReal
+DeterminantEigen(const TReal * inData)
+{
+  using RowMajor = Eigen::Matrix<TReal, VDim, VDim, Eigen::RowMajor>;
+  return Eigen::Map<const RowMajor>(inData).determinant();
+}
+
+template <typename TReal>
+TReal
+DynamicDeterminantEigen(const TReal * inData, unsigned int n)
+{
+  // Eigen's Dynamic-sized determinant always uses LU; dispatch small sizes to
+  // the compile-time direct formulas so runtime small matrices stay fast.
+  switch (n)
+  {
+    case 1:
+      return inData[0];
+    case 2:
+      return DeterminantEigen<2, TReal>(inData);
+    case 3:
+      return DeterminantEigen<3, TReal>(inData);
+    case 4:
+      return DeterminantEigen<4, TReal>(inData);
+    default:
+    {
+      using RowMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+      return Eigen::Map<const RowMajor>(inData, n, n).determinant();
+    }
+  }
+}
+} // namespace detail
+
+/** \brief Determinant of a square matrix, backed by Eigen.
+ *
+ * Eigen-backed replacement for vnl_determinant. Overloads accept an itk::Matrix,
+ * a vnl_matrix_fixed, or a runtime vnl_matrix over a zero-copy Eigen::Map; no
+ * Eigen type appears in the interface. Fixed sizes <= 4x4 -- which dominate ITK
+ * usage -- use Eigen's direct cofactor formulas; larger sizes use PartialPivLU.
+ *
+ * \ingroup ITKCommon
+ */
+template <typename TReal, unsigned int VDim>
+TReal
+Determinant(const vnl_matrix_fixed<TReal, VDim, VDim> & A)
+{
+  return detail::DeterminantEigen<VDim, TReal>(A.data_block());
+}
+
+/** Determinant of a fixed-size square itk::Matrix. */
+template <typename TReal, unsigned int VDim>
+TReal
+Determinant(const Matrix<TReal, VDim, VDim> & A)
+{
+  return Determinant<TReal, VDim>(A.GetVnlMatrix());
+}
+
+/** Determinant of a runtime-sized square vnl_matrix. */
+template <typename TReal>
+TReal
+Determinant(const vnl_matrix<TReal> & A)
+{
+  const unsigned int rows = A.rows();
+  if (rows != A.cols())
+  {
+    itkGenericExceptionMacro("itk::Math::Determinant requires a square matrix.");
+  }
+  return detail::DynamicDeterminantEigen<TReal>(A.data_block(), rows);
+}
+
+/** Non-square-typed fixed matrices (e.g. the 3 x PointDimension matrices in mesh
+ * cells) forward to the runtime path, which requires the actual matrix be
+ * square. The square overload above is more specialized and wins when the type
+ * is square. */
+template <typename TReal, unsigned int VRows, unsigned int VColumns>
+TReal
+Determinant(const vnl_matrix_fixed<TReal, VRows, VColumns> & A)
+{
+  return Determinant(A.as_ref());
+}
+
+} // namespace Math
+} // namespace itk
+
+#endif // itkMathDeterminant_h

Modules/Core/Common/include/itkMatrix.h

@@ -38,7 +38,13 @@
 #  include "vnl/algo/vnl_matrix_inverse.h"
 #  undef ITK_VNL_SVD_TRANSITIONAL_INCLUDE
 #endif
-#include "vnl/algo/vnl_determinant.h"
+#include "itkMathDeterminant.h"
+// GetInverse's singular check is Eigen-backed via itk::Math::Determinant, but
+// ITK 5.4 exposed vnl_determinant transitively through this header; keep it
+// reachable during the deprecation window so downstream code still compiles.
+#if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE)
+#  include "vnl/algo/vnl_determinant.h"
+#endif
 #include "itkMath.h"
 #include <type_traits> // For is_same.
 
@@ -323,7 +329,7 @@ class ITK_TEMPLATE_EXPORT Matrix
   [[nodiscard]] inline vnl_matrix_fixed<T, VColumns, VRows>
   GetInverse() const
   {
-    if (vnl_determinant(m_Matrix) == T{})
+    if (Math::Determinant(m_Matrix) == T{})
     {
       itkGenericExceptionMacro("Singular matrix. Determinant is 0.");
     }

Modules/Core/Common/include/itkQuadrilateralCell.hxx

@@ -18,7 +18,7 @@
 #ifndef itkQuadrilateralCell_hxx
 #define itkQuadrilateralCell_hxx
 #include "itkMath.h"
-#include "vnl/algo/vnl_determinant.h"
+#include "itkMathDeterminant.h"
 
 #include <algorithm> // For copy_n.
 
@@ -333,7 +333,7 @@ QuadrilateralCell<TCellInterface>::EvaluatePosition(CoordinateType *          x,
       mat.put(1, i, scol[i]);
     }
 
-    const double d = vnl_determinant(mat);
+    const double d = Math::Determinant(mat);
     // spell-check-disable
     // d=vtkMath::Determinant2x2(rcol,scol);
     // spell-check-enable
@@ -356,8 +356,8 @@ QuadrilateralCell<TCellInterface>::EvaluatePosition(CoordinateType *          x,
       mat2.put(1, i, fcol[i]);
     }
 
-    pcoords[0] = params[0] - vnl_determinant(mat1) / d;
-    pcoords[1] = params[1] - vnl_determinant(mat2) / d;
+    pcoords[0] = params[0] - Math::Determinant(mat1) / d;
+    pcoords[1] = params[1] - Math::Determinant(mat2) / d;
 
     if (pcoord)
     {

Modules/Core/Common/include/itkTetrahedronCell.hxx

@@ -17,7 +17,7 @@
  *=========================================================================*/
 #ifndef itkTetrahedronCell_hxx
 #define itkTetrahedronCell_hxx
-#include "vnl/algo/vnl_determinant.h"
+#include "itkMathDeterminant.h"
 
 #include <algorithm> // For copy_n.
 
@@ -106,7 +106,7 @@ TetrahedronCell<TCellInterface>::EvaluatePosition(CoordinateType *          x,
     mat.put(1, i, c2[i]);
     mat.put(2, i, c3[i]);
   }
-  const double det = vnl_determinant(mat);
+  const double det = Math::Determinant(mat);
   if (det == 0.0)
   {
     return false;
@@ -119,7 +119,7 @@ TetrahedronCell<TCellInterface>::EvaluatePosition(CoordinateType *          x,
     mat.put(2, i, c3[i]);
   }
 
-  pcoords[0] = vnl_determinant(mat) / det;
+  pcoords[0] = Math::Determinant(mat) / det;
 
   for (unsigned int i = 0; i < PointDimension; ++i)
   {
@@ -128,7 +128,7 @@ TetrahedronCell<TCellInterface>::EvaluatePosition(CoordinateType *          x,
     mat.put(2, i, c3[i]);
   }
 
-  pcoords[1] = vnl_determinant(mat) / det;
+  pcoords[1] = Math::Determinant(mat) / det;
 
   for (unsigned int i = 0; i < PointDimension; ++i)
   {
@@ -137,7 +137,7 @@ TetrahedronCell<TCellInterface>::EvaluatePosition(CoordinateType *          x,
     mat.put(2, i, rhs[i]);
   }
 
-  pcoords[2] = vnl_determinant(mat) / det;
+  pcoords[2] = Math::Determinant(mat) / det;
 
   const double p4 = 1.0 - pcoords[0] - pcoords[1] - pcoords[2];
 

Modules/Core/Common/include/itkTriangleCell.hxx

@@ -17,7 +17,6 @@
  *=========================================================================*/
 #ifndef itkTriangleCell_hxx
 #define itkTriangleCell_hxx
-#include "vnl/algo/vnl_determinant.h"
 
 #include <algorithm> // For copy_n.
 #include <cmath>     // For abs.

Modules/Core/Common/include/itkVersor.hxx

@@ -20,7 +20,7 @@
 
 #include "itkNumericTraits.h"
 #include "itkMath.h"
-#include <vnl/vnl_det.h>
+#include "itkMathDeterminant.h"
 
 namespace itk
 {
@@ -327,12 +327,12 @@ Versor<T>::Set(const MatrixType & mat)
       itk::Math::Absolute(I[2][0]) > epsilon || itk::Math::Absolute(I[2][1]) > epsilon ||
       itk::Math::Absolute(I[0][0] - itk::NumericTraits<T>::OneValue()) > epsilonDiff ||
       itk::Math::Absolute(I[1][1] - itk::NumericTraits<T>::OneValue()) > epsilonDiff ||
-      itk::Math::Absolute(I[2][2] - itk::NumericTraits<T>::OneValue()) > epsilonDiff || vnl_det(I) < 0)
+      itk::Math::Absolute(I[2][2] - itk::NumericTraits<T>::OneValue()) > epsilonDiff || Math::Determinant(I) < 0)
   {
     itkGenericExceptionMacro("The following matrix does not represent rotation to within an epsion of "
                              << epsilon << '.' << std::endl
                              << m << std::endl
-                             << "det(m * m transpose) is: " << vnl_det(I) << std::endl
+                             << "det(m * m transpose) is: " << Math::Determinant(I) << std::endl
                              << "m * m transpose is:" << std::endl
                              << I);
   }

Modules/Core/Common/test/CMakeLists.txt

@@ -1487,6 +1487,7 @@ set(
   itkMathSVDGTest.cxx
   itkMathLDLTGTest.cxx
   itkVnlCholeskyEngineGTest.cxx
+  itkMathDeterminantGTest.cxx
   itkVnlSVDEngineGTest.cxx
   itkMatrixExponentialGTest.cxx
   itkMatrixGTest.cxx