UT/UT_Matrix3.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_Matrix3_h__
#define __UT_Matrix3_h__

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

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

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

/// This class implements a 3x3 float 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_Matrix3
{
public:
    /// Construct uninitialized matrix.
    UT_Matrix3()
    {
        UT_ASSERT_COMPILETIME(sizeof(UT_Matrix3) == 9 * sizeof(float));
    }
    /// Construct identity matrix, multipled by scalar.
    explicit UT_Matrix3(fpreal 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_Matrix3(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_Matrix3(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_Matrix3(fpreal32 val00, fpreal32 val01, fpreal32 val02,
               fpreal32 val10, fpreal32 val11, fpreal32 val12,
               fpreal32 val20, fpreal32 val21, fpreal32 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_Matrix3(const UT_Matrix3 &m);
    explicit UT_Matrix3(const UT_Matrix4 &m);

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

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

    UT_Matrix3   operator-() const
                 {
                     return UT_Matrix3(-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
    inline void  addScaledMat(fpreal k, const UT_Matrix3 &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_Matrix3   &operator+=(const UT_Matrix3 &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_Matrix3   &operator-=(const UT_Matrix3 &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_Matrix3  &operator*=(const UT_Matrix3 &m);

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

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

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

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

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

    // Other matrix operations:
    // Computes the basis of a 3x3 matrix.  This is used by OBB code
    // to determine what orientation of a bounding box best represents
    // a cloud of points.  To do this, first fill this matrix with
    // the sum of p*pT for all points in the cloud.
    // The result is in this matrix itself.
    void                getBasis();

    // 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();
    
    // 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.
    static UT_Matrix3   dihedral(UT_Vector3 &a, UT_Vector3 &b,
                                     UT_Vector3 &c,int norm=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,
                           fpreal 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
    // (pscale).  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.
    // 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                orient(const UT_Vector3& v, 
                                   fpreal pscale, const UT_Vector3 *s3,
                                   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.
    fpreal      coFactor(int k, int l) const;

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

                }
    fpreal      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_Matrix3 &m)const;
    int         invertKramer();
    int         invertKramer(UT_Matrix3 &m)const;

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

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

    bool        isSymmetric(fpreal 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, fpreal theta, int norm=1);
    UT_Matrix3  rotate(UT_Vector3 &axis, fpreal theta, int norm=1) const;
    static UT_Matrix3 rotationMat(UT_Vector3 &axis, fpreal theta, int norm=1);
    void        rotate(UT_Axis3::axis a, fpreal theta);
    UT_Matrix3  rotate(UT_Axis3::axis a, fpreal theta) const;
    static UT_Matrix3 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_Matrix3  prerotate(UT_Vector3 &axis, fpreal theta, int norm=1) const;
    void        prerotate(UT_Axis3::axis a, fpreal theta);
    UT_Matrix3  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        prerotate(fpreal rx, fpreal ry, fpreal rz,
                          const UT_XformOrder &order);
    UT_Matrix3  prerotate(fpreal rx, fpreal ry, fpreal 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(fpreal rx, fpreal ry, fpreal rz,
                       const UT_XformOrder &ord);
    UT_Matrix3  rotate(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)
                {
                    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_Matrix3  scale(fpreal sx, fpreal sy, fpreal sz) const;

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

    // Multiply this matrix by a translation matrix determined by dx and dy.
    void        translate(fpreal dx, fpreal 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_Matrix3  translate(fpreal dx, fpreal dy) const;

    // Same as above, only pre-multiply this matrix by the translation matx:
    void        pretranslate(fpreal dx, fpreal 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_Matrix3  pretranslate(fpreal dx, fpreal 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;
    int         crack2D(float &r) 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);

    // Extract scales and shears from the matrix leaving the matrix as a nice
    // orthogonal rotation matrix (assuming it was only a 2D xform to begin
    // with).  The shears are:
    //                  XY
    void        extractScales2D(UT_Vector2 &scales, float *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(fpreal val);
    void        shearXZ(fpreal val);
    void        shearYZ(fpreal val);

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

    // 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(fpreal cx, fpreal cy, fpreal 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_Matrix3  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, 
                  fpreal tx=0.0f, fpreal ty=0.0f, fpreal rz=0.0f,
                  fpreal sx=1.0f, fpreal sy=1.0f,fpreal px=0.0f,fpreal py=0.0f,
                  int reverse=0);
    UT_Matrix3  xform(const UT_XformOrder &order, 
                  fpreal tx=0.0f, fpreal ty=0.0f, fpreal rz=0.0f,
                  fpreal sx=1.0f, fpreal sy=1.0f,fpreal px=0.0f,fpreal py=0.0f,
                  int reverse=0) const;

    // this version handles shears as well
    void        xform(const UT_XformOrder &order, 
                  fpreal tx, fpreal ty, fpreal rz,
                  fpreal sx, fpreal sy, fpreal s_xy,
                  fpreal px, fpreal 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, fpreal amount, int norm=1);
    UT_Matrix3  stretch(UT_Vector3 &v, fpreal 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:
    fpreal      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] : (fpreal)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.0f && matx[0][1]==0.0f &&
                            matx[0][2]==0.0f && matx[1][0]==0.0f &&
                            matx[1][1]==1.0f && matx[1][2]==0.0f &&
                            matx[2][0]==0.0f && matx[2][1]==0.0f &&
                            matx[2][2]==1.0f
                          );
                }
    
    /// Return the raw matrix data.
    // @{
    const fpreal32 *data(void) const    { return (const fpreal32 *)&matx[0][0];}
    fpreal32       *data(void)          { return (fpreal32 *)&matx[0][0]; }
    // @}

    /// Return a matrix entry. No bounds checking on subscripts.
    // @{
    inline float&operator()(unsigned row, unsigned col)
                 {
                     UT_ASSERT_P(row < 3 && col < 3);
                     return matx[row][col];
                 }
    inline float 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.
    // @{
    float       *operator()(unsigned row)
                 {
                     UT_ASSERT_P(row < 3);
                     return matx[row];
                 }
    const float *operator()(unsigned row) const
                 {
                     UT_ASSERT_P(row < 3);
                     return matx[row];
                 }
    UT_Vector3   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_Matrix3 &getIdentityMatrix();

    // I/O friends:
    friend ostream      &operator<<(ostream &os, const UT_Matrix3 &v)
                        {
                            os << className() << ' ';
                            v.outAsciiNoName(os);
                            return os;
                        }
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:
    fpreal32    matx[3][3];
};

#include "UT_Vector3.h"

inline
UT_Matrix3      &UT_Matrix3::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_Matrix3      &UT_Matrix3::operator+=(const UT_Vector3 &vec)
{
    fpreal x = vec.x(); fpreal y = vec.y(); fpreal 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_Matrix3      &UT_Matrix3::operator-=(const UT_Vector3 &vec)
{
    fpreal x = vec.x(); fpreal y = vec.y(); fpreal 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;
}

// Scaled outer product update.
// this += b * x * yT
inline
void             UT_Matrix3::outerproductUpdate(fpreal b, 
                                               const UT_Vector3 &v1,
                                               const UT_Vector3 &v2)
{
    fpreal bv1;
    bv1 = b * v1.x();
    matx[0][0]+=bv1*v2.x();
    matx[0][1]+=bv1*v2.y();
    matx[0][2]+=bv1*v2.z();
    bv1 = b * v1.y();
    matx[1][0]+=bv1*v2.x();
    matx[1][1]+=bv1*v2.y();
    matx[1][2]+=bv1*v2.z();
    bv1 = b * v1.z();
    matx[2][0]+=bv1*v2.x();
    matx[2][1]+=bv1*v2.y();
    matx[2][2]+=bv1*v2.z();
}

inline
UT_Vector3       UT_Matrix3::operator[](unsigned row) const
{
    UT_ASSERT_P(row < 3);
    return UT_Vector3(matx[row]);
}


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

inline
UT_Matrix3      operator*(fpreal sc, const UT_Matrix3 &m1)
{
    return m1*sc;
}

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

#endif

---