1 /////////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas 4 // Digital Ltd. LLC 5 // 6 // All rights reserved. 7 // 8 // Redistribution and use in source and binary forms, with or without 9 // modification, are permitted provided that the following conditions are 10 // met: 11 // * Redistributions of source code must retain the above copyright 12 // notice, this list of conditions and the following disclaimer. 13 // * Redistributions in binary form must reproduce the above 14 // copyright notice, this list of conditions and the following disclaimer 15 // in the documentation and/or other materials provided with the 16 // distribution. 17 // * Neither the name of Industrial Light & Magic nor the names of 18 // its contributors may be used to endorse or promote products derived 19 // from this software without specific prior written permission. 20 // 21 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 22 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 23 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 24 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 25 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 26 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 27 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 28 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 29 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 30 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 31 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 32 // 33 /////////////////////////////////////////////////////////////////////////// 34 35 36 37 #ifndef INCLUDED_IMATHMATRIX_H 38 #define INCLUDED_IMATHMATRIX_H 39 40 //---------------------------------------------------------------- 41 // 42 // 2D (3x3) and 3D (4x4) transformation matrix templates. 43 // 44 //---------------------------------------------------------------- 45 46 #include "ImathPlatform.h" 47 #include "ImathFun.h" 48 #include "ImathExc.h" 49 #include "ImathVec.h" 50 #include "ImathShear.h" 51 52 #include <cstring> 53 #include <iostream> 54 #include <iomanip> 55 #include <string.h> 56 57 #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER 58 // suppress exception specification warnings 59 #pragma warning(disable:4290) 60 #endif 61 62 63 namespace Imath { 64 65 enum Uninitialized {UNINITIALIZED}; 66 67 68 template <class T> class Matrix33 69 { 70 public: 71 72 //------------------- 73 // Access to elements 74 //------------------- 75 76 T x[3][3]; 77 78 T * operator [] (int i); 79 const T * operator [] (int i) const; 80 81 82 //------------- 83 // Constructors 84 //------------- 85 86 Matrix33 (Uninitialized) {} 87 88 Matrix33 (); 89 // 1 0 0 90 // 0 1 0 91 // 0 0 1 92 93 Matrix33 (T a); 94 // a a a 95 // a a a 96 // a a a 97 98 Matrix33 (const T a[3][3]); 99 // a[0][0] a[0][1] a[0][2] 100 // a[1][0] a[1][1] a[1][2] 101 // a[2][0] a[2][1] a[2][2] 102 103 Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i); 104 105 // a b c 106 // d e f 107 // g h i 108 109 110 //-------------------------------- 111 // Copy constructor and assignment 112 //-------------------------------- 113 114 Matrix33 (const Matrix33 &v); 115 template <class S> explicit Matrix33 (const Matrix33<S> &v); 116 117 const Matrix33 & operator = (const Matrix33 &v); 118 const Matrix33 & operator = (T a); 119 120 121 //---------------------- 122 // Compatibility with Sb 123 //---------------------- 124 125 T * getValue (); 126 const T * getValue () const; 127 128 template <class S> 129 void getValue (Matrix33<S> &v) const; 130 template <class S> 131 Matrix33 & setValue (const Matrix33<S> &v); 132 133 template <class S> 134 Matrix33 & setTheMatrix (const Matrix33<S> &v); 135 136 137 //--------- 138 // Identity 139 //--------- 140 141 void makeIdentity(); 142 143 144 //--------- 145 // Equality 146 //--------- 147 148 bool operator == (const Matrix33 &v) const; 149 bool operator != (const Matrix33 &v) const; 150 151 //----------------------------------------------------------------------- 152 // Compare two matrices and test if they are "approximately equal": 153 // 154 // equalWithAbsError (m, e) 155 // 156 // Returns true if the coefficients of this and m are the same with 157 // an absolute error of no more than e, i.e., for all i, j 158 // 159 // abs (this[i][j] - m[i][j]) <= e 160 // 161 // equalWithRelError (m, e) 162 // 163 // Returns true if the coefficients of this and m are the same with 164 // a relative error of no more than e, i.e., for all i, j 165 // 166 // abs (this[i] - v[i][j]) <= e * abs (this[i][j]) 167 //----------------------------------------------------------------------- 168 169 bool equalWithAbsError (const Matrix33<T> &v, T e) const; 170 bool equalWithRelError (const Matrix33<T> &v, T e) const; 171 172 173 //------------------------ 174 // Component-wise addition 175 //------------------------ 176 177 const Matrix33 & operator += (const Matrix33 &v); 178 const Matrix33 & operator += (T a); 179 Matrix33 operator + (const Matrix33 &v) const; 180 181 182 //--------------------------- 183 // Component-wise subtraction 184 //--------------------------- 185 186 const Matrix33 & operator -= (const Matrix33 &v); 187 const Matrix33 & operator -= (T a); 188 Matrix33 operator - (const Matrix33 &v) const; 189 190 191 //------------------------------------ 192 // Component-wise multiplication by -1 193 //------------------------------------ 194 195 Matrix33 operator - () const; 196 const Matrix33 & negate (); 197 198 199 //------------------------------ 200 // Component-wise multiplication 201 //------------------------------ 202 203 const Matrix33 & operator *= (T a); 204 Matrix33 operator * (T a) const; 205 206 207 //----------------------------------- 208 // Matrix-times-matrix multiplication 209 //----------------------------------- 210 211 const Matrix33 & operator *= (const Matrix33 &v); 212 Matrix33 operator * (const Matrix33 &v) const; 213 214 215 //----------------------------------------------------------------- 216 // Vector-times-matrix multiplication; see also the "operator *" 217 // functions defined below. 218 // 219 // m.multVecMatrix(src,dst) implements a homogeneous transformation 220 // by computing Vec3 (src.x, src.y, 1) * m and dividing by the 221 // result's third element. 222 // 223 // m.multDirMatrix(src,dst) multiplies src by the upper left 2x2 224 // submatrix, ignoring the rest of matrix m. 225 //----------------------------------------------------------------- 226 227 template <class S> 228 void multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const; 229 230 template <class S> 231 void multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const; 232 233 234 //------------------------ 235 // Component-wise division 236 //------------------------ 237 238 const Matrix33 & operator /= (T a); 239 Matrix33 operator / (T a) const; 240 241 242 //------------------ 243 // Transposed matrix 244 //------------------ 245 246 const Matrix33 & transpose (); 247 Matrix33 transposed () const; 248 249 250 //------------------------------------------------------------ 251 // Inverse matrix: If singExc is false, inverting a singular 252 // matrix produces an identity matrix. If singExc is true, 253 // inverting a singular matrix throws a SingMatrixExc. 254 // 255 // inverse() and invert() invert matrices using determinants; 256 // gjInverse() and gjInvert() use the Gauss-Jordan method. 257 // 258 // inverse() and invert() are significantly faster than 259 // gjInverse() and gjInvert(), but the results may be slightly 260 // less accurate. 261 // 262 //------------------------------------------------------------ 263 264 const Matrix33 & invert (bool singExc = false) 265 throw (Iex::MathExc); 266 267 Matrix33<T> inverse (bool singExc = false) const 268 throw (Iex::MathExc); 269 270 const Matrix33 & gjInvert (bool singExc = false) 271 throw (Iex::MathExc); 272 273 Matrix33<T> gjInverse (bool singExc = false) const 274 throw (Iex::MathExc); 275 276 277 //------------------------------------------------ 278 // Calculate the matrix minor of the (r,c) element 279 //------------------------------------------------ 280 281 T minorOf (const int r, const int c) const; 282 283 //--------------------------------------------------- 284 // Build a minor using the specified rows and columns 285 //--------------------------------------------------- 286 287 T fastMinor (const int r0, const int r1, 288 const int c0, const int c1) const; 289 290 //------------ 291 // Determinant 292 //------------ 293 294 T determinant() const; 295 296 //----------------------------------------- 297 // Set matrix to rotation by r (in radians) 298 //----------------------------------------- 299 300 template <class S> 301 const Matrix33 & setRotation (S r); 302 303 304 //----------------------------- 305 // Rotate the given matrix by r 306 //----------------------------- 307 308 template <class S> 309 const Matrix33 & rotate (S r); 310 311 312 //-------------------------------------------- 313 // Set matrix to scale by given uniform factor 314 //-------------------------------------------- 315 316 const Matrix33 & setScale (T s); 317 318 319 //------------------------------------ 320 // Set matrix to scale by given vector 321 //------------------------------------ 322 323 template <class S> 324 const Matrix33 & setScale (const Vec2<S> &s); 325 326 327 //---------------------- 328 // Scale the matrix by s 329 //---------------------- 330 331 template <class S> 332 const Matrix33 & scale (const Vec2<S> &s); 333 334 335 //------------------------------------------ 336 // Set matrix to translation by given vector 337 //------------------------------------------ 338 339 template <class S> 340 const Matrix33 & setTranslation (const Vec2<S> &t); 341 342 343 //----------------------------- 344 // Return translation component 345 //----------------------------- 346 347 Vec2<T> translation () const; 348 349 350 //-------------------------- 351 // Translate the matrix by t 352 //-------------------------- 353 354 template <class S> 355 const Matrix33 & translate (const Vec2<S> &t); 356 357 358 //----------------------------------------------------------- 359 // Set matrix to shear x for each y coord. by given factor xy 360 //----------------------------------------------------------- 361 362 template <class S> 363 const Matrix33 & setShear (const S &h); 364 365 366 //------------------------------------------------------------- 367 // Set matrix to shear x for each y coord. by given factor h[0] 368 // and to shear y for each x coord. by given factor h[1] 369 //------------------------------------------------------------- 370 371 template <class S> 372 const Matrix33 & setShear (const Vec2<S> &h); 373 374 375 //----------------------------------------------------------- 376 // Shear the matrix in x for each y coord. by given factor xy 377 //----------------------------------------------------------- 378 379 template <class S> 380 const Matrix33 & shear (const S &xy); 381 382 383 //----------------------------------------------------------- 384 // Shear the matrix in x for each y coord. by given factor xy 385 // and shear y for each x coord. by given factor yx 386 //----------------------------------------------------------- 387 388 template <class S> 389 const Matrix33 & shear (const Vec2<S> &h); 390 391 392 //-------------------------------------------------------- 393 // Number of the row and column dimensions, since 394 // Matrix33 is a square matrix. 395 //-------------------------------------------------------- 396 397 static unsigned int dimensions() {return 3;} 398 399 400 //------------------------------------------------- 401 // Limitations of type T (see also class limits<T>) 402 //------------------------------------------------- 403 404 static T baseTypeMin() {return limits<T>::min();} 405 static T baseTypeMax() {return limits<T>::max();} 406 static T baseTypeSmallest() {return limits<T>::smallest();} 407 static T baseTypeEpsilon() {return limits<T>::epsilon();} 408 409 typedef T BaseType; 410 typedef Vec3<T> BaseVecType; 411 412 private: 413 414 template <typename R, typename S> 415 struct isSameType 416 { 417 enum {value = 0}; 418 }; 419 420 template <typename R> 421 struct isSameType<R, R> 422 { 423 enum {value = 1}; 424 }; 425 }; 426 427 428 template <class T> class Matrix44 429 { 430 public: 431 432 //------------------- 433 // Access to elements 434 //------------------- 435 436 T x[4][4]; 437 438 T * operator [] (int i); 439 const T * operator [] (int i) const; 440 441 442 //------------- 443 // Constructors 444 //------------- 445 446 Matrix44 (Uninitialized) {} 447 448 Matrix44 (); 449 // 1 0 0 0 450 // 0 1 0 0 451 // 0 0 1 0 452 // 0 0 0 1 453 454 Matrix44 (T a); 455 // a a a a 456 // a a a a 457 // a a a a 458 // a a a a 459 460 Matrix44 (const T a[4][4]) ; 461 // a[0][0] a[0][1] a[0][2] a[0][3] 462 // a[1][0] a[1][1] a[1][2] a[1][3] 463 // a[2][0] a[2][1] a[2][2] a[2][3] 464 // a[3][0] a[3][1] a[3][2] a[3][3] 465 466 Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h, 467 T i, T j, T k, T l, T m, T n, T o, T p); 468 469 // a b c d 470 // e f g h 471 // i j k l 472 // m n o p 473 474 Matrix44 (Matrix33<T> r, Vec3<T> t); 475 // r r r 0 476 // r r r 0 477 // r r r 0 478 // t t t 1 479 480 481 //-------------------------------- 482 // Copy constructor and assignment 483 //-------------------------------- 484 485 Matrix44 (const Matrix44 &v); 486 template <class S> explicit Matrix44 (const Matrix44<S> &v); 487 488 const Matrix44 & operator = (const Matrix44 &v); 489 const Matrix44 & operator = (T a); 490 491 492 //---------------------- 493 // Compatibility with Sb 494 //---------------------- 495 496 T * getValue (); 497 const T * getValue () const; 498 499 template <class S> 500 void getValue (Matrix44<S> &v) const; 501 template <class S> 502 Matrix44 & setValue (const Matrix44<S> &v); 503 504 template <class S> 505 Matrix44 & setTheMatrix (const Matrix44<S> &v); 506 507 //--------- 508 // Identity 509 //--------- 510 511 void makeIdentity(); 512 513 514 //--------- 515 // Equality 516 //--------- 517 518 bool operator == (const Matrix44 &v) const; 519 bool operator != (const Matrix44 &v) const; 520 521 //----------------------------------------------------------------------- 522 // Compare two matrices and test if they are "approximately equal": 523 // 524 // equalWithAbsError (m, e) 525 // 526 // Returns true if the coefficients of this and m are the same with 527 // an absolute error of no more than e, i.e., for all i, j 528 // 529 // abs (this[i][j] - m[i][j]) <= e 530 // 531 // equalWithRelError (m, e) 532 // 533 // Returns true if the coefficients of this and m are the same with 534 // a relative error of no more than e, i.e., for all i, j 535 // 536 // abs (this[i] - v[i][j]) <= e * abs (this[i][j]) 537 //----------------------------------------------------------------------- 538 539 bool equalWithAbsError (const Matrix44<T> &v, T e) const; 540 bool equalWithRelError (const Matrix44<T> &v, T e) const; 541 542 543 //------------------------ 544 // Component-wise addition 545 //------------------------ 546 547 const Matrix44 & operator += (const Matrix44 &v); 548 const Matrix44 & operator += (T a); 549 Matrix44 operator + (const Matrix44 &v) const; 550 551 552 //--------------------------- 553 // Component-wise subtraction 554 //--------------------------- 555 556 const Matrix44 & operator -= (const Matrix44 &v); 557 const Matrix44 & operator -= (T a); 558 Matrix44 operator - (const Matrix44 &v) const; 559 560 561 //------------------------------------ 562 // Component-wise multiplication by -1 563 //------------------------------------ 564 565 Matrix44 operator - () const; 566 const Matrix44 & negate (); 567 568 569 //------------------------------ 570 // Component-wise multiplication 571 //------------------------------ 572 573 const Matrix44 & operator *= (T a); 574 Matrix44 operator * (T a) const; 575 576 577 //----------------------------------- 578 // Matrix-times-matrix multiplication 579 //----------------------------------- 580 581 const Matrix44 & operator *= (const Matrix44 &v); 582 Matrix44 operator * (const Matrix44 &v) const; 583 584 static void multiply (const Matrix44 &a, // assumes that 585 const Matrix44 &b, // &a != &c and 586 Matrix44 &c); // &b != &c. 587 588 589 //----------------------------------------------------------------- 590 // Vector-times-matrix multiplication; see also the "operator *" 591 // functions defined below. 592 // 593 // m.multVecMatrix(src,dst) implements a homogeneous transformation 594 // by computing Vec4 (src.x, src.y, src.z, 1) * m and dividing by 595 // the result's third element. 596 // 597 // m.multDirMatrix(src,dst) multiplies src by the upper left 3x3 598 // submatrix, ignoring the rest of matrix m. 599 //----------------------------------------------------------------- 600 601 template <class S> 602 void multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const; 603 604 template <class S> 605 void multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const; 606 607 608 //------------------------ 609 // Component-wise division 610 //------------------------ 611 612 const Matrix44 & operator /= (T a); 613 Matrix44 operator / (T a) const; 614 615 616 //------------------ 617 // Transposed matrix 618 //------------------ 619 620 const Matrix44 & transpose (); 621 Matrix44 transposed () const; 622 623 624 //------------------------------------------------------------ 625 // Inverse matrix: If singExc is false, inverting a singular 626 // matrix produces an identity matrix. If singExc is true, 627 // inverting a singular matrix throws a SingMatrixExc. 628 // 629 // inverse() and invert() invert matrices using determinants; 630 // gjInverse() and gjInvert() use the Gauss-Jordan method. 631 // 632 // inverse() and invert() are significantly faster than 633 // gjInverse() and gjInvert(), but the results may be slightly 634 // less accurate. 635 // 636 //------------------------------------------------------------ 637 638 const Matrix44 & invert (bool singExc = false) 639 throw (Iex::MathExc); 640 641 Matrix44<T> inverse (bool singExc = false) const 642 throw (Iex::MathExc); 643 644 const Matrix44 & gjInvert (bool singExc = false) 645 throw (Iex::MathExc); 646 647 Matrix44<T> gjInverse (bool singExc = false) const 648 throw (Iex::MathExc); 649 650 651 //------------------------------------------------ 652 // Calculate the matrix minor of the (r,c) element 653 //------------------------------------------------ 654 655 T minorOf (const int r, const int c) const; 656 657 //--------------------------------------------------- 658 // Build a minor using the specified rows and columns 659 //--------------------------------------------------- 660 661 T fastMinor (const int r0, const int r1, const int r2, 662 const int c0, const int c1, const int c2) const; 663 664 //------------ 665 // Determinant 666 //------------ 667 668 T determinant() const; 669 670 //-------------------------------------------------------- 671 // Set matrix to rotation by XYZ euler angles (in radians) 672 //-------------------------------------------------------- 673 674 template <class S> 675 const Matrix44 & setEulerAngles (const Vec3<S>& r); 676 677 678 //-------------------------------------------------------- 679 // Set matrix to rotation around given axis by given angle 680 //-------------------------------------------------------- 681 682 template <class S> 683 const Matrix44 & setAxisAngle (const Vec3<S>& ax, S ang); 684 685 686 //------------------------------------------- 687 // Rotate the matrix by XYZ euler angles in r 688 //------------------------------------------- 689 690 template <class S> 691 const Matrix44 & rotate (const Vec3<S> &r); 692 693 694 //-------------------------------------------- 695 // Set matrix to scale by given uniform factor 696 //-------------------------------------------- 697 698 const Matrix44 & setScale (T s); 699 700 701 //------------------------------------ 702 // Set matrix to scale by given vector 703 //------------------------------------ 704 705 template <class S> 706 const Matrix44 & setScale (const Vec3<S> &s); 707 708 709 //---------------------- 710 // Scale the matrix by s 711 //---------------------- 712 713 template <class S> 714 const Matrix44 & scale (const Vec3<S> &s); 715 716 717 //------------------------------------------ 718 // Set matrix to translation by given vector 719 //------------------------------------------ 720 721 template <class S> 722 const Matrix44 & setTranslation (const Vec3<S> &t); 723 724 725 //----------------------------- 726 // Return translation component 727 //----------------------------- 728 729 const Vec3<T> translation () const; 730 731 732 //-------------------------- 733 // Translate the matrix by t 734 //-------------------------- 735 736 template <class S> 737 const Matrix44 & translate (const Vec3<S> &t); 738 739 740 //------------------------------------------------------------- 741 // Set matrix to shear by given vector h. The resulting matrix 742 // will shear x for each y coord. by a factor of h[0] ; 743 // will shear x for each z coord. by a factor of h[1] ; 744 // will shear y for each z coord. by a factor of h[2] . 745 //------------------------------------------------------------- 746 747 template <class S> 748 const Matrix44 & setShear (const Vec3<S> &h); 749 750 751 //------------------------------------------------------------ 752 // Set matrix to shear by given factors. The resulting matrix 753 // will shear x for each y coord. by a factor of h.xy ; 754 // will shear x for each z coord. by a factor of h.xz ; 755 // will shear y for each z coord. by a factor of h.yz ; 756 // will shear y for each x coord. by a factor of h.yx ; 757 // will shear z for each x coord. by a factor of h.zx ; 758 // will shear z for each y coord. by a factor of h.zy . 759 //------------------------------------------------------------ 760 761 template <class S> 762 const Matrix44 & setShear (const Shear6<S> &h); 763 764 765 //-------------------------------------------------------- 766 // Shear the matrix by given vector. The composed matrix 767 // will be <shear> * <this>, where the shear matrix ... 768 // will shear x for each y coord. by a factor of h[0] ; 769 // will shear x for each z coord. by a factor of h[1] ; 770 // will shear y for each z coord. by a factor of h[2] . 771 //-------------------------------------------------------- 772 773 template <class S> 774 const Matrix44 & shear (const Vec3<S> &h); 775 776 //-------------------------------------------------------- 777 // Number of the row and column dimensions, since 778 // Matrix44 is a square matrix. 779 //-------------------------------------------------------- 780 781 static unsigned int dimensions() {return 4;} 782 783 784 //------------------------------------------------------------ 785 // Shear the matrix by the given factors. The composed matrix 786 // will be <shear> * <this>, where the shear matrix ... 787 // will shear x for each y coord. by a factor of h.xy ; 788 // will shear x for each z coord. by a factor of h.xz ; 789 // will shear y for each z coord. by a factor of h.yz ; 790 // will shear y for each x coord. by a factor of h.yx ; 791 // will shear z for each x coord. by a factor of h.zx ; 792 // will shear z for each y coord. by a factor of h.zy . 793 //------------------------------------------------------------ 794 795 template <class S> 796 const Matrix44 & shear (const Shear6<S> &h); 797 798 799 //------------------------------------------------- 800 // Limitations of type T (see also class limits<T>) 801 //------------------------------------------------- 802 803 static T baseTypeMin() {return limits<T>::min();} 804 static T baseTypeMax() {return limits<T>::max();} 805 static T baseTypeSmallest() {return limits<T>::smallest();} 806 static T baseTypeEpsilon() {return limits<T>::epsilon();} 807 808 typedef T BaseType; 809 typedef Vec4<T> BaseVecType; 810 811 private: 812 813 template <typename R, typename S> 814 struct isSameType 815 { 816 enum {value = 0}; 817 }; 818 819 template <typename R> 820 struct isSameType<R, R> 821 { 822 enum {value = 1}; 823 }; 824 }; 825 826 827 //-------------- 828 // Stream output 829 //-------------- 830 831 template <class T> 832 std::ostream & operator << (std::ostream & s, const Matrix33<T> &m); 833 834 template <class T> 835 std::ostream & operator << (std::ostream & s, const Matrix44<T> &m); 836 837 838 //--------------------------------------------- 839 // Vector-times-matrix multiplication operators 840 //--------------------------------------------- 841 842 template <class S, class T> 843 const Vec2<S> & operator *= (Vec2<S> &v, const Matrix33<T> &m); 844 845 template <class S, class T> 846 Vec2<S> operator * (const Vec2<S> &v, const Matrix33<T> &m); 847 848 template <class S, class T> 849 const Vec3<S> & operator *= (Vec3<S> &v, const Matrix33<T> &m); 850 851 template <class S, class T> 852 Vec3<S> operator * (const Vec3<S> &v, const Matrix33<T> &m); 853 854 template <class S, class T> 855 const Vec3<S> & operator *= (Vec3<S> &v, const Matrix44<T> &m); 856 857 template <class S, class T> 858 Vec3<S> operator * (const Vec3<S> &v, const Matrix44<T> &m); 859 860 template <class S, class T> 861 const Vec4<S> & operator *= (Vec4<S> &v, const Matrix44<T> &m); 862 863 template <class S, class T> 864 Vec4<S> operator * (const Vec4<S> &v, const Matrix44<T> &m); 865 866 //------------------------- 867 // Typedefs for convenience 868 //------------------------- 869 870 typedef Matrix33 <float> M33f; 871 typedef Matrix33 <double> M33d; 872 typedef Matrix44 <float> M44f; 873 typedef Matrix44 <double> M44d; 874 875 876 //--------------------------- 877 // Implementation of Matrix33 878 //--------------------------- 879 880 template <class T> 881 inline T * 882 Matrix33<T>::operator [] (int i) 883 { 884 return x[i]; 885 } 886 887 template <class T> 888 inline const T * 889 Matrix33<T>::operator [] (int i) const 890 { 891 return x[i]; 892 } 893 894 template <class T> 895 inline 896 Matrix33<T>::Matrix33 () 897 { 898 memset (x, 0, sizeof (x)); 899 x[0][0] = 1; 900 x[1][1] = 1; 901 x[2][2] = 1; 902 } 903 904 template <class T> 905 inline 906 Matrix33<T>::Matrix33 (T a) 907 { 908 x[0][0] = a; 909 x[0][1] = a; 910 x[0][2] = a; 911 x[1][0] = a; 912 x[1][1] = a; 913 x[1][2] = a; 914 x[2][0] = a; 915 x[2][1] = a; 916 x[2][2] = a; 917 } 918 919 template <class T> 920 inline 921 Matrix33<T>::Matrix33 (const T a[3][3]) 922 { 923 memcpy (x, a, sizeof (x)); 924 } 925 926 template <class T> 927 inline 928 Matrix33<T>::Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i) 929 { 930 x[0][0] = a; 931 x[0][1] = b; 932 x[0][2] = c; 933 x[1][0] = d; 934 x[1][1] = e; 935 x[1][2] = f; 936 x[2][0] = g; 937 x[2][1] = h; 938 x[2][2] = i; 939 } 940 941 template <class T> 942 inline 943 Matrix33<T>::Matrix33 (const Matrix33 &v) 944 { 945 memcpy (x, v.x, sizeof (x)); 946 } 947 948 template <class T> 949 template <class S> 950 inline 951 Matrix33<T>::Matrix33 (const Matrix33<S> &v) 952 { 953 x[0][0] = T (v.x[0][0]); 954 x[0][1] = T (v.x[0][1]); 955 x[0][2] = T (v.x[0][2]); 956 x[1][0] = T (v.x[1][0]); 957 x[1][1] = T (v.x[1][1]); 958 x[1][2] = T (v.x[1][2]); 959 x[2][0] = T (v.x[2][0]); 960 x[2][1] = T (v.x[2][1]); 961 x[2][2] = T (v.x[2][2]); 962 } 963 964 template <class T> 965 inline const Matrix33<T> & 966 Matrix33<T>::operator = (const Matrix33 &v) 967 { 968 memcpy (x, v.x, sizeof (x)); 969 return *this; 970 } 971 972 template <class T> 973 inline const Matrix33<T> & 974 Matrix33<T>::operator = (T a) 975 { 976 x[0][0] = a; 977 x[0][1] = a; 978 x[0][2] = a; 979 x[1][0] = a; 980 x[1][1] = a; 981 x[1][2] = a; 982 x[2][0] = a; 983 x[2][1] = a; 984 x[2][2] = a; 985 return *this; 986 } 987 988 template <class T> 989 inline T * 990 Matrix33<T>::getValue () 991 { 992 return (T *) &x[0][0]; 993 } 994 995 template <class T> 996 inline const T * 997 Matrix33<T>::getValue () const 998 { 999 return (const T *) &x[0][0]; 1000 } 1001 1002 template <class T> 1003 template <class S> 1004 inline void 1005 Matrix33<T>::getValue (Matrix33<S> &v) const 1006 { 1007 if (isSameType<S,T>::value) 1008 { 1009 memcpy (v.x, x, sizeof (x)); 1010 } 1011 else 1012 { 1013 v.x[0][0] = x[0][0]; 1014 v.x[0][1] = x[0][1]; 1015 v.x[0][2] = x[0][2]; 1016 v.x[1][0] = x[1][0]; 1017 v.x[1][1] = x[1][1]; 1018 v.x[1][2] = x[1][2]; 1019 v.x[2][0] = x[2][0]; 1020 v.x[2][1] = x[2][1]; 1021 v.x[2][2] = x[2][2]; 1022 } 1023 } 1024 1025 template <class T> 1026 template <class S> 1027 inline Matrix33<T> & 1028 Matrix33<T>::setValue (const Matrix33<S> &v) 1029 { 1030 if (isSameType<S,T>::value) 1031 { 1032 memcpy (x, v.x, sizeof (x)); 1033 } 1034 else 1035 { 1036 x[0][0] = v.x[0][0]; 1037 x[0][1] = v.x[0][1]; 1038 x[0][2] = v.x[0][2]; 1039 x[1][0] = v.x[1][0]; 1040 x[1][1] = v.x[1][1]; 1041 x[1][2] = v.x[1][2]; 1042 x[2][0] = v.x[2][0]; 1043 x[2][1] = v.x[2][1]; 1044 x[2][2] = v.x[2][2]; 1045 } 1046 1047 return *this; 1048 } 1049 1050 template <class T> 1051 template <class S> 1052 inline Matrix33<T> & 1053 Matrix33<T>::setTheMatrix (const Matrix33<S> &v) 1054 { 1055 if (isSameType<S,T>::value) 1056 { 1057 memcpy (x, v.x, sizeof (x)); 1058 } 1059 else 1060 { 1061 x[0][0] = v.x[0][0]; 1062 x[0][1] = v.x[0][1]; 1063 x[0][2] = v.x[0][2]; 1064 x[1][0] = v.x[1][0]; 1065 x[1][1] = v.x[1][1]; 1066 x[1][2] = v.x[1][2]; 1067 x[2][0] = v.x[2][0]; 1068 x[2][1] = v.x[2][1]; 1069 x[2][2] = v.x[2][2]; 1070 } 1071 1072 return *this; 1073 } 1074 1075 template <class T> 1076 inline void 1077 Matrix33<T>::makeIdentity() 1078 { 1079 memset (x, 0, sizeof (x)); 1080 x[0][0] = 1; 1081 x[1][1] = 1; 1082 x[2][2] = 1; 1083 } 1084 1085 template <class T> 1086 bool 1087 Matrix33<T>::operator == (const Matrix33 &v) const 1088 { 1089 return x[0][0] == v.x[0][0] && 1090 x[0][1] == v.x[0][1] && 1091 x[0][2] == v.x[0][2] && 1092 x[1][0] == v.x[1][0] && 1093 x[1][1] == v.x[1][1] && 1094 x[1][2] == v.x[1][2] && 1095 x[2][0] == v.x[2][0] && 1096 x[2][1] == v.x[2][1] && 1097 x[2][2] == v.x[2][2]; 1098 } 1099 1100 template <class T> 1101 bool 1102 Matrix33<T>::operator != (const Matrix33 &v) const 1103 { 1104 return x[0][0] != v.x[0][0] || 1105 x[0][1] != v.x[0][1] || 1106 x[0][2] != v.x[0][2] || 1107 x[1][0] != v.x[1][0] || 1108 x[1][1] != v.x[1][1] || 1109 x[1][2] != v.x[1][2] || 1110 x[2][0] != v.x[2][0] || 1111 x[2][1] != v.x[2][1] || 1112 x[2][2] != v.x[2][2]; 1113 } 1114 1115 template <class T> 1116 bool 1117 Matrix33<T>::equalWithAbsError (const Matrix33<T> &m, T e) const 1118 { 1119 for (int i = 0; i < 3; i++) 1120 for (int j = 0; j < 3; j++) 1121 if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e)) 1122 return false; 1123 1124 return true; 1125 } 1126 1127 template <class T> 1128 bool 1129 Matrix33<T>::equalWithRelError (const Matrix33<T> &m, T e) const 1130 { 1131 for (int i = 0; i < 3; i++) 1132 for (int j = 0; j < 3; j++) 1133 if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e)) 1134 return false; 1135 1136 return true; 1137 } 1138 1139 template <class T> 1140 const Matrix33<T> & 1141 Matrix33<T>::operator += (const Matrix33<T> &v) 1142 { 1143 x[0][0] += v.x[0][0]; 1144 x[0][1] += v.x[0][1]; 1145 x[0][2] += v.x[0][2]; 1146 x[1][0] += v.x[1][0]; 1147 x[1][1] += v.x[1][1]; 1148 x[1][2] += v.x[1][2]; 1149 x[2][0] += v.x[2][0]; 1150 x[2][1] += v.x[2][1]; 1151 x[2][2] += v.x[2][2]; 1152 1153 return *this; 1154 } 1155 1156 template <class T> 1157 const Matrix33<T> & 1158 Matrix33<T>::operator += (T a) 1159 { 1160 x[0][0] += a; 1161 x[0][1] += a; 1162 x[0][2] += a; 1163 x[1][0] += a; 1164 x[1][1] += a; 1165 x[1][2] += a; 1166 x[2][0] += a; 1167 x[2][1] += a; 1168 x[2][2] += a; 1169 1170 return *this; 1171 } 1172 1173 template <class T> 1174 Matrix33<T> 1175 Matrix33<T>::operator + (const Matrix33<T> &v) const 1176 { 1177 return Matrix33 (x[0][0] + v.x[0][0], 1178 x[0][1] + v.x[0][1], 1179 x[0][2] + v.x[0][2], 1180 x[1][0] + v.x[1][0], 1181 x[1][1] + v.x[1][1], 1182 x[1][2] + v.x[1][2], 1183 x[2][0] + v.x[2][0], 1184 x[2][1] + v.x[2][1], 1185 x[2][2] + v.x[2][2]); 1186 } 1187 1188 template <class T> 1189 const Matrix33<T> & 1190 Matrix33<T>::operator -= (const Matrix33<T> &v) 1191 { 1192 x[0][0] -= v.x[0][0]; 1193 x[0][1] -= v.x[0][1]; 1194 x[0][2] -= v.x[0][2]; 1195 x[1][0] -= v.x[1][0]; 1196 x[1][1] -= v.x[1][1]; 1197 x[1][2] -= v.x[1][2]; 1198 x[2][0] -= v.x[2][0]; 1199 x[2][1] -= v.x[2][1]; 1200 x[2][2] -= v.x[2][2]; 1201 1202 return *this; 1203 } 1204 1205 template <class T> 1206 const Matrix33<T> & 1207 Matrix33<T>::operator -= (T a) 1208 { 1209 x[0][0] -= a; 1210 x[0][1] -= a; 1211 x[0][2] -= a; 1212 x[1][0] -= a; 1213 x[1][1] -= a; 1214 x[1][2] -= a; 1215 x[2][0] -= a; 1216 x[2][1] -= a; 1217 x[2][2] -= a; 1218 1219 return *this; 1220 } 1221 1222 template <class T> 1223 Matrix33<T> 1224 Matrix33<T>::operator - (const Matrix33<T> &v) const 1225 { 1226 return Matrix33 (x[0][0] - v.x[0][0], 1227 x[0][1] - v.x[0][1], 1228 x[0][2] - v.x[0][2], 1229 x[1][0] - v.x[1][0], 1230 x[1][1] - v.x[1][1], 1231 x[1][2] - v.x[1][2], 1232 x[2][0] - v.x[2][0], 1233 x[2][1] - v.x[2][1], 1234 x[2][2] - v.x[2][2]); 1235 } 1236 1237 template <class T> 1238 Matrix33<T> 1239 Matrix33<T>::operator - () const 1240 { 1241 return Matrix33 (-x[0][0], 1242 -x[0][1], 1243 -x[0][2], 1244 -x[1][0], 1245 -x[1][1], 1246 -x[1][2], 1247 -x[2][0], 1248 -x[2][1], 1249 -x[2][2]); 1250 } 1251 1252 template <class T> 1253 const Matrix33<T> & 1254 Matrix33<T>::negate () 1255 { 1256 x[0][0] = -x[0][0]; 1257 x[0][1] = -x[0][1]; 1258 x[0][2] = -x[0][2]; 1259 x[1][0] = -x[1][0]; 1260 x[1][1] = -x[1][1]; 1261 x[1][2] = -x[1][2]; 1262 x[2][0] = -x[2][0]; 1263 x[2][1] = -x[2][1]; 1264 x[2][2] = -x[2][2]; 1265 1266 return *this; 1267 } 1268 1269 template <class T> 1270 const Matrix33<T> & 1271 Matrix33<T>::operator *= (T a) 1272 { 1273 x[0][0] *= a; 1274 x[0][1] *= a; 1275 x[0][2] *= a; 1276 x[1][0] *= a; 1277 x[1][1] *= a; 1278 x[1][2] *= a; 1279 x[2][0] *= a; 1280 x[2][1] *= a; 1281 x[2][2] *= a; 1282 1283 return *this; 1284 } 1285 1286 template <class T> 1287 Matrix33<T> 1288 Matrix33<T>::operator * (T a) const 1289 { 1290 return Matrix33 (x[0][0] * a, 1291 x[0][1] * a, 1292 x[0][2] * a, 1293 x[1][0] * a, 1294 x[1][1] * a, 1295 x[1][2] * a, 1296 x[2][0] * a, 1297 x[2][1] * a, 1298 x[2][2] * a); 1299 } 1300 1301 template <class T> 1302 inline Matrix33<T> 1303 operator * (T a, const Matrix33<T> &v) 1304 { 1305 return v * a; 1306 } 1307 1308 template <class T> 1309 const Matrix33<T> & 1310 Matrix33<T>::operator *= (const Matrix33<T> &v) 1311 { 1312 Matrix33 tmp (T (0)); 1313 1314 for (int i = 0; i < 3; i++) 1315 for (int j = 0; j < 3; j++) 1316 for (int k = 0; k < 3; k++) 1317 tmp.x[i][j] += x[i][k] * v.x[k][j]; 1318 1319 *this = tmp; 1320 return *this; 1321 } 1322 1323 template <class T> 1324 Matrix33<T> 1325 Matrix33<T>::operator * (const Matrix33<T> &v) const 1326 { 1327 Matrix33 tmp (T (0)); 1328 1329 for (int i = 0; i < 3; i++) 1330 for (int j = 0; j < 3; j++) 1331 for (int k = 0; k < 3; k++) 1332 tmp.x[i][j] += x[i][k] * v.x[k][j]; 1333 1334 return tmp; 1335 } 1336 1337 template <class T> 1338 template <class S> 1339 void 1340 Matrix33<T>::multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const 1341 { 1342 S a, b, w; 1343 1344 a = src[0] * x[0][0] + src[1] * x[1][0] + x[2][0]; 1345 b = src[0] * x[0][1] + src[1] * x[1][1] + x[2][1]; 1346 w = src[0] * x[0][2] + src[1] * x[1][2] + x[2][2]; 1347 1348 dst.x = a / w; 1349 dst.y = b / w; 1350 } 1351 1352 template <class T> 1353 template <class S> 1354 void 1355 Matrix33<T>::multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const 1356 { 1357 S a, b; 1358 1359 a = src[0] * x[0][0] + src[1] * x[1][0]; 1360 b = src[0] * x[0][1] + src[1] * x[1][1]; 1361 1362 dst.x = a; 1363 dst.y = b; 1364 } 1365 1366 template <class T> 1367 const Matrix33<T> & 1368 Matrix33<T>::operator /= (T a) 1369 { 1370 x[0][0] /= a; 1371 x[0][1] /= a; 1372 x[0][2] /= a; 1373 x[1][0] /= a; 1374 x[1][1] /= a; 1375 x[1][2] /= a; 1376 x[2][0] /= a; 1377 x[2][1] /= a; 1378 x[2][2] /= a; 1379 1380 return *this; 1381 } 1382 1383 template <class T> 1384 Matrix33<T> 1385 Matrix33<T>::operator / (T a) const 1386 { 1387 return Matrix33 (x[0][0] / a, 1388 x[0][1] / a, 1389 x[0][2] / a, 1390 x[1][0] / a, 1391 x[1][1] / a, 1392 x[1][2] / a, 1393 x[2][0] / a, 1394 x[2][1] / a, 1395 x[2][2] / a); 1396 } 1397 1398 template <class T> 1399 const Matrix33<T> & 1400 Matrix33<T>::transpose () 1401 { 1402 Matrix33 tmp (x[0][0], 1403 x[1][0], 1404 x[2][0], 1405 x[0][1], 1406 x[1][1], 1407 x[2][1], 1408 x[0][2], 1409 x[1][2], 1410 x[2][2]); 1411 *this = tmp; 1412 return *this; 1413 } 1414 1415 template <class T> 1416 Matrix33<T> 1417 Matrix33<T>::transposed () const 1418 { 1419 return Matrix33 (x[0][0], 1420 x[1][0], 1421 x[2][0], 1422 x[0][1], 1423 x[1][1], 1424 x[2][1], 1425 x[0][2], 1426 x[1][2], 1427 x[2][2]); 1428 } 1429 1430 template <class T> 1431 const Matrix33<T> & 1432 Matrix33<T>::gjInvert (bool singExc) throw (Iex::MathExc) 1433 { 1434 *this = gjInverse (singExc); 1435 return *this; 1436 } 1437 1438 template <class T> 1439 Matrix33<T> 1440 Matrix33<T>::gjInverse (bool singExc) const throw (Iex::MathExc) 1441 { 1442 int i, j, k; 1443 Matrix33 s; 1444 Matrix33 t (*this); 1445 1446 // Forward elimination 1447 1448 for (i = 0; i < 2 ; i++) 1449 { 1450 int pivot = i; 1451 1452 T pivotsize = t[i][i]; 1453 1454 if (pivotsize < 0) 1455 pivotsize = -pivotsize; 1456 1457 for (j = i + 1; j < 3; j++) 1458 { 1459 T tmp = t[j][i]; 1460 1461 if (tmp < 0) 1462 tmp = -tmp; 1463 1464 if (tmp > pivotsize) 1465 { 1466 pivot = j; 1467 pivotsize = tmp; 1468 } 1469 } 1470 1471 if (pivotsize == 0) 1472 { 1473 if (singExc) 1474 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix."); 1475 1476 return Matrix33(); 1477 } 1478 1479 if (pivot != i) 1480 { 1481 for (j = 0; j < 3; j++) 1482 { 1483 T tmp; 1484 1485 tmp = t[i][j]; 1486 t[i][j] = t[pivot][j]; 1487 t[pivot][j] = tmp; 1488 1489 tmp = s[i][j]; 1490 s[i][j] = s[pivot][j]; 1491 s[pivot][j] = tmp; 1492 } 1493 } 1494 1495 for (j = i + 1; j < 3; j++) 1496 { 1497 T f = t[j][i] / t[i][i]; 1498 1499 for (k = 0; k < 3; k++) 1500 { 1501 t[j][k] -= f * t[i][k]; 1502 s[j][k] -= f * s[i][k]; 1503 } 1504 } 1505 } 1506 1507 // Backward substitution 1508 1509 for (i = 2; i >= 0; --i) 1510 { 1511 T f; 1512 1513 if ((f = t[i][i]) == 0) 1514 { 1515 if (singExc) 1516 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix."); 1517 1518 return Matrix33(); 1519 } 1520 1521 for (j = 0; j < 3; j++) 1522 { 1523 t[i][j] /= f; 1524 s[i][j] /= f; 1525 } 1526 1527 for (j = 0; j < i; j++) 1528 { 1529 f = t[j][i]; 1530 1531 for (k = 0; k < 3; k++) 1532 { 1533 t[j][k] -= f * t[i][k]; 1534 s[j][k] -= f * s[i][k]; 1535 } 1536 } 1537 } 1538 1539 return s; 1540 } 1541 1542 template <class T> 1543 const Matrix33<T> & 1544 Matrix33<T>::invert (bool singExc) throw (Iex::MathExc) 1545 { 1546 *this = inverse (singExc); 1547 return *this; 1548 } 1549 1550 template <class T> 1551 Matrix33<T> 1552 Matrix33<T>::inverse (bool singExc) const throw (Iex::MathExc) 1553 { 1554 if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1) 1555 { 1556 Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2], 1557 x[2][1] * x[0][2] - x[0][1] * x[2][2], 1558 x[0][1] * x[1][2] - x[1][1] * x[0][2], 1559 1560 x[2][0] * x[1][2] - x[1][0] * x[2][2], 1561 x[0][0] * x[2][2] - x[2][0] * x[0][2], 1562 x[1][0] * x[0][2] - x[0][0] * x[1][2], 1563 1564 x[1][0] * x[2][1] - x[2][0] * x[1][1], 1565 x[2][0] * x[0][1] - x[0][0] * x[2][1], 1566 x[0][0] * x[1][1] - x[1][0] * x[0][1]); 1567 1568 T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0]; 1569 1570 if (Imath::abs (r) >= 1) 1571 { 1572 for (int i = 0; i < 3; ++i) 1573 { 1574 for (int j = 0; j < 3; ++j) 1575 { 1576 s[i][j] /= r; 1577 } 1578 } 1579 } 1580 else 1581 { 1582 T mr = Imath::abs (r) / limits<T>::smallest(); 1583 1584 for (int i = 0; i < 3; ++i) 1585 { 1586 for (int j = 0; j < 3; ++j) 1587 { 1588 if (mr > Imath::abs (s[i][j])) 1589 { 1590 s[i][j] /= r; 1591 } 1592 else 1593 { 1594 if (singExc) 1595 throw SingMatrixExc ("Cannot invert " 1596 "singular matrix."); 1597 return Matrix33(); 1598 } 1599 } 1600 } 1601 } 1602 1603 return s; 1604 } 1605 else 1606 { 1607 Matrix33 s ( x[1][1], 1608 -x[0][1], 1609 0, 1610 1611 -x[1][0], 1612 x[0][0], 1613 0, 1614 1615 0, 1616 0, 1617 1); 1618 1619 T r = x[0][0] * x[1][1] - x[1][0] * x[0][1]; 1620 1621 if (Imath::abs (r) >= 1) 1622 { 1623 for (int i = 0; i < 2; ++i) 1624 { 1625 for (int j = 0; j < 2; ++j) 1626 { 1627 s[i][j] /= r; 1628 } 1629 } 1630 } 1631 else 1632 { 1633 T mr = Imath::abs (r) / limits<T>::smallest(); 1634 1635 for (int i = 0; i < 2; ++i) 1636 { 1637 for (int j = 0; j < 2; ++j) 1638 { 1639 if (mr > Imath::abs (s[i][j])) 1640 { 1641 s[i][j] /= r; 1642 } 1643 else 1644 { 1645 if (singExc) 1646 throw SingMatrixExc ("Cannot invert " 1647 "singular matrix."); 1648 return Matrix33(); 1649 } 1650 } 1651 } 1652 } 1653 1654 s[2][0] = -x[2][0] * s[0][0] - x[2][1] * s[1][0]; 1655 s[2][1] = -x[2][0] * s[0][1] - x[2][1] * s[1][1]; 1656 1657 return s; 1658 } 1659 } 1660 1661 template <class T> 1662 inline T 1663 Matrix33<T>::minorOf (const int r, const int c) const 1664 { 1665 int r0 = 0 + (r < 1 ? 1 : 0); 1666 int r1 = 1 + (r < 2 ? 1 : 0); 1667 int c0 = 0 + (c < 1 ? 1 : 0); 1668 int c1 = 1 + (c < 2 ? 1 : 0); 1669 1670 return x[r0][c0]*x[r1][c1] - x[r1][c0]*x[r0][c1]; 1671 } 1672 1673 template <class T> 1674 inline T 1675 Matrix33<T>::fastMinor( const int r0, const int r1, 1676 const int c0, const int c1) const 1677 { 1678 return x[r0][c0]*x[r1][c1] - x[r0][c1]*x[r1][c0]; 1679 } 1680 1681 template <class T> 1682 inline T 1683 Matrix33<T>::determinant () const 1684 { 1685 return x[0][0]*(x[1][1]*x[2][2] - x[1][2]*x[2][1]) + 1686 x[0][1]*(x[1][2]*x[2][0] - x[1][0]*x[2][2]) + 1687 x[0][2]*(x[1][0]*x[2][1] - x[1][1]*x[2][0]); 1688 } 1689 1690 template <class T> 1691 template <class S> 1692 const Matrix33<T> & 1693 Matrix33<T>::setRotation (S r) 1694 { 1695 S cos_r, sin_r; 1696 1697 cos_r = Math<T>::cos (r); 1698 sin_r = Math<T>::sin (r); 1699 1700 x[0][0] = cos_r; 1701 x[0][1] = sin_r; 1702 x[0][2] = 0; 1703 1704 x[1][0] = -sin_r; 1705 x[1][1] = cos_r; 1706 x[1][2] = 0; 1707 1708 x[2][0] = 0; 1709 x[2][1] = 0; 1710 x[2][2] = 1; 1711 1712 return *this; 1713 } 1714 1715 template <class T> 1716 template <class S> 1717 const Matrix33<T> & 1718 Matrix33<T>::rotate (S r) 1719 { 1720 *this *= Matrix33<T>().setRotation (r); 1721 return *this; 1722 } 1723 1724 template <class T> 1725 const Matrix33<T> & 1726 Matrix33<T>::setScale (T s) 1727 { 1728 memset (x, 0, sizeof (x)); 1729 x[0][0] = s; 1730 x[1][1] = s; 1731 x[2][2] = 1; 1732 1733 return *this; 1734 } 1735 1736 template <class T> 1737 template <class S> 1738 const Matrix33<T> & 1739 Matrix33<T>::setScale (const Vec2<S> &s) 1740 { 1741 memset (x, 0, sizeof (x)); 1742 x[0][0] = s[0]; 1743 x[1][1] = s[1]; 1744 x[2][2] = 1; 1745 1746 return *this; 1747 } 1748 1749 template <class T> 1750 template <class S> 1751 const Matrix33<T> & 1752 Matrix33<T>::scale (const Vec2<S> &s) 1753 { 1754 x[0][0] *= s[0]; 1755 x[0][1] *= s[0]; 1756 x[0][2] *= s[0]; 1757 1758 x[1][0] *= s[1]; 1759 x[1][1] *= s[1]; 1760 x[1][2] *= s[1]; 1761 1762 return *this; 1763 } 1764 1765 template <class T> 1766 template <class S> 1767 const Matrix33<T> & 1768 Matrix33<T>::setTranslation (const Vec2<S> &t) 1769 { 1770 x[0][0] = 1; 1771 x[0][1] = 0; 1772 x[0][2] = 0; 1773 1774 x[1][0] = 0; 1775 x[1][1] = 1; 1776 x[1][2] = 0; 1777 1778 x[2][0] = t[0]; 1779 x[2][1] = t[1]; 1780 x[2][2] = 1; 1781 1782 return *this; 1783 } 1784 1785 template <class T> 1786 inline Vec2<T> 1787 Matrix33<T>::translation () const 1788 { 1789 return Vec2<T> (x[2][0], x[2][1]); 1790 } 1791 1792 template <class T> 1793 template <class S> 1794 const Matrix33<T> & 1795 Matrix33<T>::translate (const Vec2<S> &t) 1796 { 1797 x[2][0] += t[0] * x[0][0] + t[1] * x[1][0]; 1798 x[2][1] += t[0] * x[0][1] + t[1] * x[1][1]; 1799 x[2][2] += t[0] * x[0][2] + t[1] * x[1][2]; 1800 1801 return *this; 1802 } 1803 1804 template <class T> 1805 template <class S> 1806 const Matrix33<T> & 1807 Matrix33<T>::setShear (const S &xy) 1808 { 1809 x[0][0] = 1; 1810 x[0][1] = 0; 1811 x[0][2] = 0; 1812 1813 x[1][0] = xy; 1814 x[1][1] = 1; 1815 x[1][2] = 0; 1816 1817 x[2][0] = 0; 1818 x[2][1] = 0; 1819 x[2][2] = 1; 1820 1821 return *this; 1822 } 1823 1824 template <class T> 1825 template <class S> 1826 const Matrix33<T> & 1827 Matrix33<T>::setShear (const Vec2<S> &h) 1828 { 1829 x[0][0] = 1; 1830 x[0][1] = h[1]; 1831 x[0][2] = 0; 1832 1833 x[1][0] = h[0]; 1834 x[1][1] = 1; 1835 x[1][2] = 0; 1836 1837 x[2][0] = 0; 1838 x[2][1] = 0; 1839 x[2][2] = 1; 1840 1841 return *this; 1842 } 1843 1844 template <class T> 1845 template <class S> 1846 const Matrix33<T> & 1847 Matrix33<T>::shear (const S &xy) 1848 { 1849 // 1850 // In this case, we don't need a temp. copy of the matrix 1851 // because we never use a value on the RHS after we've 1852 // changed it on the LHS. 1853 // 1854 1855 x[1][0] += xy * x[0][0]; 1856 x[1][1] += xy * x[0][1]; 1857 x[1][2] += xy * x[0][2]; 1858 1859 return *this; 1860 } 1861 1862 template <class T> 1863 template <class S> 1864 const Matrix33<T> & 1865 Matrix33<T>::shear (const Vec2<S> &h) 1866 { 1867 Matrix33<T> P (*this); 1868 1869 x[0][0] = P[0][0] + h[1] * P[1][0]; 1870 x[0][1] = P[0][1] + h[1] * P[1][1]; 1871 x[0][2] = P[0][2] + h[1] * P[1][2]; 1872 1873 x[1][0] = P[1][0] + h[0] * P[0][0]; 1874 x[1][1] = P[1][1] + h[0] * P[0][1]; 1875 x[1][2] = P[1][2] + h[0] * P[0][2]; 1876 1877 return *this; 1878 } 1879 1880 1881 //--------------------------- 1882 // Implementation of Matrix44 1883 //--------------------------- 1884 1885 template <class T> 1886 inline T * 1887 Matrix44<T>::operator [] (int i) 1888 { 1889 return x[i]; 1890 } 1891 1892 template <class T> 1893 inline const T * 1894 Matrix44<T>::operator [] (int i) const 1895 { 1896 return x[i]; 1897 } 1898 1899 template <class T> 1900 inline 1901 Matrix44<T>::Matrix44 () 1902 { 1903 memset (x, 0, sizeof (x)); 1904 x[0][0] = 1; 1905 x[1][1] = 1; 1906 x[2][2] = 1; 1907 x[3][3] = 1; 1908 } 1909 1910 template <class T> 1911 inline 1912 Matrix44<T>::Matrix44 (T a) 1913 { 1914 x[0][0] = a; 1915 x[0][1] = a; 1916 x[0][2] = a; 1917 x[0][3] = a; 1918 x[1][0] = a; 1919 x[1][1] = a; 1920 x[1][2] = a; 1921 x[1][3] = a; 1922 x[2][0] = a; 1923 x[2][1] = a; 1924 x[2][2] = a; 1925 x[2][3] = a; 1926 x[3][0] = a; 1927 x[3][1] = a; 1928 x[3][2] = a; 1929 x[3][3] = a; 1930 } 1931 1932 template <class T> 1933 inline 1934 Matrix44<T>::Matrix44 (const T a[4][4]) 1935 { 1936 memcpy (x, a, sizeof (x)); 1937 } 1938 1939 template <class T> 1940 inline 1941 Matrix44<T>::Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h, 1942 T i, T j, T k, T l, T m, T n, T o, T p) 1943 { 1944 x[0][0] = a; 1945 x[0][1] = b; 1946 x[0][2] = c; 1947 x[0][3] = d; 1948 x[1][0] = e; 1949 x[1][1] = f; 1950 x[1][2] = g; 1951 x[1][3] = h; 1952 x[2][0] = i; 1953 x[2][1] = j; 1954 x[2][2] = k; 1955 x[2][3] = l; 1956 x[3][0] = m; 1957 x[3][1] = n; 1958 x[3][2] = o; 1959 x[3][3] = p; 1960 } 1961 1962 1963 template <class T> 1964 inline 1965 Matrix44<T>::Matrix44 (Matrix33<T> r, Vec3<T> t) 1966 { 1967 x[0][0] = r[0][0]; 1968 x[0][1] = r[0][1]; 1969 x[0][2] = r[0][2]; 1970 x[0][3] = 0; 1971 x[1][0] = r[1][0]; 1972 x[1][1] = r[1][1]; 1973 x[1][2] = r[1][2]; 1974 x[1][3] = 0; 1975 x[2][0] = r[2][0]; 1976 x[2][1] = r[2][1]; 1977 x[2][2] = r[2][2]; 1978 x[2][3] = 0; 1979 x[3][0] = t[0]; 1980 x[3][1] = t[1]; 1981 x[3][2] = t[2]; 1982 x[3][3] = 1; 1983 } 1984 1985 template <class T> 1986 inline 1987 Matrix44<T>::Matrix44 (const Matrix44 &v) 1988 { 1989 x[0][0] = v.x[0][0]; 1990 x[0][1] = v.x[0][1]; 1991 x[0][2] = v.x[0][2]; 1992 x[0][3] = v.x[0][3]; 1993 x[1][0] = v.x[1][0]; 1994 x[1][1] = v.x[1][1]; 1995 x[1][2] = v.x[1][2]; 1996 x[1][3] = v.x[1][3]; 1997 x[2][0] = v.x[2][0]; 1998 x[2][1] = v.x[2][1]; 1999 x[2][2] = v.x[2][2]; 2000 x[2][3] = v.x[2][3]; 2001 x[3][0] = v.x[3][0]; 2002 x[3][1] = v.x[3][1]; 2003 x[3][2] = v.x[3][2]; 2004 x[3][3] = v.x[3][3]; 2005 } 2006 2007 template <class T> 2008 template <class S> 2009 inline 2010 Matrix44<T>::Matrix44 (const Matrix44<S> &v) 2011 { 2012 x[0][0] = T (v.x[0][0]); 2013 x[0][1] = T (v.x[0][1]); 2014 x[0][2] = T (v.x[0][2]); 2015 x[0][3] = T (v.x[0][3]); 2016 x[1][0] = T (v.x[1][0]); 2017 x[1][1] = T (v.x[1][1]); 2018 x[1][2] = T (v.x[1][2]); 2019 x[1][3] = T (v.x[1][3]); 2020 x[2][0] = T (v.x[2][0]); 2021 x[2][1] = T (v.x[2][1]); 2022 x[2][2] = T (v.x[2][2]); 2023 x[2][3] = T (v.x[2][3]); 2024 x[3][0] = T (v.x[3][0]); 2025 x[3][1] = T (v.x[3][1]); 2026 x[3][2] = T (v.x[3][2]); 2027 x[3][3] = T (v.x[3][3]); 2028 } 2029 2030 template <class T> 2031 inline const Matrix44<T> & 2032 Matrix44<T>::operator = (const Matrix44 &v) 2033 { 2034 x[0][0] = v.x[0][0]; 2035 x[0][1] = v.x[0][1]; 2036 x[0][2] = v.x[0][2]; 2037 x[0][3] = v.x[0][3]; 2038 x[1][0] = v.x[1][0]; 2039 x[1][1] = v.x[1][1]; 2040 x[1][2] = v.x[1][2]; 2041 x[1][3] = v.x[1][3]; 2042 x[2][0] = v.x[2][0]; 2043 x[2][1] = v.x[2][1]; 2044 x[2][2] = v.x[2][2]; 2045 x[2][3] = v.x[2][3]; 2046 x[3][0] = v.x[3][0]; 2047 x[3][1] = v.x[3][1]; 2048 x[3][2] = v.x[3][2]; 2049 x[3][3] = v.x[3][3]; 2050 return *this; 2051 } 2052 2053 template <class T> 2054 inline const Matrix44<T> & 2055 Matrix44<T>::operator = (T a) 2056 { 2057 x[0][0] = a; 2058 x[0][1] = a; 2059 x[0][2] = a; 2060 x[0][3] = a; 2061 x[1][0] = a; 2062 x[1][1] = a; 2063 x[1][2] = a; 2064 x[1][3] = a; 2065 x[2][0] = a; 2066 x[2][1] = a; 2067 x[2][2] = a; 2068 x[2][3] = a; 2069 x[3][0] = a; 2070 x[3][1] = a; 2071 x[3][2] = a; 2072 x[3][3] = a; 2073 return *this; 2074 } 2075 2076 template <class T> 2077 inline T * 2078 Matrix44<T>::getValue () 2079 { 2080 return (T *) &x[0][0]; 2081 } 2082 2083 template <class T> 2084 inline const T * 2085 Matrix44<T>::getValue () const 2086 { 2087 return (const T *) &x[0][0]; 2088 } 2089 2090 template <class T> 2091 template <class S> 2092 inline void 2093 Matrix44<T>::getValue (Matrix44<S> &v) const 2094 { 2095 if (isSameType<S,T>::value) 2096 { 2097 memcpy (v.x, x, sizeof (x)); 2098 } 2099 else 2100 { 2101 v.x[0][0] = x[0][0]; 2102 v.x[0][1] = x[0][1]; 2103 v.x[0][2] = x[0][2]; 2104 v.x[0][3] = x[0][3]; 2105 v.x[1][0] = x[1][0]; 2106 v.x[1][1] = x[1][1]; 2107 v.x[1][2] = x[1][2]; 2108 v.x[1][3] = x[1][3]; 2109 v.x[2][0] = x[2][0]; 2110 v.x[2][1] = x[2][1]; 2111 v.x[2][2] = x[2][2]; 2112 v.x[2][3] = x[2][3]; 2113 v.x[3][0] = x[3][0]; 2114 v.x[3][1] = x[3][1]; 2115 v.x[3][2] = x[3][2]; 2116 v.x[3][3] = x[3][3]; 2117 } 2118 } 2119 2120 template <class T> 2121 template <class S> 2122 inline Matrix44<T> & 2123 Matrix44<T>::setValue (const Matrix44<S> &v) 2124 { 2125 if (isSameType<S,T>::value) 2126 { 2127 memcpy (x, v.x, sizeof (x)); 2128 } 2129 else 2130 { 2131 x[0][0] = v.x[0][0]; 2132 x[0][1] = v.x[0][1]; 2133 x[0][2] = v.x[0][2]; 2134 x[0][3] = v.x[0][3]; 2135 x[1][0] = v.x[1][0]; 2136 x[1][1] = v.x[1][1]; 2137 x[1][2] = v.x[1][2]; 2138 x[1][3] = v.x[1][3]; 2139 x[2][0] = v.x[2][0]; 2140 x[2][1] = v.x[2][1]; 2141 x[2][2] = v.x[2][2]; 2142 x[2][3] = v.x[2][3]; 2143 x[3][0] = v.x[3][0]; 2144 x[3][1] = v.x[3][1]; 2145 x[3][2] = v.x[3][2]; 2146 x[3][3] = v.x[3][3]; 2147 } 2148 2149 return *this; 2150 } 2151 2152 template <class T> 2153 template <class S> 2154 inline Matrix44<T> & 2155 Matrix44<T>::setTheMatrix (const Matrix44<S> &v) 2156 { 2157 if (isSameType<S,T>::value) 2158 { 2159 memcpy (x, v.x, sizeof (x)); 2160 } 2161 else 2162 { 2163 x[0][0] = v.x[0][0]; 2164 x[0][1] = v.x[0][1]; 2165 x[0][2] = v.x[0][2]; 2166 x[0][3] = v.x[0][3]; 2167 x[1][0] = v.x[1][0]; 2168 x[1][1] = v.x[1][1]; 2169 x[1][2] = v.x[1][2]; 2170 x[1][3] = v.x[1][3]; 2171 x[2][0] = v.x[2][0]; 2172 x[2][1] = v.x[2][1]; 2173 x[2][2] = v.x[2][2]; 2174 x[2][3] = v.x[2][3]; 2175 x[3][0] = v.x[3][0]; 2176 x[3][1] = v.x[3][1]; 2177 x[3][2] = v.x[3][2]; 2178 x[3][3] = v.x[3][3]; 2179 } 2180 2181 return *this; 2182 } 2183 2184 template <class T> 2185 inline void 2186 Matrix44<T>::makeIdentity() 2187 { 2188 memset (x, 0, sizeof (x)); 2189 x[0][0] = 1; 2190 x[1][1] = 1; 2191 x[2][2] = 1; 2192 x[3][3] = 1; 2193 } 2194 2195 template <class T> 2196 bool 2197 Matrix44<T>::operator == (const Matrix44 &v) const 2198 { 2199 return x[0][0] == v.x[0][0] && 2200 x[0][1] == v.x[0][1] && 2201 x[0][2] == v.x[0][2] && 2202 x[0][3] == v.x[0][3] && 2203 x[1][0] == v.x[1][0] && 2204 x[1][1] == v.x[1][1] && 2205 x[1][2] == v.x[1][2] && 2206 x[1][3] == v.x[1][3] && 2207 x[2][0] == v.x[2][0] && 2208 x[2][1] == v.x[2][1] && 2209 x[2][2] == v.x[2][2] && 2210 x[2][3] == v.x[2][3] && 2211 x[3][0] == v.x[3][0] && 2212 x[3][1] == v.x[3][1] && 2213 x[3][2] == v.x[3][2] && 2214 x[3][3] == v.x[3][3]; 2215 } 2216 2217 template <class T> 2218 bool 2219 Matrix44<T>::operator != (const Matrix44 &v) const 2220 { 2221 return x[0][0] != v.x[0][0] || 2222 x[0][1] != v.x[0][1] || 2223 x[0][2] != v.x[0][2] || 2224 x[0][3] != v.x[0][3] || 2225 x[1][0] != v.x[1][0] || 2226 x[1][1] != v.x[1][1] || 2227 x[1][2] != v.x[1][2] || 2228 x[1][3] != v.x[1][3] || 2229 x[2][0] != v.x[2][0] || 2230 x[2][1] != v.x[2][1] || 2231 x[2][2] != v.x[2][2] || 2232 x[2][3] != v.x[2][3] || 2233 x[3][0] != v.x[3][0] || 2234 x[3][1] != v.x[3][1] || 2235 x[3][2] != v.x[3][2] || 2236 x[3][3] != v.x[3][3]; 2237 } 2238 2239 template <class T> 2240 bool 2241 Matrix44<T>::equalWithAbsError (const Matrix44<T> &m, T e) const 2242 { 2243 for (int i = 0; i < 4; i++) 2244 for (int j = 0; j < 4; j++) 2245 if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e)) 2246 return false; 2247 2248 return true; 2249 } 2250 2251 template <class T> 2252 bool 2253 Matrix44<T>::equalWithRelError (const Matrix44<T> &m, T e) const 2254 { 2255 for (int i = 0; i < 4; i++) 2256 for (int j = 0; j < 4; j++) 2257 if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e)) 2258 return false; 2259 2260 return true; 2261 } 2262 2263 template <class T> 2264 const Matrix44<T> & 2265 Matrix44<T>::operator += (const Matrix44<T> &v) 2266 { 2267 x[0][0] += v.x[0][0]; 2268 x[0][1] += v.x[0][1]; 2269 x[0][2] += v.x[0][2]; 2270 x[0][3] += v.x[0][3]; 2271 x[1][0] += v.x[1][0]; 2272 x[1][1] += v.x[1][1]; 2273 x[1][2] += v.x[1][2]; 2274 x[1][3] += v.x[1][3]; 2275 x[2][0] += v.x[2][0]; 2276 x[2][1] += v.x[2][1]; 2277 x[2][2] += v.x[2][2]; 2278 x[2][3] += v.x[2][3]; 2279 x[3][0] += v.x[3][0]; 2280 x[3][1] += v.x[3][1]; 2281 x[3][2] += v.x[3][2]; 2282 x[3][3] += v.x[3][3]; 2283 2284 return *this; 2285 } 2286 2287 template <class T> 2288 const Matrix44<T> & 2289 Matrix44<T>::operator += (T a) 2290 { 2291 x[0][0] += a; 2292 x[0][1] += a; 2293 x[0][2] += a; 2294 x[0][3] += a; 2295 x[1][0] += a; 2296 x[1][1] += a; 2297 x[1][2] += a; 2298 x[1][3] += a; 2299 x[2][0] += a; 2300 x[2][1] += a; 2301 x[2][2] += a; 2302 x[2][3] += a; 2303 x[3][0] += a; 2304 x[3][1] += a; 2305 x[3][2] += a; 2306 x[3][3] += a; 2307 2308 return *this; 2309 } 2310 2311 template <class T> 2312 Matrix44<T> 2313 Matrix44<T>::operator + (const Matrix44<T> &v) const 2314 { 2315 return Matrix44 (x[0][0] + v.x[0][0], 2316 x[0][1] + v.x[0][1], 2317 x[0][2] + v.x[0][2], 2318 x[0][3] + v.x[0][3], 2319 x[1][0] + v.x[1][0], 2320 x[1][1] + v.x[1][1], 2321 x[1][2] + v.x[1][2], 2322 x[1][3] + v.x[1][3], 2323 x[2][0] + v.x[2][0], 2324 x[2][1] + v.x[2][1], 2325 x[2][2] + v.x[2][2], 2326 x[2][3] + v.x[2][3], 2327 x[3][0] + v.x[3][0], 2328 x[3][1] + v.x[3][1], 2329 x[3][2] + v.x[3][2], 2330 x[3][3] + v.x[3][3]); 2331 } 2332 2333 template <class T> 2334 const Matrix44<T> & 2335 Matrix44<T>::operator -= (const Matrix44<T> &v) 2336 { 2337 x[0][0] -= v.x[0][0]; 2338 x[0][1] -= v.x[0][1]; 2339 x[0][2] -= v.x[0][2]; 2340 x[0][3] -= v.x[0][3]; 2341 x[1][0] -= v.x[1][0]; 2342 x[1][1] -= v.x[1][1]; 2343 x[1][2] -= v.x[1][2]; 2344 x[1][3] -= v.x[1][3]; 2345 x[2][0] -= v.x[2][0]; 2346 x[2][1] -= v.x[2][1]; 2347 x[2][2] -= v.x[2][2]; 2348 x[2][3] -= v.x[2][3]; 2349 x[3][0] -= v.x[3][0]; 2350 x[3][1] -= v.x[3][1]; 2351 x[3][2] -= v.x[3][2]; 2352 x[3][3] -= v.x[3][3]; 2353 2354 return *this; 2355 } 2356 2357 template <class T> 2358 const Matrix44<T> & 2359 Matrix44<T>::operator -= (T a) 2360 { 2361 x[0][0] -= a; 2362 x[0][1] -= a; 2363 x[0][2] -= a; 2364 x[0][3] -= a; 2365 x[1][0] -= a; 2366 x[1][1] -= a; 2367 x[1][2] -= a; 2368 x[1][3] -= a; 2369 x[2][0] -= a; 2370 x[2][1] -= a; 2371 x[2][2] -= a; 2372 x[2][3] -= a; 2373 x[3][0] -= a; 2374 x[3][1] -= a; 2375 x[3][2] -= a; 2376 x[3][3] -= a; 2377 2378 return *this; 2379 } 2380 2381 template <class T> 2382 Matrix44<T> 2383 Matrix44<T>::operator - (const Matrix44<T> &v) const 2384 { 2385 return Matrix44 (x[0][0] - v.x[0][0], 2386 x[0][1] - v.x[0][1], 2387 x[0][2] - v.x[0][2], 2388 x[0][3] - v.x[0][3], 2389 x[1][0] - v.x[1][0], 2390 x[1][1] - v.x[1][1], 2391 x[1][2] - v.x[1][2], 2392 x[1][3] - v.x[1][3], 2393 x[2][0] - v.x[2][0], 2394 x[2][1] - v.x[2][1], 2395 x[2][2] - v.x[2][2], 2396 x[2][3] - v.x[2][3], 2397 x[3][0] - v.x[3][0], 2398 x[3][1] - v.x[3][1], 2399 x[3][2] - v.x[3][2], 2400 x[3][3] - v.x[3][3]); 2401 } 2402 2403 template <class T> 2404 Matrix44<T> 2405 Matrix44<T>::operator - () const 2406 { 2407 return Matrix44 (-x[0][0], 2408 -x[0][1], 2409 -x[0][2], 2410 -x[0][3], 2411 -x[1][0], 2412 -x[1][1], 2413 -x[1][2], 2414 -x[1][3], 2415 -x[2][0], 2416 -x[2][1], 2417 -x[2][2], 2418 -x[2][3], 2419 -x[3][0], 2420 -x[3][1], 2421 -x[3][2], 2422 -x[3][3]); 2423 } 2424 2425 template <class T> 2426 const Matrix44<T> & 2427 Matrix44<T>::negate () 2428 { 2429 x[0][0] = -x[0][0]; 2430 x[0][1] = -x[0][1]; 2431 x[0][2] = -x[0][2]; 2432 x[0][3] = -x[0][3]; 2433 x[1][0] = -x[1][0]; 2434 x[1][1] = -x[1][1]; 2435 x[1][2] = -x[1][2]; 2436 x[1][3] = -x[1][3]; 2437 x[2][0] = -x[2][0]; 2438 x[2][1] = -x[2][1]; 2439 x[2][2] = -x[2][2]; 2440 x[2][3] = -x[2][3]; 2441 x[3][0] = -x[3][0]; 2442 x[3][1] = -x[3][1]; 2443 x[3][2] = -x[3][2]; 2444 x[3][3] = -x[3][3]; 2445 2446 return *this; 2447 } 2448 2449 template <class T> 2450 const Matrix44<T> & 2451 Matrix44<T>::operator *= (T a) 2452 { 2453 x[0][0] *= a; 2454 x[0][1] *= a; 2455 x[0][2] *= a; 2456 x[0][3] *= a; 2457 x[1][0] *= a; 2458 x[1][1] *= a; 2459 x[1][2] *= a; 2460 x[1][3] *= a; 2461 x[2][0] *= a; 2462 x[2][1] *= a; 2463 x[2][2] *= a; 2464 x[2][3] *= a; 2465 x[3][0] *= a; 2466 x[3][1] *= a; 2467 x[3][2] *= a; 2468 x[3][3] *= a; 2469 2470 return *this; 2471 } 2472 2473 template <class T> 2474 Matrix44<T> 2475 Matrix44<T>::operator * (T a) const 2476 { 2477 return Matrix44 (x[0][0] * a, 2478 x[0][1] * a, 2479 x[0][2] * a, 2480 x[0][3] * a, 2481 x[1][0] * a, 2482 x[1][1] * a, 2483 x[1][2] * a, 2484 x[1][3] * a, 2485 x[2][0] * a, 2486 x[2][1] * a, 2487 x[2][2] * a, 2488 x[2][3] * a, 2489 x[3][0] * a, 2490 x[3][1] * a, 2491 x[3][2] * a, 2492 x[3][3] * a); 2493 } 2494 2495 template <class T> 2496 inline Matrix44<T> 2497 operator * (T a, const Matrix44<T> &v) 2498 { 2499 return v * a; 2500 } 2501 2502 template <class T> 2503 inline const Matrix44<T> & 2504 Matrix44<T>::operator *= (const Matrix44<T> &v) 2505 { 2506 Matrix44 tmp (T (0)); 2507 2508 multiply (*this, v, tmp); 2509 *this = tmp; 2510 return *this; 2511 } 2512 2513 template <class T> 2514 inline Matrix44<T> 2515 Matrix44<T>::operator * (const Matrix44<T> &v) const 2516 { 2517 Matrix44 tmp (T (0)); 2518 2519 multiply (*this, v, tmp); 2520 return tmp; 2521 } 2522 2523 template <class T> 2524 void 2525 Matrix44<T>::multiply (const Matrix44<T> &a, 2526 const Matrix44<T> &b, 2527 Matrix44<T> &c) 2528 { 2529 register const T * IMATH_RESTRICT ap = &a.x[0][0]; 2530 register const T * IMATH_RESTRICT bp = &b.x[0][0]; 2531 register T * IMATH_RESTRICT cp = &c.x[0][0]; 2532 2533 register T a0, a1, a2, a3; 2534 2535 a0 = ap[0]; 2536 a1 = ap[1]; 2537 a2 = ap[2]; 2538 a3 = ap[3]; 2539 2540 cp[0] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12]; 2541 cp[1] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13]; 2542 cp[2] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14]; 2543 cp[3] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15]; 2544 2545 a0 = ap[4]; 2546 a1 = ap[5]; 2547 a2 = ap[6]; 2548 a3 = ap[7]; 2549 2550 cp[4] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12]; 2551 cp[5] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13]; 2552 cp[6] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14]; 2553 cp[7] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15]; 2554 2555 a0 = ap[8]; 2556 a1 = ap[9]; 2557 a2 = ap[10]; 2558 a3 = ap[11]; 2559 2560 cp[8] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12]; 2561 cp[9] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13]; 2562 cp[10] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14]; 2563 cp[11] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15]; 2564 2565 a0 = ap[12]; 2566 a1 = ap[13]; 2567 a2 = ap[14]; 2568 a3 = ap[15]; 2569 2570 cp[12] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12]; 2571 cp[13] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13]; 2572 cp[14] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14]; 2573 cp[15] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15]; 2574 } 2575 2576 template <class T> template <class S> 2577 void 2578 Matrix44<T>::multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const 2579 { 2580 S a, b, c, w; 2581 2582 a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0] + x[3][0]; 2583 b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1] + x[3][1]; 2584 c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2] + x[3][2]; 2585 w = src[0] * x[0][3] + src[1] * x[1][3] + src[2] * x[2][3] + x[3][3]; 2586 2587 dst.x = a / w; 2588 dst.y = b / w; 2589 dst.z = c / w; 2590 } 2591 2592 template <class T> template <class S> 2593 void 2594 Matrix44<T>::multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const 2595 { 2596 S a, b, c; 2597 2598 a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0]; 2599 b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1]; 2600 c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2]; 2601 2602 dst.x = a; 2603 dst.y = b; 2604 dst.z = c; 2605 } 2606 2607 template <class T> 2608 const Matrix44<T> & 2609 Matrix44<T>::operator /= (T a) 2610 { 2611 x[0][0] /= a; 2612 x[0][1] /= a; 2613 x[0][2] /= a; 2614 x[0][3] /= a; 2615 x[1][0] /= a; 2616 x[1][1] /= a; 2617 x[1][2] /= a; 2618 x[1][3] /= a; 2619 x[2][0] /= a; 2620 x[2][1] /= a; 2621 x[2][2] /= a; 2622 x[2][3] /= a; 2623 x[3][0] /= a; 2624 x[3][1] /= a; 2625 x[3][2] /= a; 2626 x[3][3] /= a; 2627 2628 return *this; 2629 } 2630 2631 template <class T> 2632 Matrix44<T> 2633 Matrix44<T>::operator / (T a) const 2634 { 2635 return Matrix44 (x[0][0] / a, 2636 x[0][1] / a, 2637 x[0][2] / a, 2638 x[0][3] / a, 2639 x[1][0] / a, 2640 x[1][1] / a, 2641 x[1][2] / a, 2642 x[1][3] / a, 2643 x[2][0] / a, 2644 x[2][1] / a, 2645 x[2][2] / a, 2646 x[2][3] / a, 2647 x[3][0] / a, 2648 x[3][1] / a, 2649 x[3][2] / a, 2650 x[3][3] / a); 2651 } 2652 2653 template <class T> 2654 const Matrix44<T> & 2655 Matrix44<T>::transpose () 2656 { 2657 Matrix44 tmp (x[0][0], 2658 x[1][0], 2659 x[2][0], 2660 x[3][0], 2661 x[0][1], 2662 x[1][1], 2663 x[2][1], 2664 x[3][1], 2665 x[0][2], 2666 x[1][2], 2667 x[2][2], 2668 x[3][2], 2669 x[0][3], 2670 x[1][3], 2671 x[2][3], 2672 x[3][3]); 2673 *this = tmp; 2674 return *this; 2675 } 2676 2677 template <class T> 2678 Matrix44<T> 2679 Matrix44<T>::transposed () const 2680 { 2681 return Matrix44 (x[0][0], 2682 x[1][0], 2683 x[2][0], 2684 x[3][0], 2685 x[0][1], 2686 x[1][1], 2687 x[2][1], 2688 x[3][1], 2689 x[0][2], 2690 x[1][2], 2691 x[2][2], 2692 x[3][2], 2693 x[0][3], 2694 x[1][3], 2695 x[2][3], 2696 x[3][3]); 2697 } 2698 2699 template <class T> 2700 const Matrix44<T> & 2701 Matrix44<T>::gjInvert (bool singExc) throw (Iex::MathExc) 2702 { 2703 *this = gjInverse (singExc); 2704 return *this; 2705 } 2706 2707 template <class T> 2708 Matrix44<T> 2709 Matrix44<T>::gjInverse (bool singExc) const throw (Iex::MathExc) 2710 { 2711 int i, j, k; 2712 Matrix44 s; 2713 Matrix44 t (*this); 2714 2715 // Forward elimination 2716 2717 for (i = 0; i < 3 ; i++) 2718 { 2719 int pivot = i; 2720 2721 T pivotsize = t[i][i]; 2722 2723 if (pivotsize < 0) 2724 pivotsize = -pivotsize; 2725 2726 for (j = i + 1; j < 4; j++) 2727 { 2728 T tmp = t[j][i]; 2729 2730 if (tmp < 0) 2731 tmp = -tmp; 2732 2733 if (tmp > pivotsize) 2734 { 2735 pivot = j; 2736 pivotsize = tmp; 2737 } 2738 } 2739 2740 if (pivotsize == 0) 2741 { 2742 if (singExc) 2743 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix."); 2744 2745 return Matrix44(); 2746 } 2747 2748 if (pivot != i) 2749 { 2750 for (j = 0; j < 4; j++) 2751 { 2752 T tmp; 2753 2754 tmp = t[i][j]; 2755 t[i][j] = t[pivot][j]; 2756 t[pivot][j] = tmp; 2757 2758 tmp = s[i][j]; 2759 s[i][j] = s[pivot][j]; 2760 s[pivot][j] = tmp; 2761 } 2762 } 2763 2764 for (j = i + 1; j < 4; j++) 2765 { 2766 T f = t[j][i] / t[i][i]; 2767 2768 for (k = 0; k < 4; k++) 2769 { 2770 t[j][k] -= f * t[i][k]; 2771 s[j][k] -= f * s[i][k]; 2772 } 2773 } 2774 } 2775 2776 // Backward substitution 2777 2778 for (i = 3; i >= 0; --i) 2779 { 2780 T f; 2781 2782 if ((f = t[i][i]) == 0) 2783 { 2784 if (singExc) 2785 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix."); 2786 2787 return Matrix44(); 2788 } 2789 2790 for (j = 0; j < 4; j++) 2791 { 2792 t[i][j] /= f; 2793 s[i][j] /= f; 2794 } 2795 2796 for (j = 0; j < i; j++) 2797 { 2798 f = t[j][i]; 2799 2800 for (k = 0; k < 4; k++) 2801 { 2802 t[j][k] -= f * t[i][k]; 2803 s[j][k] -= f * s[i][k]; 2804 } 2805 } 2806 } 2807 2808 return s; 2809 } 2810 2811 template <class T> 2812 const Matrix44<T> & 2813 Matrix44<T>::invert (bool singExc) throw (Iex::MathExc) 2814 { 2815 *this = inverse (singExc); 2816 return *this; 2817 } 2818 2819 template <class T> 2820 Matrix44<T> 2821 Matrix44<T>::inverse (bool singExc) const throw (Iex::MathExc) 2822 { 2823 if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1) 2824 return gjInverse(singExc); 2825 2826 Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2], 2827 x[2][1] * x[0][2] - x[0][1] * x[2][2], 2828 x[0][1] * x[1][2] - x[1][1] * x[0][2], 2829 0, 2830 2831 x[2][0] * x[1][2] - x[1][0] * x[2][2], 2832 x[0][0] * x[2][2] - x[2][0] * x[0][2], 2833 x[1][0] * x[0][2] - x[0][0] * x[1][2], 2834 0, 2835 2836 x[1][0] * x[2][1] - x[2][0] * x[1][1], 2837 x[2][0] * x[0][1] - x[0][0] * x[2][1], 2838 x[0][0] * x[1][1] - x[1][0] * x[0][1], 2839 0, 2840 2841 0, 2842 0, 2843 0, 2844 1); 2845 2846 T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0]; 2847 2848 if (Imath::abs (r) >= 1) 2849 { 2850 for (int i = 0; i < 3; ++i) 2851 { 2852 for (int j = 0; j < 3; ++j) 2853 { 2854 s[i][j] /= r; 2855 } 2856 } 2857 } 2858 else 2859 { 2860 T mr = Imath::abs (r) / limits<T>::smallest(); 2861 2862 for (int i = 0; i < 3; ++i) 2863 { 2864 for (int j = 0; j < 3; ++j) 2865 { 2866 if (mr > Imath::abs (s[i][j])) 2867 { 2868 s[i][j] /= r; 2869 } 2870 else 2871 { 2872 if (singExc) 2873 throw SingMatrixExc ("Cannot invert singular matrix."); 2874 2875 return Matrix44(); 2876 } 2877 } 2878 } 2879 } 2880 2881 s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0]; 2882 s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1]; 2883 s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2]; 2884 2885 return s; 2886 } 2887 2888 template <class T> 2889 inline T 2890 Matrix44<T>::fastMinor( const int r0, const int r1, const int r2, 2891 const int c0, const int c1, const int c2) const 2892 { 2893 return x[r0][c0] * (x[r1][c1]*x[r2][c2] - x[r1][c2]*x[r2][c1]) 2894 + x[r0][c1] * (x[r1][c2]*x[r2][c0] - x[r1][c0]*x[r2][c2]) 2895 + x[r0][c2] * (x[r1][c0]*x[r2][c1] - x[r1][c1]*x[r2][c0]); 2896 } 2897 2898 template <class T> 2899 inline T 2900 Matrix44<T>::minorOf (const int r, const int c) const 2901 { 2902 int r0 = 0 + (r < 1 ? 1 : 0); 2903 int r1 = 1 + (r < 2 ? 1 : 0); 2904 int r2 = 2 + (r < 3 ? 1 : 0); 2905 int c0 = 0 + (c < 1 ? 1 : 0); 2906 int c1 = 1 + (c < 2 ? 1 : 0); 2907 int c2 = 2 + (c < 3 ? 1 : 0); 2908 2909 Matrix33<T> working (x[r0][c0],x[r1][c0],x[r2][c0], 2910 x[r0][c1],x[r1][c1],x[r2][c1], 2911 x[r0][c2],x[r1][c2],x[r2][c2]); 2912 2913 return working.determinant(); 2914 } 2915 2916 template <class T> 2917 inline T 2918 Matrix44<T>::determinant () const 2919 { 2920 T sum = (T)0; 2921 2922 if (x[0][3] != 0.) sum -= x[0][3] * fastMinor(1,2,3,0,1,2); 2923 if (x[1][3] != 0.) sum += x[1][3] * fastMinor(0,2,3,0,1,2); 2924 if (x[2][3] != 0.) sum -= x[2][3] * fastMinor(0,1,3,0,1,2); 2925 if (x[3][3] != 0.) sum += x[3][3] * fastMinor(0,1,2,0,1,2); 2926 2927 return sum; 2928 } 2929 2930 template <class T> 2931 template <class S> 2932 const Matrix44<T> & 2933 Matrix44<T>::setEulerAngles (const Vec3<S>& r) 2934 { 2935 S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx; 2936 2937 cos_rz = Math<T>::cos (r[2]); 2938 cos_ry = Math<T>::cos (r[1]); 2939 cos_rx = Math<T>::cos (r[0]); 2940 2941 sin_rz = Math<T>::sin (r[2]); 2942 sin_ry = Math<T>::sin (r[1]); 2943 sin_rx = Math<T>::sin (r[0]); 2944 2945 x[0][0] = cos_rz * cos_ry; 2946 x[0][1] = sin_rz * cos_ry; 2947 x[0][2] = -sin_ry; 2948 x[0][3] = 0; 2949 2950 x[1][0] = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx; 2951 x[1][1] = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx; 2952 x[1][2] = cos_ry * sin_rx; 2953 x[1][3] = 0; 2954 2955 x[2][0] = sin_rz * sin_rx + cos_rz * sin_ry * cos_rx; 2956 x[2][1] = -cos_rz * sin_rx + sin_rz * sin_ry * cos_rx; 2957 x[2][2] = cos_ry * cos_rx; 2958 x[2][3] = 0; 2959 2960 x[3][0] = 0; 2961 x[3][1] = 0; 2962 x[3][2] = 0; 2963 x[3][3] = 1; 2964 2965 return *this; 2966 } 2967 2968 template <class T> 2969 template <class S> 2970 const Matrix44<T> & 2971 Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle) 2972 { 2973 Vec3<S> unit (axis.normalized()); 2974 S sine = Math<T>::sin (angle); 2975 S cosine = Math<T>::cos (angle); 2976 2977 x[0][0] = unit[0] * unit[0] * (1 - cosine) + cosine; 2978 x[0][1] = unit[0] * unit[1] * (1 - cosine) + unit[2] * sine; 2979 x[0][2] = unit[0] * unit[2] * (1 - cosine) - unit[1] * sine; 2980 x[0][3] = 0; 2981 2982 x[1][0] = unit[0] * unit[1] * (1 - cosine) - unit[2] * sine; 2983 x[1][1] = unit[1] * unit[1] * (1 - cosine) + cosine; 2984 x[1][2] = unit[1] * unit[2] * (1 - cosine) + unit[0] * sine; 2985 x[1][3] = 0; 2986 2987 x[2][0] = unit[0] * unit[2] * (1 - cosine) + unit[1] * sine; 2988 x[2][1] = unit[1] * unit[2] * (1 - cosine) - unit[0] * sine; 2989 x[2][2] = unit[2] * unit[2] * (1 - cosine) + cosine; 2990 x[2][3] = 0; 2991 2992 x[3][0] = 0; 2993 x[3][1] = 0; 2994 x[3][2] = 0; 2995 x[3][3] = 1; 2996 2997 return *this; 2998 } 2999 3000 template <class T> 3001 template <class S> 3002 const Matrix44<T> & 3003 Matrix44<T>::rotate (const Vec3<S> &r) 3004 { 3005 S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx; 3006 S m00, m01, m02; 3007 S m10, m11, m12; 3008 S m20, m21, m22; 3009 3010 cos_rz = Math<S>::cos (r[2]); 3011 cos_ry = Math<S>::cos (r[1]); 3012 cos_rx = Math<S>::cos (r[0]); 3013 3014 sin_rz = Math<S>::sin (r[2]); 3015 sin_ry = Math<S>::sin (r[1]); 3016 sin_rx = Math<S>::sin (r[0]); 3017 3018 m00 = cos_rz * cos_ry; 3019 m01 = sin_rz * cos_ry; 3020 m02 = -sin_ry; 3021 m10 = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx; 3022 m11 = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx; 3023 m12 = cos_ry * sin_rx; 3024 m20 = -sin_rz * -sin_rx + cos_rz * sin_ry * cos_rx; 3025 m21 = cos_rz * -sin_rx + sin_rz * sin_ry * cos_rx; 3026 m22 = cos_ry * cos_rx; 3027 3028 Matrix44<T> P (*this); 3029 3030 x[0][0] = P[0][0] * m00 + P[1][0] * m01 + P[2][0] * m02; 3031 x[0][1] = P[0][1] * m00 + P[1][1] * m01 + P[2][1] * m02; 3032 x[0][2] = P[0][2] * m00 + P[1][2] * m01 + P[2][2] * m02; 3033 x[0][3] = P[0][3] * m00 + P[1][3] * m01 + P[2][3] * m02; 3034 3035 x[1][0] = P[0][0] * m10 + P[1][0] * m11 + P[2][0] * m12; 3036 x[1][1] = P[0][1] * m10 + P[1][1] * m11 + P[2][1] * m12; 3037 x[1][2] = P[0][2] * m10 + P[1][2] * m11 + P[2][2] * m12; 3038 x[1][3] = P[0][3] * m10 + P[1][3] * m11 + P[2][3] * m12; 3039 3040 x[2][0] = P[0][0] * m20 + P[1][0] * m21 + P[2][0] * m22; 3041 x[2][1] = P[0][1] * m20 + P[1][1] * m21 + P[2][1] * m22; 3042 x[2][2] = P[0][2] * m20 + P[1][2] * m21 + P[2][2] * m22; 3043 x[2][3] = P[0][3] * m20 + P[1][3] * m21 + P[2][3] * m22; 3044 3045 return *this; 3046 } 3047 3048 template <class T> 3049 const Matrix44<T> & 3050 Matrix44<T>::setScale (T s) 3051 { 3052 memset (x, 0, sizeof (x)); 3053 x[0][0] = s; 3054 x[1][1] = s; 3055 x[2][2] = s; 3056 x[3][3] = 1; 3057 3058 return *this; 3059 } 3060 3061 template <class T> 3062 template <class S> 3063 const Matrix44<T> & 3064 Matrix44<T>::setScale (const Vec3<S> &s) 3065 { 3066 memset (x, 0, sizeof (x)); 3067 x[0][0] = s[0]; 3068 x[1][1] = s[1]; 3069 x[2][2] = s[2]; 3070 x[3][3] = 1; 3071 3072 return *this; 3073 } 3074 3075 template <class T> 3076 template <class S> 3077 const Matrix44<T> & 3078 Matrix44<T>::scale (const Vec3<S> &s) 3079 { 3080 x[0][0] *= s[0]; 3081 x[0][1] *= s[0]; 3082 x[0][2] *= s[0]; 3083 x[0][3] *= s[0]; 3084 3085 x[1][0] *= s[1]; 3086 x[1][1] *= s[1]; 3087 x[1][2] *= s[1]; 3088 x[1][3] *= s[1]; 3089 3090 x[2][0] *= s[2]; 3091 x[2][1] *= s[2]; 3092 x[2][2] *= s[2]; 3093 x[2][3] *= s[2]; 3094 3095 return *this; 3096 } 3097 3098 template <class T> 3099 template <class S> 3100 const Matrix44<T> & 3101 Matrix44<T>::setTranslation (const Vec3<S> &t) 3102 { 3103 x[0][0] = 1; 3104 x[0][1] = 0; 3105 x[0][2] = 0; 3106 x[0][3] = 0; 3107 3108 x[1][0] = 0; 3109 x[1][1] = 1; 3110 x[1][2] = 0; 3111 x[1][3] = 0; 3112 3113 x[2][0] = 0; 3114 x[2][1] = 0; 3115 x[2][2] = 1; 3116 x[2][3] = 0; 3117 3118 x[3][0] = t[0]; 3119 x[3][1] = t[1]; 3120 x[3][2] = t[2]; 3121 x[3][3] = 1; 3122 3123 return *this; 3124 } 3125 3126 template <class T> 3127 inline const Vec3<T> 3128 Matrix44<T>::translation () const 3129 { 3130 return Vec3<T> (x[3][0], x[3][1], x[3][2]); 3131 } 3132 3133 template <class T> 3134 template <class S> 3135 const Matrix44<T> & 3136 Matrix44<T>::translate (const Vec3<S> &t) 3137 { 3138 x[3][0] += t[0] * x[0][0] + t[1] * x[1][0] + t[2] * x[2][0]; 3139 x[3][1] += t[0] * x[0][1] + t[1] * x[1][1] + t[2] * x[2][1]; 3140 x[3][2] += t[0] * x[0][2] + t[1] * x[1][2] + t[2] * x[2][2]; 3141 x[3][3] += t[0] * x[0][3] + t[1] * x[1][3] + t[2] * x[2][3]; 3142 3143 return *this; 3144 } 3145 3146 template <class T> 3147 template <class S> 3148 const Matrix44<T> & 3149 Matrix44<T>::setShear (const Vec3<S> &h) 3150 { 3151 x[0][0] = 1; 3152 x[0][1] = 0; 3153 x[0][2] = 0; 3154 x[0][3] = 0; 3155 3156 x[1][0] = h[0]; 3157 x[1][1] = 1; 3158 x[1][2] = 0; 3159 x[1][3] = 0; 3160 3161 x[2][0] = h[1]; 3162 x[2][1] = h[2]; 3163 x[2][2] = 1; 3164 x[2][3] = 0; 3165 3166 x[3][0] = 0; 3167 x[3][1] = 0; 3168 x[3][2] = 0; 3169 x[3][3] = 1; 3170 3171 return *this; 3172 } 3173 3174 template <class T> 3175 template <class S> 3176 const Matrix44<T> & 3177 Matrix44<T>::setShear (const Shear6<S> &h) 3178 { 3179 x[0][0] = 1; 3180 x[0][1] = h.yx; 3181 x[0][2] = h.zx; 3182 x[0][3] = 0; 3183 3184 x[1][0] = h.xy; 3185 x[1][1] = 1; 3186 x[1][2] = h.zy; 3187 x[1][3] = 0; 3188 3189 x[2][0] = h.xz; 3190 x[2][1] = h.yz; 3191 x[2][2] = 1; 3192 x[2][3] = 0; 3193 3194 x[3][0] = 0; 3195 x[3][1] = 0; 3196 x[3][2] = 0; 3197 x[3][3] = 1; 3198 3199 return *this; 3200 } 3201 3202 template <class T> 3203 template <class S> 3204 const Matrix44<T> & 3205 Matrix44<T>::shear (const Vec3<S> &h) 3206 { 3207 // 3208 // In this case, we don't need a temp. copy of the matrix 3209 // because we never use a value on the RHS after we've 3210 // changed it on the LHS. 3211 // 3212 3213 for (int i=0; i < 4; i++) 3214 { 3215 x[2][i] += h[1] * x[0][i] + h[2] * x[1][i]; 3216 x[1][i] += h[0] * x[0][i]; 3217 } 3218 3219 return *this; 3220 } 3221 3222 template <class T> 3223 template <class S> 3224 const Matrix44<T> & 3225 Matrix44<T>::shear (const Shear6<S> &h) 3226 { 3227 Matrix44<T> P (*this); 3228 3229 for (int i=0; i < 4; i++) 3230 { 3231 x[0][i] = P[0][i] + h.yx * P[1][i] + h.zx * P[2][i]; 3232 x[1][i] = h.xy * P[0][i] + P[1][i] + h.zy * P[2][i]; 3233 x[2][i] = h.xz * P[0][i] + h.yz * P[1][i] + P[2][i]; 3234 } 3235 3236 return *this; 3237 } 3238 3239 3240 //-------------------------------- 3241 // Implementation of stream output 3242 //-------------------------------- 3243 3244 template <class T> 3245 std::ostream & 3246 operator << (std::ostream &s, const Matrix33<T> &m) 3247 { 3248 std::ios_base::fmtflags oldFlags = s.flags(); 3249 int width; 3250 3251 if (s.flags() & std::ios_base::fixed) 3252 { 3253 s.setf (std::ios_base::showpoint); 3254 width = s.precision() + 5; 3255 } 3256 else 3257 { 3258 s.setf (std::ios_base::scientific); 3259 s.setf (std::ios_base::showpoint); 3260 width = s.precision() + 8; 3261 } 3262 3263 s << "(" << std::setw (width) << m[0][0] << 3264 " " << std::setw (width) << m[0][1] << 3265 " " << std::setw (width) << m[0][2] << "\n" << 3266 3267 " " << std::setw (width) << m[1][0] << 3268 " " << std::setw (width) << m[1][1] << 3269 " " << std::setw (width) << m[1][2] << "\n" << 3270 3271 " " << std::setw (width) << m[2][0] << 3272 " " << std::setw (width) << m[2][1] << 3273 " " << std::setw (width) << m[2][2] << ")\n"; 3274 3275 s.flags (oldFlags); 3276 return s; 3277 } 3278 3279 template <class T> 3280 std::ostream & 3281 operator << (std::ostream &s, const Matrix44<T> &m) 3282 { 3283 std::ios_base::fmtflags oldFlags = s.flags(); 3284 int width; 3285 3286 if (s.flags() & std::ios_base::fixed) 3287 { 3288 s.setf (std::ios_base::showpoint); 3289 width = s.precision() + 5; 3290 } 3291 else 3292 { 3293 s.setf (std::ios_base::scientific); 3294 s.setf (std::ios_base::showpoint); 3295 width = s.precision() + 8; 3296 } 3297 3298 s << "(" << std::setw (width) << m[0][0] << 3299 " " << std::setw (width) << m[0][1] << 3300 " " << std::setw (width) << m[0][2] << 3301 " " << std::setw (width) << m[0][3] << "\n" << 3302 3303 " " << std::setw (width) << m[1][0] << 3304 " " << std::setw (width) << m[1][1] << 3305 " " << std::setw (width) << m[1][2] << 3306 " " << std::setw (width) << m[1][3] << "\n" << 3307 3308 " " << std::setw (width) << m[2][0] << 3309 " " << std::setw (width) << m[2][1] << 3310 " " << std::setw (width) << m[2][2] << 3311 " " << std::setw (width) << m[2][3] << "\n" << 3312 3313 " " << std::setw (width) << m[3][0] << 3314 " " << std::setw (width) << m[3][1] << 3315 " " << std::setw (width) << m[3][2] << 3316 " " << std::setw (width) << m[3][3] << ")\n"; 3317 3318 s.flags (oldFlags); 3319 return s; 3320 } 3321 3322 3323 //--------------------------------------------------------------- 3324 // Implementation of vector-times-matrix multiplication operators 3325 //--------------------------------------------------------------- 3326 3327 template <class S, class T> 3328 inline const Vec2<S> & 3329 operator *= (Vec2<S> &v, const Matrix33<T> &m) 3330 { 3331 S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]); 3332 S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]); 3333 S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]); 3334 3335 v.x = x / w; 3336 v.y = y / w; 3337 3338 return v; 3339 } 3340 3341 template <class S, class T> 3342 inline Vec2<S> 3343 operator * (const Vec2<S> &v, const Matrix33<T> &m) 3344 { 3345 S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]); 3346 S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]); 3347 S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]); 3348 3349 return Vec2<S> (x / w, y / w); 3350 } 3351 3352 3353 template <class S, class T> 3354 inline const Vec3<S> & 3355 operator *= (Vec3<S> &v, const Matrix33<T> &m) 3356 { 3357 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]); 3358 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]); 3359 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]); 3360 3361 v.x = x; 3362 v.y = y; 3363 v.z = z; 3364 3365 return v; 3366 } 3367 3368 template <class S, class T> 3369 inline Vec3<S> 3370 operator * (const Vec3<S> &v, const Matrix33<T> &m) 3371 { 3372 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]); 3373 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]); 3374 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]); 3375 3376 return Vec3<S> (x, y, z); 3377 } 3378 3379 3380 template <class S, class T> 3381 inline const Vec3<S> & 3382 operator *= (Vec3<S> &v, const Matrix44<T> &m) 3383 { 3384 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]); 3385 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]); 3386 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]); 3387 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]); 3388 3389 v.x = x / w; 3390 v.y = y / w; 3391 v.z = z / w; 3392 3393 return v; 3394 } 3395 3396 template <class S, class T> 3397 inline Vec3<S> 3398 operator * (const Vec3<S> &v, const Matrix44<T> &m) 3399 { 3400 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]); 3401 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]); 3402 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]); 3403 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]); 3404 3405 return Vec3<S> (x / w, y / w, z / w); 3406 } 3407 3408 3409 template <class S, class T> 3410 inline const Vec4<S> & 3411 operator *= (Vec4<S> &v, const Matrix44<T> &m) 3412 { 3413 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + v.w * m[3][0]); 3414 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + v.w * m[3][1]); 3415 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + v.w * m[3][2]); 3416 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + v.w * m[3][3]); 3417 3418 v.x = x; 3419 v.y = y; 3420 v.z = z; 3421 v.w = w; 3422 3423 return v; 3424 } 3425 3426 template <class S, class T> 3427 inline Vec4<S> 3428 operator * (const Vec4<S> &v, const Matrix44<T> &m) 3429 { 3430 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + v.w * m[3][0]); 3431 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + v.w * m[3][1]); 3432 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + v.w * m[3][2]); 3433 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + v.w * m[3][3]); 3434 3435 return Vec4<S> (x, y, z, w); 3436 } 3437 3438 } // namespace Imath 3439 3440 3441 3442 #endif 3443
