PageRenderTime 249ms CodeModel.GetById 29ms app.highlight 201ms RepoModel.GetById 2ms app.codeStats 1ms

/src/FreeImage/Source/OpenEXR/Imath/ImathMatrix.h

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