text: 
  
  /* 
  A* -------------------------------------------------------------------
  B* This file contains source code for the PyMOL computer program
  C* Copyright (c) Schrodinger, LLC. 
  D* -------------------------------------------------------------------
  E* It is unlawful to modify or remove this copyright notice.
  F* -------------------------------------------------------------------
  G* Please see the accompanying LICENSE file for further information. 
  H* -------------------------------------------------------------------
  I* Additional authors of this source file include:
  -* 
  -* 
  -*
  Z* -------------------------------------------------------------------
  */
  #include"os_python.h"
  
  #include"os_predef.h"
  #include"os_std.h"
  
  #include"Isosurf.h"
  #include"Tetsurf.h"
  #include"MemoryDebug.h"
  #include"Err.h"
  #include"Crystal.h"
  #include"Vector.h"
  #include"Feedback.h"
  #include"P.h"
  
  #define Trace_OFF
  
  #ifndef true
  #define true 1
  #endif
  
  #ifndef false
  #define false 0
  #endif
  
  #define O3(field,P1,P2,P3,offs) Ffloat3(field,P1+offs[0],P2+offs[1],P3+offs[2])
  
  #define O3Ptr(field,P1,P2,P3,offs) Ffloat3p(field,P1+offs[0],P2+offs[1],P3+offs[2])
  
  #define O4(field,P1,P2,P3,P4,offs) Ffloat4(field,P1+offs[0],P2+offs[1],P3+offs[2],P4)
  
  #define O4Ptr(field,P1,P2,P3,P4,offs) Ffloat4p(field,P1+offs[0],P2+offs[1],P3+offs[2],P4)
  
  #define I3(field,P1,P2,P3) Fint3(field,P1,P2,P3)
  
  #define I3Ptr(field,P1,P2,P3) Fint3p(field,P1,P2,P3)
  
  #define I4(field,P1,P2,P3,P4) Fint4(field,P1,P2,P3,P4)
  
  #define I4Ptr(field,P1,P2,P3,P4) Fint4p(field,P1,P2,P3,P4)
  
  typedef struct {
    float Point[3];
    float Normal[3];
    int NormalFlag;
    int Link;
  } PointType;
  
  typedef struct {
    PointType *p[3];
    float n[3];
    int done;
  } TriangleType;
  
  typedef struct {
    int link;
    int tri;
  } PointLinkType;
  
  #define EdgePtPtr(field,P2,P3,P4,P5) ((PointType*)Fvoid4p(field,P2,P3,P4,P5))
  
  #define EdgePt(field,P2,P3,P4,P5) (*((PointType*)Fvoid4p(field,P2,P3,P4,P5)))
  
  struct _CTetsurf {
    PyMOLGlobals *G;
    TriangleType *Tri;
    PointLinkType *PtLink;
  
    CField *VertexCodes;
    CField *ActiveEdges;
    CField *Point;
  
    int AbsDim[3], CurDim[3], CurOff[3];
    int Max[3];
    CField *Coord, *Data, *Grad;
    float Level;
    int Edge[6020];               /* 6017 */
    int EdgeStart[256];
    int TotPrim;
  };
  
  static int TetsurfAlloc(CTetsurf * II);
  static void TetsurfPurge(CTetsurf * II);
  static int TetsurfCodeVertices(CTetsurf * II);
  
  static int TetsurfFindActiveBoxes(CTetsurf * II, int mode, int *n_strip, int n_vert,
                                    int **strip_l, float **vert,
                                    MapType * voxelmap, float *a_vert,
                                    float carvebuffer, int side);
  
  #define TetsurfSubSize		50
  
  static void copy3fn(float *v1, float *v2)
  {
    v2[0] = v1[0];
    v2[1] = v1[1];
    v2[2] = v1[2];
  }
  
  
  /*===========================================================================*/
  static int ProcessTetrahedron(int *edge, int nv, int v0, int v1, int v2, int v3,
                                int e01, int e02, int e03, int e12, int e13, int e23,
                                int reflect)
  {
    int bits = (v3 << 3) + (v2 << 2) + (v1 << 1) + v0;
    if(reflect)
      bits = 0xF - bits;
    switch (bits) {
    case 0x0:
      break;
    case 0x1:
      edge[nv++] = e01;
      edge[nv++] = e02;
      edge[nv++] = e03;
      break;
    case 0x2:
      edge[nv++] = e01;
      edge[nv++] = e13;
      edge[nv++] = e12;
      break;
    case 0x3:
      edge[nv++] = e13;
      edge[nv++] = e12;
      edge[nv++] = e02;
      edge[nv++] = e03;
      edge[nv++] = e13;
      edge[nv++] = e02;
      break;
    case 0x4:
      edge[nv++] = e12;
      edge[nv++] = e23;
      edge[nv++] = e02;
      break;
    case 0x5:
      edge[nv++] = e01;
      edge[nv++] = e12;
      edge[nv++] = e03;
      edge[nv++] = e12;
      edge[nv++] = e23;
      edge[nv++] = e03;
      break;
    case 0x6:
      edge[nv++] = e01;
      edge[nv++] = e13;
      edge[nv++] = e02;
      edge[nv++] = e13;
      edge[nv++] = e23;
      edge[nv++] = e02;
      break;
    case 0x7:
      edge[nv++] = e03;
      edge[nv++] = e13;
      edge[nv++] = e23;
      break;
    case 0x8:
      edge[nv++] = e03;
      edge[nv++] = e23;
      edge[nv++] = e13;
      break;
    case 0x9:
      edge[nv++] = e13;
      edge[nv++] = e01;
      edge[nv++] = e02;
      edge[nv++] = e02;
      edge[nv++] = e23;
      edge[nv++] = e13;
      break;
    case 0xA:
      edge[nv++] = e01;
      edge[nv++] = e03;
      edge[nv++] = e12;
      edge[nv++] = e03;
      edge[nv++] = e23;
      edge[nv++] = e12;
      break;
    case 0xB:
      edge[nv++] = e23;
      edge[nv++] = e12;
      edge[nv++] = e02;
      break;
    case 0xC:
      edge[nv++] = e13;
      edge[nv++] = e02;
      edge[nv++] = e12;
      edge[nv++] = e03;
      edge[nv++] = e02;
      edge[nv++] = e13;
      break;
    case 0xD:
      edge[nv++] = e01;
      edge[nv++] = e12;
      edge[nv++] = e13;
      break;
    case 0xE:
      edge[nv++] = e01;
      edge[nv++] = e03;
      edge[nv++] = e02;
      break;
    case 0xF:
      break;
    }
    return nv;
  }
  
  
  /*===========================================================================*/
  static CTetsurf *TetsurfNew(PyMOLGlobals * G)
  {
  
    /* there are six tetrahedrons in each cube, and there are 
       sixteen different types of tetrahedrons */
  
    /* one will encounter 2^8 different kinds of cubes. */
  
    /* bits -> tetrahedral vertices 
       zyx
       0: 000
       1: 001
       2: 010
       3: 011
       4: 100
       5: 101
       6: 110
       7: 111
     */
  
    /* edges of the tetrahedron 
  
       vertex
       3210
       0:0011
       1:0101
       2:0110
       3:1001
       4:1010
       5:1100
     */
  
    /* edges of the cube 
  
       76543210
       0:00000011 edge 000-001
       1:00000101 edge 000-010
       2:00001001 edge 000-011
       3:00010001 edge 000-100
       4:00100001 edge 000-101
       5:01000001 edge 000-110
       6:10000001 edge 000-111
       7:00001010 edge 001-011
       8:00100010 edge 001-101
       9:10000010 edge 001-111
       0xA:00001100 edge 010-011
       0xB:01000100 edge 010-110
       0xC:10000100 edge 010-111
       0xD:00110000 edge 100-101
       0xE:01010000 edge 100-110
       0xF:10010000 edge 100-111
       0x10:10001000 edge 011-111
       0x11:10100000 edge 101-111
       0x12:11000000 edge 110-111
     */
  
  #define cE_000_001 0x00
  #define cE_000_010 0x01
  #define cE_000_011 0x02
  #define cE_000_100 0x03
  #define cE_000_101 0x04
  #define cE_000_110 0x05
  #define cE_000_111 0x06
  #define cE_001_011 0x07
  #define cE_001_101 0x08
  #define cE_001_111 0x09
  #define cE_010_011 0x0A
  #define cE_010_110 0x0B
  #define cE_010_111 0x0C
  #define cE_100_101 0x0D
  #define cE_100_110 0x0E
  #define cE_100_111 0x0F
  #define cE_011_111 0x10
  #define cE_101_111 0x11
  #define cE_110_111 0x12
  
  #define cM_000_001 0x00001
  #define cM_000_010 0x00002
  #define cM_000_011 0x00004
  #define cM_000_100 0x00008
  #define cM_000_101 0x00010
  #define cM_000_110 0x00020
  #define cM_000_111 0x00040
  #define cM_001_011 0x00080
  #define cM_001_101 0x00100
  #define cM_001_111 0x00200
  #define cM_010_011 0x00400
  #define cM_010_110 0x00800
  #define cM_010_111 0x01000
  #define cM_100_101 0x02000
  #define cM_100_110 0x04000
  #define cM_100_111 0x08000
  #define cM_011_111 0x10000
  #define cM_101_111 0x20000
  #define cM_110_111 0x40000
  
    register CTetsurf *I = Calloc(CTetsurf, 1);
    int c;
    int nv = 1;
    int last_nv;
    int v000, v100, v010, v110, v001, v101, v011, v111;
  
    I->G = G;
    I->Tri = NULL;
    I->PtLink = NULL;
  
    I->VertexCodes = NULL;
    I->ActiveEdges = NULL;
    I->Point = NULL;
  
    last_nv = nv;
  
    for(c = 0; c < 256; c++) {
  
      v000 = (c & 0x01) ? 1 : 0;
      v001 = (c & 0x02) ? 1 : 0;
      v010 = (c & 0x04) ? 1 : 0;
      v011 = (c & 0x08) ? 1 : 0;
      v100 = (c & 0x10) ? 1 : 0;
      v101 = (c & 0x20) ? 1 : 0;
      v110 = (c & 0x40) ? 1 : 0;
      v111 = (c & 0x80) ? 1 : 0;
  
      /* tetrahedron 0: 000, 001, 011, 111 */
      nv = ProcessTetrahedron(I->Edge, nv, v000, v001, v011, v111,
                              cE_000_001,
                              cE_000_011,
                              cE_000_111, cE_001_011, cE_001_111, cE_011_111, 0);
      /* tetrahedron 1: 000, 001, 101, 111 */
      nv = ProcessTetrahedron(I->Edge, nv, v000, v001, v101, v111,
                              cE_000_001,
                              cE_000_101,
                              cE_000_111, cE_001_101, cE_001_111, cE_101_111, 1);
  
      /* tetrahedron 2: 000, 010, 011, 111 */
      nv = ProcessTetrahedron(I->Edge, nv, v000, v010, v011, v111,
                              cE_000_010,
                              cE_000_011,
                              cE_000_111, cE_010_011, cE_010_111, cE_011_111, 1);
      /* tetrahedron 3: 000, 010, 110, 111 */
      nv = ProcessTetrahedron(I->Edge, nv, v000, v010, v110, v111,
                              cE_000_010,
                              cE_000_110,
                              cE_000_111, cE_010_110, cE_010_111, cE_110_111, 0);
  
      /* tetrahedron 4: 000, 100, 101, 111 */
      nv = ProcessTetrahedron(I->Edge, nv, v000, v100, v101, v111,
                              cE_000_100,
                              cE_000_101,
                              cE_000_111, cE_100_101, cE_100_111, cE_101_111, 0);
  
      /* tetrahedron 5: 000, 100, 110, 111 */
      nv = ProcessTetrahedron(I->Edge, nv, v000, v100, v110, v111,
                              cE_000_100,
                              cE_000_110,
                              cE_000_111, cE_100_110, cE_100_111, cE_110_111, 1);
  
      I->Edge[nv++] = (-1);       /* sentinel */
      I->EdgeStart[c] = last_nv;
      last_nv = nv;
    }
    /*   printf("%d\n",nv);
       for(c=1;c<nv;c++) {
       if(Edge[c]<0) {
       printf("\n");
       } else {
       printf("%02X ",Edge[c]);
       }
       }
     */
    return I;
  }
  
  int TetsurfInit(PyMOLGlobals * G)
  {
    G->Tetsurf = TetsurfNew(G);
    return 1;
  }
  
  
  /*===========================================================================*/
  static void _TetsurfFree(CTetsurf * I)
  {
    FreeP(I);
  }
  
  void TetsurfFree(PyMOLGlobals * G)
  {
    _TetsurfFree(G->Tetsurf);
    G->Tetsurf = NULL;
  }
  
  
  /*===========================================================================*/
  void TetsurfGetRange(PyMOLGlobals * G,
                       Isofield * field, CCrystal * cryst, float *mn, float *mx, int *range)
  {
    float rmn[3], rmx[3];
    float imn[3], imx[3];
    float mix[24], imix[24];
    int a, b;
    PRINTFD(G, FB_Isosurface)
      " IsosurfGetRange: entered mn: %4.2f %4.2f %4.2f mx: %4.2f %4.2f %4.2f\n",
      mn[0], mn[1], mn[2], mx[0], mx[1], mx[2]
      ENDFD;
  
    for(a = 0; a < 3; a++) {
      rmn[a] = F4(field->points, 0, 0, 0, a);
      rmx[a] = F4(field->points,
                  field->dimensions[0] - 1,
                  field->dimensions[1] - 1, field->dimensions[2] - 1, a);
    }
  
    /* get min/max extents of map in fractional space */
  
    transform33f3f(cryst->RealToFrac, rmn, imn);
    transform33f3f(cryst->RealToFrac, rmx, imx);
  
    mix[0] = mn[0];
    mix[1] = mn[1];
    mix[2] = mn[2];
  
    mix[3] = mx[0];
    mix[4] = mn[1];
    mix[5] = mn[2];
  
    mix[6] = mn[0];
    mix[7] = mx[1];
    mix[8] = mn[2];
  
    mix[9] = mn[0];
    mix[10] = mn[1];
    mix[11] = mx[2];
  
    mix[12] = mx[0];
    mix[13] = mx[1];
    mix[14] = mn[2];
  
    mix[15] = mx[0];
    mix[16] = mn[1];
    mix[17] = mx[2];
  
    mix[18] = mn[0];
    mix[19] = mx[1];
    mix[20] = mx[2];
  
    mix[21] = mx[0];
    mix[22] = mx[1];
    mix[23] = mx[2];
  
    /* compute min/max of query in fractional space */
  
    for(b = 0; b < 8; b++) {
      transform33f3f(cryst->RealToFrac, mix + 3 * b, imix + 3 * b);
    }
  
    for(a = 0; a < 3; a++) {
      if(imx[a] != imn[a]) {      /* protect against div by zero */
        int b;
        int mini = 0, maxi = 0, tst_min, tst_max;
        float cur;
        for(b = 0; b < 8; b++) {
          cur =
            ((field->dimensions[a] - 1) * (imix[a + 3 * b] - imn[a]) / (imx[a] - imn[a]));
          tst_min = (int) floor(cur);
          tst_max = ((int) ceil(cur)) + 1;
  
          if(!b) {
            mini = tst_min;
            maxi = tst_max;
          } else {
            if(mini > tst_min)
              mini = tst_min;
            if(maxi <= tst_max)
              maxi = tst_max;
          }
        }
  
        range[a] = mini;
  
        range[a + 3] = maxi;
      } else {
        range[a] = 0;
        range[a + 3] = 1;
      }
      if(range[a] < 0)
        range[a] = 0;
      if(range[a] > field->dimensions[a])
        range[a] = field->dimensions[a];
      if(range[a + 3] < 0)
        range[a + 3] = 0;
      if(range[a + 3] > field->dimensions[a])
        range[a + 3] = field->dimensions[a];
    }
    PRINTFD(G, FB_Isosurface)
      " IsosurfGetRange: returning range: %d %d %d %d %d %d\n",
      range[0], range[1], range[2], range[3], range[4], range[5]
      ENDFD;
  #if 0
    float fmn[3], fmx[3];
    float rmn[3], rmx[3];
    float imn[3], imx[3];
    int a;
    transform33f3f(cryst->RealToFrac, mn, fmn);
    transform33f3f(cryst->RealToFrac, mx, fmx);
    for(a = 0; a < 3; a++) {
      rmn[a] = F4(field->points, 0, 0, 0, a);
      rmx[a] = F4(field->points,
                  field->dimensions[0] - 1,
                  field->dimensions[1] - 1, field->dimensions[2] - 1, a);
    }
  
    transform33f3f(cryst->RealToFrac, rmn, imn);
    transform33f3f(cryst->RealToFrac, rmx, imx);
  
    for(a = 0; a < 3; a++) {
      range[a] =
        floor(((field->dimensions[a] - 1) * (fmn[a] - imn[a]) / (imx[a] - imn[a])));
      range[a + 3] =
        (int) (((field->dimensions[a] - 1) * (fmx[a] - imn[a]) / (imx[a] - imn[a]))) + 1;
      if(range[a] < 0)
        range[a] = 0;
      if(range[a] > field->dimensions[a])
        range[a] = field->dimensions[a];
      if(range[a + 3] < 0)
        range[a + 3] = 0;
      if(range[a + 3] > field->dimensions[a])
        range[a + 3] = field->dimensions[a];
    }
  #endif
  
  }
  
  
  /*===========================================================================*/
  int TetsurfVolume(PyMOLGlobals * G, Isofield * field, float level, int **num,
                    float **vert, int *range, int mode, MapType * voxelmap, float *a_vert,
                    float carvebuffer, int side)
  {
  
    register CTetsurf *I;
    if(PIsGlutThread()) {
      I = G->Tetsurf;
    } else {
      I = TetsurfNew(G);
    }
  
    {
      int ok = true;
      int Steps[3];
      int c, i, j, k;
      int range_store[6];
      int n_strip = 0;
      int n_vert = 0;
      int tot_prim = 0;
  
      if(mode == 3)
        IsofieldComputeGradients(G, field);
  
      I->TotPrim = 0;
      if(range) {
        for(c = 0; c < 3; c++) {
          I->AbsDim[c] = field->dimensions[c];
          I->CurDim[c] = TetsurfSubSize + 1;
          Steps[c] = 1 + ((range[3 + c] - range[c]) - 1) / TetsurfSubSize;
        }
      } else {
        range = range_store;
        for(c = 0; c < 3; c++) {
          range[c] = 0;
          range[3 + c] = field->dimensions[c];
          I->AbsDim[c] = field->dimensions[c];
          I->CurDim[c] = TetsurfSubSize + 1;
          Steps[c] = 1 + (I->AbsDim[c] - 1) / TetsurfSubSize;
        }
      }
  
      /*   for(c=0;c<3;c++) {
         printf("range %d %d %d\n",c,range[c],range[c+3]);
         printf("steps %d\n",Steps[c]);
         }
       */
  
      I->Coord = field->points;
      I->Grad = field->gradients;
      I->Data = field->data;
      I->Level = level;
      if(ok)
        ok = TetsurfAlloc(I);
  
      if(ok) {
  
        for(i = 0; i < Steps[0]; i++)
          for(j = 0; j < Steps[1]; j++)
            for(k = 0; k < Steps[2]; k++) {
              I->CurOff[0] = TetsurfSubSize * i;
              I->CurOff[1] = TetsurfSubSize * j;
              I->CurOff[2] = TetsurfSubSize * k;
              for(c = 0; c < 3; c++)
                I->CurOff[c] += range[c];
              for(c = 0; c < 3; c++) {
                I->Max[c] = (range[3 + c] - I->CurOff[c]);
                if(I->Max[c] > (TetsurfSubSize + 1))
                  I->Max[c] = (TetsurfSubSize + 1);
              }
              /*         
                 for(c=0;c<3;c++)
                 printf(" TetsurfVolume: c: %i I->CurOff[c]: %i I->Max[c] %i\n",c,I->CurOff[c],I->Max[c]); 
               */
  
              if(ok) {
                if(TetsurfCodeVertices(I))
                  n_vert = TetsurfFindActiveBoxes(I, mode, &n_strip, n_vert, num, vert,
                                                  voxelmap, a_vert, carvebuffer, side);
              }
            }
        TetsurfPurge(I);
      }
  
      if(Feedback(G, FB_Isosurface, FB_Blather)) {
        if(mode < 2) {
          printf(" TetsurfVolume: Surface generated using %d vertices.\n", n_vert);
        } else {
          printf(" TetsurfVolume: Surface generated using %d triangles.\n", I->TotPrim);
        }
      }
  
      /* sentinel strip (0 length) */
  
      VLACheck(*num, int, n_strip);
      (*num)[n_strip] = 0;
      (n_strip)++;
  
      /* shrinks sizes for more efficient RAM usage */
  
      VLASize(*vert, float, n_vert * 3);
      VLASize(*num, int, n_strip);
  
      tot_prim = I->TotPrim;
  
      if(!PIsGlutThread()) {
        _TetsurfFree(I);
      }
      return (tot_prim);
    }
  }
  
  
  /*===========================================================================*/
  static int TetsurfAlloc(CTetsurf * II)
  {
    register CTetsurf *I = II;
    PyMOLGlobals *G = I->G;
  
    int ok = true;
    int dim4[4];
    int a;
    for(a = 0; a < 3; a++)
      dim4[a] = I->CurDim[a];
    dim4[3] = 3;
  
    I->VertexCodes = FieldNew(G, I->CurDim, 3, sizeof(int), cFieldInt);
    ErrChkPtr(G, I->VertexCodes);
    I->ActiveEdges = FieldNew(G, I->CurDim, 3, sizeof(int), cFieldInt);
    ErrChkPtr(G, I->ActiveEdges);
  
    dim4[3] = 7;                  /* seven different ways now... */
    I->Point = FieldNew(G, dim4, 4, sizeof(PointType), cFieldOther);
    ErrChkPtr(G, I->Point);
  
    I->Tri = VLAlloc(TriangleType, 50000);
    I->PtLink = VLAlloc(PointLinkType, 50000);
  
    if(!(I->VertexCodes && I->ActiveEdges && I->Point)) {
      TetsurfPurge(I);
      ok = false;
    }
  #ifdef Trace
    printf(" TetsurfAlloc: ok: %i\n", ok);
  #endif
    return (ok);
  }
  
  
  /*===========================================================================*/
  static void TetsurfPurge(CTetsurf * II)
  {
    register CTetsurf *I = II;
    if(I->Tri) {
      VLAFreeP(I->Tri);
    }
    if(I->PtLink);
    {
      VLAFreeP(I->PtLink);
    }
    if(I->VertexCodes) {
      FieldFree(I->VertexCodes);
      I->VertexCodes = NULL;
    }
    if(I->ActiveEdges) {
      FieldFree(I->ActiveEdges);
      I->ActiveEdges = NULL;
    }
    if(I->Point) {
      FieldFree(I->Point);
      I->Point = NULL;
    }
  }
  
  
  /*===========================================================================*/
  #ifdef _PYMOL_INLINE
  __inline__
  #endif
  static void TetsurfInterpolate2(float *pt, float *v0, float l0, float *v1, float l1,
                                  float level)
  {
    float ratio;
    ratio = (level - l0) / (l1 - l0);
    pt[0] = v0[0] + (v1[0] - v0[0]) * ratio;
    pt[1] = v0[1] + (v1[1] - v0[1]) * ratio;
    pt[2] = v0[2] + (v1[2] - v0[2]) * ratio;
  }
  
  
  /*===========================================================================*/
  #ifdef _PYMOL_INLINE
  __inline__
  #endif
  static void TetsurfInterpolate4(float *pt, float *v0, float l0, float *v1, float l1,
                                  float l2, float l3, float level)
  {
    float ratio;
    float v[3], l;
    average3f(v0, v1, v);
    l = (l0 + l1 + l2 + l3) * 0.25F;
    if(((l > level) && (l1 > level)) || ((l <= level) && (l0 > level))) { /* l0 vs l */
      ratio = (level - l0) / (l - l0);
      pt[0] = v0[0] + (v[0] - v0[0]) * ratio;
      pt[1] = v0[1] + (v[1] - v0[1]) * ratio;
      pt[2] = v0[2] + (v[2] - v0[2]) * ratio;
    } else {
      ratio = (level - l1) / (l - l1);
      pt[0] = v1[0] + (v[0] - v1[0]) * ratio;
      pt[1] = v1[1] + (v[1] - v1[1]) * ratio;
      pt[2] = v1[2] + (v[2] - v1[2]) * ratio;
    }
  }
  
  
  /*===========================================================================*/
  #ifdef _PYMOL_INLINE
  __inline__
  #endif
  static void TetsurfInterpolate8(float *pt, float *v0, float l0, float *v1, float l1,
                                  float l2, float l3, float l4,
                                  float l5, float l6, float l7, float level)
  {
    float ratio;
    float v[3], l;
    average3f(v0, v1, v);
    l = (l0 + l1 + l2 + l3 + l4 + l5 + l6 + l7) * 0.125F;
    if(((l > level) && (l1 > level)) || ((l <= level) && (l0 > level))) { /* l0 vs l */
      ratio = (level - l0) / (l - l0);
      pt[0] = v0[0] + (v[0] - v0[0]) * ratio;
      pt[1] = v0[1] + (v[1] - v0[1]) * ratio;
      pt[2] = v0[2] + (v[2] - v0[2]) * ratio;
    } else {
      ratio = (level - l1) / (l - l1);
      pt[0] = v1[0] + (v[0] - v1[0]) * ratio;
      pt[1] = v1[1] + (v[1] - v1[1]) * ratio;
      pt[2] = v1[2] + (v[2] - v1[2]) * ratio;
    }
  }
  
  
  /*===========================================================================*/
  static int TetsurfFindActiveBoxes(CTetsurf * II, int mode, int *n_strip, int n_vert,
                                    int **strip_l, float **vert,
                                    MapType * voxelmap, float *a_vert,
                                    float carvebuffer, int side)
  {
    register CTetsurf *I = II;
    int a, b, c, i, j, k, h, l;
  #ifdef Trace
    int ECount = 0;
  #endif
    int i000, i001, i010, i011, i100, i101, i110, i111;
    float *c000, *c001, *c010, *c011, *c100, *c101, *c110, *c111;
    float d000, d001, d010, d011, d100, d101, d110, d111;
    float *g000 = NULL, *g001 = NULL, *g010 = NULL, *g011 = NULL, *g100 = NULL, *g101 =
      NULL, *g110 = NULL, *g111 = NULL;
  
    int active;
    int n_active = 0;
    int n_start = 0;
    PointType *e[19], *p0, *p1, *p2;
    int code;
    int eidx;
    int idx;
    TriangleType *tt;
    int n_tri = 0;
    int n_link = 1;
  
    int avoid_flag = false;
    if(carvebuffer < 0.0F) {
      avoid_flag = true;
      carvebuffer = -carvebuffer;
    }
  
    FieldZero(I->Point);          /* sets initial links to zero */
    FieldZero(I->ActiveEdges);
    n_start = n_vert;
    for(i = 0; i < (I->Max[0] - 1); i++)
      for(j = 0; j < (I->Max[1] - 1); j++)
        for(k = 0; k < (I->Max[2] - 1); k++) {
          active = 0;
  
          i000 = I3(I->VertexCodes, i, j, k);
          i001 = I3(I->VertexCodes, i, j, k + 1);
          i010 = I3(I->VertexCodes, i, j + 1, k);
          i011 = I3(I->VertexCodes, i, j + 1, k + 1);
          i100 = I3(I->VertexCodes, i + 1, j, k);
          i101 = I3(I->VertexCodes, i + 1, j, k + 1);
          i110 = I3(I->VertexCodes, i + 1, j + 1, k);
          i111 = I3(I->VertexCodes, i + 1, j + 1, k + 1);
  
          if((i000 != i001) || (i001 != i010) || (i010 != i011) || (i011 != i100) || (i100 != i101) || (i101 != i110) || (i110 != i111)) {        /* this is an active box */
  
            c000 = O4Ptr(I->Coord, i, j, k, 0, I->CurOff);
            c001 = O4Ptr(I->Coord, i, j, k + 1, 0, I->CurOff);
            c010 = O4Ptr(I->Coord, i, j + 1, k, 0, I->CurOff);
            c011 = O4Ptr(I->Coord, i, j + 1, k + 1, 0, I->CurOff);
            c100 = O4Ptr(I->Coord, i + 1, j, k, 0, I->CurOff);
            c101 = O4Ptr(I->Coord, i + 1, j, k + 1, 0, I->CurOff);
            c110 = O4Ptr(I->Coord, i + 1, j + 1, k, 0, I->CurOff);
            c111 = O4Ptr(I->Coord, i + 1, j + 1, k + 1, 0, I->CurOff);
  
            if(mode == 3) {
  
              g000 = O4Ptr(I->Grad, i, j, k, 0, I->CurOff);
              g001 = O4Ptr(I->Grad, i, j, k + 1, 0, I->CurOff);
              g010 = O4Ptr(I->Grad, i, j + 1, k, 0, I->CurOff);
              g011 = O4Ptr(I->Grad, i, j + 1, k + 1, 0, I->CurOff);
              g100 = O4Ptr(I->Grad, i + 1, j, k, 0, I->CurOff);
              g101 = O4Ptr(I->Grad, i + 1, j, k + 1, 0, I->CurOff);
              g110 = O4Ptr(I->Grad, i + 1, j + 1, k, 0, I->CurOff);
              g111 = O4Ptr(I->Grad, i + 1, j + 1, k + 1, 0, I->CurOff);
            }
  
            d000 = O3(I->Data, i, j, k, I->CurOff);
            d001 = O3(I->Data, i, j, k + 1, I->CurOff);
            d010 = O3(I->Data, i, j + 1, k, I->CurOff);
            d011 = O3(I->Data, i, j + 1, k + 1, I->CurOff);
            d100 = O3(I->Data, i + 1, j, k, I->CurOff);
            d101 = O3(I->Data, i + 1, j, k + 1, I->CurOff);
            d110 = O3(I->Data, i + 1, j + 1, k, I->CurOff);
            d111 = O3(I->Data, i + 1, j + 1, k + 1, I->CurOff);
  
            e[cE_000_001] = EdgePtPtr(I->Point, i, j, k, cE_000_001);
            e[cE_000_010] = EdgePtPtr(I->Point, i, j, k, cE_000_010);
            e[cE_000_011] = EdgePtPtr(I->Point, i, j, k, cE_000_011);
            e[cE_000_100] = EdgePtPtr(I->Point, i, j, k, cE_000_100);
            e[cE_000_101] = EdgePtPtr(I->Point, i, j, k, cE_000_101);
  
            e[cE_000_110] = EdgePtPtr(I->Point, i, j, k, cE_000_110);
            e[cE_000_111] = EdgePtPtr(I->Point, i, j, k, cE_000_111);
            e[cE_001_011] = EdgePtPtr(I->Point, i, j, k + 1, cE_000_010);
            e[cE_001_101] = EdgePtPtr(I->Point, i, j, k + 1, cE_000_100);
            e[cE_001_111] = EdgePtPtr(I->Point, i, j, k + 1, cE_000_110);
  
            e[cE_010_011] = EdgePtPtr(I->Point, i, j + 1, k, cE_000_001);
            e[cE_010_110] = EdgePtPtr(I->Point, i, j + 1, k, cE_000_100);
            e[cE_010_111] = EdgePtPtr(I->Point, i, j + 1, k, cE_000_101);
            e[cE_100_101] = EdgePtPtr(I->Point, i + 1, j, k, cE_000_001);
            e[cE_100_110] = EdgePtPtr(I->Point, i + 1, j, k, cE_000_010);
  
            e[cE_100_111] = EdgePtPtr(I->Point, i + 1, j, k, cE_000_011);
            e[cE_011_111] = EdgePtPtr(I->Point, i, j + 1, k + 1, cE_000_100);
            e[cE_101_111] = EdgePtPtr(I->Point, i + 1, j, k + 1, cE_000_010);
            e[cE_110_111] = EdgePtPtr(I->Point, i + 1, j + 1, k, cE_000_001);
  
            /* Generate interpolated coordinates for all active edges */
  
            if(i000 != i001) {
              if(!(I3(I->ActiveEdges, i, j, k) & cM_000_001)) {
                I3(I->ActiveEdges, i, j, k) |= cM_000_001;
                TetsurfInterpolate2(e[cE_000_001]->Point, c000, d000, c001, d001, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_000_001]->Normal, g000, d000, g001, d001,
                                      I->Level);
              }
              active |= cM_000_001;
            }
            if(i000 != i010) {
              if(!(I3(I->ActiveEdges, i, j, k) & cM_000_010)) {
                I3(I->ActiveEdges, i, j, k) |= cM_000_010;
                TetsurfInterpolate2(e[cE_000_010]->Point, c000, d000, c010, d010, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_000_010]->Normal, g000, d000, g010, d010,
                                      I->Level);
  
              }
              active |= cM_000_010;
            }
            if(i000 != i011) {
              if(!(I3(I->ActiveEdges, i, j, k) & cM_000_011)) {
                I3(I->ActiveEdges, i, j, k) |= cM_000_011;
                TetsurfInterpolate4(e[cE_000_011]->Point, c000, d000, c011, d011, d001,
                                    d010, I->Level);
                if(mode == 3)
                  TetsurfInterpolate4(e[cE_000_011]->Normal, g000, d000, g011, d011, d001,
                                      d010, I->Level);
              }
              active |= cM_000_011;
            }
            if(i000 != i100) {
              if(!(I3(I->ActiveEdges, i, j, k) & cM_000_100)) {
                I3(I->ActiveEdges, i, j, k) |= cM_000_100;
                TetsurfInterpolate2(e[cE_000_100]->Point, c000, d000, c100, d100, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_000_100]->Normal, g000, d000, g100, d100,
                                      I->Level);
  
              }
              active |= cM_000_100;
            }
            if(i000 != i101) {
              if(!(I3(I->ActiveEdges, i, j, k) & cM_000_101)) {
                I3(I->ActiveEdges, i, j, k) |= cM_000_101;
                TetsurfInterpolate4(e[cE_000_101]->Point, c000, d000, c101, d101, d100,
                                    d001, I->Level);
                if(mode == 3)
                  TetsurfInterpolate4(e[cE_000_101]->Normal, g000, d000, g101, d101, d100,
                                      d001, I->Level);
              }
              active |= cM_000_101;
            }
            if(i000 != i110) {
              if(!(I3(I->ActiveEdges, i, j, k) & cM_000_110)) {
                I3(I->ActiveEdges, i, j, k) |= cM_000_110;
                TetsurfInterpolate4(e[cE_000_110]->Point, c000, d000, c110, d110, d100,
                                    d010, I->Level);
                if(mode == 3)
                  TetsurfInterpolate4(e[cE_000_110]->Normal, g000, d000, g110, d110, d100,
                                      d010, I->Level);
              }
              active |= cM_000_110;
            }
            if(i000 != i111) {
              if(!(I3(I->ActiveEdges, i, j, k) & cM_000_111)) {
                I3(I->ActiveEdges, i, j, k) |= cM_000_111;
                TetsurfInterpolate8(e[cE_000_111]->Point,
                                    c000, d000, c111, d111,
                                    d001, d010, d011, d100, d101, d110, I->Level);
                if(mode == 3)
                  TetsurfInterpolate8(e[cE_000_111]->Normal,
                                      g000, d000, g111, d111,
                                      d001, d010, d011, d100, d101, d110, I->Level);
              }
              active |= cM_000_111;
            }
            if(i001 != i011) {
              if(!(I3(I->ActiveEdges, i, j, k + 1) & cM_000_010)) {
                I3(I->ActiveEdges, i, j, k + 1) |= cM_000_010;
                TetsurfInterpolate2(e[cE_001_011]->Point, c001, d001, c011, d011, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_001_011]->Normal, g001, d001, g011, d011,
                                      I->Level);
              }
              active |= cM_001_011;
            }
            if(i001 != i101) {
              if(!(I3(I->ActiveEdges, i, j, k + 1) & cM_000_100)) {
                I3(I->ActiveEdges, i, j, k + 1) |= cM_000_100;
                TetsurfInterpolate2(e[cE_001_101]->Point, c001, d001, c101, d101, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_001_101]->Normal, g001, d001, g101, d101,
                                      I->Level);
              }
              active |= cM_001_101;
            }
            if(i001 != i111) {
              if(!(I3(I->ActiveEdges, i, j, k + 1) & cM_000_110)) {
                I3(I->ActiveEdges, i, j, k + 1) |= cM_000_110;
                TetsurfInterpolate4(e[cE_001_111]->Point, c001, d001, c111, d111, d101,
                                    d011, I->Level);
                if(mode == 3)
                  TetsurfInterpolate4(e[cE_001_111]->Normal, g001, d001, g111, d111, d101,
                                      d011, I->Level);
              }
              active |= cM_001_111;
            }
            if(i010 != i011) {
              if(!(I3(I->ActiveEdges, i, j + 1, k) & cM_000_001)) {
                I3(I->ActiveEdges, i, j + 1, k) |= cM_000_001;
                TetsurfInterpolate2(e[cE_010_011]->Point, c010, d010, c011, d011, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_010_011]->Normal, g010, d010, g011, d011,
                                      I->Level);
              }
              active |= cM_010_011;
            }
            if(i010 != i110) {
              if(!(I3(I->ActiveEdges, i, j + 1, k) & cM_000_100)) {
                I3(I->ActiveEdges, i, j + 1, k) |= cM_000_100;
                TetsurfInterpolate2(e[cE_010_110]->Point, c010, d010, c110, d110, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_010_110]->Normal, g010, d010, g110, d110,
                                      I->Level);
              }
              active |= cM_010_110;
            }
            if(i010 != i111) {
              if(!(I3(I->ActiveEdges, i, j + 1, k) & cM_000_101)) {
                I3(I->ActiveEdges, i, j + 1, k) |= cM_000_101;
                TetsurfInterpolate4(e[cE_010_111]->Point, c010, d010, c111, d111, d110,
                                    d011, I->Level);
                if(mode == 3)
                  TetsurfInterpolate4(e[cE_010_111]->Normal, g010, d010, g111, d111, d110,
                                      d011, I->Level);
              }
              active |= cM_010_111;
            }
            if(i100 != i101) {
              if(!(I3(I->ActiveEdges, i + 1, j, k) & cM_000_001)) {
                I3(I->ActiveEdges, i + 1, j, k) |= cM_000_001;
                TetsurfInterpolate2(e[cE_100_101]->Point, c100, d100, c101, d101, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_100_101]->Normal, g100, d100, g101, d101,
                                      I->Level);
              }
              active |= cM_100_101;
            }
            if(i100 != i110) {
              if(!(I3(I->ActiveEdges, i + 1, j, k) & cM_000_010)) {
                I3(I->ActiveEdges, i + 1, j, k) |= cM_000_010;
                TetsurfInterpolate2(e[cE_100_110]->Point, c100, d100, c110, d110, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_100_110]->Normal, g100, d100, g110, d110,
                                      I->Level);
              }
              active |= cM_100_110;
            }
            if(i100 != i111) {
              if(!(I3(I->ActiveEdges, i + 1, j, k) & cM_000_011)) {
                I3(I->ActiveEdges, i + 1, j, k) |= cM_000_011;
                TetsurfInterpolate4(e[cE_100_111]->Point, c100, d100, c111, d111, d101,
                                    d110, I->Level);
                if(mode == 3)
                  TetsurfInterpolate4(e[cE_100_111]->Normal, g100, d100, g111, d111, d101,
                                      d110, I->Level);
              }
              active |= cM_100_111;
            }
            if(i011 != i111) {
              if(!(I3(I->ActiveEdges, i, j + 1, k + 1) & cM_000_100)) {
                I3(I->ActiveEdges, i, j + 1, k + 1) |= cM_000_100;
                TetsurfInterpolate2(e[cE_011_111]->Point, c011, d011, c111, d111, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_011_111]->Normal, g011, d011, g111, d111,
                                      I->Level);
              }
              active |= cM_011_111;
            }
            if(i101 != i111) {
              if(!(I3(I->ActiveEdges, i + 1, j, k + 1) & cM_000_010)) {
                I3(I->ActiveEdges, i + 1, j, k + 1) |= cM_000_010;
                TetsurfInterpolate2(e[cE_101_111]->Point, c101, d101, c111, d111, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_101_111]->Normal, g101, d101, g111, d111,
                                      I->Level);
              }
              active |= cM_101_111;
            }
            if(i110 != i111) {
              if(!(I3(I->ActiveEdges, i + 1, j + 1, k) & cM_000_001)) {
                I3(I->ActiveEdges, i + 1, j + 1, k) |= cM_000_001;
                TetsurfInterpolate2(e[cE_110_111]->Point, c110, d110, c111, d111, I->Level);
                if(mode == 3)
                  TetsurfInterpolate2(e[cE_110_111]->Normal, g110, d110, g111, d111,
                                      I->Level);
              }
              active |= cM_110_111;
            }
  
            if(active) {
              switch (mode) {
              case 2:
              case 3:
                code =
                  (i000 ? 0x01 : 0) |
                  (i001 ? 0x02 : 0) |
                  (i010 ? 0x04 : 0) |
                  (i011 ? 0x08 : 0) |
                  (i100 ? 0x10 : 0) |
                  (i101 ? 0x20 : 0) | (i110 ? 0x40 : 0) | (i111 ? 0x80 : 0);
                eidx = I->EdgeStart[code];
                while(1) {
                  idx = I->Edge[eidx];
                  if(idx < 0)
                    break;
  
                  /* assemble a triangle from these three points */
  
                  VLACheck(I->Tri, TriangleType, n_tri);
                  tt = I->Tri + n_tri;
                  tt->p[0] = e[idx];
                  tt->p[1] = e[I->Edge[eidx + 1]];
                  tt->p[2] = e[I->Edge[eidx + 2]];
  
                  VLACheck(I->PtLink, PointLinkType, n_link + 3);
  
                  /* link this triangle into the points */
  
                  I->PtLink[n_link].tri = n_tri;
                  I->PtLink[n_link].link = tt->p[0]->Link;
                  tt->p[0]->Link = n_link;
                  n_link++;
  
                  I->PtLink[n_link].tri = n_tri;
                  I->PtLink[n_link].link = tt->p[1]->Link;
                  tt->p[1]->Link = n_link;
                  n_link++;
  
                  I->PtLink[n_link].tri = n_tri;
                  I->PtLink[n_link].link = tt->p[2]->Link;
                  tt->p[2]->Link = n_link;
                  n_link++;
                  n_tri++;
                  eidx += 3;
                }
                break;
              case 1:            /* lines */
                VLACheck(*vert, float, (n_vert * 3) + 200);
  
                code =
                  (i000 ? 0x01 : 0) |
                  (i001 ? 0x02 : 0) |
                  (i010 ? 0x04 : 0) |
                  (i011 ? 0x08 : 0) |
                  (i100 ? 0x10 : 0) |
                  (i101 ? 0x20 : 0) | (i110 ? 0x40 : 0) | (i111 ? 0x80 : 0);
                eidx = I->EdgeStart[code];
                while(1) {
                  idx = I->Edge[eidx];
                  if(idx < 0)
                    break;
                  copy3fn(e[idx]->Point, (*vert) + (n_vert * 3));
                  n_vert++;
                  copy3fn(e[I->Edge[eidx + 1]]->Point, (*vert) + (n_vert * 3));
                  n_vert++;
                  copy3fn(e[I->Edge[eidx + 1]]->Point, (*vert) + (n_vert * 3));
                  n_vert++;
                  copy3fn(e[I->Edge[eidx + 2]]->Point, (*vert) + (n_vert * 3));
                  n_vert++;
                  copy3fn(e[I->Edge[eidx + 2]]->Point, (*vert) + (n_vert * 3));
                  n_vert++;
                  copy3fn(e[idx]->Point, (*vert) + (n_vert * 3));
                  n_vert++;
                  eidx += 3;
                }
                break;
              case 0:
              default:           /* dots */
                VLACheck(*vert, float, (n_vert * 3) + 200);
  
                if(active & cM_000_001) {
                  copy3fn(e[cE_000_001]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_000_010) {
                  copy3fn(e[cE_000_010]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_000_011) {
                  copy3fn(e[cE_000_011]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_000_100) {
                  copy3fn(e[cE_000_100]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_000_101) {
                  copy3fn(e[cE_000_101]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_000_110) {
                  copy3fn(e[cE_000_110]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_000_111) {
                  copy3fn(e[cE_000_111]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
  
                if(active & cM_001_011) {
                  copy3fn(e[cE_001_011]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_001_101) {
                  copy3fn(e[cE_001_101]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_001_111) {
                  copy3fn(e[cE_001_111]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
  
                if(active & cM_010_011) {
                  copy3fn(e[cE_010_011]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_010_110) {
                  copy3fn(e[cE_010_011]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_010_111) {
                  copy3fn(e[cE_010_111]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
  
                if(active & cM_100_101) {
                  copy3fn(e[cE_100_101]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_100_110) {
                  copy3fn(e[cE_100_110]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_100_111) {
                  copy3fn(e[cE_100_111]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
  
                if(active & cM_011_111) {
                  copy3fn(e[cE_011_111]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_101_111) {
                  copy3fn(e[cE_101_111]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
                if(active & cM_110_111) {
                  copy3fn(e[cE_101_111]->Point, (*vert) + (n_vert * 3));
                  n_vert += 1;
                }
  
                break;
              }
              n_active++;
            }
          }
        }
  
    switch (mode) {
    case 2:
    case 3:
      /* compute area-weighted normal */
      for(a = 0; a < n_tri; a++) {
        float *v0, *v1, *v2;
        float vt1[3], vt2[3];
  
        tt = I->Tri + a;
        v0 = tt->p[0]->Point;
        v1 = tt->p[1]->Point;
        v2 = tt->p[2]->Point;
        tt->done = false;         /* init */
  
        subtract3f(v0, v2, vt1);
        subtract3f(v1, v2, vt2);
        if(side >= 0) {
          cross_product3f(vt2, vt1, tt->n);
        } else {
          cross_product3f(vt1, vt2, tt->n);
        }
      }
      /* compute (or lookup) normals at active points */
      for(a = 0; a < n_tri; a++) {
        float v[3];
        float dp;
        tt = I->Tri + a;
        for(b = 0; b < 3; b++) {
          if(!tt->p[b]->NormalFlag) {
            zero3f(v);
            idx = tt->p[b]->Link;
            while(idx > 0) {
              add3f(I->Tri[I->PtLink[idx].tri].n, v, v);
              idx = I->PtLink[idx].link;
            }
            if(mode == 3) {       /* gradient-based normals */
              dp = dot_product3f(v, tt->p[b]->Normal);
              if(dp < 0.0F) {
                invert3f(tt->p[b]->Normal);
                normalize3f(tt->p[b]->Normal);
              } else if(dp == 0.0F) {     /* fall back on triangle normal */
                normalize23f(v, tt->p[b]->Normal);
              } else {
                normalize3f(tt->p[b]->Normal);
              }
            } else {              /* triangle-based normals */
              normalize23f(v, tt->p[b]->Normal);
            }
            tt->p[b]->NormalFlag = true;
          }
        }
      }
      if(mode == 2) {
        /* if we're using triangle normals, then 
           do an additional averaging cycle with no weighting */
        for(a = 0; a < n_tri; a++) {
          tt = I->Tri + a;
          add3f(tt->p[0]->Normal, tt->p[1]->Normal, tt->n);
          add3f(tt->p[2]->Normal, tt->n, tt->n);
          normalize3f(tt->n);
          tt->p[0]->NormalFlag = false;
          tt->p[1]->NormalFlag = false;
          tt->p[2]->NormalFlag = false;
        }
        /* compute normals at active points */
        for(a = 0; a < n_tri; a++) {
          float v[3];
          tt = I->Tri + a;
          for(b = 0; b < 3; b++) {
            if(!tt->p[b]->NormalFlag) {
              zero3f(v);
              idx = tt->p[b]->Link;
              while(idx > 0) {
                add3f(I->Tri[I->PtLink[idx].tri].n, v, v);
                idx = I->PtLink[idx].link;
              }
              normalize23f(v, tt->p[b]->Normal);
              tt->p[b]->NormalFlag = true;
            }
          }
        }
      }
  
      /* Need to move the points now, right? */
  
      /* if we are carving, then exclude triangles outside region */
      if(voxelmap) {
        for(a = 0; a < n_tri; a++) {
          float *v;
          tt = I->Tri + a;
          c = 0;
          for(b = 0; b < 3; b++) {
            v = tt->p[b]->Point;
            MapLocus(voxelmap, v, &h, &k, &l);
            i = *(MapEStart(voxelmap, h, k, l));
            if(i) {
              j = voxelmap->EList[i++];
              while(j >= 0) {
                if(within3f(a_vert + 3 * j, v, carvebuffer)) {
                  c++;
                  break;
                }
                j = voxelmap->EList[i++];
              }
            }
          }
          if(avoid_flag) {
            if(c >= 3)
              tt->done = true;    /* exclude this triangle from the surface */
          } else {
            if(c < 3)             /* exclude this triangle from the surface */
              tt->done = true;
          }
        }
      }
  
      /* now create triangle strips (not yet optimal) */
      for(a = 0; a < n_tri; a++) {
        tt = I->Tri + a;
        n_start = n_vert;
        if(!tt->done) {
  
          VLACheck(*vert, float, (n_vert * 3) + 200);
  
          if(side >= 0) {
            /* switch order around to get "correct" triangles */
  
            copy3fn(tt->p[1]->Normal, (*vert) + (n_vert * 3));
            n_vert++;
            copy3fn(tt->p[1]->Point, (*vert) + (n_vert * 3));
            n_vert++;
  
            copy3fn(tt->p[0]->Normal, (*vert) + (n_vert * 3));
            n_vert++;
            copy3fn(tt->p[0]->Point, (*vert) + (n_vert * 3));
            n_vert++;
  
            copy3fn(tt->p[2]->Normal, (*vert) + (n_vert * 3));
            n_vert++;
            copy3fn(tt->p[2]->Point, (*vert) + (n_vert * 3));
            n_vert++;
  
            p0 = tt->p[0];
            p1 = tt->p[2];
  
          } else {
  
            copy3fn(tt->p[0]->Normal, (*vert) + (n_vert * 3));
            n_vert++;
            copy3fn(tt->p[0]->Point, (*vert) + (n_vert * 3));
            n_vert++;
  
            copy3fn(tt->p[1]->Normal, (*vert) + (n_vert * 3));
            n_vert++;
            copy3fn(tt->p[1]->Point, (*vert) + (n_vert * 3));
            n_vert++;
  
            copy3fn(tt->p[2]->Normal, (*vert) + (n_vert * 3));
            n_vert++;
            copy3fn(tt->p[2]->Point, (*vert) + (n_vert * 3));
            n_vert++;
  
            p0 = tt->p[1];
            p1 = tt->p[2];
  
          }
  
          tt->done = true;
  
          while(1) {
            p2 = NULL;
            idx = p1->Link;
            while(idx > 0) {
              tt = I->Tri + I->PtLink[idx].tri;
  
              if(!tt->done) {
                if((tt->p[0] == p0) && (tt->p[1] == p1)) {
                  p2 = tt->p[2];
                  break;
                }
  
                if((tt->p[1] == p0) && (tt->p[2] == p1)) {
                  p2 = tt->p[0];
                  break;
                }
  
                if((tt->p[2] == p0) && (tt->p[0] == p1)) {
                  p2 = tt->p[1];
                  break;
                }
  
                if((tt->p[1] == p0) && (tt->p[0] == p1)) {
                  p2 = tt->p[2];
                  break;
                }
  
                if((tt->p[2] == p0) && (tt->p[1] == p1)) {
                  p2 = tt->p[0];
                  break;
                }
  
                if((tt->p[0] == p0) && (tt->p[2] == p1)) {
                  p2 = tt->p[1];
                  break;
                }
              }
              idx = I->PtLink[idx].link;
            }
  
            if(!p2)
              break;
            tt->done = true;
            VLACheck(*vert, float, (n_vert * 3) + 200);
            copy3fn(p2->Normal, (*vert) + (n_vert * 3));
            n_vert++;
            copy3fn(p2->Point, (*vert) + (n_vert * 3));
            n_vert++;
            p0 = p1;
            p1 = p2;
          }
        }
        if(n_vert > n_start) {
          VLACheck(*strip_l, int, *n_strip);
          (*strip_l)[*n_strip] = n_vert - n_start;
          (*n_strip)++;
        }
      }
      I->TotPrim += n_tri;
      break;
  
    case 1:
      /* Need to move the points now, right? */
  
      if(n_vert > n_start) {
        VLACheck(*strip_l, int, *n_strip);
        (*strip_l)[*n_strip] = n_vert - n_start;
        (*n_strip)++;
      }
      break;
    case 0:                      /* dots */
    default:
  
      /* Need to move the points now, right? */
  
      if(n_vert > n_start) {
        VLACheck(*strip_l, int, *n_strip);
        (*strip_l)[*n_strip] = n_vert - n_start;
        (*n_strip)++;
      }
      break;
    }
    /*
       printf("n_strip %d\n",*n_strip);
       printf("n_active %d\n",n_active);
       printf("n_vert %d\n",n_vert);
       printf("mode %d\n",mode);
     */
    return (n_vert);
  }
  
  
  /*===========================================================================*/
  static int TetsurfCodeVertices(CTetsurf * II)
  {
    register CTetsurf *I = II;
    int i, j, k;
    int b0, b1;
    int flag1 = false;
    int flag2 = false;
    b0 = 1;
    if(I->Level < 0.0F)
      b0 = 0;
    b1 = 1 - b0;
  
    for(i = 0; i < I->Max[0]; i++)
      for(j = 0; j < I->Max[1]; j++)
        for(k = 0; k < I->Max[2]; k++) {
          if((O3(I->Data, i, j, k, I->CurOff) > I->Level)) {
            I3(I->VertexCodes, i, j, k) = b0;
            flag1 = true;
          } else {
            I3(I->VertexCodes, i, j, k) = b1;
            flag2 = true;
          }
        }
    return (flag1 && flag2);
  }
  