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