Reverted to rev 12673 + test for ccherett
[blender.git] / intern / elbeem / intern / isosurface.cpp
1 /******************************************************************************
2  *
3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4  * Copyright 2003-2006 Nils Thuerey
5  *
6  * Marching Cubes surface mesh generation
7  *
8  *****************************************************************************/
9
10 #include "isosurface.h"
11 #include "mcubes_tables.h"
12 #include "particletracer.h"
13 #include <algorithm>
14 #include <stdio.h>
15
16 // just use default rounding for platforms where its not available
17 #ifndef round
18 #define round(x) (x)
19 #endif
20
21 /******************************************************************************
22  * Constructor
23  *****************************************************************************/
24 IsoSurface::IsoSurface(double iso) :
25         ntlGeometryObject(),
26         mSizex(-1), mSizey(-1), mSizez(-1),
27         mpData(NULL),
28   mIsoValue( iso ), 
29         mPoints(), 
30         mUseFullEdgeArrays(false),
31         mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL),
32         mEdgeArSize(-1),
33         mIndices(),
34
35   mStart(0.0), mEnd(0.0), mDomainExtent(0.0),
36   mInitDone(false),
37         mSmoothSurface(0.0), mSmoothNormals(0.0),
38         mAcrossEdge(), mAdjacentFaces(),
39         mCutoff(-1), mCutArray(NULL), // off by default
40         mpIsoParts(NULL), mPartSize(0.), mSubdivs(0),
41         mFlagCnt(1),
42         mSCrad1(0.), mSCrad2(0.), mSCcenter(0.)
43 {
44 }
45
46
47 /******************************************************************************
48  * The real init...
49  *****************************************************************************/
50 void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent) 
51 {
52         // range 1-10 (max due to subd array in triangulate)
53         if(mSubdivs<1) mSubdivs=1;
54         if(mSubdivs>10) mSubdivs=10;
55
56         // init solver and size
57         mSizex = setx;
58         mSizey = sety;
59         if(setz == 1) {// 2D, create thin 2D surface
60                 setz = 5; 
61         }
62         mSizez = setz;
63         mDomainExtent = extent;
64
65         /* check triangulation size (for raytraing) */
66         if( ( mStart[0] >= mEnd[0] ) && ( mStart[1] >= mEnd[1] ) && ( mStart[2] >= mEnd[2] ) ){
67                 // extent was not set, use normalized one from parametrizer
68                 mStart = ntlVec3Gfx(0.0) - extent*0.5;
69                 mEnd = ntlVec3Gfx(0.0) + extent*0.5;
70         }
71
72   // init 
73         mIndices.clear();
74   mPoints.clear();
75
76         int nodes = mSizez*mSizey*mSizex;
77   mpData = new float[nodes];
78   for(int i=0;i<nodes;i++) { mpData[i] = 0.0; }
79
80   // allocate edge arrays  (last slices are never used...)
81         int initsize = -1;
82         if(mUseFullEdgeArrays) {
83                 mEdgeArSize = nodes;
84                 mpEdgeVerticesX = new int[nodes];
85                 mpEdgeVerticesY = new int[nodes];
86                 mpEdgeVerticesZ = new int[nodes];
87                 initsize = 3*nodes;
88         } else {
89                 int sliceNodes = 2*mSizex*mSizey*mSubdivs*mSubdivs;
90                 mEdgeArSize = sliceNodes;
91                 mpEdgeVerticesX = new int[sliceNodes];
92                 mpEdgeVerticesY = new int[sliceNodes];
93                 mpEdgeVerticesZ = new int[sliceNodes];
94                 initsize = 3*sliceNodes;
95         }
96   for(int i=0;i<mEdgeArSize;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
97         // WARNING - make sure this is consistent with calculateMemreqEstimate
98   
99         // marching cubes are ready 
100         mInitDone = true;
101         debMsgStd("IsoSurface::initializeIsosurface",DM_MSG,"Inited, edgenodes:"<<initsize<<" subdivs:"<<mSubdivs<<" fulledg:"<<mUseFullEdgeArrays , 10);
102 }
103
104
105
106 /*! Reset all values */
107 void IsoSurface::resetAll(gfxReal val) {
108         int nodes = mSizez*mSizey*mSizex;
109   for(int i=0;i<nodes;i++) { mpData[i] = val; }
110 }
111
112
113 /******************************************************************************
114  * Destructor
115  *****************************************************************************/
116 IsoSurface::~IsoSurface( void )
117 {
118         if(mpData) delete [] mpData;
119         if(mpEdgeVerticesX) delete [] mpEdgeVerticesX;
120         if(mpEdgeVerticesY) delete [] mpEdgeVerticesY;
121         if(mpEdgeVerticesZ) delete [] mpEdgeVerticesZ;
122 }
123
124
125
126
127
128 /******************************************************************************
129  * triangulate the scalar field given by pointer
130  *****************************************************************************/
131 void IsoSurface::triangulate( void )
132 {
133   double gsx,gsy,gsz; // grid spacing in x,y,z direction
134   double px,py,pz;    // current position in grid in x,y,z direction
135   // IsoLevelCube cubie;    // struct for a small subcube
136         myTime_t tritimestart = getTime(); 
137
138         if(!mpData) {
139                 errFatal("IsoSurface::triangulate","no LBM object, and no scalar field...!",SIMWORLD_INITERROR);
140                 return;
141         }
142
143   // get grid spacing (-2 to have same spacing as sim)
144   gsx = (mEnd[0]-mStart[0])/(double)(mSizex-2.0);
145   gsy = (mEnd[1]-mStart[1])/(double)(mSizey-2.0);
146   gsz = (mEnd[2]-mStart[2])/(double)(mSizez-2.0);
147
148   // clean up previous frame
149         mIndices.clear();
150         mPoints.clear();
151
152         // reset edge vertices
153   for(int i=0;i<mEdgeArSize;i++) {
154                 mpEdgeVerticesX[i] = -1;
155                 mpEdgeVerticesY[i] = -1;
156                 mpEdgeVerticesZ[i] = -1;
157         }
158
159         // edges between which points?
160         const int mcEdges[24] = { 
161                 0,1,  1,2,  2,3,  3,0,
162                 4,5,  5,6,  6,7,  7,4,
163                 0,4,  1,5,  2,6,  3,7 };
164
165         const int cubieOffsetX[8] = {
166                 0,1,1,0,  0,1,1,0 };
167         const int cubieOffsetY[8] = {
168                 0,0,1,1,  0,0,1,1 };
169         const int cubieOffsetZ[8] = {
170                 0,0,0,0,  1,1,1,1 };
171
172         const int coAdd=2;
173   // let the cubes march 
174         if(mSubdivs<=1) {
175
176                 pz = mStart[2]-gsz*0.5;
177                 for(int k=1;k<(mSizez-2);k++) {
178                         pz += gsz;
179                         py = mStart[1]-gsy*0.5;
180                         for(int j=1;j<(mSizey-2);j++) {
181                                 py += gsy;
182                                 px = mStart[0]-gsx*0.5;
183                                 for(int i=1;i<(mSizex-2);i++) {
184                                         px += gsx;
185                                         int cubeIndex;      // index entry of the cube 
186                                         float value[8];
187                                         int triIndices[12]; // vertex indices 
188                                         int *eVert[12];
189                                         IsoLevelVertex ilv;
190                                         
191                                         value[0] = *getData(i  ,j  ,k  );
192                                         value[1] = *getData(i+1,j  ,k  );
193                                         value[2] = *getData(i+1,j+1,k  );
194                                         value[3] = *getData(i  ,j+1,k  );
195                                         value[4] = *getData(i  ,j  ,k+1);
196                                         value[5] = *getData(i+1,j  ,k+1);
197                                         value[6] = *getData(i+1,j+1,k+1);
198                                         value[7] = *getData(i  ,j+1,k+1);
199
200                                         // check intersections of isosurface with edges, and calculate cubie index
201                                         cubeIndex = 0;
202                                         if (value[0] < mIsoValue) cubeIndex |= 1;
203                                         if (value[1] < mIsoValue) cubeIndex |= 2;
204                                         if (value[2] < mIsoValue) cubeIndex |= 4;
205                                         if (value[3] < mIsoValue) cubeIndex |= 8;
206                                         if (value[4] < mIsoValue) cubeIndex |= 16;
207                                         if (value[5] < mIsoValue) cubeIndex |= 32;
208                                         if (value[6] < mIsoValue) cubeIndex |= 64;
209                                         if (value[7] < mIsoValue) cubeIndex |= 128;
210
211                                         // No triangles to generate?
212                                         if (mcEdgeTable[cubeIndex] == 0) {
213                                                 continue;
214                                         }
215
216                                         // where to look up if this point already exists
217                                         int edgek = 0;
218                                         if(mUseFullEdgeArrays) edgek=k;
219                                         const int baseIn = ISOLEVEL_INDEX( i+0, j+0, edgek+0);
220                                         eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
221                                         eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
222                                         eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
223                                         eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
224
225                                         eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
226                                         eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, edgek+1) ];
227                                         eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+1) ];
228                                         eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
229
230                                         eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
231                                         eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, edgek+0) ];
232                                         eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, edgek+0) ];
233                                         eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
234
235                                         // grid positions
236                                         ntlVec3Gfx pos[8];
237                                         pos[0] = ntlVec3Gfx(px    ,py    ,pz);
238                                         pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
239                                         pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
240                                         pos[3] = ntlVec3Gfx(px    ,py+gsy,pz);
241                                         pos[4] = ntlVec3Gfx(px    ,py    ,pz+gsz);
242                                         pos[5] = ntlVec3Gfx(px+gsx,py    ,pz+gsz);
243                                         pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz);
244                                         pos[7] = ntlVec3Gfx(px    ,py+gsy,pz+gsz);
245
246                                         // check all edges
247                                         for(int e=0;e<12;e++) {
248                                                 if (mcEdgeTable[cubeIndex] & (1<<e)) {
249                                                         // is the vertex already calculated?
250                                                         if(*eVert[ e ] < 0) {
251                                                                 // interpolate edge
252                                                                 const int e1 = mcEdges[e*2  ];
253                                                                 const int e2 = mcEdges[e*2+1];
254                                                                 const ntlVec3Gfx p1 = pos[ e1  ];    // scalar field pos 1
255                                                                 const ntlVec3Gfx p2 = pos[ e2  ];    // scalar field pos 2
256                                                                 const float valp1  = value[ e1  ];  // scalar field val 1
257                                                                 const float valp2  = value[ e2  ];  // scalar field val 2
258                                                                 const float mu = (mIsoValue - valp1) / (valp2 - valp1);
259
260                                                                 // init isolevel vertex
261                                                                 ilv.v = p1 + (p2-p1)*mu;
262                                                                 ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) +
263                                                                                                 getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * (    mu) ;
264                                                                 mPoints.push_back( ilv );
265
266                                                                 triIndices[e] = (mPoints.size()-1);
267                                                                 // store vertex 
268                                                                 *eVert[ e ] = triIndices[e];
269                                                         }       else {
270                                                                 // retrieve  from vert array
271                                                                 triIndices[e] = *eVert[ e ];
272                                                         }
273                                                 } // along all edges 
274                                         }
275
276                                         if( (i<coAdd+mCutoff) || (j<coAdd+mCutoff) ||
277                                                         ((mCutoff>0) && (k<coAdd)) ||// bottom layer
278                                                         (i>mSizex-2-coAdd-mCutoff) ||
279                                                         (j>mSizey-2-coAdd-mCutoff) ) {
280                                                 if(mCutArray) {
281                                                         if(k < mCutArray[j*this->mSizex+i]) continue;
282                                                 } else { continue; }
283                                         }
284
285                                         // Create the triangles... 
286                                         for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
287                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
288                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] );
289                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
290                                         }
291                                         
292                                 }//i
293                         }// j
294
295                         // copy edge arrays
296                         if(!mUseFullEdgeArrays) {
297                         for(int j=0;j<(mSizey-0);j++) 
298                                 for(int i=0;i<(mSizex-0);i++) {
299                                         //int edgek = 0;
300                                         const int dst = ISOLEVEL_INDEX( i+0, j+0, 0);
301                                         const int src = ISOLEVEL_INDEX( i+0, j+0, 1);
302                                         mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ];
303                                         mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ];
304                                         mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ];
305                                         mpEdgeVerticesX[ src ]=-1;
306                                         mpEdgeVerticesY[ src ]=-1;
307                                         mpEdgeVerticesZ[ src ]=-1;
308                                 }
309                         } // */
310
311                 } // k
312
313         // precalculate normals using an approximation of the scalar field gradient 
314                 for(int ni=0;ni<(int)mPoints.size();ni++) { normalize( mPoints[ni].n ); }
315
316         } else { // subdivs
317
318 #define EDGEAR_INDEX(Ai,Aj,Ak, Bi,Bj) ((mSizex*mSizey*mSubdivs*mSubdivs*(Ak))+\
319                 (mSizex*mSubdivs*((Aj)*mSubdivs+(Bj)))+((Ai)*mSubdivs)+(Bi))
320
321 #define ISOTRILININT(fi,fj,fk) ( \
322                                 (1.-(fi))*(1.-(fj))*(1.-(fk))*orgval[0] + \
323                                 (   (fi))*(1.-(fj))*(1.-(fk))*orgval[1] + \
324                                 (   (fi))*(   (fj))*(1.-(fk))*orgval[2] + \
325                                 (1.-(fi))*(   (fj))*(1.-(fk))*orgval[3] + \
326                                 (1.-(fi))*(1.-(fj))*(   (fk))*orgval[4] + \
327                                 (   (fi))*(1.-(fj))*(   (fk))*orgval[5] + \
328                                 (   (fi))*(   (fj))*(   (fk))*orgval[6] + \
329                                 (1.-(fi))*(   (fj))*(   (fk))*orgval[7] )
330
331                 // use subdivisions
332                 gfxReal subdfac = 1./(gfxReal)(mSubdivs);
333                 gfxReal orgGsx = gsx;
334                 gfxReal orgGsy = gsy;
335                 gfxReal orgGsz = gsz;
336                 gsx *= subdfac;
337                 gsy *= subdfac;
338                 gsz *= subdfac;
339                 if(mUseFullEdgeArrays) {
340                         errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!");
341                 }
342                 
343                 ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex];
344
345                 // construct pointers
346                 // part test
347                 int pInUse = 0;
348                 int pUsedTest = 0;
349                 // reset particles
350                 // reset list array
351                 for(int k=0;k<(mSizez);k++) 
352                         for(int j=0;j<(mSizey);j++) 
353                                 for(int i=0;i<(mSizex);i++) {
354                                         arppnt[ISOLEVEL_INDEX(i,j,k)] = NULL;
355                                 }
356                 if(mpIsoParts) {
357                         for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
358                                         pit!= mpIsoParts->getParticlesEnd(); pit++) {
359                                 if( (*pit).getActive()==false ) continue;
360                                 if( (*pit).getType()!=PART_DROP) continue;
361                                 (*pit).setNext(NULL);
362                         }
363                         // build per node lists
364                         for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
365                                         pit!= mpIsoParts->getParticlesEnd(); pit++) {
366                                 if( (*pit).getActive()==false ) continue;
367                                 if( (*pit).getType()!=PART_DROP) continue;
368                                 // check lifetime ignored here
369                                 ParticleObject *p = &(*pit);
370                                 const ntlVec3Gfx ppos = p->getPos();
371                                 const int pi= (int)round(ppos[0])+0; 
372                                 const int pj= (int)round(ppos[1])+0; 
373                                 int       pk= (int)round(ppos[2])+0;// no offset necessary
374                                 // 2d should be handled by solver. if(LBMDIM==2) { pk = 0; }
375
376                                 if(pi<0) continue;
377                                 if(pj<0) continue;
378                                 if(pk<0) continue;
379                                 if(pi>mSizex-1) continue;
380                                 if(pj>mSizey-1) continue;
381                                 if(pk>mSizez-1) continue;
382                                 ParticleObject* &pnt = arppnt[ISOLEVEL_INDEX(pi,pj,pk)]; 
383                                 if(pnt) {
384                                         // append
385                                         ParticleObject* listpnt = pnt;
386                                         while(listpnt) {
387                                                 if(!listpnt->getNext()) {
388                                                         listpnt->setNext(p); listpnt = NULL;
389                                                 } else {
390                                                         listpnt = listpnt->getNext();
391                                                 }
392                                         }
393                                 } else {
394                                         // start new list
395                                         pnt = p;
396                                 }
397                                 pInUse++;
398                         }
399                 } // mpIsoParts
400
401                 debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<<pInUse<<", Subdivs:"<<mSubdivs, 9);
402                 pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz;
403                 for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) {
404                         pz += gsz;
405                         const int k = ok/mSubdivs;
406                         if(k<=0) continue; // skip zero plane
407 #pragma omp parallel for
408                         for(int j=1;j<(mSizey-2);j++) {
409                                 for(int i=1;i<(mSizex-2);i++) {
410                                         float value[8];
411                                         ntlVec3Gfx pos[8];
412                                         int cubeIndex;      // index entry of the cube 
413                                         int triIndices[12]; // vertex indices 
414                                         int *eVert[12];
415                                         IsoLevelVertex ilv;
416                                         gfxReal orgval[8];
417                                         gfxReal subdAr[2][11][11]; // max 10 subdivs!
418                                         
419                                         orgval[0] = *getData(i  ,j  ,k  );
420                                         orgval[1] = *getData(i+1,j  ,k  );
421                                         orgval[2] = *getData(i+1,j+1,k  ); // with subdivs
422                                         orgval[3] = *getData(i  ,j+1,k  );
423                                         orgval[4] = *getData(i  ,j  ,k+1);
424                                         orgval[5] = *getData(i+1,j  ,k+1);
425                                         orgval[6] = *getData(i+1,j+1,k+1); // with subdivs
426                                         orgval[7] = *getData(i  ,j+1,k+1);
427
428                                         // prebuild subsampled array slice
429                                         const int sdkOffset = ok-k*mSubdivs; 
430
431                                         for(int sdk=0; sdk<2; sdk++) 
432                                                 for(int sdj=0; sdj<mSubdivs+1; sdj++) 
433                                                         for(int sdi=0; sdi<mSubdivs+1; sdi++) {
434                                                                 subdAr[sdk][sdj][sdi] = ISOTRILININT(sdi*subdfac, sdj*subdfac, (sdkOffset+sdk)*subdfac);
435                                                         }
436
437                                         const int poDistOffset=2;
438                                         for(int pok=-poDistOffset; pok<1+poDistOffset; pok++) {
439                                                 if(k+pok<0) continue;
440                                                 if(k+pok>=mSizez-1) continue;
441                                         for(int poj=-poDistOffset; poj<1+poDistOffset; poj++) {
442                                                 if(j+poj<0) continue;
443                                                 if(j+poj>=mSizey-1) continue;
444                                         for(int poi=-poDistOffset; poi<1+poDistOffset; poi++) {
445                                                 if(i+poi<0) continue;
446                                                 if(i+poi>=mSizex-1) continue; 
447                                                 ParticleObject *p;
448                                                 p = arppnt[ISOLEVEL_INDEX(i+poi,j+poj,k+pok)];
449                                                 while(p) { // */
450                                         /*
451                                         for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
452                                                         pit!= mpIsoParts->getParticlesEnd(); pit++) { { { {
453                                                 // debug test! , full list slow!
454                                                 if(( (*pit).getActive()==false ) || ( (*pit).getType()!=PART_DROP)) continue;
455                                                 ParticleObject *p;
456                                                 p = &(*pit); // */
457
458                                                         pUsedTest++;
459                                                         ntlVec3Gfx ppos = p->getPos();
460                                                         const int spi= (int)round( (ppos[0]+1.-(gfxReal)i) *(gfxReal)mSubdivs-1.5); 
461                                                         const int spj= (int)round( (ppos[1]+1.-(gfxReal)j) *(gfxReal)mSubdivs-1.5); 
462                                                         const int spk= (int)round( (ppos[2]+1.-(gfxReal)k) *(gfxReal)mSubdivs-1.5)-sdkOffset; // why -2?
463                                                         // 2d should be handled by solver. if(LBMDIM==2) { spk = 0; }
464
465                                                         gfxReal pfLen = p->getSize()*1.5*mPartSize;  // test, was 1.1
466                                                         const gfxReal minPfLen = subdfac*0.8;
467                                                         if(pfLen<minPfLen) pfLen = minPfLen;
468                                                         //errMsg("ISOPPP"," at "<<PRINT_IJK<<"  pp"<<ppos<<"  sp"<<PRINT_VEC(spi,spj,spk)<<" pflen"<<pfLen );
469                                                         //errMsg("ISOPPP"," subdfac="<<subdfac<<" size"<<p->getSize()<<" ps"<<mPartSize );
470                                                         const int icellpsize = (int)(1.*pfLen*(gfxReal)mSubdivs)+1;
471                                                         for(int swk=-icellpsize; swk<=icellpsize; swk++) {
472                                                                 if(spk+swk<         0) { continue; }
473                                                                 if(spk+swk>         1) { continue; } // */
474                                                         for(int swj=-icellpsize; swj<=icellpsize; swj++) {
475                                                                 if(spj+swj<         0) { continue; }
476                                                                 if(spj+swj>mSubdivs+0) { continue; } // */
477                                                         for(int swi=-icellpsize; swi<=icellpsize; swi++) {
478                                                                 if(spi+swi<         0) { continue; } 
479                                                                 if(spi+swi>mSubdivs+0) { continue; } // */
480                                                                 ntlVec3Gfx cellp = ntlVec3Gfx(
481                                                                                 (1.5+(gfxReal)(spi+swi))           *subdfac + (gfxReal)(i-1),
482                                                                                 (1.5+(gfxReal)(spj+swj))           *subdfac + (gfxReal)(j-1),
483                                                                                 (1.5+(gfxReal)(spk+swk)+sdkOffset) *subdfac + (gfxReal)(k-1)
484                                                                                 );
485                                                                 //if(swi==0 && swj==0 && swk==0) subdAr[spk][spj][spi] = 1.; // DEBUG
486                                                                 // clip domain boundaries again 
487                                                                 if(cellp[0]<1.) { continue; } 
488                                                                 if(cellp[1]<1.) { continue; } 
489                                                                 if(cellp[2]<1.) { continue; } 
490                                                                 if(cellp[0]>(gfxReal)mSizex-3.) { continue; } 
491                                                                 if(cellp[1]>(gfxReal)mSizey-3.) { continue; } 
492                                                                 if(cellp[2]>(gfxReal)mSizez-3.) { continue; } 
493                                                                 gfxReal len = norm(cellp-ppos);
494                                                                 gfxReal isoadd = 0.; 
495                                                                 const gfxReal baseIsoVal = mIsoValue*1.1;
496                                                                 if(len<pfLen) { 
497                                                                         isoadd = baseIsoVal*1.;
498                                                                 } else { 
499                                                                         // falloff linear with pfLen (kernel size=2pfLen
500                                                                         isoadd = baseIsoVal*(1. - (len-pfLen)/(pfLen)); 
501                                                                 }
502                                                                 if(isoadd<0.) { continue; }
503                                                                 //errMsg("ISOPPP"," at "<<PRINT_IJK<<" sp"<<PRINT_VEC(spi+swi,spj+swj,spk+swk)<<" cellp"<<cellp<<" pp"<<ppos << " l"<< len<< " add"<< isoadd);
504                                                                 const gfxReal arval = subdAr[spk+swk][spj+swj][spi+swi];
505                                                                 if(arval>1.) { continue; }
506                                                                 subdAr[spk+swk][spj+swj][spi+swi] = arval + isoadd;
507                                                         } } }
508
509                                                         p = p->getNext();
510                                                 }
511                                         } } } // poDist loops */
512
513                                         py = mStart[1]+(((double)j-0.5)*orgGsy)-gsy;
514                                         for(int sj=0;sj<mSubdivs;sj++) {
515                                                 py += gsy;
516                                                 px = mStart[0]+(((double)i-0.5)*orgGsx)-gsx;
517                                                 for(int si=0;si<mSubdivs;si++) {
518                                                         px += gsx;
519                                                         value[0] = subdAr[0+0][sj+0][si+0]; 
520                                                         value[1] = subdAr[0+0][sj+0][si+1]; 
521                                                         value[2] = subdAr[0+0][sj+1][si+1]; 
522                                                         value[3] = subdAr[0+0][sj+1][si+0]; 
523                                                         value[4] = subdAr[0+1][sj+0][si+0]; 
524                                                         value[5] = subdAr[0+1][sj+0][si+1]; 
525                                                         value[6] = subdAr[0+1][sj+1][si+1]; 
526                                                         value[7] = subdAr[0+1][sj+1][si+0]; 
527
528                                                         // check intersections of isosurface with edges, and calculate cubie index
529                                                         cubeIndex = 0;
530                                                         if (value[0] < mIsoValue) cubeIndex |= 1;
531                                                         if (value[1] < mIsoValue) cubeIndex |= 2; // with subdivs
532                                                         if (value[2] < mIsoValue) cubeIndex |= 4;
533                                                         if (value[3] < mIsoValue) cubeIndex |= 8;
534                                                         if (value[4] < mIsoValue) cubeIndex |= 16;
535                                                         if (value[5] < mIsoValue) cubeIndex |= 32; // with subdivs
536                                                         if (value[6] < mIsoValue) cubeIndex |= 64;
537                                                         if (value[7] < mIsoValue) cubeIndex |= 128;
538
539                                                         if (mcEdgeTable[cubeIndex] >  0) {
540
541                                                         // where to look up if this point already exists
542                                                         const int edgek = 0;
543                                                         const int baseIn = EDGEAR_INDEX( i+0, j+0, edgek+0, si,sj);
544                                                         eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
545                                                         eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
546                                                         eVert[ 2] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ];
547                                                         eVert[ 3] = &mpEdgeVerticesY[ baseIn ];                             
548                                                                                                                                                                                                                                                                                                                                         
549                                                         eVert[ 4] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ];
550                                                         eVert[ 5] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+1,sj+0) ]; // with subdivs
551                                                         eVert[ 6] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+1) ];
552                                                         eVert[ 7] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ];
553                                                                                                                                                                                                                                                                                                                                         
554                                                         eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];                             
555                                                         eVert[ 9] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+0) ]; // with subdivs
556                                                         eVert[10] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+1) ];
557                                                         eVert[11] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ];
558
559                                                         // grid positions
560                                                         pos[0] = ntlVec3Gfx(px    ,py    ,pz);
561                                                         pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
562                                                         pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz); // with subdivs
563                                                         pos[3] = ntlVec3Gfx(px    ,py+gsy,pz);
564                                                         pos[4] = ntlVec3Gfx(px    ,py    ,pz+gsz);
565                                                         pos[5] = ntlVec3Gfx(px+gsx,py    ,pz+gsz);
566                                                         pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz); // with subdivs
567                                                         pos[7] = ntlVec3Gfx(px    ,py+gsy,pz+gsz);
568
569                                                         // check all edges
570                                                         for(int e=0;e<12;e++) {
571                                                                 if (mcEdgeTable[cubeIndex] & (1<<e)) {
572                                                                         // is the vertex already calculated?
573                                                                         if(*eVert[ e ] < 0) {
574                                                                                 // interpolate edge
575                                                                                 const int e1 = mcEdges[e*2  ];
576                                                                                 const int e2 = mcEdges[e*2+1];
577                                                                                 const ntlVec3Gfx p1 = pos[ e1  ];   // scalar field pos 1
578                                                                                 const ntlVec3Gfx p2 = pos[ e2  ];   // scalar field pos 2
579                                                                                 const float valp1  = value[ e1  ];  // scalar field val 1
580                                                                                 const float valp2  = value[ e2  ];  // scalar field val 2
581                                                                                 const float mu = (mIsoValue - valp1) / (valp2 - valp1);
582
583                                                                                 // init isolevel vertex
584                                                                                 ilv.v = p1 + (p2-p1)*mu; // with subdivs
585 #pragma omp critical 
586                                                                                 {
587                                                                                 mPoints.push_back( ilv );
588                                                                                 triIndices[e] = (mPoints.size()-1);
589                                                                                 }
590                                                                                 // store vertex 
591                                                                                 *eVert[ e ] = triIndices[e]; 
592                                                                         }       else {
593                                                                                 // retrieve  from vert array
594                                                                                 triIndices[e] = *eVert[ e ];
595                                                                         }
596                                                                 } // along all edges 
597                                                         }
598                                                         // removed cutoff treatment...
599 #pragma omp critical 
600                                                         {
601                                                         // Create the triangles... 
602                                                         for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
603                                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
604                                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); // with subdivs
605                                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
606                                                                 //errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() );
607                                                         }
608                                                         }
609                                                         } // triangles in edge table?
610                                                         
611                                                 }//si
612                                         }// sj
613
614                                 }//i
615                         }// j
616 #pragma omp barrier
617                         // copy edge arrays
618                         for(int j=0;j<(mSizey-0)*mSubdivs;j++) 
619                                 for(int i=0;i<(mSizex-0)*mSubdivs;i++) {
620                                         //int edgek = 0;
621                                         const int dst = EDGEAR_INDEX( 0, 0, 0, i,j);
622                                         const int src = EDGEAR_INDEX( 0, 0, 1, i,j);
623                                         mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ];
624                                         mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; // with subdivs
625                                         mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ];
626                                         mpEdgeVerticesX[ src ]=-1;
627                                         mpEdgeVerticesY[ src ]=-1; // with subdivs
628                                         mpEdgeVerticesZ[ src ]=-1;
629                                 }
630                         // */
631
632                 } // ok, k subdiv loop
633
634                 //delete [] subdAr;
635                 delete [] arppnt;
636                 computeNormals();
637         } // with subdivs
638
639         // perform smoothing
640         float smoSubdfac = 1.;
641         if(mSubdivs>0) {
642                 //smoSubdfac = 1./(float)(mSubdivs);
643                 smoSubdfac = pow(0.55,(double)mSubdivs); // slightly stronger
644         }
645         if(mSmoothSurface>0. || mSmoothNormals>0.) debMsgStd("IsoSurface::triangulate",DM_MSG,"Smoothing...",10);
646         if(mSmoothSurface>0.0) { 
647                 smoothSurface(mSmoothSurface*smoSubdfac, (mSmoothNormals<=0.0) ); 
648         }
649         if(mSmoothNormals>0.0) { 
650                 smoothNormals(mSmoothNormals*smoSubdfac); 
651         }
652
653         myTime_t tritimeend = getTime(); 
654         debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< getTimeString(tritimeend-tritimestart)<<", S("<<mSmoothSurface<<","<<mSmoothNormals<<"),"<<
655                         " verts:"<<mPoints.size()<<" tris:"<<(mIndices.size()/3)<<" subdivs:"<<mSubdivs
656                  , 10 );
657         if(mpIsoParts) debMsgStd("IsoSurface::triangulate",DM_MSG,"parts:"<<mpIsoParts->getNumParticles(), 10);
658 }
659
660
661         
662
663
664 /******************************************************************************
665  * Get triangles for rendering
666  *****************************************************************************/
667 void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles, 
668                                                                                                          vector<ntlVec3Gfx> *vertices, 
669                                                                                                          vector<ntlVec3Gfx> *normals, int objectId )
670 {
671         if(!mInitDone) {
672                 debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10);
673                 return;
674         }
675         t = 0.;
676         //return; // DEBUG
677
678   /* triangulate field */
679   triangulate();
680         //errMsg("TRIS"," "<<mIndices.size() );
681
682         // new output with vertice reuse
683         int iniVertIndex = (*vertices).size();
684         int iniNormIndex = (*normals).size();
685         if(iniVertIndex != iniNormIndex) {
686                 errFatal("getTriangles Error","For '"<<mName<<"': Vertices and normal array sizes to not match!!!",SIMWORLD_GENERICERROR);
687                 return; 
688         }
689         //errMsg("NM"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
690         //errMsg("NM"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
691
692   for(int i=0;i<(int)mPoints.size();i++) {
693                 vertices->push_back( mPoints[i].v );
694         }
695   for(int i=0;i<(int)mPoints.size();i++) {
696                 normals->push_back( mPoints[i].n );
697         }
698
699         //errMsg("N2"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
700         //errMsg("N2"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
701
702   for(int i=0;i<(int)mIndices.size();i+=3) {
703                 const int smooth = 1;
704     int t1 = mIndices[i];
705     int t2 = mIndices[i+1];
706                 int t3 = mIndices[i+2];
707                 //errMsg("NM"," tri"<<t1<<" "<<t2<<" "<<t3 );
708
709                 ntlTriangle tri;
710
711                 tri.getPoints()[0] = t1+iniVertIndex;
712                 tri.getPoints()[1] = t2+iniVertIndex;
713                 tri.getPoints()[2] = t3+iniVertIndex;
714
715                 /* init flags */
716                 int flag = 0; 
717                 if(getVisible()){ flag |= TRI_GEOMETRY; }
718                 if(getCastShadows() ) { 
719                         flag |= TRI_CASTSHADOWS; } 
720
721                 /* init geo init id */
722                 int geoiId = getGeoInitId(); 
723                 if(geoiId > 0) { 
724                         flag |= (1<< (geoiId+4)); 
725                         flag |= mGeoInitType; 
726                 } 
727
728                 tri.setFlags( flag );
729
730                 /* triangle normal missing */
731                 tri.setNormal( ntlVec3Gfx(0.0) );
732                 tri.setSmoothNormals( smooth );
733                 tri.setObjectId( objectId );
734                 triangles->push_back( tri ); 
735         }
736         //errMsg("N3"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
737         return;
738 }
739
740
741
742 inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) {
743         // WARNING - this requires a security boundary layer... 
744         ntlVec3Gfx ret(0.0);
745         ret[0] = *getData(i-1,j  ,k  ) - 
746                  *getData(i+1,j  ,k  );
747         ret[1] = *getData(i  ,j-1,k  ) - 
748                  *getData(i  ,j+1,k  );
749         ret[2] = *getData(i  ,j  ,k-1  ) - 
750                  *getData(i  ,j  ,k+1  );
751         return ret;
752 }
753
754
755
756
757 /******************************************************************************
758  * 
759  * Surface improvement, inspired by trimesh2 library
760  * (http://www.cs.princeton.edu/gfx/proj/trimesh2/)
761  * 
762  *****************************************************************************/
763
764 void IsoSurface::setSmoothRad(float radi1, float radi2, ntlVec3Gfx mscc) {
765         mSCrad1 = radi1*radi1;
766         mSCrad2 = radi2*radi2;
767         mSCcenter = mscc;
768 }
769
770 // compute normals for all generated triangles
771 void IsoSurface::computeNormals() {
772   for(int i=0;i<(int)mPoints.size();i++) {
773                 mPoints[i].n = ntlVec3Gfx(0.);
774         }
775
776   for(int i=0;i<(int)mIndices.size();i+=3) {
777     const int t1 = mIndices[i];
778     const int t2 = mIndices[i+1];
779                 const int t3 = mIndices[i+2];
780                 const ntlVec3Gfx p1 = mPoints[t1].v;
781                 const ntlVec3Gfx p2 = mPoints[t2].v;
782                 const ntlVec3Gfx p3 = mPoints[t3].v;
783                 const ntlVec3Gfx n1=p1-p2;
784                 const ntlVec3Gfx n2=p2-p3;
785                 const ntlVec3Gfx n3=p3-p1;
786                 const gfxReal len1 = normNoSqrt(n1);
787                 const gfxReal len2 = normNoSqrt(n2);
788                 const gfxReal len3 = normNoSqrt(n3);
789                 const ntlVec3Gfx norm = cross(n1,n2);
790                 mPoints[t1].n += norm * (1./(len1*len3));
791                 mPoints[t2].n += norm * (1./(len1*len2));
792                 mPoints[t3].n += norm * (1./(len2*len3));
793         }
794
795   for(int i=0;i<(int)mPoints.size();i++) {
796                 normalize(mPoints[i].n);
797         }
798 }
799
800 // Diffuse a vector field at 1 vertex, weighted by
801 // a gaussian of width 1/sqrt(invsigma2)
802 bool IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int src, float invsigma2, ntlVec3Gfx &target)
803 {
804         if((neighbors[src].size()<1) || (pointareas[src]<=0.0)) return 0;
805         const ntlVec3Gfx srcp = mPoints[src].v;
806         const ntlVec3Gfx srcn = mPoints[src].n;
807         if(mSCrad1>0.0 && mSCrad2>0.0) {
808                 ntlVec3Gfx dp = mSCcenter-srcp; dp[2] = 0.0; // only xy-plane
809                 float rd = normNoSqrt(dp);
810                 if(rd > mSCrad2) {
811                         return 0;
812                 } else if(rd > mSCrad1) {
813                         // optimize?
814                         float org = 1.0/sqrt(invsigma2);
815                         org *= (1.0- (rd-mSCrad1) / (mSCrad2-mSCrad1));
816                         invsigma2 = 1.0/(org*org);
817                         //errMsg("TRi","p"<<srcp<<" rd:"<<rd<<" r1:"<<mSCrad1<<" r2:"<<mSCrad2<<" org:"<<org<<" is:"<<invsigma2);
818                 } else {
819                 }
820         }
821         target = ntlVec3Gfx(0.0);
822         target += *(field+pointerScale*src) *pointareas[src];
823         float smstrSum = pointareas[src];
824
825         int flag = mFlagCnt; 
826         mFlagCnt++;
827         flags[src] = flag;
828         mDboundary = neighbors[src];
829         while (!mDboundary.empty()) {
830                 const int bbn = mDboundary.back();
831                 mDboundary.pop_back();
832                 if(flags[bbn]==flag) continue;
833                 flags[bbn] = flag;
834
835                 // normal check
836                 const float nvdot = dot(srcn, mPoints[bbn].n); // faster than before d2 calc?
837                 if(nvdot <= 0.0f) continue;
838
839                 // gaussian weight of width 1/sqrt(invsigma2)
840                 const float d2 = invsigma2 * normNoSqrt(mPoints[bbn].v - srcp);
841                 if(d2 >= 9.0f) continue;
842
843                 // aggressive smoothing factor
844                 float smstr = nvdot * pointareas[bbn];
845                 // Accumulate weight times field at neighbor
846                 target += *(field+pointerScale*bbn)*smstr;
847                 smstrSum += smstr;
848
849                 for(int i = 0; i < (int)neighbors[bbn].size(); i++) {
850                         const int nn = neighbors[bbn][i];
851                         if (flags[nn] == flag) continue;
852                         mDboundary.push_back(nn);
853                 }
854         }
855         target /= smstrSum;
856         return 1;
857 }
858
859         
860 // perform smoothing of the surface (and possible normals)
861 void IsoSurface::smoothSurface(float sigma, bool normSmooth)
862 {
863         int nv = mPoints.size();
864         if ((int)flags.size() != nv) flags.resize(nv);
865         int nf = mIndices.size()/3;
866
867         { // need neighbors
868                 vector<int> numneighbors(mPoints.size());
869                 int i;
870                 for (i = 0; i < (int)mIndices.size()/3; i++) {
871                         numneighbors[mIndices[i*3+0]]++;
872                         numneighbors[mIndices[i*3+1]]++;
873                         numneighbors[mIndices[i*3+2]]++;
874                 }
875
876                 neighbors.clear();
877                 neighbors.resize(mPoints.size());
878                 for (i = 0; i < (int)mPoints.size(); i++) {
879                         neighbors[i].clear();
880                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
881                 }
882
883                 for (i = 0; i < (int)mIndices.size()/3; i++) {
884                         for (int j = 0; j < 3; j++) {
885                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
886                                 int n1 =  mIndices[i*3+((j+1)%3)];
887                                 int n2 =  mIndices[i*3+((j+2)%3)];
888                                 if (std::find(me.begin(), me.end(), n1) == me.end())
889                                         me.push_back(n1);
890                                 if (std::find(me.begin(), me.end(), n2) == me.end())
891                                         me.push_back(n2);
892                         }
893                 }
894         } // need neighbor
895
896         { // need pointarea
897                 pointareas.clear();
898                 pointareas.resize(nv);
899                 cornerareas.clear();
900                 cornerareas.resize(nf);
901
902                 for (int i = 0; i < nf; i++) {
903                         // Edges
904                         ntlVec3Gfx e[3] = { 
905                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
906                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
907                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
908
909                         // Compute corner weights
910                         float area = 0.5f * norm( cross(e[0], e[1]));
911                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
912                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
913                                         l2[1] * (l2[2] + l2[0] - l2[1]),
914                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
915                         if (ew[0] <= 0.0f) {
916                                 cornerareas[i][1] = -0.25f * l2[2] * area /
917                                                                 dot(e[0] , e[2]);
918                                 cornerareas[i][2] = -0.25f * l2[1] * area /
919                                                                 dot(e[0] , e[1]);
920                                 cornerareas[i][0] = area - cornerareas[i][1] -
921                                                                 cornerareas[i][2];
922                         } else if (ew[1] <= 0.0f) {
923                                 cornerareas[i][2] = -0.25f * l2[0] * area /
924                                                                 dot(e[1] , e[0]);
925                                 cornerareas[i][0] = -0.25f * l2[2] * area /
926                                                                 dot(e[1] , e[2]);
927                                 cornerareas[i][1] = area - cornerareas[i][2] -
928                                                                 cornerareas[i][0];
929                         } else if (ew[2] <= 0.0f) {
930                                 cornerareas[i][0] = -0.25f * l2[1] * area /
931                                                                 dot(e[2] , e[1]);
932                                 cornerareas[i][1] = -0.25f * l2[0] * area /
933                                                                 dot(e[2] , e[0]);
934                                 cornerareas[i][2] = area - cornerareas[i][0] -
935                                                                 cornerareas[i][1];
936                         } else {
937                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
938                                 for (int j = 0; j < 3; j++)
939                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
940                                                                                          ew[(j+2)%3]);
941                         }
942
943                         // NT important, check this...
944 #ifndef WIN32
945                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
946                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
947                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
948 #else // WIN32
949                         // FIXME check as well...
950                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
951                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
952                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
953 #endif // WIN32
954
955                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
956                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
957                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
958                 }
959
960         } // need pointarea
961         // */
962
963         float invsigma2 = 1.0f / (sigma*sigma);
964
965         vector<ntlVec3Gfx> dflt(nv);
966         for (int i = 0; i < nv; i++) {
967                 if(diffuseVertexField( &mPoints[0].v, 2,
968                                    i, invsigma2, dflt[i])) {
969                         // Just keep the displacement
970                         dflt[i] -= mPoints[i].v;
971                 } else { dflt[i] = 0.0; } //?mPoints[i].v; }
972         }
973
974         // Slightly better small-neighborhood approximation
975         for (int i = 0; i < nf; i++) {
976                 ntlVec3Gfx c = mPoints[mIndices[i*3+0]].v +
977                           mPoints[mIndices[i*3+1]].v +
978                           mPoints[mIndices[i*3+2]].v;
979                 c /= 3.0f;
980                 for (int j = 0; j < 3; j++) {
981                         int v = mIndices[i*3+j];
982                         ntlVec3Gfx d =(c - mPoints[v].v) * 0.5f;
983                         dflt[v] += d * (cornerareas[i][j] /
984                                    pointareas[mIndices[i*3+j]] *
985                                    exp(-0.5f * invsigma2 * normNoSqrt(d)) );
986                 }
987         }
988
989         // Filter displacement field
990         vector<ntlVec3Gfx> dflt2(nv);
991         for (int i = 0; i < nv; i++) {
992                 if(diffuseVertexField( &dflt[0], 1,
993                                    i, invsigma2, dflt2[i])) { }
994                 else { /*mPoints[i].v=0.0;*/ dflt2[i] = 0.0; }//dflt2[i]; }
995         }
996
997         // Update vertex positions
998         for (int i = 0; i < nv; i++) {
999                 mPoints[i].v += dflt[i] - dflt2[i]; // second Laplacian
1000         }
1001
1002         // when normals smoothing off, this cleans up quite well
1003         // costs ca. 50% additional time though
1004         float nsFac = 1.5f;
1005         if(normSmooth) { float ninvsigma2 = 1.0f / (nsFac*nsFac*sigma*sigma);
1006                 for (int i = 0; i < nv; i++) {
1007                         if( diffuseVertexField( &mPoints[0].n, 2, i, ninvsigma2, dflt[i]) ) {
1008                                 normalize(dflt[i]);
1009                         } else {
1010                                 dflt[i] = mPoints[i].n;
1011                         }
1012                 }
1013                 for (int i = 0; i < nv; i++) {
1014                         mPoints[i].n = dflt[i];
1015                 } 
1016         } // smoothNormals copy */
1017
1018         //errMsg("SMSURF","done v:"<<sigma); // DEBUG
1019 }
1020
1021 // only smoothen the normals
1022 void IsoSurface::smoothNormals(float sigma) {
1023         // reuse from smoothSurface
1024         if(neighbors.size() != mPoints.size()) { 
1025                 // need neighbor
1026                 vector<int> numneighbors(mPoints.size());
1027                 int i;
1028                 for (i = 0; i < (int)mIndices.size()/3; i++) {
1029                         numneighbors[mIndices[i*3+0]]++;
1030                         numneighbors[mIndices[i*3+1]]++;
1031                         numneighbors[mIndices[i*3+2]]++;
1032                 }
1033
1034                 neighbors.clear();
1035                 neighbors.resize(mPoints.size());
1036                 for (i = 0; i < (int)mPoints.size(); i++) {
1037                         neighbors[i].clear();
1038                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
1039                 }
1040
1041                 for (i = 0; i < (int)mIndices.size()/3; i++) {
1042                         for (int j = 0; j < 3; j++) {
1043                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
1044                                 int n1 =  mIndices[i*3+((j+1)%3)];
1045                                 int n2 =  mIndices[i*3+((j+2)%3)];
1046                                 if (std::find(me.begin(), me.end(), n1) == me.end())
1047                                         me.push_back(n1);
1048                                 if (std::find(me.begin(), me.end(), n2) == me.end())
1049                                         me.push_back(n2);
1050                         }
1051                 }
1052         } // need neighbor
1053
1054         { // need pointarea
1055                 int nf = mIndices.size()/3, nv = mPoints.size();
1056                 pointareas.clear();
1057                 pointareas.resize(nv);
1058                 cornerareas.clear();
1059                 cornerareas.resize(nf);
1060
1061                 for (int i = 0; i < nf; i++) {
1062                         // Edges
1063                         ntlVec3Gfx e[3] = { 
1064                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
1065                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
1066                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
1067
1068                         // Compute corner weights
1069                         float area = 0.5f * norm( cross(e[0], e[1]));
1070                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
1071                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
1072                                         l2[1] * (l2[2] + l2[0] - l2[1]),
1073                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
1074                         if (ew[0] <= 0.0f) {
1075                                 cornerareas[i][1] = -0.25f * l2[2] * area /
1076                                                                 dot(e[0] , e[2]);
1077                                 cornerareas[i][2] = -0.25f * l2[1] * area /
1078                                                                 dot(e[0] , e[1]);
1079                                 cornerareas[i][0] = area - cornerareas[i][1] -
1080                                                                 cornerareas[i][2];
1081                         } else if (ew[1] <= 0.0f) {
1082                                 cornerareas[i][2] = -0.25f * l2[0] * area /
1083                                                                 dot(e[1] , e[0]);
1084                                 cornerareas[i][0] = -0.25f * l2[2] * area /
1085                                                                 dot(e[1] , e[2]);
1086                                 cornerareas[i][1] = area - cornerareas[i][2] -
1087                                                                 cornerareas[i][0];
1088                         } else if (ew[2] <= 0.0f) {
1089                                 cornerareas[i][0] = -0.25f * l2[1] * area /
1090                                                                 dot(e[2] , e[1]);
1091                                 cornerareas[i][1] = -0.25f * l2[0] * area /
1092                                                                 dot(e[2] , e[0]);
1093                                 cornerareas[i][2] = area - cornerareas[i][0] -
1094                                                                 cornerareas[i][1];
1095                         } else {
1096                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
1097                                 for (int j = 0; j < 3; j++)
1098                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
1099                                                                                          ew[(j+2)%3]);
1100                         }
1101
1102                         // NT important, check this...
1103 #ifndef WIN32
1104                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
1105                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
1106                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
1107 #else // WIN32
1108                         // FIXME check as well...
1109                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
1110                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
1111                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
1112 #endif // WIN32
1113
1114                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
1115                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
1116                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
1117                 }
1118
1119         } // need pointarea
1120
1121         int nv = mPoints.size();
1122         if ((int)flags.size() != nv) flags.resize(nv);
1123         float invsigma2 = 1.0f / (sigma*sigma);
1124
1125         vector<ntlVec3Gfx> nflt(nv);
1126         for (int i = 0; i < nv; i++) {
1127                 if(diffuseVertexField( &mPoints[0].n, 2, i, invsigma2, nflt[i])) {
1128                         normalize(nflt[i]);
1129                 } else { nflt[i]=mPoints[i].n; }
1130         }
1131
1132         // copy results
1133         for (int i = 0; i < nv; i++) { mPoints[i].n = nflt[i]; }
1134 }
1135
1136