albaMatrix3x3.h

Go to the documentation of this file.

/*=========================================================================
 
 Program: ALBA (Agile Library for Biomedical Applications)
 Module: albaMatrix3x3
 Authors: Based on vtkMath code (www.vtk.org), adapted by Marco Petrone
 
 Copyright (c) BIC
 All rights reserved. See Copyright.txt or
 
 
 This software is distributed WITHOUT ANY WARRANTY; without even
 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 PURPOSE.  See the above copyright notice for more information.
 
=========================================================================*/
 
#ifndef __albaMatrix3x3_h
#define __albaMatrix3x3_h
 
#include "albaObject.h"
#include "albaTimeStamped.h"
#include <math.h>
 
typedef double (*albaMatrix3x3Elements)[3];
 
 
class ALBA_EXPORT albaMatrix3x3: public albaObject, public albaTimeStamped
{
public:
  albaTypeMacro(albaMatrix3x3,albaObject);
  virtual void Print (std::ostream& os, const int indent=0) const;
 
  albaMatrix3x3();
  virtual ~albaMatrix3x3();
 
  albaMatrix3x3 &operator=(const albaMatrix3x3 &mat);
  albaMatrix3x3(albaMatrix3x3 &mat);
 
  void DeepCopy(albaMatrix3x3 *mat);
 
  static void MultiplyVector(const float A[3][3], const float in[3], 
                          float out[3]);
  static void MultiplyVector(const double A[3][3], const double in[3], 
                          double out[3]);
 
  void MultiplyVector(const double in[3], double out[3]) {MultiplyVector(GetElements(),in,out);}
 
  static void Multiply(const double A[3][3], const double B[3][3], 
                          double C[3][3]);
 
  static void Multiply(const albaMatrix3x3 &A,albaMatrix3x3 &B, albaMatrix3x3 &C) \
    {Multiply(A.GetElements(),B.GetElements(),C.GetElements());C.Modified();}
 
  albaMatrix3x3Elements GetElements() const {return (albaMatrix3x3Elements)m_Elements;}
 
  void SetElement(int i, int j, double value) {GetElements()[i][j]=value;}
 
  double GetElement(int i, int j) const {return GetElements()[i][j];}
 
  static void GetVersor(int axis, const albaMatrix3x3 &matrix, double versor[3]);
  void GetVersor(int axis, double versor[3]) {GetVersor(axis,*this,versor);}
 
  void Zero() { albaMatrix3x3::Zero(*GetElements()); Modified(); }
  static void Zero(double elements[9]);
 
  void Identity() { albaMatrix3x3::Identity(GetElements()); Modified();}
  static void Identity(double A[3][3]);
 
  static void Transpose(const double A[3][3], double AT[3][3]);
  void Transpose() {Transpose(GetElements(),GetElements());; Modified();}
 
  static void Invert(const double A[3][3], double AI[3][3]);
  void Invert() {albaMatrix3x3::Invert(GetElements(),GetElements());; Modified();}
  
  static void Orthogonalize(const double A[3][3], double B[3][3]);
  
  void Orthogonalize() {Orthogonalize(GetElements(),GetElements());Modified();}
 
  static void Diagonalize(const double A[3][3],double w[3],double V[3][3]);
 
  static double Determinant(double A[3][3]);
 
  double Determinant() {return Determinant(GetElements());}
 
  static inline double Determinant(const double c1[3], 
                                      const double c2[3], 
                                      const double c3[3]);
 
  static inline double Determinant(double a1, double a2, double a3, 
                                      double b1, double b2, double b3, 
                                      double c1, double c2, double c3);
  
  static void QuaternionToMatrix(const double quat[4], double A[3][3]); 
  
  void QuaternionToMatrix(const double quat[4]) {QuaternionToMatrix(quat,GetElements());}
 
  static void MatrixToQuaternion(const double A[3][3], double quat[4]);
  
  void MatrixToQuaternion(double quat[4]) {MatrixToQuaternion(GetElements(),quat);}
  
  static void SingularValueDecomposition(const double A[3][3],
                                            double U[3][3], double w[3],
                                            double VT[3][3]);
                                            
  void SingularValueDecomposition(double U[3][3], double w[3], double VT[3][3]) \
    {SingularValueDecomposition(GetElements(),U,w,VT);}
 
  /*
    Jacobi iteration for the solution of eigenvectors/eigenvalues of a 3x3
    real symmetric matrix. Square 3x3 matrix a; output eigenvalues in w;
    and output eigenvectors in v. Resulting eigenvalues/vectors are sorted
    in decreasing order; eigenvectors are normalized.*/
  static int Jacobi(const double A[3][3], double w[3], double v[3][3]);
 
  //int Jacobi(double w[3], double v[3][3]) {Jacobi(GetElements(),w,v);}
 
  static int JacobiN(double **a, int n, double *w, double **v);
 
  static void LUFactor(double A[3][3], int index[3]);
  void LUFactor(int index[3]) {LUFactor(GetElements(),index);}
 
  static void LUSolve(const double A[3][3], const int index[3], double x[3]);
  void LUSolve(const int index[3], double x[3]) \
    { LUSolve(GetElements(),index,x);}
 
  
  static double DegreesToRadians() {return 0.017453292519943295;};
  static double Pi() {return 3.1415926535897932384626;};
  static double RadiansToDegrees() {return 57.29577951308232;};
 
  double *operator[](const unsigned int i) {return &(GetElements()[i][0]);}
 
  const double *operator[](unsigned int i) const { return &(GetElements()[i][0]); }  
  
  static double Determinant2x2(double a, double b, double c, double d) {
    return (a * d - b * c);};
  static double Determinant2x2(const double c1[2], const double c2[2]) {
    return (c1[0]*c2[1] - c2[0]*c1[1]);};
 
  static double Norm(const double x[3]) {
    return sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);};
  
  static inline void Cross(const double x[3], const double y[3], double z[3]);
  static inline double Normalize(double x[3]);
  
protected:
 
  double m_Elements[3][3];  
};
 
//------------------------------------------------------------------------------
inline double albaMatrix3x3::Determinant(double A[3][3])
//------------------------------------------------------------------------------
{
  return A[0][0]*A[1][1]*A[2][2] + A[1][0]*A[2][1]*A[0][2] + 
         A[2][0]*A[0][1]*A[1][2] - A[0][0]*A[2][1]*A[1][2] - 
         A[1][0]*A[0][1]*A[2][2] - A[2][0]*A[1][1]*A[0][2];
}
//------------------------------------------------------------------------------
inline double albaMatrix3x3::Determinant(const double c1[3], 
                                      const double c2[3], 
                                      const double c3[3])
//------------------------------------------------------------------------------
{
  return c1[0]*c2[1]*c3[2] + c2[0]*c3[1]*c1[2] + c3[0]*c1[1]*c2[2] -
         c1[0]*c3[1]*c2[2] - c2[0]*c1[1]*c3[2] - c3[0]*c2[1]*c1[2];
}
//------------------------------------------------------------------------------
inline double albaMatrix3x3::Determinant(double a1, double a2, double a3, 
                                      double b1, double b2, double b3, 
                                      double c1, double c2, double c3)
//------------------------------------------------------------------------------
{
    return ( a1 * albaMatrix3x3::Determinant2x2( b2, b3, c2, c3 )
           - b1 * albaMatrix3x3::Determinant2x2( a2, a3, c2, c3 )
           + c1 * albaMatrix3x3::Determinant2x2( a2, a3, b2, b3 ) );
}
//------------------------------------------------------------------------------
// Cross product of two 3-vectors. Result vector in z[3].
inline void albaMatrix3x3::Cross(const double x[3], const double y[3], double z[3])
//------------------------------------------------------------------------------
{
  double Zx = x[1]*y[2] - x[2]*y[1]; 
  double Zy = x[2]*y[0] - x[0]*y[2];
  double Zz = x[0]*y[1] - x[1]*y[0];
  z[0] = Zx; z[1] = Zy; z[2] = Zz; 
}
//------------------------------------------------------------------------------
inline double albaMatrix3x3::Normalize(double x[3])
//------------------------------------------------------------------------------
{
  double den; 
  if ( (den = albaMatrix3x3::Norm(x)) != 0.0 )
  {
    for (int i=0; i < 3; i++)
    {
      x[i] /= den;
    }
  }
  return den;
}
 
#endif