PageRenderTime 147ms CodeModel.GetById 33ms app.highlight 101ms RepoModel.GetById 1ms 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

Large files files are truncated, but you can click here to view the full 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  

Large files files are truncated, but you can click here to view the full file