/Users/craigcornelius/Projects/SPRING Mac Release 0.2/point3d.h

Go to the documentation of this file.

// $Id: point3d.h,v 1.33 2006/05/24 16:52:41 sean Exp $
// $Copyright: (c)2001 National Biocomputation Center, Stanford University $

#ifndef Point3DBase_Header
#define Point3DBase_Header

#include <iostream>
#include <math.h>

using namespace std;

//using namespace std;

template <class T> class Point3DBase {
    
    
public:
    
    T x,y,z;
    Point3DBase() {x = y = z = 0;}
    Point3DBase(T xin, T yin, T zin) {
        x = xin; y = yin; z = zin; }
    Point3DBase(const Point3DBase<T> &pt) { 
        x = pt.x; y = pt.y; z = pt.z; }
    
    friend Point3DBase<T> operator+(const Point3DBase<T> &pt1, 
        const Point3DBase<T> &pt2) { 
        return (Point3DBase(pt1.x + pt2.x, pt1.y + pt2.y, pt1.z + pt2.z)); }
    
        /*
        friend Point3DBase<T> operator+(const Point3DBase<T> &pt, const int i) {
        return (Point3DBase<T>(pt.x + i, pt.y + i, pt.z + i)); }
        
          friend Point3DBase<T> operator+(const int i, const Point3DBase<T> &pt) {
          return (Point3DBase<T>(pt.x + i, pt.y + i, pt.z + i)); }
    */
    
    friend Point3DBase<T> operator-(const Point3DBase<T> &pt1, 
        const Point3DBase<T> &pt2) {
        return (Point3DBase<T>(pt1.x - pt2.x, pt1.y - pt2.y, pt1.z - pt2.z)); }
    
        /*
        friend Point3DBase<T> operator-(const Point3DBase<T> &pt, const int i) {
        return (Point3DBase<T>(pt.x - i, pt.y - i, pt.z - i)); }
        
          friend Point3DBase<T> operator-(const int i, const Point3DBase<T> &pt) {
          return (Point3DBase<T>(pt.x - i, pt.y - i, pt.z - i)); }
    */
    
    friend Point3DBase<T> operator-(const Point3DBase<T> &pt) {
        return (Point3DBase<T>(-pt.x, -pt.y, -pt.z)); }
    
    friend Point3DBase<T> operator*(const Point3DBase<T> &pt, const int i) {
        return (Point3DBase<T>(pt.x * i, pt.y * i, pt.z * i)); }
    
    friend Point3DBase<T> operator*(const int i, const Point3DBase<T> &pt) {
        return (Point3DBase<T>(pt.x * i, pt.y * i, pt.z * i)); }
    
    friend Point3DBase<T> operator*(const Point3DBase<T> &pt, const double d) {
        return (Point3DBase<T>((T)(pt.x * d), (T)(pt.y * d), (T)(pt.z * d))); }
    
    friend Point3DBase<T> operator*(const double d, const Point3DBase<T> &pt) {
        return (Point3DBase<T>((T)(pt.x * d), (T)(pt.y * d), (T)(pt.z * d))); }
    
    friend Point3DBase<T> operator*(const Point3DBase<T> &pt1, 
        const Point3DBase<T> &pt2) {
        return (Point3DBase<T>((T)(pt1.x * pt2.x), 
            (T)(pt1.y * pt2.y), 
            (T)(pt1.z * pt2.z))); }
    
    friend Point3DBase<T> operator/(const Point3DBase<T> &pt, const int i) {
        return (Point3DBase<T>(pt.x/i, pt.y/i, pt.z/i)); }
    
    friend Point3DBase<T> operator/(const Point3DBase<T> &pt, const double d) {
        return (Point3DBase<T>((T)(pt.x/d), (T)(pt.y/d), (T)(pt.z/d))); }
    
    friend Point3DBase<T> operator/(const Point3DBase<T> &pt1, 
        const Point3DBase<T> &pt2) {
        return (Point3DBase<T>((T)(pt1.x/pt2.x), 
            (T)(pt1.y/pt2.y), 
            (T)(pt1.z/pt2.z))); }
    
    friend ostream& operator<<(ostream &os, const Point3DBase<T> &pt) {
        return (os << pt.x << "," << pt.y << "," << pt.z); }
    
    friend istream& operator>>(istream &is, Point3DBase<T> &pt) {
        is >> pt.x; is.ignore(); is >> pt.y; is.ignore(); is >> pt.z;
        return(is); }
    
    Point3DBase<T> operator=(const Point3DBase<T> &pt) 
    { x = pt.x; y = pt.y; z = pt.z; return(*this); }
    
    Point3DBase<T> operator-=(const Point3DBase<T> &pt) 
    { x -= pt.x; y -= pt.y; z -= pt.z; return(*this); }
    
    Point3DBase<T> operator+=(const Point3DBase<T> &pt) 
    { x += pt.x; y += pt.y; z += pt.z; return(*this); }
    
    Point3DBase<T> operator*=(const int i) 
    { x *= i; y *= i; z *= i; return(*this); }
    
    Point3DBase<T> operator*=(const double d) 
    { x = (T)(double(x)*d); y = (T)(double(y)*d); z = (T)(double(z)*d); 
    return(*this); }
    
    Point3DBase<T> operator*=(const Point3DBase<T> &pt) 
    { x *= pt.x; y *= pt.y; z *= pt.z; return(*this); }
    
    Point3DBase<T> operator/=(const int i) 
    { x /= i; y /= i; z /= i; return(*this); }
    
    Point3DBase<T> operator/=(const double d) 
    { x = (T)(double(x)/d); y = (T)(double(y)/d); z = (T)(double(z)/d);
    return(*this); }
    
    Point3DBase<T> operator/=(const Point3DBase<T> &pt) 
    { x /= pt.x; y /= pt.y; z /= pt.z; return(*this); }
    
    friend int operator==(const Point3DBase<T> &pt1, const Point3DBase<T> &pt2) 
    { return((pt1.x == pt2.x) && (pt1.y == pt2.y) && (pt1.z == pt2.z)); }
    
    friend int operator!=(const Point3DBase<T> &pt1, const Point3DBase<T> &pt2) 
    { return((pt1.x != pt2.x) || (pt1.y != pt2.y) || (pt1.z != pt2.z)); }
    
    friend int operator<(const Point3DBase<T> &pt1, const Point3DBase<T> &pt2) 
    { return ((pt1.x < pt2.x) && (pt1.y < pt2.y) && (pt1.z < pt2.z)); }
    
    friend int operator>(const Point3DBase<T> &pt1, const Point3DBase<T> &pt2) 
    { return ((pt1.x > pt2.x) && (pt1.y > pt2.y) && (pt1.z > pt2.z)); }
    
    double Squared()
    { return (x*x + y*y + z*z); }
    
    double Length()
    { return (sqrt(x*x + y*y + z*z)); }
    
    Point3DBase<T> Abs()
    { return (Point3DBase<T>( fabs(x), fabs(y), fabs(z) )); }
    
    int IsZero()
    { return ( (x==0.0) && (y==0.0) && (z==0.0) ); }
    
    void Normalize()
    { 
        double len = Length();
        if (len == 0.0) {
            //printf("point3d.h::Normalize div by 0, set to 1,0,0 for now\n");
            x = T(1); y = T(0); z = T(0);
            return;
        }
        x = T(x/len); y = T(y/len); z = T(z/len);
    }
    
    Point3DBase<T> Sign()
    { return( Point3DBase<T>((x >= 0), (y >= 0), (z >= 0)) ); }
    

        // returns the dot-product of two vectors
    double Dot(const Point3DBase<T> &pt)
    {return (x*pt.x + y*pt.y + z*pt.z);}
    

        // returns the cross-product of two vectors
    Point3DBase<T> Cross(const Point3DBase<T> &pt)
    {return (Point3DBase<T>((T)(y*pt.z - z*pt.y), 
    (T)(z*pt.x - x*pt.z), 
    (T)(x*pt.y - y*pt.x))); }
    
    double Determinant(const Point3DBase<T> &pt1, const Point3DBase<T> &pt2)
    {return (x*(pt1.y*pt2.z - pt1.z*pt2.y) - y*(pt1.x*pt2.z - pt1.z*pt2.x) +
    z*(pt1.x*pt2.y - pt1.y*pt2.x)); }
    
    // BE AWARE OF OVER/UNDERFLOW!  
    // (Type should probably be double to prevent over/underflow on the subtract)
    double SquaredDist(Point3DBase<T> &pt) { 
        Point3DBase<T> diff = *this - pt;
        return (diff.Squared());
    }
    
    void interp(Point3DBase<T> u, Point3DBase<T> v, double lambda) {            
        x = u.x + lambda * (v.x - u.x);
        y = u.y + lambda * (v.y - u.y);
        z = u.z + lambda * (v.z - u.z);
    }
    

        // returns the distance
    double Dist(Point3DBase<T> &pt) { 

                return( sqrt(SquaredDist(pt)) );

        }
    

        // returns the angle in the range [0, PI] between this and the input vector

        double Angle(Point3DBase<T> &pt) {

                double dotProduct = Dot(pt);

                return acos(dotProduct);

        }


    // bounds values by min
    void ForceMin(double min) 
    { x = MAX(min, x); y = MAX(min, y); z = MAX(min, z); }
    
    // bounds values by max
    void ForceMax(double max) 
    { x = MAX(max, x); y = MAX(max, y); z = MAX(max, z); }
    

        // rotates a point given the inputted rotation matrix
    void Xform(double M[4][4]) {
        Point3DBase<T> newpos;
        newpos.x = x*M[0][0] + y*M[1][0] + z*M[2][0];
        newpos.y = x*M[0][1] + y*M[1][1] + z*M[2][1];
        newpos.z = x*M[0][2] + y*M[1][2] + z*M[2][2];
        x = newpos.x; y = newpos.y; z = newpos.z;
    }
                int insideQuadrilateral(const Point3DBase<T> p0,const Point3DBase<T> p1, 
                                                                                                const Point3DBase<T> p2,const Point3DBase<T> p3) 
{
#define EPSILON  0.00001
#define TWOPI 6.283185307179586476925287

        Point3DBase<T> verts[4];
        verts[0] = p0;
        verts[1] = p1;
        verts[2] = p2;
        verts[3] = p3; 
        
        Point3DBase<T> p4,p5;
  
  double m1,m2;
  double anglesum=0,costheta;

        for (int i=0;i<4;i++) {

        p4.x = verts[i].x - x;
        p4.y = verts[i].y - y;
        p4.z = verts[i].z - z;
        p5.x = verts[(i+1)%4].x - x;
        p5.y = verts[(i+1)%4].y - y;
        p5.z = verts[(i+1)%4].z - z;

        m1 = p4.Length();
        m2 = p5.Length();
        if (m1*m2 <= EPSILON)
         return(1); /* We are on a node, consider this inside */
        else
         costheta = (p4.x*p5.x + p4.y*p5.y + p4.z*p5.z) / (m1*m2);

        anglesum += acos(costheta);
        }
        
        double diff = anglesum - TWOPI;
        if (sqrt(diff*diff) <= EPSILON) return (1);
        return (0);
}

int insideTriangle(const Point3DBase<T> p0, const Point3DBase<T> p1,
                                                                         const Point3DBase<T> p2) 
{
#define EPSILON  0.00001
#define TWOPI 6.283185307179586476925287

        Point3DBase<T> verts[3];
        verts[0] = p0;
        verts[1] = p1;
        verts[2] = p2; 
        
        Point3DBase<T> p4,p5;
  
  double m1,m2;
  double anglesum=0.0;
        double costheta;

        for (int i=0;i<3;i++) {

        p4.x = verts[i].x - x;
        p4.y = verts[i].y - y;
        p4.z = verts[i].z - z;
        p5.x = verts[(i+1)%3].x - x;
        p5.y = verts[(i+1)%3].y - y;
        p5.z = verts[(i+1)%3].z - z;

        m1 = p4.Length();
        m2 = p5.Length();
        if (m1*m2 <= EPSILON)
         return(1); /* We are on a node, consider this inside */
        else
         costheta = (p4.x*p5.x + p4.y*p5.y + p4.z*p5.z) / (m1*m2);

        anglesum += acos(costheta);
        }
        double diff = anglesum - TWOPI;
        if (sqrt(diff*diff) <= EPSILON) return (1);
        return (0);
}
    
};


// instantiations
typedef Point3DBase<double>     Point3D;
typedef Point3DBase<int>        iPoint3D;

#endif

---