UT/UT_DMatrix3.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_DMatrix3_h__
#define __UT_DMatrix3_h__

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

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

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

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

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

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

    // Conversion constructors.
    explicit UT_DMatrix3(const UT_DMatrix4 &m);
    explicit UT_DMatrix3(const UT_Matrix3 &m);

    // Default assignment operator is fine.
    // UT_DMatrix3      &operator= (const UT_DMatrix3 &m)
    
    /// Conversion operator that returns a 3x3 from a 4x4 matrix by ignoring the
    /// last row and last column.
    // @{
    UT_DMatrix3 &operator=(const UT_Matrix4 &m);
    UT_DMatrix3 &operator=(const UT_DMatrix4 &m);
    // @}
    UT_DMatrix3 &operator=(const UT_Matrix3 &m);

    UT_DMatrix3  operator-() const
                 {
                     return UT_DMatrix3(-matx[0][0], -matx[0][1], -matx[0][2],
                                        -matx[1][0], -matx[1][1], -matx[1][2],
                                        -matx[2][0], -matx[2][1], -matx[2][2]);
                 }
    
    // Does += k * m
    void         addScaledMat(fpreal64 k, const UT_DMatrix3 &m)
                {
                    matx[0][0]+=k*m.matx[0][0];
                    matx[0][1]+=k*m.matx[0][1];
                    matx[0][2]+=k*m.matx[0][2];

                    matx[1][0]+=k*m.matx[1][0];
                    matx[1][1]+=k*m.matx[1][1];
                    matx[1][2]+=k*m.matx[1][2];

                    matx[2][0]+=k*m.matx[2][0];
                    matx[2][1]+=k*m.matx[2][1];
                    matx[2][2]+=k*m.matx[2][2];
                }
    inline
    UT_DMatrix3 &operator+=(const UT_DMatrix3 &m)
                {
                    matx[0][0]+=m.matx[0][0];
                    matx[0][1]+=m.matx[0][1];
                    matx[0][2]+=m.matx[0][2];

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

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

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

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

    void        multiply3(const UT_Matrix4 &m);
    void        multiply3(const UT_DMatrix4 &m);

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

    // Scalar operators:
    UT_DMatrix3 &operator= (fpreal64 val)
                {
                    matx[0][0] = val; matx[0][1] = 0; matx[0][2] = 0;
                    matx[1][0] = 0; matx[1][1] = val; matx[1][2] = 0;
                    matx[2][0] = 0; matx[2][1] = 0; matx[2][2] = val;
                    return *this;
                }
    UT_DMatrix3 &operator*=(fpreal64 scalar)
                {
                    matx[0][0]*=scalar; matx[0][1]*=scalar; matx[0][2]*=scalar;
                    matx[1][0]*=scalar; matx[1][1]*=scalar; matx[1][2]*=scalar;
                    matx[2][0]*=scalar; matx[2][1]*=scalar; matx[2][2]*=scalar;
                    return *this;
                }
    UT_DMatrix3 &operator/=(fpreal64 scalar)
                {
                    return operator*=( 1/scalar );
                }

    // check for equality within a tolerance level
    unsigned    isEqual( const UT_DMatrix3 &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[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[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 ) );
                }

    // Vector3 operators:
    inline
    UT_DMatrix3 &operator= (const UT_Vector3 &vec);
    inline
    UT_DMatrix3 &operator+=(const UT_Vector3 &vec);
    inline
    UT_DMatrix3 &operator-=(const UT_Vector3 &vec);

    // Scaled outer product update.
    // this += b * x * yT
    inline
    void         outerproductUpdate(fpreal64 b, 
                                    const UT_Vector3 &x, const UT_Vector3 &y);

    // Sets this matrix to the canonical 180 rotation matrix (about the y axis)
    // that is used in the orient() and dihedral() functions when given two
    // opposing vectors.
    void        arbitrary180rot();
    
    // Other matrix operations:
    // Compute the dihedral: return the matrix that computes the rotation around
    // the axes normal to both a and b, (their cross product), by the angle
    // which is between a and b. The resulting matrix maps a onto  b. If a and
    // b have the same direction, what comes out is the identity matrix. If a
    // and b are opposed, the rotation is undefined and the method returns a
    // vector c of zeroes (check with !c); if all is OK, c != 0.
    // The transformation is a symmetry around vector a, then a symmetry around
    // (norm(a) and norm(b)). If a or b needs to be normalized, pass a non-zero 
    // value in norm.
    // The second function places the result in 'm', and returns 0 if a and b
    // are not opposed; otherwise, it returns 1.
    int         dihedral(UT_Vector3 &a, UT_Vector3 &b,
                                     int norm=1);

    // Compute a lookat matrix:  These functions will compute a rotation matrix
    // (yes, with all the properties of a rotation matrix), which will provide
    // the rotates needed for "from" to look at "to".  One method uses a "roll"
    // the other method uses an "up" vector to determine the orientation.  The
    // "roll" is not as well defined as using an "up" vector.  If the two points
    // are co-incident, the result will be an identity matrix.  If norm is set,
    // then the up vector will be normalized.
    // The lookat matrix will have the -Z axis pointing at the "to" point.
    // The y axis will be pointing "up"
    void        lookat(const UT_Vector3 &from, const UT_Vector3 &to,
                           fpreal64 roll = 0);
    void        lookat(const UT_Vector3 &from, const UT_Vector3 &to,
                           const UT_Vector3 &up);


    // Compute the transform matrix to orient given a direction and an up vector
    // If the direction and up are not orthogonal, then the up vector will be
    // shifted along the direction to make it orthogonal...  If the up and
    // direction vectors are co-linear, then an identity matrix will be
    // returned.
    // The second function will return 0 to indicate this.  The inputs may
    // be modified by the orient functions.
    int         orient(UT_Vector3 &dir, UT_Vector3 &up,
                                   int norm=1);


    // Computes a transform to orient to a given direction (v) 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.
    void                orient(const UT_Vector3& v, 
                                   fpreal64 pscale,
                                   const UT_Vector3* pscale3,
                                   const UT_Vector3* up,
                                   const UT_Quaternion* q);

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

    fpreal64    determinant() 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]; }

    /// Computes the eigenvalues.  Returns number of real eigenvalues
    /// found.  Uses UT_RootFinder::cubic
    int         eigenvalues(UT_Vector3 &r, UT_Vector3 &i) const;

    /// 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_DMatrix3 &m)const;
    int         invertKramer(void);
    int         invertKramer(UT_DMatrix3 &m)const;

    // Invert *this.
    // *this should be symmetric.
    // *this does NOT have to be non-singular.  Singular matrices will be
    // inverted by filtering out the singular component.
    // This means you can reliably solve:
    // M x = b
    // for the x that is closest to satisfying b
    // Returns 0 on success, 1 if something really bad happened.
    int         safeInvertSymmetric(fpreal64 tol=1e-6f);

    // Diagonalize *this.
    // *this should be symmetric.
    // The matrix is factored into:
    // *this = Rt D R
    // The diagonalization is done with a serios of jacobi rotations.
    // *this is unchanged by the operations.
    // Returns true if successful, false if reached maximum iterations rather
    // than desired tolerance.
    // Tolerance is with respect to the maximal element of the matrix.
    bool        diagonalizeSymmetric(UT_DMatrix3 &R, UT_DMatrix3 &D, fpreal64 tol=1e-6f) const;

    // Splits a matrix into a scale & rotation component.
    // R is the rotation component, S the inverse scale.
    // If this is A, finds the inverse square root of AtA, S. 
    // We then set R = AS to give us a rotation matrix.
    void        splitRotationScale(UT_DMatrix3 &R, UT_DMatrix3 &S, fpreal64 tol=1e-6f) const;

    // Transpose this matrix or return its transpose.
    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[1][2];  matx[1][2]=matx[2][1];  matx[2][1]=tmp;
                }
    UT_DMatrix3 transpose(void) const;

    bool        isSymmetric(fpreal64 tolerance = UT_FTOLERANCE) const;

    // 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. Theta is in radians.
    void        rotate(UT_Vector3 &axis, fpreal64 theta, int norm=1);
    UT_DMatrix3 rotate(UT_Vector3 &axis, fpreal64 theta, int norm=1) const;
    static UT_DMatrix3 rotationMat(UT_Vector3 &axis, fpreal64 theta, int norm=1);
    void        rotate(UT_Axis3::axis a, fpreal64 theta);
    UT_DMatrix3 rotate(UT_Axis3::axis a, fpreal64 theta) const;
    static UT_DMatrix3 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_DMatrix3 prerotate(UT_Vector3 &axis, fpreal64 theta, int norm=1) const;
    void        prerotate(UT_Axis3::axis a, fpreal64 theta);
    UT_DMatrix3 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        prerotate(fpreal64 rx, fpreal64 ry, fpreal64 rz,
                          const UT_XformOrder &order);
    UT_DMatrix3 prerotate(fpreal64 rx, fpreal64 ry, fpreal64 rz,
                          const UT_XformOrder &order) 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_DMatrix3 rotate(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)
                {
                    matx[0][0] *= sx; matx[0][1] *= sy; matx[0][2] *= sz;
                    matx[1][0] *= sx; matx[1][1] *= sy; matx[1][2] *= sz;
                    matx[2][0] *= sx; matx[2][1] *= sy; matx[2][2] *= sz;
                }
    UT_DMatrix3 scale(fpreal64 sx, fpreal64 sy, fpreal64 sz) const;

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

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

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

    // Generate the x, y, and z values of rotation. Return 0 if successful, and
    // a non-zero value otherwise: 1 if the matrix is not a valid rotation
    // matrix, 2 if rotOrder is invalid, and 3 for other errors. rotOrder must
    // be given as a UT_XformOrder permulation.
    // WARNING: The non-const method changes the matrix!
    int         crack(UT_Vector3 &rvec, const UT_XformOrder &order,
                            int remove_scales=1);
    int         crack(UT_Vector3 &rvec, const UT_XformOrder &order,
                            int remove_scales=1) const;

    // This method will condition the matrix so that it's a valid rotation
    // matrix (i.e. crackable).  Ideally, the scales should be removed first,
    // but this method will attempt to do this as well.  It will return the
    // conditioned scales if a vector is passed in.
    void        conditionRotate(UT_Vector3 *scales = 0);

    // Extract scales and shears from the matrix leaving the matrix as a nice
    // orthogonal rotation matrix.  The shears are:
    //                  XY, XZ, and YZ
    void        extractScales(UT_Vector3 &scales, UT_Vector3 *shears=0);

    // Shear the matrix according to the values.  This is equivalent to
    // multiplying the matrix by the shear matrix with the given value.
    void        shearXY(fpreal64 val);
    void        shearXZ(fpreal64 val);
    void        shearYZ(fpreal64 val);

    // 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;
    }

    // Solve a 3x3 system of equations whose parameters are held in this matrix,
    // and whose rhs constants are cx, cy, and cz. The method returns 0 if
    // the determinant is not 0, and 1 otherwise.
    int         solve(fpreal64 cx, fpreal64 cy, fpreal64 cz, UT_Vector3 &result) 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_DMatrix3 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), degree rotation (rz), scales (sx,sy),
    // and possibly a pivot point (px,py). The second methos leaves us
    // unchanged, and returns a new (this*xform) instead. The order of the 
    // multiplies (SRT, RST, Rz, 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 
    // txy, rz, and sxy.
    void        xform(const UT_XformOrder &order, 
                      fpreal64 tx=0., fpreal64 ty=0., fpreal64 rz=0.,
                      fpreal64 sx=1., fpreal64 sy=1., fpreal64 px=0., fpreal64 py=0.,
                      int reverse=0);
    UT_DMatrix3 xform(const UT_XformOrder &order, 
                      fpreal64 tx=0., fpreal64 ty=0., fpreal64 rz=0.,
                      fpreal64 sx=1., fpreal64 sy=1., fpreal64 px=0., fpreal64 py=0.,
                      int reverse=0) const;

    // this version handles shears as well
    void        xform(const UT_XformOrder &order, 
                  fpreal64 tx, fpreal64 ty, fpreal64 rz,
                  fpreal64 sx, fpreal64 sy, fpreal64 s_xy,
                  fpreal64 px, fpreal64 py,
                  int reverse=0);

    // 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_DMatrix3 stretch(UT_Vector3 &v, fpreal64 amount, int norm=1) const;

    // This does a test to see if a vector transformed by the matrix
    // will maintain it's length (i.e. there are no scales in the matrix)
    int         isNormalized() const;

    // Return the dot product between two rows i and j:
    fpreal64    dot(unsigned i, unsigned j) const
                {
                    return (i <= 2 && j <= 2) ?
                    matx[i][0]*matx[j][0] + matx[i][1]*matx[j][1] +
                                            matx[i][2]*matx[j][2] : 0;
                }

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

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

    /// Return a matrix row. No bounds checking on subscript.
    // @{
    fpreal64    *operator()(unsigned row)
                 {
                     UT_ASSERT_P(row < 3);
                     return matx[row];
                 }
    const fpreal64*operator()(unsigned row) const
                 {
                     UT_ASSERT_P(row < 3);
                     return matx[row];
                 }
    // @}

    /// 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;
    void                 outAsciiNoName(ostream &os) const;

    static const UT_DMatrix3 &getIdentityMatrix();

    // I/O friends:
    friend ostream      &operator<<(ostream &os, const UT_DMatrix3 &v)
                        {
                            os << className() << ' ';
                            v.outAsciiNoName(os);
                            return os;
                        }
    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[1][0]==0.0 && matx[1][1]==1.0
                            && matx[1][2]==0.0 &&
                        matx[2][0]==0.0 && matx[2][1]==0.0
                            && matx[2][2]==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]; }
    // @}

    int         isZero() const
                {
                    // NB: DO NOT USE TOLERANCES!
                    return(
                        matx[0][0]==0.0 && matx[0][1]==0.0
                            && matx[0][2]==0.0 &&
                        matx[1][0]==0.0 && matx[1][1]==0.0
                            && matx[1][2]==0.0 &&
                        matx[2][0]==0.0 && matx[2][1]==0.0
                            && matx[2][2]==0.0
                        );
                }
private:
    // Check this matrix to see if it is a valid rotation matrix, and permute
    // its elements depending on the rotOrder value in crack(). checkRot() 
    // returns 0 if OK and 1 otherwise.
    int                 checkRot() const;
    void                permute(int l0, int l1, int l2);

    static const char   *className(void);

    // The matrix data:
    fpreal64 matx[3][3];
};

#include "UT_Vector3.h"

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

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

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

inline
void     UT_DMatrix3::outerproductUpdate(fpreal64 b, 
                                         const UT_Vector3 &x,
                                         const UT_Vector3 &y)
{
    fpreal64 bx;
    bx = b * x(0);
    matx[0][0]+=bx*y(0);matx[0][1]+=bx*y(1);matx[0][2]+=bx*y(2);
    bx = b * x(1);
    matx[1][0]+=bx*y(0);matx[1][1]+=bx*y(1);matx[1][2]+=bx*y(2);
    bx = b * x(2);
    matx[2][0]+=bx*y(0);matx[2][1]+=bx*y(1);matx[2][2]+=bx*y(2);
}


// Free floating functions:
inline
UT_DMatrix3      operator+(const UT_Vector3 &vec, const UT_DMatrix3 &mat)
{
    return mat+vec;
}

inline
UT_DMatrix3      operator*(fpreal64 sc, const UT_DMatrix3 &m1)
{
    return m1*sc;
}

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


#endif

---