---

Timestep.C

Go to the documentation of this file.

/*
 * Copyright (C) 2003 Eric Bohm
 *
 *This program is free software; you can redistribute it and/or modify
 *it under the terms of the GNU General Public License as published by
 *the Free Software Foundation version 2 of the License.
 * 
 *This program is distributed in the hope that it will be useful, but
 *WITHOUT ANY WARRANTY; without even the implied warranty of
 *MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *General Public License for more details.
 * 
 *You should have received a copy of the GNU General Public License along
 *with this program; if not, write to the Free Software Foundation, Inc.,
 *675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */

using namespace std;
#include <GL/glut.h>
#include <stdio.h>
#include <vector>
#include <algorithm>
#include "Timestep.h"

/* clear GL, do rotation and scaling, draw tetrahedrons */

void OutputString(GLfloat x, GLfloat y, char *string)
{
  int length;

  length = strlen(string);
  glRasterPos2f(x, y);
  for ( int i = 0; i < length; i++)
    glutBitmapCharacter(GLUT_BITMAP_9_BY_15, string[i]);
}

void OutputString3(GLfloat x, GLfloat y, GLfloat z, char *string)
{
  int length;

  length = strlen(string);
  glRasterPos3f(x, y, z);
  for ( int i = 0; i < length; i++)
    glutBitmapCharacter(GLUT_BITMAP_9_BY_15, string[i]);
}


void Timestep::draw_polygons(GLdouble ScreenWidth, GLdouble ScreenHeight,unsigned int colormap,GLdouble aspect_ratio,GLdouble scale,GLdouble xtri,GLdouble ytri,GLdouble ztri,GLdouble xtran,GLdouble ytran,GLdouble ztran,unsigned int fillpoly, bool showlabel, bool showlength, bool wire, bool pointmode, int chunk)
{


  glNewList(step, GL_COMPILE);
  glPushMatrix();
  /*  bool Texture_On=0;
  bool Alpha_Add=0;
  bool Blend_On=0;
  bool Filtering_On=0;
  */
  glEnable(GL_DEPTH_TEST);
  glEnable(GL_LINE_SMOOTH);
  glShadeModel(GL_SMOOTH);
  /*
  if (Texture_On)
    glEnable(GL_TEXTURE_2D);
  else
    glDisable(GL_TEXTURE_2D);


  if (Alpha_Add)
    glBlendFunc(GL_SRC_ALPHA,GL_ONE); 
  else
    glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);

  // If we're blending, we don't want z-buffering.
  if (Blend_On)
    glDisable(GL_DEPTH_TEST); 
  else
    glEnable(GL_DEPTH_TEST); 

  if (Filtering_On) {
    glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER,
                    GL_LINEAR_MIPMAP_LINEAR);
    glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_LINEAR);
  } else {
    glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER,
                    GL_NEAREST_MIPMAP_NEAREST);
    glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_NEAREST);
  }
  */
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glClearColor(0.0, 0.0, 0.0, 0.0);
  glLineWidth(2.0);
  gluPerspective(100.0*scale, aspect_ratio, 1.0, 100.0);
  /*  glFrustum(-radius*scale, radius*scale,
    -radius*scale, radius*scale,
    -radius*2.0, radius*scale*10.0);*/
  glMatrixMode(GL_MODELVIEW);
  gluLookAt (0.0, 0.0, radius*2.0,
             0.0,0.0,0.0,
             0.0, 1.0, 0.0);
  glMatrixMode(GL_MODELVIEW);
  glPolygonMode(GL_FRONT_AND_BACK, fillpoly);
  glRotated(xtri, 1.0, 0.0, 0.0);
  glRotated(ytri, 0.0, 1.0, 0.0);
  glRotated(ztri, 0.0, 0.0, 1.0);
  glTranslated(xtran, ytran,ztran);
  glPointSize(2.0);
  display(colormap,wire,pointmode,fillpoly,chunk);
  if (showlabel)
    showlabels();
  if (showlength)
    showlengths();
  /*  glLoadIdentity();
      glMatrixMode(GL_PROJECTION);
      glPushMatrix();
      glLoadIdentity();
      glOrtho(0,ScreenWidth,0,ScreenHeight,-1.0,1.0);
      glDisable(GL_TEXTURE_2D);
      glDisable(GL_LIGHTING);
      // We don't want depth-testing either.
      glDisable(GL_DEPTH_TEST); 
      glColor4f(1.0,1.0,1.0,1.0);
      char statstring[80];
      snprintf(statstring,80,"Timestep: %u",step);
      OutputString(1,1,statstring);
      glPopMatrix();
  */
  glPopMatrix();
  glEndList();
  glCallList(step);
}


bool Timestep::loadfile(const char *filename, GLfloat * tetcolor, unsigned int timestep,unsigned int chunk)
{
  FILE *infile;
  vector < Coordinate > coordinates;
  vector < unsigned int >vertices;
  char foo[20];
  unsigned int onevertex;
  infile = fopen(filename, "r");
  if (infile == NULL) {
    fprintf(stderr, "cannot open %s \n", filename);
    return (0);
  }
  GLdouble x, y, z;
  GLfloat r, g, b;
  step=timestep;
  while (fscanf(infile, "%lf%lf%lf%f%f%f", &x, &y, &z, &r, &g, &b)>0
         && !feof(infile))
    coordinates.push_back(Coordinate(chunk, x, y, z, r, g, b));
  // add any new coordinates to points
  // this could be more efficient, but we only do it once per file
  // so optimization is not a priority.
  for (unsigned int i = 0; i < coordinates.size(); i++) {
    vector < Coordinate >::iterator coorditer =
      find(points.begin(), points.end(), coordinates[i]);
    if (coorditer == points.end())
      points.push_back(coordinates[i]);
  }
  fscanf(infile, "%s\n", foo);
  if (strlen(foo)) {
    while (!feof(infile) && fscanf(infile, "%u", &onevertex)>0) {
      vertices.push_back(onevertex);
    }
    if(vertices.size()>0)
      {
        unsigned int chunkoffset=vertexarray.size();
        // now that we have the points and the faces add them to the drawarrays
        for (unsigned int vertex = 0; vertex < vertices.size(); vertex++) {
          if (vertices[vertex] > coordinates.size()) {
            fprintf(stderr,
                    "file format error vertex %d at %d line %d row %d is an undefined point beyond %d\n",
                    vertices[vertex], vertex, vertex / 12, vertex % 12,
                    coordinates.size());
            return (false);
          }
          Coordinate thiscoord = coordinates[vertices[vertex]];
          // put every face in the drawarray
          vertexarray.push_back(thiscoord.get_x());
          vertexarray.push_back(thiscoord.get_y());
          vertexarray.push_back(thiscoord.get_z());
          // put every point color in the pointcolor array
          pointcolor.push_back(thiscoord.get_r());
          pointcolor.push_back(thiscoord.get_g());
          pointcolor.push_back(thiscoord.get_b());
          // set the chunk color to the tetcolor parameter
          chunkcolor.push_back(tetcolor[0]);
          chunkcolor.push_back(tetcolor[1]);
          chunkcolor.push_back(tetcolor[2]);
          // put the line color in for the point
          linecolor.push_back(0.8);
          linecolor.push_back(0.8);
          linecolor.push_back(0.8);
        }
        unsigned int chunkend=vertexarray.size();
        pair <unsigned int,unsigned int> chunkoffpair;
        chunkoffpair.first= chunkoffset;
        chunkoffpair.second= chunkend-chunkoffset;
        chunkmap[chunk]=chunkoffpair;
        // each face has 3 edges which need midpoints and lengths
        for (unsigned int vertex = 0; vertex < vertices.size(); vertex += 3) {
          Coordinate c1 = coordinates[vertices[vertex]];
          Coordinate c2 = coordinates[vertices[vertex + 1]];
          Coordinate c3 = coordinates[vertices[vertex + 2]];
          Midpoint mid1_2(c1.get_x(), c2.get_x(), c1.get_y(), c2.get_y(),
                          c1.get_z(), c2.get_z());
          Midpoint mid2_3(c2.get_x(), c3.get_x(), c2.get_y(), c3.get_y(),
                          c2.get_z(), c3.get_z());
          Midpoint mid3_1(c3.get_x(), c1.get_x(), c3.get_y(), c1.get_y(),
                          c3.get_z(), c1.get_z());
          // now add any newly found midpoints
          vector < Midpoint >::iterator miter =
            find(midpoints.begin(), midpoints.end(), mid1_2);
          if (miter == midpoints.end())
            midpoints.push_back(mid1_2);
          miter = find(midpoints.begin(), midpoints.end(), mid2_3);
          if (miter == midpoints.end())
            midpoints.push_back(mid2_3);
          miter = find(midpoints.begin(), midpoints.end(), mid3_1);
          if (miter == midpoints.end())
            midpoints.push_back(mid3_1);
        }
      }
  }
  else {
    return (false);
  }
  return true;
}


void Timestep::calc_and_move()
{
  double xmin,ymin,zmin;
  double xmax,ymax,zmax;
  double xwidth,ywidth,zwidth;
  double xcenter,ycenter,zcenter;
  xmin=ymin=zmin=10000000.0;
  xmax=ymax=zmax=-10000000.0;
  for(unsigned int i=0;i<points.size();i++)
    {
      if(xmin>points[i].get_x())
        xmin= points[i].get_x();
      if(ymin>points[i].get_y())
        ymin= points[i].get_y();
      if(zmin>points[i].get_z())
        zmin= points[i].get_z();
      if(xmax<points[i].get_x())
        xmax= points[i].get_x();
      if(ymax<points[i].get_y())
        ymax= points[i].get_y();
      if(zmax<points[i].get_z())
        zmax= points[i].get_z();
    }
  xcenter=(xmax+xmin)/2;
  ycenter=(ymax+ymin)/2;
  zcenter=(zmax+zmin)/2;
  xwidth=xmax-xmin;
  ywidth=ymax-ymin;
  zwidth=zmax-zmin;
  radius=max(zwidth,max(xwidth,ywidth));
  // now move them
  for(unsigned int i=0;i<vertexarray.size()-2;i+=3)
    {
      vertexarray[i]-=xmin+xcenter;
      vertexarray[i+1]-=ymin+ycenter;
      vertexarray[i+2]-=zmin+zcenter;
    }
  for(unsigned int i=0;i<points.size();i++)
    {
      Coordinate m=points[i];
      // leave the original label
      points[i].unsafe_set(m.get_x()-xmin-xcenter,m.get_y()-ymin-ycenter,m.get_z()-zmin-zcenter);
    }
  for(unsigned int i=0;i<midpoints.size();i++)
    {
      Midpoint m=midpoints[i];
      midpoints[i].set_mx(m.get_mx()-xmin-xcenter);
      midpoints[i].set_my(m.get_my()-ymin-ycenter);
      midpoints[i].set_mz(m.get_mz()-zmin-zcenter);
    }
}



void Timestep::display(unsigned int cmap,bool wire,bool pointmode, unsigned int fillpoly, int chunk)
{
// draw it
// label a vertex
// label an edge at the midpoint
  // draw a colored polygon
  // draw a wire polygon over it
  glEnableClientState(GL_VERTEX_ARRAY);
  glEnableClientState(GL_COLOR_ARRAY);
  // cache the current mode
  if(pointmode)
    glPolygonMode(GL_FRONT_AND_BACK, GL_POINT);
  else 
    glPolygonMode(GL_FRONT_AND_BACK, fillpoly);
  switch (cmap) {
  case 0:
    glColorPointer(3, GL_FLOAT, 0, &randomcolor[0]);
    break;
  case 1:
    glColorPointer(3, GL_FLOAT, 0, &chunkcolor[0]);
    break;
  case 2:
  default:
    glColorPointer(3, GL_FLOAT, 0, &pointcolor[0]);
    break;
  }
  if(chunk>=0)
    {
      // find the offsets for the chunk
      map <unsigned int, chunk_delim>::iterator chunkoff=chunkmap.find((unsigned int) chunk);
      if(chunkoff!=chunkmap.end())
        {
          glVertexPointer(3, GL_DOUBLE, 0, &vertexarray[chunkoff->second.first]);
          switch (cmap) {
          case 0:
            glColorPointer(3, GL_FLOAT, 0, &randomcolor[chunkoff->second.first]);
            break;
          case 1:
            glColorPointer(3, GL_FLOAT, 0, &chunkcolor[chunkoff->second.first]);
            break;
          case 2:
          default:
            glColorPointer(3, GL_FLOAT, 0, &pointcolor[chunkoff->second.first]);
            break;
          }
          glDrawArrays(GL_TRIANGLES, 0, chunkoff->second.second / 3);
        }
      else
        {
          glVertexPointer(3, GL_DOUBLE, 0, &vertexarray[0]);
          glDrawArrays(GL_TRIANGLES, 0, vertexarray.size() / 3);
        }
    }
  else
    {
      glVertexPointer(3, GL_DOUBLE, 0, &vertexarray[0]);
      glDrawArrays(GL_TRIANGLES, 0, vertexarray.size() / 3);
    }
  if(wire && !pointmode)
    {
      glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
      glVertexPointer(3, GL_DOUBLE, 0, &vertexarray[0]);
      glColorPointer(3, GL_FLOAT, 0, &linecolor[0]);
      glDrawArrays(GL_TRIANGLES, 0, vertexarray.size() / 3);
    }
}



void Timestep::showlengths()
{
  glColor3f(1.0, 1.0, 1.0);
  for (unsigned int i = 0; i < midpoints.size(); i++) {
    OutputString3(midpoints[i].get_mx(), midpoints[i].get_my(),
                  midpoints[i].get_mz(), midpoints[i].get_length());
  }
}


void Timestep::showlabels()
{
  glColor3f(1.0, 1.0, 1.0);
  for (unsigned int i = 0; i < points.size(); i++) {
    OutputString3(points[i].get_x(), points[i].get_y(),
                  points[i].get_z(), points[i].get_label());

  }
}


void Timestep::set_random_color(vector <Coordinate> randompoints)
{
  //
  for(unsigned int i=0;i<vertexarray.size()-2;i+=3)
    {
      Coordinate finder(0,vertexarray[i],vertexarray[i+1],vertexarray[i+2],0.,0.,0.);
      //theoretically equal_range is more efficient.
      // but in practice it fails to find existing points
      // so its faster, but error prone.  which is unacceptable
      // so when it fails, we use find
      CoordIterPair p= equal_range(randompoints.begin(),randompoints.end(),finder);
      if(p.first !=p.second)
        {// we got it
          randomcolor.push_back(p.first->get_r());
          randomcolor.push_back(p.first->get_g());
          randomcolor.push_back(p.first->get_b());

        }
      else
        {
          vector <Coordinate>::iterator usefind= find(randompoints.begin(),randompoints.end(),finder);
          if(usefind !=randompoints.end())
            {
              randomcolor.push_back(usefind->get_r());
              randomcolor.push_back(usefind->get_g());
              randomcolor.push_back(usefind->get_b());
            }
          else
            {
              fprintf(stderr,"ERROR point not found \n");
              finder.edump();
              fprintf(stderr,"first is \n");
              if(p.first!=randompoints.end())
                p.first->edump();
              fprintf(stderr,"second is \n");
              if(p.second!=randompoints.end())
                p.second->edump();
              fprintf(stderr,"in color vector \n");
              for_each(randompoints.begin(),randompoints.end(),mem_fun_ref(&Coordinate::edump));
            }
        }
    }
}

---