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);
78
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
115
116} // end namespace detail
117
123{
124public:
132 enum class EigenValueOrder : uint8_t
133 {
137 };
138};
139// Define how to print enumeration
140extern ITKCommon_EXPORT std::ostream &
141 operator<<(std::ostream & out, const SymmetricEigenAnalysisEnums::EigenValueOrder value);
142
144
146Int2EigenValueOrderEnum(const uint8_t value)
147{
148 switch (value)
149 {
150 case 1:
151 return EigenValueOrderEnum::OrderByValue;
152 case 2:
153 return EigenValueOrderEnum::OrderByMagnitude;
154 case 3:
155 return EigenValueOrderEnum::DoNotOrder;
156 default:
157 break;
158 }
159 itkGenericExceptionMacro("Error: Invalid value for conversion.");
160}
161
162#if !defined(ITK_LEGACY_REMOVE)
164static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
165static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
166static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
167#endif
168
202
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 }
297
298
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 }
320
321
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 {
337 }
338 }
339
340
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 }
370
371
372private:
373 bool m_UseEigenLibrary{ false };
374 unsigned int m_Dimension{ 0 };
375 unsigned int m_Order{ 0 };
376 EigenValueOrderEnum m_OrderEigenValues{ EigenValueOrderEnum::OrderByValue };
377
396 void
397 ReduceToTridiagonalMatrix(double * a, double * d, double * e, double * e2) const;
398
419 void
420 ReduceToTridiagonalMatrixAndGetTransformation(const double * a, double * d, double * e, double * z) const;
421
449 unsigned int
450 ComputeEigenValuesUsingQL(double * d, double * e) const;
451
487 unsigned int
488 ComputeEigenValuesAndVectorsUsingQL(double * d, double * e, double * z) const;
489
490 /* Legacy algorithms using thread-safe netlib.
491 * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
492 */
493 unsigned int
494 ComputeEigenValuesLegacy(const TMatrix & A, TVector & D) const;
495
496 unsigned int
497 ComputeEigenValuesAndVectorsLegacy(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
498
499 /* Helper to get the matrix value type for EigenLibMatrix typename.
500 *
501 * If the TMatrix is vnl, the type is in element_type.
502 * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
503 *
504 * To use this function:
505 * using ValueType = decltype(this->GetMatrixType(true));
506 *
507 * \note The two `GetMatrixValueType` overloads have different
508 * parameter declarations (`bool` and `...`), to avoid that both
509 * functions are equally good candidates during overload resolution,
510 * in case `element_type` and `ValueType` are both nested types of
511 * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
512 */
513 template <typename QMatrix = TMatrix>
514 auto
515 GetMatrixValueType(bool) const -> typename QMatrix::element_type
516 {
517 return QMatrix::element_type();
518 }
519 template <typename QMatrix = TMatrix>
520 auto
521 GetMatrixValueType(...) const -> typename QMatrix::ValueType
522 {
523 return QMatrix::ValueType();
524 }
525
526 /* Wrapper that call the right implementation for the type of matrix. */
527 unsigned int
529 TVector & EigenValues,
530 TEigenMatrix & EigenVectors) const
531 {
532 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
533 }
534
535 /* Implementation detail using EigenLib that performs a copy of the input matrix.
536 *
537 * @param (long) implementation detail argument making this implementation less favourable
538 * to be chosen if alternatives are available.
539 *
540 * @return an unsigned int with no information value (no error code in EigenLib) */
541 template <typename QMatrix>
542 auto
544 TVector & EigenValues,
545 TEigenMatrix & EigenVectors,
546 long) const -> decltype(1U)
547 {
548 using ValueType = decltype(GetMatrixValueType(true));
549 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
550 EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
551 for (unsigned int row = 0; row < m_Dimension; ++row)
552 {
553 for (unsigned int col = 0; col < m_Dimension; ++col)
554 {
555 inputMatrix(row, col) = A(row, col);
556 }
557 }
558 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
559 const EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
560 const auto & eigenValues = solver.eigenvalues();
561 /* Column k of the returned matrix is an eigenvector corresponding to
562 * eigenvalue number $ k $ as returned by eigenvalues().
563 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
564 const auto & eigenVectors = solver.eigenvectors();
565
566 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
567 {
568 auto copyEigenValues = eigenValues;
569 auto copyEigenVectors = eigenVectors;
570 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
571 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
572
573 for (unsigned int row = 0; row < m_Dimension; ++row)
574 {
575 EigenValues[row] = copyEigenValues[row];
576 for (unsigned int col = 0; col < m_Dimension; ++col)
577 {
578 EigenVectors[row][col] = copyEigenVectors(col, row);
579 }
580 }
581 }
582 else
583 {
584 for (unsigned int row = 0; row < m_Dimension; ++row)
585 {
586 EigenValues[row] = eigenValues[row];
587 for (unsigned int col = 0; col < m_Dimension; ++col)
588 {
589 EigenVectors[row][col] = eigenVectors(col, row);
590 }
591 }
592 }
593 // No error code
594 return 1;
595 }
596
597
598 /* Implementation detail using EigenLib that do not perform a copy.
599 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
600 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
601 * should be included.
602 *
603 * @param (bool) implementation detail argument making this implementation the most favourable
604 * to be chosen from all the alternative implementations.
605 *
606 * @return an unsigned int with no information value (no error code in EigenLib) */
607 template <typename QMatrix>
608 auto
610 TVector & EigenValues,
611 TEigenMatrix & EigenVectors,
612 bool) const -> decltype(GetPointerToMatrixData(A), 1U)
613 {
614 auto pointerToData = GetPointerToMatrixData(A);
615 using PointerType = decltype(pointerToData);
616 using ValueTypeCV = std::remove_pointer_t<PointerType>;
617 using ValueType = std::remove_cv_t<ValueTypeCV>;
618 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
619 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
620 EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
621 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
622 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
623 const auto & eigenValues = solver.eigenvalues();
624 /* Column k of the returned matrix is an eigenvector corresponding to
625 * eigenvalue number $ k $ as returned by eigenvalues().
626 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
627 const auto & eigenVectors = solver.eigenvectors();
628 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
629 {
630 auto copyEigenValues = eigenValues;
631 auto copyEigenVectors = eigenVectors;
632 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
633 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
634 for (unsigned int row = 0; row < m_Dimension; ++row)
635 {
636 EigenValues[row] = copyEigenValues[row];
637 for (unsigned int col = 0; col < m_Dimension; ++col)
638 {
639 EigenVectors[row][col] = copyEigenVectors(col, row);
640 }
641 }
642 }
643 else
644 {
645 for (unsigned int row = 0; row < m_Dimension; ++row)
646 {
647 EigenValues[row] = eigenValues[row];
648 for (unsigned int col = 0; col < m_Dimension; ++col)
649 {
650 EigenVectors[row][col] = eigenVectors(col, row);
651 }
652 }
653 }
654 // No error code
655 return 1;
656 }
657
658 /* Wrapper that call the right implementation for the type of matrix. */
659 unsigned int
660 ComputeEigenValuesWithEigenLibrary(const TMatrix & A, TVector & EigenValues) const
661 {
662 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
663 }
664
665 /* Implementation detail using EigenLib that performs a copy of the input matrix.
666 *
667 * @param (long) implementation detail argument making this implementation less favourable
668 * to be chosen if alternatives are available.
669 *
670 * @return an unsigned int with no information value (no error code in EigenLib) */
671 template <typename QMatrix>
672 auto
673 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const -> decltype(1U)
674 {
675 using ValueType = decltype(GetMatrixValueType(true));
676 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
677 EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
678 for (unsigned int row = 0; row < m_Dimension; ++row)
679 {
680 for (unsigned int col = 0; col < m_Dimension; ++col)
681 {
682 inputMatrix(row, col) = A(row, col);
683 }
684 }
685 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
686 const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
687 auto eigenValues = solver.eigenvalues();
688 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
689 {
691 }
692 for (unsigned int i = 0; i < m_Dimension; ++i)
693 {
694 EigenValues[i] = eigenValues[i];
695 }
696
697 // No error code
698 return 1;
699 }
700
701 /* Implementation detail using EigenLib that do not perform a copy.
702 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
703 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
704 * should be included.
705 *
706 * @param (bool) implementation detail argument making this implementation the most favourable
707 * to be chosen from all the alternative implementations.
708 *
709 * @return an unsigned int with no information value (no error code in EigenLib) */
710 template <typename QMatrix>
711 auto
712 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
713 -> decltype(GetPointerToMatrixData(A), 1U)
714 {
715 auto pointerToData = GetPointerToMatrixData(A);
716 using PointerType = decltype(pointerToData);
717 using ValueTypeCV = std::remove_pointer_t<PointerType>;
718 using ValueType = std::remove_cv_t<ValueTypeCV>;
719 using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
720 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
721 EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
722 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
723 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
724 auto eigenValues = solver.eigenvalues();
725 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
726 {
728 }
729 for (unsigned int i = 0; i < m_Dimension; ++i)
730 {
731 EigenValues[i] = eigenValues[i];
732 }
733 // No error code
734 return 1;
735 }
736};
737
738template <typename TMatrix, typename TVector, typename TEigenMatrix>
739std::ostream &
741{
742 os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
743 os << " Dimension : " << s.GetDimension() << std::endl;
744 os << " Order : " << s.GetOrder() << std::endl;
745 os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
746 os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
747 os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
748 return os;
749}
750
751template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
753{
754public:
755#if !defined(ITK_LEGACY_REMOVE)
757 using EigenValueOrderType = EigenValueOrderEnum;
758#endif
759#if !defined(ITK_LEGACY_REMOVE)
760 // We need to expose the enum values at the class level
761 // for backwards compatibility
762 static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
763 static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
764 static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
765#endif
766
767 using MatrixType = TMatrix;
768 using EigenMatrixType = TEigenMatrix;
769 using VectorType = TVector;
770
784 unsigned int
785 ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const
786 {
787 return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
788 }
789
810 unsigned int
811 ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const
812 {
813 return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
814 }
815
816 void
818 {
819 if (b)
820 {
821 m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
822 }
823 else
824 {
825 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
826 }
827 }
828 bool
830 {
831 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
832 }
833 void
835 {
836 if (b)
837 {
838 m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
839 }
840 else
841 {
842 m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
843 }
844 }
845 bool
847 {
848 return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
849 }
850 constexpr unsigned int
851 GetOrder() const
852 {
853 return VDimension;
854 }
855 constexpr unsigned int
857 {
858 return VDimension;
859 }
860 constexpr bool
862 {
863 return true;
864 }
865
866private:
867 EigenValueOrderEnum m_OrderEigenValues{ EigenValueOrderEnum::OrderByValue };
868
869 /* Helper to get the matrix value type for EigenLibMatrix typename.
870 *
871 * If the TMatrix is vnl, the type is in element_type.
872 * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
873 *
874 * To use this function:
875 * using ValueType = decltype(this->GetMatrixType(true));
876 */
877 template <typename QMatrix = TMatrix>
878 auto
879 GetMatrixValueType(bool) const -> typename QMatrix::element_type
880 {
881 return QMatrix::element_type();
882 }
883 template <typename QMatrix = TMatrix>
884 auto
885 GetMatrixValueType(bool) const -> typename QMatrix::ValueType
886 {
887 return QMatrix::ValueType();
888 }
889
890 /* Implementation detail using EigenLib that do not perform a copy.
891 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
892 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
893 * should be included.
894 *
895 * @param (bool) implementation detail argument making this implementation the most favourable
896 * to be chosen from all the alternative implementations.
897 *
898 * @return an unsigned int with no information value (no error code in EigenLib) */
899 template <typename QMatrix>
900 auto
902 TVector & EigenValues,
903 TEigenMatrix & EigenVectors,
904 bool) const -> decltype(GetPointerToMatrixData(A), 1U)
905 {
906 auto pointerToData = GetPointerToMatrixData(A);
907 using PointerType = decltype(pointerToData);
908 using ValueTypeCV = std::remove_pointer_t<PointerType>;
909 using ValueType = std::remove_cv_t<ValueTypeCV>;
910 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
911 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
912 EigenConstMatrixMap inputMatrix(pointerToData);
913 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
914 EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
915 const auto & eigenValues = solver.eigenvalues();
916 /* Column k of the returned matrix is an eigenvector corresponding to
917 * eigenvalue number $ k $ as returned by eigenvalues().
918 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
919 const auto & eigenVectors = solver.eigenvectors();
920 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
921 {
922 auto copyEigenValues = eigenValues;
923 auto copyEigenVectors = eigenVectors;
924 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
925 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
926 for (unsigned int row = 0; row < VDimension; ++row)
927 {
928 EigenValues[row] = copyEigenValues[row];
929 for (unsigned int col = 0; col < VDimension; ++col)
930 {
931 EigenVectors[row][col] = copyEigenVectors(col, row);
932 }
933 }
934 }
935 else
936 {
937 for (unsigned int row = 0; row < VDimension; ++row)
938 {
939 EigenValues[row] = eigenValues[row];
940 for (unsigned int col = 0; col < VDimension; ++col)
941 {
942 EigenVectors[row][col] = eigenVectors(col, row);
943 }
944 }
945 }
946 // No error code
947 return 1;
948 }
949
950 /* Implementation detail using EigenLib that performs a copy of the input matrix.
951 *
952 * @param (long) implementation detail argument making this implementation less favourable
953 * to be chosen if alternatives are available.
954 *
955 * @return an unsigned int with no information value (no error code in EigenLib) */
956 template <typename QMatrix>
957 auto
959 TVector & EigenValues,
960 TEigenMatrix & EigenVectors,
961 long) const -> decltype(1U)
962 {
963 using ValueType = decltype(GetMatrixValueType(true));
964 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
965 EigenLibMatrixType inputMatrix;
966 for (unsigned int row = 0; row < VDimension; ++row)
967 {
968 for (unsigned int col = 0; col < VDimension; ++col)
969 {
970 inputMatrix(row, col) = A(row, col);
971 }
972 }
973 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
974 const EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
975 const auto & eigenValues = solver.eigenvalues();
976 /* Column k of the returned matrix is an eigenvector corresponding to
977 * eigenvalue number $ k $ as returned by eigenvalues().
978 * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
979 const auto & eigenVectors = solver.eigenvectors();
980
981 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
982 {
983 auto copyEigenValues = eigenValues;
984 auto copyEigenVectors = eigenVectors;
985 auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
986 detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
987
988 for (unsigned int row = 0; row < VDimension; ++row)
989 {
990 EigenValues[row] = copyEigenValues[row];
991 for (unsigned int col = 0; col < VDimension; ++col)
992 {
993 EigenVectors[row][col] = copyEigenVectors(col, row);
994 }
995 }
996 }
997 else
998 {
999 for (unsigned int row = 0; row < VDimension; ++row)
1000 {
1001 EigenValues[row] = eigenValues[row];
1002 for (unsigned int col = 0; col < VDimension; ++col)
1003 {
1004 EigenVectors[row][col] = eigenVectors(col, row);
1005 }
1006 }
1007 }
1008 // No error code
1009 return 1;
1010 }
1011
1012 /* Implementation detail using EigenLib that performs a copy of the input matrix.
1013 *
1014 * @param (long) implementation detail argument making this implementation less favourable
1015 * to be chosen if alternatives are available.
1016 *
1017 * @return an unsigned int with no information value (no error code in EigenLib) */
1018 template <typename QMatrix>
1019 auto
1020 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const -> decltype(1U)
1021 {
1022 using ValueType = decltype(GetMatrixValueType(true));
1023 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1024 EigenLibMatrixType inputMatrix;
1025 for (unsigned int row = 0; row < VDimension; ++row)
1026 {
1027 for (unsigned int col = 0; col < VDimension; ++col)
1028 {
1029 inputMatrix(row, col) = A(row, col);
1030 }
1031 }
1032 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1033 const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1034 auto eigenValues = solver.eigenvalues();
1035 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1036 {
1037 detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1038 }
1039 for (unsigned int i = 0; i < VDimension; ++i)
1040 {
1041 EigenValues[i] = eigenValues[i];
1042 }
1043
1044 // No error code
1045 return 1;
1046 }
1047
1048 /* Implementation detail using EigenLib that do not perform a copy.
1049 * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
1050 * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
1051 * should be included.
1052 *
1053 * @param (bool) implementation detail argument making this implementation the most favourable
1054 * to be chosen from all the alternative implementations.
1055 *
1056 * @return an unsigned int with no information value (no error code in EigenLib) */
1057 template <typename QMatrix>
1058 auto
1059 ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
1060 -> decltype(GetPointerToMatrixData(A), 1U)
1061 {
1062 auto pointerToData = GetPointerToMatrixData(A);
1063 using PointerType = decltype(pointerToData);
1064 using ValueTypeCV = std::remove_pointer_t<PointerType>;
1065 using ValueType = std::remove_cv_t<ValueTypeCV>;
1066 using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1067 using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1068 EigenConstMatrixMap inputMatrix(pointerToData);
1069 using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1070 EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1071 auto eigenValues = solver.eigenvalues();
1072 if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
1073 {
1074 detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1075 }
1076 for (unsigned int i = 0; i < VDimension; ++i)
1077 {
1078 EigenValues[i] = eigenValues[i];
1079 }
1080 // No error code
1081 return 1;
1082 }
1083};
1084
1085template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix>
1086std::ostream &
1087operator<<(std::ostream & os,
1089{
1090 os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1091 os << " Dimension : " << s.GetDimension() << std::endl;
1092 os << " Order : " << s.GetOrder() << std::endl;
1093 os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
1094 os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
1095 os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
1096 return os;
1097}
1098} // end namespace itk
1099
1100#ifndef ITK_MANUAL_INSTANTIATION
1101# include "itkSymmetricEigenAnalysis.hxx"
1102#endif
1103
1104#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:837
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)