PBRT
|
00001 00002 // 00003 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas 00004 // Digital Ltd. LLC 00005 // 00006 // All rights reserved. 00007 // 00008 // Redistribution and use in source and binary forms, with or without 00009 // modification, are permitted provided that the following conditions are 00010 // met: 00011 // * Redistributions of source code must retain the above copyright 00012 // notice, this list of conditions and the following disclaimer. 00013 // * Redistributions in binary form must reproduce the above 00014 // copyright notice, this list of conditions and the following disclaimer 00015 // in the documentation and/or other materials provided with the 00016 // distribution. 00017 // * Neither the name of Industrial Light & Magic nor the names of 00018 // its contributors may be used to endorse or promote products derived 00019 // from this software without specific prior written permission. 00020 // 00021 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 00022 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 00023 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 00024 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 00025 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 00027 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 00028 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 00029 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00030 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 00031 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00032 // 00034 00035 00036 00037 #ifndef INCLUDED_IMATHMATRIX_H 00038 #define INCLUDED_IMATHMATRIX_H 00039 00040 //---------------------------------------------------------------- 00041 // 00042 // 2D (3x3) and 3D (4x4) transformation matrix templates. 00043 // 00044 //---------------------------------------------------------------- 00045 00046 #include "ImathPlatform.h" 00047 #include "ImathFun.h" 00048 #include "ImathExc.h" 00049 #include "ImathVec.h" 00050 #include "ImathShear.h" 00051 00052 #include <iostream> 00053 #include <iomanip> 00054 00055 #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER 00056 // suppress exception specification warnings 00057 #pragma warning(disable:4290) 00058 #endif 00059 00060 00061 namespace Imath { 00062 00063 enum Uninitialized {UNINITIALIZED}; 00064 00065 00066 template <class T> class Matrix33 00067 { 00068 public: 00069 00070 //------------------- 00071 // Access to elements 00072 //------------------- 00073 00074 T x[3][3]; 00075 00076 T * operator [] (int i); 00077 const T * operator [] (int i) const; 00078 00079 00080 //------------- 00081 // Constructors 00082 //------------- 00083 00084 Matrix33 (Uninitialized) {} 00085 00086 Matrix33 (); 00087 // 1 0 0 00088 // 0 1 0 00089 // 0 0 1 00090 00091 Matrix33 (T a); 00092 // a a a 00093 // a a a 00094 // a a a 00095 00096 Matrix33 (const T a[3][3]); 00097 // a[0][0] a[0][1] a[0][2] 00098 // a[1][0] a[1][1] a[1][2] 00099 // a[2][0] a[2][1] a[2][2] 00100 00101 Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i); 00102 00103 // a b c 00104 // d e f 00105 // g h i 00106 00107 00108 //-------------------------------- 00109 // Copy constructor and assignment 00110 //-------------------------------- 00111 00112 Matrix33 (const Matrix33 &v); 00113 template <class S> explicit Matrix33 (const Matrix33<S> &v); 00114 00115 const Matrix33 & operator = (const Matrix33 &v); 00116 const Matrix33 & operator = (T a); 00117 00118 00119 //---------------------- 00120 // Compatibility with Sb 00121 //---------------------- 00122 00123 T * getValue (); 00124 const T * getValue () const; 00125 00126 template <class S> 00127 void getValue (Matrix33<S> &v) const; 00128 template <class S> 00129 Matrix33 & setValue (const Matrix33<S> &v); 00130 00131 template <class S> 00132 Matrix33 & setTheMatrix (const Matrix33<S> &v); 00133 00134 00135 //--------- 00136 // Identity 00137 //--------- 00138 00139 void makeIdentity(); 00140 00141 00142 //--------- 00143 // Equality 00144 //--------- 00145 00146 bool operator == (const Matrix33 &v) const; 00147 bool operator != (const Matrix33 &v) const; 00148 00149 //----------------------------------------------------------------------- 00150 // Compare two matrices and test if they are "approximately equal": 00151 // 00152 // equalWithAbsError (m, e) 00153 // 00154 // Returns true if the coefficients of this and m are the same with 00155 // an absolute error of no more than e, i.e., for all i, j 00156 // 00157 // abs (this[i][j] - m[i][j]) <= e 00158 // 00159 // equalWithRelError (m, e) 00160 // 00161 // Returns true if the coefficients of this and m are the same with 00162 // a relative error of no more than e, i.e., for all i, j 00163 // 00164 // abs (this[i] - v[i][j]) <= e * abs (this[i][j]) 00165 //----------------------------------------------------------------------- 00166 00167 bool equalWithAbsError (const Matrix33<T> &v, T e) const; 00168 bool equalWithRelError (const Matrix33<T> &v, T e) const; 00169 00170 00171 //------------------------ 00172 // Component-wise addition 00173 //------------------------ 00174 00175 const Matrix33 & operator += (const Matrix33 &v); 00176 const Matrix33 & operator += (T a); 00177 Matrix33 operator + (const Matrix33 &v) const; 00178 00179 00180 //--------------------------- 00181 // Component-wise subtraction 00182 //--------------------------- 00183 00184 const Matrix33 & operator -= (const Matrix33 &v); 00185 const Matrix33 & operator -= (T a); 00186 Matrix33 operator - (const Matrix33 &v) const; 00187 00188 00189 //------------------------------------ 00190 // Component-wise multiplication by -1 00191 //------------------------------------ 00192 00193 Matrix33 operator - () const; 00194 const Matrix33 & negate (); 00195 00196 00197 //------------------------------ 00198 // Component-wise multiplication 00199 //------------------------------ 00200 00201 const Matrix33 & operator *= (T a); 00202 Matrix33 operator * (T a) const; 00203 00204 00205 //----------------------------------- 00206 // Matrix-times-matrix multiplication 00207 //----------------------------------- 00208 00209 const Matrix33 & operator *= (const Matrix33 &v); 00210 Matrix33 operator * (const Matrix33 &v) const; 00211 00212 00213 //----------------------------------------------------------------- 00214 // Vector-times-matrix multiplication; see also the "operator *" 00215 // functions defined below. 00216 // 00217 // m.multVecMatrix(src,dst) implements a homogeneous transformation 00218 // by computing Vec3 (src.x, src.y, 1) * m and dividing by the 00219 // result's third element. 00220 // 00221 // m.multDirMatrix(src,dst) multiplies src by the upper left 2x2 00222 // submatrix, ignoring the rest of matrix m. 00223 //----------------------------------------------------------------- 00224 00225 template <class S> 00226 void multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const; 00227 00228 template <class S> 00229 void multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const; 00230 00231 00232 //------------------------ 00233 // Component-wise division 00234 //------------------------ 00235 00236 const Matrix33 & operator /= (T a); 00237 Matrix33 operator / (T a) const; 00238 00239 00240 //------------------ 00241 // Transposed matrix 00242 //------------------ 00243 00244 const Matrix33 & transpose (); 00245 Matrix33 transposed () const; 00246 00247 00248 //------------------------------------------------------------ 00249 // Inverse matrix: If singExc is false, inverting a singular 00250 // matrix produces an identity matrix. If singExc is true, 00251 // inverting a singular matrix throws a SingMatrixExc. 00252 // 00253 // inverse() and invert() invert matrices using determinants; 00254 // gjInverse() and gjInvert() use the Gauss-Jordan method. 00255 // 00256 // inverse() and invert() are significantly faster than 00257 // gjInverse() and gjInvert(), but the results may be slightly 00258 // less accurate. 00259 // 00260 //------------------------------------------------------------ 00261 00262 const Matrix33 & invert (bool singExc = false) 00263 throw (Iex::MathExc); 00264 00265 Matrix33<T> inverse (bool singExc = false) const 00266 throw (Iex::MathExc); 00267 00268 const Matrix33 & gjInvert (bool singExc = false) 00269 throw (Iex::MathExc); 00270 00271 Matrix33<T> gjInverse (bool singExc = false) const 00272 throw (Iex::MathExc); 00273 00274 00275 //----------------------------------------- 00276 // Set matrix to rotation by r (in radians) 00277 //----------------------------------------- 00278 00279 template <class S> 00280 const Matrix33 & setRotation (S r); 00281 00282 00283 //----------------------------- 00284 // Rotate the given matrix by r 00285 //----------------------------- 00286 00287 template <class S> 00288 const Matrix33 & rotate (S r); 00289 00290 00291 //-------------------------------------------- 00292 // Set matrix to scale by given uniform factor 00293 //-------------------------------------------- 00294 00295 const Matrix33 & setScale (T s); 00296 00297 00298 //------------------------------------ 00299 // Set matrix to scale by given vector 00300 //------------------------------------ 00301 00302 template <class S> 00303 const Matrix33 & setScale (const Vec2<S> &s); 00304 00305 00306 //---------------------- 00307 // Scale the matrix by s 00308 //---------------------- 00309 00310 template <class S> 00311 const Matrix33 & scale (const Vec2<S> &s); 00312 00313 00314 //------------------------------------------ 00315 // Set matrix to translation by given vector 00316 //------------------------------------------ 00317 00318 template <class S> 00319 const Matrix33 & setTranslation (const Vec2<S> &t); 00320 00321 00322 //----------------------------- 00323 // Return translation component 00324 //----------------------------- 00325 00326 Vec2<T> translation () const; 00327 00328 00329 //-------------------------- 00330 // Translate the matrix by t 00331 //-------------------------- 00332 00333 template <class S> 00334 const Matrix33 & translate (const Vec2<S> &t); 00335 00336 00337 //----------------------------------------------------------- 00338 // Set matrix to shear x for each y coord. by given factor xy 00339 //----------------------------------------------------------- 00340 00341 template <class S> 00342 const Matrix33 & setShear (const S &h); 00343 00344 00345 //------------------------------------------------------------- 00346 // Set matrix to shear x for each y coord. by given factor h[0] 00347 // and to shear y for each x coord. by given factor h[1] 00348 //------------------------------------------------------------- 00349 00350 template <class S> 00351 const Matrix33 & setShear (const Vec2<S> &h); 00352 00353 00354 //----------------------------------------------------------- 00355 // Shear the matrix in x for each y coord. by given factor xy 00356 //----------------------------------------------------------- 00357 00358 template <class S> 00359 const Matrix33 & shear (const S &xy); 00360 00361 00362 //----------------------------------------------------------- 00363 // Shear the matrix in x for each y coord. by given factor xy 00364 // and shear y for each x coord. by given factor yx 00365 //----------------------------------------------------------- 00366 00367 template <class S> 00368 const Matrix33 & shear (const Vec2<S> &h); 00369 00370 00371 //------------------------------------------------- 00372 // Limitations of type T (see also class limits<T>) 00373 //------------------------------------------------- 00374 00375 static T baseTypeMin() {return limits<T>::min();} 00376 static T baseTypeMax() {return limits<T>::max();} 00377 static T baseTypeSmallest() {return limits<T>::smallest();} 00378 static T baseTypeEpsilon() {return limits<T>::epsilon();} 00379 00380 private: 00381 00382 template <typename R, typename S> 00383 struct isSameType 00384 { 00385 enum {value = 0}; 00386 }; 00387 00388 template <typename R> 00389 struct isSameType<R, R> 00390 { 00391 enum {value = 1}; 00392 }; 00393 }; 00394 00395 00396 template <class T> class Matrix44 00397 { 00398 public: 00399 00400 //------------------- 00401 // Access to elements 00402 //------------------- 00403 00404 T x[4][4]; 00405 00406 T * operator [] (int i); 00407 const T * operator [] (int i) const; 00408 00409 00410 //------------- 00411 // Constructors 00412 //------------- 00413 00414 Matrix44 (Uninitialized) {} 00415 00416 Matrix44 (); 00417 // 1 0 0 0 00418 // 0 1 0 0 00419 // 0 0 1 0 00420 // 0 0 0 1 00421 00422 Matrix44 (T a); 00423 // a a a a 00424 // a a a a 00425 // a a a a 00426 // a a a a 00427 00428 Matrix44 (const T a[4][4]) ; 00429 // a[0][0] a[0][1] a[0][2] a[0][3] 00430 // a[1][0] a[1][1] a[1][2] a[1][3] 00431 // a[2][0] a[2][1] a[2][2] a[2][3] 00432 // a[3][0] a[3][1] a[3][2] a[3][3] 00433 00434 Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h, 00435 T i, T j, T k, T l, T m, T n, T o, T p); 00436 00437 // a b c d 00438 // e f g h 00439 // i j k l 00440 // m n o p 00441 00442 Matrix44 (Matrix33<T> r, Vec3<T> t); 00443 // r r r 0 00444 // r r r 0 00445 // r r r 0 00446 // t t t 1 00447 00448 00449 //-------------------------------- 00450 // Copy constructor and assignment 00451 //-------------------------------- 00452 00453 Matrix44 (const Matrix44 &v); 00454 template <class S> explicit Matrix44 (const Matrix44<S> &v); 00455 00456 const Matrix44 & operator = (const Matrix44 &v); 00457 const Matrix44 & operator = (T a); 00458 00459 00460 //---------------------- 00461 // Compatibility with Sb 00462 //---------------------- 00463 00464 T * getValue (); 00465 const T * getValue () const; 00466 00467 template <class S> 00468 void getValue (Matrix44<S> &v) const; 00469 template <class S> 00470 Matrix44 & setValue (const Matrix44<S> &v); 00471 00472 template <class S> 00473 Matrix44 & setTheMatrix (const Matrix44<S> &v); 00474 00475 //--------- 00476 // Identity 00477 //--------- 00478 00479 void makeIdentity(); 00480 00481 00482 //--------- 00483 // Equality 00484 //--------- 00485 00486 bool operator == (const Matrix44 &v) const; 00487 bool operator != (const Matrix44 &v) const; 00488 00489 //----------------------------------------------------------------------- 00490 // Compare two matrices and test if they are "approximately equal": 00491 // 00492 // equalWithAbsError (m, e) 00493 // 00494 // Returns true if the coefficients of this and m are the same with 00495 // an absolute error of no more than e, i.e., for all i, j 00496 // 00497 // abs (this[i][j] - m[i][j]) <= e 00498 // 00499 // equalWithRelError (m, e) 00500 // 00501 // Returns true if the coefficients of this and m are the same with 00502 // a relative error of no more than e, i.e., for all i, j 00503 // 00504 // abs (this[i] - v[i][j]) <= e * abs (this[i][j]) 00505 //----------------------------------------------------------------------- 00506 00507 bool equalWithAbsError (const Matrix44<T> &v, T e) const; 00508 bool equalWithRelError (const Matrix44<T> &v, T e) const; 00509 00510 00511 //------------------------ 00512 // Component-wise addition 00513 //------------------------ 00514 00515 const Matrix44 & operator += (const Matrix44 &v); 00516 const Matrix44 & operator += (T a); 00517 Matrix44 operator + (const Matrix44 &v) const; 00518 00519 00520 //--------------------------- 00521 // Component-wise subtraction 00522 //--------------------------- 00523 00524 const Matrix44 & operator -= (const Matrix44 &v); 00525 const Matrix44 & operator -= (T a); 00526 Matrix44 operator - (const Matrix44 &v) const; 00527 00528 00529 //------------------------------------ 00530 // Component-wise multiplication by -1 00531 //------------------------------------ 00532 00533 Matrix44 operator - () const; 00534 const Matrix44 & negate (); 00535 00536 00537 //------------------------------ 00538 // Component-wise multiplication 00539 //------------------------------ 00540 00541 const Matrix44 & operator *= (T a); 00542 Matrix44 operator * (T a) const; 00543 00544 00545 //----------------------------------- 00546 // Matrix-times-matrix multiplication 00547 //----------------------------------- 00548 00549 const Matrix44 & operator *= (const Matrix44 &v); 00550 Matrix44 operator * (const Matrix44 &v) const; 00551 00552 static void multiply (const Matrix44 &a, // assumes that 00553 const Matrix44 &b, // &a != &c and 00554 Matrix44 &c); // &b != &c. 00555 00556 00557 //----------------------------------------------------------------- 00558 // Vector-times-matrix multiplication; see also the "operator *" 00559 // functions defined below. 00560 // 00561 // m.multVecMatrix(src,dst) implements a homogeneous transformation 00562 // by computing Vec4 (src.x, src.y, src.z, 1) * m and dividing by 00563 // the result's third element. 00564 // 00565 // m.multDirMatrix(src,dst) multiplies src by the upper left 3x3 00566 // submatrix, ignoring the rest of matrix m. 00567 //----------------------------------------------------------------- 00568 00569 template <class S> 00570 void multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const; 00571 00572 template <class S> 00573 void multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const; 00574 00575 00576 //------------------------ 00577 // Component-wise division 00578 //------------------------ 00579 00580 const Matrix44 & operator /= (T a); 00581 Matrix44 operator / (T a) const; 00582 00583 00584 //------------------ 00585 // Transposed matrix 00586 //------------------ 00587 00588 const Matrix44 & transpose (); 00589 Matrix44 transposed () const; 00590 00591 00592 //------------------------------------------------------------ 00593 // Inverse matrix: If singExc is false, inverting a singular 00594 // matrix produces an identity matrix. If singExc is true, 00595 // inverting a singular matrix throws a SingMatrixExc. 00596 // 00597 // inverse() and invert() invert matrices using determinants; 00598 // gjInverse() and gjInvert() use the Gauss-Jordan method. 00599 // 00600 // inverse() and invert() are significantly faster than 00601 // gjInverse() and gjInvert(), but the results may be slightly 00602 // less accurate. 00603 // 00604 //------------------------------------------------------------ 00605 00606 const Matrix44 & invert (bool singExc = false) 00607 throw (Iex::MathExc); 00608 00609 Matrix44<T> inverse (bool singExc = false) const 00610 throw (Iex::MathExc); 00611 00612 const Matrix44 & gjInvert (bool singExc = false) 00613 throw (Iex::MathExc); 00614 00615 Matrix44<T> gjInverse (bool singExc = false) const 00616 throw (Iex::MathExc); 00617 00618 00619 //-------------------------------------------------------- 00620 // Set matrix to rotation by XYZ euler angles (in radians) 00621 //-------------------------------------------------------- 00622 00623 template <class S> 00624 const Matrix44 & setEulerAngles (const Vec3<S>& r); 00625 00626 00627 //-------------------------------------------------------- 00628 // Set matrix to rotation around given axis by given angle 00629 //-------------------------------------------------------- 00630 00631 template <class S> 00632 const Matrix44 & setAxisAngle (const Vec3<S>& ax, S ang); 00633 00634 00635 //------------------------------------------- 00636 // Rotate the matrix by XYZ euler angles in r 00637 //------------------------------------------- 00638 00639 template <class S> 00640 const Matrix44 & rotate (const Vec3<S> &r); 00641 00642 00643 //-------------------------------------------- 00644 // Set matrix to scale by given uniform factor 00645 //-------------------------------------------- 00646 00647 const Matrix44 & setScale (T s); 00648 00649 00650 //------------------------------------ 00651 // Set matrix to scale by given vector 00652 //------------------------------------ 00653 00654 template <class S> 00655 const Matrix44 & setScale (const Vec3<S> &s); 00656 00657 00658 //---------------------- 00659 // Scale the matrix by s 00660 //---------------------- 00661 00662 template <class S> 00663 const Matrix44 & scale (const Vec3<S> &s); 00664 00665 00666 //------------------------------------------ 00667 // Set matrix to translation by given vector 00668 //------------------------------------------ 00669 00670 template <class S> 00671 const Matrix44 & setTranslation (const Vec3<S> &t); 00672 00673 00674 //----------------------------- 00675 // Return translation component 00676 //----------------------------- 00677 00678 const Vec3<T> translation () const; 00679 00680 00681 //-------------------------- 00682 // Translate the matrix by t 00683 //-------------------------- 00684 00685 template <class S> 00686 const Matrix44 & translate (const Vec3<S> &t); 00687 00688 00689 //------------------------------------------------------------- 00690 // Set matrix to shear by given vector h. The resulting matrix 00691 // will shear x for each y coord. by a factor of h[0] ; 00692 // will shear x for each z coord. by a factor of h[1] ; 00693 // will shear y for each z coord. by a factor of h[2] . 00694 //------------------------------------------------------------- 00695 00696 template <class S> 00697 const Matrix44 & setShear (const Vec3<S> &h); 00698 00699 00700 //------------------------------------------------------------ 00701 // Set matrix to shear by given factors. The resulting matrix 00702 // will shear x for each y coord. by a factor of h.xy ; 00703 // will shear x for each z coord. by a factor of h.xz ; 00704 // will shear y for each z coord. by a factor of h.yz ; 00705 // will shear y for each x coord. by a factor of h.yx ; 00706 // will shear z for each x coord. by a factor of h.zx ; 00707 // will shear z for each y coord. by a factor of h.zy . 00708 //------------------------------------------------------------ 00709 00710 template <class S> 00711 const Matrix44 & setShear (const Shear6<S> &h); 00712 00713 00714 //-------------------------------------------------------- 00715 // Shear the matrix by given vector. The composed matrix 00716 // will be <shear> * <this>, where the shear matrix ... 00717 // will shear x for each y coord. by a factor of h[0] ; 00718 // will shear x for each z coord. by a factor of h[1] ; 00719 // will shear y for each z coord. by a factor of h[2] . 00720 //-------------------------------------------------------- 00721 00722 template <class S> 00723 const Matrix44 & shear (const Vec3<S> &h); 00724 00725 00726 //------------------------------------------------------------ 00727 // Shear the matrix by the given factors. The composed matrix 00728 // will be <shear> * <this>, where the shear matrix ... 00729 // will shear x for each y coord. by a factor of h.xy ; 00730 // will shear x for each z coord. by a factor of h.xz ; 00731 // will shear y for each z coord. by a factor of h.yz ; 00732 // will shear y for each x coord. by a factor of h.yx ; 00733 // will shear z for each x coord. by a factor of h.zx ; 00734 // will shear z for each y coord. by a factor of h.zy . 00735 //------------------------------------------------------------ 00736 00737 template <class S> 00738 const Matrix44 & shear (const Shear6<S> &h); 00739 00740 00741 //------------------------------------------------- 00742 // Limitations of type T (see also class limits<T>) 00743 //------------------------------------------------- 00744 00745 static T baseTypeMin() {return limits<T>::min();} 00746 static T baseTypeMax() {return limits<T>::max();} 00747 static T baseTypeSmallest() {return limits<T>::smallest();} 00748 static T baseTypeEpsilon() {return limits<T>::epsilon();} 00749 00750 private: 00751 00752 template <typename R, typename S> 00753 struct isSameType 00754 { 00755 enum {value = 0}; 00756 }; 00757 00758 template <typename R> 00759 struct isSameType<R, R> 00760 { 00761 enum {value = 1}; 00762 }; 00763 }; 00764 00765 00766 //-------------- 00767 // Stream output 00768 //-------------- 00769 00770 template <class T> 00771 std::ostream & operator << (std::ostream & s, const Matrix33<T> &m); 00772 00773 template <class T> 00774 std::ostream & operator << (std::ostream & s, const Matrix44<T> &m); 00775 00776 00777 //--------------------------------------------- 00778 // Vector-times-matrix multiplication operators 00779 //--------------------------------------------- 00780 00781 template <class S, class T> 00782 const Vec2<S> & operator *= (Vec2<S> &v, const Matrix33<T> &m); 00783 00784 template <class S, class T> 00785 Vec2<S> operator * (const Vec2<S> &v, const Matrix33<T> &m); 00786 00787 template <class S, class T> 00788 const Vec3<S> & operator *= (Vec3<S> &v, const Matrix33<T> &m); 00789 00790 template <class S, class T> 00791 Vec3<S> operator * (const Vec3<S> &v, const Matrix33<T> &m); 00792 00793 template <class S, class T> 00794 const Vec3<S> & operator *= (Vec3<S> &v, const Matrix44<T> &m); 00795 00796 template <class S, class T> 00797 Vec3<S> operator * (const Vec3<S> &v, const Matrix44<T> &m); 00798 00799 template <class S, class T> 00800 const Vec4<S> & operator *= (Vec4<S> &v, const Matrix44<T> &m); 00801 00802 template <class S, class T> 00803 Vec4<S> operator * (const Vec4<S> &v, const Matrix44<T> &m); 00804 00805 //------------------------- 00806 // Typedefs for convenience 00807 //------------------------- 00808 00809 typedef Matrix33 <float> M33f; 00810 typedef Matrix33 <double> M33d; 00811 typedef Matrix44 <float> M44f; 00812 typedef Matrix44 <double> M44d; 00813 00814 00815 //--------------------------- 00816 // Implementation of Matrix33 00817 //--------------------------- 00818 00819 template <class T> 00820 inline T * 00821 Matrix33<T>::operator [] (int i) 00822 { 00823 return x[i]; 00824 } 00825 00826 template <class T> 00827 inline const T * 00828 Matrix33<T>::operator [] (int i) const 00829 { 00830 return x[i]; 00831 } 00832 00833 template <class T> 00834 inline 00835 Matrix33<T>::Matrix33 () 00836 { 00837 memset (x, 0, sizeof (x)); 00838 x[0][0] = 1; 00839 x[1][1] = 1; 00840 x[2][2] = 1; 00841 } 00842 00843 template <class T> 00844 inline 00845 Matrix33<T>::Matrix33 (T a) 00846 { 00847 x[0][0] = a; 00848 x[0][1] = a; 00849 x[0][2] = a; 00850 x[1][0] = a; 00851 x[1][1] = a; 00852 x[1][2] = a; 00853 x[2][0] = a; 00854 x[2][1] = a; 00855 x[2][2] = a; 00856 } 00857 00858 template <class T> 00859 inline 00860 Matrix33<T>::Matrix33 (const T a[3][3]) 00861 { 00862 memcpy (x, a, sizeof (x)); 00863 } 00864 00865 template <class T> 00866 inline 00867 Matrix33<T>::Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i) 00868 { 00869 x[0][0] = a; 00870 x[0][1] = b; 00871 x[0][2] = c; 00872 x[1][0] = d; 00873 x[1][1] = e; 00874 x[1][2] = f; 00875 x[2][0] = g; 00876 x[2][1] = h; 00877 x[2][2] = i; 00878 } 00879 00880 template <class T> 00881 inline 00882 Matrix33<T>::Matrix33 (const Matrix33 &v) 00883 { 00884 memcpy (x, v.x, sizeof (x)); 00885 } 00886 00887 template <class T> 00888 template <class S> 00889 inline 00890 Matrix33<T>::Matrix33 (const Matrix33<S> &v) 00891 { 00892 x[0][0] = T (v.x[0][0]); 00893 x[0][1] = T (v.x[0][1]); 00894 x[0][2] = T (v.x[0][2]); 00895 x[1][0] = T (v.x[1][0]); 00896 x[1][1] = T (v.x[1][1]); 00897 x[1][2] = T (v.x[1][2]); 00898 x[2][0] = T (v.x[2][0]); 00899 x[2][1] = T (v.x[2][1]); 00900 x[2][2] = T (v.x[2][2]); 00901 } 00902 00903 template <class T> 00904 inline const Matrix33<T> & 00905 Matrix33<T>::operator = (const Matrix33 &v) 00906 { 00907 memcpy (x, v.x, sizeof (x)); 00908 return *this; 00909 } 00910 00911 template <class T> 00912 inline const Matrix33<T> & 00913 Matrix33<T>::operator = (T a) 00914 { 00915 x[0][0] = a; 00916 x[0][1] = a; 00917 x[0][2] = a; 00918 x[1][0] = a; 00919 x[1][1] = a; 00920 x[1][2] = a; 00921 x[2][0] = a; 00922 x[2][1] = a; 00923 x[2][2] = a; 00924 return *this; 00925 } 00926 00927 template <class T> 00928 inline T * 00929 Matrix33<T>::getValue () 00930 { 00931 return (T *) &x[0][0]; 00932 } 00933 00934 template <class T> 00935 inline const T * 00936 Matrix33<T>::getValue () const 00937 { 00938 return (const T *) &x[0][0]; 00939 } 00940 00941 template <class T> 00942 template <class S> 00943 inline void 00944 Matrix33<T>::getValue (Matrix33<S> &v) const 00945 { 00946 if (isSameType<S,T>::value) 00947 { 00948 memcpy (v.x, x, sizeof (x)); 00949 } 00950 else 00951 { 00952 v.x[0][0] = x[0][0]; 00953 v.x[0][1] = x[0][1]; 00954 v.x[0][2] = x[0][2]; 00955 v.x[1][0] = x[1][0]; 00956 v.x[1][1] = x[1][1]; 00957 v.x[1][2] = x[1][2]; 00958 v.x[2][0] = x[2][0]; 00959 v.x[2][1] = x[2][1]; 00960 v.x[2][2] = x[2][2]; 00961 } 00962 } 00963 00964 template <class T> 00965 template <class S> 00966 inline Matrix33<T> & 00967 Matrix33<T>::setValue (const Matrix33<S> &v) 00968 { 00969 if (isSameType<S,T>::value) 00970 { 00971 memcpy (x, v.x, sizeof (x)); 00972 } 00973 else 00974 { 00975 x[0][0] = v.x[0][0]; 00976 x[0][1] = v.x[0][1]; 00977 x[0][2] = v.x[0][2]; 00978 x[1][0] = v.x[1][0]; 00979 x[1][1] = v.x[1][1]; 00980 x[1][2] = v.x[1][2]; 00981 x[2][0] = v.x[2][0]; 00982 x[2][1] = v.x[2][1]; 00983 x[2][2] = v.x[2][2]; 00984 } 00985 00986 return *this; 00987 } 00988 00989 template <class T> 00990 template <class S> 00991 inline Matrix33<T> & 00992 Matrix33<T>::setTheMatrix (const Matrix33<S> &v) 00993 { 00994 if (isSameType<S,T>::value) 00995 { 00996 memcpy (x, v.x, sizeof (x)); 00997 } 00998 else 00999 { 01000 x[0][0] = v.x[0][0]; 01001 x[0][1] = v.x[0][1]; 01002 x[0][2] = v.x[0][2]; 01003 x[1][0] = v.x[1][0]; 01004 x[1][1] = v.x[1][1]; 01005 x[1][2] = v.x[1][2]; 01006 x[2][0] = v.x[2][0]; 01007 x[2][1] = v.x[2][1]; 01008 x[2][2] = v.x[2][2]; 01009 } 01010 01011 return *this; 01012 } 01013 01014 template <class T> 01015 inline void 01016 Matrix33<T>::makeIdentity() 01017 { 01018 memset (x, 0, sizeof (x)); 01019 x[0][0] = 1; 01020 x[1][1] = 1; 01021 x[2][2] = 1; 01022 } 01023 01024 template <class T> 01025 bool 01026 Matrix33<T>::operator == (const Matrix33 &v) const 01027 { 01028 return x[0][0] == v.x[0][0] && 01029 x[0][1] == v.x[0][1] && 01030 x[0][2] == v.x[0][2] && 01031 x[1][0] == v.x[1][0] && 01032 x[1][1] == v.x[1][1] && 01033 x[1][2] == v.x[1][2] && 01034 x[2][0] == v.x[2][0] && 01035 x[2][1] == v.x[2][1] && 01036 x[2][2] == v.x[2][2]; 01037 } 01038 01039 template <class T> 01040 bool 01041 Matrix33<T>::operator != (const Matrix33 &v) const 01042 { 01043 return x[0][0] != v.x[0][0] || 01044 x[0][1] != v.x[0][1] || 01045 x[0][2] != v.x[0][2] || 01046 x[1][0] != v.x[1][0] || 01047 x[1][1] != v.x[1][1] || 01048 x[1][2] != v.x[1][2] || 01049 x[2][0] != v.x[2][0] || 01050 x[2][1] != v.x[2][1] || 01051 x[2][2] != v.x[2][2]; 01052 } 01053 01054 template <class T> 01055 bool 01056 Matrix33<T>::equalWithAbsError (const Matrix33<T> &m, T e) const 01057 { 01058 for (int i = 0; i < 3; i++) 01059 for (int j = 0; j < 3; j++) 01060 if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e)) 01061 return false; 01062 01063 return true; 01064 } 01065 01066 template <class T> 01067 bool 01068 Matrix33<T>::equalWithRelError (const Matrix33<T> &m, T e) const 01069 { 01070 for (int i = 0; i < 3; i++) 01071 for (int j = 0; j < 3; j++) 01072 if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e)) 01073 return false; 01074 01075 return true; 01076 } 01077 01078 template <class T> 01079 const Matrix33<T> & 01080 Matrix33<T>::operator += (const Matrix33<T> &v) 01081 { 01082 x[0][0] += v.x[0][0]; 01083 x[0][1] += v.x[0][1]; 01084 x[0][2] += v.x[0][2]; 01085 x[1][0] += v.x[1][0]; 01086 x[1][1] += v.x[1][1]; 01087 x[1][2] += v.x[1][2]; 01088 x[2][0] += v.x[2][0]; 01089 x[2][1] += v.x[2][1]; 01090 x[2][2] += v.x[2][2]; 01091 01092 return *this; 01093 } 01094 01095 template <class T> 01096 const Matrix33<T> & 01097 Matrix33<T>::operator += (T a) 01098 { 01099 x[0][0] += a; 01100 x[0][1] += a; 01101 x[0][2] += a; 01102 x[1][0] += a; 01103 x[1][1] += a; 01104 x[1][2] += a; 01105 x[2][0] += a; 01106 x[2][1] += a; 01107 x[2][2] += a; 01108 01109 return *this; 01110 } 01111 01112 template <class T> 01113 Matrix33<T> 01114 Matrix33<T>::operator + (const Matrix33<T> &v) const 01115 { 01116 return Matrix33 (x[0][0] + v.x[0][0], 01117 x[0][1] + v.x[0][1], 01118 x[0][2] + v.x[0][2], 01119 x[1][0] + v.x[1][0], 01120 x[1][1] + v.x[1][1], 01121 x[1][2] + v.x[1][2], 01122 x[2][0] + v.x[2][0], 01123 x[2][1] + v.x[2][1], 01124 x[2][2] + v.x[2][2]); 01125 } 01126 01127 template <class T> 01128 const Matrix33<T> & 01129 Matrix33<T>::operator -= (const Matrix33<T> &v) 01130 { 01131 x[0][0] -= v.x[0][0]; 01132 x[0][1] -= v.x[0][1]; 01133 x[0][2] -= v.x[0][2]; 01134 x[1][0] -= v.x[1][0]; 01135 x[1][1] -= v.x[1][1]; 01136 x[1][2] -= v.x[1][2]; 01137 x[2][0] -= v.x[2][0]; 01138 x[2][1] -= v.x[2][1]; 01139 x[2][2] -= v.x[2][2]; 01140 01141 return *this; 01142 } 01143 01144 template <class T> 01145 const Matrix33<T> & 01146 Matrix33<T>::operator -= (T a) 01147 { 01148 x[0][0] -= a; 01149 x[0][1] -= a; 01150 x[0][2] -= a; 01151 x[1][0] -= a; 01152 x[1][1] -= a; 01153 x[1][2] -= a; 01154 x[2][0] -= a; 01155 x[2][1] -= a; 01156 x[2][2] -= a; 01157 01158 return *this; 01159 } 01160 01161 template <class T> 01162 Matrix33<T> 01163 Matrix33<T>::operator - (const Matrix33<T> &v) const 01164 { 01165 return Matrix33 (x[0][0] - v.x[0][0], 01166 x[0][1] - v.x[0][1], 01167 x[0][2] - v.x[0][2], 01168 x[1][0] - v.x[1][0], 01169 x[1][1] - v.x[1][1], 01170 x[1][2] - v.x[1][2], 01171 x[2][0] - v.x[2][0], 01172 x[2][1] - v.x[2][1], 01173 x[2][2] - v.x[2][2]); 01174 } 01175 01176 template <class T> 01177 Matrix33<T> 01178 Matrix33<T>::operator - () const 01179 { 01180 return Matrix33 (-x[0][0], 01181 -x[0][1], 01182 -x[0][2], 01183 -x[1][0], 01184 -x[1][1], 01185 -x[1][2], 01186 -x[2][0], 01187 -x[2][1], 01188 -x[2][2]); 01189 } 01190 01191 template <class T> 01192 const Matrix33<T> & 01193 Matrix33<T>::negate () 01194 { 01195 x[0][0] = -x[0][0]; 01196 x[0][1] = -x[0][1]; 01197 x[0][2] = -x[0][2]; 01198 x[1][0] = -x[1][0]; 01199 x[1][1] = -x[1][1]; 01200 x[1][2] = -x[1][2]; 01201 x[2][0] = -x[2][0]; 01202 x[2][1] = -x[2][1]; 01203 x[2][2] = -x[2][2]; 01204 01205 return *this; 01206 } 01207 01208 template <class T> 01209 const Matrix33<T> & 01210 Matrix33<T>::operator *= (T a) 01211 { 01212 x[0][0] *= a; 01213 x[0][1] *= a; 01214 x[0][2] *= a; 01215 x[1][0] *= a; 01216 x[1][1] *= a; 01217 x[1][2] *= a; 01218 x[2][0] *= a; 01219 x[2][1] *= a; 01220 x[2][2] *= a; 01221 01222 return *this; 01223 } 01224 01225 template <class T> 01226 Matrix33<T> 01227 Matrix33<T>::operator * (T a) const 01228 { 01229 return Matrix33 (x[0][0] * a, 01230 x[0][1] * a, 01231 x[0][2] * a, 01232 x[1][0] * a, 01233 x[1][1] * a, 01234 x[1][2] * a, 01235 x[2][0] * a, 01236 x[2][1] * a, 01237 x[2][2] * a); 01238 } 01239 01240 template <class T> 01241 inline Matrix33<T> 01242 operator * (T a, const Matrix33<T> &v) 01243 { 01244 return v * a; 01245 } 01246 01247 template <class T> 01248 const Matrix33<T> & 01249 Matrix33<T>::operator *= (const Matrix33<T> &v) 01250 { 01251 Matrix33 tmp (T (0)); 01252 01253 for (int i = 0; i < 3; i++) 01254 for (int j = 0; j < 3; j++) 01255 for (int k = 0; k < 3; k++) 01256 tmp.x[i][j] += x[i][k] * v.x[k][j]; 01257 01258 *this = tmp; 01259 return *this; 01260 } 01261 01262 template <class T> 01263 Matrix33<T> 01264 Matrix33<T>::operator * (const Matrix33<T> &v) const 01265 { 01266 Matrix33 tmp (T (0)); 01267 01268 for (int i = 0; i < 3; i++) 01269 for (int j = 0; j < 3; j++) 01270 for (int k = 0; k < 3; k++) 01271 tmp.x[i][j] += x[i][k] * v.x[k][j]; 01272 01273 return tmp; 01274 } 01275 01276 template <class T> 01277 template <class S> 01278 void 01279 Matrix33<T>::multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const 01280 { 01281 S a, b, w; 01282 01283 a = src[0] * x[0][0] + src[1] * x[1][0] + x[2][0]; 01284 b = src[0] * x[0][1] + src[1] * x[1][1] + x[2][1]; 01285 w = src[0] * x[0][2] + src[1] * x[1][2] + x[2][2]; 01286 01287 dst.x = a / w; 01288 dst.y = b / w; 01289 } 01290 01291 template <class T> 01292 template <class S> 01293 void 01294 Matrix33<T>::multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const 01295 { 01296 S a, b; 01297 01298 a = src[0] * x[0][0] + src[1] * x[1][0]; 01299 b = src[0] * x[0][1] + src[1] * x[1][1]; 01300 01301 dst.x = a; 01302 dst.y = b; 01303 } 01304 01305 template <class T> 01306 const Matrix33<T> & 01307 Matrix33<T>::operator /= (T a) 01308 { 01309 x[0][0] /= a; 01310 x[0][1] /= a; 01311 x[0][2] /= a; 01312 x[1][0] /= a; 01313 x[1][1] /= a; 01314 x[1][2] /= a; 01315 x[2][0] /= a; 01316 x[2][1] /= a; 01317 x[2][2] /= a; 01318 01319 return *this; 01320 } 01321 01322 template <class T> 01323 Matrix33<T> 01324 Matrix33<T>::operator / (T a) const 01325 { 01326 return Matrix33 (x[0][0] / a, 01327 x[0][1] / a, 01328 x[0][2] / a, 01329 x[1][0] / a, 01330 x[1][1] / a, 01331 x[1][2] / a, 01332 x[2][0] / a, 01333 x[2][1] / a, 01334 x[2][2] / a); 01335 } 01336 01337 template <class T> 01338 const Matrix33<T> & 01339 Matrix33<T>::transpose () 01340 { 01341 Matrix33 tmp (x[0][0], 01342 x[1][0], 01343 x[2][0], 01344 x[0][1], 01345 x[1][1], 01346 x[2][1], 01347 x[0][2], 01348 x[1][2], 01349 x[2][2]); 01350 *this = tmp; 01351 return *this; 01352 } 01353 01354 template <class T> 01355 Matrix33<T> 01356 Matrix33<T>::transposed () const 01357 { 01358 return Matrix33 (x[0][0], 01359 x[1][0], 01360 x[2][0], 01361 x[0][1], 01362 x[1][1], 01363 x[2][1], 01364 x[0][2], 01365 x[1][2], 01366 x[2][2]); 01367 } 01368 01369 template <class T> 01370 const Matrix33<T> & 01371 Matrix33<T>::gjInvert (bool singExc) throw (Iex::MathExc) 01372 { 01373 *this = gjInverse (singExc); 01374 return *this; 01375 } 01376 01377 template <class T> 01378 Matrix33<T> 01379 Matrix33<T>::gjInverse (bool singExc) const throw (Iex::MathExc) 01380 { 01381 int i, j, k; 01382 Matrix33 s; 01383 Matrix33 t (*this); 01384 01385 // Forward elimination 01386 01387 for (i = 0; i < 2 ; i++) 01388 { 01389 int pivot = i; 01390 01391 T pivotsize = t[i][i]; 01392 01393 if (pivotsize < 0) 01394 pivotsize = -pivotsize; 01395 01396 for (j = i + 1; j < 3; j++) 01397 { 01398 T tmp = t[j][i]; 01399 01400 if (tmp < 0) 01401 tmp = -tmp; 01402 01403 if (tmp > pivotsize) 01404 { 01405 pivot = j; 01406 pivotsize = tmp; 01407 } 01408 } 01409 01410 if (pivotsize == 0) 01411 { 01412 if (singExc) 01413 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix."); 01414 01415 return Matrix33(); 01416 } 01417 01418 if (pivot != i) 01419 { 01420 for (j = 0; j < 3; j++) 01421 { 01422 T tmp; 01423 01424 tmp = t[i][j]; 01425 t[i][j] = t[pivot][j]; 01426 t[pivot][j] = tmp; 01427 01428 tmp = s[i][j]; 01429 s[i][j] = s[pivot][j]; 01430 s[pivot][j] = tmp; 01431 } 01432 } 01433 01434 for (j = i + 1; j < 3; j++) 01435 { 01436 T f = t[j][i] / t[i][i]; 01437 01438 for (k = 0; k < 3; k++) 01439 { 01440 t[j][k] -= f * t[i][k]; 01441 s[j][k] -= f * s[i][k]; 01442 } 01443 } 01444 } 01445 01446 // Backward substitution 01447 01448 for (i = 2; i >= 0; --i) 01449 { 01450 T f; 01451 01452 if ((f = t[i][i]) == 0) 01453 { 01454 if (singExc) 01455 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix."); 01456 01457 return Matrix33(); 01458 } 01459 01460 for (j = 0; j < 3; j++) 01461 { 01462 t[i][j] /= f; 01463 s[i][j] /= f; 01464 } 01465 01466 for (j = 0; j < i; j++) 01467 { 01468 f = t[j][i]; 01469 01470 for (k = 0; k < 3; k++) 01471 { 01472 t[j][k] -= f * t[i][k]; 01473 s[j][k] -= f * s[i][k]; 01474 } 01475 } 01476 } 01477 01478 return s; 01479 } 01480 01481 template <class T> 01482 const Matrix33<T> & 01483 Matrix33<T>::invert (bool singExc) throw (Iex::MathExc) 01484 { 01485 *this = inverse (singExc); 01486 return *this; 01487 } 01488 01489 template <class T> 01490 Matrix33<T> 01491 Matrix33<T>::inverse (bool singExc) const throw (Iex::MathExc) 01492 { 01493 if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1) 01494 { 01495 Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2], 01496 x[2][1] * x[0][2] - x[0][1] * x[2][2], 01497 x[0][1] * x[1][2] - x[1][1] * x[0][2], 01498 01499 x[2][0] * x[1][2] - x[1][0] * x[2][2], 01500 x[0][0] * x[2][2] - x[2][0] * x[0][2], 01501 x[1][0] * x[0][2] - x[0][0] * x[1][2], 01502 01503 x[1][0] * x[2][1] - x[2][0] * x[1][1], 01504 x[2][0] * x[0][1] - x[0][0] * x[2][1], 01505 x[0][0] * x[1][1] - x[1][0] * x[0][1]); 01506 01507 T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0]; 01508 01509 if (Imath::abs (r) >= 1) 01510 { 01511 for (int i = 0; i < 3; ++i) 01512 { 01513 for (int j = 0; j < 3; ++j) 01514 { 01515 s[i][j] /= r; 01516 } 01517 } 01518 } 01519 else 01520 { 01521 T mr = Imath::abs (r) / limits<T>::smallest(); 01522 01523 for (int i = 0; i < 3; ++i) 01524 { 01525 for (int j = 0; j < 3; ++j) 01526 { 01527 if (mr > Imath::abs (s[i][j])) 01528 { 01529 s[i][j] /= r; 01530 } 01531 else 01532 { 01533 if (singExc) 01534 throw SingMatrixExc ("Cannot invert " 01535 "singular matrix."); 01536 return Matrix33(); 01537 } 01538 } 01539 } 01540 } 01541 01542 return s; 01543 } 01544 else 01545 { 01546 Matrix33 s ( x[1][1], 01547 -x[0][1], 01548 0, 01549 01550 -x[1][0], 01551 x[0][0], 01552 0, 01553 01554 0, 01555 0, 01556 1); 01557 01558 T r = x[0][0] * x[1][1] - x[1][0] * x[0][1]; 01559 01560 if (Imath::abs (r) >= 1) 01561 { 01562 for (int i = 0; i < 2; ++i) 01563 { 01564 for (int j = 0; j < 2; ++j) 01565 { 01566 s[i][j] /= r; 01567 } 01568 } 01569 } 01570 else 01571 { 01572 T mr = Imath::abs (r) / limits<T>::smallest(); 01573 01574 for (int i = 0; i < 2; ++i) 01575 { 01576 for (int j = 0; j < 2; ++j) 01577 { 01578 if (mr > Imath::abs (s[i][j])) 01579 { 01580 s[i][j] /= r; 01581 } 01582 else 01583 { 01584 if (singExc) 01585 throw SingMatrixExc ("Cannot invert " 01586 "singular matrix."); 01587 return Matrix33(); 01588 } 01589 } 01590 } 01591 } 01592 01593 s[2][0] = -x[2][0] * s[0][0] - x[2][1] * s[1][0]; 01594 s[2][1] = -x[2][0] * s[0][1] - x[2][1] * s[1][1]; 01595 01596 return s; 01597 } 01598 } 01599 01600 template <class T> 01601 template <class S> 01602 const Matrix33<T> & 01603 Matrix33<T>::setRotation (S r) 01604 { 01605 S cos_r, sin_r; 01606 01607 cos_r = Math<T>::cos (r); 01608 sin_r = Math<T>::sin (r); 01609 01610 x[0][0] = cos_r; 01611 x[0][1] = sin_r; 01612 x[0][2] = 0; 01613 01614 x[1][0] = -sin_r; 01615 x[1][1] = cos_r; 01616 x[1][2] = 0; 01617 01618 x[2][0] = 0; 01619 x[2][1] = 0; 01620 x[2][2] = 1; 01621 01622 return *this; 01623 } 01624 01625 template <class T> 01626 template <class S> 01627 const Matrix33<T> & 01628 Matrix33<T>::rotate (S r) 01629 { 01630 *this *= Matrix33<T>().setRotation (r); 01631 return *this; 01632 } 01633 01634 template <class T> 01635 const Matrix33<T> & 01636 Matrix33<T>::setScale (T s) 01637 { 01638 memset (x, 0, sizeof (x)); 01639 x[0][0] = s; 01640 x[1][1] = s; 01641 x[2][2] = 1; 01642 01643 return *this; 01644 } 01645 01646 template <class T> 01647 template <class S> 01648 const Matrix33<T> & 01649 Matrix33<T>::setScale (const Vec2<S> &s) 01650 { 01651 memset (x, 0, sizeof (x)); 01652 x[0][0] = s[0]; 01653 x[1][1] = s[1]; 01654 x[2][2] = 1; 01655 01656 return *this; 01657 } 01658 01659 template <class T> 01660 template <class S> 01661 const Matrix33<T> & 01662 Matrix33<T>::scale (const Vec2<S> &s) 01663 { 01664 x[0][0] *= s[0]; 01665 x[0][1] *= s[0]; 01666 x[0][2] *= s[0]; 01667 01668 x[1][0] *= s[1]; 01669 x[1][1] *= s[1]; 01670 x[1][2] *= s[1]; 01671 01672 return *this; 01673 } 01674 01675 template <class T> 01676 template <class S> 01677 const Matrix33<T> & 01678 Matrix33<T>::setTranslation (const Vec2<S> &t) 01679 { 01680 x[0][0] = 1; 01681 x[0][1] = 0; 01682 x[0][2] = 0; 01683 01684 x[1][0] = 0; 01685 x[1][1] = 1; 01686 x[1][2] = 0; 01687 01688 x[2][0] = t[0]; 01689 x[2][1] = t[1]; 01690 x[2][2] = 1; 01691 01692 return *this; 01693 } 01694 01695 template <class T> 01696 inline Vec2<T> 01697 Matrix33<T>::translation () const 01698 { 01699 return Vec2<T> (x[2][0], x[2][1]); 01700 } 01701 01702 template <class T> 01703 template <class S> 01704 const Matrix33<T> & 01705 Matrix33<T>::translate (const Vec2<S> &t) 01706 { 01707 x[2][0] += t[0] * x[0][0] + t[1] * x[1][0]; 01708 x[2][1] += t[0] * x[0][1] + t[1] * x[1][1]; 01709 x[2][2] += t[0] * x[0][2] + t[1] * x[1][2]; 01710 01711 return *this; 01712 } 01713 01714 template <class T> 01715 template <class S> 01716 const Matrix33<T> & 01717 Matrix33<T>::setShear (const S &xy) 01718 { 01719 x[0][0] = 1; 01720 x[0][1] = 0; 01721 x[0][2] = 0; 01722 01723 x[1][0] = xy; 01724 x[1][1] = 1; 01725 x[1][2] = 0; 01726 01727 x[2][0] = 0; 01728 x[2][1] = 0; 01729 x[2][2] = 1; 01730 01731 return *this; 01732 } 01733 01734 template <class T> 01735 template <class S> 01736 const Matrix33<T> & 01737 Matrix33<T>::setShear (const Vec2<S> &h) 01738 { 01739 x[0][0] = 1; 01740 x[0][1] = h[1]; 01741 x[0][2] = 0; 01742 01743 x[1][0] = h[0]; 01744 x[1][1] = 1; 01745 x[1][2] = 0; 01746 01747 x[2][0] = 0; 01748 x[2][1] = 0; 01749 x[2][2] = 1; 01750 01751 return *this; 01752 } 01753 01754 template <class T> 01755 template <class S> 01756 const Matrix33<T> & 01757 Matrix33<T>::shear (const S &xy) 01758 { 01759 // 01760 // In this case, we don't need a temp. copy of the matrix 01761 // because we never use a value on the RHS after we've 01762 // changed it on the LHS. 01763 // 01764 01765 x[1][0] += xy * x[0][0]; 01766 x[1][1] += xy * x[0][1]; 01767 x[1][2] += xy * x[0][2]; 01768 01769 return *this; 01770 } 01771 01772 template <class T> 01773 template <class S> 01774 const Matrix33<T> & 01775 Matrix33<T>::shear (const Vec2<S> &h) 01776 { 01777 Matrix33<T> P (*this); 01778 01779 x[0][0] = P[0][0] + h[1] * P[1][0]; 01780 x[0][1] = P[0][1] + h[1] * P[1][1]; 01781 x[0][2] = P[0][2] + h[1] * P[1][2]; 01782 01783 x[1][0] = P[1][0] + h[0] * P[0][0]; 01784 x[1][1] = P[1][1] + h[0] * P[0][1]; 01785 x[1][2] = P[1][2] + h[0] * P[0][2]; 01786 01787 return *this; 01788 } 01789 01790 01791 //--------------------------- 01792 // Implementation of Matrix44 01793 //--------------------------- 01794 01795 template <class T> 01796 inline T * 01797 Matrix44<T>::operator [] (int i) 01798 { 01799 return x[i]; 01800 } 01801 01802 template <class T> 01803 inline const T * 01804 Matrix44<T>::operator [] (int i) const 01805 { 01806 return x[i]; 01807 } 01808 01809 template <class T> 01810 inline 01811 Matrix44<T>::Matrix44 () 01812 { 01813 memset (x, 0, sizeof (x)); 01814 x[0][0] = 1; 01815 x[1][1] = 1; 01816 x[2][2] = 1; 01817 x[3][3] = 1; 01818 } 01819 01820 template <class T> 01821 inline 01822 Matrix44<T>::Matrix44 (T a) 01823 { 01824 x[0][0] = a; 01825 x[0][1] = a; 01826 x[0][2] = a; 01827 x[0][3] = a; 01828 x[1][0] = a; 01829 x[1][1] = a; 01830 x[1][2] = a; 01831 x[1][3] = a; 01832 x[2][0] = a; 01833 x[2][1] = a; 01834 x[2][2] = a; 01835 x[2][3] = a; 01836 x[3][0] = a; 01837 x[3][1] = a; 01838 x[3][2] = a; 01839 x[3][3] = a; 01840 } 01841 01842 template <class T> 01843 inline 01844 Matrix44<T>::Matrix44 (const T a[4][4]) 01845 { 01846 memcpy (x, a, sizeof (x)); 01847 } 01848 01849 template <class T> 01850 inline 01851 Matrix44<T>::Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h, 01852 T i, T j, T k, T l, T m, T n, T o, T p) 01853 { 01854 x[0][0] = a; 01855 x[0][1] = b; 01856 x[0][2] = c; 01857 x[0][3] = d; 01858 x[1][0] = e; 01859 x[1][1] = f; 01860 x[1][2] = g; 01861 x[1][3] = h; 01862 x[2][0] = i; 01863 x[2][1] = j; 01864 x[2][2] = k; 01865 x[2][3] = l; 01866 x[3][0] = m; 01867 x[3][1] = n; 01868 x[3][2] = o; 01869 x[3][3] = p; 01870 } 01871 01872 01873 template <class T> 01874 inline 01875 Matrix44<T>::Matrix44 (Matrix33<T> r, Vec3<T> t) 01876 { 01877 x[0][0] = r[0][0]; 01878 x[0][1] = r[0][1]; 01879 x[0][2] = r[0][2]; 01880 x[0][3] = 0; 01881 x[1][0] = r[1][0]; 01882 x[1][1] = r[1][1]; 01883 x[1][2] = r[1][2]; 01884 x[1][3] = 0; 01885 x[2][0] = r[2][0]; 01886 x[2][1] = r[2][1]; 01887 x[2][2] = r[2][2]; 01888 x[2][3] = 0; 01889 x[3][0] = t[0]; 01890 x[3][1] = t[1]; 01891 x[3][2] = t[2]; 01892 x[3][3] = 1; 01893 } 01894 01895 template <class T> 01896 inline 01897 Matrix44<T>::Matrix44 (const Matrix44 &v) 01898 { 01899 x[0][0] = v.x[0][0]; 01900 x[0][1] = v.x[0][1]; 01901 x[0][2] = v.x[0][2]; 01902 x[0][3] = v.x[0][3]; 01903 x[1][0] = v.x[1][0]; 01904 x[1][1] = v.x[1][1]; 01905 x[1][2] = v.x[1][2]; 01906 x[1][3] = v.x[1][3]; 01907 x[2][0] = v.x[2][0]; 01908 x[2][1] = v.x[2][1]; 01909 x[2][2] = v.x[2][2]; 01910 x[2][3] = v.x[2][3]; 01911 x[3][0] = v.x[3][0]; 01912 x[3][1] = v.x[3][1]; 01913 x[3][2] = v.x[3][2]; 01914 x[3][3] = v.x[3][3]; 01915 } 01916 01917 template <class T> 01918 template <class S> 01919 inline 01920 Matrix44<T>::Matrix44 (const Matrix44<S> &v) 01921 { 01922 x[0][0] = T (v.x[0][0]); 01923 x[0][1] = T (v.x[0][1]); 01924 x[0][2] = T (v.x[0][2]); 01925 x[0][3] = T (v.x[0][3]); 01926 x[1][0] = T (v.x[1][0]); 01927 x[1][1] = T (v.x[1][1]); 01928 x[1][2] = T (v.x[1][2]); 01929 x[1][3] = T (v.x[1][3]); 01930 x[2][0] = T (v.x[2][0]); 01931 x[2][1] = T (v.x[2][1]); 01932 x[2][2] = T (v.x[2][2]); 01933 x[2][3] = T (v.x[2][3]); 01934 x[3][0] = T (v.x[3][0]); 01935 x[3][1] = T (v.x[3][1]); 01936 x[3][2] = T (v.x[3][2]); 01937 x[3][3] = T (v.x[3][3]); 01938 } 01939 01940 template <class T> 01941 inline const Matrix44<T> & 01942 Matrix44<T>::operator = (const Matrix44 &v) 01943 { 01944 x[0][0] = v.x[0][0]; 01945 x[0][1] = v.x[0][1]; 01946 x[0][2] = v.x[0][2]; 01947 x[0][3] = v.x[0][3]; 01948 x[1][0] = v.x[1][0]; 01949 x[1][1] = v.x[1][1]; 01950 x[1][2] = v.x[1][2]; 01951 x[1][3] = v.x[1][3]; 01952 x[2][0] = v.x[2][0]; 01953 x[2][1] = v.x[2][1]; 01954 x[2][2] = v.x[2][2]; 01955 x[2][3] = v.x[2][3]; 01956 x[3][0] = v.x[3][0]; 01957 x[3][1] = v.x[3][1]; 01958 x[3][2] = v.x[3][2]; 01959 x[3][3] = v.x[3][3]; 01960 return *this; 01961 } 01962 01963 template <class T> 01964 inline const Matrix44<T> & 01965 Matrix44<T>::operator = (T a) 01966 { 01967 x[0][0] = a; 01968 x[0][1] = a; 01969 x[0][2] = a; 01970 x[0][3] = a; 01971 x[1][0] = a; 01972 x[1][1] = a; 01973 x[1][2] = a; 01974 x[1][3] = a; 01975 x[2][0] = a; 01976 x[2][1] = a; 01977 x[2][2] = a; 01978 x[2][3] = a; 01979 x[3][0] = a; 01980 x[3][1] = a; 01981 x[3][2] = a; 01982 x[3][3] = a; 01983 return *this; 01984 } 01985 01986 template <class T> 01987 inline T * 01988 Matrix44<T>::getValue () 01989 { 01990 return (T *) &x[0][0]; 01991 } 01992 01993 template <class T> 01994 inline const T * 01995 Matrix44<T>::getValue () const 01996 { 01997 return (const T *) &x[0][0]; 01998 } 01999 02000 template <class T> 02001 template <class S> 02002 inline void 02003 Matrix44<T>::getValue (Matrix44<S> &v) const 02004 { 02005 if (isSameType<S,T>::value) 02006 { 02007 memcpy (v.x, x, sizeof (x)); 02008 } 02009 else 02010 { 02011 v.x[0][0] = x[0][0]; 02012 v.x[0][1] = x[0][1]; 02013 v.x[0][2] = x[0][2]; 02014 v.x[0][3] = x[0][3]; 02015 v.x[1][0] = x[1][0]; 02016 v.x[1][1] = x[1][1]; 02017 v.x[1][2] = x[1][2]; 02018 v.x[1][3] = x[1][3]; 02019 v.x[2][0] = x[2][0]; 02020 v.x[2][1] = x[2][1]; 02021 v.x[2][2] = x[2][2]; 02022 v.x[2][3] = x[2][3]; 02023 v.x[3][0] = x[3][0]; 02024 v.x[3][1] = x[3][1]; 02025 v.x[3][2] = x[3][2]; 02026 v.x[3][3] = x[3][3]; 02027 } 02028 } 02029 02030 template <class T> 02031 template <class S> 02032 inline Matrix44<T> & 02033 Matrix44<T>::setValue (const Matrix44<S> &v) 02034 { 02035 if (isSameType<S,T>::value) 02036 { 02037 memcpy (x, v.x, sizeof (x)); 02038 } 02039 else 02040 { 02041 x[0][0] = v.x[0][0]; 02042 x[0][1] = v.x[0][1]; 02043 x[0][2] = v.x[0][2]; 02044 x[0][3] = v.x[0][3]; 02045 x[1][0] = v.x[1][0]; 02046 x[1][1] = v.x[1][1]; 02047 x[1][2] = v.x[1][2]; 02048 x[1][3] = v.x[1][3]; 02049 x[2][0] = v.x[2][0]; 02050 x[2][1] = v.x[2][1]; 02051 x[2][2] = v.x[2][2]; 02052 x[2][3] = v.x[2][3]; 02053 x[3][0] = v.x[3][0]; 02054 x[3][1] = v.x[3][1]; 02055 x[3][2] = v.x[3][2]; 02056 x[3][3] = v.x[3][3]; 02057 } 02058 02059 return *this; 02060 } 02061 02062 template <class T> 02063 template <class S> 02064 inline Matrix44<T> & 02065 Matrix44<T>::setTheMatrix (const Matrix44<S> &v) 02066 { 02067 if (isSameType<S,T>::value) 02068 { 02069 memcpy (x, v.x, sizeof (x)); 02070 } 02071 else 02072 { 02073 x[0][0] = v.x[0][0]; 02074 x[0][1] = v.x[0][1]; 02075 x[0][2] = v.x[0][2]; 02076 x[0][3] = v.x[0][3]; 02077 x[1][0] = v.x[1][0]; 02078 x[1][1] = v.x[1][1]; 02079 x[1][2] = v.x[1][2]; 02080 x[1][3] = v.x[1][3]; 02081 x[2][0] = v.x[2][0]; 02082 x[2][1] = v.x[2][1]; 02083 x[2][2] = v.x[2][2]; 02084 x[2][3] = v.x[2][3]; 02085 x[3][0] = v.x[3][0]; 02086 x[3][1] = v.x[3][1]; 02087 x[3][2] = v.x[3][2]; 02088 x[3][3] = v.x[3][3]; 02089 } 02090 02091 return *this; 02092 } 02093 02094 template <class T> 02095 inline void 02096 Matrix44<T>::makeIdentity() 02097 { 02098 memset (x, 0, sizeof (x)); 02099 x[0][0] = 1; 02100 x[1][1] = 1; 02101 x[2][2] = 1; 02102 x[3][3] = 1; 02103 } 02104 02105 template <class T> 02106 bool 02107 Matrix44<T>::operator == (const Matrix44 &v) const 02108 { 02109 return x[0][0] == v.x[0][0] && 02110 x[0][1] == v.x[0][1] && 02111 x[0][2] == v.x[0][2] && 02112 x[0][3] == v.x[0][3] && 02113 x[1][0] == v.x[1][0] && 02114 x[1][1] == v.x[1][1] && 02115 x[1][2] == v.x[1][2] && 02116 x[1][3] == v.x[1][3] && 02117 x[2][0] == v.x[2][0] && 02118 x[2][1] == v.x[2][1] && 02119 x[2][2] == v.x[2][2] && 02120 x[2][3] == v.x[2][3] && 02121 x[3][0] == v.x[3][0] && 02122 x[3][1] == v.x[3][1] && 02123 x[3][2] == v.x[3][2] && 02124 x[3][3] == v.x[3][3]; 02125 } 02126 02127 template <class T> 02128 bool 02129 Matrix44<T>::operator != (const Matrix44 &v) const 02130 { 02131 return x[0][0] != v.x[0][0] || 02132 x[0][1] != v.x[0][1] || 02133 x[0][2] != v.x[0][2] || 02134 x[0][3] != v.x[0][3] || 02135 x[1][0] != v.x[1][0] || 02136 x[1][1] != v.x[1][1] || 02137 x[1][2] != v.x[1][2] || 02138 x[1][3] != v.x[1][3] || 02139 x[2][0] != v.x[2][0] || 02140 x[2][1] != v.x[2][1] || 02141 x[2][2] != v.x[2][2] || 02142 x[2][3] != v.x[2][3] || 02143 x[3][0] != v.x[3][0] || 02144 x[3][1] != v.x[3][1] || 02145 x[3][2] != v.x[3][2] || 02146 x[3][3] != v.x[3][3]; 02147 } 02148 02149 template <class T> 02150 bool 02151 Matrix44<T>::equalWithAbsError (const Matrix44<T> &m, T e) const 02152 { 02153 for (int i = 0; i < 4; i++) 02154 for (int j = 0; j < 4; j++) 02155 if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e)) 02156 return false; 02157 02158 return true; 02159 } 02160 02161 template <class T> 02162 bool 02163 Matrix44<T>::equalWithRelError (const Matrix44<T> &m, T e) const 02164 { 02165 for (int i = 0; i < 4; i++) 02166 for (int j = 0; j < 4; j++) 02167 if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e)) 02168 return false; 02169 02170 return true; 02171 } 02172 02173 template <class T> 02174 const Matrix44<T> & 02175 Matrix44<T>::operator += (const Matrix44<T> &v) 02176 { 02177 x[0][0] += v.x[0][0]; 02178 x[0][1] += v.x[0][1]; 02179 x[0][2] += v.x[0][2]; 02180 x[0][3] += v.x[0][3]; 02181 x[1][0] += v.x[1][0]; 02182 x[1][1] += v.x[1][1]; 02183 x[1][2] += v.x[1][2]; 02184 x[1][3] += v.x[1][3]; 02185 x[2][0] += v.x[2][0]; 02186 x[2][1] += v.x[2][1]; 02187 x[2][2] += v.x[2][2]; 02188 x[2][3] += v.x[2][3]; 02189 x[3][0] += v.x[3][0]; 02190 x[3][1] += v.x[3][1]; 02191 x[3][2] += v.x[3][2]; 02192 x[3][3] += v.x[3][3]; 02193 02194 return *this; 02195 } 02196 02197 template <class T> 02198 const Matrix44<T> & 02199 Matrix44<T>::operator += (T a) 02200 { 02201 x[0][0] += a; 02202 x[0][1] += a; 02203 x[0][2] += a; 02204 x[0][3] += a; 02205 x[1][0] += a; 02206 x[1][1] += a; 02207 x[1][2] += a; 02208 x[1][3] += a; 02209 x[2][0] += a; 02210 x[2][1] += a; 02211 x[2][2] += a; 02212 x[2][3] += a; 02213 x[3][0] += a; 02214 x[3][1] += a; 02215 x[3][2] += a; 02216 x[3][3] += a; 02217 02218 return *this; 02219 } 02220 02221 template <class T> 02222 Matrix44<T> 02223 Matrix44<T>::operator + (const Matrix44<T> &v) const 02224 { 02225 return Matrix44 (x[0][0] + v.x[0][0], 02226 x[0][1] + v.x[0][1], 02227 x[0][2] + v.x[0][2], 02228 x[0][3] + v.x[0][3], 02229 x[1][0] + v.x[1][0], 02230 x[1][1] + v.x[1][1], 02231 x[1][2] + v.x[1][2], 02232 x[1][3] + v.x[1][3], 02233 x[2][0] + v.x[2][0], 02234 x[2][1] + v.x[2][1], 02235 x[2][2] + v.x[2][2], 02236 x[2][3] + v.x[2][3], 02237 x[3][0] + v.x[3][0], 02238 x[3][1] + v.x[3][1], 02239 x[3][2] + v.x[3][2], 02240 x[3][3] + v.x[3][3]); 02241 } 02242 02243 template <class T> 02244 const Matrix44<T> & 02245 Matrix44<T>::operator -= (const Matrix44<T> &v) 02246 { 02247 x[0][0] -= v.x[0][0]; 02248 x[0][1] -= v.x[0][1]; 02249 x[0][2] -= v.x[0][2]; 02250 x[0][3] -= v.x[0][3]; 02251 x[1][0] -= v.x[1][0]; 02252 x[1][1] -= v.x[1][1]; 02253 x[1][2] -= v.x[1][2]; 02254 x[1][3] -= v.x[1][3]; 02255 x[2][0] -= v.x[2][0]; 02256 x[2][1] -= v.x[2][1]; 02257 x[2][2] -= v.x[2][2]; 02258 x[2][3] -= v.x[2][3]; 02259 x[3][0] -= v.x[3][0]; 02260 x[3][1] -= v.x[3][1]; 02261 x[3][2] -= v.x[3][2]; 02262 x[3][3] -= v.x[3][3]; 02263 02264 return *this; 02265 } 02266 02267 template <class T> 02268 const Matrix44<T> & 02269 Matrix44<T>::operator -= (T a) 02270 { 02271 x[0][0] -= a; 02272 x[0][1] -= a; 02273 x[0][2] -= a; 02274 x[0][3] -= a; 02275 x[1][0] -= a; 02276 x[1][1] -= a; 02277 x[1][2] -= a; 02278 x[1][3] -= a; 02279 x[2][0] -= a; 02280 x[2][1] -= a; 02281 x[2][2] -= a; 02282 x[2][3] -= a; 02283 x[3][0] -= a; 02284 x[3][1] -= a; 02285 x[3][2] -= a; 02286 x[3][3] -= a; 02287 02288 return *this; 02289 } 02290 02291 template <class T> 02292 Matrix44<T> 02293 Matrix44<T>::operator - (const Matrix44<T> &v) const 02294 { 02295 return Matrix44 (x[0][0] - v.x[0][0], 02296 x[0][1] - v.x[0][1], 02297 x[0][2] - v.x[0][2], 02298 x[0][3] - v.x[0][3], 02299 x[1][0] - v.x[1][0], 02300 x[1][1] - v.x[1][1], 02301 x[1][2] - v.x[1][2], 02302 x[1][3] - v.x[1][3], 02303 x[2][0] - v.x[2][0], 02304 x[2][1] - v.x[2][1], 02305 x[2][2] - v.x[2][2], 02306 x[2][3] - v.x[2][3], 02307 x[3][0] - v.x[3][0], 02308 x[3][1] - v.x[3][1], 02309 x[3][2] - v.x[3][2], 02310 x[3][3] - v.x[3][3]); 02311 } 02312 02313 template <class T> 02314 Matrix44<T> 02315 Matrix44<T>::operator - () const 02316 { 02317 return Matrix44 (-x[0][0], 02318 -x[0][1], 02319 -x[0][2], 02320 -x[0][3], 02321 -x[1][0], 02322 -x[1][1], 02323 -x[1][2], 02324 -x[1][3], 02325 -x[2][0], 02326 -x[2][1], 02327 -x[2][2], 02328 -x[2][3], 02329 -x[3][0], 02330 -x[3][1], 02331 -x[3][2], 02332 -x[3][3]); 02333 } 02334 02335 template <class T> 02336 const Matrix44<T> & 02337 Matrix44<T>::negate () 02338 { 02339 x[0][0] = -x[0][0]; 02340 x[0][1] = -x[0][1]; 02341 x[0][2] = -x[0][2]; 02342 x[0][3] = -x[0][3]; 02343 x[1][0] = -x[1][0]; 02344 x[1][1] = -x[1][1]; 02345 x[1][2] = -x[1][2]; 02346 x[1][3] = -x[1][3]; 02347 x[2][0] = -x[2][0]; 02348 x[2][1] = -x[2][1]; 02349 x[2][2] = -x[2][2]; 02350 x[2][3] = -x[2][3]; 02351 x[3][0] = -x[3][0]; 02352 x[3][1] = -x[3][1]; 02353 x[3][2] = -x[3][2]; 02354 x[3][3] = -x[3][3]; 02355 02356 return *this; 02357 } 02358 02359 template <class T> 02360 const Matrix44<T> & 02361 Matrix44<T>::operator *= (T a) 02362 { 02363 x[0][0] *= a; 02364 x[0][1] *= a; 02365 x[0][2] *= a; 02366 x[0][3] *= a; 02367 x[1][0] *= a; 02368 x[1][1] *= a; 02369 x[1][2] *= a; 02370 x[1][3] *= a; 02371 x[2][0] *= a; 02372 x[2][1] *= a; 02373 x[2][2] *= a; 02374 x[2][3] *= a; 02375 x[3][0] *= a; 02376 x[3][1] *= a; 02377 x[3][2] *= a; 02378 x[3][3] *= a; 02379 02380 return *this; 02381 } 02382 02383 template <class T> 02384 Matrix44<T> 02385 Matrix44<T>::operator * (T a) const 02386 { 02387 return Matrix44 (x[0][0] * a, 02388 x[0][1] * a, 02389 x[0][2] * a, 02390 x[0][3] * a, 02391 x[1][0] * a, 02392 x[1][1] * a, 02393 x[1][2] * a, 02394 x[1][3] * a, 02395 x[2][0] * a, 02396 x[2][1] * a, 02397 x[2][2] * a, 02398 x[2][3] * a, 02399 x[3][0] * a, 02400 x[3][1] * a, 02401 x[3][2] * a, 02402 x[3][3] * a); 02403 } 02404 02405 template <class T> 02406 inline Matrix44<T> 02407 operator * (T a, const Matrix44<T> &v) 02408 { 02409 return v * a; 02410 } 02411 02412 template <class T> 02413 inline const Matrix44<T> & 02414 Matrix44<T>::operator *= (const Matrix44<T> &v) 02415 { 02416 Matrix44 tmp (T (0)); 02417 02418 multiply (*this, v, tmp); 02419 *this = tmp; 02420 return *this; 02421 } 02422 02423 template <class T> 02424 inline Matrix44<T> 02425 Matrix44<T>::operator * (const Matrix44<T> &v) const 02426 { 02427 Matrix44 tmp (T (0)); 02428 02429 multiply (*this, v, tmp); 02430 return tmp; 02431 } 02432 02433 template <class T> 02434 void 02435 Matrix44<T>::multiply (const Matrix44<T> &a, 02436 const Matrix44<T> &b, 02437 Matrix44<T> &c) 02438 { 02439 register const T * IMATH_RESTRICT ap = &a.x[0][0]; 02440 register const T * IMATH_RESTRICT bp = &b.x[0][0]; 02441 register T * IMATH_RESTRICT cp = &c.x[0][0]; 02442 02443 register T a0, a1, a2, a3; 02444 02445 a0 = ap[0]; 02446 a1 = ap[1]; 02447 a2 = ap[2]; 02448 a3 = ap[3]; 02449 02450 cp[0] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12]; 02451 cp[1] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13]; 02452 cp[2] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14]; 02453 cp[3] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15]; 02454 02455 a0 = ap[4]; 02456 a1 = ap[5]; 02457 a2 = ap[6]; 02458 a3 = ap[7]; 02459 02460 cp[4] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12]; 02461 cp[5] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13]; 02462 cp[6] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14]; 02463 cp[7] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15]; 02464 02465 a0 = ap[8]; 02466 a1 = ap[9]; 02467 a2 = ap[10]; 02468 a3 = ap[11]; 02469 02470 cp[8] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12]; 02471 cp[9] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13]; 02472 cp[10] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14]; 02473 cp[11] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15]; 02474 02475 a0 = ap[12]; 02476 a1 = ap[13]; 02477 a2 = ap[14]; 02478 a3 = ap[15]; 02479 02480 cp[12] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12]; 02481 cp[13] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13]; 02482 cp[14] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14]; 02483 cp[15] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15]; 02484 } 02485 02486 template <class T> template <class S> 02487 void 02488 Matrix44<T>::multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const 02489 { 02490 S a, b, c, w; 02491 02492 a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0] + x[3][0]; 02493 b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1] + x[3][1]; 02494 c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2] + x[3][2]; 02495 w = src[0] * x[0][3] + src[1] * x[1][3] + src[2] * x[2][3] + x[3][3]; 02496 02497 dst.x = a / w; 02498 dst.y = b / w; 02499 dst.z = c / w; 02500 } 02501 02502 template <class T> template <class S> 02503 void 02504 Matrix44<T>::multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const 02505 { 02506 S a, b, c; 02507 02508 a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0]; 02509 b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1]; 02510 c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2]; 02511 02512 dst.x = a; 02513 dst.y = b; 02514 dst.z = c; 02515 } 02516 02517 template <class T> 02518 const Matrix44<T> & 02519 Matrix44<T>::operator /= (T a) 02520 { 02521 x[0][0] /= a; 02522 x[0][1] /= a; 02523 x[0][2] /= a; 02524 x[0][3] /= a; 02525 x[1][0] /= a; 02526 x[1][1] /= a; 02527 x[1][2] /= a; 02528 x[1][3] /= a; 02529 x[2][0] /= a; 02530 x[2][1] /= a; 02531 x[2][2] /= a; 02532 x[2][3] /= a; 02533 x[3][0] /= a; 02534 x[3][1] /= a; 02535 x[3][2] /= a; 02536 x[3][3] /= a; 02537 02538 return *this; 02539 } 02540 02541 template <class T> 02542 Matrix44<T> 02543 Matrix44<T>::operator / (T a) const 02544 { 02545 return Matrix44 (x[0][0] / a, 02546 x[0][1] / a, 02547 x[0][2] / a, 02548 x[0][3] / a, 02549 x[1][0] / a, 02550 x[1][1] / a, 02551 x[1][2] / a, 02552 x[1][3] / a, 02553 x[2][0] / a, 02554 x[2][1] / a, 02555 x[2][2] / a, 02556 x[2][3] / a, 02557 x[3][0] / a, 02558 x[3][1] / a, 02559 x[3][2] / a, 02560 x[3][3] / a); 02561 } 02562 02563 template <class T> 02564 const Matrix44<T> & 02565 Matrix44<T>::transpose () 02566 { 02567 Matrix44 tmp (x[0][0], 02568 x[1][0], 02569 x[2][0], 02570 x[3][0], 02571 x[0][1], 02572 x[1][1], 02573 x[2][1], 02574 x[3][1], 02575 x[0][2], 02576 x[1][2], 02577 x[2][2], 02578 x[3][2], 02579 x[0][3], 02580 x[1][3], 02581 x[2][3], 02582 x[3][3]); 02583 *this = tmp; 02584 return *this; 02585 } 02586 02587 template <class T> 02588 Matrix44<T> 02589 Matrix44<T>::transposed () const 02590 { 02591 return Matrix44 (x[0][0], 02592 x[1][0], 02593 x[2][0], 02594 x[3][0], 02595 x[0][1], 02596 x[1][1], 02597 x[2][1], 02598 x[3][1], 02599 x[0][2], 02600 x[1][2], 02601 x[2][2], 02602 x[3][2], 02603 x[0][3], 02604 x[1][3], 02605 x[2][3], 02606 x[3][3]); 02607 } 02608 02609 template <class T> 02610 const Matrix44<T> & 02611 Matrix44<T>::gjInvert (bool singExc) throw (Iex::MathExc) 02612 { 02613 *this = gjInverse (singExc); 02614 return *this; 02615 } 02616 02617 template <class T> 02618 Matrix44<T> 02619 Matrix44<T>::gjInverse (bool singExc) const throw (Iex::MathExc) 02620 { 02621 int i, j, k; 02622 Matrix44 s; 02623 Matrix44 t (*this); 02624 02625 // Forward elimination 02626 02627 for (i = 0; i < 3 ; i++) 02628 { 02629 int pivot = i; 02630 02631 T pivotsize = t[i][i]; 02632 02633 if (pivotsize < 0) 02634 pivotsize = -pivotsize; 02635 02636 for (j = i + 1; j < 4; j++) 02637 { 02638 T tmp = t[j][i]; 02639 02640 if (tmp < 0) 02641 tmp = -tmp; 02642 02643 if (tmp > pivotsize) 02644 { 02645 pivot = j; 02646 pivotsize = tmp; 02647 } 02648 } 02649 02650 if (pivotsize == 0) 02651 { 02652 if (singExc) 02653 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix."); 02654 02655 return Matrix44(); 02656 } 02657 02658 if (pivot != i) 02659 { 02660 for (j = 0; j < 4; j++) 02661 { 02662 T tmp; 02663 02664 tmp = t[i][j]; 02665 t[i][j] = t[pivot][j]; 02666 t[pivot][j] = tmp; 02667 02668 tmp = s[i][j]; 02669 s[i][j] = s[pivot][j]; 02670 s[pivot][j] = tmp; 02671 } 02672 } 02673 02674 for (j = i + 1; j < 4; j++) 02675 { 02676 T f = t[j][i] / t[i][i]; 02677 02678 for (k = 0; k < 4; k++) 02679 { 02680 t[j][k] -= f * t[i][k]; 02681 s[j][k] -= f * s[i][k]; 02682 } 02683 } 02684 } 02685 02686 // Backward substitution 02687 02688 for (i = 3; i >= 0; --i) 02689 { 02690 T f; 02691 02692 if ((f = t[i][i]) == 0) 02693 { 02694 if (singExc) 02695 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix."); 02696 02697 return Matrix44(); 02698 } 02699 02700 for (j = 0; j < 4; j++) 02701 { 02702 t[i][j] /= f; 02703 s[i][j] /= f; 02704 } 02705 02706 for (j = 0; j < i; j++) 02707 { 02708 f = t[j][i]; 02709 02710 for (k = 0; k < 4; k++) 02711 { 02712 t[j][k] -= f * t[i][k]; 02713 s[j][k] -= f * s[i][k]; 02714 } 02715 } 02716 } 02717 02718 return s; 02719 } 02720 02721 template <class T> 02722 const Matrix44<T> & 02723 Matrix44<T>::invert (bool singExc) throw (Iex::MathExc) 02724 { 02725 *this = inverse (singExc); 02726 return *this; 02727 } 02728 02729 template <class T> 02730 Matrix44<T> 02731 Matrix44<T>::inverse (bool singExc) const throw (Iex::MathExc) 02732 { 02733 if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1) 02734 return gjInverse(singExc); 02735 02736 Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2], 02737 x[2][1] * x[0][2] - x[0][1] * x[2][2], 02738 x[0][1] * x[1][2] - x[1][1] * x[0][2], 02739 0, 02740 02741 x[2][0] * x[1][2] - x[1][0] * x[2][2], 02742 x[0][0] * x[2][2] - x[2][0] * x[0][2], 02743 x[1][0] * x[0][2] - x[0][0] * x[1][2], 02744 0, 02745 02746 x[1][0] * x[2][1] - x[2][0] * x[1][1], 02747 x[2][0] * x[0][1] - x[0][0] * x[2][1], 02748 x[0][0] * x[1][1] - x[1][0] * x[0][1], 02749 0, 02750 02751 0, 02752 0, 02753 0, 02754 1); 02755 02756 T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0]; 02757 02758 if (Imath::abs (r) >= 1) 02759 { 02760 for (int i = 0; i < 3; ++i) 02761 { 02762 for (int j = 0; j < 3; ++j) 02763 { 02764 s[i][j] /= r; 02765 } 02766 } 02767 } 02768 else 02769 { 02770 T mr = Imath::abs (r) / limits<T>::smallest(); 02771 02772 for (int i = 0; i < 3; ++i) 02773 { 02774 for (int j = 0; j < 3; ++j) 02775 { 02776 if (mr > Imath::abs (s[i][j])) 02777 { 02778 s[i][j] /= r; 02779 } 02780 else 02781 { 02782 if (singExc) 02783 throw SingMatrixExc ("Cannot invert singular matrix."); 02784 02785 return Matrix44(); 02786 } 02787 } 02788 } 02789 } 02790 02791 s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0]; 02792 s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1]; 02793 s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2]; 02794 02795 return s; 02796 } 02797 02798 template <class T> 02799 template <class S> 02800 const Matrix44<T> & 02801 Matrix44<T>::setEulerAngles (const Vec3<S>& r) 02802 { 02803 S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx; 02804 02805 cos_rz = Math<T>::cos (r[2]); 02806 cos_ry = Math<T>::cos (r[1]); 02807 cos_rx = Math<T>::cos (r[0]); 02808 02809 sin_rz = Math<T>::sin (r[2]); 02810 sin_ry = Math<T>::sin (r[1]); 02811 sin_rx = Math<T>::sin (r[0]); 02812 02813 x[0][0] = cos_rz * cos_ry; 02814 x[0][1] = sin_rz * cos_ry; 02815 x[0][2] = -sin_ry; 02816 x[0][3] = 0; 02817 02818 x[1][0] = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx; 02819 x[1][1] = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx; 02820 x[1][2] = cos_ry * sin_rx; 02821 x[1][3] = 0; 02822 02823 x[2][0] = sin_rz * sin_rx + cos_rz * sin_ry * cos_rx; 02824 x[2][1] = -cos_rz * sin_rx + sin_rz * sin_ry * cos_rx; 02825 x[2][2] = cos_ry * cos_rx; 02826 x[2][3] = 0; 02827 02828 x[3][0] = 0; 02829 x[3][1] = 0; 02830 x[3][2] = 0; 02831 x[3][3] = 1; 02832 02833 return *this; 02834 } 02835 02836 template <class T> 02837 template <class S> 02838 const Matrix44<T> & 02839 Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle) 02840 { 02841 Vec3<S> unit (axis.normalized()); 02842 S sine = Math<T>::sin (angle); 02843 S cosine = Math<T>::cos (angle); 02844 02845 x[0][0] = unit[0] * unit[0] * (1 - cosine) + cosine; 02846 x[0][1] = unit[0] * unit[1] * (1 - cosine) + unit[2] * sine; 02847 x[0][2] = unit[0] * unit[2] * (1 - cosine) - unit[1] * sine; 02848 x[0][3] = 0; 02849 02850 x[1][0] = unit[0] * unit[1] * (1 - cosine) - unit[2] * sine; 02851 x[1][1] = unit[1] * unit[1] * (1 - cosine) + cosine; 02852 x[1][2] = unit[1] * unit[2] * (1 - cosine) + unit[0] * sine; 02853 x[1][3] = 0; 02854 02855 x[2][0] = unit[0] * unit[2] * (1 - cosine) + unit[1] * sine; 02856 x[2][1] = unit[1] * unit[2] * (1 - cosine) - unit[0] * sine; 02857 x[2][2] = unit[2] * unit[2] * (1 - cosine) + cosine; 02858 x[2][3] = 0; 02859 02860 x[3][0] = 0; 02861 x[3][1] = 0; 02862 x[3][2] = 0; 02863 x[3][3] = 1; 02864 02865 return *this; 02866 } 02867 02868 template <class T> 02869 template <class S> 02870 const Matrix44<T> & 02871 Matrix44<T>::rotate (const Vec3<S> &r) 02872 { 02873 S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx; 02874 S m00, m01, m02; 02875 S m10, m11, m12; 02876 S m20, m21, m22; 02877 02878 cos_rz = Math<S>::cos (r[2]); 02879 cos_ry = Math<S>::cos (r[1]); 02880 cos_rx = Math<S>::cos (r[0]); 02881 02882 sin_rz = Math<S>::sin (r[2]); 02883 sin_ry = Math<S>::sin (r[1]); 02884 sin_rx = Math<S>::sin (r[0]); 02885 02886 m00 = cos_rz * cos_ry; 02887 m01 = sin_rz * cos_ry; 02888 m02 = -sin_ry; 02889 m10 = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx; 02890 m11 = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx; 02891 m12 = cos_ry * sin_rx; 02892 m20 = -sin_rz * -sin_rx + cos_rz * sin_ry * cos_rx; 02893 m21 = cos_rz * -sin_rx + sin_rz * sin_ry * cos_rx; 02894 m22 = cos_ry * cos_rx; 02895 02896 Matrix44<T> P (*this); 02897 02898 x[0][0] = P[0][0] * m00 + P[1][0] * m01 + P[2][0] * m02; 02899 x[0][1] = P[0][1] * m00 + P[1][1] * m01 + P[2][1] * m02; 02900 x[0][2] = P[0][2] * m00 + P[1][2] * m01 + P[2][2] * m02; 02901 x[0][3] = P[0][3] * m00 + P[1][3] * m01 + P[2][3] * m02; 02902 02903 x[1][0] = P[0][0] * m10 + P[1][0] * m11 + P[2][0] * m12; 02904 x[1][1] = P[0][1] * m10 + P[1][1] * m11 + P[2][1] * m12; 02905 x[1][2] = P[0][2] * m10 + P[1][2] * m11 + P[2][2] * m12; 02906 x[1][3] = P[0][3] * m10 + P[1][3] * m11 + P[2][3] * m12; 02907 02908 x[2][0] = P[0][0] * m20 + P[1][0] * m21 + P[2][0] * m22; 02909 x[2][1] = P[0][1] * m20 + P[1][1] * m21 + P[2][1] * m22; 02910 x[2][2] = P[0][2] * m20 + P[1][2] * m21 + P[2][2] * m22; 02911 x[2][3] = P[0][3] * m20 + P[1][3] * m21 + P[2][3] * m22; 02912 02913 return *this; 02914 } 02915 02916 template <class T> 02917 const Matrix44<T> & 02918 Matrix44<T>::setScale (T s) 02919 { 02920 memset (x, 0, sizeof (x)); 02921 x[0][0] = s; 02922 x[1][1] = s; 02923 x[2][2] = s; 02924 x[3][3] = 1; 02925 02926 return *this; 02927 } 02928 02929 template <class T> 02930 template <class S> 02931 const Matrix44<T> & 02932 Matrix44<T>::setScale (const Vec3<S> &s) 02933 { 02934 memset (x, 0, sizeof (x)); 02935 x[0][0] = s[0]; 02936 x[1][1] = s[1]; 02937 x[2][2] = s[2]; 02938 x[3][3] = 1; 02939 02940 return *this; 02941 } 02942 02943 template <class T> 02944 template <class S> 02945 const Matrix44<T> & 02946 Matrix44<T>::scale (const Vec3<S> &s) 02947 { 02948 x[0][0] *= s[0]; 02949 x[0][1] *= s[0]; 02950 x[0][2] *= s[0]; 02951 x[0][3] *= s[0]; 02952 02953 x[1][0] *= s[1]; 02954 x[1][1] *= s[1]; 02955 x[1][2] *= s[1]; 02956 x[1][3] *= s[1]; 02957 02958 x[2][0] *= s[2]; 02959 x[2][1] *= s[2]; 02960 x[2][2] *= s[2]; 02961 x[2][3] *= s[2]; 02962 02963 return *this; 02964 } 02965 02966 template <class T> 02967 template <class S> 02968 const Matrix44<T> & 02969 Matrix44<T>::setTranslation (const Vec3<S> &t) 02970 { 02971 x[0][0] = 1; 02972 x[0][1] = 0; 02973 x[0][2] = 0; 02974 x[0][3] = 0; 02975 02976 x[1][0] = 0; 02977 x[1][1] = 1; 02978 x[1][2] = 0; 02979 x[1][3] = 0; 02980 02981 x[2][0] = 0; 02982 x[2][1] = 0; 02983 x[2][2] = 1; 02984 x[2][3] = 0; 02985 02986 x[3][0] = t[0]; 02987 x[3][1] = t[1]; 02988 x[3][2] = t[2]; 02989 x[3][3] = 1; 02990 02991 return *this; 02992 } 02993 02994 template <class T> 02995 inline const Vec3<T> 02996 Matrix44<T>::translation () const 02997 { 02998 return Vec3<T> (x[3][0], x[3][1], x[3][2]); 02999 } 03000 03001 template <class T> 03002 template <class S> 03003 const Matrix44<T> & 03004 Matrix44<T>::translate (const Vec3<S> &t) 03005 { 03006 x[3][0] += t[0] * x[0][0] + t[1] * x[1][0] + t[2] * x[2][0]; 03007 x[3][1] += t[0] * x[0][1] + t[1] * x[1][1] + t[2] * x[2][1]; 03008 x[3][2] += t[0] * x[0][2] + t[1] * x[1][2] + t[2] * x[2][2]; 03009 x[3][3] += t[0] * x[0][3] + t[1] * x[1][3] + t[2] * x[2][3]; 03010 03011 return *this; 03012 } 03013 03014 template <class T> 03015 template <class S> 03016 const Matrix44<T> & 03017 Matrix44<T>::setShear (const Vec3<S> &h) 03018 { 03019 x[0][0] = 1; 03020 x[0][1] = 0; 03021 x[0][2] = 0; 03022 x[0][3] = 0; 03023 03024 x[1][0] = h[0]; 03025 x[1][1] = 1; 03026 x[1][2] = 0; 03027 x[1][3] = 0; 03028 03029 x[2][0] = h[1]; 03030 x[2][1] = h[2]; 03031 x[2][2] = 1; 03032 x[2][3] = 0; 03033 03034 x[3][0] = 0; 03035 x[3][1] = 0; 03036 x[3][2] = 0; 03037 x[3][3] = 1; 03038 03039 return *this; 03040 } 03041 03042 template <class T> 03043 template <class S> 03044 const Matrix44<T> & 03045 Matrix44<T>::setShear (const Shear6<S> &h) 03046 { 03047 x[0][0] = 1; 03048 x[0][1] = h.yx; 03049 x[0][2] = h.zx; 03050 x[0][3] = 0; 03051 03052 x[1][0] = h.xy; 03053 x[1][1] = 1; 03054 x[1][2] = h.zy; 03055 x[1][3] = 0; 03056 03057 x[2][0] = h.xz; 03058 x[2][1] = h.yz; 03059 x[2][2] = 1; 03060 x[2][3] = 0; 03061 03062 x[3][0] = 0; 03063 x[3][1] = 0; 03064 x[3][2] = 0; 03065 x[3][3] = 1; 03066 03067 return *this; 03068 } 03069 03070 template <class T> 03071 template <class S> 03072 const Matrix44<T> & 03073 Matrix44<T>::shear (const Vec3<S> &h) 03074 { 03075 // 03076 // In this case, we don't need a temp. copy of the matrix 03077 // because we never use a value on the RHS after we've 03078 // changed it on the LHS. 03079 // 03080 03081 for (int i=0; i < 4; i++) 03082 { 03083 x[2][i] += h[1] * x[0][i] + h[2] * x[1][i]; 03084 x[1][i] += h[0] * x[0][i]; 03085 } 03086 03087 return *this; 03088 } 03089 03090 template <class T> 03091 template <class S> 03092 const Matrix44<T> & 03093 Matrix44<T>::shear (const Shear6<S> &h) 03094 { 03095 Matrix44<T> P (*this); 03096 03097 for (int i=0; i < 4; i++) 03098 { 03099 x[0][i] = P[0][i] + h.yx * P[1][i] + h.zx * P[2][i]; 03100 x[1][i] = h.xy * P[0][i] + P[1][i] + h.zy * P[2][i]; 03101 x[2][i] = h.xz * P[0][i] + h.yz * P[1][i] + P[2][i]; 03102 } 03103 03104 return *this; 03105 } 03106 03107 03108 //-------------------------------- 03109 // Implementation of stream output 03110 //-------------------------------- 03111 03112 template <class T> 03113 std::ostream & 03114 operator << (std::ostream &s, const Matrix33<T> &m) 03115 { 03116 std::ios_base::fmtflags oldFlags = s.flags(); 03117 int width; 03118 03119 if (s.flags() & std::ios_base::fixed) 03120 { 03121 s.setf (std::ios_base::showpoint); 03122 width = s.precision() + 5; 03123 } 03124 else 03125 { 03126 s.setf (std::ios_base::scientific); 03127 s.setf (std::ios_base::showpoint); 03128 width = s.precision() + 8; 03129 } 03130 03131 s << "(" << std::setw (width) << m[0][0] << 03132 " " << std::setw (width) << m[0][1] << 03133 " " << std::setw (width) << m[0][2] << "\n" << 03134 03135 " " << std::setw (width) << m[1][0] << 03136 " " << std::setw (width) << m[1][1] << 03137 " " << std::setw (width) << m[1][2] << "\n" << 03138 03139 " " << std::setw (width) << m[2][0] << 03140 " " << std::setw (width) << m[2][1] << 03141 " " << std::setw (width) << m[2][2] << ")\n"; 03142 03143 s.flags (oldFlags); 03144 return s; 03145 } 03146 03147 template <class T> 03148 std::ostream & 03149 operator << (std::ostream &s, const Matrix44<T> &m) 03150 { 03151 std::ios_base::fmtflags oldFlags = s.flags(); 03152 int width; 03153 03154 if (s.flags() & std::ios_base::fixed) 03155 { 03156 s.setf (std::ios_base::showpoint); 03157 width = s.precision() + 5; 03158 } 03159 else 03160 { 03161 s.setf (std::ios_base::scientific); 03162 s.setf (std::ios_base::showpoint); 03163 width = s.precision() + 8; 03164 } 03165 03166 s << "(" << std::setw (width) << m[0][0] << 03167 " " << std::setw (width) << m[0][1] << 03168 " " << std::setw (width) << m[0][2] << 03169 " " << std::setw (width) << m[0][3] << "\n" << 03170 03171 " " << std::setw (width) << m[1][0] << 03172 " " << std::setw (width) << m[1][1] << 03173 " " << std::setw (width) << m[1][2] << 03174 " " << std::setw (width) << m[1][3] << "\n" << 03175 03176 " " << std::setw (width) << m[2][0] << 03177 " " << std::setw (width) << m[2][1] << 03178 " " << std::setw (width) << m[2][2] << 03179 " " << std::setw (width) << m[2][3] << "\n" << 03180 03181 " " << std::setw (width) << m[3][0] << 03182 " " << std::setw (width) << m[3][1] << 03183 " " << std::setw (width) << m[3][2] << 03184 " " << std::setw (width) << m[3][3] << ")\n"; 03185 03186 s.flags (oldFlags); 03187 return s; 03188 } 03189 03190 03191 //--------------------------------------------------------------- 03192 // Implementation of vector-times-matrix multiplication operators 03193 //--------------------------------------------------------------- 03194 03195 template <class S, class T> 03196 inline const Vec2<S> & 03197 operator *= (Vec2<S> &v, const Matrix33<T> &m) 03198 { 03199 S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]); 03200 S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]); 03201 S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]); 03202 03203 v.x = x / w; 03204 v.y = y / w; 03205 03206 return v; 03207 } 03208 03209 template <class S, class T> 03210 inline Vec2<S> 03211 operator * (const Vec2<S> &v, const Matrix33<T> &m) 03212 { 03213 S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]); 03214 S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]); 03215 S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]); 03216 03217 return Vec2<S> (x / w, y / w); 03218 } 03219 03220 03221 template <class S, class T> 03222 inline const Vec3<S> & 03223 operator *= (Vec3<S> &v, const Matrix33<T> &m) 03224 { 03225 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]); 03226 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]); 03227 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]); 03228 03229 v.x = x; 03230 v.y = y; 03231 v.z = z; 03232 03233 return v; 03234 } 03235 03236 template <class S, class T> 03237 inline Vec3<S> 03238 operator * (const Vec3<S> &v, const Matrix33<T> &m) 03239 { 03240 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]); 03241 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]); 03242 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]); 03243 03244 return Vec3<S> (x, y, z); 03245 } 03246 03247 03248 template <class S, class T> 03249 inline const Vec3<S> & 03250 operator *= (Vec3<S> &v, const Matrix44<T> &m) 03251 { 03252 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]); 03253 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]); 03254 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]); 03255 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]); 03256 03257 v.x = x / w; 03258 v.y = y / w; 03259 v.z = z / w; 03260 03261 return v; 03262 } 03263 03264 template <class S, class T> 03265 inline Vec3<S> 03266 operator * (const Vec3<S> &v, const Matrix44<T> &m) 03267 { 03268 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]); 03269 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]); 03270 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]); 03271 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]); 03272 03273 return Vec3<S> (x / w, y / w, z / w); 03274 } 03275 03276 03277 template <class S, class T> 03278 inline const Vec4<S> & 03279 operator *= (Vec4<S> &v, const Matrix44<T> &m) 03280 { 03281 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + v.w * m[3][0]); 03282 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + v.w * m[3][1]); 03283 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + v.w * m[3][2]); 03284 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + v.w * m[3][3]); 03285 03286 v.x = x; 03287 v.y = y; 03288 v.z = z; 03289 v.w = w; 03290 03291 return v; 03292 } 03293 03294 template <class S, class T> 03295 inline Vec4<S> 03296 operator * (const Vec4<S> &v, const Matrix44<T> &m) 03297 { 03298 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + v.w * m[3][0]); 03299 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + v.w * m[3][1]); 03300 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + v.w * m[3][2]); 03301 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + v.w * m[3][3]); 03302 03303 return Vec4<S> (x, y, z, w); 03304 } 03305 03306 } // namespace Imath 03307 03308 03309 03310 #endif