- removed some unecessary files & code
[blender.git] / intern / elbeem / intern / isosurface.cpp
1 /******************************************************************************
2  *
3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4  * Copyright 2005 Nils Thuerey
5  *
6  * Marching Cubes surface mesh generation
7  *
8  *****************************************************************************/
9
10 #include "isosurface.h"
11 #include "mcubes_tables.h"
12 #include "ntl_scene.h"
13
14 #include <algorithm>
15 #include <stdio.h>
16
17
18
19 /******************************************************************************
20  * Constructor
21  *****************************************************************************/
22 IsoSurface::IsoSurface(double iso, double blend) :
23         ntlGeometryObject(),
24         mSizex(-1), mSizey(-1), mSizez(-1),
25   mIsoValue( iso ), 
26         mBlendVal( blend ),
27         mPoints(), 
28         mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL),
29         mIndices(),
30
31   mStart(0.0), mEnd(0.0), mDomainExtent(0.0),
32   mInitDone(false),
33         mSmoothSurface(0.0), mSmoothNormals(0.0),
34         mAcrossEdge(), mAdjacentFaces()
35 {
36 }
37
38
39 /******************************************************************************
40  * The real init...
41  *****************************************************************************/
42 void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent) 
43 {
44         // init solver and size
45         mSizex = setx;
46         mSizey = sety;
47         if(setz == 1) {// 2D, create thin 2D surface
48                 setz = 5; 
49         }
50         mSizez = setz;
51         mDomainExtent = extent;
52
53         /* check triangulation size (for raytraing) */
54         if( ( mStart[0] >= mEnd[0] ) && ( mStart[1] >= mEnd[1] ) && ( mStart[2] >= mEnd[2] ) ){
55                 // extent was not set, use normalized one from parametrizer
56                 mStart = ntlVec3Gfx(0.0) - extent*0.5;
57                 mEnd = ntlVec3Gfx(0.0) + extent*0.5;
58         }
59
60   // init 
61         mIndices.clear();
62   mPoints.clear();
63
64         int nodes = mSizez*mSizey*mSizex;
65   mpData = new float[nodes];
66   for(int i=0;i<nodes;i++) { mpData[i] = 0.0; }
67
68   // allocate edge arrays  (last slices are never used...)
69   mpEdgeVerticesX = new int[nodes];
70   mpEdgeVerticesY = new int[nodes];
71   mpEdgeVerticesZ = new int[nodes];
72   for(int i=0;i<nodes;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
73   
74         // marching cubes are ready 
75         mInitDone = true;
76         unsigned long int memCnt = (3*sizeof(int)*nodes + sizeof(float)*nodes);
77         double memd = memCnt;
78         char *sizeStr = "";
79         const double sfac = 1000.0;
80         if(memd>sfac){ memd /= sfac; sizeStr="KB"; }
81         if(memd>sfac){ memd /= sfac; sizeStr="MB"; }
82         if(memd>sfac){ memd /= sfac; sizeStr="GB"; }
83         if(memd>sfac){ memd /= sfac; sizeStr="TB"; }
84
85         debMsgStd("IsoSurface::initializeIsosurface",DM_MSG,"Inited "<<PRINT_VEC(setx,sety,setz)<<" alloced:"<< memd<<" "<<sizeStr<<"." ,10);
86 }
87
88
89
90
91 /******************************************************************************
92  * Destructor
93  *****************************************************************************/
94 IsoSurface::~IsoSurface( void )
95 {
96
97         if(mpEdgeVerticesX) delete [] mpEdgeVerticesX;
98         if(mpEdgeVerticesY) delete [] mpEdgeVerticesY;
99         if(mpEdgeVerticesZ) delete [] mpEdgeVerticesZ;
100
101 }
102
103
104
105
106 /******************************************************************************
107  * triangulate the scalar field given by pointer
108  *****************************************************************************/
109 void IsoSurface::triangulate( void )
110 {
111   double gsx,gsy,gsz; // grid spacing in x,y,z direction
112   double px,py,pz;    // current position in grid in x,y,z direction
113   IsoLevelCube cubie;    // struct for a small subcube
114         myTime_t tritimestart = getTime(); 
115
116         if(!mpData) {
117                 errFatal("IsoSurface::triangulate","no LBM object, and no scalar field...!",SIMWORLD_INITERROR);
118                 return;
119         }
120
121   // get grid spacing
122   gsx = (mEnd[0]-mStart[0])/(double)(mSizex-1.0);
123   gsy = (mEnd[1]-mStart[1])/(double)(mSizey-1.0);
124   gsz = (mEnd[2]-mStart[2])/(double)(mSizez-1.0);
125
126   // clean up previous frame
127         mIndices.clear();
128         mPoints.clear();
129
130         // reset edge vertices
131   for(int i=0;i<(mSizex*mSizey*mSizez);i++) {
132                 mpEdgeVerticesX[i] = -1;
133                 mpEdgeVerticesY[i] = -1;
134                 mpEdgeVerticesZ[i] = -1;
135         }
136
137         ntlVec3Gfx pos[8];
138         float value[8];
139         int cubeIndex;      // index entry of the cube 
140         int triIndices[12]; // vertex indices 
141         int *eVert[12];
142         IsoLevelVertex ilv;
143
144         // edges between which points?
145         const int mcEdges[24] = { 
146                 0,1,  1,2,  2,3,  3,0,
147                 4,5,  5,6,  6,7,  7,4,
148                 0,4,  1,5,  2,6,  3,7 };
149
150         const int cubieOffsetX[8] = {
151                 0,1,1,0,  0,1,1,0 };
152         const int cubieOffsetY[8] = {
153                 0,0,1,1,  0,0,1,1 };
154         const int cubieOffsetZ[8] = {
155                 0,0,0,0,  1,1,1,1 };
156
157   // let the cubes march 
158         pz = mStart[2]-gsz;
159         for(int k=0;k<(mSizez-2);k++) {
160                 pz += gsz;
161                 py = mStart[1]-gsy;
162                 for(int j=0;j<(mSizey-2);j++) {
163       py += gsy;
164                         px = mStart[0]-gsx;
165                         for(int i=0;i<(mSizex-2);i++) {
166                         px += gsx;
167                                 int baseIn = ISOLEVEL_INDEX( i+0, j+0, k+0);
168
169                                 value[0] = *getData(i  ,j  ,k  );
170                                 value[1] = *getData(i+1,j  ,k  );
171                                 value[2] = *getData(i+1,j+1,k  );
172                                 value[3] = *getData(i  ,j+1,k  );
173                                 value[4] = *getData(i  ,j  ,k+1);
174                                 value[5] = *getData(i+1,j  ,k+1);
175                                 value[6] = *getData(i+1,j+1,k+1);
176                                 value[7] = *getData(i  ,j+1,k+1);
177                                 //errMsg("ISOT2D"," at "<<PRINT_IJK<<" "
178                                                 //<<" v0="<<value[0] <<" v1="<<value[1] <<" v2="<<value[2] <<" v3="<<value[3]
179                                                 //<<" v4="<<value[4] <<" v5="<<value[5] <<" v6="<<value[6] <<" v7="<<value[7] );
180
181                                 // check intersections of isosurface with edges, and calculate cubie index
182                                 cubeIndex = 0;
183                                 if (value[0] < mIsoValue) cubeIndex |= 1;
184                                 if (value[1] < mIsoValue) cubeIndex |= 2;
185                                 if (value[2] < mIsoValue) cubeIndex |= 4;
186                                 if (value[3] < mIsoValue) cubeIndex |= 8;
187                                 if (value[4] < mIsoValue) cubeIndex |= 16;
188                                 if (value[5] < mIsoValue) cubeIndex |= 32;
189                                 if (value[6] < mIsoValue) cubeIndex |= 64;
190                                 if (value[7] < mIsoValue) cubeIndex |= 128;
191
192                                 // No triangles to generate?
193                                 if (mcEdgeTable[cubeIndex] == 0) {
194                                         continue;
195                                 }
196
197                                 // where to look up if this point already exists
198                                 eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
199                                 eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
200                                 eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, k+0) ];
201                                 eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
202
203                                 eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, k+1) ];
204                                 eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, k+1) ];
205                                 eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, k+1) ];
206                                 eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, k+1) ];
207
208                                 eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
209                                 eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, k+0) ];
210                                 eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, k+0) ];
211                                 eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, k+0) ];
212
213                                 // grid positions
214                                 pos[0] = ntlVec3Gfx(px    ,py    ,pz);
215                                 pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
216                                 pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
217                                 pos[3] = ntlVec3Gfx(px    ,py+gsy,pz);
218                                 pos[4] = ntlVec3Gfx(px    ,py    ,pz+gsz);
219                                 pos[5] = ntlVec3Gfx(px+gsx,py    ,pz+gsz);
220                                 pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz);
221                                 pos[7] = ntlVec3Gfx(px    ,py+gsy,pz+gsz);
222
223                                 // check all edges
224                                 for(int e=0;e<12;e++) {
225
226                                         if (mcEdgeTable[cubeIndex] & (1<<e)) {
227                                                 // is the vertex already calculated?
228                                                 if(*eVert[ e ] < 0) {
229                                                         // interpolate edge
230                                                         const int e1 = mcEdges[e*2  ];
231                                                         const int e2 = mcEdges[e*2+1];
232                                                         const ntlVec3Gfx p1 = pos[ e1  ];    // scalar field pos 1
233                                                         const ntlVec3Gfx p2 = pos[ e2  ];    // scalar field pos 2
234                                                         const float valp1  = value[ e1  ];  // scalar field val 1
235                                                         const float valp2  = value[ e2  ];  // scalar field val 2
236                                                         //double mu;            // interpolation value
237                                                         //ntlVec3Gfx p;         // new point
238
239                                                         // choose if point should be calculated by interpolation,
240                                                         // or "Karolin" method 
241
242                                                         //double deltaVal = ABS(valp2-valp1);
243                                                         //if(deltaVal >-10.0) {
244                                                         // standard interpolation
245                                                         //vertInfo[e].type = 0; 
246                                                         /*if (ABS(mIsoValue-valp1) < 0.000000001) {
247                                                                 mu = 0.0;
248                                                         } else {
249                                                                 if (ABS(mIsoValue-valp2) < 0.000000001) {
250                                                                         mu = 1.0;
251                                                                 } else {
252                                                                         mu = (mIsoValue - valp1) / (valp2 - valp1);
253                                                                 }
254                                                         } */
255                                                         /*} else {
256                                                                 errorOut(" ? ");
257                                                         // use fill grade (=karo) 
258                                                         vertInfo[e].type = 1; 
259                                                         if(valp1 < valp2) { mu = 1.0- (valp1 + valp2 - mIsoValue);
260                                                         } else { mu = 0.0+ (valp1 + valp2 - mIsoValue); } 
261                                                         } */
262
263                                                         //const float mu = (mIsoValue - valp1) / (valp2 - valp1);
264                                                         float mu;
265                                                         if(valp1 < valp2) {
266                                                                 mu = 1.0 - 1.0*(valp1 + valp2 - mIsoValue);
267                                                         } else {
268                                                                 mu = 0.0 + 1.0*(valp1 + valp2 - mIsoValue);
269                                                         }
270
271                                                         float isov2 = mIsoValue;
272                                                         isov2 = (valp1+valp2)*0.5;
273                                                         mu = (isov2 - valp1) / (valp2 - valp1);
274                                                         mu = (isov2) / (valp2 - valp1);
275
276                                                         mu = (mIsoValue - valp1) / (valp2 - valp1);
277                                                         //mu *= mu;
278
279
280
281                                                         // init isolevel vertex
282                                                         ilv.v = p1 + (p2-p1)*mu;
283                                                         ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) +
284                                                                                         getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * (    mu) ;
285                                                         mPoints.push_back( ilv );
286
287                                                         triIndices[e] = (mPoints.size()-1);
288                                                         // store vertex 
289                                                         *eVert[ e ] = triIndices[e];
290                                                 }       else {
291                                                         // retrieve  from vert array
292                                                         triIndices[e] = *eVert[ e ];
293                                                 }
294                                         } // along all edges 
295
296                                 }
297
298
299                                 // Create the triangles... 
300                                 for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
301                                         //errMsg("MC","tri "<<mIndices.size() <<" "<< triIndices[ mcTriTable[cubeIndex][e+0] ]<<" "<< triIndices[ mcTriTable[cubeIndex][e+1] ]<<" "<< triIndices[ mcTriTable[cubeIndex][e+2] ] );
302                                         mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
303                                         mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] );
304                                         mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
305                                 }
306
307                                 
308       }
309     }
310
311   } // k
312   
313
314   // precalculate normals using an approximation of the scalar field gradient 
315         for(int ni=0;ni<(int)mPoints.size();ni++) {
316                 // use triangle normals?, this seems better for File-IsoSurf
317                 normalize( mPoints[ni].n );
318         }
319
320         if(mSmoothSurface>0.0) {
321                 // not needed for post normal smoothing? 
322                 // if(mSmoothNormals<=0.0) { smoothNormals(mSmoothSurface*0.5); }
323                 smoothSurface(mSmoothSurface);
324         }
325         if(mSmoothNormals>0.0) {
326                 smoothNormals(mSmoothNormals);
327         }
328         myTime_t tritimeend = getTime(); 
329         debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< ((tritimeend-tritimestart)/(double)1000.0)<<"s, ss="<<mSmoothSurface<<" sm="<<mSmoothNormals , 10 );
330 }
331
332
333         
334
335
336 /******************************************************************************
337  * Get triangles for rendering
338  *****************************************************************************/
339 void IsoSurface::getTriangles( vector<ntlTriangle> *triangles, 
340                                                                                                          vector<ntlVec3Gfx> *vertices, 
341                                                                                                          vector<ntlVec3Gfx> *normals, int objectId )
342 {
343         if(!mInitDone) {
344                 debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10);
345                 return;
346         }
347
348   /* triangulate field */
349   triangulate();
350         //errMsg("TRIS"," "<<mIndices.size() );
351
352         // new output with vertice reuse
353         int iniVertIndex = (*vertices).size();
354         int iniNormIndex = (*normals).size();
355         if(iniVertIndex != iniNormIndex) {
356                 errFatal("getTriangles Error","For '"<<mName<<"': Vertices and normal array sizes to not match!!!",SIMWORLD_GENERICERROR);
357                 return; 
358         }
359         //errMsg("NM"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
360         //errMsg("NM"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
361
362   for(int i=0;i<(int)mPoints.size();i++) {
363                 vertices->push_back( mPoints[i].v );
364         }
365   for(int i=0;i<(int)mPoints.size();i++) {
366                 normals->push_back( mPoints[i].n );
367         }
368
369         //errMsg("N2"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
370         //errMsg("N2"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
371
372   for(int i=0;i<(int)mIndices.size();i+=3) {
373                 const int smooth = 1;
374     int t1 = mIndices[i];
375     int t2 = mIndices[i+1];
376                 int t3 = mIndices[i+2];
377                 //errMsg("NM"," tri"<<t1<<" "<<t2<<" "<<t3 );
378
379                 ntlTriangle tri;
380
381                 tri.getPoints()[0] = t1+iniVertIndex;
382                 tri.getPoints()[1] = t2+iniVertIndex;
383                 tri.getPoints()[2] = t3+iniVertIndex;
384
385                 /* init flags */
386                 int flag = 0; 
387                 if(getVisible()){ flag |= TRI_GEOMETRY; }
388                 if(getCastShadows() ) { 
389                         flag |= TRI_CASTSHADOWS; } 
390                 if( (getMaterial()->getMirror()>0.0) ||  
391                                 (getMaterial()->getTransparence()>0.0) ||  
392                                 (getMaterial()->getFresnel()>0.0) ) { 
393                         flag |= TRI_MAKECAUSTICS; } 
394                 else { 
395                         flag |= TRI_NOCAUSTICS; } 
396
397                 /* init geo init id */
398                 int geoiId = getGeoInitId(); 
399                 if(geoiId > 0) { 
400                         flag |= (1<< (geoiId+4)); 
401                         flag |= mGeoInitType; 
402                 } 
403
404                 tri.setFlags( flag );
405
406                 /* triangle normal missing */
407                 tri.setNormal( ntlVec3Gfx(0.0) );
408                 tri.setSmoothNormals( smooth );
409                 tri.setObjectId( objectId );
410                 triangles->push_back( tri ); 
411         }
412         //errMsg("N3"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
413         return;
414 }
415
416
417
418 inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) {
419         // WARNING - this requires a security boundary layer... 
420         ntlVec3Gfx ret(0.0);
421         ret[0] = *getData(i-1,j  ,k  ) - 
422                  *getData(i+1,j  ,k  );
423         ret[1] = *getData(i  ,j-1,k  ) - 
424                  *getData(i  ,j+1,k  );
425         ret[2] = *getData(i  ,j  ,k-1  ) - 
426                  *getData(i  ,j  ,k+1  );
427         return ret;
428 }
429
430
431
432
433
434 /******************************************************************************
435  * 
436  * Surface improvement
437  * makes use of trimesh2 library
438  * http://www.cs.princeton.edu/gfx/proj/trimesh2/
439  *
440  * Copyright (c) 2004 Szymon Rusinkiewicz.
441  * see COPYING_trimesh2
442  * 
443  *****************************************************************************/
444
445 // Diffuse a vector field at 1 vertex, weighted by
446 // a gaussian of width 1/sqrt(invsigma2)
447 void IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int v, float invsigma2, ntlVec3Gfx &flt)
448 {
449         flt = ntlVec3Gfx(0.0);
450         flt += *(field+pointerScale*v) *pointareas[v];
451         float sum_w = pointareas[v];
452         const ntlVec3Gfx &nv = mPoints[v].n;
453
454         unsigned &flag = flag_curr;
455         flag++;
456         flags[v] = flag;
457         vector<int> boundary = neighbors[v];
458         while (!boundary.empty()) {
459                 const int bbn = boundary.back();
460                 boundary.pop_back();
461                 if (flags[bbn] == flag) continue;
462                 flags[bbn] = flag;
463
464                 // gaussian weight of width 1/sqrt(invsigma2)
465                 const float d2 = invsigma2 * normNoSqrt(mPoints[bbn].v - mPoints[v].v);
466                 if(d2 >= 9.0f) continue; // 25 also possible  , slower
467                 //float w = (d2 >=  9.0f) ? 0.0f : exp(-0.5f*d2);
468                 float w = exp(-0.5f*d2);
469                 if(dot(nv, mPoints[bbn].n) <= 0.0f) continue; // faster than before d2 calc?
470
471                 // Downweight things pointing in different directions
472                 w *= dot(nv , mPoints[bbn].n);
473                 // Surface area "belonging" to each point
474                 w *= pointareas[bbn];
475                 // Accumulate weight times field at neighbor
476                 flt += *(field+pointerScale*bbn)*w;
477
478                 sum_w += w;
479                 for (int i = 0; i < (int)neighbors[bbn].size(); i++) {
480                         int nn = neighbors[bbn][i];
481                         if (flags[nn] == flag) continue;
482                         boundary.push_back(nn);
483                 }
484         }
485         flt /= sum_w;
486 }
487
488 void IsoSurface::smoothSurface(float sigma)
489 {
490         int nv = mPoints.size();
491         if ((int)flags.size() != nv) flags.resize(nv);
492         int nf = mIndices.size()/3;
493
494         { // need neighbors
495
496                 vector<int> numneighbors(mPoints.size());
497                 int i;
498                 for (i = 0; i < (int)mIndices.size()/3; i++) {
499                         numneighbors[mIndices[i*3+0]]++;
500                         numneighbors[mIndices[i*3+1]]++;
501                         numneighbors[mIndices[i*3+2]]++;
502                 }
503
504                 neighbors.clear();
505                 neighbors.resize(mPoints.size());
506                 for (i = 0; i < (int)mPoints.size(); i++) {
507                         neighbors[i].clear();
508                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
509                 }
510
511                 for (i = 0; i < (int)mIndices.size()/3; i++) {
512                         for (int j = 0; j < 3; j++) {
513                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
514                                 int n1 =  mIndices[i*3+((j+1)%3)];
515                                 int n2 =  mIndices[i*3+((j+2)%3)];
516                                 if (std::find(me.begin(), me.end(), n1) == me.end())
517                                         me.push_back(n1);
518                                 if (std::find(me.begin(), me.end(), n2) == me.end())
519                                         me.push_back(n2);
520                         }
521                 }
522         } // need neighbor
523
524         { // need pointarea
525                 pointareas.clear();
526                 pointareas.resize(nv);
527                 cornerareas.clear();
528                 cornerareas.resize(nf);
529
530                 for (int i = 0; i < nf; i++) {
531                         // Edges
532                         ntlVec3Gfx e[3] = { 
533                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
534                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
535                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
536
537                         // Compute corner weights
538                         float area = 0.5f * norm( cross(e[0], e[1]));
539                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
540                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
541                                         l2[1] * (l2[2] + l2[0] - l2[1]),
542                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
543                         if (ew[0] <= 0.0f) {
544                                 cornerareas[i][1] = -0.25f * l2[2] * area /
545                                                                 dot(e[0] , e[2]);
546                                 cornerareas[i][2] = -0.25f * l2[1] * area /
547                                                                 dot(e[0] , e[1]);
548                                 cornerareas[i][0] = area - cornerareas[i][1] -
549                                                                 cornerareas[i][2];
550                         } else if (ew[1] <= 0.0f) {
551                                 cornerareas[i][2] = -0.25f * l2[0] * area /
552                                                                 dot(e[1] , e[0]);
553                                 cornerareas[i][0] = -0.25f * l2[2] * area /
554                                                                 dot(e[1] , e[2]);
555                                 cornerareas[i][1] = area - cornerareas[i][2] -
556                                                                 cornerareas[i][0];
557                         } else if (ew[2] <= 0.0f) {
558                                 cornerareas[i][0] = -0.25f * l2[1] * area /
559                                                                 dot(e[2] , e[1]);
560                                 cornerareas[i][1] = -0.25f * l2[0] * area /
561                                                                 dot(e[2] , e[0]);
562                                 cornerareas[i][2] = area - cornerareas[i][0] -
563                                                                 cornerareas[i][1];
564                         } else {
565                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
566                                 for (int j = 0; j < 3; j++)
567                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
568                                                                                          ew[(j+2)%3]);
569                         }
570
571                         // NT important, check this...
572 #ifndef WIN32
573                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
574                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
575                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
576 #else // WIN32
577                         // FIXME check as well...
578                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
579                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
580                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
581 #endif // WIN32
582
583                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
584                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
585                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
586                 }
587
588         } // need pointarea
589         // */
590
591         float invsigma2 = 1.0f / (sigma*sigma);
592
593         vector<ntlVec3Gfx> dflt(nv);
594         for (int i = 0; i < nv; i++) {
595                 diffuseVertexField( &mPoints[0].v, 2,
596                                    i, invsigma2, dflt[i]);
597                 // Just keep the displacement
598                 dflt[i] -= mPoints[i].v;
599         }
600
601         // Slightly better small-neighborhood approximation
602         for (int i = 0; i < nf; i++) {
603                 ntlVec3Gfx c = mPoints[mIndices[i*3+0]].v +
604                           mPoints[mIndices[i*3+1]].v +
605                           mPoints[mIndices[i*3+2]].v;
606                 c /= 3.0f;
607                 for (int j = 0; j < 3; j++) {
608                         int v = mIndices[i*3+j];
609                         ntlVec3Gfx d =(c - mPoints[v].v) * 0.5f;
610                         dflt[v] += d * (cornerareas[i][j] /
611                                    pointareas[mIndices[i*3+j]] *
612                                    exp(-0.5f * invsigma2 * normNoSqrt(d)) );
613                 }
614         }
615
616         // Filter displacement field
617         vector<ntlVec3Gfx> dflt2(nv);
618         for (int i = 0; i < nv; i++) {
619                 diffuseVertexField( &dflt[0], 1,
620                                    i, invsigma2, dflt2[i]);
621         }
622
623         // Update vertex positions
624         for (int i = 0; i < nv; i++) {
625                 mPoints[i].v += dflt[i] - dflt2[i]; // second Laplacian
626         }
627
628         // when normals smoothing off, this cleans up quite well
629         // costs ca. 50% additional time though
630         float nsFac = 1.5f;
631         { float ninvsigma2 = 1.0f / (nsFac*nsFac*sigma*sigma);
632                 for (int i = 0; i < nv; i++) {
633                         diffuseVertexField( &mPoints[0].n, 2, i, ninvsigma2, dflt[i]);
634                         normalize(dflt[i]);
635                 }
636                 for (int i = 0; i < nv; i++) {
637                         mPoints[i].n = dflt[i];
638                 } 
639         } // smoothNormals copy */
640
641         //errMsg("SMSURF","done v:"<<sigma); // DEBUG
642 }
643
644 void IsoSurface::smoothNormals(float sigma)
645 {
646         { // need neighbor
647                 vector<int> numneighbors(mPoints.size());
648                 int i;
649                 for (i = 0; i < (int)mIndices.size()/3; i++) {
650                         numneighbors[mIndices[i*3+0]]++;
651                         numneighbors[mIndices[i*3+1]]++;
652                         numneighbors[mIndices[i*3+2]]++;
653                 }
654
655                 neighbors.clear();
656                 neighbors.resize(mPoints.size());
657                 for (i = 0; i < (int)mPoints.size(); i++) {
658                         neighbors[i].clear();
659                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
660                 }
661
662                 for (i = 0; i < (int)mIndices.size()/3; i++) {
663                         for (int j = 0; j < 3; j++) {
664                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
665                                 int n1 =  mIndices[i*3+((j+1)%3)];
666                                 int n2 =  mIndices[i*3+((j+2)%3)];
667                                 if (std::find(me.begin(), me.end(), n1) == me.end())
668                                         me.push_back(n1);
669                                 if (std::find(me.begin(), me.end(), n2) == me.end())
670                                         me.push_back(n2);
671                         }
672                 }
673         } // need neighbor
674
675         { // need pointarea
676                 int nf = mIndices.size()/3, nv = mPoints.size();
677                 pointareas.clear();
678                 pointareas.resize(nv);
679                 cornerareas.clear();
680                 cornerareas.resize(nf);
681
682                 for (int i = 0; i < nf; i++) {
683                         // Edges
684                         ntlVec3Gfx e[3] = { 
685                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
686                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
687                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
688
689                         // Compute corner weights
690                         float area = 0.5f * norm( cross(e[0], e[1]));
691                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
692                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
693                                         l2[1] * (l2[2] + l2[0] - l2[1]),
694                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
695                         if (ew[0] <= 0.0f) {
696                                 cornerareas[i][1] = -0.25f * l2[2] * area /
697                                                                 dot(e[0] , e[2]);
698                                 cornerareas[i][2] = -0.25f * l2[1] * area /
699                                                                 dot(e[0] , e[1]);
700                                 cornerareas[i][0] = area - cornerareas[i][1] -
701                                                                 cornerareas[i][2];
702                         } else if (ew[1] <= 0.0f) {
703                                 cornerareas[i][2] = -0.25f * l2[0] * area /
704                                                                 dot(e[1] , e[0]);
705                                 cornerareas[i][0] = -0.25f * l2[2] * area /
706                                                                 dot(e[1] , e[2]);
707                                 cornerareas[i][1] = area - cornerareas[i][2] -
708                                                                 cornerareas[i][0];
709                         } else if (ew[2] <= 0.0f) {
710                                 cornerareas[i][0] = -0.25f * l2[1] * area /
711                                                                 dot(e[2] , e[1]);
712                                 cornerareas[i][1] = -0.25f * l2[0] * area /
713                                                                 dot(e[2] , e[0]);
714                                 cornerareas[i][2] = area - cornerareas[i][0] -
715                                                                 cornerareas[i][1];
716                         } else {
717                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
718                                 for (int j = 0; j < 3; j++)
719                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
720                                                                                          ew[(j+2)%3]);
721                         }
722
723                         // NT important, check this...
724 #ifndef WIN32
725                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
726                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
727                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
728 #else // WIN32
729                         // FIXME check as well...
730                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
731                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
732                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
733 #endif // WIN32
734
735                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
736                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
737                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
738                 }
739
740         } // need pointarea
741
742         int nv = mPoints.size();
743         if ((int)flags.size() != nv) flags.resize(nv);
744         float invsigma2 = 1.0f / (sigma*sigma);
745
746         vector<ntlVec3Gfx> nflt(nv);
747         for (int i = 0; i < nv; i++) {
748                 diffuseVertexField( &mPoints[0].n, 2, i, invsigma2, nflt[i]);
749                 normalize(nflt[i]);
750         }
751
752         for (int i = 0; i < nv; i++) {
753                 mPoints[i].n = nflt[i];
754         }
755
756         //errMsg("SMNRMLS","done v:"<<sigma); // DEBUG
757 }
758
759