Code got unreadble due to copy-paste (hint for me: revert point)
authorDaniel Genrich <daniel.genrich@gmx.net>
Mon, 26 Nov 2007 14:59:58 +0000 (14:59 +0000)
committerDaniel Genrich <daniel.genrich@gmx.net>
Mon, 26 Nov 2007 14:59:58 +0000 (14:59 +0000)
intern/elbeem/intern/isosurface.cpp

index 77530d413c7437464f72a9a1f13865d9f76dcb52..704f3fdacfed1d484676d322683b7915030d2817 100644 (file)
 #define round(x) (x)
 #endif
 
+#if PARALLEL==1        
+#include <omp.h>
+#endif
+
 /******************************************************************************
  * Constructor
  *****************************************************************************/
 IsoSurface::IsoSurface(double iso) :
-               ntlGeometryObject(),
-                                 mSizex(-1), mSizey(-1), mSizez(-1),
-                                         mpData(NULL),
-                                                         mIsoValue( iso ), 
-                                                                         mPoints(), 
-                                                                               mUseFullEdgeArrays(false),
-                                                                               mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL),
-                                                                               mEdgeArSize(-1),
-                                                                               mIndices(),
-
-                                                                               mStart(0.0), mEnd(0.0), mDomainExtent(0.0),
-                                                                               mInitDone(false),
-                                                                               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.)
+       ntlGeometryObject(),
+       mSizex(-1), mSizey(-1), mSizez(-1),
+       mpData(NULL),
+  mIsoValue( iso ), 
+       mPoints(), 
+       mUseFullEdgeArrays(false),
+       mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL),
+       mEdgeArSize(-1),
+       mIndices(),
+
+  mStart(0.0), mEnd(0.0), mDomainExtent(0.0),
+  mInitDone(false),
+       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.)
 {
 }
 
@@ -71,11 +75,11 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e
 
   // init 
        mIndices.clear();
-       mPoints.clear();
+  mPoints.clear();
 
        int nodes = mSizez*mSizey*mSizex;
-       mpData = new float[nodes];
-       for(int i=0;i<nodes;i++) { mpData[i] = 0.0; }
+  mpData = new float[nodes];
+  for(int i=0;i<nodes;i++) { mpData[i] = 0.0; }
 
   // allocate edge arrays  (last slices are never used...)
        int initsize = -1;
@@ -93,7 +97,7 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e
                mpEdgeVerticesZ = new int[sliceNodes];
                initsize = 3*sliceNodes;
        }
-       for(int i=0;i<mEdgeArSize;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
+  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 
@@ -106,7 +110,7 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e
 /*! Reset all values */
 void IsoSurface::resetAll(gfxReal val) {
        int nodes = mSizez*mSizey*mSizex;
-       for(int i=0;i<nodes;i++) { mpData[i] = val; }
+  for(int i=0;i<nodes;i++) { mpData[i] = val; }
 }
 
 
@@ -130,9 +134,9 @@ IsoSurface::~IsoSurface( void )
  *****************************************************************************/
 void IsoSurface::triangulate( void )
 {
-       double gsx,gsy,gsz; // grid spacing in x,y,z direction
-       double px,py,pz;    // current position in grid in x,y,z direction
-  // IsoLevelCube cubie;    // struct for a small subcube
+  double gsx,gsy,gsz; // grid spacing in x,y,z direction
+  double px,py,pz;    // current position in grid in x,y,z direction
+  IsoLevelCube cubie;    // struct for a small subcube
        myTime_t tritimestart = getTime(); 
 
        if(!mpData) {
@@ -141,16 +145,16 @@ void IsoSurface::triangulate( void )
        }
 
   // 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);
+  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<mEdgeArSize;i++) {
+  for(int i=0;i<mEdgeArSize;i++) {
                mpEdgeVerticesX[i] = -1;
                mpEdgeVerticesY[i] = -1;
                mpEdgeVerticesZ[i] = -1;
@@ -159,418 +163,421 @@ void IsoSurface::triangulate( void )
        // edges between which points?
        const int mcEdges[24] = { 
                0,1,  1,2,  2,3,  3,0,
-  4,5,  5,6,  6,7,  7,4,
-  0,4,  1,5,  2,6,  3,7 };
+               4,5,  5,6,  6,7,  7,4,
+               0,4,  1,5,  2,6,  3,7 };
 
-  const int cubieOffsetX[8] = {
-         0,1,1,0,  0,1,1,0 };
-         const int cubieOffsetY[8] = {
-                 0,0,1,1,  0,0,1,1 };
-                 const int cubieOffsetZ[8] = {
-                         0,0,0,0,  1,1,1,1 };
+       const int cubieOffsetX[8] = {
+               0,1,1,0,  0,1,1,0 };
+       const int cubieOffsetY[8] = {
+               0,0,1,1,  0,0,1,1 };
+       const int cubieOffsetZ[8] = {
+               0,0,0,0,  1,1,1,1 };
 
-                         const int coAdd=2;
+       const int coAdd=2;
   // let the cubes march 
-                         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;
+       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);
+                                       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;
+                                       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;
-                                                         }
+                                       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) ];
+                                       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);
+                                       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)) {
+                                       for(int e=0;e<12;e++) {
+                                               if (mcEdgeTable[cubeIndex] & (1<<e)) {
                                                        // is the vertex already calculated?
-                                                                         if(*eVert[ e ] < 0) {
+                                                       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);
+                                                               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 );
+                                                               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);
+                                                               triIndices[e] = (mPoints.size()-1);
                                                                // store vertex 
-                                                                               *eVert[ e ] = triIndices[e];
-                                                                         }     else {
+                                                               *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; }
-                                                                               }
+                                                               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] ] );
-                                                                               }
+                                       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
+                               }//i
+                       }// j
 
                        // copy edge arrays
-                                         if(!mUseFullEdgeArrays) {
-                                                 for(int j=0;j<(mSizey-0);j++) 
-                                                         for(int i=0;i<(mSizex-0);i++) {
+                       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
+                                       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 ); }
+               for(int ni=0;ni<(int)mPoints.size();ni++) { normalize( mPoints[ni].n ); }
 
-                         } else { // subdivs
+       } else { // subdivs
 
 #define EDGEAR_INDEX(Ai,Aj,Ak, Bi,Bj) ((mSizex*mSizey*mSubdivs*mSubdivs*(Ak))+\
-                                 (mSizex*mSubdivs*((Aj)*mSubdivs+(Bj)))+((Ai)*mSubdivs)+(Bi))
+               (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] )
+                               (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!");
-                                 }
+               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];
+               ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex];
 
                // construct pointers
                // part test
-                                 int pInUse = 0;
-                                 int pUsedTest = 0;
+               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);
-                                                                               }
+               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;
+                       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
+                               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) {
+                               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();
-                                                                               }
-                                                                               }
-                                                                               } else {
+                                       ParticleObject* listpnt = pnt;
+                                       while(listpnt) {
+                                               if(!listpnt->getNext()) {
+                                                       listpnt->setNext(p); listpnt = NULL;
+                                               } else {
+                                                       listpnt = listpnt->getNext();
+                                               }
+                                       }
+                               } 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
+                                       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
-                                                         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!
+#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);
+                                       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);
-                                                                               }
-
-                                                                               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) { // */
+                                       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);
+                                                       }
+
+                                       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++) { { { {
+                                       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?
+                                               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;
+                                                       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(
+                                                       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 { 
+                                                               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; }
+                                                                       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]; 
+                                                               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) {
+                                                       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 ];                             
+                                                       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[ 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) ];
+                                                       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);
+                                                       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)) {
+                                                       for(int e=0;e<12;e++) {
+                                                               if (mcEdgeTable[cubeIndex] & (1<<e)) {
                                                                        // is the vertex already calculated?
-                                                                               if(*eVert[ e ] < 0) {
+                                                                       if(*eVert[ e ] < 0) {
                                                                                // interpolate edge
                                                                                const int e1 = mcEdges[e*2  ];
                                                                                const int e2 = mcEdges[e*2+1];
@@ -582,79 +589,84 @@ void IsoSurface::triangulate( void )
 
                                                                                // init isolevel vertex
                                                                                ilv.v = p1 + (p2-p1)*mu; // with subdivs
-#pragma omp critical 
+#if PARALLEL==1        
+#pragma omp critical
+#endif
                                                                                {
                                                                                mPoints.push_back( ilv );
                                                                                triIndices[e] = (mPoints.size()-1);
                                                                                }
                                                                                // store vertex 
                                                                                *eVert[ e ] = triIndices[e]; 
-                                                                               }       else {
+                                                                       }       else {
                                                                                // retrieve  from vert array
                                                                                triIndices[e] = *eVert[ e ];
-                                                                               }
-                                                                               } // along all edges 
-                                                                               }
+                                                                       }
+                                                               } // along all edges 
+                                                       }
                                                        // removed cutoff treatment...
-#pragma omp critical 
-                                                                               {
+                                                       
                                                        // 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] ] ); // with subdivs
-                                                                               mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
+#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?
+                                                       } // triangles in edge table?
+                                                       }
+                                                       }
                                                        
-                                                                               }//si
-                                                                               }// sj
-
-                                                                 }//i
-                                                         }// j
+                                               }//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++) {
+                       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
+                                       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
 
                //delete [] subdAr;
-                                                 delete [] arppnt;
-                                                 computeNormals();
-                         } // with subdivs
+               delete [] arppnt;
+               computeNormals();
+       } // with subdivs
 
        // perform smoothing
-                         float smoSubdfac = 1.;
-                         if(mSubdivs>0) {
+       float smoSubdfac = 1.;
+       if(mSubdivs>0) {
                //smoSubdfac = 1./(float)(mSubdivs);
-                                 smoSubdfac = pow(0.55,(double)mSubdivs); // slightly stronger
-                         }
-                         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*smoSubdfac); 
-                         }
-
-                         myTime_t tritimeend = getTime(); 
-                         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);
+               smoSubdfac = pow(0.55,(double)mSubdivs); // slightly stronger
+       }
+       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*smoSubdfac); 
+       }
+
+       myTime_t tritimeend = getTime(); 
+       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);
 }
 
 
@@ -665,8 +677,8 @@ void IsoSurface::triangulate( void )
  * Get triangles for rendering
  *****************************************************************************/
 void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles, 
-                             vector<ntlVec3Gfx> *vertices, 
-        vector<ntlVec3Gfx> *normals, int objectId )
+                                                                                                        vector<ntlVec3Gfx> *vertices, 
+                                                                                                        vector<ntlVec3Gfx> *normals, int objectId )
 {
        if(!mInitDone) {
                debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10);
@@ -675,8 +687,8 @@ void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles,
        t = 0.;
        //return; // DEBUG
 
-       /* triangulate field */
-       triangulate();
+  /* triangulate field */
+  triangulate();
        //errMsg("TRIS"," "<<mIndices.size() );
 
        // new output with vertice reuse
@@ -689,20 +701,20 @@ void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles,
        //errMsg("NM"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
        //errMsg("NM"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
 
-       for(int i=0;i<(int)mPoints.size();i++) {
+  for(int i=0;i<(int)mPoints.size();i++) {
                vertices->push_back( mPoints[i].v );
        }
-       for(int i=0;i<(int)mPoints.size();i++) {
+  for(int i=0;i<(int)mPoints.size();i++) {
                normals->push_back( mPoints[i].n );
        }
 
        //errMsg("N2"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
        //errMsg("N2"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
 
-       for(int i=0;i<(int)mIndices.size();i+=3) {
+  for(int i=0;i<(int)mIndices.size();i+=3) {
                const int smooth = 1;
-               int t1 = mIndices[i];
-               int t2 = mIndices[i+1];
+    int t1 = mIndices[i];
+    int t2 = mIndices[i+1];
                int t3 = mIndices[i+2];
                //errMsg("NM"," tri"<<t1<<" "<<t2<<" "<<t3 );
 
@@ -718,20 +730,20 @@ void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles,
                if(getCastShadows() ) { 
                        flag |= TRI_CASTSHADOWS; } 
 
-                       /* init geo init id */
-                       int geoiId = getGeoInitId(); 
-                       if(geoiId > 0) { 
-                               flag |= (1<< (geoiId+4)); 
-                               flag |= mGeoInitType; 
-                       
+               /* init geo init id */
+               int geoiId = getGeoInitId(); 
+               if(geoiId > 0) { 
+                       flag |= (1<< (geoiId+4)); 
+                       flag |= mGeoInitType; 
+               } 
 
-                       tri.setFlags( flag );
+               tri.setFlags( flag );
 
-                       /* triangle normal missing */
-                       tri.setNormal( ntlVec3Gfx(0.0) );
-                       tri.setSmoothNormals( smooth );
-                       tri.setObjectId( objectId );
-                       triangles->push_back( tri ); 
+               /* triangle normal missing */
+               tri.setNormal( ntlVec3Gfx(0.0) );
+               tri.setSmoothNormals( smooth );
+               tri.setObjectId( objectId );
+               triangles->push_back( tri ); 
        }
        //errMsg("N3"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
        return;
@@ -743,11 +755,11 @@ inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) {
        // WARNING - this requires a security boundary layer... 
        ntlVec3Gfx ret(0.0);
        ret[0] = *getData(i-1,j  ,k  ) - 
-                       *getData(i+1,j  ,k  );
+                *getData(i+1,j  ,k  );
        ret[1] = *getData(i  ,j-1,k  ) - 
-                       *getData(i  ,j+1,k  );
+                *getData(i  ,j+1,k  );
        ret[2] = *getData(i  ,j  ,k-1  ) - 
-                       *getData(i  ,j  ,k+1  );
+                *getData(i  ,j  ,k+1  );
        return ret;
 }
 
@@ -769,13 +781,13 @@ void IsoSurface::setSmoothRad(float radi1, float radi2, ntlVec3Gfx mscc) {
 
 // compute normals for all generated triangles
 void IsoSurface::computeNormals() {
-       for(int i=0;i<(int)mPoints.size();i++) {
+  for(int i=0;i<(int)mPoints.size();i++) {
                mPoints[i].n = ntlVec3Gfx(0.);
        }
 
-       for(int i=0;i<(int)mIndices.size();i+=3) {
-               const int t1 = mIndices[i];
-               const int t2 = mIndices[i+1];
+  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;
@@ -792,7 +804,7 @@ void IsoSurface::computeNormals() {
                mPoints[t3].n += norm * (1./(len2*len3));
        }
 
-       for(int i=0;i<(int)mPoints.size();i++) {
+  for(int i=0;i<(int)mPoints.size();i++) {
                normalize(mPoints[i].n);
        }
 }
@@ -903,86 +915,86 @@ void IsoSurface::smoothSurface(float sigma, bool normSmooth)
                        // Edges
                        ntlVec3Gfx e[3] = { 
                                mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
-    mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
-    mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
+                               mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
+                               mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
 
                        // Compute corner weights
-    float area = 0.5f * norm( cross(e[0], e[1]));
-    float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
-    float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
-           l2[1] * (l2[2] + l2[0] - l2[1]),
-                    l2[2] * (l2[0] + l2[1] - l2[2]) };
-                    if (ew[0] <= 0.0f) {
-                            cornerareas[i][1] = -0.25f * l2[2] * area /
-                                            dot(e[0] , e[2]);
-                            cornerareas[i][2] = -0.25f * l2[1] * area /
-                                            dot(e[0] , e[1]);
-                            cornerareas[i][0] = area - cornerareas[i][1] -
-                                            cornerareas[i][2];
-                    } else if (ew[1] <= 0.0f) {
-                            cornerareas[i][2] = -0.25f * l2[0] * area /
-                                            dot(e[1] , e[0]);
-                            cornerareas[i][0] = -0.25f * l2[2] * area /
-                                            dot(e[1] , e[2]);
-                            cornerareas[i][1] = area - cornerareas[i][2] -
-                                            cornerareas[i][0];
-                    } else if (ew[2] <= 0.0f) {
-                            cornerareas[i][0] = -0.25f * l2[1] * area /
-                                            dot(e[2] , e[1]);
-                            cornerareas[i][1] = -0.25f * l2[0] * area /
-                                            dot(e[2] , e[0]);
-                            cornerareas[i][2] = area - cornerareas[i][0] -
-                                            cornerareas[i][1];
-                    } else {
-                            float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
-                            for (int j = 0; j < 3; j++)
-                                    cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
-                                                    ew[(j+2)%3]);
-                    }
+                       float area = 0.5f * norm( cross(e[0], e[1]));
+                       float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
+                       float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
+                                       l2[1] * (l2[2] + l2[0] - l2[1]),
+                                       l2[2] * (l2[0] + l2[1] - l2[2]) };
+                       if (ew[0] <= 0.0f) {
+                               cornerareas[i][1] = -0.25f * l2[2] * area /
+                                                               dot(e[0] , e[2]);
+                               cornerareas[i][2] = -0.25f * l2[1] * area /
+                                                               dot(e[0] , e[1]);
+                               cornerareas[i][0] = area - cornerareas[i][1] -
+                                                               cornerareas[i][2];
+                       } else if (ew[1] <= 0.0f) {
+                               cornerareas[i][2] = -0.25f * l2[0] * area /
+                                                               dot(e[1] , e[0]);
+                               cornerareas[i][0] = -0.25f * l2[2] * area /
+                                                               dot(e[1] , e[2]);
+                               cornerareas[i][1] = area - cornerareas[i][2] -
+                                                               cornerareas[i][0];
+                       } else if (ew[2] <= 0.0f) {
+                               cornerareas[i][0] = -0.25f * l2[1] * area /
+                                                               dot(e[2] , e[1]);
+                               cornerareas[i][1] = -0.25f * l2[0] * area /
+                                                               dot(e[2] , e[0]);
+                               cornerareas[i][2] = area - cornerareas[i][0] -
+                                                               cornerareas[i][1];
+                       } else {
+                               float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
+                               for (int j = 0; j < 3; j++)
+                                       cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
+                                                                                        ew[(j+2)%3]);
+                       }
 
                        // NT important, check this...
 #ifndef WIN32
-                    if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
-                    if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
-                    if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
+                       if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
+                       if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
+                       if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
 #else // WIN32
                        // FIXME check as well...
-                    if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
-                    if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
-                    if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
+                       if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
+                       if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
+                       if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
 #endif // WIN32
 
-                    pointareas[mIndices[i*3+0]] += cornerareas[i][0];
-                    pointareas[mIndices[i*3+1]] += cornerareas[i][1];
-                    pointareas[mIndices[i*3+2]] += cornerareas[i][2];
+                       pointareas[mIndices[i*3+0]] += cornerareas[i][0];
+                       pointareas[mIndices[i*3+1]] += cornerareas[i][1];
+                       pointareas[mIndices[i*3+2]] += cornerareas[i][2];
                }
 
        } // need pointarea
-       // */
+       // */
 
        float invsigma2 = 1.0f / (sigma*sigma);
 
        vector<ntlVec3Gfx> dflt(nv);
        for (int i = 0; i < nv; i++) {
                if(diffuseVertexField( &mPoints[0].v, 2,
-                  i, invsigma2, dflt[i])) {
+                                  i, invsigma2, dflt[i])) {
                        // Just keep the displacement
-                          dflt[i] -= mPoints[i].v;
-                  } else { dflt[i] = 0.0; } //?mPoints[i].v; }
+                       dflt[i] -= mPoints[i].v;
+               } else { dflt[i] = 0.0; } //?mPoints[i].v; }
        }
 
        // Slightly better small-neighborhood approximation
        for (int i = 0; i < nf; i++) {
                ntlVec3Gfx c = mPoints[mIndices[i*3+0]].v +
-                               mPoints[mIndices[i*3+1]].v +
-                               mPoints[mIndices[i*3+2]].v;
+                         mPoints[mIndices[i*3+1]].v +
+                         mPoints[mIndices[i*3+2]].v;
                c /= 3.0f;
                for (int j = 0; j < 3; j++) {
                        int v = mIndices[i*3+j];
                        ntlVec3Gfx d =(c - mPoints[v].v) * 0.5f;
                        dflt[v] += d * (cornerareas[i][j] /
-                                       pointareas[mIndices[i*3+j]] *
-                                       exp(-0.5f * invsigma2 * normNoSqrt(d)) );
+                                  pointareas[mIndices[i*3+j]] *
+                                  exp(-0.5f * invsigma2 * normNoSqrt(d)) );
                }
        }
 
@@ -990,7 +1002,7 @@ void IsoSurface::smoothSurface(float sigma, bool normSmooth)
        vector<ntlVec3Gfx> dflt2(nv);
        for (int i = 0; i < nv; i++) {
                if(diffuseVertexField( &dflt[0], 1,
-                  i, invsigma2, dflt2[i])) { }
+                                  i, invsigma2, dflt2[i])) { }
                else { /*mPoints[i].v=0.0;*/ dflt2[i] = 0.0; }//dflt2[i]; }
        }
 
@@ -1062,58 +1074,58 @@ void IsoSurface::smoothNormals(float sigma) {
                        // Edges
                        ntlVec3Gfx e[3] = { 
                                mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
-    mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
-    mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
+                               mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
+                               mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
 
                        // Compute corner weights
-    float area = 0.5f * norm( cross(e[0], e[1]));
-    float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
-    float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
-           l2[1] * (l2[2] + l2[0] - l2[1]),
-                    l2[2] * (l2[0] + l2[1] - l2[2]) };
-                    if (ew[0] <= 0.0f) {
-                            cornerareas[i][1] = -0.25f * l2[2] * area /
-                                            dot(e[0] , e[2]);
-                            cornerareas[i][2] = -0.25f * l2[1] * area /
-                                            dot(e[0] , e[1]);
-                            cornerareas[i][0] = area - cornerareas[i][1] -
-                                            cornerareas[i][2];
-                    } else if (ew[1] <= 0.0f) {
-                            cornerareas[i][2] = -0.25f * l2[0] * area /
-                                            dot(e[1] , e[0]);
-                            cornerareas[i][0] = -0.25f * l2[2] * area /
-                                            dot(e[1] , e[2]);
-                            cornerareas[i][1] = area - cornerareas[i][2] -
-                                            cornerareas[i][0];
-                    } else if (ew[2] <= 0.0f) {
-                            cornerareas[i][0] = -0.25f * l2[1] * area /
-                                            dot(e[2] , e[1]);
-                            cornerareas[i][1] = -0.25f * l2[0] * area /
-                                            dot(e[2] , e[0]);
-                            cornerareas[i][2] = area - cornerareas[i][0] -
-                                            cornerareas[i][1];
-                    } else {
-                            float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
-                            for (int j = 0; j < 3; j++)
-                                    cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
-                                                    ew[(j+2)%3]);
-                    }
+                       float area = 0.5f * norm( cross(e[0], e[1]));
+                       float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
+                       float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
+                                       l2[1] * (l2[2] + l2[0] - l2[1]),
+                                       l2[2] * (l2[0] + l2[1] - l2[2]) };
+                       if (ew[0] <= 0.0f) {
+                               cornerareas[i][1] = -0.25f * l2[2] * area /
+                                                               dot(e[0] , e[2]);
+                               cornerareas[i][2] = -0.25f * l2[1] * area /
+                                                               dot(e[0] , e[1]);
+                               cornerareas[i][0] = area - cornerareas[i][1] -
+                                                               cornerareas[i][2];
+                       } else if (ew[1] <= 0.0f) {
+                               cornerareas[i][2] = -0.25f * l2[0] * area /
+                                                               dot(e[1] , e[0]);
+                               cornerareas[i][0] = -0.25f * l2[2] * area /
+                                                               dot(e[1] , e[2]);
+                               cornerareas[i][1] = area - cornerareas[i][2] -
+                                                               cornerareas[i][0];
+                       } else if (ew[2] <= 0.0f) {
+                               cornerareas[i][0] = -0.25f * l2[1] * area /
+                                                               dot(e[2] , e[1]);
+                               cornerareas[i][1] = -0.25f * l2[0] * area /
+                                                               dot(e[2] , e[0]);
+                               cornerareas[i][2] = area - cornerareas[i][0] -
+                                                               cornerareas[i][1];
+                       } else {
+                               float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
+                               for (int j = 0; j < 3; j++)
+                                       cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
+                                                                                        ew[(j+2)%3]);
+                       }
 
                        // NT important, check this...
 #ifndef WIN32
-                    if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
-                    if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
-                    if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
+                       if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
+                       if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
+                       if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
 #else // WIN32
                        // FIXME check as well...
-                    if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
-                    if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
-                    if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
+                       if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
+                       if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
+                       if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
 #endif // WIN32
 
-                    pointareas[mIndices[i*3+0]] += cornerareas[i][0];
-                    pointareas[mIndices[i*3+1]] += cornerareas[i][1];
-                    pointareas[mIndices[i*3+2]] += cornerareas[i][2];
+                       pointareas[mIndices[i*3+0]] += cornerareas[i][0];
+                       pointareas[mIndices[i*3+1]] += cornerareas[i][1];
+                       pointareas[mIndices[i*3+2]] += cornerareas[i][2];
                }
 
        } // need pointarea