ITK 6.0.0
Insight Toolkit
 
Loading...
Searching...
No Matches
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);
77
78 std::sort(std::begin(indicesSortPermutations),
79 std::end(indicesSortPermutations),
80 [&eigenValues](unsigned int a, unsigned int b) {
81 return itk::Math::abs(eigenValues[a]) < itk::Math::abs(eigenValues[b]);
82 });
83 auto tmpCopy = eigenValues;
84 for (unsigned int i = 0; i < numberOfElements; ++i)
85 {
86 eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
87 }
88 return indicesSortPermutations;
89}
90
99template <typename QMatrix>
100void
101permuteColumnsWithSortIndices(QMatrix & eigenVectors, const std::vector<int> & indicesSortPermutations)
102{
103 using EigenLibPermutationMatrix = Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
104 auto numberOfElements = indicesSortPermutations.size();
105 // Creates a NxN permutation matrix copying our permutation to the matrix indices.
106 // Which holds the 1D array representation of a permutation.
107 EigenLibPermutationMatrix perm(numberOfElements);
108 perm.setIdentity();
109 std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(), perm.indices().data());
110 // Apply it
111 eigenVectors = eigenVectors * perm;
112}
113
114} // end namespace detail
115
121{
122public:
130 enum class EigenValueOrder : uint8_t
131 {
135 };
136};
137// Define how to print enumeration
138extern ITKCommon_EXPORT std::ostream &
139 operator<<(std::ostream & out, const SymmetricEigenAnalysisEnums::EigenValueOrder value);
140
142
144Int2EigenValueOrderEnum(const uint8_t value)
145{
146 switch (value)
147 {
148 case 1:
149 return EigenValueOrderEnum::OrderByValue;
150 case 2:
151 return EigenValueOrderEnum::OrderByMagnitude;
152 case 3:
153 return EigenValueOrderEnum::DoNotOrder;
154 default:
155 break;
156 }
157 itkGenericExceptionMacro("Error: Invalid value for conversion.");
158}
159
160#if !defined(ITK_LEGACY_REMOVE)
162static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
163static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
164static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
165#endif
166
200
201template <typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
202class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
203{
204public:
206#if !defined(ITK_LEGACY_REMOVE)
208 using EigenValueOrderType = EigenValueOrderEnum;
209#endif
210
212
213 SymmetricEigenAnalysis(const unsigned int dimension)
214 : m_Dimension(dimension)
215 , m_Order(dimension)
216 {}
217
219
220 using MatrixType = TMatrix;
221 using EigenMatrixType = TEigenMatrix;
222 using VectorType = TVector;
223
237 unsigned int
238 ComputeEigenValues(const TMatrix & A, TVector & D) const;
239
260 unsigned int
261 ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
262
263
265 void
266 SetOrder(const unsigned int n)
267 {
268 m_Order = n;
269 }
270
274 unsigned int
275 GetOrder() const
276 {
277 return m_Order;
278 }
279
284 void
286 {
287 if (b)
288 {
289 m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
290 }
291 else
292 {
293 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
294 }
295 }
296
297 bool
299 {
300 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
301 }
302
303
308 void
310 {
311 if (b)
312 {
313 m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
314 }
315 else
316 {
317 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
318 }
319 }
320
321 bool
323 {
324 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
325 }
326
327
330 void
331 SetDimension(const unsigned int n)
332 {
333 m_Dimension = n;
334 if (m_Order == 0)
335 {
337 }
338 }
339
342 unsigned int
344 {
345 return m_Dimension;
346 }
347
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 }
370
371private:
372 bool m_UseEigenLibrary{ false };
373 unsigned int m_Dimension{ 0 };
374 unsigned int m_Order{ 0 };
375 EigenValueOrderEnum m_OrderEigenValues{ EigenValueOrderEnum::OrderByValue };
376
395 void
396 ReduceToTridiagonalMatrix(double * a, double * d, double * e, double * e2) const;
397
418 void
419 ReduceToTridiagonalMatrixAndGetTransformation(const double * a, double * d, double * e, double * z) const;
420
448 unsigned int
449 ComputeEigenValuesUsingQL(double * d, double * e) const;
450
486 unsigned int
487 ComputeEigenValuesAndVectorsUsingQL(double * d, double * e, double * z) const;
488
489 /* Legacy algorithms using thread-safe netlib.
490 * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
491 */
492 unsigned int
493 ComputeEigenValuesLegacy(const TMatrix & A, TVector & D) const;
494
495 unsigned int
496 ComputeEigenValuesAndVectorsLegacy(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
497
498 /* Helper to get the matrix value type for EigenLibMatrix typename.
499 *
500 * If the TMatrix is vnl, the type is in element_type.
501 * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
502 *
503 * To use this function:
504 * using ValueType = decltype(this->GetMatrixType(true));
505 *
506 * \note The two `GetMatrixValueType` overloads have different
507 * parameter declarations (`bool` and `...`), to avoid that both
508 * functions are equally good candidates during overload resolution,
509 * in case `element_type` and `ValueType` are both nested types of
510 * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
511 */
512 template <typename QMatrix = TMatrix>
513 auto
514 GetMatrixValueType(bool) const -> typename QMatrix::element_type
515 {
516 return QMatrix::element_type();
517 }
518 template <typename QMatrix = TMatrix>
519 auto
520 GetMatrixValueType(...) const -> typename QMatrix::ValueType
521 {
522 return QMatrix::ValueType();
523 }
524
525 /* Wrapper that call the right implementation for the type of matrix. */
526 unsigned int
528 TVector & EigenValues,
529 TEigenMatrix & EigenVectors) const
530 {
531 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
532 }
533
534 /* Implementation detail using EigenLib that performs a copy of the input matrix.
535 *
536 * @param (long) implementation detail argument making this implementation less favourable
537 * to be chosen if alternatives are available.
538 *
539 * @return an unsigned int with no information value (no error code in EigenLib) */
540 template <typename QMatrix>
541 auto
543 TVector & EigenValues,
544 TEigenMatrix & EigenVectors,
545 long) const -> decltype(1U)
546 {
547 using ValueType = decltype(GetMatrixValueType(true));
548 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
549 EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
550 for (unsigned int row = 0; row < m_Dimension; ++row)
551 {
552 for (unsigned int col = 0; col < m_Dimension; ++col)
553 {
554 inputMatrix(row, col) = A(row, col);
555 }
556 }
557 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
558 const EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
559 const auto & eigenValues = solver.eigenvalues();
560 /* Column k of the returned matrix is an eigenvector corresponding to
561 * eigenvalue number $ k $ as returned by eigenvalues().
562 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
563 const auto & eigenVectors = solver.eigenvectors();
564
565 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
566 {
567 auto copyEigenValues = eigenValues;
568 auto copyEigenVectors = eigenVectors;
569 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
570 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
571
572 for (unsigned int row = 0; row < m_Dimension; ++row)
573 {
574 EigenValues[row] = copyEigenValues[row];
575 for (unsigned int col = 0; col < m_Dimension; ++col)
576 {
577 EigenVectors[row][col] = copyEigenVectors(col, row);
578 }
579 }
580 }
581 else
582 {
583 for (unsigned int row = 0; row < m_Dimension; ++row)
584 {
585 EigenValues[row] = eigenValues[row];
586 for (unsigned int col = 0; col < m_Dimension; ++col)
587 {
588 EigenVectors[row][col] = eigenVectors(col, row);
589 }
590 }
591 }
592 // No error code
593 return 1;
594 }
595
596
597 /* Implementation detail using EigenLib that do not perform a copy.
598 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
599 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
600 * should be included.
601 *
602 * @param (bool) implementation detail argument making this implementation the most favourable
603 * to be chosen from all the alternative implementations.
604 *
605 * @return an unsigned int with no information value (no error code in EigenLib) */
606 template <typename QMatrix>
607 auto
609 TVector & EigenValues,
610 TEigenMatrix & EigenVectors,
611 bool) const -> decltype(GetPointerToMatrixData(A), 1U)
612 {
613 auto pointerToData = GetPointerToMatrixData(A);
614 using PointerType = decltype(pointerToData);
615 using ValueTypeCV = std::remove_pointer_t<PointerType>;
616 using ValueType = std::remove_cv_t<ValueTypeCV>;
617 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
618 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
619 EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
620 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
621 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
622 const auto & eigenValues = solver.eigenvalues();
623 /* Column k of the returned matrix is an eigenvector corresponding to
624 * eigenvalue number $ k $ as returned by eigenvalues().
625 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
626 const auto & eigenVectors = solver.eigenvectors();
627 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
628 {
629 auto copyEigenValues = eigenValues;
630 auto copyEigenVectors = eigenVectors;
631 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
632 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
633 for (unsigned int row = 0; row < m_Dimension; ++row)
634 {
635 EigenValues[row] = copyEigenValues[row];
636 for (unsigned int col = 0; col < m_Dimension; ++col)
637 {
638 EigenVectors[row][col] = copyEigenVectors(col, row);
639 }
640 }
641 }
642 else
643 {
644 for (unsigned int row = 0; row < m_Dimension; ++row)
645 {
646 EigenValues[row] = eigenValues[row];
647 for (unsigned int col = 0; col < m_Dimension; ++col)
648 {
649 EigenVectors[row][col] = eigenVectors(col, row);
650 }
651 }
652 }
653 // No error code
654 return 1;
655 }
656
657 /* Wrapper that call the right implementation for the type of matrix. */
658 unsigned int
659 ComputeEigenValuesWithEigenLibrary(const TMatrix & A, TVector & EigenValues) const
660 {
661 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
662 }
663
664 /* Implementation detail using EigenLib that performs a copy of the input matrix.
665 *
666 * @param (long) implementation detail argument making this implementation less favourable
667 * to be chosen if alternatives are available.
668 *
669 * @return an unsigned int with no information value (no error code in EigenLib) */
670 template <typename QMatrix>
671 auto
672 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const -> decltype(1U)
673 {
674 using ValueType = decltype(GetMatrixValueType(true));
675 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
676 EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
677 for (unsigned int row = 0; row < m_Dimension; ++row)
678 {
679 for (unsigned int col = 0; col < m_Dimension; ++col)
680 {
681 inputMatrix(row, col) = A(row, col);
682 }
683 }
684 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
685 const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
686 auto eigenValues = solver.eigenvalues();
687 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
688 {
690 }
691 for (unsigned int i = 0; i < m_Dimension; ++i)
692 {
693 EigenValues[i] = eigenValues[i];
694 }
695
696 // No error code
697 return 1;
698 }
699
700 /* Implementation detail using EigenLib that do not perform a copy.
701 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
702 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
703 * should be included.
704 *
705 * @param (bool) implementation detail argument making this implementation the most favourable
706 * to be chosen from all the alternative implementations.
707 *
708 * @return an unsigned int with no information value (no error code in EigenLib) */
709 template <typename QMatrix>
710 auto
711 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
712 -> decltype(GetPointerToMatrixData(A), 1U)
713 {
714 auto pointerToData = GetPointerToMatrixData(A);
715 using PointerType = decltype(pointerToData);
716 using ValueTypeCV = std::remove_pointer_t<PointerType>;
717 using ValueType = std::remove_cv_t<ValueTypeCV>;
718 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
719 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
720 EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
721 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
722 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
723 auto eigenValues = solver.eigenvalues();
724 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
725 {
727 }
728 for (unsigned int i = 0; i < m_Dimension; ++i)
729 {
730 EigenValues[i] = eigenValues[i];
731 }
732 // No error code
733 return 1;
734 }
735};
736
737template <typename TMatrix, typename TVector, typename TEigenMatrix>
738std::ostream &
740{
741 os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
742 os << " Dimension : " << s.GetDimension() << std::endl;
743 os << " Order : " << s.GetOrder() << std::endl;
744 os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
745 os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
746 os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
747 return os;
748}
749
750template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
752{
753public:
754#if !defined(ITK_LEGACY_REMOVE)
756 using EigenValueOrderType = EigenValueOrderEnum;
757#endif
758#if !defined(ITK_LEGACY_REMOVE)
759 // We need to expose the enum values at the class level
760 // for backwards compatibility
761 static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
762 static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
763 static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
764#endif
765
766 using MatrixType = TMatrix;
767 using EigenMatrixType = TEigenMatrix;
768 using VectorType = TVector;
769
783 unsigned int
784 ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const
785 {
786 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
787 }
788
809 unsigned int
810 ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const
811 {
812 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
813 }
814
815 void
817 {
818 if (b)
819 {
820 m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
821 }
822 else
823 {
824 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
825 }
826 }
827 bool
829 {
830 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
831 }
832 void
834 {
835 if (b)
836 {
837 m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
838 }
839 else
840 {
841 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
842 }
843 }
844 bool
846 {
847 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
848 }
849 constexpr unsigned int
850 GetOrder() const
851 {
852 return VDimension;
853 }
854 constexpr unsigned int
856 {
857 return VDimension;
858 }
859 constexpr bool
861 {
862 return true;
863 }
864
865private:
866 EigenValueOrderEnum m_OrderEigenValues{ EigenValueOrderEnum::OrderByValue };
867
868 /* Helper to get the matrix value type for EigenLibMatrix typename.
869 *
870 * If the TMatrix is vnl, the type is in element_type.
871 * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
872 *
873 * To use this function:
874 * using ValueType = decltype(this->GetMatrixType(true));
875 */
876 template <typename QMatrix = TMatrix>
877 auto
878 GetMatrixValueType(bool) const -> typename QMatrix::element_type
879 {
880 return QMatrix::element_type();
881 }
882 template <typename QMatrix = TMatrix>
883 auto
884 GetMatrixValueType(bool) const -> typename QMatrix::ValueType
885 {
886 return QMatrix::ValueType();
887 }
888
889 /* Implementation detail using EigenLib that do not perform a copy.
890 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
891 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
892 * should be included.
893 *
894 * @param (bool) implementation detail argument making this implementation the most favourable
895 * to be chosen from all the alternative implementations.
896 *
897 * @return an unsigned int with no information value (no error code in EigenLib) */
898 template <typename QMatrix>
899 auto
901 TVector & EigenValues,
902 TEigenMatrix & EigenVectors,
903 bool) const -> decltype(GetPointerToMatrixData(A), 1U)
904 {
905 auto pointerToData = GetPointerToMatrixData(A);
906 using PointerType = decltype(pointerToData);
907 using ValueTypeCV = std::remove_pointer_t<PointerType>;
908 using ValueType = std::remove_cv_t<ValueTypeCV>;
909 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
910 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
911 EigenConstMatrixMap inputMatrix(pointerToData);
912 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
913 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
914 const auto & eigenValues = solver.eigenvalues();
915 /* Column k of the returned matrix is an eigenvector corresponding to
916 * eigenvalue number $ k $ as returned by eigenvalues().
917 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
918 const auto & eigenVectors = solver.eigenvectors();
919 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
920 {
921 auto copyEigenValues = eigenValues;
922 auto copyEigenVectors = eigenVectors;
923 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
924 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
925 for (unsigned int row = 0; row < VDimension; ++row)
926 {
927 EigenValues[row] = copyEigenValues[row];
928 for (unsigned int col = 0; col < VDimension; ++col)
929 {
930 EigenVectors[row][col] = copyEigenVectors(col, row);
931 }
932 }
933 }
934 else
935 {
936 for (unsigned int row = 0; row < VDimension; ++row)
937 {
938 EigenValues[row] = eigenValues[row];
939 for (unsigned int col = 0; col < VDimension; ++col)
940 {
941 EigenVectors[row][col] = eigenVectors(col, row);
942 }
943 }
944 }
945 // No error code
946 return 1;
947 }
948
949 /* Implementation detail using EigenLib that performs a copy of the input matrix.
950 *
951 * @param (long) implementation detail argument making this implementation less favourable
952 * to be chosen if alternatives are available.
953 *
954 * @return an unsigned int with no information value (no error code in EigenLib) */
955 template <typename QMatrix>
956 auto
958 TVector & EigenValues,
959 TEigenMatrix & EigenVectors,
960 long) const -> decltype(1U)
961 {
962 using ValueType = decltype(GetMatrixValueType(true));
963 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
964 EigenLibMatrixType inputMatrix;
965 for (unsigned int row = 0; row < VDimension; ++row)
966 {
967 for (unsigned int col = 0; col < VDimension; ++col)
968 {
969 inputMatrix(row, col) = A(row, col);
970 }
971 }
972 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
973 const EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
974 const auto & eigenValues = solver.eigenvalues();
975 /* Column k of the returned matrix is an eigenvector corresponding to
976 * eigenvalue number $ k $ as returned by eigenvalues().
977 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
978 const auto & eigenVectors = solver.eigenvectors();
979
980 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
981 {
982 auto copyEigenValues = eigenValues;
983 auto copyEigenVectors = eigenVectors;
984 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
985 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
986
987 for (unsigned int row = 0; row < VDimension; ++row)
988 {
989 EigenValues[row] = copyEigenValues[row];
990 for (unsigned int col = 0; col < VDimension; ++col)
991 {
992 EigenVectors[row][col] = copyEigenVectors(col, row);
993 }
994 }
995 }
996 else
997 {
998 for (unsigned int row = 0; row < VDimension; ++row)
999 {
1000 EigenValues[row] = eigenValues[row];
1001 for (unsigned int col = 0; col < VDimension; ++col)
1002 {
1003 EigenVectors[row][col] = eigenVectors(col, row);
1004 }
1005 }
1006 }
1007 // No error code
1008 return 1;
1009 }
1010
1011 /* Implementation detail using EigenLib that performs a copy of the input matrix.
1012 *
1013 * @param (long) implementation detail argument making this implementation less favourable
1014 * to be chosen if alternatives are available.
1015 *
1016 * @return an unsigned int with no information value (no error code in EigenLib) */
1017 template <typename QMatrix>
1018 auto
1019 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const -> decltype(1U)
1020 {
1021 using ValueType = decltype(GetMatrixValueType(true));
1022 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1023 EigenLibMatrixType inputMatrix;
1024 for (unsigned int row = 0; row < VDimension; ++row)
1025 {
1026 for (unsigned int col = 0; col < VDimension; ++col)
1027 {
1028 inputMatrix(row, col) = A(row, col);
1029 }
1030 }
1031 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1032 const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1033 auto eigenValues = solver.eigenvalues();
1034 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1035 {
1036 detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1037 }
1038 for (unsigned int i = 0; i < VDimension; ++i)
1039 {
1040 EigenValues[i] = eigenValues[i];
1041 }
1042
1043 // No error code
1044 return 1;
1045 }
1046
1047 /* Implementation detail using EigenLib that do not perform a copy.
1048 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
1049 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
1050 * should be included.
1051 *
1052 * @param (bool) implementation detail argument making this implementation the most favourable
1053 * to be chosen from all the alternative implementations.
1054 *
1055 * @return an unsigned int with no information value (no error code in EigenLib) */
1056 template <typename QMatrix>
1057 auto
1058 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
1059 -> decltype(GetPointerToMatrixData(A), 1U)
1060 {
1061 auto pointerToData = GetPointerToMatrixData(A);
1062 using PointerType = decltype(pointerToData);
1063 using ValueTypeCV = std::remove_pointer_t<PointerType>;
1064 using ValueType = std::remove_cv_t<ValueTypeCV>;
1065 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1066 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1067 EigenConstMatrixMap inputMatrix(pointerToData);
1068 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1069 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1070 auto eigenValues = solver.eigenvalues();
1071 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1072 {
1073 detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1074 }
1075 for (unsigned int i = 0; i < VDimension; ++i)
1076 {
1077 EigenValues[i] = eigenValues[i];
1078 }
1079 // No error code
1080 return 1;
1081 }
1082};
1083
1084template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix>
1085std::ostream &
1086operator<<(std::ostream & os,
1088{
1089 os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1090 os << " Dimension : " << s.GetDimension() << std::endl;
1091 os << " Order : " << s.GetOrder() << std::endl;
1092 os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
1093 os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
1094 os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
1095 return os;
1096}
1097} // end namespace itk
1098
1099#ifndef ITK_MANUAL_INSTANTIATION
1100# include "itkSymmetricEigenAnalysis.hxx"
1101#endif
1102
1103#endif
A templated class holding a M x N size Matrix.
Definition itkMatrix.h:53
InternalMatrixType & GetVnlMatrix()
Definition itkMatrix.h:214
This class contains all enum classes used by SymmetricEigenAnalysis class.
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, EigenValuesArrayType &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, EigenValuesArrayType &EigenValues, EigenVectorsMatrixType &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
itk::EigenValueOrderEnum EigenValueOrderEnum
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, TOutput &EigenValues, TInput &EigenVectors, long) const -> decltype(1U)
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TOutput &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:836
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)