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