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 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];
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)
796 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues,
true);
822 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors,
true);
859 constexpr
unsigned int
864 constexpr
unsigned int
886 template <
typename QMatrix = TMatrix>
890 return QMatrix::element_type();
892 template <
typename QMatrix = TMatrix>
896 return QMatrix::ValueType();
908 template <
typename QMatrix>
911 TVector & EigenValues,
912 TEigenMatrix & EigenVectors,
916 using PointerType = decltype(pointerToData);
917 using ValueTypeCV = std::remove_pointer_t<PointerType>;
918 using ValueType = std::remove_cv_t<ValueTypeCV>;
919 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
920 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
921 EigenConstMatrixMap inputMatrix(pointerToData);
922 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
923 EigenSolverType solver(inputMatrix);
924 const auto & eigenValues = solver.eigenvalues();
928 const auto & eigenVectors = solver.eigenvectors();
931 auto copyEigenValues = eigenValues;
932 auto copyEigenVectors = eigenVectors;
935 for (
unsigned int row = 0; row < VDimension; ++row)
937 EigenValues[row] = copyEigenValues[row];
938 for (
unsigned int col = 0; col < VDimension; ++col)
940 EigenVectors[row][col] = copyEigenVectors(col, row);
946 for (
unsigned int row = 0; row < VDimension; ++row)
948 EigenValues[row] = eigenValues[row];
949 for (
unsigned int col = 0; col < VDimension; ++col)
951 EigenVectors[row][col] = eigenVectors(col, row);
965 template <
typename QMatrix>
968 TVector & EigenValues,
969 TEigenMatrix & EigenVectors,
970 long)
const -> decltype(1U)
972 using ValueType = decltype(GetMatrixValueType(
true));
973 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
974 EigenLibMatrixType inputMatrix;
975 for (
unsigned int row = 0; row < VDimension; ++row)
977 for (
unsigned int col = 0; col < VDimension; ++col)
979 inputMatrix(row, col) = A(row, col);
982 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
983 EigenSolverType solver(inputMatrix);
984 const auto & eigenValues = solver.eigenvalues();
988 const auto & eigenVectors = solver.eigenvectors();
992 auto copyEigenValues = eigenValues;
993 auto copyEigenVectors = eigenVectors;
997 for (
unsigned int row = 0; row < VDimension; ++row)
999 EigenValues[row] = copyEigenValues[row];
1000 for (
unsigned int col = 0; col < VDimension; ++col)
1002 EigenVectors[row][col] = copyEigenVectors(col, row);
1008 for (
unsigned int row = 0; row < VDimension; ++row)
1010 EigenValues[row] = eigenValues[row];
1011 for (
unsigned int col = 0; col < VDimension; ++col)
1013 EigenVectors[row][col] = eigenVectors(col, row);
1027 template <
typename QMatrix>
1031 using ValueType = decltype(GetMatrixValueType(
true));
1032 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1033 EigenLibMatrixType inputMatrix;
1034 for (
unsigned int row = 0; row < VDimension; ++row)
1036 for (
unsigned int col = 0; col < VDimension; ++col)
1038 inputMatrix(row, col) = A(row, col);
1041 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1042 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1043 auto eigenValues = solver.eigenvalues();
1048 for (
unsigned int i = 0; i < VDimension; ++i)
1050 EigenValues[i] = eigenValues[i];
1066 template <
typename QMatrix>
1072 using PointerType = decltype(pointerToData);
1073 using ValueTypeCV = std::remove_pointer_t<PointerType>;
1074 using ValueType = std::remove_cv_t<ValueTypeCV>;
1075 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1076 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1077 EigenConstMatrixMap inputMatrix(pointerToData);
1078 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1079 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1080 auto eigenValues = solver.eigenvalues();
1085 for (
unsigned int i = 0; i < VDimension; ++i)
1087 EigenValues[i] = eigenValues[i];
1094 template <
unsigned int VDimension,
typename TMatrix,
typename TVector,
typename TEigenMatrix>
1099 os <<
"[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1100 os <<
" Dimension : " << s.
GetDimension() << std::endl;
1101 os <<
" Order : " << s.
GetOrder() << std::endl;
1109 #ifndef ITK_MANUAL_INSTANTIATION
1110 # include "itkSymmetricEigenAnalysis.hxx"