1 /*========================================================================= 2 * 3 * Copyright Insight Software Consortium 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 * http://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 itkSymmetricEigenAnalysis_h 19 #define itkSymmetricEigenAnalysis_h 20 21 #include "itkMacro.h" 22 #include "itk_eigen.h" 23 #include ITK_EIGEN(Eigenvalues) 24 #include <numeric> 25 #include <vector> 26 // For GetPointerToMatrixData 27 #include "vnl/vnl_matrix.h" 28 #include "vnl/vnl_matrix_fixed.h" 29 #include "itkMatrix.h" 30 31 namespace itk 32 { 33 namespace detail 34 { 35 /* Helper functions returning pointer to matrix data for different types. */ 36 template<typename TValueType, unsigned int NRows, unsigned int NCols> GetPointerToMatrixData(const vnl_matrix_fixed<TValueType,NRows,NCols> & inputMatrix)37 const TValueType* GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, NRows, NCols> & inputMatrix) 38 { 39 return inputMatrix.data_block(); 40 }; 41 template<typename TValueType> GetPointerToMatrixData(const vnl_matrix<TValueType> & inputMatrix)42 const TValueType* GetPointerToMatrixData(const vnl_matrix<TValueType> & inputMatrix) 43 { 44 return inputMatrix.data_block(); 45 }; 46 47 template<typename TValueType, unsigned int NRows, unsigned int NCols> GetPointerToMatrixData(const itk::Matrix<TValueType,NRows,NCols> & inputMatrix)48 const TValueType* GetPointerToMatrixData(const itk::Matrix<TValueType, NRows, NCols> & inputMatrix) 49 { 50 return inputMatrix.GetVnlMatrix().data_block(); 51 }; 52 53 /** Sort input to be ordered by magnitude, and returns container with the 54 * permutations required for the sorting. 55 * 56 * For example, if input eigenValues = {10, 0, 40}, the output would be: {2,0,1} 57 * and the eigenValues would be modified in-place: {40, 10, 0}. 58 * 59 * The permutations indices is used to order the matrix of eigenVectors. 60 * \sa permuteEigenVectorsWithSortPermutations 61 * 62 * @tparam TArray array type with operator [] 63 * @param eigenValues input array, requires operator [] 64 * @param numberOfElements size of array 65 * 66 * @return the permutations needed to sort the input array 67 */ 68 template<typename TArray> sortEigenValuesByMagnitude(TArray & eigenValues,const unsigned int numberOfElements)69 std::vector<int> sortEigenValuesByMagnitude(TArray & eigenValues, const unsigned int numberOfElements) 70 { 71 std::vector<int> indicesSortPermutations(numberOfElements, 0); 72 std::iota(std::begin(indicesSortPermutations), std::end(indicesSortPermutations), 0); 73 74 std::sort( 75 std::begin(indicesSortPermutations), 76 std::end(indicesSortPermutations), 77 [&eigenValues](unsigned int a, unsigned int b) 78 { return std::abs(eigenValues[a]) < std::abs(eigenValues[b]); } 79 ); 80 auto tmpCopy = eigenValues; 81 for (unsigned int i = 0; i < numberOfElements; ++i) 82 { 83 eigenValues[i] = tmpCopy[indicesSortPermutations[i]]; 84 } 85 return indicesSortPermutations; 86 } 87 88 /** Permute a eigenVectors matrix according to the permutation indices 89 * computed from the output of a sort function like \sa detail::sortEigenValuesByMagnitude 90 * 91 * @tparam QMatrix a Eigen3 matrix 92 * @param eigenVectors stored in columns 93 * @param indicesSortPermutations container with the permutations from the output of 94 * a sort function. 95 */ 96 template<typename QMatrix> permuteColumnsWithSortIndices(QMatrix & eigenVectors,const std::vector<int> & indicesSortPermutations)97 void permuteColumnsWithSortIndices( QMatrix & eigenVectors, 98 const std::vector<int> & indicesSortPermutations) 99 { 100 using EigenLibPermutationMatrix = 101 Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>; 102 auto numberOfElements = indicesSortPermutations.size(); 103 // Creates a NxN permutation matrix copying our permutation to the matrix indices. 104 // Which holds the 1D array representation of a permutation. 105 EigenLibPermutationMatrix perm(numberOfElements); 106 perm.setIdentity(); 107 std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(), 108 perm.indices().data()); 109 // Apply it 110 eigenVectors = eigenVectors * perm; 111 } 112 } // end namespace detail 113 114 /** \class SymmetricEigenAnalysis 115 * \brief Find Eigen values of a real 2D symmetric matrix. It 116 * serves as a thread-safe alternative to the class: 117 * vnl_symmetric_eigensystem, which uses netlib routines. 118 * 119 * The class is templated over the input matrix (which is expected to provide 120 * access to its elements with the [][] operator), matrix to store eigen 121 * values (must provide write operations on its elements with the [] operator), and 122 * EigenMatrix to store eigen vectors (must provide write access to its elements 123 * with the [][] operator). 124 * 125 * The SetOrderEigenValues() method can be used to order eigen values (and their 126 * corresponding eigen vectors if computed) in ascending order. This is the 127 * default ordering scheme. Eigen vectors and values can be obtained without 128 * ordering by calling SetOrderEigenValues(false). 129 * 130 * The SetOrderEigenMagnitudes() method can be used to order eigen values (and 131 * their corresponding eigen vectors if computed) by magnitude in ascending order. 132 * 133 * The user of this class is explicitly supposed to set the dimension of the 134 * 2D matrix using the SetDimension() method. 135 * 136 * The class contains routines taken from netlib sources (www.netlib.org). 137 * netlib/tql1.c 138 * netlib/tql2.c 139 * netlib/tred1.c 140 * netlib/tred2.c 141 * 142 * Reference: 143 * num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and 144 * wilkinson. 145 * handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). 146 * \ingroup ITKCommon 147 */ 148 149 template< typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix > 150 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis 151 { 152 public: 153 typedef enum 154 { 155 OrderByValue = 1, 156 OrderByMagnitude, 157 DoNotOrder 158 } 159 EigenValueOrderType; 160 SymmetricEigenAnalysis()161 SymmetricEigenAnalysis() : 162 m_OrderEigenValues(OrderByValue) {} 163 SymmetricEigenAnalysis(const unsigned int dimension)164 SymmetricEigenAnalysis(const unsigned int dimension): 165 m_UseEigenLibrary(false), 166 m_Dimension(dimension), 167 m_Order(dimension), 168 m_OrderEigenValues(OrderByValue) {} 169 170 ~SymmetricEigenAnalysis() = default; 171 172 using MatrixType = TMatrix; 173 using EigenMatrixType = TEigenMatrix; 174 using VectorType = TVector; 175 176 /** Compute Eigen values of A 177 * A is any type that overloads the [][] operator and contains the 178 * symmetric matrix. In practice only the upper triangle of the 179 * matrix will be accessed. (Both itk::Matrix and vnl_matrix 180 * overload [][] operator.) 181 * 182 * 'EigenValues' is any type that overloads the [][] operator and will contain 183 * the eigen values. 184 * 185 * No size checking is performed. A is expected to be a square matrix of size 186 * m_Dimension. 'EigenValues' is expected to be of length m_Dimension. 187 * The matrix is not checked to see if it is symmetric. 188 */ 189 unsigned int ComputeEigenValues( 190 const TMatrix & A, 191 TVector & EigenValues) const; 192 193 /** Compute Eigen values and vectors of A 194 * A is any type that overloads the [][] operator and contains the 195 * symmetric matrix. In practice only the upper triangle of the 196 * matrix will be accessed. (Both itk::Matrix and vnl_matrix 197 * overload [][] operator.) 198 * 199 * 'EigenValues' is any type that overloads the [] operator and will contain 200 * the eigen values. 201 * 202 * 'EigenVectors' is any type that provides access to its elements with the 203 * [][] operator. It is expected be of size m_Dimension * m_Dimension. 204 * 205 * No size checking is performed. A is expected to be a square matrix of size 206 * m_Dimension. 'EigenValues' is expected to be of length m_Dimension. 207 * The matrix is not checked to see if it is symmetric. 208 * 209 * Each row of the matrix 'EigenVectors' represents an eigen vector. (unlike MATLAB 210 * where the columns of the [EigenVectors, EigenValues] = eig(A) contains the 211 * eigenvectors). 212 */ 213 unsigned int ComputeEigenValuesAndVectors( 214 const TMatrix & A, 215 TVector & EigenValues, 216 TEigenMatrix & EigenVectors) const; 217 218 219 /** Matrix order. Defaults to matrix dimension if not set */ SetOrder(const unsigned int n)220 void SetOrder(const unsigned int n) 221 { 222 m_Order = n; 223 } 224 225 /** Get the Matrix order. Will be 0 unless explicitly set, or unless a 226 * call to SetDimension has been made in which case it will be the 227 * matrix dimension. */ GetOrder()228 unsigned int GetOrder() const { return m_Order; } 229 230 /** Set/Get methods to order the eigen values in ascending order. 231 * This is the default. ie lambda_1 < lambda_2 < .... 232 */ SetOrderEigenValues(const bool b)233 void SetOrderEigenValues(const bool b) 234 { 235 if ( b ) { m_OrderEigenValues = OrderByValue; } 236 else { m_OrderEigenValues = DoNotOrder; } 237 } 238 GetOrderEigenValues()239 bool GetOrderEigenValues() const { return ( m_OrderEigenValues == OrderByValue ); } 240 241 /** Set/Get methods to order the eigen value magnitudes in ascending order. 242 * In other words, |lambda_1| < |lambda_2| < ..... 243 */ SetOrderEigenMagnitudes(const bool b)244 void SetOrderEigenMagnitudes(const bool b) 245 { 246 if ( b ) { m_OrderEigenValues = OrderByMagnitude; } 247 else { m_OrderEigenValues = DoNotOrder; } 248 } 249 GetOrderEigenMagnitudes()250 bool GetOrderEigenMagnitudes() const { return ( m_OrderEigenValues == OrderByMagnitude ); } 251 252 /** Set the dimension of the input matrix A. A is a square matrix of 253 * size m_Dimension. */ SetDimension(const unsigned int n)254 void SetDimension(const unsigned int n) 255 { 256 m_Dimension = n; 257 if ( m_Order == 0 ) 258 { 259 m_Order = m_Dimension; 260 } 261 } 262 263 /** Get Matrix dimension, Will be 0 unless explicitly set by a 264 * call to SetDimension. */ GetDimension()265 unsigned int GetDimension() const { return m_Dimension; } 266 267 /** Set/Get to use Eigen library instead of vnl/netlib. */ SetUseEigenLibrary(const bool input)268 void SetUseEigenLibrary(const bool input) 269 { 270 m_UseEigenLibrary = input; 271 } SetUseEigenLibraryOn()272 void SetUseEigenLibraryOn() 273 { 274 m_UseEigenLibrary = true; 275 } SetUseEigenLibraryOff()276 void SetUseEigenLibraryOff() 277 { 278 m_UseEigenLibrary = false; 279 } GetUseEigenLibrary()280 bool GetUseEigenLibrary() const { return m_UseEigenLibrary; } 281 private: 282 bool m_UseEigenLibrary{false}; 283 unsigned int m_Dimension{0}; 284 unsigned int m_Order{0}; 285 EigenValueOrderType m_OrderEigenValues; 286 287 /** Reduces a real symmetric matrix to a symmetric tridiagonal matrix using 288 * orthogonal similarity transformations. 289 * 'inputMatrix' contains the real symmetric input matrix. Only the lower 290 * triangle of the matrix need be supplied. The upper triangle is unaltered. 291 * 'd' contains the diagonal elements of the tridiagonal matrix. 292 * 'e' contains the subdiagonal elements of the tridiagonal matrix in its 293 * last n-1 positions. e(1) is set to zero. 294 * 'e2' contains the squares of the corresponding elements of e. 295 * 'e2' may coincide with e if the squares are not needed. 296 * questions and comments should be directed to burton s. garbow. 297 * mathematics and computer science div, argonne national laboratory 298 * this version dated august 1983. 299 * 300 * Function adapted from netlib/tred1.c. 301 * [Changed: remove static vars, enforce const correctness. 302 * Use vnl routines as necessary]. 303 * Reference: 304 * num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson. 305 * handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). */ 306 void ReduceToTridiagonalMatrix(double *inputMatrix, double *d, 307 double *e, double *e2) const; 308 309 /** Reduces a real symmetric matrix to a symmetric tridiagonal matrix using 310 * and accumulating orthogonal similarity transformations. 311 * 'inputMatrix' contains the real symmetric input matrix. Only the lower 312 * triangle of the matrix need be supplied. The upper triangle is unaltered. 313 * 'diagonalElements' will contains the diagonal elements of the tridiagonal 314 * matrix. 315 * 'subDiagonalElements' will contain the subdiagonal elements of the tridiagonal 316 * matrix in its last n-1 positions. subDiagonalElements(1) is set to zero. 317 * 'transformMatrix' contains the orthogonal transformation matrix produced 318 * in the reduction. 319 * 320 * Questions and comments should be directed to Burton s. Garbow, 321 * Mathematics and Computer Science Div., Argonne National Laboratory. 322 * This version dated august 1983. 323 * 324 * Function adapted from netlib/tred2.c. 325 * [Changed: remove static vars, enforce const correctness. 326 * Use vnl routines as necessary]. 327 * Reference: 328 * num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson. 329 * handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). */ 330 void ReduceToTridiagonalMatrixAndGetTransformation( 331 double *inputMatrix, double *diagonalElements, 332 double *subDiagonalElements, double *transformMatrix) const; 333 334 /** Finds the eigenvalues of a symmetric tridiagonal matrix by the ql method. 335 * 336 * On input: 337 * 'd' contains the diagonal elements of the input matrix. 338 * 'e' contains the subdiagonal elements of the input matrix 339 * in its last n-1 positions. e(1) is arbitrary. 340 * On Output: 341 * 'd' contains the eigenvalues. 342 * 'e' has been destroyed. 343 * 344 * Returns: 345 * zero for normal return, 346 * j if the j-th eigenvalue has not been 347 * determined after 30 iterations. 348 * 349 * 350 * Reference 351 * This subroutine is a translation of the algol procedure tql1, 352 * num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and 353 * wilkinson. 354 * handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). 355 * 356 * Questions and comments should be directed to Burton s. Garbow, 357 * Mathematics and Computer Science Div., Argonne National Laboratory. 358 * This version dated august 1983. 359 * 360 * Function Adapted from netlib/tql1.c. 361 * [Changed: remove static vars, enforce const correctness. 362 * Use vnl routines as necessary] */ 363 unsigned int ComputeEigenValuesUsingQL(double *d, double *e) const; 364 365 /** Finds the eigenvalues and eigenvectors of a symmetric tridiagonal matrix 366 * by the ql method. 367 * 368 * On input: 369 * 'd' contains the diagonal elements of the input matrix. 370 * 'e' contains the subdiagonal elements of the input matrix 371 * in its last n-1 positions. e(1) is arbitrary. 372 * 'z' contains the transformation matrix produced in the reduction by 373 * ReduceToTridiagonalMatrixAndGetTransformation(), if performed. If the 374 * eigenvectors of the tridiagonal matrix are desired, z must contain 375 * the identity matrix. 376 377 * On Output: 378 * 'd' contains the eigenvalues. 379 * 'e' has been destroyed. 380 * 'z' contains orthonormal eigenvectors of the symmetric tridiagonal 381 * (or full) matrix. 382 * 383 * Returns: 384 * zero for normal return, 385 * j if the j-th eigenvalue has not been 386 * determined after 1000 iterations. 387 * 388 * Reference 389 * This subroutine is a translation of the algol procedure tql1, 390 * num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and 391 * wilkinson. 392 * handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). 393 * 394 * Questions and comments should be directed to Burton s. Garbow, 395 * Mathematics and Computer Science Div., Argonne National Laboratory. 396 * This version dated august 1983. 397 * 398 * Function Adapted from netlib/tql2.c. 399 * [Changed: remove static vars, enforce const correctness. 400 * Use vnl routines as necessary] 401 */ 402 unsigned int ComputeEigenValuesAndVectorsUsingQL(double *d, double *e, double *z) const; 403 404 /* Legacy algorithms using thread-safe netlib. 405 * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors 406 */ 407 unsigned int ComputeEigenValuesLegacy( 408 const TMatrix & A, 409 TVector & EigenValues) const; 410 411 unsigned int ComputeEigenValuesAndVectorsLegacy( 412 const TMatrix & A, 413 TVector & EigenValues, 414 TEigenMatrix & EigenVectors) const; 415 416 /* Helper to get the matrix value type for EigenLibMatrix typename. 417 * 418 * If the TMatrix is vnl, the type is in element_type. 419 * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType. 420 * 421 * To use this function: 422 * using ValueType = decltype(this->GetMatrixType(true)); 423 * 424 * \note The two `GetMatrixValueType` overloads have different 425 * parameter declarations (`bool` and `...`), to avoid that both 426 * functions are equally good candidates during overload resolution, 427 * in case `element_type` and `ValueType` are both nested types of 428 * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`). 429 */ 430 template<typename QMatrix = TMatrix > 431 auto GetMatrixValueType(bool) const 432 -> typename QMatrix::element_type 433 { 434 return QMatrix::element_type(); 435 } 436 template<typename QMatrix = TMatrix > 437 auto GetMatrixValueType(...) const 438 -> typename QMatrix::ValueType 439 { 440 return QMatrix::ValueType(); 441 } 442 443 /* Wrapper that call the right implementation for the type of matrix. */ ComputeEigenValuesAndVectorsWithEigenLibrary(const TMatrix & A,TVector & EigenValues,TEigenMatrix & EigenVectors)444 unsigned int ComputeEigenValuesAndVectorsWithEigenLibrary( 445 const TMatrix & A, 446 TVector & EigenValues, 447 TEigenMatrix & EigenVectors) const 448 { 449 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl( 450 A, EigenValues, EigenVectors, true); 451 } 452 453 /* Implementation detail using EigenLib that performs a copy of the input matrix. 454 * 455 * @param (long) implementation detail argument making this implementation less favourable 456 * to be chosen if alternatives are available. 457 * 458 * @return an unsigned int with no information value (no error code in EigenLib) */ 459 template<typename QMatrix > 460 auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl( 461 const QMatrix & A, 462 TVector & EigenValues, 463 TEigenMatrix & EigenVectors, 464 long) const 465 -> decltype(static_cast<unsigned int>(1)) 466 { 467 using ValueType = decltype(GetMatrixValueType(true)); 468 using EigenLibMatrixType = 469 Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; 470 EigenLibMatrixType inputMatrix( m_Dimension, m_Dimension); 471 for (unsigned int row = 0; row < m_Dimension; ++row) 472 { 473 for (unsigned int col = 0; col < m_Dimension; ++col) 474 { 475 inputMatrix(row, col) = A(row, col); 476 } 477 } 478 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>; 479 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors 480 const auto &eigenValues = solver.eigenvalues(); 481 /* Column k of the returned matrix is an eigenvector corresponding to 482 * eigenvalue number $ k $ as returned by eigenvalues(). 483 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */ 484 const auto &eigenVectors = solver.eigenvectors(); 485 486 if(m_OrderEigenValues == OrderByMagnitude) 487 { 488 auto copyEigenValues = eigenValues; 489 auto copyEigenVectors = eigenVectors; 490 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension); 491 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations); 492 493 for (unsigned int row = 0; row < m_Dimension; ++row) 494 { 495 EigenValues[row] = copyEigenValues[row]; 496 for (unsigned int col = 0; col < m_Dimension; ++col) 497 { 498 EigenVectors[row][col] = copyEigenVectors(col, row); 499 } 500 } 501 } 502 else 503 { 504 for (unsigned int row = 0; row < m_Dimension; ++row) 505 { 506 EigenValues[row] = eigenValues[row]; 507 for (unsigned int col = 0; col < m_Dimension; ++col) 508 { 509 EigenVectors[row][col] = eigenVectors(col, row); 510 } 511 } 512 } 513 // No error code 514 return 1; 515 } 516 517 518 /* Implementation detail using EigenLib that do not peform a copy. 519 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData 520 * If new types want to use this method, an appropiate overload of GetPointerToMatrixData 521 * should be included. 522 * 523 * @param (bool) implementation detail argument making this implementation the most favourable 524 * to be chosen from all the alternative implementations. 525 * 526 * @return an unsigned int with no information value (no error code in EigenLib) */ 527 template<typename QMatrix > 528 auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl( 529 const QMatrix & A, 530 TVector & EigenValues, 531 TEigenMatrix & EigenVectors, 532 bool) const 533 -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1)) 534 { 535 auto pointerToData = GetPointerToMatrixData(A); 536 using PointerType = decltype(pointerToData); 537 using ValueTypeCV = typename std::remove_pointer<PointerType>::type; 538 using ValueType = typename std::remove_cv<ValueTypeCV>::type; 539 using EigenLibMatrixType = 540 Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; 541 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>; 542 EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension); 543 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>; 544 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors 545 const auto & eigenValues = solver.eigenvalues(); 546 /* Column k of the returned matrix is an eigenvector corresponding to 547 * eigenvalue number $ k $ as returned by eigenvalues(). 548 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */ 549 const auto & eigenVectors = solver.eigenvectors(); 550 if(m_OrderEigenValues == OrderByMagnitude) 551 { 552 auto copyEigenValues = eigenValues; 553 auto copyEigenVectors = eigenVectors; 554 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension); 555 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations); 556 for (unsigned int row = 0; row < m_Dimension; ++row) 557 { 558 EigenValues[row] = copyEigenValues[row]; 559 for (unsigned int col = 0; col < m_Dimension; ++col) 560 { 561 EigenVectors[row][col] = copyEigenVectors(col, row); 562 } 563 } 564 } 565 else 566 { 567 for (unsigned int row = 0; row < m_Dimension; ++row) 568 { 569 EigenValues[row] = eigenValues[row]; 570 for (unsigned int col = 0; col < m_Dimension; ++col) 571 { 572 EigenVectors[row][col] = eigenVectors(col, row); 573 } 574 } 575 } 576 // No error code 577 return 1; 578 } 579 580 /* Wrapper that call the right implementation for the type of matrix. */ ComputeEigenValuesWithEigenLibrary(const TMatrix & A,TVector & EigenValues)581 unsigned int ComputeEigenValuesWithEigenLibrary( 582 const TMatrix & A, 583 TVector & EigenValues) const 584 { 585 return ComputeEigenValuesWithEigenLibraryImpl( 586 A, EigenValues, true); 587 } 588 589 /* Implementation detail using EigenLib that performs a copy of the input matrix. 590 * 591 * @param (long) implementation detail argument making this implementation less favourable 592 * to be chosen if alternatives are available. 593 * 594 * @return an unsigned int with no information value (no error code in EigenLib) */ 595 template<typename QMatrix > 596 auto ComputeEigenValuesWithEigenLibraryImpl( 597 const QMatrix & A, 598 TVector & EigenValues, 599 long) const 600 -> decltype(static_cast<unsigned int>(1)) 601 { 602 using ValueType = decltype(GetMatrixValueType(true)); 603 using EigenLibMatrixType = 604 Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; 605 EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension); 606 for (unsigned int row = 0; row < m_Dimension; ++row) 607 { 608 for (unsigned int col = 0; col < m_Dimension; ++col) 609 { 610 inputMatrix(row, col) = A(row, col); 611 } 612 } 613 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>; 614 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly); 615 auto eigenValues = solver.eigenvalues(); 616 if(m_OrderEigenValues == OrderByMagnitude) 617 { 618 detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension); 619 } 620 for (unsigned int i = 0; i < m_Dimension; ++i) 621 { 622 EigenValues[i] = eigenValues[i]; 623 } 624 625 // No error code 626 return 1; 627 } 628 629 /* Implementation detail using EigenLib that do not peform a copy. 630 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData 631 * If new types want to use this method, an appropiate overload of GetPointerToMatrixData 632 * should be included. 633 * 634 * @param (bool) implementation detail argument making this implementation the most favourable 635 * to be chosen from all the alternative implementations. 636 * 637 * @return an unsigned int with no information value (no error code in EigenLib) */ 638 template<typename QMatrix > 639 auto ComputeEigenValuesWithEigenLibraryImpl( 640 const QMatrix & A, 641 TVector & EigenValues, 642 bool) const 643 -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1)) 644 { 645 auto pointerToData = GetPointerToMatrixData(A); 646 using PointerType = decltype(pointerToData); 647 using ValueTypeCV = typename std::remove_pointer<PointerType>::type; 648 using ValueType = typename std::remove_cv<ValueTypeCV>::type; 649 using EigenLibMatrixType = 650 Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; 651 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>; 652 EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension); 653 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>; 654 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly); 655 auto eigenValues = solver.eigenvalues(); 656 if(m_OrderEigenValues == OrderByMagnitude) 657 { 658 detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension); 659 } 660 for (unsigned int i = 0; i < m_Dimension; ++i) 661 { 662 EigenValues[i] = eigenValues[i]; 663 } 664 // No error code 665 return 1; 666 } 667 }; 668 669 template< typename TMatrix, typename TVector, typename TEigenMatrix > 670 std::ostream & operator<<(std::ostream & os, 671 const SymmetricEigenAnalysis< TMatrix, TVector, TEigenMatrix > & s) 672 { 673 os << "[ClassType: SymmetricEigenAnalysis]" << std::endl; 674 os << " Dimension : " << s.GetDimension() << std::endl; 675 os << " Order : " << s.GetOrder() << std::endl; 676 os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl; 677 os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl; 678 os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl; 679 return os; 680 } 681 682 template< unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix > 683 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysisFixedDimension 684 { 685 public: 686 typedef enum 687 { 688 OrderByValue = 1, 689 OrderByMagnitude, 690 DoNotOrder 691 } 692 EigenValueOrderType; 693 SymmetricEigenAnalysisFixedDimension()694 SymmetricEigenAnalysisFixedDimension(): m_OrderEigenValues(OrderByValue) {} 695 ~SymmetricEigenAnalysisFixedDimension() = default; 696 697 using MatrixType = TMatrix; 698 using EigenMatrixType = TEigenMatrix; 699 using VectorType = TVector; 700 701 /** Compute Eigen values of A 702 * A is any type that overloads the [][] operator and contains the 703 * symmetric matrix. In practice only the upper triangle of the 704 * matrix will be accessed. (Both itk::Matrix and vnl_matrix 705 * overload [][] operator.) 706 * 707 * 'EigenValues' is any type that overloads the [] operator and will contain 708 * the eigen values. 709 * 710 * No size checking is performed. A is expected to be a square matrix of size 711 * VDimension. 'EigenValues' is expected to be of length VDimension. 712 * The matrix is not checked to see if it is symmetric. 713 */ ComputeEigenValues(const TMatrix & A,TVector & EigenValues)714 unsigned int ComputeEigenValues( 715 const TMatrix & A, 716 TVector & EigenValues) const 717 { 718 return ComputeEigenValuesWithEigenLibraryImpl( 719 A, EigenValues, true); 720 } 721 722 /** Compute Eigen values and vectors of A 723 * A is any type that overloads the [][] operator and contains the 724 * symmetric matrix. In practice only the upper triangle of the 725 * matrix will be accessed. (Both itk::Matrix and vnl_matrix 726 * overload [][] operator.) 727 * 728 * 'EigenValues' is any type that overloads the [] operator and will contain 729 * the eigen values. 730 * 731 * 'EigenVectors' is any type that provides access to its elements with the 732 * [][] operator. It is expected be of size VDimension * VDimension. 733 * 734 * No size checking is performed. A is expected to be a square matrix of size 735 * VDimension. 'EigenValues' is expected to be of length VDimension. 736 * The matrix is not checked to see if it is symmetric. 737 * 738 * Each row of the matrix 'EigenVectors' represents an eigen vector. (unlike MATLAB 739 * where the columns of the [EigenVectors, EigenValues] = eig(A) contains the 740 * eigenvectors). 741 */ ComputeEigenValuesAndVectors(const TMatrix & A,TVector & EigenValues,TEigenMatrix & EigenVectors)742 unsigned int ComputeEigenValuesAndVectors( 743 const TMatrix & A, 744 TVector & EigenValues, 745 TEigenMatrix & EigenVectors) const 746 { 747 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl( 748 A, EigenValues, EigenVectors, true); 749 } 750 SetOrderEigenValues(const bool b)751 void SetOrderEigenValues(const bool b) 752 { 753 if ( b ) { m_OrderEigenValues = OrderByValue; } 754 else { m_OrderEigenValues = DoNotOrder; } 755 } GetOrderEigenValues()756 bool GetOrderEigenValues() const { return ( m_OrderEigenValues == OrderByValue ); } SetOrderEigenMagnitudes(const bool b)757 void SetOrderEigenMagnitudes(const bool b) 758 { 759 if ( b ) { m_OrderEigenValues = OrderByMagnitude; } 760 else { m_OrderEigenValues = DoNotOrder; } 761 } GetOrderEigenMagnitudes()762 bool GetOrderEigenMagnitudes() const { return ( m_OrderEigenValues == OrderByMagnitude ); } GetOrder()763 constexpr unsigned int GetOrder() const { return VDimension; } GetDimension()764 constexpr unsigned int GetDimension() const { return VDimension; } GetUseEigenLibrary()765 constexpr bool GetUseEigenLibrary() const { return true; } 766 767 private: 768 EigenValueOrderType m_OrderEigenValues; 769 770 /* Helper to get the matrix value type for EigenLibMatrix typename. 771 * 772 * If the TMatrix is vnl, the type is in element_type. 773 * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType. 774 * 775 * To use this function: 776 * using ValueType = decltype(this->GetMatrixType(true)); 777 */ 778 template<typename QMatrix = TMatrix > 779 auto GetMatrixValueType(bool) const 780 -> typename QMatrix::element_type 781 { 782 return QMatrix::element_type(); 783 } 784 template<typename QMatrix = TMatrix > 785 auto GetMatrixValueType(bool) const 786 -> typename QMatrix::ValueType 787 { 788 return QMatrix::ValueType(); 789 } 790 791 /* Implementation detail using EigenLib that do not peform a copy. 792 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData 793 * If new types want to use this method, an appropiate overload of GetPointerToMatrixData 794 * should be included. 795 * 796 * @param (bool) implementation detail argument making this implementation the most favourable 797 * to be chosen from all the alternative implementations. 798 * 799 * @return an unsigned int with no information value (no error code in EigenLib) */ 800 template<typename QMatrix > 801 auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl( 802 const QMatrix & A, 803 TVector & EigenValues, 804 TEigenMatrix & EigenVectors, 805 bool) const 806 -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1)) 807 { 808 auto pointerToData = GetPointerToMatrixData(A); 809 using PointerType = decltype(pointerToData); 810 using ValueTypeCV = typename std::remove_pointer<PointerType>::type; 811 using ValueType = typename std::remove_cv<ValueTypeCV>::type; 812 using EigenLibMatrixType = Eigen::Matrix< ValueType, VDimension, VDimension, Eigen::RowMajor>; 813 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>; 814 EigenConstMatrixMap inputMatrix(pointerToData); 815 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>; 816 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors 817 const auto & eigenValues = solver.eigenvalues(); 818 /* Column k of the returned matrix is an eigenvector corresponding to 819 * eigenvalue number $ k $ as returned by eigenvalues(). 820 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */ 821 const auto & eigenVectors = solver.eigenvectors(); 822 if(m_OrderEigenValues == OrderByMagnitude) 823 { 824 auto copyEigenValues = eigenValues; 825 auto copyEigenVectors = eigenVectors; 826 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension); 827 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations); 828 for (unsigned int row = 0; row < VDimension; ++row) 829 { 830 EigenValues[row] = copyEigenValues[row]; 831 for (unsigned int col = 0; col < VDimension; ++col) 832 { 833 EigenVectors[row][col] = copyEigenVectors(col, row); 834 } 835 } 836 } 837 else 838 { 839 for (unsigned int row = 0; row < VDimension; ++row) 840 { 841 EigenValues[row] = eigenValues[row]; 842 for (unsigned int col = 0; col < VDimension; ++col) 843 { 844 EigenVectors[row][col] = eigenVectors(col, row); 845 } 846 } 847 } 848 // No error code 849 return 1; 850 } 851 852 /* Implementation detail using EigenLib that performs a copy of the input matrix. 853 * 854 * @param (long) implementation detail argument making this implementation less favourable 855 * to be chosen if alternatives are available. 856 * 857 * @return an unsigned int with no information value (no error code in EigenLib) */ 858 template<typename QMatrix > 859 auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl( 860 const QMatrix & A, 861 TVector & EigenValues, 862 TEigenMatrix & EigenVectors, 863 long) const 864 -> decltype(static_cast<unsigned int>(1)) 865 { 866 using ValueType = decltype(GetMatrixValueType(true)); 867 using EigenLibMatrixType = Eigen::Matrix< ValueType, VDimension, VDimension, Eigen::RowMajor>; 868 EigenLibMatrixType inputMatrix; 869 for (unsigned int row = 0; row < VDimension; ++row) 870 { 871 for (unsigned int col = 0; col < VDimension; ++col) 872 { 873 inputMatrix(row, col) = A(row, col); 874 } 875 } 876 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>; 877 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors 878 const auto &eigenValues = solver.eigenvalues(); 879 /* Column k of the returned matrix is an eigenvector corresponding to 880 * eigenvalue number $ k $ as returned by eigenvalues(). 881 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */ 882 const auto &eigenVectors = solver.eigenvectors(); 883 884 if(m_OrderEigenValues == OrderByMagnitude) 885 { 886 auto copyEigenValues = eigenValues; 887 auto copyEigenVectors = eigenVectors; 888 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension); 889 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations); 890 891 for (unsigned int row = 0; row < VDimension; ++row) 892 { 893 EigenValues[row] = copyEigenValues[row]; 894 for (unsigned int col = 0; col < VDimension; ++col) 895 { 896 EigenVectors[row][col] = copyEigenVectors(col, row); 897 } 898 } 899 } 900 else 901 { 902 for (unsigned int row = 0; row < VDimension; ++row) 903 { 904 EigenValues[row] = eigenValues[row]; 905 for (unsigned int col = 0; col < VDimension; ++col) 906 { 907 EigenVectors[row][col] = eigenVectors(col, row); 908 } 909 } 910 } 911 // No error code 912 return 1; 913 } 914 915 /* Implementation detail using EigenLib that performs a copy of the input matrix. 916 * 917 * @param (long) implementation detail argument making this implementation less favourable 918 * to be chosen if alternatives are available. 919 * 920 * @return an unsigned int with no information value (no error code in EigenLib) */ 921 template<typename QMatrix > 922 auto ComputeEigenValuesWithEigenLibraryImpl( 923 const QMatrix & A, 924 TVector & EigenValues, 925 long) const 926 -> decltype(static_cast<unsigned int>(1)) 927 { 928 using ValueType = decltype(GetMatrixValueType(true)); 929 using EigenLibMatrixType = Eigen::Matrix< ValueType, VDimension, VDimension, Eigen::RowMajor>; 930 EigenLibMatrixType inputMatrix; 931 for (unsigned int row = 0; row < VDimension; ++row) 932 { 933 for (unsigned int col = 0; col < VDimension; ++col) 934 { 935 inputMatrix(row, col) = A(row, col); 936 } 937 } 938 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>; 939 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly); 940 auto eigenValues = solver.eigenvalues(); 941 if(m_OrderEigenValues == OrderByMagnitude) 942 { 943 detail::sortEigenValuesByMagnitude(eigenValues, VDimension); 944 } 945 for (unsigned int i = 0; i < VDimension; ++i) 946 { 947 EigenValues[i] = eigenValues[i]; 948 } 949 950 // No error code 951 return 1; 952 } 953 954 /* Implementation detail using EigenLib that do not peform a copy. 955 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData 956 * If new types want to use this method, an appropiate overload of GetPointerToMatrixData 957 * should be included. 958 * 959 * @param (bool) implementation detail argument making this implementation the most favourable 960 * to be chosen from all the alternative implementations. 961 * 962 * @return an unsigned int with no information value (no error code in EigenLib) */ 963 template<typename QMatrix > 964 auto ComputeEigenValuesWithEigenLibraryImpl( 965 const QMatrix & A, 966 TVector & EigenValues, 967 bool) const 968 -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1)) 969 { 970 auto pointerToData = GetPointerToMatrixData(A); 971 using PointerType = decltype(pointerToData); 972 using ValueTypeCV = typename std::remove_pointer<PointerType>::type; 973 using ValueType = typename std::remove_cv<ValueTypeCV>::type; 974 using EigenLibMatrixType = Eigen::Matrix< ValueType, VDimension, VDimension, Eigen::RowMajor>; 975 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>; 976 EigenConstMatrixMap inputMatrix(pointerToData); 977 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>; 978 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly); 979 auto eigenValues = solver.eigenvalues(); 980 if(m_OrderEigenValues == OrderByMagnitude) 981 { 982 detail::sortEigenValuesByMagnitude(eigenValues, VDimension); 983 } 984 for (unsigned int i = 0; i < VDimension; ++i) 985 { 986 EigenValues[i] = eigenValues[i]; 987 } 988 // No error code 989 return 1; 990 } 991 }; 992 993 template< unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix > 994 std::ostream & operator<<(std::ostream & os, 995 const SymmetricEigenAnalysisFixedDimension<VDimension, 996 TMatrix, TVector, TEigenMatrix > & s) 997 { 998 os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl; 999 os << " Dimension : " << s.GetDimension() << std::endl; 1000 os << " Order : " << s.GetOrder() << std::endl; 1001 os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl; 1002 os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl; 1003 os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl; 1004 return os; 1005 } 1006 } // end namespace itk 1007 1008 #ifndef ITK_MANUAL_INSTANTIATION 1009 #include "itkSymmetricEigenAnalysis.hxx" 1010 #endif 1011 1012 #endif 1013