- merged latest version of fluid solver
[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                 errFatal("IsoSurface::triangulate","no LBM object, and no scalar field...!",SIMWORLD_INITERROR);
121                 return;
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         debMsgStd("IsoSurface::triangulate",DM_MSG,"Took "<< ((tritimeend-tritimestart)/(double)1000.0)<<"s) " , 10 );
336 }
337
338
339         
340
341
342 /******************************************************************************
343  * Get triangles for rendering
344  *****************************************************************************/
345 void IsoSurface::getTriangles( vector<ntlTriangle> *triangles, 
346                                                                                                          vector<ntlVec3Gfx> *vertices, 
347                                                                                                          vector<ntlVec3Gfx> *normals, int objectId )
348 {
349         if(!mInitDone) {
350                 debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10);
351                 return;
352         }
353
354   /* triangulate field */
355   triangulate();
356         //errMsg("TRIS"," "<<mIndices.size() );
357
358         // new output with vertice reuse
359         int iniVertIndex = (*vertices).size();
360         int iniNormIndex = (*normals).size();
361         if(iniVertIndex != iniNormIndex) {
362                 errFatal("getTriangles Error","For '"<<mName<<"': Vertices and normal array sizes to not match!!!",SIMWORLD_GENERICERROR);
363                 return; 
364         }
365         //errMsg("NM"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
366         //errMsg("NM"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
367
368   for(int i=0;i<(int)mPoints.size();i++) {
369                 vertices->push_back( mPoints[i].v );
370         }
371   for(int i=0;i<(int)mPoints.size();i++) {
372                 normals->push_back( mPoints[i].n );
373         }
374
375         //errMsg("N2"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
376         //errMsg("N2"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
377
378   for(int i=0;i<(int)mIndices.size();i+=3) {
379                 const int smooth = 1;
380     int t1 = mIndices[i];
381     int t2 = mIndices[i+1];
382                 int t3 = mIndices[i+2];
383                 //errMsg("NM"," tri"<<t1<<" "<<t2<<" "<<t3 );
384
385                 ntlTriangle tri;
386
387                 tri.getPoints()[0] = t1+iniVertIndex;
388                 tri.getPoints()[1] = t2+iniVertIndex;
389                 tri.getPoints()[2] = t3+iniVertIndex;
390
391                 /* init flags */
392                 int flag = 0; 
393                 if(getVisible()){ flag |= TRI_GEOMETRY; }
394                 if(getCastShadows() ) { 
395                         flag |= TRI_CASTSHADOWS; } 
396                 if( (getMaterial()->getMirror()>0.0) ||  
397                                 (getMaterial()->getTransparence()>0.0) ||  
398                                 (getMaterial()->getFresnel()>0.0) ) { 
399                         flag |= TRI_MAKECAUSTICS; } 
400                 else { 
401                         flag |= TRI_NOCAUSTICS; } 
402
403                 /* init geo init id */
404                 int geoiId = getGeoInitId(); 
405                 if(geoiId > 0) { 
406                         flag |= (1<< (geoiId+4)); 
407                         flag |= mGeoInitType; 
408                 } 
409
410                 tri.setFlags( flag );
411
412                 /* triangle normal missing */
413                 tri.setNormal( ntlVec3Gfx(0.0) );
414                 tri.setSmoothNormals( smooth );
415                 tri.setObjectId( objectId );
416                 triangles->push_back( tri ); 
417         }
418         //errMsg("N3"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
419         return;
420 }
421
422
423
424 inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) {
425         // WARNING - this assumes a security boundary layer... 
426         /*
427         if(i<=0) i=1;
428         if(j<=0) j=1;
429         if(k<=0) k=1;
430         if(i>=(mSizex-1)) i=mSizex-2;
431         if(j>=(mSizex-1)) j=mSizex-2;
432         if(k>=(mSizex-1)) k=mSizex-2; // */
433
434         ntlVec3Gfx ret(0.0);
435         ret[0] = *getData(i-1,j  ,k  ) - 
436                  *getData(i+1,j  ,k  );
437         ret[1] = *getData(i  ,j-1,k  ) - 
438                  *getData(i  ,j+1,k  );
439         ret[2] = *getData(i  ,j  ,k-1  ) - 
440                  *getData(i  ,j  ,k+1  );
441         return ret;
442 }
443
444 #define RECALCNORMALS 0
445
446 // Subdivide a mesh allways loop
447 void IsoSurface::subdivide()
448 {
449         mAdjacentFaces.clear(); 
450         mAcrossEdge.clear();
451
452         //void TriMesh::need_adjacentfaces()
453         {
454                 vector<int> numadjacentfaces(mPoints.size());
455                 //errMsg("SUBDIV ADJFA1", " "<<mPoints.size()<<" - "<<numadjacentfaces.size() );
456                 for (int i = 0; i < (int)mIndices.size()/3; i++) {
457                         numadjacentfaces[mIndices[i*3 + 0]]++;
458                         numadjacentfaces[mIndices[i*3 + 1]]++;
459                         numadjacentfaces[mIndices[i*3 + 2]]++;
460                 }
461
462                 mAdjacentFaces.resize(mPoints.size());
463                 for (int i = 0; i < (int)mPoints.size(); i++)
464                         mAdjacentFaces[i].reserve(numadjacentfaces[i]);
465
466                 for (int i = 0; i < (int)mIndices.size()/3; i++) {
467                         for (int j = 0; j < 3; j++)
468                                 mAdjacentFaces[mIndices[i*3 + j]].push_back(i);
469                 }
470
471         }
472
473         // Find the face across each edge from each other face (-1 on boundary)
474         // If topology is bad, not necessarily what one would expect...
475         //void TriMesh::need_across_edge()
476         {
477                 mAcrossEdge.resize(mIndices.size(), -1);
478
479                 for (int i = 0; i < (int)mIndices.size()/3; i++) {
480                         for (int j = 0; j < 3; j++) {
481                                 if (mAcrossEdge[i*3 + j] != -1)
482                                         continue;
483                                 int v1 = mIndices[i*3 + ((j+1)%3)];
484                                 int v2 = mIndices[i*3 + ((j+2)%3)];
485                                 const vector<int> &a1 = mAdjacentFaces[v1];
486                                 const vector<int> &a2 = mAdjacentFaces[v2];
487                                 for (int k1 = 0; k1 < (int)a1.size(); k1++) {
488                                         int other = a1[k1];
489                                         if (other == i)
490                                                 continue;
491                                         vector<int>::const_iterator it =
492                                                 std::find(a2.begin(), a2.end(), other);
493                                         if (it == a2.end())
494                                                 continue;
495
496                                         //int ind = (faces[other].indexof(v1)+1)%3;
497                                         int ind = -1;
498                                         if( mIndices[other*3+0] == (unsigned int)v1 ) ind = 0;
499                                         else if( mIndices[other*3+1] == (unsigned int)v1 ) ind = 1;
500                                         else if( mIndices[other*3+2] == (unsigned int)v1 ) ind = 2;
501                                         ind = (ind+1)%3;
502
503                                         if ( (int)mIndices[other*3 + ((ind+1)%3)] != v2)
504                                                 continue;
505                                         mAcrossEdge[i*3 + j] = other;
506                                         mAcrossEdge[other*3 + ind] = i;
507                                         break;
508                                 }
509                         }
510                 }
511
512                 //errMsg("SUBDIV ACREDG", "Done.\n");
513         }
514
515         //errMsg("SUBDIV","start");
516         // Introduce new vertices
517         int nf = (int)mIndices.size() / 3;
518
519         //vector<TriMesh::Face> newverts(nf, TriMesh::Face(-1,-1,-1));
520         vector<int> newverts(nf*3); //, TriMesh::Face(-1,-1,-1));
521         for(int j=0; j<(int)newverts.size(); j++) newverts[j] = -1;
522
523         int old_nv = (int)mPoints.size();
524         mPoints.reserve(4 * old_nv);
525         vector<int> newvert_count(old_nv + 3*nf); // wichtig...?
526         //errMsg("NC", newvert_count.size() );
527
528         for (int i = 0; i < nf; i++) {
529                 for (int j = 0; j < 3; j++) {
530                         int ae = mAcrossEdge[i*3 + j];
531                         if (newverts[i*3 + j] == -1 && ae != -1) {
532                                 if (mAcrossEdge[ae*3 + 0] == i)
533                                         newverts[i*3 + j] = newverts[ae*3 + 0];
534                                 else if (mAcrossEdge[ae*3 + 1] == i)
535                                         newverts[i*3 + j] = newverts[ae*3 + 1];
536                                 else if (mAcrossEdge[ae*3 + 2] == i)
537                                         newverts[i*3 + j] = newverts[ae*3 + 2];
538                         }
539                         if (newverts[i*3 + j] == -1) {
540                                 IsoLevelVertex ilv;
541                                 ilv.v = ntlVec3Gfx(0.0);
542                                 ilv.n = ntlVec3Gfx(0.0);
543                                 mPoints.push_back(ilv);
544                                 newverts[i*3 + j] = (int)mPoints.size() - 1;
545                                 if (ae != -1) {
546                                         if (mAcrossEdge[ae*3 + 0] == i)
547                                                 newverts[ae*3 + 0] = newverts[i*3 + j];
548                                         else if (mAcrossEdge[ae*3 + 1] == i)
549                                                 newverts[ae*3 + 1] = newverts[i*3 + j];
550                                         else if (mAcrossEdge[ae*3 + 2] == i)
551                                                 newverts[ae*3 + 2] = newverts[i*3 + j];
552                                 }
553                         }
554                         if(ae != -1) {
555                                 mPoints[newverts[i*3 + j]].v +=
556                                         mPoints[ mIndices[i*3 + ( j     )] ].v * 0.25f  + // j = 0,1,2?
557                                         mPoints[ mIndices[i*3 + ((j+1)%3)] ].v * 0.375f +
558                                         mPoints[ mIndices[i*3 + ((j+2)%3)] ].v * 0.375f;
559 #if RECALCNORMALS==0
560                                 mPoints[newverts[i*3 + j]].n +=
561                                         mPoints[ mIndices[i*3 + ( j     )] ].n * 0.25f  + // j = 0,1,2?
562                                         mPoints[ mIndices[i*3 + ((j+1)%3)] ].n * 0.375f +
563                                         mPoints[ mIndices[i*3 + ((j+2)%3)] ].n * 0.375f;
564 #endif // RECALCNORMALS==0
565                         } else {
566                                 mPoints[newverts[i*3 + j]].v +=
567                                         mPoints[ mIndices[i*3 + ((j+1)%3)] ].v * 0.5f   +
568                                         mPoints[ mIndices[i*3 + ((j+2)%3)] ].v * 0.5f  ;
569 #if RECALCNORMALS==0
570                                 mPoints[newverts[i*3 + j]].n +=
571                                         mPoints[ mIndices[i*3 + ((j+1)%3)] ].n * 0.5f   +
572                                         mPoints[ mIndices[i*3 + ((j+2)%3)] ].n * 0.5f  ;
573 #endif // RECALCNORMALS==0
574                         }
575
576                         newvert_count[newverts[i*3 + j]]++;
577                 }
578         }
579         for (int i = old_nv; i < (int)mPoints.size(); i++) {
580                 if (!newvert_count[i])
581                         continue;
582                 float scale = 1.0f / newvert_count[i];
583                 mPoints[i].v *= scale;
584
585 #if RECALCNORMALS==0
586                 //mPoints[i].n *= scale;
587                 //normalize( mPoints[i].n );
588 #endif // RECALCNORMALS==0
589         }
590
591         // Update old vertices
592         for (int i = 0; i < old_nv; i++) {
593                         ntlVec3Gfx bdyavg(0.0), nbdyavg(0.0);
594                         ntlVec3Gfx norm_bdyavg(0.0), norm_nbdyavg(0.0); // N
595                         int nbdy = 0, nnbdy = 0;
596                         int naf = (int)mAdjacentFaces[i].size();
597                         if (!naf)
598                                 continue;
599                         for (int j = 0; j < naf; j++) {
600                                 int af = mAdjacentFaces[i][j];
601
602                                 int afi = -1;
603                                 if( mIndices[af*3+0] == (unsigned int)i ) afi = 0;
604                                 else if( mIndices[af*3+1] == (unsigned int)i ) afi = 1;
605                                 else if( mIndices[af*3+2] == (unsigned int)i ) afi = 2;
606
607                                 int n1 = (afi+1) % 3;
608                                 int n2 = (afi+2) % 3;
609                                 if (mAcrossEdge[af*3 + n1] == -1) {
610                                         bdyavg += mPoints[newverts[af*3 + n1]].v;
611 #if RECALCNORMALS==0
612                                         //norm_bdyavg += mPoints[newverts[af*3 + n1]].n;
613 #endif // RECALCNORMALS==0
614                                         nbdy++;
615                                 } else {
616                                         nbdyavg += mPoints[newverts[af*3 + n1]].v;
617 #if RECALCNORMALS==0
618                                         //norm_nbdyavg += mPoints[newverts[af*3 + n1]].n;
619 #endif // RECALCNORMALS==0
620                                         nnbdy++;
621                                 }
622                                 if (mAcrossEdge[af*3 + n2] == -1) {
623                                         bdyavg += mPoints[newverts[af*3 + n2]].v;
624 #if RECALCNORMALS==0
625                                         //norm_bdyavg += mPoints[newverts[af*3 + n2]].n;
626 #endif // RECALCNORMALS==0
627                                         nbdy++;
628                                 } else {
629                                         nbdyavg += mPoints[newverts[af*3 + n2]].v;
630 #if RECALCNORMALS==0
631                                         //norm_nbdyavg += mPoints[newverts[af*3 + n2]].n;
632 #endif // RECALCNORMALS==0
633                                         nnbdy++;
634                                 }
635                         }
636
637                         float alpha;
638                         ntlVec3Gfx newpt;
639                         if (nbdy) {
640                                 newpt = bdyavg / (float) nbdy;
641                                 alpha = 0.5f;
642                         } else if (nnbdy) {
643                                 newpt = nbdyavg / (float) nnbdy;
644                                 if (nnbdy == 6)
645                                         alpha = 1.05;
646                                 else if (nnbdy == 8)
647                                         alpha = 0.86;
648                                 else if (nnbdy == 10)
649                                         alpha = 0.7;
650                                 else
651                                         alpha = 0.6;
652                         } else {
653                                 continue;
654                         }
655                         mPoints[i].v *= 1.0f - alpha;
656                         mPoints[i].v += newpt * alpha;
657
658 #if RECALCNORMALS==0
659                         //mPoints[i].n *= 1.0f - alpha;
660                         //mPoints[i].n += newpt * alpha;
661 #endif // RECALCNORMALS==0
662         }
663
664         // Insert new faces
665         mIndices.reserve(4*nf);
666         for (int i = 0; i < nf; i++) {
667                 mIndices.push_back( mIndices[i*3 + 0]);
668                 mIndices.push_back( newverts[i*3 + 2]);
669                 mIndices.push_back( newverts[i*3 + 1]);
670
671                 mIndices.push_back( mIndices[i*3 + 1]);
672                 mIndices.push_back( newverts[i*3 + 0]);
673                 mIndices.push_back( newverts[i*3 + 2]);
674
675                 mIndices.push_back( mIndices[i*3 + 2]);
676                 mIndices.push_back( newverts[i*3 + 1]);
677                 mIndices.push_back( newverts[i*3 + 0]);
678
679                 mIndices[i*3+0] = newverts[i*3+0];
680                 mIndices[i*3+1] = newverts[i*3+1];
681                 mIndices[i*3+2] = newverts[i*3+2];
682         }
683
684         // recalc normals
685 #if RECALCNORMALS==1
686         {
687                 int nf = (int)mIndices.size()/3, nv = (int)mPoints.size();
688                 for (int i = 0; i < nv; i++) {
689                         mPoints[i].n = ntlVec3Gfx(0.0);
690                 }
691                 for (int i = 0; i < nf; i++) {
692                         const ntlVec3Gfx &p0 = mPoints[mIndices[i*3+0]].v;
693                         const ntlVec3Gfx &p1 = mPoints[mIndices[i*3+1]].v;
694                         const ntlVec3Gfx &p2 = mPoints[mIndices[i*3+2]].v;
695                         ntlVec3Gfx a = p0-p1, b = p1-p2, c = p2-p0;
696                         float l2a = normNoSqrt(a), l2b = normNoSqrt(b), l2c = normNoSqrt(c);
697
698                         ntlVec3Gfx facenormal = cross(a, b);
699
700                         mPoints[mIndices[i*3+0]].n += facenormal * (1.0f / (l2a * l2c));
701                         mPoints[mIndices[i*3+1]].n += facenormal * (1.0f / (l2b * l2a));
702                         mPoints[mIndices[i*3+2]].n += facenormal * (1.0f / (l2c * l2b));
703                 }
704
705                 for (int i = 0; i < nv; i++) {
706                         normalize(mPoints[i].n);
707                 }
708         }
709 #else // RECALCNORMALS==1
710                 for (int i = 0; i < (int)mPoints.size(); i++) {
711                         normalize(mPoints[i].n);
712                 }
713 #endif // RECALCNORMALS==1
714
715         //errMsg("SUBDIV","done nv:"<<mPoints.size()<<" nf:"<<mIndices.size() );
716 }
717
718
719
720
721
722
723
724 /******************************************************************************
725  * 
726  * Surface improvement
727  * makes use of trimesh2 library
728  * http://www.cs.princeton.edu/gfx/proj/trimesh2/
729  *
730  * Copyright (c) 2004 Szymon Rusinkiewicz.
731  * see COPYING_trimesh2
732  * 
733  *****************************************************************************/
734
735
736
737 // Diffuse a vector field at 1 vertex, weighted by
738 // a gaussian of width 1/sqrt(invsigma2)
739 void IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int v, float invsigma2, ntlVec3Gfx &flt)
740 {
741         flt = ntlVec3Gfx(0.0);
742         flt += *(field+pointerScale*v) *pointareas[v];
743         float sum_w = pointareas[v];
744         const ntlVec3Gfx &nv = mPoints[v].n;
745
746         unsigned &flag = flag_curr;
747         flag++;
748         flags[v] = flag;
749         vector<int> boundary = neighbors[v];
750         while (!boundary.empty()) {
751                 const int bbn = boundary.back();
752                 boundary.pop_back();
753                 if (flags[bbn] == flag) continue;
754                 flags[bbn] = flag;
755
756                 // gaussian weight of width 1/sqrt(invsigma2)
757                 const float d2 = invsigma2 * normNoSqrt(mPoints[bbn].v - mPoints[v].v);
758                 if(d2 >= 9.0f) continue; // 25 also possible  , slower
759                 //float w = (d2 >=  9.0f) ? 0.0f : exp(-0.5f*d2);
760                 float w = exp(-0.5f*d2);
761                 if(dot(nv, mPoints[bbn].n) <= 0.0f) continue; // faster than before d2 calc?
762
763                 // Downweight things pointing in different directions
764                 w *= dot(nv , mPoints[bbn].n);
765                 // Surface area "belonging" to each point
766                 w *= pointareas[bbn];
767                 // Accumulate weight times field at neighbor
768                 flt += *(field+pointerScale*bbn)*w;
769
770                 sum_w += w;
771                 for (int i = 0; i < (int)neighbors[bbn].size(); i++) {
772                         int nn = neighbors[bbn][i];
773                         if (flags[nn] == flag) continue;
774                         boundary.push_back(nn);
775                 }
776         }
777         flt /= sum_w;
778 }
779
780
781
782 void IsoSurface::smoothSurface(float sigma)
783 {
784         int nv = mPoints.size();
785         if ((int)flags.size() != nv) flags.resize(nv);
786         int nf = mIndices.size()/3;
787
788         { // need neighbors
789
790                 vector<int> numneighbors(mPoints.size());
791                 int i;
792                 for (i = 0; i < (int)mIndices.size()/3; i++) {
793                         numneighbors[mIndices[i*3+0]]++;
794                         numneighbors[mIndices[i*3+1]]++;
795                         numneighbors[mIndices[i*3+2]]++;
796                 }
797
798                 neighbors.clear();
799                 neighbors.resize(mPoints.size());
800                 for (i = 0; i < (int)mPoints.size(); i++) {
801                         neighbors[i].clear();
802                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
803                 }
804
805                 for (i = 0; i < (int)mIndices.size()/3; i++) {
806                         for (int j = 0; j < 3; j++) {
807                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
808                                 int n1 =  mIndices[i*3+((j+1)%3)];
809                                 int n2 =  mIndices[i*3+((j+2)%3)];
810                                 if (std::find(me.begin(), me.end(), n1) == me.end())
811                                         me.push_back(n1);
812                                 if (std::find(me.begin(), me.end(), n2) == me.end())
813                                         me.push_back(n2);
814                         }
815                 }
816         } // need neighbor
817
818         { // need pointarea
819                 pointareas.clear();
820                 pointareas.resize(nv);
821                 cornerareas.clear();
822                 cornerareas.resize(nf);
823
824                 for (int i = 0; i < nf; i++) {
825                         // Edges
826                         ntlVec3Gfx e[3] = { 
827                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
828                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
829                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
830
831                         // Compute corner weights
832                         float area = 0.5f * norm( cross(e[0], e[1]));
833                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
834                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
835                                         l2[1] * (l2[2] + l2[0] - l2[1]),
836                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
837                         if (ew[0] <= 0.0f) {
838                                 cornerareas[i][1] = -0.25f * l2[2] * area /
839                                                                 dot(e[0] , e[2]);
840                                 cornerareas[i][2] = -0.25f * l2[1] * area /
841                                                                 dot(e[0] , e[1]);
842                                 cornerareas[i][0] = area - cornerareas[i][1] -
843                                                                 cornerareas[i][2];
844                         } else if (ew[1] <= 0.0f) {
845                                 cornerareas[i][2] = -0.25f * l2[0] * area /
846                                                                 dot(e[1] , e[0]);
847                                 cornerareas[i][0] = -0.25f * l2[2] * area /
848                                                                 dot(e[1] , e[2]);
849                                 cornerareas[i][1] = area - cornerareas[i][2] -
850                                                                 cornerareas[i][0];
851                         } else if (ew[2] <= 0.0f) {
852                                 cornerareas[i][0] = -0.25f * l2[1] * area /
853                                                                 dot(e[2] , e[1]);
854                                 cornerareas[i][1] = -0.25f * l2[0] * area /
855                                                                 dot(e[2] , e[0]);
856                                 cornerareas[i][2] = area - cornerareas[i][0] -
857                                                                 cornerareas[i][1];
858                         } else {
859                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
860                                 for (int j = 0; j < 3; j++)
861                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
862                                                                                          ew[(j+2)%3]);
863                         }
864
865                         // NT important, check this...
866 #ifndef WIN32
867                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
868                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
869                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
870 #else // WIN32
871                         // FIXME check as well...
872                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
873                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
874                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
875 #endif // WIN32
876
877                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
878                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
879                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
880                 }
881
882         } // need pointarea
883         // */
884
885         float invsigma2 = 1.0f / (sigma*sigma);
886
887         vector<ntlVec3Gfx> dflt(nv);
888         for (int i = 0; i < nv; i++) {
889                 diffuseVertexField( &mPoints[0].v, 2,
890                                    i, invsigma2, dflt[i]);
891                 // Just keep the displacement
892                 dflt[i] -= mPoints[i].v;
893         }
894
895         // Slightly better small-neighborhood approximation
896         for (int i = 0; i < nf; i++) {
897                 ntlVec3Gfx c = mPoints[mIndices[i*3+0]].v +
898                           mPoints[mIndices[i*3+1]].v +
899                           mPoints[mIndices[i*3+2]].v;
900                 c /= 3.0f;
901                 for (int j = 0; j < 3; j++) {
902                         int v = mIndices[i*3+j];
903                         ntlVec3Gfx d =(c - mPoints[v].v) * 0.5f;
904                         dflt[v] += d * (cornerareas[i][j] /
905                                    pointareas[mIndices[i*3+j]] *
906                                    exp(-0.5f * invsigma2 * normNoSqrt(d)) );
907                 }
908         }
909
910         // Filter displacement field
911         vector<ntlVec3Gfx> dflt2(nv);
912         for (int i = 0; i < nv; i++) {
913                 diffuseVertexField( &dflt[0], 1,
914                                    i, invsigma2, dflt2[i]);
915         }
916
917         // Update vertex positions
918         for (int i = 0; i < nv; i++) {
919                 mPoints[i].v += dflt[i] - dflt2[i]; // second Laplacian
920         }
921
922         // when normals smoothing off, this cleans up quite well
923         // costs ca. 50% additional time though
924         float nsFac = 1.5f;
925         { float ninvsigma2 = 1.0f / (nsFac*nsFac*sigma*sigma);
926                 for (int i = 0; i < nv; i++) {
927                         diffuseVertexField( &mPoints[0].n, 2, i, ninvsigma2, dflt[i]);
928                         normalize(dflt[i]);
929                 }
930                 for (int i = 0; i < nv; i++) {
931                         mPoints[i].n = dflt[i];
932                 } 
933         } // smoothNormals copy */
934
935         //errMsg("SMSURF","done v:"<<sigma); // DEBUG
936 }
937
938 void IsoSurface::smoothNormals(float sigma)
939 {
940         { // need neighbor
941                 vector<int> numneighbors(mPoints.size());
942                 int i;
943                 for (i = 0; i < (int)mIndices.size()/3; i++) {
944                         numneighbors[mIndices[i*3+0]]++;
945                         numneighbors[mIndices[i*3+1]]++;
946                         numneighbors[mIndices[i*3+2]]++;
947                 }
948
949                 neighbors.clear();
950                 neighbors.resize(mPoints.size());
951                 for (i = 0; i < (int)mPoints.size(); i++) {
952                         neighbors[i].clear();
953                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
954                 }
955
956                 for (i = 0; i < (int)mIndices.size()/3; i++) {
957                         for (int j = 0; j < 3; j++) {
958                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
959                                 int n1 =  mIndices[i*3+((j+1)%3)];
960                                 int n2 =  mIndices[i*3+((j+2)%3)];
961                                 if (std::find(me.begin(), me.end(), n1) == me.end())
962                                         me.push_back(n1);
963                                 if (std::find(me.begin(), me.end(), n2) == me.end())
964                                         me.push_back(n2);
965                         }
966                 }
967         } // need neighbor
968
969         { // need pointarea
970                 int nf = mIndices.size()/3, nv = mPoints.size();
971                 pointareas.clear();
972                 pointareas.resize(nv);
973                 cornerareas.clear();
974                 cornerareas.resize(nf);
975
976                 for (int i = 0; i < nf; i++) {
977                         // Edges
978                         ntlVec3Gfx e[3] = { 
979                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
980                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
981                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
982
983                         // Compute corner weights
984                         float area = 0.5f * norm( cross(e[0], e[1]));
985                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
986                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
987                                         l2[1] * (l2[2] + l2[0] - l2[1]),
988                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
989                         if (ew[0] <= 0.0f) {
990                                 cornerareas[i][1] = -0.25f * l2[2] * area /
991                                                                 dot(e[0] , e[2]);
992                                 cornerareas[i][2] = -0.25f * l2[1] * area /
993                                                                 dot(e[0] , e[1]);
994                                 cornerareas[i][0] = area - cornerareas[i][1] -
995                                                                 cornerareas[i][2];
996                         } else if (ew[1] <= 0.0f) {
997                                 cornerareas[i][2] = -0.25f * l2[0] * area /
998                                                                 dot(e[1] , e[0]);
999                                 cornerareas[i][0] = -0.25f * l2[2] * area /
1000                                                                 dot(e[1] , e[2]);
1001                                 cornerareas[i][1] = area - cornerareas[i][2] -
1002                                                                 cornerareas[i][0];
1003                         } else if (ew[2] <= 0.0f) {
1004                                 cornerareas[i][0] = -0.25f * l2[1] * area /
1005                                                                 dot(e[2] , e[1]);
1006                                 cornerareas[i][1] = -0.25f * l2[0] * area /
1007                                                                 dot(e[2] , e[0]);
1008                                 cornerareas[i][2] = area - cornerareas[i][0] -
1009                                                                 cornerareas[i][1];
1010                         } else {
1011                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
1012                                 for (int j = 0; j < 3; j++)
1013                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
1014                                                                                          ew[(j+2)%3]);
1015                         }
1016
1017                         // NT important, check this...
1018 #ifndef WIN32
1019                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
1020                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
1021                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
1022 #else // WIN32
1023                         // FIXME check as well...
1024                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
1025                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
1026                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
1027 #endif // WIN32
1028
1029                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
1030                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
1031                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
1032                 }
1033
1034         } // need pointarea
1035
1036         int nv = mPoints.size();
1037         if ((int)flags.size() != nv) flags.resize(nv);
1038         float invsigma2 = 1.0f / (sigma*sigma);
1039
1040         vector<ntlVec3Gfx> nflt(nv);
1041         for (int i = 0; i < nv; i++) {
1042                 diffuseVertexField( &mPoints[0].n, 2, i, invsigma2, nflt[i]);
1043                 normalize(nflt[i]);
1044         }
1045
1046         for (int i = 0; i < nv; i++) {
1047                 mPoints[i].n = nflt[i];
1048         }
1049
1050         //errMsg("SMNRMLS","done v:"<<sigma); // DEBUG
1051 }
1052
1053