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