Skip to content

Commit c420eca

Browse files
hjmjohnsonhjmjohnson
authored
Merge pull request #6499 from hjmjohnson/comp-future-legacy-vnl-qr
ENH: Extend itk::QRDecomposition and migrate in-tree vnl_qr consumers
2 parents badfa18 + ddd4db1 commit c420eca

14 files changed

Lines changed: 731 additions & 67 deletions

‎Modules/Core/Common/include/itkQRDecomposition.h‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -109,6 +109,48 @@ class QRDecomposition
109109
return out;
110110
}
111111

112+
/** Solve A X = B for a multi-column right-hand side. Each column of B is an
113+
* independent system sharing the single stored factorization, so the result
114+
* X has one solution column per column of B (X.cols() == B.cols()); column j
115+
* equals Solve(B.get_column(j)). B.rows() must match the row count and A must
116+
* have full column rank, as in the vector overload. */
117+
MatrixType
118+
Solve(const MatrixType & B) const
119+
{
120+
using RowMajor = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
121+
const unsigned int rows = m_Q.rows();
122+
const unsigned int cols = m_R.cols();
123+
const unsigned int nrhs = B.cols();
124+
itkAssertOrThrowMacro(B.rows() == rows, "QRDecomposition::Solve requires B row count to match the row count");
125+
if (cols > rows)
126+
{
127+
itkGenericExceptionMacro("QRDecomposition::Solve supports square or overdetermined systems only (rows >= cols)");
128+
}
129+
130+
Eigen::Map<const RowMajor> qMap(m_Q.data_block(), rows, rows);
131+
Eigen::Map<const RowMajor> rMap(m_R.data_block(), rows, cols);
132+
Eigen::Map<const RowMajor> bMap(B.data_block(), rows, nrhs);
133+
134+
const RowMajor qtb = qMap.adjoint() * bMap;
135+
const RowMajor x = rMap.topLeftCorner(cols, cols).template triangularView<Eigen::Upper>().solve(qtb.topRows(cols));
136+
137+
MatrixType out(cols, nrhs);
138+
Eigen::Map<RowMajor>(out.data_block(), cols, nrhs) = x;
139+
return out;
140+
}
141+
142+
/** Inverse of a square A, computed as the solution of A X = I. Requires A to
143+
* be square and full rank; a rank-deficient A yields a non-finite result. */
144+
MatrixType
145+
Inverse() const
146+
{
147+
const unsigned int rows = m_Q.rows();
148+
itkAssertOrThrowMacro(rows == m_R.cols(), "QRDecomposition::Inverse requires a square matrix");
149+
MatrixType identity(rows, rows);
150+
identity.set_identity();
151+
return this->Solve(identity);
152+
}
153+
112154
/** Determinant of a square A (zero-initialized value for non-square A). */
113155
T
114156
GetDeterminant() const

‎Modules/Core/Common/test/VNLSparseLUSolverTraitsGTest.cxx‎

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -16,13 +16,16 @@
1616
*
1717
*=========================================================================*/
1818

19-
#define ITK_LEGACY_TEST // exercises the deprecated VNLSparseLUSolverTraits on purpose
20-
#include "VNLSparseLUSolverTraits.h"
21-
#include "itkGTest.h"
22-
#include "itkMath.h" // itk::Math::Absolute
19+
#include "itkConfigure.h" // defines ITK_FUTURE_LEGACY_REMOVE before the guard below
20+
// VNLSparseLUSolverTraits is unavailable under ITK_FUTURE_LEGACY_REMOVE.
21+
#ifndef ITK_FUTURE_LEGACY_REMOVE
22+
# define ITK_LEGACY_TEST // exercises the deprecated VNLSparseLUSolverTraits on purpose
23+
# include "VNLSparseLUSolverTraits.h"
24+
# include "itkGTest.h"
25+
# include "itkMath.h" // itk::Math::Absolute
2326

24-
#include <iostream>
25-
#include <cstdlib>
27+
# include <iostream>
28+
# include <cstdlib>
2629

2730
template <class TVector>
2831
bool
@@ -234,3 +237,4 @@ DoVNLSparseLUSolverTraitsTest(int, char *[])
234237

235238

236239
TEST(VNLSparseLUSolverTraits, ConvertedLegacyTest) { EXPECT_EQ(0, DoVNLSparseLUSolverTraitsTest(0, nullptr)); }
240+
#endif

‎Modules/Core/Common/test/VnlEigenSparseLUEquivalenceGTest.cxx‎

Lines changed: 11 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -22,13 +22,17 @@
2222
// direct LU solvers for the same A x = b, so agreement to tight tolerance is
2323
// required before the Eigen algorithm may be placed behind the vnl_* API.
2424

25-
#define ITK_LEGACY_TEST // intentionally exercises the deprecated VNLSparseLUSolverTraits for equivalence
26-
#include "VNLSparseLUSolverTraits.h"
27-
#include "SparseLUSolverTraits.h"
28-
#include "itkGTest.h"
25+
#include "itkConfigure.h" // defines ITK_FUTURE_LEGACY_REMOVE before the guard below
26+
// The vnl_sparse_lu equivalence comparison is unavailable under
27+
// ITK_FUTURE_LEGACY_REMOVE, where VNLSparseLUSolverTraits is removed.
28+
#ifndef ITK_FUTURE_LEGACY_REMOVE
29+
# define ITK_LEGACY_TEST // intentionally exercises the deprecated VNLSparseLUSolverTraits for equivalence
30+
# include "VNLSparseLUSolverTraits.h"
31+
# include "SparseLUSolverTraits.h"
32+
# include "itkGTest.h"
2933

30-
#include <cmath>
31-
#include <vector>
34+
# include <cmath>
35+
# include <vector>
3236

3337
namespace
3438
{
@@ -131,3 +135,4 @@ TEST(VnlEigenSparseLUEquivalence, WellConditioned)
131135
// Ill-conditioned system: tolerance is looser because both solvers amplify the
132136
// conditioning, but they must still track each other closely.
133137
TEST(VnlEigenSparseLUEquivalence, IllConditioned) { expectEquivalent(24, 1e-9, 1e-6); }
138+
#endif

‎Modules/Core/Common/test/itkQRDecompositionGTest.cxx‎

Lines changed: 89 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,8 +20,11 @@
2020
#include "itkQRDecomposition.h"
2121

2222
// Exercise the deprecated VNL engine for the old-vs-new equivalence checks.
23-
#define ITK_LEGACY_TEST
24-
#include "vnl/algo/vnl_qr.h"
23+
// The VNL symbol is unavailable under ITK_FUTURE_LEGACY_REMOVE.
24+
#ifndef ITK_FUTURE_LEGACY_REMOVE
25+
# define ITK_LEGACY_TEST
26+
# include "vnl/algo/vnl_qr.h"
27+
#endif
2528

2629
#include <gtest/gtest.h>
2730
#include <cmath>
@@ -77,6 +80,51 @@ TEST(QRDecomposition, SolveResidual)
7780
}
7881

7982

83+
// Multi-column RHS: each column solves against the one shared factorization.
84+
TEST(QRDecomposition, SolveMatrixRHS)
85+
{
86+
const unsigned int n = 5;
87+
const unsigned int k = 3;
88+
const vnl_matrix<double> A = MakeMatrix<double>(n, n);
89+
vnl_matrix<double> B(n, k);
90+
for (unsigned int i = 0; i < n; ++i)
91+
for (unsigned int j = 0; j < k; ++j)
92+
B(i, j) = std::cos(0.3 * (i + 1) * (j + 2));
93+
94+
const itk::QRDecomposition<double> qr(A);
95+
const vnl_matrix<double> X = qr.Solve(B);
96+
ASSERT_EQ(X.rows(), n);
97+
ASSERT_EQ(X.cols(), k);
98+
EXPECT_LT((A * X - B).fro_norm() / B.fro_norm(), 1e-10);
99+
for (unsigned int j = 0; j < k; ++j)
100+
EXPECT_LT((X.get_column(j) - qr.Solve(B.get_column(j))).two_norm(), 1e-12);
101+
}
102+
103+
104+
// Overdetermined (rows > cols) multi-column solve yields the least-squares
105+
// solution, characterized by the normal equations A^T (A X - B) == 0.
106+
TEST(QRDecomposition, SolveMatrixRHSOverdetermined)
107+
{
108+
const unsigned int m = 6;
109+
const unsigned int n = 3;
110+
const unsigned int k = 2;
111+
const vnl_matrix<double> A = MakeMatrix<double>(m, n);
112+
vnl_matrix<double> B(m, k);
113+
for (unsigned int i = 0; i < m; ++i)
114+
for (unsigned int j = 0; j < k; ++j)
115+
B(i, j) = std::cos(0.35 * (i + 1) * (j + 2)) - 0.2 * (j + 1);
116+
117+
const itk::QRDecomposition<double> qr(A);
118+
const vnl_matrix<double> X = qr.Solve(B);
119+
ASSERT_EQ(X.rows(), n);
120+
ASSERT_EQ(X.cols(), k);
121+
EXPECT_LT((A.transpose() * (A * X - B)).fro_norm() / (A.transpose() * B).fro_norm(), 1e-10);
122+
for (unsigned int j = 0; j < k; ++j)
123+
EXPECT_LT((X.get_column(j) - qr.Solve(B.get_column(j))).two_norm(), 1e-12);
124+
}
125+
126+
127+
#ifndef ITK_FUTURE_LEGACY_REMOVE
80128
// itk:: solve and determinant agree with the (sign-stable) legacy vnl_qr.
81129
TEST(QRDecomposition, EquivalentToVnlQR)
82130
{
@@ -94,6 +142,45 @@ TEST(QRDecomposition, EquivalentToVnlQR)
94142
}
95143

96144

145+
// Multi-column solve matches legacy vnl_qr column-for-column. This is the exact
146+
// operation itkLandmarkBasedTransformInitializer relies on (square Q, matrix C).
147+
TEST(QRDecomposition, MatrixRHSEquivalentToVnlQR)
148+
{
149+
for (const unsigned int n : { 3u, 4u, 6u })
150+
{
151+
const unsigned int k = n - 1;
152+
const vnl_matrix<double> A = MakeMatrix<double>(n, n);
153+
vnl_matrix<double> B(n, k);
154+
for (unsigned int i = 0; i < n; ++i)
155+
for (unsigned int j = 0; j < k; ++j)
156+
B(i, j) = std::cos(0.4 * (i + 1) * (j + 2)) - 0.3 * (j + 1);
157+
158+
const vnl_matrix<double> xItk = itk::QRDecomposition<double>(A).Solve(B);
159+
const vnl_matrix<double> xVnl = vnl_qr<double>(A).solve(B);
160+
EXPECT_LT((xItk - xVnl).fro_norm() / xVnl.fro_norm(), 1e-10);
161+
}
162+
}
163+
164+
165+
// Inverse matches legacy vnl_qr.inverse(). This is the operation itk::fem::Element
166+
// (JacobianInverse / JacobianDeterminant) relies on.
167+
TEST(QRDecomposition, InverseEquivalentToVnlQR)
168+
{
169+
for (const unsigned int n : { 2u, 3u, 5u })
170+
{
171+
const vnl_matrix<double> A = MakeMatrix<double>(n, n);
172+
const vnl_matrix<double> invItk = itk::QRDecomposition<double>(A).Inverse();
173+
const vnl_matrix<double> invVnl = vnl_qr<double>(A).inverse();
174+
EXPECT_LT((invItk - invVnl).fro_norm() / invVnl.fro_norm(), 1e-10);
175+
176+
vnl_matrix<double> ident(n, n);
177+
ident.set_identity();
178+
EXPECT_LT((A * invItk - ident).fro_norm(), 1e-10);
179+
}
180+
}
181+
#endif
182+
183+
97184
// Single-precision path reconstructs.
98185
TEST(QRDecomposition, FloatReconstructs)
99186
{

‎Modules/Numerics/FEM/include/itpack.h‎

Lines changed: 0 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -525,16 +525,6 @@ extern "C"
525525
extern doublereal
526526
ddot_(integer * n, doublereal * dx, integer * incx, doublereal * dy, integer * incy);
527527

528-
/**
529-
* Subroutine that computes the determinant of a symmetric tridiagonal matrix
530-
* given by tri. det(tri - xlmda*i) = 0
531-
* \param n order of tridiagonal system
532-
* \param tri symmetric tridiagonal matrix of order n
533-
* \param xlmda argument for characteristic equation
534-
*/
535-
extern doublereal
536-
determ_(integer * n, doublereal * tri, doublereal * xlmda);
537-
538528
/**
539529
* Obtain default parameters
540530
* \param iparm array of 12 integer parameters

‎Modules/Numerics/FEM/src/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,9 @@
11
set(
22
ITKFEM_SRCS
33
dsrc2c.c
4+
itkfem_tqlrat.c
5+
itkfem_pythag.c
6+
itkfem_epslon.c
47
itkFEMElement2DC0LinearLine.cxx
58
itkFEMElement2DC0LinearLineStress.cxx
69
itkFEMElement2DC0LinearQuadrilateral.cxx

‎Modules/Numerics/FEM/src/dsrc2c.c‎

Lines changed: 2 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626

2727
#define V3P_NETLIB_SRC
2828
#include "v3p_netlib.h"
29+
#include "itkfem_eispack.h"
2930

3031
/* Modified by Peter Vanroose, Oct 2003: manual optimisation and clean-up */
3132

@@ -38,8 +39,6 @@ time(long * timer); /* #include <time.h> */
3839
extern doublereal
3940
cheby_(doublereal *, doublereal *, doublereal *, integer *, doublereal *, doublereal *);
4041
extern doublereal
41-
determ_(integer *, doublereal *, doublereal *);
42-
extern doublereal
4342
eigvns_(integer *, doublereal *, doublereal *, doublereal *, integer *);
4443
extern doublereal
4544
itpackddot_(integer *, doublereal *, integer *, doublereal *, integer *);
@@ -4433,39 +4432,6 @@ itpackddot_(integer * n, doublereal * dx, integer * incx, doublereal * dy, integ
44334432
return ret_val;
44344433
} /* itpackddot_ */
44354434

4436-
doublereal
4437-
determ_(integer * n, doublereal * tri, doublereal * xlmda)
4438-
{
4439-
/* Local variables */
4440-
static integer l;
4441-
static doublereal d1, d2, d3;
4442-
static integer icnt;
4443-
4444-
/* THIS SUBROUTINE COMPUTES THE DETERMINANT OF A SYMMETRIC */
4445-
/* TRIDIAGONAL MATRIX GIVEN BY TRI. DET(TRI - XLMDA*I) = 0 */
4446-
4447-
/* ... PARAMETER LIST */
4448-
4449-
/* N ORDER OF TRIDIAGONAL SYSTEM */
4450-
/* TRI SYMMETRIC TRIDIAGONAL MATRIX OF ORDER N */
4451-
/* XLMDA ARGUMENT FOR CHARACTERISTIC EQUATION */
4452-
4453-
d2 = tri[(*n << 1) - 2] - *xlmda;
4454-
d1 = d2 * (tri[(*n << 1) - 4] - *xlmda) - tri[(*n << 1) - 1];
4455-
if (*n == 2)
4456-
return d1;
4457-
4458-
for (icnt = 2; icnt < *n; ++icnt)
4459-
{
4460-
l = *n - icnt + 1;
4461-
d3 = d2;
4462-
d2 = d1;
4463-
d1 = (tri[((l - 1) << 1) - 2] - *xlmda) * d2 - d3 * tri[(l << 1) - 1];
4464-
}
4465-
4466-
return d1;
4467-
} /* determ_ */
4468-
44694435
/* Subroutine */
44704436
int
44714437
dfault_(integer * iparm, doublereal * rparm)
@@ -4782,7 +4748,7 @@ eigvns_(integer * n, doublereal * tri, doublereal * d, doublereal * e2, integer
47824748
}
47834749

47844750
// eqrt1s_(d, e2, n, &c__1, &c__0, ier);
4785-
v3p_netlib_tqlrat_(n, d, e2, ier);
4751+
itkfem_tqlrat_(n, d, e2, ier);
47864752

47874753
return -d[0];
47884754
} /* eigvns_ */

‎Modules/Numerics/FEM/src/itkFEMElement3DC0LinearTriangular.cxx‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@
2020

2121
#include "itkFEMElement3DC0LinearTriangular.h"
2222

23-
#include "vnl/algo/vnl_qr.h"
23+
#include "itkQRDecomposition.h"
2424

2525
namespace itk::fem
2626
{
@@ -267,7 +267,7 @@ Element3DC0LinearTriangular::JacobianInverse(const VectorType & pt, MatrixType &
267267
}
268268

269269
// invJ=vnl_svd_inverse<Float>(*pJ);
270-
invJ = vnl_qr<Float>(*pJ).inverse();
270+
invJ = itk::QRDecomposition<Float>(*pJ).Inverse();
271271

272272
/*
273273
// Note that inverse of Jacobian is not quadratic matrix

‎Modules/Numerics/FEM/src/itkFEMElementBase.cxx‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@
1717
*=========================================================================*/
1818

1919
#include "itkFEMElementBase.h"
20-
#include "vnl/algo/vnl_qr.h"
20+
#include "itkQRDecomposition.h"
2121

2222
namespace itk::fem
2323
{
@@ -305,7 +305,7 @@ Element::JacobianDeterminant(const VectorType & pt, const MatrixType * pJ) const
305305
}
306306

307307
// Float det=vnl_svd<Float>(*pJ).determinant_magnitude();
308-
Float det = vnl_qr<Float>(*pJ).determinant();
308+
Float det = itk::QRDecomposition<Float>(*pJ).GetDeterminant();
309309

310310
return det;
311311
}
@@ -324,7 +324,7 @@ Element::JacobianInverse(const VectorType & pt, MatrixType & invJ, const MatrixT
324324
}
325325

326326
// invJ=vnl_svd_inverse<Float>(*pJ);
327-
invJ = vnl_qr<Float>(*pJ).inverse();
327+
invJ = itk::QRDecomposition<Float>(*pJ).Inverse();
328328
}
329329

330330
void

‎Modules/Numerics/FEM/src/itkfem_eispack.h‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* https://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
#ifndef itkfem_eispack_h
19+
#define itkfem_eispack_h
20+
21+
/* FEM-private EISPACK eigenvalue helpers (renamed tqlrat/pythag/epslon).
22+
* v3p_netlib drops eispack under ITK_FUTURE_LEGACY_REMOVE; the itk::fem ITPACK
23+
* solver (dsrc2c.c eigvns_) needs tqlrat, so FEM owns self-contained copies and
24+
* no longer depends on the v3p_netlib eispack in any configuration. */
25+
26+
#ifdef __cplusplus
27+
extern "C"
28+
{
29+
#endif
30+
31+
int
32+
itkfem_tqlrat_(long * n, double * d, double * e2, long * ierr);
33+
double
34+
itkfem_pythag_(double * a, double * b);
35+
double
36+
itkfem_epslon_(double * x);
37+
38+
#ifdef __cplusplus
39+
}
40+
#endif
41+
42+
#endif

0 commit comments

Comments
 (0)