itkSymmetricEigenAnalysis.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkSymmetricEigenAnalysis_h
#define itkSymmetricEigenAnalysis_h
 
#include "itkMacro.h"
#include "itk_eigen.h"
#include ITK_EIGEN(Eigenvalues)
#include <numeric>
#include <vector>
// For GetPointerToMatrixData
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_matrix_fixed.h"
#include "itkMatrix.h"
 
namespace itk
{
namespace detail
{
/* Helper functions returning pointer to matrix data for different types.  */
template <typename TValueType, unsigned int VRows, unsigned int VColumns>
const TValueType *
GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, VRows, VColumns> & inputMatrix)
{
  return inputMatrix.data_block();
}
template <typename TValueType>
const TValueType *
GetPointerToMatrixData(const vnl_matrix<TValueType> & inputMatrix)
{
  return inputMatrix.data_block();
}
 
template <typename TValueType, unsigned int VRows, unsigned int VColumns>
const TValueType *
GetPointerToMatrixData(const itk::Matrix<TValueType, VRows, VColumns> & inputMatrix)
{
  return inputMatrix.GetVnlMatrix().data_block();
}

template <typename TArray>
std::vector<int>
sortEigenValuesByMagnitude(TArray & eigenValues, const unsigned int numberOfElements)
{
  std::vector<int> indicesSortPermutations(numberOfElements, 0);
  std::iota(std::begin(indicesSortPermutations), std::end(indicesSortPermutations), 0);
 
  std::sort(std::begin(indicesSortPermutations),
            std::end(indicesSortPermutations),
            [&eigenValues](unsigned int a, unsigned int b) {
              return itk::Math::abs(eigenValues[a]) < itk::Math::abs(eigenValues[b]);
            });
  auto tmpCopy = eigenValues;
  for (unsigned int i = 0; i < numberOfElements; ++i)
  {
    eigenValues[i] = tmpCopy[indicesSortPermutations[i]];
  }
  return indicesSortPermutations;
}

template <typename QMatrix>
void
permuteColumnsWithSortIndices(QMatrix & eigenVectors, const std::vector<int> & indicesSortPermutations)
{
  using EigenLibPermutationMatrix = Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
  auto numberOfElements = indicesSortPermutations.size();
  // Creates a NxN permutation matrix copying our permutation to the matrix indices.
  // Which holds the 1D array representation of a permutation.
  EigenLibPermutationMatrix perm(numberOfElements);
  perm.setIdentity();
  std::copy(indicesSortPermutations.begin(), indicesSortPermutations.end(), perm.indices().data());
  // Apply it
  eigenVectors = eigenVectors * perm;
}

 
} // end namespace detail

class SymmetricEigenAnalysisEnums
{
public:
  enum class EigenValueOrder : uint8_t
  {
    OrderByValue = 1,
    OrderByMagnitude = 2,
    DoNotOrder = 3
  };
};
// Define how to print enumeration
extern ITKCommon_EXPORT std::ostream &
                        operator<<(std::ostream & out, const SymmetricEigenAnalysisEnums::EigenValueOrder value);
 
using EigenValueOrderEnum = SymmetricEigenAnalysisEnums::EigenValueOrder;
 
inline EigenValueOrderEnum
Int2EigenValueOrderEnum(const uint8_t value)
{
  switch (value)
  {
    case 1:
      return EigenValueOrderEnum::OrderByValue;
    case 2:
      return EigenValueOrderEnum::OrderByMagnitude;
    case 3:
      return EigenValueOrderEnum::DoNotOrder;
    default:
      break;
  }
  itkGenericExceptionMacro("Error: Invalid value for conversion.");
}
 
#if !defined(ITK_LEGACY_REMOVE)
static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
#endif

 
template <typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
{
public:
  using EigenValueOrderEnum = itk::EigenValueOrderEnum;
#if !defined(ITK_LEGACY_REMOVE)
  using EigenValueOrderType = EigenValueOrderEnum;
#endif
 
  SymmetricEigenAnalysis() = default;
 
  SymmetricEigenAnalysis(const unsigned int dimension)
    : m_Dimension(dimension)
    , m_Order(dimension)
  {}
 
  ~SymmetricEigenAnalysis() = default;
 
  using MatrixType = TMatrix;
  using EigenMatrixType = TEigenMatrix;
  using VectorType = TVector;

  unsigned int
  ComputeEigenValues(const TMatrix & A, TVector & D) const;

  unsigned int
  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
 

  void
  SetOrder(const unsigned int n)
  {
    m_Order = n;
  }

  unsigned int
  GetOrder() const
  {
    return m_Order;
  }

  void
  SetOrderEigenValues(const bool b)
  {
    if (b)
    {
      m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
    }
    else
    {
      m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
    }
  }

 
  bool
  GetOrderEigenValues() const
  {
    return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
  }

  void
  SetOrderEigenMagnitudes(const bool b)
  {
    if (b)
    {
      m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
    }
    else
    {
      m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
    }
  }

 
  bool
  GetOrderEigenMagnitudes() const
  {
    return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
  }

  void
  SetDimension(const unsigned int n)
  {
    m_Dimension = n;
    if (m_Order == 0)
    {
      m_Order = m_Dimension;
    }
  }


  unsigned int
  GetDimension() const
  {
    return m_Dimension;
  }

  void
  SetUseEigenLibrary(const bool input)
  {
    m_UseEigenLibrary = input;
  }
  void
  SetUseEigenLibraryOn()
  {
    m_UseEigenLibrary = true;
  }
  void
  SetUseEigenLibraryOff()
  {
    m_UseEigenLibrary = false;
  }
  bool
  GetUseEigenLibrary() const
  {
    return m_UseEigenLibrary;
  }

 
private:
  bool                m_UseEigenLibrary{ false };
  unsigned int        m_Dimension{ 0 };
  unsigned int        m_Order{ 0 };
  EigenValueOrderEnum m_OrderEigenValues{ EigenValueOrderEnum::OrderByValue };

  void
  ReduceToTridiagonalMatrix(double * a, double * d, double * e, double * e2) const;

  void
  ReduceToTridiagonalMatrixAndGetTransformation(const double * a, double * d, double * e, double * z) const;

  unsigned int
  ComputeEigenValuesUsingQL(double * d, double * e) const;

  unsigned int
  ComputeEigenValuesAndVectorsUsingQL(double * d, double * e, double * z) const;
 
  /* Legacy algorithms using thread-safe netlib.
   * \sa ComputeEigenValues and \sa ComputeEigenValuesAndVectors
   */
  unsigned int
  ComputeEigenValuesLegacy(const TMatrix & A, TVector & D) const;
 
  unsigned int
  ComputeEigenValuesAndVectorsLegacy(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const;
 
  /* Helper to get the matrix value type for EigenLibMatrix typename.
   *
   * If the TMatrix is vnl, the type is in element_type.
   * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
   *
   * To use this function:
   * using ValueType = decltype(this->GetMatrixType(true));
   *
   * \note The two `GetMatrixValueType` overloads have different
   * parameter declarations (`bool` and `...`), to avoid that both
   * functions are equally good candidates during overload resolution,
   * in case `element_type` and `ValueType` are both nested types of
   * `TMatrix` (which is the case when `TMatrix` = `itk::Array2D`).
   */
  template <typename QMatrix = TMatrix>
  auto
  GetMatrixValueType(bool) const -> typename QMatrix::element_type
  {
    return QMatrix::element_type();
  }
  template <typename QMatrix = TMatrix>
  auto
  GetMatrixValueType(...) const -> typename QMatrix::ValueType
  {
    return QMatrix::ValueType();
  }
 
  /* Wrapper that call the right implementation for the type of matrix.  */
  unsigned int
  ComputeEigenValuesAndVectorsWithEigenLibrary(const TMatrix & A,
                                               TVector &       EigenValues,
                                               TEigenMatrix &  EigenVectors) const
  {
    return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
  }
 
  /* Implementation detail using EigenLib that performs a copy of the input matrix.
   *
   * @param (long) implementation detail argument making this implementation less favourable
   *   to be chosen if alternatives are available.
   *
   * @return an unsigned int with no information value (no error code in EigenLib) */
  template <typename QMatrix>
  auto
  ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix & A,
                                                   TVector &       EigenValues,
                                                   TEigenMatrix &  EigenVectors,
                                                   long) const -> decltype(1U)
  {
    using ValueType = decltype(GetMatrixValueType(true));
    using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
    EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
    for (unsigned int row = 0; row < m_Dimension; ++row)
    {
      for (unsigned int col = 0; col < m_Dimension; ++col)
      {
        inputMatrix(row, col) = A(row, col);
      }
    }
    using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
    const EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
    const auto &          eigenValues = solver.eigenvalues();
    /* Column  k  of the returned matrix is an eigenvector corresponding to
     * eigenvalue number $ k $ as returned by eigenvalues().
     * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
    const auto & eigenVectors = solver.eigenvectors();
 
    if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
    {
      auto copyEigenValues = eigenValues;
      auto copyEigenVectors = eigenVectors;
      auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
      detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
 
      for (unsigned int row = 0; row < m_Dimension; ++row)
      {
        EigenValues[row] = copyEigenValues[row];
        for (unsigned int col = 0; col < m_Dimension; ++col)
        {
          EigenVectors[row][col] = copyEigenVectors(col, row);
        }
      }
    }
    else
    {
      for (unsigned int row = 0; row < m_Dimension; ++row)
      {
        EigenValues[row] = eigenValues[row];
        for (unsigned int col = 0; col < m_Dimension; ++col)
        {
          EigenVectors[row][col] = eigenVectors(col, row);
        }
      }
    }
    // No error code
    return 1;
  }
 
 
  /* Implementation detail using EigenLib that do not perform a copy.
   * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
   * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
   * should be included.
   *
   * @param (bool) implementation detail argument making this implementation the most favourable
   *   to be chosen from all the alternative implementations.
   *
   * @return an unsigned int with no information value (no error code in EigenLib) */
  template <typename QMatrix>
  auto
  ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix & A,
                                                   TVector &       EigenValues,
                                                   TEigenMatrix &  EigenVectors,
                                                   bool) const -> decltype(GetPointerToMatrixData(A), 1U)
  {
    auto pointerToData = GetPointerToMatrixData(A);
    using PointerType = decltype(pointerToData);
    using ValueTypeCV = std::remove_pointer_t<PointerType>;
    using ValueType = std::remove_cv_t<ValueTypeCV>;
    using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
    using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
    EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
    using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
    EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
    const auto &    eigenValues = solver.eigenvalues();
    /* Column  k  of the returned matrix is an eigenvector corresponding to
     * eigenvalue number $ k $ as returned by eigenvalues().
     * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
    const auto & eigenVectors = solver.eigenvectors();
    if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
    {
      auto copyEigenValues = eigenValues;
      auto copyEigenVectors = eigenVectors;
      auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, m_Dimension);
      detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
      for (unsigned int row = 0; row < m_Dimension; ++row)
      {
        EigenValues[row] = copyEigenValues[row];
        for (unsigned int col = 0; col < m_Dimension; ++col)
        {
          EigenVectors[row][col] = copyEigenVectors(col, row);
        }
      }
    }
    else
    {
      for (unsigned int row = 0; row < m_Dimension; ++row)
      {
        EigenValues[row] = eigenValues[row];
        for (unsigned int col = 0; col < m_Dimension; ++col)
        {
          EigenVectors[row][col] = eigenVectors(col, row);
        }
      }
    }
    // No error code
    return 1;
  }
 
  /* Wrapper that call the right implementation for the type of matrix.  */
  unsigned int
  ComputeEigenValuesWithEigenLibrary(const TMatrix & A, TVector & EigenValues) const
  {
    return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
  }
 
  /* Implementation detail using EigenLib that performs a copy of the input matrix.
   *
   * @param (long) implementation detail argument making this implementation less favourable
   *   to be chosen if alternatives are available.
   *
   * @return an unsigned int with no information value (no error code in EigenLib) */
  template <typename QMatrix>
  auto
  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const -> decltype(1U)
  {
    using ValueType = decltype(GetMatrixValueType(true));
    using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
    EigenLibMatrixType inputMatrix(m_Dimension, m_Dimension);
    for (unsigned int row = 0; row < m_Dimension; ++row)
    {
      for (unsigned int col = 0; col < m_Dimension; ++col)
      {
        inputMatrix(row, col) = A(row, col);
      }
    }
    using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
    const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
    auto                  eigenValues = solver.eigenvalues();
    if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
    {
      detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
    }
    for (unsigned int i = 0; i < m_Dimension; ++i)
    {
      EigenValues[i] = eigenValues[i];
    }
 
    // No error code
    return 1;
  }
 
  /* Implementation detail using EigenLib that do not perform a copy.
   * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
   * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
   * should be included.
   *
   * @param (bool) implementation detail argument making this implementation the most favourable
   *   to be chosen from all the alternative implementations.
   *
   * @return an unsigned int with no information value (no error code in EigenLib) */
  template <typename QMatrix>
  auto
  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
    -> decltype(GetPointerToMatrixData(A), 1U)
  {
    auto pointerToData = GetPointerToMatrixData(A);
    using PointerType = decltype(pointerToData);
    using ValueTypeCV = std::remove_pointer_t<PointerType>;
    using ValueType = std::remove_cv_t<ValueTypeCV>;
    using EigenLibMatrixType = Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
    using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
    EigenConstMatrixMap inputMatrix(pointerToData, m_Dimension, m_Dimension);
    using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
    EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
    auto            eigenValues = solver.eigenvalues();
    if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
    {
      detail::sortEigenValuesByMagnitude(eigenValues, m_Dimension);
    }
    for (unsigned int i = 0; i < m_Dimension; ++i)
    {
      EigenValues[i] = eigenValues[i];
    }
    // No error code
    return 1;
  }
};
 
template <typename TMatrix, typename TVector, typename TEigenMatrix>
std::ostream &
operator<<(std::ostream & os, const SymmetricEigenAnalysis<TMatrix, TVector, TEigenMatrix> & s)
{
  os << "[ClassType: SymmetricEigenAnalysis]" << std::endl;
  os << "  Dimension : " << s.GetDimension() << std::endl;
  os << "  Order : " << s.GetOrder() << std::endl;
  os << "  OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
  os << "  OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
  os << "  UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
  return os;
}
 
template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix>
class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysisFixedDimension
{
public:
#if !defined(ITK_LEGACY_REMOVE)
  using EigenValueOrderType = EigenValueOrderEnum;
#endif
#if !defined(ITK_LEGACY_REMOVE)
  // We need to expose the enum values at the class level
  // for backwards compatibility
  static constexpr EigenValueOrderEnum OrderByValue = EigenValueOrderEnum::OrderByValue;
  static constexpr EigenValueOrderEnum OrderByMagnitude = EigenValueOrderEnum::OrderByMagnitude;
  static constexpr EigenValueOrderEnum DoNotOrder = EigenValueOrderEnum::DoNotOrder;
#endif
 
  using MatrixType = TMatrix;
  using EigenMatrixType = TEigenMatrix;
  using VectorType = TVector;

  unsigned int
  ComputeEigenValues(const TMatrix & A, TVector & EigenValues) const
  {
    return ComputeEigenValuesWithEigenLibraryImpl(A, EigenValues, true);
  }

  unsigned int
  ComputeEigenValuesAndVectors(const TMatrix & A, TVector & EigenValues, TEigenMatrix & EigenVectors) const
  {
    return ComputeEigenValuesAndVectorsWithEigenLibraryImpl(A, EigenValues, EigenVectors, true);
  }
 
  void
  SetOrderEigenValues(const bool b)
  {
    if (b)
    {
      m_OrderEigenValues = EigenValueOrderEnum::OrderByValue;
    }
    else
    {
      m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
    }
  }
  bool
  GetOrderEigenValues() const
  {
    return (m_OrderEigenValues == EigenValueOrderEnum::OrderByValue);
  }
  void
  SetOrderEigenMagnitudes(const bool b)
  {
    if (b)
    {
      m_OrderEigenValues = EigenValueOrderEnum::OrderByMagnitude;
    }
    else
    {
      m_OrderEigenValues = EigenValueOrderEnum::DoNotOrder;
    }
  }
  bool
  GetOrderEigenMagnitudes() const
  {
    return (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude);
  }
  constexpr unsigned int
  GetOrder() const
  {
    return VDimension;
  }
  constexpr unsigned int
  GetDimension() const
  {
    return VDimension;
  }
  constexpr bool
  GetUseEigenLibrary() const
  {
    return true;
  }
 
private:
  EigenValueOrderEnum m_OrderEigenValues{ EigenValueOrderEnum::OrderByValue };
 
  /* Helper to get the matrix value type for EigenLibMatrix typename.
   *
   * If the TMatrix is vnl, the type is in element_type.
   * In TMatrix is itk::Matrix, or any itk::FixedArray is in ValueType.
   *
   * To use this function:
   * using ValueType = decltype(this->GetMatrixType(true));
   */
  template <typename QMatrix = TMatrix>
  auto
  GetMatrixValueType(bool) const -> typename QMatrix::element_type
  {
    return QMatrix::element_type();
  }
  template <typename QMatrix = TMatrix>
  auto
  GetMatrixValueType(bool) const -> typename QMatrix::ValueType
  {
    return QMatrix::ValueType();
  }
 
  /* Implementation detail using EigenLib that do not perform a copy.
   * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
   * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
   * should be included.
   *
   * @param (bool) implementation detail argument making this implementation the most favourable
   *   to be chosen from all the alternative implementations.
   *
   * @return an unsigned int with no information value (no error code in EigenLib) */
  template <typename QMatrix>
  auto
  ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix & A,
                                                   TVector &       EigenValues,
                                                   TEigenMatrix &  EigenVectors,
                                                   bool) const -> decltype(GetPointerToMatrixData(A), 1U)
  {
    auto pointerToData = GetPointerToMatrixData(A);
    using PointerType = decltype(pointerToData);
    using ValueTypeCV = std::remove_pointer_t<PointerType>;
    using ValueType = std::remove_cv_t<ValueTypeCV>;
    using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
    using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
    EigenConstMatrixMap inputMatrix(pointerToData);
    using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
    EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
    const auto &    eigenValues = solver.eigenvalues();
    /* Column  k  of the returned matrix is an eigenvector corresponding to
     * eigenvalue number $ k $ as returned by eigenvalues().
     * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
    const auto & eigenVectors = solver.eigenvectors();
    if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
    {
      auto copyEigenValues = eigenValues;
      auto copyEigenVectors = eigenVectors;
      auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
      detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
      for (unsigned int row = 0; row < VDimension; ++row)
      {
        EigenValues[row] = copyEigenValues[row];
        for (unsigned int col = 0; col < VDimension; ++col)
        {
          EigenVectors[row][col] = copyEigenVectors(col, row);
        }
      }
    }
    else
    {
      for (unsigned int row = 0; row < VDimension; ++row)
      {
        EigenValues[row] = eigenValues[row];
        for (unsigned int col = 0; col < VDimension; ++col)
        {
          EigenVectors[row][col] = eigenVectors(col, row);
        }
      }
    }
    // No error code
    return 1;
  }
 
  /* Implementation detail using EigenLib that performs a copy of the input matrix.
   *
   * @param (long) implementation detail argument making this implementation less favourable
   *   to be chosen if alternatives are available.
   *
   * @return an unsigned int with no information value (no error code in EigenLib) */
  template <typename QMatrix>
  auto
  ComputeEigenValuesAndVectorsWithEigenLibraryImpl(const QMatrix & A,
                                                   TVector &       EigenValues,
                                                   TEigenMatrix &  EigenVectors,
                                                   long) const -> decltype(1U)
  {
    using ValueType = decltype(GetMatrixValueType(true));
    using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
    EigenLibMatrixType inputMatrix;
    for (unsigned int row = 0; row < VDimension; ++row)
    {
      for (unsigned int col = 0; col < VDimension; ++col)
      {
        inputMatrix(row, col) = A(row, col);
      }
    }
    using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
    const EigenSolverType solver(inputMatrix); // Computes EigenValues and EigenVectors
    const auto &          eigenValues = solver.eigenvalues();
    /* Column  k  of the returned matrix is an eigenvector corresponding to
     * eigenvalue number $ k $ as returned by eigenvalues().
     * The eigenvectors are normalized to have (Euclidean) norm equal to one. */
    const auto & eigenVectors = solver.eigenvectors();
 
    if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
    {
      auto copyEigenValues = eigenValues;
      auto copyEigenVectors = eigenVectors;
      auto indicesSortPermutations = detail::sortEigenValuesByMagnitude(copyEigenValues, VDimension);
      detail::permuteColumnsWithSortIndices(copyEigenVectors, indicesSortPermutations);
 
      for (unsigned int row = 0; row < VDimension; ++row)
      {
        EigenValues[row] = copyEigenValues[row];
        for (unsigned int col = 0; col < VDimension; ++col)
        {
          EigenVectors[row][col] = copyEigenVectors(col, row);
        }
      }
    }
    else
    {
      for (unsigned int row = 0; row < VDimension; ++row)
      {
        EigenValues[row] = eigenValues[row];
        for (unsigned int col = 0; col < VDimension; ++col)
        {
          EigenVectors[row][col] = eigenVectors(col, row);
        }
      }
    }
    // No error code
    return 1;
  }
 
  /* Implementation detail using EigenLib that performs a copy of the input matrix.
   *
   * @param (long) implementation detail argument making this implementation less favourable
   *   to be chosen if alternatives are available.
   *
   * @return an unsigned int with no information value (no error code in EigenLib) */
  template <typename QMatrix>
  auto
  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, long) const -> decltype(1U)
  {
    using ValueType = decltype(GetMatrixValueType(true));
    using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
    EigenLibMatrixType inputMatrix;
    for (unsigned int row = 0; row < VDimension; ++row)
    {
      for (unsigned int col = 0; col < VDimension; ++col)
      {
        inputMatrix(row, col) = A(row, col);
      }
    }
    using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
    const EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
    auto                  eigenValues = solver.eigenvalues();
    if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
    {
      detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
    }
    for (unsigned int i = 0; i < VDimension; ++i)
    {
      EigenValues[i] = eigenValues[i];
    }
 
    // No error code
    return 1;
  }
 
  /* Implementation detail using EigenLib that do not perform a copy.
   * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
   * If new types want to use this method, an appropriate overload of GetPointerToMatrixData
   * should be included.
   *
   * @param (bool) implementation detail argument making this implementation the most favourable
   *   to be chosen from all the alternative implementations.
   *
   * @return an unsigned int with no information value (no error code in EigenLib) */
  template <typename QMatrix>
  auto
  ComputeEigenValuesWithEigenLibraryImpl(const QMatrix & A, TVector & EigenValues, bool) const
    -> decltype(GetPointerToMatrixData(A), 1U)
  {
    auto pointerToData = GetPointerToMatrixData(A);
    using PointerType = decltype(pointerToData);
    using ValueTypeCV = std::remove_pointer_t<PointerType>;
    using ValueType = std::remove_cv_t<ValueTypeCV>;
    using EigenLibMatrixType = Eigen::Matrix<ValueType, VDimension, VDimension, Eigen::RowMajor>;
    using EigenConstMatrixMap = Eigen::Map<const EigenLibMatrixType>;
    EigenConstMatrixMap inputMatrix(pointerToData);
    using EigenSolverType = Eigen::SelfAdjointEigenSolver<EigenLibMatrixType>;
    EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
    auto            eigenValues = solver.eigenvalues();
    if (m_OrderEigenValues == EigenValueOrderEnum::OrderByMagnitude)
    {
      detail::sortEigenValuesByMagnitude(eigenValues, VDimension);
    }
    for (unsigned int i = 0; i < VDimension; ++i)
    {
      EigenValues[i] = eigenValues[i];
    }
    // No error code
    return 1;
  }
};
 
template <unsigned int VDimension, typename TMatrix, typename TVector, typename TEigenMatrix>
std::ostream &
operator<<(std::ostream &                                                                           os,
           const SymmetricEigenAnalysisFixedDimension<VDimension, TMatrix, TVector, TEigenMatrix> & s)
{
  os << "[ClassType: SymmetricEigenAnalysisFixedDimension]" << std::endl;
  os << "  Dimension : " << s.GetDimension() << std::endl;
  os << "  Order : " << s.GetOrder() << std::endl;
  os << "  OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
  os << "  OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
  os << "  UseEigenLibrary: " << s.GetUseEigenLibrary() << std::endl;
  return os;
}
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkSymmetricEigenAnalysis.hxx"
#endif
 
#endif