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