mlMatrix4.h

Go to the documentation of this file.

/*************************************************************************************
**
** Copyright 2007, MeVis Medical Solutions AG
**
** The user may use this file in accordance with the license agreement provided with
** the Software or, alternatively, in accordance with the terms contained in a
** written agreement between the user and MeVis Medical Solutions AG.
**
** For further information use the contact form at https://www.mevislab.de/contact
**
**************************************************************************************/
 
#ifndef ML_MATRIX4_H
#define ML_MATRIX4_H
 
 
// Include system independent file and project settings.
#include "mlLinearAlgebraSystem.h"
#include "mlLinearAlgebraDefs.h"
 
#include "mlFloatingPointMatrix.h"
 
#include "mlMatrix3.h"
#include "mlVector4.h"
#include "mlVector3.h"
 
// All declarations of this header will be in the ML_LA_START_NAMESPACE namespace.
ML_LA_START_NAMESPACE
 
//---------------------------------------------------------------
//---------------------------------------------------------------
template <class DT>
class Tmat4 : public FloatingPointMatrix<Tvec4<DT>, 4>
{
public:
 
  typedef DT ComponentType;
 
  //---------------------------------------------------------------
  // Constructors, set, and get methods.
  //---------------------------------------------------------------
  // Construct matrix from 16 zero elements.
  Tmat4();
 
  Tmat4(const DT diagValue);
 
  Tmat4(const Tvec4<DT> &row0, const Tvec4<DT> &row1, const Tvec4<DT> &row2, const Tvec4<DT> &row3);
 
  Tmat4(const Tmat4<DT> &mat);
 
  Tmat4(const float mat[16]);
 
  Tmat4(const double mat[16]);
 
  Tmat4(const double in00, const double in01, const double in02, const double in03,
        const double in10, const double in11, const double in12, const double in13,
        const double in20, const double in21, const double in22, const double in23,
        const double in30, const double in31, const double in32, const double in33);
 
  Tmat4(const Tvec3<DT> &n0, const Tvec3<DT> &n1, const Tvec3<DT> &n3, const Tvec3<DT> &t);
 
  void setValues(const float mat[16]);
 
  void getValues(float mat[16]) const;
 
  void setValues(const double mat[16]);
 
  void getValues(double mat[16]) const;
 
  void setScaleMatrix(const DT scale);
 
  inline Tvec3<DT> transformPoint(const Tvec3<DT>& sourceVec) const
  {
    return Tvec3<DT>( (*this) * sourceVec.affinePoint(), true);
  }
 
  Tmat4<DT> getRotationMatrix() const;
 
  //---------------------------------------------------------------
  // Assignment operators and other simple methods.
  //---------------------------------------------------------------
  static Tmat4<DT> getMat(const double val);
 
  void set(const DT val);
 
  bool operator<(const Tmat4<DT> &) const { return false; }
 
  const Tmat4<DT> &operator=(const Tmat4<DT> &m);
 
  const Tmat4<DT> &operator+=(const Tmat4<DT> &m);
 
  const Tmat4<DT> &operator-=(const Tmat4<DT> &m);
 
 
  const Tmat4<DT> &operator*=(const DT d);
 
  const Tmat4<DT> &operator/=(const DT d);
 
 
  //---------------------------------------------------------------
  // Special functions
  //---------------------------------------------------------------
  DT det3(DT A, DT B, DT C, DT D, DT E, DT F, DT G, DT H, DT I) const;
 
  template <class IDX_TYP>
  inline DT determinantLower3(const IDX_TYP col1, const IDX_TYP col2, const IDX_TYP col3) const
  {
    /* Determinant of 3x3 matrix with given entries */
    return det3(this->v[1][col1], this->v[1][col2], this->v[1][col3],
                this->v[2][col1], this->v[2][col2], this->v[2][col3],
                this->v[3][col1], this->v[3][col2], this->v[3][col3]);
  }
 
  DT det() const;
 
  Tmat4<DT> transpose() const;
 
  Tmat4<DT> inverse(bool* isInvertible=nullptr) const;
 
  static Tmat4<DT> getIdentity();
 
  const Tmat4<DT> &apply(MLDblFuncPtr fct);
 
}; // end of class *mat4*
 
 
//---------------------------------------------------------------
//---------------------------------------------------------------
 
template <class DT>
inline Tmat4<DT>::Tmat4()
{
  this->v[0] = this->v[1] = this->v[2] = this->v[3] = Tvec4<DT>(0);
}
 
template <class DT>
inline Tmat4<DT>::Tmat4(const DT diagValue)
{
  this->v[0][0] = this->v[1][1] = this->v[2][2] = this->v[3][3] = diagValue;
  this->v[1][0] = this->v[2][0] = this->v[3][0] = this->v[2][1] = this->v[3][1] = this->v[3][2] = 0;
  this->v[0][1] = this->v[0][2] = this->v[0][3] = this->v[1][2] = this->v[1][3] = this->v[2][3] = 0;
}
 
template <class DT>
inline Tmat4<DT>::Tmat4(const Tvec4<DT> &row0, const Tvec4<DT> &row1, const Tvec4<DT> &row2, const Tvec4<DT> &row3)
{
  this->v[0] = row0;
  this->v[1] = row1;
  this->v[2] = row2;
  this->v[3] = row3;
}
 
template <class DT>
inline Tmat4<DT>::Tmat4(const Tmat4<DT> &mat)
{
  this->v[0] = mat.v[0];
  this->v[1] = mat.v[1];
  this->v[2] = mat.v[2];
  this->v[3] = mat.v[3];
}
 
template <class DT>
inline Tmat4<DT>::Tmat4(const float mat[16])
{
  setValues(mat);
}
 
template <class DT>
inline Tmat4<DT>::Tmat4(const double mat[16])
{
  setValues(mat);
}
 
template <class DT>
inline Tmat4<DT> Tmat4<DT>::getMat(const double val)
{
  return Tmat4<DT>(val, val, val, val,
                   val, val, val, val,
                   val, val, val, val,
                   val, val, val, val);
}
 
template <class DT>
inline Tmat4<DT> Tmat4<DT>::getRotationMatrix() const
{
  Tmat4<DT> matrix = *this;
 
  // delete translation vector
  matrix[0][3] = matrix[1][3] = matrix[2][3] = 0;
 
  // normalize direction vectors
  Tvec3<DT> dirVec;
  int i, j;
  for (i = 0; i < 3; i++)
  {
    for (j = 0; j < 3; j++)
    {
      dirVec[j] = matrix[j][i];
    }
    dirVec.normalize();
    for (j = 0; j < 3; j++)
    {
      matrix[j][i] = dirVec[j];
    }
  }
  return matrix;
}
 
template <class DT>
inline void Tmat4<DT>::set(const DT val)
{
  this->v[0] = this->v[1] = this->v[2] = this->v[3] = Tvec4<DT>(val);
}
 
template <class DT>
inline const Tmat4<DT> &Tmat4<DT>::operator=(const Tmat4<DT> &m)
{
  if (&m != this){
    this->v[0] = m.v[0];
    this->v[1] = m.v[1];
    this->v[2] = m.v[2];
    this->v[3] = m.v[3];
  }
  return *this;
}
 
template <class DT>
inline const Tmat4<DT> &Tmat4<DT>::operator+=(const Tmat4<DT> &m)
{
  this->v[0] += m.v[0];
  this->v[1] += m.v[1];
  this->v[2] += m.v[2];
  this->v[3] += m.v[3];
  return *this;
}
 
template <class DT>
inline const Tmat4<DT> &Tmat4<DT>::operator-=(const Tmat4<DT> &m)
{
  this->v[0] -= m.v[0];
  this->v[1] -= m.v[1];
  this->v[2] -= m.v[2];
  this->v[3] -= m.v[3];
  return *this;
}
 
template <class DT>
inline const Tmat4<DT> &Tmat4<DT>::operator*=(const DT d)
{
  this->v[0] *= d;
  this->v[1] *= d;
  this->v[2] *= d;
  this->v[3] *= d;
  return *this;
}
 
template <class DT>
inline const Tmat4<DT> &Tmat4<DT>::operator/=(const DT d)
{
  this->v[0] /= d;
  this->v[1] /= d;
  this->v[2] /= d;
  this->v[3] /= d;
  return *this;
}
 
template <class DT>
inline const Tmat4<DT> &Tmat4<DT>::apply(MLDblFuncPtr fct)
{
  this->v[0].apply(fct);
  this->v[1].apply(fct);
  this->v[2].apply(fct);
  this->v[3].apply(fct);
  return *this;
}
 
//---------------------------------------------------------------
//---------------------------------------------------------------
#define _ML_MAT4_RC(i, j) a[i][0]*b[0][j] + a[i][1]*b[1][j] + \
a[i][2]*b[2][j] + a[i][3]*b[3][j]
 
template <class DT>
inline Tmat4<DT> operator*(const Tmat4<DT> &a, const Tmat4<DT> &b)
{
  return Tmat4<DT>(
                   Tvec4<DT>(_ML_MAT4_RC(0,0), _ML_MAT4_RC(0,1), _ML_MAT4_RC(0,2), _ML_MAT4_RC(0,3)),
                   Tvec4<DT>(_ML_MAT4_RC(1,0), _ML_MAT4_RC(1,1), _ML_MAT4_RC(1,2), _ML_MAT4_RC(1,3)),
                   Tvec4<DT>(_ML_MAT4_RC(2,0), _ML_MAT4_RC(2,1), _ML_MAT4_RC(2,2), _ML_MAT4_RC(2,3)),
                   Tvec4<DT>(_ML_MAT4_RC(3,0), _ML_MAT4_RC(3,1), _ML_MAT4_RC(3,2), _ML_MAT4_RC(3,3))
                   );
}
#undef _ML_MAT4_RC
 
 
template <class DT>
inline bool operator==(const Tmat4<DT> &a, const Tmat4<DT> &b)
{
  return ((a[0] == b[0]) &&
          (a[1] == b[1]) &&
          (a[2] == b[2]) &&
          (a[3] == b[3]));
}
 
template <class DT>
inline bool operator!=(const Tmat4<DT>& a, const Tmat4<DT>& b)
{
  return !(a == b);
}
 
 
//---------------------------------------------------------------
//---------------------------------------------------------------
template <class DT>
inline Tmat4<DT> operator-(const Tmat4<DT> &a)
{
  return Tmat4<DT>(a) *= static_cast<DT>(-1.0);
}
 
template <class DT>
inline Tmat4<DT> operator+(const Tmat4<DT> &a, const Tmat4<DT> &b)
{
  return Tmat4<DT>(a) += b;
}
 
template <class DT>
inline Tmat4<DT> operator-(const Tmat4<DT> &a, const Tmat4<DT> &b)
{
  return Tmat4<DT>(a) -= b;
}
 
template <class DT>
inline Tmat4<DT> operator*(const Tmat4<DT> &a, const DT d)
{
  return Tmat4<DT>(a) *= d;
}
 
template <class DT>
inline Tmat4<DT> operator*(const DT d, const Tmat4<DT> &a)
{
  return Tmat4<DT>(a) *= d;
}
 
template <class DT>
inline Tmat4<DT> operator/(const Tmat4<DT> &a, const DT d)
{
  return Tmat4<DT>(a) /= d;
}
 
 
//--------------------------------------------------------------------
//
//
//--------------------------------------------------------------------
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline Tmat4<DT>::Tmat4(const Tvec3<DT> &n0, const Tvec3<DT> &n1, const Tvec3<DT> &n2, const Tvec3<DT> &t)
{
  this->v[0][0]=n0[0];
  this->v[1][0]=n0[1];
  this->v[2][0]=n0[2];
 
  this->v[0][1]=n1[0];
  this->v[1][1]=n1[1];
  this->v[2][1]=n1[2];
 
  this->v[0][2]=n2[0];
  this->v[1][2]=n2[1];
  this->v[2][2]=n2[2];
 
  this->v[0][3]=t[0];
  this->v[1][3]=t[1];
  this->v[2][3]=t[2];
 
  this->v[3][0]=0;
  this->v[3][1]=0;
  this->v[3][2]=0;
  this->v[3][3]=1;
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
Tmat4<DT>::Tmat4(const double in00, const double in01, const double in02, const double in03,
                 const double in10, const double in11, const double in12, const double in13,
                 const double in20, const double in21, const double in22, const double in23,
                 const double in30, const double in31, const double in32, const double in33)
{
  this->v[0][0]=static_cast<DT>(in00); this->v[0][1]=static_cast<DT>(in01); this->v[0][2]=static_cast<DT>(in02); this->v[0][3]=static_cast<DT>(in03);
  this->v[1][0]=static_cast<DT>(in10); this->v[1][1]=static_cast<DT>(in11); this->v[1][2]=static_cast<DT>(in12); this->v[1][3]=static_cast<DT>(in13);
  this->v[2][0]=static_cast<DT>(in20); this->v[2][1]=static_cast<DT>(in21); this->v[2][2]=static_cast<DT>(in22); this->v[2][3]=static_cast<DT>(in23);
  this->v[3][0]=static_cast<DT>(in30); this->v[3][1]=static_cast<DT>(in31); this->v[3][2]=static_cast<DT>(in32); this->v[3][3]=static_cast<DT>(in33);
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
void Tmat4<DT>::setValues(const float mat[16])
{
  this->v[0][0] = mat[ 0]; this->v[0][1] = mat[ 1]; this->v[0][2] = mat[ 2]; this->v[0][3] = mat[ 3];
  this->v[1][0] = mat[ 4]; this->v[1][1] = mat[ 5]; this->v[1][2] = mat[ 6]; this->v[1][3] = mat[ 7];
  this->v[2][0] = mat[ 8]; this->v[2][1] = mat[ 9]; this->v[2][2] = mat[10]; this->v[2][3] = mat[11];
  this->v[3][0] = mat[12]; this->v[3][1] = mat[13]; this->v[3][2] = mat[14]; this->v[3][3] = mat[15];
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
void Tmat4<DT>::getValues(float mat[16]) const
{
  mat[ 0] = static_cast<float>(this->v[0][0]); mat[ 1] = static_cast<float>(this->v[0][1]); mat[ 2] = static_cast<float>(this->v[0][2]); mat[ 3] = static_cast<float>(this->v[0][3]);
  mat[ 4] = static_cast<float>(this->v[1][0]); mat[ 5] = static_cast<float>(this->v[1][1]); mat[ 6] = static_cast<float>(this->v[1][2]); mat[ 7] = static_cast<float>(this->v[1][3]);
  mat[ 8] = static_cast<float>(this->v[2][0]); mat[ 9] = static_cast<float>(this->v[2][1]); mat[10] = static_cast<float>(this->v[2][2]); mat[11] = static_cast<float>(this->v[2][3]);
  mat[12] = static_cast<float>(this->v[3][0]); mat[13] = static_cast<float>(this->v[3][1]); mat[14] = static_cast<float>(this->v[3][2]); mat[15] = static_cast<float>(this->v[3][3]);
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
void Tmat4<DT>::setValues(const double mat[16])
{
  this->v[0][0] = static_cast<DT>(mat[ 0]); this->v[0][1] = static_cast<DT>(mat[ 1]); this->v[0][2] = static_cast<DT>(mat[ 2]); this->v[0][3] = static_cast<DT>(mat[ 3]);
  this->v[1][0] = static_cast<DT>(mat[ 4]); this->v[1][1] = static_cast<DT>(mat[ 5]); this->v[1][2] = static_cast<DT>(mat[ 6]); this->v[1][3] = static_cast<DT>(mat[ 7]);
  this->v[2][0] = static_cast<DT>(mat[ 8]); this->v[2][1] = static_cast<DT>(mat[ 9]); this->v[2][2] = static_cast<DT>(mat[10]); this->v[2][3] = static_cast<DT>(mat[11]);
  this->v[3][0] = static_cast<DT>(mat[12]); this->v[3][1] = static_cast<DT>(mat[13]); this->v[3][2] = static_cast<DT>(mat[14]); this->v[3][3] = static_cast<DT>(mat[15]);
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
void Tmat4<DT>::getValues(double mat[16]) const
{
  mat[ 0] = static_cast<double>(this->v[0][0]); mat[ 1] = static_cast<double>(this->v[0][1]); mat[ 2] = static_cast<double>(this->v[0][2]); mat[ 3] = static_cast<double>(this->v[0][3]);
  mat[ 4] = static_cast<double>(this->v[1][0]); mat[ 5] = static_cast<double>(this->v[1][1]); mat[ 6] = static_cast<double>(this->v[1][2]); mat[ 7] = static_cast<double>(this->v[1][3]);
  mat[ 8] = static_cast<double>(this->v[2][0]); mat[ 9] = static_cast<double>(this->v[2][1]); mat[10] = static_cast<double>(this->v[2][2]); mat[11] = static_cast<double>(this->v[2][3]);
  mat[12] = static_cast<double>(this->v[3][0]); mat[13] = static_cast<double>(this->v[3][1]); mat[14] = static_cast<double>(this->v[3][2]); mat[15] = static_cast<double>(this->v[3][3]);
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
void Tmat4<DT>::setScaleMatrix(const DT scale)
{
  this->v[0][0] = scale; this->v[0][1] = 0;     this->v[0][2] = 0;     this->v[0][3] = 0;
  this->v[1][0] = 0;     this->v[1][1] = scale; this->v[1][2] = 0;     this->v[1][3] = 0;
  this->v[2][0] = 0;     this->v[2][1] = 0;     this->v[2][2] = scale; this->v[2][3] = 0;
  this->v[3][0] = 0;     this->v[3][1] = 0;     this->v[3][2] = 0;     this->v[3][3] = scale;
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline DT Tmat4<DT>::det3(const DT A, const DT B, const DT C, const DT D, const DT E, const DT F, const DT G, const DT H, const DT I) const
{
  /* Determinant of 3x3 matrix with given entries */
  return ((A*E*I + B*F*G + C*D*H) - (A*F*H + B*D*I + C*E*G));
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline DT Tmat4<DT>::det() const
{
  return (  this->v[0][0] * determinantLower3(1u, 2u, 3u)
          - this->v[0][1] * determinantLower3(0u, 2u, 3u)
          + this->v[0][2] * determinantLower3(0u, 1u, 3u)
          - this->v[0][3] * determinantLower3(0u, 1u, 2u));
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline Tmat4<DT> Tmat4<DT>::transpose() const
{
  return Tmat4<DT>(this->v[0][0], this->v[1][0], this->v[2][0], this->v[3][0],
                   this->v[0][1], this->v[1][1], this->v[2][1], this->v[3][1],
                   this->v[0][2], this->v[1][2], this->v[2][2], this->v[3][2],
                   this->v[0][3], this->v[1][3], this->v[2][3], this->v[3][3]);
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline Tmat4<DT> Tmat4<DT>::getIdentity()
{
  return Tmat4<DT>(1, 0, 0, 0,
                   0, 1, 0, 0,
                   0, 0, 1, 0,
                   0, 0, 0, 1);
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
Tmat4<DT> identity3D()
{
  return Tmat4<DT>::getIdentity();
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline Tmat4<DT> translation3D(const Tvec3<DT> &v)
{
  return Tmat4<DT>(Tvec4<DT>(1.0, 0.0, 0.0, v[0]),
                   Tvec4<DT>(0.0, 1.0, 0.0, v[1]),
                   Tvec4<DT>(0.0, 0.0, 1.0, v[2]),
                   Tvec4<DT>(0.0, 0.0, 0.0, 1.0));
}
 
//--------------------------------------------------------------------
//-------------------------------------------------------------------
template <class DT>
Tmat4<DT> rotation3D(Tvec3<DT> Axis, const DT angleRad)
{
  const DT c = cos(angleRad);
  const DT s = sin(angleRad);
  const DT t = 1.0 - c;
 
  Axis.normalize();
  return Tmat4<DT>(Tvec4<DT>(t * Axis[0] * Axis[0] + c,
                             t * Axis[0] * Axis[1] - s * Axis[2],
                             t * Axis[0] * Axis[2] + s * Axis[1],
                             0.0),
                   Tvec4<DT>(t * Axis[0] * Axis[1] + s * Axis[2],
                             t * Axis[1] * Axis[1] + c,
                             t * Axis[1] * Axis[2] - s * Axis[0],
                             0.0),
                   Tvec4<DT>(t * Axis[0] * Axis[2] - s * Axis[1],
                             t * Axis[1] * Axis[2] + s * Axis[0],
                             t * Axis[2] * Axis[2] + c,
                             0.0),
                   Tvec4<DT>(0.0, 0.0, 0.0, 1.0));
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline Tmat4<DT> scaling3D(const Tvec3<DT> &scaleVector)
{
  return Tmat4<DT>(Tvec4<DT>(scaleVector[0], 0.0, 0.0, 0.0),
                   Tvec4<DT>(0.0, scaleVector[1], 0.0, 0.0),
                   Tvec4<DT>(0.0, 0.0, scaleVector[2], 0.0),
                   Tvec4<DT>(0.0, 0.0, 0.0, 1.0));
}
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline Tmat4<DT> perspective3D(const DT d)
{
  Tmat4<DT> retVal;
  ML_CHECK_FLOAT_THROW(d);
 
  retVal[0][0] = 1.0;
  retVal[1][1] = 1.0;
  retVal[2][2] = 1.0;
  retVal[3][2] = 1.0/d;
  retVal[3][3] = 1.0;
  return retVal;
}
 
 
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
Tmat4<DT> Tmat4<DT>::inverse(bool* isInvertible) const
{
  // Epsilon for comparison with 0 in inversion process.
  static const DT Epsilon = static_cast<DT>(10e-14);
 
  // Use helper function from tools to invert the matrix.
  return MLInverseMatHelper(*this,
                            isInvertible,
                            Epsilon,
                            "Tmat4<DT> Tmat4<DT>::inverse(bool* isInvertible) const, matrix not invertable",
                            getIdentity(),
                            4);
}
 
 
//-----------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------
 
typedef Tmat4<MLfloat>   Matrix4f;
typedef Tmat4<MLdouble>  Matrix4d;
typedef Tmat4<MLldouble> Matrix4ld;
typedef Tmat4<MLdouble>  Matrix4;
 
 
ML_LA_END_NAMESPACE
 
namespace std
{
  //-------------------------------------------------------------------
  //-------------------------------------------------------------------
  template <class DT>
  inline std::ostream &operator<<(std::ostream &os, const ML_LA_NAMESPACE::Tmat4<DT> &m)
  {
    return os << m[0] << '\n' << m[1] << '\n' << m[2] << '\n' << m[3];
  }
 
  // Input from stream.
  template <class DT>
  inline std::istream &operator>>(std::istream &is, ML_LA_NAMESPACE::Tmat4<DT> &m)
  {
    ML_LA_NAMESPACE::Tmat4<DT> m_tmp;
    is >> m_tmp[0] >> m_tmp[1] >> m_tmp[2] >> m_tmp[3];
    if (is){ m = m_tmp; }
    return is;
  }
}
 
#endif // __mlMatrix4_H