/Users/craigcornelius/Projects/SPRING Mac Release 0.2/face.cpp

Go to the documentation of this file.

/* $Id: face.cpp,v 1.69 2006/02/01 20:18:01 craig Exp $ */
#include "util.h"
#include <assert.h>
#include "face.h"
#include "facearray.h"
#include "object.h"
#include <math.h>

#ifdef _WIN32
#include <glut.h>
#else
#include <GLUT/glut.h>
#include <OpenGL/gl.h>
#endif

#include "intersection.h"
#include "rasterfont.h"
#include "point3d.h"

const char* Face::rcsid = "@(#) $Id: face.cpp,v 1.69 2006/02/01 20:18:01 craig Exp $ $Copyright: (c)2001 National Biocomputation Center, Stanford University $";

int Face::debug = 0;

Face::Face()
{
    facearray_p = NULL;
    for (int i=0; i<3; i++) {
        node[i] = NULL;
        edge[i] = NULL;
    }
                
        splitNode[0] = splitNode[1] = NULL;
        fakeSplitNode[0] = fakeSplitNode[1] = NULL;
        splitEdge[0] = splitEdge[1] = NULL;
        num_split_edges = 0;

    marker = 0;
    boundingsphere = NULL;
}

Face::~Face()
{
    facearray_p = NULL;
    for (int i=0; i<3; i++) {
        node[i] = NULL;
        edge[i] = NULL;
    }
}

void Face::init(NodeArray* nodearray_p, EdgeArray* edgearray_p,
                FaceArray* facearray_p_in, int i0, int i1, int i2)  
{ 
    // Face knows where it is in facearray 
    facearray_p = facearray_p_in;
    
    // add nodes to this Face, add this Face to those nodes 
    node[0] = nodearray_p->getNode(i0);
    node[1] = nodearray_p->getNode(i1); 
    node[2] = nodearray_p->getNode(i2); 
    
    // call the real init
    init(nodearray_p, edgearray_p, facearray_p, node[0], node[1], node[2]);
}

void Face::init(NodeArray* nodearray_p, EdgeArray* edgearray_p,
                FaceArray* facearray_p_in, Node* n0, Node* n1, Node* n2)  
{
    // see if there's already an edge between n0 and n1, and if not, make one.
    Edge* e0 = n0->is_adjacent_to(n1);
    if (e0 == NULL)
        e0 = edgearray_p->getEdge(edgearray_p->addEdge(n0, n1));
    
    // ditto for n1 and n2
    Edge* e1 = n1->is_adjacent_to(n2);
    if (e1 == NULL)
        e1 = edgearray_p->getEdge(edgearray_p->addEdge(n1, n2));
    
    // and finally, n2 and n0
    Edge* e2 = n2->is_adjacent_to(n0);
    if (e2 == NULL)
        e2 = edgearray_p->getEdge(edgearray_p->addEdge(n2, n0));
    
    // and call the last init to setup links
    init(nodearray_p, edgearray_p, facearray_p_in, n0, n1, n2, e0, e1, e2);
}

void Face::init(NodeArray* nodearray_p, EdgeArray* edgearray_p,
                FaceArray* facearray_p_in, Node* n0, Node* n1, Node* n2,
                Edge* e0, Edge* e1, Edge* e2)  
{ 
    // set internal vars
    facearray_p = facearray_p_in;
    node[0] = n0; node[1] = n1; node[2] = n2;
    edge[0] = e0; edge[1] = e1; edge[2] = e2;
    
    // crosslink us into our node's and edge's lists
    for (int i=0; i<3; i++) {
        // add us to the node's face list
        node[i]->addFaceLink(this);
        
        // add us to the edge's face list
        edge[i]->addFaceLink(this);
    }
}

void Face::replaceNodeLink(Node* old_node, Node* new_node)
{
    for (int i=0; i<3; i++)
        if (node[i] == old_node) node[i] = new_node;
}

void Face::replaceNodeLink(int which, Node* new_node)
{
    node[which] = new_node;
}


Node *Face::getOppositeNode(Edge *e)
{
    if(hasEdgeLink(e) == 0){
        cerr << "this face and this edge are not related" << endl;
        return NULL;
    }
    for(int i = 0; i < 3; i ++){
        Node *tempN = getNodeLink(i);
        if(tempN != e->getNodeLink(0) && tempN != e->getNodeLink(1))
            return tempN;
    }
    return NULL; //actually, should not be returned!
}



void Face::replaceEdgeLink(Edge* old_edge, Edge* new_edge)
{
    for (int i=0; i<3; i++)
        if (edge[i] == old_edge) edge[i] = new_edge;
}

int Face::hasNodeLink(Node* n)
{
    for (int i=0; i<3; i++)
        if (node[i] == n) return(1);
        return(0);
}

void Face::unlink()
{int i;
    // unlink ourselves from our nodes and edges
    for ( i=0; i<3; i++) {
        // unlink from our nodes
        node[i]->deleteFaceLink(this);
        
        // unlink from edges
        edge[i]->deleteFaceLink(this);
    }
        for ( i=0; i<tetra.NumElements(); i++) {
                tetra[i]->deleteFaceLink(this);
        }
}

int Face::hasEdgeLink(Edge* e)
{
    for (int i=0; i<3; i++)
        if (edge[i] == e) return(1);
        return(0);
}

void Face::setNormal(Point3D normal_in)
{ normal = normal_in; }

void Face::SaveAsMesh(ostream& os)
{
    for (int i=0; i<3; i++)
        os << node[i]->getIndex() << " " ;
    os << endl;
}

void Face::draw(int do_texture) 
{ 
    for (int i=0; i<3; i++)
        node[i]->draw(1, 1, do_texture); // do_color, objectwide_color, do_texture
}

// make a label containing the face num
void Face::drawlabel(float offset)
{
    extern RasterFont* font_p;
    char s[20]; sprintf(s, "%d", getIndex());
    Point3D center = getCenter();
    
    Point3D ll = center + offset*normal;
    
    glBegin(GL_LINES);
    glVertex3d(center.x, center.y, center.z);
    glVertex3d(ll.x, ll.y, ll.z);
    glEnd();
    
    font_p->printString(s, ll.x,ll.y,ll.z);
}

void Face::drawnormal()
{
    Point3D p = getCenter();
    glVertex3d(p.x, p.y, p.z);
    glVertex3d(p.x+normal.x, p.y+normal.y, p.z+normal.z);
}

// want some way to mark a face as special, say colliding
void Face::drawmarker()
{
    if (!marker) return;

    // try drawing outline in purple
    glColor3f(1.0, 0.0, 1.0);
    glBegin(GL_LINE_LOOP);
    for (int i = 0; i < 3; i++)
        getNodeLink(i)->drawover(0.05);
    glEnd();
}

int Face::getIndex()
{ return(facearray_p->getIndex(this)); }

// ZVolume returns the exact volume that lies between the face
// and its projection on the X-Y plane. It may be positive or
// negative depending on the mutual orientations of the Z axis
// and the normal of the face.
double Face::ZVolume()
{
    double Zvol1 = 0, Zvol2 = 0, Zvol3 = 0, Zvolume = 0;
    int minzz = getMinZNode();
    Point3D n = getNormal();
    
    // The point with the minimum Z
    Point3D a = node[minzz]->p;
    // The two other points of the face
    Point3D b = node[(minzz+1)%3]->p;
    Point3D c = node[(minzz+2)%3]->p;
    // Projections of b and c on the XY plane
    Point3D bp, cp;
    bp.x = b.x; bp.y = b.y; bp.z = a.z;
    cp.x = c.x; cp.y = c.y; cp.z = a.z;
    
    Zvol1 = TetraVol(a,b,c,bp);
    Zvol2 = TetraVol(a,bp,cp,c);
    Zvol3 = a.z * getZShadowArea();
    
    Zvolume = Zvol1 + Zvol2 + Zvol3;
    if(n.z < 0)
        Zvolume = -Zvolume;
    
    return Zvolume;
}

// TetraVol return the volume of the tetrahedron formed
// by the four points passed in
double Face::TetraVol(Point3D a, Point3D b,
                      Point3D c, Point3D d)
{
    double vol = 0;
    vol = -a.z*b.y*c.x + a.y*b.z*c.x + a.z*b.x*c.y - a.x*b.z*c.y -
        a.y*b.x*c.z + a.x*b.y*c.z + a.z*b.y*d.x - a.y*b.z*d.x -
        a.z*c.y*d.x + b.z*c.y*d.x + a.y*c.z*d.x - b.y*c.z*d.x -
        a.z*b.x*d.y + a.x*b.z*d.y + a.z*c.x*d.y - b.z*c.x*d.y -
        a.x*c.z*d.y + b.x*c.z*d.y + a.y*b.x*d.z - a.x*b.y*d.z -
        a.y*c.x*d.z + b.y*c.x*d.z + a.x*c.y*d.z - b.x*c.y*d.z;
    vol = vol/6.0;
    return vol;
}

int Face::getMinZNode()
{
    Point3D alpha = node[0]->p, beta = node[1]->p, 
        gamma = node[2]->p;
    if(alpha.z <= beta.z)
    {
        if(alpha.z <= gamma.z) return 0;
        else return 2;
    }
    else if(beta.z <= gamma.z)
    {
        if(beta.z <= alpha.z) return 1;
        else return 0;
    }
    else if(gamma.z <= alpha.z)
    {
        if(gamma.z <= beta.z) return 2;
        else return 1;
    }
    printf("/nPb determining MinZNode, returning 0");
    return 0;
}

// This function takes in two nodes in a given face, and returns the third
Node *Face::getThirdNode(Node *n0, Node *n1)
{
    if ((n0 == node[0] && n1 == node[1]) || (n0 == node[1] && n1 == node[0]))
        return node[2];
    else if ((n0 == node[0] && n1 == node[2]) || (n0 == node[2] && n1 == node[0]))
        return node[1];
    else if ((n0 == node[1] && n1 == node[2]) || (n0 == node[2] && n1 == node[1]))
        return node[0];
    else {
        printf("Input nodes are not both in this face\n");
        return NULL;
    }
}

// returns signed distance (positive on normal side, neg on other)
double Face::PlaneDistanceToPoint(Point3D pt)
{
    if (debug) cerr << "PlaneDistanceToPoint(" << pt << ")\n";
    
    // compute plane equation for plane PQR
    // have the normal, compute the constant
    double plane_constant = normal.Dot(node[0]->p); // note, this is -d
    
    // compute the constant for the given point
    double pt_constant = normal.Dot(pt);

    // return the diff
    return (pt_constant - plane_constant);
}


//find the projection of pt on the face plane 
Point3D Face::Projection(Point3D pt)
{
  double dist = PlaneDistanceToPoint(pt);
  Point3D proj = pt - normal*(dist);
  return(proj);

  //old way, inefficient, determinePlaneEq uses normal, why recompute normal?
  // also, normal is already a unit vector
  /*
    double a, b, c, d;
    Point3D proj, normal;
    
    determinePlaneEq(&a, &b, &c, &d); //this->determinePlaneEq(&a, &b, &c, &d);
    normal = Point3D(a, b, c);
    
    proj = pt - ((normal.Dot(pt) + d) / normal.Squared()) * normal;
    return(proj);
  */
}

// should check this code ??? !!!
// computes the distance from a point to a face
double Face::FaceDistanceToPoint(Point3D pt, int early_exit, Point3D *nearpt)
{
    int i, test[3]; 
    double dist;
    Edge *e1, *e2;
    Point3D proj, normal;
    Node *p[3];
    
    normal = getNormal();
    p[0] = getNodeLink(0);
    p[1] = getNodeLink(1);
    p[2] = getNodeLink(2);
    
    proj = Projection(pt);
    
    //printf("projection = %g %g %g\n", proj.x, proj.y, proj.z);
    
    for(i = 0; i < 3; i ++)
        // tests the position of proj compared to the edge p[i]p[i+1] 
      test[i] = (((p[(i + 1) % 3]->p - p[i]->p).Cross(pt - p[i]->p)).Dot(normal)
                 >= 0) ? 0 : 1;
    
    assert(!(test[0] == 1 && test[1] == 1 && test[2] == 1));
    
    // this means that pt projects _inside_ the face, I hope
    if(test[0] == 0 && test[1] == 0 && test[2] == 0) {
      // want to return positive distance!!
      // otherwise ClosestFaceToPoint, etc. are wrong.
      if (early_exit)
        *nearpt = proj;
      return(fabs(PlaneDistanceToPoint(pt)));
    }

    if (early_exit)
      return(99999.9); // not inside face, want to return large value

    for(i = 0; i < 3; i ++){
        //printf("i = %d\n", i);
        if(test[i] == 1 && test[(i+1) % 3] == 0 && test[(i+2) % 3] == 0){
            //printf("i = %d\n", i);
            e1 = p[i]->findEdgeLink(p[(i+1) % 3]);
            dist = e1->EdgeDistanceToPoint(proj); 
            //printf("dist to edge = %g\n", dist);
            dist = sqrt(dist * dist + (pt - proj).Squared());
            return dist;
        }
    }
    
    for(i = 0; i < 3; i ++){
        //printf("i = %d\n", i);
        if(test[i] == 1 && test[(i+1) % 3] == 1 && test[(i+2) % 3] == 0){
            //printf("i = %d\n", i);
            e1 = p[i]->findEdgeLink(p[(i+1) % 3]);
            e2 = p[(i+1) % 3]->findEdgeLink(p[(i+2) % 3]);
            dist = MIN(e1->EdgeDistanceToPoint(proj),e2->EdgeDistanceToPoint(proj));
            //printf("dist to edge = %g\n", dist);
            dist = sqrt(dist * dist + (pt - proj).Squared());
            return dist;
        }
    }
    
    // if got here and nothing above hit, return large value, print error
    printf("problem in FaceDistanceToPoint\n");
    return(99999.9);
}


int Face::IsColliding(Point3D p)
{
    if (debug) cerr << "Face::IsColliding(" << p << ")\n";
    
    // first, are we close enough to the plane that the face is in?
    const double max_dist = 0.1; // mm
    double dist = fabs(PlaneDistanceToPoint(p));
    if (dist > max_dist) return(0);
    
    if (debug) 
        cerr << "  distance =  " << dist << "(max_dist = " << max_dist << ")\n";
    
    // close enough to plane, are we within the triangle bounds?
    Point3D a = node[0]->p;
    Point3D b = node[1]->p;
    Point3D c = node[2]->p;
    
    // basically, if the point is in the triangle, approximately, as point may
    // be slightly out of plane, each of these cross products
    // will point in the direction of the normal, and the length of the sum
    // will be ~ 2 instead of ~ 0.  Special case if p is collinear with one side,
    // then cross product is 0 vector, so don't normalize, the sum will have
    // length 1 and will pass the test
    
    // do the ab edge
    Point3D ab = b - a, ap = p - a;
    Point3D n1 = ab.Cross(ap);
    if (!(n1.IsZero()))
        n1.Normalize(); 
    if ((n1 + normal).Squared() < 1.0) return(0);
    
    // do the bc edge
    Point3D bc = c - b, bp = p - b;
    Point3D n2 = bc.Cross(bp);
    if (!(n2.IsZero()))
        n2.Normalize(); 
    if ((n2 + normal).Squared() < 1.0) return(0);
    
    // do the ca edge
    Point3D ca = a - c, cp = p - c;
    Point3D n3 = ca.Cross(cp);
    if (!(n3.IsZero()))
        n3.Normalize(); 
    if ((n3 + normal).Squared() < 1.0) return(0);
    
    // debugging
    if (debug) 
        cerr << "  it passed all tests- is a collision!\n";
    
    // yay!!  it worked!
    return(1);
    
#ifdef SLOW_WAY_TO_TEST_SIGNS
    // all-att-once way of doing it
    //Point3D ab = b - a, bc = c - b, ca = a - c;
    //Point3D ap = p - a, bp = p - b, cp = p - c;
    //Point3D n1 = ab.Cross(ap), n2 = bc.Cross(bp), n3 = ca.Cross(cp);
    
    n1.Normalize(); n2.Normalize(); n3.Normalize();
    if ((n1 + normal).Length() < 1.0) return(0);
    if ((n2 + normal).Length() < 1.0) return(0);
    if ((n3 + normal).Length() < 1.0) return(0);
    
#elif defined(BAD_WAY_TO_TEST_SIGNS)
    // NOTE!  This did NOT work- why???
    cerr << n1.Sign() << " " << n2.Sign() << " " << n3.Sign() << "   " 
        << normal.Sign() << endl;
    if (n1.Sign() != normal.Sign()) return(0);
    if (n2.Sign() != normal.Sign()) return(0);
    if (n3.Sign() != normal.Sign()) return(0);
#endif
}

Node* Face::getClosestNodeInFace(Point3D pt)
{
    double dist0 = node[0]->p.SquaredDist(pt);
    double dist1 = node[1]->p.SquaredDist(pt);
    double dist2 = node[2]->p.SquaredDist(pt);
    if (dist0 < dist1) {
        if (dist0 < dist2) return(node[0]);
        else return(node[2]);
    } else {
        if (dist1 < dist2) return(node[1]);
        else return(node[2]);
    }
}

// computes area of the face
// A = |p2-p1 * p3-p1|/2
double Face::getArea()
{
    Point3D p1 = node[0]->p;
    Point3D p2 = node[1]->p;
    Point3D p3 = node[2]->p;
    Point3D first_part = p2 - p1;
    Point3D last_part = p3 - p1;
    Point3D cross_part = first_part.Cross(last_part);
    
    return( cross_part.Length()/2.0 );
}

// computes area of projection of tri onto x-z plane
// for use in volume calculations
double Face::getShadowArea()
{
    Point3D p1 = node[0]->p;
    Point3D p2 = node[1]->p;
    Point3D p3 = node[2]->p;
    p1.y = p2.y = p3.y = 0;     // shadow onto the x-z plane
    Point3D first_part = p2 - p1;
    Point3D last_part = p3 - p1;
    Point3D cross_part = first_part.Cross(last_part);
    
    return( cross_part.Length()/2.0 );
}

double Face::getZShadowArea()
{
    Point3D p1 = node[0]->p;
    Point3D p2 = node[1]->p;
    Point3D p3 = node[2]->p;
    p1.z = p2.z = p3.z = 0;     // shadow onto the x-y plane
    Point3D first_part = p2 - p1;
    Point3D last_part = p3 - p1;
    Point3D cross_part = first_part.Cross(last_part);
    
    return( cross_part.Length()/2.0 );
}

void Face::getAdjacentFaces(Face** f1_p, Face** f2_p, Face** f3_p)
{
    if (!f1_p || !f2_p || !f3_p) return;
    
    { // get face on other side of edge 0
        Face* efa = edge[0]->getFaceLink(0);
        Face* efb = edge[0]->getFaceLink(1);
        if (efa == this) *f1_p = efb;
        else *f1_p = efa;
    }
    
    { // get face on other side of edge 1
        Face* efa = edge[1]->getFaceLink(0);
        Face* efb = edge[1]->getFaceLink(1);
        if (efa == this) *f2_p = efb;
        else *f2_p = efa;
    }
    
    { // get face on other side of edge 2
        Face* efa = edge[2]->getFaceLink(0);
        Face* efb = edge[2]->getFaceLink(1);
        if (efa == this) *f3_p = efb;
        else *f3_p = efa;
    }
}


ostream& operator<<(ostream& os, const Face& f_in)
{
    Face& f = (Face&)f_in;
    if (f.facearray_p)
        os << "Face#" << ((Face&)f).getIndex() << ": ";
    else os << "Face#UNLINKED: ";
    os << "Nodelinks: " ;
    { for (int i=0; i<3; i++) os << f.node[i]->getIndex() << " " ; }
    os << ", Edgelinks: ";
    { for (int i=0; i<3; i++) os << f.edge[i]->getIndex() << " " ; }
        os << ", Tetralinks: ";
  { for (int i=0; i<f.tetra.NumElements(); i++) os << f.tetra[i]->getIndex() << " " ; }

    // print out adjacent faces
    os << ", adjfaces: ";
    for (int i=0;i<3; i++) {
        Edge* cur_edge = ((Face)f).getEdgeLink(i);
        Face* other_face = cur_edge->getOtherFace((Face*)&f);
        if (other_face)
            os << other_face->getIndex() << " ";
        else os << "NULL ";
    }
    
    return(os);
}

int Face::SanityCheck(FaceArray* real_fa_p)
{
    char s[1024]; s[0] = '\0';
    char tmps[1024]; tmps[0] = '\0';
    
    // test usual stuff
    int my_index = -1;
    if (!facearray_p) strcat(s, "facearray is NULL!");
    if (facearray_p != real_fa_p) strcat(s, "facearray is BAD!");
    else {
        my_index = getIndex();
        if (my_index < 0) {
           sprintf(tmps, "getIndex() = %d\n", my_index);
           strcat(s, tmps);
        }
    }

    // is this a problem???
    //if (normal.Length() != 1.0) strcat(s, "normal is not unit-length");

        // test that we have a boundingsphere leaf if our object has a tree
        Object* obj = getFaceArray()->getObject();
        if (obj->getBoundingSphereRoot() && !boundingsphere) 
                strcat(s, "boundingsphere leaf is NULL");
    
    { // check nodelinks
        for (int i=0; i<3; i++) {
            if (!node[i]) { 
                sprintf(tmps, "nodelink[%d] is NULL", i); 
                strcat(s, tmps);
            }
            else if (!node[i]->hasFaceLink(this)) {
                sprintf(tmps, "nodelink[%d] doesn't point back to us!", i); 
                strcat(s, tmps);
            }
        }
    }
    
    { // check edgelinks
        for (int i=0; i<3; i++) {
            if (!edge[i]) {
                sprintf(tmps, "edgelink[%d] is NULL!", i); 
                strcat(s, tmps);
            }
            else if (!edge[i]->hasFaceLink(this)) {
                sprintf(tmps, "edgelink[%d] doesn't point back to us!", i); 
                strcat(s, tmps);
            }
        }
    }

    { // verify my tetra links
    for (int i=0; i<tetra.NumElements(); i++) {
        Tetra* t = tetra[i];
        if (!t) {
            sprintf(tmps, "tetralink[%d] is NULL!", i); 
            strcat(s, tmps);
        }
        else {
            if (!t->hasFaceLink(this)) {
                sprintf(tmps, "tetralink[%d] does not have a link back to us!", i); 
                strcat(s, tmps);

                // look for duplicate facelinks
                for (int j=0; j<i; j++) {
                    if (t == tetra[j]) {
                        sprintf(tmps, "tetralink[%d] is a duplicate of tetralink[%d]!", i, j); 
                        strcat(s, tmps);
                    }
                } // for
           }
        } // else
    } // for
    } // verify tetra links

    // and if there's a problem, print it out
    if (strlen(s) > 1) {
        cerr << "Face[" << my_index << "]: is INSANE:" << s << endl;
        return(0);
    }
    
    // return ok
    return(1);
}

// Determines the equation of the plane which contains the two points S 
// and E, and a vector n which is the mean of the 2 normals to the original 
// surface in S and E
void Face::determinePlaneEq(double *a_p, double *b_p, double *c_p, double *d_p)
{
    Point3D pt = getNodeLink(0)->p;
    Point3D n = getNormal();
    *a_p = n.x;
    *b_p = n.y;
    *c_p = n.z;
    *d_p = - ((*a_p) * pt.x + (*b_p) * pt.y + (*c_p) * pt.z);
}

/* Determines if this face intersects the input edge, and if so
* returns true and determines the intersection point */
int Face::intersectsEdge(Edge *e_in, Point3D *intersection_p)
{
    if (debug) 
        cerr << "Face::intersectsEdge(" << e_in << ", " << intersection_p 
        << ")\n";
    double a, b, c, d;
    int it_crosses;
    
    // find the equation of the plane containing this face
    determinePlaneEq(&a, &b, &c, &d);
    
    // and see if the input edge (line segment) crosses this plane while
    // determining point of intersection
    it_crosses = e_in->crossesPlane(a, b, c, d, intersection_p);
    
    if (!it_crosses) return 0;
    
    if (debug) cerr << "    it crosses plane!\n";
    
    // it does cross, is the intersection point within the face?
    return (IsColliding(*intersection_p));
}

// Determines if this face intersects the input face, and if so
// returns true and determines (one) point of intersection (of many?)
int Face::intersectsFace(Face *f_in, Point3D *intersection_p)
{
    if (debug) 
        cerr << "Face::intersectsFace(" << f_in << ", " << intersection_p 
        << ")\n";
    
#if (1)
    Edge *e;
    { // first check to see if this face's edges (eg the tool's edges) intersect
        // the input face (eg a face of the tube)
        // (we waste time by determining plane equation up to three times ???)
        for (int i = 0; i < 3; i++) {
            e = this->getEdgeLink(i);
            if (f_in->intersectsEdge(e, intersection_p)) return 1;
        }
    }
    
    { // now do vice versa in case input face is piercing this face
        for (int i = 0; i < 3; i++) {
            e = f_in->getEdgeLink(i);
            if (this->intersectsEdge(e, intersection_p)) return 1;
        }
    }
    
    // this face does not pierce/be-pierced-by input face
    return 0;
#else
    // try and do some smarter cutouts
    Point3D thisNormal = this->getNormal();
    double thisPlaneConst = thisNormal.Dot(this->getNodeLink(0)->p);
    // distances from f_in vertices to this face
    double q0d = thisNormal.Dot(f_in->getNodeLink(0)->p) - thisPlaneConst;
    double q1d = thisNormal.Dot(f_in->getNodeLink(1)->p) - thisPlaneConst;
    double q2d = thisNormal.Dot(f_in->getNodeLink(2)->p) - thisPlaneConst;
    
    // no intersection if f_in is completely on one side of this face
    if (((q0d > 0) && (q1d > 0) && (q2d > 0)) ||
        ((q0d < 0) && (q1d < 0) && (q2d < 0)))
        return 0;
    
    Point3D otherNormal = f_in->getNormal();
    double otherPlaneConst = otherNormal.Dot(f_in->getNodeLink(0)->p);
    // distances from this face vertices to f_in
    double p0d = otherNormal.Dot(this->getNodeLink(0)->p) - otherPlaneConst;
    double p1d = otherNormal.Dot(this->getNodeLink(1)->p) - otherPlaneConst;
    double p2d = otherNormal.Dot(this->getNodeLink(2)->p) - otherPlaneConst;
    
    // no intersection if this face is completely on one side of f_in
    if (((p0d > 0) && (p1d > 0) && (p2d > 0)) ||
        ((p0d < 0) && (p1d < 0) && (p2d < 0)))
        return 0;
    
    // now check pairs of vertices of f_in on opposite sides of this face
    if ((q0d >= 0 && q1d < 0) || (q0d < 0 && q1d >= 0)) {
        // vertices 0 and 1 -- interpolate to find plane intersection point
        intersection_p->interp(f_in->getNodeLink(0)->p, f_in->getNodeLink(1)->p,
            fabs(q0d)/(fabs(q0d)+fabs(q1d)));
        // now see if this point is actually inside this face
        if (this->IsColliding(*intersection_p))
            return 1;
    }
    if ((q0d >= 0 && q2d < 0) || (q0d < 0 && q2d >= 0)) {
        // vertices 0 and 2 -- interpolate to find plane intersection point
        intersection_p->interp(f_in->getNodeLink(0)->p, f_in->getNodeLink(2)->p,
            fabs(q0d)/(fabs(q0d)+fabs(q2d)));
        // now see if this point is actually inside this face
        if (this->IsColliding(*intersection_p))
            return 1;
    }
    if ((q1d >= 0 && q2d < 0) || (q1d < 0 && q2d >= 0)) {
        // vertices 1 and 2 -- interpolate to find plane intersection point
        intersection_p->interp(f_in->getNodeLink(1)->p, f_in->getNodeLink(2)->p,
            fabs(q1d)/(fabs(q1d)+fabs(q2d)));
        // now see if this point is actually inside this face
        if (this->IsColliding(*intersection_p))
            return 1;
    }
    
    // now check pairs of vertices of this face on opposite sides of f_in
    if ((p0d >= 0 && p1d < 0) || (p0d < 0 && p1d >= 0)) {
        // vertices 0 and 1 -- interpolate to find plane intersection point
        intersection_p->interp(this->getNodeLink(0)->p, this->getNodeLink(1)->p,
            fabs(p0d)/(fabs(p0d)+fabs(p1d)));
        // now see if this point is actually inside f_in face
        if (f_in->IsColliding(*intersection_p))
            return 1;
    }
    if ((p0d >= 0 && p2d < 0) || (p0d < 0 && p2d >= 0)) {
        // vertices 0 and 2 -- interpolate to find plane intersection point
        intersection_p->interp(this->getNodeLink(0)->p, this->getNodeLink(2)->p,
            fabs(p0d)/(fabs(p0d)+fabs(p2d)));
        // now see if this point is actually inside f_in face
        if (f_in->IsColliding(*intersection_p))
            return 1;
    }
    if ((p1d >= 0 && p2d < 0) || (p1d < 0 && p2d >= 0)) {
        // vertices 1 and 2 -- interpolate to find plane intersection point
        intersection_p->interp(this->getNodeLink(1)->p, this->getNodeLink(2)->p,
            fabs(p1d)/(fabs(p1d)+fabs(p2d)));
        // now see if this point is actually inside f_in face
        if (f_in->IsColliding(*intersection_p))
            return 1;
    }
    
    return 0;
#endif
}

// Determines if this face intersects the input edge (Chris)
bool Face::intersectsEdge(Edge *edge)
{
    double a, b, c, d;
    Point3D intersectionPoint;

    // find the equation of the plane containing this face
    determinePlaneEq(&a, &b, &c, &d);
    
    // and see if the input edge (line segment) crosses this plane while
    // determining point of intersection
    if (edge->crossesPlane(a, b, c, d, &intersectionPoint)) {
                // it does cross, is the intersection point within the face?
                return IsColliding(intersectionPoint);
        }
        else {
                return false;
        }
}

// Determines if this face intersects the input face - Chris
bool Face::intersectsFace(Face *otherFace)
{
        // supposed to be a fast intersection test function!
        Intersection* intersection = Intersection::Instance();
        return intersection->intersects(this, otherFace);
}

// Determines if this face intersects with the input box (axis aligned box) - Chris
bool Face::intersectsBox(Point3D *min, Point3D *max)
{
        // supposed to be a fast intersection test function!
        Intersection* intersection = Intersection::Instance();
        return intersection->intersects(this, min, max);
}

// i move myself
int Face::ResolveMyself(Face *f_in)
{
    if (debug) cerr << "Face[" << getIndex() << "]:ResolveMyself()\n";
    
    // first check if the faces are intersecting
    
    // if they are, then
    // project the input face onto the plane defined by my face
    for (int i=0; i<3; i++) {
        Node* movee = node[i];
        Point3D b = movee->p - f_in->getCenter();
        double vd = f_in->getNormal().Dot(b); 
        //cerr << "vd: " << vd << endl;
        
        // scale up the force to make things happen faster
        //vd *= 4.0;
        
        if (vd < 0.0) {
            movee->setFext(-vd*f_in->getNormal());      
            
            // so at the next update, the node should move
            // but make sure we remove the external force when 
            // it is not applied anymore
            
        }       
    }
    
    // return success
    return(1);
}
 void Face::addTetraLink(Tetra *t)
{
  if (debug)
    cout << "\nFace " << getIndex() << " ::addTetraLink(" << t->getIndex() << ")\n";
  tetra.Append(t);
}

 void Face::deleteTetraLink(Tetra* t)
{
  if (debug)
    cout << "\nEdge " << getIndex() << " ::deleteTetraLink(" << t->getIndex() << ")\n";

  tetra.DeleteItem(t);

}
int Face::stepOffVerts(Point3D intPt) 
{
        for (int i=0; i<3; i++) {
                if (node[i]->p.Dist(intPt) < 0.00001) { 
                        //cout << "intersection on a vert" << endl;
                        intPt = .01*(node[(i+1)%3]->p - node[(i+0)%3]->p);                              
                        return i;
                }
        }
        return -1;
}

int Face::stepOffEdges(Point3D intPt) 
{
        for (int i=0; i<3; i++) {
                if (edge[i]->EdgeDistanceToPoint(intPt) < 0.00001) { 
                 //     cout << "intersection on an edge" << endl;
                        Point3D vect = edge[i]->getNodeLink(1)->p - edge[i]->getNodeLink(0)->p;
                        intPt = .01*(normal.Cross(vect));       
                        return i;                               
                }
        }
        return -1;
}
int Face::getCase(int (&v)[9]) 
{
        //cout << "\nface: " << getIndex() << endl;
        
        int num_edges = 0;      
    {
      for (int i=0; i<3; i++) {
                if (edge[i]->getSplitNode(0) != NULL)
                        num_edges++;
      }
        }
        
        int num_faces = 0;
        if (getSplitNode(0) != NULL) {
                num_faces++;
        }
        
//      cout << "number of edges split = " << num_edges << endl;
//      cout << "number of faces split = " << num_faces << endl;

        Node* n[9];
        for (int i=0; i<9; i++) {
                n[i] = NULL;
                v[i] = -1;
        }
        
        if (num_edges == 1) {
        
                        if (num_faces == 0) {
                                cerr << "entered and exited thru the same edge " << endl;
                                return 2;
                        }
                        if (num_faces == 1) {                                                           
        
                                n[0] = getOppositeNode(splitEdge[0]);
                                int node0Idx = getNodeIndex(n[0]); cerr << "node[0] index: " << node0Idx << endl;
                                n[1] = getNodeLink((node0Idx+1)%3); cerr << "node[0] index: " << (node0Idx+1)%3 << endl;
                                n[2] = getNodeLink((node0Idx+2)%3); cerr << "node[0] index: " << (node0Idx+2)%3 << endl;
                                
                                n[3] = getSplitNode(0);

                                for (int i=0; i<4; i++) v[i] = n[i]->getIndex();

                                n[4] = splitEdge[0]->getSplitNode(0);
                                n[5] = splitEdge[0]->getSplitNode(1);

                                if (getFaceArray()->getNodeArray()->getNode(v[1])->isConnectedTo(n[5])){
                                        v[4] = n[5]->getIndex();
                                        v[5] = n[4]->getIndex();
                                }
                                else {
                                        v[4] = n[4]->getIndex();
                                        v[5] = n[5]->getIndex();                        
                                }
                                return 0;                       
                        }
                }
                if (num_edges == 2) {
                        if (num_faces == 0) {
                        
                        
                                n[0] = splitEdge[1]->getNodeLink(0);
                                n[1] = splitEdge[1]->getNodeLink(1);

                                for (int i=0; i<2; i++) {                       
                                        int numsplit = 0;
                                        for (int j=0; j< n[i]->numEdgeLinks(); j++) {
                                                Edge* myedge = n[i]->getEdgeLink(j);
                                                if (myedge->getSplitNode(0) != NULL && hasEdgeLink(myedge)) {
                                                        numsplit++;
                                                }
                                        }
                                        if (numsplit > 1) {
                                                v[0] = n[i]->getIndex();
                                                v[1] = n[(i+1)%2]->getIndex();

                                                break;
                                        }
                                }
                                n[2] = getOtherNode(n[0],n[1]);
                                v[2] = n[2]->getIndex();

                                
                                if ( getEdgeLinkIndex(splitEdge[1]) == (getEdgeLinkIndex(splitEdge[0])+1)%3 ) {
                                        n[3] = splitEdge[0]->getSplitNode(0);
                                        n[4] = splitEdge[0]->getSplitNode(1);

                                        n[5] = splitEdge[1]->getSplitNode(0);
                                        n[6] = splitEdge[1]->getSplitNode(1);
                                }
                                
                                else {
                                        n[3] = splitEdge[1]->getSplitNode(0);
                                        n[4] = splitEdge[1]->getSplitNode(1);

                                        n[5] = splitEdge[0]->getSplitNode(0);
                                        n[6] = splitEdge[0]->getSplitNode(1);           
                                        
                                        int temp = v[1];
                                        v[1] = v[2];
                                        v[2] = temp;                            
                                }

                                if (getFaceArray()->getNodeArray()->getNode(v[2])->isConnectedTo(n[4])){
                                        v[3] = n[4]->getIndex();
                                        v[4] = n[3]->getIndex();
                                }
                                else {
                                        v[3] = n[3]->getIndex();
                                        v[4] = n[4]->getIndex();                        
                                }

                                if (getFaceArray()->getNodeArray()->getNode(v[1])->isConnectedTo(n[6])){
                                        v[5] = n[6]->getIndex();
                                        v[6] = n[5]->getIndex();
                                }
                                else {
                                        v[5] = n[5]->getIndex();
                                        v[6] = n[6]->getIndex();                        
                                }
                                        
                                if (getFakeSplitNode(0) == NULL) {
                                        Point3D mid; mid.interp(n[3]->p,n[5]->p,.5);
                                        n[7] = getFaceArray()->getNodeArray()->getNode(getFaceArray()->getNodeArray()->addNode(mid));
                                        n[8] = getFaceArray()->getNodeArray()->getNode(getFaceArray()->getNodeArray()->addNode(mid));
                                        setFakeSplitNodes(n[7],n[8]);
                //                      getFaceArray()->getObject()->addToCutPath(n[7]);
                //                      getFaceArray()->getObject()->addToCutPath(n[8]);
                                }
                                n[7] = getFakeSplitNode(0);
                                n[8] = getFakeSplitNode(1);                             
                                
                                if (getFaceArray()->getNodeArray()->getNode(v[0])->isConnectedTo(n[7])) {
                                        v[7] = n[8]->getIndex();
                                        v[8] = n[7]->getIndex();
                                }
                                else {
                                        v[7] = n[7]->getIndex();
                                        v[8] = n[8]->getIndex();
                                }

                                return 1;
                        }
                }       
                
                
                // CBNOTE if you get here you have an invalid case so remove space nodes
                
        //cout << "                     has an invalid case" << endl;
        return -1;
}
Node* Face::getOtherNode(Node* n0, Node* n1)
{
        if (!hasNodeLink(n0) || !hasNodeLink(n1)) return NULL;
        
        for (int i=0; i<3; i++) {       
                if (node[i] == n0 || node[i] == n1) continue;
                return  node[i];
        }
        return NULL;
}
int Face::intersectsPlaneSweptByCuttingEdge(Edge* e_in, Point3D* intersection_p)
{

        if (intersectsEdge(e_in,intersection_p)) return 1;      


        return 0;
}

int Face::vectorCrossesMe(Point3D P1, Point3D P2, Point3D* intPt) 
{
        double a, b, c, d;
  determinePlaneEq(&a, &b, &c, &d);

        double alpha, t;

  alpha = a * (P1.x - P2.x) + b * (P1.y - P2.y) + c * (P1.z - P2.z);
  if(alpha == 0) return 0;

  t = -(a * P2.x + b * P2.y + c * P2.z + d) / alpha;
  if ((t < 0) || (t > 1)) return 0;

  intPt->x = (P1.x - P2.x) * t + P2.x;
  intPt->y = (P1.y - P2.y) * t + P2.y;
  intPt->z = (P1.z - P2.z) * t + P2.z;
        
        if (intPt->insideTriangle(node[0]->p, node[1]->p, node[2]->p)) return 1;
        return 0;

}
void Face::interpolateTextureCoords(Node* n0)
{
        double d[3]; 
        double sum_dist = 0;
        
        Point3D nt[3];
        
        n0->setTextureCoords(Point3D(0,0,0));
        
        for (int i=0; i<3; i++) {
                d[i] = node[i]->p.Dist(n0->p);
                sum_dist += d[i];
                nt[i] = node[i]->getTextureCoords();
        }

        Point3D text_coords;

        text_coords.x = (d[0]*nt[0].x + d[1]*nt[1].x + d[2]*nt[2].x)/sum_dist;
        text_coords.y = (d[0]*nt[0].y + d[1]*nt[1].y + d[2]*nt[2].y)/sum_dist;
        text_coords.z = (d[0]*nt[0].z + d[1]*nt[1].z + d[2]*nt[2].z)/sum_dist;
        
        n0->setTextureCoords(text_coords);
                        
}
int Face::getEdgeLinkIndex(Edge* e)
{
  for (int i=0; i<3; i++)
    if (edge[i] == e) return(i);
  return(-1);
}
Edge* Face::getEdge(int n0, int n1)
{
  for (int i=0; i<3; i++){
    if ( (edge[i]->getNodeLink(0)->getIndex() == n0 && edge[i]->getNodeLink(1)->getIndex() == n1) || 
                                 (edge[i]->getNodeLink(1)->getIndex() == n0 && edge[i]->getNodeLink(0)->getIndex() == n1) )
                                 return edge[i];
        }
  return(NULL);
}


//finds if faces share and edge and returns 1 if not then returns 0
Edge* Face::isEdgeAdjacent(Face *f)
{       
        for(int x=0; x<3; x++)
        for(int z=0; z<3; z++) {
        if(this->getEdgeLink(x) == f->getEdgeLink(z))
                return f->getEdgeLink(z);
        }       
        return NULL;
}

//find other point shared by two faces, meaning not p
Point3D Face::getSatelite(Face *f, Node *in)
{
        Point3D         otherp;
        Edge                    *sharededge;
        Node                    *nearnode;

        sharededge = this->isEdgeAdjacent(f);
        if(sharededge != NULL)
                otherp = sharededge->getOtherNode(in)->p;
        else //we have a gap or non-adjoining faces around a node
        {
                nearnode        = this->getClosestNodeInFace(f->getCenter(), in);
                otherp          = nearnode->p;
                //cout << "Non-adjacent faces at node ";
        }

        return otherp;
}

//anil return closest one except notme
Node* Face::getClosestNodeInFace(Point3D pt, Node* notme)
{
  double dist0 = node[0]->p.SquaredDist(pt);
  double dist1 = node[1]->p.SquaredDist(pt);
  double dist2 = node[2]->p.SquaredDist(pt);
  double sum;
        int who = 0;

        //round about way to eliminat notme to reuse old code
        sum = dist0+dist1+dist2;
        for(int i=0; i<3; i++){
                if(node[i] == notme) who = i;
        }
        if(who==0) dist0 = sum;
        if(who==1) dist1 = sum;
        if(who==2) dist2 = sum;
        
        if (dist0 < dist1) {
     if (dist0 < dist2) return(node[0]);
     else return(node[2]);
  } else {
     if (dist1 < dist2) return(node[1]);
     else return(node[2]);
  }
}

---