ITK  5.0.0
Insight Segmentation and Registration Toolkit
itkSymmetricEigenAnalysis.h
Go to the documentation of this file.
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>
37  const TValueType* GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, NRows, NCols> & inputMatrix)
38  {
39  return inputMatrix.data_block();
40  };
41  template<typename TValueType>
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>
48  const TValueType* GetPointerToMatrixData(const itk::Matrix<TValueType, NRows, NCols> & inputMatrix)
49  {
50  return inputMatrix.GetVnlMatrix().data_block();
51  };
52 
68  template<typename TArray>
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);
74 
75  std::sort(
76  std::begin(indicesSortPermutations),
77  std::end(indicesSortPermutations),
78  [&eigenValues](unsigned int a, unsigned int b)
79  { return std::abs(eigenValues[a]) < std::abs(eigenValues[b]); }
80  );
81  auto tmpCopy = eigenValues;
82  for (unsigned int i = 0; i < numberOfElements; ++i)
83  {
84  eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
85  }
86  return indicesSortPermutations;
87  }
88 
97  template<typename QMatrix>
98  void permuteColumnsWithSortIndices( QMatrix & eigenVectors,
99  const std::vector<int> & indicesSortPermutations)
100  {
101  using EigenLibPermutationMatrix =
102  Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
103  auto numberOfElements = indicesSortPermutations.size();
104  // Creates a NxN permutation matrix copying our permutation to the matrix indices.
105  // Which holds the 1D array representation of a permutation.
106  EigenLibPermutationMatrix perm(numberOfElements);
107  perm.setIdentity();
108  std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(),
109  perm.indices().data());
110  // Apply it
111  eigenVectors = eigenVectors * perm;
112  }
113 } // end namespace detail
115 
151 template< typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix >
152 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
153 {
154 public:
155  typedef enum
156  {
157  OrderByValue = 1,
159  DoNotOrder
160  }
161  EigenValueOrderType;
162 
164  m_OrderEigenValues(OrderByValue) {}
165 
166  SymmetricEigenAnalysis(const unsigned int dimension):
167  m_UseEigenLibrary(false),
168  m_Dimension(dimension),
169  m_Order(dimension),
170  m_OrderEigenValues(OrderByValue) {}
171 
172  ~SymmetricEigenAnalysis() = default;
173 
174  using MatrixType = TMatrix;
175  using EigenMatrixType = TEigenMatrix;
176  using VectorType = TVector;
177 
191  unsigned int ComputeEigenValues(
192  const TMatrix & A,
193  TVector & EigenValues) const;
194 
215  unsigned int ComputeEigenValuesAndVectors(
216  const TMatrix & A,
217  TVector & EigenValues,
218  TEigenMatrix & EigenVectors) const;
219 
220 
222  void SetOrder(const unsigned int n)
223  {
224  m_Order = n;
225  }
226 
230  unsigned int GetOrder() const { return m_Order; }
231 
235  void SetOrderEigenValues(const bool b)
236  {
237  if ( b ) { m_OrderEigenValues = OrderByValue; }
238  else { m_OrderEigenValues = DoNotOrder; }
239  }
241 
242  bool GetOrderEigenValues() const { return ( m_OrderEigenValues == OrderByValue ); }
243 
247  void SetOrderEigenMagnitudes(const bool b)
248  {
249  if ( b ) { m_OrderEigenValues = OrderByMagnitude; }
250  else { m_OrderEigenValues = DoNotOrder; }
251  }
253 
254  bool GetOrderEigenMagnitudes() const { return ( m_OrderEigenValues == OrderByMagnitude ); }
255 
258  void SetDimension(const unsigned int n)
259  {
260  m_Dimension = n;
261  if ( m_Order == 0 )
262  {
263  m_Order = m_Dimension;
264  }
265  }
267 
270  unsigned int GetDimension() const { return m_Dimension; }
271 
273  void SetUseEigenLibrary(const bool input)
274  {
275  m_UseEigenLibrary = input;
276  }
278  {
279  m_UseEigenLibrary = true;
280  }
282  {
283  m_UseEigenLibrary = false;
284  }
285  bool GetUseEigenLibrary() const { return m_UseEigenLibrary; }
286 private:
287  bool m_UseEigenLibrary{false};
288  unsigned int m_Dimension{0};
289  unsigned int m_Order{0};
292 
312  void ReduceToTridiagonalMatrix(double *inputMatrix, double *d,
313  double *e, double *e2) const;
314 
336  void ReduceToTridiagonalMatrixAndGetTransformation(
337  double *inputMatrix, double *diagonalElements,
338  double *subDiagonalElements, double *transformMatrix) const;
339 
369  unsigned int ComputeEigenValuesUsingQL(double *d, double *e) const;
370 
408  unsigned int ComputeEigenValuesAndVectorsUsingQL(double *d, double *e, double *z) const;
409 
410  /* Legacy algorithms using thread-safe netlib.
411  * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
412  */
413  unsigned int ComputeEigenValuesLegacy(
414  const TMatrix & A,
415  TVector & EigenValues) const;
416 
417  unsigned int ComputeEigenValuesAndVectorsLegacy(
418  const TMatrix & A,
419  TVector & EigenValues,
420  TEigenMatrix & EigenVectors) const;
421 
422  /* Helper to get the matrix value type for EigenLibMatrix typename.
423  *
424  * If the TMatrix is vnl, the type is in element_type.
425  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
426  *
427  * To use this function:
428  * using ValueType = decltype(this->GetMatrixType(true));
429  *
430  * \note The two `GetMatrixValueType` overloads have different
431  * parameter declarations (`bool` and `...`), to avoid that both
432  * functions are equally good candidates during overload resolution,
433  * in case `element_type` and `ValueType` are both nested types of
434  * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
435  */
436  template<typename QMatrix = TMatrix >
437  auto GetMatrixValueType(bool) const
438  -> typename QMatrix::element_type
439  {
440  return QMatrix::element_type();
441  }
442  template<typename QMatrix = TMatrix >
443  auto GetMatrixValueType(...) const
444  -> typename QMatrix::ValueType
445  {
446  return QMatrix::ValueType();
447  }
448 
449  /* Wrapper that call the right implementation for the type of matrix. */
451  const TMatrix & A,
452  TVector & EigenValues,
453  TEigenMatrix & EigenVectors) const
454  {
455  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(
456  A, EigenValues, EigenVectors, true);
457  }
458 
459  /* Implementation detail using EigenLib that performs a copy of the input matrix.
460  *
461  * @param (long) implementation detail argument making this implementation less favourable
462  * to be chosen if alternatives are available.
463  *
464  * @return an unsigned int with no information value (no error code in EigenLib) */
465  template<typename QMatrix >
467  const QMatrix & A,
468  TVector & EigenValues,
469  TEigenMatrix & EigenVectors,
470  long) const
471  -> decltype(static_cast<unsigned int>(1))
472  {
473  using ValueType = decltype(GetMatrixValueType(true));
474  using EigenLibMatrixType =
475  Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
476  EigenLibMatrixType inputMatrix( m_Dimension, m_Dimension);
477  for (unsigned int row = 0; row < m_Dimension; ++row)
478  {
479  for (unsigned int col = 0; col < m_Dimension; ++col)
480  {
481  inputMatrix(row, col) = A(row, col);
482  }
483  }
484  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
485  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
486  const auto &eigenValues = solver.eigenvalues();
487  /* Column k of the returned matrix is an eigenvector corresponding to
488  * eigenvalue number $ k $ as returned by eigenvalues().
489  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
490  const auto &eigenVectors = solver.eigenvectors();
491 
492  if(m_OrderEigenValues == OrderByMagnitude)
493  {
494  auto copyEigenValues = eigenValues;
495  auto copyEigenVectors = eigenVectors;
496  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
497  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
498 
499  for (unsigned int row = 0; row < m_Dimension; ++row)
500  {
501  EigenValues[row] = copyEigenValues[row];
502  for (unsigned int col = 0; col < m_Dimension; ++col)
503  {
504  EigenVectors[row][col] = copyEigenVectors(col, row);
505  }
506  }
507  }
508  else
509  {
510  for (unsigned int row = 0; row < m_Dimension; ++row)
511  {
512  EigenValues[row] = eigenValues[row];
513  for (unsigned int col = 0; col < m_Dimension; ++col)
514  {
515  EigenVectors[row][col] = eigenVectors(col, row);
516  }
517  }
518  }
519  // No error code
520  return 1;
521  }
522 
523 
524  /* Implementation detail using EigenLib that do not peform a copy.
525  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
526  * If new types want to use this method, an appropiate overload of GetPointerToMatrixData
527  * should be included.
528  *
529  * @param (bool) implementation detail argument making this implementation the most favourable
530  * to be chosen from all the alternative implementations.
531  *
532  * @return an unsigned int with no information value (no error code in EigenLib) */
533  template<typename QMatrix >
535  const QMatrix & A,
536  TVector & EigenValues,
537  TEigenMatrix & EigenVectors,
538  bool) const
539  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
540  {
541  auto pointerToData = GetPointerToMatrixData(A);
542  using PointerType = decltype(pointerToData);
543  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
544  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
545  using EigenLibMatrixType =
546  Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
547  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
548  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
549  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
550  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
551  const auto & eigenValues = solver.eigenvalues();
552  /* Column k of the returned matrix is an eigenvector corresponding to
553  * eigenvalue number $ k $ as returned by eigenvalues().
554  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
555  const auto & eigenVectors = solver.eigenvectors();
556  if(m_OrderEigenValues == OrderByMagnitude)
557  {
558  auto copyEigenValues = eigenValues;
559  auto copyEigenVectors = eigenVectors;
560  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
561  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
562  for (unsigned int row = 0; row < m_Dimension; ++row)
563  {
564  EigenValues[row] = copyEigenValues[row];
565  for (unsigned int col = 0; col < m_Dimension; ++col)
566  {
567  EigenVectors[row][col] = copyEigenVectors(col, row);
568  }
569  }
570  }
571  else
572  {
573  for (unsigned int row = 0; row < m_Dimension; ++row)
574  {
575  EigenValues[row] = eigenValues[row];
576  for (unsigned int col = 0; col < m_Dimension; ++col)
577  {
578  EigenVectors[row][col] = eigenVectors(col, row);
579  }
580  }
581  }
582  // No error code
583  return 1;
584  }
585 
586  /* Wrapper that call the right implementation for the type of matrix. */
588  const TMatrix & A,
589  TVector & EigenValues) const
590  {
591  return ComputeEigenValuesWithEigenLibraryImpl(
592  A, EigenValues, true);
593  }
594 
595  /* Implementation detail using EigenLib that performs a copy of the input matrix.
596  *
597  * @param (long) implementation detail argument making this implementation less favourable
598  * to be chosen if alternatives are available.
599  *
600  * @return an unsigned int with no information value (no error code in EigenLib) */
601  template<typename QMatrix >
603  const QMatrix & A,
604  TVector & EigenValues,
605  long) const
606  -> decltype(static_cast<unsigned int>(1))
607  {
608  using ValueType = decltype(GetMatrixValueType(true));
609  using EigenLibMatrixType =
610  Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
611  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
612  for (unsigned int row = 0; row < m_Dimension; ++row)
613  {
614  for (unsigned int col = 0; col < m_Dimension; ++col)
615  {
616  inputMatrix(row, col) = A(row, col);
617  }
618  }
619  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
620  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
621  auto eigenValues = solver.eigenvalues();
622  if(m_OrderEigenValues == OrderByMagnitude)
623  {
624  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
625  }
626  for (unsigned int i = 0; i < m_Dimension; ++i)
627  {
628  EigenValues[i] = eigenValues[i];
629  }
630 
631  // No error code
632  return 1;
633  }
634 
635  /* Implementation detail using EigenLib that do not peform a copy.
636  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
637  * If new types want to use this method, an appropiate overload of GetPointerToMatrixData
638  * should be included.
639  *
640  * @param (bool) implementation detail argument making this implementation the most favourable
641  * to be chosen from all the alternative implementations.
642  *
643  * @return an unsigned int with no information value (no error code in EigenLib) */
644  template<typename QMatrix >
646  const QMatrix & A,
647  TVector & EigenValues,
648  bool) const
649  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
650  {
651  auto pointerToData = GetPointerToMatrixData(A);
652  using PointerType = decltype(pointerToData);
653  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
654  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
655  using EigenLibMatrixType =
656  Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
657  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
658  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
659  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
660  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
661  auto eigenValues = solver.eigenvalues();
662  if(m_OrderEigenValues == OrderByMagnitude)
663  {
664  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
665  }
666  for (unsigned int i = 0; i < m_Dimension; ++i)
667  {
668  EigenValues[i] = eigenValues[i];
669  }
670  // No error code
671  return 1;
672  }
673 };
674 
675 template< typename TMatrix, typename TVector, typename TEigenMatrix >
676 std::ostream & operator<<(std::ostream & os,
678 {
679  os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
680  os << " Dimension : " << s.GetDimension() << std::endl;
681  os << " Order : " << s.GetOrder() << std::endl;
682  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
683  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
684  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
685  return os;
686 }
687 
688 template< unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix >
689 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysisFixedDimension
690 {
691 public:
692  typedef enum
693  {
694  OrderByValue = 1,
696  DoNotOrder
697  }
698  EigenValueOrderType;
699 
700  SymmetricEigenAnalysisFixedDimension(): m_OrderEigenValues(OrderByValue) {}
702 
703  using MatrixType = TMatrix;
704  using EigenMatrixType = TEigenMatrix;
705  using VectorType = TVector;
706 
720  unsigned int ComputeEigenValues(
721  const TMatrix & A,
722  TVector & EigenValues) const
723  {
724  return ComputeEigenValuesWithEigenLibraryImpl(
725  A, EigenValues, true);
726  }
727 
749  const TMatrix & A,
750  TVector & EigenValues,
751  TEigenMatrix & EigenVectors) const
752  {
753  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(
754  A, EigenValues, EigenVectors, true);
755  }
756 
757  void SetOrderEigenValues(const bool b)
758  {
759  if ( b ) { m_OrderEigenValues = OrderByValue; }
760  else { m_OrderEigenValues = DoNotOrder; }
761  }
762  bool GetOrderEigenValues() const { return ( m_OrderEigenValues == OrderByValue ); }
763  void SetOrderEigenMagnitudes(const bool b)
764  {
765  if ( b ) { m_OrderEigenValues = OrderByMagnitude; }
766  else { m_OrderEigenValues = DoNotOrder; }
767  }
768  bool GetOrderEigenMagnitudes() const { return ( m_OrderEigenValues == OrderByMagnitude ); }
769  constexpr unsigned int GetOrder() const { return VDimension; }
770  constexpr unsigned int GetDimension() const { return VDimension; }
771  constexpr bool GetUseEigenLibrary() const { return true; }
772 
773 private:
775 
776  /* Helper to get the matrix value type for EigenLibMatrix typename.
777  *
778  * If the TMatrix is vnl, the type is in element_type.
779  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
780  *
781  * To use this function:
782  * using ValueType = decltype(this->GetMatrixType(true));
783  */
784  template<typename QMatrix = TMatrix >
785  auto GetMatrixValueType(bool) const
786  -> typename QMatrix::element_type
787  {
788  return QMatrix::element_type();
789  }
790  template<typename QMatrix = TMatrix >
791  auto GetMatrixValueType(bool) const
792  -> typename QMatrix::ValueType
793  {
794  return QMatrix::ValueType();
795  }
796 
797  /* Implementation detail using EigenLib that do not peform a copy.
798  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
799  * If new types want to use this method, an appropiate overload of GetPointerToMatrixData
800  * should be included.
801  *
802  * @param (bool) implementation detail argument making this implementation the most favourable
803  * to be chosen from all the alternative implementations.
804  *
805  * @return an unsigned int with no information value (no error code in EigenLib) */
806  template<typename QMatrix >
808  const QMatrix & A,
809  TVector & EigenValues,
810  TEigenMatrix & EigenVectors,
811  bool) const
812  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
813  {
814  auto pointerToData = GetPointerToMatrixData(A);
815  using PointerType = decltype(pointerToData);
816  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
817  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
818  using EigenLibMatrixType = Eigen::Matrix< ValueType, VDimension, VDimension, Eigen::RowMajor>;
819  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
820  EigenConstMatrixMap inputMatrix(pointerToData);
821  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
822  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
823  const auto & eigenValues = solver.eigenvalues();
824  /* Column k of the returned matrix is an eigenvector corresponding to
825  * eigenvalue number $ k $ as returned by eigenvalues().
826  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
827  const auto & eigenVectors = solver.eigenvectors();
828  if(m_OrderEigenValues == OrderByMagnitude)
829  {
830  auto copyEigenValues = eigenValues;
831  auto copyEigenVectors = eigenVectors;
832  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
833  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
834  for (unsigned int row = 0; row < VDimension; ++row)
835  {
836  EigenValues[row] = copyEigenValues[row];
837  for (unsigned int col = 0; col < VDimension; ++col)
838  {
839  EigenVectors[row][col] = copyEigenVectors(col, row);
840  }
841  }
842  }
843  else
844  {
845  for (unsigned int row = 0; row < VDimension; ++row)
846  {
847  EigenValues[row] = eigenValues[row];
848  for (unsigned int col = 0; col < VDimension; ++col)
849  {
850  EigenVectors[row][col] = eigenVectors(col, row);
851  }
852  }
853  }
854  // No error code
855  return 1;
856  }
857 
858  /* Implementation detail using EigenLib that performs a copy of the input matrix.
859  *
860  * @param (long) implementation detail argument making this implementation less favourable
861  * to be chosen if alternatives are available.
862  *
863  * @return an unsigned int with no information value (no error code in EigenLib) */
864  template<typename QMatrix >
866  const QMatrix & A,
867  TVector & EigenValues,
868  TEigenMatrix & EigenVectors,
869  long) const
870  -> decltype(static_cast<unsigned int>(1))
871  {
872  using ValueType = decltype(GetMatrixValueType(true));
873  using EigenLibMatrixType = Eigen::Matrix< ValueType, VDimension, VDimension, Eigen::RowMajor>;
874  EigenLibMatrixType inputMatrix;
875  for (unsigned int row = 0; row < VDimension; ++row)
876  {
877  for (unsigned int col = 0; col < VDimension; ++col)
878  {
879  inputMatrix(row, col) = A(row, col);
880  }
881  }
882  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
883  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
884  const auto &eigenValues = solver.eigenvalues();
885  /* Column k of the returned matrix is an eigenvector corresponding to
886  * eigenvalue number $ k $ as returned by eigenvalues().
887  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
888  const auto &eigenVectors = solver.eigenvectors();
889 
890  if(m_OrderEigenValues == OrderByMagnitude)
891  {
892  auto copyEigenValues = eigenValues;
893  auto copyEigenVectors = eigenVectors;
894  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
895  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
896 
897  for (unsigned int row = 0; row < VDimension; ++row)
898  {
899  EigenValues[row] = copyEigenValues[row];
900  for (unsigned int col = 0; col < VDimension; ++col)
901  {
902  EigenVectors[row][col] = copyEigenVectors(col, row);
903  }
904  }
905  }
906  else
907  {
908  for (unsigned int row = 0; row < VDimension; ++row)
909  {
910  EigenValues[row] = eigenValues[row];
911  for (unsigned int col = 0; col < VDimension; ++col)
912  {
913  EigenVectors[row][col] = eigenVectors(col, row);
914  }
915  }
916  }
917  // No error code
918  return 1;
919  }
920 
921  /* Implementation detail using EigenLib that performs a copy of the input matrix.
922  *
923  * @param (long) implementation detail argument making this implementation less favourable
924  * to be chosen if alternatives are available.
925  *
926  * @return an unsigned int with no information value (no error code in EigenLib) */
927  template<typename QMatrix >
929  const QMatrix & A,
930  TVector & EigenValues,
931  long) const
932  -> decltype(static_cast<unsigned int>(1))
933  {
934  using ValueType = decltype(GetMatrixValueType(true));
935  using EigenLibMatrixType = Eigen::Matrix< ValueType, VDimension, VDimension, Eigen::RowMajor>;
936  EigenLibMatrixType inputMatrix;
937  for (unsigned int row = 0; row < VDimension; ++row)
938  {
939  for (unsigned int col = 0; col < VDimension; ++col)
940  {
941  inputMatrix(row, col) = A(row, col);
942  }
943  }
944  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
945  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
946  auto eigenValues = solver.eigenvalues();
947  if(m_OrderEigenValues == OrderByMagnitude)
948  {
949  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
950  }
951  for (unsigned int i = 0; i < VDimension; ++i)
952  {
953  EigenValues[i] = eigenValues[i];
954  }
955 
956  // No error code
957  return 1;
958  }
959 
960  /* Implementation detail using EigenLib that do not peform a copy.
961  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
962  * If new types want to use this method, an appropiate overload of GetPointerToMatrixData
963  * should be included.
964  *
965  * @param (bool) implementation detail argument making this implementation the most favourable
966  * to be chosen from all the alternative implementations.
967  *
968  * @return an unsigned int with no information value (no error code in EigenLib) */
969  template<typename QMatrix >
971  const QMatrix & A,
972  TVector & EigenValues,
973  bool) const
974  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
975  {
976  auto pointerToData = GetPointerToMatrixData(A);
977  using PointerType = decltype(pointerToData);
978  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
979  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
980  using EigenLibMatrixType = Eigen::Matrix< ValueType, VDimension, VDimension, Eigen::RowMajor>;
981  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
982  EigenConstMatrixMap inputMatrix(pointerToData);
983  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
984  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
985  auto eigenValues = solver.eigenvalues();
986  if(m_OrderEigenValues == OrderByMagnitude)
987  {
988  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
989  }
990  for (unsigned int i = 0; i < VDimension; ++i)
991  {
992  EigenValues[i] = eigenValues[i];
993  }
994  // No error code
995  return 1;
996  }
997 };
998 
999 template< unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix >
1000 std::ostream & operator<<(std::ostream & os,
1001  const SymmetricEigenAnalysisFixedDimension<VDimension,
1002  TMatrix, TVector, TEigenMatrix > & s)
1003 {
1004  os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1005  os << " Dimension : " << s.GetDimension() << std::endl;
1006  os << " Order : " << s.GetOrder() << std::endl;
1007  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
1008  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
1009  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
1010  return os;
1011 }
1012 } // end namespace itk
1013 
1014 #ifndef ITK_MANUAL_INSTANTIATION
1015 #include "itkSymmetricEigenAnalysis.hxx"
1016 #endif
1017 
1018 #endif
A templated class holding a M x N size Matrix.
Definition: itkMatrix.h:49
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(static_cast< unsigned int >(1))
Find Eigen values of a real 2D symmetric matrix. It serves as a thread-safe alternative to the class:...
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
std::ostream & operator<<(std::ostream &os, const Array< TValue > &arr)
Definition: itkArray.h:188
void SetOrder(const unsigned int n)
const TValueType * GetPointerToMatrixData(const vnl_matrix_fixed< TValueType, NRows, NCols > &inputMatrix)
unsigned int ComputeEigenValues(const TMatrix &A, TVector &EigenValues) const
void SetUseEigenLibrary(const bool input)
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(static_cast< unsigned int >(1))
auto GetMatrixValueType(bool) const -> typename QMatrix::ValueType
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(static_cast< unsigned int >(1))
auto GetMatrixValueType(...) const -> typename QMatrix::ValueType
void permuteColumnsWithSortIndices(QMatrix &eigenVectors, const std::vector< int > &indicesSortPermutations)
unsigned int ComputeEigenValuesAndVectorsWithEigenLibrary(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
std::vector< int > sortEigenValuesByMagnitude(TArray &eigenValues, const unsigned int numberOfElements)
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
void SetDimension(const unsigned int n)
InternalMatrixType & GetVnlMatrix()
Definition: itkMatrix.h:168
static constexpr double e
The base of the natural logarithm or Euler&#39;s number
Definition: itkMath.h:53
unsigned int ComputeEigenValuesWithEigenLibrary(const TMatrix &A, TVector &EigenValues) const
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(static_cast< unsigned int >(1))
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
SymmetricEigenAnalysis(const unsigned int dimension)
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))