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

Go to the documentation of this file.

// $Id: object.cpp,v 1.249 2006/05/30 20:25:19 craig Exp $ 
// $Copyright: (c)2001 National Biocomputation Center, Stanford University $


#if defined(_WIN32)                     // SYK 5/09/06 for VC++ 2005
#pragma warning (disable:4355)
#endif // _WIN32

#include "object.h"

#include <stdlib.h>     // for malloc, free, realloc, qsort
#include <stdio.h>
#include <assert.h>  // for assert functions...
#include <time.h>
#include <sys/timeb.h>

#ifdef _WIN32
#else
#include <strings.h>
#endif

#include "springCore.h"
#include "objectarray.h"
#include "tetraarray.h"
#include "facearray.h"
#include "edgearray.h"
#include "nodearray.h"

#include "node.h"
#include "edge.h"
#include "face.h"
#include "tetra.h"

#include "point3d.h"
#include "cyberware.h"
#include "rasterfont.h"
#include "geometryreplicator.h"
#include "perlinnoise.h"        // perlin noise function for fluid dynamics

#include "space.h"
#include "triangles.h"
#include "void.h"

#include "haptic_v1.h"
#include "haptic_v2.h"

#ifndef _WIN32
// If this is not known
#define _fileno fileno
#endif // _WIN32

#include "glui.h"


#ifdef __APPLE__
#include <OpenGL/gl.h>
#else
#include <GL/gl.h>
#endif // __APPLE__


// globals
const char* Object::rcsid = "@(#) $Id: object.cpp,v 1.249 2006/05/30 20:25:19 craig Exp $ $Copyright: (c)2001 National Biocomputation Center, Stanford University $";
int Object::debug = 0;

// default deformable material properties
static GLfloat default_deformable_mat_ambient[] = { 0.0f, 0.0f, 0.0f, 1.0f };
static GLfloat default_deformable_mat_diffuse[] = { 0.9f, 0.35f, 0.2f, 1.0f };
static GLfloat default_deformable_mat_specular[] = { 1.0f, 1.0f, 1.0f, 1.0f };
static GLfloat default_deformable_mat_shininess[] = { 128.0f };

// default deformable material properties
static GLfloat default_articulate_mat_ambient[] = { 0.0f, 0.0f, 0.0f, 1.0f };
static GLfloat default_articulate_mat_diffuse[] = { 0.9f, 0.35f, 0.2f, 1.0f };
static GLfloat default_articulate_mat_specular[] = { 1.0f, 1.0f, 1.0f, 1.0f };
static GLfloat default_articulate_mat_shininess[] = { 128.0f };

// default rigid body (bone) material properties
static GLfloat default_rigid_mat_ambient[] = { 0.0f, 0.0f, 0.0f, 1.0f };
static GLfloat default_rigid_mat_diffuse[] = { 0.8f, 0.65f, 0.5f, 1.0f };
static GLfloat default_rigid_mat_specular[] = { 1.0f, 1.0f, 1.0f, 1.0f };
static GLfloat default_rigid_mat_shininess[] = { 128.0f };

// default rest material properties
static GLfloat default_rest_mat_ambient[] = { 0.0f, 0.0f, 0.0f, 1.0f };
static GLfloat default_rest_mat_diffuse[] = { 0.5f, 0.5f, 0.5f, 1.0f };
static GLfloat default_rest_mat_specular[] = { 1.0f, 1.0f, 1.0f, 1.0f };
static GLfloat default_rest_mat_shininess[] = { 100.0 };

// default button dynamics material properties
static GLfloat default_button_mat_ambient[] = { 0.0f, 0.0f, 0.0f, 1.0f };
static GLfloat default_button_mat_diffuse[] = { 0.5f, 0.5f, 0.5f, 1.0f };
static GLfloat default_button_mat_specular[] = { 1.0f, 1.0f, 1.0f, 1.0f };
static GLfloat default_button_mat_shininess[] = { 100.0 };

Object::Object()
        : nodearray(this), 
                        edgearray(&nodearray), 
                        facearray(&nodearray,&edgearray),
                        tetraarray(&nodearray,&edgearray,&facearray)
{ 
    if (debug) cerr << "Object::constructor()\n";
    objectarray_p = NULL;
    type = nodes_only; 
    dynamics = deformable;
    behavior = none;
    nummethod = quasi_ordered;
    visible = 1;
    selected = 0;
    texture = NULL;
    nodearray.Init(&edgearray);
    strncpy(name, "UNINITIALIZED", NAME_LENGTH);
    sensor_p = NULL;    // not linked to any sensor
    linked_part = 0;    // entire object (class/color 0)
    tip_node = 0;
    drawover_amount = 0.0;
        being_modified = 0;
    
    IsAScreen   = 0;
    group_number = 0;
    
    closest_node = NULL;
    being_grabbed = 0;
        do_deform_only_on_grab = 0;
        do_preserve_volume = 0;
    sensor_offset = Point3D(0,0,0);     
    uses_viagra = 0;
    

        // Collision: general variables

    collision_extent = no_collisions;

        collision_force = penetration_depth_force;

        boundingBoxRoot = NULL;                 // used in BOUNDING_SPHERE_CD

        boundingSphereRoot = NULL;              // used in AABB_CD

        maxLevelCollisionCheck = (int)1e6;      // how many levels deep to check for collisions 

        ConstructCollisions();


    // labels stuff
    label_offset = 0.0;
    
    // graphics stuff
    dl_needs_refresh = 1;
    dl_id = -1;
    
    // initial params
    display_mode.mode = DisplayMode::smooth; 
    display_mode.draw_subobjcolormode = 0;
    last_display_mode = display_mode;
                
        // initialize default fluid height
        height = 0;

        // initialize hinges
    hinges.Init();

    // initialize rigid object hinges
        rigidHinge = NULL;      // default for no hinge!

    // set default material properties, based on dynamics
    set_default_material_properties();
    timer.Reset();
        
        // Initialize pointer to GLUI attribute dialog.
        current_attribute_menu = NULL;
}

Object::~Object()
{ 
    if (debug) cerr << "Object[" << name << "]::destructor()\n";
    strncpy(name, "DESTROYED", NAME_LENGTH);
    

        // deinit collision stuff
    DestructCollisions();



        // destroy bounding boxes if any


    // if we're linked to a sensor, unlink us
    if (sensor_p) sensor_p->LinkToObject(NULL);
}

void Object::Init(ObjectArray* objectarray_in, char* name_in)
{
    objectarray_p = objectarray_in;
    if (!name_in) name_in = "UNNAMED";
    strncpy(name, name_in, NAME_LENGTH);
        strncpy(path, "", PATH_LENGTH);
    objextrusions = 0;
}

void Object::set_default_material_properties()
// set default material properties based on the value of dynamics
{
    if (debug) 
        cerr << "Object[" << getName() << "]:set_default_material_properties()\n";
    
    switch (dynamics) {
                case rest:
                        memcpy(mat_ambient, default_rest_mat_ambient, 4*sizeof(float));
                        memcpy(mat_diffuse, default_rest_mat_diffuse, 4*sizeof(float));
                        memcpy(mat_specular, default_rest_mat_specular, 4*sizeof(float));
                        memcpy(mat_shininess, default_rest_mat_shininess, 1*sizeof(float));
                        break;
                case rigid:
                        memcpy(mat_ambient, default_rigid_mat_ambient, 4*sizeof(float));
                        memcpy(mat_diffuse, default_rigid_mat_diffuse, 4*sizeof(float));
                        memcpy(mat_specular, default_rigid_mat_specular, 4*sizeof(float));
                        memcpy(mat_shininess, default_rigid_mat_shininess, 1*sizeof(float));
                        break;
                case deformable:
                        memcpy(mat_ambient, default_deformable_mat_ambient, 4*sizeof(float));
                        memcpy(mat_diffuse, default_deformable_mat_diffuse, 4*sizeof(float));
                        memcpy(mat_specular, default_deformable_mat_specular, 4*sizeof(float));
                        memcpy(mat_shininess, default_deformable_mat_shininess, 1*sizeof(float));
                        break;
                case articulate:
                        memcpy(mat_ambient, default_articulate_mat_ambient, 4*sizeof(float));
                        memcpy(mat_diffuse, default_articulate_mat_diffuse, 4*sizeof(float));
                        memcpy(mat_specular, default_articulate_mat_specular, 4*sizeof(float));
                        memcpy(mat_shininess, default_articulate_mat_shininess, 1*sizeof(float));
                        break;
                case button:
                        memcpy(mat_ambient, default_button_mat_ambient, 4*sizeof(float));
                        memcpy(mat_diffuse, default_button_mat_diffuse, 4*sizeof(float));
                        memcpy(mat_specular, default_button_mat_specular, 4*sizeof(float));
                        memcpy(mat_shininess, default_button_mat_shininess, 1*sizeof(float));
                        break;
    } // end switch
    
    // by default, make the back faces the same as the front
    memcpy(back_mat_ambient, mat_ambient, 4*sizeof(float));
    memcpy(back_mat_diffuse, mat_diffuse, 4*sizeof(float));
    memcpy(back_mat_specular, mat_specular, 4*sizeof(float));
    memcpy(back_mat_shininess, mat_shininess, 1*sizeof(float));
}

// Amira format
int Object::ReadAmira(FILE* fp)
{
 if (debug) cerr << "Object[" << name << "]::ReadAmira(" 
                                                                 << int(fp) << ")\n";

        dynamics = rigid;
        set_default_material_properties();
        type = tetras_only;
        
        int num_nodes = 0;
        int num_tetra = 0;
        
        char s[MAXLINE];
  myfgets(s, MAXLINE, fp);
  int num_args = sscanf(s, "%d %d", &num_nodes, &num_tetra);
         
        if (num_args < 2) {
                cerr << "bad file" << endl;
                return 0;
        }
  // process the node lines (will also compute texture coords if needed)
  int status;
  if ((status = nodearray.getAmiraNodeData(fp,num_nodes)) == 0) 
    return(status);

  // process the face lines, this creates both faces and edges
  status = tetraarray.getAmiraTetraData(fp,num_tetra);

  // recompute normals (since surface is changing)
//  tetraarray.computeNormals();
//  nodearray.interpolateTetraNormals();
        
  // try to read in a texture (use name of object.rgb)
  setTexture(NULL);

  // return value
  return(status);
}

int Object::ReadMesh(FILE* fp)
{
    if (debug) cerr << "Object[" << name << "]::ReadMesh(" << int(fp) << ")\n";
    
    { // deal with header lines
        char s[MAXLINE];
        myfgets(s, MAXLINE, fp);    // obj# grp#
        int obj_num, group_num, dummy;
        int num_args = sscanf(s, "%d %d %d", &obj_num, &group_num, &dummy);
        if (num_args == 2) { // newer style mesh format
            if (debug) cerr << "newer style mesh file\n";
            if (group_num == 3) dynamics = rest;    // bone group number/color
            else dynamics = deformable; 
            
            myfgets(s, MAXLINE, fp);    // obj center
            myfgets(s, MAXLINE, fp);    // crap crap more_crap
        } else {        // very old style mesh format
            if (debug) cerr << "older style mesh file\n";
            dynamics = deformable; 
        }
        
        // setup material properties appropriately
        set_default_material_properties();
    }
    
    // set initial type to be just nodes
    type = nodes_only;
    
    // process the node lines (will also compute texture coords if needed)
    int status;
    if ((status = nodearray.getMeshNodeData(fp)) == 0) 
        return(status);
    
    // process the face lines, this creates both faces and edges
    status = facearray.getMeshFaceData(fp);
    
    // if we made it this far, assume we have faces
    type = faces_only;
    
#define RECALC_NORMALS_ON_READ
#ifdef RECALC_NORMALS_ON_READ
    // recompute normals (since surface is changing)
    facearray.computeNormals();
    
    // average face normals to get per-vertex normals
    nodearray.interpolateNormals();
#endif
    
    // try to read in a texture (use name of object.rgb)
    setTexture(NULL);
    
    // set some better defaults
    if (dynamics == deformable) {
        setNumMethod(Object::numerical_method(5));  // ordered quasi
        getEdgeArray()->setSpringConstants(100.0);  // stiff
        getNodeArray()->setVelocityDampConstants(100.0);
    }
    
    // return value
    return(status);
}

int Object::ReadSMF(FILE* fp)
{
    if (debug) cerr << "Object[" << name << "]::ReadSMF(" << int(fp) << ")\n";
    
    // assume bone group number/color
    dynamics = rigid;    
    
    // setup material properties appropriately
    set_default_material_properties();
    
    // set initial type to be just nodes
    type = faces_only;
    
    // process the node lines
    int status;
    if ((status = nodearray.getSMFNodeData(fp)) == 0) 
        return(status);
    
    // process the face lines, this creates both faces and edges
    status = facearray.getSMFFaceData(fp);
    
    // compute normals
    facearray.computeNormals();
    
    // average face normals to get per-vertex normals
    nodearray.interpolateNormals();
    
    // and compute the texture coords (assume cylindrical)
    nodearray.computeTextureCoords(Node::cylindrical);
    
    // try to read in a texture (use name of object.rgb)
    setTexture(NULL);
    
    // return value
    return(status);
}

// Alias/Wavefront OBJ format
int Object::ReadOBJ(FILE* fp)
{
    if (debug) cerr << "Object[" << name << "]::ReadOBJ(" << int(fp) << ")\n";
    
    // assume bone group number/color
    dynamics = rest;    
    
    // setup material properties appropriately
    set_default_material_properties();
    
    // set initial type to be just nodes
    type = faces_only;
    
    // process the node lines
    int status = nodearray.getOBJNodeData(fp);
    if (status  == 0) {

                // error condition
        return(status);

        }

        else if (status == 1) {

                // no texture coords defined, calculate cylindrical

                nodearray.computeTextureCoords(Node::cylindrical);

        }
    
    // process the face lines, this creates both faces and edges
    status = facearray.getOBJFaceData(fp);
    
    // compute normals
    facearray.computeNormals();
    
    // average face normals to get per-vertex normals
    nodearray.interpolateNormals();
    
    // try to read in a texture (use name of object.rgb)
    setTexture(NULL);
    
    // return value
    return(status);
}

// Reads an SOBJ file (=Spring OBJ) - The file contains commands
// to create the basic objects which are already defined in
// spring such as sphere and others.
int Object::ReadSOBJ(FILE* fp)
{
        // the status which is returned later
        int status = 0;

    // make a local buffer for reading/parsing
    char buffer[256];
        char object_name[80]; 
    
    // get the first line
    fgets(buffer, 256, fp);
        sscanf(buffer, "%s", object_name);

        // check what we should create
        if (!strcmp(object_name, "sphere")) {
                // parse attributes
                float layers, radius, inner;
        sscanf(buffer, "%*s %f %f %f", &layers, &radius, &inner);
                
                // and create sphere
                status = CreateSphere((int)layers, radius, inner);
        }
        else if (!strcmp(object_name, "plane")) {
                // parse attributes
                float x_length, y_length;
                int x_steps, y_steps;
                float x_offset, y_offset, z_offset;
                char orientation;
        sscanf(buffer, "%*s %f %f %d %d %f %f %f %c", &x_length, &y_length, &x_steps, &y_steps, 
                                                                  &x_offset, &y_offset, &z_offset, &orientation);
                // and create plane
                CreatePlane(x_length, y_length, x_steps, y_steps, 
                            Point3D(x_offset, y_offset, z_offset), orientation);
        }

        // return value
        return(status);
}

int Object::ReadFRD(FILE* fp)
{
    if (debug) cerr << "Object[" << name << "]::ReadFRD(" << int(fp) << ")\n";
    
    char s[MAXLINE];
    int dyn;
    float Rred, Ggreen, Bblue, Aalpha;
    
    // setup material properties appropriately
    set_default_material_properties();
    
    // set initial type to be just nodes
    type = faces_only;
    
    // process the node lines
    int status;
    if ((status = nodearray.getSMFNodeData(fp)) == 0) 
        return(status);
    
    // process the edge lines
    status = edgearray.getFRDEdgeData(fp);
    
    // process the face lines, this creates faces (and edges if needed)
    status = facearray.getSMFFaceData(fp);
    
    // process extrusion, dynamics and color data
    myfgets(s, MAXLINE, fp);
    if (s[0] == 'o')
        sscanf(s, "%*c %d", &objextrusions);
    else {
        printf("Bad read on objextrusions. objextrusions will be set to 0\n");  
        objextrusions = 0;;
    }
    
    myfgets(s, MAXLINE, fp);
    if (s[0] == 'd'){    
        sscanf(s, "%*c %d", &dyn);
        setDynamics(dynamics_type(dyn));
    }
    else {
        printf("Bad read on dynamics. Dynamics will be set to rest\n");  
        dynamics = rest;
    }
    
    myfgets(s, MAXLINE, fp);
    if (s[0] == 'r')
        sscanf(s, "%*c %e", &Rred);
    else {
        printf("Bad read on red color. Red color will be set to 0.5\n");  
        Rred = 0.5;
    }
    myfgets(s, MAXLINE, fp);
    if (s[0] == 'g')
        sscanf(s, "%*c %e", &Ggreen);
    else {
        printf("Bad read on green color. Green color will be set to 0.5\n");  
        Ggreen = 0.5;
    }
    myfgets(s, MAXLINE, fp);
    if (s[0] == 'b')
        sscanf(s, "%*c %e", &Bblue);
    else {
        printf("Bad read on blue color. Blue color will be set to 0.5\n");  
        Bblue = 0.5;
    }
    myfgets(s, MAXLINE, fp);
    if (s[0] == 'a')
        sscanf(s, "%*c %e", &Aalpha);
    else {
        printf("Bad read on alpha. Alpha will be set to 0.5\n");  
        Aalpha = 0.5;
    }
    
    setDiffuseColor(Rred, Ggreen, Bblue, Aalpha);
    
    myfgets(s, MAXLINE, fp);
    if (s[0] == 'n')
        sscanf(s, "%*c %d", &Nb_face_surf);
    else {
        printf("Bad read on the # of faces on the original mesh.\nIt will be set to the actual # of faces of this object\n");  
        Nb_face_surf = facearray.getNumFaces();;
    }
    
    double tempmass, tempspring, tempdamp;
    myfgets(s, MAXLINE, fp);
    if (s[0] == 'm')
        sscanf(s, "%*c %lf", &tempmass);
    else {
        printf("Bad read on the mass constant on the original mesh.\nIt will be set to 1\n");
        nodearray.setMasses(1.0);
    }
    nodearray.setMasses(tempmass);
    
    myfgets(s, MAXLINE, fp);
    if (s[0] == 's')
        sscanf(s, "%*c %lf", &tempspring);
    else {
        printf("Bad read on the spring constant on the original mesh.\nIt will be set to 1\n");
        edgearray.setSpringConstants(1.0);
    }
    edgearray.setSpringConstants(tempspring);
    
    myfgets(s, MAXLINE, fp);
    if (s[0] == 'p')
        sscanf(s, "%*c %lf", &tempdamp);
    else {
        printf("Bad read on the damp constant on the original mesh.\nIt will be set to 1\n");
        nodearray.setVelocityDampConstants(1.0);
    }
    nodearray.setVelocityDampConstants(tempdamp);
    
    
    // compute normals
    facearray.computeNormals();
    
    // average face normals to get per-vertex normals
    nodearray.interpolateNormals();
    
    // and compute the texture coords (assume cylindrical)
    nodearray.computeTextureCoords(Node::cylindrical);
    
    // try to read in a texture (use name of object.rgb)
    setTexture(NULL);
    
    // return value
    return(status);
}

char* get_next_token(char* buffer, int bufsize, FILE* fp)
{
    char* delimiters = " \t\n\r,";
    char* ptr = strtok(NULL, delimiters);
    
    while (!ptr && !feof(fp)) { // ran out of stuff on this line
        char* b = myfgets(buffer, bufsize, fp); //  get a new line
        ptr = strtok(b, delimiters);
    }
    return(ptr);
}

// VRML format
// Note: VRML files may have many subobjects inside, so this function may
// recursively create many objects if create_subobjects is set to true!
// (otherwise it will create one large object, but will only have one texture
// defined)
int Object::ReadVRML(FILE* fp, int create_subobjects)
{
  if (debug) cerr << "Object[" << name << "]::ReadVRML(" << int(fp) << ")\n";

  // make a line buffer
  const int bufsize = 256;
  static char buffer[bufsize]; 
  strtok(buffer, "\0"); // set up for later parsing
        
  // read the initial line of the file
  fgets(buffer, bufsize, fp);   //  get a raw first line
  strtok(buffer, "\n");         // prime strtok for later parsing
  int is_vrml2 = 0;
  if (strstr(buffer, "2.")) is_vrml2 = 1;
  if (debug) {
    cerr << "file type buffer is [" << buffer << "]\n";
    cerr << "File format is : ";
    if (is_vrml2) cerr << "VRML2.0\n";
    else cerr << "VRML1.0\n";
  }

  // create an initial transformation matrix
  Matrix initial_xform; // initialized to identity

  // read (maybe recursively) the objects in the VRML file
  int status = ReadVRMLObject(fp, initial_xform, create_subobjects, is_vrml2);

  // and return
  return(status);
}

// function to (perhaps recursively) read a VRML object
int Object::ReadVRMLObject(FILE* fp, Matrix original_xform, int create_subobjects, int is_vrml2)
{
  if (debug) cerr << "Object[" << name << "]::ReadVRMLObject(" << int(fp) << ", xform, " << create_subobjects << ", " << is_vrml2 << ")\n";
    
  int status = 1;       // assume okay
    
  // assume bone group number/color
  dynamics = rest;    
    
  // setup material properties appropriately
  set_default_material_properties();
    
  // set initial type to be faces
  type = faces_only;
    
  // note: we used to do allocateNodes, Edges, and Faces here, but now we read the items into ReallocableArrays
  // first, then allocate the right amount of space and copy them over- saves a ton of space and is much smarter...
    
  // make a line buffer
  const int bufsize = 256;
  static char buffer[bufsize]; 
  strtok(buffer, "\0"); // set up for later parsing
    
  // get some locals for parsing
  int nesting = 1;
  int last_block_of_coords_index = 0;
  int cur_part_id = -1;
  int num_defs = 0;
  ReallocableArray<float> texture_u, texture_v;
    
  // transform values
  Matrix local_xform;
  int local_trans_nest_level = 0, local_rot_nest_level = 0, 
    local_scale_nest_level = 0;
    
  // start parsing
  while (!feof(fp) && (nesting > 0)) {

    // get the next token
    char* ptr = get_next_token(buffer, bufsize, fp);
    if (!ptr) break;    // hit eof
        
    { // parse it
      //cerr << getName() << "-token: " << ptr << ", nesting = " << nesting << ", recurse=" << create_subobjects << endl;
      //getchar();
            
      // process that token
      if (strstr(ptr, "{") || strstr(ptr,"[")) nesting++;
      else if (strstr(ptr, "}") || strstr(ptr, "]")) nesting--;
      else if (!strcmp(ptr, "DEF")) {   // defines a new object
        ptr = get_next_token(buffer, bufsize, fp);      // prime for name
        setName(ptr);
        num_defs++;
        cerr << num_defs << ": DEF of [" << ptr << "]\n";
        if (create_subobjects) {
                                        
          // create a composite matrix of the original and current transforms
          Matrix composite_xform;
          composite_xform.Multiply(original_xform);
          composite_xform.Multiply(local_xform);
          //cerr << "original matrix: " << endl; original_xform.Print();
          //cerr << "local matrix: " << endl; local_xform.Print();
          //cerr << "composite matrix: " << endl; composite_xform.Print();

          // spawn a new object and tell it to go read itself
          Object* new_obj = objectarray_p->addObject(ptr);
          new_obj->ReadVRMLObject(fp, composite_xform, ++create_subobjects, is_vrml2);
        }
      }
      else if (!strcmp(ptr, "Coordinate3")) {
                                
        // now doing a new part/block of vertices
        cur_part_id++;
                
        // loop, reading in the vertices
        last_block_of_coords_index = nodearray.getNumNodes();
        ReallocableArray<Point3D> local_pts;
        int subnesting = 0;
        do {
          ptr = get_next_token(buffer, bufsize, fp);    // prime for token
          if (!ptr) break;
                                        
          if (strstr(ptr, "{") || strstr(ptr, "[")) subnesting++;
          else if (strstr(ptr, "}") || strstr(ptr, "]")) subnesting--;
          else if (!strcmp(ptr, "point")) continue;
          else {        // must be points
            float x,y,z;
            x = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  // get y
            y = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  // get z
            z = atof(ptr);
            Point3D p(x,y,z);
                                                
            // and save it
            local_pts.Append(p);
                                                
            // if last point had a } or ] at the end, have to catch that, too
            if (strstr(ptr, "}") || strstr(ptr, "]")) subnesting--;
          }
                    
          //cerr << "node is " << p << endl;
        } while (subnesting > 0);

                                // allocate space in this object for nodes and edges
        int num_read = local_pts.NumElements();
        cerr << "...read " << num_read << " nodes\n";
        nodearray.allocateNodes((int)(num_read * 1.10));
        edgearray.allocateEdges((int)(num_read * 2 * 1.10));

                                // create a composite matrix of the original and current transforms
        Matrix composite_xform;
        composite_xform.Multiply(original_xform);
        composite_xform.Multiply(local_xform);
                                //cerr << "original matrix: " << endl; original_xform.Print();
                                //cerr << "local matrix: " << endl; local_xform.Print();

                                // and copy them over into the object, transforming as we go
        for (int i=0; i<local_pts.NumElements(); i++) {
          Point3D p = local_pts[i];
          p.Xform(composite_xform.m);

          // and add it
          int index = nodearray.addNode(p);
          Node* node = nodearray.getNode(index);
          node->setPartId(cur_part_id);
        }
      }
      else if (!strcmp(ptr, "coordIndex")) {
        ReallocableArray<iPoint3D> local_faces;
        int subnesting = 0, num_read = 0;
        do {
          ptr = get_next_token(buffer, bufsize, fp);    // prime for token
                                        
          if (strstr(ptr, "{") || strstr(ptr, "[")) { subnesting++; continue; }
          if (strstr(ptr, "}") || strstr(ptr, "]")) { subnesting--; continue; }
                                        
          int i1, i2, i3;
          i1 = atoi(ptr);
          ptr = get_next_token(buffer, bufsize, fp);    // get i2
          i2 = atoi(ptr);
          ptr = get_next_token(buffer, bufsize, fp);    // get i3
          i3 = atoi(ptr);
                                        
          if ((i1 < 0) || (i2 < 0)) {
            cerr << " face[ " << facearray.getNumFaces() << "] has indices " << i1 << ", " << i2 << endl;
          }
          else { // okay values
            int is_two_params = (i3 == -1);     // if only two params, set flag
                                                
            // make relative to current block of vertices for this subobject
            i1 += last_block_of_coords_index;
            i2 += last_block_of_coords_index;
            i3 += last_block_of_coords_index;
                                                
            int nn = nodearray.getNumNodes();
            if ((i1 >= nn) || (i2 >= nn) || (i3 >= nn)) {
              cerr << " face[ " << facearray.getNumFaces() << "] has indices beyond bounds (" << nn << "): " 
                   << i1 << ", " << i2 << ", " << i3 << endl;
            }
            else if (!is_two_params) {                                                                  
              ptr = get_next_token(buffer, bufsize, fp);        // flush -1
              iPoint3D temp(i1, i2, i3);
              local_faces.Append(temp);
              num_read++;
            }
            else {  // only 2 indices- is edge-based object
              type = edges_only;
              display_mode.mode = DisplayMode::wireframe;
              edgearray.addEdge(i1,i2);
              cerr << "edge is " << i1 << " " << i2 << endl;
              num_read++;
            }
          }             
          //cerr << "face is " << i1 << " " << i2 << " " << i3 << endl;
        } while (subnesting > 0);
        cerr << "...read " << num_read << " faces\n";

                                // allocate space for faces
        facearray.allocateFaces((int)(num_read * 1.10));

                                // and copy over
        for (int i=0; i<local_faces.NumElements(); i++) {
          iPoint3D temp = local_faces[i];
          facearray.addFace(temp.x, temp.y, temp.z);    // actually i1, i2, i3
        }
      }
      else if (!strcmp(ptr, "TextureCoordinate2")) {
                                              
        // loop, reading in the vertices
        int subnesting = 0;
        do {
          ptr = get_next_token(buffer, bufsize, fp);    // prime for token
          if (!ptr) break;

          if (strstr(ptr, "{") || strstr(ptr, "[")) subnesting++;
          else if (strstr(ptr, "}") || strstr(ptr, "]")) subnesting--;
          else if (!strcmp(ptr, "point")) continue;
          else {        // must be the points
            float u, v;
            u = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  // get u
            v = atof(ptr);

            // save it
            texture_u.Append(u); texture_v.Append(v);
            //cerr << "texindex[" << texture_u.NumElements()-1 << "] = (" << u << ", " << v << ")\n";
                                                
            // if last point had a } or ] at the end, have to catch that, too
            if (strstr(ptr, "}") || strstr(ptr, "]")) subnesting--;
          }
                    
          //cerr << "node is " << p << endl;
        } while (subnesting > 0);
        cerr << "...read " << texture_u.NumElements() << " texture coords\n";
      }            
      else if (!strcmp(ptr, "textureIndex")) {  // directly gives the u,v for each vertex
        ptr = get_next_token(buffer, bufsize, fp);      // flush the '['
        // loop, reading in the vertices
        int index = last_block_of_coords_index;
        ptr = get_next_token(buffer, bufsize, fp);      // prime for u
        while (!strstr(ptr, "}") && !strstr(ptr, "]")) {
          float u, v;
          u = atof(ptr);
          ptr = get_next_token(buffer, bufsize, fp);    // get v
          v = atof(ptr);
          ptr = get_next_token(buffer, bufsize, fp);    // get -1
                    
          // set the texture coord here
          Node* n = nodearray.getNode(index);
          if (n) {
            n->setTextureCoords(u, v);
            cerr << "setting node " << index << " texture coords to " << u << ", " << v << endl;
          }

          // prime for next time
          index++;
          ptr = get_next_token(buffer, bufsize, fp);    // prime for u again
        }
      }
      else if (!strcmp(ptr, "textureCoordIndex")) {     // indirectly gives texture u,v- must index through face
        // loop, reading in the texture indices for each face
        int face_index = 0, subnesting = 0;
        do {
          ptr = get_next_token(buffer, bufsize, fp);    // get a token

          if (strstr(ptr, "{") || strstr(ptr, "[")) { subnesting++; continue; }
          if (strstr(ptr, "}") || strstr(ptr, "]")) { subnesting--; continue; }
 
          // get index to face
          Face* f = facearray.getFace(face_index);
          if (!f) { cerr << "textureCoordIndex: bad face index=" << face_index << endl; continue; }
                                        
          // loop for each vertex
          for (int i=0; i<3; i++) {
            int index = atoi(ptr);      // parse the texture array index for this node
            if ((index < 0) || (index >= texture_u.NumElements())) { 
              cerr << "textureCoordIndex: index " << index << " out of bounds\n"; 
              continue; 
            }
            float u = texture_u[index], v = texture_v[index];   // look up the values
            Node* n = f->getNodeLink(i);
            if (!n) continue;
            n->setTextureCoords(u, v);
            //cerr << "setting face#" << face_index << ".node" << i << " (node#" << nodearray.getIndex(n)
            //<< ") texture coords to texindex[" << index << "] = (" << u << ", " << v << ")\n";;
            ptr = get_next_token(buffer, bufsize, fp);  // get next index (or flush the -1 at end)
            if (strstr(ptr, "}") || strstr(ptr, "]")) { subnesting--; continue; }       // for if concatenated at end
          }

          face_index++;
        } while (subnesting > 0);
      }
      else if (!strcmp(ptr, "Normal")) {
        int subnesting = 0;
        do {
          ptr = get_next_token(buffer, bufsize, fp);    // prime for token
          if (strstr(ptr, "{") || strstr(ptr, "[")) subnesting++;
          else if (strstr(ptr, "}") || strstr(ptr, "]")) subnesting--;
          // pretty much ignore it- we'll calculate our own anyway
        } while (subnesting > 0);
      }
      else if (!strcmp(ptr, "normalIndex")) {
        do {
          ptr = get_next_token(buffer, bufsize, fp);    // get a token
        } while (!strstr(ptr, "]"));
      }
      else if (!strcmp(ptr, "MaterialBinding")) {
        do {
          ptr = get_next_token(buffer, bufsize, fp);    // get a token
        } while (!strstr(ptr, "}"));
      }
      else if (!strcmp(ptr, "Material")) {
        ptr = get_next_token(buffer, bufsize, fp);      // flush the '{'
                                
        // loop, reading in the attributes
        do {
          ptr = get_next_token(buffer, bufsize, fp);    // prime for attribute
                                        
          //cerr << "Material: " << ptr << endl;
          if (!strcmp(ptr, "ambientColor")) {
            double r,g,b;
            ptr = get_next_token(buffer, bufsize, fp);  r = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  g = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  b = atof(ptr);
            setAmbientColor(r,g,b,1.0);
          }
          else if (!strcmp(ptr, "diffuseColor")) {
            double r,g,b;
            ptr = get_next_token(buffer, bufsize, fp);  r = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  g = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  b = atof(ptr);
            setDiffuseColor(r,g,b,1.0);
          }
          else if (!strcmp(ptr, "emissiveColor")) {
            double r,g,b;
            ptr = get_next_token(buffer, bufsize, fp);  r = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  g = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  b = atof(ptr);
            // we don't really have emissive, so call it diffuse
            setDiffuseColor(r,g,b,1.0);
          }
          else if (!strcmp(ptr, "specularColor")) {
            double r,g,b;
            ptr = get_next_token(buffer, bufsize, fp);  r = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  g = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  b = atof(ptr);
            setSpecularColor(r,g,b,1.0);
          }
          else if (!strcmp(ptr, "shininess")) {
            double s;
            ptr = get_next_token(buffer, bufsize, fp);  s = atof(ptr);
            if (s <= 1.0) s *= 128;
            setShininess(s);
          }
          else if (!strcmp(ptr, "transparency")) {      // 1 means fully opaque
            // in VRML1.0, a value of 1 meant fully opaque
            // in VRML2.0, a value of 1 meant fully transparent
            // go figure
            double a;
            ptr = get_next_token(buffer, bufsize, fp);  a = atof(ptr);
            //if (is_vrml2) a = 1.0 - a;
            a = 1.0 - a;        // seems to be reversed in VRML1.0 as well?
            float* cur_diffuse = getDiffuseColor();
            setDiffuseColor(cur_diffuse[0], cur_diffuse[1], cur_diffuse[2],a);
          }
          else if (debug)
            cerr << "Material: UNKNOWN ATTRIBUTE: " << ptr << endl;
        } while (!strstr(ptr, "}"));
                                        
        // already got the '}'
      }  // end Material
      else if (!strcmp(ptr, "Texture2")) {
        ptr = get_next_token(buffer, bufsize, fp);      // flush the '{'
        ptr = get_next_token(buffer, bufsize, fp);      // flush the "filename" token

        ptr = get_next_token(buffer, bufsize, fp);      // get the filename
        if (*ptr == '\"') { ptr++; ptr[strlen(ptr)-1] = '\0'; } // get rid of the '"' chars
                                
                                // open the texture file
        if (!getTexture()) {
          cerr << getName() << ": reading texture from file [" << ptr << "]\n";
          debug = 1;
          setTexture(ptr);
          debug = 0;
        }

        if (!getTexture())
          cerr << getName() << ": could not read texture\n";

        ptr = get_next_token(buffer, bufsize, fp);      // flush the '}'
      }  // end Texture2
      else if (strstr(ptr, "Texture2Transform")) {
        do {
          ptr = get_next_token(buffer, bufsize, fp);    // get a token
        } while (!strstr(ptr, "}"));
      }
      else if (!strcmp(ptr, "Translation")) {
        local_trans_nest_level = nesting;
        ptr = get_next_token(buffer, bufsize, fp);      // flush the '{'
                
        // loop, reading in the attributes
        ptr = get_next_token(buffer, bufsize, fp);      // prime for attribute
        while (!strstr(ptr, "}") && !strstr(ptr, "]")) {
          if (!strcmp(ptr, "translation")) {
            double x,y,z;
            ptr = get_next_token(buffer, bufsize, fp);  x = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  y = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  z = atof(ptr);
            local_xform.m[3][0] = x; local_xform.m[3][1] = y; local_xform.m[3][2] = z;
          }
          ptr = get_next_token(buffer, bufsize, fp);  // prime for token again
        }

                                //ptr = get_next_token(buffer, bufsize, fp);    // flush the '}'
      }  // end Translation
      else if (!strcmp(ptr, "Rotation")) {
        local_rot_nest_level = nesting;
        ptr = get_next_token(buffer, bufsize, fp);      // flush the '{'
                
        // loop, reading in the attributes
        ptr = get_next_token(buffer, bufsize, fp);      // prime for attribute
        while (!strstr(ptr, "}") && !strstr(ptr, "]")) {
          if (!strcmp(ptr, "rotation")) {
            double x,y,z,t;
            ptr = get_next_token(buffer, bufsize, fp);  x = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  y = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  z = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  t = atof(ptr);
            // heck if I know what to do about twist
            //local_xform = Point3D(x,y,z); 
            // PUT IN THE CODE TO CONVERT THE DIRECTION+TWIST
            // TO A ROTATION MATRIX AND STORE IT IN local_xform
          }
          ptr = get_next_token(buffer, bufsize, fp);  // prime for token again
        }

                                //ptr = get_next_token(buffer, bufsize, fp);    // flush the '}'
      }  // end Rotation
      else if (!strcmp(ptr, "Scaling")) {
        local_scale_nest_level = nesting;
        ptr = get_next_token(buffer, bufsize, fp);      // flush the '{'
                
        // loop, reading in the attributes
        ptr = get_next_token(buffer, bufsize, fp);      // prime for attribute
        while (!strstr(ptr, "}") && !strstr(ptr, "]")) {
          if (!strcmp(ptr, "scaling")) {
            double x,y,z;
            ptr = get_next_token(buffer, bufsize, fp);  x = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  y = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  z = atof(ptr);
            local_xform.m[0][0] = x; local_xform.m[1][1] = y; local_xform.m[2][2] = z; 
          }
          ptr = get_next_token(buffer, bufsize, fp);  // prime for token again
        }

        ptr = get_next_token(buffer, bufsize, fp);      // flush the '}'
      }  // end Scaling
      else if (!strcmp(ptr, "MatrixTransform")) {
        local_rot_nest_level = nesting;
        ptr = get_next_token(buffer, bufsize, fp);      // flush the '{'
                
        // loop, reading in the attributes
        ptr = get_next_token(buffer, bufsize, fp);      // prime for attribute
        while (!strstr(ptr, "}")) {
          if (!strcmp(ptr, "matrix")) {
            double m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, m41, m42, m43, m44;
                        
            ptr = get_next_token(buffer, bufsize, fp);  m11 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m12 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m13 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m14 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m21 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m22 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m23 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m24 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m31 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m32 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m33 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m34 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m41 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m42 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m43 = atof(ptr);
            ptr = get_next_token(buffer, bufsize, fp);  m44 = atof(ptr);
                                                
            // assign it
            local_xform.m[0][0] = m11; local_xform.m[0][1] = m12; local_xform.m[0][2] = m13;
            local_xform.m[1][0] = m21; local_xform.m[1][1] = m22; local_xform.m[1][2] = m23;
            local_xform.m[2][0] = m31; local_xform.m[2][1] = m32; local_xform.m[2][2] = m33;
            local_xform.m[3][0] = m41; local_xform.m[3][1] = m42; local_xform.m[3][2] = m43; 
          }
          ptr = get_next_token(buffer, bufsize, fp);  // prime for token again
        }

                                // already got the '}'
      }  // end Matrix
      else {    // unknown token
        if (debug) cerr << "UNKNOWN TOKEN: " << ptr << endl;
        // flush to next one?
      }
            
      // reset local transformations if go out of scope
      if (nesting < local_rot_nest_level) local_xform.Init();
            
    } // parse it
  }  // end while
  
  // debugging
  cerr << *this << ", " << cur_part_id << " subparts\n";
  cerr << "final nesting level was " << nesting << endl;
  
  // compute stats
  nodearray.computeStats();
  
  // compute normals
  facearray.computeNormals();
  
  // average face normals to get per-vertex normals
  nodearray.interpolateNormals();

  // assume they have read in their own texture coords and texture if they wanted it
  
  // return value
  return(status);
}

GSPEC *gs = NULL;                       /* file header buffer base pointer */
int vertex_is_void(Point3D p)
{
    if ((p.x == -1) || (p.y == -1) || (p.z == -1)) return(1);
    else return(0);
}

// for now, essentially a mesh file, but later will do convert itself
// Note: Must currently convert texture image file to power of 2 dimensions
//              (assuming 512x512)
int Object::ReadCyberware(FILE* fp)
{
    if (debug) 
        cerr << "Object[" << name << "]::ReadCyberware(" << fp << ")\n";
    
    // blow away any old data in this object
    nodearray.reset();
    edgearray.reset();
    facearray.reset();
    
    // read the cyberware header
    if ((gs = cyread(NULL, _fileno(fp))) == NULL) 
        return(0);
    
    // save important info from the data
    int is_cartesian = ((gs->flags & FLAG_CARTESIAN) > 0); // else cylindrical
    int lg_size = gs->lgmax - gs->lgmin + 1;
    int lt_size = gs->ltmax - gs->ltmin + 1;
    
    // display image info for verification
    {
        fprintf(stderr, "image: %s ", gs->name);
        fprintf(stderr, "hdrsiz=%d ltincr=%ld\n", sizeof(GSPEC), gs->ltincr);
        fprintf(stderr, "ltrange=%d:%d lgrange=%d:%d\n",
            gs->ltmin, gs->ltmax, gs->lgmin, gs->lgmax);
        fprintf(stderr, "coord type=%d (0=cyl, 1=cartesian)\n", is_cartesian);
    }
    
    // allocate space
    int projected_num_nodes = lg_size * lt_size;
    int projected_num_edges = 3 * projected_num_nodes;
    int projected_num_faces = 2 * projected_num_nodes;
    nodearray.allocateNodes(projected_num_nodes);
    edgearray.allocateEdges(projected_num_edges);
    facearray.allocateFaces(projected_num_faces);
    
    // for each longitude
    long theta = gs->lgmin * gs->lgincr;
    int num_good_nodes = 0, num_void_nodes = 0;
    for (int lg = gs->lgmin; lg <= gs->lgmax; lg++, theta += gs->lgincr) {
        long y = gs->ltmin * gs->ltincr;
        
        // get some handy locals
        float cos_theta = cos(URTOR(theta));
        float sin_theta = sin(URTOR(theta));
        
        // for each lattitude
        for (int lt = gs->ltmin; lt <= gs->ltmax; lt++, y += gs->ltincr) {
            long radius = GETR(gs, lt, lg);
            
            // quick exit path
            if (radius == long(VOIDGS(gs))) {
                num_void_nodes++;
                nodearray.addNode( Point3D(-1,-1,-1) );
                continue;
            }
            
            // must be good- compute and save the value
            Point3D pt;
            if (is_cartesian) { // cartesian coords
                pt = Point3D(theta,y,radius);   // actually has the x,y,z coords
            } else {    // cylindrical coords
                long x, z;                      
                x = long(sin_theta * radius);
                z = -long(cos_theta * radius);  // neg Z is toward user
                pt = Point3D(x,y,z);
            }
            num_good_nodes++;
            
            // convert to millimeters and save it in the nodearray
            const double microns_to_millimeters = 1000.0;
            int index = nodearray.addNode( pt/microns_to_millimeters );
            
            // get the node pointer so can set texture coords
            Node* new_node = nodearray.getNode(index);
            new_node->setTextureCoords(
                double(lt-gs->ltmin)/lt_size,
                double(lg-gs->lgmin)/lg_size);
        }
    }
    cyfree(gs);
    
    // more debugging
    cerr << "Had " << num_good_nodes << " good nodes and "
        << num_void_nodes << " void nodes\n";
    
    // now center the object (otherwise is who knows where)
    cerr << "Centering object\n";
    CenterObject();
    cerr << "done Centering object\n";
    
    // now do the tessellation
    cerr << "Generating mesh\n";
    int num_polys = 0, num_void_polys = 0;
    { // create a simple polygon list
        const int ypass = (gs->ltmax - gs->ltmin + 1);
        const int xpass = (gs->lgmax - gs->lgmin + 1);
        int num_verts = nodearray.getNumNodes();
        for (int x=0; x<(num_verts/xpass)-1; x++) 
            for (int y=0; y<ypass-1; y++) {
                int index = x * ypass + y;
                int index2 = (x+1) * ypass + y;
                
                // check vertices- if any are VOID, bail
                Node* n_a = nodearray.getNode(index);
                Node* n_b = nodearray.getNode(index+1);
                Node* n_c = nodearray.getNode(index2);
                Node* n_d = nodearray.getNode(index2+1);
                if (vertex_is_void(n_a->p) || vertex_is_void(n_b->p) ||
                    vertex_is_void(n_c->p) || vertex_is_void(n_d->p)) {
                    num_void_polys++;
                    continue;
                }
                
                // First half of square
                facearray.addFace(index, index2, index+1);
                num_polys++;
                
                // Other half of square
                facearray.addFace(index+1, index2, index2+1);
                num_polys++;
            }
    }
    cerr << "done Generating mesh\n";
    
    // more debugging
    cerr << "Had " << num_polys << " good polys and " 
        << num_void_polys << " void polygons\n";
    
    // set other parameters
    dynamics = deformable;
    set_default_material_properties();
    setDiffuseColor(1.0, 1.0, 1.0, 1.0);        // make the mesh white, with texture
    type = faces_only;
    
    // read in the texture
    //setTexture(NULL);
    
#ifdef TESTING_NEW_TEXTURE_COORDS
    // compute texture coords
    if (is_cartesian)
        nodearray.computeTextureCoords(Node::linear);
    else nodearray.computeTextureCoords(Node::cylindrical);
#endif
    
    // recompute normals (since surface is new)
    facearray.computeNormals();
    
    // average face normals to get per-vertex normals
    nodearray.interpolateNormals();
    
    // return value
    return(1);
}

int Object::Read(char* dirname_part, char* filename_part)
{
        // construct the full pathname w.r.t. the world.desc file
        char filename[256];
        sprintf(filename, "%s/%s", dirname_part, filename_part);

        // store path
    setPath(dirname_part);

        // read now
        return Read(filename);
}

int Object::Read(char* filename) {
    if (debug) cerr << "Object[" << name << "]::Read(" << filename << ")\n";
    
    // open the mesh data file
    FILE* meshfile;
    if ((meshfile = fopen(filename, "r")) == NULL) {
        printf("Could not open node data file\n");
        return(0);
    }
    
    // read `em in
    char* ext = myextension(filename);
    if (!strcmp(ext, ".mesh")) ReadMesh(meshfile);
    else if (!strcmp(ext, ".smf")) ReadSMF(meshfile);
    else if (!strcmp(ext, ".obj")) ReadOBJ(meshfile);
    else if (!strcmp(ext, ".cy")) ReadCyberware(meshfile);
    else if (!strcmp(ext, ".frd")) ReadFRD(meshfile);
    else if (!strcmp(ext, ".wrl")) ReadVRML(meshfile);
    else if (!strcmp(ext, ".iv")) ReadVRML(meshfile);
    else if (!strcmp(ext, ".inp")) ReadAmira(meshfile);
        else if (!strcmp(ext, ".sobj")) ReadSOBJ(meshfile);

    // close the file
    fclose(meshfile);
    
    // setup other vars
    dl_needs_refresh = 1;
    
    // and return success
    return(0);
}

void Object::ParseAttribute(char* buffer)
{
    // locals
    int had_bambient = 0;
    int had_bdiffuse = 0;
    int had_bspecular = 0;
    
    // parse the attribute
    char attr_name[80]; sscanf(buffer, "%s", attr_name);
    cerr << " " << attr_name << " [";
    
    // parse values based on attribute keyword
    if (!strcmp(attr_name, "ambient:")) {
        float r,g,b,a;
        sscanf(buffer, "%*s %f %f %f %f", &r, &g, &b, &a);
        setAmbientColor(r,g,b,a);
        if (!had_bambient) setBackAmbientColor(r,g,b,a);
        cerr << r << "," << g << "," << b << "," << a;
    }
    else if (!strcmp(attr_name, "bambient:")) {
        float r,g,b,a;
        sscanf(buffer, "%*s %f %f %f %f", &r, &g, &b, &a);
        setBackAmbientColor(r,g,b,a);
        cerr << r << "," << g << "," << b << "," << a;
        had_bambient = 1;
    }
    else if (!strcmp(attr_name, "diffuse:")) {
        float r,g,b,a;
        sscanf(buffer, "%*s %f %f %f %f", &r, &g, &b, &a);
        setDiffuseColor(r,g,b,a);
        if (!had_bdiffuse) setBackDiffuseColor(r,g,b,a);
        cerr << r << "," << g << "," << b << "," << a;
    }
    else if (!strcmp(attr_name, "bdiffuse:")) {
        float r,g,b,a;
        sscanf(buffer, "%*s %f %f %f %f", &r, &g, &b, &a);
        setBackDiffuseColor(r,g,b,a);
        cerr << r << "," << g << "," << b << "," << a;
        had_bdiffuse = 1;
    }
    else if (!strcmp(attr_name, "specular:")) {
        float r,g,b,a;
        sscanf(buffer, "%*s %f %f %f %f", &r, &g, &b, &a);
        setSpecularColor(r,g,b,a);
        if (!had_bspecular) setBackSpecularColor(r,g,b,a);
        cerr << r << "," << g << "," << b << "," << a;
    }
    else if (!strcmp(attr_name, "bspecular:")) {
        float r,g,b,a;
        sscanf(buffer, "%*s %f %f %f %f", &r, &g, &b, &a);
        setBackSpecularColor(r,g,b,a);
        cerr << r << "," << g << "," << b << "," << a;
        had_bspecular = 1;
    }
    else if (!strcmp(attr_name, "shininess:")) {
        float s;
        sscanf(buffer, "%*s %f", &s);
        setShininess(s);
        cerr << s;
    }
    else if (!strcmp(attr_name, "bshininess:")) {
        float s;
        sscanf(buffer, "%*s %f", &s);
        setBackShininess(s);
        cerr << s;
    }
    else if (!strcmp(attr_name, "texture:")) {
        char texturefilename[256];
        sscanf(buffer, "%*s %s", texturefilename);
        setTexture(texturefilename);
        cerr << texturefilename;
    }
    else if (!strcmp(attr_name, "texturecoords:")) {
        int model;
        sscanf(buffer, "%*s %d", &model);
        setTextureCoords(Node::texture_coord_type(model));
        cerr << model;
    }
    else if (!strcmp(attr_name, "behavior:")) {
        char str[80];
        sscanf(buffer, "%*s %s", str);
        behavior_type b = none;
        if (!strcmp(str, "none")) b = none;
        else if (!strcmp(str, "grabbing")) b = grabbing;
        else if (!strcmp(str, "cutting")) b = cutting;
        else if (!strcmp(str, "deleting")) b = deleting;
        else if (!strcmp(str, "poking")) b = poking;
        else if (!strcmp(str, "ablating")) b = ablating;
        else if (!strcmp(str, "cauterizing")) b = cauterizing;
        else if (!strcmp(str, "dilating")) b = dilating;
        else if (!strcmp(str, "coloring")) b = coloring;
        else if (!strcmp(str, "quisar_test")) b = quisar_test;
        else if (!strcmp(str, "beating")) b = beating;
        else if (!strcmp(str, "fluid")) b = fluid;
        else {  // old format-uses a number, so use old numbering
            int val = atoi(str);
            switch (val) {
            case int(none): b = none; break;
            case int(grabbing): b = grabbing; break;
            case int(cutting): b = cutting; break;
            case int(deleting): b = deleting; break;
            case int(poking): b = poking; break;
            case int(ablating): b = ablating; break;
            case int(cauterizing): b = cauterizing; break;
            case int(dilating): b = dilating; break;
            case int(coloring): b = coloring; break;
            case int(quisar_test): b = quisar_test; break;
            case int(beating): b = beating; break;
            case int(fluid): b = fluid; break;
            }
        }
        setBehavior(b);
        cerr << b;
    }
    else if (!strcmp(attr_name, "dynamics:")) {
        char str[80];
        sscanf(buffer, "%*s %s", str);
        dynamics_type d = rest;
        if (!strcmp(str, "rest")) d = rest;
        else if (!strcmp(str, "rigid")) d = rigid;
        else if (!strcmp(str, "deformable")) d = deformable;
        else if (!strcmp(str, "articulate")) d = articulate;
        else if (!strcmp(str, "button")) d = button;
        else {  // old format-uses a number, so use old numbering
            int val = atoi(str);
            switch (val) {
                                case 0: d = rest; break;
                                case 1: d = rigid; break;
                                case 2: d = deformable; break;
                                case 3: d = articulate; break;
                                case 4: d = button; break;
            }
        }
        setDynamics(d);
        cerr << d;
    }
    else if (!strcmp(attr_name, "mass:")) {
        float m;
        sscanf(buffer, "%*s %f", &m);
        getNodeArray()->setMasses(m);
        cerr << m;
        }
    else if (!strcmp(attr_name, "deformonlyongrab:")) {
        int i;
        sscanf(buffer, "%*s %d", &i);
                do_deform_only_on_grab = i;
        cerr << i;
    }
    else if (!strcmp(attr_name, "springconstant:")) {
        float sc;
        sscanf(buffer, "%*s %f", &sc);
        getEdgeArray()->setSpringConstants(sc);
        cerr << sc;
    }
    else if (!strcmp(attr_name, "dampconstant:")) {
        float dc;
        sscanf(buffer, "%*s %f", &dc);
        getNodeArray()->setVelocityDampConstants(dc);
        cerr << dc;
    }
    else if (!strcmp(attr_name, "nummethod:")) {
        int nm;
        sscanf(buffer, "%*s %d", &nm);
        setNumMethod(Object::numerical_method(nm));
        cerr << nm;
    }
    else if (!strcmp(attr_name, "visible:")) {
        int on;
        sscanf(buffer, "%*s %d", &on);
        SetVisible(on);
        cerr << on;
    }
        else if (!strcmp(attr_name, "translate:")) {
        float x, y, z;
        sscanf(buffer, "%*s %f %f %f", &x, &y, &z);
                Point3D vector(x,y,z);
        Move(vector);
        cerr << x << " " << y << " " << z;
    }
        else if (!strcmp(attr_name, "scale:")) {
        float scale_factor;
        sscanf(buffer, "%*s %f", &scale_factor);
        Scale(scale_factor);
        cerr << scale_factor;
    }
        else if (!strcmp(attr_name, "rotate:")) {
        char axis; double angle;
        sscanf(buffer, "%*s %c %lf", &axis, &angle);
        Rotate(axis, angle);
        cerr << axis << " " << angle;
    }
    else if (!strcmp(attr_name, "collisionextent:")) {
        int ce;
        sscanf(buffer, "%*s %d", &ce);
        setCollisionExtent(Object::collision_extent_type(ce));
        cerr << ce;
    }
    else if (!strcmp(attr_name, "linkedpart:")) {
        int lp;
        sscanf(buffer, "%*s %d", &lp);
        setLinkedPart(lp);
        cerr << lp;
    }
    else if (!strcmp(attr_name, "hingenode:")) {
                        int i; 
                        int si;
                        float x,y,z;
                        double maxopen, maxclosed;
                        int col0 = -1;
                        int col1 = -1;
                        int col2 = -1;

                        sscanf(buffer, "%*s %d %d %f %f %f %lf %lf %d %d %d", 
                                &si, &i, &x, &y, &z, &maxopen, &maxclosed, 
                                &col0, &col1, &col2);
                        Point3D r(x,y,z);
                        this->addToHinges(si,i,r,maxopen,maxclosed,col0,col1,col2);

                        cerr << si << " " << i << " " << r << " " << maxopen << " " 
                                 << maxclosed << " " << col0 << " " << col1 << " " << col2;
        }
    else if (!strcmp(attr_name, "sensoroffset:")) {
        float x,y,z;
        sscanf(buffer, "%*s %f %f %f", &x, &y, &z);
        Point3D p(x,y,z);
        setSensorOffset(p);
        cerr << p;
    }
    else if (!strcmp(attr_name, "tipnode:")) {
        int tn;
        sscanf(buffer, "%*s %d", &tn);
        setTip(tn);
        cerr << tn;
    }
    else if (!strcmp(attr_name, "collisionfaces:")) {

        // make some locals

        char* delimiters = " \t";

        char* ptr = strtok(buffer, delimiters); // flush attr_name

        

        // parse the face indices on the line and add them to the collisionfaces list

        while ((ptr = strtok(NULL, delimiters)) != NULL ) {

            int face_index = atoi(ptr);

            Face* face = getFaceArray()->getFace(face_index);

            if (!face) {

                cerr << "Bad collision face index: " << face_index << endl;

                        }

            else {

                // add to collision face list

                collision_faces.Append(face);

            }

        }  

    }   
        else if (!strcmp(attr_name, "cuttingfaces:")) {
                char* delimiters = " \t\n";
                char* ptr = strtok(buffer, delimiters);
                while(ptr) {
                        ptr = strtok(NULL,delimiters);
                        if (ptr == NULL) break;
                        int f = atoi(ptr);
                        cerr << f << " " ;
                        addToCuttingFaces(getFaceArray()->getFace(f));
                }
        }
        else if (!strcmp(attr_name, "cuttingedges:")) {
                char* delimiters  = " \t\n";
                char* ptr = strtok(buffer, delimiters);
                while(ptr) {
                        ptr = strtok(NULL,delimiters);
                        if (ptr == NULL) break;
                        int e = atoi(ptr);
                        cerr << e << " " ;
                        addToCuttingEdges(getEdgeArray()->getEdge(e));
                }
        }
    else {    // unknown
        cerr << "UNKNOWN ATTRIBUTE";
    }
    
    // end the line
    cerr << "]\n";
}

void Object::addToHinges(int si, int ai, Point3D r, double maxopen, double maxclosed, 
                                                                                                 int col0, int col1, int col2)
{ 
        hinge h;
  h.sensor_index = si;
        h.node_index = ai; 
        h.rotation_vector = r;
        h.max_open = maxopen;
        h.max_closed = maxclosed;

        int posnum = 0;
        int negnum = 0;
        
        for (int i=0; i<3; i++) {
                h.colors_pos[i] = -1;
                h.colors_neg[i] = -1;
        }


        if (col0 > -1) {
                h.colors_pos[posnum] = col0;
                posnum++;
        }
        else {
                h.colors_neg[negnum] = col0;
                negnum++;
        }

        if (col1 > -1) {
                h.colors_pos[posnum] = col1;
                posnum++;
        }
        else {
                h.colors_neg[negnum] = col1;
                negnum++;
        }

        if (col2 > -1) {
                h.colors_pos[posnum] = col2;
                posnum++;
        }
        else {
                h.colors_neg[negnum] = col2;
                negnum++;
        }

        hinges.Append(h);
        return;
}
  
void Object::setAmbientColor(float red, float green, float blue, float alpha)
{ 
    // set the values
    mat_ambient[0] = red; 
    mat_ambient[1] = green; 
    mat_ambient[2] = blue;
    mat_ambient[3] = alpha;
    
    // update display list
    dl_needs_refresh = 1; 
}

void Object::setDiffuseColor(float red, float green, float blue, float alpha)
{ 
#ifdef SUN
    // Hack: There is a bug in the Solaris version of OpenGL where 
    // alpha-blended textures aren't handled correctly and produce an
    // artifact.  The workaround is to have the r,g,b values the same
    // as the alpha if you're doing this
    if (texture) {
        red = green = blue = alpha;
    }
#endif
    
    // set the values
    mat_diffuse[0] = red; 
    mat_diffuse[1] = green; 
    mat_diffuse[2] = blue;
    mat_diffuse[3] = alpha;
    
    // update display list
    dl_needs_refresh = 1; 
}

void Object::setSpecularColor(float red, float green, float blue, float alpha)
{ 
    // set the values
    mat_specular[0] = red; 
    mat_specular[1] = green; 
    mat_specular[2] = blue;
    mat_specular[3] = alpha;
    
    // update display list
    dl_needs_refresh = 1; 
}

void Object::setShininess(float s)
{ 
    // set the values
    mat_shininess[0] = s; 
    
    // update display list
    dl_needs_refresh = 1; 
}

void Object::setBackAmbientColor(float red, float green, float blue, float alpha)
{ 
    // set the values
    back_mat_ambient[0] = red; 
    back_mat_ambient[1] = green; 
    back_mat_ambient[2] = blue;
    back_mat_ambient[3] = alpha;
    
    // update display list
    dl_needs_refresh = 1; 
}

void Object::setBackDiffuseColor(float red, float green, float blue, float alpha)
{ 
    // Hack: There is a bug in the Solaris version of OpenGL where 
    // alpha-blended textures aren't handled correctly and produce an
    // artifact.  The workaround is to have the r,g,b values the same
    // as the alpha if you're doing this
    if (texture) {
        red = green = blue = alpha;
    }
    
    // set the values
    back_mat_diffuse[0] = red; 
    back_mat_diffuse[1] = green; 
    back_mat_diffuse[2] = blue;
    back_mat_diffuse[3] = alpha;
    
    // update display list
    dl_needs_refresh = 1; 
}

void Object::setBackSpecularColor(float red, float green, float blue, float alpha)
{ 
    // set the values
    back_mat_specular[0] = red; 
    back_mat_specular[1] = green; 
    back_mat_specular[2] = blue;
    back_mat_specular[3] = alpha;
    
    // update display list
    dl_needs_refresh = 1; 
}

void Object::setBackShininess(float s)
{ 
    // set the values
    back_mat_shininess[0] = s; 
    
    // update display list
    dl_needs_refresh = 1; 
}


// Note: caller should also do a nodearray.computeTextureCoords()
void Object::setTexture(char* filename)
{
    char texture_filename[256];
    if (filename)     // use the name we were given
        strncpy(texture_filename, filename, 256);
    else {    // create name from name of object
                char regular_filename[256];
                char* filePath = getPath();
                
                if (strcmp(getPath(), "")) {
                        sprintf(regular_filename, "%s/%s", getPath(), getName());
                }
                else {
                        sprintf(regular_filename, "%s", getName());
                }
        strncpy(texture_filename, regular_filename, 256);        // copy over name
        char* ptr = myrindex(texture_filename, '.');      // strip off extension
        if (ptr) {
                        sprintf(ptr, ".rgb");                    // add on a .rgb
                }
        else {
                        sprintf(texture_filename, "%s.rgb", texture_filename);
                }
    }

    // create the texture
    // NOTE: width and height must be powers of 2 for texture! 
    // (maybe next release of OGL, but was still case in 1.2!)
    if (debug) 
        cerr << "Reading texture from file [" << texture_filename << "]\n";
    
    // read the texture
    texture = read_texture(texture_filename, 
        &texture_size.x, &texture_size.y, &texture_size.z);
    
    // debugging
    if (debug) {
        if (texture) 
            cerr << " size = " << texture_size.x << "x" << texture_size.y 
            << " @ " << texture_size.z << endl;
        else cerr << " texture [" << texture_filename << "] not found!\n";
    }
    
    // setup other vars
    dl_needs_refresh = 1;
}

void Object::setTextureCoords(Node::texture_coord_type t)
{
    nodearray.computeTextureCoords(t);
    dl_needs_refresh = 1;
}

void Object::moveTexture(Point3D translation, Point3D scaling)
{
    nodearray.moveTexture(translation, scaling);
    dl_needs_refresh = 1;
}

void Object::tipUpdate()
{
        if (nodearray.getNumNodes() <= 0) return;

    tipdelta = getTip() - oldtip;
    oldtip = getTip();
}

void Object::Update()
{
    if (debug) cerr << "Object[" << name << "]::Update()\n";

        // tell other threads that we're modifying this one now
        // (poor man's mutex- just want to indicate, not lock)
        being_modified = 1;




    // handle sensors/motion links here
    if (sensor_p) {    
        
        // Note: I'd ifdefed this out again to *quickly* provide support for 
        // multihinged objects for Cynthia's articulating hand- the right thing 
        // to do is to put it into the more general xformObject framework.  But 
        // either way, either method needs to provide support for xforms of 
        // "rest'(static) objects, hinged objects, multihinged objects, 
        // "plunger" objects (syringes), etc
#define NEW_XFORM_OBJECT_WAY
#ifdef NEW_XFORM_OBJECT_WAY
                // xform hinged object (hinged objects are still a special case 
                // until we generalize it into the xformObject() method
        if (linked_part == -1)
           xformHingedObject();
        else // update us (handles parts, everything, etc)
           xformObject();
#else
        // if we're linked to a sensor, update us
        // this is kinda like an external force acting on the
        // whole object (maybe it should be a force, now it's
        // just a 'move here' kinda thing)
        if (linked_part == -1) {                        // xform hinged object
           xformHingedObject();
        }
        else if (linked_part == 0) {            // xform whole thing
            nodearray.setToInitialPosition();   // xform from initial position
            nodearray.xformNodesOfPart(-1, sensor_p->m, sensor_offset, 
                sensor_p->pos);
        }
        else {  // just a part of object moved- do translation only
            // get closest node in our mesh
            //Node* closest_node = findClosestNode(sensor_p->pos);
            //Point3D motion_vector = sensor_p->pos - closest_node->p;
            Point3D motion_vector = sensor_p->pos - sensor_p->init_pos;
            nodearray.moveNodesOfPart(linked_part, motion_vector, 1);
        }
#endif
        

                // indication that nodes were modified
        addEverythingToUpdateList();


        // no matter what, we've moved
        dl_needs_refresh = 1;
    }
    
    // Note: we don't reset the nodearray's collision flags here because we only
    //        do collision testing on the first update cycle now
    switch (dynamics) {
                case rest: {    
                        // for objects that have no dynamics of their own
                        // Note: they CANNOT be picked up and DON'T deform (deformable) and 
                        //        DON'T undergo rigid body dynamics
                        break;
                }

                case button: {
                        // for objects that have no dynamics of their own
                        // Note: they CAN be picked up and DON'T deform (deformable) and 
                        //        DON'T undergo rigid body dynamics
                        break;
                }

                case rigid: {
                        // rigid bodies don't spontaneously move (sez Newton)
                        // but maybe they should handle a per-node velocity, etc
        
                        // if we're not attached to a sensor, then process rigid body dynamics
                        if (!sensor_p) {
            
                                // Note: There is a problem if we do the current code, then grab the 
                                // object with a sensor.  In the current updateNodesRigid() code, it has
                                // to reset the init_pos of the nodes as it goes.  The problem is
                                // that then, later, if we grab the tool, our sense of the initial
                                // position is wrong, as is the sensor_offset, hinge_point, etc


                                // update node positions- rigid method
                                int nodes_moved = nodearray.updateNodesRigid();


                                // if have updated, assume we need to refresh graphics
                                if (nodes_moved) {
                                        // indication that nodes were modified

                                        addEverythingToUpdateList();


                                        dl_needs_refresh = 1;
                                }


                                // this is for our pole and ring edge objects, so that the quads
                                // redraw as front-facing as the camera moves
                                if (type == edges_only) {
                                        dl_needs_refresh = 1;

                                }

                        }
                        break;
                }

                case deformable: {
                        // useful mode for some applications
                        if (do_deform_only_on_grab && !being_grabbed) {
                                if (nodearray.orderedNodesFlag)  {
                                        edgearray.computeNewRestLengths();
                                        nodearray.orderedNodesFlag = 0;
                                }
                                break;
                        }

                        // update node positions
                        int nodes_moved = nodearray.updateNodes();
        
                        // if have updated and moved something, we need to update graphics
                        if (nodes_moved) {
                                // indication that nodes were modified

                                addEverythingToUpdateList();


                                if (do_preserve_volume) { //attempt to maintain volume
                                        PreserveVolume();
                                }

                                dl_needs_refresh = 1;
                        }               
                        // cbnote and add moved faces to bs tree

        
                        break;
                }
        
                case articulate:  {
                        // update myself
                        extern game_enum game_mode;
                        if (game_mode == knot_game || game_mode == vessel_knot_game) {
                                nodearray.updateThread();

                        }
                        else {

                                // for vessel_game or others
                                nodearray.updateArticulateEdge();
                        }


                        // indication that nodes were modified

                        addEverythingToUpdateList();



                        // if have updated, assume we need to refresh graphics
                        dl_needs_refresh = 1;
        
                        break;
                }
    }  // end switch




        // no longer in use
        being_modified = 0;
}

// print out position and angles of the given plane tool
void Object::getPosAngle(Point3D* pos_p, Point3D* rot_p)
{
    if (strncmp(name, "plane #", 7))
        return;
    // node 0 is the center, normal of node 0 is the direction
    Node *node0 = nodearray.getNode(0);
    Point3D pos = node0->p;
    Point3D normal = node0->getNormal();
    double dotx = normal.Dot(Point3D(1,0,0));
    double doty = normal.Dot(Point3D(0,1,0));
    double dotz = normal.Dot(Point3D(0,0,1));
    double thetax = (acos(dotx))*(180.0/PI);
    double thetay = (acos(doty))*(180.0/PI);
    double thetaz = (acos(dotz))*(180.0/PI);
    
    *pos_p = pos;
    rot_p->x = thetax;
    rot_p->y = thetay;
    rot_p->z = thetaz;
}

// split an object in two with respect to the given planar tool
void Object::SplitColor(Object* tool)
{
    if (debug) cerr << "Object::SplitColor(" << tool << ")\n";
    
    if (!strstr(tool->name, "grabby") && !strstr(tool->name, "plane"))
        return;    // only split about the planar tools
    
    // the 0th face of the tool is the correct plane
    Face *splitting_face = tool->facearray.getFace(0);
    nodearray.SplitColor(splitting_face);
    
    // update graphics required in case viewing
    dl_needs_refresh = 1;
}

void Object::DeleteNodesOfColor(int color)
{
    if (debug) cerr << "Object::DeleteNodeOfColor(" << color << ")\n";
    facearray.EraseFacesOfColor(color); // get rid of the faces...
    edgearray.DeleteFacelessEdges();
    //nodearray.DeleteNodesOfColor(color);      // and the nodes that love them
    dl_needs_refresh = 1;
}

double Object::distanceToBoundingBox(Point3D pt)
{
    // get bounding box of object
    Point3D cur_min, cur_max;
    getBoundingBox(&cur_min, &cur_max);
    
    // if we're inside, forget it
    if ((pt.x > cur_min.x) && (pt.y > cur_min.y) && (pt.z > cur_min.z) &&
        (pt.x < cur_max.x) && (pt.y < cur_max.y) && (pt.z < cur_max.z))
        return(0);      // are inside the bbox
    
    // otherwise, calc distance to it
    // this is totally wrongo:
    //Point3D dist_vect(MIN(cur_min.x-pt.x,pt.x-cur_max.x),
    //    MIN(cur_min.y-pt.y,pt.y-cur_max.y),
    //    MIN(cur_min.z-pt.z,pt.z-cur_max.z));
    
    // let's try again:
    double xdist, ydist, zdist;
    xdist = ydist = zdist = 0.0;
    if (pt.x > cur_max.x)
      xdist = pt.x - cur_max.x;
    else if (pt.x < cur_min.x)
      xdist = cur_min.x - pt.x;
    if (pt.y > cur_max.y)
      ydist = pt.y - cur_max.y;
    else if (pt.y < cur_min.y)
      ydist = cur_min.y - pt.y;
    if (pt.z > cur_max.z)
      zdist = pt.z - cur_max.z;
    else if (pt.z < cur_min.z)
      zdist = cur_min.z - pt.z;

    // and return the distance
    return(sqrt(xdist*xdist + ydist*ydist + zdist*zdist));
}

// this object is a tool (ie has a behavior)
void Object::HandleBehavior()
{
  if (debug) cerr << "Object[" << getName() << "]::HandleBehavior()\n";

  // update the notion of the tip delta movement
  tipUpdate();
    
  // handle both interactive and spontaneous behaviors
  HandleInteractiveBehavior();
  HandleSpontaneousBehavior();
}

// Handles interaction type behavior, i.e. a tool interacting with the object. (Chris, HI)
void Object::HandleInteractiveBehavior() 
{
  // Check if sensor is available, otherwise there is no interactive behavior to process 
  if (!sensor_p) {    
    closest_node = NULL;    // reset things
    return;
  }   

  // Local constants
  // Note: setting this value this high certainly seems odd, but if it is
  //             set lower, then grabbing is somewhat problematic.  It seems
  //             that we all too eaily loose our grasp on our closest_node and
  //             break our grab and do it in such a way that it's still fixed 
  //             in position.  So, for now, we put it big- later, let's solve this.
  const double max_distance_to_affect = 250.5;    // mm
    
  // Vectors for force calculation
  Point3D haptic_force(0,0,0);
  float twist_force = 0, grip_force = 0;
        
  //
  // INTERACTION BEHAVIORS
  //
  // The rationale is this- a tool affects an object that it's in contact
  // (i.e., colliding) with.  Therefore, this is a special case of a
  // collision (for example, a drill collides with a skull and, in so
  // doing, makes a hole).
  // 
  // However, once in contact (meaning once the closest_node is set),
  // we remain in contact.  This takes care of a few things: sometimes
  // when in contact, the tool can barely move out of the bbox of the
  // affected object and it would release.  Also, this way is faster
  // because we don't do a search for what object we're in contact with
  // if we've already got one.
  //
  // Also note- the way this is structured right now, on first contact,
  // we can only affect one object at a time (the first one we find that 
  // intersects our bounding sphere/box).
  //
  switch (behavior) {          
  case quisar_test: {
    // if not the v2 sensor, forget it
    if (sensor_p->type != Sensor::haptic_v2_sensor) break;
      
    // cast away!
    Haptic_v2* haptic_p = (Haptic_v2*)sensor_p;
      
    // This thing is supposed to handle the first contact: that's all.
    Object* other_obj = NULL;
      
    // we drop a grabbed node if sensor becomes inactive
    if (closest_node && (haptic_p->active[0] < 0.5)) {
      closest_node->setFixedType(Node::NOT);  
      being_grabbed = 0;
      closest_node->getNodeArray()->getObject()->being_grabbed = 0;
      closest_node = NULL;
    }
      
    // now we find a closest node if we don't have one from prior grab
    if (!closest_node)
      closest_node = objectarray_p->FindClosestUnlinkedNode(getTip());
      
    // and find the object it belongs to
    if (closest_node)
      other_obj = closest_node->getNodeArray()->getObject();
      
    // we may not have found any closest node if we're outside bboxes
    if (!closest_node) {
      haptic_p->my_space = haptic_p->my_space_v;
      haptic_p->points_to_save->erase_all();
      break;
    }
      
    Point3D vect = getTip() - closest_node->p;
      
    if ((haptic_p->active[0] < 0.5) || 
        (vect.Length() > max_distance_to_affect)) {
      // We handle only the triangles right now.
      haptic_p->my_space = haptic_p->my_space_t;
      haptic_p->points_to_save->erase_all();
      haptic_p->my_space_t->erase_all();
      triangle *t;
      Face *f;
      for(int i=0; i< closest_node->numFaceLinks(); i++) {
        f  = closest_node->getFaceLink(i);
        t = new triangle(f->getNodeLink(0)->p, f->getNodeLink(1)->p, 
                         f->getNodeLink(2)->p, f->getNodeLink(0)->getNormal(),
                         f->getNodeLink(1)->getNormal(), 
                         f->getNodeLink(2)->getNormal(), (char) 127);
        haptic_p->my_space_t->add_new_triangle(t);
      }
      break;
    }
      
    // now we have a closest node, and he is close enough => we
    // need to send the matrix of points
    if (!being_grabbed) { // this is our first grab, maybe do special

                if (other_obj->getNumMethod() == quasi_ordered ||

                        other_obj->getNumMethod() == euler_ordered) {      
                        // order nodes
                        other_obj->getNodeArray()->OrderNodes(closest_node);
                }
                being_grabbed = 1;
                other_obj->being_grabbed = 1;
    }
      
    // set the position of the node to our tip
    // if (use_fext) closest_node->setFext(50.0*vect); // using forces
    closest_node->p = getTip();
      
    // make it stop moving, since we've got it
    closest_node->v = Point3D(0,0,0);
      
    // we need to fix a node, otherwise it can snap back if we have
    // multiple Update calls between calls to HandleBehavior
    closest_node->setFixedType(Node::POS); 
            
    // Now we must find which points we don't have the force.
    // Finding the force for those, putting what we need in the
    // space_i we have, and it will be sent later.
    Point3D pttmp;
    point *mypoint;
    double xmin = ceil((closest_node->p.x - RAYON) / MATRIX) * MATRIX;
    double ymin = ceil((closest_node->p.y - RAYON) / MATRIX) * MATRIX;
    double zmin = ceil((closest_node->p.z - RAYON) / MATRIX) * MATRIX;
            
    haptic_p->my_space_i->erase_all();
            
    for(double x=xmin; x <= (closest_node->p.x + RAYON); x+=MATRIX) {
      for(double y=ymin; y <= (closest_node->p.y + RAYON); y+=MATRIX) {
        for(double z=zmin; z <= (closest_node->p.z + RAYON); z+=MATRIX) {
          //             cerr << x << ":" << y << ":" << z << endl;
          pttmp.x = x; pttmp.y = y; pttmp.z = z;
          if((closest_node->p - pttmp).Squared() <= RAYON * RAYON) {
            mypoint = haptic_p->points_to_save->isthere(&pttmp);
            if(mypoint) {
              point *pointcopy;
              pointcopy = new point(mypoint);
              haptic_p->my_space_i->add_new_point(pointcopy);
            }
            else {
              point *pointcopy;
              mypoint = new point(&pttmp, what_force_if(closest_node, &pttmp));
              pointcopy = new point(mypoint);
              haptic_p->points_to_save->add_new_point(pointcopy);
              haptic_p->my_space_i->add_new_point(mypoint);
            }
          }
        } // for z
      } // for y
    } // for x
    haptic_p->my_space = haptic_p->my_space_i;
  }  // case quisar test
    
  case grabbing: {
    // need to find closest node, its object, and grab it:
    Object* other_obj = NULL;
    float sensor_activation = sensor_p->active[0];      // activation status of sensor
    float prev_sensor_activation = prev_grabbing_activation;
    prev_grabbing_activation = sensor_activation;

    // check if we were already grabbing before, but didn't catch anything
    // -> if so, we won't try grabbing again, because the tool has to be opened 
    //    again first, to allow grabbing again
    if ((prev_sensor_activation >= 0.5) && (being_grabbed == 0)) {
                                // let's just stop right here!
      break;
    }

    // we drop a grabbed node if sensor becomes inactive
    // this also sets the closes_node to NULL even if we didn't grab previously!  
    if (closest_node && (sensor_activation < 0.5)) {
      closest_node->setFixedType(Node::NOT);  
            
#ifdef JOEL
      int hack_index = closest_node->getIndex();
                                //------------------------- HACK here ------------------------------
      Object *other_vessel = objectarray_p->findObject("other vessel");
      if (other_vessel != NULL) {
        Object *left_vessel = objectarray_p->findObject("left vessel");
        Object *right_vessel = objectarray_p->findObject("right vessel");
                
        Node* n = left_vessel->nodearray.getNode(hack_index);
        n->setFixedType(Node::NOT);
        n = right_vessel->nodearray.getNode(hack_index);
        n->setFixedType(Node::NOT);
        n = other_vessel->nodearray.getNode(hack_index);
        n->setFixedType(Node::NOT);
      }
                                // -----------------------------------------------------------------
#endif
                                // Only resets grabbing, if we were really grabbing before
                                // Chris: This fixes a bug if a second laparascope comes close to the
                                //      object. The moving in of the second laparoscope made the object
                                //      automatically get detached/ungrabbed from the first laparoscope!
      if (being_grabbed == 1) {
        being_grabbed = 0;
        closest_node->getNodeArray()->getObject()->being_grabbed = 0;
      }

      closest_node = NULL;
    }    // if sensor inactive
        
    // now we find a closest node if we don't have one from prior grab
    if (!closest_node)
      closest_node = objectarray_p->FindClosestUnlinkedNode(getTip());
        
    // and find the object it belongs to
    if (closest_node)
      other_obj = closest_node->getNodeArray()->getObject();
        
    // we may not have found any closest node if we're outside bboxes
    if (!closest_node) break;
        
    // don't grab tied nodes
    if (closest_node->getTieNode())  {
      closest_node = NULL;
      break;
    }
        
    // we may be inactive
    if (sensor_activation < 0.5) break; // keep the closest node for fun
        
    // if it's too far away, forget it
    Point3D vect = getTip() - closest_node->p;
    if (vect.Length() > max_distance_to_affect) {
      closest_node = NULL;
      if (debug) cerr << "Too far away to grab\n";
      break;
    }
        
    // now we actually grab something
    if (!being_grabbed) { // this is our first grab, maybe do special

                if (other_obj->getNumMethod() == quasi_ordered ||

                        other_obj->getNumMethod() == euler_ordered) {
                        other_obj->getNodeArray()->OrderNodes(closest_node);
                }
                being_grabbed = 1;
                other_obj->being_grabbed = 1;
    }
        
    // check if button
    if (closest_node->getNodeArray()->getObject()->getDynamics() == button) {
                                // break, because there are no forces rendered for button
      break;
    }

    // set the old position of the grabbed node -- this allows us, for
    // example, to break a grab into steps when pulling something too
    // far, as in updateThread
    closest_node->p_old = closest_node->p;
    // set the position of the node to our tip
    // if (use_fext) closest_node->setFext(50.0*vect); // using forces
    closest_node->p = getTip();
        
    // make it stop moving, since we've got it
    closest_node->v = Point3D(0,0,0);
        
    // we need to fix a node, otherwise it can snap back if we have
    // multiple Update calls between calls to HandleBehavior
    closest_node->setFixedType(Node::POS); 
        
#ifdef JOEL
    // ----------------------------and HACK here ---------------------
    Object *other_vessel = objectarray_p->findObject("other vessel");
    if (other_vessel != NULL) {
      Object *left_vessel = objectarray_p->findObject("left vessel");
      Object *right_vessel = objectarray_p->findObject("right vessel");
            
      Point3D disp =closest_node->p - closest_node->init_pos;
      int hack_index = closest_node->getIndex();
      Node* n = left_vessel->nodearray.getNode(hack_index);
      n->p = n->init_pos + disp;
      n->v = Point3D(0.0,0.0,0.0);
      n->setFixedType(Node::POS);
      n = right_vessel->nodearray.getNode(hack_index);
      n->p = n->init_pos + disp;
      n->v = Point3D(0.0,0.0,0.0);
      n->setFixedType(Node::POS);
      n = other_vessel->nodearray.getNode(hack_index);
      n->p = n->init_pos + disp;
      n->v = Point3D(0.0,0.0,0.0);
      n->setFixedType(Node::POS);
    }
    // ----------------------------------------------------------------
#endif
        

        // calculate haptic force
    if (closest_node->getNodeArray()->getObject()->getRigidHinge() == NULL) {
      // regular haptics force
      haptic_force = (closest_node->getF());
    }
    else {
      // Chris: haptics force if we pull on a hinged object.
      // The farther away you go, the stronger you are pulled back
      // to the hinged object
      Point3D tip_pt = getTip();
      float distance = closest_node->p_old.Dist(tip_pt);
      haptic_force = closest_node->p_old - tip_pt;
      haptic_force *= 10 * distance * distance;
    }
                        
    // if we're grabbing something, give us some force = hack
    grip_force = (sensor_activation - 0.5)*2 * 100;     
    twist_force = 0;    // but no twist for now
                        
    if (debug) 
      cerr << "haptic_force: " << haptic_force 
           << ",twist=" << twist_force << ", grip=" << grip_force 
           << endl;
    break;
  } // case grabbing
                
  case deleting: {
    Object* other_obj = NULL;
      
#ifdef DELETE_CURRENT_OBJECT_ONLY
    other_obj = current_object;
#else
    // find the closest node 
    closest_node = objectarray_p->FindClosestUnlinkedNode(getTip());
      
    // and find the object it belongs to
    if (closest_node) {
      other_obj = closest_node->getNodeArray()->getObject();
    }
#endif
      
    // call Jeremie's cutting code here
    if (sensor_p->active[0] && other_obj) {
      other_obj->CuttingByDeleting(this);
    }
    break;
  }  // case deleting
      
  case coloring: {
    if (sensor_p->active[0] < 0.5) break;
      
    // get the object
    closest_node = objectarray_p->FindClosestUnlinkedNode(getTip());
    if (!closest_node) break;
      
    // find the object it belongs to
    Object* other_obj = closest_node->getNodeArray()->getObject();
    cerr << "Coloring on object [" << other_obj->getName() << "]\n";
    int color = 10;
      
    // get our region
    Point3D min, max;
    getBoundingBox(&min, &max);
      
    // and color it in our region
    other_obj->ColorNodesWithinBoundingBox(min, max, color);
      
    break;
  } // case coloring 
      
  case poking: {
    // Note: This is handled by the ResolveCollisions() framework now
    break;
  } // case poking
      
  case ablating:
  case cauterizing:

    //#define SOUND
#if defined(SOUND) && defined(SUN)
    // if the tool is on, make a sound
    if (sensor_p->num_active > 1) {
      if (sensor_p->active[1] > 0.5) {
        system("/usr/bin/audioplay fastbusy.au &");
      }
    }
#endif
    break;

  case none: 
    // no behavior defined, do nothing
    break;
                
  default: 
    // not handled here
    break;
                        
  } // switch behavior
    
  // render the haptic force 
  // (used to have sensor must be active=kill switch)
  if (sensor_p->type == Sensor::haptic_v1_sensor) {
    Haptic_v1* haptic_p = (Haptic_v1*)sensor_p;
    haptic_p->force = haptic_force;
    if (haptic_p->GetGridMode() == false) {
      haptic_p->twist_force = twist_force; 
      haptic_p->grip_force = grip_force; 
    }
    else { 
      // the haptic controller uses a forcegrid...
      if (closest_node && (sensor_p->active[0] > 0.5)) {
        haptic_p->SetActive(true);
      }
      else {
        haptic_p->SetActive(false);
      }
    }     
    //if (haptic_force.z != 0)
    //cout << "\nForce sent to haptic_v1: " << haptic_force.x << " " 
    //<< haptic_force.y << " " << haptic_force.z << endl;
    if (debug) {
      cerr << "[" << getName() << "]- haptic force = " << haptic_force 
           << ", " << haptic_force.Length() << endl;
    }
  }
}

// Handles spontaneous behavior such as beating or fluid type dynamics which happen
// without outside influence. (Chris, HI)
void Object::HandleSpontaneousBehavior() 
{  
  //
  // SPONTANEOUS BEHAVIORS
  //
  // Spontaneous behaviors move on their own without any outside influence.
  //
  switch (behavior) {          
  case beating: {
    const double repeat_interval = 1.0;  // seconds
    const double magnitude = 0.2;
    extern Timer global_timer;    // global system timer
    double t = global_timer.GetElapsedTime();
    t /= repeat_interval;   // interval- every n seconds
    t = t - int(t);   // get just the fractional part of a second
    t *= magnitude;

    // scale it in place
    nodearray.setToInitialPosition();
    nodearray.ScaleByNormal(t);
    //edgearray.scaleRestLengths(value);
    dl_needs_refresh = 1;
    addEverythingToUpdateList();
    break;
  } // case beating

  case fluid: {
    // use perlin noise function to create fluid effect in y (height) direction (Chris: HI)
    extern int iteration;
    // factor to reduce and save cpu processing power for fluids 
    const int iterations_per_update = 5;     
    // bigger value = slower liquid movement = more viscosity 
    // -> The spring constant donates the viscosity level!
    float viscosity = getEdgeArray()->getEdge(0)->getSpringConstant();
    // get perlin noise function
    PerlinNoise* noise = PerlinNoise::Instance();
    if ((iteration % iterations_per_update) == 0) {
      NodeArray* fluid_nodes = getNodeArray();
      int numnodes = fluid_nodes->getNumNodes();
      for (int i = 0; i < numnodes; i++) {
        Node* node = fluid_nodes->getNode(i);
        double y = noise->pnoise(node->p.x, node->p.z, 
                                 ((float)iteration / iterations_per_update / viscosity));
        node->p.y = getHeight() + y / 10;
      }
                                // to be called, otherwise nothing gets updated and shown on the screen
      addEverythingToUpdateList();
      dl_needs_refresh = 1;
    }
    break;
  } // case fluid

  default: 
    // not handled here
    break;
                                
  } // switch behavior
}


// top level draw function
void Object::Draw()
{
  if (visible) {
    DrawGeometry();
  }
  if (selected)
    DrawBoundingBox();
}

// this draws the edges and/or faces and/or nodes of the actual object
void Object::DrawGeometry()
{ 
  if (debug) cerr << "Object[" << name << "]::Draw(" << &display_mode << ")\n";
    
  // find out if we're selecting/picking
  GLint rendermode;
  glGetIntegerv(GL_RENDER_MODE, &rendermode);
    
  // if have updated, assume we need to refresh graphics
  extern int do_display_lists;
  if (do_display_lists && !dl_needs_refresh && 
      (display_mode == last_display_mode) && (rendermode == GL_RENDER)) {
    glCallList(dl_id);
    return;
  }
    
  // if never here before, create a new list
  if (dl_id == -1) {

          dl_id = glGenLists(1);

  }
  else {

          glDeleteLists(dl_id, 1);

          // dl_id = glGenLists(1); Chris: Is this needed?

  }


  // setup for DL list refresh

  if (do_display_lists) {
        glNewList(dl_id, GL_COMPILE_AND_EXECUTE);

  }
  dl_needs_refresh = 0;
  last_display_mode = display_mode;
  if (debug) cerr << "Object[" << getName() << "- regen display-list\n";
    
  // since have updated node positions, recalculate our stats here
  // (we used to do it on every iteration, but that was overkill)
  Point3D mmin, mmax; getBoundingBox(&mmin, &mmax);
  nodearray.computeStats();
    
  // if doing selection, load up the object number on the name stack
  if (rendermode == GL_SELECT)
    glPushName(objectarray_p->getIndex(this));
    
  // draw special things- labels, normals, etc
  {
    glDisable(GL_LIGHTING);
    glShadeModel(GL_FLAT);
        
    // if to draw normals on nodes
    if (display_mode.draw_nodenormals) 
      { nodearray.drawnormals(); }
        
    // if to draw normals on faces
    if (display_mode.draw_facenormals) 
      { facearray.drawnormals(); }
        
    // if to draw mesh over polys
    // if to draw normals on faces
    if (display_mode.draw_tetranormals) 
      { tetraarray.drawnormals(); }

    // if to draw mesh over polys
    if ((type == faces_only) && display_mode.draw_wireover) { 
      glColor3d(0,0, 0); 
      edgearray.drawover(drawover_amount); 
    }
                
    // if to draw mesh over polys
    if ((type == tetras_only) && display_mode.draw_wireover) { 
      glColor3d(0, 0, 0); 
                                                
      for (int i=0; i<facearray.getNumFaces(); i++) {
        Face* myface = facearray.getFace(i);
        if (myface->getTetraLink(1) == NULL) {  
          for (int j=0; j<3; j++) {
            Edge* myedge = myface->getEdgeLink(j);
            if (myedge) myedge->drawover(drawover_amount);
          }
        }
      } 
    }        
    // if to draw initial position
    if (display_mode.draw_initial) 
      { edgearray.drawinitial(); nodearray.drawinitial(); }
        
    // if to draw initial position
    if (display_mode.draw_collisions) 
      { facearray.drawmarkers(); }
        
    // if to draw labels for nodes
    if (display_mode.draw_nodelabels) 
      { glColor3d(0.5, 1.0, 1.0); nodearray.drawlabels(label_offset); }
        
    // if to draw labels for edges
    if (display_mode.draw_edgelabels) 
      { glColor3d(1.0, 1.0, 0.5); edgearray.drawlabels(label_offset); }
        
    // if to draw labels for faces
    if (display_mode.draw_facelabels) 
      { glColor3d(0.5, 1.0, 0.5); facearray.drawlabels(label_offset); } 
    
    // if to draw labels for faces
    if (display_mode.draw_tetralabels) 
      { glColor3d(0.0, 1.0, 1.0); tetraarray.drawlabels(label_offset); } 
  
    // if to draw labels for objects
    if (display_mode.draw_objectlabels) 
      { glColor3d(1.0, 0.0, 0.0); DrawLabel(); }
  }
    
  // do textures if supposed to
  extern int do_textures;
  if (do_textures && texture && (type == faces_only) &&
      ((display_mode.mode == DisplayMode::flat) || 
       (display_mode.mode == DisplayMode::smooth))){
    // if we have a texture, deal with it here
    glEnable(GL_TEXTURE_2D);
        
    // 4-byte-wise alignment (word-align) for speed
    glPixelStorei(GL_UNPACK_ALIGNMENT, 4);
        
    // setup texturing parameters
    //MODULATE,BLEND,DECAL,REPLACE
    if (texture_mode == Object::texture_env_mode(0))
      glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
    else if (texture_mode == texture_env_mode(1))
      glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_BLEND);
    else if (texture_mode == Object::texture_env_mode(2))
      glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL);
    else if (texture_mode == Object::texture_env_mode(3))
      glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
        
    if (texture_wrap == Object::clamp) {        
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
    }
    else {
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
    }
        
    if (texture_filter == Object::linear) {
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    }
    else {
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    }
        
    // debugging if we need it
    if (debug) 
      cerr << "applying " << texture_size.x << "x" << texture_size.y 
           << " texture...\n";
        
    /*********************************************/
    // if texture is not a power of 2, crop out the nearest part that is
    int crop_w = 1,crop_h = 1;
    int grab_w = 256,grab_h = 256;
    if ((texture_size.x == 64) || (texture_size.x == 128) || 
                (texture_size.x == 256) || (texture_size.x == 512) || (texture_size.x == 1024))
        {
                grab_w = texture_size.x;
                crop_w = 0;
        }
    if ((texture_size.y == 64) || (texture_size.y == 128) || 
                (texture_size.y == 256) || (texture_size.y == 512) || (texture_size.y == 1024))
        {
                grab_h = texture_size.y;
                crop_h = 0;
        }
        
    unsigned char* texture2 = NULL;
    int d = texture_size.z;
    if (crop_w*crop_h) {
                cerr << "cropping texture to be power of 2...\n";
                if (!texture2) {
                        texture2 = (unsigned char*)malloc(grab_w*grab_h*d*sizeof(unsigned char));
                }
                
                int texWidth = (grab_w > texture_size.x ? texture_size.x : grab_w);
                int texHeight = (grab_h > texture_size.y ? texture_size.y : grab_h);
                
                for (int count = 0; count < texHeight; count++)
                        memcpy((char *)texture2 + count*grab_w*d, (char *)texture + (count)*texture_size.x*d, texWidth*d);
    }
    else 
                texture2 = (unsigned char*)texture;
        
    switch (d) {
    case 1:
      glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE, grab_w, 
                   grab_h, 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, texture2);
      break;
    case 3:
      glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, grab_w, 
                   grab_h, 0, GL_RGB, GL_UNSIGNED_BYTE, texture2);
      break;
    case 4:
      glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, grab_w, 
                   grab_h, 0, GL_RGBA, GL_UNSIGNED_BYTE, texture2);
      break;
    default:
      break;
    }
    //  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, grab_w, 
    //          grab_h, 0, GL_RGBA, GL_UNSIGNED_BYTE, texture);
        
        
  }
  /********************************************************************/
  // if we should just display nodes, do so
  switch (type) {
  case nodes_only: {
    glDisable(GL_LIGHTING);
    nodearray.DrawNodes();
    break;
  }
  case edges_only: {
    // here we add texture to the thread object -- must figure out
    // why this interferes with plane texture
    if (do_textures && texture) {
      glEnable(GL_TEXTURE_2D);
      glPixelStorei(GL_UNPACK_ALIGNMENT, 4);
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
      glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
      glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, texture_size.x, 
                   texture_size.y,0, GL_RGBA, GL_UNSIGNED_BYTE, texture);

      // should think about mipmaps, like the following but use
      // glbindtextures ???
      //glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, 
      //                GL_NEAREST_MIPMAP_NEAREST);
      //gluBuild2DMipmaps(GL_TEXTURE_2D,GL_RGBA,texture_size.x,
      //                  texture_size.y,GL_RGBA,GL_UNSIGNED_BYTE,texture);
      glDisable(GL_LIGHTING);
    }
    else { // do colors/material properties
      // set front face properties
      glEnable(GL_LIGHTING);
      glShadeModel(GL_FLAT);
      glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient);
      glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
      glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
      glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
      // set back face properties
      glMaterialfv(GL_BACK, GL_AMBIENT, back_mat_ambient);
      glMaterialfv(GL_BACK, GL_DIFFUSE, back_mat_diffuse);
      glMaterialfv(GL_BACK, GL_SPECULAR, back_mat_specular);
      glMaterialfv(GL_BACK, GL_SHININESS, back_mat_shininess);  
    }
    edgearray.DrawEdges();
    glDisable(GL_TEXTURE_2D);
    if (display_mode.mode == DisplayMode::wireframe)
      nodearray.DrawNodes();
    break;
  }
  case faces_only: {
    // recompute normals (since surface is changing)
    facearray.computeNormals();
            
    // average face normals to get per-vertex normals
    nodearray.interpolateNormals();
            
    // depending on the display mode
    if (display_mode.mode == DisplayMode::wireframe) {
      glDisable(GL_LIGHTING);
      edgearray.DrawEdges();
      nodearray.DrawNodes();
    }
            
    // for flat or solid shaded modes
    if ((display_mode.mode == DisplayMode::flat) || 
        (display_mode.mode == DisplayMode::smooth)) {
      // turn on lights
      glEnable(GL_LIGHTING);
                
      // setup shading
      if (display_mode.mode == DisplayMode::flat) glShadeModel(GL_FLAT);
      else glShadeModel(GL_SMOOTH);
                
      { // do colors/material properties
        // set front face properties
        glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient);
        glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
        glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
        glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
                    
        // set back face properties
        glMaterialfv(GL_BACK, GL_AMBIENT, back_mat_ambient);
        glMaterialfv(GL_BACK, GL_DIFFUSE, back_mat_diffuse);
        glMaterialfv(GL_BACK, GL_SPECULAR, back_mat_specular);
        glMaterialfv(GL_BACK, GL_SHININESS, back_mat_shininess);        
      }
                
      // do the draw
      facearray.DrawFaces(int(do_textures && texture));
    }
    break;
  }
  case tetras_only: {
    if (dynamics == deformable) {
      // recompute normals (since surface is changing)
      facearray.computeNormals();

      // average face normals to get per-vertex normals
      nodearray.interpolateTetraNormals();
    }

    // depending on the display mode
    if (display_mode.mode == DisplayMode::wireframe) {
      glDisable(GL_LIGHTING);
      edgearray.DrawEdges();
      nodearray.DrawNodes();
    }

    // for flat or solid shaded modes
    if ((display_mode.mode == DisplayMode::flat) || 
        (display_mode.mode == DisplayMode::smooth)) {
      // turn on lights
      glEnable(GL_LIGHTING);

      // setup shading
      if (display_mode.mode == DisplayMode::flat) glShadeModel(GL_FLAT);
      else glShadeModel(GL_SMOOTH);

      { // do colors/material properties
        glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient);
        glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
        glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
        glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);

        glColor4f(mat_ambient[0],mat_ambient[1],mat_ambient[2],mat_ambient[3]);

        glMaterialfv(GL_BACK, GL_AMBIENT, back_mat_ambient);
        glMaterialfv(GL_BACK, GL_DIFFUSE, back_mat_diffuse);
        glMaterialfv(GL_BACK, GL_SPECULAR, back_mat_specular);
        glMaterialfv(GL_BACK, GL_SHININESS, back_mat_shininess);                                                
      }

      // do the draw
      tetraarray.DrawTetras(int(do_textures && texture));
    }
    break;
  }
  } // end switch(type)
        
  // finish up - other objects may not be textured, so turn off
  glDisable(GL_TEXTURE_2D);    
        
  { 

        // draw tie nodes if set.
        if (display_mode.draw_tienodes) {

          nodearray.drawtienodes();

        }

    // if to draw bounding sphere
    if (display_mode.draw_boundingspheres) {

                // check what type of collision bounds to draw

                extern springCore* springApp;

                switch (springApp->getCollisionAlgorithm()) {



                        // BOUNDING_SPHERE type collision detection

                        case springCore::BOUNDING_SPHERE_CD: {

                                if (boundingSphereRoot) {
                                        boundingSphereRoot->drawSubtree();

                                }

                                break;

                        }



                        // AABB type collision detection (Axis-Aligned Bounding-Boxes)

                        case springCore::AABB_CD: {

                                if (boundingBoxRoot) {

                                        boundingBoxRoot->draw();

                                }

                                break;

                        }

                }
        }


    //draw principal axis
    if (display_mode.draw_principalaxis) {
      getNodeArray()->PrincipalAxes();
      DrawPrincipalAxes();
    }
  }
        
  // if doing selection, pop the object number off the name stack
  if (rendermode == GL_SELECT)
    glPopName();
        
  // end of display list

  if (do_display_lists) {
        glEndList();

  }
}

// this draws a wireframe bounding box around the object
void Object::DrawBoundingBox()
{
    Point3D min, max;
    getBoundingBox(&min, &max);
    glDisable(GL_LIGHTING);
    glColor3f(1.0, 0.0, 0.0);  // red for now
    // front face
    glBegin(GL_LINE_LOOP);
    glVertex3d(min.x, min.y, min.z);
    glVertex3d(max.x, min.y, min.z);
    glVertex3d(max.x, max.y, min.z);
    glVertex3d(min.x, max.y, min.z);
    glEnd();
    // back face
    glBegin(GL_LINE_LOOP);
    glVertex3d(min.x, min.y, max.z);
    glVertex3d(max.x, min.y, max.z);
    glVertex3d(max.x, max.y, max.z);
    glVertex3d(min.x, max.y, max.z);
    glEnd();
    // connecting edges
    glBegin(GL_LINES);
    glVertex3d(min.x, min.y, min.z);
    glVertex3d(min.x, min.y, max.z);
    glVertex3d(max.x, min.y, min.z);
    glVertex3d(max.x, min.y, max.z);
    glVertex3d(max.x, max.y, min.z);
    glVertex3d(max.x, max.y, max.z);
    glVertex3d(min.x, max.y, min.z);
    glVertex3d(min.x, max.y, max.z);
    glEnd();
}

void Object::DrawLabel()
{
    // draw object name in center of object
    extern RasterFont* font_p;
    Point3D center = getCenter();
    font_p->printString(getName(), center.x,center.y,center.z);
}

void Object::DrawPrincipalAxes()
{
        //draw princial axes.  Note that higher values of sig[] contain the 
        //relative importance of principal axis, ie the bigger one. could draw 
        //in this order if desired but not set up to do so now
  
        NodeArray       *na = getNodeArray(); 
        Point3D bmin, bmax;
        getBoundingBox(&bmin, &bmax);
        double temp1 = MAX3(fabs(bmin.x), fabs(bmin.y), fabs(bmin.z));
        double temp2 = MAX3(fabs(bmax.x), fabs(bmax.y), fabs(bmax.z));
        double max = MAX(temp1, temp2);  // size of axes
        double xx, xy, xz, yx, yy, yz, zx, zy, zz;

        double U[3][3];  //matrix contains the eigen vectors that correspond to P.A.
        na->getUmatrix(U);

        //assign these vectors to axis directions
        xx=U[0][0];
        xy=U[1][0];
        xz=U[2][0];
        yx=U[0][1];
        yy=U[1][1];
        yz=U[2][1];
        zx=U[0][2];
        zy=U[1][2];
        zz=U[2][2];

  glDisable(GL_LIGHTING);
  glBegin(GL_LINES);
  glColor3f(1.0, 0.0, 0.0);
  glVertex3d(-max*xx, -max*xy, -max*xz);
  glVertex3d(max*xx, max*xy, max*xz);
  glColor3f(1.0, 1.0, 0.0);
  glVertex3d(-max*yx, -max*yy, -max*yz);
  glVertex3d(max*yx, max*yy, max*yz);
  glColor3f(0.0, 1.0, 0.0);
  glVertex3d(-max*zx, -max*zy, -max*zz);
  glVertex3d(max*zx, max*zy, max*zz);
  glEnd();
}

void Object::GeometryReplicate(GeometryReplicator* gr)
{
  if (debug) 
          cerr << "Object[" << getName() << "]::GeometryReplicate(" << gr << ")\n";

  // sanity check
  if (!gr) return; 

  // if nothing new, then nothing to update!
  if (!dl_needs_refresh) return;

  // if currently being modified, don't send- 
  // let them use the geometry they've got
  if (being_modified) return;

  // message buffer + locals
  const int buf_size = 1024;
  char buffer[buf_size];
  Timer timer;

  { // first, send our object ID
    int s = sprintf(buffer, "%c%c%s", GR_CODE, GR_OBJECT, getName());
    gr->SendMessage(buffer, s);
  }

  { // next, send the object attributes 
    // (visibility, ambient, diffuse, specular, shininess)
        unsigned char* ubuffer = (unsigned char*)buffer;
    int n = 0;
    ubuffer[n++] = GR_CODE; ubuffer[n++] = GR_OBJ_ATTR;
        ubuffer[n++] = isVisible();
    { for (int i=0; i<4; i++) ubuffer[n++] = int(mat_ambient[i] * 255); }
    { for (int i=0; i<4; i++) ubuffer[n++] = int(mat_diffuse[i] * 255); }
    { for (int i=0; i<4; i++) ubuffer[n++] = int(mat_specular[i] * 255); }
    { for (int i=0; i<4; i++) ubuffer[n++] = int(mat_ambient[i] * 255); }
    ubuffer[n++] = int(mat_shininess[0] * 255);

    // (do we need to send the number of nodes/faces as an attribute to make 
    // it easier for the other side?)

    gr->SendMessage((char*)ubuffer, n);
  }

  // next send our texture (if we have one)
  // other things to add- texture modes- clamp/repeat, etc
  if (debug) { 
          cerr << "\ttime to send attributes = " << timer.GetElapsedTime() << endl;
          timer.Reset();
  }
  if (texture) {
        unsigned char* ubuffer = (unsigned char*)buffer;
    int n = 0;
        // encode the texture header
    ubuffer[n++] = GR_CODE; ubuffer[n++] = GR_TEXTURE;
        ubuffer[n++] = texture_size.x / 255; ubuffer[n++] = texture_size.x % 255;
        ubuffer[n++] = texture_size.y / 255; ubuffer[n++] = texture_size.y % 255;
        ubuffer[n++] = texture_size.z / 255; ubuffer[n++] = texture_size.z % 255;

        // and the texture data itself
    int size = texture_size.x * texture_size.y * texture_size.z;
        memcpy(&(buffer[n]), texture, size);
   gr->SendMessage(buffer, n+size);
  }
  if (debug) {
          cerr << "\ttime to send texture = " << timer.GetElapsedTime() << endl;
          timer.Reset();
  }

  // then send our nodes
  nodearray.GeometryReplicate(gr);
  if (debug) {
          cerr << "\ttime to send nodes = " << timer.GetElapsedTime() << endl;
          timer.Reset();
  }

  // get the nodearray to cache the index for each node to make it's index
  // lookup fast for the facearray function
  nodearray.CacheIndices();

  // and send our faces
  facearray.GeometryReplicate(gr);
  if (debug) {
          cerr << "\ttime to send faces = " << timer.GetElapsedTime() << endl;
          timer.Reset();
  }
}

// tests if this object and input object have intersecting bounding spheres
// (faster test than bounding box test)
int Object::intersectsBoundingSphere(Object* test_obj)
{
    Point3D our_center = getCenter();
    Point3D test_obj_center = test_obj->getCenter();
    double dist_between_centers = our_center.Dist(test_obj_center);
    double dist_between_bspheres = getRadius() + test_obj->getRadius();
    int retval = ( (dist_between_centers < dist_between_bspheres) );
    
    if (debug && retval)
        cerr << "[" << getName() << "] is inside bsphere of object ["
        << test_obj->getName() << "] dist = "
        << dist_between_centers << ", radii = "
        << getRadius() << " + " << test_obj->getRadius() 
        << "=" << dist_between_bspheres << endl;
    
    return(retval);
}

// tests if this object and input object have intersecting bounding boxes
// can grow each bounding box by some error e (defaults to 0)
int Object::intersectsBoundingBox(Object *obj, double e)
{
    Point3D thismin, thismax;
    Point3D objmin, objmax;
    this->getBoundingBox(&thismin, &thismax);
    obj->getBoundingBox(&objmin, &objmax);
    if (e > 0.0) { // grow the bounding boxes
        Point3D p(e, e, e);
        thismin -= p;
        thismax += p;
        objmin -= p;
        objmax += p;
    }
    if (((thismax.x < objmin.x) || (thismin.x > objmax.x)) || 
        ((thismax.y < objmin.y) || (thismin.y > objmax.y)) || 
        ((thismax.z < objmin.z) || (thismin.z > objmax.z)))
        return(0);  // they can't intersect if one of these conditions is true
    else
        return(1);
}

Node* Object::findClosestNode(Point3D pt)
{
    if (debug) cerr << "Object[" << name << "]::findClosestNode(" << pt << ")\n";
    
    if (nodearray.getNumNodes() == 0) 
        return(NULL);
    
    // ask the nodearray to do it for us- they have the info
//#define USE_SBBS
#ifdef USE_SBBS
    int closest_node_index = nodearray.findCollisionNode(pt);
#else
    int closest_node_index = nodearray.findClosestNode(pt);
#endif
    if (closest_node_index < 0) return(NULL);
    
    // return the closest one
    return( nodearray.getNode(closest_node_index) );
}

// Takes another object and returns the closest node in our mesh to them,
// the closest node in their mesh to us, and the distance
// Uses dumb, but correct brute force O(n^2) algorithm
double Object::findClosestNodesBetweenObjects(Object* other_obj, 
                                              Node** our_node_p, Node** other_node_p)
{
    if (debug) 
        cerr << "Object[" << name << "]::findClosestNodesBetweenObjects(" 
        << other_obj << "," << our_node_p << "," << other_node_p << ")\n";
    
    // ask nodearray to do it
    double dist = nodearray.findClosestNodesBetweenNodeArrays(
        other_obj->getNodeArray(), our_node_p, other_node_p);
    
    // finally, return the distance
    return( dist );
}

/* test with a 1-D object */
void Object::CreateTestLine(int num_segs, double length, Point3D offset,
                            int fix_end)
{
    type = edges_only;
    dynamics = deformable;
    behavior = none;
    uses_viagra = 0;
    set_default_material_properties();
    nodearray.reset();
    edgearray.reset();
    facearray.reset(); // don't need to allocate faces
    nodearray.allocateNodes(num_segs + 1);
    edgearray.allocateEdges(num_segs);
    
    double seg_length = length/num_segs;
    for (int i = 0; i < (num_segs + 1); i++) {
        Point3D pt(offset.x + (double)i*seg_length, offset.y, offset.z);
        nodearray.addNode(pt);
        if (i != 0)
            edgearray.addEdge(i, i-1);
    }
    if (fix_end)
        nodearray.getNode(num_segs)->setFixedType(Node::POS);
    
    nodearray.computeStats();
}

// create a cloud of randomly placed nodes
void Object::CreateNodeCloud(int num_nodes, int radius, Point3D offset)
{
    type = nodes_only;
    dynamics = deformable;
    behavior = none;
    uses_viagra = 0;
    set_default_material_properties();
    nodearray.reset();
    edgearray.reset();
    facearray.reset(); // don't need to allocate faces
    nodearray.allocateNodes(num_nodes);
    
    for (int i = 0; i < num_nodes; i++) {
        Point3D pt(rand() % radius, rand() % radius, rand() % radius);
        pt += offset;
        nodearray.addNode(pt);
    }
    
    nodearray.computeStats();
}

// triangular equilateral mesh
void Object::CreateEquiMesh(int numlevels, double edgelength)
{
    // set the type, dynamics, and material of object
    type = faces_only;
    dynamics = deformable;
    behavior = none;
    uses_viagra = 0;
    set_default_material_properties();
    
    if (numlevels < 1) numlevels = 1;
    // allocate enough nodes, edges, and faces
    int projected_num_nodes = (numlevels + 1)*(numlevels + 2)/2;
    nodearray.reset();
    edgearray.reset();
    facearray.reset();
    projected_num_nodes *= 2; // need safety factor for adding nodes, etc.
    nodearray.allocateNodes(projected_num_nodes);
    edgearray.allocateEdges(projected_num_nodes * 3);
    facearray.allocateFaces(projected_num_nodes * 2);
    
    double cos30 = cos(PI/6.0);
    int node_index = 0;
    for (int i = 0; i <= numlevels; i++) {
        for (int j = 0; j <= i; j++) {
            // add the nodes
            Point3D pt(i*edgelength*cos30, i*edgelength*0.5 - j*edgelength, 0.0);
            nodearray.addNode(pt);
            // fix if in last level
            if (i == numlevels)
                nodearray.getNode(node_index)->setFixedType(Node::POS);
            // and add faces when their 3 nodes all exist
            if (i > 0) {
                if (j > 0 && j < i) {
                    facearray.addFace(node_index, node_index - i - 1, node_index - 1);
                    facearray.addFace(node_index, node_index - i, node_index - i - 1);
                }
                else if (j == i) {
                    facearray.addFace(node_index, node_index - i - 1, node_index - 1);
                }
            }
            node_index++;
        }
    }
}

/*
* orientation: 'x' for length going in direction of x axis, 'y' for y, etc
* minimal tube will have rings of edges, connected with straight edges to 
* adjacent rings.  
* add one diagonal between rings to get triangles (diag1)
* add another to have symmetry (diag2)
* add a smaller tube inside with straight edges to first tube (inner)
* add cross links between the 2 layers in adjacent rings (cross1)
* add cross links between the 2 layers in the same ring (cross2)
* rl = -1, restLength default, sc1 (radial) sc2 (lateral) springConstants 
* rl_out = resting length of the springs on the outer cylinder
* rl_in = resting length of the springs on the inner cylinder
* rl_betw = resting length of the springs linking the 2 cylinders
*/
void Object::CreateTube(char orientation, int num_radii, int num_rings, double radius, 
                        double length, Point3D offset, int diag1, int diag2,
                        int fix_num, double inner_radius, 
                        int cross1, int cross2, int capped,
                        double sc1, double sc2, double rl_out,
                        double rl_in, double rl_betw)
{
    if (debug) cerr << "Object[" << name << "]::CreateTube(" << num_radii 
        << ", " << radius << ", " << inner_radius << ", "
        << length << ", " << offset << ")\n";
    
    // set the type, dynamics, and material of object
    type = faces_only;
    dynamics = deformable;
    behavior = none;
    uses_viagra = 0;
    set_default_material_properties();
    
    // allocate enough nodes, edges, and faces
    int projected_num_nodes = num_rings * num_radii * 2; // *2 in case inner hull
    nodearray.reset();
    edgearray.reset();
    facearray.reset();
    projected_num_nodes *= 2; // need safety factor for adding nodes, etc.
    nodearray.allocateNodes(projected_num_nodes);
    edgearray.allocateEdges(projected_num_nodes * 6);
    facearray.allocateFaces(projected_num_nodes * 2);
    
    int node_index = 0;
    
    // avoid division by 0:
    int xdivide = ((num_rings == 1) ? 1 : (num_rings - 1));
    // outer hull
    {
        for (int ring = 0; ring < num_rings; ring++) { // for each ring
            for (int spoke = 0; spoke < num_radii; spoke++) { // for each spoke
                // compute a node
                double angle = spoke * 2*PI/num_radii;
                Point3D pt;
                
                // use orientation to figure out how to map the values
                switch (orientation) {
                default:
                case 'x' : pt = Point3D(ring*length/(xdivide),
                               radius*cos(angle),-radius*sin(angle));
                    break;
                case 'y' : pt = Point3D(radius*cos(angle),
                               ring*length/(xdivide),-radius*sin(angle));
                    break;
                case 'z' : pt = Point3D(radius*cos(angle),
                               -radius*sin(angle),ring*length/(xdivide));
                    break;
                } // end switch
                
                // add in the offset and create the node at this point
                pt += offset;
                nodearray.addNode(pt);
                
                // caveat: make sure nodes exist before adding edges 
                // and edges exist before adding faces
                
                // add edges in a ring
                if (spoke != 0)
                    edgearray.addEdge(node_index, (node_index - 1), rl_out, sc1);
                if (spoke == (num_radii - 1))
                    edgearray.addEdge(node_index, (node_index - num_radii + 1), rl_out, sc1);
                
                // add straight edges between adjacent rings
                if (ring != 0)
                    edgearray.addEdge(node_index, (node_index - num_radii), rl_out, sc2);
                
                // add one diagonal between adjacent rings
                if (diag1 && ring != 0) {
                    if (spoke != 0)
                        edgearray.addEdge(node_index,(node_index - num_radii - 1),rl_out,sc2);
                    else
                        edgearray.addEdge(node_index, (node_index - 1), rl_out, sc2);
                }
                
                // add other diagonal between adjacent rings
                if (diag2 && ring != 0) {
                    if (spoke != (num_radii - 1))
                        edgearray.addEdge(node_index,(node_index - num_radii + 1),rl_out,sc2);
                    else
                        edgearray.addEdge(node_index,(node_index -2*num_radii + 1),rl_out,sc2);
                }
                
                // add faces using the first diagonal, if we have it
                if (diag1 && ring != 0) {
                    int other = node_index - num_radii; // neighbor in previous ring
                    if (spoke != 0) {
                        facearray.addFace(node_index, other, (other - 1));
                        facearray.addFace(node_index, (other - 1), (node_index - 1));
                    }
                    if (spoke == 0)
                        facearray.addFace(node_index, other, (node_index - 1));
                    if (spoke == (num_radii - 1))
                        facearray.addFace(node_index, (other + 1), other);
                }
                
                // fix some rings on the end of the tube
                if ((((num_rings - ring) <= fix_num) && (offset.x > 0)) ||
                    ((ring < fix_num) && (offset.x < 0))) {
                    nodearray.getNode(node_index)->setFixedType(Node::POS);
                }
                
                node_index++;
            }
        }
    }
    
    // inner hull
    if (inner_radius > 0.0) {
        int outer_num = node_index;  // number of nodes in outer hull
        radius = inner_radius;
        
        for (int ring = 0; ring < num_rings; ring++) { // for each ring
            for (int spoke = 0; spoke < num_radii; spoke++) { // for each spoke
                // compute a node
                double angle = spoke * 2*PI/num_radii;
                Point3D pt;
                
                // use orientation to figure out how to map the values
                switch (orientation) {
                default:
                case 'x' : pt = Point3D(ring*length/(xdivide),
                               radius*cos(angle),-radius*sin(angle));
                    break;
                case 'y' : pt = Point3D(radius*cos(angle),
                               ring*length/(xdivide),-radius*sin(angle));
                    break;
                case 'z' : pt = Point3D(radius*cos(angle),
                               -radius*sin(angle),ring*length/(xdivide));
                    break;
                } // end switch
                
                // add in the offset and create the node at this point
                pt += offset;
                nodearray.addNode(pt);
                
                // caveat: make sure nodes exist before adding edges 
                // and edges exist before adding faces
                
                // add edges in a ring
                if (spoke != 0)
                    edgearray.addEdge(node_index, (node_index - 1), rl_in, sc1);
                if (spoke == (num_radii - 1))
                    edgearray.addEdge(node_index, (node_index - num_radii + 1), rl_in, sc1);
                
                // add straight edges between adjacent rings
                if (ring != 0)
                    edgearray.addEdge(node_index, (node_index - num_radii), rl_in, sc2);
                
                // add one diagonal between adjacent rings
                if (diag1 && ring != 0) {
                    if (spoke != 0)
                        edgearray.addEdge(node_index,(node_index - num_radii - 1),rl_in,sc2);
                    else
                        edgearray.addEdge(node_index, (node_index - 1), rl_in, sc2);
                }
                
                // add other diagonal between adjacent rings
                if (diag2 && ring != 0) {
                    if (spoke != (num_radii - 1))
                        edgearray.addEdge(node_index,(node_index - num_radii + 1),rl_in,sc2);
                    else
                        edgearray.addEdge(node_index,(node_index -2*num_radii + 1),rl_in,sc2);
                }
                
                // add radial edge connecting inner and outer hull
                edgearray.addEdge(node_index, node_index - outer_num, rl_betw, sc1);
                
                // add cross members between corresponding inner and outer nodes of
                // adjacent rings
                if (cross1) {
                    int other = node_index - outer_num;
                    if (ring != (num_rings - 1))
                        edgearray.addEdge(node_index, (other + num_radii), rl_betw, sc2);
                    if (ring != 0)
                        edgearray.addEdge(node_index, (other - num_radii), rl_betw, sc2);
                }
                
                // add cross members between inner and outer nodes of the same ring
                if (cross2) {
                    int other = node_index - outer_num;
                    if (spoke != 0)
                        edgearray.addEdge(node_index, (other - 1), rl_betw, sc1);
                    else
                        edgearray.addEdge(node_index, (other + num_radii - 1), rl_betw, sc1);
                    if (spoke != (num_radii - 1))
                        edgearray.addEdge(node_index, (other + 1), rl_betw, sc1);
                    else
                        edgearray.addEdge(node_index, (other - num_radii + 1), rl_betw, sc1);
                }
                
                // ??? could eliminate faces on inner tube
                // add faces using the first diagonal, if we have it
                /*if (diag1 && ring != 0) {
                int other = node_index - num_radii; // neighbor in previous ring
                if (spoke != 0) {
                facearray.addFace(node_index, other, (other - 1));
                facearray.addFace(node_index, (other - 1), (node_index - 1));
                }
                if (spoke == 0)
                facearray.addFace(node_index, other, (node_index - 1));
                if (spoke == (num_radii - 1))
                facearray.addFace(node_index, (other + 1), other);
            }*/
                
                // add faces to cap walls at the ends
                // NB: this will add (or use) one of the cross2 links
                if (capped) {
                    int other = node_index - outer_num;
                    if (ring == 0) {
                        if (spoke != 0) {
                            facearray.addFace(node_index, (node_index - 1), (other - 1));
                            facearray.addFace(node_index, (other - 1), other);
                        }
                        if (spoke == 0)
                            facearray.addFace(node_index, (other + num_radii - 1), other);
                        if (spoke == (num_radii - 1))
                            facearray.addFace(node_index, other,(node_index - num_radii +1));
                    }
                    if (ring == (num_rings -1)) {  // make normals go the other way
                        if (spoke != 0) {
                            facearray.addFace(node_index, (other - 1), (node_index - 1));
                            facearray.addFace(node_index, other, (other - 1));
                        }
                        if (spoke == 0)
                            facearray.addFace(node_index, other, (other + num_radii - 1));
                        if (spoke == (num_radii - 1))
                            facearray.addFace(node_index,(node_index - num_radii +1),other);
                    }
                }
                
                // fix some rings on the end of the tube
                if ((((num_rings - ring) <= fix_num) && (offset.x > 0)) ||
                    ((ring < fix_num) && (offset.x < 0))) {
                    nodearray.getNode(node_index)->setFixedType(Node::POS);
                }
                
                node_index++;
      }
    }
  }
  
  // cap the ends (with many triangles all meeting in the center) even 
  // without an inner cylinder
  if (inner_radius == 0.0 && capped) {
      int outer_num = node_index;
      Point3D center0 = offset;
      Point3D center1 = offset + (nodearray.getNode(outer_num - num_radii)->p - 
          nodearray.getNode(0)->p);
      
      int i;
      nodearray.addNode(center0);
      for (i = 0; i < num_radii; i++) {
          facearray.addFace(node_index, i, (i+1)%num_radii);
      }
      node_index++;
      
      if (num_rings > 1) {
          nodearray.addNode(center1);
          for (i = 0; i < num_radii; i++) {
              facearray.addFace(node_index, (i+1)%num_radii + outer_num - num_radii,
                  i + outer_num - num_radii);
          }
          node_index++;
      }
  }
  
  // tell nodearray to recalculate bounding boxes, etc
  nodearray.computeStats();
  
  // recompute surface normals (since is a new surface)
  facearray.computeNormals();
  
  // average face normals to get per-vertex normals
  nodearray.interpolateNormals();
}

// Create a triangulated sphere with given radius
int Object::CreateSphere(int numlayers, double radius, double inner)
{
    // set the type, dynamics, and material of object
    type = faces_only;
    dynamics = deformable;
    behavior = none;
    uses_viagra = 0;
    set_default_material_properties();
    
    if (numlayers < 2) numlayers = 2;
    // allocate enough nodes, edges, and faces
    int projected_num_nodes = 2000;
    nodearray.reset();
    edgearray.reset();
    facearray.reset();
    nodearray.allocateNodes(projected_num_nodes);
    edgearray.allocateEdges(projected_num_nodes * 3);
    facearray.allocateFaces(projected_num_nodes * 2);
    
    Point3D pt;
    int shells = 1;
    int top = 0;
    int outer = 0;
    if (inner > 0.0)
        shells = 2;
    for (int k = 0; k < shells; k++) {
        if (k == 1)
            radius = inner;
        //
        // first add all the nodes
        //
        { // northern hemisphere
            for (int layer = 0; layer < numlayers; layer++) { 
                if (layer == 0) { // do north pole
                    pt = Point3D(0, radius, 0);
                    nodearray.addNode(pt);
                }
                else {
                    double h = radius*sin((PI/2)*(numlayers-1-layer)/(numlayers-1));
                    double r = sqrt(radius*radius - h*h);
                    int n = 5*(int)(pow((double) 2, (layer-1)));        // SYK - 5/09/06 for VC2005
                    for (int i = 0; i < n; i++) {
                        double angle = i*2*PI/n;
                        pt = Point3D(r*sin(angle), h, r*cos(angle));
                        nodearray.addNode(pt);
                    }
                }
            }
        }
        
        if (k == 0)
            top = nodearray.getNumNodes(); 
        
        {  // southern hemisphere
            for (int layer = 0; layer < (numlayers - 1); layer ++) { 
                if (layer == 0) { // south pole (fix it)
                    pt = Point3D(0, -radius, 0);
                    nodearray.addNode(pt);
                    nodearray.getNode(top)->setFixedType(Node::POS);
                }
                else {
                    double h = radius*sin((PI/2)*(numlayers-1-layer)/(numlayers-1));
                    double r = sqrt(radius*radius - h*h);
                    int n = 5*(int)(pow((double) 2, (layer-1)));  // SYK - 5/09/06 for VC2005
                    for (int i = 0; i < n; i++) {
                        double angle = i*2*PI/n;
                        pt = Point3D(r*sin(angle), -h, r*cos(angle));
                        nodearray.addNode(pt);
                    }
                }
            }
        }
        
        // and add some edges between inner/outer sphere
        if (k == 1) {
            for (int m = 0; m < outer; m++)
                edgearray.addEdge(m, outer+m);
        }
        // now add all the faces
        // 5 triangles for poles
        for (int i = 0; i < 5; i++) {
            if (numlayers > 2) {
                facearray.addFace(outer,outer+i+1, outer+(i+1)%5 + 1);
                facearray.addFace(outer+top, outer+(i+1)%5 +1+top, outer+i+1 +top);
            }
            else { // if this is the final layer
                facearray.addFace(outer, outer+i+1, outer+(i+1)%5 + 1);
                facearray.addFace(outer+top, outer+(i+1)%5 + 1, outer+i+1);
            }
        }
        int min = outer; int min2;
        int max = outer; int max2;
        int fl;
        for (int layer = 1; layer < (numlayers - 1); layer++) {
            min = max + 1;                   // min node in this layer
            max = max + 5*(int)(pow((double) 2,(layer-1)));  // max node in this layer // SYK - 5/09/06 for VC2005
            min2 = max + 1;                  // min node in next layer
            max2 = max + 5*(int)(pow((double) 2,layer));     // max node in next layer // SYK - 5/09/06 for VC2005
            if (layer == (numlayers-2)) // flag for final layer of south
                fl = 0;
            else
                fl = top;
            int n = 5*(int)(pow((double) 2, (layer-1)));        // SYK - 5/09/06 for VC2005
            for (int i = 0; i < n; i++) {
                // need to add 3 triangles for each node in interior layer
                facearray.addFace(min + i, min2 + 2*i, min2 + 2*i +1);
                if (i == (n - 1))
                    facearray.addFace(min + i, min2 + 2*i + 1, min);
                else
                    facearray.addFace(min + i, min2 + 2*i + 1, min +i+1);
                if (i == 0)
                    facearray.addFace(min + i, max2, min2 + 2*i);
                else
                    facearray.addFace(min + i, min2 + 2*i - 1, min2 + 2*i);
                // and again for southern hemisphere
                facearray.addFace(min +i+top, min2 + 2*i +1+fl, min2 + 2*i +fl);
                if (i == (n - 1))
                    facearray.addFace(min +i+top, min+top, min2 + 2*i + 1+fl);
                else
                    facearray.addFace(min +i+top, min +i+1+top, min2 + 2*i + 1+fl);
                if (i == 0)
                    facearray.addFace(min +i+top, min2 + 2*i +fl, max2 +fl);
                else
                    facearray.addFace(min +i+top, min2 + 2*i +fl, min2 + 2*i -1+fl);
            }
        }
        outer = nodearray.getNumNodes(); printf("%d %d\n", top, outer);
  } // for k -- inner shell
  // tell nodearray to recalculate bounding boxes, etc
  nodearray.computeStats();
  
  // recompute surface normals (since is a new surface)
  facearray.computeNormals();
  
  // average face normals to get per-vertex normals
  nodearray.interpolateNormals();
  return (top); // might want to fix the south pole
}

// creates an articulate thread of given length (in mm), with nodes spaced
// X mm apart, starting at given offset (orientation: x=0,y=1,z=2)
void Object::CreateThread(double length, double node_spacing, int orientation,
                          Point3D offset, double edgethresh_in, int first_seg)
{
    if (debug) cerr << "Object[" << name << "]::CreateThread(" << length << 
        ", " << node_spacing << ", " << offset << ")\n";
    
    // set the type of object
    type = edges_only;
    // set the thickness of the edges
    edgethresh = edgethresh_in;
    dynamics = articulate;
    //dynamics = deformable;
    nummethod = quasi;
    behavior = none;
    uses_viagra = 0;
    set_default_material_properties(); // setup material properties appropriately
    
    // will turn during knot tie (iff we have >=2 objects colliding)
    //do_internal_collisions = 1;                
    
    int num_nodes = (int)(length/node_spacing) + 1;

    // calculate how many nodes and edges there will be and allocate enough
    nodearray.reset();
    edgearray.reset();
    facearray.reset();
    nodearray.allocateNodes(2*num_nodes);  // safety factor
    edgearray.allocateEdges(2*num_nodes);  // #edges = #nodes (well, -1)
    facearray.allocateFaces(10);           // there aren't any faces
    
    // create the thread
    int last_node_index = -1;
    for (int i=0; i < num_nodes; i++) {
        Point3D pt = offset;
        switch (orientation) {
        case 0:  pt.x += i*node_spacing; break; // x-direction
        case 1:  pt.y += i*node_spacing; break; // y-direction
        case 2:  pt.z += i*node_spacing; break; // z-direction
        }
        int node_index = nodearray.addNode(pt);
        //Node* node = nodearray.getNode(node_index);
        //node->setMass(100.0);
        
        if (last_node_index == -1)
            last_node_index = node_index;
        else {
            // add in this edge
            int edge_index = 
                edgearray.addEdge(last_node_index, node_index, -1, 20, 40);
            edgearray.getEdge(edge_index)->setSpringConstant(2.0);
            
            // save the back pointer
            last_node_index = node_index;
        }
        if (i == 0)
          i += first_seg;
        //if (i == 0|| i == 3)
        //i += 2; // make the first 2 edges 3 times as long, they're the needle!
    }
    
    // tell nodearray to recalculate bounding boxes, etc
    nodearray.computeStats();
    
    // recompute surface normals (since is a new surface)
    facearray.computeNormals();
    
    // average face normals to get per-vertex normals
    nodearray.interpolateNormals();
}

// creates an edge object in ring shape, for wrapping thread around
void Object::CreateEdgeRing(Point3D center, double radius, int numedges,
                            double edgethresh_in)
{
    if (debug) cerr << "Object[" << name << "]::CreateEdgeRing(" << center << 
        ", " << radius << ", " << numedges << ")\n";
    
    // set the type of object
    type = edges_only;
    // set the thickness of the edges
    edgethresh = edgethresh_in;
    dynamics = rigid;
    nummethod = quasi;
    behavior = none;
    uses_viagra = 0;
    set_default_material_properties(); // setup material properties appropriately
    
    int num_nodes = numedges;

    // calculate how many nodes and edges there will be and allocate enough
    nodearray.reset();