← Back to branch summary

~paparazzi-uav/paparazzi/v5.0-manual

« back to all changes in this revision

Viewing changes to sw/ext/opencv_bebop/opencv/3rdparty/openexr/Imath/ImathMatrix.h

  • Committer: Paparazzi buildbot
  • Date: 2016-05-18 15:00:29 UTC
  • Revision ID: felix.ruess+docbot@gmail.com-20160518150029-e8lgzi5kvb4p7un9
Manual import commit 4b8bbb730080dac23cf816b98908dacfabe2a8ec from v5.0 branch.

Show diffs side-by-side

added added

removed removed

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