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

Go to the documentation of this file.

// $Id: facearray.cpp,v 1.72 2006/05/25 22:23:47 craig Exp $
#include "face.h"
#include "facearray.h"
#include "object.h"
#include "objectarray.h"
#include <stdio.h>
#include <stdlib.h>
#ifdef _WIN32
#include <glut.h>
#else
#include <GLUT/glut.h>
#include <OpenGL/gl.h>
#endif

#include <string>

using namespace std;

const char* FaceArray::rcsid = "@(#) $Id: facearray.cpp,v 1.72 2006/05/25 22:23:47 craig Exp $ $Copyright: (c)2001 National Biocomputation Center, Stanford University $";
int FaceArray::debug = 0;

FaceArray::FaceArray(NodeArray* nodearray_p_in, EdgeArray* edgearray_p_in) 
: nodearray_p(nodearray_p_in), edgearray_p(edgearray_p_in)
{
    faces = NULL;
    numfaces = maxfaces = 0;
}

FaceArray::~FaceArray() 
{
    // go through and free memory
    reset();
    
    // and set some noticable values
    numfaces = maxfaces = -1;
}

void FaceArray::reset(void) 
{
    if (faces) {
        // go through and delete each face
        for (int i=0; i<numfaces; i++)
            delete faces[i];
        
        // free up the array itself
        delete [] faces;
    }
    
    // and reset values
    numfaces = maxfaces = 0;
}

// the same face, or if the node is part of the edge, it returns NULL.
Face* FaceArray::getFace(Edge *e, Node *n)
{
    // check that n is not one of the node which defines e!!!
    if( (n == e->getNodeLink(0)) || (n == e->getNodeLink(1))) {
        cerr << "The node cannot be one of the edge link!!" << endl;
        return NULL;
    }
    
    for(int i = 0; i < e->numFaceLinks(); i ++){
        Face *curFace = e->getFaceLink(i);
        for(int j = 0; j < 3; j ++)
            if(n == curFace->getNodeLink(j))
                return curFace;
    }
    
    cerr << "The node is not 'adjacent' to the edge!" << endl;
    return NULL; //in case the node is not "adjacent" to the edge...
}

int FaceArray::getIndex(Face* f)
{
    // loop through, looking for it
    for (int i=0; i<numfaces; i++)
        if (faces[i] == f)
            return(i);
        
        // didn't find it
        return(-1);
}

void FaceArray::allocateFaces(int maxfaces_in) 
{ 
    maxfaces = maxfaces_in;
    faces = new Face*[maxfaces]; 
    numfaces=0; 
    if (!faces) printf("Could not allocate space for %d faces\n", maxfaces);
}

int FaceArray::addFace()
{
    faces[numfaces] = new Face;
    if (numfaces >= maxfaces) {
        cerr << "FaceArray::addFace()- blew bounds of array (" << maxfaces 
            << ")!\n";
        return(numfaces);
    }
    
    // return index of this guy and bump counter up
    return(numfaces++); 
}


int FaceArray::addFace(int i1, int i2, int i3)
{ 
    int index = addFace();
    faces[index]->init(nodearray_p, edgearray_p, this, i1,i2,i3); 
    return(index);
}


/* function to open the mesh data file and extract the data into globals
* faces[] and numfaces */
int FaceArray::getMeshFaceData(FILE *meshfile)
{
    char s[MAXLINE];
    int numfaces_given, numedges_given;
    
    // get the number of faces from the first line
    if (myfgets(s, MAXLINE, meshfile) == NULL) {
        printf("No lines in face data file\n");
        return(0);
    }
    if (sscanf(s, "%d", &numfaces_given) != 1) {
        printf("Could not read numfaces\n");
        return(0);
    }
    printf("%d faces\n", numfaces_given);
    
    // allocate space for Face objects + 100% for new faces
    allocateFaces(numfaces_given * 2);
    
    // allocate space for Edge objects
    numedges_given = 3 * numfaces_given;    // usually 2xnumfaces = numedges
    printf("%d maxedges\n", numedges_given);
    edgearray_p->allocateEdges(numedges_given);
    
    // loop through and get data for each face
    for (int i=0; i<numfaces_given; i++) {
        if (myfgets(s, MAXLINE, meshfile) == NULL) {
            printf("Not enough faces in data file\n");
            return(0);
        }
        
        { // initialize face with current data line
            int i1, i2, i3;
            int numread = sscanf(s, "%d %d %d", &i1, &i2, &i3);
            
            if (numread != 3) {
                printf("Bad read on face %d- '%d %d %d' - IGNORED\n", i, i1,i2,i3);
                continue;
            }
            
#ifdef DO_READ_SANITY_CHECKS
            // good read: set face info (this will also take care of edges!!)
            int new_index = addFace(i1,i2,i3);
            
            // do a sanity check to make sure
            Face* new_face = getFace(new_index);
            int ok = new_face->SanityCheck(this);
            
            // if there was a problem (adding > 2 faces to a edge, etc), delete
            if (!ok) {
                new_face->unlink();
                DeleteFace(new_face);
            }
#else
            // good read: set face info (this will also take care of edges!!)
            addFace(i1,i2,i3);
#endif
        }
        
    }
    
    // return success
    return(1);
}

/* function to open the mesh data file and extract the data into globals
* faces[] and numfaces */
int FaceArray::getSMFFaceData(FILE *meshfile)
{
    char s[MAXLINE];
    int i, numfaces_given, numedges_given;
    
    // read file once to get the number of verts and polys
    numfaces_given = 0;
    while (!feof(meshfile)) {
        char s[256];
        fgets(s, 256, meshfile);
        if (feof(meshfile)) break;
        switch (s[0]) {
        case 'f' : numfaces_given++; break;
        case 'e' : //edge
        case 'v' : // vertex 
        case '#' : // comment
        default  : // who knows?
            break;
        }
    }
    
    // reset input file pointer
    fseek(meshfile, 0, SEEK_SET);
    
    printf("%d faces\n", numfaces_given);
    
    // allocate space for Face objects + 50% for new faces
    allocateFaces((int)(numfaces_given * 1.5));
    
    // allocate space for Edge objects
    if(edgearray_p->getNumEdges()==0)
    {
        numedges_given = 3 * numfaces_given;    // usually 2xnumfaces = numedges
        printf("%d maxedges\n", numedges_given);
        edgearray_p->allocateEdges(numedges_given);
    }
    
    // loop through and get data for each face
    for (i = 0; i < numfaces_given; i++) {
        do {
            if (myfgets(s, MAXLINE, meshfile) == NULL) {
                printf("Not enough faces in data file\n");
                return(0);
            }
        } while (s[0] != 'f');
        
        { // initialize face with current data line
            int i1, i2, i3;
            int numread = sscanf(s, "%*c %d %d %d", &i1, &i2, &i3);
            
            if (debug) cerr << "Face #" << i << ": " << i1 << " " 
                << i2 << " " << i3 << endl;
            
            if (numread != 3) {
                printf("Bad read on face %d- '%d %d %d' - IGNORED\n", i, i1,i2,i3);
                continue;
            }
            
            // SMF is 1-based, we are 0-index-based
            i1--; i2--; i3--;
            
#ifdef DO_READ_SANITY_CHECKS
            // good read: set face info (this will also take care of edges!!)
            int new_index = addFace(i1,i2,i3);
            
            // do a sanity check to make sure
            Face* new_face = getFace(new_index);
            int ok = new_face->SanityCheck(this);
            
            // if there was a problem (adding > 2 faces to a edge, etc), delete
            if (!ok) {
                new_face->unlink();
                DeleteFace(new_face);
            }
#else
            // good read: set face info (this will also take care of edges!!)
            addFace(i1,i2,i3);
#endif
        }
        
    }
    
    // return success
    return(1);
}

int FaceArray::getOBJFaceData(FILE *meshfile)
{
    char s[MAXLINE];
    int i, numfaces_given, numedges_given;
    
    // read file once to get the number of verts and polys
    numfaces_given = 0;
    while (!feof(meshfile)) {
        char s[256];
        fgets(s, 256, meshfile);
        if (feof(meshfile)) break;
        switch (s[0]) {
        case 'f' : numfaces_given++; break;
        case 'g' : // group
        case 'v' : // vertex 
        case '#' : // comment
        default  : // who knows?
            break;
        }
    }
    
    // reset input file pointer
    fseek(meshfile, 0, SEEK_SET);
    
    printf("%d faces\n", numfaces_given);
    
    // allocate space for Face objects + 50% for new faces
    // (could have rects given, so is 2x num face lines
    allocateFaces((int)(numfaces_given * 2.5));
    
    // allocate space for Edge objects
    if(edgearray_p->getNumEdges()==0)
    {
        numedges_given = 6 * numfaces_given;    // usually 2xnumfaces = numedges
        printf("%d maxedges\n", numedges_given);
        edgearray_p->allocateEdges(numedges_given);
    }
    
    // loop through and get data for each face
    for (i = 0; i < numfaces_given; i++) {
        do {
            if (myfgets(s, MAXLINE, meshfile) == NULL) {
                printf("Not enough faces in data file\n");
                return(0);
            }
        } while (s[0] != 'f');
        
        { // initialize face with current data line
            int i1, i2, i3, i4;
            char arg1[20], arg2[20], arg3[20], arg4[20];        
            // currently only handle triangles
            int numread = sscanf(s, "%*s %s %s %s %s", arg1, arg2, arg3, arg4);
            i1 = atoi(arg1); i2 = atoi(arg2); i3 = atoi(arg3); i4 = atoi(arg4);
            
            if (debug) {
                cerr << "Face #" << i << ": " << i1 << " " 
                    << i2 << " " << i3 << endl;
                if (numread > 3) 
                    cerr << "Face #" << i << ": " << i2 << " " 
                    << i3 << " " << i4 << endl;
            }
            
            // OBJ is 1-based, we are 0-index-based
            i1--; i2--; i3--; i4--;
            
            // good read: set face info (this will also take care of edges!!)
            addFace(i1,i2,i3);
            if (numread > 3) addFace(i1,i3,i4);
        }
        
    }
    
    // return success
    return(1);
}

void FaceArray::SaveMeshFaceData(ostream& os)
{
    os << numfaces << endl;
    for (int i=0; i<numfaces; i++)
        faces[i]->SaveAsMesh(os);
}

void FaceArray::SaveSMFFaceData(ostream& os)
{
    nodearray_p->CacheIndices();
    for (int i=0; i<numfaces; i++) {
        Face* f = getFace(i);
        os << "f";
        for (int j=0; j<3; j++) {
            Node* n = f->getNodeLink(j);
            os << " " << n->getCachedIndex() + 1;       // SMF is 1-index based
        }
        os << endl;
    }
}

void FaceArray::SaveVRMLFaceData(ostream& os)
{
    nodearray_p->CacheIndices();
    for (int i=0; i<numfaces; i++) {
        Face* f = getFace(i);
        for (int j=0; j<3; j++) {
            Node* n = f->getNodeLink(j);
            os << "    " << n->getCachedIndex() << ", ";
        }
        os << "-1,\n";
    }
}

void FaceArray::SaveSTLFaceData(ostream& os)
{
    for (int i=0; i<numfaces; i++) {
        Face* f = getFace(i);
        
        // give it the normal
        Point3D normal = f->getNormal();
        os << "FACET NORMAL " << normal << "\n";
        
        // give it the vertex info
        os << "OUTER LOOP\n";
        for (int j=0; j<3; j++) {
            Node* n = f->getNodeLink(j);
            os << "VERTEX " << n->p.x << " " << n->p.y << " " << n->p.z << "\n";
        }
        os << "ENDLOOP\n";
        os << "ENDFACET\n";
    }
}

/* Function for drawing faces */
void FaceArray::DrawFaces(int do_texture)
{
    // find out if we're in select (picking) mode.  We do this because
    // the glLoadName slows things down quite a bit- a new token per triangle!
    // (and even if we're just doing a callobj()!)
    // We should only do this if in picking mode...
    
    GLint rendermode;
    glGetIntegerv(GL_RENDER_MODE, &rendermode);
    if (debug) cerr << "DrawFaces(): rendermode = " << rendermode << endl;
    
    if (rendermode == GL_RENDER) {
        { // draw the triangles
            glBegin(GL_TRIANGLES);
            for (int i = 0; i < numfaces; i++)
                faces[i]->draw(do_texture);
            glEnd();
        }
    }
    else {      // assume GL_SELECT
        for (int i = 0; i < numfaces; i++) {
            // must be before glBegin for it to pick it up
            glPushName(i / 65535);      // push on most significant word
            glPushName(i % 65535);      // push on least significant word
            glBegin(GL_TRIANGLES);
            
            // do the draw (no texture needed)
            faces[i]->draw(0);
            
            // pop `em off
            glEnd();
            glPopName();        // pop off least significant word
            glPopName();        // pop off most significant word
        }
    }
}

void FaceArray::drawlabels(float offset)
{
    for (int i=0; i<numfaces; i++)
        faces[i]->drawlabel(offset);
}

void FaceArray::drawmarkers()
{
    for (int i=0; i<numfaces; i++)
        faces[i]->drawmarker();
}

void FaceArray::drawnormals()
{
    glBegin(GL_LINES);
    glColor3f(1.0, 0.5, 1.0); // draw the normals in light purple
    for (int i=0; i<numfaces; i++)
        faces[i]->drawnormal();
    glEnd();
}

// function to compute the normals of each face 
void FaceArray::computeNormals(void)
{
    for (int i = 0; i < numfaces; i++)
        faces[i]->computenormal();
}

void FaceArray::computeIsobarys(void)
{
    for (int i = 0; i < numfaces; i++)
        faces[i]->computeisobary();
}

void FaceArray::DeleteFace(int index)
{
    if (debug) cerr << "FaceArray::DeleteFace(" << index << ")\n";
    
    // sanity check
    if (index < 0) return;
    
    // delete the one we've got
    delete faces[index];
    
    // pack the rest all down
    for (int i=index; i<numfaces-1; i++) 
        faces[i] = faces[i+1];
    
    // and decrement number of faces
    numfaces--;
}

void FaceArray::DeleteFace(Face* f)
{
    if (debug) cerr << "FaceArray::DeleteFace(" << f << ")\n";
    
    // get index of face
    int index = getIndex(f);
    
    // do the delete
    DeleteFace(index);
}

void FaceArray::EraseFace(Face * f)
{
    Face *tempFace;
    
    tempFace = f;
    f->unlink();
    DeleteFace( tempFace );
}

void FaceArray::EraseFacesOfColor(int color)
{
    if (debug) cerr << "EraseFacesOfColor()\n";
    
    int open_slot = 0;
    for (int i = 0; i < numfaces; i++) {
        Face* face = getFace(i);
        
        // assume did vertex coloring and other nodes have same color
        Node* n = face->getNodeLink(0); 
        if (n->getPartId() == color) {  
            face->unlink();     // unlink it from its friends
        }
        else { // save it
            *faces[open_slot] = *faces[i];
            open_slot++;
        }
    }
    
    // reset the number of faces left
    numfaces = open_slot;
}

void FaceArray::DeleteAllFaces()
{
    if (debug) cerr << "DeleteAllFaces()\n";
    
    // turns out, all we have to do is to reset the number (is static allocation)
    numfaces = 0;
}

void FaceArray::WriteDXF(ostream& os)
{
    for (int i=0; i<numfaces; i++) {
        Face* face = getFace(i);
        os << "0\n3DFACE\n";
        Node* node1 = face->getNodeLink(0);
        Node* node2 = face->getNodeLink(1);
        Node* node3 = face->getNodeLink(2);
        os << "10\n" << node1->p.x << "\n20\n" << node1->p.y << "\n30\n" 
            << node1->p.z << endl;
        os << "11\n" << node2->p.x << "\n21\n" << node2->p.y << "\n31\n" 
            << node2->p.z << endl;
        os << "12\n" << node3->p.x << "\n22\n" << node3->p.y << "\n32\n" 
            << node3->p.z << endl;
        os << "13\n" << node3->p.x << "\n23\n" << node3->p.y << "\n33\n" 
            << node3->p.z << endl;
    }
}

double FaceArray::getSurfaceArea()
{
    if (debug) cerr << "FaceArray::getSurfaceArea()\n";
    
    // loop over all faces, accumulating surface area
    double area_sum = 0.0;
    for (int i=0; i<numfaces; i++) {
        Face* face = getFace(i);
        area_sum += face->getArea();
    }
    
    // and return it
    return(area_sum);
}

double FaceArray::getVolume()
{
    if (debug) cerr << "FaceArray::getVolume()\n";
    
    // get the minimum value in the y direction (min.y of bbox)
    Point3D min, max;
    nodearray_p->getBoundingBox(&min, &max);
    
    // loop over all faces, accumulating volume
    double volume_sum = 0.0;
    for (int i=0; i<numfaces; i++) {
        Face* face = getFace(i);
        Point3D n = face->getNormal();
        Point3D center = face->getCenter();
        
        double height = center.y - min.y;
        double volume = height * face->getShadowArea();
        
        // see which way the face is pointing
        if (n.y < 0) // pointing in opposite direction as x-z plane normal=neg vol 
            volume = -volume;
        
        // accumulate the volume of this face
        volume_sum += volume;
    }
    
    // return the value
    return(volume_sum);
}

double FaceArray::getExactZVolume(int Num_Face_Surface)
{
    if (debug) cerr << "FaceArray::getExactZVolume()\n";
    // Note that I use Num_Face_Surface and not numfaces on purpose, 
    // because if the object has been extruded, NumFaces != numfaces
    // and we really want to use NumFaces which should be the number
    // of faces of the object before it was extruded.
    
    // loop over all faces, accumulating volume
    double volume_sum = 0.0;
    for (int i=0; i<Num_Face_Surface; i++) 
        volume_sum += getFace(i)->ZVolume();
    return volume_sum;
}

// this function creates a new node at the given S, then calls the
// function below.  We created it because sometimes we need to have
// the new node created, and other times (like when we're splitting
// two adjacent faces), we already have the split node
Node* FaceArray::subdivideFace(Face *initial_p, Point3D S, Face* *t1_p, 
                               Face* *t2_p, Face* *t3_p) 
{
    cerr << "FaceArray::subdivideFace(" << initial_p << ", " << S << ", "
        << t1_p << "," << t2_p << "," << t3_p << "," << ")\n";
    
    // S becomes a new Node 
    cerr << "initial_p is face # " << initial_p->getIndex() << endl;
    int index_S = nodearray_p->addNode(S);
    Node* node_S = nodearray_p->getNode(index_S);
    cerr << " Created new node S of index " << index_S <<": "<< *node_S << endl;
    
    // and call the other subdivide
    Node* ret_ptr = subdivideFace(initial_p, node_S, t1_p, t2_p, t3_p);
    
    // test if we were unsuccessful (if so, kill our created node)
    if (ret_ptr == NULL) nodearray_p->DeleteNode(node_S);
    
    // and return their status/new_node
    return(ret_ptr);
}

// subdivides a face in 2 or 3 faces in function of the position of the point
// in the triangle. If the point is outside, it returns NULL, otherwise it
// returns a pointer to the new node created due to the split and the faces
// that were created
// notice: this function assumes the point S is in the plane of the triangle.
Node* FaceArray::subdivideFace(Face *initial_p, Node* node_S, Face* *t1_p, 
                               Face* *t2_p, Face* *t3_p) 
{
    cerr << "FaceArray::subdivideFace(" << initial_p << ", " << node_S << ", "
        << t1_p << "," << t2_p << "," << t3_p << "," << ")\n";
    cerr << "SUBDIVIDING FACE: " << *initial_p << endl;

    const double precision = 0.001;
    double dist, dot;
    Point3D center = initial_p->getCenter();
    Node *p[3];
    p[0] = initial_p->getNodeLink(0);
    p[1] = initial_p->getNodeLink(1);
    p[2] = initial_p->getNodeLink(2);
    
    Point3D p_S[3];
    Point3D p_p[3];
    Point3D p_center[3];
    Point3D cross1, cross2;
    Face *t1, *t2, *t3;

    for (int i = 0; i < 3; i++){
        p_S[i] = node_S->p - p[i]->p;
        
        // ie, if S is too close to one of the vertices, we do not subdivide
        if (p_S[i].Length() < precision) { 
            cerr << "too close to vertex (length = " << p_S[i].Length() 
                << "<" << precision << "\n";
            return(NULL); 
        }
        
        p_p[i] = p[(i + 1) % 3]->p - p[i]->p;
        p_center[i] = center - p[i]->p;
        cross1 = p_S[i].Cross(p_p[i]);
        cross2 = p_center[i].Cross(p_p[i]);
        dot = cross1.Dot(cross2);
        
        if (dot < -precision) { 
            cerr << "dot < 0 (=" << dot << ")\n"; 
            return(NULL); 
        }
        
        dist = cross1.Length() / p_p[i].Length();
        
        //case when S in on an edge of the initial face so we subdivide in 2 faces
        if(dist <= precision) {
            cerr << " split into 2 faces\n";
            
            //should destroy the old sphere tree before mucking with subdivide
            // avoid coredumps from other objects colliding with this one
            getObject()->getObjectArray()->ResetCollisions();
            getObject()->DestroySphereTree();
            
            // unlink the parent face from its nodes and edges
            // this way, the edges will be open for the new faces to bind to
            initial_p->unlink();
            
            // and delete the face from our list
            DeleteFace(initial_p);
            
            // create 2 new faces
            int new_index1 = addFace(); int new_index2 = addFace();
            t1 = getFace(new_index1); t2 = getFace(new_index2);
            
            // initialize our new faces and link in to new edges
            t1->init(nodearray_p, edgearray_p, this, node_S, p[(i+1) % 3],
                p[(i+2)%3]);
            t2->init(nodearray_p, edgearray_p, this, node_S, p[(i+2) % 3],p[i]);
            cerr << "split_node = " << *node_S << endl;
            cerr << "subtri t1: " << *t1 << endl;
            cerr << "subtri t2: " << *t2 << endl;
            
            // tell new faces to compute their normals
            t1->computenormal();
            t2->computenormal();
            
            // set up faces for return
            if (t1_p) *t1_p = t1;
            if (t2_p) *t2_p = t2;
            if (t3_p) *t3_p = NULL;
            
            getObject()->CreateSphereTree();
            return(node_S);
        }
    }
    
    /* case when S inside the initial face so we subdivide in 3 faces */
    cerr << " split into 3 faces\n";
    
    //should destroy the old sphere tree before mucking with subdivide
    // avoid coredumps from other objects colliding with this one
    getObject()->getObjectArray()->ResetCollisions();
    getObject()->DestroySphereTree();
    
    // create the 3 new faces
    int new_index1 = addFace();
    int new_index2 = addFace(); 
    int new_index3 = addFace();
    t1 = getFace(new_index1); 
    t2 = getFace(new_index2); 
    t3 = getFace(new_index3);
    
    // unlink the parent face from its nodes and edges
    // this way, the edges will be open for the new faces to bind to
    initial_p->unlink();
    
    // delete it from our list
    DeleteFace(initial_p);
    
    // initialize our new faces and link in to new edges
    t1->init(nodearray_p, edgearray_p, this, node_S, p[1], p[2]);
    t2->init(nodearray_p, edgearray_p, this, node_S, p[2], p[0]);
    t3->init(nodearray_p, edgearray_p, this, node_S, p[0], p[1]);
    
    // tell new faces to compute their normals
    t1->computenormal();
    t2->computenormal();
    t3->computenormal();
    
    cerr << "split_node = " << *node_S << endl;
    cerr << "subtri t1: " << *t1 << endl;
    cerr << "subtri t2: " << *t2 << endl;
    cerr << "subtri t3: " << *t3 << endl;
    
    if (t1_p) *t1_p = t1;
    if (t2_p) *t2_p = t2;
    if (t3_p) *t3_p = t3;
    
    getObject()->CreateSphereTree();
    return(node_S);
}

/* This function undoes the changes done by subdivideFace -- given a node, it
* deletes this node, incident faces and edges, and creates a new face which
* is the union of the deleted faces */
Face* FaceArray::unSubdivideFace(Node *n_in)
{
    printf("unsubdividing ...\n");
    // destroy old sphere tree before changes
    // avoid coredumps from other objects colliding with this one
    (getObject()->getObjectArray())->ResetCollisions();
    getObject()->DestroySphereTree();
    int i;
    int n_facedegree = n_in->numFaceLinks();
    int n_edgedegree = n_in->numEdgeLinks();
    Face* f[3];
    Edge* e[3];
    Node* n[3];
    // make sure n is connected to 3 edges and 2 or 3 faces (2 cases of subdiv)
    if (n_facedegree < 2 || n_facedegree > 3 || n_edgedegree != 3) {
        printf("node not created by subdivideFace, or resubdivided fd %d ed %d\n",
            n_facedegree, n_edgedegree);
        // if we subdivided and then resubdivided one of the new faces, we must
        // unsubdivide in reverse order, or else it will bail out here!!
        cerr << "more node info: " << *n_in << endl;
        return(NULL);
    }
    //get the incident faces and edges, and the 3 nodes connected to these edges
    for (i = 0; i < n_facedegree; i++) {
        f[i] = n_in->getFaceLink(i);
        cerr << "old face " << *f[i] << endl;
    }
    for (i = 0; i < n_edgedegree; i++) {
        e[i] = n_in->getEdgeLink(i);
        n[i] = e[i]->getOtherNode(n_in);
    }
    // now unlink the faces from adjacent nodes and edges, and delete them
    for (i = 0; i < n_facedegree; i++) {
        f[i]->unlink();
        DeleteFace(f[i]);
    }
    // now unlink the edges from adjacent nodes and delete them
    for (i = 0; i < n_edgedegree; i++) {
        n[i]->deleteEdgeLink(e[i]);
        n_in->deleteEdgeLink(e[i]);  //unnecessary since we're about to delete it
        edgearray_p->DeleteEdge(e[i]);
    }
    // now delete the given node
    nodearray_p->DeleteNode(n_in);
    // now create a new face from the 3 nodes
    i = addFace();
    Face* f_out = getFace(i);
    f_out->init(nodearray_p, edgearray_p, this, n[0], n[1], n[2]);
    f_out->computenormal(); // why?
    cerr << "new face " << *f_out << endl;
    getObject()->CreateSphereTree();
    return(f_out);
}

void FaceArray::SanityCheck(NodeArray* real_na_p, EdgeArray* real_ea_p)
{
    string s;
    
    if (!nodearray_p) s += "nodearray is NULL!";
    if (nodearray_p != real_na_p) s += "nodearray is BAD!";
    if (!edgearray_p) s += "edgearray is NULL!";
    if (edgearray_p != real_ea_p) s += "edgearray is BAD!";
    if (numfaces < 0) s += "numfaces < 0!";
    if (numfaces > maxfaces) s += "numfaces > maxfaces!";
    if (maxfaces < 0) s += "maxfaces < 0!";
    
    if (s.size() > 1) 
        cerr << "FaceArray is INSANE: " << s.c_str() << endl;
    
    // check all our faces
    for (int i=0; i<numfaces; i++)
        faces[i]->SanityCheck(this);
}

// Eliminates any faces that are insane
void FaceArray::Cleanup()
{
    // check all our faces
    for (int i=0; i<numfaces; i++) {
        Face* f = faces[i];
        if (!f->SanityCheck(this)) {    // is insane- kill it
            f->unlink();        // unlink it from it's edges/nodes
            DeleteFace(i);      // and kill it
            
            i--;        // backup the counter because the DeleteFace() shifted all down
        }
    }
}

double FaceArray::DistanceToNearestFace(Point3D p)
{
    double min_dist = 1e20;
    for (int i=0; i<numfaces; i++) {
        Face* f = getFace(i);
        
        // NOTE: THIS IS NOT THE RIGHT FUNCTION!  WE NEED A
        // DistanceToPoint(p) FUNCTION THAT TELLS THE DISTANCE
        // FROM THE GIVEN POINT TO THE FACE (NOT JUST IT'S PLANE)
        //double cur_dist = fabs(f->PlaneDistanceToPoint(p));
        
        // for now, approximate by distance to center
        Point3D face_center = f->getCenter();
        double cur_dist = p.Dist(face_center);
        
        if (cur_dist < min_dist)
            min_dist = cur_dist;
    }
    
    return(min_dist);
}


//finds the closest face of the object to the given pt point.
Face* FaceArray::ClosestFaceToPoint(Point3D pt)
{
    int numF = getNumFaces();
    Face* current = getFace(0);
    Face* closest = current;
    double minDist = current->FaceDistanceToPoint(pt); 
    for (int i = 1; i < numF; i ++){
        current = getFace(i);
        double currDist = current->FaceDistanceToPoint(pt);
        if( minDist > currDist ){
            minDist = currDist;
            closest = current;
        }
    }
    
    return(closest);
}

// this function reorders the facearray based on a viewing matrix
// in this way, we can then display the triangles from farthest Z-value
// to nearest Z-value (w.r.t. the viewer) to make transparency look good
// so, this function gives us the array as farthest to nearest faces
void FaceArray::ReorderFacesBasedOnViewingMatrix(double* model, double* proj, 
                                                 GLint* view)
{
    // get some local storage
    double* depths = new double[getNumFaces()];
    double min = 1e20, max = -1e20;
    
    { // for each face, project it to the viewer and save the resulting depth
        for (int i=0; i<getNumFaces(); i++) {
            // get the face, then its center
            Face* f = getFace(i);
            Point3D center = f->getCenter();
            
            // project it
            double projx, projy, depth;
            gluProject(center.x,center.y,center.z, model,proj,view, 
                &projx,&projy,&depth);
            
            // save the value
            depths[i] = depth;
            if (depth > max) max = depth;
            if (depth < min) min = depth;
        }
    }
    
    // debugging
    if (debug) cerr << "Max = " << max << ", " << "Min = " << min << endl;
    
    { // now reorder all the triangles in our list based on their depth
        // implements a very dumb O(n^2) sort
        for (int i=0; i<getNumFaces(); i++) {
            for (int j=i+1; j<getNumFaces(); j++) {
                
                // if face i is farther than face j, swap `em
                if (depths[i] < depths[j]) {
                    Face* tmp = faces[i];
                    faces[i] = faces[j];
                    faces[j] = tmp;
                    
                    double d_tmp = depths[i];
                    depths[i] = depths[j];
                    depths[j] = d_tmp;
                }
            }
        }
    }
    
    // clean up
    delete [] depths;
}

// Note: the Node::getIndex() call is slow, unless we first call 
// NodeArray::CacheIndices(), which we're assuming we have
void FaceArray::GeometryReplicate(GeometryReplicator* gr)
{
  if (debug) cerr << "FaceArray::GeometryReplicate(" << gr << ")\n";

  // locals
  const int buf_size = 1490;    //1500 is native packet side - need some padding
  char buffer[buf_size];
/*
  // loop, sending each face's info
  for (int i=0; i<getNumFaces(); i++) {
        Face* face = getFace(i);
        Node* i0 = face->getNodeLink(0);
        Node* i1 = face->getNodeLink(1);
        Node* i2 = face->getNodeLink(2);
        int s = sprintf(buffer, "%c%c%d %d %d %d", GR_CODE, GR_FACE, i, 
                i0->getIndex(), i1->getIndex(), i2->getIndex());
        gr->SendMessage(buffer, s);
  }
*/

        int numNodesInBuffer = 0;
        int faceEntry[3];
        u_long* netFaceEntry = (u_long*)faceEntry;

        int nodeEntrySize = sizeof(int) + 3*sizeof(int);

        u_long netvalue;

    int s = sprintf(buffer, "%c%c", GR_CODE, GR_FACE);
        netvalue = htonl(numNodesInBuffer);
        memcpy(&buffer[s],&netvalue,sizeof(int));
        s += sizeof(int);
    
        for (int i=0; i<getNumFaces(); i++) 
        {
                if (s + nodeEntrySize < buf_size)
                {
                        Face* face = getFace(i);
                        Node* i0 = face->getNodeLink(0);
                        Node* i1 = face->getNodeLink(1);
                        Node* i2 = face->getNodeLink(2);

                        faceEntry[0] = i0->getIndex();
                        faceEntry[1] = i1->getIndex();
                        faceEntry[2] = i2->getIndex();

                        netFaceEntry[0] = htonl(netFaceEntry[0]);
                        netFaceEntry[1] = htonl(netFaceEntry[1]);
                        netFaceEntry[2] = htonl(netFaceEntry[2]);


                        //place the index in the entry
                        netvalue = htonl(i);

                        memcpy(&buffer[s],&netvalue,sizeof(int));
                        s += sizeof(int);
                        //send all the node data
                        memcpy(&buffer[s],netFaceEntry,3*sizeof(int));
                        s += 3*sizeof(int);

                        numNodesInBuffer++;

                }

                if ((s + nodeEntrySize >= buf_size) || (i == getNumFaces() - 1) )
                {
                        //copy in the actual number of nodes
                        netvalue = htonl(numNodesInBuffer);
                        memcpy(&buffer[2],&netvalue,sizeof(int));
                        gr->SendMessage(buffer, s);

                        //reload the buffer with the header and reset s
                        s = sprintf(buffer, "%c%c", GR_CODE, GR_FACE);
                        numNodesInBuffer = 0;
                        netvalue = htonl(numNodesInBuffer);
                        memcpy(&buffer[2],&netvalue,sizeof(int));
                        s += sizeof(int);
                }
    }


}

ostream& operator<<(ostream& os, const FaceArray& fa)
{
    os << fa.numfaces << endl;
    for (int i=0; i<fa.numfaces; i++) {
        const FaceArray* fa_p = &fa;
        Face* cur_face = ((FaceArray*)fa_p)->getFace(i);
        os << (*cur_face) << endl;
    }
    return(os);
}

---