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;
134 OrderByMagnitude = 2,
139extern ITKCommon_EXPORT std::ostream &
158 itkGenericExceptionMacro(
"Error: Invalid value for conversion.");
161#if !defined(ITK_LEGACY_REMOVE)
203template <
typename TMatrix,
typename TVector,
typename TEigenMatrix = TMatrix>
208#if !defined(ITK_LEGACY_REMOVE)
216 : m_Dimension(dimension)
336 m_Order = m_Dimension;
353 m_UseEigenLibrary = input;
358 m_UseEigenLibrary =
true;
363 m_UseEigenLibrary =
false;
368 return m_UseEigenLibrary;
373 bool m_UseEigenLibrary{
false };
374 unsigned int m_Dimension{ 0 };
375 unsigned int m_Order{ 0 };
519 template <
typename QMatrix = TMatrix>
523 return QMatrix::element_type();
525 template <
typename QMatrix = TMatrix>
529 return QMatrix::ValueType();
535 TVector & EigenValues,
536 TEigenMatrix & EigenVectors)
const
538 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors,
true);
547 template <
typename QMatrix>
550 TVector & EigenValues,
551 TEigenMatrix & EigenVectors,
552 long)
const ->
decltype(1U)
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)
559 for (
unsigned int col = 0; col < m_Dimension; ++col)
561 inputMatrix(row, col) = A(row, col);
564 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
565 EigenSolverType solver(inputMatrix);
566 const auto & eigenValues = solver.eigenvalues();
570 const auto & eigenVectors = solver.eigenvectors();
574 auto copyEigenValues = eigenValues;
575 auto copyEigenVectors = eigenVectors;
579 for (
unsigned int row = 0; row < m_Dimension; ++row)
581 EigenValues[row] = copyEigenValues[row];
582 for (
unsigned int col = 0; col < m_Dimension; ++col)
584 EigenVectors[row][col] = copyEigenVectors(col, row);
590 for (
unsigned int row = 0; row < m_Dimension; ++row)
592 EigenValues[row] = eigenValues[row];
593 for (
unsigned int col = 0; col < m_Dimension; ++col)
595 EigenVectors[row][col] = eigenVectors(col, row);
613 template <
typename QMatrix>
616 TVector & EigenValues,
617 TEigenMatrix & EigenVectors,
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);
629 const auto & eigenValues = solver.eigenvalues();
633 const auto & eigenVectors = solver.eigenvectors();
636 auto copyEigenValues = eigenValues;
637 auto copyEigenVectors = eigenVectors;
640 for (
unsigned int row = 0; row < m_Dimension; ++row)
642 EigenValues[row] = copyEigenValues[row];
643 for (
unsigned int col = 0; col < m_Dimension; ++col)
645 EigenVectors[row][col] = copyEigenVectors(col, row);
651 for (
unsigned int row = 0; row < m_Dimension; ++row)
653 EigenValues[row] = eigenValues[row];
654 for (
unsigned int col = 0; col < m_Dimension; ++col)
656 EigenVectors[row][col] = eigenVectors(col, row);
668 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues,
true);
677 template <
typename QMatrix>
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)
686 for (
unsigned int col = 0; col < m_Dimension; ++col)
688 inputMatrix(row, col) = A(row, col);
691 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
692 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
693 auto eigenValues = solver.eigenvalues();
698 for (
unsigned int i = 0; i < m_Dimension; ++i)
700 EigenValues[i] = eigenValues[i];
716 template <
typename QMatrix>
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();
735 for (
unsigned int i = 0; i < m_Dimension; ++i)
737 EigenValues[i] = eigenValues[i];
744template <
typename TMatrix,
typename TVector,
typename TEigenMatrix>
748 os <<
"[ClassType: SymmetricEigenAnalysis]" << std::endl;
750 os <<
" Order : " << s.
GetOrder() << std::endl;
757template <
unsigned int VDimension,
typename TMatrix,
typename TVector,
typename TEigenMatrix = TMatrix>
761#if !defined(ITK_LEGACY_REMOVE)
765#if !defined(ITK_LEGACY_REMOVE)
793 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues,
true);
819 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors,
true);
856 constexpr unsigned int
861 constexpr unsigned int
883 template <
typename QMatrix = TMatrix>
887 return QMatrix::element_type();
889 template <
typename QMatrix = TMatrix>
893 return QMatrix::ValueType();
905 template <
typename QMatrix>
908 TVector & EigenValues,
909 TEigenMatrix & EigenVectors,
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);
921 const auto & eigenValues = solver.eigenvalues();
925 const auto & eigenVectors = solver.eigenvectors();
928 auto copyEigenValues = eigenValues;
929 auto copyEigenVectors = eigenVectors;
932 for (
unsigned int row = 0; row < VDimension; ++row)
934 EigenValues[row] = copyEigenValues[row];
935 for (
unsigned int col = 0; col < VDimension; ++col)
937 EigenVectors[row][col] = copyEigenVectors(col, row);
943 for (
unsigned int row = 0; row < VDimension; ++row)
945 EigenValues[row] = eigenValues[row];
946 for (
unsigned int col = 0; col < VDimension; ++col)
948 EigenVectors[row][col] = eigenVectors(col, row);
962 template <
typename QMatrix>
965 TVector & EigenValues,
966 TEigenMatrix & EigenVectors,
967 long)
const ->
decltype(1U)
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)
974 for (
unsigned int col = 0; col < VDimension; ++col)
976 inputMatrix(row, col) = A(row, col);
979 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
980 EigenSolverType solver(inputMatrix);
981 const auto & eigenValues = solver.eigenvalues();
985 const auto & eigenVectors = solver.eigenvectors();
989 auto copyEigenValues = eigenValues;
990 auto copyEigenVectors = eigenVectors;
994 for (
unsigned int row = 0; row < VDimension; ++row)
996 EigenValues[row] = copyEigenValues[row];
997 for (
unsigned int col = 0; col < VDimension; ++col)
999 EigenVectors[row][col] = copyEigenVectors(col, row);
1005 for (
unsigned int row = 0; row < VDimension; ++row)
1007 EigenValues[row] = eigenValues[row];
1008 for (
unsigned int col = 0; col < VDimension; ++col)
1010 EigenVectors[row][col] = eigenVectors(col, row);
1024 template <
typename QMatrix>
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)
1033 for (
unsigned int col = 0; col < VDimension; ++col)
1035 inputMatrix(row, col) = A(row, col);
1038 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1039 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1040 auto eigenValues = solver.eigenvalues();
1045 for (
unsigned int i = 0; i < VDimension; ++i)
1047 EigenValues[i] = eigenValues[i];
1063 template <
typename QMatrix>
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();
1082 for (
unsigned int i = 0; i < VDimension; ++i)
1084 EigenValues[i] = eigenValues[i];
1091template <
unsigned int VDimension,
typename TMatrix,
typename TVector,
typename TEigenMatrix>
1096 os <<
"[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1097 os <<
" Dimension : " << s.
GetDimension() << std::endl;
1098 os <<
" Order : " << s.
GetOrder() << std::endl;
1106#ifndef ITK_MANUAL_INSTANTIATION
1107# 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, TVector &EigenValues, long) const -> decltype(1U)
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
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, TVector &EigenValues, TEigenMatrix &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
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)
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, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(1U)
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &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()
static constexpr double e
const TValueType * GetPointerToMatrixData(const itk::Matrix< TValueType, VRows, VColumns > &inputMatrix)
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)