source: NonGTP/OpenEXR/include/Imath/ImathMatrix.h @ 855

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