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