/vol/vipdata/irtk/contrib++/include/irtkMatrix3x3.h

00001 /*=========================================================================
00002 
00003   Library   : Image Registration Toolkit (IRTK)
00004   Module    : $Id: irtkMatrix3x3.h 2 2008-12-23 12:40:14Z dr $
00005   Copyright : Imperial College, Department of Computing
00006               Visual Information Processing (VIP), 2008 onwards
00007   Date      : $Date: 2008-12-23 12:40:14 +0000 (Tue, 23 Dec 2008) $
00008   Version   : $Revision: 2 $
00009   Changes   : $Author: dr $
00010 
00011 =========================================================================*/
00012 
00013 #ifndef IRTKMATRIX3X3_H
00014 
00015 #define IRTKMATRIX3X3_H
00016 
00017 #include <irtkVector3.h>
00018 
00019 // NOTE.  The (x,y,z) coordinate system is assumed to be right-handed.
00020 // Coordinate axis rotation matrices are of the form
00021 //   RX =    1       0       0
00022 //           0     cos(t) -sin(t)
00023 //           0     sin(t)  cos(t)
00024 // where t > 0 indicates a counterclockwise rotation in the yz-plane
00025 //   RY =  cos(t)    0     sin(t)
00026 //           0       1       0
00027 //        -sin(t)    0     cos(t)
00028 // where t > 0 indicates a counterclockwise rotation in the zx-plane
00029 //   RZ =  cos(t) -sin(t)    0
00030 //         sin(t)  cos(t)    0
00031 //           0       0       1
00032 // where t > 0 indicates a counterclockwise rotation in the xy-plane.
00033 
00034 
00035 class irtkMatrix3x3
00036 {
00037 public:
00038   // construction
00039   irtkMatrix3x3 ();
00040   irtkMatrix3x3 (const double aafEntry[3][3]);
00041   irtkMatrix3x3 (const irtkMatrix3x3& rkMatrix);
00042   irtkMatrix3x3 (double fEntry00, double fEntry01, double fEntry02,
00043                  double fEntry10, double fEntry11, double fEntry12,
00044                  double fEntry20, double fEntry21, double fEntry22);
00045 
00046   // member access, allows use of construct mat[r][c]
00047   double* operator[] (int iRow) const;
00048   operator double* ();
00049   irtkVector3 GetColumn (int iCol) const;
00050 
00051   // assignment and comparison
00052   irtkMatrix3x3& operator= (const irtkMatrix3x3& rkMatrix);
00053   bool operator== (const irtkMatrix3x3& rkMatrix) const;
00054   bool operator!= (const irtkMatrix3x3& rkMatrix) const;
00055 
00056   // arithmetic operations
00057   irtkMatrix3x3 operator+ (const irtkMatrix3x3& rkMatrix) const;
00058   irtkMatrix3x3 operator- (const irtkMatrix3x3& rkMatrix) const;
00059   irtkMatrix3x3 operator* (const irtkMatrix3x3& rkMatrix) const;
00060   irtkMatrix3x3 operator- () const;
00061 
00062   // matrix * vector [3x3 * 3x1 = 3x1]
00063   irtkVector3 operator* (const irtkVector3& rkVector) const;
00064 
00065   // vector * matrix [1x3 * 3x3 = 1x3]
00066   friend irtkVector3 operator* (const irtkVector3& rkVector,
00067                                 const irtkMatrix3x3& rkMatrix);
00068 
00069   // matrix * scalar
00070   irtkMatrix3x3 operator* (double fScalar) const;
00071 
00072   // scalar * matrix
00073   friend irtkMatrix3x3 operator* (double fScalar, const irtkMatrix3x3& rkMatrix);
00074 
00075   // utilities
00076   irtkMatrix3x3 Transpose () const;
00077   bool Inverse (irtkMatrix3x3& rkInverse, double fTolerance = 1e-06) const;
00078   irtkMatrix3x3 Inverse (double fTolerance = 1e-06) const;
00079   double Determinant () const;
00080 
00081   // singular value decomposition
00082   void SingularValueDecomposition (irtkMatrix3x3& rkL, irtkVector3& rkS,
00083                                    irtkMatrix3x3& rkR) const;
00084   void SingularValueComposition (const irtkMatrix3x3& rkL,
00085                                  const irtkVector3& rkS, const irtkMatrix3x3& rkR);
00086 
00087   // Gram-Schmidt orthonormalization (applied to columns of rotation matrix)
00088   void Orthonormalize ();
00089 
00090   // orthogonal Q, diagonal D, upper triangular U stored as (u01,u02,u12)
00091   void QDUDecomposition (irtkMatrix3x3& rkQ, irtkVector3& rkD,
00092                          irtkVector3& rkU) const;
00093 
00094   double SpectralNorm () const;
00095 
00096   // matrix must be orthonormal
00097   void ToAxisAngle (irtkVector3& rkAxis, double& rfRadians) const;
00098   void FromAxisAngle (const irtkVector3& rkAxis, double fRadians);
00099 
00100   // The matrix must be orthonormal.  The decomposition is yaw*pitch*roll
00101   // where yaw is rotation about the Up vector, pitch is rotation about the
00102   // Right axis, and roll is rotation about the Direction axis.
00103   bool ToEulerAnglesXYZ (float& rfYAngle, float& rfPAngle,
00104                          float& rfRAngle) const;
00105   bool ToEulerAnglesXZY (float& rfYAngle, float& rfPAngle,
00106                          float& rfRAngle) const;
00107   bool ToEulerAnglesYXZ (float& rfYAngle, float& rfPAngle,
00108                          float& rfRAngle) const;
00109   bool ToEulerAnglesYZX (float& rfYAngle, float& rfPAngle,
00110                          float& rfRAngle) const;
00111   bool ToEulerAnglesZXY (float& rfYAngle, float& rfPAngle,
00112                          float& rfRAngle) const;
00113   bool ToEulerAnglesZYX (float& rfYAngle, float& rfPAngle,
00114                          float& rfRAngle) const;
00115   void FromEulerAnglesXYZ (float fYAngle, float fPAngle, float fRAngle);
00116   void FromEulerAnglesXZY (float fYAngle, float fPAngle, float fRAngle);
00117   void FromEulerAnglesYXZ (float fYAngle, float fPAngle, float fRAngle);
00118   void FromEulerAnglesYZX (float fYAngle, float fPAngle, float fRAngle);
00119   void FromEulerAnglesZXY (float fYAngle, float fPAngle, float fRAngle);
00120   void FromEulerAnglesZYX (float fYAngle, float fPAngle, float fRAngle);
00121 
00122   // eigensolver, matrix must be symmetric
00123   void EigenSolveSymmetric (double afEigenvalue[3],
00124                             irtkVector3 akEigenvector[3]) const;
00125 
00126   static void TensorProduct (const irtkVector3& rkU, const irtkVector3& rkV,
00127                              irtkMatrix3x3& rkProduct);
00128 
00129   static const double EPSILON;
00130   static const irtkMatrix3x3 ZERO;
00131   static const irtkMatrix3x3 IDENTITY;
00132 
00133 protected:
00134   // support for eigensolver
00135   void Tridiagonal (double afDiag[3], double afSubDiag[3]);
00136   bool QLAlgorithm (double afDiag[3], double afSubDiag[3]);
00137 
00138   // support for singular value decomposition
00139   static const double ms_fSvdEpsilon;
00140   static const int ms_iSvdMaxIterations;
00141   static void Bidiagonalize (irtkMatrix3x3& kA, irtkMatrix3x3& kL,
00142                              irtkMatrix3x3& kR);
00143   static void GolubKahanStep (irtkMatrix3x3& kA, irtkMatrix3x3& kL,
00144                               irtkMatrix3x3& kR);
00145 
00146   // support for spectral norm
00147   static double MaxCubicRoot (double afCoeff[3]);
00148 
00149   double m_aafEntry[3][3];
00150 };
00151 
00152 #endif