initial commit of the fluid simulator.
[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 // IMPL ----------------------------------------------------------------------------------------
11 #include "isosurface.h"
12 #include "mcubes_tables.h"
13 #include "ntl_scene.h"
14
15 #include <algorithm>
16 #include <stdio.h>
17 #define MCUBES_MAXPOLNUM  10000
18 #define MCUBES_MAXVERTNUM 30000
19
20
21
22 /******************************************************************************
23  * Constructor
24  *****************************************************************************/
25 IsoSurface::IsoSurface(double iso, double blend) :
26         ntlGeometryObject(),
27         mSizex(-1), mSizey(-1), mSizez(-1),
28   mIsoValue( iso ), 
29         mBlendVal( blend ),
30         mPoints(), 
31         mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL),
32         mIndices(),
33
34   mStart(0.0), mEnd(0.0), mDomainExtent(0.0),
35   mInitDone(false),
36         mLoopSubdivs(0), mSmoothSurface(0.0), mSmoothNormals(0.0),
37         mAcrossEdge(), mAdjacentFaces()
38 {
39 }
40
41
42 /******************************************************************************
43  * The real init...
44  *****************************************************************************/
45 void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent) 
46 {
47         // init solver and size
48         mSizex = setx;
49         mSizey = sety;
50         if(setz == 1) {// 2D, create thin 2D surface
51                 setz = 5; 
52         }
53         mSizez = setz;
54         mDomainExtent = extent;
55
56         /* check triangulation size (for raytraing) */
57         if( ( mStart[0] >= mEnd[0] ) && ( mStart[1] >= mEnd[1] ) && ( mStart[2] >= mEnd[2] ) ){
58                 // extent was not set, use normalized one from parametrizer
59                 mStart = ntlVec3Gfx(0.0) - extent*0.5;
60                 mEnd = ntlVec3Gfx(0.0) + extent*0.5;
61         }
62
63   // init 
64         mIndices.clear();
65   mPoints.clear();
66
67         int nodes = mSizez*mSizey*mSizex;
68   mpData = new float[nodes];
69   for(int i=0;i<nodes;i++) { mpData[i] = 0.0; }
70
71   // allocate edge arrays  (last slices are never used...)
72   mpEdgeVerticesX = new int[nodes];
73   mpEdgeVerticesY = new int[nodes];
74   mpEdgeVerticesZ = new int[nodes];
75   for(int i=0;i<nodes;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
76   
77         // marching cubes are ready 
78         mInitDone = true;
79         unsigned long int memCnt = (3*sizeof(int)*nodes + sizeof(float)*nodes);
80         double memd = memCnt;
81         char *sizeStr = "";
82         const double sfac = 1000.0;
83         if(memd>sfac){ memd /= sfac; sizeStr="KB"; }
84         if(memd>sfac){ memd /= sfac; sizeStr="MB"; }
85         if(memd>sfac){ memd /= sfac; sizeStr="GB"; }
86         if(memd>sfac){ memd /= sfac; sizeStr="TB"; }
87
88         debMsgStd("IsoSurface::initializeIsosurface",DM_MSG,"Inited "<<PRINT_VEC(setx,sety,setz)<<" alloced:"<< memd<<" "<<sizeStr<<"." ,10);
89 }
90
91
92
93
94 /******************************************************************************
95  * Destructor
96  *****************************************************************************/
97 IsoSurface::~IsoSurface( void )
98 {
99
100         if(mpEdgeVerticesX) delete [] mpEdgeVerticesX;
101         if(mpEdgeVerticesY) delete [] mpEdgeVerticesY;
102         if(mpEdgeVerticesZ) delete [] mpEdgeVerticesZ;
103
104 }
105
106
107
108
109 /******************************************************************************
110  * triangulate the scalar field given by pointer
111  *****************************************************************************/
112 void IsoSurface::triangulate( void )
113 {
114   double gsx,gsy,gsz; // grid spacing in x,y,z direction
115   double px,py,pz;    // current position in grid in x,y,z direction
116   IsoLevelCube cubie;    // struct for a small subcube
117         myTime_t tritimestart = getTime(); 
118
119         if(!mpData) {
120                 errorOut("IsoSurface::triangulate fatal error: no LBM object, and no scalar field...!");
121                 exit( -1 );
122         }
123
124   // get grid spacing
125   gsx = (mEnd[0]-mStart[0])/(double)(mSizex-1.0);
126   gsy = (mEnd[1]-mStart[1])/(double)(mSizey-1.0);
127   gsz = (mEnd[2]-mStart[2])/(double)(mSizez-1.0);
128
129   // clean up previous frame
130         mIndices.clear();
131         mPoints.clear();
132
133         // reset edge vertices
134   for(int i=0;i<(mSizex*mSizey*mSizez);i++) {
135                 mpEdgeVerticesX[i] = -1;
136                 mpEdgeVerticesY[i] = -1;
137                 mpEdgeVerticesZ[i] = -1;
138         }
139
140         ntlVec3Gfx pos[8];
141         float value[8];
142         int cubeIndex;      // index entry of the cube 
143         int triIndices[12]; // vertex indices 
144         int *eVert[12];
145         IsoLevelVertex ilv;
146
147         // edges between which points?
148         const int mcEdges[24] = { 
149                 0,1,  1,2,  2,3,  3,0,
150                 4,5,  5,6,  6,7,  7,4,
151                 0,4,  1,5,  2,6,  3,7 };
152
153         const int cubieOffsetX[8] = {
154                 0,1,1,0,  0,1,1,0 };
155         const int cubieOffsetY[8] = {
156                 0,0,1,1,  0,0,1,1 };
157         const int cubieOffsetZ[8] = {
158                 0,0,0,0,  1,1,1,1 };
159
160   // let the cubes march 
161         pz = mStart[2]-gsz;
162         for(int k=0;k<(mSizez-2);k++) {
163                 pz += gsz;
164                 py = mStart[1]-gsy;
165                 for(int j=0;j<(mSizey-2);j++) {
166       py += gsy;
167                         px = mStart[0]-gsx;
168                         for(int i=0;i<(mSizex-2);i++) {
169                         px += gsx;
170                                 int baseIn = ISOLEVEL_INDEX( i+0, j+0, k+0);
171
172                                 value[0] = *getData(i  ,j  ,k  );
173                                 value[1] = *getData(i+1,j  ,k  );
174                                 value[2] = *getData(i+1,j+1,k  );
175                                 value[3] = *getData(i  ,j+1,k  );
176                                 value[4] = *getData(i  ,j  ,k+1);
177                                 value[5] = *getData(i+1,j  ,k+1);
178                                 value[6] = *getData(i+1,j+1,k+1);
179                                 value[7] = *getData(i  ,j+1,k+1);
180                                 //errMsg("ISOT2D"," at "<<PRINT_IJK<<" "
181                                                 //<<" v0="<<value[0] <<" v1="<<value[1] <<" v2="<<value[2] <<" v3="<<value[3]
182                                                 //<<" v4="<<value[4] <<" v5="<<value[5] <<" v6="<<value[6] <<" v7="<<value[7] );
183
184                                 // check intersections of isosurface with edges, and calculate cubie index
185                                 cubeIndex = 0;
186                                 if (value[0] < mIsoValue) cubeIndex |= 1;
187                                 if (value[1] < mIsoValue) cubeIndex |= 2;
188                                 if (value[2] < mIsoValue) cubeIndex |= 4;
189                                 if (value[3] < mIsoValue) cubeIndex |= 8;
190                                 if (value[4] < mIsoValue) cubeIndex |= 16;
191                                 if (value[5] < mIsoValue) cubeIndex |= 32;
192                                 if (value[6] < mIsoValue) cubeIndex |= 64;
193                                 if (value[7] < mIsoValue) cubeIndex |= 128;
194
195                                 // No triangles to generate?
196                                 if (mcEdgeTable[cubeIndex] == 0) {
197                                         continue;
198                                 }
199
200                                 // where to look up if this point already exists
201                                 eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
202                                 eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
203                                 eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, k+0) ];
204                                 eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
205
206                                 eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, k+1) ];
207                                 eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, k+1) ];
208                                 eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, k+1) ];
209                                 eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, k+1) ];
210
211                                 eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
212                                 eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, k+0) ];
213                                 eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, k+0) ];
214                                 eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, k+0) ];
215
216                                 // grid positions
217                                 pos[0] = ntlVec3Gfx(px    ,py    ,pz);
218                                 pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
219                                 pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
220                                 pos[3] = ntlVec3Gfx(px    ,py+gsy,pz);
221                                 pos[4] = ntlVec3Gfx(px    ,py    ,pz+gsz);
222                                 pos[5] = ntlVec3Gfx(px+gsx,py    ,pz+gsz);
223                                 pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz);
224                                 pos[7] = ntlVec3Gfx(px    ,py+gsy,pz+gsz);
225
226                                 // check all edges
227                                 for(int e=0;e<12;e++) {
228
229                                         if (mcEdgeTable[cubeIndex] & (1<<e)) {
230                                                 // is the vertex already calculated?
231                                                 if(*eVert[ e ] < 0) {
232                                                         // interpolate edge
233                                                         const int e1 = mcEdges[e*2  ];
234                                                         const int e2 = mcEdges[e*2+1];
235                                                         const ntlVec3Gfx p1 = pos[ e1  ];    // scalar field pos 1
236                                                         const ntlVec3Gfx p2 = pos[ e2  ];    // scalar field pos 2
237                                                         const float valp1  = value[ e1  ];  // scalar field val 1
238                                                         const float valp2  = value[ e2  ];  // scalar field val 2
239                                                         //double mu;            // interpolation value
240                                                         //ntlVec3Gfx p;         // new point
241
242                                                         // choose if point should be calculated by interpolation,
243                                                         // or "Karolin" method 
244
245                                                         //double deltaVal = ABS(valp2-valp1);
246                                                         //if(deltaVal >-10.0) {
247                                                         // standard interpolation
248                                                         //vertInfo[e].type = 0; 
249                                                         /*if (ABS(mIsoValue-valp1) < 0.000000001) {
250                                                                 mu = 0.0;
251                                                         } else {
252                                                                 if (ABS(mIsoValue-valp2) < 0.000000001) {
253                                                                         mu = 1.0;
254                                                                 } else {
255                                                                         mu = (mIsoValue - valp1) / (valp2 - valp1);
256                                                                 }
257                                                         } */
258                                                         /*} else {
259                                                                 errorOut(" ? ");
260                                                         // use fill grade (=karo) 
261                                                         vertInfo[e].type = 1; 
262                                                         if(valp1 < valp2) { mu = 1.0- (valp1 + valp2 - mIsoValue);
263                                                         } else { mu = 0.0+ (valp1 + valp2 - mIsoValue); } 
264                                                         } */
265
266                                                         //const float mu = (mIsoValue - valp1) / (valp2 - valp1);
267                                                         float mu;
268                                                         if(valp1 < valp2) {
269                                                                 mu = 1.0 - 1.0*(valp1 + valp2 - mIsoValue);
270                                                         } else {
271                                                                 mu = 0.0 + 1.0*(valp1 + valp2 - mIsoValue);
272                                                         }
273
274                                                         float isov2 = mIsoValue;
275                                                         isov2 = (valp1+valp2)*0.5;
276                                                         mu = (isov2 - valp1) / (valp2 - valp1);
277                                                         mu = (isov2) / (valp2 - valp1);
278
279                                                         mu = (mIsoValue - valp1) / (valp2 - valp1);
280                                                         //mu *= mu;
281
282
283
284                                                         // init isolevel vertex
285                                                         ilv.v = p1 + (p2-p1)*mu;
286                                                         ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) +
287                                                                                         getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * (    mu) ;
288                                                         mPoints.push_back( ilv );
289
290                                                         triIndices[e] = (mPoints.size()-1);
291                                                         // store vertex 
292                                                         *eVert[ e ] = triIndices[e];
293                                                 }       else {
294                                                         // retrieve  from vert array
295                                                         triIndices[e] = *eVert[ e ];
296                                                 }
297                                         } // along all edges 
298
299                                 }
300
301
302                                 // Create the triangles... 
303                                 for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
304                                         //errMsg("MC","tri "<<mIndices.size() <<" "<< triIndices[ mcTriTable[cubeIndex][e+0] ]<<" "<< triIndices[ mcTriTable[cubeIndex][e+1] ]<<" "<< triIndices[ mcTriTable[cubeIndex][e+2] ] );
305                                         mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
306                                         mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] );
307                                         mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
308                                 }
309
310                                 
311       }
312     }
313
314   } // k
315   
316
317   // precalculate normals using an approximation of the scalar field gradient 
318         for(int ni=0;ni<(int)mPoints.size();ni++) {
319                 // use triangle normals?, this seems better for File-IsoSurf
320                 normalize( mPoints[ni].n );
321         }
322
323         if(mSmoothSurface>0.0) {
324                 // not needed for post normal smoothing? 
325                 // if(mSmoothNormals<=0.0) { smoothNormals(mSmoothSurface*0.5); }
326                 smoothSurface(mSmoothSurface);
327         }
328         for(int i=0; i<mLoopSubdivs; i++) {
329                 subdivide();
330         }
331         if(mSmoothNormals>0.0) {
332                 smoothNormals(mSmoothNormals);
333         }
334         myTime_t tritimeend = getTime(); 
335 #if ELBEEM_BLENDER!=1
336         debMsgStd("IsoSurface::triangulate",DM_MSG,"Took "<< ((tritimeend-tritimestart)/(double)1000.0)<<"s) " , 10 );
337 #endif // ELBEEM_BLENDER!=1
338 }
339
340
341         
342
343
344 /******************************************************************************
345  * Get triangles for rendering
346  *****************************************************************************/
347 void IsoSurface::getTriangles( vector<ntlTriangle> *triangles, 
348                                                                                                          vector<ntlVec3Gfx> *vertices, 
349                                                                                                          vector<ntlVec3Gfx> *normals, int objectId )
350 {
351         if(!mInitDone) {
352                 debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10);
353                 return;
354         }
355
356   /* triangulate field */
357   triangulate();
358         //errMsg("TRIS"," "<<mIndices.size() );
359
360         // new output with vertice reuse
361         int iniVertIndex = (*vertices).size();
362         int iniNormIndex = (*normals).size();
363         if(iniVertIndex != iniNormIndex) {
364                 errorOut("getTriangles Error for '"<<mName<<"': Vertices and normal array sizes to not match!!!");
365                 exit(1); }
366         //errMsg("NM"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
367         //errMsg("NM"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
368
369   for(int i=0;i<(int)mPoints.size();i++) {
370                 vertices->push_back( mPoints[i].v );
371         }
372   for(int i=0;i<(int)mPoints.size();i++) {
373                 normals->push_back( mPoints[i].n );
374         }
375
376         //errMsg("N2"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
377         //errMsg("N2"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
378
379   for(int i=0;i<(int)mIndices.size();i+=3) {
380                 const int smooth = 1;
381     int t1 = mIndices[i];
382     int t2 = mIndices[i+1];
383                 int t3 = mIndices[i+2];
384                 //errMsg("NM"," tri"<<t1<<" "<<t2<<" "<<t3 );
385
386                 ntlTriangle tri;
387
388                 tri.getPoints()[0] = t1+iniVertIndex;
389                 tri.getPoints()[1] = t2+iniVertIndex;
390                 tri.getPoints()[2] = t3+iniVertIndex;
391
392                 /* init flags */
393                 int flag = 0; 
394                 if(getVisible()){ flag |= TRI_GEOMETRY; }
395                 if(getCastShadows() ) { 
396                         flag |= TRI_CASTSHADOWS; } 
397                 if( (getMaterial()->getMirror()>0.0) ||  
398                                 (getMaterial()->getTransparence()>0.0) ||  
399                                 (getMaterial()->getFresnel()>0.0) ) { 
400                         flag |= TRI_MAKECAUSTICS; } 
401                 else { 
402                         flag |= TRI_NOCAUSTICS; } 
403
404                 /* init geo init id */
405                 int geoiId = getGeoInitId(); 
406                 if(geoiId > 0) { 
407                         flag |= (1<< (geoiId+4)); 
408                         flag |= mGeoInitType; 
409                 } 
410
411                 tri.setFlags( flag );
412
413                 /* triangle normal missing */
414                 tri.setNormal( ntlVec3Gfx(0.0) );
415                 tri.setSmoothNormals( smooth );
416                 tri.setObjectId( objectId );
417                 triangles->push_back( tri ); 
418         }
419         //errMsg("N3"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
420         return;
421 }
422
423
424
425 inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) {
426         // WARNING - this assumes a security boundary layer... 
427         /*
428         if(i<=0) i=1;
429         if(j<=0) j=1;
430         if(k<=0) k=1;
431         if(i>=(mSizex-1)) i=mSizex-2;
432         if(j>=(mSizex-1)) j=mSizex-2;
433         if(k>=(mSizex-1)) k=mSizex-2; // */
434
435         ntlVec3Gfx ret(0.0);
436         ret[0] = *getData(i-1,j  ,k  ) - 
437                  *getData(i+1,j  ,k  );
438         ret[1] = *getData(i  ,j-1,k  ) - 
439                  *getData(i  ,j+1,k  );
440         ret[2] = *getData(i  ,j  ,k-1  ) - 
441                  *getData(i  ,j  ,k+1  );
442         return ret;
443 }
444
445 #define RECALCNORMALS 0
446
447 // Subdivide a mesh allways loop
448 void IsoSurface::subdivide()
449 {
450         int i;
451
452         mAdjacentFaces.clear(); 
453         mAcrossEdge.clear();
454
455         //void TriMesh::need_adjacentfaces()
456         {
457                 vector<int> numadjacentfaces(mPoints.size());
458                 //errMsg("SUBDIV ADJFA1", " "<<mPoints.size()<<" - "<<numadjacentfaces.size() );
459                 int i;
460                 for (i = 0; i < (int)mIndices.size()/3; i++) {
461                         numadjacentfaces[mIndices[i*3 + 0]]++;
462                         numadjacentfaces[mIndices[i*3 + 1]]++;
463                         numadjacentfaces[mIndices[i*3 + 2]]++;
464                 }
465
466                 mAdjacentFaces.resize(mPoints.size());
467                 for (i = 0; i < (int)mPoints.size(); i++)
468                         mAdjacentFaces[i].reserve(numadjacentfaces[i]);
469
470                 for (i = 0; i < (int)mIndices.size()/3; i++) {
471                         for (int j = 0; j < 3; j++)
472                                 mAdjacentFaces[mIndices[i*3 + j]].push_back(i);
473                 }
474
475         }
476
477         // Find the face across each edge from each other face (-1 on boundary)
478         // If topology is bad, not necessarily what one would expect...
479         //void TriMesh::need_across_edge()
480         {
481                 mAcrossEdge.resize(mIndices.size(), -1);
482
483                 for (int i = 0; i < (int)mIndices.size()/3; i++) {
484                         for (int j = 0; j < 3; j++) {
485                                 if (mAcrossEdge[i*3 + j] != -1)
486                                         continue;
487                                 int v1 = mIndices[i*3 + ((j+1)%3)];
488                                 int v2 = mIndices[i*3 + ((j+2)%3)];
489                                 const vector<int> &a1 = mAdjacentFaces[v1];
490                                 const vector<int> &a2 = mAdjacentFaces[v2];
491                                 for (int k1 = 0; k1 < (int)a1.size(); k1++) {
492                                         int other = a1[k1];
493                                         if (other == i)
494                                                 continue;
495                                         vector<int>::const_iterator it =
496                                                 std::find(a2.begin(), a2.end(), other);
497                                         if (it == a2.end())
498                                                 continue;
499
500                                         //int ind = (faces[other].indexof(v1)+1)%3;
501                                         int ind = -1;
502                                         if( mIndices[other*3+0] == (unsigned int)v1 ) ind = 0;
503                                         else if( mIndices[other*3+1] == (unsigned int)v1 ) ind = 1;
504                                         else if( mIndices[other*3+2] == (unsigned int)v1 ) ind = 2;
505                                         ind = (ind+1)%3;
506
507                                         if ( (int)mIndices[other*3 + ((ind+1)%3)] != v2)
508                                                 continue;
509                                         mAcrossEdge[i*3 + j] = other;
510                                         mAcrossEdge[other*3 + ind] = i;
511                                         break;
512                                 }
513                         }
514                 }
515
516                 //errMsg("SUBDIV ACREDG", "Done.\n");
517         }
518
519         //errMsg("SUBDIV","start");
520         // Introduce new vertices
521         int nf = (int)mIndices.size() / 3;
522
523         //vector<TriMesh::Face> newverts(nf, TriMesh::Face(-1,-1,-1));
524         vector<int> newverts(nf*3); //, TriMesh::Face(-1,-1,-1));
525         for(int j=0; j<(int)newverts.size(); j++) newverts[j] = -1;
526
527         int old_nv = (int)mPoints.size();
528         mPoints.reserve(4 * old_nv);
529         vector<int> newvert_count(old_nv + 3*nf); // wichtig...?
530         //errMsg("NC", newvert_count.size() );
531
532         for (i = 0; i < nf; i++) {
533                 for (int j = 0; j < 3; j++) {
534                         int ae = mAcrossEdge[i*3 + j];
535                         if (newverts[i*3 + j] == -1 && ae != -1) {
536                                 if (mAcrossEdge[ae*3 + 0] == i)
537                                         newverts[i*3 + j] = newverts[ae*3 + 0];
538                                 else if (mAcrossEdge[ae*3 + 1] == i)
539                                         newverts[i*3 + j] = newverts[ae*3 + 1];
540                                 else if (mAcrossEdge[ae*3 + 2] == i)
541                                         newverts[i*3 + j] = newverts[ae*3 + 2];
542                         }
543                         if (newverts[i*3 + j] == -1) {
544                                 IsoLevelVertex ilv;
545                                 ilv.v = ntlVec3Gfx(0.0);
546                                 ilv.n = ntlVec3Gfx(0.0);
547                                 mPoints.push_back(ilv);
548                                 newverts[i*3 + j] = (int)mPoints.size() - 1;
549                                 if (ae != -1) {
550                                         if (mAcrossEdge[ae*3 + 0] == i)
551                                                 newverts[ae*3 + 0] = newverts[i*3 + j];
552                                         else if (mAcrossEdge[ae*3 + 1] == i)
553                                                 newverts[ae*3 + 1] = newverts[i*3 + j];
554                                         else if (mAcrossEdge[ae*3 + 2] == i)
555                                                 newverts[ae*3 + 2] = newverts[i*3 + j];
556                                 }
557                         }
558                         if(ae != -1) {
559                                 mPoints[newverts[i*3 + j]].v +=
560                                         mPoints[ mIndices[i*3 + ( j     )] ].v * 0.25f  + // j = 0,1,2?
561                                         mPoints[ mIndices[i*3 + ((j+1)%3)] ].v * 0.375f +
562                                         mPoints[ mIndices[i*3 + ((j+2)%3)] ].v * 0.375f;
563 #if RECALCNORMALS==0
564                                 mPoints[newverts[i*3 + j]].n +=
565                                         mPoints[ mIndices[i*3 + ( j     )] ].n * 0.25f  + // j = 0,1,2?
566                                         mPoints[ mIndices[i*3 + ((j+1)%3)] ].n * 0.375f +
567                                         mPoints[ mIndices[i*3 + ((j+2)%3)] ].n * 0.375f;
568 #endif // RECALCNORMALS==0
569                         } else {
570                                 mPoints[newverts[i*3 + j]].v +=
571                                         mPoints[ mIndices[i*3 + ((j+1)%3)] ].v * 0.5f   +
572                                         mPoints[ mIndices[i*3 + ((j+2)%3)] ].v * 0.5f  ;
573 #if RECALCNORMALS==0
574                                 mPoints[newverts[i*3 + j]].n +=
575                                         mPoints[ mIndices[i*3 + ((j+1)%3)] ].n * 0.5f   +
576                                         mPoints[ mIndices[i*3 + ((j+2)%3)] ].n * 0.5f  ;
577 #endif // RECALCNORMALS==0
578                         }
579
580                         newvert_count[newverts[i*3 + j]]++;
581                 }
582         }
583         for (i = old_nv; i < (int)mPoints.size(); i++) {
584                 if (!newvert_count[i])
585                         continue;
586                 float scale = 1.0f / newvert_count[i];
587                 mPoints[i].v *= scale;
588
589 #if RECALCNORMALS==0
590                 //mPoints[i].n *= scale;
591                 //normalize( mPoints[i].n );
592 #endif // RECALCNORMALS==0
593         }
594
595         // Update old vertices
596         for (i = 0; i < old_nv; i++) {
597                         ntlVec3Gfx bdyavg(0.0), nbdyavg(0.0);
598                         ntlVec3Gfx norm_bdyavg(0.0), norm_nbdyavg(0.0); // N
599                         int nbdy = 0, nnbdy = 0;
600                         int naf = (int)mAdjacentFaces[i].size();
601                         if (!naf)
602                                 continue;
603                         for (int j = 0; j < naf; j++) {
604                                 int af = mAdjacentFaces[i][j];
605
606                                 int afi = -1;
607                                 if( mIndices[af*3+0] == (unsigned int)i ) afi = 0;
608                                 else if( mIndices[af*3+1] == (unsigned int)i ) afi = 1;
609                                 else if( mIndices[af*3+2] == (unsigned int)i ) afi = 2;
610
611                                 int n1 = (afi+1) % 3;
612                                 int n2 = (afi+2) % 3;
613                                 if (mAcrossEdge[af*3 + n1] == -1) {
614                                         bdyavg += mPoints[newverts[af*3 + n1]].v;
615 #if RECALCNORMALS==0
616                                         //norm_bdyavg += mPoints[newverts[af*3 + n1]].n;
617 #endif // RECALCNORMALS==0
618                                         nbdy++;
619                                 } else {
620                                         nbdyavg += mPoints[newverts[af*3 + n1]].v;
621 #if RECALCNORMALS==0
622                                         //norm_nbdyavg += mPoints[newverts[af*3 + n1]].n;
623 #endif // RECALCNORMALS==0
624                                         nnbdy++;
625                                 }
626                                 if (mAcrossEdge[af*3 + n2] == -1) {
627                                         bdyavg += mPoints[newverts[af*3 + n2]].v;
628 #if RECALCNORMALS==0
629                                         //norm_bdyavg += mPoints[newverts[af*3 + n2]].n;
630 #endif // RECALCNORMALS==0
631                                         nbdy++;
632                                 } else {
633                                         nbdyavg += mPoints[newverts[af*3 + n2]].v;
634 #if RECALCNORMALS==0
635                                         //norm_nbdyavg += mPoints[newverts[af*3 + n2]].n;
636 #endif // RECALCNORMALS==0
637                                         nnbdy++;
638                                 }
639                         }
640
641                         float alpha;
642                         ntlVec3Gfx newpt;
643                         if (nbdy) {
644                                 newpt = bdyavg / (float) nbdy;
645                                 alpha = 0.5f;
646                         } else if (nnbdy) {
647                                 newpt = nbdyavg / (float) nnbdy;
648                                 if (nnbdy == 6)
649                                         alpha = 1.05;
650                                 else if (nnbdy == 8)
651                                         alpha = 0.86;
652                                 else if (nnbdy == 10)
653                                         alpha = 0.7;
654                                 else
655                                         alpha = 0.6;
656                         } else {
657                                 continue;
658                         }
659                         mPoints[i].v *= 1.0f - alpha;
660                         mPoints[i].v += newpt * alpha;
661
662 #if RECALCNORMALS==0
663                         //mPoints[i].n *= 1.0f - alpha;
664                         //mPoints[i].n += newpt * alpha;
665 #endif // RECALCNORMALS==0
666         }
667
668         // Insert new faces
669         mIndices.reserve(4*nf);
670         for (i = 0; i < nf; i++) {
671                 mIndices.push_back( mIndices[i*3 + 0]);
672                 mIndices.push_back( newverts[i*3 + 2]);
673                 mIndices.push_back( newverts[i*3 + 1]);
674
675                 mIndices.push_back( mIndices[i*3 + 1]);
676                 mIndices.push_back( newverts[i*3 + 0]);
677                 mIndices.push_back( newverts[i*3 + 2]);
678
679                 mIndices.push_back( mIndices[i*3 + 2]);
680                 mIndices.push_back( newverts[i*3 + 1]);
681                 mIndices.push_back( newverts[i*3 + 0]);
682
683                 mIndices[i*3+0] = newverts[i*3+0];
684                 mIndices[i*3+1] = newverts[i*3+1];
685                 mIndices[i*3+2] = newverts[i*3+2];
686         }
687
688         // recalc normals
689 #if RECALCNORMALS==1
690         {
691                 int nf = (int)mIndices.size()/3, nv = (int)mPoints.size(), i;
692                 for (i = 0; i < nv; i++) {
693                         mPoints[i].n = ntlVec3Gfx(0.0);
694                 }
695                 for (i = 0; i < nf; i++) {
696                         const ntlVec3Gfx &p0 = mPoints[mIndices[i*3+0]].v;
697                         const ntlVec3Gfx &p1 = mPoints[mIndices[i*3+1]].v;
698                         const ntlVec3Gfx &p2 = mPoints[mIndices[i*3+2]].v;
699                         ntlVec3Gfx a = p0-p1, b = p1-p2, c = p2-p0;
700                         float l2a = normNoSqrt(a), l2b = normNoSqrt(b), l2c = normNoSqrt(c);
701
702                         ntlVec3Gfx facenormal = cross(a, b);
703
704                         mPoints[mIndices[i*3+0]].n += facenormal * (1.0f / (l2a * l2c));
705                         mPoints[mIndices[i*3+1]].n += facenormal * (1.0f / (l2b * l2a));
706                         mPoints[mIndices[i*3+2]].n += facenormal * (1.0f / (l2c * l2b));
707                 }
708
709                 for (i = 0; i < nv; i++) {
710                         normalize(mPoints[i].n);
711                 }
712         }
713 #else // RECALCNORMALS==1
714                 for (i = 0; i < (int)mPoints.size(); i++) {
715                         normalize(mPoints[i].n);
716                 }
717 #endif // RECALCNORMALS==1
718
719         //errMsg("SUBDIV","done nv:"<<mPoints.size()<<" nf:"<<mIndices.size() );
720 }
721
722
723
724
725
726
727
728 /******************************************************************************
729  * 
730  * Surface improvement
731  * makes use of trimesh2 library
732  * http://www.cs.princeton.edu/gfx/proj/trimesh2/
733  *
734  * Copyright (c) 2004 Szymon Rusinkiewicz.
735  * see COPYING_trimesh2
736  * 
737  *****************************************************************************/
738
739
740
741 // Diffuse a vector field at 1 vertex, weighted by
742 // a gaussian of width 1/sqrt(invsigma2)
743 void IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int v, float invsigma2, ntlVec3Gfx &flt)
744 {
745         flt = ntlVec3Gfx(0.0);
746         flt += *(field+pointerScale*v) *pointareas[v];
747         float sum_w = pointareas[v];
748         const ntlVec3Gfx &nv = mPoints[v].n;
749
750         unsigned &flag = flag_curr;
751         flag++;
752         flags[v] = flag;
753         vector<int> boundary = neighbors[v];
754         while (!boundary.empty()) {
755                 const int bbn = boundary.back();
756                 boundary.pop_back();
757                 if (flags[bbn] == flag) continue;
758                 flags[bbn] = flag;
759
760                 // gaussian weight of width 1/sqrt(invsigma2)
761                 const float d2 = invsigma2 * normNoSqrt(mPoints[bbn].v - mPoints[v].v);
762                 if(d2 >= 9.0f) continue; // 25 also possible  , slower
763                 //float w = (d2 >=  9.0f) ? 0.0f : exp(-0.5f*d2);
764                 float w = exp(-0.5f*d2);
765                 if(dot(nv, mPoints[bbn].n) <= 0.0f) continue; // faster than before d2 calc?
766
767                 // Downweight things pointing in different directions
768                 w *= dot(nv , mPoints[bbn].n);
769                 // Surface area "belonging" to each point
770                 w *= pointareas[bbn];
771                 // Accumulate weight times field at neighbor
772                 flt += *(field+pointerScale*bbn)*w;
773
774                 sum_w += w;
775                 for (int i = 0; i < (int)neighbors[bbn].size(); i++) {
776                         int nn = neighbors[bbn][i];
777                         if (flags[nn] == flag) continue;
778                         boundary.push_back(nn);
779                 }
780         }
781         flt /= sum_w;
782 }
783
784
785
786 void IsoSurface::smoothSurface(float sigma)
787 {
788         int nv = mPoints.size();
789         if ((int)flags.size() != nv) flags.resize(nv);
790         int nf = mIndices.size()/3;
791
792         { // need neighbors
793
794                 vector<int> numneighbors(mPoints.size());
795                 int i;
796                 for (i = 0; i < (int)mIndices.size()/3; i++) {
797                         numneighbors[mIndices[i*3+0]]++;
798                         numneighbors[mIndices[i*3+1]]++;
799                         numneighbors[mIndices[i*3+2]]++;
800                 }
801
802                 neighbors.clear();
803                 neighbors.resize(mPoints.size());
804                 for (i = 0; i < (int)mPoints.size(); i++) {
805                         neighbors[i].clear();
806                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
807                 }
808
809                 for (i = 0; i < (int)mIndices.size()/3; i++) {
810                         for (int j = 0; j < 3; j++) {
811                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
812                                 int n1 =  mIndices[i*3+((j+1)%3)];
813                                 int n2 =  mIndices[i*3+((j+2)%3)];
814                                 if (std::find(me.begin(), me.end(), n1) == me.end())
815                                         me.push_back(n1);
816                                 if (std::find(me.begin(), me.end(), n2) == me.end())
817                                         me.push_back(n2);
818                         }
819                 }
820         } // need neighbor
821
822         { // need pointarea
823                 pointareas.clear();
824                 pointareas.resize(nv);
825                 cornerareas.clear();
826                 cornerareas.resize(nf);
827
828                 for (int i = 0; i < nf; i++) {
829                         // Edges
830                         ntlVec3Gfx e[3] = { 
831                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
832                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
833                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
834
835                         // Compute corner weights
836                         float area = 0.5f * norm( cross(e[0], e[1]));
837                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
838                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
839                                         l2[1] * (l2[2] + l2[0] - l2[1]),
840                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
841                         if (ew[0] <= 0.0f) {
842                                 cornerareas[i][1] = -0.25f * l2[2] * area /
843                                                                 dot(e[0] , e[2]);
844                                 cornerareas[i][2] = -0.25f * l2[1] * area /
845                                                                 dot(e[0] , e[1]);
846                                 cornerareas[i][0] = area - cornerareas[i][1] -
847                                                                 cornerareas[i][2];
848                         } else if (ew[1] <= 0.0f) {
849                                 cornerareas[i][2] = -0.25f * l2[0] * area /
850                                                                 dot(e[1] , e[0]);
851                                 cornerareas[i][0] = -0.25f * l2[2] * area /
852                                                                 dot(e[1] , e[2]);
853                                 cornerareas[i][1] = area - cornerareas[i][2] -
854                                                                 cornerareas[i][0];
855                         } else if (ew[2] <= 0.0f) {
856                                 cornerareas[i][0] = -0.25f * l2[1] * area /
857                                                                 dot(e[2] , e[1]);
858                                 cornerareas[i][1] = -0.25f * l2[0] * area /
859                                                                 dot(e[2] , e[0]);
860                                 cornerareas[i][2] = area - cornerareas[i][0] -
861                                                                 cornerareas[i][1];
862                         } else {
863                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
864                                 for (int j = 0; j < 3; j++)
865                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
866                                                                                          ew[(j+2)%3]);
867                         }
868
869                         // NT important, check this...
870 #ifndef WIN32
871                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
872                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
873                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
874 #else // WIN32
875                         // FIXME check as well...
876                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
877                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
878                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
879 #endif // WIN32
880
881                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
882                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
883                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
884                 }
885
886         } // need pointarea
887         // */
888
889         float invsigma2 = 1.0f / (sigma*sigma);
890
891         vector<ntlVec3Gfx> dflt(nv);
892         for (int i = 0; i < nv; i++) {
893                 diffuseVertexField( &mPoints[0].v, 2,
894                                    i, invsigma2, dflt[i]);
895                 // Just keep the displacement
896                 dflt[i] -= mPoints[i].v;
897         }
898
899         // Slightly better small-neighborhood approximation
900         for (int i = 0; i < nf; i++) {
901                 ntlVec3Gfx c = mPoints[mIndices[i*3+0]].v +
902                           mPoints[mIndices[i*3+1]].v +
903                           mPoints[mIndices[i*3+2]].v;
904                 c /= 3.0f;
905                 for (int j = 0; j < 3; j++) {
906                         int v = mIndices[i*3+j];
907                         ntlVec3Gfx d =(c - mPoints[v].v) * 0.5f;
908                         dflt[v] += d * (cornerareas[i][j] /
909                                    pointareas[mIndices[i*3+j]] *
910                                    exp(-0.5f * invsigma2 * normNoSqrt(d)) );
911                 }
912         }
913
914         // Filter displacement field
915         vector<ntlVec3Gfx> dflt2(nv);
916         for (int i = 0; i < nv; i++) {
917                 diffuseVertexField( &dflt[0], 1,
918                                    i, invsigma2, dflt2[i]);
919         }
920
921         // Update vertex positions
922         for (int i = 0; i < nv; i++) {
923                 mPoints[i].v += dflt[i] - dflt2[i]; // second Laplacian
924         }
925
926         // when normals smoothing off, this cleans up quite well
927         // costs ca. 50% additional time though
928         float nsFac = 1.5f;
929         { float ninvsigma2 = 1.0f / (nsFac*nsFac*sigma*sigma);
930                 for (int i = 0; i < nv; i++) {
931                         diffuseVertexField( &mPoints[0].n, 2, i, ninvsigma2, dflt[i]);
932                         normalize(dflt[i]);
933                 }
934                 for (int i = 0; i < nv; i++) {
935                         mPoints[i].n = dflt[i];
936                 } 
937         } // smoothNormals copy */
938
939         //errMsg("SMSURF","done v:"<<sigma); // DEBUG
940 }
941
942 void IsoSurface::smoothNormals(float sigma)
943 {
944         { // need neighbor
945                 vector<int> numneighbors(mPoints.size());
946                 int i;
947                 for (i = 0; i < (int)mIndices.size()/3; i++) {
948                         numneighbors[mIndices[i*3+0]]++;
949                         numneighbors[mIndices[i*3+1]]++;
950                         numneighbors[mIndices[i*3+2]]++;
951                 }
952
953                 neighbors.clear();
954                 neighbors.resize(mPoints.size());
955                 for (i = 0; i < (int)mPoints.size(); i++) {
956                         neighbors[i].clear();
957                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
958                 }
959
960                 for (i = 0; i < (int)mIndices.size()/3; i++) {
961                         for (int j = 0; j < 3; j++) {
962                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
963                                 int n1 =  mIndices[i*3+((j+1)%3)];
964                                 int n2 =  mIndices[i*3+((j+2)%3)];
965                                 if (std::find(me.begin(), me.end(), n1) == me.end())
966                                         me.push_back(n1);
967                                 if (std::find(me.begin(), me.end(), n2) == me.end())
968                                         me.push_back(n2);
969                         }
970                 }
971         } // need neighbor
972
973         { // need pointarea
974                 int nf = mIndices.size()/3, nv = mPoints.size();
975                 pointareas.clear();
976                 pointareas.resize(nv);
977                 cornerareas.clear();
978                 cornerareas.resize(nf);
979
980                 for (int i = 0; i < nf; i++) {
981                         // Edges
982                         ntlVec3Gfx e[3] = { 
983                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
984                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
985                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
986
987                         // Compute corner weights
988                         float area = 0.5f * norm( cross(e[0], e[1]));
989                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
990                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
991                                         l2[1] * (l2[2] + l2[0] - l2[1]),
992                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
993                         if (ew[0] <= 0.0f) {
994                                 cornerareas[i][1] = -0.25f * l2[2] * area /
995                                                                 dot(e[0] , e[2]);
996                                 cornerareas[i][2] = -0.25f * l2[1] * area /
997                                                                 dot(e[0] , e[1]);
998                                 cornerareas[i][0] = area - cornerareas[i][1] -
999                                                                 cornerareas[i][2];
1000                         } else if (ew[1] <= 0.0f) {
1001                                 cornerareas[i][2] = -0.25f * l2[0] * area /
1002                                                                 dot(e[1] , e[0]);
1003                                 cornerareas[i][0] = -0.25f * l2[2] * area /
1004                                                                 dot(e[1] , e[2]);
1005                                 cornerareas[i][1] = area - cornerareas[i][2] -
1006                                                                 cornerareas[i][0];
1007                         } else if (ew[2] <= 0.0f) {
1008                                 cornerareas[i][0] = -0.25f * l2[1] * area /
1009                                                                 dot(e[2] , e[1]);
1010                                 cornerareas[i][1] = -0.25f * l2[0] * area /
1011                                                                 dot(e[2] , e[0]);
1012                                 cornerareas[i][2] = area - cornerareas[i][0] -
1013                                                                 cornerareas[i][1];
1014                         } else {
1015                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
1016                                 for (int j = 0; j < 3; j++)
1017                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
1018                                                                                          ew[(j+2)%3]);
1019                         }
1020
1021                         // NT important, check this...
1022 #ifndef WIN32
1023                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
1024                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
1025                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
1026 #else // WIN32
1027                         // FIXME check as well...
1028                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
1029                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
1030                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
1031 #endif // WIN32
1032
1033                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
1034                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
1035                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
1036                 }
1037
1038         } // need pointarea
1039
1040         int nv = mPoints.size();
1041         if ((int)flags.size() != nv) flags.resize(nv);
1042         float invsigma2 = 1.0f / (sigma*sigma);
1043
1044         vector<ntlVec3Gfx> nflt(nv);
1045         for (int i = 0; i < nv; i++) {
1046                 diffuseVertexField( &mPoints[0].n, 2, i, invsigma2, nflt[i]);
1047                 normalize(nflt[i]);
1048         }
1049
1050         for (int i = 0; i < nv; i++) {
1051                 mPoints[i].n = nflt[i];
1052         }
1053
1054         //errMsg("SMNRMLS","done v:"<<sigma); // DEBUG
1055 }
1056
1057