18#ifndef itkSymmetricEigenAnalysis_h
19#define itkSymmetricEigenAnalysis_h
23#include ITK_EIGEN(Eigenvalues)
27#include "vnl/vnl_matrix.h"
28#include "vnl/vnl_matrix_fixed.h"
36template <
typename TValueType,
unsigned int VRows,
unsigned int VColumns>
40 return inputMatrix.data_block();
42template <
typename TValueType>
46 return inputMatrix.data_block();
49template <
typename TValueType,
unsigned int VRows,
unsigned int VColumns>
71template <
typename TArray>
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) {
84 auto tmpCopy = eigenValues;
85 for (
unsigned int i = 0; i < numberOfElements; ++i)
87 eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
89 return indicesSortPermutations;
100template <
typename QMatrix>
104 using EigenLibPermutationMatrix = Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
105 auto numberOfElements = indicesSortPermutations.size();
108 EigenLibPermutationMatrix perm(numberOfElements);
110 std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(), perm.indices().data());
112 eigenVectors = eigenVectors * perm;
140extern ITKCommon_EXPORT std::ostream &
151 return EigenValueOrderEnum::OrderByValue;
153 return EigenValueOrderEnum::OrderByMagnitude;
155 return EigenValueOrderEnum::DoNotOrder;
159 itkGenericExceptionMacro(
"Error: Invalid value for conversion.");
162#if !defined(ITK_LEGACY_REMOVE)
165static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
203template <
typename TMatrix,
typename TVector,
typename TEigenMatrix = TMatrix>
208#if !defined(ITK_LEGACY_REMOVE)
513 template <
typename QMatrix = TMatrix>
517 return QMatrix::element_type();
519 template <
typename QMatrix = TMatrix>
523 return QMatrix::ValueType();
529 TVector & EigenValues,
530 TEigenMatrix & EigenVectors)
const
541 template <
typename QMatrix>
544 TVector & EigenValues,
545 TEigenMatrix & EigenVectors,
546 long)
const ->
decltype(1U)
549 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
551 for (
unsigned int row = 0; row <
m_Dimension; ++row)
553 for (
unsigned int col = 0; col <
m_Dimension; ++col)
555 inputMatrix(row, col) = A(row, col);
558 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
559 const EigenSolverType solver(inputMatrix);
560 const auto & eigenValues = solver.eigenvalues();
564 const auto & eigenVectors = solver.eigenvectors();
568 auto copyEigenValues = eigenValues;
569 auto copyEigenVectors = eigenVectors;
573 for (
unsigned int row = 0; row <
m_Dimension; ++row)
575 EigenValues[row] = copyEigenValues[row];
576 for (
unsigned int col = 0; col <
m_Dimension; ++col)
578 EigenVectors[row][col] = copyEigenVectors(col, row);
584 for (
unsigned int row = 0; row <
m_Dimension; ++row)
586 EigenValues[row] = eigenValues[row];
587 for (
unsigned int col = 0; col <
m_Dimension; ++col)
589 EigenVectors[row][col] = eigenVectors(col, row);
607 template <
typename QMatrix>
610 TVector & EigenValues,
611 TEigenMatrix & EigenVectors,
612 bool)
const ->
decltype(GetPointerToMatrixData(A), 1U)
614 auto pointerToData = GetPointerToMatrixData(A);
615 using PointerType =
decltype(pointerToData);
616 using ValueTypeCV = std::remove_pointer_t<PointerType>;
617 using ValueType = std::remove_cv_t<ValueTypeCV>;
618 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
619 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
621 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
622 EigenSolverType solver(inputMatrix);
623 const auto & eigenValues = solver.eigenvalues();
627 const auto & eigenVectors = solver.eigenvectors();
630 auto copyEigenValues = eigenValues;
631 auto copyEigenVectors = eigenVectors;
634 for (
unsigned int row = 0; row <
m_Dimension; ++row)
636 EigenValues[row] = copyEigenValues[row];
637 for (
unsigned int col = 0; col <
m_Dimension; ++col)
639 EigenVectors[row][col] = copyEigenVectors(col, row);
645 for (
unsigned int row = 0; row <
m_Dimension; ++row)
647 EigenValues[row] = eigenValues[row];
648 for (
unsigned int col = 0; col <
m_Dimension; ++col)
650 EigenVectors[row][col] = eigenVectors(col, row);
671 template <
typename QMatrix>
676 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
678 for (
unsigned int row = 0; row <
m_Dimension; ++row)
680 for (
unsigned int col = 0; col <
m_Dimension; ++col)
682 inputMatrix(row, col) = A(row, col);
685 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
686 const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
687 auto eigenValues = solver.eigenvalues();
694 EigenValues[i] = eigenValues[i];
710 template <
typename QMatrix>
713 ->
decltype(GetPointerToMatrixData(A), 1U)
715 auto pointerToData = GetPointerToMatrixData(A);
716 using PointerType =
decltype(pointerToData);
717 using ValueTypeCV = std::remove_pointer_t<PointerType>;
718 using ValueType = std::remove_cv_t<ValueTypeCV>;
719 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
720 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
722 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
723 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
724 auto eigenValues = solver.eigenvalues();
731 EigenValues[i] = eigenValues[i];
738template <
typename TMatrix,
typename TVector,
typename TEigenMatrix>
742 os <<
"[ClassType: SymmetricEigenAnalysis]" << std::endl;
744 os <<
" Order : " << s.
GetOrder() << std::endl;
751template <
unsigned int VDimension,
typename TMatrix,
typename TVector,
typename TEigenMatrix = TMatrix>
755#if !defined(ITK_LEGACY_REMOVE)
759#if !defined(ITK_LEGACY_REMOVE)
763 static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
850 constexpr unsigned int
855 constexpr unsigned int
877 template <
typename QMatrix = TMatrix>
881 return QMatrix::element_type();
883 template <
typename QMatrix = TMatrix>
887 return QMatrix::ValueType();
899 template <
typename QMatrix>
902 TVector & EigenValues,
903 TEigenMatrix & EigenVectors,
904 bool)
const ->
decltype(GetPointerToMatrixData(A), 1U)
906 auto pointerToData = GetPointerToMatrixData(A);
907 using PointerType =
decltype(pointerToData);
908 using ValueTypeCV = std::remove_pointer_t<PointerType>;
909 using ValueType = std::remove_cv_t<ValueTypeCV>;
910 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
911 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
912 EigenConstMatrixMap inputMatrix(pointerToData);
913 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
914 EigenSolverType solver(inputMatrix);
915 const auto & eigenValues = solver.eigenvalues();
919 const auto & eigenVectors = solver.eigenvectors();
922 auto copyEigenValues = eigenValues;
923 auto copyEigenVectors = eigenVectors;
926 for (
unsigned int row = 0; row < VDimension; ++row)
928 EigenValues[row] = copyEigenValues[row];
929 for (
unsigned int col = 0; col < VDimension; ++col)
931 EigenVectors[row][col] = copyEigenVectors(col, row);
937 for (
unsigned int row = 0; row < VDimension; ++row)
939 EigenValues[row] = eigenValues[row];
940 for (
unsigned int col = 0; col < VDimension; ++col)
942 EigenVectors[row][col] = eigenVectors(col, row);
956 template <
typename QMatrix>
959 TVector & EigenValues,
960 TEigenMatrix & EigenVectors,
961 long)
const ->
decltype(1U)
964 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
965 EigenLibMatrixType inputMatrix;
966 for (
unsigned int row = 0; row < VDimension; ++row)
968 for (
unsigned int col = 0; col < VDimension; ++col)
970 inputMatrix(row, col) = A(row, col);
973 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
974 const EigenSolverType solver(inputMatrix);
975 const auto & eigenValues = solver.eigenvalues();
979 const auto & eigenVectors = solver.eigenvectors();
983 auto copyEigenValues = eigenValues;
984 auto copyEigenVectors = eigenVectors;
988 for (
unsigned int row = 0; row < VDimension; ++row)
990 EigenValues[row] = copyEigenValues[row];
991 for (
unsigned int col = 0; col < VDimension; ++col)
993 EigenVectors[row][col] = copyEigenVectors(col, row);
999 for (
unsigned int row = 0; row < VDimension; ++row)
1001 EigenValues[row] = eigenValues[row];
1002 for (
unsigned int col = 0; col < VDimension; ++col)
1004 EigenVectors[row][col] = eigenVectors(col, row);
1018 template <
typename QMatrix>
1023 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1024 EigenLibMatrixType inputMatrix;
1025 for (
unsigned int row = 0; row < VDimension; ++row)
1027 for (
unsigned int col = 0; col < VDimension; ++col)
1029 inputMatrix(row, col) = A(row, col);
1032 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1033 const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1034 auto eigenValues = solver.eigenvalues();
1039 for (
unsigned int i = 0; i < VDimension; ++i)
1041 EigenValues[i] = eigenValues[i];
1057 template <
typename QMatrix>
1060 ->
decltype(GetPointerToMatrixData(A), 1U)
1062 auto pointerToData = GetPointerToMatrixData(A);
1063 using PointerType =
decltype(pointerToData);
1064 using ValueTypeCV = std::remove_pointer_t<PointerType>;
1065 using ValueType = std::remove_cv_t<ValueTypeCV>;
1066 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1067 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1068 EigenConstMatrixMap inputMatrix(pointerToData);
1069 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1070 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1071 auto eigenValues = solver.eigenvalues();
1076 for (
unsigned int i = 0; i < VDimension; ++i)
1078 EigenValues[i] = eigenValues[i];
1085template <
unsigned int VDimension,
typename TMatrix,
typename TVector,
typename TEigenMatrix>
1090 os <<
"[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1091 os <<
" Dimension : " << s.
GetDimension() << std::endl;
1092 os <<
" Order : " << s.
GetOrder() << std::endl;
1100#ifndef ITK_MANUAL_INSTANTIATION
1101# include "itkSymmetricEigenAnalysis.hxx"
A templated class holding a M x N size Matrix.
InternalMatrixType & GetVnlMatrix()
This class contains all enum classes used by SymmetricEigenAnalysis class.
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
bool GetOrderEigenMagnitudes() const
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, EigenValuesArrayType &EigenValues, long) const -> decltype(1U)
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
EigenValueOrderEnum m_OrderEigenValues
constexpr bool GetUseEigenLibrary() const
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
constexpr unsigned int GetOrder() const
void SetOrderEigenValues(const bool b)
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, EigenValuesArrayType &EigenValues, EigenVectorsMatrixType &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
unsigned int ComputeEigenValues(const TMatrix &A, TVector &EigenValues) const
bool GetOrderEigenValues() const
constexpr unsigned int GetDimension() const
void SetOrderEigenMagnitudes(const bool b)
TEigenMatrix EigenMatrixType
auto GetMatrixValueType(bool) const -> typename QMatrix::ValueType
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(1U)
Find Eigen values of a real 2D symmetric matrix. It serves as a thread-safe alternative to the class:...
void SetOrderEigenMagnitudes(const bool b)
bool GetOrderEigenValues() const
void SetUseEigenLibraryOn()
void SetDimension(const unsigned int n)
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
itk::EigenValueOrderEnum EigenValueOrderEnum
unsigned int GetDimension() const
unsigned int ComputeEigenValuesLegacy(const TMatrix &A, TVector &D) const
unsigned int ComputeEigenValuesAndVectorsWithEigenLibrary(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
unsigned int ComputeEigenValuesWithEigenLibrary(const TMatrix &A, TVector &EigenValues) const
unsigned int GetOrder() const
bool GetUseEigenLibrary() const
unsigned int ComputeEigenValuesAndVectorsUsingQL(double *d, double *e, double *z) const
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
EigenValueOrderEnum m_OrderEigenValues
unsigned int ComputeEigenValuesAndVectorsLegacy(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
auto GetMatrixValueType(...) const -> typename QMatrix::ValueType
SymmetricEigenAnalysis(const unsigned int dimension)
bool GetOrderEigenMagnitudes() const
void SetOrder(const unsigned int n)
void SetOrderEigenValues(const bool b)
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
void SetUseEigenLibrary(const bool input)
unsigned int ComputeEigenValuesUsingQL(double *d, double *e) const
SymmetricEigenAnalysis()=default
void ReduceToTridiagonalMatrixAndGetTransformation(const double *a, double *d, double *e, double *z) const
unsigned int ComputeEigenValues(const TMatrix &A, TVector &D) const
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TOutput &EigenValues, TInput &EigenVectors, long) const -> decltype(1U)
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TOutput &EigenValues, long) const -> decltype(1U)
void ReduceToTridiagonalMatrix(double *a, double *d, double *e, double *e2) const
~SymmetricEigenAnalysis()=default
TEigenMatrix EigenMatrixType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
void SetUseEigenLibraryOff()
const TValueType * GetPointerToMatrixData(const vnl_matrix_fixed< TValueType, VRows, VColumns > &inputMatrix)
void permuteColumnsWithSortIndices(QMatrix &eigenVectors, const std::vector< int > &indicesSortPermutations)
std::vector< int > sortEigenValuesByMagnitude(TArray &eigenValues, const unsigned int numberOfElements)
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
ITKCommon_EXPORT std::ostream & operator<<(std::ostream &out, typename AnatomicalOrientation::CoordinateEnum value)
SymmetricEigenAnalysisEnums::EigenValueOrder EigenValueOrderEnum
EigenValueOrderEnum Int2EigenValueOrderEnum(const uint8_t value)