1
///////////////////////////////////////////////////////////////////////////
3
// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
6
// All rights reserved.
8
// Redistribution and use in source and binary forms, with or without
9
// modification, are permitted provided that the following conditions are
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
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.
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.
33
///////////////////////////////////////////////////////////////////////////
37
#ifndef INCLUDED_IMATHMATRIX_H
38
#define INCLUDED_IMATHMATRIX_H
40
//----------------------------------------------------------------
42
// 2D (3x3) and 3D (4x4) transformation matrix templates.
44
//----------------------------------------------------------------
46
#include "ImathPlatform.h"
50
#include "ImathShear.h"
55
#if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
56
// suppress exception specification warnings
57
#pragma warning(disable:4290)
64
template <class T> class Matrix33
74
T * operator [] (int i);
75
const T * operator [] (int i) const;
92
Matrix33 (const T a[3][3]);
93
// a[0][0] a[0][1] a[0][2]
94
// a[1][0] a[1][1] a[1][2]
95
// a[2][0] a[2][1] a[2][2]
97
Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i);
104
//--------------------------------
105
// Copy constructor and assignment
106
//--------------------------------
108
Matrix33 (const Matrix33 &v);
110
const Matrix33 & operator = (const Matrix33 &v);
111
const Matrix33 & operator = (T a);
114
//----------------------
115
// Compatibility with Sb
116
//----------------------
119
const T * getValue () const;
122
void getValue (Matrix33<S> &v) const;
124
Matrix33 & setValue (const Matrix33<S> &v);
127
Matrix33 & setTheMatrix (const Matrix33<S> &v);
141
bool operator == (const Matrix33 &v) const;
142
bool operator != (const Matrix33 &v) const;
144
//-----------------------------------------------------------------------
145
// Compare two matrices and test if they are "approximately equal":
147
// equalWithAbsError (m, e)
149
// Returns true if the coefficients of this and m are the same with
150
// an absolute error of no more than e, i.e., for all i, j
152
// abs (this[i][j] - m[i][j]) <= e
154
// equalWithRelError (m, e)
156
// Returns true if the coefficients of this and m are the same with
157
// a relative error of no more than e, i.e., for all i, j
159
// abs (this[i] - v[i][j]) <= e * abs (this[i][j])
160
//-----------------------------------------------------------------------
162
bool equalWithAbsError (const Matrix33<T> &v, T e) const;
163
bool equalWithRelError (const Matrix33<T> &v, T e) const;
166
//------------------------
167
// Component-wise addition
168
//------------------------
170
const Matrix33 & operator += (const Matrix33 &v);
171
const Matrix33 & operator += (T a);
172
Matrix33 operator + (const Matrix33 &v) const;
175
//---------------------------
176
// Component-wise subtraction
177
//---------------------------
179
const Matrix33 & operator -= (const Matrix33 &v);
180
const Matrix33 & operator -= (T a);
181
Matrix33 operator - (const Matrix33 &v) const;
184
//------------------------------------
185
// Component-wise multiplication by -1
186
//------------------------------------
188
Matrix33 operator - () const;
189
const Matrix33 & negate ();
192
//------------------------------
193
// Component-wise multiplication
194
//------------------------------
196
const Matrix33 & operator *= (T a);
197
Matrix33 operator * (T a) const;
200
//-----------------------------------
201
// Matrix-times-matrix multiplication
202
//-----------------------------------
204
const Matrix33 & operator *= (const Matrix33 &v);
205
Matrix33 operator * (const Matrix33 &v) const;
208
//---------------------------------------------
209
// Vector-times-matrix multiplication; see also
210
// the "operator *" functions defined below.
211
//---------------------------------------------
214
void multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const;
217
void multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const;
220
//------------------------
221
// Component-wise division
222
//------------------------
224
const Matrix33 & operator /= (T a);
225
Matrix33 operator / (T a) const;
232
const Matrix33 & transpose ();
233
Matrix33 transposed () const;
236
//------------------------------------------------------------
237
// Inverse matrix: If singExc is false, inverting a singular
238
// matrix produces an identity matrix. If singExc is true,
239
// inverting a singular matrix throws a SingMatrixExc.
241
// inverse() and invert() invert matrices using determinants;
242
// gjInverse() and gjInvert() use the Gauss-Jordan method.
244
// inverse() and invert() are significantly faster than
245
// gjInverse() and gjInvert(), but the results may be slightly
248
//------------------------------------------------------------
250
const Matrix33 & invert (bool singExc = false)
251
throw (Iex::MathExc);
253
Matrix33<T> inverse (bool singExc = false) const
254
throw (Iex::MathExc);
256
const Matrix33 & gjInvert (bool singExc = false)
257
throw (Iex::MathExc);
259
Matrix33<T> gjInverse (bool singExc = false) const
260
throw (Iex::MathExc);
263
//-----------------------------------------
264
// Set matrix to rotation by r (in radians)
265
//-----------------------------------------
268
const Matrix33 & setRotation (S r);
271
//-----------------------------
272
// Rotate the given matrix by r
273
//-----------------------------
276
const Matrix33 & rotate (S r);
279
//--------------------------------------------
280
// Set matrix to scale by given uniform factor
281
//--------------------------------------------
283
const Matrix33 & setScale (T s);
286
//------------------------------------
287
// Set matrix to scale by given vector
288
//------------------------------------
291
const Matrix33 & setScale (const Vec2<S> &s);
294
//----------------------
295
// Scale the matrix by s
296
//----------------------
299
const Matrix33 & scale (const Vec2<S> &s);
302
//------------------------------------------
303
// Set matrix to translation by given vector
304
//------------------------------------------
307
const Matrix33 & setTranslation (const Vec2<S> &t);
310
//-----------------------------
311
// Return translation component
312
//-----------------------------
314
Vec2<T> translation () const;
317
//--------------------------
318
// Translate the matrix by t
319
//--------------------------
322
const Matrix33 & translate (const Vec2<S> &t);
325
//-----------------------------------------------------------
326
// Set matrix to shear x for each y coord. by given factor xy
327
//-----------------------------------------------------------
330
const Matrix33 & setShear (const S &h);
333
//-------------------------------------------------------------
334
// Set matrix to shear x for each y coord. by given factor h[0]
335
// and to shear y for each x coord. by given factor h[1]
336
//-------------------------------------------------------------
339
const Matrix33 & setShear (const Vec2<S> &h);
342
//-----------------------------------------------------------
343
// Shear the matrix in x for each y coord. by given factor xy
344
//-----------------------------------------------------------
347
const Matrix33 & shear (const S &xy);
350
//-----------------------------------------------------------
351
// Shear the matrix in x for each y coord. by given factor xy
352
// and shear y for each x coord. by given factor yx
353
//-----------------------------------------------------------
356
const Matrix33 & shear (const Vec2<S> &h);
359
//-------------------------------------------------
360
// Limitations of type T (see also class limits<T>)
361
//-------------------------------------------------
363
static T baseTypeMin() {return limits<T>::min();}
364
static T baseTypeMax() {return limits<T>::max();}
365
static T baseTypeSmallest() {return limits<T>::smallest();}
366
static T baseTypeEpsilon() {return limits<T>::epsilon();}
370
template <class T> class Matrix44
374
//-------------------
375
// Access to elements
376
//-------------------
380
T * operator [] (int i);
381
const T * operator [] (int i) const;
400
Matrix44 (const T a[4][4]) ;
401
// a[0][0] a[0][1] a[0][2] a[0][3]
402
// a[1][0] a[1][1] a[1][2] a[1][3]
403
// a[2][0] a[2][1] a[2][2] a[2][3]
404
// a[3][0] a[3][1] a[3][2] a[3][3]
406
Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,
407
T i, T j, T k, T l, T m, T n, T o, T p);
414
Matrix44 (Matrix33<T> r, Vec3<T> t);
421
//--------------------------------
422
// Copy constructor and assignment
423
//--------------------------------
425
Matrix44 (const Matrix44 &v);
427
const Matrix44 & operator = (const Matrix44 &v);
428
const Matrix44 & operator = (T a);
431
//----------------------
432
// Compatibility with Sb
433
//----------------------
436
const T * getValue () const;
439
void getValue (Matrix44<S> &v) const;
441
Matrix44 & setValue (const Matrix44<S> &v);
444
Matrix44 & setTheMatrix (const Matrix44<S> &v);
457
bool operator == (const Matrix44 &v) const;
458
bool operator != (const Matrix44 &v) const;
460
//-----------------------------------------------------------------------
461
// Compare two matrices and test if they are "approximately equal":
463
// equalWithAbsError (m, e)
465
// Returns true if the coefficients of this and m are the same with
466
// an absolute error of no more than e, i.e., for all i, j
468
// abs (this[i][j] - m[i][j]) <= e
470
// equalWithRelError (m, e)
472
// Returns true if the coefficients of this and m are the same with
473
// a relative error of no more than e, i.e., for all i, j
475
// abs (this[i] - v[i][j]) <= e * abs (this[i][j])
476
//-----------------------------------------------------------------------
478
bool equalWithAbsError (const Matrix44<T> &v, T e) const;
479
bool equalWithRelError (const Matrix44<T> &v, T e) const;
482
//------------------------
483
// Component-wise addition
484
//------------------------
486
const Matrix44 & operator += (const Matrix44 &v);
487
const Matrix44 & operator += (T a);
488
Matrix44 operator + (const Matrix44 &v) const;
491
//---------------------------
492
// Component-wise subtraction
493
//---------------------------
495
const Matrix44 & operator -= (const Matrix44 &v);
496
const Matrix44 & operator -= (T a);
497
Matrix44 operator - (const Matrix44 &v) const;
500
//------------------------------------
501
// Component-wise multiplication by -1
502
//------------------------------------
504
Matrix44 operator - () const;
505
const Matrix44 & negate ();
508
//------------------------------
509
// Component-wise multiplication
510
//------------------------------
512
const Matrix44 & operator *= (T a);
513
Matrix44 operator * (T a) const;
516
//-----------------------------------
517
// Matrix-times-matrix multiplication
518
//-----------------------------------
520
const Matrix44 & operator *= (const Matrix44 &v);
521
Matrix44 operator * (const Matrix44 &v) const;
523
static void multiply (const Matrix44 &a, // assumes that
524
const Matrix44 &b, // &a != &c and
525
Matrix44 &c); // &b != &c.
528
//---------------------------------------------
529
// Vector-times-matrix multiplication; see also
530
// the "operator *" functions defined below.
531
//---------------------------------------------
534
void multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const;
537
void multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const;
540
//------------------------
541
// Component-wise division
542
//------------------------
544
const Matrix44 & operator /= (T a);
545
Matrix44 operator / (T a) const;
552
const Matrix44 & transpose ();
553
Matrix44 transposed () const;
556
//------------------------------------------------------------
557
// Inverse matrix: If singExc is false, inverting a singular
558
// matrix produces an identity matrix. If singExc is true,
559
// inverting a singular matrix throws a SingMatrixExc.
561
// inverse() and invert() invert matrices using determinants;
562
// gjInverse() and gjInvert() use the Gauss-Jordan method.
564
// inverse() and invert() are significantly faster than
565
// gjInverse() and gjInvert(), but the results may be slightly
568
//------------------------------------------------------------
570
const Matrix44 & invert (bool singExc = false)
571
throw (Iex::MathExc);
573
Matrix44<T> inverse (bool singExc = false) const
574
throw (Iex::MathExc);
576
const Matrix44 & gjInvert (bool singExc = false)
577
throw (Iex::MathExc);
579
Matrix44<T> gjInverse (bool singExc = false) const
580
throw (Iex::MathExc);
583
//--------------------------------------------------------
584
// Set matrix to rotation by XYZ euler angles (in radians)
585
//--------------------------------------------------------
588
const Matrix44 & setEulerAngles (const Vec3<S>& r);
591
//--------------------------------------------------------
592
// Set matrix to rotation around given axis by given angle
593
//--------------------------------------------------------
596
const Matrix44 & setAxisAngle (const Vec3<S>& ax, S ang);
599
//-------------------------------------------
600
// Rotate the matrix by XYZ euler angles in r
601
//-------------------------------------------
604
const Matrix44 & rotate (const Vec3<S> &r);
607
//--------------------------------------------
608
// Set matrix to scale by given uniform factor
609
//--------------------------------------------
611
const Matrix44 & setScale (T s);
614
//------------------------------------
615
// Set matrix to scale by given vector
616
//------------------------------------
619
const Matrix44 & setScale (const Vec3<S> &s);
622
//----------------------
623
// Scale the matrix by s
624
//----------------------
627
const Matrix44 & scale (const Vec3<S> &s);
630
//------------------------------------------
631
// Set matrix to translation by given vector
632
//------------------------------------------
635
const Matrix44 & setTranslation (const Vec3<S> &t);
638
//-----------------------------
639
// Return translation component
640
//-----------------------------
642
const Vec3<T> translation () const;
645
//--------------------------
646
// Translate the matrix by t
647
//--------------------------
650
const Matrix44 & translate (const Vec3<S> &t);
653
//-------------------------------------------------------------
654
// Set matrix to shear by given vector h. The resulting matrix
655
// will shear x for each y coord. by a factor of h[0] ;
656
// will shear x for each z coord. by a factor of h[1] ;
657
// will shear y for each z coord. by a factor of h[2] .
658
//-------------------------------------------------------------
661
const Matrix44 & setShear (const Vec3<S> &h);
664
//------------------------------------------------------------
665
// Set matrix to shear by given factors. The resulting matrix
666
// will shear x for each y coord. by a factor of h.xy ;
667
// will shear x for each z coord. by a factor of h.xz ;
668
// will shear y for each z coord. by a factor of h.yz ;
669
// will shear y for each x coord. by a factor of h.yx ;
670
// will shear z for each x coord. by a factor of h.zx ;
671
// will shear z for each y coord. by a factor of h.zy .
672
//------------------------------------------------------------
675
const Matrix44 & setShear (const Shear6<S> &h);
678
//--------------------------------------------------------
679
// Shear the matrix by given vector. The composed matrix
680
// will be <shear> * <this>, where the shear matrix ...
681
// will shear x for each y coord. by a factor of h[0] ;
682
// will shear x for each z coord. by a factor of h[1] ;
683
// will shear y for each z coord. by a factor of h[2] .
684
//--------------------------------------------------------
687
const Matrix44 & shear (const Vec3<S> &h);
690
//------------------------------------------------------------
691
// Shear the matrix by the given factors. The composed matrix
692
// will be <shear> * <this>, where the shear matrix ...
693
// will shear x for each y coord. by a factor of h.xy ;
694
// will shear x for each z coord. by a factor of h.xz ;
695
// will shear y for each z coord. by a factor of h.yz ;
696
// will shear y for each x coord. by a factor of h.yx ;
697
// will shear z for each x coord. by a factor of h.zx ;
698
// will shear z for each y coord. by a factor of h.zy .
699
//------------------------------------------------------------
702
const Matrix44 & shear (const Shear6<S> &h);
705
//-------------------------------------------------
706
// Limitations of type T (see also class limits<T>)
707
//-------------------------------------------------
709
static T baseTypeMin() {return limits<T>::min();}
710
static T baseTypeMax() {return limits<T>::max();}
711
static T baseTypeSmallest() {return limits<T>::smallest();}
712
static T baseTypeEpsilon() {return limits<T>::epsilon();}
721
std::ostream & operator << (std::ostream & s, const Matrix33<T> &m);
724
std::ostream & operator << (std::ostream & s, const Matrix44<T> &m);
727
//---------------------------------------------
728
// Vector-times-matrix multiplication operators
729
//---------------------------------------------
731
template <class S, class T>
732
const Vec2<S> & operator *= (Vec2<S> &v, const Matrix33<T> &m);
734
template <class S, class T>
735
Vec2<S> operator * (const Vec2<S> &v, const Matrix33<T> &m);
737
template <class S, class T>
738
const Vec3<S> & operator *= (Vec3<S> &v, const Matrix33<T> &m);
740
template <class S, class T>
741
Vec3<S> operator * (const Vec3<S> &v, const Matrix33<T> &m);
743
template <class S, class T>
744
const Vec3<S> & operator *= (Vec3<S> &v, const Matrix44<T> &m);
746
template <class S, class T>
747
Vec3<S> operator * (const Vec3<S> &v, const Matrix44<T> &m);
750
//-------------------------
751
// Typedefs for convenience
752
//-------------------------
754
typedef Matrix33 <float> M33f;
755
typedef Matrix33 <double> M33d;
756
typedef Matrix44 <float> M44f;
757
typedef Matrix44 <double> M44d;
760
//---------------------------
761
// Implementation of Matrix33
762
//---------------------------
766
Matrix33<T>::operator [] (int i)
773
Matrix33<T>::operator [] (int i) const
780
Matrix33<T>::Matrix33 ()
795
Matrix33<T>::Matrix33 (T a)
810
Matrix33<T>::Matrix33 (const T a[3][3])
825
Matrix33<T>::Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i)
840
Matrix33<T>::Matrix33 (const Matrix33 &v)
854
inline const Matrix33<T> &
855
Matrix33<T>::operator = (const Matrix33 &v)
870
inline const Matrix33<T> &
871
Matrix33<T>::operator = (T a)
887
Matrix33<T>::getValue ()
889
return (T *) &x[0][0];
894
Matrix33<T>::getValue () const
896
return (const T *) &x[0][0];
902
Matrix33<T>::getValue (Matrix33<S> &v) const
918
Matrix33<T>::setValue (const Matrix33<S> &v)
935
Matrix33<T>::setTheMatrix (const Matrix33<S> &v)
951
Matrix33<T>::makeIdentity()
966
Matrix33<T>::operator == (const Matrix33 &v) const
968
return x[0][0] == v.x[0][0] &&
969
x[0][1] == v.x[0][1] &&
970
x[0][2] == v.x[0][2] &&
971
x[1][0] == v.x[1][0] &&
972
x[1][1] == v.x[1][1] &&
973
x[1][2] == v.x[1][2] &&
974
x[2][0] == v.x[2][0] &&
975
x[2][1] == v.x[2][1] &&
976
x[2][2] == v.x[2][2];
981
Matrix33<T>::operator != (const Matrix33 &v) const
983
return x[0][0] != v.x[0][0] ||
984
x[0][1] != v.x[0][1] ||
985
x[0][2] != v.x[0][2] ||
986
x[1][0] != v.x[1][0] ||
987
x[1][1] != v.x[1][1] ||
988
x[1][2] != v.x[1][2] ||
989
x[2][0] != v.x[2][0] ||
990
x[2][1] != v.x[2][1] ||
991
x[2][2] != v.x[2][2];
996
Matrix33<T>::equalWithAbsError (const Matrix33<T> &m, T e) const
998
for (int i = 0; i < 3; i++)
999
for (int j = 0; j < 3; j++)
1000
if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e))
1008
Matrix33<T>::equalWithRelError (const Matrix33<T> &m, T e) const
1010
for (int i = 0; i < 3; i++)
1011
for (int j = 0; j < 3; j++)
1012
if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e))
1020
Matrix33<T>::operator += (const Matrix33<T> &v)
1022
x[0][0] += v.x[0][0];
1023
x[0][1] += v.x[0][1];
1024
x[0][2] += v.x[0][2];
1025
x[1][0] += v.x[1][0];
1026
x[1][1] += v.x[1][1];
1027
x[1][2] += v.x[1][2];
1028
x[2][0] += v.x[2][0];
1029
x[2][1] += v.x[2][1];
1030
x[2][2] += v.x[2][2];
1037
Matrix33<T>::operator += (T a)
1054
Matrix33<T>::operator + (const Matrix33<T> &v) const
1056
return Matrix33 (x[0][0] + v.x[0][0],
1057
x[0][1] + v.x[0][1],
1058
x[0][2] + v.x[0][2],
1059
x[1][0] + v.x[1][0],
1060
x[1][1] + v.x[1][1],
1061
x[1][2] + v.x[1][2],
1062
x[2][0] + v.x[2][0],
1063
x[2][1] + v.x[2][1],
1064
x[2][2] + v.x[2][2]);
1069
Matrix33<T>::operator -= (const Matrix33<T> &v)
1071
x[0][0] -= v.x[0][0];
1072
x[0][1] -= v.x[0][1];
1073
x[0][2] -= v.x[0][2];
1074
x[1][0] -= v.x[1][0];
1075
x[1][1] -= v.x[1][1];
1076
x[1][2] -= v.x[1][2];
1077
x[2][0] -= v.x[2][0];
1078
x[2][1] -= v.x[2][1];
1079
x[2][2] -= v.x[2][2];
1086
Matrix33<T>::operator -= (T a)
1103
Matrix33<T>::operator - (const Matrix33<T> &v) const
1105
return Matrix33 (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]);
1118
Matrix33<T>::operator - () const
1120
return Matrix33 (-x[0][0],
1133
Matrix33<T>::negate ()
1150
Matrix33<T>::operator *= (T a)
1167
Matrix33<T>::operator * (T a) const
1169
return Matrix33 (x[0][0] * a,
1182
operator * (T a, const Matrix33<T> &v)
1189
Matrix33<T>::operator *= (const Matrix33<T> &v)
1191
Matrix33 tmp (T (0));
1193
for (int i = 0; i < 3; i++)
1194
for (int j = 0; j < 3; j++)
1195
for (int k = 0; k < 3; k++)
1196
tmp.x[i][j] += x[i][k] * v.x[k][j];
1204
Matrix33<T>::operator * (const Matrix33<T> &v) const
1206
Matrix33 tmp (T (0));
1208
for (int i = 0; i < 3; i++)
1209
for (int j = 0; j < 3; j++)
1210
for (int k = 0; k < 3; k++)
1211
tmp.x[i][j] += x[i][k] * v.x[k][j];
1219
Matrix33<T>::multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const
1223
a = src[0] * x[0][0] + src[1] * x[1][0] + x[2][0];
1224
b = src[0] * x[0][1] + src[1] * x[1][1] + x[2][1];
1225
w = src[0] * x[0][2] + src[1] * x[1][2] + x[2][2];
1234
Matrix33<T>::multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const
1238
a = src[0] * x[0][0] + src[1] * x[1][0];
1239
b = src[0] * x[0][1] + src[1] * x[1][1];
1247
Matrix33<T>::operator /= (T a)
1264
Matrix33<T>::operator / (T a) const
1266
return Matrix33 (x[0][0] / a,
1279
Matrix33<T>::transpose ()
1281
Matrix33 tmp (x[0][0],
1296
Matrix33<T>::transposed () const
1298
return Matrix33 (x[0][0],
1311
Matrix33<T>::gjInvert (bool singExc) throw (Iex::MathExc)
1313
*this = gjInverse (singExc);
1319
Matrix33<T>::gjInverse (bool singExc) const throw (Iex::MathExc)
1325
// Forward elimination
1327
for (i = 0; i < 2 ; i++)
1331
T pivotsize = t[i][i];
1334
pivotsize = -pivotsize;
1336
for (j = i + 1; j < 3; j++)
1343
if (tmp > pivotsize)
1353
throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
1360
for (j = 0; j < 3; j++)
1365
t[i][j] = t[pivot][j];
1369
s[i][j] = s[pivot][j];
1374
for (j = i + 1; j < 3; j++)
1376
T f = t[j][i] / t[i][i];
1378
for (k = 0; k < 3; k++)
1380
t[j][k] -= f * t[i][k];
1381
s[j][k] -= f * s[i][k];
1386
// Backward substitution
1388
for (i = 2; i >= 0; --i)
1392
if ((f = t[i][i]) == 0)
1395
throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
1400
for (j = 0; j < 3; j++)
1406
for (j = 0; j < i; j++)
1410
for (k = 0; k < 3; k++)
1412
t[j][k] -= f * t[i][k];
1413
s[j][k] -= f * s[i][k];
1423
Matrix33<T>::invert (bool singExc) throw (Iex::MathExc)
1425
*this = inverse (singExc);
1431
Matrix33<T>::inverse (bool singExc) const throw (Iex::MathExc)
1433
if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1)
1435
Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
1436
x[2][1] * x[0][2] - x[0][1] * x[2][2],
1437
x[0][1] * x[1][2] - x[1][1] * x[0][2],
1439
x[2][0] * x[1][2] - x[1][0] * x[2][2],
1440
x[0][0] * x[2][2] - x[2][0] * x[0][2],
1441
x[1][0] * x[0][2] - x[0][0] * x[1][2],
1443
x[1][0] * x[2][1] - x[2][0] * x[1][1],
1444
x[2][0] * x[0][1] - x[0][0] * x[2][1],
1445
x[0][0] * x[1][1] - x[1][0] * x[0][1]);
1447
T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
1449
if (Imath::abs (r) >= 1)
1451
for (int i = 0; i < 3; ++i)
1453
for (int j = 0; j < 3; ++j)
1461
T mr = Imath::abs (r) / limits<T>::smallest();
1463
for (int i = 0; i < 3; ++i)
1465
for (int j = 0; j < 3; ++j)
1467
if (mr > Imath::abs (s[i][j]))
1474
throw SingMatrixExc ("Cannot invert "
1475
"singular matrix.");
1486
Matrix33 s ( x[1][1],
1498
T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
1500
if (Imath::abs (r) >= 1)
1502
for (int i = 0; i < 2; ++i)
1504
for (int j = 0; j < 2; ++j)
1512
T mr = Imath::abs (r) / limits<T>::smallest();
1514
for (int i = 0; i < 2; ++i)
1516
for (int j = 0; j < 2; ++j)
1518
if (mr > Imath::abs (s[i][j]))
1525
throw SingMatrixExc ("Cannot invert "
1526
"singular matrix.");
1533
s[2][0] = -x[2][0] * s[0][0] - x[2][1] * s[1][0];
1534
s[2][1] = -x[2][0] * s[0][1] - x[2][1] * s[1][1];
1543
Matrix33<T>::setRotation (S r)
1547
cos_r = Math<T>::cos (r);
1548
sin_r = Math<T>::sin (r);
1568
Matrix33<T>::rotate (S r)
1570
*this *= Matrix33<T>().setRotation (r);
1576
Matrix33<T>::setScale (T s)
1596
Matrix33<T>::setScale (const Vec2<S> &s)
1616
Matrix33<T>::scale (const Vec2<S> &s)
1632
Matrix33<T>::setTranslation (const Vec2<S> &t)
1651
Matrix33<T>::translation () const
1653
return Vec2<T> (x[2][0], x[2][1]);
1659
Matrix33<T>::translate (const Vec2<S> &t)
1661
x[2][0] += t[0] * x[0][0] + t[1] * x[1][0];
1662
x[2][1] += t[0] * x[0][1] + t[1] * x[1][1];
1663
x[2][2] += t[0] * x[0][2] + t[1] * x[1][2];
1671
Matrix33<T>::setShear (const S &xy)
1691
Matrix33<T>::setShear (const Vec2<S> &h)
1711
Matrix33<T>::shear (const S &xy)
1714
// In this case, we don't need a temp. copy of the matrix
1715
// because we never use a value on the RHS after we've
1716
// changed it on the LHS.
1719
x[1][0] += xy * x[0][0];
1720
x[1][1] += xy * x[0][1];
1721
x[1][2] += xy * x[0][2];
1729
Matrix33<T>::shear (const Vec2<S> &h)
1731
Matrix33<T> P (*this);
1733
x[0][0] = P[0][0] + h[1] * P[1][0];
1734
x[0][1] = P[0][1] + h[1] * P[1][1];
1735
x[0][2] = P[0][2] + h[1] * P[1][2];
1737
x[1][0] = P[1][0] + h[0] * P[0][0];
1738
x[1][1] = P[1][1] + h[0] * P[0][1];
1739
x[1][2] = P[1][2] + h[0] * P[0][2];
1745
//---------------------------
1746
// Implementation of Matrix44
1747
//---------------------------
1751
Matrix44<T>::operator [] (int i)
1758
Matrix44<T>::operator [] (int i) const
1765
Matrix44<T>::Matrix44 ()
1787
Matrix44<T>::Matrix44 (T a)
1809
Matrix44<T>::Matrix44 (const T a[4][4])
1831
Matrix44<T>::Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,
1832
T i, T j, T k, T l, T m, T n, T o, T p)
1855
Matrix44<T>::Matrix44 (Matrix33<T> r, Vec3<T> t)
1877
Matrix44<T>::Matrix44 (const Matrix44 &v)
1879
x[0][0] = v.x[0][0];
1880
x[0][1] = v.x[0][1];
1881
x[0][2] = v.x[0][2];
1882
x[0][3] = v.x[0][3];
1883
x[1][0] = v.x[1][0];
1884
x[1][1] = v.x[1][1];
1885
x[1][2] = v.x[1][2];
1886
x[1][3] = v.x[1][3];
1887
x[2][0] = v.x[2][0];
1888
x[2][1] = v.x[2][1];
1889
x[2][2] = v.x[2][2];
1890
x[2][3] = v.x[2][3];
1891
x[3][0] = v.x[3][0];
1892
x[3][1] = v.x[3][1];
1893
x[3][2] = v.x[3][2];
1894
x[3][3] = v.x[3][3];
1898
inline const Matrix44<T> &
1899
Matrix44<T>::operator = (const Matrix44 &v)
1901
x[0][0] = v.x[0][0];
1902
x[0][1] = v.x[0][1];
1903
x[0][2] = v.x[0][2];
1904
x[0][3] = v.x[0][3];
1905
x[1][0] = v.x[1][0];
1906
x[1][1] = v.x[1][1];
1907
x[1][2] = v.x[1][2];
1908
x[1][3] = v.x[1][3];
1909
x[2][0] = v.x[2][0];
1910
x[2][1] = v.x[2][1];
1911
x[2][2] = v.x[2][2];
1912
x[2][3] = v.x[2][3];
1913
x[3][0] = v.x[3][0];
1914
x[3][1] = v.x[3][1];
1915
x[3][2] = v.x[3][2];
1916
x[3][3] = v.x[3][3];
1921
inline const Matrix44<T> &
1922
Matrix44<T>::operator = (T a)
1945
Matrix44<T>::getValue ()
1947
return (T *) &x[0][0];
1952
Matrix44<T>::getValue () const
1954
return (const T *) &x[0][0];
1960
Matrix44<T>::getValue (Matrix44<S> &v) const
1962
v.x[0][0] = x[0][0];
1963
v.x[0][1] = x[0][1];
1964
v.x[0][2] = x[0][2];
1965
v.x[0][3] = x[0][3];
1966
v.x[1][0] = x[1][0];
1967
v.x[1][1] = x[1][1];
1968
v.x[1][2] = x[1][2];
1969
v.x[1][3] = x[1][3];
1970
v.x[2][0] = x[2][0];
1971
v.x[2][1] = x[2][1];
1972
v.x[2][2] = x[2][2];
1973
v.x[2][3] = x[2][3];
1974
v.x[3][0] = x[3][0];
1975
v.x[3][1] = x[3][1];
1976
v.x[3][2] = x[3][2];
1977
v.x[3][3] = x[3][3];
1982
inline Matrix44<T> &
1983
Matrix44<T>::setValue (const Matrix44<S> &v)
1985
x[0][0] = v.x[0][0];
1986
x[0][1] = v.x[0][1];
1987
x[0][2] = v.x[0][2];
1988
x[0][3] = v.x[0][3];
1989
x[1][0] = v.x[1][0];
1990
x[1][1] = v.x[1][1];
1991
x[1][2] = v.x[1][2];
1992
x[1][3] = v.x[1][3];
1993
x[2][0] = v.x[2][0];
1994
x[2][1] = v.x[2][1];
1995
x[2][2] = v.x[2][2];
1996
x[2][3] = v.x[2][3];
1997
x[3][0] = v.x[3][0];
1998
x[3][1] = v.x[3][1];
1999
x[3][2] = v.x[3][2];
2000
x[3][3] = v.x[3][3];
2006
inline Matrix44<T> &
2007
Matrix44<T>::setTheMatrix (const Matrix44<S> &v)
2009
x[0][0] = v.x[0][0];
2010
x[0][1] = v.x[0][1];
2011
x[0][2] = v.x[0][2];
2012
x[0][3] = v.x[0][3];
2013
x[1][0] = v.x[1][0];
2014
x[1][1] = v.x[1][1];
2015
x[1][2] = v.x[1][2];
2016
x[1][3] = v.x[1][3];
2017
x[2][0] = v.x[2][0];
2018
x[2][1] = v.x[2][1];
2019
x[2][2] = v.x[2][2];
2020
x[2][3] = v.x[2][3];
2021
x[3][0] = v.x[3][0];
2022
x[3][1] = v.x[3][1];
2023
x[3][2] = v.x[3][2];
2024
x[3][3] = v.x[3][3];
2030
Matrix44<T>::makeIdentity()
2052
Matrix44<T>::operator == (const Matrix44 &v) const
2054
return x[0][0] == v.x[0][0] &&
2055
x[0][1] == v.x[0][1] &&
2056
x[0][2] == v.x[0][2] &&
2057
x[0][3] == v.x[0][3] &&
2058
x[1][0] == v.x[1][0] &&
2059
x[1][1] == v.x[1][1] &&
2060
x[1][2] == v.x[1][2] &&
2061
x[1][3] == v.x[1][3] &&
2062
x[2][0] == v.x[2][0] &&
2063
x[2][1] == v.x[2][1] &&
2064
x[2][2] == v.x[2][2] &&
2065
x[2][3] == v.x[2][3] &&
2066
x[3][0] == v.x[3][0] &&
2067
x[3][1] == v.x[3][1] &&
2068
x[3][2] == v.x[3][2] &&
2069
x[3][3] == v.x[3][3];
2074
Matrix44<T>::operator != (const Matrix44 &v) const
2076
return x[0][0] != v.x[0][0] ||
2077
x[0][1] != v.x[0][1] ||
2078
x[0][2] != v.x[0][2] ||
2079
x[0][3] != v.x[0][3] ||
2080
x[1][0] != v.x[1][0] ||
2081
x[1][1] != v.x[1][1] ||
2082
x[1][2] != v.x[1][2] ||
2083
x[1][3] != v.x[1][3] ||
2084
x[2][0] != v.x[2][0] ||
2085
x[2][1] != v.x[2][1] ||
2086
x[2][2] != v.x[2][2] ||
2087
x[2][3] != v.x[2][3] ||
2088
x[3][0] != v.x[3][0] ||
2089
x[3][1] != v.x[3][1] ||
2090
x[3][2] != v.x[3][2] ||
2091
x[3][3] != v.x[3][3];
2096
Matrix44<T>::equalWithAbsError (const Matrix44<T> &m, T e) const
2098
for (int i = 0; i < 4; i++)
2099
for (int j = 0; j < 4; j++)
2100
if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e))
2108
Matrix44<T>::equalWithRelError (const Matrix44<T> &m, T e) const
2110
for (int i = 0; i < 4; i++)
2111
for (int j = 0; j < 4; j++)
2112
if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e))
2120
Matrix44<T>::operator += (const Matrix44<T> &v)
2122
x[0][0] += v.x[0][0];
2123
x[0][1] += v.x[0][1];
2124
x[0][2] += v.x[0][2];
2125
x[0][3] += v.x[0][3];
2126
x[1][0] += v.x[1][0];
2127
x[1][1] += v.x[1][1];
2128
x[1][2] += v.x[1][2];
2129
x[1][3] += v.x[1][3];
2130
x[2][0] += v.x[2][0];
2131
x[2][1] += v.x[2][1];
2132
x[2][2] += v.x[2][2];
2133
x[2][3] += v.x[2][3];
2134
x[3][0] += v.x[3][0];
2135
x[3][1] += v.x[3][1];
2136
x[3][2] += v.x[3][2];
2137
x[3][3] += v.x[3][3];
2144
Matrix44<T>::operator += (T a)
2168
Matrix44<T>::operator + (const Matrix44<T> &v) const
2170
return Matrix44 (x[0][0] + v.x[0][0],
2171
x[0][1] + v.x[0][1],
2172
x[0][2] + v.x[0][2],
2173
x[0][3] + v.x[0][3],
2174
x[1][0] + v.x[1][0],
2175
x[1][1] + v.x[1][1],
2176
x[1][2] + v.x[1][2],
2177
x[1][3] + v.x[1][3],
2178
x[2][0] + v.x[2][0],
2179
x[2][1] + v.x[2][1],
2180
x[2][2] + v.x[2][2],
2181
x[2][3] + v.x[2][3],
2182
x[3][0] + v.x[3][0],
2183
x[3][1] + v.x[3][1],
2184
x[3][2] + v.x[3][2],
2185
x[3][3] + v.x[3][3]);
2190
Matrix44<T>::operator -= (const Matrix44<T> &v)
2192
x[0][0] -= v.x[0][0];
2193
x[0][1] -= v.x[0][1];
2194
x[0][2] -= v.x[0][2];
2195
x[0][3] -= v.x[0][3];
2196
x[1][0] -= v.x[1][0];
2197
x[1][1] -= v.x[1][1];
2198
x[1][2] -= v.x[1][2];
2199
x[1][3] -= v.x[1][3];
2200
x[2][0] -= v.x[2][0];
2201
x[2][1] -= v.x[2][1];
2202
x[2][2] -= v.x[2][2];
2203
x[2][3] -= v.x[2][3];
2204
x[3][0] -= v.x[3][0];
2205
x[3][1] -= v.x[3][1];
2206
x[3][2] -= v.x[3][2];
2207
x[3][3] -= v.x[3][3];
2214
Matrix44<T>::operator -= (T a)
2238
Matrix44<T>::operator - (const Matrix44<T> &v) const
2240
return Matrix44 (x[0][0] - v.x[0][0],
2241
x[0][1] - v.x[0][1],
2242
x[0][2] - v.x[0][2],
2243
x[0][3] - v.x[0][3],
2244
x[1][0] - v.x[1][0],
2245
x[1][1] - v.x[1][1],
2246
x[1][2] - v.x[1][2],
2247
x[1][3] - v.x[1][3],
2248
x[2][0] - v.x[2][0],
2249
x[2][1] - v.x[2][1],
2250
x[2][2] - v.x[2][2],
2251
x[2][3] - v.x[2][3],
2252
x[3][0] - v.x[3][0],
2253
x[3][1] - v.x[3][1],
2254
x[3][2] - v.x[3][2],
2255
x[3][3] - v.x[3][3]);
2260
Matrix44<T>::operator - () const
2262
return Matrix44 (-x[0][0],
2282
Matrix44<T>::negate ()
2306
Matrix44<T>::operator *= (T a)
2330
Matrix44<T>::operator * (T a) const
2332
return Matrix44 (x[0][0] * a,
2352
operator * (T a, const Matrix44<T> &v)
2358
inline const Matrix44<T> &
2359
Matrix44<T>::operator *= (const Matrix44<T> &v)
2361
Matrix44 tmp (T (0));
2363
multiply (*this, v, tmp);
2370
Matrix44<T>::operator * (const Matrix44<T> &v) const
2372
Matrix44 tmp (T (0));
2374
multiply (*this, v, tmp);
2380
Matrix44<T>::multiply (const Matrix44<T> &a,
2381
const Matrix44<T> &b,
2384
register const T * restrict ap = &a.x[0][0];
2385
register const T * restrict bp = &b.x[0][0];
2386
register T * restrict cp = &c.x[0][0];
2388
register T a0, a1, a2, a3;
2395
cp[0] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2396
cp[1] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2397
cp[2] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2398
cp[3] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2405
cp[4] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2406
cp[5] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2407
cp[6] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2408
cp[7] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2415
cp[8] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2416
cp[9] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2417
cp[10] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2418
cp[11] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2425
cp[12] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2426
cp[13] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2427
cp[14] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2428
cp[15] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2431
template <class T> template <class S>
2433
Matrix44<T>::multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const
2437
a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0] + x[3][0];
2438
b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1] + x[3][1];
2439
c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2] + x[3][2];
2440
w = src[0] * x[0][3] + src[1] * x[1][3] + src[2] * x[2][3] + x[3][3];
2447
template <class T> template <class S>
2449
Matrix44<T>::multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const
2453
a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0];
2454
b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1];
2455
c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2];
2464
Matrix44<T>::operator /= (T a)
2488
Matrix44<T>::operator / (T a) const
2490
return Matrix44 (x[0][0] / a,
2510
Matrix44<T>::transpose ()
2512
Matrix44 tmp (x[0][0],
2534
Matrix44<T>::transposed () const
2536
return Matrix44 (x[0][0],
2556
Matrix44<T>::gjInvert (bool singExc) throw (Iex::MathExc)
2558
*this = gjInverse (singExc);
2564
Matrix44<T>::gjInverse (bool singExc) const throw (Iex::MathExc)
2570
// Forward elimination
2572
for (i = 0; i < 3 ; i++)
2576
T pivotsize = t[i][i];
2579
pivotsize = -pivotsize;
2581
for (j = i + 1; j < 4; j++)
2588
if (tmp > pivotsize)
2598
throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
2605
for (j = 0; j < 4; j++)
2610
t[i][j] = t[pivot][j];
2614
s[i][j] = s[pivot][j];
2619
for (j = i + 1; j < 4; j++)
2621
T f = t[j][i] / t[i][i];
2623
for (k = 0; k < 4; k++)
2625
t[j][k] -= f * t[i][k];
2626
s[j][k] -= f * s[i][k];
2631
// Backward substitution
2633
for (i = 3; i >= 0; --i)
2637
if ((f = t[i][i]) == 0)
2640
throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
2645
for (j = 0; j < 4; j++)
2651
for (j = 0; j < i; j++)
2655
for (k = 0; k < 4; k++)
2657
t[j][k] -= f * t[i][k];
2658
s[j][k] -= f * s[i][k];
2668
Matrix44<T>::invert (bool singExc) throw (Iex::MathExc)
2670
*this = inverse (singExc);
2676
Matrix44<T>::inverse (bool singExc) const throw (Iex::MathExc)
2678
if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)
2679
return gjInverse(singExc);
2681
Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
2682
x[2][1] * x[0][2] - x[0][1] * x[2][2],
2683
x[0][1] * x[1][2] - x[1][1] * x[0][2],
2686
x[2][0] * x[1][2] - x[1][0] * x[2][2],
2687
x[0][0] * x[2][2] - x[2][0] * x[0][2],
2688
x[1][0] * x[0][2] - x[0][0] * x[1][2],
2691
x[1][0] * x[2][1] - x[2][0] * x[1][1],
2692
x[2][0] * x[0][1] - x[0][0] * x[2][1],
2693
x[0][0] * x[1][1] - x[1][0] * x[0][1],
2701
T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
2703
if (Imath::abs (r) >= 1)
2705
for (int i = 0; i < 3; ++i)
2707
for (int j = 0; j < 3; ++j)
2715
T mr = Imath::abs (r) / limits<T>::smallest();
2717
for (int i = 0; i < 3; ++i)
2719
for (int j = 0; j < 3; ++j)
2721
if (mr > Imath::abs (s[i][j]))
2728
throw SingMatrixExc ("Cannot invert singular matrix.");
2736
s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0];
2737
s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1];
2738
s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2];
2746
Matrix44<T>::setEulerAngles (const Vec3<S>& r)
2748
S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
2750
cos_rz = Math<T>::cos (r[2]);
2751
cos_ry = Math<T>::cos (r[1]);
2752
cos_rx = Math<T>::cos (r[0]);
2754
sin_rz = Math<T>::sin (r[2]);
2755
sin_ry = Math<T>::sin (r[1]);
2756
sin_rx = Math<T>::sin (r[0]);
2758
x[0][0] = cos_rz * cos_ry;
2759
x[0][1] = sin_rz * cos_ry;
2763
x[1][0] = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
2764
x[1][1] = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
2765
x[1][2] = cos_ry * sin_rx;
2768
x[2][0] = sin_rz * sin_rx + cos_rz * sin_ry * cos_rx;
2769
x[2][1] = -cos_rz * sin_rx + sin_rz * sin_ry * cos_rx;
2770
x[2][2] = cos_ry * cos_rx;
2784
Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle)
2786
Vec3<S> unit (axis.normalized());
2787
S sine = Math<T>::sin (angle);
2788
S cosine = Math<T>::cos (angle);
2790
x[0][0] = unit[0] * unit[0] * (1 - cosine) + cosine;
2791
x[0][1] = unit[0] * unit[1] * (1 - cosine) + unit[2] * sine;
2792
x[0][2] = unit[0] * unit[2] * (1 - cosine) - unit[1] * sine;
2795
x[1][0] = unit[0] * unit[1] * (1 - cosine) - unit[2] * sine;
2796
x[1][1] = unit[1] * unit[1] * (1 - cosine) + cosine;
2797
x[1][2] = unit[1] * unit[2] * (1 - cosine) + unit[0] * sine;
2800
x[2][0] = unit[0] * unit[2] * (1 - cosine) + unit[1] * sine;
2801
x[2][1] = unit[1] * unit[2] * (1 - cosine) - unit[0] * sine;
2802
x[2][2] = unit[2] * unit[2] * (1 - cosine) + cosine;
2816
Matrix44<T>::rotate (const Vec3<S> &r)
2818
S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
2823
cos_rz = Math<S>::cos (r[2]);
2824
cos_ry = Math<S>::cos (r[1]);
2825
cos_rx = Math<S>::cos (r[0]);
2827
sin_rz = Math<S>::sin (r[2]);
2828
sin_ry = Math<S>::sin (r[1]);
2829
sin_rx = Math<S>::sin (r[0]);
2831
m00 = cos_rz * cos_ry;
2832
m01 = sin_rz * cos_ry;
2834
m10 = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
2835
m11 = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
2836
m12 = cos_ry * sin_rx;
2837
m20 = -sin_rz * -sin_rx + cos_rz * sin_ry * cos_rx;
2838
m21 = cos_rz * -sin_rx + sin_rz * sin_ry * cos_rx;
2839
m22 = cos_ry * cos_rx;
2841
Matrix44<T> P (*this);
2843
x[0][0] = P[0][0] * m00 + P[1][0] * m01 + P[2][0] * m02;
2844
x[0][1] = P[0][1] * m00 + P[1][1] * m01 + P[2][1] * m02;
2845
x[0][2] = P[0][2] * m00 + P[1][2] * m01 + P[2][2] * m02;
2846
x[0][3] = P[0][3] * m00 + P[1][3] * m01 + P[2][3] * m02;
2848
x[1][0] = P[0][0] * m10 + P[1][0] * m11 + P[2][0] * m12;
2849
x[1][1] = P[0][1] * m10 + P[1][1] * m11 + P[2][1] * m12;
2850
x[1][2] = P[0][2] * m10 + P[1][2] * m11 + P[2][2] * m12;
2851
x[1][3] = P[0][3] * m10 + P[1][3] * m11 + P[2][3] * m12;
2853
x[2][0] = P[0][0] * m20 + P[1][0] * m21 + P[2][0] * m22;
2854
x[2][1] = P[0][1] * m20 + P[1][1] * m21 + P[2][1] * m22;
2855
x[2][2] = P[0][2] * m20 + P[1][2] * m21 + P[2][2] * m22;
2856
x[2][3] = P[0][3] * m20 + P[1][3] * m21 + P[2][3] * m22;
2863
Matrix44<T>::setScale (T s)
2891
Matrix44<T>::setScale (const Vec3<S> &s)
2919
Matrix44<T>::scale (const Vec3<S> &s)
2942
Matrix44<T>::setTranslation (const Vec3<S> &t)
2968
inline const Vec3<T>
2969
Matrix44<T>::translation () const
2971
return Vec3<T> (x[3][0], x[3][1], x[3][2]);
2977
Matrix44<T>::translate (const Vec3<S> &t)
2979
x[3][0] += t[0] * x[0][0] + t[1] * x[1][0] + t[2] * x[2][0];
2980
x[3][1] += t[0] * x[0][1] + t[1] * x[1][1] + t[2] * x[2][1];
2981
x[3][2] += t[0] * x[0][2] + t[1] * x[1][2] + t[2] * x[2][2];
2982
x[3][3] += t[0] * x[0][3] + t[1] * x[1][3] + t[2] * x[2][3];
2990
Matrix44<T>::setShear (const Vec3<S> &h)
3018
Matrix44<T>::setShear (const Shear6<S> &h)
3046
Matrix44<T>::shear (const Vec3<S> &h)
3049
// In this case, we don't need a temp. copy of the matrix
3050
// because we never use a value on the RHS after we've
3051
// changed it on the LHS.
3054
for (int i=0; i < 4; i++)
3056
x[2][i] += h[1] * x[0][i] + h[2] * x[1][i];
3057
x[1][i] += h[0] * x[0][i];
3066
Matrix44<T>::shear (const Shear6<S> &h)
3068
Matrix44<T> P (*this);
3070
for (int i=0; i < 4; i++)
3072
x[0][i] = P[0][i] + h.yx * P[1][i] + h.zx * P[2][i];
3073
x[1][i] = h.xy * P[0][i] + P[1][i] + h.zy * P[2][i];
3074
x[2][i] = h.xz * P[0][i] + h.yz * P[1][i] + P[2][i];
3081
//--------------------------------
3082
// Implementation of stream output
3083
//--------------------------------
3087
operator << (std::ostream &s, const Matrix33<T> &m)
3089
std::ios_base::fmtflags oldFlags = s.flags();
3092
if (s.flags() & std::ios_base::fixed)
3094
s.setf (std::ios_base::showpoint);
3095
width = s.precision() + 5;
3099
s.setf (std::ios_base::scientific);
3100
s.setf (std::ios_base::showpoint);
3101
width = s.precision() + 8;
3104
s << "(" << std::setw (width) << m[0][0] <<
3105
" " << std::setw (width) << m[0][1] <<
3106
" " << std::setw (width) << m[0][2] << "\n" <<
3108
" " << std::setw (width) << m[1][0] <<
3109
" " << std::setw (width) << m[1][1] <<
3110
" " << std::setw (width) << m[1][2] << "\n" <<
3112
" " << std::setw (width) << m[2][0] <<
3113
" " << std::setw (width) << m[2][1] <<
3114
" " << std::setw (width) << m[2][2] << ")\n";
3122
operator << (std::ostream &s, const Matrix44<T> &m)
3124
std::ios_base::fmtflags oldFlags = s.flags();
3127
if (s.flags() & std::ios_base::fixed)
3129
s.setf (std::ios_base::showpoint);
3130
width = s.precision() + 5;
3134
s.setf (std::ios_base::scientific);
3135
s.setf (std::ios_base::showpoint);
3136
width = s.precision() + 8;
3139
s << "(" << std::setw (width) << m[0][0] <<
3140
" " << std::setw (width) << m[0][1] <<
3141
" " << std::setw (width) << m[0][2] <<
3142
" " << std::setw (width) << m[0][3] << "\n" <<
3144
" " << std::setw (width) << m[1][0] <<
3145
" " << std::setw (width) << m[1][1] <<
3146
" " << std::setw (width) << m[1][2] <<
3147
" " << std::setw (width) << m[1][3] << "\n" <<
3149
" " << std::setw (width) << m[2][0] <<
3150
" " << std::setw (width) << m[2][1] <<
3151
" " << std::setw (width) << m[2][2] <<
3152
" " << std::setw (width) << m[2][3] << "\n" <<
3154
" " << std::setw (width) << m[3][0] <<
3155
" " << std::setw (width) << m[3][1] <<
3156
" " << std::setw (width) << m[3][2] <<
3157
" " << std::setw (width) << m[3][3] << ")\n";
3164
//---------------------------------------------------------------
3165
// Implementation of vector-times-matrix multiplication operators
3166
//---------------------------------------------------------------
3168
template <class S, class T>
3169
inline const Vec2<S> &
3170
operator *= (Vec2<S> &v, const Matrix33<T> &m)
3172
S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
3173
S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
3174
S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
3182
template <class S, class T>
3184
operator * (const Vec2<S> &v, const Matrix33<T> &m)
3186
S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
3187
S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
3188
S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
3190
return Vec2<S> (x / w, y / w);
3194
template <class S, class T>
3195
inline const Vec3<S> &
3196
operator *= (Vec3<S> &v, const Matrix33<T> &m)
3198
S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
3199
S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
3200
S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
3210
template <class S, class T>
3212
operator * (const Vec3<S> &v, const Matrix33<T> &m)
3214
S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
3215
S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
3216
S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
3218
return Vec3<S> (x, y, z);
3222
template <class S, class T>
3223
inline const Vec3<S> &
3224
operator *= (Vec3<S> &v, const Matrix44<T> &m)
3226
S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
3227
S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
3228
S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
3229
S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
3238
template <class S, class T>
3240
operator * (const Vec3<S> &v, const Matrix44<T> &m)
3242
S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
3243
S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
3244
S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
3245
S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
3247
return Vec3<S> (x / w, y / w, z / w);
3250
} // namespace Imath
1
///////////////////////////////////////////////////////////////////////////
3
// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
6
// All rights reserved.
8
// Redistribution and use in source and binary forms, with or without
9
// modification, are permitted provided that the following conditions are
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
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.
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.
33
///////////////////////////////////////////////////////////////////////////
37
#ifndef INCLUDED_IMATHMATRIX_H
38
#define INCLUDED_IMATHMATRIX_H
40
//----------------------------------------------------------------
42
// 2D (3x3) and 3D (4x4) transformation matrix templates.
44
//----------------------------------------------------------------
46
#include "ImathPlatform.h"
50
#include "ImathShear.h"
55
#if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
56
// suppress exception specification warnings
57
#pragma warning(disable:4290)
64
template <class T> class Matrix33
74
T * operator [] (int i);
75
const T * operator [] (int i) const;
92
Matrix33 (const T a[3][3]);
93
// a[0][0] a[0][1] a[0][2]
94
// a[1][0] a[1][1] a[1][2]
95
// a[2][0] a[2][1] a[2][2]
97
Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i);
104
//--------------------------------
105
// Copy constructor and assignment
106
//--------------------------------
108
Matrix33 (const Matrix33 &v);
110
const Matrix33 & operator = (const Matrix33 &v);
111
const Matrix33 & operator = (T a);
114
//----------------------
115
// Compatibility with Sb
116
//----------------------
119
const T * getValue () const;
122
void getValue (Matrix33<S> &v) const;
124
Matrix33 & setValue (const Matrix33<S> &v);
127
Matrix33 & setTheMatrix (const Matrix33<S> &v);
141
bool operator == (const Matrix33 &v) const;
142
bool operator != (const Matrix33 &v) const;
144
//-----------------------------------------------------------------------
145
// Compare two matrices and test if they are "approximately equal":
147
// equalWithAbsError (m, e)
149
// Returns true if the coefficients of this and m are the same with
150
// an absolute error of no more than e, i.e., for all i, j
152
// abs (this[i][j] - m[i][j]) <= e
154
// equalWithRelError (m, e)
156
// Returns true if the coefficients of this and m are the same with
157
// a relative error of no more than e, i.e., for all i, j
159
// abs (this[i] - v[i][j]) <= e * abs (this[i][j])
160
//-----------------------------------------------------------------------
162
bool equalWithAbsError (const Matrix33<T> &v, T e) const;
163
bool equalWithRelError (const Matrix33<T> &v, T e) const;
166
//------------------------
167
// Component-wise addition
168
//------------------------
170
const Matrix33 & operator += (const Matrix33 &v);
171
const Matrix33 & operator += (T a);
172
Matrix33 operator + (const Matrix33 &v) const;
175
//---------------------------
176
// Component-wise subtraction
177
//---------------------------
179
const Matrix33 & operator -= (const Matrix33 &v);
180
const Matrix33 & operator -= (T a);
181
Matrix33 operator - (const Matrix33 &v) const;
184
//------------------------------------
185
// Component-wise multiplication by -1
186
//------------------------------------
188
Matrix33 operator - () const;
189
const Matrix33 & negate ();
192
//------------------------------
193
// Component-wise multiplication
194
//------------------------------
196
const Matrix33 & operator *= (T a);
197
Matrix33 operator * (T a) const;
200
//-----------------------------------
201
// Matrix-times-matrix multiplication
202
//-----------------------------------
204
const Matrix33 & operator *= (const Matrix33 &v);
205
Matrix33 operator * (const Matrix33 &v) const;
208
//---------------------------------------------
209
// Vector-times-matrix multiplication; see also
210
// the "operator *" functions defined below.
211
//---------------------------------------------
214
void multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const;
217
void multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const;
220
//------------------------
221
// Component-wise division
222
//------------------------
224
const Matrix33 & operator /= (T a);
225
Matrix33 operator / (T a) const;
232
const Matrix33 & transpose ();
233
Matrix33 transposed () const;
236
//------------------------------------------------------------
237
// Inverse matrix: If singExc is false, inverting a singular
238
// matrix produces an identity matrix. If singExc is true,
239
// inverting a singular matrix throws a SingMatrixExc.
241
// inverse() and invert() invert matrices using determinants;
242
// gjInverse() and gjInvert() use the Gauss-Jordan method.
244
// inverse() and invert() are significantly faster than
245
// gjInverse() and gjInvert(), but the results may be slightly
248
//------------------------------------------------------------
250
const Matrix33 & invert (bool singExc = false)
251
throw (Iex::MathExc);
253
Matrix33<T> inverse (bool singExc = false) const
254
throw (Iex::MathExc);
256
const Matrix33 & gjInvert (bool singExc = false)
257
throw (Iex::MathExc);
259
Matrix33<T> gjInverse (bool singExc = false) const
260
throw (Iex::MathExc);
263
//-----------------------------------------
264
// Set matrix to rotation by r (in radians)
265
//-----------------------------------------
268
const Matrix33 & setRotation (S r);
271
//-----------------------------
272
// Rotate the given matrix by r
273
//-----------------------------
276
const Matrix33 & rotate (S r);
279
//--------------------------------------------
280
// Set matrix to scale by given uniform factor
281
//--------------------------------------------
283
const Matrix33 & setScale (T s);
286
//------------------------------------
287
// Set matrix to scale by given vector
288
//------------------------------------
291
const Matrix33 & setScale (const Vec2<S> &s);
294
//----------------------
295
// Scale the matrix by s
296
//----------------------
299
const Matrix33 & scale (const Vec2<S> &s);
302
//------------------------------------------
303
// Set matrix to translation by given vector
304
//------------------------------------------
307
const Matrix33 & setTranslation (const Vec2<S> &t);
310
//-----------------------------
311
// Return translation component
312
//-----------------------------
314
Vec2<T> translation () const;
317
//--------------------------
318
// Translate the matrix by t
319
//--------------------------
322
const Matrix33 & translate (const Vec2<S> &t);
325
//-----------------------------------------------------------
326
// Set matrix to shear x for each y coord. by given factor xy
327
//-----------------------------------------------------------
330
const Matrix33 & setShear (const S &h);
333
//-------------------------------------------------------------
334
// Set matrix to shear x for each y coord. by given factor h[0]
335
// and to shear y for each x coord. by given factor h[1]
336
//-------------------------------------------------------------
339
const Matrix33 & setShear (const Vec2<S> &h);
342
//-----------------------------------------------------------
343
// Shear the matrix in x for each y coord. by given factor xy
344
//-----------------------------------------------------------
347
const Matrix33 & shear (const S &xy);
350
//-----------------------------------------------------------
351
// Shear the matrix in x for each y coord. by given factor xy
352
// and shear y for each x coord. by given factor yx
353
//-----------------------------------------------------------
356
const Matrix33 & shear (const Vec2<S> &h);
359
//-------------------------------------------------
360
// Limitations of type T (see also class limits<T>)
361
//-------------------------------------------------
363
static T baseTypeMin() {return limits<T>::min();}
364
static T baseTypeMax() {return limits<T>::max();}
365
static T baseTypeSmallest() {return limits<T>::smallest();}
366
static T baseTypeEpsilon() {return limits<T>::epsilon();}
370
template <class T> class Matrix44
374
//-------------------
375
// Access to elements
376
//-------------------
380
T * operator [] (int i);
381
const T * operator [] (int i) const;
400
Matrix44 (const T a[4][4]) ;
401
// a[0][0] a[0][1] a[0][2] a[0][3]
402
// a[1][0] a[1][1] a[1][2] a[1][3]
403
// a[2][0] a[2][1] a[2][2] a[2][3]
404
// a[3][0] a[3][1] a[3][2] a[3][3]
406
Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,
407
T i, T j, T k, T l, T m, T n, T o, T p);
414
Matrix44 (Matrix33<T> r, Vec3<T> t);
421
//--------------------------------
422
// Copy constructor and assignment
423
//--------------------------------
425
Matrix44 (const Matrix44 &v);
427
const Matrix44 & operator = (const Matrix44 &v);
428
const Matrix44 & operator = (T a);
431
//----------------------
432
// Compatibility with Sb
433
//----------------------
436
const T * getValue () const;
439
void getValue (Matrix44<S> &v) const;
441
Matrix44 & setValue (const Matrix44<S> &v);
444
Matrix44 & setTheMatrix (const Matrix44<S> &v);
457
bool operator == (const Matrix44 &v) const;
458
bool operator != (const Matrix44 &v) const;
460
//-----------------------------------------------------------------------
461
// Compare two matrices and test if they are "approximately equal":
463
// equalWithAbsError (m, e)
465
// Returns true if the coefficients of this and m are the same with
466
// an absolute error of no more than e, i.e., for all i, j
468
// abs (this[i][j] - m[i][j]) <= e
470
// equalWithRelError (m, e)
472
// Returns true if the coefficients of this and m are the same with
473
// a relative error of no more than e, i.e., for all i, j
475
// abs (this[i] - v[i][j]) <= e * abs (this[i][j])
476
//-----------------------------------------------------------------------
478
bool equalWithAbsError (const Matrix44<T> &v, T e) const;
479
bool equalWithRelError (const Matrix44<T> &v, T e) const;
482
//------------------------
483
// Component-wise addition
484
//------------------------
486
const Matrix44 & operator += (const Matrix44 &v);
487
const Matrix44 & operator += (T a);
488
Matrix44 operator + (const Matrix44 &v) const;
491
//---------------------------
492
// Component-wise subtraction
493
//---------------------------
495
const Matrix44 & operator -= (const Matrix44 &v);
496
const Matrix44 & operator -= (T a);
497
Matrix44 operator - (const Matrix44 &v) const;
500
//------------------------------------
501
// Component-wise multiplication by -1
502
//------------------------------------
504
Matrix44 operator - () const;
505
const Matrix44 & negate ();
508
//------------------------------
509
// Component-wise multiplication
510
//------------------------------
512
const Matrix44 & operator *= (T a);
513
Matrix44 operator * (T a) const;
516
//-----------------------------------
517
// Matrix-times-matrix multiplication
518
//-----------------------------------
520
const Matrix44 & operator *= (const Matrix44 &v);
521
Matrix44 operator * (const Matrix44 &v) const;
523
static void multiply (const Matrix44 &a, // assumes that
524
const Matrix44 &b, // &a != &c and
525
Matrix44 &c); // &b != &c.
528
//---------------------------------------------
529
// Vector-times-matrix multiplication; see also
530
// the "operator *" functions defined below.
531
//---------------------------------------------
534
void multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const;
537
void multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const;
540
//------------------------
541
// Component-wise division
542
//------------------------
544
const Matrix44 & operator /= (T a);
545
Matrix44 operator / (T a) const;
552
const Matrix44 & transpose ();
553
Matrix44 transposed () const;
556
//------------------------------------------------------------
557
// Inverse matrix: If singExc is false, inverting a singular
558
// matrix produces an identity matrix. If singExc is true,
559
// inverting a singular matrix throws a SingMatrixExc.
561
// inverse() and invert() invert matrices using determinants;
562
// gjInverse() and gjInvert() use the Gauss-Jordan method.
564
// inverse() and invert() are significantly faster than
565
// gjInverse() and gjInvert(), but the results may be slightly
568
//------------------------------------------------------------
570
const Matrix44 & invert (bool singExc = false)
571
throw (Iex::MathExc);
573
Matrix44<T> inverse (bool singExc = false) const
574
throw (Iex::MathExc);
576
const Matrix44 & gjInvert (bool singExc = false)
577
throw (Iex::MathExc);
579
Matrix44<T> gjInverse (bool singExc = false) const
580
throw (Iex::MathExc);
583
//--------------------------------------------------------
584
// Set matrix to rotation by XYZ euler angles (in radians)
585
//--------------------------------------------------------
588
const Matrix44 & setEulerAngles (const Vec3<S>& r);
591
//--------------------------------------------------------
592
// Set matrix to rotation around given axis by given angle
593
//--------------------------------------------------------
596
const Matrix44 & setAxisAngle (const Vec3<S>& ax, S ang);
599
//-------------------------------------------
600
// Rotate the matrix by XYZ euler angles in r
601
//-------------------------------------------
604
const Matrix44 & rotate (const Vec3<S> &r);
607
//--------------------------------------------
608
// Set matrix to scale by given uniform factor
609
//--------------------------------------------
611
const Matrix44 & setScale (T s);
614
//------------------------------------
615
// Set matrix to scale by given vector
616
//------------------------------------
619
const Matrix44 & setScale (const Vec3<S> &s);
622
//----------------------
623
// Scale the matrix by s
624
//----------------------
627
const Matrix44 & scale (const Vec3<S> &s);
630
//------------------------------------------
631
// Set matrix to translation by given vector
632
//------------------------------------------
635
const Matrix44 & setTranslation (const Vec3<S> &t);
638
//-----------------------------
639
// Return translation component
640
//-----------------------------
642
const Vec3<T> translation () const;
645
//--------------------------
646
// Translate the matrix by t
647
//--------------------------
650
const Matrix44 & translate (const Vec3<S> &t);
653
//-------------------------------------------------------------
654
// Set matrix to shear by given vector h. The resulting matrix
655
// will shear x for each y coord. by a factor of h[0] ;
656
// will shear x for each z coord. by a factor of h[1] ;
657
// will shear y for each z coord. by a factor of h[2] .
658
//-------------------------------------------------------------
661
const Matrix44 & setShear (const Vec3<S> &h);
664
//------------------------------------------------------------
665
// Set matrix to shear by given factors. The resulting matrix
666
// will shear x for each y coord. by a factor of h.xy ;
667
// will shear x for each z coord. by a factor of h.xz ;
668
// will shear y for each z coord. by a factor of h.yz ;
669
// will shear y for each x coord. by a factor of h.yx ;
670
// will shear z for each x coord. by a factor of h.zx ;
671
// will shear z for each y coord. by a factor of h.zy .
672
//------------------------------------------------------------
675
const Matrix44 & setShear (const Shear6<S> &h);
678
//--------------------------------------------------------
679
// Shear the matrix by given vector. The composed matrix
680
// will be <shear> * <this>, where the shear matrix ...
681
// will shear x for each y coord. by a factor of h[0] ;
682
// will shear x for each z coord. by a factor of h[1] ;
683
// will shear y for each z coord. by a factor of h[2] .
684
//--------------------------------------------------------
687
const Matrix44 & shear (const Vec3<S> &h);
690
//------------------------------------------------------------
691
// Shear the matrix by the given factors. The composed matrix
692
// will be <shear> * <this>, where the shear matrix ...
693
// will shear x for each y coord. by a factor of h.xy ;
694
// will shear x for each z coord. by a factor of h.xz ;
695
// will shear y for each z coord. by a factor of h.yz ;
696
// will shear y for each x coord. by a factor of h.yx ;
697
// will shear z for each x coord. by a factor of h.zx ;
698
// will shear z for each y coord. by a factor of h.zy .
699
//------------------------------------------------------------
702
const Matrix44 & shear (const Shear6<S> &h);
705
//-------------------------------------------------
706
// Limitations of type T (see also class limits<T>)
707
//-------------------------------------------------
709
static T baseTypeMin() {return limits<T>::min();}
710
static T baseTypeMax() {return limits<T>::max();}
711
static T baseTypeSmallest() {return limits<T>::smallest();}
712
static T baseTypeEpsilon() {return limits<T>::epsilon();}
721
std::ostream & operator << (std::ostream & s, const Matrix33<T> &m);
724
std::ostream & operator << (std::ostream & s, const Matrix44<T> &m);
727
//---------------------------------------------
728
// Vector-times-matrix multiplication operators
729
//---------------------------------------------
731
template <class S, class T>
732
const Vec2<S> & operator *= (Vec2<S> &v, const Matrix33<T> &m);
734
template <class S, class T>
735
Vec2<S> operator * (const Vec2<S> &v, const Matrix33<T> &m);
737
template <class S, class T>
738
const Vec3<S> & operator *= (Vec3<S> &v, const Matrix33<T> &m);
740
template <class S, class T>
741
Vec3<S> operator * (const Vec3<S> &v, const Matrix33<T> &m);
743
template <class S, class T>
744
const Vec3<S> & operator *= (Vec3<S> &v, const Matrix44<T> &m);
746
template <class S, class T>
747
Vec3<S> operator * (const Vec3<S> &v, const Matrix44<T> &m);
750
//-------------------------
751
// Typedefs for convenience
752
//-------------------------
754
typedef Matrix33 <float> M33f;
755
typedef Matrix33 <double> M33d;
756
typedef Matrix44 <float> M44f;
757
typedef Matrix44 <double> M44d;
760
//---------------------------
761
// Implementation of Matrix33
762
//---------------------------
766
Matrix33<T>::operator [] (int i)
773
Matrix33<T>::operator [] (int i) const
780
Matrix33<T>::Matrix33 ()
795
Matrix33<T>::Matrix33 (T a)
810
Matrix33<T>::Matrix33 (const T a[3][3])
825
Matrix33<T>::Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i)
840
Matrix33<T>::Matrix33 (const Matrix33 &v)
854
inline const Matrix33<T> &
855
Matrix33<T>::operator = (const Matrix33 &v)
870
inline const Matrix33<T> &
871
Matrix33<T>::operator = (T a)
887
Matrix33<T>::getValue ()
889
return (T *) &x[0][0];
894
Matrix33<T>::getValue () const
896
return (const T *) &x[0][0];
902
Matrix33<T>::getValue (Matrix33<S> &v) const
918
Matrix33<T>::setValue (const Matrix33<S> &v)
935
Matrix33<T>::setTheMatrix (const Matrix33<S> &v)
951
Matrix33<T>::makeIdentity()
966
Matrix33<T>::operator == (const Matrix33 &v) const
968
return x[0][0] == v.x[0][0] &&
969
x[0][1] == v.x[0][1] &&
970
x[0][2] == v.x[0][2] &&
971
x[1][0] == v.x[1][0] &&
972
x[1][1] == v.x[1][1] &&
973
x[1][2] == v.x[1][2] &&
974
x[2][0] == v.x[2][0] &&
975
x[2][1] == v.x[2][1] &&
976
x[2][2] == v.x[2][2];
981
Matrix33<T>::operator != (const Matrix33 &v) const
983
return x[0][0] != v.x[0][0] ||
984
x[0][1] != v.x[0][1] ||
985
x[0][2] != v.x[0][2] ||
986
x[1][0] != v.x[1][0] ||
987
x[1][1] != v.x[1][1] ||
988
x[1][2] != v.x[1][2] ||
989
x[2][0] != v.x[2][0] ||
990
x[2][1] != v.x[2][1] ||
991
x[2][2] != v.x[2][2];
996
Matrix33<T>::equalWithAbsError (const Matrix33<T> &m, T e) const
998
for (int i = 0; i < 3; i++)
999
for (int j = 0; j < 3; j++)
1000
if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e))
1008
Matrix33<T>::equalWithRelError (const Matrix33<T> &m, T e) const
1010
for (int i = 0; i < 3; i++)
1011
for (int j = 0; j < 3; j++)
1012
if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e))
1020
Matrix33<T>::operator += (const Matrix33<T> &v)
1022
x[0][0] += v.x[0][0];
1023
x[0][1] += v.x[0][1];
1024
x[0][2] += v.x[0][2];
1025
x[1][0] += v.x[1][0];
1026
x[1][1] += v.x[1][1];
1027
x[1][2] += v.x[1][2];
1028
x[2][0] += v.x[2][0];
1029
x[2][1] += v.x[2][1];
1030
x[2][2] += v.x[2][2];
1037
Matrix33<T>::operator += (T a)
1054
Matrix33<T>::operator + (const Matrix33<T> &v) const
1056
return Matrix33 (x[0][0] + v.x[0][0],
1057
x[0][1] + v.x[0][1],
1058
x[0][2] + v.x[0][2],
1059
x[1][0] + v.x[1][0],
1060
x[1][1] + v.x[1][1],
1061
x[1][2] + v.x[1][2],
1062
x[2][0] + v.x[2][0],
1063
x[2][1] + v.x[2][1],
1064
x[2][2] + v.x[2][2]);
1069
Matrix33<T>::operator -= (const Matrix33<T> &v)
1071
x[0][0] -= v.x[0][0];
1072
x[0][1] -= v.x[0][1];
1073
x[0][2] -= v.x[0][2];
1074
x[1][0] -= v.x[1][0];
1075
x[1][1] -= v.x[1][1];
1076
x[1][2] -= v.x[1][2];
1077
x[2][0] -= v.x[2][0];
1078
x[2][1] -= v.x[2][1];
1079
x[2][2] -= v.x[2][2];
1086
Matrix33<T>::operator -= (T a)
1103
Matrix33<T>::operator - (const Matrix33<T> &v) const
1105
return Matrix33 (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]);
1118
Matrix33<T>::operator - () const
1120
return Matrix33 (-x[0][0],
1133
Matrix33<T>::negate ()
1150
Matrix33<T>::operator *= (T a)
1167
Matrix33<T>::operator * (T a) const
1169
return Matrix33 (x[0][0] * a,
1182
operator * (T a, const Matrix33<T> &v)
1189
Matrix33<T>::operator *= (const Matrix33<T> &v)
1191
Matrix33 tmp (T (0));
1193
for (int i = 0; i < 3; i++)
1194
for (int j = 0; j < 3; j++)
1195
for (int k = 0; k < 3; k++)
1196
tmp.x[i][j] += x[i][k] * v.x[k][j];
1204
Matrix33<T>::operator * (const Matrix33<T> &v) const
1206
Matrix33 tmp (T (0));
1208
for (int i = 0; i < 3; i++)
1209
for (int j = 0; j < 3; j++)
1210
for (int k = 0; k < 3; k++)
1211
tmp.x[i][j] += x[i][k] * v.x[k][j];
1219
Matrix33<T>::multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const
1223
a = src[0] * x[0][0] + src[1] * x[1][0] + x[2][0];
1224
b = src[0] * x[0][1] + src[1] * x[1][1] + x[2][1];
1225
w = src[0] * x[0][2] + src[1] * x[1][2] + x[2][2];
1234
Matrix33<T>::multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const
1238
a = src[0] * x[0][0] + src[1] * x[1][0];
1239
b = src[0] * x[0][1] + src[1] * x[1][1];
1247
Matrix33<T>::operator /= (T a)
1264
Matrix33<T>::operator / (T a) const
1266
return Matrix33 (x[0][0] / a,
1279
Matrix33<T>::transpose ()
1281
Matrix33 tmp (x[0][0],
1296
Matrix33<T>::transposed () const
1298
return Matrix33 (x[0][0],
1311
Matrix33<T>::gjInvert (bool singExc) throw (Iex::MathExc)
1313
*this = gjInverse (singExc);
1319
Matrix33<T>::gjInverse (bool singExc) const throw (Iex::MathExc)
1325
// Forward elimination
1327
for (i = 0; i < 2 ; i++)
1331
T pivotsize = t[i][i];
1334
pivotsize = -pivotsize;
1336
for (j = i + 1; j < 3; j++)
1343
if (tmp > pivotsize)
1353
throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
1360
for (j = 0; j < 3; j++)
1365
t[i][j] = t[pivot][j];
1369
s[i][j] = s[pivot][j];
1374
for (j = i + 1; j < 3; j++)
1376
T f = t[j][i] / t[i][i];
1378
for (k = 0; k < 3; k++)
1380
t[j][k] -= f * t[i][k];
1381
s[j][k] -= f * s[i][k];
1386
// Backward substitution
1388
for (i = 2; i >= 0; --i)
1392
if ((f = t[i][i]) == 0)
1395
throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
1400
for (j = 0; j < 3; j++)
1406
for (j = 0; j < i; j++)
1410
for (k = 0; k < 3; k++)
1412
t[j][k] -= f * t[i][k];
1413
s[j][k] -= f * s[i][k];
1423
Matrix33<T>::invert (bool singExc) throw (Iex::MathExc)
1425
*this = inverse (singExc);
1431
Matrix33<T>::inverse (bool singExc) const throw (Iex::MathExc)
1433
if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1)
1435
Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
1436
x[2][1] * x[0][2] - x[0][1] * x[2][2],
1437
x[0][1] * x[1][2] - x[1][1] * x[0][2],
1439
x[2][0] * x[1][2] - x[1][0] * x[2][2],
1440
x[0][0] * x[2][2] - x[2][0] * x[0][2],
1441
x[1][0] * x[0][2] - x[0][0] * x[1][2],
1443
x[1][0] * x[2][1] - x[2][0] * x[1][1],
1444
x[2][0] * x[0][1] - x[0][0] * x[2][1],
1445
x[0][0] * x[1][1] - x[1][0] * x[0][1]);
1447
T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
1449
if (Imath::abs (r) >= 1)
1451
for (int i = 0; i < 3; ++i)
1453
for (int j = 0; j < 3; ++j)
1461
T mr = Imath::abs (r) / limits<T>::smallest();
1463
for (int i = 0; i < 3; ++i)
1465
for (int j = 0; j < 3; ++j)
1467
if (mr > Imath::abs (s[i][j]))
1474
throw SingMatrixExc ("Cannot invert "
1475
"singular matrix.");
1486
Matrix33 s ( x[1][1],
1498
T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
1500
if (Imath::abs (r) >= 1)
1502
for (int i = 0; i < 2; ++i)
1504
for (int j = 0; j < 2; ++j)
1512
T mr = Imath::abs (r) / limits<T>::smallest();
1514
for (int i = 0; i < 2; ++i)
1516
for (int j = 0; j < 2; ++j)
1518
if (mr > Imath::abs (s[i][j]))
1525
throw SingMatrixExc ("Cannot invert "
1526
"singular matrix.");
1533
s[2][0] = -x[2][0] * s[0][0] - x[2][1] * s[1][0];
1534
s[2][1] = -x[2][0] * s[0][1] - x[2][1] * s[1][1];
1543
Matrix33<T>::setRotation (S r)
1547
cos_r = Math<T>::cos (r);
1548
sin_r = Math<T>::sin (r);
1568
Matrix33<T>::rotate (S r)
1570
*this *= Matrix33<T>().setRotation (r);
1576
Matrix33<T>::setScale (T s)
1596
Matrix33<T>::setScale (const Vec2<S> &s)
1616
Matrix33<T>::scale (const Vec2<S> &s)
1632
Matrix33<T>::setTranslation (const Vec2<S> &t)
1651
Matrix33<T>::translation () const
1653
return Vec2<T> (x[2][0], x[2][1]);
1659
Matrix33<T>::translate (const Vec2<S> &t)
1661
x[2][0] += t[0] * x[0][0] + t[1] * x[1][0];
1662
x[2][1] += t[0] * x[0][1] + t[1] * x[1][1];
1663
x[2][2] += t[0] * x[0][2] + t[1] * x[1][2];
1671
Matrix33<T>::setShear (const S &xy)
1691
Matrix33<T>::setShear (const Vec2<S> &h)
1711
Matrix33<T>::shear (const S &xy)
1714
// In this case, we don't need a temp. copy of the matrix
1715
// because we never use a value on the RHS after we've
1716
// changed it on the LHS.
1719
x[1][0] += xy * x[0][0];
1720
x[1][1] += xy * x[0][1];
1721
x[1][2] += xy * x[0][2];
1729
Matrix33<T>::shear (const Vec2<S> &h)
1731
Matrix33<T> P (*this);
1733
x[0][0] = P[0][0] + h[1] * P[1][0];
1734
x[0][1] = P[0][1] + h[1] * P[1][1];
1735
x[0][2] = P[0][2] + h[1] * P[1][2];
1737
x[1][0] = P[1][0] + h[0] * P[0][0];
1738
x[1][1] = P[1][1] + h[0] * P[0][1];
1739
x[1][2] = P[1][2] + h[0] * P[0][2];
1745
//---------------------------
1746
// Implementation of Matrix44
1747
//---------------------------
1751
Matrix44<T>::operator [] (int i)
1758
Matrix44<T>::operator [] (int i) const
1765
Matrix44<T>::Matrix44 ()
1787
Matrix44<T>::Matrix44 (T a)
1809
Matrix44<T>::Matrix44 (const T a[4][4])
1831
Matrix44<T>::Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,
1832
T i, T j, T k, T l, T m, T n, T o, T p)
1855
Matrix44<T>::Matrix44 (Matrix33<T> r, Vec3<T> t)
1877
Matrix44<T>::Matrix44 (const Matrix44 &v)
1879
x[0][0] = v.x[0][0];
1880
x[0][1] = v.x[0][1];
1881
x[0][2] = v.x[0][2];
1882
x[0][3] = v.x[0][3];
1883
x[1][0] = v.x[1][0];
1884
x[1][1] = v.x[1][1];
1885
x[1][2] = v.x[1][2];
1886
x[1][3] = v.x[1][3];
1887
x[2][0] = v.x[2][0];
1888
x[2][1] = v.x[2][1];
1889
x[2][2] = v.x[2][2];
1890
x[2][3] = v.x[2][3];
1891
x[3][0] = v.x[3][0];
1892
x[3][1] = v.x[3][1];
1893
x[3][2] = v.x[3][2];
1894
x[3][3] = v.x[3][3];
1898
inline const Matrix44<T> &
1899
Matrix44<T>::operator = (const Matrix44 &v)
1901
x[0][0] = v.x[0][0];
1902
x[0][1] = v.x[0][1];
1903
x[0][2] = v.x[0][2];
1904
x[0][3] = v.x[0][3];
1905
x[1][0] = v.x[1][0];
1906
x[1][1] = v.x[1][1];
1907
x[1][2] = v.x[1][2];
1908
x[1][3] = v.x[1][3];
1909
x[2][0] = v.x[2][0];
1910
x[2][1] = v.x[2][1];
1911
x[2][2] = v.x[2][2];
1912
x[2][3] = v.x[2][3];
1913
x[3][0] = v.x[3][0];
1914
x[3][1] = v.x[3][1];
1915
x[3][2] = v.x[3][2];
1916
x[3][3] = v.x[3][3];
1921
inline const Matrix44<T> &
1922
Matrix44<T>::operator = (T a)
1945
Matrix44<T>::getValue ()
1947
return (T *) &x[0][0];
1952
Matrix44<T>::getValue () const
1954
return (const T *) &x[0][0];
1960
Matrix44<T>::getValue (Matrix44<S> &v) const
1962
v.x[0][0] = x[0][0];
1963
v.x[0][1] = x[0][1];
1964
v.x[0][2] = x[0][2];
1965
v.x[0][3] = x[0][3];
1966
v.x[1][0] = x[1][0];
1967
v.x[1][1] = x[1][1];
1968
v.x[1][2] = x[1][2];
1969
v.x[1][3] = x[1][3];
1970
v.x[2][0] = x[2][0];
1971
v.x[2][1] = x[2][1];
1972
v.x[2][2] = x[2][2];
1973
v.x[2][3] = x[2][3];
1974
v.x[3][0] = x[3][0];
1975
v.x[3][1] = x[3][1];
1976
v.x[3][2] = x[3][2];
1977
v.x[3][3] = x[3][3];
1982
inline Matrix44<T> &
1983
Matrix44<T>::setValue (const Matrix44<S> &v)
1985
x[0][0] = v.x[0][0];
1986
x[0][1] = v.x[0][1];
1987
x[0][2] = v.x[0][2];
1988
x[0][3] = v.x[0][3];
1989
x[1][0] = v.x[1][0];
1990
x[1][1] = v.x[1][1];
1991
x[1][2] = v.x[1][2];
1992
x[1][3] = v.x[1][3];
1993
x[2][0] = v.x[2][0];
1994
x[2][1] = v.x[2][1];
1995
x[2][2] = v.x[2][2];
1996
x[2][3] = v.x[2][3];
1997
x[3][0] = v.x[3][0];
1998
x[3][1] = v.x[3][1];
1999
x[3][2] = v.x[3][2];
2000
x[3][3] = v.x[3][3];
2006
inline Matrix44<T> &
2007
Matrix44<T>::setTheMatrix (const Matrix44<S> &v)
2009
x[0][0] = v.x[0][0];
2010
x[0][1] = v.x[0][1];
2011
x[0][2] = v.x[0][2];
2012
x[0][3] = v.x[0][3];
2013
x[1][0] = v.x[1][0];
2014
x[1][1] = v.x[1][1];
2015
x[1][2] = v.x[1][2];
2016
x[1][3] = v.x[1][3];
2017
x[2][0] = v.x[2][0];
2018
x[2][1] = v.x[2][1];
2019
x[2][2] = v.x[2][2];
2020
x[2][3] = v.x[2][3];
2021
x[3][0] = v.x[3][0];
2022
x[3][1] = v.x[3][1];
2023
x[3][2] = v.x[3][2];
2024
x[3][3] = v.x[3][3];
2030
Matrix44<T>::makeIdentity()
2052
Matrix44<T>::operator == (const Matrix44 &v) const
2054
return x[0][0] == v.x[0][0] &&
2055
x[0][1] == v.x[0][1] &&
2056
x[0][2] == v.x[0][2] &&
2057
x[0][3] == v.x[0][3] &&
2058
x[1][0] == v.x[1][0] &&
2059
x[1][1] == v.x[1][1] &&
2060
x[1][2] == v.x[1][2] &&
2061
x[1][3] == v.x[1][3] &&
2062
x[2][0] == v.x[2][0] &&
2063
x[2][1] == v.x[2][1] &&
2064
x[2][2] == v.x[2][2] &&
2065
x[2][3] == v.x[2][3] &&
2066
x[3][0] == v.x[3][0] &&
2067
x[3][1] == v.x[3][1] &&
2068
x[3][2] == v.x[3][2] &&
2069
x[3][3] == v.x[3][3];
2074
Matrix44<T>::operator != (const Matrix44 &v) const
2076
return x[0][0] != v.x[0][0] ||
2077
x[0][1] != v.x[0][1] ||
2078
x[0][2] != v.x[0][2] ||
2079
x[0][3] != v.x[0][3] ||
2080
x[1][0] != v.x[1][0] ||
2081
x[1][1] != v.x[1][1] ||
2082
x[1][2] != v.x[1][2] ||
2083
x[1][3] != v.x[1][3] ||
2084
x[2][0] != v.x[2][0] ||
2085
x[2][1] != v.x[2][1] ||
2086
x[2][2] != v.x[2][2] ||
2087
x[2][3] != v.x[2][3] ||
2088
x[3][0] != v.x[3][0] ||
2089
x[3][1] != v.x[3][1] ||
2090
x[3][2] != v.x[3][2] ||
2091
x[3][3] != v.x[3][3];
2096
Matrix44<T>::equalWithAbsError (const Matrix44<T> &m, T e) const
2098
for (int i = 0; i < 4; i++)
2099
for (int j = 0; j < 4; j++)
2100
if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e))
2108
Matrix44<T>::equalWithRelError (const Matrix44<T> &m, T e) const
2110
for (int i = 0; i < 4; i++)
2111
for (int j = 0; j < 4; j++)
2112
if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e))
2120
Matrix44<T>::operator += (const Matrix44<T> &v)
2122
x[0][0] += v.x[0][0];
2123
x[0][1] += v.x[0][1];
2124
x[0][2] += v.x[0][2];
2125
x[0][3] += v.x[0][3];
2126
x[1][0] += v.x[1][0];
2127
x[1][1] += v.x[1][1];
2128
x[1][2] += v.x[1][2];
2129
x[1][3] += v.x[1][3];
2130
x[2][0] += v.x[2][0];
2131
x[2][1] += v.x[2][1];
2132
x[2][2] += v.x[2][2];
2133
x[2][3] += v.x[2][3];
2134
x[3][0] += v.x[3][0];
2135
x[3][1] += v.x[3][1];
2136
x[3][2] += v.x[3][2];
2137
x[3][3] += v.x[3][3];
2144
Matrix44<T>::operator += (T a)
2168
Matrix44<T>::operator + (const Matrix44<T> &v) const
2170
return Matrix44 (x[0][0] + v.x[0][0],
2171
x[0][1] + v.x[0][1],
2172
x[0][2] + v.x[0][2],
2173
x[0][3] + v.x[0][3],
2174
x[1][0] + v.x[1][0],
2175
x[1][1] + v.x[1][1],
2176
x[1][2] + v.x[1][2],
2177
x[1][3] + v.x[1][3],
2178
x[2][0] + v.x[2][0],
2179
x[2][1] + v.x[2][1],
2180
x[2][2] + v.x[2][2],
2181
x[2][3] + v.x[2][3],
2182
x[3][0] + v.x[3][0],
2183
x[3][1] + v.x[3][1],
2184
x[3][2] + v.x[3][2],
2185
x[3][3] + v.x[3][3]);
2190
Matrix44<T>::operator -= (const Matrix44<T> &v)
2192
x[0][0] -= v.x[0][0];
2193
x[0][1] -= v.x[0][1];
2194
x[0][2] -= v.x[0][2];
2195
x[0][3] -= v.x[0][3];
2196
x[1][0] -= v.x[1][0];
2197
x[1][1] -= v.x[1][1];
2198
x[1][2] -= v.x[1][2];
2199
x[1][3] -= v.x[1][3];
2200
x[2][0] -= v.x[2][0];
2201
x[2][1] -= v.x[2][1];
2202
x[2][2] -= v.x[2][2];
2203
x[2][3] -= v.x[2][3];
2204
x[3][0] -= v.x[3][0];
2205
x[3][1] -= v.x[3][1];
2206
x[3][2] -= v.x[3][2];
2207
x[3][3] -= v.x[3][3];
2214
Matrix44<T>::operator -= (T a)
2238
Matrix44<T>::operator - (const Matrix44<T> &v) const
2240
return Matrix44 (x[0][0] - v.x[0][0],
2241
x[0][1] - v.x[0][1],
2242
x[0][2] - v.x[0][2],
2243
x[0][3] - v.x[0][3],
2244
x[1][0] - v.x[1][0],
2245
x[1][1] - v.x[1][1],
2246
x[1][2] - v.x[1][2],
2247
x[1][3] - v.x[1][3],
2248
x[2][0] - v.x[2][0],
2249
x[2][1] - v.x[2][1],
2250
x[2][2] - v.x[2][2],
2251
x[2][3] - v.x[2][3],
2252
x[3][0] - v.x[3][0],
2253
x[3][1] - v.x[3][1],
2254
x[3][2] - v.x[3][2],
2255
x[3][3] - v.x[3][3]);
2260
Matrix44<T>::operator - () const
2262
return Matrix44 (-x[0][0],
2282
Matrix44<T>::negate ()
2306
Matrix44<T>::operator *= (T a)
2330
Matrix44<T>::operator * (T a) const
2332
return Matrix44 (x[0][0] * a,
2352
operator * (T a, const Matrix44<T> &v)
2358
inline const Matrix44<T> &
2359
Matrix44<T>::operator *= (const Matrix44<T> &v)
2361
Matrix44 tmp (T (0));
2363
multiply (*this, v, tmp);
2370
Matrix44<T>::operator * (const Matrix44<T> &v) const
2372
Matrix44 tmp (T (0));
2374
multiply (*this, v, tmp);
2380
Matrix44<T>::multiply (const Matrix44<T> &a,
2381
const Matrix44<T> &b,
2384
register const T * restrict ap = &a.x[0][0];
2385
register const T * restrict bp = &b.x[0][0];
2386
register T * restrict cp = &c.x[0][0];
2388
register T a0, a1, a2, a3;
2395
cp[0] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2396
cp[1] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2397
cp[2] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2398
cp[3] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2405
cp[4] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2406
cp[5] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2407
cp[6] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2408
cp[7] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2415
cp[8] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2416
cp[9] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2417
cp[10] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2418
cp[11] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2425
cp[12] = a0 * bp[0] + a1 * bp[4] + a2 * bp[8] + a3 * bp[12];
2426
cp[13] = a0 * bp[1] + a1 * bp[5] + a2 * bp[9] + a3 * bp[13];
2427
cp[14] = a0 * bp[2] + a1 * bp[6] + a2 * bp[10] + a3 * bp[14];
2428
cp[15] = a0 * bp[3] + a1 * bp[7] + a2 * bp[11] + a3 * bp[15];
2431
template <class T> template <class S>
2433
Matrix44<T>::multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const
2437
a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0] + x[3][0];
2438
b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1] + x[3][1];
2439
c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2] + x[3][2];
2440
w = src[0] * x[0][3] + src[1] * x[1][3] + src[2] * x[2][3] + x[3][3];
2447
template <class T> template <class S>
2449
Matrix44<T>::multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const
2453
a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0];
2454
b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1];
2455
c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2];
2464
Matrix44<T>::operator /= (T a)
2488
Matrix44<T>::operator / (T a) const
2490
return Matrix44 (x[0][0] / a,
2510
Matrix44<T>::transpose ()
2512
Matrix44 tmp (x[0][0],
2534
Matrix44<T>::transposed () const
2536
return Matrix44 (x[0][0],
2556
Matrix44<T>::gjInvert (bool singExc) throw (Iex::MathExc)
2558
*this = gjInverse (singExc);
2564
Matrix44<T>::gjInverse (bool singExc) const throw (Iex::MathExc)
2570
// Forward elimination
2572
for (i = 0; i < 3 ; i++)
2576
T pivotsize = t[i][i];
2579
pivotsize = -pivotsize;
2581
for (j = i + 1; j < 4; j++)
2588
if (tmp > pivotsize)
2598
throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
2605
for (j = 0; j < 4; j++)
2610
t[i][j] = t[pivot][j];
2614
s[i][j] = s[pivot][j];
2619
for (j = i + 1; j < 4; j++)
2621
T f = t[j][i] / t[i][i];
2623
for (k = 0; k < 4; k++)
2625
t[j][k] -= f * t[i][k];
2626
s[j][k] -= f * s[i][k];
2631
// Backward substitution
2633
for (i = 3; i >= 0; --i)
2637
if ((f = t[i][i]) == 0)
2640
throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
2645
for (j = 0; j < 4; j++)
2651
for (j = 0; j < i; j++)
2655
for (k = 0; k < 4; k++)
2657
t[j][k] -= f * t[i][k];
2658
s[j][k] -= f * s[i][k];
2668
Matrix44<T>::invert (bool singExc) throw (Iex::MathExc)
2670
*this = inverse (singExc);
2676
Matrix44<T>::inverse (bool singExc) const throw (Iex::MathExc)
2678
if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)
2679
return gjInverse(singExc);
2681
Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
2682
x[2][1] * x[0][2] - x[0][1] * x[2][2],
2683
x[0][1] * x[1][2] - x[1][1] * x[0][2],
2686
x[2][0] * x[1][2] - x[1][0] * x[2][2],
2687
x[0][0] * x[2][2] - x[2][0] * x[0][2],
2688
x[1][0] * x[0][2] - x[0][0] * x[1][2],
2691
x[1][0] * x[2][1] - x[2][0] * x[1][1],
2692
x[2][0] * x[0][1] - x[0][0] * x[2][1],
2693
x[0][0] * x[1][1] - x[1][0] * x[0][1],
2701
T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
2703
if (Imath::abs (r) >= 1)
2705
for (int i = 0; i < 3; ++i)
2707
for (int j = 0; j < 3; ++j)
2715
T mr = Imath::abs (r) / limits<T>::smallest();
2717
for (int i = 0; i < 3; ++i)
2719
for (int j = 0; j < 3; ++j)
2721
if (mr > Imath::abs (s[i][j]))
2728
throw SingMatrixExc ("Cannot invert singular matrix.");
2736
s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0];
2737
s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1];
2738
s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2];
2746
Matrix44<T>::setEulerAngles (const Vec3<S>& r)
2748
S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
2750
cos_rz = Math<T>::cos (r[2]);
2751
cos_ry = Math<T>::cos (r[1]);
2752
cos_rx = Math<T>::cos (r[0]);
2754
sin_rz = Math<T>::sin (r[2]);
2755
sin_ry = Math<T>::sin (r[1]);
2756
sin_rx = Math<T>::sin (r[0]);
2758
x[0][0] = cos_rz * cos_ry;
2759
x[0][1] = sin_rz * cos_ry;
2763
x[1][0] = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
2764
x[1][1] = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
2765
x[1][2] = cos_ry * sin_rx;
2768
x[2][0] = sin_rz * sin_rx + cos_rz * sin_ry * cos_rx;
2769
x[2][1] = -cos_rz * sin_rx + sin_rz * sin_ry * cos_rx;
2770
x[2][2] = cos_ry * cos_rx;
2784
Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle)
2786
Vec3<S> unit (axis.normalized());
2787
S sine = Math<T>::sin (angle);
2788
S cosine = Math<T>::cos (angle);
2790
x[0][0] = unit[0] * unit[0] * (1 - cosine) + cosine;
2791
x[0][1] = unit[0] * unit[1] * (1 - cosine) + unit[2] * sine;
2792
x[0][2] = unit[0] * unit[2] * (1 - cosine) - unit[1] * sine;
2795
x[1][0] = unit[0] * unit[1] * (1 - cosine) - unit[2] * sine;
2796
x[1][1] = unit[1] * unit[1] * (1 - cosine) + cosine;
2797
x[1][2] = unit[1] * unit[2] * (1 - cosine) + unit[0] * sine;
2800
x[2][0] = unit[0] * unit[2] * (1 - cosine) + unit[1] * sine;
2801
x[2][1] = unit[1] * unit[2] * (1 - cosine) - unit[0] * sine;
2802
x[2][2] = unit[2] * unit[2] * (1 - cosine) + cosine;
2816
Matrix44<T>::rotate (const Vec3<S> &r)
2818
S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
2823
cos_rz = Math<S>::cos (r[2]);
2824
cos_ry = Math<S>::cos (r[1]);
2825
cos_rx = Math<S>::cos (r[0]);
2827
sin_rz = Math<S>::sin (r[2]);
2828
sin_ry = Math<S>::sin (r[1]);
2829
sin_rx = Math<S>::sin (r[0]);
2831
m00 = cos_rz * cos_ry;
2832
m01 = sin_rz * cos_ry;
2834
m10 = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
2835
m11 = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
2836
m12 = cos_ry * sin_rx;
2837
m20 = -sin_rz * -sin_rx + cos_rz * sin_ry * cos_rx;
2838
m21 = cos_rz * -sin_rx + sin_rz * sin_ry * cos_rx;
2839
m22 = cos_ry * cos_rx;
2841
Matrix44<T> P (*this);
2843
x[0][0] = P[0][0] * m00 + P[1][0] * m01 + P[2][0] * m02;
2844
x[0][1] = P[0][1] * m00 + P[1][1] * m01 + P[2][1] * m02;
2845
x[0][2] = P[0][2] * m00 + P[1][2] * m01 + P[2][2] * m02;
2846
x[0][3] = P[0][3] * m00 + P[1][3] * m01 + P[2][3] * m02;
2848
x[1][0] = P[0][0] * m10 + P[1][0] * m11 + P[2][0] * m12;
2849
x[1][1] = P[0][1] * m10 + P[1][1] * m11 + P[2][1] * m12;
2850
x[1][2] = P[0][2] * m10 + P[1][2] * m11 + P[2][2] * m12;
2851
x[1][3] = P[0][3] * m10 + P[1][3] * m11 + P[2][3] * m12;
2853
x[2][0] = P[0][0] * m20 + P[1][0] * m21 + P[2][0] * m22;
2854
x[2][1] = P[0][1] * m20 + P[1][1] * m21 + P[2][1] * m22;
2855
x[2][2] = P[0][2] * m20 + P[1][2] * m21 + P[2][2] * m22;
2856
x[2][3] = P[0][3] * m20 + P[1][3] * m21 + P[2][3] * m22;
2863
Matrix44<T>::setScale (T s)
2891
Matrix44<T>::setScale (const Vec3<S> &s)
2919
Matrix44<T>::scale (const Vec3<S> &s)
2942
Matrix44<T>::setTranslation (const Vec3<S> &t)
2968
inline const Vec3<T>
2969
Matrix44<T>::translation () const
2971
return Vec3<T> (x[3][0], x[3][1], x[3][2]);
2977
Matrix44<T>::translate (const Vec3<S> &t)
2979
x[3][0] += t[0] * x[0][0] + t[1] * x[1][0] + t[2] * x[2][0];
2980
x[3][1] += t[0] * x[0][1] + t[1] * x[1][1] + t[2] * x[2][1];
2981
x[3][2] += t[0] * x[0][2] + t[1] * x[1][2] + t[2] * x[2][2];
2982
x[3][3] += t[0] * x[0][3] + t[1] * x[1][3] + t[2] * x[2][3];
2990
Matrix44<T>::setShear (const Vec3<S> &h)
3018
Matrix44<T>::setShear (const Shear6<S> &h)
3046
Matrix44<T>::shear (const Vec3<S> &h)
3049
// In this case, we don't need a temp. copy of the matrix
3050
// because we never use a value on the RHS after we've
3051
// changed it on the LHS.
3054
for (int i=0; i < 4; i++)
3056
x[2][i] += h[1] * x[0][i] + h[2] * x[1][i];
3057
x[1][i] += h[0] * x[0][i];
3066
Matrix44<T>::shear (const Shear6<S> &h)
3068
Matrix44<T> P (*this);
3070
for (int i=0; i < 4; i++)
3072
x[0][i] = P[0][i] + h.yx * P[1][i] + h.zx * P[2][i];
3073
x[1][i] = h.xy * P[0][i] + P[1][i] + h.zy * P[2][i];
3074
x[2][i] = h.xz * P[0][i] + h.yz * P[1][i] + P[2][i];
3081
//--------------------------------
3082
// Implementation of stream output
3083
//--------------------------------
3087
operator << (std::ostream &s, const Matrix33<T> &m)
3089
std::ios_base::fmtflags oldFlags = s.flags();
3092
if (s.flags() & std::ios_base::fixed)
3094
s.setf (std::ios_base::showpoint);
3095
width = s.precision() + 5;
3099
s.setf (std::ios_base::scientific);
3100
s.setf (std::ios_base::showpoint);
3101
width = s.precision() + 8;
3104
s << "(" << std::setw (width) << m[0][0] <<
3105
" " << std::setw (width) << m[0][1] <<
3106
" " << std::setw (width) << m[0][2] << "\n" <<
3108
" " << std::setw (width) << m[1][0] <<
3109
" " << std::setw (width) << m[1][1] <<
3110
" " << std::setw (width) << m[1][2] << "\n" <<
3112
" " << std::setw (width) << m[2][0] <<
3113
" " << std::setw (width) << m[2][1] <<
3114
" " << std::setw (width) << m[2][2] << ")\n";
3122
operator << (std::ostream &s, const Matrix44<T> &m)
3124
std::ios_base::fmtflags oldFlags = s.flags();
3127
if (s.flags() & std::ios_base::fixed)
3129
s.setf (std::ios_base::showpoint);
3130
width = s.precision() + 5;
3134
s.setf (std::ios_base::scientific);
3135
s.setf (std::ios_base::showpoint);
3136
width = s.precision() + 8;
3139
s << "(" << std::setw (width) << m[0][0] <<
3140
" " << std::setw (width) << m[0][1] <<
3141
" " << std::setw (width) << m[0][2] <<
3142
" " << std::setw (width) << m[0][3] << "\n" <<
3144
" " << std::setw (width) << m[1][0] <<
3145
" " << std::setw (width) << m[1][1] <<
3146
" " << std::setw (width) << m[1][2] <<
3147
" " << std::setw (width) << m[1][3] << "\n" <<
3149
" " << std::setw (width) << m[2][0] <<
3150
" " << std::setw (width) << m[2][1] <<
3151
" " << std::setw (width) << m[2][2] <<
3152
" " << std::setw (width) << m[2][3] << "\n" <<
3154
" " << std::setw (width) << m[3][0] <<
3155
" " << std::setw (width) << m[3][1] <<
3156
" " << std::setw (width) << m[3][2] <<
3157
" " << std::setw (width) << m[3][3] << ")\n";
3164
//---------------------------------------------------------------
3165
// Implementation of vector-times-matrix multiplication operators
3166
//---------------------------------------------------------------
3168
template <class S, class T>
3169
inline const Vec2<S> &
3170
operator *= (Vec2<S> &v, const Matrix33<T> &m)
3172
S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
3173
S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
3174
S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
3182
template <class S, class T>
3184
operator * (const Vec2<S> &v, const Matrix33<T> &m)
3186
S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
3187
S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
3188
S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
3190
return Vec2<S> (x / w, y / w);
3194
template <class S, class T>
3195
inline const Vec3<S> &
3196
operator *= (Vec3<S> &v, const Matrix33<T> &m)
3198
S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
3199
S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
3200
S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
3210
template <class S, class T>
3212
operator * (const Vec3<S> &v, const Matrix33<T> &m)
3214
S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
3215
S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
3216
S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
3218
return Vec3<S> (x, y, z);
3222
template <class S, class T>
3223
inline const Vec3<S> &
3224
operator *= (Vec3<S> &v, const Matrix44<T> &m)
3226
S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
3227
S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
3228
S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
3229
S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
3238
template <class S, class T>
3240
operator * (const Vec3<S> &v, const Matrix44<T> &m)
3242
S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
3243
S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
3244
S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
3245
S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
3247
return Vec3<S> (x / w, y / w, z / w);
3250
} // namespace Imath
added
removed
Source/OpenEXR/Imath/ImathMatrix.h
