/src/freeimage-svn/FreeImage/trunk/Source/OpenEXR/Imath/ImathMatrix.h
Line  | Count  | Source  | 
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  |  |  | 
64  |  | IMATH_INTERNAL_NAMESPACE_HEADER_ENTER  | 
65  |  |  | 
66  |  | enum Uninitialized {UNINITIALIZED}; | 
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  |  |  | 
87  |  |     Matrix33 (Uninitialized) {} | 
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();} | 
407  |  |     static T            baseTypeSmallest()      {return limits<T>::smallest();} | 
408  |  |     static T            baseTypeEpsilon()       {return limits<T>::epsilon();} | 
409  |  |  | 
410  |  |     typedef T   BaseType;  | 
411  |  |     typedef Vec3<T> BaseVecType;  | 
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  |  |  | 
447  |  |     Matrix44 (Uninitialized) {} | 
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  |  |  | 
475  |  |     Matrix44 (Matrix33<T> r, Vec3<T> t);  | 
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();} | 
806  |  |     static T            baseTypeSmallest()      {return limits<T>::smallest();} | 
807  |  |     static T            baseTypeEpsilon()       {return limits<T>::epsilon();} | 
808  |  |  | 
809  |  |     typedef T   BaseType;  | 
810  |  |     typedef Vec4<T> BaseVecType;  | 
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  |  |  | 
871  |  | typedef Matrix33 <float>  M33f;  | 
872  |  | typedef Matrix33 <double> M33d;  | 
873  |  | typedef Matrix44 <float>  M44f;  | 
874  |  | typedef Matrix44 <double> M44d;  | 
875  |  |  | 
876  |  |  | 
877  |  | //---------------------------  | 
878  |  | // Implementation of Matrix33  | 
879  |  | //---------------------------  | 
880  |  |  | 
881  |  | template <class T>  | 
882  |  | inline T *  | 
883  |  | Matrix33<T>::operator [] (int i)  | 
884  | 0  | { | 
885  | 0  |     return x[i];  | 
886  | 0  | } Unexecuted instantiation: Imath_2_2::Matrix33<float>::operator[](int) Unexecuted instantiation: Imath_2_2::Matrix33<double>::operator[](int)  | 
887  |  |  | 
888  |  | template <class T>  | 
889  |  | inline const T *  | 
890  |  | Matrix33<T>::operator [] (int i) const  | 
891  | 0  | { | 
892  | 0  |     return x[i];  | 
893  | 0  | } Unexecuted instantiation: Imath_2_2::Matrix33<float>::operator[](int) const Unexecuted instantiation: Imath_2_2::Matrix33<double>::operator[](int) const  | 
894  |  |  | 
895  |  | template <class T>  | 
896  |  | inline  | 
897  |  | Matrix33<T>::Matrix33 ()  | 
898  | 0  | { | 
899  | 0  |     memset (x, 0, sizeof (x));  | 
900  | 0  |     x[0][0] = 1;  | 
901  | 0  |     x[1][1] = 1;  | 
902  | 0  |     x[2][2] = 1;  | 
903  | 0  | } Unexecuted instantiation: Imath_2_2::Matrix33<double>::Matrix33() Unexecuted instantiation: Imath_2_2::Matrix33<float>::Matrix33()  | 
904  |  |  | 
905  |  | template <class T>  | 
906  |  | inline  | 
907  |  | Matrix33<T>::Matrix33 (T a)  | 
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  | 
944  |  | Matrix33<T>::Matrix33 (const Matrix33 &v)  | 
945  |  | { | 
946  |  |     memcpy (x, v.x, sizeof (x));  | 
947  |  | }  | 
948  |  |  | 
949  |  | template <class T>  | 
950  |  | template <class S>  | 
951  |  | inline  | 
952  |  | Matrix33<T>::Matrix33 (const Matrix33<S> &v)  | 
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> &  | 
967  |  | Matrix33<T>::operator = (const Matrix33 &v)  | 
968  | 0  | { | 
969  | 0  |     memcpy (x, v.x, sizeof (x));  | 
970  | 0  |     return *this;  | 
971  | 0  | } Unexecuted instantiation: Imath_2_2::Matrix33<double>::operator=(Imath_2_2::Matrix33<double> const&) Unexecuted instantiation: Imath_2_2::Matrix33<float>::operator=(Imath_2_2::Matrix33<float> const&)  | 
972  |  |  | 
973  |  | template <class T>  | 
974  |  | inline const Matrix33<T> &  | 
975  |  | Matrix33<T>::operator = (T a)  | 
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 *  | 
991  |  | Matrix33<T>::getValue ()  | 
992  |  | { | 
993  |  |     return (T *) &x[0][0];  | 
994  |  | }  | 
995  |  |  | 
996  |  | template <class T>  | 
997  |  | inline const T *  | 
998  |  | Matrix33<T>::getValue () const  | 
999  |  | { | 
1000  |  |     return (const T *) &x[0][0];  | 
1001  |  | }  | 
1002  |  |  | 
1003  |  | template <class T>  | 
1004  |  | template <class S>  | 
1005  |  | inline void  | 
1006  |  | Matrix33<T>::getValue (Matrix33<S> &v) const  | 
1007  |  | { | 
1008  |  |     if (isSameType<S,T>::value)  | 
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> &  | 
1029  |  | Matrix33<T>::setValue (const Matrix33<S> &v)  | 
1030  |  | { | 
1031  |  |     if (isSameType<S,T>::value)  | 
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> &  | 
1054  |  | Matrix33<T>::setTheMatrix (const Matrix33<S> &v)  | 
1055  |  | { | 
1056  |  |     if (isSameType<S,T>::value)  | 
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  | 
1078  |  | Matrix33<T>::makeIdentity()  | 
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  | 
1088  |  | Matrix33<T>::operator == (const Matrix33 &v) const  | 
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  | 
1103  |  | Matrix33<T>::operator != (const Matrix33 &v) const  | 
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  | 
1118  |  | Matrix33<T>::equalWithAbsError (const Matrix33<T> &m, T e) const  | 
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  | 
1130  |  | Matrix33<T>::equalWithRelError (const Matrix33<T> &m, T e) const  | 
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> &  | 
1142  |  | Matrix33<T>::operator += (const Matrix33<T> &v)  | 
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> &  | 
1159  |  | Matrix33<T>::operator += (T a)  | 
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>  | 
1175  |  | Matrix33<T>  | 
1176  |  | Matrix33<T>::operator + (const Matrix33<T> &v) const  | 
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> &  | 
1191  |  | Matrix33<T>::operator -= (const Matrix33<T> &v)  | 
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> &  | 
1208  |  | Matrix33<T>::operator -= (T a)  | 
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>  | 
1224  |  | Matrix33<T>  | 
1225  |  | Matrix33<T>::operator - (const Matrix33<T> &v) const  | 
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>  | 
1239  |  | Matrix33<T>  | 
1240  |  | Matrix33<T>::operator - () const  | 
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> &  | 
1255  |  | Matrix33<T>::negate ()  | 
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> &  | 
1272  |  | Matrix33<T>::operator *= (T a)  | 
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>  | 
1288  |  | Matrix33<T>  | 
1289  |  | Matrix33<T>::operator * (T a) const  | 
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>  | 
1304  |  | operator * (T a, const Matrix33<T> &v)  | 
1305  |  | { | 
1306  |  |     return v * a;  | 
1307  |  | }  | 
1308  |  |  | 
1309  |  | template <class T>  | 
1310  |  | const Matrix33<T> &  | 
1311  |  | Matrix33<T>::operator *= (const Matrix33<T> &v)  | 
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>  | 
1325  |  | Matrix33<T>  | 
1326  |  | Matrix33<T>::operator * (const Matrix33<T> &v) const  | 
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  | 
1341  |  | Matrix33<T>::multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const  | 
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  | 
1356  |  | Matrix33<T>::multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const  | 
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> &  | 
1369  |  | Matrix33<T>::operator /= (T a)  | 
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>  | 
1385  |  | Matrix33<T>  | 
1386  |  | Matrix33<T>::operator / (T a) const  | 
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> &  | 
1401  |  | Matrix33<T>::transpose ()  | 
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>  | 
1417  |  | Matrix33<T>  | 
1418  |  | Matrix33<T>::transposed () const  | 
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>  | 
1440  |  | Matrix33<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>  | 
1552  |  | Matrix33<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  |  |         { | 
1583  |  |             T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / limits<T>::smallest();  | 
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  |  |         { | 
1634  |  |             T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / limits<T>::smallest();  | 
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  | 
1684  |  | Matrix33<T>::determinant () const  | 
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> &  | 
1694  |  | Matrix33<T>::setRotation (S r)  | 
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> &  | 
1719  |  | Matrix33<T>::rotate (S r)  | 
1720  |  | { | 
1721  |  |     *this *= Matrix33<T>().setRotation (r);  | 
1722  |  |     return *this;  | 
1723  |  | }  | 
1724  |  |  | 
1725  |  | template <class T>  | 
1726  |  | const Matrix33<T> &  | 
1727  |  | Matrix33<T>::setScale (T s)  | 
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> &  | 
1740  |  | Matrix33<T>::setScale (const Vec2<S> &s)  | 
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> &  | 
1753  |  | Matrix33<T>::scale (const Vec2<S> &s)  | 
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> &  | 
1769  |  | Matrix33<T>::setTranslation (const Vec2<S> &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>   | 
1788  |  | Matrix33<T>::translation () const  | 
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> &  | 
1796  |  | Matrix33<T>::translate (const Vec2<S> &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> &  | 
1808  |  | Matrix33<T>::setShear (const S &xy)  | 
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> &  | 
1828  |  | Matrix33<T>::setShear (const Vec2<S> &h)  | 
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> &  | 
1866  |  | Matrix33<T>::shear (const Vec2<S> &h)  | 
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 *  | 
1888  |  | Matrix44<T>::operator [] (int i)  | 
1889  | 0  | { | 
1890  | 0  |     return x[i];  | 
1891  | 0  | } Unexecuted instantiation: Imath_2_2::Matrix44<float>::operator[](int) Unexecuted instantiation: Imath_2_2::Matrix44<double>::operator[](int)  | 
1892  |  |  | 
1893  |  | template <class T>  | 
1894  |  | inline const T *  | 
1895  |  | Matrix44<T>::operator [] (int i) const  | 
1896  | 0  | { | 
1897  | 0  |     return x[i];  | 
1898  | 0  | } Unexecuted instantiation: Imath_2_2::Matrix44<float>::operator[](int) const Unexecuted instantiation: Imath_2_2::Matrix44<double>::operator[](int) const  | 
1899  |  |  | 
1900  |  | template <class T>  | 
1901  |  | inline  | 
1902  |  | Matrix44<T>::Matrix44 ()  | 
1903  | 0  | { | 
1904  | 0  |     memset (x, 0, sizeof (x));  | 
1905  | 0  |     x[0][0] = 1;  | 
1906  | 0  |     x[1][1] = 1;  | 
1907  | 0  |     x[2][2] = 1;  | 
1908  | 0  |     x[3][3] = 1;  | 
1909  | 0  | } Unexecuted instantiation: Imath_2_2::Matrix44<double>::Matrix44() Unexecuted instantiation: Imath_2_2::Matrix44<float>::Matrix44()  | 
1910  |  |  | 
1911  |  | template <class T>  | 
1912  |  | inline  | 
1913  |  | Matrix44<T>::Matrix44 (T a)  | 
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  | 0  | { | 
1945  | 0  |     x[0][0] = a;  | 
1946  | 0  |     x[0][1] = b;  | 
1947  | 0  |     x[0][2] = c;  | 
1948  | 0  |     x[0][3] = d;  | 
1949  | 0  |     x[1][0] = e;  | 
1950  | 0  |     x[1][1] = f;  | 
1951  | 0  |     x[1][2] = g;  | 
1952  | 0  |     x[1][3] = h;  | 
1953  | 0  |     x[2][0] = i;  | 
1954  | 0  |     x[2][1] = j;  | 
1955  | 0  |     x[2][2] = k;  | 
1956  | 0  |     x[2][3] = l;  | 
1957  | 0  |     x[3][0] = m;  | 
1958  | 0  |     x[3][1] = n;  | 
1959  | 0  |     x[3][2] = o;  | 
1960  | 0  |     x[3][3] = p;  | 
1961  | 0  | }  | 
1962  |  |  | 
1963  |  |  | 
1964  |  | template <class T>  | 
1965  |  | inline  | 
1966  |  | Matrix44<T>::Matrix44 (Matrix33<T> r, Vec3<T> t)  | 
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  | 
1988  |  | Matrix44<T>::Matrix44 (const Matrix44 &v)  | 
1989  | 0  | { | 
1990  | 0  |     x[0][0] = v.x[0][0];  | 
1991  | 0  |     x[0][1] = v.x[0][1];  | 
1992  | 0  |     x[0][2] = v.x[0][2];  | 
1993  | 0  |     x[0][3] = v.x[0][3];  | 
1994  | 0  |     x[1][0] = v.x[1][0];  | 
1995  | 0  |     x[1][1] = v.x[1][1];  | 
1996  | 0  |     x[1][2] = v.x[1][2];  | 
1997  | 0  |     x[1][3] = v.x[1][3];  | 
1998  | 0  |     x[2][0] = v.x[2][0];  | 
1999  | 0  |     x[2][1] = v.x[2][1];  | 
2000  | 0  |     x[2][2] = v.x[2][2];  | 
2001  | 0  |     x[2][3] = v.x[2][3];  | 
2002  | 0  |     x[3][0] = v.x[3][0];  | 
2003  | 0  |     x[3][1] = v.x[3][1];  | 
2004  | 0  |     x[3][2] = v.x[3][2];  | 
2005  | 0  |     x[3][3] = v.x[3][3];  | 
2006  | 0  | }  | 
2007  |  |  | 
2008  |  | template <class T>  | 
2009  |  | template <class S>  | 
2010  |  | inline  | 
2011  |  | Matrix44<T>::Matrix44 (const Matrix44<S> &v)  | 
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> &  | 
2033  |  | Matrix44<T>::operator = (const Matrix44 &v)  | 
2034  | 0  | { | 
2035  | 0  |     x[0][0] = v.x[0][0];  | 
2036  | 0  |     x[0][1] = v.x[0][1];  | 
2037  | 0  |     x[0][2] = v.x[0][2];  | 
2038  | 0  |     x[0][3] = v.x[0][3];  | 
2039  | 0  |     x[1][0] = v.x[1][0];  | 
2040  | 0  |     x[1][1] = v.x[1][1];  | 
2041  | 0  |     x[1][2] = v.x[1][2];  | 
2042  | 0  |     x[1][3] = v.x[1][3];  | 
2043  | 0  |     x[2][0] = v.x[2][0];  | 
2044  | 0  |     x[2][1] = v.x[2][1];  | 
2045  | 0  |     x[2][2] = v.x[2][2];  | 
2046  | 0  |     x[2][3] = v.x[2][3];  | 
2047  | 0  |     x[3][0] = v.x[3][0];  | 
2048  | 0  |     x[3][1] = v.x[3][1];  | 
2049  | 0  |     x[3][2] = v.x[3][2];  | 
2050  | 0  |     x[3][3] = v.x[3][3];  | 
2051  | 0  |     return *this;  | 
2052  | 0  | } Unexecuted instantiation: Imath_2_2::Matrix44<double>::operator=(Imath_2_2::Matrix44<double> const&) Unexecuted instantiation: Imath_2_2::Matrix44<float>::operator=(Imath_2_2::Matrix44<float> const&)  | 
2053  |  |  | 
2054  |  | template <class T>  | 
2055  |  | inline const Matrix44<T> &  | 
2056  |  | Matrix44<T>::operator = (T a)  | 
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 *  | 
2079  |  | Matrix44<T>::getValue ()  | 
2080  |  | { | 
2081  |  |     return (T *) &x[0][0];  | 
2082  |  | }  | 
2083  |  |  | 
2084  |  | template <class T>  | 
2085  |  | inline const T *  | 
2086  |  | Matrix44<T>::getValue () const  | 
2087  |  | { | 
2088  |  |     return (const T *) &x[0][0];  | 
2089  |  | }  | 
2090  |  |  | 
2091  |  | template <class T>  | 
2092  |  | template <class S>  | 
2093  |  | inline void  | 
2094  |  | Matrix44<T>::getValue (Matrix44<S> &v) const  | 
2095  |  | { | 
2096  |  |     if (isSameType<S,T>::value)  | 
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> &  | 
2124  |  | Matrix44<T>::setValue (const Matrix44<S> &v)  | 
2125  |  | { | 
2126  |  |     if (isSameType<S,T>::value)  | 
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> &  | 
2156  |  | Matrix44<T>::setTheMatrix (const Matrix44<S> &v)  | 
2157  |  | { | 
2158  |  |     if (isSameType<S,T>::value)  | 
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  | 
2187  |  | Matrix44<T>::makeIdentity()  | 
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  | 
2198  |  | Matrix44<T>::operator == (const Matrix44 &v) const  | 
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  | 
2220  |  | Matrix44<T>::operator != (const Matrix44 &v) const  | 
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  | 
2242  |  | Matrix44<T>::equalWithAbsError (const Matrix44<T> &m, T e) const  | 
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  | 
2254  |  | Matrix44<T>::equalWithRelError (const Matrix44<T> &m, T e) const  | 
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> &  | 
2266  |  | Matrix44<T>::operator += (const Matrix44<T> &v)  | 
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> &  | 
2290  |  | Matrix44<T>::operator += (T a)  | 
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>  | 
2313  |  | Matrix44<T>  | 
2314  |  | Matrix44<T>::operator + (const Matrix44<T> &v) const  | 
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> &  | 
2336  |  | Matrix44<T>::operator -= (const Matrix44<T> &v)  | 
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> &  | 
2360  |  | Matrix44<T>::operator -= (T a)  | 
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>  | 
2383  |  | Matrix44<T>  | 
2384  |  | Matrix44<T>::operator - (const Matrix44<T> &v) const  | 
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>  | 
2405  |  | Matrix44<T>  | 
2406  |  | Matrix44<T>::operator - () const  | 
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> &  | 
2428  |  | Matrix44<T>::negate ()  | 
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> &  | 
2452  |  | Matrix44<T>::operator *= (T a)  | 
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>  | 
2475  |  | Matrix44<T>  | 
2476  |  | Matrix44<T>::operator * (T a) const  | 
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>  | 
2498  |  | operator * (T a, const Matrix44<T> &v)  | 
2499  |  | { | 
2500  |  |     return v * a;  | 
2501  |  | }  | 
2502  |  |  | 
2503  |  | template <class T>  | 
2504  |  | inline const Matrix44<T> &  | 
2505  |  | Matrix44<T>::operator *= (const Matrix44<T> &v)  | 
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>  | 
2516  |  | Matrix44<T>::operator * (const Matrix44<T> &v) const  | 
2517  |  | { | 
2518  |  |     Matrix44 tmp (T (0));  | 
2519  |  |  | 
2520  |  |     multiply (*this, v, tmp);  | 
2521  |  |     return tmp;  | 
2522  |  | }  | 
2523  |  |  | 
2524  |  | template <class T>  | 
2525  |  | void  | 
2526  |  | Matrix44<T>::multiply (const Matrix44<T> &a,  | 
2527  |  |                        const Matrix44<T> &b,  | 
2528  |  |                        Matrix44<T> &c)  | 
2529  |  | { | 
2530  |  |     register const T * IMATH_RESTRICT ap = &a.x[0][0];  | 
2531  |  |     register const T * IMATH_RESTRICT bp = &b.x[0][0];  | 
2532  |  |     register       T * IMATH_RESTRICT cp = &c.x[0][0];  | 
2533  |  |  | 
2534  |  |     register 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  | 
2579  |  | Matrix44<T>::multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const  | 
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  | 
2595  |  | Matrix44<T>::multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const  | 
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> &  | 
2610  |  | Matrix44<T>::operator /= (T a)  | 
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>  | 
2633  |  | Matrix44<T>  | 
2634  |  | Matrix44<T>::operator / (T a) const  | 
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> &  | 
2656  |  | Matrix44<T>::transpose ()  | 
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>  | 
2679  |  | Matrix44<T>  | 
2680  |  | Matrix44<T>::transposed () const  | 
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>  | 
2709  |  | Matrix44<T>  | 
2710  |  | Matrix44<T>::gjInverse (bool singExc) const throw (IEX_NAMESPACE::MathExc)  | 
2711  | 0  | { | 
2712  | 0  |     int i, j, k;  | 
2713  | 0  |     Matrix44 s;  | 
2714  | 0  |     Matrix44 t (*this);  | 
2715  |  |  | 
2716  |  |     // Forward elimination  | 
2717  |  | 
  | 
2718  | 0  |     for (i = 0; i < 3 ; i++)  | 
2719  | 0  |     { | 
2720  | 0  |         int pivot = i;  | 
2721  |  | 
  | 
2722  | 0  |         T pivotsize = t[i][i];  | 
2723  |  | 
  | 
2724  | 0  |         if (pivotsize < 0)  | 
2725  | 0  |             pivotsize = -pivotsize;  | 
2726  |  | 
  | 
2727  | 0  |         for (j = i + 1; j < 4; j++)  | 
2728  | 0  |         { | 
2729  | 0  |             T tmp = t[j][i];  | 
2730  |  | 
  | 
2731  | 0  |             if (tmp < 0)  | 
2732  | 0  |                 tmp = -tmp;  | 
2733  |  | 
  | 
2734  | 0  |             if (tmp > pivotsize)  | 
2735  | 0  |             { | 
2736  | 0  |                 pivot = j;  | 
2737  | 0  |                 pivotsize = tmp;  | 
2738  | 0  |             }  | 
2739  | 0  |         }  | 
2740  |  | 
  | 
2741  | 0  |         if (pivotsize == 0)  | 
2742  | 0  |         { | 
2743  | 0  |             if (singExc)  | 
2744  | 0  |                 throw ::IMATH_INTERNAL_NAMESPACE::SingMatrixExc ("Cannot invert singular matrix."); | 
2745  |  |  | 
2746  | 0  |             return Matrix44();  | 
2747  | 0  |         }  | 
2748  |  |  | 
2749  | 0  |         if (pivot != i)  | 
2750  | 0  |         { | 
2751  | 0  |             for (j = 0; j < 4; j++)  | 
2752  | 0  |             { | 
2753  | 0  |                 T tmp;  | 
2754  |  | 
  | 
2755  | 0  |                 tmp = t[i][j];  | 
2756  | 0  |                 t[i][j] = t[pivot][j];  | 
2757  | 0  |                 t[pivot][j] = tmp;  | 
2758  |  | 
  | 
2759  | 0  |                 tmp = s[i][j];  | 
2760  | 0  |                 s[i][j] = s[pivot][j];  | 
2761  | 0  |                 s[pivot][j] = tmp;  | 
2762  | 0  |             }  | 
2763  | 0  |         }  | 
2764  |  | 
  | 
2765  | 0  |         for (j = i + 1; j < 4; j++)  | 
2766  | 0  |         { | 
2767  | 0  |             T f = t[j][i] / t[i][i];  | 
2768  |  | 
  | 
2769  | 0  |             for (k = 0; k < 4; k++)  | 
2770  | 0  |             { | 
2771  | 0  |                 t[j][k] -= f * t[i][k];  | 
2772  | 0  |                 s[j][k] -= f * s[i][k];  | 
2773  | 0  |             }  | 
2774  | 0  |         }  | 
2775  | 0  |     }  | 
2776  |  |  | 
2777  |  |     // Backward substitution  | 
2778  |  |  | 
2779  | 0  |     for (i = 3; i >= 0; --i)  | 
2780  | 0  |     { | 
2781  | 0  |         T f;  | 
2782  |  | 
  | 
2783  | 0  |         if ((f = t[i][i]) == 0)  | 
2784  | 0  |         { | 
2785  | 0  |             if (singExc)  | 
2786  | 0  |                 throw ::IMATH_INTERNAL_NAMESPACE::SingMatrixExc ("Cannot invert singular matrix."); | 
2787  |  |  | 
2788  | 0  |             return Matrix44();  | 
2789  | 0  |         }  | 
2790  |  |  | 
2791  | 0  |         for (j = 0; j < 4; j++)  | 
2792  | 0  |         { | 
2793  | 0  |             t[i][j] /= f;  | 
2794  | 0  |             s[i][j] /= f;  | 
2795  | 0  |         }  | 
2796  |  | 
  | 
2797  | 0  |         for (j = 0; j < i; j++)  | 
2798  | 0  |         { | 
2799  | 0  |             f = t[j][i];  | 
2800  |  | 
  | 
2801  | 0  |             for (k = 0; k < 4; k++)  | 
2802  | 0  |             { | 
2803  | 0  |                 t[j][k] -= f * t[i][k];  | 
2804  | 0  |                 s[j][k] -= f * s[i][k];  | 
2805  | 0  |             }  | 
2806  | 0  |         }  | 
2807  | 0  |     }  | 
2808  |  |  | 
2809  | 0  |     return s;  | 
2810  | 0  | }  | 
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>  | 
2821  |  | Matrix44<T>  | 
2822  |  | Matrix44<T>::inverse (bool singExc) const throw (IEX_NAMESPACE::MathExc)  | 
2823  | 0  | { | 
2824  | 0  |     if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)  | 
2825  | 0  |         return gjInverse(singExc);  | 
2826  |  |  | 
2827  | 0  |     Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],  | 
2828  | 0  |                 x[2][1] * x[0][2] - x[0][1] * x[2][2],  | 
2829  | 0  |                 x[0][1] * x[1][2] - x[1][1] * x[0][2],  | 
2830  | 0  |                 0,  | 
2831  |  | 
  | 
2832  | 0  |                 x[2][0] * x[1][2] - x[1][0] * x[2][2],  | 
2833  | 0  |                 x[0][0] * x[2][2] - x[2][0] * x[0][2],  | 
2834  | 0  |                 x[1][0] * x[0][2] - x[0][0] * x[1][2],  | 
2835  | 0  |                 0,  | 
2836  |  | 
  | 
2837  | 0  |                 x[1][0] * x[2][1] - x[2][0] * x[1][1],  | 
2838  | 0  |                 x[2][0] * x[0][1] - x[0][0] * x[2][1],  | 
2839  | 0  |                 x[0][0] * x[1][1] - x[1][0] * x[0][1],  | 
2840  | 0  |                 0,  | 
2841  |  | 
  | 
2842  | 0  |                 0,  | 
2843  | 0  |                 0,  | 
2844  | 0  |                 0,  | 
2845  | 0  |                 1);  | 
2846  |  | 
  | 
2847  | 0  |     T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];  | 
2848  |  | 
  | 
2849  | 0  |     if (IMATH_INTERNAL_NAMESPACE::abs (r) >= 1)  | 
2850  | 0  |     { | 
2851  | 0  |         for (int i = 0; i < 3; ++i)  | 
2852  | 0  |         { | 
2853  | 0  |             for (int j = 0; j < 3; ++j)  | 
2854  | 0  |             { | 
2855  | 0  |                 s[i][j] /= r;  | 
2856  | 0  |             }  | 
2857  | 0  |         }  | 
2858  | 0  |     }  | 
2859  | 0  |     else  | 
2860  | 0  |     { | 
2861  | 0  |         T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / limits<T>::smallest();  | 
2862  |  | 
  | 
2863  | 0  |         for (int i = 0; i < 3; ++i)  | 
2864  | 0  |         { | 
2865  | 0  |             for (int j = 0; j < 3; ++j)  | 
2866  | 0  |             { | 
2867  | 0  |                 if (mr > IMATH_INTERNAL_NAMESPACE::abs (s[i][j]))  | 
2868  | 0  |                 { | 
2869  | 0  |                     s[i][j] /= r;  | 
2870  | 0  |                 }  | 
2871  | 0  |                 else  | 
2872  | 0  |                 { | 
2873  | 0  |                     if (singExc)  | 
2874  | 0  |                         throw SingMatrixExc ("Cannot invert singular matrix."); | 
2875  |  |  | 
2876  | 0  |                     return Matrix44();  | 
2877  | 0  |                 }  | 
2878  | 0  |             }  | 
2879  | 0  |         }  | 
2880  | 0  |     }  | 
2881  |  |  | 
2882  | 0  |     s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0];  | 
2883  | 0  |     s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1];  | 
2884  | 0  |     s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2];  | 
2885  |  | 
  | 
2886  | 0  |     return s;  | 
2887  | 0  | }  | 
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  | 
2919  |  | Matrix44<T>::determinant () const  | 
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> &  | 
2934  |  | Matrix44<T>::setEulerAngles (const Vec3<S>& r)  | 
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> &  | 
2972  |  | Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle)  | 
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> &  | 
3004  |  | Matrix44<T>::rotate (const Vec3<S> &r)  | 
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> &  | 
3051  |  | Matrix44<T>::setScale (T s)  | 
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> &  | 
3065  |  | Matrix44<T>::setScale (const Vec3<S> &s)  | 
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> &  | 
3079  |  | Matrix44<T>::scale (const Vec3<S> &s)  | 
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> &  | 
3102  |  | Matrix44<T>::setTranslation (const Vec3<S> &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>  | 
3129  |  | Matrix44<T>::translation () const  | 
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> &  | 
3137  |  | Matrix44<T>::translate (const Vec3<S> &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> &  | 
3150  |  | Matrix44<T>::setShear (const Vec3<S> &h)  | 
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> &  | 
3178  |  | Matrix44<T>::setShear (const Shear6<S> &h)  | 
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> &  | 
3206  |  | Matrix44<T>::shear (const Vec3<S> &h)  | 
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> &  | 
3226  |  | Matrix44<T>::shear (const Shear6<S> &h)  | 
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> &  | 
3330  |  | operator *= (Vec2<S> &v, const Matrix33<T> &m)  | 
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>  | 
3344  |  | operator * (const Vec2<S> &v, const Matrix33<T> &m)  | 
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> &  | 
3356  |  | operator *= (Vec3<S> &v, const Matrix33<T> &m)  | 
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>  | 
3371  |  | operator * (const Vec3<S> &v, const Matrix33<T> &m)  | 
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> &  | 
3383  |  | operator *= (Vec3<S> &v, const Matrix44<T> &m)  | 
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>  | 
3399  |  | operator * (const Vec3<S> &v, const Matrix44<T> &m)  | 
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> &  | 
3412  |  | operator *= (Vec4<S> &v, const Matrix44<T> &m)  | 
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>  | 
3429  |  | operator * (const Vec4<S> &v, const Matrix44<T> &m)  | 
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  |  |  | 
3439  |  | IMATH_INTERNAL_NAMESPACE_HEADER_EXIT  | 
3440  |  |  | 
3441  |  | #endif // INCLUDED_IMATHMATRIX_H  |