From c8fedf00ca7af244d3a9b6ae15971fc1dae7e3ba Mon Sep 17 00:00:00 2001 From: "Hans J. Johnson" Date: Wed, 8 Jul 2026 07:26:33 -0500 Subject: [PATCH 1/4] ENH: Extend itk::Math::SVD for the vnl_svd consumer idioms Extend the Eigen-backed itk::Math::SVD helpers so every vnl_svd / vnl_matrix_inverse call site in ITK has a direct replacement: pseudo-inverse, solve, rank, recompose, null vector, and the WellCondition and DeterminantMagnitude predicates. Route the eigenvector sign convention through itkEigenDecompositionSignConvention so results are backend-stable. --- .../itkEigenDecompositionSignConvention.h | 11 +- Modules/Core/Common/include/itkMathSVD.h | 277 +++++++++++++++++- Modules/Core/Common/src/CMakeLists.txt | 1 + Modules/Core/Common/src/itkMathSVD.cxx | 50 ++++ 4 files changed, 320 insertions(+), 19 deletions(-) create mode 100644 Modules/Core/Common/src/itkMathSVD.cxx diff --git a/Modules/Core/Common/include/itkEigenDecompositionSignConvention.h b/Modules/Core/Common/include/itkEigenDecompositionSignConvention.h index 0acc983e192..c4d5c7158e7 100644 --- a/Modules/Core/Common/include/itkEigenDecompositionSignConvention.h +++ b/Modules/Core/Common/include/itkEigenDecompositionSignConvention.h @@ -54,12 +54,13 @@ CanonicalizeEigenvectorColumnSigns(vnl_matrix & V) /** Same canonicalization as \c CanonicalizeEigenvectorColumnSigns, applied to \a u, * with the identical per-column flip mirrored onto \a paired so a factor pair (e.g. * the U and V of an SVD) stays - * consistent. Templated on the matrix type so it serves vnl_matrix and - * vnl_matrix_fixed; the sign is well-defined only when the leading per-column + * consistent. Templated on the matrix types so it serves vnl_matrix and + * vnl_matrix_fixed, including pairs of distinct fixed types (a rectangular SVD's + * U and V); the sign is well-defined only when the leading per-column * magnitude is unambiguous (distinct singular values). */ -template +template void -CanonicalizeColumnSignsPaired(TMatrix & u, TMatrix & paired) +CanonicalizeColumnSignsPaired(TMatrixU & u, TMatrixV & paired) { // u and paired share a column count but may differ in row count (an SVD's // thin U is rows x k, V is cols x k), so each is flipped over its own rows. @@ -75,7 +76,7 @@ CanonicalizeColumnSignsPaired(TMatrix & u, TMatrix & paired) pivot = i; } } - if (u(pivot, j) < typename TMatrix::element_type{ 0 }) + if (u(pivot, j) < typename TMatrixU::element_type{ 0 }) { for (unsigned int i = 0; i < uRows; ++i) { diff --git a/Modules/Core/Common/include/itkMathSVD.h b/Modules/Core/Common/include/itkMathSVD.h index cbc0610186f..f320e24d413 100644 --- a/Modules/Core/Common/include/itkMathSVD.h +++ b/Modules/Core/Common/include/itkMathSVD.h @@ -19,7 +19,6 @@ #define itkMathSVD_h #include "itkMacro.h" -#include "itkMatrix.h" #include "itkEigenDecompositionSignConvention.h" #include "vnl/vnl_matrix.h" #include "vnl/vnl_matrix_fixed.h" @@ -30,9 +29,16 @@ #include ITK_EIGEN(Dense) #include +#include namespace itk { +// Forward-declared to break the itkMatrix.h <-> itkMathSVD.h include cycle; the +// itk::Matrix SVD overload needs the complete type only at instantiation, where +// the caller has already included itkMatrix.h. +template +class Matrix; + namespace Math { @@ -47,15 +53,15 @@ ResolveRcond(TReal rcond, unsigned int n) } // Moore-Penrose pseudo-inverse V diag(1/w) U^T, with singular values at or below -// rcond*max(w) treated as zero. Works for vnl_matrix and vnl_matrix_fixed. -template -TMatrix -PseudoInverse(const TMatrix & U, const TVector & W, const TMatrix & V, TReal rcond) +// rcond*max(w) treated as zero. U and V may be distinct types (rectangular fixed). +template +auto +PseudoInverse(const TMatrixU & U, const TVector & W, const TMatrixV & V, TReal rcond) { const unsigned int k = W.size(); // W is sorted descending (Eigen guarantee), so W[0] is the max without an O(k) scan. const TReal tol = ResolveRcond(rcond, k) * W[0]; - TMatrix scaledV = V; + TMatrixV scaledV = V; for (unsigned int col = 0; col < k; ++col) { const TReal s = (W[col] > tol) ? TReal{ 1 } / W[col] : TReal{ 0 }; @@ -68,13 +74,13 @@ PseudoInverse(const TMatrix & U, const TVector & W, const TMatrix & V, TReal rco } // Least-squares / minimum-norm solution x = V diag(1/w) U^T b of A x = b. -template -TVector -SolveLinear(const TMatrix & U, const TVector & W, const TMatrix & V, const TVector & b, TReal rcond) +template +auto +SolveLinear(const TMatrixU & U, const TVector & W, const TMatrixV & V, const TVectorB & b, TReal rcond) { const unsigned int n = W.size(); const TReal tol = ResolveRcond(rcond, n) * W[0]; // W sorted descending; W[0] == max - TVector utb = U.transpose() * b; + auto utb = U.transpose() * b; for (unsigned int k = 0; k < n; ++k) { utb[k] = (W[k] > tol) ? utb[k] / W[k] : TReal{ 0 }; @@ -101,13 +107,13 @@ NumericalRank(const TVector & W, TReal rcond) // Reconstruct U diag(w) V^T, treating singular values at or below rcond*max(w) as // zero (a truncated, rank-reduced reconstruction of A). -template -TMatrix -Recompose(const TMatrix & U, const TVector & W, const TMatrix & V, TReal rcond) +template +auto +Recompose(const TMatrixU & U, const TVector & W, const TMatrixV & V, TReal rcond) { const unsigned int k = W.size(); const TReal tol = ResolveRcond(rcond, k) * W[0]; // W sorted descending; W[0] == max - TMatrix scaledU = U; + TMatrixU scaledU = U; for (unsigned int col = 0; col < k; ++col) { const TReal s = (W[col] > tol) ? W[col] : TReal{ 0 }; @@ -118,6 +124,49 @@ Recompose(const TMatrix & U, const TVector & W, const TMatrix & V, TReal rcond) } return scaledU * V.transpose(); } + +// Reconstruct U diag(w') V^T with caller-supplied singular values (no truncation); +// serves the modify-W-then-recompose idiom (value clamping/inversion/zeroing). +template +auto +RecomposeWith(const TMatrixU & U, const TVector & modifiedW, const TMatrixV & V) +{ + TMatrixU scaledU = U; + for (unsigned int col = 0; col < modifiedW.size(); ++col) + { + for (unsigned int i = 0; i < scaledU.rows(); ++i) + { + scaledU(i, col) *= modifiedW[col]; + } + } + return scaledU * V.transpose(); +} + +// Reciprocal condition number sigma_min/sigma_max; 0 when the matrix is singular. +template +auto +WellCondition(const TVector & W) -> std::decay_t +{ + const unsigned int k = W.size(); + if (k == 0 || W[0] == 0) + { + return 0; + } + return W[k - 1] / W[0]; // W sorted descending +} + +// Product of the singular values (|det A| for a square A). +template +auto +DeterminantMagnitude(const TVector & W) -> std::decay_t +{ + std::decay_t product{ 1 }; + for (unsigned int k = 0; k < W.size(); ++k) + { + product *= W[k]; + } + return product; +} } // namespace detail /** Result of a fixed-size square SVD: A == U * diag(W) * V^T, W descending. @@ -156,6 +205,35 @@ struct FixedSquareSVDResult { return detail::Recompose(U, W, V, rcond); } + + /** Reconstruct U diag(modifiedW) V^T with caller-supplied singular values. */ + vnl_matrix_fixed + RecomposeWith(const vnl_vector_fixed & modifiedW) const + { + return detail::RecomposeWith(U, modifiedW, V); + } + + /** Singular vector for the smallest singular value (last column of V); spans + * the nullspace when A is rank-deficient. */ + vnl_vector_fixed + NullVector() const + { + return V.get_column(VDim - 1); + } + + /** Reciprocal condition number sigma_min/sigma_max (0 = singular, 1 = perfect). */ + TReal + WellCondition() const + { + return detail::WellCondition(W); + } + + /** Product of the singular values (magnitude of the determinant for a square A). */ + TReal + DeterminantMagnitude() const + { + return detail::DeterminantMagnitude(W); + } }; /** Result of a runtime-sized SVD (square or rectangular): A == U diag(W) V^T, @@ -195,6 +273,114 @@ struct SVDResult { return detail::Recompose(U, W, V, rcond); } + + /** Reconstruct U diag(modifiedW) V^T with caller-supplied singular values. */ + vnl_matrix + RecomposeWith(const vnl_vector & modifiedW) const + { + return detail::RecomposeWith(U, modifiedW, V); + } + + /** Singular vector for the smallest singular value (last column of V); spans + * the nullspace when A is rank-deficient. Defined only for rows >= cols: an + * underdetermined input's thin V (cols x min(rows,cols)) cannot span the + * nullspace, so this throws rather than return a misleading direction. */ + vnl_vector + NullVector() const + { + if (U.rows() < V.rows()) + { + itkGenericExceptionMacro( + "NullVector() requires rows >= cols; the thin V of an underdetermined input does not span the nullspace."); + } + return V.get_column(V.cols() - 1); + } + + /** Reciprocal condition number sigma_min/sigma_max (0 = singular, 1 = perfect). */ + TReal + WellCondition() const + { + return detail::WellCondition(W); + } + + /** Product of the singular values (magnitude of the determinant for a square A). */ + TReal + DeterminantMagnitude() const + { + return detail::DeterminantMagnitude(W); + } +}; + +/** Result of a fixed-size rectangular SVD with thin factors: A (VRows x VCols) + * == U diag(W) V^T, where U is VRows x K, V is VCols x K and K = min(VRows, + * VCols); W is descending. \a rcond < 0 auto-selects a K*epsilon threshold. */ +template +struct FixedRectangularSVDResult +{ + static constexpr unsigned int K = (VRows < VCols) ? VRows : VCols; + + vnl_matrix_fixed U{}; + vnl_vector_fixed W{}; + vnl_matrix_fixed V{}; + + /** Moore-Penrose pseudo-inverse A^+ (VCols x VRows). */ + vnl_matrix_fixed + PseudoInverse(TReal rcond = TReal{ -1 }) const + { + return detail::PseudoInverse(U, W, V, rcond); + } + + /** Least-squares / minimum-norm solution of A x = b. */ + vnl_vector_fixed + Solve(const vnl_vector_fixed & b, TReal rcond = TReal{ -1 }) const + { + return detail::SolveLinear(U, W, V, b, rcond); + } + + /** Numerical rank (count of singular values above rcond*max(w)). */ + unsigned int + Rank(TReal rcond = TReal{ -1 }) const + { + return detail::NumericalRank(W, rcond); + } + + /** Reconstruct A with singular values at or below rcond*max(w) zeroed. */ + vnl_matrix_fixed + Recompose(TReal rcond = TReal{ -1 }) const + { + return detail::Recompose(U, W, V, rcond); + } + + /** Reconstruct U diag(modifiedW) V^T with caller-supplied singular values. */ + vnl_matrix_fixed + RecomposeWith(const vnl_vector_fixed & modifiedW) const + { + return detail::RecomposeWith(U, modifiedW, V); + } + + /** Singular vector for the smallest singular value (last column of V); spans + * the nullspace when A is rank-deficient. Overdetermined shapes only: the thin + * V of an underdetermined input does not span the nullspace. */ + vnl_vector_fixed + NullVector() const + { + static_assert(VRows >= VCols, "NullVector() requires VRows >= VCols (thin V cannot span the nullspace)."); + return V.get_column(K - 1); + } + + /** Reciprocal condition number sigma_min/sigma_max (0 = singular, 1 = perfect). */ + TReal + WellCondition() const + { + return detail::WellCondition(W); + } + + /** Product of the singular values (magnitude of the determinant for a square A). */ + TReal + DeterminantMagnitude() const + { + return detail::DeterminantMagnitude(W); + } }; namespace detail @@ -315,6 +501,44 @@ RectangularSVDEigen(const TReal * inData, } } // namespace detail +// The Eigen JacobiSVD/BDCSVD instantiations behind the detail engines dominate +// compile memory: without pre-instantiation, GCC needs ~2 GB per translation unit +// that inverts an itk::Matrix. Declare the common specializations extern and define +// them once in itkMathSVD.cxx so consumers link instead of re-instantiating them. +#if defined(ITKCommon_EXPORTS) +# define ITKCommon_EXPORT_EXPLICIT ITK_TEMPLATE_EXPORT +#else +# define ITKCommon_EXPORT_EXPLICIT ITKCommon_EXPORT +#endif + +#define ITK_MATH_SVD_FIXED_DIMS(F) F(1) F(2) F(3) F(4) F(5) F(6) + +namespace detail +{ +ITK_GCC_PRAGMA_DIAG_PUSH() +ITK_GCC_PRAGMA_DIAG(ignored "-Wattributes") + +#define ITK_MATH_SVD_EXTERN_FIXED(D) \ + extern template ITKCommon_EXPORT_EXPLICIT void SquareSVDEigen(const float *, float *, float *, float *); \ + extern template ITKCommon_EXPORT_EXPLICIT void SquareSVDEigen( \ + const double *, double *, double *, double *); +ITK_MATH_SVD_FIXED_DIMS(ITK_MATH_SVD_EXTERN_FIXED) +#undef ITK_MATH_SVD_EXTERN_FIXED + +extern template ITKCommon_EXPORT_EXPLICIT void +DynamicSquareSVDEigen(const float *, unsigned int, float *, float *, float *); +extern template ITKCommon_EXPORT_EXPLICIT void +DynamicSquareSVDEigen(const double *, unsigned int, double *, double *, double *); +extern template ITKCommon_EXPORT_EXPLICIT void +RectangularSVDEigen(const float *, unsigned int, unsigned int, float *, float *, float *); +extern template ITKCommon_EXPORT_EXPLICIT void +RectangularSVDEigen(const double *, unsigned int, unsigned int, double *, double *, double *); + +ITK_GCC_PRAGMA_DIAG_POP() +} // namespace detail + +#undef ITKCommon_EXPORT_EXPLICIT + /** \brief Singular value decomposition A = U diag(W) V^T, backed by Eigen. * * Opt-in Eigen-backed alternative to vnl_svd: the optimized Eigen path is @@ -366,6 +590,31 @@ SVD(const Matrix & A, bool canonicalizeSigns = true) return SVD(A.GetVnlMatrix(), canonicalizeSigns); } +/** SVD of a fixed-size rectangular vnl_matrix_fixed (thin U/V factors); serves + * fixed non-square consumers such as an OutputDim x InputDim transform Jacobian. + * Square fixed inputs take the dedicated square overload. */ +template > +FixedRectangularSVDResult +SVD(const vnl_matrix_fixed & A, bool canonicalizeSigns = true) +{ + FixedRectangularSVDResult result; + detail::RectangularSVDEigen( + A.data_block(), VRows, VCols, result.U.data_block(), result.W.data_block(), result.V.data_block()); + if (canonicalizeSigns) + { + itk::detail::CanonicalizeColumnSignsPaired(result.U, result.V); + } + return result; +} + +/** SVD of a fixed-size rectangular itk::Matrix (thin U/V factors). */ +template > +FixedRectangularSVDResult +SVD(const Matrix & A, bool canonicalizeSigns = true) +{ + return SVD(A.GetVnlMatrix(), canonicalizeSigns); +} + /** SVD of a runtime-sized vnl_matrix (square or rectangular). Square inputs take * the fast square path (JacobiSVD+NoQRPreconditioner / BDCSVD); rectangular inputs * use BDCSVD with thin U/V. The dispatch is a single dimension comparison. */ diff --git a/Modules/Core/Common/src/CMakeLists.txt b/Modules/Core/Common/src/CMakeLists.txt index 1de6835aba2..d7da6aa50ba 100644 --- a/Modules/Core/Common/src/CMakeLists.txt +++ b/Modules/Core/Common/src/CMakeLists.txt @@ -95,6 +95,7 @@ set( itkLoggerOutput.cxx itkLoggerThreadWrapper.cxx itkLogOutput.cxx + itkMathSVD.cxx itkMemoryProbe.cxx itkMemoryProbesCollectorBase.cxx itkMemoryUsageObserver.cxx diff --git a/Modules/Core/Common/src/itkMathSVD.cxx b/Modules/Core/Common/src/itkMathSVD.cxx new file mode 100644 index 00000000000..4e69de406f0 --- /dev/null +++ b/Modules/Core/Common/src/itkMathSVD.cxx @@ -0,0 +1,50 @@ +/*========================================================================= + * + * 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 "itkMathSVD.h" + +// Define once here the Eigen JacobiSVD/BDCSVD specializations that itkMathSVD.h +// declares extern, so consumer translation units link instead of paying the +// per-TU Eigen instantiation cost. +namespace itk +{ +namespace Math +{ +namespace detail +{ +ITK_GCC_PRAGMA_DIAG_PUSH() +ITK_GCC_PRAGMA_DIAG(ignored "-Wattributes") + +#define ITK_MATH_SVD_INST_FIXED(D) \ + template ITKCommon_EXPORT void SquareSVDEigen(const float *, float *, float *, float *); \ + template ITKCommon_EXPORT void SquareSVDEigen(const double *, double *, double *, double *); +ITK_MATH_SVD_FIXED_DIMS(ITK_MATH_SVD_INST_FIXED) +#undef ITK_MATH_SVD_INST_FIXED + +template ITKCommon_EXPORT void +DynamicSquareSVDEigen(const float *, unsigned int, float *, float *, float *); +template ITKCommon_EXPORT void +DynamicSquareSVDEigen(const double *, unsigned int, double *, double *, double *); +template ITKCommon_EXPORT void +RectangularSVDEigen(const float *, unsigned int, unsigned int, float *, float *, float *); +template ITKCommon_EXPORT void +RectangularSVDEigen(const double *, unsigned int, unsigned int, double *, double *, double *); + +ITK_GCC_PRAGMA_DIAG_POP() +} // namespace detail +} // namespace Math +} // namespace itk From ea3249dff5c8f67ae16b238454ba226b30cbce1c Mon Sep 17 00:00:00 2001 From: "Hans J. Johnson" Date: Wed, 8 Jul 2026 07:26:36 -0500 Subject: [PATCH 2/4] ENH: Test itk::Math::SVD against vnl across the consumer idioms Exercise the extended itk::Math::SVD surface (pseudo-inverse, solve, rank, WellCondition, DeterminantMagnitude) and cross-check it against the vnl_svd result on well- and ill-conditioned inputs. --- Modules/Core/Common/test/itkMathSVDGTest.cxx | 172 ++++++++++++++++-- Modules/Core/Common/test/itkNumericsGTest.cxx | 10 +- .../test/itkSymmetricEigenAnalysisTest.cxx | 1 + 3 files changed, 170 insertions(+), 13 deletions(-) diff --git a/Modules/Core/Common/test/itkMathSVDGTest.cxx b/Modules/Core/Common/test/itkMathSVDGTest.cxx index 8a51b2882e8..055e88612c5 100644 --- a/Modules/Core/Common/test/itkMathSVDGTest.cxx +++ b/Modules/Core/Common/test/itkMathSVDGTest.cxx @@ -17,12 +17,16 @@ *=========================================================================*/ // First include the header file to be tested: +#include "itkMatrix.h" #include "itkMathSVD.h" // This test cross-checks against the vnl_svd reference engine; mark it so it keeps -// compiling once vnl_svd is deprecated under ITK_LEGACY_REMOVE. -#define ITK_LEGACY_TEST -#include "vnl/algo/vnl_svd.h" +// compiling once vnl_svd is deprecated under ITK_LEGACY_REMOVE. The oracle is +// unavailable under ITK_FUTURE_LEGACY_REMOVE. +#ifndef ITK_FUTURE_LEGACY_REMOVE +# define ITK_LEGACY_TEST +# include "vnl/algo/vnl_svd.h" +#endif #include #include @@ -76,6 +80,7 @@ TEST(MathSVD, FixedReconstructs) } } +#ifndef ITK_FUTURE_LEGACY_REMOVE // Singular values agree with vnl_svd (the engine being supplemented). TEST(MathSVD, SingularValuesMatchVnl) { @@ -85,6 +90,7 @@ TEST(MathSVD, SingularValuesMatchVnl) for (unsigned int i = 0; i < 6; ++i) EXPECT_NEAR(r.W[i], ref.W()(i, i), 1e-12); } +#endif // ITK_FUTURE_LEGACY_REMOVE // Default canonicalization makes the largest-magnitude element of each U column // positive, giving a deterministic, build-independent sign convention. @@ -106,6 +112,7 @@ TEST(MathSVD, SignsCanonical) } } +#ifndef ITK_FUTURE_LEGACY_REMOVE // U V^T is invariant to the SVD sign/basis ambiguity: it matches vnl_svd even // for a degenerate (repeated singular value) matrix. This is the property the // geometry call sites (orthogonalization) rely on. @@ -143,6 +150,7 @@ TEST(MathSVD, UVtransposeInvariantOnDegenerate) diff = std::max(diff, std::abs(rotEig(i, j) - rotVnl(i, j))); EXPECT_LT(diff, 1e-12); } +#endif // ITK_FUTURE_LEGACY_REMOVE // The itk::Matrix overload forwards to the same computation. TEST(MathSVD, ItkMatrixOverload) @@ -186,9 +194,11 @@ TEST(MathSVD, DynamicReconstructs) } EXPECT_LT(err, 1e-11) << "n=" << n; +#ifndef ITK_FUTURE_LEGACY_REMOVE const vnl_svd ref(A); for (unsigned int i = 0; i < n; ++i) EXPECT_NEAR(r.W[i], ref.W()(i, i), 1e-11) << "n=" << n << " i=" << i; +#endif } } @@ -202,14 +212,16 @@ TEST(MathSVD, DynamicRejectsEmpty) // PseudoInverse() agrees with vnl_svd and satisfies the Moore-Penrose identity. TEST(MathSVD, PseudoInverseMatchesVnl) { - const auto A = MakeFixed(); - const auto pinv = itk::Math::SVD(A).PseudoInverse(); + const auto A = MakeFixed(); + const auto pinv = itk::Math::SVD(A).PseudoInverse(); +#ifndef ITK_FUTURE_LEGACY_REMOVE const vnl_matrix pinvVnl = vnl_svd(A.as_matrix()).pinverse(); double dInv = 0.0; for (unsigned int i = 0; i < 6; ++i) for (unsigned int j = 0; j < 6; ++j) dInv = std::max(dInv, std::abs(pinv(i, j) - pinvVnl(i, j))); EXPECT_LT(dInv, 1e-10); +#endif // A * A^+ * A == A const vnl_matrix_fixed recon = A * pinv * A; @@ -234,11 +246,33 @@ TEST(MathSVD, SolveMatchesVnl) const vnl_vector_fixed residual = A * x - b; EXPECT_LT(residual.inf_norm(), 1e-10); +#ifndef ITK_FUTURE_LEGACY_REMOVE // agrees with vnl_svd's solve const vnl_vector xVnl = vnl_svd(A.as_matrix()).solve(b.as_vector()); EXPECT_LT((x.as_vector() - xVnl).inf_norm(), 1e-10); +#endif } +#ifndef ITK_FUTURE_LEGACY_REMOVE +// WellCondition() (sigma_min/sigma_max) agrees with vnl_svd::well_condition(), +// the accessor the NIfTI direction-cosine check migrated onto. +TEST(MathSVD, WellConditionMatchesVnl) +{ + const auto A = MakeFixed(); + EXPECT_NEAR(itk::Math::SVD(A).WellCondition(), vnl_svd(A.as_matrix()).well_condition(), 1e-12); +} + +// DeterminantMagnitude() (product of singular values) agrees with +// vnl_svd::determinant_magnitude(), the accessor the Mahalanobis check migrated onto. +TEST(MathSVD, DeterminantMagnitudeMatchesVnl) +{ + const auto A = MakeFixed(); + const double itkDet = itk::Math::SVD(A).DeterminantMagnitude(); + const double vnlDet = vnl_svd(A.as_matrix()).determinant_magnitude(); + EXPECT_NEAR(itkDet, vnlDet, 1e-10 * std::abs(vnlDet)); +} +#endif // ITK_FUTURE_LEGACY_REMOVE + // rank() reports full rank for a well-conditioned matrix and the reduced rank of // a deliberately rank-deficient one. TEST(MathSVD, Rank) @@ -270,13 +304,15 @@ TEST(MathSVD, IllConditionedMatchesVnl) for (unsigned int j = 0; j < N; ++j) A(i, j) = 1.0 / (i + j + 1.0); // Hilbert - const auto r = itk::Math::SVD(A); - const vnl_svd ref(A); + const auto r = itk::Math::SVD(A); +#ifndef ITK_FUTURE_LEGACY_REMOVE // singular values agree relative to the largest - const double wmax = r.W[0]; + const vnl_svd ref(A); + const double wmax = r.W[0]; for (unsigned int i = 0; i < N; ++i) EXPECT_NEAR(r.W[i], ref.W()(i, i), 1e-9 * wmax) << "i=" << i; +#endif // reconstruction A == U diag(W) V^T is well-conditioned and stays ~machine eps double err = 0.0; @@ -439,15 +475,17 @@ TEST(MathSVD, Rectangular) EXPECT_LT(reconErr, 1e-12) << "m=" << m << " n=" << n; // PseudoInverse is n x m and matches vnl_svd - const auto pinv = r.PseudoInverse(); - const vnl_matrix pinvVnl = vnl_svd(A).pinverse(); + const auto pinv = r.PseudoInverse(); EXPECT_EQ(pinv.rows(), n); EXPECT_EQ(pinv.cols(), m); - double pinvErr = 0.0; +#ifndef ITK_FUTURE_LEGACY_REMOVE + const vnl_matrix pinvVnl = vnl_svd(A).pinverse(); + double pinvErr = 0.0; for (unsigned int i = 0; i < n; ++i) for (unsigned int j = 0; j < m; ++j) pinvErr = std::max(pinvErr, std::abs(pinv(i, j) - pinvVnl(i, j))); EXPECT_LT(pinvErr, 1e-10) << "m=" << m << " n=" << n; +#endif // Moore-Penrose: A A^+ A == A const vnl_matrix recon = A * pinv * A; @@ -471,7 +509,117 @@ TEST(MathSVD, RectangularSolveMatchesVnl) b[i] = static_cast(i) - 1.5; const vnl_vector x = itk::Math::SVD(A).Solve(b); - const vnl_vector xVnl = vnl_svd(A).solve(b); EXPECT_EQ(x.size(), 3u); +#ifndef ITK_FUTURE_LEGACY_REMOVE + const vnl_vector xVnl = vnl_svd(A).solve(b); EXPECT_LT((x - xVnl).inf_norm(), 1e-10); +#endif +} + +// NullVector() matches vnl_svd::nullvector() up to sign, and lies in the nullspace. +TEST(MathSVD, NullVector) +{ + // Rank-2 square 3x3: third row = row0 + row1. + vnl_matrix A(3, 3); + A(0, 0) = 1.0; + A(0, 1) = 2.0; + A(0, 2) = 3.0; + A(1, 0) = 4.0; + A(1, 1) = 5.0; + A(1, 2) = 6.0; + for (unsigned int j = 0; j < 3; ++j) + A(2, j) = A(0, j) + A(1, j); + + const vnl_vector nv = itk::Math::SVD(A).NullVector(); + EXPECT_LT((A * nv).inf_norm(), 1e-12); + +#ifndef ITK_FUTURE_LEGACY_REMOVE + const vnl_vector nvVnl = vnl_svd(A).nullvector(); + const double align = std::abs(dot_product(nv, nvVnl)) / (nv.magnitude() * nvVnl.magnitude()); + EXPECT_NEAR(align, 1.0, 1e-12); +#endif + + // Fixed-size overload agrees. + vnl_matrix_fixed Af; + Af.copy_in(A.data_block()); + const vnl_vector_fixed nvf = itk::Math::SVD(Af).NullVector(); + EXPECT_LT((A * nvf.as_vector()).inf_norm(), 1e-12); +} + +// NullVector() refuses an underdetermined input whose thin V cannot span the nullspace. +TEST(MathSVD, NullVectorRejectsUnderdetermined) +{ + vnl_matrix A(2, 4); + A.fill(1.0); + A(0, 1) = 2.0; + A(1, 3) = 3.0; + EXPECT_THROW(itk::Math::SVD(A).NullVector(), itk::ExceptionObject); +} + +// RecomposeWith() reproduces A for unmodified W and honors hand-edited values. +TEST(MathSVD, RecomposeWith) +{ + const auto Af = MakeFixed(); + const auto r = itk::Math::SVD(Af); + + const vnl_matrix_fixed same = r.RecomposeWith(r.W); + double err = 0.0; + for (unsigned int i = 0; i < 4; ++i) + for (unsigned int j = 0; j < 4; ++j) + err = std::max(err, std::abs(same(i, j) - Af(i, j))); + EXPECT_LT(err, 1e-12); + + // Zeroing the tail value reproduces the rank-truncated Recompose(). + vnl_vector_fixed wMod = r.W; + wMod[3] = 0.0; + const vnl_matrix_fixed truncated = r.RecomposeWith(wMod); + // Truncate via rcond just above the smallest normalized singular value. + const double rcond = (r.W[3] / r.W[0]) * (1.0 + 1e-6); + const vnl_matrix_fixed viaRcond = r.Recompose(rcond); + for (unsigned int i = 0; i < 4; ++i) + for (unsigned int j = 0; j < 4; ++j) + EXPECT_NEAR(truncated(i, j), viaRcond(i, j), 1e-12); + + // Inverted values give the transposed pseudo-inverse: U diag(1/w) V^T == (A^+)^T. + vnl_vector_fixed wInv; + for (unsigned int k = 0; k < 4; ++k) + wInv[k] = 1.0 / r.W[k]; + const vnl_matrix_fixed recomposedInv = r.RecomposeWith(wInv); + const vnl_matrix_fixed pinvT = r.PseudoInverse(0.0).transpose(); + for (unsigned int i = 0; i < 4; ++i) + for (unsigned int j = 0; j < 4; ++j) + EXPECT_NEAR(recomposedInv(i, j), pinvT(i, j), 1e-10); +} + +// Fixed rectangular overload: factors match the dynamic path and the +// pseudo-inverse matches the vnl_svd reference (the itkTransform Jacobian use). +TEST(MathSVD, FixedRectangular) +{ + constexpr unsigned int rows = 5; + constexpr unsigned int cols = 3; + vnl_matrix_fixed Af; + for (unsigned int i = 0; i < rows; ++i) + for (unsigned int j = 0; j < cols; ++j) + Af(i, j) = std::cos(0.4 * i + 0.9 * j) + (i == j ? 2.0 : 0.0); + + const auto rf = itk::Math::SVD(Af); + const auto rd = itk::Math::SVD(vnl_matrix(Af.as_matrix())); + for (unsigned int k = 0; k < 3; ++k) + EXPECT_NEAR(rf.W[k], rd.W[k], 1e-12); + +#ifndef ITK_FUTURE_LEGACY_REMOVE + const vnl_matrix_fixed pinv = rf.PseudoInverse(); + const vnl_matrix pinvVnl = vnl_svd(Af.as_matrix()).inverse(); + for (unsigned int i = 0; i < cols; ++i) + for (unsigned int j = 0; j < rows; ++j) + EXPECT_NEAR(pinv(i, j), pinvVnl(i, j), 1e-10); +#endif + + // NullVector on an overdetermined fixed shape lies in the nullspace of a + // rank-deficient input (third column = sum of the first two). + vnl_matrix_fixed Adef = Af; + for (unsigned int i = 0; i < rows; ++i) + Adef(i, 2) = Adef(i, 0) + Adef(i, 1); + const vnl_vector_fixed nv = itk::Math::SVD(Adef).NullVector(); + EXPECT_LT((Adef.as_matrix() * nv.as_vector()).inf_norm(), 1e-12); } diff --git a/Modules/Core/Common/test/itkNumericsGTest.cxx b/Modules/Core/Common/test/itkNumericsGTest.cxx index 15e662f633c..a3814ddb90e 100644 --- a/Modules/Core/Common/test/itkNumericsGTest.cxx +++ b/Modules/Core/Common/test/itkNumericsGTest.cxx @@ -15,10 +15,16 @@ * limitations under the License. * *=========================================================================*/ +#define ITK_LEGACY_TEST #include "itkGTest.h" +#include "itkConfigure.h" #include -#include "vnl/algo/vnl_svd.h" +// This converted legacy test exercises vnl_svd, which is unavailable under +// ITK_FUTURE_LEGACY_REMOVE. +#ifndef ITK_FUTURE_LEGACY_REMOVE + +# include "vnl/algo/vnl_svd.h" namespace { @@ -77,3 +83,5 @@ TEST(Numerics, ConvertedLegacyTest) EXPECT_NEAR(result(1, 0), 7.0 / 3.0, 1e-6); EXPECT_NEAR(result(2, 0), -7.0 / 6.0, 1e-6); } + +#endif // ITK_FUTURE_LEGACY_REMOVE diff --git a/Modules/Core/Common/test/itkSymmetricEigenAnalysisTest.cxx b/Modules/Core/Common/test/itkSymmetricEigenAnalysisTest.cxx index 52b4fb57013..72fa8a87a65 100644 --- a/Modules/Core/Common/test/itkSymmetricEigenAnalysisTest.cxx +++ b/Modules/Core/Common/test/itkSymmetricEigenAnalysisTest.cxx @@ -24,6 +24,7 @@ #include #include #include "itkTestingMacros.h" +#include "vnl/vnl_diag_matrix.h" // Test template instantiations for various supported template arguments: From fc7ab6bb55033821c10f33395e7977c57c4ea58d Mon Sep 17 00:00:00 2001 From: "Hans J. Johnson" Date: Wed, 8 Jul 2026 07:27:45 -0500 Subject: [PATCH 3/4] ENH: Migrate ITK consumers off vnl_svd/vnl_matrix_inverse to itk::Math::SVD Replace direct vnl_svd, vnl_svd_fixed and vnl_matrix_inverse use across the transforms, diffusion-tensor reconstruction, image moments, label-map filters, membership functions, Mahalanobis metrics, the Quasi-Newton optimizer, NIfTI and MINC IO and the classifier estimators with the Eigen-backed itk::Math::SVD API, including itk::Matrix::GetInverse. Where a widely included header (itkTransform, itkMatrix, the affine transforms) previously pulled a vnl_svd* header transitively, keep that include reachable under #if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE) so downstream code that relied on the transitive include still compiles during the deprecation window. --- .../RegistrationITKv4/ImageRegistration9.cxx | 9 +++++---- Modules/Core/Common/include/itkMatrix.h | 17 ++++++++++++++--- .../Common/include/itkVariableSizeMatrix.h | 7 +++++-- .../Transform/include/itkAffineTransform.hxx | 4 +++- .../include/itkCenteredAffineTransform.hxx | 4 +++- ...itkFixedCenterOfRotationAffineTransform.hxx | 4 +++- .../Transform/include/itkKernelTransform.h | 4 +++- .../include/itkMatrixOffsetTransformBase.hxx | 4 +++- .../include/itkScalableAffineTransform.hxx | 4 +++- .../Core/Transform/include/itkTransform.hxx | 18 +++++++++++++----- .../Transform/test/itkRigid3DTransformTest.cxx | 4 ++-- ...iffusionTensor3DReconstructionImageFilter.h | 14 ++++---------- ...fusionTensor3DReconstructionImageFilter.hxx | 3 ++- .../include/itkImageMomentsCalculator.hxx | 1 + .../include/itkImagePCAShapeModelEstimator.h | 4 +++- .../include/itkShapeLabelMapFilter.hxx | 1 + .../include/itkStatisticsLabelMapFilter.hxx | 1 + Modules/IO/MINC/src/itkMINCImageIO.cxx | 3 ++- Modules/IO/NIFTI/src/itkNiftiImageIO.cxx | 9 ++++----- .../include/itkQuasiNewtonOptimizerv4.h | 5 ++++- .../include/itkQuasiNewtonOptimizerv4.hxx | 4 +--- .../include/itkGaussianMembershipFunction.hxx | 18 +++++++----------- ...tkMahalanobisDistanceMembershipFunction.hxx | 11 +++++++---- .../include/itkMahalanobisDistanceMetric.h | 5 ++++- .../include/itkMahalanobisDistanceMetric.hxx | 2 +- .../RegistrationITKv3/ImageRegistration9.cxx | 13 +++++++------ .../include/itkImageGaussianModelEstimator.h | 4 +++- .../include/itkImageKmeansModelEstimator.h | 4 +++- 28 files changed, 112 insertions(+), 69 deletions(-) diff --git a/Examples/RegistrationITKv4/ImageRegistration9.cxx b/Examples/RegistrationITKv4/ImageRegistration9.cxx index 988e75dc6b9..83fcf1828eb 100644 --- a/Examples/RegistrationITKv4/ImageRegistration9.cxx +++ b/Examples/RegistrationITKv4/ImageRegistration9.cxx @@ -73,6 +73,7 @@ // that will monitor the evolution of the registration process. // #include "itkCommand.h" +#include "itkMathSVD.h" class CommandIterationUpdate : public itk::Command { public: @@ -112,9 +113,9 @@ class CommandIterationUpdate : public itk::Command p[0][1] = static_cast(optimizer->GetCurrentPosition()[1]); p[1][0] = static_cast(optimizer->GetCurrentPosition()[2]); p[1][1] = static_cast(optimizer->GetCurrentPosition()[3]); - vnl_svd svd(p); + const auto svd = itk::Math::SVD(p); vnl_matrix r(2, 2); - r = svd.U() * vnl_transpose(svd.V()); + r = svd.U * svd.V.transpose(); const double angle = std::asin(r[1][0]); std::cout << " AffineAngle: " << angle * 180.0 / itk::Math::pi << std::endl; @@ -411,9 +412,9 @@ main(int argc, char * argv[]) p[0][1] = static_cast(finalParameters[1]); p[1][0] = static_cast(finalParameters[2]); p[1][1] = static_cast(finalParameters[3]); - vnl_svd svd(p); + const auto svd = itk::Math::SVD(p); vnl_matrix r(2, 2); - r = svd.U() * vnl_transpose(svd.V()); + r = svd.U * svd.V.transpose(); const double angle = std::asin(r[1][0]); const double angleInDegrees = angle * 180.0 / itk::Math::pi; diff --git a/Modules/Core/Common/include/itkMatrix.h b/Modules/Core/Common/include/itkMatrix.h index 570c879b4e7..84931fafee0 100644 --- a/Modules/Core/Common/include/itkMatrix.h +++ b/Modules/Core/Common/include/itkMatrix.h @@ -24,8 +24,20 @@ #include #include "vnl/vnl_matrix_fixed.hxx" // Get the templates #include "vnl/vnl_transpose.h" -#include "vnl/algo/vnl_matrix_inverse.h" +#include "itkMathSVD.h" #include "vnl/vnl_matrix.h" + +// GetInverse is Eigen-backed via itk::Math::SVD, but ITK 5.4 exposed +// vnl_matrix_inverse (and through it vnl_svd / vnl_diag_matrix) transitively to +// every itkMatrix.h consumer. Retain that surface for release-5.4 source +// compatibility in default builds; ITK proper includes what it uses, so +// ITK_LEGACY_REMOVE drops the leak. The transitional guard keeps the +// vnl_matrix_inverse deprecation warning off itkMatrix consumers. +#if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE) +# define ITK_VNL_SVD_TRANSITIONAL_INCLUDE +# include "vnl/algo/vnl_matrix_inverse.h" +# undef ITK_VNL_SVD_TRANSITIONAL_INCLUDE +#endif #include "vnl/algo/vnl_determinant.h" #include "itkMath.h" #include // For is_same. @@ -315,8 +327,7 @@ class ITK_TEMPLATE_EXPORT Matrix { itkGenericExceptionMacro("Singular matrix. Determinant is 0."); } - const vnl_matrix_inverse inverse(m_Matrix.as_ref()); - return vnl_matrix_fixed{ inverse.as_matrix() }; + return vnl_matrix_fixed{ Math::SVD(m_Matrix.as_ref()).PseudoInverse() }; } /** Return the transposed matrix. */ diff --git a/Modules/Core/Common/include/itkVariableSizeMatrix.h b/Modules/Core/Common/include/itkVariableSizeMatrix.h index 293b9e700c6..cc421acbb48 100644 --- a/Modules/Core/Common/include/itkVariableSizeMatrix.h +++ b/Modules/Core/Common/include/itkVariableSizeMatrix.h @@ -21,7 +21,10 @@ #include "itkPoint.h" #include "itkCovariantVector.h" #include "vnl/vnl_matrix_fixed.h" -#include "vnl/algo/vnl_matrix_inverse.h" +#include "itkMathSVD.h" +#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/vnl_transpose.h" #include "vnl/vnl_matrix.h" #include "itkArray.h" @@ -209,7 +212,7 @@ class ITK_TEMPLATE_EXPORT VariableSizeMatrix [[nodiscard]] inline vnl_matrix GetInverse() const { - return vnl_matrix_inverse(m_Matrix).as_matrix(); + return itk::Math::SVD(m_Matrix).PseudoInverse(); } /** Return the transposed matrix. */ diff --git a/Modules/Core/Transform/include/itkAffineTransform.hxx b/Modules/Core/Transform/include/itkAffineTransform.hxx index d31c52d02e0..195e151ddfa 100644 --- a/Modules/Core/Transform/include/itkAffineTransform.hxx +++ b/Modules/Core/Transform/include/itkAffineTransform.hxx @@ -19,7 +19,9 @@ #define itkAffineTransform_hxx #include "itkNumericTraits.h" -#include "vnl/algo/vnl_matrix_inverse.h" +#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 namespace itk { diff --git a/Modules/Core/Transform/include/itkCenteredAffineTransform.hxx b/Modules/Core/Transform/include/itkCenteredAffineTransform.hxx index 485bcbcd326..169e9215045 100644 --- a/Modules/Core/Transform/include/itkCenteredAffineTransform.hxx +++ b/Modules/Core/Transform/include/itkCenteredAffineTransform.hxx @@ -19,7 +19,9 @@ #define itkCenteredAffineTransform_hxx #include "itkNumericTraits.h" -#include "vnl/algo/vnl_matrix_inverse.h" +#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 namespace itk { diff --git a/Modules/Core/Transform/include/itkFixedCenterOfRotationAffineTransform.hxx b/Modules/Core/Transform/include/itkFixedCenterOfRotationAffineTransform.hxx index 381785f268a..c54ddaf5c39 100644 --- a/Modules/Core/Transform/include/itkFixedCenterOfRotationAffineTransform.hxx +++ b/Modules/Core/Transform/include/itkFixedCenterOfRotationAffineTransform.hxx @@ -19,8 +19,10 @@ #define itkFixedCenterOfRotationAffineTransform_hxx #include "itkNumericTraits.h" +#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 "itkAffineTransform.h" -#include "vnl/algo/vnl_matrix_inverse.h" namespace itk { diff --git a/Modules/Core/Transform/include/itkKernelTransform.h b/Modules/Core/Transform/include/itkKernelTransform.h index 30343348b3e..219b3588358 100644 --- a/Modules/Core/Transform/include/itkKernelTransform.h +++ b/Modules/Core/Transform/include/itkKernelTransform.h @@ -19,6 +19,9 @@ #define itkKernelTransform_h #include "itkTransform.h" +#if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE) +# include "vnl/algo/vnl_svd.h" // transitional transitive include; dropped on ITK legacy removal +#endif #include "itkMatrix.h" #include "itkPointSet.h" #include @@ -27,7 +30,6 @@ #include "vnl/vnl_matrix.h" #include "vnl/vnl_vector.h" #include "vnl/vnl_vector_fixed.h" -#include "vnl/algo/vnl_svd.h" namespace itk { diff --git a/Modules/Core/Transform/include/itkMatrixOffsetTransformBase.hxx b/Modules/Core/Transform/include/itkMatrixOffsetTransformBase.hxx index 4699e908c1d..ec6c449bdd3 100644 --- a/Modules/Core/Transform/include/itkMatrixOffsetTransformBase.hxx +++ b/Modules/Core/Transform/include/itkMatrixOffsetTransformBase.hxx @@ -19,7 +19,9 @@ #define itkMatrixOffsetTransformBase_hxx #include "itkNumericTraits.h" -#include "vnl/algo/vnl_matrix_inverse.h" +#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 "itkMath.h" #include "itkCrossHelper.h" diff --git a/Modules/Core/Transform/include/itkScalableAffineTransform.hxx b/Modules/Core/Transform/include/itkScalableAffineTransform.hxx index 5d3334dee01..31ab2004eaf 100644 --- a/Modules/Core/Transform/include/itkScalableAffineTransform.hxx +++ b/Modules/Core/Transform/include/itkScalableAffineTransform.hxx @@ -19,9 +19,11 @@ #define itkScalableAffineTransform_hxx #include "itkMath.h" +#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 "itkNumericTraits.h" #include "itkPrintHelper.h" -#include "vnl/algo/vnl_matrix_inverse.h" namespace itk { diff --git a/Modules/Core/Transform/include/itkTransform.hxx b/Modules/Core/Transform/include/itkTransform.hxx index beb1f7f837a..e2a1ac6cfb5 100644 --- a/Modules/Core/Transform/include/itkTransform.hxx +++ b/Modules/Core/Transform/include/itkTransform.hxx @@ -19,7 +19,10 @@ #define itkTransform_hxx #include "itkCrossHelper.h" -#include "vnl/algo/vnl_svd_fixed.h" +#include "itkMathSVD.h" +#if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE) +# include "vnl/algo/vnl_svd_fixed.h" // transitional transitive include; dropped on ITK legacy removal +#endif #include "itkMetaProgrammingLibrary.h" namespace itk @@ -532,10 +535,15 @@ Transform::ComputeInver JacobianPositionType forward_jacobian; this->ComputeJacobianWithRespectToPosition(pnt, forward_jacobian); - using SVDAlgorithmType = - vnl_svd_fixed; - const SVDAlgorithmType svd(forward_jacobian); - jacobian.set(svd.inverse().data_block()); + if constexpr (VOutputDimension == VInputDimension) + { + // Square: fixed-size pseudo-inverse avoids the dynamic as_matrix() round-trip. + jacobian = Math::SVD(forward_jacobian).PseudoInverse(); + } + else + { + jacobian.set(Math::SVD(forward_jacobian.as_matrix()).PseudoInverse().data_block()); + } } template diff --git a/Modules/Core/Transform/test/itkRigid3DTransformTest.cxx b/Modules/Core/Transform/test/itkRigid3DTransformTest.cxx index 9d0ccda584a..2e65b27fd16 100644 --- a/Modules/Core/Transform/test/itkRigid3DTransformTest.cxx +++ b/Modules/Core/Transform/test/itkRigid3DTransformTest.cxx @@ -709,7 +709,7 @@ itkRigid3DTransformTest(int, char *[]) } TransformType::MatrixType expectedInverse{}; expectedInverse.SetIdentity(); - ITK_TEST_EXPECT_EQUAL(inverse->GetMatrix(), expectedInverse); + ITK_TEST_EXPECT_TRUE(inverse->GetMatrix().GetVnlMatrix().is_equal(expectedInverse.GetVnlMatrix(), 1e-10)); // An orthogonal transform should have an inverse: use a rotation matrix transform->SetIdentity(); @@ -743,7 +743,7 @@ itkRigid3DTransformTest(int, char *[]) expectedInverse[0][1] = -sinth; expectedInverse[1][0] = sinth; expectedInverse[1][1] = costh; - ITK_TEST_EXPECT_EQUAL(inverse->GetMatrix(), expectedInverse); + ITK_TEST_EXPECT_TRUE(inverse->GetMatrix().GetVnlMatrix().is_equal(expectedInverse.GetVnlMatrix(), 1e-10)); // Cannot test a singular matrix (i.e. not having an inverse) since ITK does not allow to set a // non-orthogonal rotation matrix to the transform diff --git a/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.h b/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.h index 3bd5427de8b..4e6966dc72e 100644 --- a/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.h +++ b/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.h @@ -24,7 +24,10 @@ #include "vnl/vnl_matrix.h" #include "vnl/vnl_vector_fixed.h" #include "vnl/vnl_matrix_fixed.h" -#include "vnl/algo/vnl_svd.h" +#include "itkMathSVD.h" +#if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE) +# include "vnl/algo/vnl_svd.h" // transitional transitive include; dropped on ITK legacy removal +#endif #include "itkVectorContainer.h" #include "itkVectorImage.h" #include "ITKDiffusionTensorImageExport.h" @@ -106,15 +109,6 @@ operator<<(std::ostream & out, const DiffusionTensor3DReconstructionImageFilterE * * For additional details see \cite westin2002 and \cite westin2002a. * - * \warning - * Although this filter has been written to support multiple threads, please - * set the number of threads to 1. - \code - filter->SetNumberOfWorkUnits(1); - \endcode - * This is due to buggy code in netlib/dsvdc, that is called by vnl_svd. - * (used to compute the pseudo-inverse to find the dual tensor basis). - * * \author Xiaodong Tao, GE, for contributing parts of this class. * \author Casey Goodlet, UNC for patches to support multiple baseline * images and other improvements. diff --git a/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.hxx b/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.hxx index 010b4be8fbf..118147e497e 100644 --- a/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.hxx +++ b/Modules/Filtering/DiffusionTensorImage/include/itkDiffusionTensor3DReconstructionImageFilter.hxx @@ -25,6 +25,7 @@ #include "itkImageMaskSpatialObject.h" #include "vnl/vnl_vector.h" #include "itkTotalProgressReporter.h" +#include "itkMathSVD.h" namespace itk { @@ -435,7 +436,7 @@ DiffusionTensor3DReconstructionImageFilter{ m_TensorBasis.as_matrix() }.pinverse(); + m_TensorBasisInverse = itk::Math::SVD(m_TensorBasis.as_matrix()).PseudoInverse(); m_BMatrix.inplace_transpose(); } diff --git a/Modules/Filtering/ImageStatistics/include/itkImageMomentsCalculator.hxx b/Modules/Filtering/ImageStatistics/include/itkImageMomentsCalculator.hxx index 5dff4b62f6e..5b188aa2c49 100644 --- a/Modules/Filtering/ImageStatistics/include/itkImageMomentsCalculator.hxx +++ b/Modules/Filtering/ImageStatistics/include/itkImageMomentsCalculator.hxx @@ -21,6 +21,7 @@ #include "itkRealEigenDecomposition.h" #include "itkSymmetricEigenDecomposition.h" #include "itkImageRegionConstIteratorWithIndex.h" +#include "vnl/vnl_diag_matrix.h" namespace itk { diff --git a/Modules/Filtering/ImageStatistics/include/itkImagePCAShapeModelEstimator.h b/Modules/Filtering/ImageStatistics/include/itkImagePCAShapeModelEstimator.h index e28016916e5..15751ccd588 100644 --- a/Modules/Filtering/ImageStatistics/include/itkImagePCAShapeModelEstimator.h +++ b/Modules/Filtering/ImageStatistics/include/itkImagePCAShapeModelEstimator.h @@ -19,13 +19,15 @@ #define itkImagePCAShapeModelEstimator_h #include +#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 #include #include "vnl/vnl_vector.h" #include "vnl/vnl_matrix.h" #include "itkMath.h" -#include "vnl/algo/vnl_matrix_inverse.h" #include "itkImageRegionIterator.h" #include "itkMacro.h" diff --git a/Modules/Filtering/LabelMap/include/itkShapeLabelMapFilter.hxx b/Modules/Filtering/LabelMap/include/itkShapeLabelMapFilter.hxx index ef70df144bd..1ee8e14bdda 100644 --- a/Modules/Filtering/LabelMap/include/itkShapeLabelMapFilter.hxx +++ b/Modules/Filtering/LabelMap/include/itkShapeLabelMapFilter.hxx @@ -30,6 +30,7 @@ #include "itkMath.h" #include "itkLexicographicCompare.h" #include +#include "vnl/vnl_diag_matrix.h" #include namespace itk diff --git a/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx b/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx index 0e99d00ee2f..0fc63e69ac0 100644 --- a/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx +++ b/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx @@ -23,6 +23,7 @@ #include "itkProgressReporter.h" #include "itkRealEigenDecomposition.h" #include "itkSymmetricEigenDecomposition.h" +#include "vnl/vnl_diag_matrix.h" namespace itk { diff --git a/Modules/IO/MINC/src/itkMINCImageIO.cxx b/Modules/IO/MINC/src/itkMINCImageIO.cxx index 07c5c5bf0bb..667adab852c 100644 --- a/Modules/IO/MINC/src/itkMINCImageIO.cxx +++ b/Modules/IO/MINC/src/itkMINCImageIO.cxx @@ -22,6 +22,7 @@ #include "vnl/vnl_vector.h" #include "itkMetaDataObject.h" #include "itkArray.h" +#include "itkMathSVD.h" #include "itkPrintHelper.h" #include "itkMakeUniqueForOverwrite.h" @@ -983,7 +984,7 @@ MINCImageIO::WriteImageInformation() origin[i] = this->GetOrigin(i); } - const vnl_matrix inverseDirectionCosines{ vnl_matrix_inverse(directionCosineMatrix).as_matrix() }; + const vnl_matrix inverseDirectionCosines{ itk::Math::SVD(directionCosineMatrix).PseudoInverse() }; origin *= inverseDirectionCosines; // transform to minc convention diff --git a/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx b/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx index e21d616ebc8..81c163c460b 100644 --- a/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx +++ b/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx @@ -1555,7 +1555,7 @@ IsAffine(const mat44 & nifti_mat) } } - const double condition = vnl_matrix_inverse(mat.as_matrix()).well_condition(); + const double condition = itk::Math::SVD(mat.as_matrix()).WellCondition(); // Check matrix is invertible by testing condition number of inverse if (!(condition > std::numeric_limits::epsilon())) { @@ -1564,15 +1564,14 @@ IsAffine(const mat44 & nifti_mat) // Calculate the inverse and separate the inverse translation component // and the top 3x3 part of the inverse matrix - const vnl_matrix_fixed inv4x4Matrix = vnl_matrix_inverse(mat.as_matrix()).as_matrix(); + const vnl_matrix_fixed inv4x4Matrix = itk::Math::SVD(mat.as_matrix()).PseudoInverse(); const vnl_vector_fixed inv4x4Translation(inv4x4Matrix[0][3], inv4x4Matrix[1][3], inv4x4Matrix[2][3]); const vnl_matrix_fixed inv4x4Top3x3 = inv4x4Matrix.extract(3, 3, 0, 0); // Grab just the top 3x3 matrix const vnl_matrix_fixed top3x3Matrix = mat.extract(3, 3, 0, 0); - const vnl_matrix_fixed invTop3x3Matrix = - vnl_matrix_inverse(top3x3Matrix.as_matrix()).as_matrix(); - const vnl_vector_fixed inv3x3Translation = -(invTop3x3Matrix * mat.get_column(3).extract(3)); + const vnl_matrix_fixed invTop3x3Matrix = itk::Math::SVD(top3x3Matrix.as_matrix()).PseudoInverse(); + const vnl_vector_fixed inv3x3Translation = -(invTop3x3Matrix * mat.get_column(3).extract(3)); // Make sure we adhere to the conditions of a 4x4 invertible affine transform matrix const double diff_matrix_array_one_norm = (inv4x4Top3x3 - invTop3x3Matrix).array_one_norm(); diff --git a/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.h b/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.h index db3eb67b45e..4fec9f50ab2 100644 --- a/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.h +++ b/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.h @@ -22,7 +22,10 @@ #include "itkBooleanStdVector.h" #include "itkGradientDescentOptimizerv4.h" -#include "vnl/algo/vnl_matrix_inverse.h" +#include "itkMathSVD.h" +#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" namespace itk diff --git a/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.hxx b/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.hxx index b1750f41ca3..9489b0928ee 100644 --- a/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.hxx +++ b/Modules/Numerics/Optimizersv4/include/itkQuasiNewtonOptimizerv4.hxx @@ -399,9 +399,7 @@ QuasiNewtonOptimizerv4Template::ComputeHessianAnd gradient[p] = this->m_Gradient[offset + p]; } - const vnl_matrix hessianInverse{ - vnl_matrix_inverse(newHessian).as_matrix() - }; + const vnl_matrix hessianInverse{ itk::Math::SVD(newHessian).PseudoInverse() }; // gradient is already negated const DerivativeType newtonStep{ hessianInverse * gradient }; for (SizeValueType p = 0; p < numLocalPara; ++p) diff --git a/Modules/Numerics/Statistics/include/itkGaussianMembershipFunction.hxx b/Modules/Numerics/Statistics/include/itkGaussianMembershipFunction.hxx index 69f14bc7477..aec0aa57871 100644 --- a/Modules/Numerics/Statistics/include/itkGaussianMembershipFunction.hxx +++ b/Modules/Numerics/Statistics/include/itkGaussianMembershipFunction.hxx @@ -102,16 +102,12 @@ GaussianMembershipFunction::SetCovariance(const CovarianceMa } // the inverse of the covariance matrix is first computed by SVD - const vnl_matrix_inverse inv_cov(cov.GetVnlMatrix()); - - // Compute the *signed* determinant of the covariance matrix. - // vnl_matrix_inverse::determinant_magnitude() (which used to be used - // here) returns the product of singular values and is always - // non-negative, so it cannot detect a non-positive-definite matrix. - // Use vnl_determinant on the original covariance instead; this is - // O(n^3) but n is the measurement-vector dimension (typically very - // small) so the cost is negligible compared to the SVD already - // computed for the inverse. + 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()); if (det <= 0.0) @@ -131,7 +127,7 @@ GaussianMembershipFunction::SetCovariance(const CovarianceMa if (m_CovarianceNonsingular) { // allocate the memory for m_InverseCovariance matrix - m_InverseCovariance.GetVnlMatrix() = inv_cov.inverse(); + m_InverseCovariance.GetVnlMatrix() = inv_cov.PseudoInverse(); // calculate coefficient C of multivariate gaussian m_PreFactor = 1.0 / (std::sqrt(det) * std::pow(std::sqrt(2.0 * itk::Math::pi), diff --git a/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMembershipFunction.hxx b/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMembershipFunction.hxx index 5df0726d321..00d9ad6995a 100644 --- a/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMembershipFunction.hxx +++ b/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMembershipFunction.hxx @@ -21,7 +21,10 @@ #include "vnl/vnl_vector.h" #include "vnl/vnl_matrix.h" -#include "vnl/algo/vnl_matrix_inverse.h" +#include "itkMathSVD.h" +#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 namespace itk::Statistics { @@ -92,10 +95,10 @@ MahalanobisDistanceMembershipFunction::SetCovariance(const CovarianceMa m_Covariance = cov; // the inverse of the covariance matrix is first computed by SVD - const vnl_matrix_inverse inv_cov(m_Covariance.GetVnlMatrix()); + const auto inv_cov = itk::Math::SVD(m_Covariance.GetVnlMatrix()); // the determinant is then costless this way - const double det = inv_cov.determinant_magnitude(); + const double det = inv_cov.DeterminantMagnitude(); if (det < 0.) { @@ -109,7 +112,7 @@ MahalanobisDistanceMembershipFunction::SetCovariance(const CovarianceMa if (m_CovarianceNonsingular) { // allocate the memory for m_InverseCovariance matrix - m_InverseCovariance.GetVnlMatrix() = inv_cov.inverse(); + m_InverseCovariance.GetVnlMatrix() = inv_cov.PseudoInverse(); } else { diff --git a/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.h b/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.h index 459530c8d3b..cc0471b1fb8 100644 --- a/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.h +++ b/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.h @@ -22,7 +22,10 @@ #include "vnl/vnl_vector_ref.h" #include "vnl/vnl_transpose.h" #include "vnl/vnl_matrix.h" -#include "vnl/algo/vnl_matrix_inverse.h" +#include "itkMathSVD.h" +#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" diff --git a/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.hxx b/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.hxx index e6f7ab76b39..b35e1391ff6 100644 --- a/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.hxx +++ b/Modules/Numerics/Statistics/include/itkMahalanobisDistanceMetric.hxx @@ -128,7 +128,7 @@ MahalanobisDistanceMetric::CalculateInverseCovariance() } else { - m_InverseCovariance = vnl_matrix_inverse(m_Covariance).as_matrix(); + m_InverseCovariance = itk::Math::SVD(m_Covariance).PseudoInverse(); } } // end inverse calculations } diff --git a/Modules/Registration/Common/test/RegistrationITKv3/ImageRegistration9.cxx b/Modules/Registration/Common/test/RegistrationITKv3/ImageRegistration9.cxx index 3f9b51281a5..2c207acc998 100644 --- a/Modules/Registration/Common/test/RegistrationITKv3/ImageRegistration9.cxx +++ b/Modules/Registration/Common/test/RegistrationITKv3/ImageRegistration9.cxx @@ -65,6 +65,7 @@ // that will monitor the evolution of the registration process. // #include "itkCommand.h" +#include "itkMathSVD.h" class CommandIterationUpdate : public itk::Command { public: @@ -104,9 +105,9 @@ class CommandIterationUpdate : public itk::Command p[0][1] = static_cast(optimizer->GetCurrentPosition()[1]); p[1][0] = static_cast(optimizer->GetCurrentPosition()[2]); p[1][1] = static_cast(optimizer->GetCurrentPosition()[3]); - vnl_svd svd(p); + const auto svd = itk::Math::SVD(p); vnl_matrix r(2, 2); - r = svd.U() * vnl_transpose(svd.V()); + r = svd.U * svd.V.transpose(); const double angle = std::asin(r[1][0]); std::cout << " AffineAngle: " << angle * 180.0 / itk::Math::pi << std::endl; } @@ -359,15 +360,15 @@ main(int argc, char * argv[]) p[0][1] = static_cast(finalParameters[1]); p[1][0] = static_cast(finalParameters[2]); p[1][1] = static_cast(finalParameters[3]); - vnl_svd svd(p); + const auto svd = itk::Math::SVD(p); vnl_matrix r(2, 2); - r = svd.U() * vnl_transpose(svd.V()); + r = svd.U * svd.V.transpose(); const double angle = std::asin(r[1][0]); const double angleInDegrees = angle * 180.0 / itk::Math::pi; - std::cout << " Scale 1 = " << svd.W(0) << std::endl; - std::cout << " Scale 2 = " << svd.W(1) << std::endl; + std::cout << " Scale 1 = " << svd.W[0] << std::endl; + std::cout << " Scale 2 = " << svd.W[1] << std::endl; std::cout << " Angle (degrees) = " << angleInDegrees << std::endl; diff --git a/Modules/Segmentation/Classifiers/include/itkImageGaussianModelEstimator.h b/Modules/Segmentation/Classifiers/include/itkImageGaussianModelEstimator.h index ae809b49d1d..d5dce722306 100644 --- a/Modules/Segmentation/Classifiers/include/itkImageGaussianModelEstimator.h +++ b/Modules/Segmentation/Classifiers/include/itkImageGaussianModelEstimator.h @@ -19,6 +19,9 @@ #define itkImageGaussianModelEstimator_h #include +#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 #include // For unique_ptr. @@ -26,7 +29,6 @@ #include "vnl/vnl_matrix.h" #include "vnl/vnl_matrix_fixed.h" #include "itkMath.h" -#include "vnl/algo/vnl_matrix_inverse.h" #include "itkImageRegionIterator.h" #include "itkMacro.h" diff --git a/Modules/Segmentation/Classifiers/include/itkImageKmeansModelEstimator.h b/Modules/Segmentation/Classifiers/include/itkImageKmeansModelEstimator.h index 4411e67c0c2..88fdd59a5a6 100644 --- a/Modules/Segmentation/Classifiers/include/itkImageKmeansModelEstimator.h +++ b/Modules/Segmentation/Classifiers/include/itkImageKmeansModelEstimator.h @@ -19,13 +19,15 @@ #define itkImageKmeansModelEstimator_h #include +#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 #include #include "vnl/vnl_vector.h" #include "vnl/vnl_matrix.h" #include "itkMath.h" -#include "vnl/algo/vnl_matrix_inverse.h" #include "itkImageRegionIterator.h" #include "itkMacro.h" From 232a2ad2d448cf5fa303ea1e416fb2742af2884d Mon Sep 17 00:00:00 2001 From: "Hans J. Johnson" Date: Wed, 8 Jul 2026 07:27:50 -0500 Subject: [PATCH 4/4] ENH: Deprecate the vnl_svd family and build-gate the LINPACK svdc engine Add the three-state deprecation guard to vnl_svd, vnl_svd_fixed, vnl_svd_economy, vnl_matrix_inverse and vnl_solve_qp, pointing callers at the Eigen-backed itk::Math::SVD: silent by default, a #warning under ITK_LEGACY_REMOVE and a hard #error under ITK_FUTURE_LEGACY_REMOVE. Build-gate the vnl_svd family and the LINPACK svdc netlib sources out of the ITK_FUTURE_LEGACY_REMOVE build so the dead engine no longer compiles once the APIs are gone, and move the vnl algo tests off the deprecated spellings. --- .../Core/Common/test/itkVnlSVDEngineGTest.cxx | 17 +++++++---- .../VNL/src/vxl/core/vnl/algo/CMakeLists.txt | 29 +++++++++++++++---- .../src/vxl/core/vnl/algo/tests/test_algo.cxx | 7 ++++- .../core/vnl/algo/tests/test_complex_algo.cxx | 7 ++++- .../vxl/core/vnl/algo/vnl_matrix_inverse.h | 17 +++++++++++ .../VNL/src/vxl/core/vnl/algo/vnl_solve_qp.h | 19 ++++++++++++ .../VNL/src/vxl/core/vnl/algo/vnl_svd.h | 21 ++++++++++++++ .../src/vxl/core/vnl/algo/vnl_svd_economy.h | 19 ++++++++++++ .../VNL/src/vxl/core/vnl/algo/vnl_svd_fixed.h | 20 +++++++++++++ .../VNL/src/vxl/v3p/netlib/CMakeLists.txt | 5 ++++ 10 files changed, 149 insertions(+), 12 deletions(-) diff --git a/Modules/Core/Common/test/itkVnlSVDEngineGTest.cxx b/Modules/Core/Common/test/itkVnlSVDEngineGTest.cxx index 2f536036983..5ec0060facb 100644 --- a/Modules/Core/Common/test/itkVnlSVDEngineGTest.cxx +++ b/Modules/Core/Common/test/itkVnlSVDEngineGTest.cxx @@ -21,12 +21,17 @@ // core/vnl/algo/tests/test_svd.cxx so the coverage runs in ITK CI and guards // any future change of the underlying SVD engine. -#define ITK_LEGACY_TEST -#include "vnl/algo/vnl_svd.h" -#include "vnl/vnl_random.h" +#include "itkConfigure.h" -#include -#include +// The deprecated engine under test is unavailable under ITK_FUTURE_LEGACY_REMOVE. +#ifndef ITK_FUTURE_LEGACY_REMOVE + +# define ITK_LEGACY_TEST +# include "vnl/algo/vnl_svd.h" +# include "vnl/vnl_random.h" + +# include +# include namespace { @@ -186,3 +191,5 @@ TEST(VnlSVDEngine, Nullvector) nullvector>(1e-5, rng); nullvector>(1e-12, rng); } + +#endif // ITK_FUTURE_LEGACY_REMOVE diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/CMakeLists.txt b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/CMakeLists.txt index 7243f3cf737..cf350d28aa1 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/CMakeLists.txt +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/CMakeLists.txt @@ -10,10 +10,6 @@ set( vnl_algo_sources vnl_netlib.h # matrix decompositions - vnl_svd.hxx vnl_svd.h - vnl_svd_economy.hxx vnl_svd_economy.h - vnl_svd_fixed.hxx vnl_svd_fixed.h - vnl_matrix_inverse.hxx vnl_matrix_inverse.h vnl_qr.hxx vnl_qr.h vnl_cholesky.cxx vnl_cholesky.h vnl_ldl_cholesky.cxx vnl_ldl_cholesky.h @@ -29,7 +25,6 @@ set( vnl_algo_sources vnl_powell.cxx vnl_powell.h vnl_brent.cxx vnl_brent.h vnl_lsqr.cxx vnl_lsqr.h - vnl_solve_qp.cxx vnl_solve_qp.h vnl_bracket_minimum.cxx vnl_bracket_minimum.h vnl_brent_minimizer.cxx vnl_brent_minimizer.h @@ -59,11 +54,35 @@ if(NOT ITK_FUTURE_LEGACY_REMOVE) ) endif() +# vnl_matrix_inverse is a thin vnl_svd wrapper deprecated for itk::Math::SVD; +# drop it and its template instantiations under ITK_FUTURE_LEGACY_REMOVE. +if(NOT ITK_FUTURE_LEGACY_REMOVE) + list(APPEND vnl_algo_sources + vnl_matrix_inverse.hxx vnl_matrix_inverse.h + ) +endif() + +# LINPACK-svdc-backed legacy: vnl_svd{,_economy,_fixed} and vnl_solve_qp (its +# only client outside the deprecated eigensystem/matrix_inverse wrappers) are +# the consumers of the netlib svdc engine. They are deprecated for itk::Math::SVD; +# drop them and their template instantiations once the APIs are removed, leaving +# no svdc dependency to build. +if(NOT ITK_FUTURE_LEGACY_REMOVE) + list(APPEND vnl_algo_sources + vnl_svd.hxx vnl_svd.h + vnl_svd_economy.hxx vnl_svd_economy.h + vnl_svd_fixed.hxx vnl_svd_fixed.h + vnl_solve_qp.cxx vnl_solve_qp.h + ) +endif() + aux_source_directory(Templates vnl_algo_sources) if(ITK_FUTURE_LEGACY_REMOVE) list(FILTER vnl_algo_sources EXCLUDE REGEX "vnl_symmetric_eigensystem\\+") list(FILTER vnl_algo_sources EXCLUDE REGEX "vnl_scatter_3x3\\+") + list(FILTER vnl_algo_sources EXCLUDE REGEX "vnl_matrix_inverse\\+") + list(FILTER vnl_algo_sources EXCLUDE REGEX "vnl_svd.*\\+") endif() # If VXL_INSTALL_INCLUDE_DIR is the default value diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_algo.cxx b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_algo.cxx index 5a79cb43de7..dbfc5521faa 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_algo.cxx +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_algo.cxx @@ -24,7 +24,10 @@ #include #include #include -#include +#ifndef ITK_FUTURE_LEGACY_REMOVE +# define ITK_LEGACY_TEST // exercising the deprecated vnl_matrix_inverse +# include +#endif #include #include #include @@ -42,10 +45,12 @@ test_matrix_inverse() vnl_matrix V0 = svd.V(); TEST_NEAR("vnl_svd_economy", V[0][1], V0[0][1], 1e-6); +#ifndef ITK_FUTURE_LEGACY_REMOVE const vnl_matrix inv{ vnl_matrix_inverse(m).as_matrix() }; vnl_matrix identity(4, 4); identity.set_identity(); TEST_NEAR("vnl_matrix_inverse", (m * inv - identity).array_inf_norm(), 0, 1e-6); +#endif } class F_test_powell : public vnl_cost_function diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_complex_algo.cxx b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_complex_algo.cxx index 6def6279695..c1af00b8702 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_complex_algo.cxx +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_complex_algo.cxx @@ -14,7 +14,10 @@ #include "testlib/testlib_test.h" #include #include -#include +#ifndef ITK_FUTURE_LEGACY_REMOVE +# define ITK_LEGACY_TEST // exercising the deprecated vnl_matrix_inverse +# include +#endif static void test_matrix_inverse() @@ -42,10 +45,12 @@ test_matrix_inverse() vnl_matrix> V0 = svd.V(); TEST_NEAR("complex vnl_svd_economy", V[0][1], V0[0][1], 1e-6); +#ifndef ITK_FUTURE_LEGACY_REMOVE const vnl_matrix> inv{ vnl_matrix_inverse>(m).as_matrix() }; vnl_matrix> identity(4, 4); identity.set_identity(); TEST_NEAR("complex vnl_matrix_inverse", (m * inv - identity).array_inf_norm(), 0, 1e-6); +#endif } void diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_matrix_inverse.h b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_matrix_inverse.h index 9a14a81e0b5..072e77f3a0d 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_matrix_inverse.h +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_matrix_inverse.h @@ -1,6 +1,23 @@ // This is core/vnl/algo/vnl_matrix_inverse.h #ifndef vnl_matrix_inverse_h_ #define vnl_matrix_inverse_h_ + +// ITK deprecation shim: vnl_matrix_inverse is a thin vnl_svd wrapper being +// retired; itk::Math::SVD (PseudoInverse/Solve, Eigen-backed) is the supported +// replacement. The guard is active only when itkConfigure.h is reachable (an +// ITK consumer), so ITK's own VXL build is unaffected. +#if __has_include() +# include +# if defined(ITK_FUTURE_LEGACY_REMOVE) +# error "vnl/algo/vnl_matrix_inverse.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed)." +# elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST) +# if defined(_MSC_VER) +# pragma message("vnl/algo/vnl_matrix_inverse.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed).") +# else +# warning "vnl/algo/vnl_matrix_inverse.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed)." +# endif +# endif +#endif //: // \file // \brief Calculates inverse of a matrix (wrapper around vnl_svd) diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_solve_qp.h b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_solve_qp.h index f25b44ffdfb..f0e98cef43f 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_solve_qp.h +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_solve_qp.h @@ -2,6 +2,25 @@ #ifndef vnl_solve_qp_h_ #define vnl_solve_qp_h_ +// ITK deprecation shim: vnl_solve_qp is the last non-eigensystem consumer of the +// LINPACK-svdc engine (through vnl_svd); it has no ITK client and no itk::Math +// replacement. The guard is active only when itkConfigure.h is reachable (an ITK +// consumer), so ITK's own VXL build is unaffected. +#if __has_include() +# include +# if defined(ITK_FUTURE_LEGACY_REMOVE) +# error \ + "vnl/algo/vnl_solve_qp.h is deprecated and removed with the vnl_svd/LINPACK-svdc engine; no itk::Math replacement." +# elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST) +# if defined(_MSC_VER) +# pragma message( \ + "vnl/algo/vnl_solve_qp.h is deprecated and removed with the vnl_svd/LINPACK-svdc engine; no itk::Math replacement.") +# else +# warning \ + "vnl/algo/vnl_solve_qp.h is deprecated and removed with the vnl_svd/LINPACK-svdc engine; no itk::Math replacement." +# endif +# endif +#endif //: // \file // \brief Functions to solve various forms of constrained quadratic programming diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd.h b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd.h index ed2fff68f39..c3a4aa57bbd 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd.h +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd.h @@ -1,6 +1,27 @@ // This is core/vnl/algo/vnl_svd.h #ifndef vnl_svd_h_ #define vnl_svd_h_ + +// ITK deprecation shim: the Eigen-backed itk::Math::SVD (itkMathSVD.h) is the +// supported replacement for vnl_svd, offering PseudoInverse()/Solve()/Rank()/ +// Recompose() with fixed- and runtime-sized overloads. The guard is active only +// when itkConfigure.h is reachable (an ITK consumer), so ITK's own VXL build is +// unaffected. +#if __has_include() +# include +# if defined(ITK_FUTURE_LEGACY_REMOVE) +# error \ + "vnl/algo/vnl_svd.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed)." +# elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST) +# if defined(_MSC_VER) +# pragma message( \ + "vnl/algo/vnl_svd.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed).") +# else +# warning \ + "vnl/algo/vnl_svd.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed)." +# endif +# endif +#endif //: // \file // \brief Holds the singular value decomposition of a vnl_matrix. diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd_economy.h b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd_economy.h index 3de6278987f..44fc99820c5 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd_economy.h +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd_economy.h @@ -1,6 +1,25 @@ // This is core/vnl/algo/vnl_svd_economy.h #ifndef vnl_svd_economy_h_ #define vnl_svd_economy_h_ + +// ITK deprecation shim: the Eigen-backed itk::Math::SVD (itkMathSVD.h) is the +// supported replacement. The guard is active only when itkConfigure.h is +// reachable (an ITK consumer); ITK's own VXL build is unaffected. +#if __has_include() +# include +# if defined(ITK_FUTURE_LEGACY_REMOVE) +# error \ + "vnl/algo/vnl_svd_economy.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed)." +# elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST) +# if defined(_MSC_VER) +# pragma message( \ + "vnl/algo/vnl_svd_economy.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed).") +# else +# warning \ + "vnl/algo/vnl_svd_economy.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed)." +# endif +# endif +#endif //: // \file // \brief SVD wrapper that doesn't compute the left singular vectors, U. diff --git a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd_fixed.h b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd_fixed.h index 5df8b92bfa6..1b1a7a4460e 100644 --- a/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd_fixed.h +++ b/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd_fixed.h @@ -1,6 +1,26 @@ // This is core/vnl/algo/vnl_svd_fixed.h #ifndef vnl_svd_fixed_h_ #define vnl_svd_fixed_h_ + +// ITK deprecation shim: the Eigen-backed itk::Math::SVD (itkMathSVD.h) is the +// supported replacement (fixed-size square overload). The guard is active only +// when itkConfigure.h is reachable (an ITK consumer); ITK's own VXL build is +// unaffected. +#if __has_include() +# include +# if defined(ITK_FUTURE_LEGACY_REMOVE) +# error \ + "vnl/algo/vnl_svd_fixed.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed)." +# elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST) +# if defined(_MSC_VER) +# pragma message( \ + "vnl/algo/vnl_svd_fixed.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed).") +# else +# warning \ + "vnl/algo/vnl_svd_fixed.h is deprecated; migrate to itk::Math::SVD (itkMathSVD.h, Eigen-backed)." +# endif +# endif +#endif //: // \file // \brief Holds the singular value decomposition of a vnl_matrix_fixed. diff --git a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/CMakeLists.txt b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/CMakeLists.txt index e9db4a81ad0..838060afc7e 100644 --- a/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/CMakeLists.txt +++ b/Modules/ThirdParty/VNL/src/vxl/v3p/netlib/CMakeLists.txt @@ -116,6 +116,11 @@ set(V3P_NETLIB_eispack_SOURCES if(ITK_FUTURE_LEGACY_REMOVE) set(V3P_NETLIB_eispack_SOURCES "") endif() +# linpack svdc ({s,d,c,z}svdc) backs only the deprecated vnl_svd family and +# vnl_solve_qp. Drop the dead engine once those APIs are removed. +if(ITK_FUTURE_LEGACY_REMOVE) + set(V3P_NETLIB_linpack_SOURCES "") +endif() set(V3P_NETLIB_napack_SOURCES napack/cg.c napack/cg.h )