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