ITK  6.0.0
Insight Toolkit
itkSymmetricEigenAnalysis.h
Go to the documentation of this file.
1/*=========================================================================
2 *
3 * Copyright NumFOCUS
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18#ifndef itkSymmetricEigenAnalysis_h
19#define itkSymmetricEigenAnalysis_h
20
21#include "itkMacro.h"
22#include "itk_eigen.h"
23#include ITK_EIGEN(Eigenvalues)
24#include <numeric>
25#include <vector>
26// For GetPointerToMatrixData
27#include "vnl/vnl_matrix.h"
28#include "vnl/vnl_matrix_fixed.h"
29#include "itkMatrix.h"
30
31namespace itk
32{
33namespace detail
34{
35/* Helper functions returning pointer to matrix data for different types. */
36template <typename TValueType, unsigned int VRows, unsigned int VColumns>
37const TValueType *
38GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, VRows, VColumns> & inputMatrix)
39{
40 return inputMatrix.data_block();
41}
42template <typename TValueType>
43const TValueType *
44GetPointerToMatrixData(const vnl_matrix<TValueType> & inputMatrix)
45{
46 return inputMatrix.data_block();
47}
48
49template <typename TValueType, unsigned int VRows, unsigned int VColumns>
50const TValueType *
52{
53 return inputMatrix.GetVnlMatrix().data_block();
54}
55
71template <typename TArray>
72std::vector<int>
73sortEigenValuesByMagnitude(TArray & eigenValues, const unsigned int numberOfElements)
74{
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) {
82 return itk::Math::abs(eigenValues[a]) < itk::Math::abs(eigenValues[b]);
83 });
84 auto tmpCopy = eigenValues;
85 for (unsigned int i = 0; i < numberOfElements; ++i)
86 {
87 eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
88 }
89 return indicesSortPermutations;
90}
91
100template <typename QMatrix>
101void
102permuteColumnsWithSortIndices(QMatrix & eigenVectors, const std::vector<int> & indicesSortPermutations)
103{
104 using EigenLibPermutationMatrix = Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
105 auto numberOfElements = indicesSortPermutations.size();
106 // Creates a NxN permutation matrix copying our permutation to the matrix indices.
107 // Which holds the 1D array representation of a permutation.
108 EigenLibPermutationMatrix perm(numberOfElements);
109 perm.setIdentity();
110 std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(), perm.indices().data());
111 // Apply it
112 eigenVectors = eigenVectors * perm;
113}
114} // end namespace detail
122{
123public:
131 enum class EigenValueOrder : uint8_t
132 {
133 OrderByValue = 1,
134 OrderByMagnitude = 2,
135 DoNotOrder = 3
136 };
137};
138// Define how to print enumeration
139extern ITKCommon_EXPORT std::ostream &
140 operator<<(std::ostream & out, const SymmetricEigenAnalysisEnums::EigenValueOrder value);
141
143
145Int2EigenValueOrderEnum(const uint8_t value)
146{
147 switch (value)
148 {
149 case 1:
151 case 2:
153 case 3:
155 default:
156 break;
157 }
158 itkGenericExceptionMacro("Error: Invalid value for conversion.");
159}
160
161#if !defined(ITK_LEGACY_REMOVE)
163static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
164static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
165static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
166#endif
167
203template <typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
204class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
205{
206public:
208#if !defined(ITK_LEGACY_REMOVE)
210 using EigenValueOrderType = EigenValueOrderEnum;
211#endif
212
214
215 SymmetricEigenAnalysis(const unsigned int dimension)
216 : m_Dimension(dimension)
217 , m_Order(dimension)
218 {}
219
221
222 using MatrixType = TMatrix;
223 using EigenMatrixType = TEigenMatrix;
224 using VectorType = TVector;
225
239 unsigned int
240 ComputeEigenValues(const TMatrix & A, TVector & D) const;
241
262 unsigned int
263 ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
264
265
267 void
268 SetOrder(const unsigned int n)
269 {
270 m_Order = n;
271 }
272
276 unsigned int
277 GetOrder() const
278 {
279 return m_Order;
280 }
281
285 void
287 {
288 if (b)
289 {
290 m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
291 }
292 else
293 {
294 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
295 }
296 }
299 bool
301 {
302 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
303 }
304
308 void
310 {
311 if (b)
312 {
313 m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
314 }
315 else
316 {
317 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
318 }
319 }
322 bool
324 {
325 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
326 }
327
330 void
331 SetDimension(const unsigned int n)
332 {
333 m_Dimension = n;
334 if (m_Order == 0)
335 {
336 m_Order = m_Dimension;
337 }
338 }
343 unsigned int
345 {
346 return m_Dimension;
347 }
348
350 void
351 SetUseEigenLibrary(const bool input)
352 {
353 m_UseEigenLibrary = input;
354 }
355 void
357 {
358 m_UseEigenLibrary = true;
359 }
360 void
362 {
363 m_UseEigenLibrary = false;
364 }
365 bool
367 {
368 return m_UseEigenLibrary;
369 }
372private:
373 bool m_UseEigenLibrary{ false };
374 unsigned int m_Dimension{ 0 };
375 unsigned int m_Order{ 0 };
377
397 void
398 ReduceToTridiagonalMatrix(double * a, double * d, double * e, double * e2) const;
399
421 void
422 ReduceToTridiagonalMatrixAndGetTransformation(const double * a, double * d, double * e, double * z) const;
423
453 unsigned int
454 ComputeEigenValuesUsingQL(double * d, double * e) const;
455
493 unsigned int
494 ComputeEigenValuesAndVectorsUsingQL(double * d, double * e, double * z) const;
495
496 /* Legacy algorithms using thread-safe netlib.
497 * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
498 */
499 unsigned int
500 ComputeEigenValuesLegacy(const TMatrix & A, TVector & D) const;
501
502 unsigned int
503 ComputeEigenValuesAndVectorsLegacy(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
504
505 /* Helper to get the matrix value type for EigenLibMatrix typename.
506 *
507 * If the TMatrix is vnl, the type is in element_type.
508 * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
509 *
510 * To use this function:
511 * using ValueType = decltype(this->GetMatrixType(true));
512 *
513 * \note The two `GetMatrixValueType` overloads have different
514 * parameter declarations (`bool` and `...`), to avoid that both
515 * functions are equally good candidates during overload resolution,
516 * in case `element_type` and `ValueType` are both nested types of
517 * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
518 */
519 template <typename QMatrix = TMatrix>
520 auto
521 GetMatrixValueType(bool) const -> typename QMatrix::element_type
522 {
523 return QMatrix::element_type();
524 }
525 template <typename QMatrix = TMatrix>
526 auto
527 GetMatrixValueType(...) const -> typename QMatrix::ValueType
528 {
529 return QMatrix::ValueType();
530 }
531
532 /* Wrapper that call the right implementation for the type of matrix. */
533 unsigned int
535 TVector & EigenValues,
536 TEigenMatrix & EigenVectors) const
537 {
538 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
539 }
540
541 /* Implementation detail using EigenLib that performs a copy of the input matrix.
542 *
543 * @param (long) implementation detail argument making this implementation less favourable
544 * to be chosen if alternatives are available.
545 *
546 * @return an unsigned int with no information value (no error code in EigenLib) */
547 template <typename QMatrix>
548 auto
550 TVector & EigenValues,
551 TEigenMatrix & EigenVectors,
552 long) const -> decltype(1U)
553 {
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)
558 {
559 for (unsigned int col = 0; col < m_Dimension; ++col)
560 {
561 inputMatrix(row, col) = A(row, col);
562 }
563 }
564 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
565 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
566 const auto & eigenValues = solver.eigenvalues();
567 /* Column k of the returned matrix is an eigenvector corresponding to
568 * eigenvalue number $ k $ as returned by eigenvalues().
569 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
570 const auto & eigenVectors = solver.eigenvectors();
571
572 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
573 {
574 auto copyEigenValues = eigenValues;
575 auto copyEigenVectors = eigenVectors;
576 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
577 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
578
579 for (unsigned int row = 0; row < m_Dimension; ++row)
580 {
581 EigenValues[row] = copyEigenValues[row];
582 for (unsigned int col = 0; col < m_Dimension; ++col)
583 {
584 EigenVectors[row][col] = copyEigenVectors(col, row);
585 }
586 }
587 }
588 else
589 {
590 for (unsigned int row = 0; row < m_Dimension; ++row)
591 {
592 EigenValues[row] = eigenValues[row];
593 for (unsigned int col = 0; col < m_Dimension; ++col)
594 {
595 EigenVectors[row][col] = eigenVectors(col, row);
596 }
597 }
598 }
599 // No error code
600 return 1;
601 }
602
603
604 /* Implementation detail using EigenLib that do not perform a copy.
605 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
606 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
607 * should be included.
608 *
609 * @param (bool) implementation detail argument making this implementation the most favourable
610 * to be chosen from all the alternative implementations.
611 *
612 * @return an unsigned int with no information value (no error code in EigenLib) */
613 template <typename QMatrix>
614 auto
616 TVector & EigenValues,
617 TEigenMatrix & EigenVectors,
618 bool) const -> decltype(GetPointerToMatrixData(A), 1U)
619 {
620 auto pointerToData = GetPointerToMatrixData(A);
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); // Computes EigenValues and EigenVectors
629 const auto & eigenValues = solver.eigenvalues();
630 /* Column k of the returned matrix is an eigenvector corresponding to
631 * eigenvalue number $ k $ as returned by eigenvalues().
632 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
633 const auto & eigenVectors = solver.eigenvectors();
634 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
635 {
636 auto copyEigenValues = eigenValues;
637 auto copyEigenVectors = eigenVectors;
638 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
639 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
640 for (unsigned int row = 0; row < m_Dimension; ++row)
641 {
642 EigenValues[row] = copyEigenValues[row];
643 for (unsigned int col = 0; col < m_Dimension; ++col)
644 {
645 EigenVectors[row][col] = copyEigenVectors(col, row);
646 }
647 }
648 }
649 else
650 {
651 for (unsigned int row = 0; row < m_Dimension; ++row)
652 {
653 EigenValues[row] = eigenValues[row];
654 for (unsigned int col = 0; col < m_Dimension; ++col)
655 {
656 EigenVectors[row][col] = eigenVectors(col, row);
657 }
658 }
659 }
660 // No error code
661 return 1;
662 }
663
664 /* Wrapper that call the right implementation for the type of matrix. */
665 unsigned int
666 ComputeEigenValuesWithEigenLibrary(const TMatrix & A, TVector & EigenValues) const
667 {
668 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
669 }
670
671 /* Implementation detail using EigenLib that performs a copy of the input matrix.
672 *
673 * @param (long) implementation detail argument making this implementation less favourable
674 * to be chosen if alternatives are available.
675 *
676 * @return an unsigned int with no information value (no error code in EigenLib) */
677 template <typename QMatrix>
678 auto
679 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const -> decltype(1U)
680 {
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)
685 {
686 for (unsigned int col = 0; col < m_Dimension; ++col)
687 {
688 inputMatrix(row, col) = A(row, col);
689 }
690 }
691 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
692 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
693 auto eigenValues = solver.eigenvalues();
694 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
695 {
696 detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
697 }
698 for (unsigned int i = 0; i < m_Dimension; ++i)
699 {
700 EigenValues[i] = eigenValues[i];
701 }
702
703 // No error code
704 return 1;
705 }
706
707 /* Implementation detail using EigenLib that do not perform a copy.
708 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
709 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
710 * should be included.
711 *
712 * @param (bool) implementation detail argument making this implementation the most favourable
713 * to be chosen from all the alternative implementations.
714 *
715 * @return an unsigned int with no information value (no error code in EigenLib) */
716 template <typename QMatrix>
717 auto
718 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
719 -> decltype(GetPointerToMatrixData(A), 1U)
720 {
721 auto pointerToData = GetPointerToMatrixData(A);
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();
731 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
732 {
733 detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
734 }
735 for (unsigned int i = 0; i < m_Dimension; ++i)
736 {
737 EigenValues[i] = eigenValues[i];
738 }
739 // No error code
740 return 1;
741 }
742};
743
744template <typename TMatrix, typename TVector, typename TEigenMatrix>
745std::ostream &
747{
748 os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
749 os << " Dimension : " << s.GetDimension() << std::endl;
750 os << " Order : " << s.GetOrder() << std::endl;
751 os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
752 os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
753 os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
754 return os;
755}
756
757template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
759{
760public:
761#if !defined(ITK_LEGACY_REMOVE)
763 using EigenValueOrderType = EigenValueOrderEnum;
764#endif
765#if !defined(ITK_LEGACY_REMOVE)
766 // We need to expose the enum values at the class level
767 // for backwards compatibility
768 static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
769 static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
770 static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
771#endif
772
773 using MatrixType = TMatrix;
774 using EigenMatrixType = TEigenMatrix;
775 using VectorType = TVector;
776
790 unsigned int
791 ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const
792 {
793 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
794 }
795
816 unsigned int
817 ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const
818 {
819 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
820 }
821
822 void
824 {
825 if (b)
826 {
827 m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
828 }
829 else
830 {
831 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
832 }
833 }
834 bool
836 {
837 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
838 }
839 void
841 {
842 if (b)
843 {
844 m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
845 }
846 else
847 {
848 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
849 }
850 }
851 bool
853 {
854 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
855 }
856 constexpr unsigned int
857 GetOrder() const
858 {
859 return VDimension;
860 }
861 constexpr unsigned int
863 {
864 return VDimension;
865 }
866 constexpr bool
868 {
869 return true;
870 }
871
872private:
874
875 /* Helper to get the matrix value type for EigenLibMatrix typename.
876 *
877 * If the TMatrix is vnl, the type is in element_type.
878 * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
879 *
880 * To use this function:
881 * using ValueType = decltype(this->GetMatrixType(true));
882 */
883 template <typename QMatrix = TMatrix>
884 auto
885 GetMatrixValueType(bool) const -> typename QMatrix::element_type
886 {
887 return QMatrix::element_type();
888 }
889 template <typename QMatrix = TMatrix>
890 auto
891 GetMatrixValueType(bool) const -> typename QMatrix::ValueType
892 {
893 return QMatrix::ValueType();
894 }
895
896 /* Implementation detail using EigenLib that do not perform a copy.
897 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
898 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
899 * should be included.
900 *
901 * @param (bool) implementation detail argument making this implementation the most favourable
902 * to be chosen from all the alternative implementations.
903 *
904 * @return an unsigned int with no information value (no error code in EigenLib) */
905 template <typename QMatrix>
906 auto
908 TVector & EigenValues,
909 TEigenMatrix & EigenVectors,
910 bool) const -> decltype(GetPointerToMatrixData(A), 1U)
911 {
912 auto pointerToData = GetPointerToMatrixData(A);
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); // Computes EigenValues and EigenVectors
921 const auto & eigenValues = solver.eigenvalues();
922 /* Column k of the returned matrix is an eigenvector corresponding to
923 * eigenvalue number $ k $ as returned by eigenvalues().
924 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
925 const auto & eigenVectors = solver.eigenvectors();
926 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
927 {
928 auto copyEigenValues = eigenValues;
929 auto copyEigenVectors = eigenVectors;
930 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
931 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
932 for (unsigned int row = 0; row < VDimension; ++row)
933 {
934 EigenValues[row] = copyEigenValues[row];
935 for (unsigned int col = 0; col < VDimension; ++col)
936 {
937 EigenVectors[row][col] = copyEigenVectors(col, row);
938 }
939 }
940 }
941 else
942 {
943 for (unsigned int row = 0; row < VDimension; ++row)
944 {
945 EigenValues[row] = eigenValues[row];
946 for (unsigned int col = 0; col < VDimension; ++col)
947 {
948 EigenVectors[row][col] = eigenVectors(col, row);
949 }
950 }
951 }
952 // No error code
953 return 1;
954 }
955
956 /* Implementation detail using EigenLib that performs a copy of the input matrix.
957 *
958 * @param (long) implementation detail argument making this implementation less favourable
959 * to be chosen if alternatives are available.
960 *
961 * @return an unsigned int with no information value (no error code in EigenLib) */
962 template <typename QMatrix>
963 auto
965 TVector & EigenValues,
966 TEigenMatrix & EigenVectors,
967 long) const -> decltype(1U)
968 {
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)
973 {
974 for (unsigned int col = 0; col < VDimension; ++col)
975 {
976 inputMatrix(row, col) = A(row, col);
977 }
978 }
979 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
980 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
981 const auto & eigenValues = solver.eigenvalues();
982 /* Column k of the returned matrix is an eigenvector corresponding to
983 * eigenvalue number $ k $ as returned by eigenvalues().
984 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
985 const auto & eigenVectors = solver.eigenvectors();
986
987 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
988 {
989 auto copyEigenValues = eigenValues;
990 auto copyEigenVectors = eigenVectors;
991 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
992 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
993
994 for (unsigned int row = 0; row < VDimension; ++row)
995 {
996 EigenValues[row] = copyEigenValues[row];
997 for (unsigned int col = 0; col < VDimension; ++col)
998 {
999 EigenVectors[row][col] = copyEigenVectors(col, row);
1000 }
1001 }
1002 }
1003 else
1004 {
1005 for (unsigned int row = 0; row < VDimension; ++row)
1006 {
1007 EigenValues[row] = eigenValues[row];
1008 for (unsigned int col = 0; col < VDimension; ++col)
1009 {
1010 EigenVectors[row][col] = eigenVectors(col, row);
1011 }
1012 }
1013 }
1014 // No error code
1015 return 1;
1016 }
1017
1018 /* Implementation detail using EigenLib that performs a copy of the input matrix.
1019 *
1020 * @param (long) implementation detail argument making this implementation less favourable
1021 * to be chosen if alternatives are available.
1022 *
1023 * @return an unsigned int with no information value (no error code in EigenLib) */
1024 template <typename QMatrix>
1025 auto
1026 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const -> decltype(1U)
1027 {
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)
1032 {
1033 for (unsigned int col = 0; col < VDimension; ++col)
1034 {
1035 inputMatrix(row, col) = A(row, col);
1036 }
1037 }
1038 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1039 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1040 auto eigenValues = solver.eigenvalues();
1041 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1042 {
1043 detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1044 }
1045 for (unsigned int i = 0; i < VDimension; ++i)
1046 {
1047 EigenValues[i] = eigenValues[i];
1048 }
1049
1050 // No error code
1051 return 1;
1052 }
1053
1054 /* Implementation detail using EigenLib that do not perform a copy.
1055 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
1056 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
1057 * should be included.
1058 *
1059 * @param (bool) implementation detail argument making this implementation the most favourable
1060 * to be chosen from all the alternative implementations.
1061 *
1062 * @return an unsigned int with no information value (no error code in EigenLib) */
1063 template <typename QMatrix>
1064 auto
1065 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
1066 -> decltype(GetPointerToMatrixData(A), 1U)
1067 {
1068 auto pointerToData = GetPointerToMatrixData(A);
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();
1078 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1079 {
1080 detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1081 }
1082 for (unsigned int i = 0; i < VDimension; ++i)
1083 {
1084 EigenValues[i] = eigenValues[i];
1085 }
1086 // No error code
1087 return 1;
1088 }
1089};
1090
1091template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix>
1092std::ostream &
1093operator<<(std::ostream & os,
1095{
1096 os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1097 os << " Dimension : " << s.GetDimension() << std::endl;
1098 os << " Order : " << s.GetOrder() << std::endl;
1099 os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
1100 os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
1101 os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
1102 return os;
1103}
1104} // end namespace itk
1105
1106#ifndef ITK_MANUAL_INSTANTIATION
1107# include "itkSymmetricEigenAnalysis.hxx"
1108#endif
1109
1110#endif
A templated class holding a M x N size Matrix.
Definition: itkMatrix.h:53
InternalMatrixType & GetVnlMatrix()
Definition: itkMatrix.h:197
This class contains all enum classes used by SymmetricEigenAnalysis class.
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(1U)
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), 1U)
unsigned int ComputeEigenValues(const TMatrix &A, TVector &EigenValues) const
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 SetDimension(const unsigned int n)
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) 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 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)
void SetOrder(const unsigned int n)
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
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
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
bool abs(bool x)
Definition: itkMath.h:840
static constexpr double e
Definition: itkMath.h:56
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)