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