ImathMatrix.h source code [include/OpenEXR/ImathMatrix.h]

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	{
885	return x[i];
886	}
887
888	template <class T>
889	inline const T *
890	Matrix33<T>::operator [] (int i) const
891	{
892	return x[i];
893	}
894
895	template <class T>
896	inline
897	Matrix33<T>::Matrix33 ()
898	{
899	memset (x, `0`, sizeof (x));
900	x[`0`][`0`] = `1`;
901	x[`1`][`1`] = `1`;
902	x[`2`][`2`] = `1`;
903	}
904
905	template <class T>
906	inline
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	{
969	memcpy (x, v.x, sizeof (x));
970	return *this;
971	}
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	{
1890	return x[i];
1891	}
1892
1893	template <class T>
1894	inline const T *
1895	Matrix44<T>::operator [] (int i) const
1896	{
1897	return x[i];
1898	}
1899
1900	template <class T>
1901	inline
1902	Matrix44<T>::Matrix44 ()
1903	{
1904	memset (x, `0`, sizeof (x));
1905	x[`0`][`0`] = `1`;
1906	x[`1`][`1`] = `1`;
1907	x[`2`][`2`] = `1`;
1908	x[`3`][`3`] = `1`;
1909	}
1910
1911	template <class T>
1912	inline
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	{
1945	x[`0`][`0`] = a;
1946	x[`0`][`1`] = b;
1947	x[`0`][`2`] = c;
1948	x[`0`][`3`] = d;
1949	x[`1`][`0`] = e;
1950	x[`1`][`1`] = f;
1951	x[`1`][`2`] = g;
1952	x[`1`][`3`] = h;
1953	x[`2`][`0`] = i;
1954	x[`2`][`1`] = j;
1955	x[`2`][`2`] = k;
1956	x[`2`][`3`] = l;
1957	x[`3`][`0`] = m;
1958	x[`3`][`1`] = n;
1959	x[`3`][`2`] = o;
1960	x[`3`][`3`] = p;
1961	}
1962
1963
1964	template <class T>
1965	inline
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	{
1990	x[`0`][`0`] = v.x[`0`][`0`];
1991	x[`0`][`1`] = v.x[`0`][`1`];
1992	x[`0`][`2`] = v.x[`0`][`2`];
1993	x[`0`][`3`] = v.x[`0`][`3`];
1994	x[`1`][`0`] = v.x[`1`][`0`];
1995	x[`1`][`1`] = v.x[`1`][`1`];
1996	x[`1`][`2`] = v.x[`1`][`2`];
1997	x[`1`][`3`] = v.x[`1`][`3`];
1998	x[`2`][`0`] = v.x[`2`][`0`];
1999	x[`2`][`1`] = v.x[`2`][`1`];
2000	x[`2`][`2`] = v.x[`2`][`2`];
2001	x[`2`][`3`] = v.x[`2`][`3`];
2002	x[`3`][`0`] = v.x[`3`][`0`];
2003	x[`3`][`1`] = v.x[`3`][`1`];
2004	x[`3`][`2`] = v.x[`3`][`2`];
2005	x[`3`][`3`] = v.x[`3`][`3`];
2006	}
2007
2008	template <class T>
2009	template <class S>
2010	inline
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	{
2035	x[`0`][`0`] = v.x[`0`][`0`];
2036	x[`0`][`1`] = v.x[`0`][`1`];
2037	x[`0`][`2`] = v.x[`0`][`2`];
2038	x[`0`][`3`] = v.x[`0`][`3`];
2039	x[`1`][`0`] = v.x[`1`][`0`];
2040	x[`1`][`1`] = v.x[`1`][`1`];
2041	x[`1`][`2`] = v.x[`1`][`2`];
2042	x[`1`][`3`] = v.x[`1`][`3`];
2043	x[`2`][`0`] = v.x[`2`][`0`];
2044	x[`2`][`1`] = v.x[`2`][`1`];
2045	x[`2`][`2`] = v.x[`2`][`2`];
2046	x[`2`][`3`] = v.x[`2`][`3`];
2047	x[`3`][`0`] = v.x[`3`][`0`];
2048	x[`3`][`1`] = v.x[`3`][`1`];
2049	x[`3`][`2`] = v.x[`3`][`2`];
2050	x[`3`][`3`] = v.x[`3`][`3`];
2051	return *this;
2052	}
2053
2054	template <class T>
2055	inline const Matrix44<T> &
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	{
2712	int i, j, k;
2713	Matrix44 s;
2714	Matrix44 t (*this);
2715
2716	// Forward elimination
2717
2718	for (i = `0`; i < `3` ; i++)
2719	{
2720	int pivot = i;
2721
2722	T pivotsize = t[i][i];
2723
2724	if (pivotsize < `0`)
2725	pivotsize = -pivotsize;
2726
2727	for (j = i + `1`; j < `4`; j++)
2728	{
2729	T tmp = t[j][i];
2730
2731	if (tmp < `0`)
2732	tmp = -tmp;
2733
2734	if (tmp > pivotsize)
2735	{
2736	pivot = j;
2737	pivotsize = tmp;
2738	}
2739	}
2740
2741	if (pivotsize == `0`)
2742	{
2743	if (singExc)
2744	throw ::IMATH_INTERNAL_NAMESPACE::SingMatrixExc ("Cannot invert singular matrix.");
2745
2746	return Matrix44();
2747	}
2748
2749	if (pivot != i)
2750	{
2751	for (j = `0`; j < `4`; j++)
2752	{
2753	T tmp;
2754
2755	tmp = t[i][j];
2756	t[i][j] = t[pivot][j];
2757	t[pivot][j] = tmp;
2758
2759	tmp = s[i][j];
2760	s[i][j] = s[pivot][j];
2761	s[pivot][j] = tmp;
2762	}
2763	}
2764
2765	for (j = i + `1`; j < `4`; j++)
2766	{
2767	T f = t[j][i] / t[i][i];
2768
2769	for (k = `0`; k < `4`; k++)
2770	{
2771	t[j][k] -= f * t[i][k];
2772	s[j][k] -= f * s[i][k];
2773	}
2774	}
2775	}
2776
2777	// Backward substitution
2778
2779	for (i = `3`; i >= `0`; --i)
2780	{
2781	T f;
2782
2783	if ((f = t[i][i]) == `0`)
2784	{
2785	if (singExc)
2786	throw ::IMATH_INTERNAL_NAMESPACE::SingMatrixExc ("Cannot invert singular matrix.");
2787
2788	return Matrix44();
2789	}
2790
2791	for (j = `0`; j < `4`; j++)
2792	{
2793	t[i][j] /= f;
2794	s[i][j] /= f;
2795	}
2796
2797	for (j = `0`; j < i; j++)
2798	{
2799	f = t[j][i];
2800
2801	for (k = `0`; k < `4`; k++)
2802	{
2803	t[j][k] -= f * t[i][k];
2804	s[j][k] -= f * s[i][k];
2805	}
2806	}
2807	}
2808
2809	return s;
2810	}
2811
2812	template <class T>
2813	const Matrix44<T> &
2814	Matrix44<T>::invert (bool singExc) throw (IEX_NAMESPACE::MathExc)
2815	{
2816	*this = inverse (singExc);
2817	return *this;
2818	}
2819
2820	template <class T>
2821	Matrix44<T>
2822	Matrix44<T>::inverse (bool singExc) const throw (IEX_NAMESPACE::MathExc)
2823	{
2824	if (x[`0`][`3`] != `0` \|\| x[`1`][`3`] != `0` \|\| x[`2`][`3`] != `0` \|\| x[`3`][`3`] != `1`)
2825	return gjInverse(singExc);
2826
2827	Matrix44 s (x[`1`][`1`] * x[`2`][`2`] - x[`2`][`1`] * x[`1`][`2`],
2828	x[`2`][`1`] * x[`0`][`2`] - x[`0`][`1`] * x[`2`][`2`],
2829	x[`0`][`1`] * x[`1`][`2`] - x[`1`][`1`] * x[`0`][`2`],
2830	`0`,
2831
2832	x[`2`][`0`] * x[`1`][`2`] - x[`1`][`0`] * x[`2`][`2`],
2833	x[`0`][`0`] * x[`2`][`2`] - x[`2`][`0`] * x[`0`][`2`],
2834	x[`1`][`0`] * x[`0`][`2`] - x[`0`][`0`] * x[`1`][`2`],
2835	`0`,
2836
2837	x[`1`][`0`] * x[`2`][`1`] - x[`2`][`0`] * x[`1`][`1`],
2838	x[`2`][`0`] * x[`0`][`1`] - x[`0`][`0`] * x[`2`][`1`],
2839	x[`0`][`0`] * x[`1`][`1`] - x[`1`][`0`] * x[`0`][`1`],
2840	`0`,
2841
2842	`0`,
2843	`0`,
2844	`0`,
2845	`1`);
2846
2847	T r = x[`0`][`0`] * s[`0`][`0`] + x[`0`][`1`] * s[`1`][`0`] + x[`0`][`2`] * s[`2`][`0`];
2848
2849	if (IMATH_INTERNAL_NAMESPACE::abs (r) >= `1`)
2850	{
2851	for (int i = `0`; i < `3`; ++i)
2852	{
2853	for (int j = `0`; j < `3`; ++j)
2854	{
2855	s[i][j] /= r;
2856	}
2857	}
2858	}
2859	else
2860	{
2861	T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / limits<T>::smallest();
2862
2863	for (int i = `0`; i < `3`; ++i)
2864	{
2865	for (int j = `0`; j < `3`; ++j)
2866	{
2867	if (mr > IMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
2868	{
2869	s[i][j] /= r;
2870	}
2871	else
2872	{
2873	if (singExc)
2874	throw SingMatrixExc ("Cannot invert singular matrix.");
2875
2876	return Matrix44();
2877	}
2878	}
2879	}
2880	}
2881
2882	s[`3`][`0`] = -x[`3`][`0`] * s[`0`][`0`] - x[`3`][`1`] * s[`1`][`0`] - x[`3`][`2`] * s[`2`][`0`];
2883	s[`3`][`1`] = -x[`3`][`0`] * s[`0`][`1`] - x[`3`][`1`] * s[`1`][`1`] - x[`3`][`2`] * s[`2`][`1`];
2884	s[`3`][`2`] = -x[`3`][`0`] * s[`0`][`2`] - x[`3`][`1`] * s[`1`][`2`] - x[`3`][`2`] * s[`2`][`2`];
2885
2886	return s;
2887	}
2888
2889	template <class T>
2890	inline T
2891	Matrix44<T>::fastMinor( const int r0, const int r1, const int r2,
2892	const int c0, const int c1, const int c2) const
2893	{
2894	return x[r0][c0] * (x[r1][c1]x[r2][c2] - x[r1][c2]x[r2][c1])
2895	+ x[r0][c1] * (x[r1][c2]x[r2][c0] - x[r1][c0]x[r2][c2])
2896	+ x[r0][c2] * (x[r1][c0]x[r2][c1] - x[r1][c1]x[r2][c0]);
2897	}
2898
2899	template <class T>
2900	inline T
2901	Matrix44<T>::minorOf (const int r, const int c) const
2902	{
2903	int r0 = `0` + (r < `1` ? `1` : `0`);
2904	int r1 = `1` + (r < `2` ? `1` : `0`);
2905	int r2 = `2` + (r < `3` ? `1` : `0`);
2906	int c0 = `0` + (c < `1` ? `1` : `0`);
2907	int c1 = `1` + (c < `2` ? `1` : `0`);
2908	int c2 = `2` + (c < `3` ? `1` : `0`);
2909
2910	Matrix33<T> working (x[r0][c0],x[r1][c0],x[r2][c0],
2911	x[r0][c1],x[r1][c1],x[r2][c1],
2912	x[r0][c2],x[r1][c2],x[r2][c2]);
2913
2914	return working.determinant();
2915	}
2916
2917	template <class T>
2918	inline T
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
3442

Browse the source code of include/OpenEXR/ImathMatrix.h