diff --git a/Modules/Core/Common/include/itkExtractImageFilter.hxx b/Modules/Core/Common/include/itkExtractImageFilter.hxx index f546cfef284..c84087f98e5 100644 --- a/Modules/Core/Common/include/itkExtractImageFilter.hxx +++ b/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,7 +194,7 @@ ExtractImageFilter::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."); } @@ -201,7 +202,7 @@ ExtractImageFilter::GenerateOutputInformation() break; case DirectionCollapseStrategyEnum::DIRECTIONCOLLAPSETOGUESS: { - if (vnl_determinant(outputDirection.GetVnlMatrix()) == 0.0) + if (Math::Determinant(outputDirection.GetVnlMatrix()) == 0.0) { outputDirection.SetIdentity(); } diff --git a/Modules/Core/Common/include/itkHexahedronCell.hxx b/Modules/Core/Common/include/itkHexahedronCell.hxx index 2d6fe6121ea..f1ad184b6a7 100644 --- a/Modules/Core/Common/include/itkHexahedronCell.hxx +++ b/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 // For copy_n. @@ -398,7 +398,7 @@ HexahedronCell::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::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) { diff --git a/Modules/Core/Common/include/itkImageBase.hxx b/Modules/Core/Common/include/itkImageBase.hxx index d0428454da9..19b480950bd 100644 --- a/Modules/Core/Common/include/itkImageBase.hxx +++ b/Modules/Core/Common/include/itkImageBase.hxx @@ -34,6 +34,7 @@ #include "itkSpatialOrientation.h" #include #include "itkMath.h" +#include "itkMathDeterminant.h" namespace itk { @@ -133,7 +134,7 @@ ImageBase::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); diff --git a/Modules/Core/Common/include/itkMathDeterminant.h b/Modules/Core/Common/include/itkMathDeterminant.h new file mode 100644 index 00000000000..83db69055bd --- /dev/null +++ b/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 +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 +TReal +DeterminantEigen(const TReal * inData) +{ + using RowMajor = Eigen::Matrix; + return Eigen::Map(inData).determinant(); +} + +template +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; + return Eigen::Map(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 +TReal +Determinant(const vnl_matrix_fixed & A) +{ + return detail::DeterminantEigen(A.data_block()); +} + +/** Determinant of a fixed-size square itk::Matrix. */ +template +TReal +Determinant(const Matrix & A) +{ + return Determinant(A.GetVnlMatrix()); +} + +/** Determinant of a runtime-sized square vnl_matrix. */ +template +TReal +Determinant(const vnl_matrix & A) +{ + const unsigned int rows = A.rows(); + if (rows != A.cols()) + { + itkGenericExceptionMacro("itk::Math::Determinant requires a square matrix."); + } + return detail::DynamicDeterminantEigen(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 +TReal +Determinant(const vnl_matrix_fixed & A) +{ + return Determinant(A.as_ref()); +} + +} // namespace Math +} // namespace itk + +#endif // itkMathDeterminant_h diff --git a/Modules/Core/Common/include/itkMatrix.h b/Modules/Core/Common/include/itkMatrix.h index 84931fafee0..0115acd2b1c 100644 --- a/Modules/Core/Common/include/itkMatrix.h +++ b/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 // For is_same. @@ -323,7 +329,7 @@ class ITK_TEMPLATE_EXPORT Matrix [[nodiscard]] inline vnl_matrix_fixed GetInverse() const { - if (vnl_determinant(m_Matrix) == T{}) + if (Math::Determinant(m_Matrix) == T{}) { itkGenericExceptionMacro("Singular matrix. Determinant is 0."); } diff --git a/Modules/Core/Common/include/itkQuadrilateralCell.hxx b/Modules/Core/Common/include/itkQuadrilateralCell.hxx index 32b8845cc10..f81eec486d1 100644 --- a/Modules/Core/Common/include/itkQuadrilateralCell.hxx +++ b/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 // For copy_n. @@ -333,7 +333,7 @@ QuadrilateralCell::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::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) { diff --git a/Modules/Core/Common/include/itkTetrahedronCell.hxx b/Modules/Core/Common/include/itkTetrahedronCell.hxx index f8021181b4a..816c515b2ae 100644 --- a/Modules/Core/Common/include/itkTetrahedronCell.hxx +++ b/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 // For copy_n. @@ -106,7 +106,7 @@ TetrahedronCell::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::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::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::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]; diff --git a/Modules/Core/Common/include/itkTriangleCell.hxx b/Modules/Core/Common/include/itkTriangleCell.hxx index af1518cf29c..f0f08608421 100644 --- a/Modules/Core/Common/include/itkTriangleCell.hxx +++ b/Modules/Core/Common/include/itkTriangleCell.hxx @@ -17,7 +17,6 @@ *=========================================================================*/ #ifndef itkTriangleCell_hxx #define itkTriangleCell_hxx -#include "vnl/algo/vnl_determinant.h" #include // For copy_n. #include // For abs. diff --git a/Modules/Core/Common/include/itkVersor.hxx b/Modules/Core/Common/include/itkVersor.hxx index dfe60f51fa6..f9f633e9631 100644 --- a/Modules/Core/Common/include/itkVersor.hxx +++ b/Modules/Core/Common/include/itkVersor.hxx @@ -20,7 +20,7 @@ #include "itkNumericTraits.h" #include "itkMath.h" -#include +#include "itkMathDeterminant.h" namespace itk { @@ -327,12 +327,12 @@ Versor::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::OneValue()) > epsilonDiff || itk::Math::Absolute(I[1][1] - itk::NumericTraits::OneValue()) > epsilonDiff || - itk::Math::Absolute(I[2][2] - itk::NumericTraits::OneValue()) > epsilonDiff || vnl_det(I) < 0) + itk::Math::Absolute(I[2][2] - itk::NumericTraits::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); } diff --git a/Modules/Core/Common/test/CMakeLists.txt b/Modules/Core/Common/test/CMakeLists.txt index 1593d8456f6..668dece3c95 100644 --- a/Modules/Core/Common/test/CMakeLists.txt +++ b/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 diff --git a/Modules/Core/Common/test/itkMathDeterminantGTest.cxx b/Modules/Core/Common/test/itkMathDeterminantGTest.cxx new file mode 100644 index 00000000000..e2f777a9d02 --- /dev/null +++ b/Modules/Core/Common/test/itkMathDeterminantGTest.cxx @@ -0,0 +1,179 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// First include the header file to be tested: +#include "itkMathDeterminant.h" + +#include "itkMatrix.h" + +// This test cross-checks against the vnl reference engines; mark it so it keeps +// compiling once vnl_determinant/vnl_det are deprecated under ITK_LEGACY_REMOVE. +#define ITK_LEGACY_TEST +#include "vnl/algo/vnl_determinant.h" +#include "vnl/vnl_det.h" + +#include +#include + +namespace +{ +// The analytic 2x2/3x3/4x4 fixtures from vnl's test_determinant, whose values +// are verified below by an independent closed-form expansion. +const double kM2[2][2] = { { 0.60684258354179, 0.89129896614890 }, { 0.48598246870930, 0.76209683302739 } }; +const double kM3[3][3] = { { 0.45646766516834, 0.44470336435319, 0.92181297074480 }, + { 0.01850364324822, 0.61543234810009, 0.73820724581067 }, + { 0.82140716429525, 0.79193703742704, 0.17626614449462 } }; + +template +itk::Matrix +FillMatrix(const double src[VDim][VDim]) +{ + itk::Matrix A; + for (unsigned int i = 0; i < VDim; ++i) + for (unsigned int j = 0; j < VDim; ++j) + A(i, j) = static_cast(src[i][j]); + return A; +} +} // namespace + +// Closed-form 2x2 reference: ad - bc. +TEST(MathDeterminant, TwoByTwoAnalytic) +{ + const auto A = FillMatrix(kM2); + const double expected = kM2[0][0] * kM2[1][1] - kM2[0][1] * kM2[1][0]; + EXPECT_NEAR(itk::Math::Determinant(A), expected, 1e-14); +} + +// Closed-form 3x3 reference: cofactor expansion along the first row. +TEST(MathDeterminant, ThreeByThreeAnalytic) +{ + const auto A = FillMatrix(kM3); + const double expected = kM3[0][0] * (kM3[1][1] * kM3[2][2] - kM3[1][2] * kM3[2][1]) - + kM3[0][1] * (kM3[1][0] * kM3[2][2] - kM3[1][2] * kM3[2][0]) + + kM3[0][2] * (kM3[1][0] * kM3[2][1] - kM3[1][1] * kM3[2][0]); + EXPECT_NEAR(itk::Math::Determinant(A), expected, 1e-14); +} + +// The three overloads (itk::Matrix, vnl_matrix_fixed, vnl_matrix) return the +// same value for the same data. +TEST(MathDeterminant, OverloadConsistency) +{ + const auto A = FillMatrix(kM3); + const vnl_matrix_fixed fixedA = A.GetVnlMatrix(); + const vnl_matrix dynA = fixedA.as_matrix(); + + const double d = itk::Math::Determinant(A); + EXPECT_DOUBLE_EQ(itk::Math::Determinant(fixedA), d); + EXPECT_NEAR(itk::Math::Determinant(dynA), d, 1e-14); +} + +// Identity has unit determinant; a diagonal matrix yields the product of its +// diagonal; both stress the direct fixed-size path. +TEST(MathDeterminant, IdentityAndDiagonal) +{ + EXPECT_DOUBLE_EQ(itk::Math::Determinant(itk::Matrix::GetIdentity()), 1.0); + + itk::Matrix D; + D.SetIdentity(); + D(0, 0) = 2.0; + D(1, 1) = -3.0; + D(2, 2) = 0.5; + D(3, 3) = 4.0; + EXPECT_DOUBLE_EQ(itk::Math::Determinant(D), 2.0 * -3.0 * 0.5 * 4.0); +} + +// A rank-deficient matrix (two identical rows) has a zero determinant. +TEST(MathDeterminant, SingularIsZero) +{ + itk::Matrix A; + A(0, 0) = 1.0; + A(0, 1) = 2.0; + A(0, 2) = 3.0; + A(1, 0) = 1.0; + A(1, 1) = 2.0; + A(1, 2) = 3.0; + A(2, 0) = 4.0; + A(2, 1) = 5.0; + A(2, 2) = 6.0; + EXPECT_NEAR(itk::Math::Determinant(A), 0.0, 1e-14); +} + +// Upper-triangular determinant is the product of the diagonal; exercises the +// runtime PartialPivLU path at a size beyond the direct formulas. +TEST(MathDeterminant, LargeTriangular) +{ + constexpr unsigned int N = 6; + vnl_matrix A(N, N, 0.0); + double expected = 1.0; + for (unsigned int i = 0; i < N; ++i) + { + A(i, i) = static_cast(i) - 2.5; // nonzero diagonal, mixed sign + expected *= A(i, i); + for (unsigned int j = i + 1; j < N; ++j) + { + A(i, j) = 0.3 * (i + 1) + 0.7 * j; // upper triangle does not affect det + } + } + EXPECT_NEAR(itk::Math::Determinant(A), expected, 1e-10 * std::abs(expected)); +} + +// det(cA) = c^n det(A) for an n x n matrix. +TEST(MathDeterminant, ScalingLaw) +{ + const auto A = FillMatrix(kM3); + const double base = itk::Math::Determinant(A); + itk::Matrix scaled = A * 2.0; + EXPECT_NEAR(itk::Math::Determinant(scaled), 8.0 * base, 1e-13); +} + +// float instantiation computes with float precision. +TEST(MathDeterminant, FloatType) +{ + const auto A = FillMatrix(kM2); + const float expected = static_cast(kM2[0][0] * kM2[1][1] - kM2[0][1] * kM2[1][0]); + EXPECT_NEAR(itk::Math::Determinant(A), expected, 1e-6f); +} + +// A non-square dynamic matrix is a usage error. +TEST(MathDeterminant, NonSquareThrows) +{ + const vnl_matrix rect(2, 3, 1.0); + EXPECT_THROW(itk::Math::Determinant(rect), itk::ExceptionObject); +} + +// Equivalence with the vnl engines being supplemented, on well- and +// ill-conditioned inputs, by relative agreement. +TEST(MathDeterminant, MatchesVnl) +{ + const auto A3 = FillMatrix(kM3); + const vnl_matrix_fixed f3 = A3.GetVnlMatrix(); + EXPECT_NEAR(itk::Math::Determinant(A3), vnl_det(f3), 1e-13); + + const vnl_matrix d3 = f3.as_matrix(); + const double ref = vnl_determinant(d3); + EXPECT_NEAR(itk::Math::Determinant(d3), ref, 1e-13 * std::max(1.0, std::abs(ref))); + + // Ill-conditioned: a scaled Hilbert-like matrix. + constexpr unsigned int N = 8; + vnl_matrix H(N, N); + for (unsigned int i = 0; i < N; ++i) + for (unsigned int j = 0; j < N; ++j) + H(i, j) = 1.0 / (i + j + 1.0); + const double refH = vnl_determinant(H); + EXPECT_NEAR(itk::Math::Determinant(H), refH, 1e-9 * std::max(1.0, std::abs(refH))); +} diff --git a/Modules/Core/TestKernel/include/itkTestingExtractSliceImageFilter.hxx b/Modules/Core/TestKernel/include/itkTestingExtractSliceImageFilter.hxx index 4249e0d5c50..dfcd9b3dc31 100644 --- a/Modules/Core/TestKernel/include/itkTestingExtractSliceImageFilter.hxx +++ b/Modules/Core/TestKernel/include/itkTestingExtractSliceImageFilter.hxx @@ -19,6 +19,7 @@ #define itkTestingExtractSliceImageFilter_hxx #include "itkImageRegionIterator.h" +#include "itkMathDeterminant.h" #include "itkObjectFactory.h" #include "itkTotalProgressReporter.h" @@ -200,7 +201,7 @@ ExtractSliceImageFilter::GenerateOutputInformation() break; case TestExtractSliceImageFilterCollapseStrategyEnum::DIRECTIONCOLLAPSETOSUBMATRIX: { - if (vnl_determinant(outputDirection.GetVnlMatrix()) == 0.0) + if (Math::Determinant(outputDirection.GetVnlMatrix()) == 0.0) { itkExceptionStringMacro("Invalid submatrix extracted for collapsed direction."); } @@ -208,7 +209,7 @@ ExtractSliceImageFilter::GenerateOutputInformation() break; case TestExtractSliceImageFilterCollapseStrategyEnum::DIRECTIONCOLLAPSETOGUESS: { - if (vnl_determinant(outputDirection.GetVnlMatrix()) == 0.0) + if (Math::Determinant(outputDirection.GetVnlMatrix()) == 0.0) { outputDirection.SetIdentity(); } diff --git a/Modules/Core/Transform/include/itkComposeScaleSkewVersor3DTransform.hxx b/Modules/Core/Transform/include/itkComposeScaleSkewVersor3DTransform.hxx index 40065d47dcc..b84e74bb07c 100644 --- a/Modules/Core/Transform/include/itkComposeScaleSkewVersor3DTransform.hxx +++ b/Modules/Core/Transform/include/itkComposeScaleSkewVersor3DTransform.hxx @@ -19,6 +19,7 @@ #define itkComposeScaleSkewVersor3DTransform_hxx #include "itkMath.h" +#include "itkMathDeterminant.h" #include "vnl/vnl_inverse.h" @@ -280,7 +281,7 @@ ComposeScaleSkewVersor3DTransform::ComputeMatrixParameters M(2, 2) /= m_Scale[2]; m_Skew[1] = ortho0 / m_Scale[0]; m_Skew[2] = ortho1 / m_Scale[1]; - if (vnl_determinant(M.GetVnlMatrix()) < 0) + if (Math::Determinant(M.GetVnlMatrix()) < 0) { m_Scale[0] *= -1; M(0, 0) *= -1; diff --git a/Modules/Core/Transform/include/itkSimilarity3DTransform.hxx b/Modules/Core/Transform/include/itkSimilarity3DTransform.hxx index 7ebd97e10cb..a5d98b9122d 100644 --- a/Modules/Core/Transform/include/itkSimilarity3DTransform.hxx +++ b/Modules/Core/Transform/include/itkSimilarity3DTransform.hxx @@ -19,7 +19,7 @@ #define itkSimilarity3DTransform_hxx #include "itkMath.h" -#include "vnl/vnl_det.h" +#include "itkMathDeterminant.h" #include "itkPrintHelper.h" namespace itk @@ -85,7 +85,7 @@ Similarity3DTransform::SetMatrix(const MatrixType & matrix // multiplied by the scale factor, then its determinant // must be equal to the cube of the scale factor. // - const double det = vnl_det(matrix.GetVnlMatrix()); + const double det = Math::Determinant(matrix.GetVnlMatrix()); if (det == 0.0) { @@ -290,7 +290,7 @@ Similarity3DTransform::ComputeMatrixParameters() { MatrixType matrix = this->GetMatrix(); - m_Scale = itk::Math::cbrt(vnl_det(matrix.GetVnlMatrix())); + m_Scale = itk::Math::cbrt(Math::Determinant(matrix.GetVnlMatrix())); matrix /= m_Scale; diff --git a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.h b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.h index 61a4aef91b4..ca301db6cc4 100644 --- a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.h +++ b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.h @@ -23,7 +23,6 @@ #include "itkVector.h" #include "vnl/vnl_matrix_fixed.h" #include "vnl/vnl_vector_fixed.h" -#include "vnl/vnl_det.h" namespace itk { @@ -221,7 +220,7 @@ class ITK_TEMPLATE_EXPORT DisplacementFieldJacobianDeterminantFilter BeforeThreadedGenerateData() override; /** DisplacementFieldJacobianDeterminantFilter can be implemented as a - * multithreaded filter (we're only using vnl_det(), which is trivially + * multithreaded filter (the determinant computation is trivially * thread safe). Therefore, this implementation provides a * DynamicThreadedGenerateData() routine which is called for each * processing thread. The output image data is allocated diff --git a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx index b75e10e7b5c..65e11859c7b 100644 --- a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx +++ b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx @@ -27,7 +27,7 @@ #include "itkMath.h" #include "itkPrintHelper.h" -#include "vnl/algo/vnl_determinant.h" +#include "itkMathDeterminant.h" namespace itk { @@ -251,7 +251,7 @@ DisplacementFieldJacobianDeterminantFilter } // Return determinant of physical Jacobian: - return vnl_determinant(physicalGrad); + return Math::Determinant(physicalGrad); } template diff --git a/Modules/IO/IOFDF/src/itkFDFImageIO.cxx b/Modules/IO/IOFDF/src/itkFDFImageIO.cxx index c66ff0cd3f8..c95758fe0db 100644 --- a/Modules/IO/IOFDF/src/itkFDFImageIO.cxx +++ b/Modules/IO/IOFDF/src/itkFDFImageIO.cxx @@ -17,6 +17,7 @@ *=========================================================================*/ #include "itkFDFImageIO.h" #include "itkFDFCommonImageIO.h" +#include "itkMathDeterminant.h" #include "itkByteSwapper.h" #include "itkRGBPixel.h" @@ -175,7 +176,7 @@ FDFImageIO::ReadImageInformation() // direction matrix in the file is 3x3. // if direction matrix is degenerate, punt and set // directions to identity - if (vnl_determinant(testDirections) == 0) + if (itk::Math::Determinant(testDirections) == 0) { for (unsigned int i = 0; i < numDim; i++) { diff --git a/Modules/IO/ImageBase/include/itkImageSeriesWriter.hxx b/Modules/IO/ImageBase/include/itkImageSeriesWriter.hxx index 87e6779991c..c5a29c299ed 100644 --- a/Modules/IO/ImageBase/include/itkImageSeriesWriter.hxx +++ b/Modules/IO/ImageBase/include/itkImageSeriesWriter.hxx @@ -25,7 +25,7 @@ #include "itkImageAlgorithm.h" #include "itkMetaDataObject.h" #include "itkArray.h" -#include "vnl/algo/vnl_determinant.h" +#include "itkMathDeterminant.h" #include namespace itk @@ -214,7 +214,7 @@ ImageSeriesWriter::WriteFiles() // therefore, replacing the orientation with an identity // is as arbitrary as any other choice. // - if (vnl_determinant(direction.GetVnlMatrix()) == 0.0) + if (Math::Determinant(direction.GetVnlMatrix()) == 0.0) { direction.SetIdentity(); } diff --git a/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h b/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h index 32917bcdc09..b736ac2d6ca 100644 --- a/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h +++ b/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h @@ -24,7 +24,6 @@ #include "itkCovariantVector.h" #include "vnl/vnl_matrix_fixed.h" #include "vnl/vnl_vector_fixed.h" -#include "vnl/vnl_det.h" #include "vnl/vnl_vector_fixed_ref.h" #include "vnl/vnl_vector.h" #include "itkTransform.h" diff --git a/Modules/Nonunit/Review/include/itkLabelGeometryImageFilter.h b/Modules/Nonunit/Review/include/itkLabelGeometryImageFilter.h index 24b85ce54a6..d33ab7377f3 100644 --- a/Modules/Nonunit/Review/include/itkLabelGeometryImageFilter.h +++ b/Modules/Nonunit/Review/include/itkLabelGeometryImageFilter.h @@ -25,7 +25,7 @@ #include #include #include "itkSymmetricEigenDecomposition.h" -#include "vnl/vnl_det.h" +#include "itkMathDeterminant.h" #include "itkMath.h" namespace itk diff --git a/Modules/Nonunit/Review/include/itkLabelGeometryImageFilter.hxx b/Modules/Nonunit/Review/include/itkLabelGeometryImageFilter.hxx index da5c6826800..32493b819d8 100644 --- a/Modules/Nonunit/Review/include/itkLabelGeometryImageFilter.hxx +++ b/Modules/Nonunit/Review/include/itkLabelGeometryImageFilter.hxx @@ -60,13 +60,9 @@ CalculateRotationMatrix(const itk::SymmetricEigenDecomposition & eig) // can fix this by making one of them negative. Make the last // eigenvector (with smallest eigenvalue) negative. float matrixDet = NAN; - if constexpr (VDimension == 2) + if constexpr (VDimension == 2 || VDimension == 3) { - matrixDet = vnl_det(rotationMatrix[0], rotationMatrix[1]); - } - else if constexpr (VDimension == 3) - { - matrixDet = vnl_det(rotationMatrix[0], rotationMatrix[1], rotationMatrix[2]); + matrixDet = static_cast(itk::Math::Determinant(rotationMatrix)); } else { diff --git a/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.h b/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.h index 4fec9f50ab2..987db761366 100644 --- a/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.h +++ b/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.h @@ -26,7 +26,7 @@ #if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE) # include "vnl/algo/vnl_matrix_inverse.h" // transitional transitive include; dropped on ITK legacy removal #endif -#include "vnl/algo/vnl_determinant.h" +#include "itkMathDeterminant.h" namespace itk { diff --git a/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.hxx b/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.hxx index ddded032df5..1eb4003ce3d 100644 --- a/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.hxx +++ b/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.hxx @@ -386,7 +386,7 @@ QuasiNewtonOptimizerv4Template::ComputeHessianAnd TInternalComputationValueType threshold = NumericTraits::epsilon(); - if (itk::Math::Absolute(vnl_determinant(newHessian)) <= threshold) + if (itk::Math::Absolute(Math::Determinant(newHessian)) <= threshold) { return false; } diff --git a/Modules/Numerics/Statistics/include/itkGaussianMembershipFunction.hxx b/Modules/Numerics/Statistics/include/itkGaussianMembershipFunction.hxx index aec0aa57871..ffc2dd10b4a 100644 --- a/Modules/Numerics/Statistics/include/itkGaussianMembershipFunction.hxx +++ b/Modules/Numerics/Statistics/include/itkGaussianMembershipFunction.hxx @@ -18,7 +18,7 @@ #ifndef itkGaussianMembershipFunction_hxx #define itkGaussianMembershipFunction_hxx -#include "vnl/algo/vnl_determinant.h" +#include "itkMathDeterminant.h" namespace itk::Statistics { @@ -104,11 +104,9 @@ GaussianMembershipFunction::SetCovariance(const CovarianceMa // the inverse of the covariance matrix is first computed by SVD const auto inv_cov = itk::Math::SVD(cov.GetVnlMatrix()); - // Signed determinant: the SVD's singular-value product is always - // non-negative and cannot detect a non-positive-definite matrix, so use - // vnl_determinant on the covariance. O(n^3) but n is the (tiny) - // measurement-vector dimension. - const double det = vnl_determinant(cov.GetVnlMatrix()); + // Signed determinant of the covariance: a magnitude-only determinant is + // always non-negative and cannot detect a non-positive-definite matrix. + const double det = itk::Math::Determinant(cov.GetVnlMatrix()); if (det <= 0.0) { diff --git a/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.h b/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.h index cc0471b1fb8..ca2e5f0bc0a 100644 --- a/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.h +++ b/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.h @@ -26,7 +26,6 @@ #if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE) # include "vnl/algo/vnl_matrix_inverse.h" // transitional transitive include; dropped on ITK legacy removal #endif -#include "vnl/algo/vnl_determinant.h" #include "itkArray.h" #include "itkDistanceMetric.h" diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_determinant.h b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_determinant.h index c1e6d0ff508..6baddc51b84 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_determinant.h +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_determinant.h @@ -1,6 +1,19 @@ // This is core/vnl/algo/vnl_determinant.h #ifndef vnl_algo_determinant_h_ #define vnl_algo_determinant_h_ + +#if __has_include() +# include +# if defined(ITK_FUTURE_LEGACY_REMOVE) +# error "vnl_determinant was removed; migrate to itk::Math::Determinant (itkMathDeterminant.h, Eigen-backed)." +# elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST) +# if defined(_MSC_VER) +# pragma message("vnl_determinant is deprecated; migrate to itk::Math::Determinant.") +# else +# warning "vnl_determinant is deprecated; migrate to itk::Math::Determinant." +# endif +# endif +#endif //: // \file // \brief calculates the determinant of a matrix diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/vnl_det.h b/Modules/ThirdParty/VNL/src/vxl/core/vnl/vnl_det.h index bd2751331c0..97c8fb01d06 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/vnl_det.h +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/vnl_det.h @@ -1,6 +1,19 @@ // This is core/vnl/vnl_det.h #ifndef vnl_det_h_ #define vnl_det_h_ + +#if __has_include() +# include +# if defined(ITK_FUTURE_LEGACY_REMOVE) +# error "vnl_det was removed; migrate to itk::Math::Determinant (itkMathDeterminant.h, Eigen-backed)." +# elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST) +# if defined(_MSC_VER) +# pragma message("vnl_det is deprecated; migrate to itk::Math::Determinant.") +# else +# warning "vnl_det is deprecated; migrate to itk::Math::Determinant." +# endif +# endif +#endif //: // \file // \brief Direct evaluation of 2x2, 3x3 and 4x4 determinants. diff --git a/Modules/Video/Core/include/itkVideoStream.hxx b/Modules/Video/Core/include/itkVideoStream.hxx index db69b988e67..aeb60e00906 100644 --- a/Modules/Video/Core/include/itkVideoStream.hxx +++ b/Modules/Video/Core/include/itkVideoStream.hxx @@ -18,6 +18,8 @@ #ifndef itkVideoStream_hxx #define itkVideoStream_hxx +#include "itkMathDeterminant.h" + namespace itk { @@ -174,7 +176,7 @@ void VideoStream::SetFrameDirection(SizeValueType frameNumber, typename TFrameType::DirectionType direction) { // Determinant is non-zero - if (itk::Math::Absolute(vnl_determinant(direction.GetVnlMatrix())) <= itk::Math::eps) + if (itk::Math::Absolute(Math::Determinant(direction.GetVnlMatrix())) <= itk::Math::eps) { itkExceptionMacro("Bad direction, determinant is 0. Direction is " << direction); }