itkSymmetricEigenAnalysis.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  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
 *
 *         http://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 NRows, unsigned int NCols>
    const TValueType* GetPointerToMatrixData(const vnl_matrix_fixed<TValueType, NRows, NCols> & 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 NRows, unsigned int NCols>
    const TValueType* GetPointerToMatrixData(const itk::Matrix<TValueType, NRows, NCols> & 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 std::abs(eigenValues[a]) < std::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

template< typename TMatrix, typename TVector, typename TEigenMatrix = TMatrix >
class ITK_TEMPLATE_EXPORT SymmetricEigenAnalysis
{
public:
  typedef enum
  {
    OrderByValue = 1,
    OrderByMagnitude,
    DoNotOrder
  }
  EigenValueOrderType;

  SymmetricEigenAnalysis() :
  m_OrderEigenValues(OrderByValue) {}

  SymmetricEigenAnalysis(const unsigned int dimension):
  m_UseEigenLibrary(false),
  m_Dimension(dimension),
  m_Order(dimension),
  m_OrderEigenValues(OrderByValue) {}

  ~SymmetricEigenAnalysis() = default;

  using MatrixType = TMatrix;
  using EigenMatrixType = TEigenMatrix;
  using VectorType = TVector;

  unsigned int ComputeEigenValues(
    const TMatrix  & A,
    TVector        & EigenValues) 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 = OrderByValue;     }
    else   { m_OrderEigenValues = DoNotOrder;       }
  }

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

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

  bool GetOrderEigenMagnitudes() const { return ( m_OrderEigenValues == 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};
  EigenValueOrderType m_OrderEigenValues;

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

  void ReduceToTridiagonalMatrixAndGetTransformation(
    double *inputMatrix, double *diagonalElements,
    double *subDiagonalElements, double *transformMatrix) 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        & EigenValues) 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(static_cast<unsigned int>(1))
  {
  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>;
  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 == 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 peform a copy.
   * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
   * If new types want to use this method, an appropiate 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), static_cast<unsigned int>(1))
    {
    auto pointerToData = GetPointerToMatrixData(A);
    using PointerType = decltype(pointerToData);
    using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
    using ValueType = typename std::remove_cv<ValueTypeCV>::type;
    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 == 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(static_cast<unsigned int>(1))
    {
    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>;
    EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
    auto eigenValues = solver.eigenvalues();
    if(m_OrderEigenValues == 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 peform a copy.
   * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
   * If new types want to use this method, an appropiate 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), static_cast<unsigned int>(1))
    {
    auto pointerToData = GetPointerToMatrixData(A);
    using PointerType = decltype(pointerToData);
    using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
    using ValueType = typename std::remove_cv<ValueTypeCV>::type;
    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 == 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:
  typedef enum
  {
    OrderByValue = 1,
    OrderByMagnitude,
    DoNotOrder
  }
  EigenValueOrderType;

  SymmetricEigenAnalysisFixedDimension(): m_OrderEigenValues(OrderByValue) {}
  ~SymmetricEigenAnalysisFixedDimension() = default;

  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 = OrderByValue;     }
    else   { m_OrderEigenValues = DoNotOrder;       }
  }
  bool GetOrderEigenValues() const { return ( m_OrderEigenValues == OrderByValue ); }
  void SetOrderEigenMagnitudes(const bool b)
  {
    if ( b ) { m_OrderEigenValues = OrderByMagnitude; }
    else   { m_OrderEigenValues = DoNotOrder;       }
  }
  bool GetOrderEigenMagnitudes() const { return ( m_OrderEigenValues == OrderByMagnitude ); }
  constexpr unsigned int GetOrder() const { return VDimension; }
  constexpr unsigned int GetDimension() const { return VDimension; }
  constexpr bool GetUseEigenLibrary() const { return true; }

private:
  EigenValueOrderType m_OrderEigenValues;

  /* 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 peform a copy.
   * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
   * If new types want to use this method, an appropiate 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), static_cast<unsigned int>(1))
    {
    auto pointerToData = GetPointerToMatrixData(A);
    using PointerType = decltype(pointerToData);
    using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
    using ValueType = typename std::remove_cv<ValueTypeCV>::type;
    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 == 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(static_cast<unsigned int>(1))
  {
  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>;
  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 == 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(static_cast<unsigned int>(1))
    {
    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>;
    EigenSolverType solver(inputMatrix, Eigen::EigenvaluesOnly);
    auto eigenValues = solver.eigenvalues();
    if(m_OrderEigenValues == 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 peform a copy.
   * It needs the existence of a pointer to matrix data. \sa GetPointerToMatrixData
   * If new types want to use this method, an appropiate 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), static_cast<unsigned int>(1))
    {
    auto pointerToData = GetPointerToMatrixData(A);
    using PointerType = decltype(pointerToData);
    using ValueTypeCV = typename std::remove_pointer<PointerType>::type;
    using ValueType = typename std::remove_cv<ValueTypeCV>::type;
    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 == 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