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"

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

  SymmetricEigenAnalysis():
  m_Dimension(0),
  m_Order(0),
  m_OrderEigenValues(OrderByValue) {}

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

  ~SymmetricEigenAnalysis() {}

  typedef TMatrix      MatrixType;
  typedef TEigenMatrix EigenMatrixType;
  typedef TVector      VectorType;

  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; }

private:
  unsigned int        m_Dimension;
  unsigned int        m_Order;
  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;
};

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;
  return os;
}
} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkSymmetricEigenAnalysis.hxx"
#endif

#endif