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