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 NRows,
unsigned int NCols>
40 return inputMatrix.data_block();
42 template <
typename TValueType>
46 return inputMatrix.data_block();
49 template <
typename TValueType,
unsigned int NRows,
unsigned int NCols>
71 template <
typename TArray>
75 std::vector<int> indicesSortPermutations(numberOfElements, 0);
76 std::iota(std::begin(indicesSortPermutations), std::end(indicesSortPermutations), 0);
80 std::begin(indicesSortPermutations),
81 std::end(indicesSortPermutations),
82 [&eigenValues](
unsigned int a,
unsigned int b) {
return std::abs(eigenValues[a]) < std::abs(eigenValues[b]); });
83 auto tmpCopy = eigenValues;
84 for (
unsigned int i = 0; i < numberOfElements; ++i)
86 eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
88 return indicesSortPermutations;
99 template <
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;
138 extern ITKCommon_EXPORT std::ostream &
157 itkGenericExceptionMacro(<<
"Error: Invalid value for conversion.");
160 #if !defined(ITK_LEGACY_REMOVE)
202 template <
typename TMatrix,
typename TVector,
typename TEigenMatrix = TMatrix>
207 #if !defined(ITK_LEGACY_REMOVE)
215 : m_Dimension(dimension)
239 ComputeEigenValues(
const TMatrix & A, TVector & EigenValues)
const;
262 ComputeEigenValuesAndVectors(
const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors)
const;
335 m_Order = m_Dimension;
352 m_UseEigenLibrary = input;
357 m_UseEigenLibrary =
true;
362 m_UseEigenLibrary =
false;
367 return m_UseEigenLibrary;
372 bool m_UseEigenLibrary{
false };
373 unsigned int m_Dimension{ 0 };
374 unsigned int m_Order{ 0 };
397 ReduceToTridiagonalMatrix(
double * inputMatrix,
double * d,
double *
e,
double * e2)
const;
421 ReduceToTridiagonalMatrixAndGetTransformation(
double * inputMatrix,
422 double * diagonalElements,
423 double * subDiagonalElements,
424 double * transformMatrix)
const;
456 ComputeEigenValuesUsingQL(
double * d,
double *
e)
const;
496 ComputeEigenValuesAndVectorsUsingQL(
double * d,
double *
e,
double * z)
const;
502 ComputeEigenValuesLegacy(
const TMatrix & A, TVector & EigenValues)
const;
505 ComputeEigenValuesAndVectorsLegacy(
const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors)
const;
521 template <
typename QMatrix = TMatrix>
525 return QMatrix::element_type();
527 template <
typename QMatrix = TMatrix>
531 return QMatrix::ValueType();
537 TVector & EigenValues,
538 TEigenMatrix & EigenVectors)
const
540 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors,
true);
549 template <
typename QMatrix>
552 TVector & EigenValues,
553 TEigenMatrix & EigenVectors,
554 long)
const -> decltype(static_cast<unsigned int>(1))
556 using ValueType = decltype(GetMatrixValueType(
true));
557 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
558 EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
559 for (
unsigned int row = 0; row < m_Dimension; ++row)
561 for (
unsigned int col = 0; col < m_Dimension; ++col)
563 inputMatrix(row, col) = A(row, col);
566 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
567 EigenSolverType solver(inputMatrix);
568 const auto & eigenValues = solver.eigenvalues();
572 const auto & eigenVectors = solver.eigenvectors();
576 auto copyEigenValues = eigenValues;
577 auto copyEigenVectors = eigenVectors;
581 for (
unsigned int row = 0; row < m_Dimension; ++row)
583 EigenValues[row] = copyEigenValues[row];
584 for (
unsigned int col = 0; col < m_Dimension; ++col)
586 EigenVectors[row][col] = copyEigenVectors(col, row);
592 for (
unsigned int row = 0; row < m_Dimension; ++row)
594 EigenValues[row] = eigenValues[row];
595 for (
unsigned int col = 0; col < m_Dimension; ++col)
597 EigenVectors[row][col] = eigenVectors(col, row);
615 template <
typename QMatrix>
618 TVector & EigenValues,
619 TEigenMatrix & EigenVectors,
624 using PointerType = decltype(pointerToData);
625 using ValueTypeCV =
typename std::remove_pointer<PointerType>::type;
626 using ValueType =
typename std::remove_cv<ValueTypeCV>::type;
627 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
628 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
629 EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
630 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
631 EigenSolverType solver(inputMatrix);
632 const auto & eigenValues = solver.eigenvalues();
636 const auto & eigenVectors = solver.eigenvectors();
639 auto copyEigenValues = eigenValues;
640 auto copyEigenVectors = eigenVectors;
643 for (
unsigned int row = 0; row < m_Dimension; ++row)
645 EigenValues[row] = copyEigenValues[row];
646 for (
unsigned int col = 0; col < m_Dimension; ++col)
648 EigenVectors[row][col] = copyEigenVectors(col, row);
654 for (
unsigned int row = 0; row < m_Dimension; ++row)
656 EigenValues[row] = eigenValues[row];
657 for (
unsigned int col = 0; col < m_Dimension; ++col)
659 EigenVectors[row][col] = eigenVectors(col, row);
671 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues,
true);
680 template <
typename QMatrix>
683 -> decltype(static_cast<unsigned int>(1))
685 using ValueType = decltype(GetMatrixValueType(
true));
686 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
687 EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
688 for (
unsigned int row = 0; row < m_Dimension; ++row)
690 for (
unsigned int col = 0; col < m_Dimension; ++col)
692 inputMatrix(row, col) = A(row, col);
695 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
696 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
697 auto eigenValues = solver.eigenvalues();
702 for (
unsigned int i = 0; i < m_Dimension; ++i)
704 EigenValues[i] = eigenValues[i];
720 template <
typename QMatrix>
726 using PointerType = decltype(pointerToData);
727 using ValueTypeCV =
typename std::remove_pointer<PointerType>::type;
728 using ValueType =
typename std::remove_cv<ValueTypeCV>::type;
729 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
730 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
731 EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
732 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
733 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
734 auto eigenValues = solver.eigenvalues();
739 for (
unsigned int i = 0; i < m_Dimension; ++i)
741 EigenValues[i] = eigenValues[i];
748 template <
typename TMatrix,
typename TVector,
typename TEigenMatrix>
752 os <<
"[ClassType: SymmetricEigenAnalysis]" << std::endl;
754 os <<
" Order : " << s.
GetOrder() << std::endl;
761 template <
unsigned int VDimension,
typename TMatrix,
typename TVector,
typename TEigenMatrix = TMatrix>
765 #if !defined(ITK_LEGACY_REMOVE)
769 #if !defined(ITK_LEGACY_REMOVE)
800 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues,
true);
826 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors,
true);
863 constexpr
unsigned int
868 constexpr
unsigned int
890 template <
typename QMatrix = TMatrix>
894 return QMatrix::element_type();
896 template <
typename QMatrix = TMatrix>
900 return QMatrix::ValueType();
912 template <
typename QMatrix>
915 TVector & EigenValues,
916 TEigenMatrix & EigenVectors,
921 using PointerType = decltype(pointerToData);
922 using ValueTypeCV =
typename std::remove_pointer<PointerType>::type;
923 using ValueType =
typename std::remove_cv<ValueTypeCV>::type;
924 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
925 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
926 EigenConstMatrixMap inputMatrix(pointerToData);
927 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
928 EigenSolverType solver(inputMatrix);
929 const auto & eigenValues = solver.eigenvalues();
933 const auto & eigenVectors = solver.eigenvectors();
936 auto copyEigenValues = eigenValues;
937 auto copyEigenVectors = eigenVectors;
940 for (
unsigned int row = 0; row < VDimension; ++row)
942 EigenValues[row] = copyEigenValues[row];
943 for (
unsigned int col = 0; col < VDimension; ++col)
945 EigenVectors[row][col] = copyEigenVectors(col, row);
951 for (
unsigned int row = 0; row < VDimension; ++row)
953 EigenValues[row] = eigenValues[row];
954 for (
unsigned int col = 0; col < VDimension; ++col)
956 EigenVectors[row][col] = eigenVectors(col, row);
970 template <
typename QMatrix>
973 TVector & EigenValues,
974 TEigenMatrix & EigenVectors,
975 long)
const -> decltype(static_cast<unsigned int>(1))
977 using ValueType = decltype(GetMatrixValueType(
true));
978 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
979 EigenLibMatrixType inputMatrix;
980 for (
unsigned int row = 0; row < VDimension; ++row)
982 for (
unsigned int col = 0; col < VDimension; ++col)
984 inputMatrix(row, col) = A(row, col);
987 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
988 EigenSolverType solver(inputMatrix);
989 const auto & eigenValues = solver.eigenvalues();
993 const auto & eigenVectors = solver.eigenvectors();
997 auto copyEigenValues = eigenValues;
998 auto copyEigenVectors = eigenVectors;
1002 for (
unsigned int row = 0; row < VDimension; ++row)
1004 EigenValues[row] = copyEigenValues[row];
1005 for (
unsigned int col = 0; col < VDimension; ++col)
1007 EigenVectors[row][col] = copyEigenVectors(col, row);
1013 for (
unsigned int row = 0; row < VDimension; ++row)
1015 EigenValues[row] = eigenValues[row];
1016 for (
unsigned int col = 0; col < VDimension; ++col)
1018 EigenVectors[row][col] = eigenVectors(col, row);
1032 template <
typename QMatrix>
1035 -> decltype(static_cast<unsigned int>(1))
1037 using ValueType = decltype(GetMatrixValueType(
true));
1038 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1039 EigenLibMatrixType inputMatrix;
1040 for (
unsigned int row = 0; row < VDimension; ++row)
1042 for (
unsigned int col = 0; col < VDimension; ++col)
1044 inputMatrix(row, col) = A(row, col);
1047 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1048 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1049 auto eigenValues = solver.eigenvalues();
1054 for (
unsigned int i = 0; i < VDimension; ++i)
1056 EigenValues[i] = eigenValues[i];
1072 template <
typename QMatrix>
1078 using PointerType = decltype(pointerToData);
1079 using ValueTypeCV =
typename std::remove_pointer<PointerType>::type;
1080 using ValueType =
typename std::remove_cv<ValueTypeCV>::type;
1081 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1082 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1083 EigenConstMatrixMap inputMatrix(pointerToData);
1084 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1085 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1086 auto eigenValues = solver.eigenvalues();
1091 for (
unsigned int i = 0; i < VDimension; ++i)
1093 EigenValues[i] = eigenValues[i];
1100 template <
unsigned int VDimension,
typename TMatrix,
typename TVector,
typename TEigenMatrix>
1105 os <<
"[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1106 os <<
" Dimension : " << s.
GetDimension() << std::endl;
1107 os <<
" Order : " << s.
GetOrder() << std::endl;
1115 #ifndef ITK_MANUAL_INSTANTIATION
1116 # include "itkSymmetricEigenAnalysis.hxx"