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

Go to the documentation of this file.

/* $Id: edgearray.cpp,v 1.54 2006/05/25 22:23:11 craig Exp $ */
#include "edgearray.h"
#include "object.h"
#include "nodearray.h"
#include "edge.h"
#include "node.h"
#include "object.h"
#include "knot.h"
#ifdef _WIN32
#include <glut.h>
#else
#include <GLUT/glut.h>
#include <OpenGL/gl.h>
#endif
#include <stdlib.h>


#include <string>

using namespace std;

const char* EdgeArray::rcsid = "@(#) $Id: edgearray.cpp,v 1.54 2006/05/25 22:23:11 craig Exp $ $Copyright: (c)2001 National Biocomputation Center, Stanford University $";

int EdgeArray::debug = 0;
int default_maxedges = 1000;

EdgeArray::EdgeArray(NodeArray* nodearray_p_in)
: nodearray_p(nodearray_p_in)
{ 
    edges = NULL;
    numedges = maxedges = 0; 
// !!! TEMPORARY
    do_nonlinear = 1;
    strcpy(edgefilename, "edge.data");
}

EdgeArray::~EdgeArray()
{ 
    reset();
    numedges = maxedges = -1;
}

void EdgeArray::reset(void)
{ 
    if (edges) {
        // delete all the edges
        for (int i=0; i<numedges; i++)
            delete edges[i];
        
        // and delete the array itself
        delete [] edges; 
    }
    
    numedges = maxedges = 0;
}

void EdgeArray::scaleRestLengths(double scalefactor)
{ 
    for (int i=0; i<numedges; i++)
        edges[i]->scaleRestLength(scalefactor);
}

void EdgeArray::computeNewRestLengths()
{ 
    for (int i=0; i<numedges; i++) 
        edges[i]->computeNewRestLength();
}

int EdgeArray::getIndex(Edge* e)
{
    // go look for it
    for (int i=0; i<numedges; i++)
        if (edges[i] == e) 
            return(i);
        
        // didn't find it
        return(-1);
}

void EdgeArray::setEdgeFilename(char *s)
{
    strcpy(edgefilename, s);
}

void EdgeArray::allocateEdges(int maxedges_in)
{ 
    numedges=0;
    maxedges = maxedges_in;
    default_maxedges = maxedges;    // save as new default
    edges = new Edge*[maxedges]; 
    if (edges == NULL) {
        cerr << "EdgeArray: Could not allocate memory for " << maxedges 
            << " edges!\n";
    }
}

int EdgeArray::addEdge(int i0, int i1, double restLength,
                       double springConstant, double dampConstant) 
{ 
    //  get the real nodes
    Node* n0 = nodearray_p->getNode(i0);
    Node* n1 = nodearray_p->getNode(i1);
    
    // and call the other function
    int s = addEdge(n0, n1, restLength, springConstant, 
        dampConstant);
    
    // and report their value
    return(s);
}

int EdgeArray::addEdge(Node* n0, Node* n1, double restLength,
                       double springConstant, 
                       double dampConstant)
{ 
    // if no space, go and get some
    if (maxedges == 0) allocateEdges(default_maxedges);
    
    // do it
    edges[numedges] = new Edge;
    edges[numedges]->init(nodearray_p,this, n0,n1,
        restLength, springConstant,dampConstant);
    
    // if beyond end, report error
    if (numedges >= maxedges) {
        cerr << "EdgeArray: blew bounds of memory! (" << maxedges << ")\n";
        return(numedges);
    }
    else
        return(numedges++);
}

/* function to open the edge data file and extract the data into
* edges[] and numedges */
int EdgeArray::getEdgeData() {
    FILE *edgefile;
    char s[MAXLINE];
    int numedges_given;
    
    // open the edge data file
    if ((edgefile = fopen(edgefilename, "r")) == NULL) {
        printf("Could not open edge data file\n");
        return 0;
    }
    // get the number of edges from the first line
    if (myfgets(s, MAXLINE, edgefile) == NULL) {
        printf("No lines in edge data file\n");
        return 0;
    }
    if (sscanf(s, "%d", &numedges_given) != 1) {
        printf("Could not read numedges\n");
        return 0;
    }
    printf("%d edges\n", numedges_given);
    printf("Edge# Node1# Node2#\n");
    
    // allocate space for Edge objects + 50%
    allocateEdges((int)(numedges_given * 1.5));
    
    // loop through and get data for each edge
    for (int i = 0; i < numedges_given; i++) {
        if (myfgets(s, MAXLINE, edgefile) == NULL) {
            printf("Not enough edges in data file\n");
            return 0;
        }
        
        // initialize one edge with current data line
        int index = addEdge(0,0,0,0,0);    // create a dummy one
        if (edges[index]->init(nodearray_p, this, s) == 0) {    // load it up
            printf("edges[%d].init failed\n", i);
            return 0;
        }
        if (debug)
            printf("%d %d %d\n", edges[i]->getIndex(),
            nodearray_p->getIndex(edges[i]->getNodeLink(0)),
            nodearray_p->getIndex(edges[i]->getNodeLink(1)));
    }
    if (fclose(edgefile) == EOF) {
        printf("Could not close edge data file\n");
        return 0;
    }
    
    return 1;
}

void EdgeArray::SaveFrdEdgeData(ostream& os)
{
    nodearray_p->CacheIndices();
    for (int i=0; i<numedges; i++) {
        Edge* ed = getEdge(i);
        os << "e";
        for (int j=0; j<2; j++) {
            Node* n = ed->getNodeLink(j);
            os << " " << n->getCachedIndex() + 1;    // Frd is 1-index based (as SMF...)
        }
        os << endl;
    }
}

int EdgeArray::getFRDEdgeData(FILE *meshfile)
{
    char s[MAXLINE];
    int i, numedges_given;
    
    // read file once to get the number of verts and polys
    numedges_given = 0;
    while (!feof(meshfile)) {
        char s[256];
        fgets(s, 256, meshfile);
        if (feof(meshfile)) break;
        switch (s[0]) {
        case 'e' : numedges_given++; break;
        case 'f' : // faces
        case 'v' : // vertex 
        case '#' : // comment
        default  : // who knows?
            break;
        }
    }
    
    // reset input file pointer
    fseek(meshfile, 0, SEEK_SET);
    
    printf("%d edges\n", numedges_given);
    
    // allocate space for Edge objects + 50% for new edges
    allocateEdges((int)(numedges_given * 1.5));
    
    // loop through and get data for each edge
    for (i = 0; i < numedges_given; i++) 
    {
        do 
        {
            if (myfgets(s, MAXLINE, meshfile) == NULL) 
            {
                printf("Not enough edges in data file\n");
                return(0);
            }
        } while (s[0] != 'e');
        
        { // initialize Edge with current data line
            int i1, i2;
            int numread = sscanf(s, "%*c %d %d", &i1, &i2);
            
            if (debug) cerr << "Edge #" << i << ": " << i1 << " " 
                << i2 << endl;
            
            if (numread != 2) 
            {
                printf("Bad read on edge %d- '%d %d' - IGNORED\n", i, i1,i2);
                continue;
            }
            
            // FRD is 1-based, we are 0-index-based (as SMF...)
            i1--; i2--;
            
#ifdef DO_READ_SANITY_CHECKS
            // good read: set edge info
            int new_index = addEdge(i1,i2);
            // do a sanity check to make sure
            Edge* new_edge = getEdge(new_index);
            int ok = new_edge->SanityCheck(this);
            
            // if there was a problem delete
            if (!ok) 
            {
                //new_edge->unlink();
                DeleteEdge(new_edge);
            }
#else
            // good read: set edge info
            addEdge(i1,i2);
            
#endif
        }
        
    }
    
    // return success
    return(1);
}

/* function to calculate the torsion force between two faces for every edge,
* and apply this to the exterior nodes of these faces */
void EdgeArray::addTorsions(void)
{
    // loop through all the edges
    for (int i = 0; i < numedges; i++) {
        edges[i]->addTorsion();
    }
}

// Calculates the force vector for each edge based on spring length, spring
// constant, and rest length.
// Requires the position of the end nodes.
// Note that the force is directional, depending if the spring is compressed or streteched.
// Using this avoids computing the force for each end node.
void EdgeArray::ComputeForceVectors(int dampingType)
{
        int i;

        switch (dampingType) {
        case FULL:
                // Loop over all edges.
                for (i = 0; i < numedges; i++) {
                        edges[i]->ComputeSpringForceFullDamp();
                }
                break;
        case REL:
                // Loop over all edges.
                for (i = 0; i < numedges; i++) {
                        edges[i]->ComputeSpringForceRelDamp();
                }
                break;
        default:
                // ??? What about ABS - I think this will be computed in the nodes - more complex.
                for (i = 0; i < numedges; i++) {
                        edges[i]->ComputeSpringForce(dampingType);
                }
                break;
        }

}

// Calculates the force vector for each edge based on spring length, spring3
// constant, and rest length.
// Important: This version increase the force if the distance is too far from the rest length!
// Requires the position of the end nodes.
// Note that the force is directional, depending if the spring is compressed or streteched.
// Using this avoids computing the force for each end node.
void EdgeArray::ComputeForceVectorsNonlinear(int dampingType)
{
        // Loop over all edges, computing the non-linear spring force for each.
        int i;
        for (i = 0; i < numedges; i++) {
                edges[i]->ComputeSpringForceNonlinear(dampingType);
        }
}

// this function projects all the edges of a thread type object onto a given
// plane (the z-plane), and finds all the crossings, as a precursor to 
// finding a knot -- then it finds the knot
void EdgeArray::KnotIdentify(int testing)
{
  extern char knotstring[100]; // keeps old global value

  double zdowker[MAXCROSS];
  int zorder[MAXCROSS];
  int zdk[MAXCROSS/2];
  int zcount = 0;
  int i, j;

  // one dummy edge, connects first and last nodes to close the knot
  Edge edummy;
  edummy.init(NULL, NULL, getEdge(numedges-1)->getNodeLink(1),
              getEdge(0)->getNodeLink(0));

  // this part actually finds all the 2D crossings
  for (i = 0; i < numedges; i++) {
    // get info for the first line segment, n1 + s*(n2-n1)
    Edge *ea = getEdge(i);
    Node *n1 = ea->getNodeLink(0);
    Node *n2 = ea->getNodeLink(1);
    double x1 = n1->p.x; double y1 = n1->p.y; double z1 = n1->p.z;
    double x2 = n2->p.x; double y2 = n2->p.y; double z2 = n2->p.z;

    for (j = i+2; j < numedges+1; j++) {
      // get info for the second line segment, n3 + t*(n4-n3)
      Edge *eb;
      if (j == numedges)
        eb = &edummy;
      else
        eb = getEdge(j);
      Node *n3 = eb->getNodeLink(0);
      Node *n4 = eb->getNodeLink(1);
      double x3 = n3->p.x; double y3 = n3->p.y; double z3 = n3->p.z;
      double x4 = n4->p.x; double y4 = n4->p.y; double z4 = n4->p.z;

      // Solve for intersection point when projected to z = 0 plane:
      // (see http://astronomy.swin.edu.au/pbourke/geometry/lineline2d/ )
      double x2minusx1 = x2 - x1;
      double x3minusx1 = x3 - x1;
      double x4minusx3 = x4 - x3;
      double y2minusy1 = y2 - y1;
      double y3minusy1 = y3 - y1;
      double y4minusy3 = y4 - y3;
      double denom = x2minusx1*y4minusy3 - x4minusx3*y2minusy1;
      // check for zero denom (and possibly zero num)
      if (denom == 0.0)
        ; // skip for now
      else {
        double snum = x3minusx1*y4minusy3 - y3minusy1*x4minusx3;
        double tnum = x3minusx1*y2minusy1 - y3minusy1*x2minusx1;
        double s = snum/denom;
        double t = tnum/denom;
        if (s >= 0.0 && s < 1.0 && t >= 0.0 && t < 1.0) {
          if (zcount >= MAXCROSS) {
            // too many crossings, bail out
            strcpy(knotstring, "too many crossings");
            return;
          }
          // add s,t to make these exact, segment can be in 2 crossings
          zdowker[zcount++] = i + s;
          if (zdowker[zcount-1] == 0.0) zdowker[zcount-1] = 0.01;
          zdowker[zcount++] = j + t;
          if (testing)
            printf("z Segment %5.2f crosses ", i+s);
          // Unproject to get over/under/through
          double zs = z1 + s*(z2 - z1);
          double zt = z3 + t*(z4 - z3);
          if (testing) {
            if (zs > zt)
              printf("over");
            else if (zs < zt)
              printf("under");
            else // zs == zt
              printf("through");
            if (t > 0.99) t = 0.99; // only affects this printf
            printf(" segment %5.2f \n", j+t);
          }
          // negate the over segment
          if (zs > zt)
            zdowker[zcount - 2] = -zdowker[zcount - 2];
          else if (zs < zt)
            zdowker[zcount - 1] = -zdowker[zcount - 1];
        }
      }
    } // for j
  } // for i

  // the first crossnumber will always be 1, by the way we find the
  // intersections.  They will also either all be in even/odd pairs, or it
  // is definitely not a knot without some extra closing segment.

  // first pair off and number all the crossing segments from 1 to 2*crossnum
  for (i = 0; i < zcount; i++) {
    double crosssegment = fabs(zdowker[i]);
    int crossnumber = 0;
    for (j = 0; j < zcount; j++)
      if (crosssegment >= fabs(zdowker[j]))
        crossnumber++;
    if (zdowker[i] > 0.0)
      zorder[i] = crossnumber;
    else
      zorder[i] = -crossnumber;
    if (testing)
      printf("%d ", zorder[i]);
  }
  if (testing)
    printf("\n");

  if (zcount == 0) { // no crossings
    strcpy(knotstring, "    unknot");
    return;
  }

  // next convert to Dowker-Thistlethwaite notation, list of even numbers which
  // are paired with 1,3,5,... (+ for under, - for over)
  for (i = 0; i < zcount; i = i+2) { // zcount is multiple of 2
    if (abs(zorder[i])%2 == 0 && abs(zorder[i+1])%2 == 1)
      zdk[abs(zorder[i+1])/2] = zorder[i];
    else if (abs(zorder[i])%2 == 1 && abs(zorder[i+1])%2 == 0)
      zdk[abs(zorder[i])/2] = zorder[i+1];
    else { // uh oh, a pair of two odd or two even crossings
      if (testing)
        printf("both odd or both even, problem\n");
      strcpy(knotstring, "odd/even problem");
      return;
    }
  }

  // ok, now create a knot, and send it to our knot finder
  int knots[10][MAXCROSS];  //up to 10 knots, should really be 2+MAXCROSS/2?
  int ret;
  Knot k;
  char* knotname;

  knots[0][0] = zcount/2;
  knots[0][1] = 1;
  if (testing)
    printf("%d  %d  ", zcount/2, 1);
  for (i = 0; i < zcount/2; i++) {
    knots[0][i+2] = zdk[i];
    if (testing)
      printf("%d ", zdk[i]);
  }
  if (testing)
    printf("\n");
    
  // send in our knot, get back array of prime factors
  ret = k.find_init(knots[0], &knots[0]);
  if (testing)
    printf("found %d knots\n", ret);
  if (ret > 0) {
    strcpy(knotstring, "");
    for (i = 0 ; i < ret; i++) {
      knotname = k.get_knotid(knots[i]);
      if (testing)
        printf("knot %d: %s : ", i, knotname);
      if (i > 0)
        strcat(knotstring, "+");
      strcat(knotstring, knotname);
      if (testing) {
        for (j = 0; j < knots[i][0] + 2; j++)
          printf("%d ", knots[i][j]);
        printf("\n");
      }
    }
  }
  else
    strcpy(knotstring, "    unknot");

}



/* function to draw a line for every edge
* using energy to decide color */
void EdgeArray::DrawEdges()
{
  /* We really don't care about this energy stuff
    double totalenergy, energypercent;
    
    totalenergy = getEnergy() + nodearray_p->getEnergy();
    // loop through all the edges
    for (int i = 0; i < numedges; i++) {
        // use energy to decide color of each edge
        if (totalenergy == 0.0)
            energypercent = 0.0;
        else
            energypercent = (edges[i]->getEnergy())/totalenergy;

        // for if drawing the compression/stretching of edge
        // ranges from white (no energy) to red (100% of system's energy)
        glColor3d(1.0, (1.0 - energypercent), (1.0 - energypercent));

                // hack for suture/microsurgery application
        if ((getNodeArray()->getObject()->getType() == Object::edges_only) &&
            (i < RIGIDNEEDLE))
           glColor3d(0.5, 0.5, 0.5);

        // just do the draw
        edges[i]->draw();
    }
  */

  // if we're texturing and drawing a thread, let's add some rope texture
  // first, draw quads.  perpendicular lines should be EDGETHRESH in length,
  // centered on the edge, bisecting two adjacent edges, and ideally, facing
  // the camera.
  extern int do_textures;
  // even if texture is NULL, the glTexCoord commands don't complain
  // so draw the quads for any edge object, if we're doing textures
  if (do_textures && getObject()->getType() == Object::edges_only) {
    extern int random_flag;
    extern int testing_flag;
    double EDGETHRESH = getObject()->getEdgeThresh();

    // this does no smoothing, but links up the first and last quads
    if (strstr(getObject()->getName(), "ring")) {
    Point3D cam_vec, cam_pos;
    cam_pos = getCameraPos();
    Point3D perp, perp0, perp1, oldperp, v0, v1, pt, u0, u1;
    Edge *e0; Edge *e1;
    double widthpct = 0.0;  // for the u (width) texture coordinate
    glColor4f(1.0, 1.0, 1.0, (getObject()->getDiffuseColor())[3]);
    glBegin(GL_QUADS);
    // head -- the beginning of the first quad
    e1 = getEdge(0);
    u1 = e1->UnitDirectionVector();
    pt = e1->getNodeLink(0)->p;
    // this is perpendicular to edge and vector from camera to edge
    cam_vec = pt - cam_pos;
    perp1 = u1.Cross(cam_vec);
    perp1.Normalize();
    // HACK for ring object:
    e0 = getEdge(numedges-1);
    u0 = e0->UnitDirectionVector();
    perp0 = u0.Cross(e0->getNodeLink(0)->p - cam_pos);
    perp0.Normalize();
    perp = perp0 + perp1;
    perp.Normalize();
    oldperp = perp;
    // end HACK
    v0 = pt + 0.5*EDGETHRESH*perp;
    v1 = pt - 0.5*EDGETHRESH*perp;
    glTexCoord2f(widthpct, 0.0);
    glVertex3d(v0.x, v0.y, v0.z);
    glTexCoord2f(widthpct, 1.0);
    glVertex3d(v1.x, v1.y, v1.z);
    //interior
    for (int i = 1; i < numedges; i++) {
      perp0 = perp1;
      e0 = e1;
      u0 = u1;
      e1 = getEdge(i);
      u1 = e1->UnitDirectionVector();
      pt = e1->getNodeLink(0)->p;
      perp1 = u1.Cross(cam_vec);
      perp1.Normalize();
      // average of perpendiculars to each edge
      perp = perp0 + perp1;
      perp.Normalize();
      v0 = pt + 0.5*EDGETHRESH*perp;
      v1 = pt - 0.5*EDGETHRESH*perp;
      // since our texture is 128 x 32, only use 1/4 of it per quad:
      widthpct += 0.25;
      // end previous quad
      glTexCoord2f(widthpct, 1.0);
      glVertex3d(v1.x, v1.y, v1.z);
      glTexCoord2f(widthpct, 0.0);
      glVertex3d(v0.x, v0.y, v0.z);
      // and begin next quad
      //if (widthpct > 0.99) // eg = 1.0
      //  widthpct = 0.0;    // perhaps unnecessary because of GL_REPEAT
      glTexCoord2f(widthpct, 0.0);
      glVertex3d(v0.x, v0.y, v0.z);
      glTexCoord2f(widthpct, 1.0);
      glVertex3d(v1.x, v1.y, v1.z);
    }
    //tail
    pt = e1->getNodeLink(1)->p;
    perp = oldperp; // HACK for ring object
    v0 = pt + 0.5*EDGETHRESH*perp;
    v1 = pt - 0.5*EDGETHRESH*perp;
    widthpct += 0.25;
    glTexCoord2f(widthpct, 1.0);
    glVertex3d(v1.x, v1.y, v1.z);
    glTexCoord2f(widthpct, 0.0);
    glVertex3d(v0.x, v0.y, v0.z);
    glEnd();
    } // if not smoothing, but matching up first and last quads
    // this does parabolic smoothing of adjacent quads
    else {
    Point3D cam_vec, cam_pos;
    cam_pos = getCameraPos();
    Point3D perp, oldperp, v0, v1, pt, oldpt, u0, u1;
    Edge *e0; Edge *e1;
    double widthpct = 0.0;  // for the u (width) texture coordinate
    glColor4f(1.0, 1.0, 1.0, (getObject()->getDiffuseColor())[3]);
    glBegin(GL_QUADS);
    // head -- the beginning of the first quad
    e1 = getEdge(0);
    u1 = e1->UnitDirectionVector();
    pt = e1->getNodeLink(0)->p;
    // this is perpendicular to edge and vector from camera to edge
    cam_vec = pt - cam_pos;
    perp = u1.Cross(cam_vec);
    perp.Normalize();
    v0 = pt + 0.5*EDGETHRESH*perp;
    v1 = pt - 0.5*EDGETHRESH*perp;
    extern int needle;
    if (needle) { // this draws a pointy needle
      //glColor4f(0.5, 0.5, 0.5, (getObject()->getDiffuseColor())[3]);
      v0 = pt;
      v1 = pt;
    }
    if (!needle)
      glTexCoord2f(widthpct, 0.0);
    glVertex3d(v0.x, v0.y, v0.z);
    if (!needle)
      glTexCoord2f(widthpct, 1.0);
    glVertex3d(v1.x, v1.y, v1.z);
    // interior
    for (int i = 1; i < numedges; i++) {
      e0 = e1;
      u0 = u1;
      e1 = getEdge(i);
      u1 = e1->UnitDirectionVector();
      // now find the angle between the two edge vectors
      double dot = u0.Dot(u1);
      double theta = (acos(dot))*(180.0/PI);  // theta between 0 and 180
      int k = (int)theta/15; //add an extra quad (per side) per 15 degrees
      k = MIN(k, 5);
      if ((i == 1) && needle)
        k = 0;
      for (int j = -k; j <= k; j++) {
        oldpt = pt;
        oldperp = perp;
        double param;
        // ok, param is a number in the interval (1-sqrt(2)/2, sqrt(2)/2),
        // where .414 is the size of that interval (.707 - .293)
        // Why?  Because parametrization A(1-t)(1-t) + B(t)(1-t) + C(t)(t)
        // is A/2 + blah for t = .293, and is blah + C/2 for t = .707
        // and is A/4 + B/2 + C/4 for param = .5
        param = 0.5 + (0.414*j)/(2*k+1);
        widthpct = 0.25*e1->getIndex() + 0.25*j/(2*k+1);
        pt = (1.0-param)*(1.0-param)*e0->getNodeLink(0)->p +
          2.0*param*(1.0-param)*e1->getNodeLink(0)->p +
          param*param*e1->getNodeLink(1)->p;
        if ((i == 1) && needle) // b/c needle is different length
          pt = e1->getNodeLink(0)->p;
        cam_vec = pt - cam_pos;
        perp = (pt - oldpt).Cross(cam_vec);
        perp.Normalize();
        v0 = pt + 0.5*EDGETHRESH*perp;
        v1 = pt - 0.5*EDGETHRESH*perp;
        // if perp flipped direction, unflip it for previous quad:
        dot = perp.Dot(oldperp);
        theta = (acos(dot))*(180.0/PI);  // theta between 0 and 180
        //if (theta > 135.0)   //test to see where we flipped
        //glColor3f(1.0,0.0,0.0);
        //else
        //glColor3f(1.0,1.0,1.0);
        if (theta > 135.0) {
          v0 = pt - 0.5*EDGETHRESH*perp;
          v1 = pt + 0.5*EDGETHRESH*perp;
        }
        // end previous quad
        if (i > 1 || !needle)
          glTexCoord2f(widthpct, 1.0);
        glVertex3d(v1.x, v1.y, v1.z);
        if (i > 1 || !needle)
          glTexCoord2f(widthpct, 0.0);
        glVertex3d(v0.x, v0.y, v0.z);
        // and now reflip perp for next quad if we unflipped it
        if (theta > 135.0) {
          v0 = pt + 0.5*EDGETHRESH*perp;
          v1 = pt - 0.5*EDGETHRESH*perp;
        }
        // and begin next quad
        if (needle && i == 1) // after needle, get rope color correct
          glColor4f(1.0, 1.0, 1.0, (getObject()->getDiffuseColor())[3]);
        glTexCoord2f(widthpct, 0.0);
        glVertex3d(v0.x, v0.y, v0.z);
        glTexCoord2f(widthpct, 1.0);
        glVertex3d(v1.x, v1.y, v1.z);
      } // for j (extra quads for large angles between edges)
    } // for i (quad per edge)
    // tail -- the end of the last quad
    pt = e1->getNodeLink(1)->p;
    cam_vec = pt - cam_pos;
    perp = u1.Cross(cam_vec);
    perp.Normalize();
    v0 = pt + 0.5*EDGETHRESH*perp;
    v1 = pt - 0.5*EDGETHRESH*perp;
    widthpct = 0.25*(e1->getIndex() + 1);
    glTexCoord2f(widthpct, 1.0);
    glVertex3d(v1.x, v1.y, v1.z);
    glTexCoord2f(widthpct, 0.0);
    glVertex3d(v0.x, v0.y, v0.z);
    glEnd();
    } // if smoothing
  } // if texturing
  else { // no textures, just draw each edge in white or needle color
    int microhack = 0;
    if (strstr(getObject()->getName(), "thread"))
      microhack = 1;
    for (int i = 0; i < numedges; i++) {
      glColor3d(1.0, 1.0, 1.0);
      // hack for suture/microsurgery application
      if (microhack && (i < RIGIDNEEDLE))
        glColor3d(0.5, 0.5, 0.5);
      edges[i]->draw();
    }
  }
}

void EdgeArray::drawover(double scale)
{
    // loop through all the edges
    for (int i = 0; i < numedges; i++) {
        edges[i]->drawover(scale);
    }
}

void EdgeArray::drawinitial()
{
    // loop through all the edges
    for (int i = 0; i < numedges; i++) {
        // draw the initial configuration in gray 
        edges[i]->drawInitial();
    }
}

void EdgeArray::drawlabels(float offset)
{
    // loop through all the edges
    for (int i = 0; i < numedges; i++) {
        // draw the labels
        edges[i]->drawlabel(offset);
    }
}

// this sets the spring constant for every edge to be value
void EdgeArray::setSpringConstants(double value, int restLengthScale)
{
  for (int i = 0; i < numedges; i++) {
    edges[i]->setSpringConstant(value, restLengthScale);
  }
}

// this sets the damp constant for every edge to be value
// NOTE: not using dampConstant usually, using velocityDampConstant in node
void EdgeArray::setDampConstants(double value)
{
    for (int i = 0; i < numedges; i++) {
        edges[i]->setDampConstant(value);
    }
}

// function to calculate total energy of the edges
double EdgeArray::getEnergy(void)
{
    int i;
    double total = 0.0;
    for (i = 0; i < numedges; i++) 
        total += edges[i]->getEnergy();
    
    return total;
}

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

void EdgeArray::DeleteEdge(Edge* e)
{
    if (debug) cerr << "EdgeArray::DeleteEdge(" << e << ")\n";
    
    // get index of edge
    int index = getIndex(e);
    
    // and delete it
    DeleteEdge(index);
}

void EdgeArray::EraseEdge(Edge* e)
{
    Edge *tempEdge = e;
    Node *node0 = e->getNodeLink(0);
    
    node0->deleteEdgeLink( e );
    e->getOtherNode(node0)->deleteEdgeLink( e );
    DeleteEdge( tempEdge );
}

void EdgeArray::DeleteFacelessEdges()
{
    int open_slot = 0;
    for (int i=0; i<numedges; i++) {
        if (edges[i]->numFaceLinks() <= 0) 
            delete edges[i];
        else edges[open_slot++] = edges[i];
    }
    numedges = open_slot;
}

void EdgeArray::DeleteAllEdges()
{
    // turns out all we have to do is reset the number
    numedges = 0;
}


void EdgeArray::SanityCheck()
{
    string s("");
    
    // check our stuff
    if (numedges < 0) s += "numedges < 0!";
    if (numedges > maxedges) s += "numedges > maxedges!";
    if (maxedges < 0) s += "maxedges < 0!";
    if (!edges) s += "edges is NULL!";
    if (!nodearray_p) s += "nodearray_p is NULL!";
    
    // if there's a problem, print it out
    //s << ends;
    if (s.size() > 0)
        cerr << "EdgeArray: INSANE: " << s.c_str() << endl;
    
    // check all the edges
    for (int i=0; i<numedges; i++)
        edges[i]->SanityCheck(this);
}

void EdgeArray::Cleanup()
{
    // check all the edges
    for (int i=0; i<numedges; i++) {
        Edge* e = getEdge(i);
        int s = e->SanityCheck(this);
        if (!s) EraseEdge(e);
    }
}

int EdgeArray::CompressDuplicateEdges()
{
    int num_dup_edges = 0;
    
    for (int i=0; i<getNumEdges(); i++) {
        Edge* e = getEdge(i);
        Node* n0 = e->getNodeLink(0);
        Node* n1 = e->getNodeLink(1);
        
        for (int j=0; j<getNumEdges(); j++) {
            Edge* other_e = getEdge(j);
            
            // if is a duplicate (between same nodes as edge e)
            if (other_e->hasNodeLink(n0) && other_e->hasNodeLink(n1)) {
                // move our faces to them
                // (note: if we have a really wacky geometry, then each edge could
                // have duplicate faces and we'd be moving too many faces over,
                // but then the mesh is so bad that this is okay to do anyway)
                for (int k=0; k<other_e->numFaceLinks(); k++) {
                    Face* f = other_e->getFaceLink(k);
                    if (!e->hasFaceLink(f)) e->addFaceLink(f);
                }
                
                // and delete it
                DeleteEdge(j);
                
                // and keep counter
                num_dup_edges++;
            }
        }
    }
    
    // and return how many we squashed out
    return(num_dup_edges);
}

---