ITK  5.1.0
Insight Segmentation and Registration Toolkit
itkSymmetricEigenAnalysis.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
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  * http://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 
31 namespace itk
32 {
33 namespace detail
34 {
35 /* Helper functions returning pointer to matrix data for different types. */
36 template <typename TValueType, unsigned int NRows, unsigned int NCols>
37 const TValueType *
38 GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, NRows, NCols> & inputMatrix)
39 {
40  return inputMatrix.data_block();
41 };
42 template <typename TValueType>
43 const TValueType *
44 GetPointerToMatrixData(const vnl_matrix<TValueType> & inputMatrix)
45 {
46  return inputMatrix.data_block();
47 };
48 
49 template <typename TValueType, unsigned int NRows, unsigned int NCols>
50 const TValueType *
52 {
53  return inputMatrix.GetVnlMatrix().data_block();
54 };
55 
71 template <typename TArray>
72 std::vector<int>
73 sortEigenValuesByMagnitude(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(
80  std::begin(indicesSortPermutations),
81  std::end(indicesSortPermutations),
82  [&eigenValues](unsigned int a, unsigned int b) { return std::abs(eigenValues[a]) < std::abs(eigenValues[b]); });
83  auto tmpCopy = eigenValues;
84  for (unsigned int i = 0; i < numberOfElements; ++i)
85  {
86  eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
87  }
88  return indicesSortPermutations;
89 }
90 
99 template <typename QMatrix>
100 void
101 permuteColumnsWithSortIndices(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 } // end namespace detail
115 
120 enum class OrderType : uint8_t
121 {
122  OrderByValue = 1,
124  DoNotOrder
125 };
126 
162 template <typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
163 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
164 {
165 public:
168 #if !defined(ITK_LEGACY_REMOVE)
169  // We need to expose the enum values at the class level
170  // for backwards compatibility
174 #endif
175 
177  : m_OrderEigenValues(OrderType::OrderByValue)
178  {}
179 
180  SymmetricEigenAnalysis(const unsigned int dimension)
181  :
182 
183  m_Dimension(dimension)
184  , m_Order(dimension)
185  , m_OrderEigenValues(OrderType::OrderByValue)
186  {}
187 
188  ~SymmetricEigenAnalysis() = default;
189 
190  using MatrixType = TMatrix;
191  using EigenMatrixType = TEigenMatrix;
192  using VectorType = TVector;
193 
207  unsigned int
208  ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const;
209 
230  unsigned int
231  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
232 
233 
235  void
236  SetOrder(const unsigned int n)
237  {
238  m_Order = n;
239  }
240 
244  unsigned int
245  GetOrder() const
246  {
247  return m_Order;
248  }
249 
253  void
254  SetOrderEigenValues(const bool b)
255  {
256  if (b)
257  {
258  m_OrderEigenValues = OrderType::OrderByValue;
259  }
260  else
261  {
262  m_OrderEigenValues = OrderType::DoNotOrder;
263  }
264  }
266 
267  bool
269  {
270  return (m_OrderEigenValues == OrderType::OrderByValue);
271  }
272 
276  void
278  {
279  if (b)
280  {
281  m_OrderEigenValues = OrderType::OrderByMagnitude;
282  }
283  else
284  {
285  m_OrderEigenValues = OrderType::DoNotOrder;
286  }
287  }
289 
290  bool
292  {
293  return (m_OrderEigenValues == OrderType::OrderByMagnitude);
294  }
295 
298  void
299  SetDimension(const unsigned int n)
300  {
301  m_Dimension = n;
302  if (m_Order == 0)
303  {
304  m_Order = m_Dimension;
305  }
306  }
308 
311  unsigned int
312  GetDimension() const
313  {
314  return m_Dimension;
315  }
316 
318  void
319  SetUseEigenLibrary(const bool input)
320  {
321  m_UseEigenLibrary = input;
322  }
323  void
325  {
326  m_UseEigenLibrary = true;
327  }
328  void
330  {
331  m_UseEigenLibrary = false;
332  }
333  bool
335  {
336  return m_UseEigenLibrary;
337  }
339 
340 private:
341  bool m_UseEigenLibrary{ false };
342  unsigned int m_Dimension{ 0 };
343  unsigned int m_Order{ 0 };
345 
365  void
366  ReduceToTridiagonalMatrix(double * inputMatrix, double * d, double * e, double * e2) const;
367 
389  void
390  ReduceToTridiagonalMatrixAndGetTransformation(double * inputMatrix,
391  double * diagonalElements,
392  double * subDiagonalElements,
393  double * transformMatrix) const;
394 
424  unsigned int
425  ComputeEigenValuesUsingQL(double * d, double * e) const;
426 
464  unsigned int
465  ComputeEigenValuesAndVectorsUsingQL(double * d, double * e, double * z) const;
466 
467  /* Legacy algorithms using thread-safe netlib.
468  * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
469  */
470  unsigned int
471  ComputeEigenValuesLegacy(const TMatrix & A, TVector & EigenValues) const;
472 
473  unsigned int
474  ComputeEigenValuesAndVectorsLegacy(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
475 
476  /* Helper to get the matrix value type for EigenLibMatrix typename.
477  *
478  * If the TMatrix is vnl, the type is in element_type.
479  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
480  *
481  * To use this function:
482  * using ValueType = decltype(this->GetMatrixType(true));
483  *
484  * \note The two `GetMatrixValueType` overloads have different
485  * parameter declarations (`bool` and `...`), to avoid that both
486  * functions are equally good candidates during overload resolution,
487  * in case `element_type` and `ValueType` are both nested types of
488  * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
489  */
490  template <typename QMatrix = TMatrix>
491  auto
492  GetMatrixValueType(bool) const -> typename QMatrix::element_type
493  {
494  return QMatrix::element_type();
495  }
496  template <typename QMatrix = TMatrix>
497  auto
498  GetMatrixValueType(...) const -> typename QMatrix::ValueType
499  {
500  return QMatrix::ValueType();
501  }
502 
503  /* Wrapper that call the right implementation for the type of matrix. */
504  unsigned int
506  TVector & EigenValues,
507  TEigenMatrix & EigenVectors) const
508  {
509  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
510  }
511 
512  /* Implementation detail using EigenLib that performs a copy of the input matrix.
513  *
514  * @param (long) implementation detail argument making this implementation less favourable
515  * to be chosen if alternatives are available.
516  *
517  * @return an unsigned int with no information value (no error code in EigenLib) */
518  template <typename QMatrix>
519  auto
521  TVector & EigenValues,
522  TEigenMatrix & EigenVectors,
523  long) const -> decltype(static_cast<unsigned int>(1))
524  {
525  using ValueType = decltype(GetMatrixValueType(true));
526  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
527  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
528  for (unsigned int row = 0; row < m_Dimension; ++row)
529  {
530  for (unsigned int col = 0; col < m_Dimension; ++col)
531  {
532  inputMatrix(row, col) = A(row, col);
533  }
534  }
535  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
536  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
537  const auto & eigenValues = solver.eigenvalues();
538  /* Column k of the returned matrix is an eigenvector corresponding to
539  * eigenvalue number $ k $ as returned by eigenvalues().
540  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
541  const auto & eigenVectors = solver.eigenvectors();
542 
543  if (m_OrderEigenValues == OrderType::OrderByMagnitude)
544  {
545  auto copyEigenValues = eigenValues;
546  auto copyEigenVectors = eigenVectors;
547  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
548  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
549 
550  for (unsigned int row = 0; row < m_Dimension; ++row)
551  {
552  EigenValues[row] = copyEigenValues[row];
553  for (unsigned int col = 0; col < m_Dimension; ++col)
554  {
555  EigenVectors[row][col] = copyEigenVectors(col, row);
556  }
557  }
558  }
559  else
560  {
561  for (unsigned int row = 0; row < m_Dimension; ++row)
562  {
563  EigenValues[row] = eigenValues[row];
564  for (unsigned int col = 0; col < m_Dimension; ++col)
565  {
566  EigenVectors[row][col] = eigenVectors(col, row);
567  }
568  }
569  }
570  // No error code
571  return 1;
572  }
573 
574 
575  /* Implementation detail using EigenLib that do not peform a copy.
576  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
577  * If new types want to use this method, an appropiate overload of GetPointerToMatrixData
578  * should be included.
579  *
580  * @param (bool) implementation detail argument making this implementation the most favourable
581  * to be chosen from all the alternative implementations.
582  *
583  * @return an unsigned int with no information value (no error code in EigenLib) */
584  template <typename QMatrix>
585  auto
587  TVector & EigenValues,
588  TEigenMatrix & EigenVectors,
589  bool) const
590  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
591  {
592  auto pointerToData = GetPointerToMatrixData(A);
593  using PointerType = decltype(pointerToData);
594  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
595  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
596  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
597  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
598  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
599  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
600  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
601  const auto & eigenValues = solver.eigenvalues();
602  /* Column k of the returned matrix is an eigenvector corresponding to
603  * eigenvalue number $ k $ as returned by eigenvalues().
604  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
605  const auto & eigenVectors = solver.eigenvectors();
606  if (m_OrderEigenValues == OrderType::OrderByMagnitude)
607  {
608  auto copyEigenValues = eigenValues;
609  auto copyEigenVectors = eigenVectors;
610  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
611  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
612  for (unsigned int row = 0; row < m_Dimension; ++row)
613  {
614  EigenValues[row] = copyEigenValues[row];
615  for (unsigned int col = 0; col < m_Dimension; ++col)
616  {
617  EigenVectors[row][col] = copyEigenVectors(col, row);
618  }
619  }
620  }
621  else
622  {
623  for (unsigned int row = 0; row < m_Dimension; ++row)
624  {
625  EigenValues[row] = eigenValues[row];
626  for (unsigned int col = 0; col < m_Dimension; ++col)
627  {
628  EigenVectors[row][col] = eigenVectors(col, row);
629  }
630  }
631  }
632  // No error code
633  return 1;
634  }
635 
636  /* Wrapper that call the right implementation for the type of matrix. */
637  unsigned int
638  ComputeEigenValuesWithEigenLibrary(const TMatrix & A, TVector & EigenValues) const
639  {
640  return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
641  }
642 
643  /* Implementation detail using EigenLib that performs a copy of the input matrix.
644  *
645  * @param (long) implementation detail argument making this implementation less favourable
646  * to be chosen if alternatives are available.
647  *
648  * @return an unsigned int with no information value (no error code in EigenLib) */
649  template <typename QMatrix>
650  auto
651  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const
652  -> decltype(static_cast<unsigned int>(1))
653  {
654  using ValueType = decltype(GetMatrixValueType(true));
655  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
656  EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
657  for (unsigned int row = 0; row < m_Dimension; ++row)
658  {
659  for (unsigned int col = 0; col < m_Dimension; ++col)
660  {
661  inputMatrix(row, col) = A(row, col);
662  }
663  }
664  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
665  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
666  auto eigenValues = solver.eigenvalues();
667  if (m_OrderEigenValues == OrderType::OrderByMagnitude)
668  {
669  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
670  }
671  for (unsigned int i = 0; i < m_Dimension; ++i)
672  {
673  EigenValues[i] = eigenValues[i];
674  }
675 
676  // No error code
677  return 1;
678  }
679 
680  /* Implementation detail using EigenLib that do not peform a copy.
681  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
682  * If new types want to use this method, an appropiate overload of GetPointerToMatrixData
683  * should be included.
684  *
685  * @param (bool) implementation detail argument making this implementation the most favourable
686  * to be chosen from all the alternative implementations.
687  *
688  * @return an unsigned int with no information value (no error code in EigenLib) */
689  template <typename QMatrix>
690  auto
691  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
692  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
693  {
694  auto pointerToData = GetPointerToMatrixData(A);
695  using PointerType = decltype(pointerToData);
696  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
697  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
698  using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
699  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
700  EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
701  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
702  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
703  auto eigenValues = solver.eigenvalues();
704  if (m_OrderEigenValues == OrderType::OrderByMagnitude)
705  {
706  detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
707  }
708  for (unsigned int i = 0; i < m_Dimension; ++i)
709  {
710  EigenValues[i] = eigenValues[i];
711  }
712  // No error code
713  return 1;
714  }
715 };
716 
717 template <typename TMatrix, typename TVector, typename TEigenMatrix>
718 std::ostream &
720 {
721  os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
722  os << " Dimension : " << s.GetDimension() << std::endl;
723  os << " Order : " << s.GetOrder() << std::endl;
724  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
725  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
726  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
727  return os;
728 }
729 
730 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
731 class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysisFixedDimension
732 {
733 public:
736 #if !defined(ITK_LEGACY_REMOVE)
737  // We need to expose the enum values at the class level
738  // for backwards compatibility
742 #endif
743 
745  : m_OrderEigenValues(OrderType::OrderByValue)
746  {}
748 
749  using MatrixType = TMatrix;
750  using EigenMatrixType = TEigenMatrix;
751  using VectorType = TVector;
752 
766  unsigned int
767  ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const
768  {
769  return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
770  }
771 
792  unsigned int
793  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const
794  {
795  return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
796  }
797 
798  void
799  SetOrderEigenValues(const bool b)
800  {
801  if (b)
802  {
803  m_OrderEigenValues = OrderType::OrderByValue;
804  }
805  else
806  {
807  m_OrderEigenValues = OrderType::DoNotOrder;
808  }
809  }
810  bool
812  {
813  return (m_OrderEigenValues == OrderType::OrderByValue);
814  }
815  void
817  {
818  if (b)
819  {
820  m_OrderEigenValues = OrderType::OrderByMagnitude;
821  }
822  else
823  {
824  m_OrderEigenValues = OrderType::DoNotOrder;
825  }
826  }
827  bool
829  {
830  return (m_OrderEigenValues == OrderType::OrderByMagnitude);
831  }
832  constexpr unsigned int
833  GetOrder() const
834  {
835  return VDimension;
836  }
837  constexpr unsigned int
838  GetDimension() const
839  {
840  return VDimension;
841  }
842  constexpr bool
844  {
845  return true;
846  }
847 
848 private:
850 
851  /* Helper to get the matrix value type for EigenLibMatrix typename.
852  *
853  * If the TMatrix is vnl, the type is in element_type.
854  * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
855  *
856  * To use this function:
857  * using ValueType = decltype(this->GetMatrixType(true));
858  */
859  template <typename QMatrix = TMatrix>
860  auto
861  GetMatrixValueType(bool) const -> typename QMatrix::element_type
862  {
863  return QMatrix::element_type();
864  }
865  template <typename QMatrix = TMatrix>
866  auto
867  GetMatrixValueType(bool) const -> typename QMatrix::ValueType
868  {
869  return QMatrix::ValueType();
870  }
871 
872  /* Implementation detail using EigenLib that do not peform a copy.
873  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
874  * If new types want to use this method, an appropiate overload of GetPointerToMatrixData
875  * should be included.
876  *
877  * @param (bool) implementation detail argument making this implementation the most favourable
878  * to be chosen from all the alternative implementations.
879  *
880  * @return an unsigned int with no information value (no error code in EigenLib) */
881  template <typename QMatrix>
882  auto
884  TVector & EigenValues,
885  TEigenMatrix & EigenVectors,
886  bool) const
887  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
888  {
889  auto pointerToData = GetPointerToMatrixData(A);
890  using PointerType = decltype(pointerToData);
891  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
892  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
893  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
894  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
895  EigenConstMatrixMap inputMatrix(pointerToData);
896  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
897  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
898  const auto & eigenValues = solver.eigenvalues();
899  /* Column k of the returned matrix is an eigenvector corresponding to
900  * eigenvalue number $ k $ as returned by eigenvalues().
901  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
902  const auto & eigenVectors = solver.eigenvectors();
903  if (m_OrderEigenValues == OrderType::OrderByMagnitude)
904  {
905  auto copyEigenValues = eigenValues;
906  auto copyEigenVectors = eigenVectors;
907  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
908  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
909  for (unsigned int row = 0; row < VDimension; ++row)
910  {
911  EigenValues[row] = copyEigenValues[row];
912  for (unsigned int col = 0; col < VDimension; ++col)
913  {
914  EigenVectors[row][col] = copyEigenVectors(col, row);
915  }
916  }
917  }
918  else
919  {
920  for (unsigned int row = 0; row < VDimension; ++row)
921  {
922  EigenValues[row] = eigenValues[row];
923  for (unsigned int col = 0; col < VDimension; ++col)
924  {
925  EigenVectors[row][col] = eigenVectors(col, row);
926  }
927  }
928  }
929  // No error code
930  return 1;
931  }
932 
933  /* Implementation detail using EigenLib that performs a copy of the input matrix.
934  *
935  * @param (long) implementation detail argument making this implementation less favourable
936  * to be chosen if alternatives are available.
937  *
938  * @return an unsigned int with no information value (no error code in EigenLib) */
939  template <typename QMatrix>
940  auto
942  TVector & EigenValues,
943  TEigenMatrix & EigenVectors,
944  long) const -> decltype(static_cast<unsigned int>(1))
945  {
946  using ValueType = decltype(GetMatrixValueType(true));
947  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
948  EigenLibMatrixType inputMatrix;
949  for (unsigned int row = 0; row < VDimension; ++row)
950  {
951  for (unsigned int col = 0; col < VDimension; ++col)
952  {
953  inputMatrix(row, col) = A(row, col);
954  }
955  }
956  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
957  EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
958  const auto & eigenValues = solver.eigenvalues();
959  /* Column k of the returned matrix is an eigenvector corresponding to
960  * eigenvalue number $ k $ as returned by eigenvalues().
961  * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
962  const auto & eigenVectors = solver.eigenvectors();
963 
964  if (m_OrderEigenValues == OrderType::OrderByMagnitude)
965  {
966  auto copyEigenValues = eigenValues;
967  auto copyEigenVectors = eigenVectors;
968  auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
969  detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
970 
971  for (unsigned int row = 0; row < VDimension; ++row)
972  {
973  EigenValues[row] = copyEigenValues[row];
974  for (unsigned int col = 0; col < VDimension; ++col)
975  {
976  EigenVectors[row][col] = copyEigenVectors(col, row);
977  }
978  }
979  }
980  else
981  {
982  for (unsigned int row = 0; row < VDimension; ++row)
983  {
984  EigenValues[row] = eigenValues[row];
985  for (unsigned int col = 0; col < VDimension; ++col)
986  {
987  EigenVectors[row][col] = eigenVectors(col, row);
988  }
989  }
990  }
991  // No error code
992  return 1;
993  }
994 
995  /* Implementation detail using EigenLib that performs a copy of the input matrix.
996  *
997  * @param (long) implementation detail argument making this implementation less favourable
998  * to be chosen if alternatives are available.
999  *
1000  * @return an unsigned int with no information value (no error code in EigenLib) */
1001  template <typename QMatrix>
1002  auto
1003  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const
1004  -> decltype(static_cast<unsigned int>(1))
1005  {
1006  using ValueType = decltype(GetMatrixValueType(true));
1007  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1008  EigenLibMatrixType inputMatrix;
1009  for (unsigned int row = 0; row < VDimension; ++row)
1010  {
1011  for (unsigned int col = 0; col < VDimension; ++col)
1012  {
1013  inputMatrix(row, col) = A(row, col);
1014  }
1015  }
1016  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1017  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1018  auto eigenValues = solver.eigenvalues();
1019  if (m_OrderEigenValues == OrderType::OrderByMagnitude)
1020  {
1021  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1022  }
1023  for (unsigned int i = 0; i < VDimension; ++i)
1024  {
1025  EigenValues[i] = eigenValues[i];
1026  }
1027 
1028  // No error code
1029  return 1;
1030  }
1031 
1032  /* Implementation detail using EigenLib that do not peform a copy.
1033  * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
1034  * If new types want to use this method, an appropiate overload of GetPointerToMatrixData
1035  * should be included.
1036  *
1037  * @param (bool) implementation detail argument making this implementation the most favourable
1038  * to be chosen from all the alternative implementations.
1039  *
1040  * @return an unsigned int with no information value (no error code in EigenLib) */
1041  template <typename QMatrix>
1042  auto
1043  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
1044  -> decltype(GetPointerToMatrixData(A), static_cast<unsigned int>(1))
1045  {
1046  auto pointerToData = GetPointerToMatrixData(A);
1047  using PointerType = decltype(pointerToData);
1048  using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
1049  using ValueType = typename std::remove_cv<ValueTypeCV>::type;
1050  using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
1051  using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
1052  EigenConstMatrixMap inputMatrix(pointerToData);
1053  using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
1054  EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
1055  auto eigenValues = solver.eigenvalues();
1056  if (m_OrderEigenValues == OrderType::OrderByMagnitude)
1057  {
1058  detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
1059  }
1060  for (unsigned int i = 0; i < VDimension; ++i)
1061  {
1062  EigenValues[i] = eigenValues[i];
1063  }
1064  // No error code
1065  return 1;
1066  }
1067 };
1068 
1069 template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix>
1070 std::ostream &
1071 operator<<(std::ostream & os,
1073 {
1074  os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
1075  os << " Dimension : " << s.GetDimension() << std::endl;
1076  os << " Order : " << s.GetOrder() << std::endl;
1077  os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
1078  os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
1079  os << " UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
1080  return os;
1081 }
1082 } // end namespace itk
1083 
1084 #ifndef ITK_MANUAL_INSTANTIATION
1085 # include "itkSymmetricEigenAnalysis.hxx"
1086 #endif
1087 
1088 #endif
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::VectorType
TOutput VectorType
Definition: itkSymmetricEigenAnalysis.h:751
itk::uint8_t
::uint8_t uint8_t
Definition: itkIntTypes.h:29
itk::SymmetricEigenAnalysis::SymmetricEigenAnalysis
SymmetricEigenAnalysis()
Definition: itkSymmetricEigenAnalysis.h:176
itk::SymmetricEigenAnalysisFixedDimension::m_OrderEigenValues
EigenValueOrderType m_OrderEigenValues
Definition: itkSymmetricEigenAnalysis.h:849
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectors
unsigned int ComputeEigenValuesAndVectors(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:793
itk::SymmetricEigenAnalysis::SetDimension
void SetDimension(const unsigned int n)
Definition: itkSymmetricEigenAnalysis.h:299
itk::SymmetricEigenAnalysis::GetOrder
unsigned int GetOrder() const
Definition: itkSymmetricEigenAnalysis.h:245
itk::Matrix::GetVnlMatrix
InternalMatrixType & GetVnlMatrix()
Definition: itkMatrix.h:177
itkMatrix.h
itk::operator<<
std::ostream & operator<<(std::ostream &os, const Array< TValue > &arr)
Definition: itkArray.h:208
itk::SymmetricEigenAnalysis::GetUseEigenLibrary
bool GetUseEigenLibrary() const
Definition: itkSymmetricEigenAnalysis.h:334
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:520
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:586
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:691
itk::SymmetricEigenAnalysis::SetUseEigenLibraryOn
void SetUseEigenLibraryOn()
Definition: itkSymmetricEigenAnalysis.h:324
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::VectorType
TOutputImage::PixelType VectorType
Definition: itkSymmetricEigenAnalysis.h:192
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:883
itk::SymmetricEigenAnalysis::SymmetricEigenAnalysis
SymmetricEigenAnalysis(const unsigned int dimension)
Definition: itkSymmetricEigenAnalysis.h:180
itk::SymmetricEigenAnalysisFixedDimension::SetOrderEigenValues
void SetOrderEigenValues(const bool b)
Definition: itkSymmetricEigenAnalysis.h:799
itk::detail::permuteColumnsWithSortIndices
void permuteColumnsWithSortIndices(QMatrix &eigenVectors, const std::vector< int > &indicesSortPermutations)
Definition: itkSymmetricEigenAnalysis.h:101
EigenValueOrderType
itk::SymmetricEigenAnalysis::GetOrderEigenValues
bool GetOrderEigenValues() const
Definition: itkSymmetricEigenAnalysis.h:268
itk::SymmetricEigenAnalysisFixedDimension::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::ValueType
Definition: itkSymmetricEigenAnalysis.h:867
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
Definition: itkSymmetricEigenAnalysis.h:492
itk::SymmetricEigenAnalysis::ComputeEigenValuesAndVectorsWithEigenLibrary
unsigned int ComputeEigenValuesAndVectorsWithEigenLibrary(const TMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors) const
Definition: itkSymmetricEigenAnalysis.h:505
itk::SymmetricEigenAnalysisFixedDimension::SymmetricEigenAnalysisFixedDimension
SymmetricEigenAnalysisFixedDimension()
Definition: itkSymmetricEigenAnalysis.h:744
itk::detail::GetPointerToMatrixData
const TValueType * GetPointerToMatrixData(const vnl_matrix_fixed< TValueType, NRows, NCols > &inputMatrix)
Definition: itkSymmetricEigenAnalysis.h:38
itk::SymmetricEigenAnalysis::SetOrder
void SetOrder(const unsigned int n)
Definition: itkSymmetricEigenAnalysis.h:236
itkMacro.h
itk::SymmetricEigenAnalysis::GetDimension
unsigned int GetDimension() const
Definition: itkSymmetricEigenAnalysis.h:312
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:1003
itk::SymmetricEigenAnalysis::GetMatrixValueType
auto GetMatrixValueType(...) const -> typename QMatrix::ValueType
Definition: itkSymmetricEigenAnalysis.h:498
itk::detail::sortEigenValuesByMagnitude
std::vector< int > sortEigenValuesByMagnitude(TArray &eigenValues, const unsigned int numberOfElements)
Definition: itkSymmetricEigenAnalysis.h:73
itk::SymmetricEigenAnalysisFixedDimension::SetOrderEigenMagnitudes
void SetOrderEigenMagnitudes(const bool b)
Definition: itkSymmetricEigenAnalysis.h:816
itk::SymmetricEigenAnalysisFixedDimension::GetOrderEigenMagnitudes
bool GetOrderEigenMagnitudes() const
Definition: itkSymmetricEigenAnalysis.h:828
itk::SymmetricEigenAnalysis::SetOrderEigenMagnitudes
void SetOrderEigenMagnitudes(const bool b)
Definition: itkSymmetricEigenAnalysis.h:277
itk::SymmetricEigenAnalysis::SetUseEigenLibraryOff
void SetUseEigenLibraryOff()
Definition: itkSymmetricEigenAnalysis.h:329
itk::OrderType::OrderByValue
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::EigenMatrixType
TInput EigenMatrixType
Definition: itkSymmetricEigenAnalysis.h:750
itk::Matrix
A templated class holding a M x N size Matrix.
Definition: itkMatrix.h:50
itk::OrderType::DoNotOrder
itk::SymmetricEigenAnalysisFixedDimension::GetMatrixValueType
auto GetMatrixValueType(bool) const -> typename QMatrix::element_type
Definition: itkSymmetricEigenAnalysis.h:861
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, bool) const -> decltype(GetPointerToMatrixData(A), static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:1043
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::EigenMatrixType
TInputImage::PixelType EigenMatrixType
Definition: itkSymmetricEigenAnalysis.h:191
itk::SymmetricEigenAnalysis::SetUseEigenLibrary
void SetUseEigenLibrary(const bool input)
Definition: itkSymmetricEigenAnalysis.h:319
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValuesAndVectorsWithEigenLibraryImpl
auto ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, TEigenMatrix &EigenVectors, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:941
itk::SymmetricEigenAnalysisFixedDimension::GetDimension
constexpr unsigned int GetDimension() const
Definition: itkSymmetricEigenAnalysis.h:838
itk::SymmetricEigenAnalysisFixedDimension
Definition: itkSymmetricEigenAnalysis.h:731
itk::OrderType
OrderType
Definition: itkSymmetricEigenAnalysis.h:120
itk::SymmetricEigenAnalysisFixedDimension::GetUseEigenLibrary
constexpr bool GetUseEigenLibrary() const
Definition: itkSymmetricEigenAnalysis.h:843
itk::SymmetricEigenAnalysis::SetOrderEigenValues
void SetOrderEigenValues(const bool b)
Definition: itkSymmetricEigenAnalysis.h:254
itk
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
Definition: itkArray.h:26
itk::SymmetricEigenAnalysisFixedDimension::GetOrderEigenValues
bool GetOrderEigenValues() const
Definition: itkSymmetricEigenAnalysis.h:811
itk::SymmetricEigenAnalysis::m_OrderEigenValues
EigenValueOrderType m_OrderEigenValues
Definition: itkSymmetricEigenAnalysis.h:344
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibraryImpl
auto ComputeEigenValuesWithEigenLibraryImpl(const QMatrix &A, TVector &EigenValues, long) const -> decltype(static_cast< unsigned int >(1))
Definition: itkSymmetricEigenAnalysis.h:651
itk::OrderType::OrderByMagnitude
itk::Math::e
static constexpr double e
Definition: itkMath.h:53
itk::SymmetricEigenAnalysis::ComputeEigenValuesWithEigenLibrary
unsigned int ComputeEigenValuesWithEigenLibrary(const TMatrix &A, TVector &EigenValues) const
Definition: itkSymmetricEigenAnalysis.h:638
itk::SymmetricEigenAnalysisFixedDimension< TMatrixDimension, TInput, TOutput >::MatrixType
TInput MatrixType
Definition: itkSymmetricEigenAnalysis.h:749
itk::SymmetricEigenAnalysisFixedDimension::GetOrder
constexpr unsigned int GetOrder() const
Definition: itkSymmetricEigenAnalysis.h:833
itk::SymmetricEigenAnalysis< TInputImage::PixelType, TOutputImage::PixelType >::MatrixType
TInputImage::PixelType MatrixType
Definition: itkSymmetricEigenAnalysis.h:190
itk::detail::GetPointerToMatrixData
const TValueType * GetPointerToMatrixData(const itk::Matrix< TValueType, NRows, NCols > &inputMatrix)
Definition: itkSymmetricEigenAnalysis.h:51
itk::SymmetricEigenAnalysis::GetOrderEigenMagnitudes
bool GetOrderEigenMagnitudes() const
Definition: itkSymmetricEigenAnalysis.h:291
itk::SymmetricEigenAnalysisFixedDimension::ComputeEigenValues
unsigned int ComputeEigenValues(const TMatrix &A, TVector &EigenValues) const
Definition: itkSymmetricEigenAnalysis.h:767
OrderType
itk::SymmetricEigenAnalysis
Find Eigen values of a real 2D symmetric matrix. It serves as a thread-safe alternative to the class:...
Definition: itkSymmetricEigenAnalysis.h:163