itkSymmetricEigenAnalysis.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkSymmetricEigenAnalysis.h,v $
  Language:  C++
  Date:      $Date: 2009-03-03 15:09:43 $
  Version:   $Revision: 1.6 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#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, VectorType &d,
                                 double *e, double *e2) const;

  void ReduceToTridiagonalMatrixAndGetTransformation(
                  double *inputMatrix, VectorType &diagonalElements,
                  double *subDiagonalElements, double *transformMatrix) const;

  /* Finds the eigenvalues of a symmetric tridiagonal matrix by the ql method.
   *
   * On input:
   * 'd' contains the diagonal elements of the input matrix.
   * 'e' contains the subdiagonal elements of the input matrix
   * in its last n-1 positions.  e(1) is arbitrary.
   * On Output:
   * 'd' contains the eigenvalues.
   * 'e' has been destroyed.
   *
   * Returns:
   *          zero       for normal return,
   *          j          if the j-th eigenvalue has not been
   *                     determined after 30 iterations.
   *
   *
   * Reference
   *     this subroutine is a translation of the algol procedure tql1,
   *     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
   *     wilkinson.
   *     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
   *
   *     Questions and comments should be directed to burton s. garbow,
   *     mathematics and computer science div, argonne national laboratory
   *     this version dated august 1983.
   *
   *  Function Adapted from netlib/tql1.c.
   *  [Changed: remove static vars, enforce const correctness.
   *            Use vnl routines as necessary]                      */
  unsigned int ComputeEigenValuesUsingQL(
                         VectorType &d, double *e) const;


  /* Finds the eigenvalues and eigenvectors of a symmetric tridiagonal matrix
   * by the ql method.
   *
   * On input:
   * 'd' contains the diagonal elements of the input matrix.
   * 'e' contains the subdiagonal elements of the input matrix
   * in its last n-1 positions.  e(1) is arbitrary.
   * 'z' contains the transformation matrix produced in the reduction by
   * ReduceToTridiagonalMatrixAndGetTransformation(), if performed. If the
   * eigenvectors of the tridiagonal matrix are desired, z must contain
   * the identity matrix.

   * On Output:
   * 'd' contains the eigenvalues.
   * 'e' has been destroyed.
   * 'z' contains orthonormal eigenvectors of the symmetric tridiagonal
   * (or full) matrix.
   *
   * Returns:
   *          zero       for normal return,
   *          j          if the j-th eigenvalue has not been
   *                     determined after 1000 iterations.
   *
   * Reference
   *     this subroutine is a translation of the algol procedure tql1,
   *     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
   *     wilkinson.
   *     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
   *
   *     Questions and comments should be directed to burton s. garbow,
   *     mathematics and computer science div, argonne national laboratory
   *     this version dated august 1983.
   *
   *  Function Adapted from netlib/tql2.c.
   *  [Changed: remove static vars, enforce const correctness.
   *            Use vnl routines as necessary]
   */
  unsigned int ComputeEigenValuesAndVectorsUsingQL(
                                   VectorType &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.txx"
#endif

#endif