ITK  5.3.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 & D) 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 * a, double * d, double * e, double * e2) const;
398 
420  void
421  ReduceToTridiagonalMatrixAndGetTransformation(const double * a, double * d, double * e, double * z) const;
422 
452  unsigned int
453  ComputeEigenValuesUsingQL(double * d, double * e) const;
454 
492  unsigned int
493  ComputeEigenValuesAndVectorsUsingQL(double * d, double * e, double * z) const;
494 
495  /* Legacy algorithms using thread-safe netlib.
496  * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
497  */
498  unsigned int
499  ComputeEigenValuesLegacy(const TMatrix & A, TVector & D) const;
500 
501  unsigned int
502  ComputeEigenValuesAndVectorsLegacy(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
503 
504  /* Helper to get the matrix value type for EigenLibMatrix typename.
505  *
506  * If the TMatrix is vnl, the type is in element_type.
507  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
508  *
509  * To use this function:
510  * using ValueType = decltype(this->GetMatrixType(true));
511  *
512  * \note The two `GetMatrixValueType` overloads have different
513  * parameter declarations (`bool` and `...`), to avoid that both
514  * functions are equally good candidates during overload resolution,
515  * in case `element_type` and `ValueType` are both nested types of
516  * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
517  */
518  template <typename QMatrix = TMatrix>
519  auto
520  GetMatrixValueType(bool) const -> typename QMatrix::element_type
521  {
522  return QMatrix::element_type();
523  }
524  template <typename QMatrix = TMatrix>
525  auto
526  GetMatrixValueType(...) const -> typename QMatrix::ValueType
527  {
528  return QMatrix::ValueType();
529  }
530 
531  /* Wrapper that call the right implementation for the type of matrix. */
532  unsigned int
534  TVector & EigenValues,
535  TEigenMatrix & EigenVectors) const
536  {
537  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
538  }
539 
540  /* Implementation detail using EigenLib that performs a copy of the input matrix.
541  *
542  * @param (long) implementation detail argument making this implementation less favourable
543  * to be chosen if alternatives are available.
544  *
545  * @return an unsigned int with no information value (no error code in EigenLib) */
546  template <typename QMatrix>
547  auto
549  TVector & EigenValues,
550  TEigenMatrix & EigenVectors,
551  long) const -> decltype(static_cast<unsigned int>(1))
552  {
553  using ValueType = decltype(GetMatrixValueType(true));
554  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
555  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
556  for (unsigned int row = 0; row < m_Dimension; ++row)
557  {
558  for (unsigned int col = 0; col < m_Dimension; ++col)
559  {
560  inputMatrix(row, col) = A(row, col);
561  }
562  }
563  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
564  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
565  const auto & eigenValues = solver.eigenvalues();
566  /* Column k of the returned matrix is an eigenvector corresponding to
567  * eigenvalue number $ k $ as returned by eigenvalues().
568  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
569  const auto & eigenVectors = solver.eigenvectors();
570 
571  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
572  {
573  auto copyEigenValues = eigenValues;
574  auto copyEigenVectors = eigenVectors;
575  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
576  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
577 
578  for (unsigned int row = 0; row < m_Dimension; ++row)
579  {
580  EigenValues[row] = copyEigenValues[row];
581  for (unsigned int col = 0; col < m_Dimension; ++col)
582  {
583  EigenVectors[row][col] = copyEigenVectors(col, row);
584  }
585  }
586  }
587  else
588  {
589  for (unsigned int row = 0; row < m_Dimension; ++row)
590  {
591  EigenValues[row] = eigenValues[row];
592  for (unsigned int col = 0; col < m_Dimension; ++col)
593  {
594  EigenVectors[row][col] = eigenVectors(col, row);
595  }
596  }
597  }
598  // No error code
599  return 1;
600  }
601 
602 
603  /* Implementation detail using EigenLib that do not peform a copy.
604  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
605  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
606  * should be included.
607  *
608  * @param (bool) implementation detail argument making this implementation the most favourable
609  * to be chosen from all the alternative implementations.
610  *
611  * @return an unsigned int with no information value (no error code in EigenLib) */
612  template <typename QMatrix>
613  auto
615  TVector & EigenValues,
616  TEigenMatrix & EigenVectors,
617  bool) const
618  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
619  {
620  auto pointerToData = GetPointerToMatrixData(A);
621  using PointerType = decltype(pointerToData);
622  using ValueTypeCV = std::remove_pointer_t<PointerType>;
623  using ValueType = std::remove_cv_t<ValueTypeCV>;
624  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
625  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
626  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
627  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
628  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
629  const auto & eigenValues = solver.eigenvalues();
630  /* Column k of the returned matrix is an eigenvector corresponding to
631  * eigenvalue number $ k $ as returned by eigenvalues().
632  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
633  const auto & eigenVectors = solver.eigenvectors();
634  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
635  {
636  auto copyEigenValues = eigenValues;
637  auto copyEigenVectors = eigenVectors;
638  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
639  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
640  for (unsigned int row = 0; row < m_Dimension; ++row)
641  {
642  EigenValues[row] = copyEigenValues[row];
643  for (unsigned int col = 0; col < m_Dimension; ++col)
644  {
645  EigenVectors[row][col] = copyEigenVectors(col, row);
646  }
647  }
648  }
649  else
650  {
651  for (unsigned int row = 0; row < m_Dimension; ++row)
652  {
653  EigenValues[row] = eigenValues[row];
654  for (unsigned int col = 0; col < m_Dimension; ++col)
655  {
656  EigenVectors[row][col] = eigenVectors(col, row);
657  }
658  }
659  }
660  // No error code
661  return 1;
662  }
663 
664  /* Wrapper that call the right implementation for the type of matrix. */
665  unsigned int
666  ComputeEigenValuesWithEigenLibrary(const TMatrix & A, TVector & EigenValues) const
667  {
668  return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
669  }
670 
671  /* Implementation detail using EigenLib that performs a copy of the input matrix.
672  *
673  * @param (long) implementation detail argument making this implementation less favourable
674  * to be chosen if alternatives are available.
675  *
676  * @return an unsigned int with no information value (no error code in EigenLib) */
677  template <typename QMatrix>
678  auto
679  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const
680  -> decltype(static_cast<unsigned int>(1))
681  {
682  using ValueType = decltype(GetMatrixValueType(true));
683  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
684  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
685  for (unsigned int row = 0; row < m_Dimension; ++row)
686  {
687  for (unsigned int col = 0; col < m_Dimension; ++col)
688  {
689  inputMatrix(row, col) = A(row, col);
690  }
691  }
692  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
693  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
694  auto eigenValues = solver.eigenvalues();
695  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
696  {
697  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
698  }
699  for (unsigned int i = 0; i < m_Dimension; ++i)
700  {
701  EigenValues[i] = eigenValues[i];
702  }
703 
704  // No error code
705  return 1;
706  }
707 
708  /* Implementation detail using EigenLib that do not peform a copy.
709  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
710  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
711  * should be included.
712  *
713  * @param (bool) implementation detail argument making this implementation the most favourable
714  * to be chosen from all the alternative implementations.
715  *
716  * @return an unsigned int with no information value (no error code in EigenLib) */
717  template <typename QMatrix>
718  auto
719  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
720  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
721  {
722  auto pointerToData = GetPointerToMatrixData(A);
723  using PointerType = decltype(pointerToData);
724  using ValueTypeCV = std::remove_pointer_t<PointerType>;
725  using ValueType = std::remove_cv_t<ValueTypeCV>;
726  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
727  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
728  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
729  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
730  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
731  auto eigenValues = solver.eigenvalues();
732  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
733  {
734  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
735  }
736  for (unsigned int i = 0; i < m_Dimension; ++i)
737  {
738  EigenValues[i] = eigenValues[i];
739  }
740  // No error code
741  return 1;
742  }
743 };
744 
745 template <typename TMatrix, typename TVector, typename TEigenMatrix>
746 std::ostream &
748 {
749  os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
750  os << " Dimension : " << s.GetDimension() << std::endl;
751  os << " Order : " << s.GetOrder() << std::endl;
752  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
753  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
754  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
755  return os;
756 }
757 
758 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
759 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysisFixedDimension
760 {
761 public:
762 #if !defined(ITK_LEGACY_REMOVE)
763 
764  using EigenValueOrderType = EigenValueOrderEnum;
765 #endif
766 #if !defined(ITK_LEGACY_REMOVE)
767  // We need to expose the enum values at the class level
768  // for backwards compatibility
769  static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
770  static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
771  static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
772 #endif
773 
776 
777  using MatrixType = TMatrix;
778  using EigenMatrixType = TEigenMatrix;
779  using VectorType = TVector;
780 
794  unsigned int
795  ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const
796  {
797  return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
798  }
799 
820  unsigned int
821  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const
822  {
823  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
824  }
825 
826  void
827  SetOrderEigenValues(const bool b)
828  {
829  if (b)
830  {
831  m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
832  }
833  else
834  {
835  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
836  }
837  }
838  bool
840  {
841  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
842  }
843  void
845  {
846  if (b)
847  {
848  m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
849  }
850  else
851  {
852  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
853  }
854  }
855  bool
857  {
858  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
859  }
860  constexpr unsigned int
861  GetOrder() const
862  {
863  return VDimension;
864  }
865  constexpr unsigned int
866  GetDimension() const
867  {
868  return VDimension;
869  }
870  constexpr bool
872  {
873  return true;
874  }
875 
876 private:
878 
879  /* Helper to get the matrix value type for EigenLibMatrix typename.
880  *
881  * If the TMatrix is vnl, the type is in element_type.
882  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
883  *
884  * To use this function:
885  * using ValueType = decltype(this->GetMatrixType(true));
886  */
887  template <typename QMatrix = TMatrix>
888  auto
889  GetMatrixValueType(bool) const -> typename QMatrix::element_type
890  {
891  return QMatrix::element_type();
892  }
893  template <typename QMatrix = TMatrix>
894  auto
895  GetMatrixValueType(bool) const -> typename QMatrix::ValueType
896  {
897  return QMatrix::ValueType();
898  }
899 
900  /* Implementation detail using EigenLib that do not peform a copy.
901  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
902  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
903  * should be included.
904  *
905  * @param (bool) implementation detail argument making this implementation the most favourable
906  * to be chosen from all the alternative implementations.
907  *
908  * @return an unsigned int with no information value (no error code in EigenLib) */
909  template <typename QMatrix>
910  auto
912  TVector & EigenValues,
913  TEigenMatrix & EigenVectors,
914  bool) const
915  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
916  {
917  auto pointerToData = GetPointerToMatrixData(A);
918  using PointerType = decltype(pointerToData);
919  using ValueTypeCV = std::remove_pointer_t<PointerType>;
920  using ValueType = std::remove_cv_t<ValueTypeCV>;
921  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
922  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
923  EigenConstMatrixMap inputMatrix(pointerToData);
924  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
925  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
926  const auto & eigenValues = solver.eigenvalues();
927  /* Column k of the returned matrix is an eigenvector corresponding to
928  * eigenvalue number $ k $ as returned by eigenvalues().
929  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
930  const auto & eigenVectors = solver.eigenvectors();
931  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
932  {
933  auto copyEigenValues = eigenValues;
934  auto copyEigenVectors = eigenVectors;
935  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
936  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
937  for (unsigned int row = 0; row < VDimension; ++row)
938  {
939  EigenValues[row] = copyEigenValues[row];
940  for (unsigned int col = 0; col < VDimension; ++col)
941  {
942  EigenVectors[row][col] = copyEigenVectors(col, row);
943  }
944  }
945  }
946  else
947  {
948  for (unsigned int row = 0; row < VDimension; ++row)
949  {
950  EigenValues[row] = eigenValues[row];
951  for (unsigned int col = 0; col < VDimension; ++col)
952  {
953  EigenVectors[row][col] = eigenVectors(col, row);
954  }
955  }
956  }
957  // No error code
958  return 1;
959  }
960 
961  /* Implementation detail using EigenLib that performs a copy of the input matrix.
962  *
963  * @param (long) implementation detail argument making this implementation less favourable
964  * to be chosen if alternatives are available.
965  *
966  * @return an unsigned int with no information value (no error code in EigenLib) */
967  template <typename QMatrix>
968  auto
970  TVector & EigenValues,
971  TEigenMatrix & EigenVectors,
972  long) const -> decltype(static_cast<unsigned int>(1))
973  {
974  using ValueType = decltype(GetMatrixValueType(true));
975  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
976  EigenLibMatrixType inputMatrix;
977  for (unsigned int row = 0; row < VDimension; ++row)
978  {
979  for (unsigned int col = 0; col < VDimension; ++col)
980  {
981  inputMatrix(row, col) = A(row, col);
982  }
983  }
984  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
985  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
986  const auto & eigenValues = solver.eigenvalues();
987  /* Column k of the returned matrix is an eigenvector corresponding to
988  * eigenvalue number $ k $ as returned by eigenvalues().
989  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
990  const auto & eigenVectors = solver.eigenvectors();
991 
992  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
993  {
994  auto copyEigenValues = eigenValues;
995  auto copyEigenVectors = eigenVectors;
996  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
997  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
998 
999  for (unsigned int row = 0; row < VDimension; ++row)
1000  {
1001  EigenValues[row] = copyEigenValues[row];
1002  for (unsigned int col = 0; col < VDimension; ++col)
1003  {
1004  EigenVectors[row][col] = copyEigenVectors(col, row);
1005  }
1006  }
1007  }
1008  else
1009  {
1010  for (unsigned int row = 0; row < VDimension; ++row)
1011  {
1012  EigenValues[row] = eigenValues[row];
1013  for (unsigned int col = 0; col < VDimension; ++col)
1014  {
1015  EigenVectors[row][col] = eigenVectors(col, row);
1016  }
1017  }
1018  }
1019  // No error code
1020  return 1;
1021  }
1022 
1023  /* Implementation detail using EigenLib that performs a copy of the input matrix.
1024  *
1025  * @param (long) implementation detail argument making this implementation less favourable
1026  * to be chosen if alternatives are available.
1027  *
1028  * @return an unsigned int with no information value (no error code in EigenLib) */
1029  template <typename QMatrix>
1030  auto
1031  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const
1032  -> decltype(static_cast<unsigned int>(1))
1033  {
1034  using ValueType = decltype(GetMatrixValueType(true));
1035  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1036  EigenLibMatrixType inputMatrix;
1037  for (unsigned int row = 0; row < VDimension; ++row)
1038  {
1039  for (unsigned int col = 0; col < VDimension; ++col)
1040  {
1041  inputMatrix(row, col) = A(row, col);
1042  }
1043  }
1044  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1045  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1046  auto eigenValues = solver.eigenvalues();
1047  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1048  {
1049  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1050  }
1051  for (unsigned int i = 0; i < VDimension; ++i)
1052  {
1053  EigenValues[i] = eigenValues[i];
1054  }
1055 
1056  // No error code
1057  return 1;
1058  }
1059 
1060  /* Implementation detail using EigenLib that do not peform a copy.
1061  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
1062  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
1063  * should be included.
1064  *
1065  * @param (bool) implementation detail argument making this implementation the most favourable
1066  * to be chosen from all the alternative implementations.
1067  *
1068  * @return an unsigned int with no information value (no error code in EigenLib) */
1069  template <typename QMatrix>
1070  auto
1071  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
1072  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
1073  {
1074  auto pointerToData = GetPointerToMatrixData(A);
1075  using PointerType = decltype(pointerToData);
1076  using ValueTypeCV = std::remove_pointer_t<PointerType>;
1077  using ValueType = std::remove_cv_t<ValueTypeCV>;
1078  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1079  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1080  EigenConstMatrixMap inputMatrix(pointerToData);
1081  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1082  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1083  auto eigenValues = solver.eigenvalues();
1084  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1085  {
1086  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1087  }
1088  for (unsigned int i = 0; i < VDimension; ++i)
1089  {
1090  EigenValues[i] = eigenValues[i];
1091  }
1092  // No error code
1093  return 1;
1094  }
1095 };
1096 
1097 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix>
1098 std::ostream &
1099 operator<<(std::ostream & os,
1101 {
1102  os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1103  os << " Dimension : " << s.GetDimension() << std::endl;
1104  os << " Order : " << s.GetOrder() << std::endl;
1105  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
1106  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
1107  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
1108  return os;
1109 }
1110 } // end namespace itk
1111 
1112 #ifndef ITK_MANUAL_INSTANTIATION
1113 # include "itkSymmetricEigenAnalysis.hxx"
1114 #endif
1115 
1116 #endif
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::VectorType
TOutput VectorType
Definition: itkSymmetricEigenAnalysis.h:779
EigenValueOrder
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectors
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:821
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:216
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:548
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:614
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:719
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:911
itk::SymmetricEigenAnalysis::SymmetricEigenAnalysis
SymmetricEigenAnalysis(const unsigned int dimension)
Definition: itkSymmetricEigenAnalysis.h:214
itk::SymmetricEigenAnalysisFixedDimension::SetOrderEigenValues
void SetOrderEigenValues(const bool b)
Definition: itkSymmetricEigenAnalysis.h:827
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:895
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
Definition: itkSymmetricEigenAnalysis.h:520
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibrary
unsigned int ComputeEigenValuesAndVectorsWithEigenLibrary(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:533
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:1031
itk::Int2EigenValueOrderEnum
EigenValueOrderEnum Int2EigenValueOrderEnum(const uint8_t value)
Definition: itkSymmetricEigenAnalysis.h:144
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(...) const -> typename QMatrix::ValueType
Definition: itkSymmetricEigenAnalysis.h:526
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:844
itk::SymmetricEigenAnalysisFixedDimension::GetOrderEigenMagnitudes
bool GetOrderEigenMagnitudes() const
Definition: itkSymmetricEigenAnalysis.h:856
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:778
itk::uint8_t
std::uint8_t uint8_t
Definition: itkIntTypes.h:29
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:889
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:1071
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:969
itk::SymmetricEigenAnalysisFixedDimension::GetDimension
constexpr unsigned int GetDimension() const
Definition: itkSymmetricEigenAnalysis.h:866
itk::SymmetricEigenAnalysisFixedDimension
Definition: itkSymmetricEigenAnalysis.h:759
itk::SymmetricEigenAnalysisFixedDimension::GetUseEigenLibrary
constexpr bool GetUseEigenLibrary() const
Definition: itkSymmetricEigenAnalysis.h:871
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:839
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:679
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:666
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::MatrixType
TInput MatrixType
Definition: itkSymmetricEigenAnalysis.h:777
itk::SymmetricEigenAnalysisFixedDimension::GetOrder
constexpr unsigned int GetOrder() const
Definition: itkSymmetricEigenAnalysis.h:861
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:795
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