UT/UT_Matrix4.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_Matrix4_h__
#define __UT_Matrix4_h__

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

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

// Free floating operators that return a UT_Matrix4 object. 
UT_API UT_Matrix4       operator+(const UT_Matrix4 &m1,  const UT_Matrix4 &m2);
UT_API UT_Matrix4       operator-(const UT_Matrix4 &m1,  const UT_Matrix4 &m2);
UT_API UT_Matrix4       operator*(const UT_Matrix4 &m1,  const UT_Matrix4 &m2);
UT_API UT_Matrix4       operator+(const UT_Matrix4 &mat, const UT_Vector4 &vec);
UT_API UT_Matrix4       operator+(const UT_Vector4 &vec, const UT_Matrix4 &mat);
UT_API UT_Matrix4       operator-(const UT_Matrix4 &mat, const UT_Vector4 &vec);
UT_API UT_Matrix4       operator-(const UT_Vector4 &vec, const UT_Matrix4 &mat);
UT_API UT_Matrix4       operator+(const UT_Matrix4 &mat, fpreal sc);
inline UT_Matrix4       operator-(const UT_Matrix4 &mat, fpreal sc);
UT_API UT_Matrix4       operator*(const UT_Matrix4 &mat, fpreal sc);
inline UT_Matrix4       operator/(const UT_Matrix4 &mat, fpreal sc);
inline UT_Matrix4       operator+(fpreal sc, const UT_Matrix4 &mat);
UT_API UT_Matrix4       operator-(fpreal sc, const UT_Matrix4 &mat);
inline UT_Matrix4       operator*(fpreal sc, const UT_Matrix4 &mat);
UT_API UT_Matrix4       operator/(fpreal sc, const UT_Matrix4 &mat);

/// This class implements a 4x4 fpreal 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_Matrix4
{
public:
    /// Construct uninitialized matrix.
    UT_Matrix4()
    {
        UT_ASSERT_COMPILETIME(sizeof(UT_Matrix4) == 16 * sizeof(float));
    }
    /// Construct identity matrix, multipled by scalar.
    explicit UT_Matrix4(fpreal 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_Matrix4(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_Matrix4(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_Matrix4(fpreal32 val00, fpreal32 val01, fpreal32 val02, fpreal32 val03,
               fpreal32 val10, fpreal32 val11, fpreal32 val12, fpreal32 val13,
               fpreal32 val20, fpreal32 val21, fpreal32 val22, fpreal32 val23,
               fpreal32 val30, fpreal32 val31, fpreal32 val32, fpreal32 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_Matrix4(const UT_Matrix4 &m);

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

    // Default assignment operator is fine.
    //UT_Matrix4        &operator= (const UT_Matrix4 &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_Matrix4  &operator=(const UT_Matrix3 &m);
    UT_Matrix4  &operator=(const UT_DMatrix3 &m);
    // @}
    UT_Matrix4  &operator=(const UT_DMatrix4 &m);

    /// Conversion from a symmetric to a non symmetric matrix
    UT_Matrix4  &operator=(const UT_SymMatrix4 &m)
    {
        matx[0][0] = m.q00;     matx[0][1] = m.q01;     matx[0][2] = m.q02;
        matx[0][3] = m.q03;     matx[1][0] = m.q01;     matx[1][1] = m.q11;
        matx[1][2] = m.q12;     matx[1][3] = m.q13;     matx[2][0] = m.q02;
        matx[2][1] = m.q12;     matx[2][2] = m.q22;     matx[2][3] = m.q23;
        matx[3][0] = m.q03;     matx[3][1] = m.q13;     matx[3][2] = m.q23;
        matx[3][3] = m.q33;
        return *this;
    }

    UT_Matrix4   operator-() const
                 {
                     return UT_Matrix4(
                        -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_Matrix4   &operator+=(const UT_Matrix4 &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_Matrix4   &operator-=(const UT_Matrix4 &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_Matrix4  &operator*=(const UT_Matrix4 &m);

    // test for exact floating point equality.  
    // for equality within a threshold, see isEqual()
    unsigned    operator==(const UT_Matrix4 &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_Matrix4 &m) const
                {
                    return !(*this == m);
                }       

    // Scalar operators:
    UT_Matrix4  &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;
                }
    inline UT_Matrix4  &operator*=(fpreal 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;
                }
    inline UT_Matrix4  &operator/=(fpreal scalar)
                {
                    return operator*=( 1.0f/scalar );
                }

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

    // Other matrix operations:
    // left multiply the matrix by a matrix 
    // (operator*= does right-multiplication)
    void        leftMult( const UT_Matrix4 &m );
    void        preMultiply(const UT_Matrix4 &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.
    fpreal      coFactor(int k, int l) const;

    fpreal      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
    fpreal      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]) );

                }
    fpreal      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(fpreal tol = 0.0F);
    int         invertDouble();
    // @}

    /// 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_Matrix4 &m) const;
    int         invertDouble(UT_Matrix4 &m) const;
    // @}
    int         invertKramer(void);
    int         invertKramer(UT_Matrix4 &m)const;

    // Computes a transform to orient to a given direction (v) at a given
    // position (p) and with a scale (s). 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 (orient) 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.
    //
    // The matrix is scaled non-uniformly about each axis using s3, if s3
    // is non-zero.  A uniform scale of pscale is applied regardless, so if
    // s3 is non-zero, the x axis will be scaled by pscale * s3->x().
    void instance(const UT_Vector3& p, const UT_Vector3& v, fpreal s,
                  const UT_Vector3* s3,
                  const UT_Vector3* up, const UT_Quaternion* q,
                  const UT_Vector3* tr, const UT_Quaternion* orient);

    // Transpose this matrix or return its transpose.
    void        transpose(void)
                {
                    fpreal 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_Matrix4  transpose(void) const;

    // check for equality within a tolerance level
    unsigned    isEqual( const UT_Matrix4 &m, 
                         fpreal tolerance=UT_FTOLERANCE ) 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, fpreal theta, int norm=1);
    UT_Matrix4  rotate(UT_Vector3 &axis, fpreal theta, int norm=1) const;
    static UT_Matrix4 rotationMat(UT_Vector3 &axis, fpreal theta, int norm=1);
    void        rotate(UT_Axis3::axis a, fpreal theta);
    UT_Matrix4  rotate(UT_Axis3::axis a, fpreal theta) const;
    static UT_Matrix4 rotationMat(UT_Axis3::axis a, fpreal theta);

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

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

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

    // Multiply this matrix by a scale matrix with diagonal (sx, sy, sz):
    void        scale(fpreal sx, fpreal sy, fpreal sz, fpreal sw = 1.0f)
                {
                    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_Matrix4  scale(fpreal sx, fpreal sy, fpreal sz, fpreal sw = 1.0f) const;

    // Same as above, only pre-multiply this matrix by the scale matx:
    void        prescale(fpreal sx, fpreal sy, fpreal sz, fpreal sw = 1.0f)
                {
                    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_Matrix4  prescale(fpreal sx, fpreal sy, fpreal sz, fpreal sw = 1.0f) const;

    // post-multiply this matrix by the shear matrix formed by (sxy, sxz, syz)
    void        shear(fpreal s_xy, fpreal s_xz, fpreal 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(fpreal dx, fpreal dy, fpreal dz = 0.0f)
                {
                    fpreal 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_Matrix4  translate(fpreal dx, fpreal dy, fpreal dz = 0.0f) const;

    // Same as above, only pre-multiply this matrix by the translation matx:
    void        pretranslate(fpreal dx, fpreal dy, fpreal dz = 0.0f)
                {
                    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_Matrix4  pretranslate(fpreal dx, fpreal dy, fpreal dz = 0.0f) 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_Matrix4  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, 
                      fpreal tx=0.0f, fpreal ty=0.0f, fpreal tz=0.0f,
                      fpreal rx=0.0f, fpreal ry=0.0f, fpreal rz=0.0f,
                      fpreal sx=1.0f, fpreal sy=1.0f, fpreal sz=1.0f,
                      fpreal px=0.0f, fpreal py=0.0f, fpreal pz=0.0f,
                      int reverse=0);
    UT_Matrix4  xform(const UT_XformOrder &order, 
                      fpreal tx=0.0f, fpreal ty=0.0f, fpreal tz=0.0f,
                      fpreal rx=0.0f, fpreal ry=0.0f, fpreal rz=0.0f,
                      fpreal sx=1.0f, fpreal sy=1.0f, fpreal sz=1.0f,
                      fpreal px=0.0f, fpreal py=0.0f, fpreal pz=0.0f,
                      int reverse=0) const;

    // This version handles shears as well
    void        xform(const UT_XformOrder &order, 
                      fpreal tx, fpreal ty, fpreal tz,
                      fpreal rx, fpreal ry, fpreal rz,
                      fpreal sx, fpreal sy, fpreal sz,
                      fpreal s_xy, fpreal s_xz, fpreal s_yz,
                      fpreal px, fpreal py, fpreal 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,
                               fpreal tx, fpreal ty, fpreal tz,
                               fpreal rx, fpreal ry, fpreal rz,
                               fpreal sx, fpreal sy, fpreal sz,
                               fpreal px, fpreal py, fpreal pz);

    void                 rotate(const UT_XformOrder &order,
                                applyType type, char limit,
                                fpreal rx, fpreal ry, fpreal 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;

    // These versions treat the matrix as only containing a 2D 
    // transformation, that is in x, y, with a rotate around z.
    int         explode2D(const UT_XformOrder &order,
                          float &r, UT_Vector2 &s, UT_Vector2 &t,
                          float *shears = 0) const;

    int         explode2D(const UT_XformOrder &order,
                        float &r, UT_Vector2 &s,
                        UT_Vector2 &t, const UT_Vector2 &p,
                        float *shears = 0) const;

    void        extractRotate(UT_Matrix3 &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, fpreal amount, int norm=1);
    UT_Matrix4  stretch(UT_Vector3 &v, fpreal amount, int norm=1) const;

    fpreal      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] : (fpreal)0;
                }

    // Matrix += b * v1 * v2T
    void        outerproductUpdate(fpreal b, 
                             const UT_Vector4 &v1, const UT_Vector4 &v2);

    /// 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.0f && matx[0][1]==0.0f 
                            && matx[0][2]==0.0f && matx[0][3]==0.0f &&
                        matx[1][0]==0.0f && matx[1][1]==1.0f 
                            && matx[1][2]==0.0f && matx[1][3]==0.0f && 
                        matx[2][0]==0.0f && matx[2][1]==0.0f 
                            && matx[2][2]==1.0f && matx[2][3]==0.0f &&
                        matx[3][0]==0.0f && matx[3][1]==0.0f 
                            && matx[3][2]==0.0f && matx[3][3]==1.0f 
                        );
                }


    /// Return the raw matrix data.
    // @{
    const float *data(void) const       { return (const float *)&matx[0][0]; }
    float       *data(void)             { return (float *)&matx[0][0];       }
    // @}

    /// Return a matrix entry. No bounds checking on subscripts.
    // @{
    float       &operator()(unsigned row, unsigned col)
                 {
                     UT_ASSERT_P(row < 4 && col < 4);
                     return matx[row][col];
                 }
    float        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.
    // @{
    float       *operator()(unsigned row)
                 {
                     UT_ASSERT_P(row < 4);
                     return matx[row];
                 }
    const float *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))
    fpreal       getEuclideanNorm() const
                 { return SYSsqrt(getEuclideanNorm2()); }
    /// Euclidean norm squared.
    fpreal       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;

    static const UT_Matrix4 &getIdentityMatrix();

    // I/O friends:
    friend ostream      &operator<<(ostream &os, const UT_Matrix4 &v)
                        {
                            os << className() << ' ';
                            v.outAsciiNoName(os);
                            return os;
                        }

    // The matrix data:
    float matx[4][4];

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

#include "UT_Vector4.h"

// Vector4 operators:

inline
UT_Matrix4  &UT_Matrix4::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_Matrix4  &UT_Matrix4::operator+=(const UT_Vector4 &vec)
{
    fpreal x = vec.x(); fpreal y = vec.y();
    fpreal z = vec.z(); fpreal 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_Matrix4  &UT_Matrix4::operator-=(const UT_Vector4 &vec)
{
    fpreal x = vec.x(); fpreal y = vec.y();
    fpreal z = vec.z(); fpreal 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
void        UT_Matrix4::getTranslates(UT_Vector3 &translates) const
{
    translates.x() = matx[3][0];
    translates.y() = matx[3][1];
    translates.z() = matx[3][2];
}

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

inline
UT_Vector4  UT_Matrix4::operator[](unsigned row) const
{
    UT_ASSERT_P(row < 4);
    return UT_Vector4(matx[row]);
}

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

inline  
UT_Matrix4      operator/(const UT_Matrix4 &m1, fpreal scalar)
{
    return (m1 * (1.0f/scalar));
}

#endif

---