nnewvertices=0;  for(i=1;i<nvertices+1;i++){    iim1=ordered_closestnodes[i-1];    if(i!=nvertices)ii=ordered_closestnodes[i];    inverse_vertexmap[iim1] = nnewvertices;    nn++;    vx = vertices[3*iim1  ]; sumx += vx;    vy = vertices[3*iim1+1]; sumy += vy;    vz = vertices[3*iim1+2]; sumz += vz;    if(surface->dataflag==1)sumt += tvertices[iim1];    flag=0;    if(i!=nvertices&&closestnodes[ii]!=closestnodes[iim1])flag=1;    if(i==nvertices||flag==1){      vertexmap[nnewvertices++] = iim1;      vertices[3*iim1  ]=sumx/nn;      vertices[3*iim1+1]=sumy/nn;      vertices[3*iim1+2]=sumz/nn;      if(surface->dataflag==1)tvertices[iim1]=sumt/nn;      nn=0;       sumx = 0; sumy = 0;sumz = 0; sumt = 0;    }  }  if(NewMemory((void **)&newvertices,3*nnewvertices*sizeof(unsigned short))==0){    return 1;  }  if(surface->dataflag==1){    if(NewMemory((void **)&newtvertices,nnewvertices*sizeof(unsigned short))==0){      return 1;    }  }  for(i=0;i<nnewvertices;i++){    ii = vertexmap[i];    newvertices[3*i  ] = vertices[3*ii  ];    newvertices[3*i+1] = vertices[3*ii+1];    newvertices[3*i+2] = vertices[3*ii+2];    if(surface->dataflag==1)newtvertices[i]=tvertices[ii];  }  FREEMEMORY(vertices);  surface->vertices = newvertices;  surface->nvertices = nnewvertices;  if(surface->dataflag==1){    FREEMEMORY(tvertices);    surface->tvertices=newtvertices;  }  triangles = surface->triangles;  for(i=0;i<ntriangles;i++){    triangles[i] = inverse_vertexmap[triangles[i]];  }  FREEMEMORY(ordered_closestnodes);  FREEMEMORY(vertexmap);  FREEMEMORY(inverse_vertexmap);  return 0;}/* ------------------ UpdateIsosurface ------------------------ */int UpdateIsosurface(isosurface *surface,                       float *xvert, float *yvert, float *zvert, float *tvert,                       int *closestnodes, int nvert, int *triangles, int ntriangles){  int n,ns, noldvert, *is;  float *xs=NULL, *ys=NULL, *zs=NULL, *ts=NULL;  int *cn=NULL;  if(ResizeSurface(surface,nvert,ntriangles,0)!=0)return 1;  /* copy vertex data */  if(nvert>0){    noldvert = surface->nvertices;    xs = surface->xvert + noldvert;    ys = surface->yvert + noldvert;    zs = surface->zvert + noldvert;    if(tvert!=NULL)ts = surface->tvert + noldvert;    cn = surface->closestnodes + noldvert;    for(n=0;n<nvert;n++){      xs[n] = xvert[n];      ys[n] = yvert[n];      zs[n] = zvert[n];      if(tvert!=NULL)ts[n]=tvert[n];      cn[n] = closestnodes[n];    }    surface->nvertices = noldvert + nvert;  }  /* copy triangle path data */  if(ntriangles>0){    ns = surface->ntriangles;    is = surface->triangles + ns;    for(n=0;n<ntriangles;n++){is[n] = triangles[n]+noldvert;}    surface->ntriangles = ns + ntriangles;  }  return 0;}/* ------------------ ResizeSurface ------------------------ */int ResizeSurface(isosurface *surfacedata, int incvert, int inctriangles, int incnorm){  int maxnum, *itemp=NULL;  float *temp=NULL;  /* resize vertex data if necessary */  maxnum = surfacedata->nvertices+incvert;  if(maxnum>surfacedata->maxvertices){    maxnum += INCPOINTS;    surfacedata->maxvertices = maxnum;    temp = surfacedata->xvert;    if(temp==NULL){      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    else{      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    surfacedata->xvert = temp;    temp = surfacedata->yvert;    if(temp==NULL){      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    else{      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    surfacedata->yvert = temp;    temp = surfacedata->zvert;    if(temp==NULL){      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    else{      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    surfacedata->zvert = temp;    if(surfacedata->dataflag==1){      temp = surfacedata->tvert;      if(temp==NULL){        if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;      }      else{        if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;      }      surfacedata->tvert = temp;    }    itemp = surfacedata->closestnodes;    if(itemp==NULL){      if(NewMemory((void **)&itemp,maxnum*sizeof(int))==0)return 1;    }    else{      if(ResizeMemory((void **)&itemp,maxnum*sizeof(int))==0)return 1;    }    surfacedata->closestnodes = itemp;  }  /* resize norm data if necessary */  maxnum = surfacedata->nnorm+incnorm;  if(maxnum>surfacedata->maxnorm){    maxnum += INCPOINTS;    surfacedata->maxnorm = maxnum;    temp = surfacedata->xnorm;    if(temp==NULL){      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    else{      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    surfacedata->xnorm = temp;    temp = surfacedata->ynorm;    if(temp==NULL){      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    else{      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    surfacedata->ynorm = temp;    temp = surfacedata->znorm;    if(temp==NULL){      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    else{      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;    }    surfacedata->znorm = temp;  }  /* resize triangles data if necessary */  maxnum = surfacedata->ntriangles+inctriangles;  if(maxnum>surfacedata->maxtriangles){    maxnum += INCPOINTS;    surfacedata->maxtriangles = maxnum;    itemp = surfacedata->triangles;    if(itemp==NULL){      if(NewMemory((void **)&itemp,maxnum*sizeof(int))==0)return 1;    }    else{      if(ResizeMemory((void **)&itemp,maxnum*sizeof(int))==0)return 1;    }    surfacedata->triangles = itemp;  }  return 0;}/* ------------------ InitIsosurface ------------------------ */void InitIsosurface(isosurface *surfacedata, float level, float *color,int colorindex){  surfacedata->level = level;  surfacedata->color = color;  surfacedata->colorindex=colorindex;  surfacedata->defined=1;  surfacedata->nvertices=0;  surfacedata->maxvertices=-1;  surfacedata->ntriangles=0;  surfacedata->maxtriangles=-1;  surfacedata->nnorm=0;  surfacedata->maxnorm=-1;  surfacedata->normtype=0;  surfacedata->xvert=NULL;  surfacedata->yvert=NULL;  surfacedata->zvert=NULL;  surfacedata->xnorm=NULL;  surfacedata->ynorm=NULL;  surfacedata->znorm=NULL;  surfacedata->vertices=NULL;  surfacedata->tvertices=NULL;  surfacedata->rank=NULL;  surfacedata->sortedlist=NULL;  surfacedata->triangles=NULL;  surfacedata->plottype = 0;  surfacedata->closestnodes=NULL;  surfacedata->tvert=NULL;}/* ------------------ freesurface ------------------------ */void freesurface(isosurface *surfacedata){  if(surfacedata->defined==0)return;  FREEMEMORY(surfacedata->xvert);  FREEMEMORY(surfacedata->yvert);  FREEMEMORY(surfacedata->zvert);  FREEMEMORY(surfacedata->xnorm);  FREEMEMORY(surfacedata->ynorm);  FREEMEMORY(surfacedata->znorm);  FREEMEMORY(surfacedata->triangles);  FREEMEMORY(surfacedata->vertices);  FREEMEMORY(surfacedata->sortedlist);  FREEMEMORY(surfacedata->rank);  FREEMEMORY(surfacedata->closestnodes);  FREEMEMORY(surfacedata->tvert);  FREEMEMORY(surfacedata->tvertices);  surfacedata->defined=0;}void SmoothIsoSurface(isosurface *surfacedata){  int *triangles_i;  unsigned short *vertices_i,*v1,*v2,*v3;  int n,i1,i2,i3;  int ntriangles_i,nvertices_i;  float out[3],area;  float *xyznorm;  short *norm,*vertexnorm;  triangles_i = surfacedata->triangles;  vertices_i = surfacedata->vertices;  ntriangles_i = surfacedata->ntriangles;  nvertices_i = surfacedata->nvertices;  if(nvertices_i==0||ntriangles_i==0)return;  NewMemory((void **)&xyznorm,3*nvertices_i*sizeof(float));  NewMemory((void **)&norm,ntriangles_i*sizeof(short));  NewMemory((void **)&vertexnorm,3*nvertices_i*sizeof(short));  surfacedata->norm=norm;  surfacedata->vertexnorm=vertexnorm;  for(n=0;n<3*nvertices_i;n++){xyznorm[n]=0.0;}  for(n=0;n<ntriangles_i/3;n++){    i1=3*triangles_i[3*n];    i2=3*triangles_i[3*n+1];    i3=3*triangles_i[3*n+2]; /*   for(k=0;k<3;k++){ //     vv1[k]=xyzmin[k]+xyzmaxdiff*v1[k]/65535.; //     vv2[k]=xyzmin[k]+xyzmaxdiff*v2[k]/65535.; //     vv3[k]=xyzmin[k]+xyzmaxdiff*v3[k]/65535.; //   } */    v1=vertices_i+i1;    v2=vertices_i+i2;    v3=vertices_i+i3;    calcNormal2(v1,v2,v3,out,&area);    norm[3*n  ]=(short)(out[0]*32767);    norm[3*n+1]=(short)(out[1]*32767);    norm[3*n+2]=(short)(out[2]*32767);	  xyznorm[i1  ] += out[0]*area;    xyznorm[i1+1] += out[1]*area;    xyznorm[i1+2] += out[2]*area;    xyznorm[i2  ] += out[0]*area;    xyznorm[i2+1] += out[1]*area;    xyznorm[i2+2] += out[2]*area;    xyznorm[i3  ] += out[0]*area;    xyznorm[i3+1] += out[1]*area;    xyznorm[i3+2] += out[2]*area;  }  for(n=0;n<nvertices_i;n++){    ReduceToUnit(xyznorm+3*n);    vertexnorm[3*n  ]=(short)(xyznorm[3*n  ]*32767);    vertexnorm[3*n+1]=(short)(xyznorm[3*n+1]*32767);    vertexnorm[3*n+2]=(short)(xyznorm[3*n+2]*32767);  }  FREEMEMORY(xyznorm);}/* ------------------ GetNormalSurface ------------------------ */int GetNormalSurface(isosurface *surfacedata){  float out[3];  float vertx[3], verty[3], vertz[3];  float *xnorm=NULL, *ynorm=NULL, *znorm=NULL;  float *xvert=NULL, *yvert=NULL, *zvert=NULL;  int *triangles=NULL;  int ntriangles;  int nn, n, index;    ntriangles = (surfacedata->ntriangles)/3;  if(ntriangles==0)return 0;  xvert = surfacedata->xvert;  yvert = surfacedata->yvert;  zvert = surfacedata->zvert;  triangles = surfacedata->triangles;  FREEMEMORY(surfacedata->xnorm);  FREEMEMORY(surfacedata->ynorm);  FREEMEMORY(surfacedata->znorm);  if(NewMemory((void **)&surfacedata->xnorm,ntriangles*sizeof(int))==0||     NewMemory((void **)&surfacedata->ynorm,ntriangles*sizeof(int))==0||     NewMemory((void **)&surfacedata->znorm,ntriangles*sizeof(int))==0){    freesurface(surfacedata);    return 1;  }  xnorm = surfacedata->xnorm;  ynorm = surfacedata->ynorm;  znorm = surfacedata->znorm;  nn = 0;  for(n=0;n<ntriangles;n++){    index = triangles[nn++];    vertx[0] = xvert[index]; verty[0] = yvert[index]; vertz[0] = zvert[index];    index = triangles[nn++];    vertx[1] = xvert[index]; verty[1] = yvert[index]; vertz[1] = zvert[index];    index = triangles[nn++];    vertx[2] = xvert[index]; verty[2] = yvert[index]; vertz[2] = zvert[index];    calcNormal(vertx,verty,vertz,out);    xnorm[n]=out[0];    ynorm[n]=out[1];    znorm[n]=out[2];  }  return 0;}/* ------------------ DrawIsosurface ------------------------ */#ifdef pp_DRAWISOvoid DrawIsosurface(isosurface *surfacedata){  float shiny=10.;  GLfloat *color=NULL;  GLfloat specular[]={0.5,0.5,0.5,1.0};  float *xvert, *yvert, *zvert, *xnorm, *ynorm, *znorm;  float x0, y0, z0;  float x, y, z;  int *triangles, ntriangles;  int index,n,nn,nvert;  xvert = surfacedata->xvert;  yvert = surfacedata->yvert;  zvert = surfacedata->zvert;  nvert = surfacedata->nvertices;  xnorm = surfacedata->xnorm;  ynorm = surfacedata->ynorm;  znorm = surfacedata->znorm;  triangles = surfacedata->triangles;  ntriangles = surfacedata->ntriangles;  nn = 0;  color = surfacedata->color;  if(surfacedata->plottype==0){    glPushAttrib(GL_LIGHTING_BIT);#ifdef ISO_DEBUG    glLineWidth(1.0);    glBegin(GL_LINES);    nn=0;    for(n=0;n<ntriangles/3;n++){      index = triangles[nn];      x = xvert[index]; y = yvert[index]; z = zvert[index];      glVertex3f(x,y,z);      glVertex3f(x+xnorm[n],y+ynorm[n],z+znorm[n]);      nn +=3;    }    glEnd(); #endif    nn = 0;    glEnable(GL_LIGHTING);    glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,shiny);    glBegin(GL_TRIANGLES);    glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,color);    glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,specular);    for (n = 0; n < ntriangles/3; n++) {      glNormal3f(xnorm[n],ynorm[n],znorm[n]);      index = triangles[nn++];      x = xvert[index]; y = yvert[index]; z = zvert[index];      glVertex3f(x,y,z);      index = triangles[nn++];      x = xvert[index]; y = yvert[index]; z = zvert[index];      glVertex3f(x,y,z);      index = triangles[nn++];      x = xvert[index]; y = yvert[index]; z = zvert[index];      glVertex3f(x,y,z);    }    glEnd();    glPopAttrib();  }  if(surfacedata->plottype==1){    glColor3fv(color);    glBegin(GL_LINES);    for (n = 0; n < ntriangles/3; n++) {      index = triangles[nn++];      x = xvert[index]; y = yvert[index]; z = zvert[index];      glVertex3f(x,y,z);      x0 = x; y0 = y; z0 = z;      index = triangles[nn++];      x = xvert[index]; y = yvert[index]; z = zvert[index];      glVertex3f(x,y,z);      glVertex3f(x,y,z);