itkMatrix.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 itkMatrix_h
#define itkMatrix_h

#include "itkPoint.h"
#include "itkCovariantVector.h"

#include <vxl_version.h>
#include "vnl/vnl_matrix_fixed.hxx" // Get the templates
#include "vnl/vnl_transpose.h"
#include "vnl/algo/vnl_matrix_inverse.h"
#include "vnl/vnl_matrix.h"
#include "vnl/algo/vnl_determinant.h"
#include "itkMath.h"

namespace itk
{
template< typename T, unsigned int NRows = 3, unsigned int NColumns = 3 >
class ITK_TEMPLATE_EXPORT Matrix
{
public:
  using Self = Matrix;

  using ValueType = T;
  using ComponentType = T;

  static constexpr unsigned int RowDimensions = NRows;
  static constexpr unsigned int ColumnDimensions = NColumns;

  using InternalMatrixType = vnl_matrix_fixed< T, NRows, NColumns >;

  using CompatibleSquareMatrixType = Matrix<T, NColumns, NColumns>;

  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;

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

  Self operator *(const CompatibleSquareMatrixType & matrix) const;

  template<unsigned int OuterDim>
   Matrix<T, NRows, OuterDim>  operator*(const  vnl_matrix_fixed< T, NRows, OuterDim > & matrix) const
      {
      const Matrix<T, NRows, OuterDim> result ( m_Matrix * matrix );
      return result;
      }

  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 CompatibleSquareMatrixType & 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) const
  {
    Self result(*this);

    result *= value;
    return result;
  }

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

  Self operator/(const T & value) const
  {
    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()
  {
    return m_Matrix;
  }

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

  inline void SetIdentity()
  {
    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 ( Math::NotExactlyEquals(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() const
  {
    if ( vnl_determinant(m_Matrix) == NumericTraits<T>::ZeroValue() )
      {
      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() const
  {
    return m_Matrix.transpose();
  }

  Matrix():m_Matrix(NumericTraits< T >::ZeroValue()) {}

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

  void swap(Self &other)
    {
      // allow for augment dependent look up
      using std::swap;
      swap(this->m_Matrix, other.m_Matrix);
    }
private:
  InternalMatrixType m_Matrix;
};

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


template< typename T, unsigned int NRows, unsigned int NColumns >
inline void swap( const Matrix< T, NRows, NColumns > &a, const Matrix< T, NRows, NColumns >&b )
{
  a.swap(b);
}

} // end namespace itk

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

#endif