18 #ifndef itkSymmetricEigenAnalysis_h
19 #define itkSymmetricEigenAnalysis_h
22 #include "itk_eigen.h"
23 #include ITK_EIGEN(Eigenvalues)
27 #include "vnl/vnl_matrix.h"
28 #include "vnl/vnl_matrix_fixed.h"
36 template <
typename TValueType,
unsigned int VRows,
unsigned int VColumns>
40 return inputMatrix.data_block();
42 template <
typename TValueType>
46 return inputMatrix.data_block();
49 template <
typename TValueType,
unsigned int VRows,
unsigned int VColumns>
71 template <
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;
100 template <
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;
139 extern ITKCommon_EXPORT std::ostream &
158 itkGenericExceptionMacro(
"Error: Invalid value for conversion.");
161 #if !defined(ITK_LEGACY_REMOVE)
203 template <
typename TMatrix,
typename TVector,
typename TEigenMatrix = TMatrix>
208 #if !defined(ITK_LEGACY_REMOVE)
216 : m_Dimension(dimension)
240 ComputeEigenValues(
const TMatrix & A, TVector & D)
const;
263 ComputeEigenValuesAndVectors(
const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors)
const;
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 };
398 ReduceToTridiagonalMatrix(
double * a,
double * d,
double *
e,
double * e2)
const;
422 ReduceToTridiagonalMatrixAndGetTransformation(
const double * a,
double * d,
double *
e,
double * z)
const;
454 ComputeEigenValuesUsingQL(
double * d,
double *
e)
const;
494 ComputeEigenValuesAndVectorsUsingQL(
double * d,
double *
e,
double * z)
const;
500 ComputeEigenValuesLegacy(
const TMatrix & A, TVector & D)
const;
503 ComputeEigenValuesAndVectorsLegacy(
const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors)
const;
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 const 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 const 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];
744 template <
typename TMatrix,
typename TVector,
typename TEigenMatrix>
748 os <<
"[ClassType: SymmetricEigenAnalysis]" << std::endl;
750 os <<
" Order : " << s.
GetOrder() << std::endl;
757 template <
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 const 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 const 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];
1091 template <
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"