ITK  5.1.0
Insight Toolkit
itkSymmetricEigenAnalysis.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef itkSymmetricEigenAnalysis_h
19 #define itkSymmetricEigenAnalysis_h
20 
21 #include "itkMacro.h"
22 #include "itk_eigen.h"
23 #include ITK_EIGEN(Eigenvalues)
24 #include <numeric>
25 #include <vector>
26 // For GetPointerToMatrixData
27 #include "vnl/vnl_matrix.h"
28 #include "vnl/vnl_matrix_fixed.h"
29 #include "itkMatrix.h"
30 
31 namespace itk
32 {
33 namespace detail
34 {
35 /* Helper functions returning pointer to matrix data for different types. */
36 template <typename TValueType, unsigned int NRows, unsigned int NCols>
37 const TValueType *
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 
124 {
125  OrderByValue = 1,
127  DoNotOrder
128 };
129 
165 template <typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
166 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
167 {
168 public:
171 #if !defined(ITK_LEGACY_REMOVE)
172  // We need to expose the enum values at the class level
173  // for backwards compatibility
177 #endif
178 
180  : m_OrderEigenValues(EigenValueOrderEnum::OrderByValue)
181  {}
182 
183  SymmetricEigenAnalysis(const unsigned int dimension)
184  :
185 
186  m_Dimension(dimension)
187  , m_Order(dimension)
188  , m_OrderEigenValues(EigenValueOrderEnum::OrderByValue)
189  {}
190 
191  ~SymmetricEigenAnalysis() = default;
192 
193  using MatrixType = TMatrix;
194  using EigenMatrixType = TEigenMatrix;
195  using VectorType = TVector;
196 
210  unsigned int
211  ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const;
212 
233  unsigned int
234  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
235 
236 
238  void
239  SetOrder(const unsigned int n)
240  {
241  m_Order = n;
242  }
243 
247  unsigned int
248  GetOrder() const
249  {
250  return m_Order;
251  }
252 
256  void
257  SetOrderEigenValues(const bool b)
258  {
259  if (b)
260  {
261  m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
262  }
263  else
264  {
265  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
266  }
267  }
269 
270  bool
272  {
273  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
274  }
275 
279  void
281  {
282  if (b)
283  {
284  m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
285  }
286  else
287  {
288  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
289  }
290  }
292 
293  bool
295  {
296  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
297  }
298 
301  void
302  SetDimension(const unsigned int n)
303  {
304  m_Dimension = n;
305  if (m_Order == 0)
306  {
307  m_Order = m_Dimension;
308  }
309  }
311 
314  unsigned int
315  GetDimension() const
316  {
317  return m_Dimension;
318  }
319 
321  void
322  SetUseEigenLibrary(const bool input)
323  {
324  m_UseEigenLibrary = input;
325  }
326  void
328  {
329  m_UseEigenLibrary = true;
330  }
331  void
333  {
334  m_UseEigenLibrary = false;
335  }
336  bool
338  {
339  return m_UseEigenLibrary;
340  }
342 
343 private:
344  bool m_UseEigenLibrary{ false };
345  unsigned int m_Dimension{ 0 };
346  unsigned int m_Order{ 0 };
348 
368  void
369  ReduceToTridiagonalMatrix(double * inputMatrix, double * d, double * e, double * e2) const;
370 
392  void
393  ReduceToTridiagonalMatrixAndGetTransformation(double * inputMatrix,
394  double * diagonalElements,
395  double * subDiagonalElements,
396  double * transformMatrix) const;
397 
427  unsigned int
428  ComputeEigenValuesUsingQL(double * d, double * e) const;
429 
467  unsigned int
468  ComputeEigenValuesAndVectorsUsingQL(double * d, double * e, double * z) const;
469 
470  /* Legacy algorithms using thread-safe netlib.
471  * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
472  */
473  unsigned int
474  ComputeEigenValuesLegacy(const TMatrix & A, TVector & EigenValues) const;
475 
476  unsigned int
477  ComputeEigenValuesAndVectorsLegacy(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
478 
479  /* Helper to get the matrix value type for EigenLibMatrix typename.
480  *
481  * If the TMatrix is vnl, the type is in element_type.
482  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
483  *
484  * To use this function:
485  * using ValueType = decltype(this->GetMatrixType(true));
486  *
487  * \note The two `GetMatrixValueType` overloads have different
488  * parameter declarations (`bool` and `...`), to avoid that both
489  * functions are equally good candidates during overload resolution,
490  * in case `element_type` and `ValueType` are both nested types of
491  * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
492  */
493  template <typename QMatrix = TMatrix>
494  auto
495  GetMatrixValueType(bool) const -> typename QMatrix::element_type
496  {
497  return QMatrix::element_type();
498  }
499  template <typename QMatrix = TMatrix>
500  auto
501  GetMatrixValueType(...) const -> typename QMatrix::ValueType
502  {
503  return QMatrix::ValueType();
504  }
505 
506  /* Wrapper that call the right implementation for the type of matrix. */
507  unsigned int
509  TVector & EigenValues,
510  TEigenMatrix & EigenVectors) const
511  {
512  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
513  }
514 
515  /* Implementation detail using EigenLib that performs a copy of the input matrix.
516  *
517  * @param (long) implementation detail argument making this implementation less favourable
518  * to be chosen if alternatives are available.
519  *
520  * @return an unsigned int with no information value (no error code in EigenLib) */
521  template <typename QMatrix>
522  auto
524  TVector & EigenValues,
525  TEigenMatrix & EigenVectors,
526  long) const -> decltype(static_cast<unsigned int>(1))
527  {
528  using ValueType = decltype(GetMatrixValueType(true));
529  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
530  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
531  for (unsigned int row = 0; row < m_Dimension; ++row)
532  {
533  for (unsigned int col = 0; col < m_Dimension; ++col)
534  {
535  inputMatrix(row, col) = A(row, col);
536  }
537  }
538  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
539  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
540  const auto & eigenValues = solver.eigenvalues();
541  /* Column k of the returned matrix is an eigenvector corresponding to
542  * eigenvalue number $ k $ as returned by eigenvalues().
543  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
544  const auto & eigenVectors = solver.eigenvectors();
545 
546  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
547  {
548  auto copyEigenValues = eigenValues;
549  auto copyEigenVectors = eigenVectors;
550  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
551  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
552 
553  for (unsigned int row = 0; row < m_Dimension; ++row)
554  {
555  EigenValues[row] = copyEigenValues[row];
556  for (unsigned int col = 0; col < m_Dimension; ++col)
557  {
558  EigenVectors[row][col] = copyEigenVectors(col, row);
559  }
560  }
561  }
562  else
563  {
564  for (unsigned int row = 0; row < m_Dimension; ++row)
565  {
566  EigenValues[row] = eigenValues[row];
567  for (unsigned int col = 0; col < m_Dimension; ++col)
568  {
569  EigenVectors[row][col] = eigenVectors(col, row);
570  }
571  }
572  }
573  // No error code
574  return 1;
575  }
576 
577 
578  /* Implementation detail using EigenLib that do not peform a copy.
579  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
580  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
581  * should be included.
582  *
583  * @param (bool) implementation detail argument making this implementation the most favourable
584  * to be chosen from all the alternative implementations.
585  *
586  * @return an unsigned int with no information value (no error code in EigenLib) */
587  template <typename QMatrix>
588  auto
590  TVector & EigenValues,
591  TEigenMatrix & EigenVectors,
592  bool) const
593  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
594  {
595  auto pointerToData = GetPointerToMatrixData(A);
596  using PointerType = decltype(pointerToData);
597  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
598  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
599  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
600  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
601  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
602  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
603  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
604  const auto & eigenValues = solver.eigenvalues();
605  /* Column k of the returned matrix is an eigenvector corresponding to
606  * eigenvalue number $ k $ as returned by eigenvalues().
607  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
608  const auto & eigenVectors = solver.eigenvectors();
609  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
610  {
611  auto copyEigenValues = eigenValues;
612  auto copyEigenVectors = eigenVectors;
613  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
614  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
615  for (unsigned int row = 0; row < m_Dimension; ++row)
616  {
617  EigenValues[row] = copyEigenValues[row];
618  for (unsigned int col = 0; col < m_Dimension; ++col)
619  {
620  EigenVectors[row][col] = copyEigenVectors(col, row);
621  }
622  }
623  }
624  else
625  {
626  for (unsigned int row = 0; row < m_Dimension; ++row)
627  {
628  EigenValues[row] = eigenValues[row];
629  for (unsigned int col = 0; col < m_Dimension; ++col)
630  {
631  EigenVectors[row][col] = eigenVectors(col, row);
632  }
633  }
634  }
635  // No error code
636  return 1;
637  }
638 
639  /* Wrapper that call the right implementation for the type of matrix. */
640  unsigned int
641  ComputeEigenValuesWithEigenLibrary(const TMatrix & A, TVector & EigenValues) const
642  {
643  return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
644  }
645 
646  /* Implementation detail using EigenLib that performs a copy of the input matrix.
647  *
648  * @param (long) implementation detail argument making this implementation less favourable
649  * to be chosen if alternatives are available.
650  *
651  * @return an unsigned int with no information value (no error code in EigenLib) */
652  template <typename QMatrix>
653  auto
654  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const
655  -> decltype(static_cast<unsigned int>(1))
656  {
657  using ValueType = decltype(GetMatrixValueType(true));
658  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
659  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
660  for (unsigned int row = 0; row < m_Dimension; ++row)
661  {
662  for (unsigned int col = 0; col < m_Dimension; ++col)
663  {
664  inputMatrix(row, col) = A(row, col);
665  }
666  }
667  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
668  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
669  auto eigenValues = solver.eigenvalues();
670  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
671  {
672  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
673  }
674  for (unsigned int i = 0; i < m_Dimension; ++i)
675  {
676  EigenValues[i] = eigenValues[i];
677  }
678 
679  // No error code
680  return 1;
681  }
682 
683  /* Implementation detail using EigenLib that do not peform a copy.
684  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
685  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
686  * should be included.
687  *
688  * @param (bool) implementation detail argument making this implementation the most favourable
689  * to be chosen from all the alternative implementations.
690  *
691  * @return an unsigned int with no information value (no error code in EigenLib) */
692  template <typename QMatrix>
693  auto
694  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
695  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
696  {
697  auto pointerToData = GetPointerToMatrixData(A);
698  using PointerType = decltype(pointerToData);
699  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
700  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
701  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
702  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
703  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
704  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
705  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
706  auto eigenValues = solver.eigenvalues();
707  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
708  {
709  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
710  }
711  for (unsigned int i = 0; i < m_Dimension; ++i)
712  {
713  EigenValues[i] = eigenValues[i];
714  }
715  // No error code
716  return 1;
717  }
718 };
719 
720 template <typename TMatrix, typename TVector, typename TEigenMatrix>
721 std::ostream &
723 {
724  os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
725  os << " Dimension : " << s.GetDimension() << std::endl;
726  os << " Order : " << s.GetOrder() << std::endl;
727  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
728  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
729  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
730  return os;
731 }
732 
733 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
734 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysisFixedDimension
735 {
736 public:
739 #if !defined(ITK_LEGACY_REMOVE)
740  // We need to expose the enum values at the class level
741  // for backwards compatibility
745 #endif
746 
748  : m_OrderEigenValues(EigenValueOrderEnum::OrderByValue)
749  {}
751 
752  using MatrixType = TMatrix;
753  using EigenMatrixType = TEigenMatrix;
754  using VectorType = TVector;
755 
769  unsigned int
770  ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const
771  {
772  return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
773  }
774 
795  unsigned int
796  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const
797  {
798  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
799  }
800 
801  void
802  SetOrderEigenValues(const bool b)
803  {
804  if (b)
805  {
806  m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
807  }
808  else
809  {
810  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
811  }
812  }
813  bool
815  {
816  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
817  }
818  void
820  {
821  if (b)
822  {
823  m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
824  }
825  else
826  {
827  m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
828  }
829  }
830  bool
832  {
833  return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
834  }
835  constexpr unsigned int
836  GetOrder() const
837  {
838  return VDimension;
839  }
840  constexpr unsigned int
841  GetDimension() const
842  {
843  return VDimension;
844  }
845  constexpr bool
847  {
848  return true;
849  }
850 
851 private:
853 
854  /* Helper to get the matrix value type for EigenLibMatrix typename.
855  *
856  * If the TMatrix is vnl, the type is in element_type.
857  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
858  *
859  * To use this function:
860  * using ValueType = decltype(this->GetMatrixType(true));
861  */
862  template <typename QMatrix = TMatrix>
863  auto
864  GetMatrixValueType(bool) const -> typename QMatrix::element_type
865  {
866  return QMatrix::element_type();
867  }
868  template <typename QMatrix = TMatrix>
869  auto
870  GetMatrixValueType(bool) const -> typename QMatrix::ValueType
871  {
872  return QMatrix::ValueType();
873  }
874 
875  /* Implementation detail using EigenLib that do not peform a copy.
876  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
877  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
878  * should be included.
879  *
880  * @param (bool) implementation detail argument making this implementation the most favourable
881  * to be chosen from all the alternative implementations.
882  *
883  * @return an unsigned int with no information value (no error code in EigenLib) */
884  template <typename QMatrix>
885  auto
887  TVector & EigenValues,
888  TEigenMatrix & EigenVectors,
889  bool) const
890  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
891  {
892  auto pointerToData = GetPointerToMatrixData(A);
893  using PointerType = decltype(pointerToData);
894  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
895  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
896  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
897  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
898  EigenConstMatrixMap inputMatrix(pointerToData);
899  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
900  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
901  const auto & eigenValues = solver.eigenvalues();
902  /* Column k of the returned matrix is an eigenvector corresponding to
903  * eigenvalue number $ k $ as returned by eigenvalues().
904  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
905  const auto & eigenVectors = solver.eigenvectors();
906  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
907  {
908  auto copyEigenValues = eigenValues;
909  auto copyEigenVectors = eigenVectors;
910  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
911  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
912  for (unsigned int row = 0; row < VDimension; ++row)
913  {
914  EigenValues[row] = copyEigenValues[row];
915  for (unsigned int col = 0; col < VDimension; ++col)
916  {
917  EigenVectors[row][col] = copyEigenVectors(col, row);
918  }
919  }
920  }
921  else
922  {
923  for (unsigned int row = 0; row < VDimension; ++row)
924  {
925  EigenValues[row] = eigenValues[row];
926  for (unsigned int col = 0; col < VDimension; ++col)
927  {
928  EigenVectors[row][col] = eigenVectors(col, row);
929  }
930  }
931  }
932  // No error code
933  return 1;
934  }
935 
936  /* Implementation detail using EigenLib that performs a copy of the input matrix.
937  *
938  * @param (long) implementation detail argument making this implementation less favourable
939  * to be chosen if alternatives are available.
940  *
941  * @return an unsigned int with no information value (no error code in EigenLib) */
942  template <typename QMatrix>
943  auto
945  TVector & EigenValues,
946  TEigenMatrix & EigenVectors,
947  long) const -> decltype(static_cast<unsigned int>(1))
948  {
949  using ValueType = decltype(GetMatrixValueType(true));
950  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
951  EigenLibMatrixType inputMatrix;
952  for (unsigned int row = 0; row < VDimension; ++row)
953  {
954  for (unsigned int col = 0; col < VDimension; ++col)
955  {
956  inputMatrix(row, col) = A(row, col);
957  }
958  }
959  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
960  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
961  const auto & eigenValues = solver.eigenvalues();
962  /* Column k of the returned matrix is an eigenvector corresponding to
963  * eigenvalue number $ k $ as returned by eigenvalues().
964  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
965  const auto & eigenVectors = solver.eigenvectors();
966 
967  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
968  {
969  auto copyEigenValues = eigenValues;
970  auto copyEigenVectors = eigenVectors;
971  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
972  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
973 
974  for (unsigned int row = 0; row < VDimension; ++row)
975  {
976  EigenValues[row] = copyEigenValues[row];
977  for (unsigned int col = 0; col < VDimension; ++col)
978  {
979  EigenVectors[row][col] = copyEigenVectors(col, row);
980  }
981  }
982  }
983  else
984  {
985  for (unsigned int row = 0; row < VDimension; ++row)
986  {
987  EigenValues[row] = eigenValues[row];
988  for (unsigned int col = 0; col < VDimension; ++col)
989  {
990  EigenVectors[row][col] = eigenVectors(col, row);
991  }
992  }
993  }
994  // No error code
995  return 1;
996  }
997 
998  /* Implementation detail using EigenLib that performs a copy of the input matrix.
999  *
1000  * @param (long) implementation detail argument making this implementation less favourable
1001  * to be chosen if alternatives are available.
1002  *
1003  * @return an unsigned int with no information value (no error code in EigenLib) */
1004  template <typename QMatrix>
1005  auto
1006  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const
1007  -> decltype(static_cast<unsigned int>(1))
1008  {
1009  using ValueType = decltype(GetMatrixValueType(true));
1010  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1011  EigenLibMatrixType inputMatrix;
1012  for (unsigned int row = 0; row < VDimension; ++row)
1013  {
1014  for (unsigned int col = 0; col < VDimension; ++col)
1015  {
1016  inputMatrix(row, col) = A(row, col);
1017  }
1018  }
1019  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1020  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1021  auto eigenValues = solver.eigenvalues();
1022  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1023  {
1024  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1025  }
1026  for (unsigned int i = 0; i < VDimension; ++i)
1027  {
1028  EigenValues[i] = eigenValues[i];
1029  }
1030 
1031  // No error code
1032  return 1;
1033  }
1034 
1035  /* Implementation detail using EigenLib that do not peform a copy.
1036  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
1037  * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
1038  * should be included.
1039  *
1040  * @param (bool) implementation detail argument making this implementation the most favourable
1041  * to be chosen from all the alternative implementations.
1042  *
1043  * @return an unsigned int with no information value (no error code in EigenLib) */
1044  template <typename QMatrix>
1045  auto
1046  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
1047  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
1048  {
1049  auto pointerToData = GetPointerToMatrixData(A);
1050  using PointerType = decltype(pointerToData);
1051  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
1052  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
1053  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1054  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1055  EigenConstMatrixMap inputMatrix(pointerToData);
1056  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1057  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1058  auto eigenValues = solver.eigenvalues();
1059  if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1060  {
1061  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1062  }
1063  for (unsigned int i = 0; i < VDimension; ++i)
1064  {
1065  EigenValues[i] = eigenValues[i];
1066  }
1067  // No error code
1068  return 1;
1069  }
1070 };
1071 
1072 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix>
1073 std::ostream &
1074 operator<<(std::ostream & os,
1076 {
1077  os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1078  os << " Dimension : " << s.GetDimension() << std::endl;
1079  os << " Order : " << s.GetOrder() << std::endl;
1080  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
1081  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
1082  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
1083  return os;
1084 }
1085 } // end namespace itk
1086 
1087 #ifndef ITK_MANUAL_INSTANTIATION
1088 # include "itkSymmetricEigenAnalysis.hxx"
1089 #endif
1090 
1091 #endif
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::VectorType
TOutput VectorType
Definition: itkSymmetricEigenAnalysis.h:754
itk::uint8_t
::uint8_t uint8_t
Definition: itkIntTypes.h:29
itk::SymmetricEigenAnalysis::SymmetricEigenAnalysis
SymmetricEigenAnalysis()
Definition: itkSymmetricEigenAnalysis.h:179
EigenValueOrderEnum
itk::SymmetricEigenAnalysisFixedDimension::m_OrderEigenValues
EigenValueOrderType m_OrderEigenValues
Definition: itkSymmetricEigenAnalysis.h:852
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectors
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:796
itk::SymmetricEigenAnalysis::SetDimension
void SetDimension(const unsigned int n)
Definition: itkSymmetricEigenAnalysis.h:302
itk::SymmetricEigenAnalysis::GetOrder
unsigned int GetOrder() const
Definition: itkSymmetricEigenAnalysis.h:248
itk::Matrix::GetVnlMatrix
InternalMatrixType & GetVnlMatrix()
Definition: itkMatrix.h:177
itkMatrix.h
itk::operator<<
std::ostream & operator<<(std::ostream &os, const Array< TValue > &arr)
Definition: itkArray.h:211
itk::SymmetricEigenAnalysis::GetUseEigenLibrary
bool GetUseEigenLibrary() const
Definition: itkSymmetricEigenAnalysis.h:337
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:523
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:589
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:694
itk::SymmetricEigenAnalysis::SetUseEigenLibraryOn
void SetUseEigenLibraryOn()
Definition: itkSymmetricEigenAnalysis.h:327
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::VectorType
TOutputImage::PixelType VectorType
Definition: itkSymmetricEigenAnalysis.h:195
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:886
itk::SymmetricEigenAnalysis::SymmetricEigenAnalysis
SymmetricEigenAnalysis(const unsigned int dimension)
Definition: itkSymmetricEigenAnalysis.h:183
itk::SymmetricEigenAnalysisFixedDimension::SetOrderEigenValues
void SetOrderEigenValues(const bool b)
Definition: itkSymmetricEigenAnalysis.h:802
itk::detail::permuteColumnsWithSortIndices
void permuteColumnsWithSortIndices(QMatrix &eigenVectors, const std::vector< int > &indicesSortPermutations)
Definition: itkSymmetricEigenAnalysis.h:101
EigenValueOrderType
itk::SymmetricEigenAnalysis::GetOrderEigenValues
bool GetOrderEigenValues() const
Definition: itkSymmetricEigenAnalysis.h:271
itk::SymmetricEigenAnalysisFixedDimension::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::ValueType
Definition: itkSymmetricEigenAnalysis.h:870
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
Definition: itkSymmetricEigenAnalysis.h:495
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibrary
unsigned int ComputeEigenValuesAndVectorsWithEigenLibrary(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:508
itk::EigenValueOrderEnum::OrderByValue
itk::SymmetricEigenAnalysisFixedDimension::SymmetricEigenAnalysisFixedDimension
SymmetricEigenAnalysisFixedDimension()
Definition: itkSymmetricEigenAnalysis.h:747
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:239
itkMacro.h
itk::SymmetricEigenAnalysis::GetDimension
unsigned int GetDimension() const
Definition: itkSymmetricEigenAnalysis.h:315
itk::EigenValueOrderEnum
EigenValueOrderEnum
Definition: itkSymmetricEigenAnalysis.h:123
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:1006
itk::EigenValueOrderEnum::DoNotOrder
itk::EigenValueOrderEnum::OrderByMagnitude
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(...) const -> typename QMatrix::ValueType
Definition: itkSymmetricEigenAnalysis.h:501
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:819
itk::SymmetricEigenAnalysisFixedDimension::GetOrderEigenMagnitudes
bool GetOrderEigenMagnitudes() const
Definition: itkSymmetricEigenAnalysis.h:831
itk::SymmetricEigenAnalysis::SetOrderEigenMagnitudes
void SetOrderEigenMagnitudes(const bool b)
Definition: itkSymmetricEigenAnalysis.h:280
itk::SymmetricEigenAnalysis::SetUseEigenLibraryOff
void SetUseEigenLibraryOff()
Definition: itkSymmetricEigenAnalysis.h:332
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::EigenMatrixType
TInput EigenMatrixType
Definition: itkSymmetricEigenAnalysis.h:753
itk::Matrix
A templated class holding a M x N size Matrix.
Definition: itkMatrix.h:50
itk::SymmetricEigenAnalysisFixedDimension::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
Definition: itkSymmetricEigenAnalysis.h:864
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:1046
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::EigenMatrixType
TInputImage::PixelType EigenMatrixType
Definition: itkSymmetricEigenAnalysis.h:194
itk::SymmetricEigenAnalysis::SetUseEigenLibrary
void SetUseEigenLibrary(const bool input)
Definition: itkSymmetricEigenAnalysis.h:322
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:944
itk::SymmetricEigenAnalysisFixedDimension::GetDimension
constexpr unsigned int GetDimension() const
Definition: itkSymmetricEigenAnalysis.h:841
itk::SymmetricEigenAnalysisFixedDimension
Definition: itkSymmetricEigenAnalysis.h:734
itk::SymmetricEigenAnalysisFixedDimension::GetUseEigenLibrary
constexpr bool GetUseEigenLibrary() const
Definition: itkSymmetricEigenAnalysis.h:846
itk::SymmetricEigenAnalysis::SetOrderEigenValues
void SetOrderEigenValues(const bool b)
Definition: itkSymmetricEigenAnalysis.h:257
itk
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
Definition: itkArray.h:26
itk::SymmetricEigenAnalysisFixedDimension::GetOrderEigenValues
bool GetOrderEigenValues() const
Definition: itkSymmetricEigenAnalysis.h:814
itk::SymmetricEigenAnalysis::m_OrderEigenValues
EigenValueOrderType m_OrderEigenValues
Definition: itkSymmetricEigenAnalysis.h:347
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:654
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:641
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::MatrixType
TInput MatrixType
Definition: itkSymmetricEigenAnalysis.h:752
itk::SymmetricEigenAnalysisFixedDimension::GetOrder
constexpr unsigned int GetOrder() const
Definition: itkSymmetricEigenAnalysis.h:836
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::MatrixType
TInputImage::PixelType MatrixType
Definition: itkSymmetricEigenAnalysis.h:193
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:294
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValues
unsigned int ComputeEigenValues(const TMatrix &A, TVector &EigenValues) const
Definition: itkSymmetricEigenAnalysis.h:770
itk::SymmetricEigenAnalysis
Find Eigen values of a real 2D symmetric matrix. It serves as a thread-safe alternative to the class:...
Definition: itkSymmetricEigenAnalysis.h:166