Coverage Report

Created: 2023-12-08 06:53

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