1 ///////////////////////////////////////////////////////////////////////////
3 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
6 // All rights reserved.
8 // Redistribution and use in source and binary forms, with or without
9 // modification, are permitted provided that the following conditions are
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
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.
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.
33 ///////////////////////////////////////////////////////////////////////////
37 #ifndef INCLUDED_IMATHMATRIX_H
38 #define INCLUDED_IMATHMATRIX_H
40 //----------------------------------------------------------------
42 // 2D (3x3) and 3D (4x4) transformation matrix templates.
44 //----------------------------------------------------------------
46 #include "ImathPlatform.h"
50 #include "ImathShear.h"
59 template <class T> class Matrix33
69 T * operator [] (int i);
70 const T * operator [] (int i) const;
87 Matrix33 (const T a[3][3]);
88 // a[0][0] a[0][1] a[0][2]
89 // a[1][0] a[1][1] a[1][2]
90 // a[2][0] a[2][1] a[2][2]
92 Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i);
99 //--------------------------------
100 // Copy constructor and assignment
101 //--------------------------------
103 Matrix33 (const Matrix33 &v);
105 const Matrix33 & operator = (const Matrix33 &v);
106 const Matrix33 & operator = (T a);
109 //----------------------
110 // Compatibility with Sb
111 //----------------------
114 const T * getValue () const;
117 void getValue (Matrix33<S> &v) const;
119 Matrix33 & setValue (const Matrix33<S> &v);
122 Matrix33 & setTheMatrix (const Matrix33<S> &v);
136 bool operator == (const Matrix33 &v) const;
137 bool operator != (const Matrix33 &v) const;
139 //-----------------------------------------------------------------------
140 // Compare two matrices and test if they are "approximately equal":
142 // equalWithAbsError (m, e)
144 // Returns true if the coefficients of this and m are the same with
145 // an absolute error of no more than e, i.e., for all i, j
147 // abs (this[i][j] - m[i][j]) <= e
149 // equalWithRelError (m, e)
151 // Returns true if the coefficients of this and m are the same with
152 // a relative error of no more than e, i.e., for all i, j
154 // abs (this[i] - v[i][j]) <= e * abs (this[i][j])
155 //-----------------------------------------------------------------------
157 bool equalWithAbsError (const Matrix33<T> &v, T e) const;
158 bool equalWithRelError (const Matrix33<T> &v, T e) const;
161 //------------------------
162 // Component-wise addition
163 //------------------------
165 const Matrix33 & operator += (const Matrix33 &v);
166 const Matrix33 & operator += (T a);
167 Matrix33 operator + (const Matrix33 &v) const;
170 //---------------------------
171 // Component-wise subtraction
172 //---------------------------
174 const Matrix33 & operator -= (const Matrix33 &v);
175 const Matrix33 & operator -= (T a);
176 Matrix33 operator - (const Matrix33 &v) const;
179 //------------------------------------
180 // Component-wise multiplication by -1
181 //------------------------------------
183 Matrix33 operator - () const;
184 const Matrix33 & negate ();
187 //------------------------------
188 // Component-wise multiplication
189 //------------------------------
191 const Matrix33 & operator *= (T a);
192 Matrix33 operator * (T a) const;
195 //-----------------------------------
196 // Matrix-times-matrix multiplication
197 //-----------------------------------
199 const Matrix33 & operator *= (const Matrix33 &v);
200 Matrix33 operator * (const Matrix33 &v) const;
203 //---------------------------------------------
204 // Vector-times-matrix multiplication; see also
205 // the "operator *" functions defined below.
206 //---------------------------------------------
209 void multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const;
212 void multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const;
215 //------------------------
216 // Component-wise division
217 //------------------------
219 const Matrix33 & operator /= (T a);
220 Matrix33 operator / (T a) const;
227 const Matrix33 & transpose ();
228 Matrix33 transposed () const;
231 //------------------------------------------------------------
232 // Inverse matrix: If singExc is false, inverting a singular
233 // matrix produces an identity matrix. If singExc is true,
234 // inverting a singular matrix throws a SingMatrixExc.
236 // inverse() and invert() invert matrices using determinants;
237 // gjInverse() and gjInvert() use the Gauss-Jordan method.
239 // inverse() and invert() are significantly faster than
240 // gjInverse() and gjInvert(), but the results may be slightly
243 //------------------------------------------------------------
245 const Matrix33 & invert (bool singExc = false)
246 throw (Iex::MathExc);
248 Matrix33<T> inverse (bool singExc = false) const
249 throw (Iex::MathExc);
251 const Matrix33 & gjInvert (bool singExc = false)
252 throw (Iex::MathExc);
254 Matrix33<T> gjInverse (bool singExc = false) const
255 throw (Iex::MathExc);
258 //-----------------------------------------
259 // Set matrix to rotation by r (in radians)
260 //-----------------------------------------
263 const Matrix33 & setRotation (S r);
266 //-----------------------------
267 // Rotate the given matrix by r
268 //-----------------------------
271 const Matrix33 & rotate (S r);
274 //--------------------------------------------
275 // Set matrix to scale by given uniform factor
276 //--------------------------------------------
278 const Matrix33 & setScale (T s);
281 //------------------------------------
282 // Set matrix to scale by given vector
283 //------------------------------------
286 const Matrix33 & setScale (const Vec2<S> &s);
289 //----------------------
290 // Scale the matrix by s
291 //----------------------
294 const Matrix33 & scale (const Vec2<S> &s);
297 //------------------------------------------
298 // Set matrix to translation by given vector
299 //------------------------------------------
302 const Matrix33 & setTranslation (const Vec2<S> &t);
305 //-----------------------------
306 // Return translation component
307 //-----------------------------
309 Vec2<T> translation () const;
312 //--------------------------
313 // Translate the matrix by t
314 //--------------------------
317 const Matrix33 & translate (const Vec2<S> &t);
320 //-----------------------------------------------------------
321 // Set matrix to shear x for each y coord. by given factor xy
322 //-----------------------------------------------------------
325 const Matrix33 & setShear (const S &h);
328 //-------------------------------------------------------------
329 // Set matrix to shear x for each y coord. by given factor h[0]
330 // and to shear y for each x coord. by given factor h[1]
331 //-------------------------------------------------------------
334 const Matrix33 & setShear (const Vec2<S> &h);
337 //-----------------------------------------------------------
338 // Shear the matrix in x for each y coord. by given factor xy
339 //-----------------------------------------------------------
342 const Matrix33 & shear (const S &xy);
345 //-----------------------------------------------------------
346 // Shear the matrix in x for each y coord. by given factor xy
347 // and shear y for each x coord. by given factor yx
348 //-----------------------------------------------------------
351 const Matrix33 & shear (const Vec2<S> &h);
354 //-------------------------------------------------
355 // Limitations of type T (see also class limits<T>)
356 //-------------------------------------------------
358 static T baseTypeMin() {return limits<T>::min();}
359 static T baseTypeMax() {return limits<T>::max();}
360 static T baseTypeSmallest() {return limits<T>::smallest();}
361 static T baseTypeEpsilon() {return limits<T>::epsilon();}
365 template <class T> class Matrix44
369 //-------------------
370 // Access to elements
371 //-------------------
375 T * operator [] (int i);
376 const T * operator [] (int i) const;
395 Matrix44 (const T a[4][4]) ;
396 // a[0][0] a[0][1] a[0][2] a[0][3]
397 // a[1][0] a[1][1] a[1][2] a[1][3]
398 // a[2][0] a[2][1] a[2][2] a[2][3]
399 // a[3][0] a[3][1] a[3][2] a[3][3]
401 Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,
402 T i, T j, T k, T l, T m, T n, T o, T p);
409 Matrix44 (Matrix33<T> r, Vec3<T> t);
416 //--------------------------------
417 // Copy constructor and assignment
418 //--------------------------------
420 Matrix44 (const Matrix44 &v);
422 const Matrix44 & operator = (const Matrix44 &v);
423 const Matrix44 & operator = (T a);
426 //----------------------
427 // Compatibility with Sb
428 //----------------------
431 const T * getValue () const;
434 void getValue (Matrix44<S> &v) const;
436 Matrix44 & setValue (const Matrix44<S> &v);
439 Matrix44 & setTheMatrix (const Matrix44<S> &v);
452 bool operator == (const Matrix44 &v) const;
453 bool operator != (const Matrix44 &v) const;
455 //-----------------------------------------------------------------------
456 // Compare two matrices and test if they are "approximately equal":
458 // equalWithAbsError (m, e)
460 // Returns true if the coefficients of this and m are the same with
461 // an absolute error of no more than e, i.e., for all i, j
463 // abs (this[i][j] - m[i][j]) <= e
465 // equalWithRelError (m, e)
467 // Returns true if the coefficients of this and m are the same with
468 // a relative error of no more than e, i.e., for all i, j
470 // abs (this[i] - v[i][j]) <= e * abs (this[i][j])
471 //-----------------------------------------------------------------------
473 bool equalWithAbsError (const Matrix44<T> &v, T e) const;
474 bool equalWithRelError (const Matrix44<T> &v, T e) const;
477 //------------------------
478 // Component-wise addition
479 //------------------------
481 const Matrix44 & operator += (const Matrix44 &v);
482 const Matrix44 & operator += (T a);
483 Matrix44 operator + (const Matrix44 &v) const;
486 //---------------------------
487 // Component-wise subtraction
488 //---------------------------
490 const Matrix44 & operator -= (const Matrix44 &v);
491 const Matrix44 & operator -= (T a);
492 Matrix44 operator - (const Matrix44 &v) const;
495 //------------------------------------
496 // Component-wise multiplication by -1
497 //------------------------------------
499 Matrix44 operator - () const;
500 const Matrix44 & negate ();
503 //------------------------------
504 // Component-wise multiplication
505 //------------------------------
507 const Matrix44 & operator *= (T a);
508 Matrix44 operator * (T a) const;
511 //-----------------------------------
512 // Matrix-times-matrix multiplication
513 //-----------------------------------
515 const Matrix44 & operator *= (const Matrix44 &v);
516 Matrix44 operator * (const Matrix44 &v) const;
518 static void multiply (const Matrix44 &a, // assumes that
519 const Matrix44 &b, // &a != &c and
520 Matrix44 &c); // &b != &c.
523 //---------------------------------------------
524 // Vector-times-matrix multiplication; see also
525 // the "operator *" functions defined below.
526 //---------------------------------------------
529 void multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const;
532 void multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const;
535 //------------------------
536 // Component-wise division
537 //------------------------
539 const Matrix44 & operator /= (T a);
540 Matrix44 operator / (T a) const;
547 const Matrix44 & transpose ();
548 Matrix44 transposed () const;
551 //------------------------------------------------------------
552 // Inverse matrix: If singExc is false, inverting a singular
553 // matrix produces an identity matrix. If singExc is true,
554 // inverting a singular matrix throws a SingMatrixExc.
556 // inverse() and invert() invert matrices using determinants;
557 // gjInverse() and gjInvert() use the Gauss-Jordan method.
559 // inverse() and invert() are significantly faster than
560 // gjInverse() and gjInvert(), but the results may be slightly
563 //------------------------------------------------------------
565 const Matrix44 & invert (bool singExc = false)
566 throw (Iex::MathExc);
568 Matrix44<T> inverse (bool singExc = false) const
569 throw (Iex::MathExc);
571 const Matrix44 & gjInvert (bool singExc = false)
572 throw (Iex::MathExc);
574 Matrix44<T> gjInverse (bool singExc = false) const
575 throw (Iex::MathExc);
578 //--------------------------------------------------------
579 // Set matrix to rotation by XYZ euler angles (in radians)
580 //--------------------------------------------------------
583 const Matrix44 & setEulerAngles (const Vec3<S>& r);
586 //--------------------------------------------------------
587 // Set matrix to rotation around given axis by given angle
588 //--------------------------------------------------------
591 const Matrix44 & setAxisAngle (const Vec3<S>& ax, S ang);
594 //-------------------------------------------
595 // Rotate the matrix by XYZ euler angles in r
596 //-------------------------------------------
599 const Matrix44 & rotate (const Vec3<S> &r);
602 //--------------------------------------------
603 // Set matrix to scale by given uniform factor
604 //--------------------------------------------
606 const Matrix44 & setScale (T s);
609 //------------------------------------
610 // Set matrix to scale by given vector
611 //------------------------------------
614 const Matrix44 & setScale (const Vec3<S> &s);
617 //----------------------
618 // Scale the matrix by s
619 //----------------------
622 const Matrix44 & scale (const Vec3<S> &s);
625 //------------------------------------------
626 // Set matrix to translation by given vector
627 //------------------------------------------
630 const Matrix44 & setTranslation (const Vec3<S> &t);
633 //-----------------------------
634 // Return translation component
635 //-----------------------------
637 const Vec3<T> translation () const;
640 //--------------------------
641 // Translate the matrix by t
642 //--------------------------
645 const Matrix44 & translate (const Vec3<S> &t);
648 //-------------------------------------------------------------
649 // Set matrix to shear by given vector h. The resulting matrix
650 // will shear x for each y coord. by a factor of h[0] ;
651 // will shear x for each z coord. by a factor of h[1] ;
652 // will shear y for each z coord. by a factor of h[2] .
653 //-------------------------------------------------------------
656 const Matrix44 & setShear (const Vec3<S> &h);
659 //------------------------------------------------------------
660 // Set matrix to shear by given factors. The resulting matrix
661 // will shear x for each y coord. by a factor of h.xy ;
662 // will shear x for each z coord. by a factor of h.xz ;
663 // will shear y for each z coord. by a factor of h.yz ;
664 // will shear y for each x coord. by a factor of h.yx ;
665 // will shear z for each x coord. by a factor of h.zx ;
666 // will shear z for each y coord. by a factor of h.zy .
667 //------------------------------------------------------------
670 const Matrix44 & setShear (const Shear6<S> &h);
673 //--------------------------------------------------------
674 // Shear the matrix by given vector. The composed matrix
675 // will be <shear> * <this>, where the shear matrix ...
676 // will shear x for each y coord. by a factor of h[0] ;
677 // will shear x for each z coord. by a factor of h[1] ;
678 // will shear y for each z coord. by a factor of h[2] .
679 //--------------------------------------------------------
682 const Matrix44 & shear (const Vec3<S> &h);
685 //------------------------------------------------------------
686 // Shear the matrix by the given factors. The composed matrix
687 // will be <shear> * <this>, where the shear matrix ...
688 // will shear x for each y coord. by a factor of h.xy ;
689 // will shear x for each z coord. by a factor of h.xz ;
690 // will shear y for each z coord. by a factor of h.yz ;
691 // will shear y for each x coord. by a factor of h.yx ;
692 // will shear z for each x coord. by a factor of h.zx ;
693 // will shear z for each y coord. by a factor of h.zy .
694 //------------------------------------------------------------
697 const Matrix44 & shear (const Shear6<S> &h);
700 //-------------------------------------------------
701 // Limitations of type T (see also class limits<T>)
702 //-------------------------------------------------
704 static T baseTypeMin() {return limits<T>::min();}
705 static T baseTypeMax() {return limits<T>::max();}
706 static T baseTypeSmallest() {return limits<T>::smallest();}
707 static T baseTypeEpsilon() {return limits<T>::epsilon();}
716 std::ostream & operator << (std::ostream & s, const Matrix33<T> &m);
719 std::ostream & operator << (std::ostream & s, const Matrix44<T> &m);
722 //---------------------------------------------
723 // Vector-times-matrix multiplication operators
724 //---------------------------------------------
726 template <class S, class T>
727 const Vec2<S> & operator *= (Vec2<S> &v, const Matrix33<T> &m);
729 template <class S, class T>
730 Vec2<S> operator * (const Vec2<S> &v, const Matrix33<T> &m);
732 template <class S, class T>
733 const Vec3<S> & operator *= (Vec3<S> &v, const Matrix33<T> &m);
735 template <class S, class T>
736 Vec3<S> operator * (const Vec3<S> &v, const Matrix33<T> &m);
738 template <class S, class T>
739 const Vec3<S> & operator *= (Vec3<S> &v, const Matrix44<T> &m);
741 template <class S, class T>
742 Vec3<S> operator * (const Vec3<S> &v, const Matrix44<T> &m);
745 //-------------------------
746 // Typedefs for convenience
747 //-------------------------
749 typedef Matrix33 <float> M33f;
750 typedef Matrix33 <double> M33d;
751 typedef Matrix44 <float> M44f;
752 typedef Matrix44 <double> M44d;
755 //---------------------------
756 // Implementation of Matrix33
757 //---------------------------
761 Matrix33<T>::operator [] (int i)
768 Matrix33<T>::operator [] (int i) const
775 Matrix33<T>::Matrix33 ()
790 Matrix33<T>::Matrix33 (T a)
805 Matrix33<T>::Matrix33 (const T a[3][3])
820 Matrix33<T>::Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i)
835 Matrix33<T>::Matrix33 (const Matrix33 &v)
849 inline const Matrix33<T> &
850 Matrix33<T>::operator = (const Matrix33 &v)
865 inline const Matrix33<T> &
866 Matrix33<T>::operator = (T a)
882 Matrix33<T>::getValue ()
884 return (T *) &x[0][0];
889 Matrix33<T>::getValue () const
891 return (const T *) &x[0][0];
897 Matrix33<T>::getValue (Matrix33<S> &v) const
913 Matrix33<T>::setValue (const Matrix33<S> &v)
930 Matrix33<T>::setTheMatrix (const Matrix33<S> &v)
946 Matrix33<T>::makeIdentity()
961 Matrix33<T>::operator == (const Matrix33 &v) const
963 return x[0][0] == v.x[0][0] &&
964 x[0][1] == v.x[0][1] &&
965 x[0][2] == v.x[0][2] &&
966 x[1][0] == v.x[1][0] &&
967 x[1][1] == v.x[1][1] &&
968 x[1][2] == v.x[1][2] &&
969 x[2][0] == v.x[2][0] &&
970 x[2][1] == v.x[2][1] &&
971 x[2][2] == v.x[2][2];
976 Matrix33<T>::operator != (const Matrix33 &v) const
978 return x[0][0] != v.x[0][0] ||
979 x[0][1] != v.x[0][1] ||
980 x[0][2] != v.x[0][2] ||
981 x[1][0] != v.x[1][0] ||
982 x[1][1] != v.x[1][1] ||
983 x[1][2] != v.x[1][2] ||
984 x[2][0] != v.x[2][0] ||
985 x[2][1] != v.x[2][1] ||
986 x[2][2] != v.x[2][2];
991 Matrix33<T>::equalWithAbsError (const Matrix33<T> &m, T e) const
993 for (int i = 0; i < 3; i++)
994 for (int j = 0; j < 3; j++)
995 if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e))
1003 Matrix33<T>::equalWithRelError (const Matrix33<T> &m, T e) const
1005 for (int i = 0; i < 3; i++)
1006 for (int j = 0; j < 3; j++)
1007 if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e))
1015 Matrix33<T>::operator += (const Matrix33<T> &v)
1017 x[0][0] += v.x[0][0];
1018 x[0][1] += v.x[0][1];
1019 x[0][2] += v.x[0][2];
1020 x[1][0] += v.x[1][0];
1021 x[1][1] += v.x[1][1];
1022 x[1][2] += v.x[1][2];
1023 x[2][0] += v.x[2][0];
1024 x[2][1] += v.x[2][1];
1025 x[2][2] += v.x[2][2];
1032 Matrix33<T>::operator += (T a)
1049 Matrix33<T>::operator + (const Matrix33<T> &v) const
1051 return Matrix33 (x[0][0] + v.x[0][0],
1052 x[0][1] + v.x[0][1],
1053 x[0][2] + v.x[0][2],
1054 x[1][0] + v.x[1][0],
1055 x[1][1] + v.x[1][1],
1056 x[1][2] + v.x[1][2],
1057 x[2][0] + v.x[2][0],
1058 x[2][1] + v.x[2][1],
1059 x[2][2] + v.x[2][2]);
1064 Matrix33<T>::operator -= (const Matrix33<T> &v)
1066 x[0][0] -= v.x[0][0];
1067 x[0][1] -= v.x[0][1];
1068 x[0][2] -= v.x[0][2];
1069 x[1][0] -= v.x[1][0];
1070 x[1][1] -= v.x[1][1];
1071 x[1][2] -= v.x[1][2];
1072 x[2][0] -= v.x[2][0];
1073 x[2][1] -= v.x[2][1];
1074 x[2][2] -= v.x[2][2];
1081 Matrix33<T>::operator -= (T a)
1098 Matrix33<T>::operator - (const Matrix33<T> &v) const
1100 return Matrix33 (x[0][0] - v.x[0][0],
1101 x[0][1] - v.x[0][1],
1102 x[0][2] - v.x[0][2],
1103 x[1][0] - v.x[1][0],
1104 x[1][1] - v.x[1][1],
1105 x[1][2] - v.x[1][2],
1106 x[2][0] - v.x[2][0],
1107 x[2][1] - v.x[2][1],
1108 x[2][2] - v.x[2][2]);
1113 Matrix33<T>::operator - () const
1115 return Matrix33 (-x[0][0],
1128 Matrix33<T>::negate ()
1145 Matrix33<T>::operator *= (T a)
1162 Matrix33<T>::operator * (T a) const
1164 return Matrix33 (x[0][0] * a,
1177 operator * (T a, const Matrix33<T> &v)
1184 Matrix33<T>::operator *= (const Matrix33<T> &v)
1186 Matrix33 tmp (T (0));
1188 for (int i = 0; i < 3; i++)
1189 for (int j = 0; j < 3; j++)
1190 for (int k = 0; k < 3; k++)
1191 tmp.x[i][j] += x[i][k] * v.x[k][j];
1199 Matrix33<T>::operator * (const Matrix33<T> &v) const
1201 Matrix33 tmp (T (0));
1203 for (int i = 0; i < 3; i++)
1204 for (int j = 0; j < 3; j++)
1205 for (int k = 0; k < 3; k++)
1206 tmp.x[i][j] += x[i][k] * v.x[k][j];
1214 Matrix33<T>::multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const
1218 a = src[0] * x[0][0] + src[1] * x[1][0] + x[2][0];
1219 b = src[0] * x[0][1] + src[1] * x[1][1] + x[2][1];
1220 w = src[0] * x[0][2] + src[1] * x[1][2] + x[2][2];
1229 Matrix33<T>::multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const
1233 a = src[0] * x[0][0] + src[1] * x[1][0];
1234 b = src[0] * x[0][1] + src[1] * x[1][1];
1242 Matrix33<T>::operator /= (T a)
1259 Matrix33<T>::operator / (T a) const
1261 return Matrix33 (x[0][0] / a,
1274 Matrix33<T>::transpose ()
1276 Matrix33 tmp (x[0][0],
1291 Matrix33<T>::transposed () const
1293 return Matrix33 (x[0][0],
1306 Matrix33<T>::gjInvert (bool singExc) throw (Iex::MathExc)
1308 *this = gjInverse (singExc);
1314 Matrix33<T>::gjInverse (bool singExc) const throw (Iex::MathExc)
1320 // Forward elimination
1322 for (i = 0; i < 2 ; i++)
1326 T pivotsize = t[i][i];
1329 pivotsize = -pivotsize;
1331 for (j = i + 1; j < 3; j++)
1338 if (tmp > pivotsize)
1348 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
1355 for (j = 0; j < 3; j++)
1360 t[i][j] = t[pivot][j];
1364 s[i][j] = s[pivot][j];
1369 for (j = i + 1; j < 3; j++)
1371 T f = t[j][i] / t[i][i];
1373 for (k = 0; k < 3; k++)
1375 t[j][k] -= f * t[i][k];
1376 s[j][k] -= f * s[i][k];
1381 // Backward substitution
1383 for (i = 2; i >= 0; --i)
1387 if ((f = t[i][i]) == 0)
1390 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
1395 for (j = 0; j < 3; j++)
1401 for (j = 0; j < i; j++)
1405 for (k = 0; k < 3; k++)
1407 t[j][k] -= f * t[i][k];
1408 s[j][k] -= f * s[i][k];
1418 Matrix33<T>::invert (bool singExc) throw (Iex::MathExc)
1420 *this = inverse (singExc);
1426 Matrix33<T>::inverse (bool singExc) const throw (Iex::MathExc)
1428 if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1)
1430 Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
1431 x[2][1] * x[0][2] - x[0][1] * x[2][2],
1432 x[0][1] * x[1][2] - x[1][1] * x[0][2],
1434 x[2][0] * x[1][2] - x[1][0] * x[2][2],
1435 x[0][0] * x[2][2] - x[2][0] * x[0][2],
1436 x[1][0] * x[0][2] - x[0][0] * x[1][2],
1438 x[1][0] * x[2][1] - x[2][0] * x[1][1],
1439 x[2][0] * x[0][1] - x[0][0] * x[2][1],
1440 x[0][0] * x[1][1] - x[1][0] * x[0][1]);
1442 T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
1444 if (Imath::abs (r) >= 1)
1446 for (int i = 0; i < 3; ++i)
1448 for (int j = 0; j < 3; ++j)
1456 T mr = Imath::abs (r) / limits<T>::smallest();
1458 for (int i = 0; i < 3; ++i)
1460 for (int j = 0; j < 3; ++j)
1462 if (mr > Imath::abs (s[i][j]))
1469 throw SingMatrixExc ("Cannot invert "
1470 "singular matrix.");
1481 Matrix33 s ( x[1][1],
1493 T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
1495 if (Imath::abs (r) >= 1)
1497 for (int i = 0; i < 2; ++i)
1499 for (int j = 0; j < 2; ++j)
1507 T mr = Imath::abs (r) / limits<T>::smallest();
1509 for (int i = 0; i < 2; ++i)
1511 for (int j = 0; j < 2; ++j)
1513 if (mr > Imath::abs (s[i][j]))
1520 throw SingMatrixExc ("Cannot invert "
1521 "singular matrix.");
1528 s[2][0] = -x[2][0] * s[0][0] - x[2][1] * s[1][0];
1529 s[2][1] = -x[2][0] * s[0][1] - x[2][1] * s[1][1];
1538 Matrix33<T>::setRotation (S r)
1542 cos_r = Math<T>::cos (r);
1543 sin_r = Math<T>::sin (r);
1563 Matrix33<T>::rotate (S r)
1565 *this *= Matrix33<T>().setRotation (r);
1571 Matrix33<T>::setScale (T s)
1591 Matrix33<T>::setScale (const Vec2<S> &s)
1611 Matrix33<T>::scale (const Vec2<S> &s)
1627 Matrix33<T>::setTranslation (const Vec2<S> &t)
1646 Matrix33<T>::translation () const
1648 return Vec2<T> (x[2][0], x[2][1]);
1654 Matrix33<T>::translate (const Vec2<S> &t)
1656 x[2][0] += t[0] * x[0][0] + t[1] * x[1][0];
1657 x[2][1] += t[0] * x[0][1] + t[1] * x[1][1];
1658 x[2][2] += t[0] * x[0][2] + t[1] * x[1][2];
1666 Matrix33<T>::setShear (const S &xy)
1686 Matrix33<T>::setShear (const Vec2<S> &h)
1706 Matrix33<T>::shear (const S &xy)
1709 // In this case, we don't need a temp. copy of the matrix
1710 // because we never use a value on the RHS after we've
1711 // changed it on the LHS.
1714 x[1][0] += xy * x[0][0];
1715 x[1][1] += xy * x[0][1];
1716 x[1][2] += xy * x[0][2];
1724 Matrix33<T>::shear (const Vec2<S> &h)
1726 Matrix33<T> P (*this);
1728 x[0][0] = P[0][0] + h[1] * P[1][0];
1729 x[0][1] = P[0][1] + h[1] * P[1][1];
1730 x[0][2] = P[0][2] + h[1] * P[1][2];
1732 x[1][0] = P[1][0] + h[0] * P[0][0];
1733 x[1][1] = P[1][1] + h[0] * P[0][1];
1734 x[1][2] = P[1][2] + h[0] * P[0][2];
1740 //---------------------------
1741 // Implementation of Matrix44
1742 //---------------------------
1746 Matrix44<T>::operator [] (int i)
1753 Matrix44<T>::operator [] (int i) const
1760 Matrix44<T>::Matrix44 ()
1782 Matrix44<T>::Matrix44 (T a)
1804 Matrix44<T>::Matrix44 (const T a[4][4])
1826 Matrix44<T>::Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,
1827 T i, T j, T k, T l, T m, T n, T o, T p)
1850 Matrix44<T>::Matrix44 (Matrix33<T> r, Vec3<T> t)
1872 Matrix44<T>::Matrix44 (const Matrix44 &v)
1874 x[0][0] = v.x[0][0];
1875 x[0][1] = v.x[0][1];
1876 x[0][2] = v.x[0][2];
1877 x[0][3] = v.x[0][3];
1878 x[1][0] = v.x[1][0];
1879 x[1][1] = v.x[1][1];
1880 x[1][2] = v.x[1][2];
1881 x[1][3] = v.x[1][3];
1882 x[2][0] = v.x[2][0];
1883 x[2][1] = v.x[2][1];
1884 x[2][2] = v.x[2][2];
1885 x[2][3] = v.x[2][3];
1886 x[3][0] = v.x[3][0];
1887 x[3][1] = v.x[3][1];
1888 x[3][2] = v.x[3][2];
1889 x[3][3] = v.x[3][3];
1893 inline const Matrix44<T> &
1894 Matrix44<T>::operator = (const Matrix44 &v)
1896 x[0][0] = v.x[0][0];
1897 x[0][1] = v.x[0][1];
1898 x[0][2] = v.x[0][2];
1899 x[0][3] = v.x[0][3];
1900 x[1][0] = v.x[1][0];
1901 x[1][1] = v.x[1][1];
1902 x[1][2] = v.x[1][2];
1903 x[1][3] = v.x[1][3];
1904 x[2][0] = v.x[2][0];
1905 x[2][1] = v.x[2][1];
1906 x[2][2] = v.x[2][2];
1907 x[2][3] = v.x[2][3];
1908 x[3][0] = v.x[3][0];
1909 x[3][1] = v.x[3][1];
1910 x[3][2] = v.x[3][2];
1911 x[3][3] = v.x[3][3];
1916 inline const Matrix44<T> &
1917 Matrix44<T>::operator = (T a)
1940 Matrix44<T>::getValue ()
1942 return (T *) &x[0][0];
1947 Matrix44<T>::getValue () const
1949 return (const T *) &x[0][0];
1955 Matrix44<T>::getValue (Matrix44<S> &v) const
1957 v.x[0][0] = x[0][0];
1958 v.x[0][1] = x[0][1];
1959 v.x[0][2] = x[0][2];
1960 v.x[0][3] = x[0][3];
1961 v.x[1][0] = x[1][0];
1962 v.x[1][1] = x[1][1];
1963 v.x[1][2] = x[1][2];
1964 v.x[1][3] = x[1][3];
1965 v.x[2][0] = x[2][0];
1966 v.x[2][1] = x[2][1];
1967 v.x[2][2] = x[2][2];
1968 v.x[2][3] = x[2][3];
1969 v.x[3][0] = x[3][0];
1970 v.x[3][1] = x[3][1];
1971 v.x[3][2] = x[3][2];
1972 v.x[3][3] = x[3][3];
1977 inline Matrix44<T> &
1978 Matrix44<T>::setValue (const Matrix44<S> &v)
1980 x[0][0] = v.x[0][0];
1981 x[0][1] = v.x[0][1];
1982 x[0][2] = v.x[0][2];
1983 x[0][3] = v.x[0][3];
1984 x[1][0] = v.x[1][0];
1985 x[1][1] = v.x[1][1];
1986 x[1][2] = v.x[1][2];
1987 x[1][3] = v.x[1][3];
1988 x[2][0] = v.x[2][0];
1989 x[2][1] = v.x[2][1];
1990 x[2][2] = v.x[2][2];
1991 x[2][3] = v.x[2][3];
1992 x[3][0] = v.x[3][0];
1993 x[3][1] = v.x[3][1];
1994 x[3][2] = v.x[3][2];
1995 x[3][3] = v.x[3][3];
2001 inline Matrix44<T> &
2002 Matrix44<T>::setTheMatrix (const Matrix44<S> &v)
2004 x[0][0] = v.x[0][0];
2005 x[0][1] = v.x[0][1];
2006 x[0][2] = v.x[0][2];
2007 x[0][3] = v.x[0][3];
2008 x[1][0] = v.x[1][0];
2009 x[1][1] = v.x[1][1];
2010 x[1][2] = v.x[1][2];
2011 x[1][3] = v.x[1][3];
2012 x[2][0] = v.x[2][0];
2013 x[2][1] = v.x[2][1];
2014 x[2][2] = v.x[2][2];
2015 x[2][3] = v.x[2][3];
2016 x[3][0] = v.x[3][0];
2017 x[3][1] = v.x[3][1];
2018 x[3][2] = v.x[3][2];
2019 x[3][3] = v.x[3][3];
2025 Matrix44<T>::makeIdentity()
2047 Matrix44<T>::operator == (const Matrix44 &v) const
2049 return x[0][0] == v.x[0][0] &&
2050 x[0][1] == v.x[0][1] &&
2051 x[0][2] == v.x[0][2] &&
2052 x[0][3] == v.x[0][3] &&
2053 x[1][0] == v.x[1][0] &&
2054 x[1][1] == v.x[1][1] &&
2055 x[1][2] == v.x[1][2] &&
2056 x[1][3] == v.x[1][3] &&
2057 x[2][0] == v.x[2][0] &&
2058 x[2][1] == v.x[2][1] &&
2059 x[2][2] == v.x[2][2] &&
2060 x[2][3] == v.x[2][3] &&
2061 x[3][0] == v.x[3][0] &&
2062 x[3][1] == v.x[3][1] &&
2063 x[3][2] == v.x[3][2] &&
2064 x[3][3] == v.x[3][3];
2069 Matrix44<T>::operator != (const Matrix44 &v) const
2071 return x[0][0] != v.x[0][0] ||
2072 x[0][1] != v.x[0][1] ||
2073 x[0][2] != v.x[0][2] ||
2074 x[0][3] != v.x[0][3] ||
2075 x[1][0] != v.x[1][0] ||
2076 x[1][1] != v.x[1][1] ||
2077 x[1][2] != v.x[1][2] ||
2078 x[1][3] != v.x[1][3] ||
2079 x[2][0] != v.x[2][0] ||
2080 x[2][1] != v.x[2][1] ||
2081 x[2][2] != v.x[2][2] ||
2082 x[2][3] != v.x[2][3] ||
2083 x[3][0] != v.x[3][0] ||
2084 x[3][1] != v.x[3][1] ||
2085 x[3][2] != v.x[3][2] ||
2086 x[3][3] != v.x[3][3];
2091 Matrix44<T>::equalWithAbsError (const Matrix44<T> &m, T e) const
2093 for (int i = 0; i < 4; i++)
2094 for (int j = 0; j < 4; j++)
2095 if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e))
2103 Matrix44<T>::equalWithRelError (const Matrix44<T> &m, T e) const
2105 for (int i = 0; i < 4; i++)
2106 for (int j = 0; j < 4; j++)
2107 if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e))
2115 Matrix44<T>::operator += (const Matrix44<T> &v)
2117 x[0][0] += v.x[0][0];
2118 x[0][1] += v.x[0][1];
2119 x[0][2] += v.x[0][2];
2120 x[0][3] += v.x[0][3];
2121 x[1][0] += v.x[1][0];
2122 x[1][1] += v.x[1][1];
2123 x[1][2] += v.x[1][2];
2124 x[1][3] += v.x[1][3];
2125 x[2][0] += v.x[2][0];
2126 x[2][1] += v.x[2][1];
2127 x[2][2] += v.x[2][2];
2128 x[2][3] += v.x[2][3];
2129 x[3][0] += v.x[3][0];
2130 x[3][1] += v.x[3][1];
2131 x[3][2] += v.x[3][2];
2132 x[3][3] += v.x[3][3];
2139 Matrix44<T>::operator += (T a)
2163 Matrix44<T>::operator + (const Matrix44<T> &v) const
2165 return Matrix44 (x[0][0] + v.x[0][0],
2166 x[0][1] + v.x[0][1],
2167 x[0][2] + v.x[0][2],
2168 x[0][3] + v.x[0][3],
2169 x[1][0] + v.x[1][0],
2170 x[1][1] + v.x[1][1],
2171 x[1][2] + v.x[1][2],
2172 x[1][3] + v.x[1][3],
2173 x[2][0] + v.x[2][0],
2174 x[2][1] + v.x[2][1],
2175 x[2][2] + v.x[2][2],
2176 x[2][3] + v.x[2][3],
2177 x[3][0] + v.x[3][0],
2178 x[3][1] + v.x[3][1],
2179 x[3][2] + v.x[3][2],
2180 x[3][3] + v.x[3][3]);
2185 Matrix44<T>::operator -= (const Matrix44<T> &v)
2187 x[0][0] -= v.x[0][0];
2188 x[0][1] -= v.x[0][1];
2189 x[0][2] -= v.x[0][2];
2190 x[0][3] -= v.x[0][3];
2191 x[1][0] -= v.x[1][0];
2192 x[1][1] -= v.x[1][1];
2193 x[1][2] -= v.x[1][2];
2194 x[1][3] -= v.x[1][3];
2195 x[2][0] -= v.x[2][0];
2196 x[2][1] -= v.x[2][1];
2197 x[2][2] -= v.x[2][2];
2198 x[2][3] -= v.x[2][3];
2199 x[3][0] -= v.x[3][0];
2200 x[3][1] -= v.x[3][1];
2201 x[3][2] -= v.x[3][2];
2202 x[3][3] -= v.x[3][3];
2209 Matrix44<T>::operator -= (T a)
2233 Matrix44<T>::operator - (const Matrix44<T> &v) const
2235 return Matrix44 (x[0][0] - v.x[0][0],
2236 x[0][1] - v.x[0][1],
2237 x[0][2] - v.x[0][2],
2238 x[0][3] - v.x[0][3],
2239 x[1][0] - v.x[1][0],
2240 x[1][1] - v.x[1][1],
2241 x[1][2] - v.x[1][2],
2242 x[1][3] - v.x[1][3],
2243 x[2][0] - v.x[2][0],
2244 x[2][1] - v.x[2][1],
2245 x[2][2] - v.x[2][2],
2246 x[2][3] - v.x[2][3],
2247 x[3][0] - v.x[3][0],
2248 x[3][1] - v.x[3][1],
2249 x[3][2] - v.x[3][2],
2250 x[3][3] - v.x[3][3]);
2255 Matrix44<T>::operator - () const
2257 return Matrix44 (-x[0][0],
2277 Matrix44<T>::negate ()
2301 Matrix44<T>::operator *= (T a)
2325 Matrix44<T>::operator * (T a) const
2327 return Matrix44 (x[0][0] * a,
2347 operator * (T a, const Matrix44<T> &v)
2353 inline const Matrix44<T> &
2354 Matrix44<T>::operator *= (const Matrix44<T> &v)
2356 Matrix44 tmp (T (0));
2358 multiply (*this, v, tmp);
2365 Matrix44<T>::operator * (const Matrix44<T> &v) const
2367 Matrix44 tmp (T (0));
2369 multiply (*this, v, tmp);
2375 Matrix44<T>::multiply (const Matrix44<T> &a,
2376 const Matrix44<T> &b,
2379 register const T * restrict ap = &a.x[0][0];
2380 register const T * restrict bp = &b.x[0][0];
2381 register T * restrict cp = &c.x[0][0];
2383 register T a0, a1, a2, a3;
2390 cp[0] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2391 cp[1] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2392 cp[2] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2393 cp[3] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2400 cp[4] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2401 cp[5] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2402 cp[6] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2403 cp[7] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2410 cp[8] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2411 cp[9] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2412 cp[10] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2413 cp[11] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2420 cp[12] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2421 cp[13] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2422 cp[14] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2423 cp[15] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2426 template <class T> template <class S>
2428 Matrix44<T>::multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const
2432 a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0] + x[3][0];
2433 b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1] + x[3][1];
2434 c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2] + x[3][2];
2435 w = src[0] * x[0][3] + src[1] * x[1][3] + src[2] * x[2][3] + x[3][3];
2442 template <class T> template <class S>
2444 Matrix44<T>::multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const
2448 a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0];
2449 b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1];
2450 c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2];
2459 Matrix44<T>::operator /= (T a)
2483 Matrix44<T>::operator / (T a) const
2485 return Matrix44 (x[0][0] / a,
2505 Matrix44<T>::transpose ()
2507 Matrix44 tmp (x[0][0],
2529 Matrix44<T>::transposed () const
2531 return Matrix44 (x[0][0],
2551 Matrix44<T>::gjInvert (bool singExc) throw (Iex::MathExc)
2553 *this = gjInverse (singExc);
2559 Matrix44<T>::gjInverse (bool singExc) const throw (Iex::MathExc)
2565 // Forward elimination
2567 for (i = 0; i < 3 ; i++)
2571 T pivotsize = t[i][i];
2574 pivotsize = -pivotsize;
2576 for (j = i + 1; j < 4; j++)
2583 if (tmp > pivotsize)
2593 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
2600 for (j = 0; j < 4; j++)
2605 t[i][j] = t[pivot][j];
2609 s[i][j] = s[pivot][j];
2614 for (j = i + 1; j < 4; j++)
2616 T f = t[j][i] / t[i][i];
2618 for (k = 0; k < 4; k++)
2620 t[j][k] -= f * t[i][k];
2621 s[j][k] -= f * s[i][k];
2626 // Backward substitution
2628 for (i = 3; i >= 0; --i)
2632 if ((f = t[i][i]) == 0)
2635 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
2640 for (j = 0; j < 4; j++)
2646 for (j = 0; j < i; j++)
2650 for (k = 0; k < 4; k++)
2652 t[j][k] -= f * t[i][k];
2653 s[j][k] -= f * s[i][k];
2663 Matrix44<T>::invert (bool singExc) throw (Iex::MathExc)
2665 *this = inverse (singExc);
2671 Matrix44<T>::inverse (bool singExc) const throw (Iex::MathExc)
2673 if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)
2674 return gjInverse(singExc);
2676 Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
2677 x[2][1] * x[0][2] - x[0][1] * x[2][2],
2678 x[0][1] * x[1][2] - x[1][1] * x[0][2],
2681 x[2][0] * x[1][2] - x[1][0] * x[2][2],
2682 x[0][0] * x[2][2] - x[2][0] * x[0][2],
2683 x[1][0] * x[0][2] - x[0][0] * x[1][2],
2686 x[1][0] * x[2][1] - x[2][0] * x[1][1],
2687 x[2][0] * x[0][1] - x[0][0] * x[2][1],
2688 x[0][0] * x[1][1] - x[1][0] * x[0][1],
2696 T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
2698 if (Imath::abs (r) >= 1)
2700 for (int i = 0; i < 3; ++i)
2702 for (int j = 0; j < 3; ++j)
2710 T mr = Imath::abs (r) / limits<T>::smallest();
2712 for (int i = 0; i < 3; ++i)
2714 for (int j = 0; j < 3; ++j)
2716 if (mr > Imath::abs (s[i][j]))
2723 throw SingMatrixExc ("Cannot invert singular matrix.");
2731 s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0];
2732 s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1];
2733 s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2];
2741 Matrix44<T>::setEulerAngles (const Vec3<S>& r)
2743 S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
2745 cos_rz = Math<T>::cos (r[2]);
2746 cos_ry = Math<T>::cos (r[1]);
2747 cos_rx = Math<T>::cos (r[0]);
2749 sin_rz = Math<T>::sin (r[2]);
2750 sin_ry = Math<T>::sin (r[1]);
2751 sin_rx = Math<T>::sin (r[0]);
2753 x[0][0] = cos_rz * cos_ry;
2754 x[0][1] = sin_rz * cos_ry;
2758 x[1][0] = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
2759 x[1][1] = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
2760 x[1][2] = cos_ry * sin_rx;
2763 x[2][0] = sin_rz * sin_rx + cos_rz * sin_ry * cos_rx;
2764 x[2][1] = -cos_rz * sin_rx + sin_rz * sin_ry * cos_rx;
2765 x[2][2] = cos_ry * cos_rx;
2779 Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle)
2781 Vec3<S> unit (axis.normalized());
2782 S sine = Math<T>::sin (angle);
2783 S cosine = Math<T>::cos (angle);
2785 x[0][0] = unit[0] * unit[0] * (1 - cosine) + cosine;
2786 x[0][1] = unit[0] * unit[1] * (1 - cosine) + unit[2] * sine;
2787 x[0][2] = unit[0] * unit[2] * (1 - cosine) - unit[1] * sine;
2790 x[1][0] = unit[0] * unit[1] * (1 - cosine) - unit[2] * sine;
2791 x[1][1] = unit[1] * unit[1] * (1 - cosine) + cosine;
2792 x[1][2] = unit[1] * unit[2] * (1 - cosine) + unit[0] * sine;
2795 x[2][0] = unit[0] * unit[2] * (1 - cosine) + unit[1] * sine;
2796 x[2][1] = unit[1] * unit[2] * (1 - cosine) - unit[0] * sine;
2797 x[2][2] = unit[2] * unit[2] * (1 - cosine) + cosine;
2811 Matrix44<T>::rotate (const Vec3<S> &r)
2813 S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
2818 cos_rz = Math<S>::cos (r[2]);
2819 cos_ry = Math<S>::cos (r[1]);
2820 cos_rx = Math<S>::cos (r[0]);
2822 sin_rz = Math<S>::sin (r[2]);
2823 sin_ry = Math<S>::sin (r[1]);
2824 sin_rx = Math<S>::sin (r[0]);
2826 m00 = cos_rz * cos_ry;
2827 m01 = sin_rz * cos_ry;
2829 m10 = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
2830 m11 = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
2831 m12 = cos_ry * sin_rx;
2832 m20 = -sin_rz * -sin_rx + cos_rz * sin_ry * cos_rx;
2833 m21 = cos_rz * -sin_rx + sin_rz * sin_ry * cos_rx;
2834 m22 = cos_ry * cos_rx;
2836 Matrix44<T> P (*this);
2838 x[0][0] = P[0][0] * m00 + P[1][0] * m01 + P[2][0] * m02;
2839 x[0][1] = P[0][1] * m00 + P[1][1] * m01 + P[2][1] * m02;
2840 x[0][2] = P[0][2] * m00 + P[1][2] * m01 + P[2][2] * m02;
2841 x[0][3] = P[0][3] * m00 + P[1][3] * m01 + P[2][3] * m02;
2843 x[1][0] = P[0][0] * m10 + P[1][0] * m11 + P[2][0] * m12;
2844 x[1][1] = P[0][1] * m10 + P[1][1] * m11 + P[2][1] * m12;
2845 x[1][2] = P[0][2] * m10 + P[1][2] * m11 + P[2][2] * m12;
2846 x[1][3] = P[0][3] * m10 + P[1][3] * m11 + P[2][3] * m12;
2848 x[2][0] = P[0][0] * m20 + P[1][0] * m21 + P[2][0] * m22;
2849 x[2][1] = P[0][1] * m20 + P[1][1] * m21 + P[2][1] * m22;
2850 x[2][2] = P[0][2] * m20 + P[1][2] * m21 + P[2][2] * m22;
2851 x[2][3] = P[0][3] * m20 + P[1][3] * m21 + P[2][3] * m22;
2858 Matrix44<T>::setScale (T s)
2886 Matrix44<T>::setScale (const Vec3<S> &s)
2914 Matrix44<T>::scale (const Vec3<S> &s)
2937 Matrix44<T>::setTranslation (const Vec3<S> &t)
2963 inline const Vec3<T>
2964 Matrix44<T>::translation () const
2966 return Vec3<T> (x[3][0], x[3][1], x[3][2]);
2972 Matrix44<T>::translate (const Vec3<S> &t)
2974 x[3][0] += t[0] * x[0][0] + t[1] * x[1][0] + t[2] * x[2][0];
2975 x[3][1] += t[0] * x[0][1] + t[1] * x[1][1] + t[2] * x[2][1];
2976 x[3][2] += t[0] * x[0][2] + t[1] * x[1][2] + t[2] * x[2][2];
2977 x[3][3] += t[0] * x[0][3] + t[1] * x[1][3] + t[2] * x[2][3];
2985 Matrix44<T>::setShear (const Vec3<S> &h)
3013 Matrix44<T>::setShear (const Shear6<S> &h)
3041 Matrix44<T>::shear (const Vec3<S> &h)
3044 // In this case, we don't need a temp. copy of the matrix
3045 // because we never use a value on the RHS after we've
3046 // changed it on the LHS.
3049 for (int i=0; i < 4; i++)
3051 x[2][i] += h[1] * x[0][i] + h[2] * x[1][i];
3052 x[1][i] += h[0] * x[0][i];
3061 Matrix44<T>::shear (const Shear6<S> &h)
3063 Matrix44<T> P (*this);
3065 for (int i=0; i < 4; i++)
3067 x[0][i] = P[0][i] + h.yx * P[1][i] + h.zx * P[2][i];
3068 x[1][i] = h.xy * P[0][i] + P[1][i] + h.zy * P[2][i];
3069 x[2][i] = h.xz * P[0][i] + h.yz * P[1][i] + P[2][i];
3076 //--------------------------------
3077 // Implementation of stream output
3078 //--------------------------------
3082 operator << (std::ostream &s, const Matrix33<T> &m)
3084 std::ios_base::fmtflags oldFlags = s.flags();
3087 if (s.flags() & std::ios_base::fixed)
3089 s.setf (std::ios_base::showpoint);
3090 width = s.precision() + 5;
3094 s.setf (std::ios_base::scientific);
3095 s.setf (std::ios_base::showpoint);
3096 width = s.precision() + 8;
3099 s << "(" << std::setw (width) << m[0][0] <<
3100 " " << std::setw (width) << m[0][1] <<
3101 " " << std::setw (width) << m[0][2] << "\n" <<
3103 " " << std::setw (width) << m[1][0] <<
3104 " " << std::setw (width) << m[1][1] <<
3105 " " << std::setw (width) << m[1][2] << "\n" <<
3107 " " << std::setw (width) << m[2][0] <<
3108 " " << std::setw (width) << m[2][1] <<
3109 " " << std::setw (width) << m[2][2] << ")\n";
3117 operator << (std::ostream &s, const Matrix44<T> &m)
3119 std::ios_base::fmtflags oldFlags = s.flags();
3122 if (s.flags() & std::ios_base::fixed)
3124 s.setf (std::ios_base::showpoint);
3125 width = s.precision() + 5;
3129 s.setf (std::ios_base::scientific);
3130 s.setf (std::ios_base::showpoint);
3131 width = s.precision() + 8;
3134 s << "(" << std::setw (width) << m[0][0] <<
3135 " " << std::setw (width) << m[0][1] <<
3136 " " << std::setw (width) << m[0][2] <<
3137 " " << std::setw (width) << m[0][3] << "\n" <<
3139 " " << std::setw (width) << m[1][0] <<
3140 " " << std::setw (width) << m[1][1] <<
3141 " " << std::setw (width) << m[1][2] <<
3142 " " << std::setw (width) << m[1][3] << "\n" <<
3144 " " << std::setw (width) << m[2][0] <<
3145 " " << std::setw (width) << m[2][1] <<
3146 " " << std::setw (width) << m[2][2] <<
3147 " " << std::setw (width) << m[2][3] << "\n" <<
3149 " " << std::setw (width) << m[3][0] <<
3150 " " << std::setw (width) << m[3][1] <<
3151 " " << std::setw (width) << m[3][2] <<
3152 " " << std::setw (width) << m[3][3] << ")\n";
3159 //---------------------------------------------------------------
3160 // Implementation of vector-times-matrix multiplication operators
3161 //---------------------------------------------------------------
3163 template <class S, class T>
3164 inline const Vec2<S> &
3165 operator *= (Vec2<S> &v, const Matrix33<T> &m)
3167 S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
3168 S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
3169 S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
3177 template <class S, class T>
3179 operator * (const Vec2<S> &v, const Matrix33<T> &m)
3181 S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
3182 S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
3183 S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
3185 return Vec2<S> (x / w, y / w);
3189 template <class S, class T>
3190 inline const Vec3<S> &
3191 operator *= (Vec3<S> &v, const Matrix33<T> &m)
3193 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
3194 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
3195 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
3205 template <class S, class T>
3207 operator * (const Vec3<S> &v, const Matrix33<T> &m)
3209 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
3210 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
3211 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
3213 return Vec3<S> (x, y, z);
3217 template <class S, class T>
3218 inline const Vec3<S> &
3219 operator *= (Vec3<S> &v, const Matrix44<T> &m)
3221 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
3222 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
3223 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
3224 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
3233 template <class S, class T>
3235 operator * (const Vec3<S> &v, const Matrix44<T> &m)
3237 S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
3238 S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
3239 S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
3240 S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
3242 return Vec3<S> (x / w, y / w, z / w);
3245 } // namespace Imath