svn merge -r 12716:12856 https://svn.blender.org/svnroot/bf-blender/trunk/blender
[blender.git] / intern / elbeem / intern / isosurface.cpp
index 92ed8c63eaba51106c1c2d0f64592313e5e7e6a9..14aa22d06b59d4f4c1640004bc73c1a10baec028 100644 (file)
@@ -1,40 +1,54 @@
 /******************************************************************************
  *
  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
- * Copyright 2005 Nils Thuerey
+ * Copyright 2003-2006 Nils Thuerey
  *
  * Marching Cubes surface mesh generation
  *
  *****************************************************************************/
 
-// IMPL ----------------------------------------------------------------------------------------
 #include "isosurface.h"
 #include "mcubes_tables.h"
-#include "ntl_scene.h"
-
+#include "particletracer.h"
 #include <algorithm>
 #include <stdio.h>
-#define MCUBES_MAXPOLNUM  10000
-#define MCUBES_MAXVERTNUM 30000
 
+#if !defined(linux) && (defined (__sparc) || defined (__sparc__))
+#include <ieeefp.h>
+#endif
+
+
+// just use default rounding for platforms where its not available
+#ifndef round
+#define round(x) (x)
+#endif
 
+#if PARALLEL==1        
+#include <omp.h>
+#endif
 
 /******************************************************************************
  * Constructor
  *****************************************************************************/
-IsoSurface::IsoSurface(double iso, double blend) :
+IsoSurface::IsoSurface(double iso) :
        ntlGeometryObject(),
        mSizex(-1), mSizey(-1), mSizez(-1),
+       mpData(NULL),
   mIsoValue( iso ), 
-       mBlendVal( blend ),
        mPoints(), 
+       mUseFullEdgeArrays(false),
        mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL),
+       mEdgeArSize(-1),
        mIndices(),
 
   mStart(0.0), mEnd(0.0), mDomainExtent(0.0),
   mInitDone(false),
-       mLoopSubdivs(0), mSmoothSurface(0.0), mSmoothNormals(0.0),
-       mAcrossEdge(), mAdjacentFaces()
+       mSmoothSurface(0.0), mSmoothNormals(0.0),
+       mAcrossEdge(), mAdjacentFaces(),
+       mCutoff(-1), mCutArray(NULL), // off by default
+       mpIsoParts(NULL), mPartSize(0.), mSubdivs(0),
+       mFlagCnt(1),
+       mSCrad1(0.), mSCrad2(0.), mSCcenter(0.)
 {
 }
 
@@ -44,6 +58,10 @@ IsoSurface::IsoSurface(double iso, double blend) :
  *****************************************************************************/
 void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent) 
 {
+       // range 1-10 (max due to subd array in triangulate)
+       if(mSubdivs<1) mSubdivs=1;
+       if(mSubdivs>10) mSubdivs=10;
+
        // init solver and size
        mSizex = setx;
        mSizey = sety;
@@ -69,43 +87,53 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e
   for(int i=0;i<nodes;i++) { mpData[i] = 0.0; }
 
   // allocate edge arrays  (last slices are never used...)
-  mpEdgeVerticesX = new int[nodes];
-  mpEdgeVerticesY = new int[nodes];
-  mpEdgeVerticesZ = new int[nodes];
-  for(int i=0;i<nodes;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
+       int initsize = -1;
+       if(mUseFullEdgeArrays) {
+               mEdgeArSize = nodes;
+               mpEdgeVerticesX = new int[nodes];
+               mpEdgeVerticesY = new int[nodes];
+               mpEdgeVerticesZ = new int[nodes];
+               initsize = 3*nodes;
+       } else {
+               int sliceNodes = 2*mSizex*mSizey*mSubdivs*mSubdivs;
+               mEdgeArSize = sliceNodes;
+               mpEdgeVerticesX = new int[sliceNodes];
+               mpEdgeVerticesY = new int[sliceNodes];
+               mpEdgeVerticesZ = new int[sliceNodes];
+               initsize = 3*sliceNodes;
+       }
+  for(int i=0;i<mEdgeArSize;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
+       // WARNING - make sure this is consistent with calculateMemreqEstimate
   
        // marching cubes are ready 
        mInitDone = true;
-       unsigned long int memCnt = (3*sizeof(int)*nodes + sizeof(float)*nodes);
-       double memd = memCnt;
-       char *sizeStr = "";
-       const double sfac = 1000.0;
-       if(memd>sfac){ memd /= sfac; sizeStr="KB"; }
-       if(memd>sfac){ memd /= sfac; sizeStr="MB"; }
-       if(memd>sfac){ memd /= sfac; sizeStr="GB"; }
-       if(memd>sfac){ memd /= sfac; sizeStr="TB"; }
-
-       debMsgStd("IsoSurface::initializeIsosurface",DM_MSG,"Inited "<<PRINT_VEC(setx,sety,setz)<<" alloced:"<< memd<<" "<<sizeStr<<"." ,10);
+       debMsgStd("IsoSurface::initializeIsosurface",DM_MSG,"Inited, edgenodes:"<<initsize<<" subdivs:"<<mSubdivs<<" fulledg:"<<mUseFullEdgeArrays , 10);
 }
 
 
 
+/*! Reset all values */
+void IsoSurface::resetAll(gfxReal val) {
+       int nodes = mSizez*mSizey*mSizex;
+  for(int i=0;i<nodes;i++) { mpData[i] = val; }
+}
+
 
 /******************************************************************************
  * Destructor
  *****************************************************************************/
 IsoSurface::~IsoSurface( void )
 {
-
+       if(mpData) delete [] mpData;
        if(mpEdgeVerticesX) delete [] mpEdgeVerticesX;
        if(mpEdgeVerticesY) delete [] mpEdgeVerticesY;
        if(mpEdgeVerticesZ) delete [] mpEdgeVerticesZ;
-
 }
 
 
 
 
+
 /******************************************************************************
  * triangulate the scalar field given by pointer
  *****************************************************************************/
@@ -121,29 +149,22 @@ void IsoSurface::triangulate( void )
                return;
        }
 
-  // get grid spacing
-  gsx = (mEnd[0]-mStart[0])/(double)(mSizex-1.0);
-  gsy = (mEnd[1]-mStart[1])/(double)(mSizey-1.0);
-  gsz = (mEnd[2]-mStart[2])/(double)(mSizez-1.0);
+  // get grid spacing (-2 to have same spacing as sim)
+  gsx = (mEnd[0]-mStart[0])/(double)(mSizex-2.0);
+  gsy = (mEnd[1]-mStart[1])/(double)(mSizey-2.0);
+  gsz = (mEnd[2]-mStart[2])/(double)(mSizez-2.0);
 
   // clean up previous frame
        mIndices.clear();
        mPoints.clear();
 
        // reset edge vertices
-  for(int i=0;i<(mSizex*mSizey*mSizez);i++) {
+  for(int i=0;i<mEdgeArSize;i++) {
                mpEdgeVerticesX[i] = -1;
                mpEdgeVerticesY[i] = -1;
                mpEdgeVerticesZ[i] = -1;
        }
 
-       ntlVec3Gfx pos[8];
-       float value[8];
-       int cubeIndex;      // index entry of the cube 
-       int triIndices[12]; // vertex indices 
-       int *eVert[12];
-       IsoLevelVertex ilv;
-
        // edges between which points?
        const int mcEdges[24] = { 
                0,1,  1,2,  2,3,  3,0,
@@ -157,182 +178,500 @@ void IsoSurface::triangulate( void )
        const int cubieOffsetZ[8] = {
                0,0,0,0,  1,1,1,1 };
 
+       const int coAdd=2;
   // let the cubes march 
-       pz = mStart[2]-gsz;
-       for(int k=0;k<(mSizez-2);k++) {
-               pz += gsz;
-               py = mStart[1]-gsy;
-               for(int j=0;j<(mSizey-2);j++) {
-      py += gsy;
-                       px = mStart[0]-gsx;
-                       for(int i=0;i<(mSizex-2);i++) {
-                       px += gsx;
-                               int baseIn = ISOLEVEL_INDEX( i+0, j+0, k+0);
-
-                               value[0] = *getData(i  ,j  ,k  );
-                               value[1] = *getData(i+1,j  ,k  );
-                               value[2] = *getData(i+1,j+1,k  );
-                               value[3] = *getData(i  ,j+1,k  );
-                               value[4] = *getData(i  ,j  ,k+1);
-                               value[5] = *getData(i+1,j  ,k+1);
-                               value[6] = *getData(i+1,j+1,k+1);
-                               value[7] = *getData(i  ,j+1,k+1);
-                               //errMsg("ISOT2D"," at "<<PRINT_IJK<<" "
-                                               //<<" v0="<<value[0] <<" v1="<<value[1] <<" v2="<<value[2] <<" v3="<<value[3]
-                                               //<<" v4="<<value[4] <<" v5="<<value[5] <<" v6="<<value[6] <<" v7="<<value[7] );
-
-                               // check intersections of isosurface with edges, and calculate cubie index
-                               cubeIndex = 0;
-                               if (value[0] < mIsoValue) cubeIndex |= 1;
-                               if (value[1] < mIsoValue) cubeIndex |= 2;
-                               if (value[2] < mIsoValue) cubeIndex |= 4;
-                               if (value[3] < mIsoValue) cubeIndex |= 8;
-                               if (value[4] < mIsoValue) cubeIndex |= 16;
-                               if (value[5] < mIsoValue) cubeIndex |= 32;
-                               if (value[6] < mIsoValue) cubeIndex |= 64;
-                               if (value[7] < mIsoValue) cubeIndex |= 128;
-
-                               // No triangles to generate?
-                               if (mcEdgeTable[cubeIndex] == 0) {
-                                       continue;
-                               }
-
-                               // where to look up if this point already exists
-                               eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
-                               eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
-                               eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, k+0) ];
-                               eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
-
-                               eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, k+1) ];
-                               eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, k+1) ];
-                               eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, k+1) ];
-                               eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, k+1) ];
-
-                               eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
-                               eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, k+0) ];
-                               eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, k+0) ];
-                               eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, k+0) ];
-
-                               // grid positions
-                               pos[0] = ntlVec3Gfx(px    ,py    ,pz);
-                               pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
-                               pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
-                               pos[3] = ntlVec3Gfx(px    ,py+gsy,pz);
-                               pos[4] = ntlVec3Gfx(px    ,py    ,pz+gsz);
-                               pos[5] = ntlVec3Gfx(px+gsx,py    ,pz+gsz);
-                               pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz);
-                               pos[7] = ntlVec3Gfx(px    ,py+gsy,pz+gsz);
-
-                               // check all edges
-                               for(int e=0;e<12;e++) {
-
-                                       if (mcEdgeTable[cubeIndex] & (1<<e)) {
-                                               // is the vertex already calculated?
-                                               if(*eVert[ e ] < 0) {
-                                                       // interpolate edge
-                                                       const int e1 = mcEdges[e*2  ];
-                                                       const int e2 = mcEdges[e*2+1];
-                                                       const ntlVec3Gfx p1 = pos[ e1  ];    // scalar field pos 1
-                                                       const ntlVec3Gfx p2 = pos[ e2  ];    // scalar field pos 2
-                                                       const float valp1  = value[ e1  ];  // scalar field val 1
-                                                       const float valp2  = value[ e2  ];  // scalar field val 2
-                                                       //double mu;            // interpolation value
-                                                       //ntlVec3Gfx p;         // new point
-
-                                                       // choose if point should be calculated by interpolation,
-                                                       // or "Karolin" method 
-
-                                                       //double deltaVal = ABS(valp2-valp1);
-                                                       //if(deltaVal >-10.0) {
-                                                       // standard interpolation
-                                                       //vertInfo[e].type = 0; 
-                                                       /*if (ABS(mIsoValue-valp1) < 0.000000001) {
-                                                               mu = 0.0;
-                                                       } else {
-                                                               if (ABS(mIsoValue-valp2) < 0.000000001) {
-                                                                       mu = 1.0;
-                                                               } else {
-                                                                       mu = (mIsoValue - valp1) / (valp2 - valp1);
-                                                               }
-                                                       } */
-                                                       /*} else {
-                                                               errorOut(" ? ");
-                                                       // use fill grade (=karo) 
-                                                       vertInfo[e].type = 1; 
-                                                       if(valp1 < valp2) { mu = 1.0- (valp1 + valp2 - mIsoValue);
-                                                       } else { mu = 0.0+ (valp1 + valp2 - mIsoValue); } 
-                                                       } */
-
-                                                       //const float mu = (mIsoValue - valp1) / (valp2 - valp1);
-                                                       float mu;
-                                                       if(valp1 < valp2) {
-                                                               mu = 1.0 - 1.0*(valp1 + valp2 - mIsoValue);
-                                                       } else {
-                                                               mu = 0.0 + 1.0*(valp1 + valp2 - mIsoValue);
+       if(mSubdivs<=1) {
+
+               pz = mStart[2]-gsz*0.5;
+               for(int k=1;k<(mSizez-2);k++) {
+                       pz += gsz;
+                       py = mStart[1]-gsy*0.5;
+                       for(int j=1;j<(mSizey-2);j++) {
+                               py += gsy;
+                               px = mStart[0]-gsx*0.5;
+                               for(int i=1;i<(mSizex-2);i++) {
+                                       px += gsx;
+                                       int cubeIndex;      // index entry of the cube 
+                                       float value[8];
+                                       int triIndices[12]; // vertex indices 
+                                       int *eVert[12];
+                                       IsoLevelVertex ilv;
+                                       
+                                       value[0] = *getData(i  ,j  ,k  );
+                                       value[1] = *getData(i+1,j  ,k  );
+                                       value[2] = *getData(i+1,j+1,k  );
+                                       value[3] = *getData(i  ,j+1,k  );
+                                       value[4] = *getData(i  ,j  ,k+1);
+                                       value[5] = *getData(i+1,j  ,k+1);
+                                       value[6] = *getData(i+1,j+1,k+1);
+                                       value[7] = *getData(i  ,j+1,k+1);
+
+                                       // check intersections of isosurface with edges, and calculate cubie index
+                                       cubeIndex = 0;
+                                       if (value[0] < mIsoValue) cubeIndex |= 1;
+                                       if (value[1] < mIsoValue) cubeIndex |= 2;
+                                       if (value[2] < mIsoValue) cubeIndex |= 4;
+                                       if (value[3] < mIsoValue) cubeIndex |= 8;
+                                       if (value[4] < mIsoValue) cubeIndex |= 16;
+                                       if (value[5] < mIsoValue) cubeIndex |= 32;
+                                       if (value[6] < mIsoValue) cubeIndex |= 64;
+                                       if (value[7] < mIsoValue) cubeIndex |= 128;
+
+                                       // No triangles to generate?
+                                       if (mcEdgeTable[cubeIndex] == 0) {
+                                               continue;
+                                       }
+
+                                       // where to look up if this point already exists
+                                       int edgek = 0;
+                                       if(mUseFullEdgeArrays) edgek=k;
+                                       const int baseIn = ISOLEVEL_INDEX( i+0, j+0, edgek+0);
+                                       eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
+                                       eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
+                                       eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
+                                       eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
+
+                                       eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
+                                       eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, edgek+1) ];
+                                       eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+1) ];
+                                       eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
+
+                                       eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
+                                       eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, edgek+0) ];
+                                       eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, edgek+0) ];
+                                       eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
+
+                                       // grid positions
+                                       ntlVec3Gfx pos[8];
+                                       pos[0] = ntlVec3Gfx(px    ,py    ,pz);
+                                       pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
+                                       pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
+                                       pos[3] = ntlVec3Gfx(px    ,py+gsy,pz);
+                                       pos[4] = ntlVec3Gfx(px    ,py    ,pz+gsz);
+                                       pos[5] = ntlVec3Gfx(px+gsx,py    ,pz+gsz);
+                                       pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz);
+                                       pos[7] = ntlVec3Gfx(px    ,py+gsy,pz+gsz);
+
+                                       // check all edges
+                                       for(int e=0;e<12;e++) {
+                                               if (mcEdgeTable[cubeIndex] & (1<<e)) {
+                                                       // is the vertex already calculated?
+                                                       if(*eVert[ e ] < 0) {
+                                                               // interpolate edge
+                                                               const int e1 = mcEdges[e*2  ];
+                                                               const int e2 = mcEdges[e*2+1];
+                                                               const ntlVec3Gfx p1 = pos[ e1  ];    // scalar field pos 1
+                                                               const ntlVec3Gfx p2 = pos[ e2  ];    // scalar field pos 2
+                                                               const float valp1  = value[ e1  ];  // scalar field val 1
+                                                               const float valp2  = value[ e2  ];  // scalar field val 2
+                                                               const float mu = (mIsoValue - valp1) / (valp2 - valp1);
+
+                                                               // init isolevel vertex
+                                                               ilv.v = p1 + (p2-p1)*mu;
+                                                               ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) +
+                                                                                               getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * (    mu) ;
+                                                               mPoints.push_back( ilv );
+
+                                                               triIndices[e] = (mPoints.size()-1);
+                                                               // store vertex 
+                                                               *eVert[ e ] = triIndices[e];
+                                                       }       else {
+                                                               // retrieve  from vert array
+                                                               triIndices[e] = *eVert[ e ];
                                                        }
-
-                                                       float isov2 = mIsoValue;
-                                                       isov2 = (valp1+valp2)*0.5;
-                                                       mu = (isov2 - valp1) / (valp2 - valp1);
-                                                       mu = (isov2) / (valp2 - valp1);
-
-                                                       mu = (mIsoValue - valp1) / (valp2 - valp1);
-                                                       //mu *= mu;
-
-
-
-                                                       // init isolevel vertex
-                                                       ilv.v = p1 + (p2-p1)*mu;
-                                                       ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) +
-                                                                                       getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * (    mu) ;
-                                                       mPoints.push_back( ilv );
-
-                                                       triIndices[e] = (mPoints.size()-1);
-                                                       // store vertex 
-                                                       *eVert[ e ] = triIndices[e];
-                                               }       else {
-                                                       // retrieve  from vert array
-                                                       triIndices[e] = *eVert[ e ];
+                                               } // along all edges 
+                                       }
+
+                                       if( (i<coAdd+mCutoff) || (j<coAdd+mCutoff) ||
+                                                       ((mCutoff>0) && (k<coAdd)) ||// bottom layer
+                                                       (i>mSizex-2-coAdd-mCutoff) ||
+                                                       (j>mSizey-2-coAdd-mCutoff) ) {
+                                               if(mCutArray) {
+                                                       if(k < mCutArray[j*this->mSizex+i]) continue;
+                                               } else { continue; }
+                                       }
+
+                                       // Create the triangles... 
+                                       for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
+                                               mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
+                                               mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] );
+                                               mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
+                                       }
+                                       
+                               }//i
+                       }// j
+
+                       // copy edge arrays
+                       if(!mUseFullEdgeArrays) {
+                       for(int j=0;j<(mSizey-0);j++) 
+                               for(int i=0;i<(mSizex-0);i++) {
+                                       //int edgek = 0;
+                                       const int dst = ISOLEVEL_INDEX( i+0, j+0, 0);
+                                       const int src = ISOLEVEL_INDEX( i+0, j+0, 1);
+                                       mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ];
+                                       mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ];
+                                       mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ];
+                                       mpEdgeVerticesX[ src ]=-1;
+                                       mpEdgeVerticesY[ src ]=-1;
+                                       mpEdgeVerticesZ[ src ]=-1;
+                               }
+                       } // */
+
+               } // k
+
+       // precalculate normals using an approximation of the scalar field gradient 
+               for(int ni=0;ni<(int)mPoints.size();ni++) { normalize( mPoints[ni].n ); }
+
+       } else { // subdivs
+
+#define EDGEAR_INDEX(Ai,Aj,Ak, Bi,Bj) ((mSizex*mSizey*mSubdivs*mSubdivs*(Ak))+\
+               (mSizex*mSubdivs*((Aj)*mSubdivs+(Bj)))+((Ai)*mSubdivs)+(Bi))
+
+#define ISOTRILININT(fi,fj,fk) ( \
+                               (1.-(fi))*(1.-(fj))*(1.-(fk))*orgval[0] + \
+                               (   (fi))*(1.-(fj))*(1.-(fk))*orgval[1] + \
+                               (   (fi))*(   (fj))*(1.-(fk))*orgval[2] + \
+                               (1.-(fi))*(   (fj))*(1.-(fk))*orgval[3] + \
+                               (1.-(fi))*(1.-(fj))*(   (fk))*orgval[4] + \
+                               (   (fi))*(1.-(fj))*(   (fk))*orgval[5] + \
+                               (   (fi))*(   (fj))*(   (fk))*orgval[6] + \
+                               (1.-(fi))*(   (fj))*(   (fk))*orgval[7] )
+
+               // use subdivisions
+               gfxReal subdfac = 1./(gfxReal)(mSubdivs);
+               gfxReal orgGsx = gsx;
+               gfxReal orgGsy = gsy;
+               gfxReal orgGsz = gsz;
+               gsx *= subdfac;
+               gsy *= subdfac;
+               gsz *= subdfac;
+               if(mUseFullEdgeArrays) {
+                       errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!");
+               }
+               
+               ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex];
+
+               // construct pointers
+               // part test
+               int pInUse = 0;
+               int pUsedTest = 0;
+               // reset particles
+               // reset list array
+               for(int k=0;k<(mSizez);k++) 
+                       for(int j=0;j<(mSizey);j++) 
+                               for(int i=0;i<(mSizex);i++) {
+                                       arppnt[ISOLEVEL_INDEX(i,j,k)] = NULL;
+                               }
+               if(mpIsoParts) {
+                       for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
+                                       pit!= mpIsoParts->getParticlesEnd(); pit++) {
+                               if( (*pit).getActive()==false ) continue;
+                               if( (*pit).getType()!=PART_DROP) continue;
+                               (*pit).setNext(NULL);
+                       }
+                       // build per node lists
+                       for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
+                                       pit!= mpIsoParts->getParticlesEnd(); pit++) {
+                               if( (*pit).getActive()==false ) continue;
+                               if( (*pit).getType()!=PART_DROP) continue;
+                               // check lifetime ignored here
+                               ParticleObject *p = &(*pit);
+                               const ntlVec3Gfx ppos = p->getPos();
+                               const int pi= (int)round(ppos[0])+0; 
+                               const int pj= (int)round(ppos[1])+0; 
+                               int       pk= (int)round(ppos[2])+0;// no offset necessary
+                               // 2d should be handled by solver. if(LBMDIM==2) { pk = 0; }
+
+                               if(pi<0) continue;
+                               if(pj<0) continue;
+                               if(pk<0) continue;
+                               if(pi>mSizex-1) continue;
+                               if(pj>mSizey-1) continue;
+                               if(pk>mSizez-1) continue;
+                               ParticleObject* &pnt = arppnt[ISOLEVEL_INDEX(pi,pj,pk)]; 
+                               if(pnt) {
+                                       // append
+                                       ParticleObject* listpnt = pnt;
+                                       while(listpnt) {
+                                               if(!listpnt->getNext()) {
+                                                       listpnt->setNext(p); listpnt = NULL;
+                                               } else {
+                                                       listpnt = listpnt->getNext();
                                                }
-                                       } // along all edges 
-
+                                       }
+                               } else {
+                                       // start new list
+                                       pnt = p;
                                }
+                               pInUse++;
+                       }
+               } // mpIsoParts
 
+               debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<<pInUse<<", Subdivs:"<<mSubdivs, 9);
+               pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz;
+       
+               for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) {
+                       pz += gsz;
+                       const int k = ok/mSubdivs;
+                       if(k<=0) continue; // skip zero plane
+#if PARALLEL==1        
+#pragma omp parallel for
+#endif 
+                       for(int j=1;j<(mSizey-2);j++) {
+                               for(int i=1;i<(mSizex-2);i++) {
+                                       float value[8];
+                                       ntlVec3Gfx pos[8];
+                                       int cubeIndex;      // index entry of the cube 
+                                       int triIndices[12]; // vertex indices 
+                                       int *eVert[12];
+                                       IsoLevelVertex ilv;
+                                       gfxReal orgval[8];
+                                       gfxReal subdAr[2][11][11]; // max 10 subdivs!
+                                       
+                                       orgval[0] = *getData(i  ,j  ,k  );
+                                       orgval[1] = *getData(i+1,j  ,k  );
+                                       orgval[2] = *getData(i+1,j+1,k  ); // with subdivs
+                                       orgval[3] = *getData(i  ,j+1,k  );
+                                       orgval[4] = *getData(i  ,j  ,k+1);
+                                       orgval[5] = *getData(i+1,j  ,k+1);
+                                       orgval[6] = *getData(i+1,j+1,k+1); // with subdivs
+                                       orgval[7] = *getData(i  ,j+1,k+1);
+
+                                       // prebuild subsampled array slice
+                                       const int sdkOffset = ok-k*mSubdivs; 
+
+                                       for(int sdk=0; sdk<2; sdk++) 
+                                               for(int sdj=0; sdj<mSubdivs+1; sdj++) 
+                                                       for(int sdi=0; sdi<mSubdivs+1; sdi++) {
+                                                               subdAr[sdk][sdj][sdi] = ISOTRILININT(sdi*subdfac, sdj*subdfac, (sdkOffset+sdk)*subdfac);
+                                                       }
 
-                               // Create the triangles... 
-                               for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
-                                       //errMsg("MC","tri "<<mIndices.size() <<" "<< triIndices[ mcTriTable[cubeIndex][e+0] ]<<" "<< triIndices[ mcTriTable[cubeIndex][e+1] ]<<" "<< triIndices[ mcTriTable[cubeIndex][e+2] ] );
-                                       mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
-                                       mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] );
-                                       mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
+                                       const int poDistOffset=2;
+                                       for(int pok=-poDistOffset; pok<1+poDistOffset; pok++) {
+                                               if(k+pok<0) continue;
+                                               if(k+pok>=mSizez-1) continue;
+                                       for(int poj=-poDistOffset; poj<1+poDistOffset; poj++) {
+                                               if(j+poj<0) continue;
+                                               if(j+poj>=mSizey-1) continue;
+                                       for(int poi=-poDistOffset; poi<1+poDistOffset; poi++) {
+                                               if(i+poi<0) continue;
+                                               if(i+poi>=mSizex-1) continue; 
+                                               ParticleObject *p;
+                                               p = arppnt[ISOLEVEL_INDEX(i+poi,j+poj,k+pok)];
+                                               while(p) { // */
+                                       /*
+                                       for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
+                                                       pit!= mpIsoParts->getParticlesEnd(); pit++) { { { {
+                                               // debug test! , full list slow!
+                                               if(( (*pit).getActive()==false ) || ( (*pit).getType()!=PART_DROP)) continue;
+                                               ParticleObject *p;
+                                               p = &(*pit); // */
+
+                                                       pUsedTest++;
+                                                       ntlVec3Gfx ppos = p->getPos();
+                                                       const int spi= (int)round( (ppos[0]+1.-(gfxReal)i) *(gfxReal)mSubdivs-1.5); 
+                                                       const int spj= (int)round( (ppos[1]+1.-(gfxReal)j) *(gfxReal)mSubdivs-1.5); 
+                                                       const int spk= (int)round( (ppos[2]+1.-(gfxReal)k) *(gfxReal)mSubdivs-1.5)-sdkOffset; // why -2?
+                                                       // 2d should be handled by solver. if(LBMDIM==2) { spk = 0; }
+
+                                                       gfxReal pfLen = p->getSize()*1.5*mPartSize;  // test, was 1.1
+                                                       const gfxReal minPfLen = subdfac*0.8;
+                                                       if(pfLen<minPfLen) pfLen = minPfLen;
+                                                       //errMsg("ISOPPP"," at "<<PRINT_IJK<<"  pp"<<ppos<<"  sp"<<PRINT_VEC(spi,spj,spk)<<" pflen"<<pfLen );
+                                                       //errMsg("ISOPPP"," subdfac="<<subdfac<<" size"<<p->getSize()<<" ps"<<mPartSize );
+                                                       const int icellpsize = (int)(1.*pfLen*(gfxReal)mSubdivs)+1;
+                                                       for(int swk=-icellpsize; swk<=icellpsize; swk++) {
+                                                               if(spk+swk<         0) { continue; }
+                                                               if(spk+swk>         1) { continue; } // */
+                                                       for(int swj=-icellpsize; swj<=icellpsize; swj++) {
+                                                               if(spj+swj<         0) { continue; }
+                                                               if(spj+swj>mSubdivs+0) { continue; } // */
+                                                       for(int swi=-icellpsize; swi<=icellpsize; swi++) {
+                                                               if(spi+swi<         0) { continue; } 
+                                                               if(spi+swi>mSubdivs+0) { continue; } // */
+                                                               ntlVec3Gfx cellp = ntlVec3Gfx(
+                                                                               (1.5+(gfxReal)(spi+swi))           *subdfac + (gfxReal)(i-1),
+                                                                               (1.5+(gfxReal)(spj+swj))           *subdfac + (gfxReal)(j-1),
+                                                                               (1.5+(gfxReal)(spk+swk)+sdkOffset) *subdfac + (gfxReal)(k-1)
+                                                                               );
+                                                               //if(swi==0 && swj==0 && swk==0) subdAr[spk][spj][spi] = 1.; // DEBUG
+                                                               // clip domain boundaries again 
+                                                               if(cellp[0]<1.) { continue; } 
+                                                               if(cellp[1]<1.) { continue; } 
+                                                               if(cellp[2]<1.) { continue; } 
+                                                               if(cellp[0]>(gfxReal)mSizex-3.) { continue; } 
+                                                               if(cellp[1]>(gfxReal)mSizey-3.) { continue; } 
+                                                               if(cellp[2]>(gfxReal)mSizez-3.) { continue; } 
+                                                               gfxReal len = norm(cellp-ppos);
+                                                               gfxReal isoadd = 0.; 
+                                                               const gfxReal baseIsoVal = mIsoValue*1.1;
+                                                               if(len<pfLen) { 
+                                                                       isoadd = baseIsoVal*1.;
+                                                               } else { 
+                                                                       // falloff linear with pfLen (kernel size=2pfLen
+                                                                       isoadd = baseIsoVal*(1. - (len-pfLen)/(pfLen)); 
+                                                               }
+                                                               if(isoadd<0.) { continue; }
+                                                               //errMsg("ISOPPP"," at "<<PRINT_IJK<<" sp"<<PRINT_VEC(spi+swi,spj+swj,spk+swk)<<" cellp"<<cellp<<" pp"<<ppos << " l"<< len<< " add"<< isoadd);
+                                                               const gfxReal arval = subdAr[spk+swk][spj+swj][spi+swi];
+                                                               if(arval>1.) { continue; }
+                                                               subdAr[spk+swk][spj+swj][spi+swi] = arval + isoadd;
+                                                       } } }
+
+                                                       p = p->getNext();
+                                               }
+                                       } } } // poDist loops */
+
+                                       py = mStart[1]+(((double)j-0.5)*orgGsy)-gsy;
+                                       for(int sj=0;sj<mSubdivs;sj++) {
+                                               py += gsy;
+                                               px = mStart[0]+(((double)i-0.5)*orgGsx)-gsx;
+                                               for(int si=0;si<mSubdivs;si++) {
+                                                       px += gsx;
+                                                       value[0] = subdAr[0+0][sj+0][si+0]; 
+                                                       value[1] = subdAr[0+0][sj+0][si+1]; 
+                                                       value[2] = subdAr[0+0][sj+1][si+1]; 
+                                                       value[3] = subdAr[0+0][sj+1][si+0]; 
+                                                       value[4] = subdAr[0+1][sj+0][si+0]; 
+                                                       value[5] = subdAr[0+1][sj+0][si+1]; 
+                                                       value[6] = subdAr[0+1][sj+1][si+1]; 
+                                                       value[7] = subdAr[0+1][sj+1][si+0]; 
+
+                                                       // check intersections of isosurface with edges, and calculate cubie index
+                                                       cubeIndex = 0;
+                                                       if (value[0] < mIsoValue) cubeIndex |= 1;
+                                                       if (value[1] < mIsoValue) cubeIndex |= 2; // with subdivs
+                                                       if (value[2] < mIsoValue) cubeIndex |= 4;
+                                                       if (value[3] < mIsoValue) cubeIndex |= 8;
+                                                       if (value[4] < mIsoValue) cubeIndex |= 16;
+                                                       if (value[5] < mIsoValue) cubeIndex |= 32; // with subdivs
+                                                       if (value[6] < mIsoValue) cubeIndex |= 64;
+                                                       if (value[7] < mIsoValue) cubeIndex |= 128;
+
+                                                       if (mcEdgeTable[cubeIndex] >  0) {
+
+                                                       // where to look up if this point already exists
+                                                       const int edgek = 0;
+                                                       const int baseIn = EDGEAR_INDEX( i+0, j+0, edgek+0, si,sj);
+                                                       eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
+                                                       eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
+                                                       eVert[ 2] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ];
+                                                       eVert[ 3] = &mpEdgeVerticesY[ baseIn ];                             
+                                                                                                                                                                                                                                                                                                                                       
+                                                       eVert[ 4] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ];
+                                                       eVert[ 5] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+1,sj+0) ]; // with subdivs
+                                                       eVert[ 6] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+1) ];
+                                                       eVert[ 7] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ];
+                                                                                                                                                                                                                                                                                                                                       
+                                                       eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];                             
+                                                       eVert[ 9] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+0) ]; // with subdivs
+                                                       eVert[10] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+1) ];
+                                                       eVert[11] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ];
+
+                                                       // grid positions
+                                                       pos[0] = ntlVec3Gfx(px    ,py    ,pz);
+                                                       pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
+                                                       pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz); // with subdivs
+                                                       pos[3] = ntlVec3Gfx(px    ,py+gsy,pz);
+                                                       pos[4] = ntlVec3Gfx(px    ,py    ,pz+gsz);
+                                                       pos[5] = ntlVec3Gfx(px+gsx,py    ,pz+gsz);
+                                                       pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz); // with subdivs
+                                                       pos[7] = ntlVec3Gfx(px    ,py+gsy,pz+gsz);
+
+                                                       // check all edges
+                                                       for(int e=0;e<12;e++) {
+                                                               if (mcEdgeTable[cubeIndex] & (1<<e)) {
+                                                                       // is the vertex already calculated?
+                                                                       if(*eVert[ e ] < 0) {
+                                                                               // interpolate edge
+                                                                               const int e1 = mcEdges[e*2  ];
+                                                                               const int e2 = mcEdges[e*2+1];
+                                                                               const ntlVec3Gfx p1 = pos[ e1  ];   // scalar field pos 1
+                                                                               const ntlVec3Gfx p2 = pos[ e2  ];   // scalar field pos 2
+                                                                               const float valp1  = value[ e1  ];  // scalar field val 1
+                                                                               const float valp2  = value[ e2  ];  // scalar field val 2
+                                                                               const float mu = (mIsoValue - valp1) / (valp2 - valp1);
+
+                                                                               // init isolevel vertex
+                                                                               ilv.v = p1 + (p2-p1)*mu; // with subdivs
+#if PARALLEL==1        
+#pragma omp critical
+#endif
+                                                                               {
+                                                                               mPoints.push_back( ilv );
+                                                                               triIndices[e] = (mPoints.size()-1);
+                                                                               }
+                                                                               // store vertex 
+                                                                               *eVert[ e ] = triIndices[e]; 
+                                                                       }       else {
+                                                                               // retrieve  from vert array
+                                                                               triIndices[e] = *eVert[ e ];
+                                                                       }
+                                                               } // along all edges 
+                                                       }
+                                                       // removed cutoff treatment...
+                                                       
+                                                       // Create the triangles... 
+#if PARALLEL==1        
+#pragma omp critical
+#endif
+                                                       {
+                                                       for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
+                                                               mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
+                                                               mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); // with subdivs
+                                                               mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
+                                                               //errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() );
+                                                       } // triangles in edge table?
+                                                       }
+                                                       }
+                                                       
+                                               }//si
+                                       }// sj
+
+                               }//i
+                       }// j
+                       
+                       // copy edge arrays
+                       for(int j=0;j<(mSizey-0)*mSubdivs;j++) 
+                               for(int i=0;i<(mSizex-0)*mSubdivs;i++) {
+                                       //int edgek = 0;
+                                       const int dst = EDGEAR_INDEX( 0, 0, 0, i,j);
+                                       const int src = EDGEAR_INDEX( 0, 0, 1, i,j);
+                                       mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ];
+                                       mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; // with subdivs
+                                       mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ];
+                                       mpEdgeVerticesX[ src ]=-1;
+                                       mpEdgeVerticesY[ src ]=-1; // with subdivs
+                                       mpEdgeVerticesZ[ src ]=-1;
                                }
+                       // */
 
-                               
-      }
-    }
+               } // ok, k subdiv loop
 
-  } // k
-  
-
-  // precalculate normals using an approximation of the scalar field gradient 
-       for(int ni=0;ni<(int)mPoints.size();ni++) {
-               // use triangle normals?, this seems better for File-IsoSurf
-               normalize( mPoints[ni].n );
-       }
+               //delete [] subdAr;
+               delete [] arppnt;
+               computeNormals();
+       } // with subdivs
 
-       if(mSmoothSurface>0.0) {
-               // not needed for post normal smoothing? 
-               // if(mSmoothNormals<=0.0) { smoothNormals(mSmoothSurface*0.5); }
-               smoothSurface(mSmoothSurface);
+       // perform smoothing
+       float smoSubdfac = 1.;
+       if(mSubdivs>0) {
+               //smoSubdfac = 1./(float)(mSubdivs);
+               smoSubdfac = pow(0.55,(double)mSubdivs); // slightly stronger
        }
-       for(int i=0; i<mLoopSubdivs; i++) {
-               subdivide();
+       if(mSmoothSurface>0. || mSmoothNormals>0.) debMsgStd("IsoSurface::triangulate",DM_MSG,"Smoothing...",10);
+       if(mSmoothSurface>0.0) { 
+               smoothSurface(mSmoothSurface*smoSubdfac, (mSmoothNormals<=0.0) ); 
        }
-       if(mSmoothNormals>0.0) {
-               smoothNormals(mSmoothNormals);
+       if(mSmoothNormals>0.0) { 
+               smoothNormals(mSmoothNormals*smoSubdfac); 
        }
+
        myTime_t tritimeend = getTime(); 
-       debMsgStd("IsoSurface::triangulate",DM_MSG,"Took "<< ((tritimeend-tritimestart)/(double)1000.0)<<"s) " , 10 );
+       debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< getTimeString(tritimeend-tritimestart)<<", S("<<mSmoothSurface<<","<<mSmoothNormals<<"),"<<
+                       " verts:"<<mPoints.size()<<" tris:"<<(mIndices.size()/3)<<" subdivs:"<<mSubdivs
+                , 10 );
+       if(mpIsoParts) debMsgStd("IsoSurface::triangulate",DM_MSG,"parts:"<<mpIsoParts->getNumParticles(), 10);
 }
 
 
@@ -342,7 +681,7 @@ void IsoSurface::triangulate( void )
 /******************************************************************************
  * Get triangles for rendering
  *****************************************************************************/
-void IsoSurface::getTriangles( vector<ntlTriangle> *triangles, 
+void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles, 
                                                                                                         vector<ntlVec3Gfx> *vertices, 
                                                                                                         vector<ntlVec3Gfx> *normals, int objectId )
 {
@@ -350,6 +689,8 @@ void IsoSurface::getTriangles( vector<ntlTriangle> *triangles,
                debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10);
                return;
        }
+       t = 0.;
+       //return; // DEBUG
 
   /* triangulate field */
   triangulate();
@@ -393,12 +734,6 @@ void IsoSurface::getTriangles( vector<ntlTriangle> *triangles,
                if(getVisible()){ flag |= TRI_GEOMETRY; }
                if(getCastShadows() ) { 
                        flag |= TRI_CASTSHADOWS; } 
-               if( (getMaterial()->getMirror()>0.0) ||  
-                               (getMaterial()->getTransparence()>0.0) ||  
-                               (getMaterial()->getFresnel()>0.0) ) { 
-                       flag |= TRI_MAKECAUSTICS; } 
-               else { 
-                       flag |= TRI_NOCAUSTICS; } 
 
                /* init geo init id */
                int geoiId = getGeoInitId(); 
@@ -422,15 +757,7 @@ void IsoSurface::getTriangles( vector<ntlTriangle> *triangles,
 
 
 inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) {
-       // WARNING - this assumes a security boundary layer... 
-       /*
-       if(i<=0) i=1;
-       if(j<=0) j=1;
-       if(k<=0) k=1;
-       if(i>=(mSizex-1)) i=mSizex-2;
-       if(j>=(mSizex-1)) j=mSizex-2;
-       if(k>=(mSizex-1)) k=mSizex-2; // */
-
+       // WARNING - this requires a security boundary layer... 
        ntlVec3Gfx ret(0.0);
        ret[0] = *getData(i-1,j  ,k  ) - 
                 *getData(i+1,j  ,k  );
@@ -441,352 +768,120 @@ inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) {
        return ret;
 }
 
-#define RECALCNORMALS 0
 
-// Subdivide a mesh allways loop
-void IsoSurface::subdivide()
-{
-       mAdjacentFaces.clear(); 
-       mAcrossEdge.clear();
-
-       //void TriMesh::need_adjacentfaces()
-       {
-               vector<int> numadjacentfaces(mPoints.size());
-               //errMsg("SUBDIV ADJFA1", " "<<mPoints.size()<<" - "<<numadjacentfaces.size() );
-               for (int i = 0; i < (int)mIndices.size()/3; i++) {
-                       numadjacentfaces[mIndices[i*3 + 0]]++;
-                       numadjacentfaces[mIndices[i*3 + 1]]++;
-                       numadjacentfaces[mIndices[i*3 + 2]]++;
-               }
 
-               mAdjacentFaces.resize(mPoints.size());
-               for (int i = 0; i < (int)mPoints.size(); i++)
-                       mAdjacentFaces[i].reserve(numadjacentfaces[i]);
 
-               for (int i = 0; i < (int)mIndices.size()/3; i++) {
-                       for (int j = 0; j < 3; j++)
-                               mAdjacentFaces[mIndices[i*3 + j]].push_back(i);
-               }
-
-       }
-
-       // Find the face across each edge from each other face (-1 on boundary)
-       // If topology is bad, not necessarily what one would expect...
-       //void TriMesh::need_across_edge()
-       {
-               mAcrossEdge.resize(mIndices.size(), -1);
-
-               for (int i = 0; i < (int)mIndices.size()/3; i++) {
-                       for (int j = 0; j < 3; j++) {
-                               if (mAcrossEdge[i*3 + j] != -1)
-                                       continue;
-                               int v1 = mIndices[i*3 + ((j+1)%3)];
-                               int v2 = mIndices[i*3 + ((j+2)%3)];
-                               const vector<int> &a1 = mAdjacentFaces[v1];
-                               const vector<int> &a2 = mAdjacentFaces[v2];
-                               for (int k1 = 0; k1 < (int)a1.size(); k1++) {
-                                       int other = a1[k1];
-                                       if (other == i)
-                                               continue;
-                                       vector<int>::const_iterator it =
-                                               std::find(a2.begin(), a2.end(), other);
-                                       if (it == a2.end())
-                                               continue;
-
-                                       //int ind = (faces[other].indexof(v1)+1)%3;
-                                       int ind = -1;
-                                       if( mIndices[other*3+0] == (unsigned int)v1 ) ind = 0;
-                                       else if( mIndices[other*3+1] == (unsigned int)v1 ) ind = 1;
-                                       else if( mIndices[other*3+2] == (unsigned int)v1 ) ind = 2;
-                                       ind = (ind+1)%3;
-
-                                       if ( (int)mIndices[other*3 + ((ind+1)%3)] != v2)
-                                               continue;
-                                       mAcrossEdge[i*3 + j] = other;
-                                       mAcrossEdge[other*3 + ind] = i;
-                                       break;
-                               }
-                       }
-               }
-
-               //errMsg("SUBDIV ACREDG", "Done.\n");
-       }
-
-       //errMsg("SUBDIV","start");
-       // Introduce new vertices
-       int nf = (int)mIndices.size() / 3;
-
-       //vector<TriMesh::Face> newverts(nf, TriMesh::Face(-1,-1,-1));
-       vector<int> newverts(nf*3); //, TriMesh::Face(-1,-1,-1));
-       for(int j=0; j<(int)newverts.size(); j++) newverts[j] = -1;
-
-       int old_nv = (int)mPoints.size();
-       mPoints.reserve(4 * old_nv);
-       vector<int> newvert_count(old_nv + 3*nf); // wichtig...?
-       //errMsg("NC", newvert_count.size() );
-
-       for (int i = 0; i < nf; i++) {
-               for (int j = 0; j < 3; j++) {
-                       int ae = mAcrossEdge[i*3 + j];
-                       if (newverts[i*3 + j] == -1 && ae != -1) {
-                               if (mAcrossEdge[ae*3 + 0] == i)
-                                       newverts[i*3 + j] = newverts[ae*3 + 0];
-                               else if (mAcrossEdge[ae*3 + 1] == i)
-                                       newverts[i*3 + j] = newverts[ae*3 + 1];
-                               else if (mAcrossEdge[ae*3 + 2] == i)
-                                       newverts[i*3 + j] = newverts[ae*3 + 2];
-                       }
-                       if (newverts[i*3 + j] == -1) {
-                               IsoLevelVertex ilv;
-                               ilv.v = ntlVec3Gfx(0.0);
-                               ilv.n = ntlVec3Gfx(0.0);
-                               mPoints.push_back(ilv);
-                               newverts[i*3 + j] = (int)mPoints.size() - 1;
-                               if (ae != -1) {
-                                       if (mAcrossEdge[ae*3 + 0] == i)
-                                               newverts[ae*3 + 0] = newverts[i*3 + j];
-                                       else if (mAcrossEdge[ae*3 + 1] == i)
-                                               newverts[ae*3 + 1] = newverts[i*3 + j];
-                                       else if (mAcrossEdge[ae*3 + 2] == i)
-                                               newverts[ae*3 + 2] = newverts[i*3 + j];
-                               }
-                       }
-                       if(ae != -1) {
-                               mPoints[newverts[i*3 + j]].v +=
-                                       mPoints[ mIndices[i*3 + ( j     )] ].v * 0.25f  + // j = 0,1,2?
-                                       mPoints[ mIndices[i*3 + ((j+1)%3)] ].v * 0.375f +
-                                       mPoints[ mIndices[i*3 + ((j+2)%3)] ].v * 0.375f;
-#if RECALCNORMALS==0
-                               mPoints[newverts[i*3 + j]].n +=
-                                       mPoints[ mIndices[i*3 + ( j     )] ].n * 0.25f  + // j = 0,1,2?
-                                       mPoints[ mIndices[i*3 + ((j+1)%3)] ].n * 0.375f +
-                                       mPoints[ mIndices[i*3 + ((j+2)%3)] ].n * 0.375f;
-#endif // RECALCNORMALS==0
-                       } else {
-                               mPoints[newverts[i*3 + j]].v +=
-                                       mPoints[ mIndices[i*3 + ((j+1)%3)] ].v * 0.5f   +
-                                       mPoints[ mIndices[i*3 + ((j+2)%3)] ].v * 0.5f  ;
-#if RECALCNORMALS==0
-                               mPoints[newverts[i*3 + j]].n +=
-                                       mPoints[ mIndices[i*3 + ((j+1)%3)] ].n * 0.5f   +
-                                       mPoints[ mIndices[i*3 + ((j+2)%3)] ].n * 0.5f  ;
-#endif // RECALCNORMALS==0
-                       }
-
-                       newvert_count[newverts[i*3 + j]]++;
-               }
-       }
-       for (int i = old_nv; i < (int)mPoints.size(); i++) {
-               if (!newvert_count[i])
-                       continue;
-               float scale = 1.0f / newvert_count[i];
-               mPoints[i].v *= scale;
-
-#if RECALCNORMALS==0
-               //mPoints[i].n *= scale;
-               //normalize( mPoints[i].n );
-#endif // RECALCNORMALS==0
-       }
-
-       // Update old vertices
-       for (int i = 0; i < old_nv; i++) {
-                       ntlVec3Gfx bdyavg(0.0), nbdyavg(0.0);
-                       ntlVec3Gfx norm_bdyavg(0.0), norm_nbdyavg(0.0); // N
-                       int nbdy = 0, nnbdy = 0;
-                       int naf = (int)mAdjacentFaces[i].size();
-                       if (!naf)
-                               continue;
-                       for (int j = 0; j < naf; j++) {
-                               int af = mAdjacentFaces[i][j];
-
-                               int afi = -1;
-                               if( mIndices[af*3+0] == (unsigned int)i ) afi = 0;
-                               else if( mIndices[af*3+1] == (unsigned int)i ) afi = 1;
-                               else if( mIndices[af*3+2] == (unsigned int)i ) afi = 2;
-
-                               int n1 = (afi+1) % 3;
-                               int n2 = (afi+2) % 3;
-                               if (mAcrossEdge[af*3 + n1] == -1) {
-                                       bdyavg += mPoints[newverts[af*3 + n1]].v;
-#if RECALCNORMALS==0
-                                       //norm_bdyavg += mPoints[newverts[af*3 + n1]].n;
-#endif // RECALCNORMALS==0
-                                       nbdy++;
-                               } else {
-                                       nbdyavg += mPoints[newverts[af*3 + n1]].v;
-#if RECALCNORMALS==0
-                                       //norm_nbdyavg += mPoints[newverts[af*3 + n1]].n;
-#endif // RECALCNORMALS==0
-                                       nnbdy++;
-                               }
-                               if (mAcrossEdge[af*3 + n2] == -1) {
-                                       bdyavg += mPoints[newverts[af*3 + n2]].v;
-#if RECALCNORMALS==0
-                                       //norm_bdyavg += mPoints[newverts[af*3 + n2]].n;
-#endif // RECALCNORMALS==0
-                                       nbdy++;
-                               } else {
-                                       nbdyavg += mPoints[newverts[af*3 + n2]].v;
-#if RECALCNORMALS==0
-                                       //norm_nbdyavg += mPoints[newverts[af*3 + n2]].n;
-#endif // RECALCNORMALS==0
-                                       nnbdy++;
-                               }
-                       }
+/******************************************************************************
+ * 
+ * Surface improvement, inspired by trimesh2 library
+ * (http://www.cs.princeton.edu/gfx/proj/trimesh2/)
+ * 
+ *****************************************************************************/
 
-                       float alpha;
-                       ntlVec3Gfx newpt;
-                       if (nbdy) {
-                               newpt = bdyavg / (float) nbdy;
-                               alpha = 0.5f;
-                       } else if (nnbdy) {
-                               newpt = nbdyavg / (float) nnbdy;
-                               if (nnbdy == 6)
-                                       alpha = 1.05;
-                               else if (nnbdy == 8)
-                                       alpha = 0.86;
-                               else if (nnbdy == 10)
-                                       alpha = 0.7;
-                               else
-                                       alpha = 0.6;
-                       } else {
-                               continue;
-                       }
-                       mPoints[i].v *= 1.0f - alpha;
-                       mPoints[i].v += newpt * alpha;
+void IsoSurface::setSmoothRad(float radi1, float radi2, ntlVec3Gfx mscc) {
+       mSCrad1 = radi1*radi1;
+       mSCrad2 = radi2*radi2;
+       mSCcenter = mscc;
+}
 
-#if RECALCNORMALS==0
-                       //mPoints[i].n *= 1.0f - alpha;
-                       //mPoints[i].n += newpt * alpha;
-#endif // RECALCNORMALS==0
+// compute normals for all generated triangles
+void IsoSurface::computeNormals() {
+  for(int i=0;i<(int)mPoints.size();i++) {
+               mPoints[i].n = ntlVec3Gfx(0.);
        }
 
-       // Insert new faces
-       mIndices.reserve(4*nf);
-       for (int i = 0; i < nf; i++) {
-               mIndices.push_back( mIndices[i*3 + 0]);
-               mIndices.push_back( newverts[i*3 + 2]);
-               mIndices.push_back( newverts[i*3 + 1]);
-
-               mIndices.push_back( mIndices[i*3 + 1]);
-               mIndices.push_back( newverts[i*3 + 0]);
-               mIndices.push_back( newverts[i*3 + 2]);
-
-               mIndices.push_back( mIndices[i*3 + 2]);
-               mIndices.push_back( newverts[i*3 + 1]);
-               mIndices.push_back( newverts[i*3 + 0]);
-
-               mIndices[i*3+0] = newverts[i*3+0];
-               mIndices[i*3+1] = newverts[i*3+1];
-               mIndices[i*3+2] = newverts[i*3+2];
+  for(int i=0;i<(int)mIndices.size();i+=3) {
+    const int t1 = mIndices[i];
+    const int t2 = mIndices[i+1];
+               const int t3 = mIndices[i+2];
+               const ntlVec3Gfx p1 = mPoints[t1].v;
+               const ntlVec3Gfx p2 = mPoints[t2].v;
+               const ntlVec3Gfx p3 = mPoints[t3].v;
+               const ntlVec3Gfx n1=p1-p2;
+               const ntlVec3Gfx n2=p2-p3;
+               const ntlVec3Gfx n3=p3-p1;
+               const gfxReal len1 = normNoSqrt(n1);
+               const gfxReal len2 = normNoSqrt(n2);
+               const gfxReal len3 = normNoSqrt(n3);
+               const ntlVec3Gfx norm = cross(n1,n2);
+               mPoints[t1].n += norm * (1./(len1*len3));
+               mPoints[t2].n += norm * (1./(len1*len2));
+               mPoints[t3].n += norm * (1./(len2*len3));
        }
 
-       // recalc normals
-#if RECALCNORMALS==1
-       {
-               int nf = (int)mIndices.size()/3, nv = (int)mPoints.size();
-               for (int i = 0; i < nv; i++) {
-                       mPoints[i].n = ntlVec3Gfx(0.0);
-               }
-               for (int i = 0; i < nf; i++) {
-                       const ntlVec3Gfx &p0 = mPoints[mIndices[i*3+0]].v;
-                       const ntlVec3Gfx &p1 = mPoints[mIndices[i*3+1]].v;
-                       const ntlVec3Gfx &p2 = mPoints[mIndices[i*3+2]].v;
-                       ntlVec3Gfx a = p0-p1, b = p1-p2, c = p2-p0;
-                       float l2a = normNoSqrt(a), l2b = normNoSqrt(b), l2c = normNoSqrt(c);
-
-                       ntlVec3Gfx facenormal = cross(a, b);
-
-                       mPoints[mIndices[i*3+0]].n += facenormal * (1.0f / (l2a * l2c));
-                       mPoints[mIndices[i*3+1]].n += facenormal * (1.0f / (l2b * l2a));
-                       mPoints[mIndices[i*3+2]].n += facenormal * (1.0f / (l2c * l2b));
-               }
-
-               for (int i = 0; i < nv; i++) {
-                       normalize(mPoints[i].n);
-               }
+  for(int i=0;i<(int)mPoints.size();i++) {
+               normalize(mPoints[i].n);
        }
-#else // RECALCNORMALS==1
-               for (int i = 0; i < (int)mPoints.size(); i++) {
-                       normalize(mPoints[i].n);
-               }
-#endif // RECALCNORMALS==1
-
-       //errMsg("SUBDIV","done nv:"<<mPoints.size()<<" nf:"<<mIndices.size() );
 }
 
-
-
-
-
-
-
-/******************************************************************************
- * 
- * Surface improvement
- * makes use of trimesh2 library
- * http://www.cs.princeton.edu/gfx/proj/trimesh2/
- *
- * Copyright (c) 2004 Szymon Rusinkiewicz.
- * see COPYING_trimesh2
- * 
- *****************************************************************************/
-
-
-
 // Diffuse a vector field at 1 vertex, weighted by
 // a gaussian of width 1/sqrt(invsigma2)
-void IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int v, float invsigma2, ntlVec3Gfx &flt)
+bool IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int src, float invsigma2, ntlVec3Gfx &target)
 {
-       flt = ntlVec3Gfx(0.0);
-       flt += *(field+pointerScale*v) *pointareas[v];
-       float sum_w = pointareas[v];
-       const ntlVec3Gfx &nv = mPoints[v].n;
-
-       unsigned &flag = flag_curr;
-       flag++;
-       flags[v] = flag;
-       vector<int> boundary = neighbors[v];
-       while (!boundary.empty()) {
-               const int bbn = boundary.back();
-               boundary.pop_back();
-               if (flags[bbn] == flag) continue;
+       if((neighbors[src].size()<1) || (pointareas[src]<=0.0)) return 0;
+       const ntlVec3Gfx srcp = mPoints[src].v;
+       const ntlVec3Gfx srcn = mPoints[src].n;
+       if(mSCrad1>0.0 && mSCrad2>0.0) {
+               ntlVec3Gfx dp = mSCcenter-srcp; dp[2] = 0.0; // only xy-plane
+               float rd = normNoSqrt(dp);
+               if(rd > mSCrad2) {
+                       return 0;
+               } else if(rd > mSCrad1) {
+                       // optimize?
+                       float org = 1.0/sqrt(invsigma2);
+                       org *= (1.0- (rd-mSCrad1) / (mSCrad2-mSCrad1));
+                       invsigma2 = 1.0/(org*org);
+                       //errMsg("TRi","p"<<srcp<<" rd:"<<rd<<" r1:"<<mSCrad1<<" r2:"<<mSCrad2<<" org:"<<org<<" is:"<<invsigma2);
+               } else {
+               }
+       }
+       target = ntlVec3Gfx(0.0);
+       target += *(field+pointerScale*src) *pointareas[src];
+       float smstrSum = pointareas[src];
+
+       int flag = mFlagCnt; 
+       mFlagCnt++;
+       flags[src] = flag;
+       mDboundary = neighbors[src];
+       while (!mDboundary.empty()) {
+               const int bbn = mDboundary.back();
+               mDboundary.pop_back();
+               if(flags[bbn]==flag) continue;
                flags[bbn] = flag;
 
+               // normal check
+               const float nvdot = dot(srcn, mPoints[bbn].n); // faster than before d2 calc?
+               if(nvdot <= 0.0f) continue;
+
                // gaussian weight of width 1/sqrt(invsigma2)
-               const float d2 = invsigma2 * normNoSqrt(mPoints[bbn].v - mPoints[v].v);
-               if(d2 >= 9.0f) continue; // 25 also possible  , slower
-               //float w = (d2 >=  9.0f) ? 0.0f : exp(-0.5f*d2);
-               float w = exp(-0.5f*d2);
-               if(dot(nv, mPoints[bbn].n) <= 0.0f) continue; // faster than before d2 calc?
-
-               // Downweight things pointing in different directions
-               w *= dot(nv , mPoints[bbn].n);
-               // Surface area "belonging" to each point
-               w *= pointareas[bbn];
+               const float d2 = invsigma2 * normNoSqrt(mPoints[bbn].v - srcp);
+               if(d2 >= 9.0f) continue;
+
+               // aggressive smoothing factor
+               float smstr = nvdot * pointareas[bbn];
                // Accumulate weight times field at neighbor
-               flt += *(field+pointerScale*bbn)*w;
+               target += *(field+pointerScale*bbn)*smstr;
+               smstrSum += smstr;
 
-               sum_w += w;
-               for (int i = 0; i < (int)neighbors[bbn].size(); i++) {
-                       int nn = neighbors[bbn][i];
+               for(int i = 0; i < (int)neighbors[bbn].size(); i++) {
+                       const int nn = neighbors[bbn][i];
                        if (flags[nn] == flag) continue;
-                       boundary.push_back(nn);
+                       mDboundary.push_back(nn);
                }
        }
-       flt /= sum_w;
+       target /= smstrSum;
+       return 1;
 }
 
-
-
-void IsoSurface::smoothSurface(float sigma)
+       
+// perform smoothing of the surface (and possible normals)
+void IsoSurface::smoothSurface(float sigma, bool normSmooth)
 {
        int nv = mPoints.size();
        if ((int)flags.size() != nv) flags.resize(nv);
        int nf = mIndices.size()/3;
 
        { // need neighbors
-
                vector<int> numneighbors(mPoints.size());
                int i;
                for (i = 0; i < (int)mIndices.size()/3; i++) {
@@ -886,10 +981,11 @@ void IsoSurface::smoothSurface(float sigma)
 
        vector<ntlVec3Gfx> dflt(nv);
        for (int i = 0; i < nv; i++) {
-               diffuseVertexField( &mPoints[0].v, 2,
-                                  i, invsigma2, dflt[i]);
-               // Just keep the displacement
-               dflt[i] -= mPoints[i].v;
+               if(diffuseVertexField( &mPoints[0].v, 2,
+                                  i, invsigma2, dflt[i])) {
+                       // Just keep the displacement
+                       dflt[i] -= mPoints[i].v;
+               } else { dflt[i] = 0.0; } //?mPoints[i].v; }
        }
 
        // Slightly better small-neighborhood approximation
@@ -910,8 +1006,9 @@ void IsoSurface::smoothSurface(float sigma)
        // Filter displacement field
        vector<ntlVec3Gfx> dflt2(nv);
        for (int i = 0; i < nv; i++) {
-               diffuseVertexField( &dflt[0], 1,
-                                  i, invsigma2, dflt2[i]);
+               if(diffuseVertexField( &dflt[0], 1,
+                                  i, invsigma2, dflt2[i])) { }
+               else { /*mPoints[i].v=0.0;*/ dflt2[i] = 0.0; }//dflt2[i]; }
        }
 
        // Update vertex positions
@@ -922,10 +1019,13 @@ void IsoSurface::smoothSurface(float sigma)
        // when normals smoothing off, this cleans up quite well
        // costs ca. 50% additional time though
        float nsFac = 1.5f;
-       { float ninvsigma2 = 1.0f / (nsFac*nsFac*sigma*sigma);
+       if(normSmooth) { float ninvsigma2 = 1.0f / (nsFac*nsFac*sigma*sigma);
                for (int i = 0; i < nv; i++) {
-                       diffuseVertexField( &mPoints[0].n, 2, i, ninvsigma2, dflt[i]);
-                       normalize(dflt[i]);
+                       if( diffuseVertexField( &mPoints[0].n, 2, i, ninvsigma2, dflt[i]) ) {
+                               normalize(dflt[i]);
+                       } else {
+                               dflt[i] = mPoints[i].n;
+                       }
                }
                for (int i = 0; i < nv; i++) {
                        mPoints[i].n = dflt[i];
@@ -935,9 +1035,11 @@ void IsoSurface::smoothSurface(float sigma)
        //errMsg("SMSURF","done v:"<<sigma); // DEBUG
 }
 
-void IsoSurface::smoothNormals(float sigma)
-{
-       { // need neighbor
+// only smoothen the normals
+void IsoSurface::smoothNormals(float sigma) {
+       // reuse from smoothSurface
+       if(neighbors.size() != mPoints.size()) { 
+               // need neighbor
                vector<int> numneighbors(mPoints.size());
                int i;
                for (i = 0; i < (int)mIndices.size()/3; i++) {
@@ -1039,15 +1141,13 @@ void IsoSurface::smoothNormals(float sigma)
 
        vector<ntlVec3Gfx> nflt(nv);
        for (int i = 0; i < nv; i++) {
-               diffuseVertexField( &mPoints[0].n, 2, i, invsigma2, nflt[i]);
-               normalize(nflt[i]);
-       }
-
-       for (int i = 0; i < nv; i++) {
-               mPoints[i].n = nflt[i];
+               if(diffuseVertexField( &mPoints[0].n, 2, i, invsigma2, nflt[i])) {
+                       normalize(nflt[i]);
+               } else { nflt[i]=mPoints[i].n; }
        }
 
-       //errMsg("SMNRMLS","done v:"<<sigma); // DEBUG
+       // copy results
+       for (int i = 0; i < nv; i++) { mPoints[i].n = nflt[i]; }
 }