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);
78 std::sort(std::begin(indicesSortPermutations),
79 std::end(indicesSortPermutations),
80 [&eigenValues](
unsigned int a,
unsigned int b) {
83 auto tmpCopy = eigenValues;
84 for (
unsigned int i = 0; i < numberOfElements; ++i)
86 eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
88 return indicesSortPermutations;
99template <
typename QMatrix>
103 using EigenLibPermutationMatrix = Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
104 auto numberOfElements = indicesSortPermutations.size();
107 EigenLibPermutationMatrix perm(numberOfElements);
109 std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(), perm.indices().data());
111 eigenVectors = eigenVectors * perm;
138extern ITKCommon_EXPORT std::ostream &
149 return EigenValueOrderEnum::OrderByValue;
151 return EigenValueOrderEnum::OrderByMagnitude;
153 return EigenValueOrderEnum::DoNotOrder;
157 itkGenericExceptionMacro(
"Error: Invalid value for conversion.");
160#if !defined(ITK_LEGACY_REMOVE)
163static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
201template <
typename TMatrix,
typename TVector,
typename TEigenMatrix = TMatrix>
206#if !defined(ITK_LEGACY_REMOVE)
512 template <
typename QMatrix = TMatrix>
516 return QMatrix::element_type();
518 template <
typename QMatrix = TMatrix>
522 return QMatrix::ValueType();
528 TVector & EigenValues,
529 TEigenMatrix & EigenVectors)
const
540 template <
typename QMatrix>
543 TVector & EigenValues,
544 TEigenMatrix & EigenVectors,
545 long)
const ->
decltype(1U)
548 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
550 for (
unsigned int row = 0; row <
m_Dimension; ++row)
552 for (
unsigned int col = 0; col <
m_Dimension; ++col)
554 inputMatrix(row, col) = A(row, col);
557 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
558 const EigenSolverType solver(inputMatrix);
559 const auto & eigenValues = solver.eigenvalues();
563 const auto & eigenVectors = solver.eigenvectors();
567 auto copyEigenValues = eigenValues;
568 auto copyEigenVectors = eigenVectors;
572 for (
unsigned int row = 0; row <
m_Dimension; ++row)
574 EigenValues[row] = copyEigenValues[row];
575 for (
unsigned int col = 0; col <
m_Dimension; ++col)
577 EigenVectors[row][col] = copyEigenVectors(col, row);
583 for (
unsigned int row = 0; row <
m_Dimension; ++row)
585 EigenValues[row] = eigenValues[row];
586 for (
unsigned int col = 0; col <
m_Dimension; ++col)
588 EigenVectors[row][col] = eigenVectors(col, row);
606 template <
typename QMatrix>
609 TVector & EigenValues,
610 TEigenMatrix & EigenVectors,
611 bool)
const ->
decltype(GetPointerToMatrixData(A), 1U)
613 auto pointerToData = GetPointerToMatrixData(A);
614 using PointerType =
decltype(pointerToData);
615 using ValueTypeCV = std::remove_pointer_t<PointerType>;
616 using ValueType = std::remove_cv_t<ValueTypeCV>;
617 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
618 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
620 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
621 EigenSolverType solver(inputMatrix);
622 const auto & eigenValues = solver.eigenvalues();
626 const auto & eigenVectors = solver.eigenvectors();
629 auto copyEigenValues = eigenValues;
630 auto copyEigenVectors = eigenVectors;
633 for (
unsigned int row = 0; row <
m_Dimension; ++row)
635 EigenValues[row] = copyEigenValues[row];
636 for (
unsigned int col = 0; col <
m_Dimension; ++col)
638 EigenVectors[row][col] = copyEigenVectors(col, row);
644 for (
unsigned int row = 0; row <
m_Dimension; ++row)
646 EigenValues[row] = eigenValues[row];
647 for (
unsigned int col = 0; col <
m_Dimension; ++col)
649 EigenVectors[row][col] = eigenVectors(col, row);
670 template <
typename QMatrix>
675 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
677 for (
unsigned int row = 0; row <
m_Dimension; ++row)
679 for (
unsigned int col = 0; col <
m_Dimension; ++col)
681 inputMatrix(row, col) = A(row, col);
684 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
685 const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
686 auto eigenValues = solver.eigenvalues();
693 EigenValues[i] = eigenValues[i];
709 template <
typename QMatrix>
712 ->
decltype(GetPointerToMatrixData(A), 1U)
714 auto pointerToData = GetPointerToMatrixData(A);
715 using PointerType =
decltype(pointerToData);
716 using ValueTypeCV = std::remove_pointer_t<PointerType>;
717 using ValueType = std::remove_cv_t<ValueTypeCV>;
718 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
719 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
721 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
722 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
723 auto eigenValues = solver.eigenvalues();
730 EigenValues[i] = eigenValues[i];
737template <
typename TMatrix,
typename TVector,
typename TEigenMatrix>
741 os <<
"[ClassType: SymmetricEigenAnalysis]" << std::endl;
743 os <<
" Order : " << s.
GetOrder() << std::endl;
750template <
unsigned int VDimension,
typename TMatrix,
typename TVector,
typename TEigenMatrix = TMatrix>
754#if !defined(ITK_LEGACY_REMOVE)
758#if !defined(ITK_LEGACY_REMOVE)
762 static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
849 constexpr unsigned int
854 constexpr unsigned int
876 template <
typename QMatrix = TMatrix>
880 return QMatrix::element_type();
882 template <
typename QMatrix = TMatrix>
886 return QMatrix::ValueType();
898 template <
typename QMatrix>
901 TVector & EigenValues,
902 TEigenMatrix & EigenVectors,
903 bool)
const ->
decltype(GetPointerToMatrixData(A), 1U)
905 auto pointerToData = GetPointerToMatrixData(A);
906 using PointerType =
decltype(pointerToData);
907 using ValueTypeCV = std::remove_pointer_t<PointerType>;
908 using ValueType = std::remove_cv_t<ValueTypeCV>;
909 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
910 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
911 EigenConstMatrixMap inputMatrix(pointerToData);
912 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
913 EigenSolverType solver(inputMatrix);
914 const auto & eigenValues = solver.eigenvalues();
918 const auto & eigenVectors = solver.eigenvectors();
921 auto copyEigenValues = eigenValues;
922 auto copyEigenVectors = eigenVectors;
925 for (
unsigned int row = 0; row < VDimension; ++row)
927 EigenValues[row] = copyEigenValues[row];
928 for (
unsigned int col = 0; col < VDimension; ++col)
930 EigenVectors[row][col] = copyEigenVectors(col, row);
936 for (
unsigned int row = 0; row < VDimension; ++row)
938 EigenValues[row] = eigenValues[row];
939 for (
unsigned int col = 0; col < VDimension; ++col)
941 EigenVectors[row][col] = eigenVectors(col, row);
955 template <
typename QMatrix>
958 TVector & EigenValues,
959 TEigenMatrix & EigenVectors,
960 long)
const ->
decltype(1U)
963 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
964 EigenLibMatrixType inputMatrix;
965 for (
unsigned int row = 0; row < VDimension; ++row)
967 for (
unsigned int col = 0; col < VDimension; ++col)
969 inputMatrix(row, col) = A(row, col);
972 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
973 const EigenSolverType solver(inputMatrix);
974 const auto & eigenValues = solver.eigenvalues();
978 const auto & eigenVectors = solver.eigenvectors();
982 auto copyEigenValues = eigenValues;
983 auto copyEigenVectors = eigenVectors;
987 for (
unsigned int row = 0; row < VDimension; ++row)
989 EigenValues[row] = copyEigenValues[row];
990 for (
unsigned int col = 0; col < VDimension; ++col)
992 EigenVectors[row][col] = copyEigenVectors(col, row);
998 for (
unsigned int row = 0; row < VDimension; ++row)
1000 EigenValues[row] = eigenValues[row];
1001 for (
unsigned int col = 0; col < VDimension; ++col)
1003 EigenVectors[row][col] = eigenVectors(col, row);
1017 template <
typename QMatrix>
1022 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1023 EigenLibMatrixType inputMatrix;
1024 for (
unsigned int row = 0; row < VDimension; ++row)
1026 for (
unsigned int col = 0; col < VDimension; ++col)
1028 inputMatrix(row, col) = A(row, col);
1031 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1032 const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1033 auto eigenValues = solver.eigenvalues();
1038 for (
unsigned int i = 0; i < VDimension; ++i)
1040 EigenValues[i] = eigenValues[i];
1056 template <
typename QMatrix>
1059 ->
decltype(GetPointerToMatrixData(A), 1U)
1061 auto pointerToData = GetPointerToMatrixData(A);
1062 using PointerType =
decltype(pointerToData);
1063 using ValueTypeCV = std::remove_pointer_t<PointerType>;
1064 using ValueType = std::remove_cv_t<ValueTypeCV>;
1065 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1066 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1067 EigenConstMatrixMap inputMatrix(pointerToData);
1068 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1069 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1070 auto eigenValues = solver.eigenvalues();
1075 for (
unsigned int i = 0; i < VDimension; ++i)
1077 EigenValues[i] = eigenValues[i];
1084template <
unsigned int VDimension,
typename TMatrix,
typename TVector,
typename TEigenMatrix>
1089 os <<
"[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1090 os <<
" Dimension : " << s.
GetDimension() << std::endl;
1091 os <<
" Order : " << s.
GetOrder() << std::endl;
1099#ifndef ITK_MANUAL_INSTANTIATION
1100# 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)