UT/UT_DMatrix4.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * Produced by:
 *      Cristin Barghiel
 *      Side Effects Software Inc.
 *      123 Front Street West, Suite 1401
 *      Toronto, Ontario
 *      Canada   M5J 2M2
 *      416-366-4607
 */

#ifndef __UT_DMatrix4_h__
#define __UT_DMatrix4_h__

#include "UT_API.h"
#include <iostream.h>
#include "UT_Assert.h"
#include "UT_Axis.h"
#include "UT_Math.h"

class UT_IStream;
class UT_Matrix3;
class UT_Matrix4;
class UT_DMatrix3;
class UT_DMatrix4;
class UT_XformOrder;
class UT_Quaternion;
class UT_Vector3;
class UT_Vector4;

// Free floating operators that return a UT_DMatrix4 object. 
UT_API UT_DMatrix4      operator+(const UT_DMatrix4 &m1,const UT_DMatrix4 &m2);
UT_API UT_DMatrix4      operator-(const UT_DMatrix4 &m1,const UT_DMatrix4 &m2);
UT_API UT_DMatrix4      operator*(const UT_DMatrix4 &m1,const UT_DMatrix4 &m2);
UT_API UT_DMatrix4      operator+(const UT_DMatrix4 &m, const UT_Vector4 &v);
UT_API UT_DMatrix4      operator+(const UT_Vector4 &v, const UT_DMatrix4 &m);
UT_API UT_DMatrix4      operator-(const UT_DMatrix4 &m, const UT_Vector4 &v);
UT_API UT_DMatrix4      operator-(const UT_Vector4 &v, const UT_DMatrix4 &m);
UT_API UT_DMatrix4      operator+(const UT_DMatrix4 &mat, fpreal64 sc);
inline UT_DMatrix4      operator-(const UT_DMatrix4 &mat, fpreal64 sc);
UT_API UT_DMatrix4      operator*(const UT_DMatrix4 &mat, fpreal64 sc);
inline UT_DMatrix4      operator/(const UT_DMatrix4 &mat, fpreal64 sc);
inline UT_DMatrix4      operator+(fpreal64 sc, const UT_DMatrix4 &mat);
UT_API UT_DMatrix4      operator-(fpreal64 sc, const UT_DMatrix4 &mat);
inline UT_DMatrix4      operator*(fpreal64 sc, const UT_DMatrix4 &mat);
UT_API UT_DMatrix4      operator/(fpreal64 sc, const UT_DMatrix4 &mat);

/// This class implements a 4x4 double matrix in row-major order.
///
/// Most of Houdini operates with row vectors that are left-multiplied with
/// matrices. e.g.,   z = v * M
///   As a result, translation data is in row 3 of the matrix, rather than
/// column 3.
///
/// This convention, combined with row-major order, is directly compatible
/// with OpenGL matrix requirements.
class UT_API UT_DMatrix4
{
public:
    /// Construct uninitialized matrix.
    UT_DMatrix4()
    { }
    /// Construct identity matrix, multipled by scalar.
    explicit UT_DMatrix4(fpreal64 val)
    {
        matx[0][0] = val; matx[0][1] = 0; matx[0][2] = 0; matx[0][3] = 0;
        matx[1][0] = 0; matx[1][1] = val; matx[1][2] = 0; matx[1][3] = 0;
        matx[2][0] = 0; matx[2][1] = 0; matx[2][2] = val; matx[2][3] = 0;
        matx[3][0] = 0; matx[3][1] = 0; matx[3][2] = 0; matx[3][3] = val;
    }
    /// Construct a deep copy of the input row-major data.
    // @{
    explicit UT_DMatrix4(const fpreal32 m[4][4])
    {
        matx[0][0]=m[0][0]; matx[0][1]=m[0][1]; 
        matx[0][2]=m[0][2]; matx[0][3]=m[0][3];

        matx[1][0]=m[1][0]; matx[1][1]=m[1][1]; 
        matx[1][2]=m[1][2]; matx[1][3]=m[1][3];

        matx[2][0]=m[2][0]; matx[2][1]=m[2][1]; 
        matx[2][2]=m[2][2]; matx[2][3]=m[2][3];

        matx[3][0]=m[3][0]; matx[3][1]=m[3][1]; 
        matx[3][2]=m[3][2]; matx[3][3]=m[3][3];
    }
    explicit UT_DMatrix4(const fpreal64 m[4][4])
    {
        matx[0][0]=m[0][0]; matx[0][1]=m[0][1]; 
        matx[0][2]=m[0][2]; matx[0][3]=m[0][3];

        matx[1][0]=m[1][0]; matx[1][1]=m[1][1]; 
        matx[1][2]=m[1][2]; matx[1][3]=m[1][3];

        matx[2][0]=m[2][0]; matx[2][1]=m[2][1]; 
        matx[2][2]=m[2][2]; matx[2][3]=m[2][3];

        matx[3][0]=m[3][0]; matx[3][1]=m[3][1]; 
        matx[3][2]=m[3][2]; matx[3][3]=m[3][3];
    }
    // @}

    /// This constructor is for convenience; in many situations, it's less
    /// efficient than the array-based constructors.
    UT_DMatrix4(fpreal64 val00, fpreal64 val01, fpreal64 val02, fpreal64 val03,
                fpreal64 val10, fpreal64 val11, fpreal64 val12, fpreal64 val13,
                fpreal64 val20, fpreal64 val21, fpreal64 val22, fpreal64 val23,
                fpreal64 val30, fpreal64 val31, fpreal64 val32, fpreal64 val33)
    {
        matx[0][0] = val00; matx[0][1] = val01; matx[0][2] = val02;
                                                            matx[0][3] = val03;
        matx[1][0] = val10; matx[1][1] = val11; matx[1][2] = val12;
                                                            matx[1][3] = val13;
        matx[2][0] = val20; matx[2][1] = val21; matx[2][2] = val22;
                                                            matx[2][3] = val23;
        matx[3][0] = val30; matx[3][1] = val31; matx[3][2] = val32;
                                                            matx[3][3] = val33;
    }

    // Default copy constructor is fine.
    // UT_DMatrix4(const UT_DMatrix4 &m);

    /// Conversion constructor.
    explicit UT_DMatrix4(const UT_Matrix4 &m);

    // Default assignment operator is fine.
    //UT_DMatrix4       &operator= (const UT_DMatrix4 &m);

    /// Conversion operator that expands a 3x3 into a 4x4 matrix by adding a
    /// row and column of zeroes, except the diagonal element which is 1.
    // @{
    UT_DMatrix4 &operator=(const UT_Matrix3 &m);
    UT_DMatrix4 &operator=(const UT_Matrix4 &m);
    UT_DMatrix4 &operator=(const UT_DMatrix3 &m);
    // @}

    UT_DMatrix4  operator-() const
                 {
                     return UT_DMatrix4(
                        -matx[0][0], -matx[0][1], -matx[0][2], -matx[0][3],
                        -matx[1][0], -matx[1][1], -matx[1][2], -matx[1][3],
                        -matx[2][0], -matx[2][1], -matx[2][2], -matx[2][3],
                        -matx[3][0], -matx[3][1], -matx[3][2], -matx[3][3]);
                 }
    
    inline
    UT_DMatrix4 &operator+=(const UT_DMatrix4 &m)
                {
                    matx[0][0]+=m.matx[0][0];   matx[0][1]+=m.matx[0][1]; 
                    matx[0][2]+=m.matx[0][2]; matx[0][3]+=m.matx[0][3];

                    matx[1][0]+=m.matx[1][0];   matx[1][1]+=m.matx[1][1]; 
                    matx[1][2]+=m.matx[1][2]; matx[1][3]+=m.matx[1][3];

                    matx[2][0]+=m.matx[2][0];   matx[2][1]+=m.matx[2][1]; 
                    matx[2][2]+=m.matx[2][2]; matx[2][3]+=m.matx[2][3];

                    matx[3][0]+=m.matx[3][0];   matx[3][1]+=m.matx[3][1];
                    matx[3][2]+=m.matx[3][2]; matx[3][3]+=m.matx[3][3];
                    return *this;
                }
    inline
    UT_DMatrix4 &operator-=(const UT_DMatrix4 &m)
                {
                    matx[0][0]-=m.matx[0][0];   matx[0][1]-=m.matx[0][1]; 
                    matx[0][2]-=m.matx[0][2]; matx[0][3]-=m.matx[0][3];

                    matx[1][0]-=m.matx[1][0];   matx[1][1]-=m.matx[1][1]; 
                    matx[1][2]-=m.matx[1][2]; matx[1][3]-=m.matx[1][3];

                    matx[2][0]-=m.matx[2][0];   matx[2][1]-=m.matx[2][1]; 
                    matx[2][2]-=m.matx[2][2]; matx[2][3]-=m.matx[2][3];

                    matx[3][0]-=m.matx[3][0];   matx[3][1]-=m.matx[3][1];
                    matx[3][2]-=m.matx[3][2]; matx[3][3]-=m.matx[3][3];
                    return *this;
                }
    UT_DMatrix4 &operator*=(const UT_DMatrix4 &m);

    // test for exact floating point equality.  
    // for equality within a threshold, see isEqual()
    unsigned    operator==(const UT_DMatrix4 &m) const
                {
                    return (&m == this) ? 1U : (
                    matx[0][0]==m.matx[0][0] && matx[0][1]==m.matx[0][1] &&
                    matx[0][2]==m.matx[0][2] && matx[0][3]==m.matx[0][3] &&

                    matx[1][0]==m.matx[1][0] && matx[1][1]==m.matx[1][1] &&
                    matx[1][2]==m.matx[1][2] && matx[1][3]==m.matx[1][3] &&

                    matx[2][0]==m.matx[2][0] && matx[2][1]==m.matx[2][1] &&
                    matx[2][2]==m.matx[2][2] && matx[2][3]==m.matx[2][3] &&

                    matx[3][0]==m.matx[3][0] && matx[3][1]==m.matx[3][1] &&
                    matx[3][2]==m.matx[3][2] && matx[3][3]==m.matx[3][3] );
                }

    unsigned    operator!=(const UT_DMatrix4 &m) const
                {
                    return !(*this == m);
                }       

    // Scalar operators:
    UT_DMatrix4  &operator= (fpreal v)
                {
                    matx[0][0]= v; matx[0][1]= 0; matx[0][2]= 0; matx[0][3]= 0;
                    matx[1][0]= 0; matx[1][1]= v; matx[1][2]= 0; matx[1][3]= 0;
                    matx[2][0]= 0; matx[2][1]= 0; matx[2][2]= v; matx[2][3]= 0;
                    matx[3][0]= 0; matx[3][1]= 0; matx[3][2]= 0; matx[3][3]= v;
                    return *this;
                }
    UT_DMatrix4  &operator*=(fpreal64 scalar)
                {
                    matx[0][0]*=scalar; matx[0][1]*=scalar; 
                    matx[0][2]*=scalar; matx[0][3]*=scalar;

                    matx[1][0]*=scalar; matx[1][1]*=scalar; 
                    matx[1][2]*=scalar; matx[1][3]*=scalar;

                    matx[2][0]*=scalar; matx[2][1]*=scalar; 
                    matx[2][2]*=scalar; matx[2][3]*=scalar;

                    matx[3][0]*=scalar; matx[3][1]*=scalar; 
                    matx[3][2]*=scalar; matx[3][3]*=scalar;
                    return *this;
                }
    UT_DMatrix4  &operator/=(fpreal64 scalar)
                {
                    return operator*=( 1/scalar );
                }

    // Vector4 operators:
    inline
    UT_DMatrix4 &operator=(const UT_Vector4 &vec);
    inline
    UT_DMatrix4 &operator+=(const UT_Vector4 &vec);
    inline
    UT_DMatrix4 &operator-=(const UT_Vector4 &vec);

    // Other matrix operations:
    // left multiply the matrix by a matrix 
    // (operator*= does right-multiplication)
    void        leftMult( const UT_DMatrix4 &m );
    void        preMultiply(const UT_DMatrix4 &m)       { leftMult(m); }

    // Return the cofactor of the matrix, ie the determinant of the 3x3
    // submatrix that results from removing row 'k' and column 'l' from the 
    // 4x4.
    fpreal64    coFactor(int k, int l) const;

    fpreal64    determinant() const
                {
                    return(matx[0][0]*coFactor(0,0) +
                           matx[0][1]*coFactor(0,1) +
                           matx[0][2]*coFactor(0,2) +
                           matx[0][3]*coFactor(0,3) );
                }
    /// Compute determinant of the upper-left 3x3 sub-matrix
    fpreal64    determinant3() const
                {
                    return(matx[0][0]*
                               (matx[1][1]*matx[2][2]-matx[1][2]*matx[2][1]) +
                           matx[0][1]*
                               (matx[1][2]*matx[2][0]-matx[1][0]*matx[2][2]) +
                           matx[0][2]*
                               (matx[1][0]*matx[2][1]-matx[1][1]*matx[2][0]) );

                }
    fpreal64    trace() const
                { return matx[0][0]+matx[1][1]+matx[2][2]+matx[3][3]; }

    /// Invert this matrix and return 0 if OK, 1 if singular.
    int         invert();
    /// Invert the matrix and return 0 if OK, 1 if singular.
    /// Puts the inverted matrix in m, and leaves this matrix unchanged.
    int         invert(UT_DMatrix4 &m)const;
    int         invertKramer();
    int         invertKramer(UT_DMatrix4 &m)const;


    // Computes a transform to orient to a given direction (v) at a given
    // position (p) and with a scale (s). A secondary scale (s3) can be used to
    // specify the scale in x, y and z.  The up vector (up) is optional and
    // will orient the matrix to this up vector. If no up vector is given, the
    // z axis will be oriented to point in the v direction. If a quaternion (q)
    // is specified, the orientation will be additionally transformed by the
    // rotation specified by the quaternion. If a translation (tr) is
    // specified, the entire frame of reference will be moved by this
    // translation (unaffected by the scale or rotation).  If an orientation
    // (or) is specified, the orientation (using the velocity and up vector)
    // will not be performed and this orientation will instead be used to
    // define an original orientation.

    void instance(const UT_Vector3& p, const UT_Vector3& v, fpreal64 s,
                  const UT_Vector3* s3,
                  const UT_Vector3* up, const UT_Quaternion* q,
                  const UT_Vector3* tr, const UT_Quaternion* orient);

    // Transpose this matri or return its traanspose.
    void        transpose(void)
                {
                    fpreal64 tmp;
                    tmp=matx[0][1];  matx[0][1]=matx[1][0];  matx[1][0]=tmp;
                    tmp=matx[0][2];  matx[0][2]=matx[2][0];  matx[2][0]=tmp;
                    tmp=matx[0][3];  matx[0][3]=matx[3][0];  matx[3][0]=tmp;
                    tmp=matx[1][2];  matx[1][2]=matx[2][1];  matx[2][1]=tmp;
                    tmp=matx[1][3];  matx[1][3]=matx[3][1];  matx[3][1]=tmp;
                    tmp=matx[2][3];  matx[2][3]=matx[3][2];  matx[3][2]=tmp;
                }
    UT_DMatrix4 transpose(void) const;

    // check for equality within a tolerance level
    unsigned    isEqual( const UT_DMatrix4 &m, 
                         fpreal64 tolerance=UT_DTOLERANCE ) const
                {
                    return (&m == this) ? 1U : (
                    UTisEqual( matx[0][0], m.matx[0][0], tolerance ) &&
                    UTisEqual( matx[0][1], m.matx[0][1], tolerance ) &&
                    UTisEqual( matx[0][2], m.matx[0][2], tolerance ) &&
                    UTisEqual( matx[0][3], m.matx[0][3], tolerance ) &&

                    UTisEqual( matx[1][0], m.matx[1][0], tolerance ) &&
                    UTisEqual( matx[1][1], m.matx[1][1], tolerance ) &&
                    UTisEqual( matx[1][2], m.matx[1][2], tolerance ) &&
                    UTisEqual( matx[1][3], m.matx[1][3], tolerance ) &&

                    UTisEqual( matx[2][0], m.matx[2][0], tolerance ) &&
                    UTisEqual( matx[2][1], m.matx[2][1], tolerance ) &&
                    UTisEqual( matx[2][2], m.matx[2][2], tolerance ) &&
                    UTisEqual( matx[2][3], m.matx[2][3], tolerance ) &&

                    UTisEqual( matx[3][0], m.matx[3][0], tolerance ) &&
                    UTisEqual( matx[3][1], m.matx[3][1], tolerance ) &&
                    UTisEqual( matx[3][2], m.matx[3][2], tolerance ) &&
                    UTisEqual( matx[3][3], m.matx[3][3], tolerance ) );
                }

    // Multiply this matrix by a 3x3 rotation matrix determined by the axis 
    // and angle of rotation. If 'norm' is not 0, the axis vector is 
    // normalized before computing the rotation matrix. rotationMat() returns
    // a rotation matrix, and could as well be defined as a free floating
    // function.
    void        rotate(UT_Vector3 &axis, fpreal64 theta, int norm=1);
    UT_DMatrix4 rotate(UT_Vector3 &axis, fpreal64 theta, int norm=1) const;
    static UT_DMatrix4 rotationMat(UT_Vector3 &axis, fpreal64 theta, int norm=1);
    void        rotate(UT_Axis3::axis a, fpreal64 theta);
    UT_DMatrix4 rotate(UT_Axis3::axis a, fpreal64 theta) const;
    static UT_DMatrix4 rotationMat(UT_Axis3::axis a, fpreal64 theta);

    // Same as above, only pre-multiply this matrix by the rotation matx:
    void        prerotate(UT_Vector3 &axis, fpreal64 theta, int norm=1);
    UT_DMatrix4 prerotate(UT_Vector3 &axis, fpreal64 theta, int norm=1) const;
    void        prerotate(UT_Axis3::axis a, fpreal64 theta);
    UT_DMatrix4 prerotate(UT_Axis3::axis a, fpreal64 theta) const;

    // Rotate by rx, ry, rz radians around the three basic axes in the order
    // given by UT_XformOrder.
    void        rotate(fpreal64 rx, fpreal64 ry, fpreal64 rz,
                       const UT_XformOrder &ord);
    UT_DMatrix4 rotate(fpreal64 rx, fpreal64 ry, fpreal64 rz, 
                       const UT_XformOrder &ord) const;

    // Same as above, but pre-multiply
    void        prerotate(fpreal64 rx, fpreal64 ry, fpreal64 rz, 
                          const UT_XformOrder &ord);
    UT_DMatrix4 prerotate(fpreal64 rx, fpreal64 ry, fpreal64 rz, 
                          const UT_XformOrder &ord) const;

    // Multiply this matrix by a scale matrix with diagonal (sx, sy, sz):
    void        scale(fpreal64 sx, fpreal64 sy, fpreal64 sz, fpreal64 sw = 1.)
                {
                    matx[0][0] *= sx; matx[0][1] *= sy; 
                    matx[0][2] *= sz; matx[0][3] *= sw;

                    matx[1][0] *= sx; matx[1][1] *= sy; 
                    matx[1][2] *= sz; matx[1][3] *= sw;

                    matx[2][0] *= sx; matx[2][1] *= sy; 
                    matx[2][2] *= sz; matx[2][3] *= sw;

                    matx[3][0] *= sx; matx[3][1] *= sy;
                    matx[3][2] *= sz; matx[3][3] *= sw;
                }
    UT_DMatrix4 scale(fpreal64 sx, fpreal64 sy, fpreal64 sz, fpreal64 sw = 1.) const;

    // Same as above, only pre-multiply this matrix by the scale matx:
    void        prescale(fpreal64 sx, fpreal64 sy, fpreal64 sz, fpreal64 sw = 1.)
                {
                    matx[0][0] *= sx; matx[1][0] *= sy; 
                    matx[2][0] *= sz; matx[3][0] *= sw;

                    matx[0][1] *= sx; matx[1][1] *= sy; 
                    matx[2][1] *= sz; matx[3][1] *= sw;

                    matx[0][2] *= sx; matx[1][2] *= sy; 
                    matx[2][2] *= sz; matx[3][2] *= sw;

                    matx[0][3] *= sx; matx[1][3] *= sy; 
                    matx[2][3] *= sz; matx[3][3] *= sw;
                }
    UT_DMatrix4 prescale(fpreal64 sx, fpreal64 sy, fpreal64 sz, fpreal64 sw = 1.) const;

    // post-multiply this matrix by the shear matrix formed by (sxy, sxz, syz)
    void        shear(fpreal64 s_xy, fpreal64 s_xz, fpreal64 s_yz)
    {
        matx[0][0] += matx[0][1]*s_xy + matx[0][2]*s_xz;
        matx[0][1] += matx[0][2]*s_yz;

        matx[1][0] += matx[1][1]*s_xy + matx[1][2]*s_xz;
        matx[1][1] += matx[1][2]*s_yz;

        matx[2][0] += matx[2][1]*s_xy + matx[2][2]*s_xz;
        matx[2][1] += matx[2][2]*s_yz;

        matx[3][0] += matx[3][1]*s_xy + matx[3][2]*s_xz;
        matx[3][1] += matx[3][2]*s_yz;
    }

    // Multiply this matrix by a translation matrix determined by dx, dy, dz.
    void        translate(fpreal64 dx, fpreal64 dy, fpreal64 dz = 0.)
                {
                    fpreal64 a;
                    a = matx[0][3];
                    matx[0][0] += a*dx; matx[0][1] += a*dy; matx[0][2] += a*dz;
                    a = matx[1][3];
                    matx[1][0] += a*dx; matx[1][1] += a*dy; matx[1][2] += a*dz;
                    a = matx[2][3];
                    matx[2][0] += a*dx; matx[2][1] += a*dy; matx[2][2] += a*dz;
                    a = matx[3][3];
                    matx[3][0] += a*dx; matx[3][1] += a*dy; matx[3][2] += a*dz;
                }
    UT_DMatrix4 translate(fpreal64 dx, fpreal64 dy, fpreal64 dz = 0.) const;

    // Same as above, only pre-multiply this matrix by the translation matx:
    void        pretranslate(fpreal64 dx, fpreal64 dy, fpreal64 dz = 0.)
                {
                    matx[3][0] += matx[0][0]*dx + matx[1][0]*dy + matx[2][0]*dz;
                    matx[3][1] += matx[0][1]*dx + matx[1][1]*dy + matx[2][1]*dz;
                    matx[3][2] += matx[0][2]*dx + matx[1][2]*dy + matx[2][2]*dz;
                    matx[3][3] += matx[0][3]*dx + matx[1][3]*dy + matx[2][3]*dz;
                }
    UT_DMatrix4 pretranslate(fpreal64 dx, fpreal64 dy, fpreal64 dz = 0.) const;

    // Space change operation: right multiply this matrix by the 3x3 matrix
    // of the transformation which moves the vector space defined by
    // (iSrc, jSrc, cross(iSrc,jSrc)) into the space defined by
    // (iDest, jDest, cross(iDest,jDest)). iSrc, jSrc, iDest, and jDest will
    // be normalized before the operation if norm is 1. This matrix transforms
    // iSrc into iDest, and jSrc into jDest.
    void        changeSpace(UT_Vector3 &iSrc, UT_Vector3 &jSrc,
                            UT_Vector3 &iDest,UT_Vector3 &jDest,
                            int norm=1);
    UT_DMatrix4 changeSpace(UT_Vector3 &iSrc, UT_Vector3 &jSrc,
                            UT_Vector3 &iDest,UT_Vector3 &jDest,
                            int norm=1) const;

    // Multiply this matrix by the general transform matrix built from 
    // translations (tx,ty,tz), degree rotations (rx,ry,rz), scales (sx,sy,sz),
    // and possibly a pivot point (px,py,pz). The second methos leaves us
    // unchanged, and returns a new (this*xform) instead. The order of the 
    // multiplies (SRT, RST, RxRyRz, etc) is stored in 'order'. Normally you
    // will ignore the 'reverse' parameter, which tells us to build the
    // matrix from last to first xform, and to apply some inverses to 
    // txyz, rxyz, and sxyz.
    void        xform(const UT_XformOrder &order, 
                      fpreal64 tx=0., fpreal64 ty=0., fpreal64 tz=0.,
                      fpreal64 rx=0., fpreal64 ry=0., fpreal64 rz=0.,
                      fpreal64 sx=1., fpreal64 sy=1., fpreal64 sz=1.,
                      fpreal64 px=0., fpreal64 py=0., fpreal64 pz=0.,
                      int reverse=0);
    UT_DMatrix4 xform(const UT_XformOrder &order, 
                      fpreal64 tx=0., fpreal64 ty=0., fpreal64 tz=0.,
                      fpreal64 rx=0., fpreal64 ry=0., fpreal64 rz=0.,
                      fpreal64 sx=1., fpreal64 sy=1., fpreal64 sz=1.,
                      fpreal64 px=0., fpreal64 py=0., fpreal64 pz=0.,
                      int reverse=0) const;

    // This version handles shears as well
    void        xform(const UT_XformOrder &order, 
                      fpreal64 tx, fpreal64 ty, fpreal64 tz,
                      fpreal64 rx, fpreal64 ry, fpreal64 rz,
                      fpreal64 sx, fpreal64 sy, fpreal64 sz,
                      fpreal64 s_xy, fpreal64 s_xz, fpreal64 s_yz,
                      fpreal64 px, fpreal64 py, fpreal64 pz,
                      int reverse=0);

    // These versions of xform and rotate can be used to apply selected parts
    // of a transform. The applyType controls how far into the xform
    // it should go. For example: applyXform(order, BEFORE_EQUAL, 'T', ...)
    // will apply all the components of the transform up to and equal to the
    // translates (depending on UT_XformOrder)
    //
    // For xform char can be T, R, S, P, or p (P is for pivot)
    // For rotate char can be X, Y, or Z
    //
    // TODO add a reverse option

    enum applyType { BEFORE=1, EQUAL=2, AFTER=4, BEFORE_EQUAL=4, AFTER_EQUAL=6};
    void                 xform(const UT_XformOrder &order,
                               applyType type, char limit,
                               fpreal64 tx, fpreal64 ty, fpreal64 tz,
                               fpreal64 rx, fpreal64 ry, fpreal64 rz,
                               fpreal64 sx, fpreal64 sy, fpreal64 sz,
                               fpreal64 px, fpreal64 py, fpreal64 pz);

    void                 rotate(const UT_XformOrder &order,
                                applyType type, char limit,
                                fpreal64 rx, fpreal64 ry, fpreal64 rz);

    // extract only the translates from the matrix;
    inline
    void        getTranslates(UT_Vector3 &translates) const;
    inline
    void        setTranslates(const UT_Vector3 &translates);

    // This is a super-crack that returns the translation, scale, and radian
    // rotation vectors given a transformation order and a valid xform matrix.
    // It returns 0 if succesful, and non-zero otherwise: 1 if the embedded
    // rotation matrix is invalid, 2 if the rotation order is invalid, 
    // and 3 for other problems. If any of the scaling values is 0, the method 
    // returns all zeroes, and a 0 return value.
    int         explode(const UT_XformOrder &order,
                        UT_Vector3 &r, UT_Vector3 &s, UT_Vector3 &t,
                        UT_Vector3 *shears = 0) const;

    // This version of explode returns the t, r, and s, given a pivot value.
    int         explode(const UT_XformOrder &order,
                        UT_Vector3 &r, UT_Vector3 &s,
                        UT_Vector3 &t, const UT_Vector3 &p,
                        UT_Vector3 *shears = 0) const;

    void        extractRotate(UT_Matrix3 &dst) const;
    void        extractRotate(UT_DMatrix3 &dst) const;

    // Right multiply this matrix by a 3x3 matrix which scales by a given 
    // amount along the direction of vector v. When applied to a vector w,
    // the stretched matrix (*this) stretches w in v by the amount given.
    // If norm is non-zero, v will be normalized prior to the operation.
    void        stretch(UT_Vector3 &v, fpreal64 amount, int norm=1);
    UT_DMatrix4 stretch(UT_Vector3 &v, fpreal64 amount, int norm=1) const;

    fpreal64    dot(unsigned i, unsigned j) const
                {
                    return (i <= 3 && j <= 3) ?
                    matx[i][0]*matx[j][0] + matx[i][1]*matx[j][1] +
                    matx[i][2]*matx[j][2] + matx[i][3]*matx[j][3] : 0;
                }

    /// Set the matrix to identity
    void        identity()      { *this = 1; }
    /// Set the matrix to zero
    void        zero()          { *this = 0; }

    int         isIdentity() const
                {
                    // NB: DO NOT USE TOLERANCES!
                    return( 
                        matx[0][0]==1.0 && matx[0][1]==0.0 
                            && matx[0][2]==0.0 && matx[0][3]==0.0 &&
                        matx[1][0]==0.0 && matx[1][1]==1.0 
                            && matx[1][2]==0.0 && matx[1][3]==0.0 && 
                        matx[2][0]==0.0 && matx[2][1]==0.0 
                            && matx[2][2]==1.0 && matx[2][3]==0.0 &&
                        matx[3][0]==0.0 && matx[3][1]==0.0 
                            && matx[3][2]==0.0 && matx[3][3]==1.0 
                        );
                }

    /// Return the raw matrix data.
    const fpreal64      *data(void) const
                         { return (const fpreal64 *)&matx[0][0]; }
    fpreal64            *data(void)
                         { return (fpreal64 *)&matx[0][0]; }
    /// Return a matrix entry. No bounds checking on subscripts.
    // @{
    inline
    fpreal64    &operator()(unsigned row, unsigned col)
                 {
                     UT_ASSERT_P(row < 4 && col < 4);
                     return matx[row][col];
                 }
    inline
    fpreal64     operator()(unsigned row, unsigned col) const
                 {
                     UT_ASSERT_P(row < 4 && col < 4);
                     return matx[row][col];
                 }
    // @}

    /// Return a matrix row. No bounds checking on subscript.
    // @{
    fpreal64    *operator()(unsigned row)
                 {
                     UT_ASSERT_P(row < 4);
                     return matx[row];
                 }
    const fpreal64 *operator()(unsigned row) const
                 {
                     UT_ASSERT_P(row < 4);
                     return matx[row];
                 }
    inline
    UT_Vector4  operator[](unsigned row) const;
    // @}
    
    /// Euclidean or Frobenius norm of a matrix.
    /// Does sqrt(sum(a_ij ^2))
    fpreal64     getEuclideanNorm() const
                 { return SYSsqrt(getEuclideanNorm2()); }
    /// Euclidean norm squared.
    fpreal64     getEuclideanNorm2() const;


    // I/O methods.  Return 0 if read/write successful, -1 if unsuccessful.
    int                  save(ostream &os, int binary) const;
    bool                 load(UT_IStream &is);

    void                 outAsciiNoName(ostream &os) const;
    bool                 inAsciiNoName(UT_IStream &is);

    static const UT_DMatrix4 &getIdentityMatrix();

    // I/O friends:
    friend ostream      &operator<<(ostream &os, const UT_DMatrix4 &v)
                        {
                            os << className() << ' ';
                            v.outAsciiNoName(os);
                            return os;
                        }
    // The matrix data:
    fpreal64 matx[4][4];

private:
    static const char   *className(void);
};

#include "UT_Vector4.h"

inline
UT_DMatrix4 &UT_DMatrix4::operator=(const UT_Vector4 &vec)
{
    matx[0][0] = matx[0][1] = matx[0][2] = matx[0][3] = vec.x();
    matx[1][0] = matx[1][1] = matx[1][2] = matx[1][3] = vec.y();
    matx[2][0] = matx[2][1] = matx[2][2] = matx[2][3] = vec.z();
    matx[3][0] = matx[3][1] = matx[3][2] = matx[3][3] = vec.w();
    return *this;
}

inline
UT_DMatrix4 &UT_DMatrix4::operator+=(const UT_Vector4 &vec)
{
    fpreal64 x = vec.x(); fpreal64 y = vec.y();
    fpreal64 z = vec.z(); fpreal64 w = vec.w();
    matx[0][0]+=x; matx[0][1]+=x; matx[0][2]+=x; matx[0][3]+=x;
    matx[1][0]+=y; matx[1][1]+=y; matx[1][2]+=y; matx[1][3]+=y;
    matx[2][0]+=z; matx[2][1]+=z; matx[2][2]+=z; matx[2][3]+=z;
    matx[3][0]+=w; matx[3][1]+=w; matx[3][2]+=w; matx[3][3]+=w;
    return *this;
}

inline
UT_DMatrix4 &UT_DMatrix4::operator-=(const UT_Vector4 &vec)
{
    fpreal64 x = vec.x(); fpreal64 y = vec.y();
    fpreal64 z = vec.z(); fpreal64 w = vec.w();
    matx[0][0]-=x; matx[0][1]-=x; matx[0][2]-=x; matx[0][3]-=x;
    matx[1][0]-=y; matx[1][1]-=y; matx[1][2]-=y; matx[1][3]-=y;
    matx[2][0]-=z; matx[2][1]-=z; matx[2][2]-=z; matx[2][3]-=z;
    matx[3][0]-=w; matx[3][1]-=w; matx[3][2]-=w; matx[3][3]-=w;
    return *this;
}

inline
UT_Vector4  UT_DMatrix4::operator[](unsigned row) const
{
    return UT_Vector4(matx[row]);
}

inline void UT_DMatrix4::getTranslates(UT_Vector3 &translates) const
{
    translates.x() = (fpreal) matx[3][0];
    translates.y() = (fpreal) matx[3][1];
    translates.z() = (fpreal) matx[3][2];
}

inline void UT_DMatrix4::setTranslates(const UT_Vector3 &translates)
{
    matx[3][0] = translates.x();
    matx[3][1] = translates.y();
    matx[3][2] = translates.z();
}


// Free floating functions:
inline  
UT_DMatrix4      operator*(fpreal64 sc, const UT_DMatrix4 &m1)
{
    return m1*sc;
}

inline  
UT_DMatrix4      operator/(const UT_DMatrix4 &m1, fpreal64 scalar)
{
    return m1 * (1.0/scalar);
}


#endif

---