PBRT
/home/felix/UBC/projects/AdaptiveLightfieldSampling/pbrt_v2/src/3rdparty/ilmbase-1.0.2/ImathMatrix.h
00001 
00002 //
00003 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
00004 // Digital Ltd. LLC
00005 // 
00006 // All rights reserved.
00007 // 
00008 // Redistribution and use in source and binary forms, with or without
00009 // modification, are permitted provided that the following conditions are
00010 // met:
00011 // *       Redistributions of source code must retain the above copyright
00012 // notice, this list of conditions and the following disclaimer.
00013 // *       Redistributions in binary form must reproduce the above
00014 // copyright notice, this list of conditions and the following disclaimer
00015 // in the documentation and/or other materials provided with the
00016 // distribution.
00017 // *       Neither the name of Industrial Light & Magic nor the names of
00018 // its contributors may be used to endorse or promote products derived
00019 // from this software without specific prior written permission. 
00020 // 
00021 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00022 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00023 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00024 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
00025 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
00027 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
00028 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
00029 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
00030 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
00031 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00032 //
00034 
00035 
00036 
00037 #ifndef INCLUDED_IMATHMATRIX_H
00038 #define INCLUDED_IMATHMATRIX_H
00039 
00040 //----------------------------------------------------------------
00041 //
00042 //      2D (3x3) and 3D (4x4) transformation matrix templates.
00043 //
00044 //----------------------------------------------------------------
00045 
00046 #include "ImathPlatform.h"
00047 #include "ImathFun.h"
00048 #include "ImathExc.h"
00049 #include "ImathVec.h"
00050 #include "ImathShear.h"
00051 
00052 #include <iostream>
00053 #include <iomanip>
00054 
00055 #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
00056 // suppress exception specification warnings
00057 #pragma warning(disable:4290)
00058 #endif
00059 
00060 
00061 namespace Imath {
00062 
00063 enum Uninitialized {UNINITIALIZED};
00064 
00065 
00066 template <class T> class Matrix33
00067 {
00068   public:
00069 
00070     //-------------------
00071     // Access to elements
00072     //-------------------
00073 
00074     T           x[3][3];
00075 
00076     T *         operator [] (int i);
00077     const T *   operator [] (int i) const;
00078 
00079 
00080     //-------------
00081     // Constructors
00082     //-------------
00083 
00084     Matrix33 (Uninitialized) {}
00085 
00086     Matrix33 ();
00087                                 // 1 0 0
00088                                 // 0 1 0
00089                                 // 0 0 1
00090 
00091     Matrix33 (T a);
00092                                 // a a a
00093                                 // a a a
00094                                 // a a a
00095 
00096     Matrix33 (const T a[3][3]);
00097                                 // a[0][0] a[0][1] a[0][2]
00098                                 // a[1][0] a[1][1] a[1][2]
00099                                 // a[2][0] a[2][1] a[2][2]
00100 
00101     Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i);
00102 
00103                                 // a b c
00104                                 // d e f
00105                                 // g h i
00106 
00107 
00108     //--------------------------------
00109     // Copy constructor and assignment
00110     //--------------------------------
00111 
00112     Matrix33 (const Matrix33 &v);
00113     template <class S> explicit Matrix33 (const Matrix33<S> &v);
00114 
00115     const Matrix33 &    operator = (const Matrix33 &v);
00116     const Matrix33 &    operator = (T a);
00117 
00118 
00119     //----------------------
00120     // Compatibility with Sb
00121     //----------------------
00122     
00123     T *                 getValue ();
00124     const T *           getValue () const;
00125 
00126     template <class S>
00127     void                getValue (Matrix33<S> &v) const;
00128     template <class S>
00129     Matrix33 &          setValue (const Matrix33<S> &v);
00130 
00131     template <class S>
00132     Matrix33 &          setTheMatrix (const Matrix33<S> &v);
00133 
00134 
00135     //---------
00136     // Identity
00137     //---------
00138 
00139     void                makeIdentity();
00140 
00141 
00142     //---------
00143     // Equality
00144     //---------
00145 
00146     bool                operator == (const Matrix33 &v) const;
00147     bool                operator != (const Matrix33 &v) const;
00148 
00149     //-----------------------------------------------------------------------
00150     // Compare two matrices and test if they are "approximately equal":
00151     //
00152     // equalWithAbsError (m, e)
00153     //
00154     //      Returns true if the coefficients of this and m are the same with
00155     //      an absolute error of no more than e, i.e., for all i, j
00156     //
00157     //      abs (this[i][j] - m[i][j]) <= e
00158     //
00159     // equalWithRelError (m, e)
00160     //
00161     //      Returns true if the coefficients of this and m are the same with
00162     //      a relative error of no more than e, i.e., for all i, j
00163     //
00164     //      abs (this[i] - v[i][j]) <= e * abs (this[i][j])
00165     //-----------------------------------------------------------------------
00166 
00167     bool                equalWithAbsError (const Matrix33<T> &v, T e) const;
00168     bool                equalWithRelError (const Matrix33<T> &v, T e) const;
00169 
00170 
00171     //------------------------
00172     // Component-wise addition
00173     //------------------------
00174 
00175     const Matrix33 &    operator += (const Matrix33 &v);
00176     const Matrix33 &    operator += (T a);
00177     Matrix33            operator + (const Matrix33 &v) const;
00178 
00179 
00180     //---------------------------
00181     // Component-wise subtraction
00182     //---------------------------
00183 
00184     const Matrix33 &    operator -= (const Matrix33 &v);
00185     const Matrix33 &    operator -= (T a);
00186     Matrix33            operator - (const Matrix33 &v) const;
00187 
00188 
00189     //------------------------------------
00190     // Component-wise multiplication by -1
00191     //------------------------------------
00192 
00193     Matrix33            operator - () const;
00194     const Matrix33 &    negate ();
00195 
00196 
00197     //------------------------------
00198     // Component-wise multiplication
00199     //------------------------------
00200 
00201     const Matrix33 &    operator *= (T a);
00202     Matrix33            operator * (T a) const;
00203 
00204 
00205     //-----------------------------------
00206     // Matrix-times-matrix multiplication
00207     //-----------------------------------
00208 
00209     const Matrix33 &    operator *= (const Matrix33 &v);
00210     Matrix33            operator * (const Matrix33 &v) const;
00211 
00212 
00213     //-----------------------------------------------------------------
00214     // Vector-times-matrix multiplication; see also the "operator *"
00215     // functions defined below.
00216     //
00217     // m.multVecMatrix(src,dst) implements a homogeneous transformation
00218     // by computing Vec3 (src.x, src.y, 1) * m and dividing by the
00219     // result's third element.
00220     //
00221     // m.multDirMatrix(src,dst) multiplies src by the upper left 2x2
00222     // submatrix, ignoring the rest of matrix m.
00223     //-----------------------------------------------------------------
00224 
00225     template <class S>
00226     void                multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const;
00227 
00228     template <class S>
00229     void                multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const;
00230 
00231 
00232     //------------------------
00233     // Component-wise division
00234     //------------------------
00235 
00236     const Matrix33 &    operator /= (T a);
00237     Matrix33            operator / (T a) const;
00238 
00239 
00240     //------------------
00241     // Transposed matrix
00242     //------------------
00243 
00244     const Matrix33 &    transpose ();
00245     Matrix33            transposed () const;
00246 
00247 
00248     //------------------------------------------------------------
00249     // Inverse matrix: If singExc is false, inverting a singular
00250     // matrix produces an identity matrix.  If singExc is true,
00251     // inverting a singular matrix throws a SingMatrixExc.
00252     //
00253     // inverse() and invert() invert matrices using determinants;
00254     // gjInverse() and gjInvert() use the Gauss-Jordan method.
00255     //
00256     // inverse() and invert() are significantly faster than
00257     // gjInverse() and gjInvert(), but the results may be slightly
00258     // less accurate.
00259     // 
00260     //------------------------------------------------------------
00261 
00262     const Matrix33 &    invert (bool singExc = false)
00263                         throw (Iex::MathExc);
00264 
00265     Matrix33<T>         inverse (bool singExc = false) const
00266                         throw (Iex::MathExc);
00267 
00268     const Matrix33 &    gjInvert (bool singExc = false)
00269                         throw (Iex::MathExc);
00270 
00271     Matrix33<T>         gjInverse (bool singExc = false) const
00272                         throw (Iex::MathExc);
00273 
00274 
00275     //-----------------------------------------
00276     // Set matrix to rotation by r (in radians)
00277     //-----------------------------------------
00278 
00279     template <class S>
00280     const Matrix33 &    setRotation (S r);
00281 
00282 
00283     //-----------------------------
00284     // Rotate the given matrix by r
00285     //-----------------------------
00286 
00287     template <class S>
00288     const Matrix33 &    rotate (S r);
00289 
00290 
00291     //--------------------------------------------
00292     // Set matrix to scale by given uniform factor
00293     //--------------------------------------------
00294 
00295     const Matrix33 &    setScale (T s);
00296 
00297 
00298     //------------------------------------
00299     // Set matrix to scale by given vector
00300     //------------------------------------
00301 
00302     template <class S>
00303     const Matrix33 &    setScale (const Vec2<S> &s);
00304 
00305 
00306     //----------------------
00307     // Scale the matrix by s
00308     //----------------------
00309 
00310     template <class S>
00311     const Matrix33 &    scale (const Vec2<S> &s);
00312 
00313 
00314     //------------------------------------------
00315     // Set matrix to translation by given vector
00316     //------------------------------------------
00317 
00318     template <class S>
00319     const Matrix33 &    setTranslation (const Vec2<S> &t);
00320 
00321 
00322     //-----------------------------
00323     // Return translation component
00324     //-----------------------------
00325 
00326     Vec2<T>             translation () const;
00327 
00328 
00329     //--------------------------
00330     // Translate the matrix by t
00331     //--------------------------
00332 
00333     template <class S>
00334     const Matrix33 &    translate (const Vec2<S> &t);
00335 
00336 
00337     //-----------------------------------------------------------
00338     // Set matrix to shear x for each y coord. by given factor xy
00339     //-----------------------------------------------------------
00340 
00341     template <class S>
00342     const Matrix33 &    setShear (const S &h);
00343 
00344 
00345     //-------------------------------------------------------------
00346     // Set matrix to shear x for each y coord. by given factor h[0]
00347     // and to shear y for each x coord. by given factor h[1]
00348     //-------------------------------------------------------------
00349 
00350     template <class S>
00351     const Matrix33 &    setShear (const Vec2<S> &h);
00352 
00353 
00354     //-----------------------------------------------------------
00355     // Shear the matrix in x for each y coord. by given factor xy
00356     //-----------------------------------------------------------
00357 
00358     template <class S>
00359     const Matrix33 &    shear (const S &xy);
00360 
00361 
00362     //-----------------------------------------------------------
00363     // Shear the matrix in x for each y coord. by given factor xy
00364     // and shear y for each x coord. by given factor yx
00365     //-----------------------------------------------------------
00366 
00367     template <class S>
00368     const Matrix33 &    shear (const Vec2<S> &h);
00369 
00370 
00371     //-------------------------------------------------
00372     // Limitations of type T (see also class limits<T>)
00373     //-------------------------------------------------
00374 
00375     static T            baseTypeMin()           {return limits<T>::min();}
00376     static T            baseTypeMax()           {return limits<T>::max();}
00377     static T            baseTypeSmallest()      {return limits<T>::smallest();}
00378     static T            baseTypeEpsilon()       {return limits<T>::epsilon();}
00379 
00380   private:
00381 
00382     template <typename R, typename S>
00383     struct isSameType
00384     {
00385         enum {value = 0};
00386     };
00387 
00388     template <typename R>
00389     struct isSameType<R, R>
00390     {
00391         enum {value = 1};
00392     };
00393 };
00394 
00395 
00396 template <class T> class Matrix44
00397 {
00398   public:
00399 
00400     //-------------------
00401     // Access to elements
00402     //-------------------
00403 
00404     T           x[4][4];
00405 
00406     T *         operator [] (int i);
00407     const T *   operator [] (int i) const;
00408 
00409 
00410     //-------------
00411     // Constructors
00412     //-------------
00413 
00414     Matrix44 (Uninitialized) {}
00415 
00416     Matrix44 ();
00417                                 // 1 0 0 0
00418                                 // 0 1 0 0
00419                                 // 0 0 1 0
00420                                 // 0 0 0 1
00421 
00422     Matrix44 (T a);
00423                                 // a a a a
00424                                 // a a a a
00425                                 // a a a a
00426                                 // a a a a
00427 
00428     Matrix44 (const T a[4][4]) ;
00429                                 // a[0][0] a[0][1] a[0][2] a[0][3]
00430                                 // a[1][0] a[1][1] a[1][2] a[1][3]
00431                                 // a[2][0] a[2][1] a[2][2] a[2][3]
00432                                 // a[3][0] a[3][1] a[3][2] a[3][3]
00433 
00434     Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,
00435               T i, T j, T k, T l, T m, T n, T o, T p);
00436 
00437                                 // a b c d
00438                                 // e f g h
00439                                 // i j k l
00440                                 // m n o p
00441 
00442     Matrix44 (Matrix33<T> r, Vec3<T> t);
00443                                 // r r r 0
00444                                 // r r r 0
00445                                 // r r r 0
00446                                 // t t t 1
00447 
00448 
00449     //--------------------------------
00450     // Copy constructor and assignment
00451     //--------------------------------
00452 
00453     Matrix44 (const Matrix44 &v);
00454     template <class S> explicit Matrix44 (const Matrix44<S> &v);
00455 
00456     const Matrix44 &    operator = (const Matrix44 &v);
00457     const Matrix44 &    operator = (T a);
00458 
00459 
00460     //----------------------
00461     // Compatibility with Sb
00462     //----------------------
00463     
00464     T *                 getValue ();
00465     const T *           getValue () const;
00466 
00467     template <class S>
00468     void                getValue (Matrix44<S> &v) const;
00469     template <class S>
00470     Matrix44 &          setValue (const Matrix44<S> &v);
00471 
00472     template <class S>
00473     Matrix44 &          setTheMatrix (const Matrix44<S> &v);
00474 
00475     //---------
00476     // Identity
00477     //---------
00478 
00479     void                makeIdentity();
00480 
00481 
00482     //---------
00483     // Equality
00484     //---------
00485 
00486     bool                operator == (const Matrix44 &v) const;
00487     bool                operator != (const Matrix44 &v) const;
00488 
00489     //-----------------------------------------------------------------------
00490     // Compare two matrices and test if they are "approximately equal":
00491     //
00492     // equalWithAbsError (m, e)
00493     //
00494     //      Returns true if the coefficients of this and m are the same with
00495     //      an absolute error of no more than e, i.e., for all i, j
00496     //
00497     //      abs (this[i][j] - m[i][j]) <= e
00498     //
00499     // equalWithRelError (m, e)
00500     //
00501     //      Returns true if the coefficients of this and m are the same with
00502     //      a relative error of no more than e, i.e., for all i, j
00503     //
00504     //      abs (this[i] - v[i][j]) <= e * abs (this[i][j])
00505     //-----------------------------------------------------------------------
00506 
00507     bool                equalWithAbsError (const Matrix44<T> &v, T e) const;
00508     bool                equalWithRelError (const Matrix44<T> &v, T e) const;
00509 
00510 
00511     //------------------------
00512     // Component-wise addition
00513     //------------------------
00514 
00515     const Matrix44 &    operator += (const Matrix44 &v);
00516     const Matrix44 &    operator += (T a);
00517     Matrix44            operator + (const Matrix44 &v) const;
00518 
00519 
00520     //---------------------------
00521     // Component-wise subtraction
00522     //---------------------------
00523 
00524     const Matrix44 &    operator -= (const Matrix44 &v);
00525     const Matrix44 &    operator -= (T a);
00526     Matrix44            operator - (const Matrix44 &v) const;
00527 
00528 
00529     //------------------------------------
00530     // Component-wise multiplication by -1
00531     //------------------------------------
00532 
00533     Matrix44            operator - () const;
00534     const Matrix44 &    negate ();
00535 
00536 
00537     //------------------------------
00538     // Component-wise multiplication
00539     //------------------------------
00540 
00541     const Matrix44 &    operator *= (T a);
00542     Matrix44            operator * (T a) const;
00543 
00544 
00545     //-----------------------------------
00546     // Matrix-times-matrix multiplication
00547     //-----------------------------------
00548 
00549     const Matrix44 &    operator *= (const Matrix44 &v);
00550     Matrix44            operator * (const Matrix44 &v) const;
00551 
00552     static void         multiply (const Matrix44 &a,    // assumes that
00553                                   const Matrix44 &b,    // &a != &c and
00554                                   Matrix44 &c);         // &b != &c.
00555 
00556 
00557     //-----------------------------------------------------------------
00558     // Vector-times-matrix multiplication; see also the "operator *"
00559     // functions defined below.
00560     //
00561     // m.multVecMatrix(src,dst) implements a homogeneous transformation
00562     // by computing Vec4 (src.x, src.y, src.z, 1) * m and dividing by
00563     // the result's third element.
00564     //
00565     // m.multDirMatrix(src,dst) multiplies src by the upper left 3x3
00566     // submatrix, ignoring the rest of matrix m.
00567     //-----------------------------------------------------------------
00568 
00569     template <class S>
00570     void                multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const;
00571 
00572     template <class S>
00573     void                multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const;
00574 
00575 
00576     //------------------------
00577     // Component-wise division
00578     //------------------------
00579 
00580     const Matrix44 &    operator /= (T a);
00581     Matrix44            operator / (T a) const;
00582 
00583 
00584     //------------------
00585     // Transposed matrix
00586     //------------------
00587 
00588     const Matrix44 &    transpose ();
00589     Matrix44            transposed () const;
00590 
00591 
00592     //------------------------------------------------------------
00593     // Inverse matrix: If singExc is false, inverting a singular
00594     // matrix produces an identity matrix.  If singExc is true,
00595     // inverting a singular matrix throws a SingMatrixExc.
00596     //
00597     // inverse() and invert() invert matrices using determinants;
00598     // gjInverse() and gjInvert() use the Gauss-Jordan method.
00599     //
00600     // inverse() and invert() are significantly faster than
00601     // gjInverse() and gjInvert(), but the results may be slightly
00602     // less accurate.
00603     // 
00604     //------------------------------------------------------------
00605 
00606     const Matrix44 &    invert (bool singExc = false)
00607                         throw (Iex::MathExc);
00608 
00609     Matrix44<T>         inverse (bool singExc = false) const
00610                         throw (Iex::MathExc);
00611 
00612     const Matrix44 &    gjInvert (bool singExc = false)
00613                         throw (Iex::MathExc);
00614 
00615     Matrix44<T>         gjInverse (bool singExc = false) const
00616                         throw (Iex::MathExc);
00617 
00618 
00619     //--------------------------------------------------------
00620     // Set matrix to rotation by XYZ euler angles (in radians)
00621     //--------------------------------------------------------
00622 
00623     template <class S>
00624     const Matrix44 &    setEulerAngles (const Vec3<S>& r);
00625 
00626 
00627     //--------------------------------------------------------
00628     // Set matrix to rotation around given axis by given angle
00629     //--------------------------------------------------------
00630 
00631     template <class S>
00632     const Matrix44 &    setAxisAngle (const Vec3<S>& ax, S ang);
00633 
00634 
00635     //-------------------------------------------
00636     // Rotate the matrix by XYZ euler angles in r
00637     //-------------------------------------------
00638 
00639     template <class S>
00640     const Matrix44 &    rotate (const Vec3<S> &r);
00641 
00642 
00643     //--------------------------------------------
00644     // Set matrix to scale by given uniform factor
00645     //--------------------------------------------
00646 
00647     const Matrix44 &    setScale (T s);
00648 
00649 
00650     //------------------------------------
00651     // Set matrix to scale by given vector
00652     //------------------------------------
00653 
00654     template <class S>
00655     const Matrix44 &    setScale (const Vec3<S> &s);
00656 
00657 
00658     //----------------------
00659     // Scale the matrix by s
00660     //----------------------
00661 
00662     template <class S>
00663     const Matrix44 &    scale (const Vec3<S> &s);
00664 
00665 
00666     //------------------------------------------
00667     // Set matrix to translation by given vector
00668     //------------------------------------------
00669 
00670     template <class S>
00671     const Matrix44 &    setTranslation (const Vec3<S> &t);
00672 
00673 
00674     //-----------------------------
00675     // Return translation component
00676     //-----------------------------
00677 
00678     const Vec3<T>       translation () const;
00679 
00680 
00681     //--------------------------
00682     // Translate the matrix by t
00683     //--------------------------
00684 
00685     template <class S>
00686     const Matrix44 &    translate (const Vec3<S> &t);
00687 
00688 
00689     //-------------------------------------------------------------
00690     // Set matrix to shear by given vector h.  The resulting matrix
00691     //    will shear x for each y coord. by a factor of h[0] ;
00692     //    will shear x for each z coord. by a factor of h[1] ;
00693     //    will shear y for each z coord. by a factor of h[2] .
00694     //-------------------------------------------------------------
00695 
00696     template <class S>
00697     const Matrix44 &    setShear (const Vec3<S> &h);
00698 
00699 
00700     //------------------------------------------------------------
00701     // Set matrix to shear by given factors.  The resulting matrix
00702     //    will shear x for each y coord. by a factor of h.xy ;
00703     //    will shear x for each z coord. by a factor of h.xz ;
00704     //    will shear y for each z coord. by a factor of h.yz ; 
00705     //    will shear y for each x coord. by a factor of h.yx ;
00706     //    will shear z for each x coord. by a factor of h.zx ;
00707     //    will shear z for each y coord. by a factor of h.zy .
00708     //------------------------------------------------------------
00709 
00710     template <class S>
00711     const Matrix44 &    setShear (const Shear6<S> &h);
00712 
00713 
00714     //--------------------------------------------------------
00715     // Shear the matrix by given vector.  The composed matrix 
00716     // will be <shear> * <this>, where the shear matrix ...
00717     //    will shear x for each y coord. by a factor of h[0] ;
00718     //    will shear x for each z coord. by a factor of h[1] ;
00719     //    will shear y for each z coord. by a factor of h[2] .
00720     //--------------------------------------------------------
00721 
00722     template <class S>
00723     const Matrix44 &    shear (const Vec3<S> &h);
00724 
00725 
00726     //------------------------------------------------------------
00727     // Shear the matrix by the given factors.  The composed matrix 
00728     // will be <shear> * <this>, where the shear matrix ...
00729     //    will shear x for each y coord. by a factor of h.xy ;
00730     //    will shear x for each z coord. by a factor of h.xz ;
00731     //    will shear y for each z coord. by a factor of h.yz ;
00732     //    will shear y for each x coord. by a factor of h.yx ;
00733     //    will shear z for each x coord. by a factor of h.zx ;
00734     //    will shear z for each y coord. by a factor of h.zy .
00735     //------------------------------------------------------------
00736 
00737     template <class S>
00738     const Matrix44 &    shear (const Shear6<S> &h);
00739 
00740 
00741     //-------------------------------------------------
00742     // Limitations of type T (see also class limits<T>)
00743     //-------------------------------------------------
00744 
00745     static T            baseTypeMin()           {return limits<T>::min();}
00746     static T            baseTypeMax()           {return limits<T>::max();}
00747     static T            baseTypeSmallest()      {return limits<T>::smallest();}
00748     static T            baseTypeEpsilon()       {return limits<T>::epsilon();}
00749 
00750   private:
00751 
00752     template <typename R, typename S>
00753     struct isSameType
00754     {
00755         enum {value = 0};
00756     };
00757 
00758     template <typename R>
00759     struct isSameType<R, R>
00760     {
00761         enum {value = 1};
00762     };
00763 };
00764 
00765 
00766 //--------------
00767 // Stream output
00768 //--------------
00769 
00770 template <class T>
00771 std::ostream &  operator << (std::ostream & s, const Matrix33<T> &m); 
00772 
00773 template <class T>
00774 std::ostream &  operator << (std::ostream & s, const Matrix44<T> &m); 
00775 
00776 
00777 //---------------------------------------------
00778 // Vector-times-matrix multiplication operators
00779 //---------------------------------------------
00780 
00781 template <class S, class T>
00782 const Vec2<S> &            operator *= (Vec2<S> &v, const Matrix33<T> &m);
00783 
00784 template <class S, class T>
00785 Vec2<S>                    operator * (const Vec2<S> &v, const Matrix33<T> &m);
00786 
00787 template <class S, class T>
00788 const Vec3<S> &            operator *= (Vec3<S> &v, const Matrix33<T> &m);
00789 
00790 template <class S, class T>
00791 Vec3<S>                    operator * (const Vec3<S> &v, const Matrix33<T> &m);
00792 
00793 template <class S, class T>
00794 const Vec3<S> &            operator *= (Vec3<S> &v, const Matrix44<T> &m);
00795 
00796 template <class S, class T>
00797 Vec3<S>                    operator * (const Vec3<S> &v, const Matrix44<T> &m);
00798 
00799 template <class S, class T>
00800 const Vec4<S> &            operator *= (Vec4<S> &v, const Matrix44<T> &m);
00801 
00802 template <class S, class T>
00803 Vec4<S>                    operator * (const Vec4<S> &v, const Matrix44<T> &m);
00804 
00805 //-------------------------
00806 // Typedefs for convenience
00807 //-------------------------
00808 
00809 typedef Matrix33 <float>  M33f;
00810 typedef Matrix33 <double> M33d;
00811 typedef Matrix44 <float>  M44f;
00812 typedef Matrix44 <double> M44d;
00813 
00814 
00815 //---------------------------
00816 // Implementation of Matrix33
00817 //---------------------------
00818 
00819 template <class T>
00820 inline T *
00821 Matrix33<T>::operator [] (int i)
00822 {
00823     return x[i];
00824 }
00825 
00826 template <class T>
00827 inline const T *
00828 Matrix33<T>::operator [] (int i) const
00829 {
00830     return x[i];
00831 }
00832 
00833 template <class T>
00834 inline
00835 Matrix33<T>::Matrix33 ()
00836 {
00837     memset (x, 0, sizeof (x));
00838     x[0][0] = 1;
00839     x[1][1] = 1;
00840     x[2][2] = 1;
00841 }
00842 
00843 template <class T>
00844 inline
00845 Matrix33<T>::Matrix33 (T a)
00846 {
00847     x[0][0] = a;
00848     x[0][1] = a;
00849     x[0][2] = a;
00850     x[1][0] = a;
00851     x[1][1] = a;
00852     x[1][2] = a;
00853     x[2][0] = a;
00854     x[2][1] = a;
00855     x[2][2] = a;
00856 }
00857 
00858 template <class T>
00859 inline
00860 Matrix33<T>::Matrix33 (const T a[3][3]) 
00861 {
00862     memcpy (x, a, sizeof (x));
00863 }
00864 
00865 template <class T>
00866 inline
00867 Matrix33<T>::Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i)
00868 {
00869     x[0][0] = a;
00870     x[0][1] = b;
00871     x[0][2] = c;
00872     x[1][0] = d;
00873     x[1][1] = e;
00874     x[1][2] = f;
00875     x[2][0] = g;
00876     x[2][1] = h;
00877     x[2][2] = i;
00878 }
00879 
00880 template <class T>
00881 inline
00882 Matrix33<T>::Matrix33 (const Matrix33 &v)
00883 {
00884     memcpy (x, v.x, sizeof (x));
00885 }
00886 
00887 template <class T>
00888 template <class S>
00889 inline
00890 Matrix33<T>::Matrix33 (const Matrix33<S> &v)
00891 {
00892     x[0][0] = T (v.x[0][0]);
00893     x[0][1] = T (v.x[0][1]);
00894     x[0][2] = T (v.x[0][2]);
00895     x[1][0] = T (v.x[1][0]);
00896     x[1][1] = T (v.x[1][1]);
00897     x[1][2] = T (v.x[1][2]);
00898     x[2][0] = T (v.x[2][0]);
00899     x[2][1] = T (v.x[2][1]);
00900     x[2][2] = T (v.x[2][2]);
00901 }
00902 
00903 template <class T>
00904 inline const Matrix33<T> &
00905 Matrix33<T>::operator = (const Matrix33 &v)
00906 {
00907     memcpy (x, v.x, sizeof (x));
00908     return *this;
00909 }
00910 
00911 template <class T>
00912 inline const Matrix33<T> &
00913 Matrix33<T>::operator = (T a)
00914 {
00915     x[0][0] = a;
00916     x[0][1] = a;
00917     x[0][2] = a;
00918     x[1][0] = a;
00919     x[1][1] = a;
00920     x[1][2] = a;
00921     x[2][0] = a;
00922     x[2][1] = a;
00923     x[2][2] = a;
00924     return *this;
00925 }
00926 
00927 template <class T>
00928 inline T *
00929 Matrix33<T>::getValue ()
00930 {
00931     return (T *) &x[0][0];
00932 }
00933 
00934 template <class T>
00935 inline const T *
00936 Matrix33<T>::getValue () const
00937 {
00938     return (const T *) &x[0][0];
00939 }
00940 
00941 template <class T>
00942 template <class S>
00943 inline void
00944 Matrix33<T>::getValue (Matrix33<S> &v) const
00945 {
00946     if (isSameType<S,T>::value)
00947     {
00948         memcpy (v.x, x, sizeof (x));
00949     }
00950     else
00951     {
00952         v.x[0][0] = x[0][0];
00953         v.x[0][1] = x[0][1];
00954         v.x[0][2] = x[0][2];
00955         v.x[1][0] = x[1][0];
00956         v.x[1][1] = x[1][1];
00957         v.x[1][2] = x[1][2];
00958         v.x[2][0] = x[2][0];
00959         v.x[2][1] = x[2][1];
00960         v.x[2][2] = x[2][2];
00961     }
00962 }
00963 
00964 template <class T>
00965 template <class S>
00966 inline Matrix33<T> &
00967 Matrix33<T>::setValue (const Matrix33<S> &v)
00968 {
00969     if (isSameType<S,T>::value)
00970     {
00971         memcpy (x, v.x, sizeof (x));
00972     }
00973     else
00974     {
00975         x[0][0] = v.x[0][0];
00976         x[0][1] = v.x[0][1];
00977         x[0][2] = v.x[0][2];
00978         x[1][0] = v.x[1][0];
00979         x[1][1] = v.x[1][1];
00980         x[1][2] = v.x[1][2];
00981         x[2][0] = v.x[2][0];
00982         x[2][1] = v.x[2][1];
00983         x[2][2] = v.x[2][2];
00984     }
00985 
00986     return *this;
00987 }
00988 
00989 template <class T>
00990 template <class S>
00991 inline Matrix33<T> &
00992 Matrix33<T>::setTheMatrix (const Matrix33<S> &v)
00993 {
00994     if (isSameType<S,T>::value)
00995     {
00996         memcpy (x, v.x, sizeof (x));
00997     }
00998     else
00999     {
01000         x[0][0] = v.x[0][0];
01001         x[0][1] = v.x[0][1];
01002         x[0][2] = v.x[0][2];
01003         x[1][0] = v.x[1][0];
01004         x[1][1] = v.x[1][1];
01005         x[1][2] = v.x[1][2];
01006         x[2][0] = v.x[2][0];
01007         x[2][1] = v.x[2][1];
01008         x[2][2] = v.x[2][2];
01009     }
01010 
01011     return *this;
01012 }
01013 
01014 template <class T>
01015 inline void
01016 Matrix33<T>::makeIdentity()
01017 {
01018     memset (x, 0, sizeof (x));
01019     x[0][0] = 1;
01020     x[1][1] = 1;
01021     x[2][2] = 1;
01022 }
01023 
01024 template <class T>
01025 bool
01026 Matrix33<T>::operator == (const Matrix33 &v) const
01027 {
01028     return x[0][0] == v.x[0][0] &&
01029            x[0][1] == v.x[0][1] &&
01030            x[0][2] == v.x[0][2] &&
01031            x[1][0] == v.x[1][0] &&
01032            x[1][1] == v.x[1][1] &&
01033            x[1][2] == v.x[1][2] &&
01034            x[2][0] == v.x[2][0] &&
01035            x[2][1] == v.x[2][1] &&
01036            x[2][2] == v.x[2][2];
01037 }
01038 
01039 template <class T>
01040 bool
01041 Matrix33<T>::operator != (const Matrix33 &v) const
01042 {
01043     return x[0][0] != v.x[0][0] ||
01044            x[0][1] != v.x[0][1] ||
01045            x[0][2] != v.x[0][2] ||
01046            x[1][0] != v.x[1][0] ||
01047            x[1][1] != v.x[1][1] ||
01048            x[1][2] != v.x[1][2] ||
01049            x[2][0] != v.x[2][0] ||
01050            x[2][1] != v.x[2][1] ||
01051            x[2][2] != v.x[2][2];
01052 }
01053 
01054 template <class T>
01055 bool
01056 Matrix33<T>::equalWithAbsError (const Matrix33<T> &m, T e) const
01057 {
01058     for (int i = 0; i < 3; i++)
01059         for (int j = 0; j < 3; j++)
01060             if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e))
01061                 return false;
01062 
01063     return true;
01064 }
01065 
01066 template <class T>
01067 bool
01068 Matrix33<T>::equalWithRelError (const Matrix33<T> &m, T e) const
01069 {
01070     for (int i = 0; i < 3; i++)
01071         for (int j = 0; j < 3; j++)
01072             if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e))
01073                 return false;
01074 
01075     return true;
01076 }
01077 
01078 template <class T>
01079 const Matrix33<T> &
01080 Matrix33<T>::operator += (const Matrix33<T> &v)
01081 {
01082     x[0][0] += v.x[0][0];
01083     x[0][1] += v.x[0][1];
01084     x[0][2] += v.x[0][2];
01085     x[1][0] += v.x[1][0];
01086     x[1][1] += v.x[1][1];
01087     x[1][2] += v.x[1][2];
01088     x[2][0] += v.x[2][0];
01089     x[2][1] += v.x[2][1];
01090     x[2][2] += v.x[2][2];
01091 
01092     return *this;
01093 }
01094 
01095 template <class T>
01096 const Matrix33<T> &
01097 Matrix33<T>::operator += (T a)
01098 {
01099     x[0][0] += a;
01100     x[0][1] += a;
01101     x[0][2] += a;
01102     x[1][0] += a;
01103     x[1][1] += a;
01104     x[1][2] += a;
01105     x[2][0] += a;
01106     x[2][1] += a;
01107     x[2][2] += a;
01108   
01109     return *this;
01110 }
01111 
01112 template <class T>
01113 Matrix33<T>
01114 Matrix33<T>::operator + (const Matrix33<T> &v) const
01115 {
01116     return Matrix33 (x[0][0] + v.x[0][0],
01117                      x[0][1] + v.x[0][1],
01118                      x[0][2] + v.x[0][2],
01119                      x[1][0] + v.x[1][0],
01120                      x[1][1] + v.x[1][1],
01121                      x[1][2] + v.x[1][2],
01122                      x[2][0] + v.x[2][0],
01123                      x[2][1] + v.x[2][1],
01124                      x[2][2] + v.x[2][2]);
01125 }
01126 
01127 template <class T>
01128 const Matrix33<T> &
01129 Matrix33<T>::operator -= (const Matrix33<T> &v)
01130 {
01131     x[0][0] -= v.x[0][0];
01132     x[0][1] -= v.x[0][1];
01133     x[0][2] -= v.x[0][2];
01134     x[1][0] -= v.x[1][0];
01135     x[1][1] -= v.x[1][1];
01136     x[1][2] -= v.x[1][2];
01137     x[2][0] -= v.x[2][0];
01138     x[2][1] -= v.x[2][1];
01139     x[2][2] -= v.x[2][2];
01140   
01141     return *this;
01142 }
01143 
01144 template <class T>
01145 const Matrix33<T> &
01146 Matrix33<T>::operator -= (T a)
01147 {
01148     x[0][0] -= a;
01149     x[0][1] -= a;
01150     x[0][2] -= a;
01151     x[1][0] -= a;
01152     x[1][1] -= a;
01153     x[1][2] -= a;
01154     x[2][0] -= a;
01155     x[2][1] -= a;
01156     x[2][2] -= a;
01157   
01158     return *this;
01159 }
01160 
01161 template <class T>
01162 Matrix33<T>
01163 Matrix33<T>::operator - (const Matrix33<T> &v) const
01164 {
01165     return Matrix33 (x[0][0] - v.x[0][0],
01166                      x[0][1] - v.x[0][1],
01167                      x[0][2] - v.x[0][2],
01168                      x[1][0] - v.x[1][0],
01169                      x[1][1] - v.x[1][1],
01170                      x[1][2] - v.x[1][2],
01171                      x[2][0] - v.x[2][0],
01172                      x[2][1] - v.x[2][1],
01173                      x[2][2] - v.x[2][2]);
01174 }
01175 
01176 template <class T>
01177 Matrix33<T>
01178 Matrix33<T>::operator - () const
01179 {
01180     return Matrix33 (-x[0][0],
01181                      -x[0][1],
01182                      -x[0][2],
01183                      -x[1][0],
01184                      -x[1][1],
01185                      -x[1][2],
01186                      -x[2][0],
01187                      -x[2][1],
01188                      -x[2][2]);
01189 }
01190 
01191 template <class T>
01192 const Matrix33<T> &
01193 Matrix33<T>::negate ()
01194 {
01195     x[0][0] = -x[0][0];
01196     x[0][1] = -x[0][1];
01197     x[0][2] = -x[0][2];
01198     x[1][0] = -x[1][0];
01199     x[1][1] = -x[1][1];
01200     x[1][2] = -x[1][2];
01201     x[2][0] = -x[2][0];
01202     x[2][1] = -x[2][1];
01203     x[2][2] = -x[2][2];
01204 
01205     return *this;
01206 }
01207 
01208 template <class T>
01209 const Matrix33<T> &
01210 Matrix33<T>::operator *= (T a)
01211 {
01212     x[0][0] *= a;
01213     x[0][1] *= a;
01214     x[0][2] *= a;
01215     x[1][0] *= a;
01216     x[1][1] *= a;
01217     x[1][2] *= a;
01218     x[2][0] *= a;
01219     x[2][1] *= a;
01220     x[2][2] *= a;
01221   
01222     return *this;
01223 }
01224 
01225 template <class T>
01226 Matrix33<T>
01227 Matrix33<T>::operator * (T a) const
01228 {
01229     return Matrix33 (x[0][0] * a,
01230                      x[0][1] * a,
01231                      x[0][2] * a,
01232                      x[1][0] * a,
01233                      x[1][1] * a,
01234                      x[1][2] * a,
01235                      x[2][0] * a,
01236                      x[2][1] * a,
01237                      x[2][2] * a);
01238 }
01239 
01240 template <class T>
01241 inline Matrix33<T>
01242 operator * (T a, const Matrix33<T> &v)
01243 {
01244     return v * a;
01245 }
01246 
01247 template <class T>
01248 const Matrix33<T> &
01249 Matrix33<T>::operator *= (const Matrix33<T> &v)
01250 {
01251     Matrix33 tmp (T (0));
01252 
01253     for (int i = 0; i < 3; i++)
01254         for (int j = 0; j < 3; j++)
01255             for (int k = 0; k < 3; k++)
01256                 tmp.x[i][j] += x[i][k] * v.x[k][j];
01257 
01258     *this = tmp;
01259     return *this;
01260 }
01261 
01262 template <class T>
01263 Matrix33<T>
01264 Matrix33<T>::operator * (const Matrix33<T> &v) const
01265 {
01266     Matrix33 tmp (T (0));
01267 
01268     for (int i = 0; i < 3; i++)
01269         for (int j = 0; j < 3; j++)
01270             for (int k = 0; k < 3; k++)
01271                 tmp.x[i][j] += x[i][k] * v.x[k][j];
01272 
01273     return tmp;
01274 }
01275 
01276 template <class T>
01277 template <class S>
01278 void
01279 Matrix33<T>::multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const
01280 {
01281     S a, b, w;
01282 
01283     a = src[0] * x[0][0] + src[1] * x[1][0] + x[2][0];
01284     b = src[0] * x[0][1] + src[1] * x[1][1] + x[2][1];
01285     w = src[0] * x[0][2] + src[1] * x[1][2] + x[2][2];
01286 
01287     dst.x = a / w;
01288     dst.y = b / w;
01289 }
01290 
01291 template <class T>
01292 template <class S>
01293 void
01294 Matrix33<T>::multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const
01295 {
01296     S a, b;
01297 
01298     a = src[0] * x[0][0] + src[1] * x[1][0];
01299     b = src[0] * x[0][1] + src[1] * x[1][1];
01300 
01301     dst.x = a;
01302     dst.y = b;
01303 }
01304 
01305 template <class T>
01306 const Matrix33<T> &
01307 Matrix33<T>::operator /= (T a)
01308 {
01309     x[0][0] /= a;
01310     x[0][1] /= a;
01311     x[0][2] /= a;
01312     x[1][0] /= a;
01313     x[1][1] /= a;
01314     x[1][2] /= a;
01315     x[2][0] /= a;
01316     x[2][1] /= a;
01317     x[2][2] /= a;
01318   
01319     return *this;
01320 }
01321 
01322 template <class T>
01323 Matrix33<T>
01324 Matrix33<T>::operator / (T a) const
01325 {
01326     return Matrix33 (x[0][0] / a,
01327                      x[0][1] / a,
01328                      x[0][2] / a,
01329                      x[1][0] / a,
01330                      x[1][1] / a,
01331                      x[1][2] / a,
01332                      x[2][0] / a,
01333                      x[2][1] / a,
01334                      x[2][2] / a);
01335 }
01336 
01337 template <class T>
01338 const Matrix33<T> &
01339 Matrix33<T>::transpose ()
01340 {
01341     Matrix33 tmp (x[0][0],
01342                   x[1][0],
01343                   x[2][0],
01344                   x[0][1],
01345                   x[1][1],
01346                   x[2][1],
01347                   x[0][2],
01348                   x[1][2],
01349                   x[2][2]);
01350     *this = tmp;
01351     return *this;
01352 }
01353 
01354 template <class T>
01355 Matrix33<T>
01356 Matrix33<T>::transposed () const
01357 {
01358     return Matrix33 (x[0][0],
01359                      x[1][0],
01360                      x[2][0],
01361                      x[0][1],
01362                      x[1][1],
01363                      x[2][1],
01364                      x[0][2],
01365                      x[1][2],
01366                      x[2][2]);
01367 }
01368 
01369 template <class T>
01370 const Matrix33<T> &
01371 Matrix33<T>::gjInvert (bool singExc) throw (Iex::MathExc)
01372 {
01373     *this = gjInverse (singExc);
01374     return *this;
01375 }
01376 
01377 template <class T>
01378 Matrix33<T>
01379 Matrix33<T>::gjInverse (bool singExc) const throw (Iex::MathExc)
01380 {
01381     int i, j, k;
01382     Matrix33 s;
01383     Matrix33 t (*this);
01384 
01385     // Forward elimination
01386 
01387     for (i = 0; i < 2 ; i++)
01388     {
01389         int pivot = i;
01390 
01391         T pivotsize = t[i][i];
01392 
01393         if (pivotsize < 0)
01394             pivotsize = -pivotsize;
01395 
01396         for (j = i + 1; j < 3; j++)
01397         {
01398             T tmp = t[j][i];
01399 
01400             if (tmp < 0)
01401                 tmp = -tmp;
01402 
01403             if (tmp > pivotsize)
01404             {
01405                 pivot = j;
01406                 pivotsize = tmp;
01407             }
01408         }
01409 
01410         if (pivotsize == 0)
01411         {
01412             if (singExc)
01413                 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
01414 
01415             return Matrix33();
01416         }
01417 
01418         if (pivot != i)
01419         {
01420             for (j = 0; j < 3; j++)
01421             {
01422                 T tmp;
01423 
01424                 tmp = t[i][j];
01425                 t[i][j] = t[pivot][j];
01426                 t[pivot][j] = tmp;
01427 
01428                 tmp = s[i][j];
01429                 s[i][j] = s[pivot][j];
01430                 s[pivot][j] = tmp;
01431             }
01432         }
01433 
01434         for (j = i + 1; j < 3; j++)
01435         {
01436             T f = t[j][i] / t[i][i];
01437 
01438             for (k = 0; k < 3; k++)
01439             {
01440                 t[j][k] -= f * t[i][k];
01441                 s[j][k] -= f * s[i][k];
01442             }
01443         }
01444     }
01445 
01446     // Backward substitution
01447 
01448     for (i = 2; i >= 0; --i)
01449     {
01450         T f;
01451 
01452         if ((f = t[i][i]) == 0)
01453         {
01454             if (singExc)
01455                 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
01456 
01457             return Matrix33();
01458         }
01459 
01460         for (j = 0; j < 3; j++)
01461         {
01462             t[i][j] /= f;
01463             s[i][j] /= f;
01464         }
01465 
01466         for (j = 0; j < i; j++)
01467         {
01468             f = t[j][i];
01469 
01470             for (k = 0; k < 3; k++)
01471             {
01472                 t[j][k] -= f * t[i][k];
01473                 s[j][k] -= f * s[i][k];
01474             }
01475         }
01476     }
01477 
01478     return s;
01479 }
01480 
01481 template <class T>
01482 const Matrix33<T> &
01483 Matrix33<T>::invert (bool singExc) throw (Iex::MathExc)
01484 {
01485     *this = inverse (singExc);
01486     return *this;
01487 }
01488 
01489 template <class T>
01490 Matrix33<T>
01491 Matrix33<T>::inverse (bool singExc) const throw (Iex::MathExc)
01492 {
01493     if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1)
01494     {
01495         Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
01496                     x[2][1] * x[0][2] - x[0][1] * x[2][2],
01497                     x[0][1] * x[1][2] - x[1][1] * x[0][2],
01498 
01499                     x[2][0] * x[1][2] - x[1][0] * x[2][2],
01500                     x[0][0] * x[2][2] - x[2][0] * x[0][2],
01501                     x[1][0] * x[0][2] - x[0][0] * x[1][2],
01502 
01503                     x[1][0] * x[2][1] - x[2][0] * x[1][1],
01504                     x[2][0] * x[0][1] - x[0][0] * x[2][1],
01505                     x[0][0] * x[1][1] - x[1][0] * x[0][1]);
01506 
01507         T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
01508 
01509         if (Imath::abs (r) >= 1)
01510         {
01511             for (int i = 0; i < 3; ++i)
01512             {
01513                 for (int j = 0; j < 3; ++j)
01514                 {
01515                     s[i][j] /= r;
01516                 }
01517             }
01518         }
01519         else
01520         {
01521             T mr = Imath::abs (r) / limits<T>::smallest();
01522 
01523             for (int i = 0; i < 3; ++i)
01524             {
01525                 for (int j = 0; j < 3; ++j)
01526                 {
01527                     if (mr > Imath::abs (s[i][j]))
01528                     {
01529                         s[i][j] /= r;
01530                     }
01531                     else
01532                     {
01533                         if (singExc)
01534                             throw SingMatrixExc ("Cannot invert "
01535                                                  "singular matrix.");
01536                         return Matrix33();
01537                     }
01538                 }
01539             }
01540         }
01541 
01542         return s;
01543     }
01544     else
01545     {
01546         Matrix33 s ( x[1][1],
01547                     -x[0][1],
01548                      0, 
01549 
01550                     -x[1][0],
01551                      x[0][0],
01552                      0,
01553 
01554                      0,
01555                      0,
01556                      1);
01557 
01558         T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
01559 
01560         if (Imath::abs (r) >= 1)
01561         {
01562             for (int i = 0; i < 2; ++i)
01563             {
01564                 for (int j = 0; j < 2; ++j)
01565                 {
01566                     s[i][j] /= r;
01567                 }
01568             }
01569         }
01570         else
01571         {
01572             T mr = Imath::abs (r) / limits<T>::smallest();
01573 
01574             for (int i = 0; i < 2; ++i)
01575             {
01576                 for (int j = 0; j < 2; ++j)
01577                 {
01578                     if (mr > Imath::abs (s[i][j]))
01579                     {
01580                         s[i][j] /= r;
01581                     }
01582                     else
01583                     {
01584                         if (singExc)
01585                             throw SingMatrixExc ("Cannot invert "
01586                                                  "singular matrix.");
01587                         return Matrix33();
01588                     }
01589                 }
01590             }
01591         }
01592 
01593         s[2][0] = -x[2][0] * s[0][0] - x[2][1] * s[1][0];
01594         s[2][1] = -x[2][0] * s[0][1] - x[2][1] * s[1][1];
01595 
01596         return s;
01597     }
01598 }
01599 
01600 template <class T>
01601 template <class S>
01602 const Matrix33<T> &
01603 Matrix33<T>::setRotation (S r)
01604 {
01605     S cos_r, sin_r;
01606 
01607     cos_r = Math<T>::cos (r);
01608     sin_r = Math<T>::sin (r);
01609 
01610     x[0][0] =  cos_r;
01611     x[0][1] =  sin_r;
01612     x[0][2] =  0;
01613 
01614     x[1][0] =  -sin_r;
01615     x[1][1] =  cos_r;
01616     x[1][2] =  0;
01617 
01618     x[2][0] =  0;
01619     x[2][1] =  0;
01620     x[2][2] =  1;
01621 
01622     return *this;
01623 }
01624 
01625 template <class T>
01626 template <class S>
01627 const Matrix33<T> &
01628 Matrix33<T>::rotate (S r)
01629 {
01630     *this *= Matrix33<T>().setRotation (r);
01631     return *this;
01632 }
01633 
01634 template <class T>
01635 const Matrix33<T> &
01636 Matrix33<T>::setScale (T s)
01637 {
01638     memset (x, 0, sizeof (x));
01639     x[0][0] = s;
01640     x[1][1] = s;
01641     x[2][2] = 1;
01642 
01643     return *this;
01644 }
01645 
01646 template <class T>
01647 template <class S>
01648 const Matrix33<T> &
01649 Matrix33<T>::setScale (const Vec2<S> &s)
01650 {
01651     memset (x, 0, sizeof (x));
01652     x[0][0] = s[0];
01653     x[1][1] = s[1];
01654     x[2][2] = 1;
01655 
01656     return *this;
01657 }
01658 
01659 template <class T>
01660 template <class S>
01661 const Matrix33<T> &
01662 Matrix33<T>::scale (const Vec2<S> &s)
01663 {
01664     x[0][0] *= s[0];
01665     x[0][1] *= s[0];
01666     x[0][2] *= s[0];
01667 
01668     x[1][0] *= s[1];
01669     x[1][1] *= s[1];
01670     x[1][2] *= s[1];
01671 
01672     return *this;
01673 }
01674 
01675 template <class T>
01676 template <class S>
01677 const Matrix33<T> &
01678 Matrix33<T>::setTranslation (const Vec2<S> &t)
01679 {
01680     x[0][0] = 1;
01681     x[0][1] = 0;
01682     x[0][2] = 0;
01683 
01684     x[1][0] = 0;
01685     x[1][1] = 1;
01686     x[1][2] = 0;
01687 
01688     x[2][0] = t[0];
01689     x[2][1] = t[1];
01690     x[2][2] = 1;
01691 
01692     return *this;
01693 }
01694 
01695 template <class T>
01696 inline Vec2<T> 
01697 Matrix33<T>::translation () const
01698 {
01699     return Vec2<T> (x[2][0], x[2][1]);
01700 }
01701 
01702 template <class T>
01703 template <class S>
01704 const Matrix33<T> &
01705 Matrix33<T>::translate (const Vec2<S> &t)
01706 {
01707     x[2][0] += t[0] * x[0][0] + t[1] * x[1][0];
01708     x[2][1] += t[0] * x[0][1] + t[1] * x[1][1];
01709     x[2][2] += t[0] * x[0][2] + t[1] * x[1][2];
01710 
01711     return *this;
01712 }
01713 
01714 template <class T>
01715 template <class S>
01716 const Matrix33<T> &
01717 Matrix33<T>::setShear (const S &xy)
01718 {
01719     x[0][0] = 1;
01720     x[0][1] = 0;
01721     x[0][2] = 0;
01722 
01723     x[1][0] = xy;
01724     x[1][1] = 1;
01725     x[1][2] = 0;
01726 
01727     x[2][0] = 0;
01728     x[2][1] = 0;
01729     x[2][2] = 1;
01730 
01731     return *this;
01732 }
01733 
01734 template <class T>
01735 template <class S>
01736 const Matrix33<T> &
01737 Matrix33<T>::setShear (const Vec2<S> &h)
01738 {
01739     x[0][0] = 1;
01740     x[0][1] = h[1];
01741     x[0][2] = 0;
01742 
01743     x[1][0] = h[0];
01744     x[1][1] = 1;
01745     x[1][2] = 0;
01746 
01747     x[2][0] = 0;
01748     x[2][1] = 0;
01749     x[2][2] = 1;
01750 
01751     return *this;
01752 }
01753 
01754 template <class T>
01755 template <class S>
01756 const Matrix33<T> &
01757 Matrix33<T>::shear (const S &xy)
01758 {
01759     //
01760     // In this case, we don't need a temp. copy of the matrix 
01761     // because we never use a value on the RHS after we've 
01762     // changed it on the LHS.
01763     // 
01764 
01765     x[1][0] += xy * x[0][0];
01766     x[1][1] += xy * x[0][1];
01767     x[1][2] += xy * x[0][2];
01768 
01769     return *this;
01770 }
01771 
01772 template <class T>
01773 template <class S>
01774 const Matrix33<T> &
01775 Matrix33<T>::shear (const Vec2<S> &h)
01776 {
01777     Matrix33<T> P (*this);
01778     
01779     x[0][0] = P[0][0] + h[1] * P[1][0];
01780     x[0][1] = P[0][1] + h[1] * P[1][1];
01781     x[0][2] = P[0][2] + h[1] * P[1][2];
01782     
01783     x[1][0] = P[1][0] + h[0] * P[0][0];
01784     x[1][1] = P[1][1] + h[0] * P[0][1];
01785     x[1][2] = P[1][2] + h[0] * P[0][2];
01786 
01787     return *this;
01788 }
01789 
01790 
01791 //---------------------------
01792 // Implementation of Matrix44
01793 //---------------------------
01794 
01795 template <class T>
01796 inline T *
01797 Matrix44<T>::operator [] (int i)
01798 {
01799     return x[i];
01800 }
01801 
01802 template <class T>
01803 inline const T *
01804 Matrix44<T>::operator [] (int i) const
01805 {
01806     return x[i];
01807 }
01808 
01809 template <class T>
01810 inline
01811 Matrix44<T>::Matrix44 ()
01812 {
01813     memset (x, 0, sizeof (x));
01814     x[0][0] = 1;
01815     x[1][1] = 1;
01816     x[2][2] = 1;
01817     x[3][3] = 1;
01818 }
01819 
01820 template <class T>
01821 inline
01822 Matrix44<T>::Matrix44 (T a)
01823 {
01824     x[0][0] = a;
01825     x[0][1] = a;
01826     x[0][2] = a;
01827     x[0][3] = a;
01828     x[1][0] = a;
01829     x[1][1] = a;
01830     x[1][2] = a;
01831     x[1][3] = a;
01832     x[2][0] = a;
01833     x[2][1] = a;
01834     x[2][2] = a;
01835     x[2][3] = a;
01836     x[3][0] = a;
01837     x[3][1] = a;
01838     x[3][2] = a;
01839     x[3][3] = a;
01840 }
01841 
01842 template <class T>
01843 inline
01844 Matrix44<T>::Matrix44 (const T a[4][4]) 
01845 {
01846     memcpy (x, a, sizeof (x));
01847 }
01848 
01849 template <class T>
01850 inline
01851 Matrix44<T>::Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,
01852                        T i, T j, T k, T l, T m, T n, T o, T p)
01853 {
01854     x[0][0] = a;
01855     x[0][1] = b;
01856     x[0][2] = c;
01857     x[0][3] = d;
01858     x[1][0] = e;
01859     x[1][1] = f;
01860     x[1][2] = g;
01861     x[1][3] = h;
01862     x[2][0] = i;
01863     x[2][1] = j;
01864     x[2][2] = k;
01865     x[2][3] = l;
01866     x[3][0] = m;
01867     x[3][1] = n;
01868     x[3][2] = o;
01869     x[3][3] = p;
01870 }
01871 
01872 
01873 template <class T>
01874 inline
01875 Matrix44<T>::Matrix44 (Matrix33<T> r, Vec3<T> t)
01876 {
01877     x[0][0] = r[0][0];
01878     x[0][1] = r[0][1];
01879     x[0][2] = r[0][2];
01880     x[0][3] = 0;
01881     x[1][0] = r[1][0];
01882     x[1][1] = r[1][1];
01883     x[1][2] = r[1][2];
01884     x[1][3] = 0;
01885     x[2][0] = r[2][0];
01886     x[2][1] = r[2][1];
01887     x[2][2] = r[2][2];
01888     x[2][3] = 0;
01889     x[3][0] = t[0];
01890     x[3][1] = t[1];
01891     x[3][2] = t[2];
01892     x[3][3] = 1;
01893 }
01894 
01895 template <class T>
01896 inline
01897 Matrix44<T>::Matrix44 (const Matrix44 &v)
01898 {
01899     x[0][0] = v.x[0][0];
01900     x[0][1] = v.x[0][1];
01901     x[0][2] = v.x[0][2];
01902     x[0][3] = v.x[0][3];
01903     x[1][0] = v.x[1][0];
01904     x[1][1] = v.x[1][1];
01905     x[1][2] = v.x[1][2];
01906     x[1][3] = v.x[1][3];
01907     x[2][0] = v.x[2][0];
01908     x[2][1] = v.x[2][1];
01909     x[2][2] = v.x[2][2];
01910     x[2][3] = v.x[2][3];
01911     x[3][0] = v.x[3][0];
01912     x[3][1] = v.x[3][1];
01913     x[3][2] = v.x[3][2];
01914     x[3][3] = v.x[3][3];
01915 }
01916 
01917 template <class T>
01918 template <class S>
01919 inline
01920 Matrix44<T>::Matrix44 (const Matrix44<S> &v)
01921 {
01922     x[0][0] = T (v.x[0][0]);
01923     x[0][1] = T (v.x[0][1]);
01924     x[0][2] = T (v.x[0][2]);
01925     x[0][3] = T (v.x[0][3]);
01926     x[1][0] = T (v.x[1][0]);
01927     x[1][1] = T (v.x[1][1]);
01928     x[1][2] = T (v.x[1][2]);
01929     x[1][3] = T (v.x[1][3]);
01930     x[2][0] = T (v.x[2][0]);
01931     x[2][1] = T (v.x[2][1]);
01932     x[2][2] = T (v.x[2][2]);
01933     x[2][3] = T (v.x[2][3]);
01934     x[3][0] = T (v.x[3][0]);
01935     x[3][1] = T (v.x[3][1]);
01936     x[3][2] = T (v.x[3][2]);
01937     x[3][3] = T (v.x[3][3]);
01938 }
01939 
01940 template <class T>
01941 inline const Matrix44<T> &
01942 Matrix44<T>::operator = (const Matrix44 &v)
01943 {
01944     x[0][0] = v.x[0][0];
01945     x[0][1] = v.x[0][1];
01946     x[0][2] = v.x[0][2];
01947     x[0][3] = v.x[0][3];
01948     x[1][0] = v.x[1][0];
01949     x[1][1] = v.x[1][1];
01950     x[1][2] = v.x[1][2];
01951     x[1][3] = v.x[1][3];
01952     x[2][0] = v.x[2][0];
01953     x[2][1] = v.x[2][1];
01954     x[2][2] = v.x[2][2];
01955     x[2][3] = v.x[2][3];
01956     x[3][0] = v.x[3][0];
01957     x[3][1] = v.x[3][1];
01958     x[3][2] = v.x[3][2];
01959     x[3][3] = v.x[3][3];
01960     return *this;
01961 }
01962 
01963 template <class T>
01964 inline const Matrix44<T> &
01965 Matrix44<T>::operator = (T a)
01966 {
01967     x[0][0] = a;
01968     x[0][1] = a;
01969     x[0][2] = a;
01970     x[0][3] = a;
01971     x[1][0] = a;
01972     x[1][1] = a;
01973     x[1][2] = a;
01974     x[1][3] = a;
01975     x[2][0] = a;
01976     x[2][1] = a;
01977     x[2][2] = a;
01978     x[2][3] = a;
01979     x[3][0] = a;
01980     x[3][1] = a;
01981     x[3][2] = a;
01982     x[3][3] = a;
01983     return *this;
01984 }
01985 
01986 template <class T>
01987 inline T *
01988 Matrix44<T>::getValue ()
01989 {
01990     return (T *) &x[0][0];
01991 }
01992 
01993 template <class T>
01994 inline const T *
01995 Matrix44<T>::getValue () const
01996 {
01997     return (const T *) &x[0][0];
01998 }
01999 
02000 template <class T>
02001 template <class S>
02002 inline void
02003 Matrix44<T>::getValue (Matrix44<S> &v) const
02004 {
02005     if (isSameType<S,T>::value)
02006     {
02007         memcpy (v.x, x, sizeof (x));
02008     }
02009     else
02010     {
02011         v.x[0][0] = x[0][0];
02012         v.x[0][1] = x[0][1];
02013         v.x[0][2] = x[0][2];
02014         v.x[0][3] = x[0][3];
02015         v.x[1][0] = x[1][0];
02016         v.x[1][1] = x[1][1];
02017         v.x[1][2] = x[1][2];
02018         v.x[1][3] = x[1][3];
02019         v.x[2][0] = x[2][0];
02020         v.x[2][1] = x[2][1];
02021         v.x[2][2] = x[2][2];
02022         v.x[2][3] = x[2][3];
02023         v.x[3][0] = x[3][0];
02024         v.x[3][1] = x[3][1];
02025         v.x[3][2] = x[3][2];
02026         v.x[3][3] = x[3][3];
02027     }
02028 }
02029 
02030 template <class T>
02031 template <class S>
02032 inline Matrix44<T> &
02033 Matrix44<T>::setValue (const Matrix44<S> &v)
02034 {
02035     if (isSameType<S,T>::value)
02036     {
02037         memcpy (x, v.x, sizeof (x));
02038     }
02039     else
02040     {
02041         x[0][0] = v.x[0][0];
02042         x[0][1] = v.x[0][1];
02043         x[0][2] = v.x[0][2];
02044         x[0][3] = v.x[0][3];
02045         x[1][0] = v.x[1][0];
02046         x[1][1] = v.x[1][1];
02047         x[1][2] = v.x[1][2];
02048         x[1][3] = v.x[1][3];
02049         x[2][0] = v.x[2][0];
02050         x[2][1] = v.x[2][1];
02051         x[2][2] = v.x[2][2];
02052         x[2][3] = v.x[2][3];
02053         x[3][0] = v.x[3][0];
02054         x[3][1] = v.x[3][1];
02055         x[3][2] = v.x[3][2];
02056         x[3][3] = v.x[3][3];
02057     }
02058 
02059     return *this;
02060 }
02061 
02062 template <class T>
02063 template <class S>
02064 inline Matrix44<T> &
02065 Matrix44<T>::setTheMatrix (const Matrix44<S> &v)
02066 {
02067     if (isSameType<S,T>::value)
02068     {
02069         memcpy (x, v.x, sizeof (x));
02070     }
02071     else
02072     {
02073         x[0][0] = v.x[0][0];
02074         x[0][1] = v.x[0][1];
02075         x[0][2] = v.x[0][2];
02076         x[0][3] = v.x[0][3];
02077         x[1][0] = v.x[1][0];
02078         x[1][1] = v.x[1][1];
02079         x[1][2] = v.x[1][2];
02080         x[1][3] = v.x[1][3];
02081         x[2][0] = v.x[2][0];
02082         x[2][1] = v.x[2][1];
02083         x[2][2] = v.x[2][2];
02084         x[2][3] = v.x[2][3];
02085         x[3][0] = v.x[3][0];
02086         x[3][1] = v.x[3][1];
02087         x[3][2] = v.x[3][2];
02088         x[3][3] = v.x[3][3];
02089     }
02090 
02091     return *this;
02092 }
02093 
02094 template <class T>
02095 inline void
02096 Matrix44<T>::makeIdentity()
02097 {
02098     memset (x, 0, sizeof (x));
02099     x[0][0] = 1;
02100     x[1][1] = 1;
02101     x[2][2] = 1;
02102     x[3][3] = 1;
02103 }
02104 
02105 template <class T>
02106 bool
02107 Matrix44<T>::operator == (const Matrix44 &v) const
02108 {
02109     return x[0][0] == v.x[0][0] &&
02110            x[0][1] == v.x[0][1] &&
02111            x[0][2] == v.x[0][2] &&
02112            x[0][3] == v.x[0][3] &&
02113            x[1][0] == v.x[1][0] &&
02114            x[1][1] == v.x[1][1] &&
02115            x[1][2] == v.x[1][2] &&
02116            x[1][3] == v.x[1][3] &&
02117            x[2][0] == v.x[2][0] &&
02118            x[2][1] == v.x[2][1] &&
02119            x[2][2] == v.x[2][2] &&
02120            x[2][3] == v.x[2][3] &&
02121            x[3][0] == v.x[3][0] &&
02122            x[3][1] == v.x[3][1] &&
02123            x[3][2] == v.x[3][2] &&
02124            x[3][3] == v.x[3][3];
02125 }
02126 
02127 template <class T>
02128 bool
02129 Matrix44<T>::operator != (const Matrix44 &v) const
02130 {
02131     return x[0][0] != v.x[0][0] ||
02132            x[0][1] != v.x[0][1] ||
02133            x[0][2] != v.x[0][2] ||
02134            x[0][3] != v.x[0][3] ||
02135            x[1][0] != v.x[1][0] ||
02136            x[1][1] != v.x[1][1] ||
02137            x[1][2] != v.x[1][2] ||
02138            x[1][3] != v.x[1][3] ||
02139            x[2][0] != v.x[2][0] ||
02140            x[2][1] != v.x[2][1] ||
02141            x[2][2] != v.x[2][2] ||
02142            x[2][3] != v.x[2][3] ||
02143            x[3][0] != v.x[3][0] ||
02144            x[3][1] != v.x[3][1] ||
02145            x[3][2] != v.x[3][2] ||
02146            x[3][3] != v.x[3][3];
02147 }
02148 
02149 template <class T>
02150 bool
02151 Matrix44<T>::equalWithAbsError (const Matrix44<T> &m, T e) const
02152 {
02153     for (int i = 0; i < 4; i++)
02154         for (int j = 0; j < 4; j++)
02155             if (!Imath::equalWithAbsError ((*this)[i][j], m[i][j], e))
02156                 return false;
02157 
02158     return true;
02159 }
02160 
02161 template <class T>
02162 bool
02163 Matrix44<T>::equalWithRelError (const Matrix44<T> &m, T e) const
02164 {
02165     for (int i = 0; i < 4; i++)
02166         for (int j = 0; j < 4; j++)
02167             if (!Imath::equalWithRelError ((*this)[i][j], m[i][j], e))
02168                 return false;
02169 
02170     return true;
02171 }
02172 
02173 template <class T>
02174 const Matrix44<T> &
02175 Matrix44<T>::operator += (const Matrix44<T> &v)
02176 {
02177     x[0][0] += v.x[0][0];
02178     x[0][1] += v.x[0][1];
02179     x[0][2] += v.x[0][2];
02180     x[0][3] += v.x[0][3];
02181     x[1][0] += v.x[1][0];
02182     x[1][1] += v.x[1][1];
02183     x[1][2] += v.x[1][2];
02184     x[1][3] += v.x[1][3];
02185     x[2][0] += v.x[2][0];
02186     x[2][1] += v.x[2][1];
02187     x[2][2] += v.x[2][2];
02188     x[2][3] += v.x[2][3];
02189     x[3][0] += v.x[3][0];
02190     x[3][1] += v.x[3][1];
02191     x[3][2] += v.x[3][2];
02192     x[3][3] += v.x[3][3];
02193 
02194     return *this;
02195 }
02196 
02197 template <class T>
02198 const Matrix44<T> &
02199 Matrix44<T>::operator += (T a)
02200 {
02201     x[0][0] += a;
02202     x[0][1] += a;
02203     x[0][2] += a;
02204     x[0][3] += a;
02205     x[1][0] += a;
02206     x[1][1] += a;
02207     x[1][2] += a;
02208     x[1][3] += a;
02209     x[2][0] += a;
02210     x[2][1] += a;
02211     x[2][2] += a;
02212     x[2][3] += a;
02213     x[3][0] += a;
02214     x[3][1] += a;
02215     x[3][2] += a;
02216     x[3][3] += a;
02217 
02218     return *this;
02219 }
02220 
02221 template <class T>
02222 Matrix44<T>
02223 Matrix44<T>::operator + (const Matrix44<T> &v) const
02224 {
02225     return Matrix44 (x[0][0] + v.x[0][0],
02226                      x[0][1] + v.x[0][1],
02227                      x[0][2] + v.x[0][2],
02228                      x[0][3] + v.x[0][3],
02229                      x[1][0] + v.x[1][0],
02230                      x[1][1] + v.x[1][1],
02231                      x[1][2] + v.x[1][2],
02232                      x[1][3] + v.x[1][3],
02233                      x[2][0] + v.x[2][0],
02234                      x[2][1] + v.x[2][1],
02235                      x[2][2] + v.x[2][2],
02236                      x[2][3] + v.x[2][3],
02237                      x[3][0] + v.x[3][0],
02238                      x[3][1] + v.x[3][1],
02239                      x[3][2] + v.x[3][2],
02240                      x[3][3] + v.x[3][3]);
02241 }
02242 
02243 template <class T>
02244 const Matrix44<T> &
02245 Matrix44<T>::operator -= (const Matrix44<T> &v)
02246 {
02247     x[0][0] -= v.x[0][0];
02248     x[0][1] -= v.x[0][1];
02249     x[0][2] -= v.x[0][2];
02250     x[0][3] -= v.x[0][3];
02251     x[1][0] -= v.x[1][0];
02252     x[1][1] -= v.x[1][1];
02253     x[1][2] -= v.x[1][2];
02254     x[1][3] -= v.x[1][3];
02255     x[2][0] -= v.x[2][0];
02256     x[2][1] -= v.x[2][1];
02257     x[2][2] -= v.x[2][2];
02258     x[2][3] -= v.x[2][3];
02259     x[3][0] -= v.x[3][0];
02260     x[3][1] -= v.x[3][1];
02261     x[3][2] -= v.x[3][2];
02262     x[3][3] -= v.x[3][3];
02263 
02264     return *this;
02265 }
02266 
02267 template <class T>
02268 const Matrix44<T> &
02269 Matrix44<T>::operator -= (T a)
02270 {
02271     x[0][0] -= a;
02272     x[0][1] -= a;
02273     x[0][2] -= a;
02274     x[0][3] -= a;
02275     x[1][0] -= a;
02276     x[1][1] -= a;
02277     x[1][2] -= a;
02278     x[1][3] -= a;
02279     x[2][0] -= a;
02280     x[2][1] -= a;
02281     x[2][2] -= a;
02282     x[2][3] -= a;
02283     x[3][0] -= a;
02284     x[3][1] -= a;
02285     x[3][2] -= a;
02286     x[3][3] -= a;
02287 
02288     return *this;
02289 }
02290 
02291 template <class T>
02292 Matrix44<T>
02293 Matrix44<T>::operator - (const Matrix44<T> &v) const
02294 {
02295     return Matrix44 (x[0][0] - v.x[0][0],
02296                      x[0][1] - v.x[0][1],
02297                      x[0][2] - v.x[0][2],
02298                      x[0][3] - v.x[0][3],
02299                      x[1][0] - v.x[1][0],
02300                      x[1][1] - v.x[1][1],
02301                      x[1][2] - v.x[1][2],
02302                      x[1][3] - v.x[1][3],
02303                      x[2][0] - v.x[2][0],
02304                      x[2][1] - v.x[2][1],
02305                      x[2][2] - v.x[2][2],
02306                      x[2][3] - v.x[2][3],
02307                      x[3][0] - v.x[3][0],
02308                      x[3][1] - v.x[3][1],
02309                      x[3][2] - v.x[3][2],
02310                      x[3][3] - v.x[3][3]);
02311 }
02312 
02313 template <class T>
02314 Matrix44<T>
02315 Matrix44<T>::operator - () const
02316 {
02317     return Matrix44 (-x[0][0],
02318                      -x[0][1],
02319                      -x[0][2],
02320                      -x[0][3],
02321                      -x[1][0],
02322                      -x[1][1],
02323                      -x[1][2],
02324                      -x[1][3],
02325                      -x[2][0],
02326                      -x[2][1],
02327                      -x[2][2],
02328                      -x[2][3],
02329                      -x[3][0],
02330                      -x[3][1],
02331                      -x[3][2],
02332                      -x[3][3]);
02333 }
02334 
02335 template <class T>
02336 const Matrix44<T> &
02337 Matrix44<T>::negate ()
02338 {
02339     x[0][0] = -x[0][0];
02340     x[0][1] = -x[0][1];
02341     x[0][2] = -x[0][2];
02342     x[0][3] = -x[0][3];
02343     x[1][0] = -x[1][0];
02344     x[1][1] = -x[1][1];
02345     x[1][2] = -x[1][2];
02346     x[1][3] = -x[1][3];
02347     x[2][0] = -x[2][0];
02348     x[2][1] = -x[2][1];
02349     x[2][2] = -x[2][2];
02350     x[2][3] = -x[2][3];
02351     x[3][0] = -x[3][0];
02352     x[3][1] = -x[3][1];
02353     x[3][2] = -x[3][2];
02354     x[3][3] = -x[3][3];
02355 
02356     return *this;
02357 }
02358 
02359 template <class T>
02360 const Matrix44<T> &
02361 Matrix44<T>::operator *= (T a)
02362 {
02363     x[0][0] *= a;
02364     x[0][1] *= a;
02365     x[0][2] *= a;
02366     x[0][3] *= a;
02367     x[1][0] *= a;
02368     x[1][1] *= a;
02369     x[1][2] *= a;
02370     x[1][3] *= a;
02371     x[2][0] *= a;
02372     x[2][1] *= a;
02373     x[2][2] *= a;
02374     x[2][3] *= a;
02375     x[3][0] *= a;
02376     x[3][1] *= a;
02377     x[3][2] *= a;
02378     x[3][3] *= a;
02379 
02380     return *this;
02381 }
02382 
02383 template <class T>
02384 Matrix44<T>
02385 Matrix44<T>::operator * (T a) const
02386 {
02387     return Matrix44 (x[0][0] * a,
02388                      x[0][1] * a,
02389                      x[0][2] * a,
02390                      x[0][3] * a,
02391                      x[1][0] * a,
02392                      x[1][1] * a,
02393                      x[1][2] * a,
02394                      x[1][3] * a,
02395                      x[2][0] * a,
02396                      x[2][1] * a,
02397                      x[2][2] * a,
02398                      x[2][3] * a,
02399                      x[3][0] * a,
02400                      x[3][1] * a,
02401                      x[3][2] * a,
02402                      x[3][3] * a);
02403 }
02404 
02405 template <class T>
02406 inline Matrix44<T>
02407 operator * (T a, const Matrix44<T> &v)
02408 {
02409     return v * a;
02410 }
02411 
02412 template <class T>
02413 inline const Matrix44<T> &
02414 Matrix44<T>::operator *= (const Matrix44<T> &v)
02415 {
02416     Matrix44 tmp (T (0));
02417 
02418     multiply (*this, v, tmp);
02419     *this = tmp;
02420     return *this;
02421 }
02422 
02423 template <class T>
02424 inline Matrix44<T>
02425 Matrix44<T>::operator * (const Matrix44<T> &v) const
02426 {
02427     Matrix44 tmp (T (0));
02428 
02429     multiply (*this, v, tmp);
02430     return tmp;
02431 }
02432 
02433 template <class T>
02434 void
02435 Matrix44<T>::multiply (const Matrix44<T> &a,
02436                        const Matrix44<T> &b,
02437                        Matrix44<T> &c)
02438 {
02439         register const T * IMATH_RESTRICT ap = &a.x[0][0];
02440         register const T * IMATH_RESTRICT bp = &b.x[0][0];
02441         register       T * IMATH_RESTRICT cp = &c.x[0][0];
02442 
02443     register T a0, a1, a2, a3;
02444 
02445     a0 = ap[0];
02446     a1 = ap[1];
02447     a2 = ap[2];
02448     a3 = ap[3];
02449 
02450     cp[0]  = a0 * bp[0]  + a1 * bp[4]  + a2 * bp[8]  + a3 * bp[12];
02451     cp[1]  = a0 * bp[1]  + a1 * bp[5]  + a2 * bp[9]  + a3 * bp[13];
02452     cp[2]  = a0 * bp[2]  + a1 * bp[6]  + a2 * bp[10] + a3 * bp[14];
02453     cp[3]  = a0 * bp[3]  + a1 * bp[7]  + a2 * bp[11] + a3 * bp[15];
02454 
02455     a0 = ap[4];
02456     a1 = ap[5];
02457     a2 = ap[6];
02458     a3 = ap[7];
02459 
02460     cp[4]  = a0 * bp[0]  + a1 * bp[4]  + a2 * bp[8]  + a3 * bp[12];
02461     cp[5]  = a0 * bp[1]  + a1 * bp[5]  + a2 * bp[9]  + a3 * bp[13];
02462     cp[6]  = a0 * bp[2]  + a1 * bp[6]  + a2 * bp[10] + a3 * bp[14];
02463     cp[7]  = a0 * bp[3]  + a1 * bp[7]  + a2 * bp[11] + a3 * bp[15];
02464 
02465     a0 = ap[8];
02466     a1 = ap[9];
02467     a2 = ap[10];
02468     a3 = ap[11];
02469 
02470     cp[8]  = a0 * bp[0]  + a1 * bp[4]  + a2 * bp[8]  + a3 * bp[12];
02471     cp[9]  = a0 * bp[1]  + a1 * bp[5]  + a2 * bp[9]  + a3 * bp[13];
02472     cp[10] = a0 * bp[2]  + a1 * bp[6]  + a2 * bp[10] + a3 * bp[14];
02473     cp[11] = a0 * bp[3]  + a1 * bp[7]  + a2 * bp[11] + a3 * bp[15];
02474 
02475     a0 = ap[12];
02476     a1 = ap[13];
02477     a2 = ap[14];
02478     a3 = ap[15];
02479 
02480     cp[12] = a0 * bp[0]  + a1 * bp[4]  + a2 * bp[8]  + a3 * bp[12];
02481     cp[13] = a0 * bp[1]  + a1 * bp[5]  + a2 * bp[9]  + a3 * bp[13];
02482     cp[14] = a0 * bp[2]  + a1 * bp[6]  + a2 * bp[10] + a3 * bp[14];
02483     cp[15] = a0 * bp[3]  + a1 * bp[7]  + a2 * bp[11] + a3 * bp[15];
02484 }
02485 
02486 template <class T> template <class S>
02487 void
02488 Matrix44<T>::multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const
02489 {
02490     S a, b, c, w;
02491 
02492     a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0] + x[3][0];
02493     b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1] + x[3][1];
02494     c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2] + x[3][2];
02495     w = src[0] * x[0][3] + src[1] * x[1][3] + src[2] * x[2][3] + x[3][3];
02496 
02497     dst.x = a / w;
02498     dst.y = b / w;
02499     dst.z = c / w;
02500 }
02501 
02502 template <class T> template <class S>
02503 void
02504 Matrix44<T>::multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const
02505 {
02506     S a, b, c;
02507 
02508     a = src[0] * x[0][0] + src[1] * x[1][0] + src[2] * x[2][0];
02509     b = src[0] * x[0][1] + src[1] * x[1][1] + src[2] * x[2][1];
02510     c = src[0] * x[0][2] + src[1] * x[1][2] + src[2] * x[2][2];
02511 
02512     dst.x = a;
02513     dst.y = b;
02514     dst.z = c;
02515 }
02516 
02517 template <class T>
02518 const Matrix44<T> &
02519 Matrix44<T>::operator /= (T a)
02520 {
02521     x[0][0] /= a;
02522     x[0][1] /= a;
02523     x[0][2] /= a;
02524     x[0][3] /= a;
02525     x[1][0] /= a;
02526     x[1][1] /= a;
02527     x[1][2] /= a;
02528     x[1][3] /= a;
02529     x[2][0] /= a;
02530     x[2][1] /= a;
02531     x[2][2] /= a;
02532     x[2][3] /= a;
02533     x[3][0] /= a;
02534     x[3][1] /= a;
02535     x[3][2] /= a;
02536     x[3][3] /= a;
02537 
02538     return *this;
02539 }
02540 
02541 template <class T>
02542 Matrix44<T>
02543 Matrix44<T>::operator / (T a) const
02544 {
02545     return Matrix44 (x[0][0] / a,
02546                      x[0][1] / a,
02547                      x[0][2] / a,
02548                      x[0][3] / a,
02549                      x[1][0] / a,
02550                      x[1][1] / a,
02551                      x[1][2] / a,
02552                      x[1][3] / a,
02553                      x[2][0] / a,
02554                      x[2][1] / a,
02555                      x[2][2] / a,
02556                      x[2][3] / a,
02557                      x[3][0] / a,
02558                      x[3][1] / a,
02559                      x[3][2] / a,
02560                      x[3][3] / a);
02561 }
02562 
02563 template <class T>
02564 const Matrix44<T> &
02565 Matrix44<T>::transpose ()
02566 {
02567     Matrix44 tmp (x[0][0],
02568                   x[1][0],
02569                   x[2][0],
02570                   x[3][0],
02571                   x[0][1],
02572                   x[1][1],
02573                   x[2][1],
02574                   x[3][1],
02575                   x[0][2],
02576                   x[1][2],
02577                   x[2][2],
02578                   x[3][2],
02579                   x[0][3],
02580                   x[1][3],
02581                   x[2][3],
02582                   x[3][3]);
02583     *this = tmp;
02584     return *this;
02585 }
02586 
02587 template <class T>
02588 Matrix44<T>
02589 Matrix44<T>::transposed () const
02590 {
02591     return Matrix44 (x[0][0],
02592                      x[1][0],
02593                      x[2][0],
02594                      x[3][0],
02595                      x[0][1],
02596                      x[1][1],
02597                      x[2][1],
02598                      x[3][1],
02599                      x[0][2],
02600                      x[1][2],
02601                      x[2][2],
02602                      x[3][2],
02603                      x[0][3],
02604                      x[1][3],
02605                      x[2][3],
02606                      x[3][3]);
02607 }
02608 
02609 template <class T>
02610 const Matrix44<T> &
02611 Matrix44<T>::gjInvert (bool singExc) throw (Iex::MathExc)
02612 {
02613     *this = gjInverse (singExc);
02614     return *this;
02615 }
02616 
02617 template <class T>
02618 Matrix44<T>
02619 Matrix44<T>::gjInverse (bool singExc) const throw (Iex::MathExc)
02620 {
02621     int i, j, k;
02622     Matrix44 s;
02623     Matrix44 t (*this);
02624 
02625     // Forward elimination
02626 
02627     for (i = 0; i < 3 ; i++)
02628     {
02629         int pivot = i;
02630 
02631         T pivotsize = t[i][i];
02632 
02633         if (pivotsize < 0)
02634             pivotsize = -pivotsize;
02635 
02636         for (j = i + 1; j < 4; j++)
02637         {
02638             T tmp = t[j][i];
02639 
02640             if (tmp < 0)
02641                 tmp = -tmp;
02642 
02643             if (tmp > pivotsize)
02644             {
02645                 pivot = j;
02646                 pivotsize = tmp;
02647             }
02648         }
02649 
02650         if (pivotsize == 0)
02651         {
02652             if (singExc)
02653                 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
02654 
02655             return Matrix44();
02656         }
02657 
02658         if (pivot != i)
02659         {
02660             for (j = 0; j < 4; j++)
02661             {
02662                 T tmp;
02663 
02664                 tmp = t[i][j];
02665                 t[i][j] = t[pivot][j];
02666                 t[pivot][j] = tmp;
02667 
02668                 tmp = s[i][j];
02669                 s[i][j] = s[pivot][j];
02670                 s[pivot][j] = tmp;
02671             }
02672         }
02673 
02674         for (j = i + 1; j < 4; j++)
02675         {
02676             T f = t[j][i] / t[i][i];
02677 
02678             for (k = 0; k < 4; k++)
02679             {
02680                 t[j][k] -= f * t[i][k];
02681                 s[j][k] -= f * s[i][k];
02682             }
02683         }
02684     }
02685 
02686     // Backward substitution
02687 
02688     for (i = 3; i >= 0; --i)
02689     {
02690         T f;
02691 
02692         if ((f = t[i][i]) == 0)
02693         {
02694             if (singExc)
02695                 throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");
02696 
02697             return Matrix44();
02698         }
02699 
02700         for (j = 0; j < 4; j++)
02701         {
02702             t[i][j] /= f;
02703             s[i][j] /= f;
02704         }
02705 
02706         for (j = 0; j < i; j++)
02707         {
02708             f = t[j][i];
02709 
02710             for (k = 0; k < 4; k++)
02711             {
02712                 t[j][k] -= f * t[i][k];
02713                 s[j][k] -= f * s[i][k];
02714             }
02715         }
02716     }
02717 
02718     return s;
02719 }
02720 
02721 template <class T>
02722 const Matrix44<T> &
02723 Matrix44<T>::invert (bool singExc) throw (Iex::MathExc)
02724 {
02725     *this = inverse (singExc);
02726     return *this;
02727 }
02728 
02729 template <class T>
02730 Matrix44<T>
02731 Matrix44<T>::inverse (bool singExc) const throw (Iex::MathExc)
02732 {
02733     if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)
02734         return gjInverse(singExc);
02735 
02736     Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
02737                 x[2][1] * x[0][2] - x[0][1] * x[2][2],
02738                 x[0][1] * x[1][2] - x[1][1] * x[0][2],
02739                 0,
02740 
02741                 x[2][0] * x[1][2] - x[1][0] * x[2][2],
02742                 x[0][0] * x[2][2] - x[2][0] * x[0][2],
02743                 x[1][0] * x[0][2] - x[0][0] * x[1][2],
02744                 0,
02745 
02746                 x[1][0] * x[2][1] - x[2][0] * x[1][1],
02747                 x[2][0] * x[0][1] - x[0][0] * x[2][1],
02748                 x[0][0] * x[1][1] - x[1][0] * x[0][1],
02749                 0,
02750 
02751                 0,
02752                 0,
02753                 0,
02754                 1);
02755 
02756     T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
02757 
02758     if (Imath::abs (r) >= 1)
02759     {
02760         for (int i = 0; i < 3; ++i)
02761         {
02762             for (int j = 0; j < 3; ++j)
02763             {
02764                 s[i][j] /= r;
02765             }
02766         }
02767     }
02768     else
02769     {
02770         T mr = Imath::abs (r) / limits<T>::smallest();
02771 
02772         for (int i = 0; i < 3; ++i)
02773         {
02774             for (int j = 0; j < 3; ++j)
02775             {
02776                 if (mr > Imath::abs (s[i][j]))
02777                 {
02778                     s[i][j] /= r;
02779                 }
02780                 else
02781                 {
02782                     if (singExc)
02783                         throw SingMatrixExc ("Cannot invert singular matrix.");
02784 
02785                     return Matrix44();
02786                 }
02787             }
02788         }
02789     }
02790 
02791     s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0];
02792     s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1];
02793     s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2];
02794 
02795     return s;
02796 }
02797 
02798 template <class T>
02799 template <class S>
02800 const Matrix44<T> &
02801 Matrix44<T>::setEulerAngles (const Vec3<S>& r)
02802 {
02803     S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
02804     
02805     cos_rz = Math<T>::cos (r[2]);
02806     cos_ry = Math<T>::cos (r[1]);
02807     cos_rx = Math<T>::cos (r[0]);
02808     
02809     sin_rz = Math<T>::sin (r[2]);
02810     sin_ry = Math<T>::sin (r[1]);
02811     sin_rx = Math<T>::sin (r[0]);
02812     
02813     x[0][0] =  cos_rz * cos_ry;
02814     x[0][1] =  sin_rz * cos_ry;
02815     x[0][2] = -sin_ry;
02816     x[0][3] =  0;
02817     
02818     x[1][0] = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
02819     x[1][1] =  cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
02820     x[1][2] =  cos_ry * sin_rx;
02821     x[1][3] =  0;
02822     
02823     x[2][0] =  sin_rz * sin_rx + cos_rz * sin_ry * cos_rx;
02824     x[2][1] = -cos_rz * sin_rx + sin_rz * sin_ry * cos_rx;
02825     x[2][2] =  cos_ry * cos_rx;
02826     x[2][3] =  0;
02827 
02828     x[3][0] =  0;
02829     x[3][1] =  0;
02830     x[3][2] =  0;
02831     x[3][3] =  1;
02832 
02833     return *this;
02834 }
02835 
02836 template <class T>
02837 template <class S>
02838 const Matrix44<T> &
02839 Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle)
02840 {
02841     Vec3<S> unit (axis.normalized());
02842     S sine   = Math<T>::sin (angle);
02843     S cosine = Math<T>::cos (angle);
02844 
02845     x[0][0] = unit[0] * unit[0] * (1 - cosine) + cosine;
02846     x[0][1] = unit[0] * unit[1] * (1 - cosine) + unit[2] * sine;
02847     x[0][2] = unit[0] * unit[2] * (1 - cosine) - unit[1] * sine;
02848     x[0][3] = 0;
02849 
02850     x[1][0] = unit[0] * unit[1] * (1 - cosine) - unit[2] * sine;
02851     x[1][1] = unit[1] * unit[1] * (1 - cosine) + cosine;
02852     x[1][2] = unit[1] * unit[2] * (1 - cosine) + unit[0] * sine;
02853     x[1][3] = 0;
02854 
02855     x[2][0] = unit[0] * unit[2] * (1 - cosine) + unit[1] * sine;
02856     x[2][1] = unit[1] * unit[2] * (1 - cosine) - unit[0] * sine;
02857     x[2][2] = unit[2] * unit[2] * (1 - cosine) + cosine;
02858     x[2][3] = 0;
02859 
02860     x[3][0] = 0;
02861     x[3][1] = 0;
02862     x[3][2] = 0;
02863     x[3][3] = 1;
02864 
02865     return *this;
02866 }
02867 
02868 template <class T>
02869 template <class S>
02870 const Matrix44<T> &
02871 Matrix44<T>::rotate (const Vec3<S> &r)
02872 {
02873     S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
02874     S m00, m01, m02;
02875     S m10, m11, m12;
02876     S m20, m21, m22;
02877 
02878     cos_rz = Math<S>::cos (r[2]);
02879     cos_ry = Math<S>::cos (r[1]);
02880     cos_rx = Math<S>::cos (r[0]);
02881     
02882     sin_rz = Math<S>::sin (r[2]);
02883     sin_ry = Math<S>::sin (r[1]);
02884     sin_rx = Math<S>::sin (r[0]);
02885 
02886     m00 =  cos_rz *  cos_ry;
02887     m01 =  sin_rz *  cos_ry;
02888     m02 = -sin_ry;
02889     m10 = -sin_rz *  cos_rx + cos_rz * sin_ry * sin_rx;
02890     m11 =  cos_rz *  cos_rx + sin_rz * sin_ry * sin_rx;
02891     m12 =  cos_ry *  sin_rx;
02892     m20 = -sin_rz * -sin_rx + cos_rz * sin_ry * cos_rx;
02893     m21 =  cos_rz * -sin_rx + sin_rz * sin_ry * cos_rx;
02894     m22 =  cos_ry *  cos_rx;
02895 
02896     Matrix44<T> P (*this);
02897 
02898     x[0][0] = P[0][0] * m00 + P[1][0] * m01 + P[2][0] * m02;
02899     x[0][1] = P[0][1] * m00 + P[1][1] * m01 + P[2][1] * m02;
02900     x[0][2] = P[0][2] * m00 + P[1][2] * m01 + P[2][2] * m02;
02901     x[0][3] = P[0][3] * m00 + P[1][3] * m01 + P[2][3] * m02;
02902 
02903     x[1][0] = P[0][0] * m10 + P[1][0] * m11 + P[2][0] * m12;
02904     x[1][1] = P[0][1] * m10 + P[1][1] * m11 + P[2][1] * m12;
02905     x[1][2] = P[0][2] * m10 + P[1][2] * m11 + P[2][2] * m12;
02906     x[1][3] = P[0][3] * m10 + P[1][3] * m11 + P[2][3] * m12;
02907 
02908     x[2][0] = P[0][0] * m20 + P[1][0] * m21 + P[2][0] * m22;
02909     x[2][1] = P[0][1] * m20 + P[1][1] * m21 + P[2][1] * m22;
02910     x[2][2] = P[0][2] * m20 + P[1][2] * m21 + P[2][2] * m22;
02911     x[2][3] = P[0][3] * m20 + P[1][3] * m21 + P[2][3] * m22;
02912 
02913     return *this;
02914 }
02915 
02916 template <class T>
02917 const Matrix44<T> &
02918 Matrix44<T>::setScale (T s)
02919 {
02920     memset (x, 0, sizeof (x));
02921     x[0][0] = s;
02922     x[1][1] = s;
02923     x[2][2] = s;
02924     x[3][3] = 1;
02925 
02926     return *this;
02927 }
02928 
02929 template <class T>
02930 template <class S>
02931 const Matrix44<T> &
02932 Matrix44<T>::setScale (const Vec3<S> &s)
02933 {
02934     memset (x, 0, sizeof (x));
02935     x[0][0] = s[0];
02936     x[1][1] = s[1];
02937     x[2][2] = s[2];
02938     x[3][3] = 1;
02939 
02940     return *this;
02941 }
02942 
02943 template <class T>
02944 template <class S>
02945 const Matrix44<T> &
02946 Matrix44<T>::scale (const Vec3<S> &s)
02947 {
02948     x[0][0] *= s[0];
02949     x[0][1] *= s[0];
02950     x[0][2] *= s[0];
02951     x[0][3] *= s[0];
02952 
02953     x[1][0] *= s[1];
02954     x[1][1] *= s[1];
02955     x[1][2] *= s[1];
02956     x[1][3] *= s[1];
02957 
02958     x[2][0] *= s[2];
02959     x[2][1] *= s[2];
02960     x[2][2] *= s[2];
02961     x[2][3] *= s[2];
02962 
02963     return *this;
02964 }
02965 
02966 template <class T>
02967 template <class S>
02968 const Matrix44<T> &
02969 Matrix44<T>::setTranslation (const Vec3<S> &t)
02970 {
02971     x[0][0] = 1;
02972     x[0][1] = 0;
02973     x[0][2] = 0;
02974     x[0][3] = 0;
02975 
02976     x[1][0] = 0;
02977     x[1][1] = 1;
02978     x[1][2] = 0;
02979     x[1][3] = 0;
02980 
02981     x[2][0] = 0;
02982     x[2][1] = 0;
02983     x[2][2] = 1;
02984     x[2][3] = 0;
02985 
02986     x[3][0] = t[0];
02987     x[3][1] = t[1];
02988     x[3][2] = t[2];
02989     x[3][3] = 1;
02990 
02991     return *this;
02992 }
02993 
02994 template <class T>
02995 inline const Vec3<T>
02996 Matrix44<T>::translation () const
02997 {
02998     return Vec3<T> (x[3][0], x[3][1], x[3][2]);
02999 }
03000 
03001 template <class T>
03002 template <class S>
03003 const Matrix44<T> &
03004 Matrix44<T>::translate (const Vec3<S> &t)
03005 {
03006     x[3][0] += t[0] * x[0][0] + t[1] * x[1][0] + t[2] * x[2][0];
03007     x[3][1] += t[0] * x[0][1] + t[1] * x[1][1] + t[2] * x[2][1];
03008     x[3][2] += t[0] * x[0][2] + t[1] * x[1][2] + t[2] * x[2][2];
03009     x[3][3] += t[0] * x[0][3] + t[1] * x[1][3] + t[2] * x[2][3];
03010 
03011     return *this;
03012 }
03013 
03014 template <class T>
03015 template <class S>
03016 const Matrix44<T> &
03017 Matrix44<T>::setShear (const Vec3<S> &h)
03018 {
03019     x[0][0] = 1;
03020     x[0][1] = 0;
03021     x[0][2] = 0;
03022     x[0][3] = 0;
03023 
03024     x[1][0] = h[0];
03025     x[1][1] = 1;
03026     x[1][2] = 0;
03027     x[1][3] = 0;
03028 
03029     x[2][0] = h[1];
03030     x[2][1] = h[2];
03031     x[2][2] = 1;
03032     x[2][3] = 0;
03033 
03034     x[3][0] = 0;
03035     x[3][1] = 0;
03036     x[3][2] = 0;
03037     x[3][3] = 1;
03038 
03039     return *this;
03040 }
03041 
03042 template <class T>
03043 template <class S>
03044 const Matrix44<T> &
03045 Matrix44<T>::setShear (const Shear6<S> &h)
03046 {
03047     x[0][0] = 1;
03048     x[0][1] = h.yx;
03049     x[0][2] = h.zx;
03050     x[0][3] = 0;
03051 
03052     x[1][0] = h.xy;
03053     x[1][1] = 1;
03054     x[1][2] = h.zy;
03055     x[1][3] = 0;
03056 
03057     x[2][0] = h.xz;
03058     x[2][1] = h.yz;
03059     x[2][2] = 1;
03060     x[2][3] = 0;
03061 
03062     x[3][0] = 0;
03063     x[3][1] = 0;
03064     x[3][2] = 0;
03065     x[3][3] = 1;
03066 
03067     return *this;
03068 }
03069 
03070 template <class T>
03071 template <class S>
03072 const Matrix44<T> &
03073 Matrix44<T>::shear (const Vec3<S> &h)
03074 {
03075     //
03076     // In this case, we don't need a temp. copy of the matrix 
03077     // because we never use a value on the RHS after we've 
03078     // changed it on the LHS.
03079     // 
03080 
03081     for (int i=0; i < 4; i++)
03082     {
03083         x[2][i] += h[1] * x[0][i] + h[2] * x[1][i];
03084         x[1][i] += h[0] * x[0][i];
03085     }
03086 
03087     return *this;
03088 }
03089 
03090 template <class T>
03091 template <class S>
03092 const Matrix44<T> &
03093 Matrix44<T>::shear (const Shear6<S> &h)
03094 {
03095     Matrix44<T> P (*this);
03096 
03097     for (int i=0; i < 4; i++)
03098     {
03099         x[0][i] = P[0][i] + h.yx * P[1][i] + h.zx * P[2][i];
03100         x[1][i] = h.xy * P[0][i] + P[1][i] + h.zy * P[2][i];
03101         x[2][i] = h.xz * P[0][i] + h.yz * P[1][i] + P[2][i];
03102     }
03103 
03104     return *this;
03105 }
03106 
03107 
03108 //--------------------------------
03109 // Implementation of stream output
03110 //--------------------------------
03111 
03112 template <class T>
03113 std::ostream &
03114 operator << (std::ostream &s, const Matrix33<T> &m)
03115 {
03116     std::ios_base::fmtflags oldFlags = s.flags();
03117     int width;
03118 
03119     if (s.flags() & std::ios_base::fixed)
03120     {
03121         s.setf (std::ios_base::showpoint);
03122         width = s.precision() + 5;
03123     }
03124     else
03125     {
03126         s.setf (std::ios_base::scientific);
03127         s.setf (std::ios_base::showpoint);
03128         width = s.precision() + 8;
03129     }
03130 
03131     s << "(" << std::setw (width) << m[0][0] <<
03132          " " << std::setw (width) << m[0][1] <<
03133          " " << std::setw (width) << m[0][2] << "\n" <<
03134 
03135          " " << std::setw (width) << m[1][0] <<
03136          " " << std::setw (width) << m[1][1] <<
03137          " " << std::setw (width) << m[1][2] << "\n" <<
03138 
03139          " " << std::setw (width) << m[2][0] <<
03140          " " << std::setw (width) << m[2][1] <<
03141          " " << std::setw (width) << m[2][2] << ")\n";
03142 
03143     s.flags (oldFlags);
03144     return s;
03145 }
03146 
03147 template <class T>
03148 std::ostream &
03149 operator << (std::ostream &s, const Matrix44<T> &m)
03150 {
03151     std::ios_base::fmtflags oldFlags = s.flags();
03152     int width;
03153 
03154     if (s.flags() & std::ios_base::fixed)
03155     {
03156         s.setf (std::ios_base::showpoint);
03157         width = s.precision() + 5;
03158     }
03159     else
03160     {
03161         s.setf (std::ios_base::scientific);
03162         s.setf (std::ios_base::showpoint);
03163         width = s.precision() + 8;
03164     }
03165 
03166     s << "(" << std::setw (width) << m[0][0] <<
03167          " " << std::setw (width) << m[0][1] <<
03168          " " << std::setw (width) << m[0][2] <<
03169          " " << std::setw (width) << m[0][3] << "\n" <<
03170 
03171          " " << std::setw (width) << m[1][0] <<
03172          " " << std::setw (width) << m[1][1] <<
03173          " " << std::setw (width) << m[1][2] <<
03174          " " << std::setw (width) << m[1][3] << "\n" <<
03175 
03176          " " << std::setw (width) << m[2][0] <<
03177          " " << std::setw (width) << m[2][1] <<
03178          " " << std::setw (width) << m[2][2] <<
03179          " " << std::setw (width) << m[2][3] << "\n" <<
03180 
03181          " " << std::setw (width) << m[3][0] <<
03182          " " << std::setw (width) << m[3][1] <<
03183          " " << std::setw (width) << m[3][2] <<
03184          " " << std::setw (width) << m[3][3] << ")\n";
03185 
03186     s.flags (oldFlags);
03187     return s;
03188 }
03189 
03190 
03191 //---------------------------------------------------------------
03192 // Implementation of vector-times-matrix multiplication operators
03193 //---------------------------------------------------------------
03194 
03195 template <class S, class T>
03196 inline const Vec2<S> &
03197 operator *= (Vec2<S> &v, const Matrix33<T> &m)
03198 {
03199     S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
03200     S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
03201     S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
03202 
03203     v.x = x / w;
03204     v.y = y / w;
03205 
03206     return v;
03207 }
03208 
03209 template <class S, class T>
03210 inline Vec2<S>
03211 operator * (const Vec2<S> &v, const Matrix33<T> &m)
03212 {
03213     S x = S(v.x * m[0][0] + v.y * m[1][0] + m[2][0]);
03214     S y = S(v.x * m[0][1] + v.y * m[1][1] + m[2][1]);
03215     S w = S(v.x * m[0][2] + v.y * m[1][2] + m[2][2]);
03216 
03217     return Vec2<S> (x / w, y / w);
03218 }
03219 
03220 
03221 template <class S, class T>
03222 inline const Vec3<S> &
03223 operator *= (Vec3<S> &v, const Matrix33<T> &m)
03224 {
03225     S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
03226     S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
03227     S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
03228 
03229     v.x = x;
03230     v.y = y;
03231     v.z = z;
03232 
03233     return v;
03234 }
03235 
03236 template <class S, class T>
03237 inline Vec3<S>
03238 operator * (const Vec3<S> &v, const Matrix33<T> &m)
03239 {
03240     S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0]);
03241     S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1]);
03242     S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
03243 
03244     return Vec3<S> (x, y, z);
03245 }
03246 
03247 
03248 template <class S, class T>
03249 inline const Vec3<S> &
03250 operator *= (Vec3<S> &v, const Matrix44<T> &m)
03251 {
03252     S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
03253     S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
03254     S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
03255     S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
03256 
03257     v.x = x / w;
03258     v.y = y / w;
03259     v.z = z / w;
03260 
03261     return v;
03262 }
03263 
03264 template <class S, class T>
03265 inline Vec3<S>
03266 operator * (const Vec3<S> &v, const Matrix44<T> &m)
03267 {
03268     S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]);
03269     S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]);
03270     S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]);
03271     S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + m[3][3]);
03272 
03273     return Vec3<S> (x / w, y / w, z / w);
03274 }
03275 
03276 
03277 template <class S, class T>
03278 inline const Vec4<S> &
03279 operator *= (Vec4<S> &v, const Matrix44<T> &m)
03280 {
03281     S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + v.w * m[3][0]);
03282     S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + v.w * m[3][1]);
03283     S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + v.w * m[3][2]);
03284     S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + v.w * m[3][3]);
03285 
03286     v.x = x;
03287     v.y = y;
03288     v.z = z;
03289     v.w = w;
03290 
03291     return v;
03292 }
03293 
03294 template <class S, class T>
03295 inline Vec4<S>
03296 operator * (const Vec4<S> &v, const Matrix44<T> &m)
03297 {
03298     S x = S(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + v.w * m[3][0]);
03299     S y = S(v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + v.w * m[3][1]);
03300     S z = S(v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + v.w * m[3][2]);
03301     S w = S(v.x * m[0][3] + v.y * m[1][3] + v.z * m[2][3] + v.w * m[3][3]);
03302 
03303     return Vec4<S> (x, y, z, w);
03304 }
03305 
03306 } // namespace Imath
03307 
03308 
03309 
03310 #endif