MeVis/Foundation/Sources/MLLinearAlgebra/mlMatrix4.h

Go to the documentation of this file.

// **InsertLicense** code 
//=====================================================================================

//=====================================================================================

#ifndef __mlMatrix4_H
#define __mlMatrix4_H

// Include system independency file and project settings.
#ifndef __mlLinearAlgebraSystem_H
#include "mlLinearAlgebraSystem.h"
#endif
#ifndef __mlLinearAlgebraDefs_H
#include "mlLinearAlgebraDefs.h"
#endif

#ifndef __mlFloatingPointMatrix_H
#include "mlFloatingPointMatrix.h"
#endif

#ifndef __mlMatrix3_H
#include "mlMatrix3.h"
#endif
#ifndef __mlVector4_H
#include "mlVector4.h"
#endif

// 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();
  
  // Build matrix that has the argument at the diagonal values, zero otherwise
  Tmat4(const DT diagValue);

  // Composes a matrix from the four row vectors row0, row1, row2 and row3.
  Tmat4(const Tvec4<DT> &row0, const Tvec4<DT> &row1, const Tvec4<DT> &row2, const Tvec4<DT> &row3);
  
  // Copy constructor from the Tmat4 mat.
  Tmat4(const Tmat4<DT> &mat);
  
  // Constructor from 16 floating point values in an array given by mat, row by row.
  Tmat4(const float mat[16]);
  
  // Constructor from 16 double values in an array given by mat, row by row.
  Tmat4(const double mat[16]);
  
  // Initialize all matrix elements explicitly with scalars,
  // filling it row by row. 
  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);
  
  // Construct matrix from three base vectors n0, ... n3 and a translation t,
  // all given as COLUMN vectors.
  Tmat4(const Tvec3<DT> &n0, const Tvec3<DT> &n1, const Tvec3<DT> &n3, const Tvec3<DT> &t);
  
  // Copy contents from float array mat into *this, row by row.
  void setValues(const float mat[16]);
  
  // Copy contents of *this into floating point matrix mat, row by row.
  // Note that range and precision of the float values may not be
  // sufficient for double or long double matrix contents.
  void getValues(float mat[16]) const;
  
  // Copy contents from double array mat into *this, row by row.
  void setValues(const double mat[16]);
  
  // Copy contents of *this into into double matrix mat, row by row.
  // Note that range and precision of the double values may not be
  // sufficient for long double matrix contents.
  void getValues(double mat[16]) const;
  
  // Set diagonal matrix with scale on diagonal.
  void setScaleMatrix(const DT scale);

  inline Tvec3<DT> transformPoint(const Tvec3<DT>& sourceVec) const
  {
    return Tvec3<DT>( (*this) * sourceVec.affinePoint(), true);
  }

#ifdef ML_DEPRECATED




  inline ML_DEPRECATED void TransformPoint(const Tvec3<DT>& sourceVec, Tvec3<DT>& targetVec) const
  {
    targetVec = transformPoint(sourceVec);
  }
#endif  
  
  //---------------------------------------------------------------
  // Assignment operators and other simple stuff.
  //---------------------------------------------------------------
  // Return a matrix filled with values val.
  static Tmat4<DT> getMat(const double val);
  
  // Set all values to val.
  void set(const double val);
  
  bool operator<(const Tmat4<DT> &) const { return false; }
  
  // Assignment from a Tmat4.
  const Tmat4<DT> &operator=(const Tmat4<DT> &m);
  
  // Incrementation by a Tmat4.
  const Tmat4<DT> &operator+=(const Tmat4<DT> &m);
  
  // Decrement by a Tmat4.
  const Tmat4<DT> &operator-=(const Tmat4<DT> &m);
  
  
  // Multiplication by a constant d.
  const Tmat4<DT> &operator*=(const DT d);
  
  // Division by a constant d. Division by zero is not handled and must be avoided by caller.
  const Tmat4<DT> &operator/=(const DT d);
  
  
  //---------------------------------------------------------------
  // Special functions
  //---------------------------------------------------------------
  // Determine the determinant of a 3x3 matrix given by A,B,C,D,E,F,G,H,I.
  DT det3(DT A, DT B, DT C, DT D, DT E, DT F, DT G, DT H, DT I) const;
  
  // Determine the (sub)determinant of columns given by col1, col2 and col3.
  DT determinantLower3(const int col1, const int col2, const int col3) const;
  
  // Returns the determinant of this.
  DT det() const;
  
  // Return the transposed *this.
  Tmat4<DT> transpose() const;
  
  // Returns the inverse. Gauss-Jordan elimination with partial pivoting.
  // If a non-NULL Boolean pointer is passed to isInvertible 
  // then true is returned in *isInvertible in the case of a 
  // successful inversion or false if the inversion is not possible
  // (function return is the identity then).
  // If a NULL pointer is passed as isInvertible the matrix must
  // be invertible, otherwise errors will occur.
  Tmat4<DT> inverse(bool* isInvertible=NULL) const;
  
  // Return identity matrix.
  static Tmat4<DT> getIdentity();
  
  // Apply the function fct to each component.
  const Tmat4<DT> &apply(MLDblFuncPtr fct);
  
}; // end of class *mat4*

//---------------------------------------------------------------

//---------------------------------------------------------------

template <class DT>
inline Tmat4<DT>::Tmat4()
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4<DT>::Tmat4()");
  
  this->v[0] = this->v[1] = this->v[2] = this->v[3] = Tvec4<DT>(0);
}

// Construct matrix that has the argument at the diagonal values, zero otherwise
template <class DT>
inline Tmat4<DT>::Tmat4(const DT diagValue)
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4<DT>::Tmat4(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)
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4<DT>::Tmat4(const Tvec4<DT> &row0, ..., 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)
{
  ML_TRACE_IN_TIME_CRITICAL("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])
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4<DT>::Tmat4(const float mat[16])");
  
  setValues(mat);
}

template <class DT>
inline Tmat4<DT>::Tmat4(const double mat[16])
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4<DT>::Tmat4(const double mat[16])");
  
  setValues(mat);
}

template <class DT>
inline Tmat4<DT> Tmat4<DT>::getMat(const double val)
{
  ML_TRACE_IN_TIME_CRITICAL("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 void Tmat4<DT>::set(const double val)
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4<DT>::set(const double 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)
{
  ML_TRACE_IN_TIME_CRITICAL("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)
{
  ML_TRACE_IN_TIME_CRITICAL("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)
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4<DT>::operator-=(const Tmat4 &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)
{
  ML_TRACE_IN_TIME_CRITICAL("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)
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4<DT>::operator/=(const DT d)");
  ML_TRY
  {
    this->v[0] /= d; 
    this->v[1] /= d; 
    this->v[2] /= d; 
    this->v[3] /= d;
  }
  ML_CATCH_RETHROW;
  return *this;
}

template <class DT>
inline const Tmat4<DT> &Tmat4<DT>::apply(MLDblFuncPtr fct)
{
  ML_TRACE_IN_TIME_CRITICAL("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)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: operator!=(const Tmat4<DT> &a, const Tmat4<DT> &b)");
  
  return !(a == b);
}


//---------------------------------------------------------------

//---------------------------------------------------------------
template <class DT>
inline Tmat4<DT> operator-(const Tmat4<DT> &a)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN("Tmat4<DT>::Tmat4(const Tvec3<DT> &n0, ..., 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)
{
  ML_TRACE_IN("Tmat4<DT>::Tmat4(const double in00, ..., 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])
{
  ML_TRACE_IN( "Tmat4::setValues( )" );
  ML_TRY
  {
    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];
  }
  ML_CATCH_RETHROW;
}

//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
void Tmat4<DT>::getValues(float mat[16]) const
{
  ML_TRACE_IN( "Tmat4::setValues( )" );
  ML_TRY
  {
    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]);
  }
  ML_CATCH_RETHROW;
}

//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
void Tmat4<DT>::setValues(const double mat[16])
{
  ML_TRACE_IN( "Tmat4::setValues( )" );
  ML_TRY
  {
    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]);
  }
  ML_CATCH_RETHROW;
}

//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
void Tmat4<DT>::getValues(double mat[16]) const
{
  ML_TRACE_IN( "Tmat4::setValues( )" );
  ML_TRY
  {
    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]);
  }
  ML_CATCH_RETHROW;
}

//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
void Tmat4<DT>::setScaleMatrix(const DT scale)
{
  ML_TRACE_IN( "Tmat4::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
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4::det3()");
  
  /* 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>
DT Tmat4<DT>::determinantLower3(const int col1, const int col2, const int col3) const
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4::determinantLower3()");
  
  /* 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]);
}

//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline DT Tmat4<DT>::det() const
{
  ML_TRACE_IN("Tmat4::det()");
  
  DT determ=0;
  ML_TRY
  {
    determ = (  this->v[0][0] * determinantLower3(1, 2, 3)
              - this->v[0][1] * determinantLower3(0, 2, 3)
              + this->v[0][2] * determinantLower3(0, 1, 3)
              - this->v[0][3] * determinantLower3(0, 1, 2));
  }
  ML_CATCH_RETHROW;
  return determ;
}

//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline Tmat4<DT> Tmat4<DT>::transpose() const
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4::transpose()");
  
  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()
{
  ML_TRACE_IN_TIME_CRITICAL("Tmat4::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()
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: identity3D()");
  
  return Tmat4<DT>::getIdentity();
}

//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline Tmat4<DT> translation3D(const Tvec3<DT> &v)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN("mlMatrix4.h: rotation3D(Tvec3<DT> Axis, const DT angleRad)");
  
  Tmat4<DT> retVal;
  ML_TRY
  {
    const DT c = cos(angleRad);
    const DT s = sin(angleRad);
    const DT t = 1.0 - c;
    
    Axis.normalize();
    retVal = 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));
  }
  ML_CATCH_RETHROW;
  return retVal;
}

//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
inline Tmat4<DT> scaling3D(const Tvec3<DT> &scaleVector)
{
  ML_TRACE_IN_TIME_CRITICAL("mlMatrix4.h: 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)
{
  ML_TRACE_IN( "perspective3D( )" );
  
  Tmat4<DT> retVal;
  ML_TRY
  {
    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;
  }
  ML_CATCH_RETHROW;
  return retVal;
}


//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class DT>
Tmat4<DT> Tmat4<DT>::inverse(bool* isInvertible) const
{
  ML_TRACE_IN( "Tmat4<DT>::inverse( )" );
  
  Tmat4<DT> retVal;
  ML_TRY
  {
    // 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.
    retVal = MLInverseMatHelper(*this,
                                isInvertible, 
                                Epsilon,
                                "Tmat4<DT> Tmat4<DT>::inverse(bool* isInvertible) const, matrix not invertable",
                                getIdentity(),
                                4);
  }
  ML_CATCH_RETHROW;
  return retVal;
}


//-----------------------------------------------------------------------------------

//-----------------------------------------------------------------------------------

typedef Tmat4<MLfloat>   Matrix4f;
typedef Tmat4<MLdouble>  Matrix4d;
typedef Tmat4<MLldouble> Matrix4ld;
typedef Tmat4<MLdouble>  Matrix4;


#ifdef ML_DEPRECATED



ML_DEPRECATED typedef Tmat4<MLfloat>   matf4;
ML_DEPRECATED typedef Tmat4<MLdouble>  matd4;
ML_DEPRECATED typedef Tmat4<MLldouble> matld4;
ML_DEPRECATED typedef Tmat4<MLdouble> mat4;


#endif // ML_DEPRECATED



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