ITK  5.2.0
Insight Toolkit
itkSymmetricEigenAnalysis.h
Go to the documentation of this file.
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  * 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 *
38 GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, NRows, NCols> & inputMatrix)
39 {
40  return inputMatrix.data_block();
41 };
42 template <typename TValueType>
43 const TValueType *
44 GetPointerToMatrixData(const vnl_matrix<TValueType> & inputMatrix)
45 {
46  return inputMatrix.data_block();
47 };
48 
49 template <typename TValueType, unsigned int NRows, unsigned int NCols>
50 const TValueType *
52 {
53  return inputMatrix.GetVnlMatrix().data_block();
54 };
55 
71 template <typename TArray>
72 std::vector<int>
73 sortEigenValuesByMagnitude(TArray & eigenValues, const unsigned int numberOfElements)
74 {
75  std::vector<int> indicesSortPermutations(numberOfElements, 0);
76  std::iota(std::begin(indicesSortPermutations), std::end(indicesSortPermutations), 0);
78 
79  std::sort(
80  std::begin(indicesSortPermutations),
81  std::end(indicesSortPermutations),
82  [&eigenValues](unsigned int a, unsigned int b) { return std::abs(eigenValues[a]) < std::abs(eigenValues[b]); });
83  auto tmpCopy = eigenValues;
84  for (unsigned int i = 0; i < numberOfElements; ++i)
85  {
86  eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
87  }
88  return indicesSortPermutations;
89 }
90 
99 template <typename QMatrix>
100 void
101 permuteColumnsWithSortIndices(QMatrix & eigenVectors, const std::vector<int> & indicesSortPermutations)
102 {
103  using EigenLibPermutationMatrix = Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
104  auto numberOfElements = indicesSortPermutations.size();
105  // Creates a NxN permutation matrix copying our permutation to the matrix indices.
106  // Which holds the 1D array representation of a permutation.
107  EigenLibPermutationMatrix perm(numberOfElements);
108  perm.setIdentity();
109  std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(), perm.indices().data());
110  // Apply it
111  eigenVectors = eigenVectors * perm;
112 }
113 } // end namespace detail
115 
121 {
122 public:
131  {
132  OrderByValue = 1,
133  OrderByMagnitude = 2,
134  DoNotOrder = 3
135  };
136 };
137 // Define how to print enumeration
138 extern ITKCommon_EXPORT std::ostream &
139  operator<<(std::ostream & out, const SymmetricEigenAnalysisEnums::EigenValueOrder value);
140 
142 
143 inline EigenValueOrderEnum
145 {
146  switch (value)
147  {
148  case 1:
150  case 2:
152  case 3:
154  default:
155  break;
156  }
157  itkGenericExceptionMacro(<< "Error: Invalid value for conversion.");
158 }
159 
160 #if !defined(ITK_LEGACY_REMOVE)
161 
162 static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
163 static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
164 static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
165 #endif
166 
202 template <typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
203 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
204 {
205 public:
207 #if !defined(ITK_LEGACY_REMOVE)
208 
209  using EigenValueOrderType = EigenValueOrderEnum;
210 #endif
211 
212  SymmetricEigenAnalysis() = default;
213 
214  SymmetricEigenAnalysis(const unsigned int dimension)
215  : m_Dimension(dimension)
216  , m_Order(dimension)
217  {}
218 
219  ~SymmetricEigenAnalysis() = default;
220 
221  using MatrixType = TMatrix;
222  using EigenMatrixType = TEigenMatrix;
223  using VectorType = TVector;
224 
238  unsigned int
239  ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const;
240 
261  unsigned int
262  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
263 
264 
266  void
267  SetOrder(const unsigned int n)
268  {
269  m_Order = n;
270  }
271 
275  unsigned int
276  GetOrder() const
277  {
278  return m_Order;
279  }
280 
284  void
285  SetOrderEigenValues(const bool b)
286  {
287  if (b)
288  {
289  m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
290  }
291  else
292  {
293  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
294  }
295  }
297 
298  bool
300  {
301  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
302  }
303 
307  void
309  {
310  if (b)
311  {
312  m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
313  }
314  else
315  {
316  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
317  }
318  }
320 
321  bool
323  {
324  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
325  }
326 
329  void
330  SetDimension(const unsigned int n)
331  {
332  m_Dimension = n;
333  if (m_Order == 0)
334  {
335  m_Order = m_Dimension;
336  }
337  }
339 
342  unsigned int
343  GetDimension() const
344  {
345  return m_Dimension;
346  }
347 
349  void
350  SetUseEigenLibrary(const bool input)
351  {
352  m_UseEigenLibrary = input;
353  }
354  void
356  {
357  m_UseEigenLibrary = true;
358  }
359  void
361  {
362  m_UseEigenLibrary = false;
363  }
364  bool
366  {
367  return m_UseEigenLibrary;
368  }
370 
371 private:
372  bool m_UseEigenLibrary{ false };
373  unsigned int m_Dimension{ 0 };
374  unsigned int m_Order{ 0 };
376 
396  void
397  ReduceToTridiagonalMatrix(double * inputMatrix, double * d, double * e, double * e2) const;
398 
420  void
421  ReduceToTridiagonalMatrixAndGetTransformation(double * inputMatrix,
422  double * diagonalElements,
423  double * subDiagonalElements,
424  double * transformMatrix) const;
425 
455  unsigned int
456  ComputeEigenValuesUsingQL(double * d, double * e) const;
457 
495  unsigned int
496  ComputeEigenValuesAndVectorsUsingQL(double * d, double * e, double * z) const;
497 
498  /* Legacy algorithms using thread-safe netlib.
499  * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
500  */
501  unsigned int
502  ComputeEigenValuesLegacy(const TMatrix & A, TVector & EigenValues) const;
503 
504  unsigned int
505  ComputeEigenValuesAndVectorsLegacy(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
506 
507  /* Helper to get the matrix value type for EigenLibMatrix typename.
508  *
509  * If the TMatrix is vnl, the type is in element_type.
510  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
511  *
512  * To use this function:
513  * using ValueType = decltype(this->GetMatrixType(true));
514  *
515  * \note The two `GetMatrixValueType` overloads have different
516  * parameter declarations (`bool` and `...`), to avoid that both
517  * functions are equally good candidates during overload resolution,
518  * in case `element_type` and `ValueType` are both nested types of
519  * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
520  */
521  template <typename QMatrix = TMatrix>
522  auto
523  GetMatrixValueType(bool) const -> typename QMatrix::element_type
524  {
525  return QMatrix::element_type();
526  }
527  template <typename QMatrix = TMatrix>
528  auto
529  GetMatrixValueType(...) const -> typename QMatrix::ValueType
530  {
531  return QMatrix::ValueType();
532  }
533 
534  /* Wrapper that call the right implementation for the type of matrix. */
535  unsigned int
537  TVector & EigenValues,
538  TEigenMatrix & EigenVectors) const
539  {
540  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
541  }
542 
543  /* Implementation detail using EigenLib that performs a copy of the input matrix.
544  *
545  * @param (long) implementation detail argument making this implementation less favourable
546  * to be chosen if alternatives are available.
547  *
548  * @return an unsigned int with no information value (no error code in EigenLib) */
549  template <typename QMatrix>
550  auto
552  TVector & EigenValues,
553  TEigenMatrix & EigenVectors,
554  long) const -> decltype(static_cast<unsigned int>(1))
555  {
556  using ValueType = decltype(GetMatrixValueType(true));
557  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
558  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
559  for (unsigned int row = 0; row < m_Dimension; ++row)
560  {
561  for (unsigned int col = 0; col < m_Dimension; ++col)
562  {
563  inputMatrix(row, col) = A(row, col);
564  }
565  }
566  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
567  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
568  const auto & eigenValues = solver.eigenvalues();
569  /* Column k of the returned matrix is an eigenvector corresponding to
570  * eigenvalue number $ k $ as returned by eigenvalues().
571  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
572  const auto & eigenVectors = solver.eigenvectors();
573 
574  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
575  {
576  auto copyEigenValues = eigenValues;
577  auto copyEigenVectors = eigenVectors;
578  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
579  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
580 
581  for (unsigned int row = 0; row < m_Dimension; ++row)
582  {
583  EigenValues[row] = copyEigenValues[row];
584  for (unsigned int col = 0; col < m_Dimension; ++col)
585  {
586  EigenVectors[row][col] = copyEigenVectors(col, row);
587  }
588  }
589  }
590  else
591  {
592  for (unsigned int row = 0; row < m_Dimension; ++row)
593  {
594  EigenValues[row] = eigenValues[row];
595  for (unsigned int col = 0; col < m_Dimension; ++col)
596  {
597  EigenVectors[row][col] = eigenVectors(col, row);
598  }
599  }
600  }
601  // No error code
602  return 1;
603  }
604 
605 
606  /* Implementation detail using EigenLib that do not peform a copy.
607  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
608  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
609  * should be included.
610  *
611  * @param (bool) implementation detail argument making this implementation the most favourable
612  * to be chosen from all the alternative implementations.
613  *
614  * @return an unsigned int with no information value (no error code in EigenLib) */
615  template <typename QMatrix>
616  auto
618  TVector & EigenValues,
619  TEigenMatrix & EigenVectors,
620  bool) const
621  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
622  {
623  auto pointerToData = GetPointerToMatrixData(A);
624  using PointerType = decltype(pointerToData);
625  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
626  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
627  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
628  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
629  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
630  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
631  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
632  const auto & eigenValues = solver.eigenvalues();
633  /* Column k of the returned matrix is an eigenvector corresponding to
634  * eigenvalue number $ k $ as returned by eigenvalues().
635  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
636  const auto & eigenVectors = solver.eigenvectors();
637  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
638  {
639  auto copyEigenValues = eigenValues;
640  auto copyEigenVectors = eigenVectors;
641  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
642  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
643  for (unsigned int row = 0; row < m_Dimension; ++row)
644  {
645  EigenValues[row] = copyEigenValues[row];
646  for (unsigned int col = 0; col < m_Dimension; ++col)
647  {
648  EigenVectors[row][col] = copyEigenVectors(col, row);
649  }
650  }
651  }
652  else
653  {
654  for (unsigned int row = 0; row < m_Dimension; ++row)
655  {
656  EigenValues[row] = eigenValues[row];
657  for (unsigned int col = 0; col < m_Dimension; ++col)
658  {
659  EigenVectors[row][col] = eigenVectors(col, row);
660  }
661  }
662  }
663  // No error code
664  return 1;
665  }
666 
667  /* Wrapper that call the right implementation for the type of matrix. */
668  unsigned int
669  ComputeEigenValuesWithEigenLibrary(const TMatrix & A, TVector & EigenValues) const
670  {
671  return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
672  }
673 
674  /* Implementation detail using EigenLib that performs a copy of the input matrix.
675  *
676  * @param (long) implementation detail argument making this implementation less favourable
677  * to be chosen if alternatives are available.
678  *
679  * @return an unsigned int with no information value (no error code in EigenLib) */
680  template <typename QMatrix>
681  auto
682  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const
683  -> decltype(static_cast<unsigned int>(1))
684  {
685  using ValueType = decltype(GetMatrixValueType(true));
686  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
687  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
688  for (unsigned int row = 0; row < m_Dimension; ++row)
689  {
690  for (unsigned int col = 0; col < m_Dimension; ++col)
691  {
692  inputMatrix(row, col) = A(row, col);
693  }
694  }
695  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
696  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
697  auto eigenValues = solver.eigenvalues();
698  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
699  {
700  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
701  }
702  for (unsigned int i = 0; i < m_Dimension; ++i)
703  {
704  EigenValues[i] = eigenValues[i];
705  }
706 
707  // No error code
708  return 1;
709  }
710 
711  /* Implementation detail using EigenLib that do not peform a copy.
712  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
713  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
714  * should be included.
715  *
716  * @param (bool) implementation detail argument making this implementation the most favourable
717  * to be chosen from all the alternative implementations.
718  *
719  * @return an unsigned int with no information value (no error code in EigenLib) */
720  template <typename QMatrix>
721  auto
722  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
723  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
724  {
725  auto pointerToData = GetPointerToMatrixData(A);
726  using PointerType = decltype(pointerToData);
727  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
728  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
729  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
730  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
731  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
732  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
733  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
734  auto eigenValues = solver.eigenvalues();
735  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
736  {
737  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
738  }
739  for (unsigned int i = 0; i < m_Dimension; ++i)
740  {
741  EigenValues[i] = eigenValues[i];
742  }
743  // No error code
744  return 1;
745  }
746 };
747 
748 template <typename TMatrix, typename TVector, typename TEigenMatrix>
749 std::ostream &
751 {
752  os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
753  os << " Dimension : " << s.GetDimension() << std::endl;
754  os << " Order : " << s.GetOrder() << std::endl;
755  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
756  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
757  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
758  return os;
759 }
760 
761 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
762 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysisFixedDimension
763 {
764 public:
765 #if !defined(ITK_LEGACY_REMOVE)
766 
767  using EigenValueOrderType = EigenValueOrderEnum;
768 #endif
769 #if !defined(ITK_LEGACY_REMOVE)
770  // We need to expose the enum values at the class level
771  // for backwards compatibility
772  static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
773  static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
774  static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
775 #endif
776 
779 
780  using MatrixType = TMatrix;
781  using EigenMatrixType = TEigenMatrix;
782  using VectorType = TVector;
783 
797  unsigned int
798  ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const
799  {
800  return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
801  }
802 
823  unsigned int
824  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const
825  {
826  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
827  }
828 
829  void
830  SetOrderEigenValues(const bool b)
831  {
832  if (b)
833  {
834  m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
835  }
836  else
837  {
838  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
839  }
840  }
841  bool
843  {
844  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
845  }
846  void
848  {
849  if (b)
850  {
851  m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
852  }
853  else
854  {
855  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
856  }
857  }
858  bool
860  {
861  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
862  }
863  constexpr unsigned int
864  GetOrder() const
865  {
866  return VDimension;
867  }
868  constexpr unsigned int
869  GetDimension() const
870  {
871  return VDimension;
872  }
873  constexpr bool
875  {
876  return true;
877  }
878 
879 private:
881 
882  /* Helper to get the matrix value type for EigenLibMatrix typename.
883  *
884  * If the TMatrix is vnl, the type is in element_type.
885  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
886  *
887  * To use this function:
888  * using ValueType = decltype(this->GetMatrixType(true));
889  */
890  template <typename QMatrix = TMatrix>
891  auto
892  GetMatrixValueType(bool) const -> typename QMatrix::element_type
893  {
894  return QMatrix::element_type();
895  }
896  template <typename QMatrix = TMatrix>
897  auto
898  GetMatrixValueType(bool) const -> typename QMatrix::ValueType
899  {
900  return QMatrix::ValueType();
901  }
902 
903  /* Implementation detail using EigenLib that do not peform a copy.
904  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
905  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
906  * should be included.
907  *
908  * @param (bool) implementation detail argument making this implementation the most favourable
909  * to be chosen from all the alternative implementations.
910  *
911  * @return an unsigned int with no information value (no error code in EigenLib) */
912  template <typename QMatrix>
913  auto
915  TVector & EigenValues,
916  TEigenMatrix & EigenVectors,
917  bool) const
918  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
919  {
920  auto pointerToData = GetPointerToMatrixData(A);
921  using PointerType = decltype(pointerToData);
922  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
923  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
924  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
925  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
926  EigenConstMatrixMap inputMatrix(pointerToData);
927  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
928  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
929  const auto & eigenValues = solver.eigenvalues();
930  /* Column k of the returned matrix is an eigenvector corresponding to
931  * eigenvalue number $ k $ as returned by eigenvalues().
932  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
933  const auto & eigenVectors = solver.eigenvectors();
934  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
935  {
936  auto copyEigenValues = eigenValues;
937  auto copyEigenVectors = eigenVectors;
938  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
939  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
940  for (unsigned int row = 0; row < VDimension; ++row)
941  {
942  EigenValues[row] = copyEigenValues[row];
943  for (unsigned int col = 0; col < VDimension; ++col)
944  {
945  EigenVectors[row][col] = copyEigenVectors(col, row);
946  }
947  }
948  }
949  else
950  {
951  for (unsigned int row = 0; row < VDimension; ++row)
952  {
953  EigenValues[row] = eigenValues[row];
954  for (unsigned int col = 0; col < VDimension; ++col)
955  {
956  EigenVectors[row][col] = eigenVectors(col, row);
957  }
958  }
959  }
960  // No error code
961  return 1;
962  }
963 
964  /* Implementation detail using EigenLib that performs a copy of the input matrix.
965  *
966  * @param (long) implementation detail argument making this implementation less favourable
967  * to be chosen if alternatives are available.
968  *
969  * @return an unsigned int with no information value (no error code in EigenLib) */
970  template <typename QMatrix>
971  auto
973  TVector & EigenValues,
974  TEigenMatrix & EigenVectors,
975  long) const -> decltype(static_cast<unsigned int>(1))
976  {
977  using ValueType = decltype(GetMatrixValueType(true));
978  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
979  EigenLibMatrixType inputMatrix;
980  for (unsigned int row = 0; row < VDimension; ++row)
981  {
982  for (unsigned int col = 0; col < VDimension; ++col)
983  {
984  inputMatrix(row, col) = A(row, col);
985  }
986  }
987  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
988  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
989  const auto & eigenValues = solver.eigenvalues();
990  /* Column k of the returned matrix is an eigenvector corresponding to
991  * eigenvalue number $ k $ as returned by eigenvalues().
992  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
993  const auto & eigenVectors = solver.eigenvectors();
994 
995  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
996  {
997  auto copyEigenValues = eigenValues;
998  auto copyEigenVectors = eigenVectors;
999  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
1000  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
1001 
1002  for (unsigned int row = 0; row < VDimension; ++row)
1003  {
1004  EigenValues[row] = copyEigenValues[row];
1005  for (unsigned int col = 0; col < VDimension; ++col)
1006  {
1007  EigenVectors[row][col] = copyEigenVectors(col, row);
1008  }
1009  }
1010  }
1011  else
1012  {
1013  for (unsigned int row = 0; row < VDimension; ++row)
1014  {
1015  EigenValues[row] = eigenValues[row];
1016  for (unsigned int col = 0; col < VDimension; ++col)
1017  {
1018  EigenVectors[row][col] = eigenVectors(col, row);
1019  }
1020  }
1021  }
1022  // No error code
1023  return 1;
1024  }
1025 
1026  /* Implementation detail using EigenLib that performs a copy of the input matrix.
1027  *
1028  * @param (long) implementation detail argument making this implementation less favourable
1029  * to be chosen if alternatives are available.
1030  *
1031  * @return an unsigned int with no information value (no error code in EigenLib) */
1032  template <typename QMatrix>
1033  auto
1034  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const
1035  -> decltype(static_cast<unsigned int>(1))
1036  {
1037  using ValueType = decltype(GetMatrixValueType(true));
1038  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1039  EigenLibMatrixType inputMatrix;
1040  for (unsigned int row = 0; row < VDimension; ++row)
1041  {
1042  for (unsigned int col = 0; col < VDimension; ++col)
1043  {
1044  inputMatrix(row, col) = A(row, col);
1045  }
1046  }
1047  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1048  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1049  auto eigenValues = solver.eigenvalues();
1050  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1051  {
1052  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1053  }
1054  for (unsigned int i = 0; i < VDimension; ++i)
1055  {
1056  EigenValues[i] = eigenValues[i];
1057  }
1058 
1059  // No error code
1060  return 1;
1061  }
1062 
1063  /* Implementation detail using EigenLib that do not peform a copy.
1064  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
1065  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
1066  * should be included.
1067  *
1068  * @param (bool) implementation detail argument making this implementation the most favourable
1069  * to be chosen from all the alternative implementations.
1070  *
1071  * @return an unsigned int with no information value (no error code in EigenLib) */
1072  template <typename QMatrix>
1073  auto
1074  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
1075  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
1076  {
1077  auto pointerToData = GetPointerToMatrixData(A);
1078  using PointerType = decltype(pointerToData);
1079  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
1080  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
1081  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1082  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1083  EigenConstMatrixMap inputMatrix(pointerToData);
1084  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1085  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1086  auto eigenValues = solver.eigenvalues();
1087  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1088  {
1089  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1090  }
1091  for (unsigned int i = 0; i < VDimension; ++i)
1092  {
1093  EigenValues[i] = eigenValues[i];
1094  }
1095  // No error code
1096  return 1;
1097  }
1098 };
1099 
1100 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix>
1101 std::ostream &
1102 operator<<(std::ostream & os,
1104 {
1105  os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1106  os << " Dimension : " << s.GetDimension() << std::endl;
1107  os << " Order : " << s.GetOrder() << std::endl;
1108  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
1109  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
1110  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
1111  return os;
1112 }
1113 } // end namespace itk
1114 
1115 #ifndef ITK_MANUAL_INSTANTIATION
1116 # include "itkSymmetricEigenAnalysis.hxx"
1117 #endif
1118 
1119 #endif
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::VectorType
TOutput VectorType
Definition: itkSymmetricEigenAnalysis.h:782
itk::uint8_t
::uint8_t uint8_t
Definition: itkIntTypes.h:29
EigenValueOrder
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectors
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:824
itk::SymmetricEigenAnalysis::SetDimension
void SetDimension(const unsigned int n)
Definition: itkSymmetricEigenAnalysis.h:330
itk::SymmetricEigenAnalysis::GetOrder
unsigned int GetOrder() const
Definition: itkSymmetricEigenAnalysis.h:276
itk::Matrix::GetVnlMatrix
InternalMatrixType & GetVnlMatrix()
Definition: itkMatrix.h:178
itk::SymmetricEigenAnalysisEnums::EigenValueOrder::OrderByValue
itkMatrix.h
itk::operator<<
std::ostream & operator<<(std::ostream &os, const Array< TValue > &arr)
Definition: itkArray.h:213
itk::SymmetricEigenAnalysis::GetUseEigenLibrary
bool GetUseEigenLibrary() const
Definition: itkSymmetricEigenAnalysis.h:365
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:551
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:617
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:722
itk::SymmetricEigenAnalysis::SetUseEigenLibraryOn
void SetUseEigenLibraryOn()
Definition: itkSymmetricEigenAnalysis.h:355
itk::SymmetricEigenAnalysisEnums::EigenValueOrder::OrderByMagnitude
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::VectorType
TOutputImage::PixelType VectorType
Definition: itkSymmetricEigenAnalysis.h:223
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:914
itk::SymmetricEigenAnalysis::SymmetricEigenAnalysis
SymmetricEigenAnalysis(const unsigned int dimension)
Definition: itkSymmetricEigenAnalysis.h:214
itk::SymmetricEigenAnalysisFixedDimension::SetOrderEigenValues
void SetOrderEigenValues(const bool b)
Definition: itkSymmetricEigenAnalysis.h:830
itk::detail::permuteColumnsWithSortIndices
void permuteColumnsWithSortIndices(QMatrix &eigenVectors, const std::vector< int > &indicesSortPermutations)
Definition: itkSymmetricEigenAnalysis.h:101
itk::SymmetricEigenAnalysis::GetOrderEigenValues
bool GetOrderEigenValues() const
Definition: itkSymmetricEigenAnalysis.h:299
itk::SymmetricEigenAnalysisFixedDimension::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::ValueType
Definition: itkSymmetricEigenAnalysis.h:898
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
Definition: itkSymmetricEigenAnalysis.h:523
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibrary
unsigned int ComputeEigenValuesAndVectorsWithEigenLibrary(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:536
itk::EigenValueOrderEnum
SymmetricEigenAnalysisEnums::EigenValueOrder EigenValueOrderEnum
Definition: itkSymmetricEigenAnalysis.h:141
itk::detail::GetPointerToMatrixData
const TValueType * GetPointerToMatrixData(const vnl_matrix_fixed< TValueType, NRows, NCols > &inputMatrix)
Definition: itkSymmetricEigenAnalysis.h:38
itk::SymmetricEigenAnalysis::SetOrder
void SetOrder(const unsigned int n)
Definition: itkSymmetricEigenAnalysis.h:267
itkMacro.h
itk::SymmetricEigenAnalysis::GetDimension
unsigned int GetDimension() const
Definition: itkSymmetricEigenAnalysis.h:343
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:1034
itk::Int2EigenValueOrderEnum
EigenValueOrderEnum Int2EigenValueOrderEnum(const uint8_t value)
Definition: itkSymmetricEigenAnalysis.h:144
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(...) const -> typename QMatrix::ValueType
Definition: itkSymmetricEigenAnalysis.h:529
itk::detail::sortEigenValuesByMagnitude
std::vector< int > sortEigenValuesByMagnitude(TArray &eigenValues, const unsigned int numberOfElements)
Definition: itkSymmetricEigenAnalysis.h:73
itk::SymmetricEigenAnalysisFixedDimension::SetOrderEigenMagnitudes
void SetOrderEigenMagnitudes(const bool b)
Definition: itkSymmetricEigenAnalysis.h:847
itk::SymmetricEigenAnalysisFixedDimension::GetOrderEigenMagnitudes
bool GetOrderEigenMagnitudes() const
Definition: itkSymmetricEigenAnalysis.h:859
itk::SymmetricEigenAnalysis::SetOrderEigenMagnitudes
void SetOrderEigenMagnitudes(const bool b)
Definition: itkSymmetricEigenAnalysis.h:308
itk::SymmetricEigenAnalysisEnums::EigenValueOrder::DoNotOrder
itk::SymmetricEigenAnalysis::SetUseEigenLibraryOff
void SetUseEigenLibraryOff()
Definition: itkSymmetricEigenAnalysis.h:360
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::EigenMatrixType
TInput EigenMatrixType
Definition: itkSymmetricEigenAnalysis.h:781
itk::Matrix
A templated class holding a M x N size Matrix.
Definition: itkMatrix.h:51
itk::SymmetricEigenAnalysisFixedDimension::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
Definition: itkSymmetricEigenAnalysis.h:892
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:1074
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::EigenMatrixType
TInputImage::PixelType EigenMatrixType
Definition: itkSymmetricEigenAnalysis.h:222
itk::SymmetricEigenAnalysis::SetUseEigenLibrary
void SetUseEigenLibrary(const bool input)
Definition: itkSymmetricEigenAnalysis.h:350
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:972
itk::SymmetricEigenAnalysisFixedDimension::GetDimension
constexpr unsigned int GetDimension() const
Definition: itkSymmetricEigenAnalysis.h:869
itk::SymmetricEigenAnalysisFixedDimension
Definition: itkSymmetricEigenAnalysis.h:762
itk::SymmetricEigenAnalysisFixedDimension::GetUseEigenLibrary
constexpr bool GetUseEigenLibrary() const
Definition: itkSymmetricEigenAnalysis.h:874
itk::SymmetricEigenAnalysis::SetOrderEigenValues
void SetOrderEigenValues(const bool b)
Definition: itkSymmetricEigenAnalysis.h:285
itk
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
Definition: itkAnnulusOperator.h:24
itk::SymmetricEigenAnalysisFixedDimension::GetOrderEigenValues
bool GetOrderEigenValues() const
Definition: itkSymmetricEigenAnalysis.h:842
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:682
itk::Math::e
static constexpr double e
Definition: itkMath.h:54
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibrary
unsigned int ComputeEigenValuesWithEigenLibrary(const TMatrix &A, TVector &EigenValues) const
Definition: itkSymmetricEigenAnalysis.h:669
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::MatrixType
TInput MatrixType
Definition: itkSymmetricEigenAnalysis.h:780
itk::SymmetricEigenAnalysisFixedDimension::GetOrder
constexpr unsigned int GetOrder() const
Definition: itkSymmetricEigenAnalysis.h:864
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::MatrixType
TInputImage::PixelType MatrixType
Definition: itkSymmetricEigenAnalysis.h:221
itk::SymmetricEigenAnalysisEnums
This class contains all enum classes used by SymmetricEigenAnalysis class.
Definition: itkSymmetricEigenAnalysis.h:120
itk::detail::GetPointerToMatrixData
const TValueType * GetPointerToMatrixData(const itk::Matrix< TValueType, NRows, NCols > &inputMatrix)
Definition: itkSymmetricEigenAnalysis.h:51
itk::SymmetricEigenAnalysis::GetOrderEigenMagnitudes
bool GetOrderEigenMagnitudes() const
Definition: itkSymmetricEigenAnalysis.h:322
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValues
unsigned int ComputeEigenValues(const TMatrix &A, TVector &EigenValues) const
Definition: itkSymmetricEigenAnalysis.h:798
itk::SymmetricEigenAnalysisEnums::EigenValueOrder
EigenValueOrder
Definition: itkSymmetricEigenAnalysis.h:130
itk::SymmetricEigenAnalysis
Find Eigen values of a real 2D symmetric matrix. It serves as a thread-safe alternative to the class:...
Definition: itkSymmetricEigenAnalysis.h:203