Wm4Matrix3.h

Go to the documentation of this file.
00001 // Wild Magic Source Code
00002 // David Eberly
00003 // http://www.geometrictools.com
00004 // Copyright (c) 1998-2008
00005 //
00006 // This library is free software; you can redistribute it and/or modify it
00007 // under the terms of the GNU Lesser General Public License as published by
00008 // the Free Software Foundation; either version 2.1 of the License, or (at
00009 // your option) any later version.  The license is available for reading at
00010 // either of the locations:
00011 //     http://www.gnu.org/copyleft/lgpl.html
00012 //     http://www.geometrictools.com/License/WildMagicLicense.pdf
00013 //
00014 // Version: 4.0.5 (2007/08/31)
00015 
00016 #ifndef WM4MATRIX3_H
00017 #define WM4MATRIX3_H
00018 
00019 // Matrix operations are applied on the left.  For example, given a matrix M
00020 // and a vector V, matrix-times-vector is M*V.  That is, V is treated as a
00021 // column vector.  Some graphics APIs use V*M where V is treated as a row
00022 // vector.  In this context the "M" matrix is really a transpose of the M as
00023 // represented in Wild Magic.  Similarly, to apply two matrix operations M0
00024 // and M1, in that order, you compute M1*M0 so that the transform of a vector
00025 // is (M1*M0)*V = M1*(M0*V).  Some graphics APIs use M0*M1, but again these
00026 // matrices are the transpose of those as represented in Wild Magic.  You
00027 // must therefore be careful about how you interface the transformation code
00028 // with graphics APIS.
00029 //
00030 // For memory organization it might seem natural to chose Real[N][N] for the
00031 // matrix storage, but this can be a problem on a platform/console that
00032 // chooses to store the data in column-major rather than row-major format.
00033 // To avoid potential portability problems, the matrix is stored as Real[N*N]
00034 // and organized in row-major order.  That is, the entry of the matrix in row
00035 // r (0 <= r < N) and column c (0 <= c < N) is stored at index i = c+N*r
00036 // (0 <= i < N*N).
00037 
00038 // The (x,y,z) coordinate system is assumed to be right-handed.  Coordinate
00039 // axis rotation matrices are of the form
00040 //   RX =    1       0       0
00041 //           0     cos(t) -sin(t)
00042 //           0     sin(t)  cos(t)
00043 // where t > 0 indicates a counterclockwise rotation in the yz-plane
00044 //   RY =  cos(t)    0     sin(t)
00045 //           0       1       0
00046 //        -sin(t)    0     cos(t)
00047 // where t > 0 indicates a counterclockwise rotation in the zx-plane
00048 //   RZ =  cos(t) -sin(t)    0
00049 //         sin(t)  cos(t)    0
00050 //           0       0       1
00051 // where t > 0 indicates a counterclockwise rotation in the xy-plane.
00052 
00053 #include "Wm4FoundationLIB.h"
00054 #include "Wm4Vector3.h"
00055 
00056 namespace Wm4
00057 {
00058 
00059 template <class Real>
00060 class Matrix3
00061 {
00062 public:
00063     // If bZero is true, create the zero matrix.  Otherwise, create the
00064     // identity matrix.
00065     Matrix3 (bool bZero = true);
00066 
00067     // copy constructor
00068     Matrix3 (const Matrix3& rkM);
00069 
00070     // input Mrc is in row r, column c.
00071     Matrix3 (Real fM00, Real fM01, Real fM02,
00072              Real fM10, Real fM11, Real fM12,
00073              Real fM20, Real fM21, Real fM22);
00074 
00075     // Create a matrix from an array of numbers.  The input array is
00076     // interpreted based on the Boolean input as
00077     //   true:  entry[0..8]={m00,m01,m02,m10,m11,m12,m20,m21,m22} [row major]
00078     //   false: entry[0..8]={m00,m10,m20,m01,m11,m21,m02,m12,m22} [col major]
00079     Matrix3 (const Real afEntry[9], bool bRowMajor);
00080 
00081     // Create matrices based on vector input.  The Boolean is interpreted as
00082     //   true: vectors are columns of the matrix
00083     //   false: vectors are rows of the matrix
00084     Matrix3 (const Vector3<Real>& rkU, const Vector3<Real>& rkV,
00085         const Vector3<Real>& rkW, bool bColumns);
00086     Matrix3 (const Vector3<Real>* akV, bool bColumns);
00087 
00088     // create a diagonal matrix
00089     Matrix3 (Real fM00, Real fM11, Real fM22);
00090 
00091     // Create rotation matrices (positive angle - counterclockwise).  The
00092     // angle must be in radians, not degrees.
00093     Matrix3 (const Vector3<Real>& rkAxis, Real fAngle);
00094 
00095     // create a tensor product U*V^T
00096     Matrix3 (const Vector3<Real>& rkU, const Vector3<Real>& rkV);
00097 
00098     // create various matrices
00099     Matrix3& MakeZero ();
00100     Matrix3& MakeIdentity ();
00101     Matrix3& MakeDiagonal (Real fM00, Real fM11, Real fM22);
00102     Matrix3& FromAxisAngle (const Vector3<Real>& rkAxis, Real fAngle);
00103     Matrix3& MakeTensorProduct (const Vector3<Real>& rkU,
00104         const Vector3<Real>& rkV);
00105 
00106     // member access
00107     inline operator const Real* () const;
00108     inline operator Real* ();
00109     inline const Real* operator[] (int iRow) const;
00110     inline Real* operator[] (int iRow);
00111     inline Real operator() (int iRow, int iCol) const;
00112     inline Real& operator() (int iRow, int iCol);
00113     void SetRow (int iRow, const Vector3<Real>& rkV);
00114     Vector3<Real> GetRow (int iRow) const;
00115     void SetColumn (int iCol, const Vector3<Real>& rkV);
00116     Vector3<Real> GetColumn (int iCol) const;
00117     void GetColumnMajor (Real* afCMajor) const;
00118 
00119     // assignment
00120     inline Matrix3& operator= (const Matrix3& rkM);
00121 
00122     // comparison
00123     bool operator== (const Matrix3& rkM) const;
00124     bool operator!= (const Matrix3& rkM) const;
00125     bool operator<  (const Matrix3& rkM) const;
00126     bool operator<= (const Matrix3& rkM) const;
00127     bool operator>  (const Matrix3& rkM) const;
00128     bool operator>= (const Matrix3& rkM) const;
00129 
00130     // arithmetic operations
00131     inline Matrix3 operator+ (const Matrix3& rkM) const;
00132     inline Matrix3 operator- (const Matrix3& rkM) const;
00133     inline Matrix3 operator* (const Matrix3& rkM) const;
00134     inline Matrix3 operator* (Real fScalar) const;
00135     inline Matrix3 operator/ (Real fScalar) const;
00136     inline Matrix3 operator- () const;
00137 
00138     // arithmetic updates
00139     inline Matrix3& operator+= (const Matrix3& rkM);
00140     inline Matrix3& operator-= (const Matrix3& rkM);
00141     inline Matrix3& operator*= (Real fScalar);
00142     inline Matrix3& operator/= (Real fScalar);
00143 
00144     // matrix times vector
00145     inline Vector3<Real> operator* (const Vector3<Real>& rkV) const;  // M * v
00146 
00147     // other operations
00148     Matrix3 Transpose () const;  // M^T
00149     Matrix3 TransposeTimes (const Matrix3& rkM) const;  // this^T * M
00150     Matrix3 TimesTranspose (const Matrix3& rkM) const;  // this * M^T
00151     Matrix3 Inverse () const;
00152     Matrix3 Adjoint () const;
00153     Real Determinant () const;
00154     Real QForm (const Vector3<Real>& rkU,
00155         const Vector3<Real>& rkV) const;  // u^T*M*v
00156     Matrix3 TimesDiagonal (const Vector3<Real>& rkDiag) const;  // M*D
00157     Matrix3 DiagonalTimes (const Vector3<Real>& rkDiag) const;  // D*M
00158 
00159     // The matrix must be a rotation for these functions to be valid.  The
00160     // last function uses Gram-Schmidt orthonormalization applied to the
00161     // columns of the rotation matrix.  The angle must be in radians, not
00162     // degrees.
00163     void ToAxisAngle (Vector3<Real>& rkAxis, Real& rfAngle) const;
00164     void Orthonormalize ();
00165 
00166     // The matrix must be symmetric.  Factor M = R * D * R^T where
00167     // R = [u0|u1|u2] is a rotation matrix with columns u0, u1, and u2 and
00168     // D = diag(d0,d1,d2) is a diagonal matrix whose diagonal entries are d0,
00169     // d1, and d2.  The eigenvector u[i] corresponds to eigenvector d[i].
00170     // The eigenvalues are ordered as d0 <= d1 <= d2.
00171     void EigenDecomposition (Matrix3& rkRot, Matrix3& rkDiag) const;
00172 
00173     // Create rotation matrices from Euler angles.
00174     Matrix3& FromEulerAnglesXYZ (Real fXAngle, Real fYAngle, Real fZAngle);
00175     Matrix3& FromEulerAnglesXZY (Real fXAngle, Real fZAngle, Real fYAngle);
00176     Matrix3& FromEulerAnglesYXZ (Real fYAngle, Real fXAngle, Real fZAngle);
00177     Matrix3& FromEulerAnglesYZX (Real fYAngle, Real fZAngle, Real fXAngle);
00178     Matrix3& FromEulerAnglesZXY (Real fZAngle, Real fXAngle, Real fYAngle);
00179     Matrix3& FromEulerAnglesZYX (Real fZAngle, Real fYAngle, Real fXAngle);
00180 
00181     // TODO.  From EulerAnglesUVU for all combinations of U and V.
00182 
00183     // Extract Euler angles from rotation matrices.
00184     enum EulerResult
00185     {
00186         // The solution is unique.
00187         EA_UNIQUE,
00188 
00189         // The solution is not unique.  A sum of angles is constant.
00190         EA_NOT_UNIQUE_SUM,
00191 
00192         // The solution is not unique.  A difference of angles is constant.
00193         EA_NOT_UNIQUE_DIF
00194     };
00195 
00196     // The return values are in the specified ranges:
00197     //   xAngle in [-pi,pi], yAngle in [-pi/2,pi/2], zAngle in [-pi,pi]
00198     // When the solution is not unique, zAngle = 0 is returned.  Generally,
00199     // the set of solutions is
00200     //   EA_NOT_UNIQUE_SUM:  zAngle + xAngle = c
00201     //   EA_NOT_UNIQUE_DIF:  zAngle - xAngle = c
00202     // for some angle c.
00203     EulerResult ToEulerAnglesXYZ (Real& rfXAngle, Real& rfYAngle,
00204         Real& rfZAngle) const;
00205 
00206     // The return values are in the specified ranges:
00207     //   xAngle in [-pi,pi], zAngle in [-pi/2,pi/2], yAngle in [-pi,pi]
00208     // When the solution is not unique, yAngle = 0 is returned.  Generally,
00209     // the set of solutions is
00210     //   EA_NOT_UNIQUE_SUM:  yAngle + xAngle = c
00211     //   EA_NOT_UNIQUE_DIF:  yAngle - xAngle = c
00212     // for some angle c.
00213     EulerResult ToEulerAnglesXZY (Real& rfXAngle, Real& rfZAngle,
00214         Real& rfYAngle) const;
00215 
00216     // The return values are in the specified ranges:
00217     //   yAngle in [-pi,pi], xAngle in [-pi/2,pi/2], zAngle in [-pi,pi]
00218     // When the solution is not unique, zAngle = 0 is returned.  Generally,
00219     // the set of solutions is
00220     //   EA_NOT_UNIQUE_SUM:  zAngle + yAngle = c
00221     //   EA_NOT_UNIQUE_DIF:  zAngle - yAngle = c
00222     // for some angle c.
00223     EulerResult ToEulerAnglesYXZ (Real& rfYAngle, Real& rfXAngle,
00224         Real& rfZAngle) const;
00225 
00226     // The return values are in the specified ranges:
00227     //   yAngle in [-pi,pi], zAngle in [-pi/2,pi/2], xAngle in [-pi,pi]
00228     // When the solution is not unique, xAngle = 0 is returned.  Generally,
00229     // the set of solutions is
00230     //   EA_NOT_UNIQUE_SUM:  xAngle + yAngle = c
00231     //   EA_NOT_UNIQUE_DIF:  xAngle - yAngle = c
00232     // for some angle c.
00233     EulerResult ToEulerAnglesYZX (Real& rfYAngle, Real& rfZAngle,
00234         Real& rfXAngle) const;
00235 
00236     // The return values are in the specified ranges:
00237     //   zAngle in [-pi,pi], xAngle in [-pi/2,pi/2], yAngle in [-pi,pi]
00238     // When the solution is not unique, yAngle = 0 is returned.  Generally,
00239     // the set of solutions is
00240     //   EA_NOT_UNIQUE_SUM:  yAngle + zAngle = c
00241     //   EA_NOT_UNIQUE_DIF:  yAngle - zAngle = c
00242     // for some angle c.
00243     EulerResult ToEulerAnglesZXY (Real& rfZAngle, Real& rfXAngle,
00244         Real& rfYAngle) const;
00245 
00246     // The return values are in the specified ranges:
00247     //   zAngle in [-pi,pi], yAngle in [-pi/2,pi/2], xAngle in [-pi,pi]
00248     // When the solution is not unique, xAngle = 0 is/ returned.  Generally,
00249     // the set of solutions is
00250     //   EA_NOT_UNIQUE_SUM:  xAngle + zAngle = c
00251     //   EA_NOT_UNIQUE_DIF:  xAngle - zAngle = c
00252     // for some angle c.
00253     EulerResult ToEulerAnglesZYX (Real& rfZAngle, Real& rfYAngle,
00254         Real& rfXAngle) const;
00255 
00256     // The return values are in the specified ranges:
00257     //   x0Angle in [-pi,pi], yAngle in [0,pi], x1Angle in [-pi,pi]
00258     // When the solution is not unique, x1Angle = 0 is returned.  Generally,
00259     // the set of solutions is
00260     //   EA_NOT_UNIQUE_SUM:  x1Angle + x0Angle = c
00261     //   EA_NOT_UNIQUE_DIF:  x1Angle - x0Angle = c
00262     // for some angle c.
00263     EulerResult ToEulerAnglesXYX (Real& rfX0Angle, Real& rfYAngle,
00264         Real& rfX1Angle) const;
00265 
00266     // The return values are in the specified ranges:
00267     //   x0Angle in [-pi,pi], zAngle in [0,pi], x1Angle in [-pi,pi]
00268     // When the solution is not unique, x1Angle = 0 is returned.  Generally,
00269     // the set of solutions is
00270     //   EA_NOT_UNIQUE_SUM:  x1Angle + x0Angle = c
00271     //   EA_NOT_UNIQUE_DIF:  x1Angle - x0Angle = c
00272     // for some angle c.
00273     EulerResult ToEulerAnglesXZX (Real& rfX0Angle, Real& rfZAngle,
00274         Real& rfX1Angle) const;
00275 
00276     // The return values are in the specified ranges:
00277     //   y0Angle in [-pi,pi], xAngle in [0,pi], y1Angle in [-pi,pi]
00278     // When the solution is not unique, y1Angle = 0 is returned.  Generally,
00279     // the set of solutions is
00280     //   EA_NOT_UNIQUE_SUM:  y1Angle + y0Angle = c
00281     //   EA_NOT_UNIQUE_DIF:  y1Angle - y0Angle = c
00282     // for some angle c.
00283     EulerResult ToEulerAnglesYXY (Real& rfY0Angle, Real& rfXAngle,
00284         Real& rfY1Angle) const;
00285 
00286     // The return values are in the specified ranges:
00287     //   y0Angle in [-pi,pi], zAngle in [0,pi], y1Angle in [-pi,pi]
00288     // When the solution is not unique, y1Angle = 0 is returned.  Generally,
00289     // the set of solutions is
00290     //   EA_NOT_UNIQUE_SUM:  y1Angle + y0Angle = c
00291     //   EA_NOT_UNIQUE_DIF:  y1Angle - y0Angle = c
00292     // for some angle c.
00293     EulerResult ToEulerAnglesYZY (Real& rfY0Angle, Real& rfZAngle,
00294         Real& rfY1Angle) const;
00295 
00296     // The return values are in the specified ranges:
00297     //   z0Angle in [-pi,pi], xAngle in [0,pi], z1Angle in [-pi,pi]
00298     // When the solution is not unique, z1Angle = 0 is returned.  Generally,
00299     // the set of solutions is
00300     //   EA_NOT_UNIQUE_SUM:  z1Angle + z0Angle = c
00301     //   EA_NOT_UNIQUE_DIF:  z1Angle - z0Angle = c
00302     // for some angle c.
00303     EulerResult ToEulerAnglesZXZ (Real& rfZ0Angle, Real& rfXAngle,
00304         Real& rfZ1Angle) const;
00305 
00306     // The return values are in the specified ranges:
00307     //   z0Angle in [-pi,pi], yAngle in [0,pi], z1Angle in [-pi,pi]
00308     // When the solution is not unique, z1Angle = 0 is returned.  Generally,
00309     // the set of solutions is
00310     //   EA_NOT_UNIQUE_SUM:  z1Angle + z0Angle = c
00311     //   EA_NOT_UNIQUE_DIF:  z1Angle - z0Angle = c
00312     // for some angle c.
00313     EulerResult ToEulerAnglesZYZ (Real& rfZ0Angle, Real& rfYAngle,
00314         Real& rfZ1Angle) const;
00315 
00316     // SLERP (spherical linear interpolation) without quaternions.  Computes
00317     // R(t) = R0*(Transpose(R0)*R1)^t.  If Q is a rotation matrix with
00318     // unit-length axis U and angle A, then Q^t is a rotation matrix with
00319     // unit-length axis U and rotation angle t*A.
00320     Matrix3& Slerp (Real fT, const Matrix3& rkR0, const Matrix3& rkR1);
00321 
00322     // Singular value decomposition, M = L*D*Transpose(R), where L and R are
00323     // orthogonal and D is a diagonal matrix whose diagonal entries are
00324     // nonnegative.
00325     void SingularValueDecomposition (Matrix3& rkL, Matrix3& rkD,
00326         Matrix3& rkRTranspose) const;
00327 
00328     // For debugging purposes.  Take the output of SingularValueDecomposition
00329     // and multiply to see if you get M.
00330     void SingularValueComposition (const Matrix3& rkL, const Matrix3& rkD,
00331         const Matrix3& rkRTranspose);
00332 
00333     // Polar decomposition, M = Q*S, where Q is orthogonal and S is symmetric.
00334     // This uses the singular value decomposition:
00335     //   M = L*D*Transpose(R) = (L*Transpose(R))*(R*D*Transpose(R)) = Q*S
00336     // where Q = L*Transpose(R) and S = R*D*Transpose(R).
00337     void PolarDecomposition (Matrix3& rkQ, Matrix3& rkS);
00338 
00339     // Factor M = Q*D*U with orthogonal Q, diagonal D, upper triangular U.
00340     void QDUDecomposition (Matrix3& rkQ, Matrix3& rkD, Matrix3& rkU) const;
00341 
00342     // special matrices
00343     WM4_FOUNDATION_ITEM static const Matrix3 ZERO;
00344     WM4_FOUNDATION_ITEM static const Matrix3 IDENTITY;
00345 
00346 private:
00347     // Support for eigendecomposition.  The Tridiagonalize function applies
00348     // a Householder transformation to the matrix.  If that transformation
00349     // is the identity (the matrix is already tridiagonal), then the return
00350     // value is 'false'.  Otherwise, the transformation is a reflection and
00351     // the return value is 'true'.  The QLAlgorithm returns 'true' iff the
00352     // QL iteration scheme converged.
00353     bool Tridiagonalize (Real afDiag[3], Real afSubd[2]);
00354     bool QLAlgorithm (Real afDiag[3], Real afSubd[2]);
00355 
00356     // support for singular value decomposition
00357     static void Bidiagonalize (Matrix3& rkA, Matrix3& rkL, Matrix3& rkR);
00358     static void GolubKahanStep (Matrix3& rkA, Matrix3& rkL, Matrix3& rkR);
00359 
00360     // support for comparisons
00361     int CompareArrays (const Matrix3& rkM) const;
00362 
00363     Real m_afEntry[9];
00364 };
00365 
00366 // c * M
00367 template <class Real>
00368 inline Matrix3<Real> operator* (Real fScalar, const Matrix3<Real>& rkM);
00369 
00370 // v^T * M
00371 template <class Real>
00372 inline Vector3<Real> operator* (const Vector3<Real>& rkV,
00373     const Matrix3<Real>& rkM);
00374 
00375 #include "Wm4Matrix3.inl"
00376 
00377 typedef Matrix3<float> Matrix3f;
00378 typedef Matrix3<double> Matrix3d;
00379 
00380 }
00381 
00382 #endif

Generated on Fri Feb 13 13:58:10 2009 for meshmorph by  doxygen 1.5.1