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>
#include "UT_Assert.h"
#include "UT_Math.h"

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

#include "UT_VectorTypes.h"
class UT_IStream;
class UT_JSONWriter;
class UT_JSONValue;
class UT_JSONParser;
// 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.

template <typename T, typename S>
inline  UT_Vector3T<T>  operator*(const UT_Vector3T<T> &v, const UT_Matrix3T<S>  &m);
template <typename T, typename S>
inline  UT_Vector3T<T>  operator*(const UT_Vector3T<T> &v, const UT_Matrix4T<S>  &m);
template <typename T>
inline  UT_Vector3T<T>  operator+(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2);
template <typename T>
inline  UT_Vector3T<T>  operator-(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2);
template <typename T, typename S>
inline  UT_Vector3T<T>  operator+(const UT_Vector3T<T> &v, S scalar);
template <typename T, typename S>
inline  UT_Vector3T<T>  operator-(const UT_Vector3T<T> &v, S scalar);
template <typename T, typename S>
inline  UT_Vector3T<T>  operator*(const UT_Vector3T<T> &v, S scalar);
template <typename T, typename S>
inline  UT_Vector3T<T>  operator/(const UT_Vector3T<T> &v, S scalar);
template <typename T, typename S>
inline  UT_Vector3T<T>  operator+(S scalar, const UT_Vector3T<T> &v);
template <typename T, typename S>
inline  UT_Vector3T<T>  operator-(S scalar, const UT_Vector3T<T> &v);
template <typename T, typename S>
inline  UT_Vector3T<T>  operator*(S scalar, const UT_Vector3T<T> &v);
template <typename T, typename S>
inline  UT_Vector3T<T>  operator/(S scalar, const UT_Vector3T<T> &v);

/// The dot and cross products between two vectors (see operator*() too)
// @{
template <typename T>
inline  T               dot      (const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2);
template <typename T>
inline  UT_Vector3T<T>  cross    (const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2);
// @}

/// The orthogonal projection of a vector u onto a vector v
template <typename T>
inline  UT_Vector3T<T>  project  (const UT_Vector3T<T> &u, const UT_Vector3T<T> &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)

template <typename T, typename S>
inline UT_Vector3T<T>   rowVecMult(const UT_Vector3T<T> &v, const UT_Matrix3T<S> &m);
template <typename T, typename S>
inline UT_Vector3T<T>   colVecMult(const UT_Matrix3T<S> &m, const UT_Vector3T<T> &v);
template <typename T, typename S>
inline UT_Vector3T<T>   rowVecMult(const UT_Vector3T<T> &v, const UT_Matrix4T<S> &m);
template <typename T, typename S>
inline UT_Vector3T<T>   colVecMult(const UT_Matrix4T<S> &m, const UT_Vector3T<T> &v);

template <typename T, typename S>
inline UT_Vector3T<T>   rowVecMult3(const UT_Vector3T<T> &v, const UT_Matrix4T<S> &m);
template <typename T, typename S>
inline UT_Vector3T<T>   colVecMult3(const UT_Matrix4T<S> &m, const UT_Vector3T<T> &v);
// @}

// TODO: make UT_Vector4 versions of these:

/// Compute the distance between two points
template <typename T>
inline  T               distance3d(const UT_Vector3T<T> &p1, const UT_Vector3T<T> &p2);
/// Compute the distance squared
template <typename T>
inline  T               distance2(const UT_Vector3T<T> &p1, const UT_Vector3T<T> &p2);
template <typename T>
inline  T               segmentPointDist2( const UT_Vector3T<T> &pos, 
                                           const UT_Vector3T<T> &pt1, 
                                           const UT_Vector3T<T> &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.
template <typename T>
UT_API double           intersectLines(const UT_Vector3T<T> &p1,
                                       const UT_Vector3T<T> &v1,
                                       const UT_Vector3T<T> &p2,
                                       const UT_Vector3T<T> &v2,
                                       T &t1, T &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.
template <typename T>
UT_API bool             intersectSegments(const UT_Vector3T<T> &p0, 
                                          const UT_Vector3T<T> &p1,
                                          const UT_Vector3T<T> &a, 
                                          const UT_Vector3T<T> &b, T &t);

/// Returns the U coordinates of the closest points on each of the two 
/// parallel line segments
template <typename T>
UT_API UT_Vector2T<T>   segmentClosestParallel(const UT_Vector3T<T> &p0, 
                                         const UT_Vector3T<T> &p1,
                                         const UT_Vector3T<T> &a, 
                                         const UT_Vector3T<T> &b);
/// Returns the U coordinates of the closest points on each of the two 
/// line segments
template <typename T>
UT_API UT_Vector2T<T>   segmentClosest(const UT_Vector3T<T> &p0, 
                                         const UT_Vector3T<T> &p1,
                                         const UT_Vector3T<T> &a, 
                                         const UT_Vector3T<T> &b);

/// 3D Vector class.
template <typename T>
class UT_API UT_Vector3T
{
public:

    typedef T           value_type;
    static const int    tuple_size = 3;

    /// Default constructor.
    /// No data is initialized! Use it for extra speed.
    inline       UT_Vector3T(void)
                 {
                 }
    inline       UT_Vector3T(T vx, T vy, T vz)
                 {
                    vec[0] = vx; vec[1] = vy; vec[2] = vz;
                 }
    inline       UT_Vector3T(const fpreal16 v[tuple_size])
                 {
                    vec[0] = v[0]; vec[1] = v[1]; vec[2] = v[2];
                 }
    inline       UT_Vector3T(const fpreal32 v[tuple_size])
                 {
                    vec[0] = v[0]; vec[1] = v[1]; vec[2] = v[2];
                 }
    inline       UT_Vector3T(const fpreal64 v[tuple_size])
                 {
                    vec[0] = v[0]; vec[1] = v[1]; vec[2] = v[2];
                 }
    // TODO: make this explicit.
    inline       UT_Vector3T(const UT_Vector4T<T> &v);
    // Default copy constructor is fine.
    template <typename S> 
    inline UT_Vector3T(const UT_Vector3T<S> &v)
    { vec[0] = v.x(); vec[1] = v.y(); vec[2] = v.z(); }

    // Default assignment operator is fine.
    template <typename S> 
    inline UT_Vector3T<T> &operator=(const UT_Vector3T<S> &v)
    { vec[0] = v.x(); vec[1] = v.y(); vec[2] = v.z(); return *this; }

    /// 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_Vector3T<T>      &operator=(const UT_Vector4T<T> &v);

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

    int         equalZero(T 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_Vector3T<T> &vect, T 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));
                }

    bool        isNan() const
                { return SYSisNan(vec[0]) || SYSisNan(vec[1]) || SYSisNan(vec[2]); }

    void        clampZero(T 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_Vector3T<T> &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.
    // @{
    template <typename S>
    inline UT_Vector3T<T>       &operator*=(const UT_Matrix3T<S> &m)
                 { rowVecMult(m); return *this; }
    template <typename S>
    inline UT_Vector3T<T>       &operator*=(const UT_Matrix4T<S> &m)
                 { rowVecMult(m); return *this; }

    template <typename S>
    inline void  multiply3(const UT_Matrix4T<S> &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).
    // @{
    template <typename S>
    inline void  multiplyT(const UT_Matrix3T<S> &mat)
                 { colVecMult(mat); }
    template <typename S>
    inline void  multiply3T(const UT_Matrix4T<S> &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".
    // @{
    template <typename S>
    inline void  multiply3(UT_Vector3T<T> &dest, const UT_Matrix4T<S> &mat) const
                 {
                    dest = ::rowVecMult3(*this, mat);
                 }
    template <typename S>
    inline void  multiplyT(UT_Vector3T<T> &dest, const UT_Matrix3T<S> &mat) const
                 {
                    dest = ::colVecMult(mat, *this);
                 }
    template <typename S>
    inline void  multiply3T(UT_Vector3T<T> &dest, const UT_Matrix4T<S> &mat) const
                 {
                    dest = ::colVecMult3(mat, *this);
                 }
    template <typename S>
    inline void  multiply(UT_Vector3T<T> &dest, const UT_Matrix4T<S> &mat) const
                 {
                    dest = ::rowVecMult(*this, mat);
                 }
    template <typename S>
    inline void  multiply(UT_Vector3T<T> &dest, const UT_Matrix3T<S> &mat) const
                 {
                    dest = ::rowVecMult(*this, mat);
                 }
    // @}

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

    inline void  cross(const UT_Vector3T<T> &v)
                 {
                    operator=(::cross(*this, v));
                 }
    inline
    T    dot(const UT_Vector3T<T> &v) const
                 {
                    return ::dot(*this, v);
                 }
    T    normalize (void)
                 {
                    T dn   = vec[0]*vec[0] + vec[1]*vec[1] +
                                 vec[2]*vec[2];
                    if (dn > std::numeric_limits<T>::min() && dn != 1.0F)
                    {
                        dn = SYSsqrt(dn);
                        (*this) /= dn;
                    }

                    return dn;
                 }

    void         normal(const UT_Vector3T<T> &va, const UT_Vector3T<T> &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_Vector4T<T> &va, const UT_Vector4T<T> &vb);

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

    /// Find the maximum component
    T    maxComponent() const
                 {
                    return SYSmax(vec[0], vec[1], vec[2]);
                 }
    T    minComponent() const
                 {
                    return SYSmin(vec[0], vec[1], vec[2]);
                 }
    T    avgComponent() const
                 {
                    return SYSavg(vec[0], vec[1], vec[2]);
                 }

    /// These allow you to find out what indices to use for different axes
    // @{
    int          findMinAbsAxis() const
                 {
                    T   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
                 {
                    T   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_Vector3T<T> &X, UT_Vector3T<T> &Y) const
                 {
                         if (SYSabs(x()) < 0.6F) Y = UT_Vector3T<T>(1, 0, 0);
                    else if (SYSabs(z()) < 0.6F) Y = UT_Vector3T<T>(0, 1, 0);
                    else                        Y = UT_Vector3T<T>(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
    T    length(void) const { return SYSsqrt(dot(*this)); }
    /// The vector length squared.
    inline
    T    length2(void) const { return dot(*this); }

    /// Calculates the orthogonal projection of a vector u on the *this vector 
    UT_Vector3T<T>       project(const UT_Vector3T<T> &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.
    template <typename S>
    UT_Matrix3T<S> project(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_Vector3T<T>       projection(const UT_Vector3T<T> &p,
                            const UT_Vector3T<T> &v) const;

    /// Projects this onto the line segement [a,b].  The returned point
    /// will lie between a and b.
    UT_Vector3T<T>       projectOnSegment(const UT_Vector3T<T> &va,
                                  const UT_Vector3T<T> &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_Vector3T<T>       projectOnSegment(const UT_Vector3T<T> &va, const UT_Vector3T<T> &vb,
                                  T &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_Vector3T<T> &p1, const UT_Vector3T<T> &v1,
                              const UT_Vector3T<T> &p2, const UT_Vector3T<T> &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_Vector3T<T> &pa,  const UT_Vector3T<T> &pb,
                                 const UT_Vector3T<T> &p2, const UT_Vector3T<T> &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
    /// std::numeric_limits<T>::max() and true is returned.
    bool        areCollinear(const UT_Vector3T<T> &p0, const UT_Vector3T<T> &p1,
                             T *t = 0, T 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_Vector3T<T>      getBary(const UT_Vector3T<T> &t0, const UT_Vector3T<T> &t1,
                        const UT_Vector3T<T> &t2, bool *degen = NULL) const;


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

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

    /// Compute a hash
    unsigned    hash() const    { return SYSvector_hash(data(), tuple_size); }

    // TODO: eliminate these methods. They're redundant, given good inline
    // constructors.
    /// Set the values of the vector components
    void         assign(T xx = 0.0f, T yy = 0.0f, T zz = 0.0f)
                 {
                     vec[0] = xx; vec[1] = yy; vec[2] = zz;
                 }
    /// Set the values of the vector components
    void         assign(const T *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)
                 {
                     if (vec[2] != 0)
                     {
                         T 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_Vector3T<T> &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_Vector3T<T> &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)
    template <typename S>
    void         getDual(UT_Matrix3T<S> &dual) const;

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

    /// @{
    /// Methods to serialize to a JSON stream.  The vector is stored as an
    /// array of 3 reals.
    bool                save(UT_JSONWriter &w) const;
    bool                save(UT_JSONValue &v) const;
    bool                load(UT_JSONParser &p);
    /// @}

    /// Returns the vector size
    static int          entries()       { return tuple_size; }

    /// The data.
    T    vec[tuple_size];

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

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


template <typename T>
inline UT_Vector3T<T>::UT_Vector3T(const UT_Vector4T<T> &v)
{
    vec[0] = v.x();
    vec[1] = v.y();
    vec[2] = v.z();
}

// Free floating functions:
template <typename T>
inline
UT_Vector3T<T>  operator+(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2)
{
    return UT_Vector3T<T>(v1.x()+v2.x(), v1.y()+v2.y(), v1.z()+v2.z());
}

template <typename T>
inline
UT_Vector3T<T>  operator-(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2)
{
    return UT_Vector3T<T>(v1.x()-v2.x(), v1.y()-v2.y(), v1.z()-v2.z());
}

template <typename T, typename S>
inline
UT_Vector3T<T>  operator+(const UT_Vector3T<T> &v, S scalar)
{
    return UT_Vector3T<T>(v.x()+scalar, v.y()+scalar, v.z()+scalar);
}

template <typename T>
inline
UT_Vector3T<T>  operator*(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2)
{
    return UT_Vector3T<T>(v1.x()*v2.x(), v1.y()*v2.y(), v1.z()*v2.z());
}

template <typename T>
inline
UT_Vector3T<T>  operator/(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2)
{
    return UT_Vector3T<T>(v1.x()/v2.x(), v1.y()/v2.y(), v1.z()/v2.z());
}

template <typename T, typename S>
inline
UT_Vector3T<T>  operator+(S scalar, const UT_Vector3T<T> &v)
{
    return UT_Vector3T<T>(v.x()+scalar, v.y()+scalar, v.z()+scalar);
}

template <typename T, typename S>
inline
UT_Vector3T<T>  operator-(const UT_Vector3T<T> &v, S scalar)
{
    return UT_Vector3T<T>(v.x()-scalar, v.y()-scalar, v.z()-scalar);
}

template <typename T, typename S>
inline
UT_Vector3T<T>  operator-(S scalar, const UT_Vector3T<T> &v)
{
    return UT_Vector3T<T>(scalar-v.x(), scalar-v.y(), scalar-v.z());
}

template <typename T, typename S>
inline
UT_Vector3T<T>  operator*(const UT_Vector3T<T> &v, S scalar)
{
    return UT_Vector3T<T>(v.x()*scalar, v.y()*scalar, v.z()*scalar);
}

template <typename T, typename S>
inline
UT_Vector3T<T>  operator*(S scalar, const UT_Vector3T<T> &v)
{
    return UT_Vector3T<T>(v.x()*scalar, v.y()*scalar, v.z()*scalar);
}

template <typename T, typename S>
inline
UT_Vector3T<T>  operator/(const UT_Vector3T<T> &v, S scalar)
{
    // This has to be T because S may be int for "v = v/2" code
    // For the same reason we must cast the 1
    T           inv = ((T)1) / scalar;
    return UT_Vector3T<T>(v.x()*inv, v.y()*inv, v.z()*inv);
}

template <typename T, typename S>
inline
UT_Vector3T<T>  operator/(S scalar, const UT_Vector3T<T> &v)
{
    return UT_Vector3T<T>(scalar/v.x(), scalar/v.y(), scalar/v.z());
}

template <typename T>
inline
T               dot(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2)
{
    return v1.x()*v2.x() + v1.y()*v2.y() + v1.z()*v2.z();
}

template <typename T>
inline
UT_Vector3T<T>  cross(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2)
{
    // compute the cross product:
    return UT_Vector3T<T>(
        v1.y()*v2.z() - v1.z()*v2.y(),
        v1.z()*v2.x() - v1.x()*v2.z(),
        v1.x()*v2.y() - v1.y()*v2.x()
    );
}

template <typename T>
inline
UT_Vector3T<T>  project(const UT_Vector3T<T> &u, const UT_Vector3T<T> &v)
{
    return dot(u, v) / v.length2() * v;
}

template <typename T, typename S>
inline
UT_Vector3T<T>  rowVecMult(const UT_Vector3T<T> &v, const UT_Matrix3T<S> &m)
{
    return UT_Vector3T<T>(
        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)
    );
}

template <typename T, typename S>
inline
UT_Vector3T<T>  colVecMult(const UT_Matrix3T<S> &m, const UT_Vector3T<T> &v)
{
    return UT_Vector3T<T>(
        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)
    );
}

template <typename T, typename S>
inline UT_Vector3T<T>
rowVecMult(const UT_Vector3T<T> &v, const UT_Matrix4T<S> &m)
{
    return UT_Vector3T<T>(
        v.x()*m(0,0) + v.y()*m(1,0) + v.z()*m(2,0) + m(3,0),
        v.x()*m(0,1) + v.y()*m(1,1) + v.z()*m(2,1) + m(3,1),
        v.x()*m(0,2) + v.y()*m(1,2) + v.z()*m(2,2) + m(3,2)
    );
}

template <typename T, typename S>
inline UT_Vector3T<T>
rowVecMult3(const UT_Vector3T<T> &v, const UT_Matrix4T<S> &m)
{
    return UT_Vector3T<T>(
        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)
    );
}

template <typename T, typename S>
inline UT_Vector3T<T>
colVecMult(const UT_Matrix4T<S> &m, const UT_Vector3T<T> &v)
{
    return UT_Vector3T<T>(
        v.x()*m(0,0) + v.y()*m(0,1) + v.z()*m(0,2) + m(0,3),
        v.x()*m(1,0) + v.y()*m(1,1) + v.z()*m(1,2) + m(1,3),
        v.x()*m(2,0) + v.y()*m(2,1) + v.z()*m(2,2) + m(2,3)
    );
}

template <typename T, typename S>
inline UT_Vector3T<T>
colVecMult3(const UT_Matrix4T<S> &m, const UT_Vector3T<T> &v)
{
    return UT_Vector3T<T>(
        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)
    );
}


template <typename T, typename S>
inline
UT_Vector3T<T>  operator*(const UT_Vector3T<T> &v, const UT_Matrix3T<S> &m)
{
    return rowVecMult(v, m);
}

template <typename T, typename S>
inline
UT_Vector3T<T>  operator*(const UT_Vector3T<T> &v, const UT_Matrix4T<S> &m)
{
    return rowVecMult(v, m);
}

template <typename T>
inline
T       distance3d(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2)
{
    return (v1 - v2).length();
}
template <typename T>
inline
T       distance2(const UT_Vector3T<T> &v1, const UT_Vector3T<T> &v2)
{
    return (v1 - v2).length2();
}

// calculate distance squared of pos to the line segment defined by pt1 to pt2
template <typename T>
inline
T       segmentPointDist2(const UT_Vector3T<T> &pos, 
                          const UT_Vector3T<T> &pt1, const UT_Vector3T<T> &pt2 )
{
    UT_Vector3T<T>      vec;
    T   proj_t;
    T   veclen2;

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

    if( proj_t <= (T)0.0 )
    {
        // in bottom cap region, calculate distance from pt1
        vec = pos - pt1;
    }
    else if( proj_t >= veclen2 )
    {
        // in top cap region, calculate distance from pt2
        vec = pos - pt2;
    }
    else
    {
        // middle region, calculate distance from projected pt
        proj_t /= veclen2;
        vec = (pt1 + (proj_t * vec)) - pos;
    }

    return dot(vec, vec);
}

#endif

---