UT/UT_Vector3.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
 *
 * NAME:        Utility Library (C++)
 *
 * COMMENTS:
 *      This class handles fpreal vectors of dimension 3.
 *
 * WARNING:
 *      This class should NOT contain any virtual methods, nor should it
 *      define more member data. The size of UT_Vector3 must always be 
 *      12 bytes (3 floats).
 */

#ifndef __UT_Vector3_h__
#define __UT_Vector3_h__

#include "UT_API.h"
#include <math.h>
#include <limits.h>
#include "UT_Assert.h"
#include "UT_Math.h"

#include <iostream.h>
#include <vector>

class UT_XformOrder;
class UT_Matrix3;
class UT_Matrix4;
class UT_DMatrix3;
class UT_DMatrix4;
class UT_Vector2;
class UT_Vector3;
class UT_Vector4;
class UT_IStream;
// Free floating functions:

// Right-multiply operators (M*v) have been removed. They had previously
// been defined to return v*M, which was too counterintuitive. Once HDK
// etc. users have a chance to update their code (post 7.0) we could
// reintroduce a right-multiply operator that does a colVecMult.

inline  UT_Vector3      operator*(const UT_Vector3 &v, const UT_Matrix3  &m);
inline  UT_Vector3      operator*(const UT_Vector3 &v, const UT_DMatrix3 &m);
inline  UT_Vector3      operator*(const UT_Vector3 &v, const UT_Matrix4  &m);
inline  UT_Vector3      operator*(const UT_Vector3 &v, const UT_DMatrix4 &m);
inline  UT_Vector3      operator+(const UT_Vector3 &v1, const UT_Vector3 &v2);
inline  UT_Vector3      operator-(const UT_Vector3 &v1, const UT_Vector3 &v2);
inline  UT_Vector3      operator+(const UT_Vector3 &v, fpreal scalar);
inline  UT_Vector3      operator-(const UT_Vector3 &v, fpreal scalar);
inline  UT_Vector3      operator*(const UT_Vector3 &v, fpreal scalar);
inline  UT_Vector3      operator/(const UT_Vector3 &v, fpreal scalar);
inline  UT_Vector3      operator+(fpreal scalar, const UT_Vector3 &v);
inline  UT_Vector3      operator-(fpreal scalar, const UT_Vector3 &v);
inline  UT_Vector3      operator*(fpreal scalar, const UT_Vector3 &v);
inline  UT_Vector3      operator/(fpreal scalar, const UT_Vector3 &v);

/// The dot and cross products between two vectors (see operator*() too)
// @{
inline  fpreal          dot      (const UT_Vector3 &v1, const UT_Vector3 &v2);
inline  UT_Vector3      cross    (const UT_Vector3 &v1, const UT_Vector3 &v2);
// @}

/// The orthogonal projection of a vector u onto a vector v
inline  UT_Vector3      project  (const UT_Vector3 &u, const UT_Vector3 &v);

/// Multiplication of a row or column vector by a matrix (ie. right vs. left
/// multiplication respectively). The operator*() declared above is an alias 
/// for rowVecMult(). The functions that take a 4x4 matrix first extend
/// the vector to 4D (with a trailing 1.0 element).
//
// @{
// Notes on optimisation of matrix/vector multiplies:
// - multiply(dest, mat) functions have been deprecated in favour of
// rowVecMult/colVecMult routines, which produce temporaries. For these to
// offer comparable performance, the compiler has to optimize away the
// temporary, but most modern compilers can do this. Performance tests with
// gcc3.3 indicate that this is a realistic expectation for modern
// compilers.
// - since matrix/vector multiplies cannot be done without temporary data,
// the "primary" functions are the global matrix/vector
// rowVecMult/colVecMult routines, rather than the member functions.
// - inlining is explicitly requested only for non-deprecated functions
// involving the native types (UT_Vector3 and UT_Matrix3)

inline UT_Vector3       rowVecMult(const UT_Vector3 &v, const UT_Matrix3 &m);
inline UT_Vector3       colVecMult(const UT_Matrix3 &m, const UT_Vector3 &v);
UT_API  UT_Vector3      rowVecMult(const UT_Vector3 &v, const UT_DMatrix3 &m);
UT_API  UT_Vector3      rowVecMult(const UT_Vector3 &v, const UT_Matrix4 &m);
UT_API  UT_Vector3      rowVecMult(const UT_Vector3 &v, const UT_DMatrix4 &m);
UT_API  UT_Vector3      colVecMult(const UT_DMatrix3 &m, const UT_Vector3 &v);
UT_API  UT_Vector3      colVecMult(const UT_Matrix4 &m, const UT_Vector3 &v);
UT_API  UT_Vector3      colVecMult(const UT_DMatrix4 &m, const UT_Vector3 &v);

UT_API  UT_Vector3      rowVecMult3(const UT_Vector3 &v, const UT_Matrix4 &m);
UT_API  UT_Vector3      rowVecMult3(const UT_Vector3 &v, const UT_DMatrix4 &m);
UT_API  UT_Vector3      colVecMult3(const UT_Matrix4 &m, const UT_Vector3 &v);
UT_API  UT_Vector3      colVecMult3(const UT_DMatrix4 &m, const UT_Vector3 &v);
// @}

// TODO: make UT_Vector4 versions of these:

/// Compute the distance between two points
inline  fpreal          distance3d(const UT_Vector3 &p1, const UT_Vector3 &p2);
/// Compute the distance squared
inline  fpreal          distance2(const UT_Vector3 &p1, const UT_Vector3 &p2);
inline  fpreal          segmentPointDist2( const UT_Vector3 &pos, 
                                           const UT_Vector3 &pt1, 
                                           const UT_Vector3 &pt2 );

/// Intersect the lines p1 + v1 * t1 and p2 + v2 * t2.
/// t1 and t2 are set so that the lines intersect when
/// projected to the plane defined by the two lines.
/// This function returns a value which indicates how close
/// to being parallel the lines are. Closer to zero means
/// more parallel. This is done so that the user of this
/// function can decide what epsilon they want to use.
UT_API double           intersectLines(const UT_Vector3 &p1,
                                       const UT_Vector3 &v1,
                                       const UT_Vector3 &p2,
                                       const UT_Vector3 &v2,
                                       float &t1, float &t2);

/// Returns true if the segments from p0 to p1 and from a to b intersect, and
/// t will contain the parametric value of the intersection on the segment a-b.
/// Otherwise returns false.  Parallel segments will return false.  T is
/// guaranteed to be between 0.0 and 1.0 if this function returns true.
UT_API bool             intersectSegments(const UT_Vector3 &p0, 
                                          const UT_Vector3 &p1,
                                          const UT_Vector3 &a, 
                                          const UT_Vector3 &b, float &t);

/// Returns the U coordinates of the closest points on each of the two 
/// parallel line segments
UT_API UT_Vector2       segmentClosestParallel(const UT_Vector3 &p0, 
                                         const UT_Vector3 &p1,
                                         const UT_Vector3 &a, 
                                         const UT_Vector3 &b);
/// Returns the U coordinates of the closest points on each of the two 
/// line segments
UT_API UT_Vector2       segmentClosest(const UT_Vector3 &p0, 
                                         const UT_Vector3 &p1,
                                         const UT_Vector3 &a, 
                                         const UT_Vector3 &b);

/// 3D Vector class.
class UT_API UT_Vector3 {
public:
    /// Default constructor.
    /// No data is initialized! Use it for extra speed.
    inline       UT_Vector3(void)
                 {
                    UT_ASSERT_COMPILETIME(sizeof(UT_Vector3) ==
                                          3 * sizeof(float));
                 }
    inline       UT_Vector3(fpreal vx, fpreal vy, fpreal vz)
                 {
                    vec[0] = vx; vec[1] = vy; vec[2] = vz;
                 }
                 UT_Vector3(const fpreal32 v[3])
                 {
                    vec[0] = v[0]; vec[1] = v[1]; vec[2] = v[2];
                 }
                 UT_Vector3(const fpreal64 v[3])
                 {
                    vec[0] = v[0]; vec[1] = v[1]; vec[2] = v[2];
                 }
    // TODO: make this explicit.
    inline       UT_Vector3(const UT_Vector4 &v);
    // Default copy constructor is fine.
    //           UT_Vector3(const UT_Vector3 &v);
                ~UT_Vector3(void)               {}

    // Default assignment operator is fine.
    //UT_Vector3&operator=(const UT_Vector3 &v);

    /// Assignment operator that truncates a V4 to a V3. 
    // TODO: remove this. This should require an explicit UT_Vector3()
    // construction, since it's unsafe.
    UT_Vector3  &operator=(const UT_Vector4 &v);

    UT_Vector3  operator-() const
                {
                    return UT_Vector3(-vec[0], -vec[1], -vec[2]);
                }
    inline
    UT_Vector3  &operator+=(const UT_Vector3 &v)
                {
                    vec[0] += v.vec[0];
                    vec[1] += v.vec[1];
                    vec[2] += v.vec[2];
                    return *this;
                }
    inline
    UT_Vector3  &operator-=(const UT_Vector3 &v)
                {
                    vec[0] -= v.vec[0];
                    vec[1] -= v.vec[1];
                    vec[2] -= v.vec[2];
                    return *this;
                }
    unsigned    operator==(const UT_Vector3 &v) const
                {
                    return (vec[0] == v.vec[0] && vec[1] == v.vec[1] &&
                                                  vec[2] == v.vec[2] );
                    
                }
    unsigned    operator!=(const UT_Vector3 &v) const { return !(*this == v); } 

    int         equalZero(fpreal tol = 0.00001f) const
                {
                    return (vec[0] >= -tol && vec[0] <= tol) &&
                           (vec[1] >= -tol && vec[1] <= tol) &&
                           (vec[2] >= -tol && vec[2] <= tol);
                }

    int         isEqual(const UT_Vector3 &vect, fpreal tol = 0.00001f) const
                {
                    return ((vec[0] >= vect.vec[0]-tol) &&
                            (vec[0] <= vect.vec[0]+tol) &&
                            (vec[1] >= vect.vec[1]-tol) &&
                            (vec[1] <= vect.vec[1]+tol) &&
                            (vec[2] >= vect.vec[2]-tol) &&
                            (vec[2] <= vect.vec[2]+tol));
                }

    void        clampZero(fpreal tol = 0.00001f) 
                {
                    if (vec[0] >= -tol && vec[0] <= tol) vec[0] = 0;
                    if (vec[1] >= -tol && vec[1] <= tol) vec[1] = 0;
                    if (vec[2] >= -tol && vec[2] <= tol) vec[2] = 0;
                }


    void         negate() { operator=(operator-()); }


    void         multiplyComponents(const UT_Vector3 &v)
                 {
                     vec[0] *= v.vec[0];
                     vec[1] *= v.vec[1];
                     vec[2] *= v.vec[2];
                 }

    /// If you need a multiplication operator that left multiplies the vector
    /// by a matrix (M * v), use the following colVecMult() functions. If
    /// you'd rather not use operator*=() for right-multiplications (v * M),
    /// use the following rowVecMult() functions. The methods that take a 4x4
    /// matrix first extend this vector to 4D by adding an element equal to 1.0.
    // @{
    inline void  rowVecMult(const UT_Matrix3 &m)
                 { 
                    operator=(::rowVecMult(*this, m));
                 }
    inline void  rowVecMult(const UT_Matrix4 &m)
                 { 
                    operator=(::rowVecMult(*this, m));
                 }
    inline void  rowVecMult(const UT_DMatrix3 &m)
                 { 
                    operator=(::rowVecMult(*this, m));
                 }
    inline void  rowVecMult(const UT_DMatrix4 &m)
                 { 
                    operator=(::rowVecMult(*this, m));
                 }
    inline void  colVecMult(const UT_Matrix3 &m)
                 { 
                    operator=(::colVecMult(m, *this));
                 }
    inline void  colVecMult(const UT_Matrix4 &m)
                 { 
                    operator=(::colVecMult(m, *this));
                 }
    inline void  colVecMult(const UT_DMatrix3 &m)
                 { 
                    operator=(::colVecMult(m, *this));
                 }
    inline void  colVecMult(const UT_DMatrix4 &m)
                 { 
                    operator=(::colVecMult(m, *this));
                 }
    // @}


    /// This multiply will not extend the vector by adding a fourth element.
    /// Instead, it converts the Matrix4 to a Matrix3.  This means that
    /// the translate component of the matrix is not applied to the vector
    // @{
    inline void  rowVecMult3(const UT_Matrix4 &m)
                 { 
                    operator=(::rowVecMult3(*this, m));
                 }
    inline void  rowVecMult3(const UT_DMatrix4 &m)
                 { 
                    operator=(::rowVecMult3(*this, m));
                 }
    inline void  colVecMult3(const UT_Matrix4 &m)
                 { 
                    operator=(::colVecMult3(m, *this));
                 }
    inline void  colVecMult3(const UT_DMatrix4 &m)
                 { 
                    operator=(::colVecMult3(m, *this));
                 }
    // @}



    /// The *=, multiply, multiply3 and multiplyT routines are provided for
    /// legacy reasons. They all assume that *this is a row vector. Generally,
    /// the rowVecMult and colVecMult methods are preferred, since they're
    /// more explicit about the row vector assumption.
    // @{
    inline
    UT_Vector3  &operator*=(const UT_Matrix3 &m)
                 { rowVecMult(m); return *this; }
    UT_Vector3  &operator*=(const UT_DMatrix3 &m)
                 { rowVecMult(m); return *this; }
    UT_Vector3  &operator*=(const UT_Matrix4 &m)
                 { rowVecMult(m); return *this; }
    UT_Vector3  &operator*=(const UT_DMatrix4 &m)
                 { rowVecMult(m); return *this; }

    inline void  multiply3(const UT_Matrix4 &mat)
                 {
                    rowVecMult3(mat);
                 }
    inline void  multiply3(const UT_DMatrix4 &mat)
                 {
                    rowVecMult3(mat);
                 }
    // @}

    /// This multiply will multiply the (row) vector by the transpose of the
    /// matrix instead of the matrix itself.  This is faster than
    /// transposing the matrix, then multiplying (as well there's potentially
    /// less storage requirements).
    // @{
    inline void  multiplyT(const UT_Matrix3 &mat)
                 {
                    colVecMult(mat);
                 }
    inline void  multiplyT(const UT_DMatrix3 &mat)
                 {
                    colVecMult(mat);
                 }
    inline void  multiply3T(const UT_Matrix4 &mat)
                 {
                    colVecMult3(mat);
                 }
    inline void  multiply3T(const UT_DMatrix4 &mat)
                 {
                    colVecMult3(mat);
                 }
    // @}

    /// The following methods implement multiplies (row) vector by a matrix,
    /// however, the resulting vector is specified by the dest parameter
    /// These operations are safe even if "dest" is the same as "this".
    // @{
    inline void  multiply3(UT_Vector3 &dest, const UT_Matrix4 &mat) const
                 {
                    dest = ::rowVecMult3(*this, mat);
                 }
    inline void  multiplyT(UT_Vector3 &dest, const UT_Matrix3 &mat) const
                 {
                    dest = ::colVecMult(mat, *this);
                 }
    inline void  multiply3T(UT_Vector3 &dest, const UT_Matrix4 &mat) const
                 {
                    dest = ::colVecMult3(mat, *this);
                 }
    inline void  multiply(UT_Vector3 &dest, const UT_Matrix4 &mat) const
                 {
                    dest = ::rowVecMult(*this, mat);
                 }
    inline void  multiply(UT_Vector3 &dest, const UT_Matrix3 &mat) const
                 {
                    dest = ::rowVecMult(*this, mat);
                 }

    inline void  multiply3(UT_Vector3 &dest, const UT_DMatrix4 &mat) const
                 {
                    dest = ::rowVecMult3(*this, mat);
                 }
    inline void  multiplyT(UT_Vector3 &dest, const UT_DMatrix3 &mat) const
                 {
                    dest = ::colVecMult(mat, *this);
                 }
    inline void  multiply3T(UT_Vector3 &dest, const UT_DMatrix4 &mat) const
                 {
                    dest = ::colVecMult3(mat, *this);
                 }
    inline void  multiply(UT_Vector3 &dest, const UT_DMatrix4 &mat) const
                 {
                    dest = ::rowVecMult(*this, mat);
                 }
    inline void  multiply(UT_Vector3 &dest, const UT_DMatrix3 &mat) const
                 {
                    dest = ::rowVecMult(*this, mat);
                 }
    // @}

    UT_Vector3  &operator=(fpreal scalar)
                 {
                    vec[0] = vec[1] = vec[2] = scalar;
                    return *this;
                 }
    UT_Vector3  &operator+=(fpreal scalar)
                 {
                    vec[0] += scalar;
                    vec[1] += scalar;
                    vec[2] += scalar;
                    return *this;
                 }
    UT_Vector3  &operator-=(fpreal scalar)
                 {
                    return operator+=(-scalar);
                 }
    inline
    UT_Vector3  &operator*=(fpreal scalar)
                 {
                    vec[0] *= scalar;
                    vec[1] *= scalar;
                    vec[2] *= scalar;
                    return *this;
                 }
    inline
    UT_Vector3  &operator*=(const UT_Vector3 &v)
                 {
                    vec[0] *= v.vec[0];
                    vec[1] *= v.vec[1];
                    vec[2] *= v.vec[2];
                    return *this;
                 }
    inline
    UT_Vector3  &operator/=(fpreal scalar)
                 {
                    return operator*=( 1.0f/scalar );
                 }
    inline
    UT_Vector3  &operator/=(const UT_Vector3 &v)
                 {
                    vec[0] /= v.vec[0];
                    vec[1] /= v.vec[1];
                    vec[2] /= v.vec[2];
                    return *this;
                 }

    inline void  cross(const UT_Vector3 &v)
                 {
                    operator=(::cross(*this, v));
                 }
    inline
    fpreal       dot(const UT_Vector3 &v) const
                 {
                    return ::dot(*this, v);
                 }
    fpreal       normalize (void)
                 {
                    fpreal dn   = vec[0]*vec[0] + vec[1]*vec[1] +
                                 vec[2]*vec[2];
                    if (dn > FLT_MIN && dn != 1.0F)
                    {
                        dn = SYSsqrt(dn);
                        (*this) /= dn;
                    }

                    return dn;
                 }

    void         normal(const UT_Vector3 &va, const UT_Vector3 &vb)
                 {
                    vec[0] += (va.vec[2]+vb.vec[2])*(vb.vec[1]-va.vec[1]);
                    vec[1] += (va.vec[0]+vb.vec[0])*(vb.vec[2]-va.vec[2]);
                    vec[2] += (va.vec[1]+vb.vec[1])*(vb.vec[0]-va.vec[0]);
                 }

    void         normal(const UT_Vector4 &va, const UT_Vector4 &vb);

    /// Finds an arbitrary perpendicular to v, and sets this to it.
    void         arbitraryPerp(const UT_Vector3 &v);
    /// Makes this orthogonal to the given vector.  If they are colinear,
    /// does an arbitrary perp
    void         makeOrthonormal(const UT_Vector3 &v);

    /// Find the maximum component
    fpreal       maxComponent() const
                 {
                    fpreal      m = vec[0];
                    if (vec[1] > m) m = vec[1];
                    if (vec[2] > m) m = vec[2];
                    return m;
                 }
    fpreal       minComponent() const
                 {
                    fpreal      m = vec[0];
                    if (vec[1] < m) m = vec[1];
                    if (vec[2] < m) m = vec[2];
                    return m;
                 }
    fpreal       avgComponent() const
                 {
                    return (vec[0]+vec[1]+vec[2])*(1.0F/3.0F);
                 }

    /// These allow you to find out what indices to use for different axes
    // @{
    int          findMinAbsAxis() const
                 {
                    fpreal      ax = SYSabs(x()), ay = SYSabs(y());
                    if (ax < ay)
                        return (SYSabs(z()) < ax) ? 2 : 0;
                    else
                        return (SYSabs(z()) < ay) ? 2 : 1;
                 }
    int          findMaxAbsAxis() const
                 {
                    fpreal      ax = SYSabs(x()), ay = SYSabs(y());
                    if (ax >= ay)
                        return (SYSabs(z()) >= ax) ? 2 : 0;
                    else
                        return (SYSabs(z()) >= ay) ? 2 : 1;
                 }
    // @}

    /// Given this vector as the z-axis, get a frame of reference such that the
    /// X and Y vectors are orthonormal to the vector.  This vector should be
    /// normalized.
    void         getFrameOfReference(UT_Vector3 &X, UT_Vector3 &Y) const
                 {
                         if (SYSabs(x()) < 0.6F) Y = UT_Vector3(1, 0, 0);
                    else if (SYSabs(z()) < 0.6F) Y = UT_Vector3(0, 1, 0);
                    else                        Y = UT_Vector3(0, 0, 1);
                    X = ::cross(Y, *this);
                    X.normalize();
                    Y = ::cross(*this, X);
                 }

    /// The vector length (not to be confused with the vector dimension).
    inline
    fpreal       length(void) const { return SYSsqrt(dot(*this)); }
    /// The vector length squared.
    inline
    fpreal       length2(void) const { return dot(*this); }

    /// Calculates the orthogonal projection of a vector u on the *this vector 
    UT_Vector3   project(const UT_Vector3 &u) const;

    /// Create a matrix of projection onto this vector: the matrix transforms
    /// a vector v into its projection on the direction of (*this) vector,
    /// ie. dot(*this, v) * this->normalize();
    /// If we need to be normalized, set norm to a non-zero value.
    UT_Matrix3   project(int norm=1);
    UT_DMatrix3  projectD(int norm=1);

    /// Vector p (representing a point in 3-space) and vector v define
    /// a line. This member returns the projection of "this" onto the
    /// line (the point on the line that is closest to this point).
    UT_Vector3   projection(const UT_Vector3 &p,
                            const UT_Vector3 &v) const;

    /// Projects this onto the line segement [a,b].  The returned point
    /// will lie between a and b.
    UT_Vector3   projectOnSegment(const UT_Vector3 &va,
                                  const UT_Vector3 &vb) const;
    /// Projects this onto the line segment [a, b].  The fpreal t is set
    /// to the parametric position of intersection, a being 0 and b being 1.
    UT_Vector3   projectOnSegment(const UT_Vector3 &va, const UT_Vector3 &vb,
                                  fpreal &t) const;

    /// Create a matrix of symmetry around this vector: the matrix transforms 
    /// a vector v into its symmetry around (*this), ie. two times the 
    /// projection of v onto (*this) minus v.
    /// If we need to be normalized, set norm to a non-zero value.
    UT_Matrix3   symmetry(int norm=1);

    /// This method stores in (*this) the intersection between two 3D lines,
    /// p1+t*v1 and p2+u*v2. If the two lines do not actually intersect, we
    /// shift the 2nd line along the perpendicular on both lines (along the
    /// line of min distance) and return the shifted intersection point; this
    /// point thus lies on the 1st line.
    /// If we find an intersection point (shifted or not) we return 0; if
    /// the two lines are parallel we return -1; and if they intersect 
    /// behind our back we return -2. When we return -2 there still is a
    /// valid intersection point in (*this).
    int         lineIntersect(const UT_Vector3 &p1, const UT_Vector3 &v1,
                              const UT_Vector3 &p2, const UT_Vector3 &v2);

    /// Compute the intersection of vector p2+t*v2 and the line segment between
    /// points pa and pb. If the two lines do not intersect we shift the 
    /// (p2, v2) line along the line of min distance and return the point 
    /// where it intersects the segment. If we find an intersection point 
    /// along the stretch between pa and pb, we return 0. If the lines are
    /// parallel we return -1. If they intersect before pa we return -2, and 
    /// if after pb, we return -3. The intersection point is valid with 
    /// return codes 0,-2,-3.
    int         segLineIntersect(const UT_Vector3 &pa,  const UT_Vector3 &pb,
                                 const UT_Vector3 &p2, const UT_Vector3 &v2);

    /// Determines whether or not the points p0, p1 and "this" are collinear.
    /// If they are t contains the parametric value of where "this" is found
    /// on the segment from p0 to p1 and returns true.  Otherwise returns
    /// false.  If p0 and p1 are equal, t is set to FLT_MAX and true is 
    /// returned.
    bool        areCollinear(const UT_Vector3 &p0, const UT_Vector3 &p1,
                             float *t = 0, fpreal tol = 1e-5) const;

    /// Compute (homogeneous) barycentric co-ordinates of this point
    /// relative to the triangle defined by t0, t1 and t2. (The point is
    /// projected into the triangle's plane.)
    UT_Vector3  getBary(const UT_Vector3 &t0, const UT_Vector3 &t1,
                        const UT_Vector3 &t2, bool *degen = NULL) const;


    /// Compute the signed distance from us to a line.
    fpreal      distance(const UT_Vector3 &p1, const UT_Vector3 &v1) const;
    /// Compute the signed distance between two lines.
    fpreal      distance(const UT_Vector3 &p1, const UT_Vector3 &v1,
                         const UT_Vector3 &p2, const UT_Vector3 &v2) const;

    /// Return the components of the vector. The () operator does NOT check
    /// for the boundary condition.
    // @{
    const float *data(void) const       { return &vec[0]; }
    float       *data(void)             { return &vec[0]; }
    inline float&x(void)                { return vec[0]; } 
    inline float x(void) const          { return vec[0]; } 
    inline float&y(void)                { return vec[1]; } 
    inline float y(void) const          { return vec[1]; } 
    inline float&z(void)                { return vec[2]; } 
    inline float z(void) const          { return vec[2]; } 
    inline float&r(void)                { return vec[0]; } 
    inline float r(void) const          { return vec[0]; } 
    inline float&g(void)                { return vec[1]; } 
    inline float g(void) const          { return vec[1]; } 
    inline float&b(void)                { return vec[2]; } 
    inline float b(void) const          { return vec[2]; } 
    inline float&operator()(unsigned i)
                 {
                     UT_ASSERT_P(i < 3);
                     return vec[i];
                 }
    inline float operator()(unsigned i) const
                 {
                     UT_ASSERT_P(i < 3);
                     return vec[i];
                 }
    inline float&operator[](unsigned i)
                 {
                     UT_ASSERT_P(i < 3);
                     return vec[i];
                 }
    inline float operator[](unsigned i) const
                 {
                     UT_ASSERT_P(i < 3);
                     return vec[i];
                 }
    std::vector<float> asStdVector() const;
    // @}

    // TODO: eliminate these methods. They're redundant, given good inline
    // constructors.
    /// Set the values of the vector components
    void         assign(fpreal xx = 0.0f, fpreal yy = 0.0f, fpreal zz = 0.0f)
                 {
                     vec[0] = xx; vec[1] = yy; vec[2] = zz;
                 }
    /// Set the values of the vector components
    void         assign(const float *v)
                 {
                    vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2];
                 }

    /// Express the point in homogeneous coordinates or vice-versa
    // @{
    void         homogenize(void)
                 {
                     vec[0] *= vec[2];
                     vec[1] *= vec[2];
                 }
    void         dehomogenize(void)
                 {
                     fpreal denom = 1.0f / vec[2];
                     vec[0] *= denom;
                     vec[1] *= denom;
                 }
    // @}

    /// assuming that "this" is a rotation (in radians, of course), the
    /// equivalent set of rotations which are closest to the "base" rotation
    /// are found. The equivalent rotations are the same as the original
    /// rotations +2*n*PI
    void         roundAngles(const UT_Vector3 &base);

    /// conversion between degrees and radians
    // @{
    void         degToRad();
    void         radToDeg();
    // @}

    /// It seems that given any rotation matrix and transform order, 
    /// there are two distinct triples of rotations that will result in
    /// the same overall rotation. This method will find the closest of
    /// the two after finding the closest using the above method.
    void         roundAngles(const UT_Vector3 &b, const UT_XformOrder &o);

    /// Return the dual of the vector
    /// The dual is a matrix which acts like the cross product when
    /// multiplied by other vectors.
    /// The following are equivalent:
    /// a.getDual(A);  c = colVecMult(A, b)
    /// c = cross(a, b)
    void         getDual(UT_Matrix3 &dual) const;
    void         getDual(UT_DMatrix3 &dual) const;

    /// Protected I/O methods
    // @{
    void         save(ostream &os, int binary = 0) const;
    bool         load(UT_IStream &is);
    // @}


    /// The data.
    float vec[3];

private:
    /// I/O friends
    // @{
    friend ostream      &operator<<(ostream &os, const UT_Vector3 &v)
                        {
                            v.save(os);
                            return os;
                        }
    // @}
};

// Required for inline rowVecMult and colVecMult routines
#include "UT_Matrix3.h"
// Required for constructor
#include "UT_Vector4.h"


inline UT_Vector3::UT_Vector3(const UT_Vector4 &v)
{
    vec[0] = v.x();
    vec[1] = v.y();
    vec[2] = v.z();
}

// Free floating functions:
inline
UT_Vector3      operator+(const UT_Vector3 &v1, const UT_Vector3 &v2)
{
    return UT_Vector3(v1.x()+v2.x(), v1.y()+v2.y(), v1.z()+v2.z());
}

inline
UT_Vector3      operator-(const UT_Vector3 &v1, const UT_Vector3 &v2)
{
    return UT_Vector3(v1.x()-v2.x(), v1.y()-v2.y(), v1.z()-v2.z());
}

inline
UT_Vector3      operator+(const UT_Vector3 &v, fpreal scalar)
{
    return UT_Vector3(v.x()+scalar, v.y()+scalar, v.z()+scalar);
}

inline
UT_Vector3      operator*(const UT_Vector3 &v1, const UT_Vector3 &v2)
{
    return UT_Vector3(v1.x()*v2.x(), v1.y()*v2.y(), v1.z()*v2.z());
}

inline
UT_Vector3      operator/(const UT_Vector3 &v1, const UT_Vector3 &v2)
{
    return UT_Vector3(v1.x()/v2.x(), v1.y()/v2.y(), v1.z()/v2.z());
}

inline
UT_Vector3      operator+(fpreal scalar, const UT_Vector3 &v)
{
    return UT_Vector3(v.x()+scalar, v.y()+scalar, v.z()+scalar);
}

inline
UT_Vector3      operator-(const UT_Vector3 &v, fpreal scalar)
{
    return UT_Vector3(v.x()-scalar, v.y()-scalar, v.z()-scalar);
}

inline
UT_Vector3      operator-(fpreal scalar, const UT_Vector3 &v)
{
    return UT_Vector3(scalar-v.x(), scalar-v.y(), scalar-v.z());
}

inline
UT_Vector3      operator*(const UT_Vector3 &v, fpreal scalar)
{
    return UT_Vector3(v.x()*scalar, v.y()*scalar, v.z()*scalar);
}

inline
UT_Vector3      operator*(fpreal scalar, const UT_Vector3 &v)
{
    return UT_Vector3(v.x()*scalar, v.y()*scalar, v.z()*scalar);
}

inline
UT_Vector3      operator/(const UT_Vector3 &v, fpreal scalar)
{
    return UT_Vector3(v.x()/scalar, v.y()/scalar, v.z()/scalar);
}

inline
UT_Vector3      operator/(fpreal scalar, const UT_Vector3 &v)
{
    return UT_Vector3(scalar/v.x(), scalar/v.y(), scalar/v.z());
}

inline
fpreal          dot(const UT_Vector3 &v1, const UT_Vector3 &v2)
{
    return v1.x()*v2.x() + v1.y()*v2.y() + v1.z()*v2.z();
}

inline
UT_Vector3      cross(const UT_Vector3 &v1, const UT_Vector3 &v2)
{
    // compute the cross product:
    return UT_Vector3(
        v1.y()*v2.z() - v1.z()*v2.y(),
        v1.z()*v2.x() - v1.x()*v2.z(),
        v1.x()*v2.y() - v1.y()*v2.x()
    );
}

inline
UT_Vector3      project(const UT_Vector3 &u, const UT_Vector3 &v)
{
    return dot(u, v) / v.length2() * v;
}

inline
UT_Vector3      rowVecMult(const UT_Vector3 &v, const UT_Matrix3 &m)
{
    return UT_Vector3(
        v.x()*m(0,0) + v.y()*m(1,0) + v.z()*m(2,0),
        v.x()*m(0,1) + v.y()*m(1,1) + v.z()*m(2,1),
        v.x()*m(0,2) + v.y()*m(1,2) + v.z()*m(2,2)
    );
}

inline
UT_Vector3      colVecMult(const UT_Matrix3 &m, const UT_Vector3 &v)
{
    return UT_Vector3(
        v.x()*m(0,0) + v.y()*m(0,1) + v.z()*m(0,2),
        v.x()*m(1,0) + v.y()*m(1,1) + v.z()*m(1,2),
        v.x()*m(2,0) + v.y()*m(2,1) + v.z()*m(2,2)
    );
}

inline
UT_Vector3      operator*(const UT_Vector3 &v, const UT_Matrix3 &m)
{
    return rowVecMult(v, m);
}

inline
UT_Vector3      operator*(const UT_Vector3 &v, const UT_DMatrix3 &m)
{
    return rowVecMult(v, m);
}

inline
UT_Vector3      operator*(const UT_Vector3 &v, const UT_Matrix4 &m)
{
    return rowVecMult(v, m);
}

inline
UT_Vector3      operator*(const UT_Vector3 &v, const UT_DMatrix4 &m)
{
    return rowVecMult(v, m);
}

inline
fpreal  distance3d(const UT_Vector3 &v1, const UT_Vector3 &v2)
{
    return (v1 - v2).length();
}
inline
fpreal  distance2(const UT_Vector3 &v1, const UT_Vector3 &v2)
{
    return (v1 - v2).length2();
}

// calculate distance squared of pos to the line segment defined by pt1 to pt2
inline
fpreal  segmentPointDist2(const UT_Vector3 &pos, 
                          const UT_Vector3 &pt1, const UT_Vector3 &pt2 )
{
    UT_Vector3  vec;
    fpreal      proj_t;

    vec = pt2 - pt1;
    proj_t = vec.dot( pos - pt1 ) / vec.length2();

    if( UTisLess( proj_t, (fpreal)0.0 ) )
    {
        // in bottom cap region, calculate distance from pt1
        vec = pos - pt1;
    }
    else if( UTisGreater( proj_t, (fpreal)1.0 ) )
    {
        // in top cap region, calculate distance from pt2
        vec = pos - pt2;
    }
    else
    {
        // middle region, calculate distance from projected pt
        vec = (pt1 + (proj_t * vec)) - pos;
    }

    return dot(vec, vec);
}

#endif

---