itkSymmetricEigenDecomposition.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 itkSymmetricEigenDecomposition_h
#define itkSymmetricEigenDecomposition_h
 
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_vector.h"
#include "vnl/vnl_diag_matrix.h"
#include "itkEigenDecompositionSignConvention.h"
#include "itkEigenDecompositionSolverInfo.h"
#include "itkMacro.h"
#include "itk_eigen.h"
#include ITK_EIGEN(Dense)
 
namespace itk
{

template <typename T>
class SymmetricEigenDecomposition
{
public:
  explicit SymmetricEigenDecomposition(const vnl_matrix<T> & M, bool canonicalizeSigns = true)
    : V(M.rows(), M.cols())
    , D(M.rows())
  {
    using RowMajor = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
    using Vector = Eigen::Matrix<T, Eigen::Dynamic, 1>;
    const unsigned int n = M.rows();
 
    Eigen::Map<const RowMajor>                    mMap(M.data_block(), n, n);
    const Eigen::SelfAdjointEigenSolver<RowMajor> solver(mMap);
 
    if (solver.info() != Eigen::Success)
    {
      itkGenericExceptionMacro(<< "SymmetricEigenDecomposition: Eigen SelfAdjointEigenSolver failed: "
                               << detail::EigenComputationInfoString(solver.info()) << '.');
    }
 
    Eigen::Map<RowMajor>(V.data_block(), n, n) = solver.eigenvectors();
    Eigen::Map<Vector>(D.data_block(), n) = solver.eigenvalues();
 
    if (canonicalizeSigns)
    {
      detail::CanonicalizeEigenvectorColumnSigns(V);
    }
  }

  vnl_matrix<T> V;

  vnl_diag_matrix<T> D;

  vnl_vector<T>
  get_eigenvector(int i) const
  {
    return V.get_column(i);
  }

  T
  get_eigenvalue(int i) const
  {
    return D(i, i);
  }

  vnl_vector<T>
  nullvector() const
  {
    return V.get_column(0);
  }

  vnl_matrix<T>
  recompose() const
  {
    return V * D * V.transpose();
  }
};
 
} // namespace itk
 
#endif // itkSymmetricEigenDecomposition_h