ITK  6.0.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  * https://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef 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 VRows, unsigned int VColumns>
37 const TValueType *
38 GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, VRows, VColumns> & 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 VRows, unsigned int VColumns>
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);
79  std::sort(std::begin(indicesSortPermutations),
80  std::end(indicesSortPermutations),
81  [&eigenValues](unsigned int a, unsigned int b) {
82  return itk::Math::abs(eigenValues[a]) < itk::Math::abs(eigenValues[b]);
83  });
84  auto tmpCopy = eigenValues;
85  for (unsigned int i = 0; i < numberOfElements; ++i)
86  {
87  eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
88  }
89  return indicesSortPermutations;
90 }
91 
100 template <typename QMatrix>
101 void
102 permuteColumnsWithSortIndices(QMatrix & eigenVectors, const std::vector<int> & indicesSortPermutations)
103 {
104  using EigenLibPermutationMatrix = Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
105  auto numberOfElements = indicesSortPermutations.size();
106  // Creates a NxN permutation matrix copying our permutation to the matrix indices.
107  // Which holds the 1D array representation of a permutation.
108  EigenLibPermutationMatrix perm(numberOfElements);
109  perm.setIdentity();
110  std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(), perm.indices().data());
111  // Apply it
112  eigenVectors = eigenVectors * perm;
113 }
114 } // end namespace detail
122 {
123 public:
131  enum class EigenValueOrder : uint8_t
132  {
133  OrderByValue = 1,
134  OrderByMagnitude = 2,
135  DoNotOrder = 3
136  };
137 };
138 // Define how to print enumeration
139 extern ITKCommon_EXPORT std::ostream &
140  operator<<(std::ostream & out, const SymmetricEigenAnalysisEnums::EigenValueOrder value);
141 
143 
144 inline EigenValueOrderEnum
145 Int2EigenValueOrderEnum(const uint8_t value)
146 {
147  switch (value)
148  {
149  case 1:
151  case 2:
153  case 3:
155  default:
156  break;
157  }
158  itkGenericExceptionMacro("Error: Invalid value for conversion.");
159 }
160 
161 #if !defined(ITK_LEGACY_REMOVE)
162 
163 static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
164 static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
165 static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
166 #endif
167 
203 template <typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
204 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
205 {
206 public:
208 #if !defined(ITK_LEGACY_REMOVE)
209 
210  using EigenValueOrderType = EigenValueOrderEnum;
211 #endif
212 
213  SymmetricEigenAnalysis() = default;
214 
215  SymmetricEigenAnalysis(const unsigned int dimension)
216  : m_Dimension(dimension)
217  , m_Order(dimension)
218  {}
219 
220  ~SymmetricEigenAnalysis() = default;
221 
222  using MatrixType = TMatrix;
223  using EigenMatrixType = TEigenMatrix;
224  using VectorType = TVector;
225 
239  unsigned int
240  ComputeEigenValues(const TMatrix & A, TVector & D) const;
241 
262  unsigned int
263  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
264 
265 
267  void
268  SetOrder(const unsigned int n)
269  {
270  m_Order = n;
271  }
272 
276  unsigned int
277  GetOrder() const
278  {
279  return m_Order;
280  }
281 
285  void
286  SetOrderEigenValues(const bool b)
287  {
288  if (b)
289  {
290  m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
291  }
292  else
293  {
294  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
295  }
296  }
299  bool
301  {
302  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
303  }
304 
308  void
310  {
311  if (b)
312  {
313  m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
314  }
315  else
316  {
317  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
318  }
319  }
322  bool
324  {
325  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
326  }
327 
330  void
331  SetDimension(const unsigned int n)
332  {
333  m_Dimension = n;
334  if (m_Order == 0)
335  {
336  m_Order = m_Dimension;
337  }
338  }
343  unsigned int
344  GetDimension() const
345  {
346  return m_Dimension;
347  }
348 
350  void
351  SetUseEigenLibrary(const bool input)
352  {
353  m_UseEigenLibrary = input;
354  }
355  void
357  {
358  m_UseEigenLibrary = true;
359  }
360  void
362  {
363  m_UseEigenLibrary = false;
364  }
365  bool
367  {
368  return m_UseEigenLibrary;
369  }
372 private:
373  bool m_UseEigenLibrary{ false };
374  unsigned int m_Dimension{ 0 };
375  unsigned int m_Order{ 0 };
377 
397  void
398  ReduceToTridiagonalMatrix(double * a, double * d, double * e, double * e2) const;
399 
421  void
422  ReduceToTridiagonalMatrixAndGetTransformation(const double * a, double * d, double * e, double * z) const;
423 
453  unsigned int
454  ComputeEigenValuesUsingQL(double * d, double * e) const;
455 
493  unsigned int
494  ComputeEigenValuesAndVectorsUsingQL(double * d, double * e, double * z) const;
495 
496  /* Legacy algorithms using thread-safe netlib.
497  * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
498  */
499  unsigned int
500  ComputeEigenValuesLegacy(const TMatrix & A, TVector & D) const;
501 
502  unsigned int
503  ComputeEigenValuesAndVectorsLegacy(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
504 
505  /* Helper to get the matrix value type for EigenLibMatrix typename.
506  *
507  * If the TMatrix is vnl, the type is in element_type.
508  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
509  *
510  * To use this function:
511  * using ValueType = decltype(this->GetMatrixType(true));
512  *
513  * \note The two `GetMatrixValueType` overloads have different
514  * parameter declarations (`bool` and `...`), to avoid that both
515  * functions are equally good candidates during overload resolution,
516  * in case `element_type` and `ValueType` are both nested types of
517  * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
518  */
519  template <typename QMatrix = TMatrix>
520  auto
521  GetMatrixValueType(bool) const -> typename QMatrix::element_type
522  {
523  return QMatrix::element_type();
524  }
525  template <typename QMatrix = TMatrix>
526  auto
527  GetMatrixValueType(...) const -> typename QMatrix::ValueType
528  {
529  return QMatrix::ValueType();
530  }
531 
532  /* Wrapper that call the right implementation for the type of matrix. */
533  unsigned int
535  TVector & EigenValues,
536  TEigenMatrix & EigenVectors) const
537  {
538  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
539  }
540 
541  /* Implementation detail using EigenLib that performs a copy of the input matrix.
542  *
543  * @param (long) implementation detail argument making this implementation less favourable
544  * to be chosen if alternatives are available.
545  *
546  * @return an unsigned int with no information value (no error code in EigenLib) */
547  template <typename QMatrix>
548  auto
550  TVector & EigenValues,
551  TEigenMatrix & EigenVectors,
552  long) const -> decltype(1U)
553  {
554  using ValueType = decltype(GetMatrixValueType(true));
555  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
556  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
557  for (unsigned int row = 0; row < m_Dimension; ++row)
558  {
559  for (unsigned int col = 0; col < m_Dimension; ++col)
560  {
561  inputMatrix(row, col) = A(row, col);
562  }
563  }
564  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
565  const EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
566  const auto & eigenValues = solver.eigenvalues();
567  /* Column k of the returned matrix is an eigenvector corresponding to
568  * eigenvalue number $ k $ as returned by eigenvalues().
569  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
570  const auto & eigenVectors = solver.eigenvectors();
571 
572  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
573  {
574  auto copyEigenValues = eigenValues;
575  auto copyEigenVectors = eigenVectors;
576  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
577  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
578 
579  for (unsigned int row = 0; row < m_Dimension; ++row)
580  {
581  EigenValues[row] = copyEigenValues[row];
582  for (unsigned int col = 0; col < m_Dimension; ++col)
583  {
584  EigenVectors[row][col] = copyEigenVectors(col, row);
585  }
586  }
587  }
588  else
589  {
590  for (unsigned int row = 0; row < m_Dimension; ++row)
591  {
592  EigenValues[row] = eigenValues[row];
593  for (unsigned int col = 0; col < m_Dimension; ++col)
594  {
595  EigenVectors[row][col] = eigenVectors(col, row);
596  }
597  }
598  }
599  // No error code
600  return 1;
601  }
602 
603 
604  /* Implementation detail using EigenLib that do not perform a copy.
605  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
606  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
607  * should be included.
608  *
609  * @param (bool) implementation detail argument making this implementation the most favourable
610  * to be chosen from all the alternative implementations.
611  *
612  * @return an unsigned int with no information value (no error code in EigenLib) */
613  template <typename QMatrix>
614  auto
616  TVector & EigenValues,
617  TEigenMatrix & EigenVectors,
618  bool) const -> decltype(GetPointerToMatrixData(A), 1U)
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 -> decltype(1U)
680  {
681  using ValueType = decltype(GetMatrixValueType(true));
682  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
683  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
684  for (unsigned int row = 0; row < m_Dimension; ++row)
685  {
686  for (unsigned int col = 0; col < m_Dimension; ++col)
687  {
688  inputMatrix(row, col) = A(row, col);
689  }
690  }
691  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
692  const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
693  auto eigenValues = solver.eigenvalues();
694  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
695  {
696  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
697  }
698  for (unsigned int i = 0; i < m_Dimension; ++i)
699  {
700  EigenValues[i] = eigenValues[i];
701  }
702 
703  // No error code
704  return 1;
705  }
706 
707  /* Implementation detail using EigenLib that do not perform a copy.
708  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
709  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
710  * should be included.
711  *
712  * @param (bool) implementation detail argument making this implementation the most favourable
713  * to be chosen from all the alternative implementations.
714  *
715  * @return an unsigned int with no information value (no error code in EigenLib) */
716  template <typename QMatrix>
717  auto
718  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
719  -> decltype(GetPointerToMatrixData(A), 1U)
720  {
721  auto pointerToData = GetPointerToMatrixData(A);
722  using PointerType = decltype(pointerToData);
723  using ValueTypeCV = std::remove_pointer_t<PointerType>;
724  using ValueType = std::remove_cv_t<ValueTypeCV>;
725  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
726  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
727  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
728  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
729  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
730  auto eigenValues = solver.eigenvalues();
731  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
732  {
733  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
734  }
735  for (unsigned int i = 0; i < m_Dimension; ++i)
736  {
737  EigenValues[i] = eigenValues[i];
738  }
739  // No error code
740  return 1;
741  }
742 };
743 
744 template <typename TMatrix, typename TVector, typename TEigenMatrix>
745 std::ostream &
747 {
748  os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
749  os << " Dimension : " << s.GetDimension() << std::endl;
750  os << " Order : " << s.GetOrder() << std::endl;
751  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
752  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
753  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
754  return os;
755 }
756 
757 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
758 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysisFixedDimension
759 {
760 public:
761 #if !defined(ITK_LEGACY_REMOVE)
762 
763  using EigenValueOrderType = EigenValueOrderEnum;
764 #endif
765 #if !defined(ITK_LEGACY_REMOVE)
766  // We need to expose the enum values at the class level
767  // for backwards compatibility
768  static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
769  static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
770  static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
771 #endif
772 
773  using MatrixType = TMatrix;
774  using EigenMatrixType = TEigenMatrix;
775  using VectorType = TVector;
776 
790  unsigned int
791  ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const
792  {
793  return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
794  }
795 
816  unsigned int
817  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const
818  {
819  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
820  }
821 
822  void
823  SetOrderEigenValues(const bool b)
824  {
825  if (b)
826  {
827  m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
828  }
829  else
830  {
831  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
832  }
833  }
834  bool
836  {
837  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
838  }
839  void
841  {
842  if (b)
843  {
844  m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
845  }
846  else
847  {
848  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
849  }
850  }
851  bool
853  {
854  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
855  }
856  constexpr unsigned int
857  GetOrder() const
858  {
859  return VDimension;
860  }
861  constexpr unsigned int
862  GetDimension() const
863  {
864  return VDimension;
865  }
866  constexpr bool
868  {
869  return true;
870  }
871 
872 private:
874 
875  /* Helper to get the matrix value type for EigenLibMatrix typename.
876  *
877  * If the TMatrix is vnl, the type is in element_type.
878  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
879  *
880  * To use this function:
881  * using ValueType = decltype(this->GetMatrixType(true));
882  */
883  template <typename QMatrix = TMatrix>
884  auto
885  GetMatrixValueType(bool) const -> typename QMatrix::element_type
886  {
887  return QMatrix::element_type();
888  }
889  template <typename QMatrix = TMatrix>
890  auto
891  GetMatrixValueType(bool) const -> typename QMatrix::ValueType
892  {
893  return QMatrix::ValueType();
894  }
895 
896  /* Implementation detail using EigenLib that do not perform a copy.
897  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
898  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
899  * should be included.
900  *
901  * @param (bool) implementation detail argument making this implementation the most favourable
902  * to be chosen from all the alternative implementations.
903  *
904  * @return an unsigned int with no information value (no error code in EigenLib) */
905  template <typename QMatrix>
906  auto
908  TVector & EigenValues,
909  TEigenMatrix & EigenVectors,
910  bool) const -> decltype(GetPointerToMatrixData(A), 1U)
911  {
912  auto pointerToData = GetPointerToMatrixData(A);
913  using PointerType = decltype(pointerToData);
914  using ValueTypeCV = std::remove_pointer_t<PointerType>;
915  using ValueType = std::remove_cv_t<ValueTypeCV>;
916  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
917  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
918  EigenConstMatrixMap inputMatrix(pointerToData);
919  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
920  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
921  const auto & eigenValues = solver.eigenvalues();
922  /* Column k of the returned matrix is an eigenvector corresponding to
923  * eigenvalue number $ k $ as returned by eigenvalues().
924  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
925  const auto & eigenVectors = solver.eigenvectors();
926  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
927  {
928  auto copyEigenValues = eigenValues;
929  auto copyEigenVectors = eigenVectors;
930  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
931  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
932  for (unsigned int row = 0; row < VDimension; ++row)
933  {
934  EigenValues[row] = copyEigenValues[row];
935  for (unsigned int col = 0; col < VDimension; ++col)
936  {
937  EigenVectors[row][col] = copyEigenVectors(col, row);
938  }
939  }
940  }
941  else
942  {
943  for (unsigned int row = 0; row < VDimension; ++row)
944  {
945  EigenValues[row] = eigenValues[row];
946  for (unsigned int col = 0; col < VDimension; ++col)
947  {
948  EigenVectors[row][col] = eigenVectors(col, row);
949  }
950  }
951  }
952  // No error code
953  return 1;
954  }
955 
956  /* Implementation detail using EigenLib that performs a copy of the input matrix.
957  *
958  * @param (long) implementation detail argument making this implementation less favourable
959  * to be chosen if alternatives are available.
960  *
961  * @return an unsigned int with no information value (no error code in EigenLib) */
962  template <typename QMatrix>
963  auto
965  TVector & EigenValues,
966  TEigenMatrix & EigenVectors,
967  long) const -> decltype(1U)
968  {
969  using ValueType = decltype(GetMatrixValueType(true));
970  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
971  EigenLibMatrixType inputMatrix;
972  for (unsigned int row = 0; row < VDimension; ++row)
973  {
974  for (unsigned int col = 0; col < VDimension; ++col)
975  {
976  inputMatrix(row, col) = A(row, col);
977  }
978  }
979  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
980  const EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
981  const auto & eigenValues = solver.eigenvalues();
982  /* Column k of the returned matrix is an eigenvector corresponding to
983  * eigenvalue number $ k $ as returned by eigenvalues().
984  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
985  const auto & eigenVectors = solver.eigenvectors();
986 
987  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
988  {
989  auto copyEigenValues = eigenValues;
990  auto copyEigenVectors = eigenVectors;
991  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
992  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
993 
994  for (unsigned int row = 0; row < VDimension; ++row)
995  {
996  EigenValues[row] = copyEigenValues[row];
997  for (unsigned int col = 0; col < VDimension; ++col)
998  {
999  EigenVectors[row][col] = copyEigenVectors(col, row);
1000  }
1001  }
1002  }
1003  else
1004  {
1005  for (unsigned int row = 0; row < VDimension; ++row)
1006  {
1007  EigenValues[row] = eigenValues[row];
1008  for (unsigned int col = 0; col < VDimension; ++col)
1009  {
1010  EigenVectors[row][col] = eigenVectors(col, row);
1011  }
1012  }
1013  }
1014  // No error code
1015  return 1;
1016  }
1017 
1018  /* Implementation detail using EigenLib that performs a copy of the input matrix.
1019  *
1020  * @param (long) implementation detail argument making this implementation less favourable
1021  * to be chosen if alternatives are available.
1022  *
1023  * @return an unsigned int with no information value (no error code in EigenLib) */
1024  template <typename QMatrix>
1025  auto
1026  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const -> decltype(1U)
1027  {
1028  using ValueType = decltype(GetMatrixValueType(true));
1029  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1030  EigenLibMatrixType inputMatrix;
1031  for (unsigned int row = 0; row < VDimension; ++row)
1032  {
1033  for (unsigned int col = 0; col < VDimension; ++col)
1034  {
1035  inputMatrix(row, col) = A(row, col);
1036  }
1037  }
1038  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1039  const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1040  auto eigenValues = solver.eigenvalues();
1041  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1042  {
1043  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1044  }
1045  for (unsigned int i = 0; i < VDimension; ++i)
1046  {
1047  EigenValues[i] = eigenValues[i];
1048  }
1049 
1050  // No error code
1051  return 1;
1052  }
1053 
1054  /* Implementation detail using EigenLib that do not perform a copy.
1055  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
1056  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
1057  * should be included.
1058  *
1059  * @param (bool) implementation detail argument making this implementation the most favourable
1060  * to be chosen from all the alternative implementations.
1061  *
1062  * @return an unsigned int with no information value (no error code in EigenLib) */
1063  template <typename QMatrix>
1064  auto
1065  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
1066  -> decltype(GetPointerToMatrixData(A), 1U)
1067  {
1068  auto pointerToData = GetPointerToMatrixData(A);
1069  using PointerType = decltype(pointerToData);
1070  using ValueTypeCV = std::remove_pointer_t<PointerType>;
1071  using ValueType = std::remove_cv_t<ValueTypeCV>;
1072  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1073  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1074  EigenConstMatrixMap inputMatrix(pointerToData);
1075  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1076  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1077  auto eigenValues = solver.eigenvalues();
1078  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1079  {
1080  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1081  }
1082  for (unsigned int i = 0; i < VDimension; ++i)
1083  {
1084  EigenValues[i] = eigenValues[i];
1085  }
1086  // No error code
1087  return 1;
1088  }
1089 };
1090 
1091 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix>
1092 std::ostream &
1093 operator<<(std::ostream & os,
1095 {
1096  os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1097  os << " Dimension : " << s.GetDimension() << std::endl;
1098  os << " Order : " << s.GetOrder() << std::endl;
1099  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
1100  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
1101  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
1102  return os;
1103 }
1104 } // end namespace itk
1105 
1106 #ifndef ITK_MANUAL_INSTANTIATION
1107 # include "itkSymmetricEigenAnalysis.hxx"
1108 #endif
1109 
1110 #endif
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::VectorType
TOutput VectorType
Definition: itkSymmetricEigenAnalysis.h:775
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectors
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:817
itk::SymmetricEigenAnalysis::SetDimension
void SetDimension(const unsigned int n)
Definition: itkSymmetricEigenAnalysis.h:331
itk::SymmetricEigenAnalysis::GetOrder
unsigned int GetOrder() const
Definition: itkSymmetricEigenAnalysis.h:277
itkMatrix.h
itk::SymmetricEigenAnalysis::GetUseEigenLibrary
bool GetUseEigenLibrary() const
Definition: itkSymmetricEigenAnalysis.h:366
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
Definition: itkSymmetricEigenAnalysis.h:907
itk::SymmetricEigenAnalysisEnums::EigenValueOrder::OrderByMagnitude
itk::SymmetricEigenAnalysis::SetUseEigenLibraryOn
void SetUseEigenLibraryOn()
Definition: itkSymmetricEigenAnalysis.h:356
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
Definition: itkSymmetricEigenAnalysis.h:1065
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::VectorType
TOutputImage::PixelType VectorType
Definition: itkSymmetricEigenAnalysis.h:224
itk::SymmetricEigenAnalysis::SymmetricEigenAnalysis
SymmetricEigenAnalysis(const unsigned int dimension)
Definition: itkSymmetricEigenAnalysis.h:215
itk::SymmetricEigenAnalysisFixedDimension::SetOrderEigenValues
void SetOrderEigenValues(const bool b)
Definition: itkSymmetricEigenAnalysis.h:823
itk::SymmetricEigenAnalysisEnums::EigenValueOrder::OrderByValue
itk::detail::permuteColumnsWithSortIndices
void permuteColumnsWithSortIndices(QMatrix &eigenVectors, const std::vector< int > &indicesSortPermutations)
Definition: itkSymmetricEigenAnalysis.h:102
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(1U)
Definition: itkSymmetricEigenAnalysis.h:1026
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(1U)
Definition: itkSymmetricEigenAnalysis.h:964
itk::operator<<
ITKCommon_EXPORT std::ostream & operator<<(std::ostream &out, typename AnatomicalOrientation::CoordinateEnum value)
itk::SymmetricEigenAnalysis::GetOrderEigenValues
bool GetOrderEigenValues() const
Definition: itkSymmetricEigenAnalysis.h:300
itk::SymmetricEigenAnalysisFixedDimension::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::ValueType
Definition: itkSymmetricEigenAnalysis.h:891
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
Definition: itkSymmetricEigenAnalysis.h:521
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibrary
unsigned int ComputeEigenValuesAndVectorsWithEigenLibrary(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:534
itk::Math::abs
bool abs(bool x)
Definition: itkMath.h:839
itk::EigenValueOrderEnum
SymmetricEigenAnalysisEnums::EigenValueOrder EigenValueOrderEnum
Definition: itkSymmetricEigenAnalysis.h:142
itk::SymmetricEigenAnalysis::SetOrder
void SetOrder(const unsigned int n)
Definition: itkSymmetricEigenAnalysis.h:268
itk::SymmetricEigenAnalysisEnums::EigenValueOrder
EigenValueOrder
Definition: itkSymmetricEigenAnalysis.h:131
itkMacro.h
itk::SymmetricEigenAnalysisEnums::EigenValueOrder::DoNotOrder
itk::SymmetricEigenAnalysis::GetDimension
unsigned int GetDimension() const
Definition: itkSymmetricEigenAnalysis.h:344
itk::Matrix::GetVnlMatrix
InternalMatrixType & GetVnlMatrix()
Definition: itkMatrix.h:214
itk::detail::GetPointerToMatrixData
const TValueType * GetPointerToMatrixData(const vnl_matrix_fixed< TValueType, VRows, VColumns > &inputMatrix)
Definition: itkSymmetricEigenAnalysis.h:38
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(1U)
Definition: itkSymmetricEigenAnalysis.h:679
itk::Int2EigenValueOrderEnum
EigenValueOrderEnum Int2EigenValueOrderEnum(const uint8_t value)
Definition: itkSymmetricEigenAnalysis.h:145
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(...) const -> typename QMatrix::ValueType
Definition: itkSymmetricEigenAnalysis.h:527
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:840
itk::SymmetricEigenAnalysisFixedDimension::GetOrderEigenMagnitudes
bool GetOrderEigenMagnitudes() const
Definition: itkSymmetricEigenAnalysis.h:852
itk::SymmetricEigenAnalysis::SetOrderEigenMagnitudes
void SetOrderEigenMagnitudes(const bool b)
Definition: itkSymmetricEigenAnalysis.h:309
itk::SymmetricEigenAnalysis::SetUseEigenLibraryOff
void SetUseEigenLibraryOff()
Definition: itkSymmetricEigenAnalysis.h:361
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::EigenMatrixType
TInput EigenMatrixType
Definition: itkSymmetricEigenAnalysis.h:774
itk::Matrix
A templated class holding a M x N size Matrix.
Definition: itkMatrix.h:52
itk::SymmetricEigenAnalysisFixedDimension::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
Definition: itkSymmetricEigenAnalysis.h:885
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::EigenMatrixType
TInputImage::PixelType EigenMatrixType
Definition: itkSymmetricEigenAnalysis.h:223
itk::SymmetricEigenAnalysis::SetUseEigenLibrary
void SetUseEigenLibrary(const bool input)
Definition: itkSymmetricEigenAnalysis.h:351
itk::SymmetricEigenAnalysisFixedDimension::GetDimension
constexpr unsigned int GetDimension() const
Definition: itkSymmetricEigenAnalysis.h:862
itk::SymmetricEigenAnalysisFixedDimension
Definition: itkSymmetricEigenAnalysis.h:758
itk::SymmetricEigenAnalysisFixedDimension::GetUseEigenLibrary
constexpr bool GetUseEigenLibrary() const
Definition: itkSymmetricEigenAnalysis.h:867
itk::SymmetricEigenAnalysis::SetOrderEigenValues
void SetOrderEigenValues(const bool b)
Definition: itkSymmetricEigenAnalysis.h:286
itk
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
Definition: itkAnatomicalOrientation.h:29
itk::SymmetricEigenAnalysisFixedDimension::GetOrderEigenValues
bool GetOrderEigenValues() const
Definition: itkSymmetricEigenAnalysis.h:835
itk::Math::e
static constexpr double e
Definition: itkMath.h:56
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:773
itk::SymmetricEigenAnalysisFixedDimension::GetOrder
constexpr unsigned int GetOrder() const
Definition: itkSymmetricEigenAnalysis.h:857
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(1U)
Definition: itkSymmetricEigenAnalysis.h:549
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::MatrixType
TInputImage::PixelType MatrixType
Definition: itkSymmetricEigenAnalysis.h:222
itk::SymmetricEigenAnalysisEnums
This class contains all enum classes used by SymmetricEigenAnalysis class.
Definition: itkSymmetricEigenAnalysis.h:121
itk::detail::GetPointerToMatrixData
const TValueType * GetPointerToMatrixData(const itk::Matrix< TValueType, VRows, VColumns > &inputMatrix)
Definition: itkSymmetricEigenAnalysis.h:51
itk::SymmetricEigenAnalysis::GetOrderEigenMagnitudes
bool GetOrderEigenMagnitudes() const
Definition: itkSymmetricEigenAnalysis.h:323
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValues
unsigned int ComputeEigenValues(const TMatrix &A, TVector &EigenValues) const
Definition: itkSymmetricEigenAnalysis.h:791
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
Definition: itkSymmetricEigenAnalysis.h:615
itk::SymmetricEigenAnalysis
Find Eigen values of a real 2D symmetric matrix. It serves as a thread-safe alternative to the class:...
Definition: itkSymmetricEigenAnalysis.h:204
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
Definition: itkSymmetricEigenAnalysis.h:718