edge = edges[n];    v1 = case2[edge2vertex[edge][0]];    v2 = case2[edge2vertex[edge][1]];    val1 = vals[v1]-level; val2 = vals[v2]-level;    denom = val2 - val1;    factor = 0.5;    if(denom!=0.0)factor = -val1/denom;    if(factor<0.5){closestnodes[n]=nodeindexes[v1];}    else{closestnodes[n]=nodeindexes[v2];}    if(factor>1.0){/*factor=1;*/outofbounds=1;}    if(factor<0.0){/*factor=0.0;*/outofbounds=1;}    xx = xxval[v1]*(1.0-factor) + xxval[v2]*factor;    yy = yyval[v1]*(1.0-factor) + yyval[v2]*factor;    zz = zzval[v1]*(1.0-factor) + zzval[v2]*factor;    xvert[n] = xx;    yvert[n] = yy;    zvert[n] = zz;    if(tvert!=NULL)tvert[n]= tvals[v1]*(1.0-factor) + tvals[v2]*factor;  }  if(outofbounds==1){    printf("*** warning - computed isosurface vertices are out of bounds for :\n");    printf("case number=%i level=%f\n",casenum,level);    printf("values=");    for(n=0;n<8;n++){printf("%f ",vals[n]);}    printf("\n");    printf("x=%f %f y=%f %f z=%f %f\n\n",x[0],x[1],y[0],y[1],z[1],z[2]);  }  /* copy coordinates to output array */  *nvert = nedges;  *ntriangles = npath;  for(n=0;n<npath;n++){triangles[n] = path[n];}}/* ------------------ calcNormal2 ------------------------ */void calcNormal2(unsigned short *v1, unsigned short *v2, unsigned short *v3, float *out, float *area){  float u[3], v[3];  int i;  for(i=0;i<3;i++){    u[i]=v2[i]-v1[i];    v[i]=v3[i]-v1[i];  }  out[0] = u[1]*v[2] - u[2]*v[1];  out[1] = u[2]*v[0] - u[0]*v[2];  out[2] = u[0]*v[1] - u[1]*v[0];  *area = sqrt(out[0]*out[0]+out[1]*out[1]+out[2]*out[2])/2.0;  ReduceToUnit(out);}/* ------------------ calcNormal ------------------------ */void calcNormal(float *xx, float *yy, float *zz, float *out){  float u[3],v[3], p1[3], p2[3], p3[3];  static const int x = 0;  static const int y = 1;  static const int z = 2;  p1[x]=xx[0]; p1[y]=yy[0]; p1[z]=zz[0];  p2[x]=xx[1]; p2[y]=yy[1]; p2[z]=zz[1];  p3[x]=xx[2]; p3[y]=yy[2]; p3[z]=zz[2];  u[x] = p2[x] - p1[x];  u[y] = p2[y] - p1[y];  u[z] = p2[z] - p1[z];  v[x] = p3[x] - p1[x];  v[y] = p3[y] - p1[y];  v[z] = p3[z] - p1[z];  out[x] = u[y]*v[z] - u[z]*v[y];  out[y] = u[z]*v[x] - u[x]*v[z];  out[z] = u[x]*v[y] - u[y]*v[x];  ReduceToUnit(out);}/* ------------------ ReduceToUnit ------------------------ */void ReduceToUnit(float *v){  float length;  length = (float)sqrt((v[0]*v[0]+v[1]*v[1]+v[2]*v[2]));  if(length==0.0f)length=1.0f;  v[0] /= length;  v[1] /= length;  v[2] /= length;}/* ------------------ GetIsoSurface ------------------------ */int GetIsosurface(isosurface *surface, float *data, float *tdata, int *iblank_cell, float level,                   float *xplt, int nx, float *yplt, int ny, float *zplt, int nz,                   int isooffset                   ){#define ijk(i,j,k) ((i)+(j)*nx+(k)*nx*ny)#define ijkcell(i,j,k) ((i)+(j)*ibar+(k)*ijbar)#define ij(i,j) ((i)+(j)*nx)  int ibar,ijbar;  float xvert[12], yvert[12], zvert[12], tvert[12], *tvertptr=NULL;  int triangles[18];  int nvert, ntriangles;  int i, j, k;  float xxx[2], yyy[2], zzz[2], vals[8], tvals[8], *tvalsptr=NULL;  int nodeindexes[8], closestnodes[18];  float *xx, *yy, *zz;  int ijkbase,ip1jk,ijkp1,ip1jkp1,ijp1k,ip1jp1k,ijp1kp1,ip1jp1kp1;  int ijbase;  int nxy;    ibar = nx-1;  ijbar = (nx-1)*(ny-1);  xx = xxx;  yy = yyy;  zz = zzz;  nxy = nx*ny;  if(surface->defined==0)return 0;  if(tdata!=NULL){    tvalsptr=tvals;    tvertptr=tvert;  }  for(i=0;i<nx-isooffset;){    xx[0]=xplt[i];    xx[1]=xplt[i+isooffset];    for(j=0;j<ny-isooffset;){      yy[0]=yplt[j];      yy[1]=yplt[j+isooffset];      ijbase = ij(i,j);      for(k=0;k<nz-isooffset;){        ijkbase = ijbase + k*nxy;        ip1jk = ijkbase + isooffset;        ijkp1 = ijkbase + isooffset*nxy;        ip1jkp1 = ijkbase + isooffset*(1+nxy);        ijp1k = ijkbase + isooffset*nx;        ip1jp1k = ijkbase + isooffset*(1+nx);        ijp1kp1 = ijkbase + isooffset*(nx+nxy);        ip1jp1kp1 = ijkbase + isooffset*(1+nx+nxy);        if(isooffset!=1||iblank_cell==NULL||iblank_cell[ijkcell(i,j,k)]!=0){          zz[0]=zplt[k];          zz[1]=zplt[k+isooffset];          nodeindexes[0]=ijkbase;          nodeindexes[1]=ijp1k;          nodeindexes[2]=ip1jp1k;          nodeindexes[3]=ip1jk;          nodeindexes[4]=ijkp1;          nodeindexes[5]=ijp1kp1;          nodeindexes[6]=ip1jp1kp1;          nodeindexes[7]=ip1jkp1;          vals[0]=data[ijkbase];          vals[1]=data[ijp1k];          vals[2]=data[ip1jp1k];          vals[3]=data[ip1jk];          vals[4]=data[ijkp1];          vals[5]=data[ijp1kp1];          vals[6]=data[ip1jp1kp1];          vals[7]=data[ip1jkp1];          if(tdata!=NULL){            tvals[0]=tdata[ijkbase];            tvals[1]=tdata[ijp1k];            tvals[2]=tdata[ip1jp1k];            tvals[3]=tdata[ip1jk];            tvals[4]=tdata[ijkp1];            tvals[5]=tdata[ijp1kp1];            tvals[6]=tdata[ip1jp1kp1];            tvals[7]=tdata[ip1jkp1];          }          GetIsobox(xx, yy, zz, vals, tvalsptr, nodeindexes, level,                    xvert, yvert, zvert, tvertptr, closestnodes, &nvert, triangles, &ntriangles);          if(UpdateIsosurface(surface, xvert, yvert, zvert, tvertptr,                               closestnodes, nvert, triangles, ntriangles)!=0)return 1;        }        k+=isooffset;      }      j+=isooffset;    }    i+=isooffset;  }  surface->defined=1;  return 0;}/* ------------------ compareisonodes ------------------------ */int compareisonodes( const void *arg1, const void *arg2 ){  int i, j;  i=*(int *)arg1;  j=*(int *)arg2;  i *= 3;  j *= 3;  if(vertices[i]<vertices[j])return -1;  if(vertices[i]>vertices[j])return 1;  i++; j++;  if(vertices[i]<vertices[j])return -1;  if(vertices[i]>vertices[j])return 1;  i++; j++;  if(vertices[i]<vertices[j])return -1;  if(vertices[i]>vertices[j])return 1;  return 0;}/* ------------------ computerank ------------------------ */int computerank( const void *arg1, const void *arg2 ){  int i, j;  i=*(int *)arg1;  j=*(int *)arg2;  if(sortedlist[i]<sortedlist[j])return -1;  if(sortedlist[i]>sortedlist[j])return 1;  return 0;}int order_closestnodes( const void *arg1, const void *arg2 ){  int i, j;  int ii,jj;  i=*(int *)arg1;  j=*(int *)arg2;  ii = closestnodes[i];  jj = closestnodes[j];  if(ii<jj)return -1;  if(ii>jj)return 1;  return 0;}/* ------------------ CompressIsosurface ------------------------ */int CompressIsosurface(isosurface *surface, int reduce_triangles,                        float xmin, float xmax,                         float ymin, float ymax,                         float zmin, float zmax){  int i,j,nn;  float *x=NULL, *y=NULL, *z=NULL, *t=NULL;  int *map=NULL,*map2=NULL,nmap2,*vertexmap=NULL,*inverse_vertexmap=NULL;  int nvertices;  unsigned short *newvertices=NULL;  int *triangles=NULL,*newtriangles=NULL,*newtriangles2=NULL,ntriangles;  float xyzmaxdiff;  int *cs=NULL, *ordered_closestnodes=NULL;  int v1, v2, v3;  int ii,iim1;  unsigned short vx, vy, vz;  int sumx, sumy, sumz, sumt;  int nnewvertices;  float tmin, tmax, tmaxmin;  int flag;  unsigned short *tvertices,*newtvertices;  nvertices=surface->nvertices;  if(nvertices==0)return 0;  vertices=surface->vertices;  x=surface->xvert;   y=surface->yvert;   z=surface->zvert;   if(surface->tvert!=NULL)t=surface->tvert;  sortedlist=surface->sortedlist;  rank=surface->rank;  triangles=surface->triangles;  ntriangles=surface->ntriangles;    if(surface->dataflag==1){    tvertices=surface->tvertices;    FREEMEMORY(tvertices);    if(NewMemory((void **)&tvertices,nvertices*sizeof(unsigned short))==0){      FREEMEMORY(tvertices);      return 1;    }    surface->tvertices=tvertices;  }  FREEMEMORY(vertices); FREEMEMORY(rank);  if(NewMemory((void **)&vertices,3*nvertices*sizeof(unsigned short))==0||     NewMemory((void **)&rank,nvertices*sizeof(int))==0||     NewMemory((void **)&map,nvertices*sizeof(int))==0||     NewMemory((void **)&map2,nvertices*sizeof(int))==0||     NewMemory((void **)&sortedlist,nvertices*sizeof(int))==0){    FREEMEMORY(vertices);    FREEMEMORY(rank);    FREEMEMORY(map);    FREEMEMORY(map2);    FREEMEMORY(sortedlist);    return 1;  }  surface->vertices=vertices;  surface->rank=rank;  xyzmaxdiff = xmax - xmin;  if(ymax-ymin>xyzmaxdiff)xyzmaxdiff=ymax-ymin;  if(zmax-zmin>xyzmaxdiff)xyzmaxdiff=zmax-zmin;  surface->xmin=xmin;  surface->ymin=ymin;  surface->zmin=zmin;  surface->xyzmaxdiff=xyzmaxdiff;  if(surface->dataflag==1){    tmin=t[0];    tmax=tmin;    for(i=1;i<nvertices;i++){      if(t[i]<tmin)tmin=t[i];      if(t[i]>tmax)tmax=t[i];    }    surface->tmin=tmin;    surface->tmax=tmax;    tmaxmin=tmax-tmin;    if(tmaxmin>0.0){      for(i=0;i<nvertices;i++){        tvertices[i]=65535*(t[i]-tmin)/tmaxmin;      }    }    else{      for(i=0;i<nvertices;i++){        tvertices[i]=0;      }    }  }  for(i=0;i<nvertices;i++){    vertices[3*i+0]=65535*(x[i]-xmin)/xyzmaxdiff;    vertices[3*i+1]=65535*(y[i]-ymin)/xyzmaxdiff;    vertices[3*i+2]=65535*(z[i]-zmin)/xyzmaxdiff;    sortedlist[i]=i;    map[i]=i;  	rank[i]=i;  }  qsort((int *)sortedlist,(size_t)nvertices,sizeof(int),compareisonodes);  qsort((int *)rank,(size_t)nvertices,sizeof(int),computerank);  j=0;  map[0]=0;  map2[0]=0;  nmap2=1;  for(i=1;i<nvertices;i++){    if(compareisonodes(sortedlist+i-1,sortedlist+i)!=0){  	  j++;  	  map2[j]=i;  	  nmap2++;    }    map[i]=j;  }  nvertices=nmap2;  if(NewMemory((void **)&newvertices,3*nvertices*sizeof(unsigned short))==0||     NewMemory((void **)&cs,nvertices*sizeof(int))==0){    FREEMEMORY(newvertices);    FREEMEMORY(cs);    return 1;  }  if(surface->dataflag==1){    if(NewMemory((void **)&newtvertices,nvertices*sizeof(unsigned short))==0){      FREEMEMORY(newtvertices);      return 1;    }  }  closestnodes=surface->closestnodes;  for(i=0;i<nvertices;i++){	  j=sortedlist[map2[i]];	  newvertices[3*i+0]=vertices[3*j+0];	  newvertices[3*i+1]=vertices[3*j+1];	  newvertices[3*i+2]=vertices[3*j+2];    cs[i] = closestnodes[j];  }  if(surface->dataflag==1){    for(i=0;i<nvertices;i++){	    j=sortedlist[map2[i]];	    newtvertices[i]=tvertices[j];    }  }  FREEMEMORY(closestnodes);  surface->closestnodes = cs;  if(NewMemory((void **)&newtriangles,ntriangles*sizeof(int))==0){    return 1;  }  for(i=0;i<ntriangles;i++){newtriangles[i]=map[rank[triangles[i]]];}  surface->triangles=newtriangles;  surface->vertices=newvertices;  surface->nvertices=nvertices;  FREEMEMORY(vertices);  if(surface->dataflag==1){    surface->tvertices=newtvertices;    FREEMEMORY(tvertices);  }  FREEMEMORY(triangles);  FREEMEMORY(map);FREEMEMORY(map2);FREEMEMORY(sortedlist);  if(reduce_triangles!=1)return 0;  /* phase II compression, reduce the number of triangles */  /* first, eliminate triangles whose nodes (2 or more) are closest to the same grid point */  if(NewMemory((void **)&newtriangles2,ntriangles*sizeof(int))==0){    return 1;  }  nn=0;  for(i=0;i<ntriangles/3;i++){    v1=newtriangles[3*i];    v2=newtriangles[3*i+1];    v3=newtriangles[3*i+2];    if(cs[v1]!=cs[v2]&&cs[v1]!=cs[v3]&&cs[v2]!=cs[v3]){      newtriangles2[nn++]=v1;      newtriangles2[nn++]=v2;      newtriangles2[nn++]=v3;    }  }  FREEMEMORY(newtriangles);  surface->triangles=newtriangles2;  ntriangles=nn;  surface->ntriangles=ntriangles;  /* sort the closestnodes list */  closestnodes=surface->closestnodes;  if(NewMemory((void **)&ordered_closestnodes,nvertices*sizeof(int))==0){    return 1;  }  for(i=0;i<nvertices;i++){ordered_closestnodes[i]=i;}  qsort((int *)ordered_closestnodes,(size_t)nvertices,sizeof(int),order_closestnodes);  if(NewMemory((void **)&vertexmap,nvertices*sizeof(int))==0||     NewMemory((void **)&inverse_vertexmap,nvertices*sizeof(int))==0){    FREEMEMORY(vertexmap);    FREEMEMORY(inverse_vertexmap);    return 1;  }  for(i=0;i<nvertices;i++){vertexmap[i]=i;}  nn=0;  sumx=0; sumy = 0; sumz = 0; sumt=0;  /* average nodes */  nn=0;  vertices = surface->vertices;  if(surface->dataflag==1)tvertices = surface->tvertices;