itkMatrix.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkMatrix.h,v $
  Language:  C++
  Date:      $Date: 2009-02-07 22:38:33 $
  Version:   $Revision: 1.31 $

  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 __itkMatrix_h
#define __itkMatrix_h


#include "itkPoint.h"
#include "itkVector.h"
#include "itkCovariantVector.h"
#include "vnl/vnl_matrix_fixed.h"
#include "vnl/algo/vnl_matrix_inverse.h"
#include "vnl/vnl_transpose.h"
#include "vnl/vnl_matrix.h"
#include <vnl/algo/vnl_determinant.h>

namespace itk
{

template<class T, unsigned int NRows=3, unsigned int NColumns=3>
class Matrix
{
public:
  typedef Matrix  Self;

  typedef T ValueType;
  typedef T ComponentType;

  itkStaticConstMacro(RowDimensions,    unsigned int, NRows);
  itkStaticConstMacro(ColumnDimensions, unsigned int, NColumns);

  typedef vnl_matrix_fixed<T,NRows,NColumns>  InternalMatrixType;

  Vector<T,NRows> operator*(const Vector<T,NColumns> & vector) const;

  Point<T,NRows> operator*(const Point<T,NColumns> & vector) const;

  CovariantVector<T,NRows>
  operator*(const CovariantVector<T,NColumns> & vector) const;

  Self operator*(const Self & matrix) const;

  Self operator+(const Self & matrix) const;
  const Self & operator+=(const Self & matrix );

  Self operator-(const Self & matrix) const;
  const Self & operator-=(const Self & matrix );

  vnl_matrix<T> operator*(const vnl_matrix<T> & matrix) const;

  void operator*=(const Self & matrix);

  void operator*=(const vnl_matrix<T> & matrix);

  vnl_vector<T> operator*(const vnl_vector<T> & matrix) const;

  void operator*=(const T & value)
    { m_Matrix *= value; }

  Self operator*(const T & value)
    {
    Self result( *this );
    result *= value;
    return result;
    }

  void operator /= (const T & value)
    {
    m_Matrix /= value;
    }

  Self operator/(const T & value)
    {
    Self result( *this );
    result /= value;
    return result;
    }

  inline T & operator()( unsigned int row, unsigned int col )
    {
    return m_Matrix(row,col);
    }

  inline const T & operator()( unsigned int row, unsigned int col ) const
    { return m_Matrix(row,col); }

  inline T * operator[]( unsigned int i )
    { return m_Matrix[i]; }

  inline const T * operator[]( unsigned int i ) const
    { return m_Matrix[i]; }

  inline InternalMatrixType & GetVnlMatrix( void )
    { return m_Matrix; }

  inline const InternalMatrixType & GetVnlMatrix( void ) const
    { return m_Matrix; }

  inline void SetIdentity( void )
    { m_Matrix.set_identity(); }

  inline void Fill( const T & value )
    { m_Matrix.fill( value ); }

  inline const Self & operator=( const vnl_matrix<T> & matrix)
    {
    m_Matrix = matrix;
    return *this;
    }

  inline Matrix(const vnl_matrix<T> & matrix)
    {
    this->operator=(matrix);
    }

  inline bool operator==( const Self & matrix) const
    {
    bool equal = true;
    for( unsigned int r=0; r<NRows; r++)
      {
      for( unsigned int c=0; c<NColumns; c++ )
        {
        if (m_Matrix(r,c) != matrix.m_Matrix(r,c))
          {
          equal = false;
          break;
          }
        }
      }
    return equal;
    }

  inline bool operator!=( const Self & matrix) const
    {
    return !this->operator==(matrix);
    }

  inline const Self & operator=( const InternalMatrixType & matrix )
    {
    this->m_Matrix = matrix;
    return *this;
    }

  inline explicit Matrix(const InternalMatrixType & matrix)
    {
    this->operator=(matrix);
    }

  inline const Self & operator=( const Self & matrix)
    {
    m_Matrix = matrix.m_Matrix;
    return *this;
    }

  inline vnl_matrix_fixed<T,NColumns,NRows> GetInverse( void ) const
    {
    if (vnl_determinant(m_Matrix) == 0.0)
      {
      itkGenericExceptionMacro(<< "Singular matrix. Determinant is 0.");
      }
    vnl_matrix<T> temp = vnl_matrix_inverse<T>( m_Matrix );
    return temp;
    }

  inline vnl_matrix_fixed<T,NColumns,NRows> GetTranspose( void ) const
    {
    return m_Matrix.transpose();
    }

  Matrix() : m_Matrix(NumericTraits<T>::Zero) {};

  Matrix(const Self & matrix) : m_Matrix( matrix.m_Matrix ) {};

private:
  InternalMatrixType     m_Matrix;

};

template< class T, unsigned int NRows, unsigned int NColumns >
ITK_EXPORT std::ostream& operator<<( std::ostream& os, const Matrix<T,NRows,NColumns> & v)
{
  os << v.GetVnlMatrix();
  return os;
}


} // end namespace itk

// Define instantiation macro for this template.
#define ITK_TEMPLATE_Matrix(_, EXPORT, x, y) namespace itk { \
  _(3(class EXPORT Matrix< ITK_TEMPLATE_3 x >)) \
  namespace Templates { typedef Matrix< ITK_TEMPLATE_3 x > Matrix##y; } \
  }

#if ITK_TEMPLATE_EXPLICIT
# include "Templates/itkMatrix+-.h"
#endif

#if ITK_TEMPLATE_TXX
# include "itkMatrix.txx"
#endif

#endif