itkMatrix.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 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 "itkMathSVD.h"
#include "vnl/vnl_matrix.h"
 
// GetInverse is Eigen-backed via itk::Math::SVD, but ITK 5.4 exposed
// vnl_matrix_inverse (and through it vnl_svd / vnl_diag_matrix) transitively to
// every itkMatrix.h consumer. Retain that surface for release-5.4 source
// compatibility in default builds; ITK proper includes what it uses, so
// ITK_LEGACY_REMOVE drops the leak. The transitional guard keeps the
// vnl_matrix_inverse deprecation warning off itkMatrix consumers.
#if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE)
#  define ITK_VNL_SVD_TRANSITIONAL_INCLUDE
#  include "vnl/algo/vnl_matrix_inverse.h"
#  undef ITK_VNL_SVD_TRANSITIONAL_INCLUDE
#endif
#include "itkMathDeterminant.h"
// GetInverse's singular check is Eigen-backed via itk::Math::Determinant, but
// ITK 5.4 exposed vnl_determinant transitively through this header; keep it
// reachable during the deprecation window so downstream code still compiles.
#if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE)
#  include "vnl/algo/vnl_determinant.h"
#endif
#include "itkMath.h"
#include <type_traits> // For is_same.
 
namespace itk
{
 
template <typename T, unsigned int VRows = 3, unsigned int VColumns = 3>
class ITK_TEMPLATE_EXPORT Matrix
{
public:
  using Self = Matrix;

  using ValueType = T;
  using ComponentType = T;

  using pointer = ValueType *;

  using const_pointer = const ValueType *;

  using reference = ValueType &;

  using const_reference = const ValueType &;

  using iterator = ValueType *;

  using const_iterator = const ValueType *;

  static constexpr unsigned int RowDimensions = VRows;
  static constexpr unsigned int ColumnDimensions = VColumns;

  using InternalMatrixType = vnl_matrix_fixed<T, VRows, VColumns>;

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

  Vector<T, VRows>
  operator*(const Vector<T, VColumns> & vect) const;

  Point<T, VRows>
  operator*(const Point<T, VColumns> & pnt) const;

  CovariantVector<T, VRows>
  operator*(const CovariantVector<T, VColumns> & covect) const;

  vnl_vector_fixed<T, VRows>
  operator*(const vnl_vector_fixed<T, VColumns> & inVNLvect) const;

  Self
  operator*(const CompatibleSquareMatrixType & matrix) const;
 
  template <unsigned int OuterDim>
  Matrix<T, VRows, OuterDim>
  operator*(const vnl_matrix_fixed<T, VRows, OuterDim> & matrix) const
  {
    const Matrix<T, VRows, 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> & vc) 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;
  }

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

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

  static Self
  GetIdentity()
  {
    InternalMatrixType internalMatrix;
    internalMatrix.set_identity();
    return Self{ internalMatrix };
  }

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

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

  inline explicit Matrix(const vnl_matrix<T> & matrix)
    : m_Matrix(matrix)
  {}

  template <typename TElement>
  explicit Matrix(const TElement (&elements)[VRows][VColumns])
    : m_Matrix(&elements[0][0])
  {
    static_assert(std::is_same_v<TElement, T>,
                  "The type of an element should correspond with this itk::Matrix instantiation.");
  }

  inline bool
  operator==(const Self & matrix) const
  {
    bool equal = true;
 
    for (unsigned int r = 0; r < VRows; ++r)
    {
      for (unsigned int c = 0; c < VColumns; ++c)
      {
        if (Math::NotExactlyEquals(m_Matrix(r, c), matrix.m_Matrix(r, c)))
        {
          equal = false;
          break;
        }
      }
    }
    return equal;
  }
 
  ITK_UNEQUAL_OPERATOR_MEMBER_FUNCTION(Self);
 
 
  inline Self &
  operator=(const InternalMatrixType & matrix)
  {
    this->m_Matrix = matrix;
    return *this;
  }

  inline Matrix(const InternalMatrixType & matrix)
    : m_Matrix(matrix)
  {}

  [[nodiscard]] inline vnl_matrix_fixed<T, VColumns, VRows>
  GetInverse() const
  {
    if (Math::Determinant(m_Matrix) == T{})
    {
      itkGenericExceptionMacro("Singular matrix. Determinant is 0.");
    }
    return vnl_matrix_fixed<T, VColumns, VRows>{ Math::SVD(m_Matrix.as_ref()).PseudoInverse() };
  }

  [[nodiscard]] inline vnl_matrix_fixed<T, VColumns, VRows>
  GetTranspose() const
  {
    return vnl_matrix_fixed<T, VColumns, VRows>{ m_Matrix.transpose().as_matrix() };
  }

  Matrix() = default;

  [[nodiscard]] constexpr unsigned int
  size() const
  {
    return m_Matrix.size();
  }

  [[nodiscard]] iterator
  begin()
  {
    return m_Matrix.begin();
  }

  [[nodiscard]] iterator
  end()
  {
    return m_Matrix.end();
  }

  [[nodiscard]] const_iterator
  begin() const
  {
    return m_Matrix.begin();
  }

  [[nodiscard]] const_iterator
  end() const
  {
    return m_Matrix.end();
  }

  [[nodiscard]] const_iterator
  cbegin() const
  {
    return m_Matrix.begin();
  }

  [[nodiscard]] const_iterator
  cend() const
  {
    return m_Matrix.end();
  }
 
  void
  swap(Self & other) noexcept
  {
    // allow for augment dependent look up
    using std::swap;
    swap(this->m_Matrix, other.m_Matrix);
  }
 
  inline void
  PrintSelf(std::ostream & os, Indent indent) const
  {
    os << indent << "Matrix (" << VRows << "x" << VColumns << ")\n";
    for (unsigned int r = 0; r < VRows; ++r)
    {
      os << indent << "  ";
      for (unsigned int c = 0; c < VColumns; ++c)
      {
        os << m_Matrix(r, c) << " ";
      }
      os << "\n";
    }
  }
 
private:
  InternalMatrixType m_Matrix{};
};
 
template <typename T, unsigned int VRows, unsigned int VColumns>
std::ostream &
operator<<(std::ostream & os, const Matrix<T, VRows, VColumns> & v)
{
  os << v.GetVnlMatrix();
  return os;
}
 
 
template <typename T, unsigned int VRows, unsigned int VColumns>
inline void
swap(const Matrix<T, VRows, VColumns> & a, const Matrix<T, VRows, VColumns> & b) noexcept
{
  a.swap(b);
}
 
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkMatrix.hxx"
#endif
 
#endif