intern/elbeem/intern/isosurface.cpp

   1 /******************************************************************************
   2  *
   3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
   4  * Copyright 2003-2006 Nils Thuerey
   5  *
   6  * Marching Cubes surface mesh generation
   7  *
   8  *****************************************************************************/
   9
  10 #include "isosurface.h"
  11 #include "mcubes_tables.h"
  12 #include "particletracer.h"
  13 #include <algorithm>
  14 #include <stdio.h>
  15
  16 // just use default rounding for platforms where its not available
  17 #ifndef round
  18 #define round(x) (x)
  19 #endif
  20
  21 /******************************************************************************
  22  * Constructor
  23  *****************************************************************************/
  24 IsoSurface::IsoSurface(double iso) :
  25         ntlGeometryObject(),
  26         mSizex(-1), mSizey(-1), mSizez(-1),
  27         mpData(NULL),
  28   mIsoValue( iso ),
  29         mPoints(),
  30         mUseFullEdgeArrays(false),
  31         mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL),
  32         mEdgeArSize(-1),
  33         mIndices(),
  34
  35   mStart(0.0), mEnd(0.0), mDomainExtent(0.0),
  36   mInitDone(false),
  37         mSmoothSurface(0.0), mSmoothNormals(0.0),
  38         mAcrossEdge(), mAdjacentFaces(),
  39         mCutoff(-1), mCutArray(NULL), // off by default
  40         mpIsoParts(NULL), mPartSize(0.), mSubdivs(0),
  41         mFlagCnt(1),
  42         mSCrad1(0.), mSCrad2(0.), mSCcenter(0.)
  43 {
  44 }
  45
  46
  47 /******************************************************************************
  48  * The real init...
  49  *****************************************************************************/
  50 void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent)
  51 {
  52         // range 1-10 (max due to subd array in triangulate)
  53         if(mSubdivs<1) mSubdivs=1;
  54         if(mSubdivs>10) mSubdivs=10;
  55
  56         // init solver and size
  57         mSizex = setx;
  58         mSizey = sety;
  59         if(setz == 1) {// 2D, create thin 2D surface
  60                 setz = 5;
  61         }
  62         mSizez = setz;
  63         mDomainExtent = extent;
  64
  65         /* check triangulation size (for raytraing) */
  66         if( ( mStart[0] >= mEnd[0] ) && ( mStart[1] >= mEnd[1] ) && ( mStart[2] >= mEnd[2] ) ){
  67                 // extent was not set, use normalized one from parametrizer
  68                 mStart = ntlVec3Gfx(0.0) - extent*0.5;
  69                 mEnd = ntlVec3Gfx(0.0) + extent*0.5;
  70         }
  71
  72   // init
  73         mIndices.clear();
  74   mPoints.clear();
  75
  76         int nodes = mSizez*mSizey*mSizex;
  77   mpData = new float[nodes];
  78   for(int i=0;i<nodes;i++) { mpData[i] = 0.0; }
  79
  80   // allocate edge arrays  (last slices are never used...)
  81         int initsize = -1;
  82         if(mUseFullEdgeArrays) {
  83                 mEdgeArSize = nodes;
  84                 mpEdgeVerticesX = new int[nodes];
  85                 mpEdgeVerticesY = new int[nodes];
  86                 mpEdgeVerticesZ = new int[nodes];
  87                 initsize = 3*nodes;
  88         } else {
  89                 int sliceNodes = 2*mSizex*mSizey*mSubdivs*mSubdivs;
  90                 mEdgeArSize = sliceNodes;
  91                 mpEdgeVerticesX = new int[sliceNodes];
  92                 mpEdgeVerticesY = new int[sliceNodes];
  93                 mpEdgeVerticesZ = new int[sliceNodes];
  94                 initsize = 3*sliceNodes;
  95         }
  96   for(int i=0;i<mEdgeArSize;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
  97         // WARNING - make sure this is consistent with calculateMemreqEstimate
  98
  99         // marching cubes are ready
 100         mInitDone = true;
 101         debMsgStd("IsoSurface::initializeIsosurface",DM_MSG,"Inited, edgenodes:"<<initsize<<" subdivs:"<<mSubdivs<<" fulledg:"<<mUseFullEdgeArrays , 10);
 102 }
 103
 104
 105
 106 /*! Reset all values */
 107 void IsoSurface::resetAll(gfxReal val) {
 108         int nodes = mSizez*mSizey*mSizex;
 109   for(int i=0;i<nodes;i++) { mpData[i] = val; }
 110 }
 111
 112
 113 /******************************************************************************
 114  * Destructor
 115  *****************************************************************************/
 116 IsoSurface::~IsoSurface( void )
 117 {
 118         if(mpData) delete [] mpData;
 119         if(mpEdgeVerticesX) delete [] mpEdgeVerticesX;
 120         if(mpEdgeVerticesY) delete [] mpEdgeVerticesY;
 121         if(mpEdgeVerticesZ) delete [] mpEdgeVerticesZ;
 122 }
 123
 124
 125
 126
 127
 128 /******************************************************************************
 129  * triangulate the scalar field given by pointer
 130  *****************************************************************************/
 131 void IsoSurface::triangulate( void )
 132 {
 133   double gsx,gsy,gsz; // grid spacing in x,y,z direction
 134   double px,py,pz;    // current position in grid in x,y,z direction
 135   IsoLevelCube cubie;    // struct for a small subcube
 136         myTime_t tritimestart = getTime();
 137
 138         if(!mpData) {
 139                 errFatal("IsoSurface::triangulate","no LBM object, and no scalar field...!",SIMWORLD_INITERROR);
 140                 return;
 141         }
 142
 143   // get grid spacing (-2 to have same spacing as sim)
 144   gsx = (mEnd[0]-mStart[0])/(double)(mSizex-2.0);
 145   gsy = (mEnd[1]-mStart[1])/(double)(mSizey-2.0);
 146   gsz = (mEnd[2]-mStart[2])/(double)(mSizez-2.0);
 147
 148   // clean up previous frame
 149         mIndices.clear();
 150         mPoints.clear();
 151
 152         // reset edge vertices
 153   for(int i=0;i<mEdgeArSize;i++) {
 154                 mpEdgeVerticesX[i] = -1;
 155                 mpEdgeVerticesY[i] = -1;
 156                 mpEdgeVerticesZ[i] = -1;
 157         }
 158
 159         ntlVec3Gfx pos[8];
 160         float value[8];
 161         int cubeIndex;      // index entry of the cube
 162         int triIndices[12]; // vertex indices
 163         int *eVert[12];
 164         IsoLevelVertex ilv;
 165
 166         // edges between which points?
 167         const int mcEdges[24] = {
 168                 0,1,  1,2,  2,3,  3,0,
 169                 4,5,  5,6,  6,7,  7,4,
 170                 0,4,  1,5,  2,6,  3,7 };
 171
 172         const int cubieOffsetX[8] = {
 173                 0,1,1,0,  0,1,1,0 };
 174         const int cubieOffsetY[8] = {
 175                 0,0,1,1,  0,0,1,1 };
 176         const int cubieOffsetZ[8] = {
 177                 0,0,0,0,  1,1,1,1 };
 178
 179         const int coAdd=2;
 180   // let the cubes march
 181         if(mSubdivs<=1) {
 182
 183                 pz = mStart[2]-gsz*0.5;
 184                 for(int k=1;k<(mSizez-2);k++) {
 185                         pz += gsz;
 186                         py = mStart[1]-gsy*0.5;
 187                         for(int j=1;j<(mSizey-2);j++) {
 188                                 py += gsy;
 189                                 px = mStart[0]-gsx*0.5;
 190                                 for(int i=1;i<(mSizex-2);i++) {
 191                                         px += gsx;
 192
 193                                         value[0] = *getData(i  ,j  ,k  );
 194                                         value[1] = *getData(i+1,j  ,k  );
 195                                         value[2] = *getData(i+1,j+1,k  );
 196                                         value[3] = *getData(i  ,j+1,k  );
 197                                         value[4] = *getData(i  ,j  ,k+1);
 198                                         value[5] = *getData(i+1,j  ,k+1);
 199                                         value[6] = *getData(i+1,j+1,k+1);
 200                                         value[7] = *getData(i  ,j+1,k+1);
 201
 202                                         // check intersections of isosurface with edges, and calculate cubie index
 203                                         cubeIndex = 0;
 204                                         if (value[0] < mIsoValue) cubeIndex |= 1;
 205                                         if (value[1] < mIsoValue) cubeIndex |= 2;
 206                                         if (value[2] < mIsoValue) cubeIndex |= 4;
 207                                         if (value[3] < mIsoValue) cubeIndex |= 8;
 208                                         if (value[4] < mIsoValue) cubeIndex |= 16;
 209                                         if (value[5] < mIsoValue) cubeIndex |= 32;
 210                                         if (value[6] < mIsoValue) cubeIndex |= 64;
 211                                         if (value[7] < mIsoValue) cubeIndex |= 128;
 212
 213                                         // No triangles to generate?
 214                                         if (mcEdgeTable[cubeIndex] == 0) {
 215                                                 continue;
 216                                         }
 217
 218                                         // where to look up if this point already exists
 219                                         int edgek = 0;
 220                                         if(mUseFullEdgeArrays) edgek=k;
 221                                         const int baseIn = ISOLEVEL_INDEX( i+0, j+0, edgek+0);
 222                                         eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
 223                                         eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
 224                                         eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
 225                                         eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
 226
 227                                         eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
 228                                         eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, edgek+1) ];
 229                                         eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+1) ];
 230                                         eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
 231
 232                                         eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
 233                                         eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, edgek+0) ];
 234                                         eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, edgek+0) ];
 235                                         eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
 236
 237                                         // grid positions
 238                                         pos[0] = ntlVec3Gfx(px    ,py    ,pz);
 239                                         pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
 240                                         pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
 241                                         pos[3] = ntlVec3Gfx(px    ,py+gsy,pz);
 242                                         pos[4] = ntlVec3Gfx(px    ,py    ,pz+gsz);
 243                                         pos[5] = ntlVec3Gfx(px+gsx,py    ,pz+gsz);
 244                                         pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz);
 245                                         pos[7] = ntlVec3Gfx(px    ,py+gsy,pz+gsz);
 246
 247                                         // check all edges
 248                                         for(int e=0;e<12;e++) {
 249                                                 if (mcEdgeTable[cubeIndex] & (1<<e)) {
 250                                                         // is the vertex already calculated?
 251                                                         if(*eVert[ e ] < 0) {
 252                                                                 // interpolate edge
 253                                                                 const int e1 = mcEdges[e*2  ];
 254                                                                 const int e2 = mcEdges[e*2+1];
 255                                                                 const ntlVec3Gfx p1 = pos[ e1  ];    // scalar field pos 1
 256                                                                 const ntlVec3Gfx p2 = pos[ e2  ];    // scalar field pos 2
 257                                                                 const float valp1  = value[ e1  ];  // scalar field val 1
 258                                                                 const float valp2  = value[ e2  ];  // scalar field val 2
 259                                                                 const float 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                                         if( (i<coAdd+mCutoff) || (j<coAdd+mCutoff) ||
 278                                                         ((mCutoff>0) && (k<coAdd)) ||// bottom layer
 279                                                         (i>mSizex-2-coAdd-mCutoff) ||
 280                                                         (j>mSizey-2-coAdd-mCutoff) ) {
 281                                                 if(mCutArray) {
 282                                                         if(k < mCutArray[j*this->mSizex+i]) continue;
 283                                                 } else { continue; }
 284                                         }
 285
 286                                         // Create the triangles...
 287                                         for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
 288                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
 289                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] );
 290                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
 291                                         }
 292
 293                                 }//i
 294                         }// j
 295
 296                         // copy edge arrays
 297                         if(!mUseFullEdgeArrays) {
 298                         for(int j=0;j<(mSizey-0);j++)
 299                                 for(int i=0;i<(mSizex-0);i++) {
 300                                         //int edgek = 0;
 301                                         const int dst = ISOLEVEL_INDEX( i+0, j+0, 0);
 302                                         const int src = ISOLEVEL_INDEX( i+0, j+0, 1);
 303                                         mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ];
 304                                         mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ];
 305                                         mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ];
 306                                         mpEdgeVerticesX[ src ]=-1;
 307                                         mpEdgeVerticesY[ src ]=-1;
 308                                         mpEdgeVerticesZ[ src ]=-1;
 309                                 }
 310                         } // */
 311
 312                 } // k
 313
 314         // precalculate normals using an approximation of the scalar field gradient
 315                 for(int ni=0;ni<(int)mPoints.size();ni++) { normalize( mPoints[ni].n ); }
 316
 317         } else { // subdivs
 318
 319 #define EDGEAR_INDEX(Ai,Aj,Ak, Bi,Bj) ((mSizex*mSizey*mSubdivs*mSubdivs*(Ak))+\
 320                 (mSizex*mSubdivs*((Aj)*mSubdivs+(Bj)))+((Ai)*mSubdivs)+(Bi))
 321
 322 #define ISOTRILININT(fi,fj,fk) ( \
 323                                 (1.-(fi))*(1.-(fj))*(1.-(fk))*orgval[0] + \
 324                                 (   (fi))*(1.-(fj))*(1.-(fk))*orgval[1] + \
 325                                 (   (fi))*(   (fj))*(1.-(fk))*orgval[2] + \
 326                                 (1.-(fi))*(   (fj))*(1.-(fk))*orgval[3] + \
 327                                 (1.-(fi))*(1.-(fj))*(   (fk))*orgval[4] + \
 328                                 (   (fi))*(1.-(fj))*(   (fk))*orgval[5] + \
 329                                 (   (fi))*(   (fj))*(   (fk))*orgval[6] + \
 330                                 (1.-(fi))*(   (fj))*(   (fk))*orgval[7] )
 331
 332                 // use subdivisions
 333                 gfxReal subdfac = 1./(gfxReal)(mSubdivs);
 334                 gfxReal orgGsx = gsx;
 335                 gfxReal orgGsy = gsy;
 336                 gfxReal orgGsz = gsz;
 337                 gsx *= subdfac;
 338                 gsy *= subdfac;
 339                 gsz *= subdfac;
 340                 if(mUseFullEdgeArrays) {
 341                         errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!");
 342                 }
 343
 344                 // subdiv local arrays
 345                 gfxReal orgval[8];
 346                 gfxReal subdAr[2][11][11]; // max 10 subdivs!
 347                 ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex];
 348
 349                 // construct pointers
 350                 // part test
 351                 int pInUse = 0;
 352                 int pUsedTest = 0;
 353                 // reset particles
 354                 // reset list array
 355                 for(int k=0;k<(mSizez);k++)
 356                         for(int j=0;j<(mSizey);j++)
 357                                 for(int i=0;i<(mSizex);i++) {
 358                                         arppnt[ISOLEVEL_INDEX(i,j,k)] = NULL;
 359                                 }
 360                 if(mpIsoParts) {
 361                         for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
 362                                         pit!= mpIsoParts->getParticlesEnd(); pit++) {
 363                                 if( (*pit).getActive()==false ) continue;
 364                                 if( (*pit).getType()!=PART_DROP) continue;
 365                                 (*pit).setNext(NULL);
 366                         }
 367                         // build per node lists
 368                         for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
 369                                         pit!= mpIsoParts->getParticlesEnd(); pit++) {
 370                                 if( (*pit).getActive()==false ) continue;
 371                                 if( (*pit).getType()!=PART_DROP) continue;
 372                                 // check lifetime ignored here
 373                                 ParticleObject *p = &(*pit);
 374                                 const ntlVec3Gfx ppos = p->getPos();
 375                                 const int pi= (int)round(ppos[0])+0;
 376                                 const int pj= (int)round(ppos[1])+0;
 377                                 int       pk= (int)round(ppos[2])+0;// no offset necessary
 378                                 // 2d should be handled by solver. if(LBMDIM==2) { pk = 0; }
 379
 380                                 if(pi<0) continue;
 381                                 if(pj<0) continue;
 382                                 if(pk<0) continue;
 383                                 if(pi>mSizex-1) continue;
 384                                 if(pj>mSizey-1) continue;
 385                                 if(pk>mSizez-1) continue;
 386                                 ParticleObject* &pnt = arppnt[ISOLEVEL_INDEX(pi,pj,pk)];
 387                                 if(pnt) {
 388                                         // append
 389                                         ParticleObject* listpnt = pnt;
 390                                         while(listpnt) {
 391                                                 if(!listpnt->getNext()) {
 392                                                         listpnt->setNext(p); listpnt = NULL;
 393                                                 } else {
 394                                                         listpnt = listpnt->getNext();
 395                                                 }
 396                                         }
 397                                 } else {
 398                                         // start new list
 399                                         pnt = p;
 400                                 }
 401                                 pInUse++;
 402                         }
 403                 } // mpIsoParts
 404
 405                 debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<<pInUse<<", Subdivs:"<<mSubdivs, 9);
 406                 pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz;
 407                 for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) {
 408                         pz += gsz;
 409                         const int k = ok/mSubdivs;
 410                         if(k<=0) continue; // skip zero plane
 411                         for(int j=1;j<(mSizey-2);j++) {
 412                                 for(int i=1;i<(mSizex-2);i++) {
 413
 414                                         orgval[0] = *getData(i  ,j  ,k  );
 415                                         orgval[1] = *getData(i+1,j  ,k  );
 416                                         orgval[2] = *getData(i+1,j+1,k  ); // with subdivs
 417                                         orgval[3] = *getData(i  ,j+1,k  );
 418                                         orgval[4] = *getData(i  ,j  ,k+1);
 419                                         orgval[5] = *getData(i+1,j  ,k+1);
 420                                         orgval[6] = *getData(i+1,j+1,k+1); // with subdivs
 421                                         orgval[7] = *getData(i  ,j+1,k+1);
 422
 423                                         // prebuild subsampled array slice
 424                                         const int sdkOffset = ok-k*mSubdivs;
 425                                         for(int sdk=0; sdk<2; sdk++)
 426                                                 for(int sdj=0; sdj<mSubdivs+1; sdj++)
 427                                                         for(int sdi=0; sdi<mSubdivs+1; sdi++) {
 428                                                                 subdAr[sdk][sdj][sdi] = ISOTRILININT(sdi*subdfac, sdj*subdfac, (sdkOffset+sdk)*subdfac);
 429                                                         }
 430
 431                                         const int poDistOffset=2;
 432                                         for(int pok=-poDistOffset; pok<1+poDistOffset; pok++) {
 433                                                 if(k+pok<0) continue;
 434                                                 if(k+pok>=mSizez-1) continue;
 435                                         for(int poj=-poDistOffset; poj<1+poDistOffset; poj++) {
 436                                                 if(j+poj<0) continue;
 437                                                 if(j+poj>=mSizey-1) continue;
 438                                         for(int poi=-poDistOffset; poi<1+poDistOffset; poi++) {
 439                                                 if(i+poi<0) continue;
 440                                                 if(i+poi>=mSizex-1) continue;
 441                                                 ParticleObject *p;
 442                                                 p = arppnt[ISOLEVEL_INDEX(i+poi,j+poj,k+pok)];
 443                                                 while(p) { // */
 444                                         /*
 445                                         for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
 446                                                         pit!= mpIsoParts->getParticlesEnd(); pit++) { { { {
 447                                                 // debug test! , full list slow!
 448                                                 if(( (*pit).getActive()==false ) || ( (*pit).getType()!=PART_DROP)) continue;
 449                                                 ParticleObject *p;
 450                                                 p = &(*pit); // */
 451
 452                                                         pUsedTest++;
 453                                                         ntlVec3Gfx ppos = p->getPos();
 454                                                         const int spi= (int)round( (ppos[0]+1.-(gfxReal)i) *(gfxReal)mSubdivs-1.5);
 455                                                         const int spj= (int)round( (ppos[1]+1.-(gfxReal)j) *(gfxReal)mSubdivs-1.5);
 456                                                         const int spk= (int)round( (ppos[2]+1.-(gfxReal)k) *(gfxReal)mSubdivs-1.5)-sdkOffset; // why -2?
 457                                                         // 2d should be handled by solver. if(LBMDIM==2) { spk = 0; }
 458
 459                                                         gfxReal pfLen = p->getSize()*1.5*mPartSize;  // test, was 1.1
 460                                                         const gfxReal minPfLen = subdfac*0.8;
 461                                                         if(pfLen<minPfLen) pfLen = minPfLen;
 462                                                         //errMsg("ISOPPP"," at "<<PRINT_IJK<<"  pp"<<ppos<<"  sp"<<PRINT_VEC(spi,spj,spk)<<" pflen"<<pfLen );
 463                                                         //errMsg("ISOPPP"," subdfac="<<subdfac<<" size"<<p->getSize()<<" ps"<<mPartSize );
 464                                                         const int icellpsize = (int)(1.*pfLen*(gfxReal)mSubdivs)+1;
 465                                                         for(int swk=-icellpsize; swk<=icellpsize; swk++) {
 466                                                                 if(spk+swk<         0) { continue; }
 467                                                                 if(spk+swk>         1) { continue; } // */
 468                                                         for(int swj=-icellpsize; swj<=icellpsize; swj++) {
 469                                                                 if(spj+swj<         0) { continue; }
 470                                                                 if(spj+swj>mSubdivs+0) { continue; } // */
 471                                                         for(int swi=-icellpsize; swi<=icellpsize; swi++) {
 472                                                                 if(spi+swi<         0) { continue; }
 473                                                                 if(spi+swi>mSubdivs+0) { continue; } // */
 474                                                                 ntlVec3Gfx cellp = ntlVec3Gfx(
 475                                                                                 (1.5+(gfxReal)(spi+swi))           *subdfac + (gfxReal)(i-1),
 476                                                                                 (1.5+(gfxReal)(spj+swj))           *subdfac + (gfxReal)(j-1),
 477                                                                                 (1.5+(gfxReal)(spk+swk)+sdkOffset) *subdfac + (gfxReal)(k-1)
 478                                                                                 );
 479                                                                 //if(swi==0 && swj==0 && swk==0) subdAr[spk][spj][spi] = 1.; // DEBUG
 480                                                                 // clip domain boundaries again
 481                                                                 if(cellp[0]<1.) { continue; }
 482                                                                 if(cellp[1]<1.) { continue; }
 483                                                                 if(cellp[2]<1.) { continue; }
 484                                                                 if(cellp[0]>(gfxReal)mSizex-3.) { continue; }
 485                                                                 if(cellp[1]>(gfxReal)mSizey-3.) { continue; }
 486                                                                 if(cellp[2]>(gfxReal)mSizez-3.) { continue; }
 487                                                                 gfxReal len = norm(cellp-ppos);
 488                                                                 gfxReal isoadd = 0.;
 489                                                                 const gfxReal baseIsoVal = mIsoValue*1.1;
 490                                                                 if(len<pfLen) {
 491                                                                         isoadd = baseIsoVal*1.;
 492                                                                 } else {
 493                                                                         // falloff linear with pfLen (kernel size=2pfLen
 494                                                                         isoadd = baseIsoVal*(1. - (len-pfLen)/(pfLen));
 495                                                                 }
 496                                                                 if(isoadd<0.) { continue; }
 497                                                                 //errMsg("ISOPPP"," at "<<PRINT_IJK<<" sp"<<PRINT_VEC(spi+swi,spj+swj,spk+swk)<<" cellp"<<cellp<<" pp"<<ppos << " l"<< len<< " add"<< isoadd);
 498                                                                 const gfxReal arval = subdAr[spk+swk][spj+swj][spi+swi];
 499                                                                 if(arval>1.) { continue; }
 500                                                                 subdAr[spk+swk][spj+swj][spi+swi] = arval + isoadd;
 501                                                         } } }
 502
 503                                                         p = p->getNext();
 504                                                 }
 505                                         } } } // poDist loops */
 506
 507                                         py = mStart[1]+(((double)j-0.5)*orgGsy)-gsy;
 508                                         for(int sj=0;sj<mSubdivs;sj++) {
 509                                                 py += gsy;
 510                                                 px = mStart[0]+(((double)i-0.5)*orgGsx)-gsx;
 511                                                 for(int si=0;si<mSubdivs;si++) {
 512                                                         px += gsx;
 513                                                         value[0] = subdAr[0+0][sj+0][si+0];
 514                                                         value[1] = subdAr[0+0][sj+0][si+1];
 515                                                         value[2] = subdAr[0+0][sj+1][si+1];
 516                                                         value[3] = subdAr[0+0][sj+1][si+0];
 517                                                         value[4] = subdAr[0+1][sj+0][si+0];
 518                                                         value[5] = subdAr[0+1][sj+0][si+1];
 519                                                         value[6] = subdAr[0+1][sj+1][si+1];
 520                                                         value[7] = subdAr[0+1][sj+1][si+0];
 521
 522                                                         // check intersections of isosurface with edges, and calculate cubie index
 523                                                         cubeIndex = 0;
 524                                                         if (value[0] < mIsoValue) cubeIndex |= 1;
 525                                                         if (value[1] < mIsoValue) cubeIndex |= 2; // with subdivs
 526                                                         if (value[2] < mIsoValue) cubeIndex |= 4;
 527                                                         if (value[3] < mIsoValue) cubeIndex |= 8;
 528                                                         if (value[4] < mIsoValue) cubeIndex |= 16;
 529                                                         if (value[5] < mIsoValue) cubeIndex |= 32; // with subdivs
 530                                                         if (value[6] < mIsoValue) cubeIndex |= 64;
 531                                                         if (value[7] < mIsoValue) cubeIndex |= 128;
 532
 533                                                         if (mcEdgeTable[cubeIndex] >  0) {
 534
 535                                                         // where to look up if this point already exists
 536                                                         const int edgek = 0;
 537                                                         const int baseIn = EDGEAR_INDEX( i+0, j+0, edgek+0, si,sj);
 538                                                         eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
 539                                                         eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
 540                                                         eVert[ 2] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ];
 541                                                         eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
 542
 543                                                         eVert[ 4] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ];
 544                                                         eVert[ 5] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+1,sj+0) ]; // with subdivs
 545                                                         eVert[ 6] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+1) ];
 546                                                         eVert[ 7] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ];
 547
 548                                                         eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
 549                                                         eVert[ 9] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+0) ]; // with subdivs
 550                                                         eVert[10] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+1) ];
 551                                                         eVert[11] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ];
 552
 553                                                         // grid positions
 554                                                         pos[0] = ntlVec3Gfx(px    ,py    ,pz);
 555                                                         pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
 556                                                         pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz); // with subdivs
 557                                                         pos[3] = ntlVec3Gfx(px    ,py+gsy,pz);
 558                                                         pos[4] = ntlVec3Gfx(px    ,py    ,pz+gsz);
 559                                                         pos[5] = ntlVec3Gfx(px+gsx,py    ,pz+gsz);
 560                                                         pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz); // with subdivs
 561                                                         pos[7] = ntlVec3Gfx(px    ,py+gsy,pz+gsz);
 562
 563                                                         // check all edges
 564                                                         for(int e=0;e<12;e++) {
 565                                                                 if (mcEdgeTable[cubeIndex] & (1<<e)) {
 566                                                                         // is the vertex already calculated?
 567                                                                         if(*eVert[ e ] < 0) {
 568                                                                                 // interpolate edge
 569                                                                                 const int e1 = mcEdges[e*2  ];
 570                                                                                 const int e2 = mcEdges[e*2+1];
 571                                                                                 const ntlVec3Gfx p1 = pos[ e1  ];   // scalar field pos 1
 572                                                                                 const ntlVec3Gfx p2 = pos[ e2  ];   // scalar field pos 2
 573                                                                                 const float valp1  = value[ e1  ];  // scalar field val 1
 574                                                                                 const float valp2  = value[ e2  ];  // scalar field val 2
 575                                                                                 const float mu = (mIsoValue - valp1) / (valp2 - valp1);
 576
 577                                                                                 // init isolevel vertex
 578                                                                                 ilv.v = p1 + (p2-p1)*mu; // with subdivs
 579                                                                                 mPoints.push_back( ilv );
 580                                                                                 triIndices[e] = (mPoints.size()-1);
 581                                                                                 // store vertex
 582                                                                                 *eVert[ e ] = triIndices[e];
 583                                                                         }       else {
 584                                                                                 // retrieve  from vert array
 585                                                                                 triIndices[e] = *eVert[ e ];
 586                                                                         }
 587                                                                 } // along all edges
 588                                                         }
 589                                                         // removed cutoff treatment...
 590
 591                                                         // Create the triangles...
 592                                                         for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
 593                                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
 594                                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); // with subdivs
 595                                                                 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
 596                                                                 //errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() );
 597                                                         }
 598
 599                                                         } // triangles in edge table?
 600
 601                                                 }//si
 602                                         }// sj
 603
 604                                 }//i
 605                         }// j
 606
 607                         // copy edge arrays
 608                         for(int j=0;j<(mSizey-0)*mSubdivs;j++)
 609                                 for(int i=0;i<(mSizex-0)*mSubdivs;i++) {
 610                                         //int edgek = 0;
 611                                         const int dst = EDGEAR_INDEX( 0, 0, 0, i,j);
 612                                         const int src = EDGEAR_INDEX( 0, 0, 1, i,j);
 613                                         mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ];
 614                                         mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; // with subdivs
 615                                         mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ];
 616                                         mpEdgeVerticesX[ src ]=-1;
 617                                         mpEdgeVerticesY[ src ]=-1; // with subdivs
 618                                         mpEdgeVerticesZ[ src ]=-1;
 619                                 }
 620                         // */
 621
 622                 } // ok, k subdiv loop
 623
 624                 //delete [] subdAr;
 625                 delete [] arppnt;
 626                 computeNormals();
 627         } // with subdivs
 628
 629         // perform smoothing
 630         float smoSubdfac = 1.;
 631         if(mSubdivs>0) {
 632                 //smoSubdfac = 1./(float)(mSubdivs);
 633                 smoSubdfac = pow(0.55,(double)mSubdivs); // slightly stronger
 634         }
 635         if(mSmoothSurface>0. || mSmoothNormals>0.) debMsgStd("IsoSurface::triangulate",DM_MSG,"Smoothing...",10);
 636         if(mSmoothSurface>0.0) {
 637                 smoothSurface(mSmoothSurface*smoSubdfac, (mSmoothNormals<=0.0) );
 638         }
 639         if(mSmoothNormals>0.0) {
 640                 smoothNormals(mSmoothNormals*smoSubdfac);
 641         }
 642
 643         myTime_t tritimeend = getTime();
 644         debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< getTimeString(tritimeend-tritimestart)<<", S("<<mSmoothSurface<<","<<mSmoothNormals<<"),"<<
 645                         " verts:"<<mPoints.size()<<" tris:"<<(mIndices.size()/3)<<" subdivs:"<<mSubdivs
 646                  , 10 );
 647         if(mpIsoParts) debMsgStd("IsoSurface::triangulate",DM_MSG,"parts:"<<mpIsoParts->getNumParticles(), 10);
 648 }
 649
 650
 651
 652
 653
 654 /******************************************************************************
 655  * Get triangles for rendering
 656  *****************************************************************************/
 657 void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles,
 658                                                                                                          vector<ntlVec3Gfx> *vertices,
 659                                                                                                          vector<ntlVec3Gfx> *normals, int objectId )
 660 {
 661         if(!mInitDone) {
 662                 debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10);
 663                 return;
 664         }
 665         t = 0.;
 666         //return; // DEBUG
 667
 668   /* triangulate field */
 669   triangulate();
 670         //errMsg("TRIS"," "<<mIndices.size() );
 671
 672         // new output with vertice reuse
 673         int iniVertIndex = (*vertices).size();
 674         int iniNormIndex = (*normals).size();
 675         if(iniVertIndex != iniNormIndex) {
 676                 errFatal("getTriangles Error","For '"<<mName<<"': Vertices and normal array sizes to not match!!!",SIMWORLD_GENERICERROR);
 677                 return;
 678         }
 679         //errMsg("NM"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
 680         //errMsg("NM"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
 681
 682   for(int i=0;i<(int)mPoints.size();i++) {
 683                 vertices->push_back( mPoints[i].v );
 684         }
 685   for(int i=0;i<(int)mPoints.size();i++) {
 686                 normals->push_back( mPoints[i].n );
 687         }
 688
 689         //errMsg("N2"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
 690         //errMsg("N2"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
 691
 692   for(int i=0;i<(int)mIndices.size();i+=3) {
 693                 const int smooth = 1;
 694     int t1 = mIndices[i];
 695     int t2 = mIndices[i+1];
 696                 int t3 = mIndices[i+2];
 697                 //errMsg("NM"," tri"<<t1<<" "<<t2<<" "<<t3 );
 698
 699                 ntlTriangle tri;
 700
 701                 tri.getPoints()[0] = t1+iniVertIndex;
 702                 tri.getPoints()[1] = t2+iniVertIndex;
 703                 tri.getPoints()[2] = t3+iniVertIndex;
 704
 705                 /* init flags */
 706                 int flag = 0;
 707                 if(getVisible()){ flag |= TRI_GEOMETRY; }
 708                 if(getCastShadows() ) {
 709                         flag |= TRI_CASTSHADOWS; }
 710
 711                 /* init geo init id */
 712                 int geoiId = getGeoInitId();
 713                 if(geoiId > 0) {
 714                         flag |= (1<< (geoiId+4));
 715                         flag |= mGeoInitType;
 716                 }
 717
 718                 tri.setFlags( flag );
 719
 720                 /* triangle normal missing */
 721                 tri.setNormal( ntlVec3Gfx(0.0) );
 722                 tri.setSmoothNormals( smooth );
 723                 tri.setObjectId( objectId );
 724                 triangles->push_back( tri );
 725         }
 726         //errMsg("N3"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
 727         return;
 728 }
 729
 730
 731
 732 inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) {
 733         // WARNING - this requires a security boundary layer...
 734         ntlVec3Gfx ret(0.0);
 735         ret[0] = *getData(i-1,j  ,k  ) -
 736                  *getData(i+1,j  ,k  );
 737         ret[1] = *getData(i  ,j-1,k  ) -
 738                  *getData(i  ,j+1,k  );
 739         ret[2] = *getData(i  ,j  ,k-1  ) -
 740                  *getData(i  ,j  ,k+1  );
 741         return ret;
 742 }
 743
 744
 745
 746
 747 /******************************************************************************
 748  *
 749  * Surface improvement, inspired by trimesh2 library
 750  * (http://www.cs.princeton.edu/gfx/proj/trimesh2/)
 751  *
 752  *****************************************************************************/
 753
 754 void IsoSurface::setSmoothRad(float radi1, float radi2, ntlVec3Gfx mscc) {
 755         mSCrad1 = radi1*radi1;
 756         mSCrad2 = radi2*radi2;
 757         mSCcenter = mscc;
 758 }
 759
 760 // compute normals for all generated triangles
 761 void IsoSurface::computeNormals() {
 762   for(int i=0;i<(int)mPoints.size();i++) {
 763                 mPoints[i].n = ntlVec3Gfx(0.);
 764         }
 765
 766   for(int i=0;i<(int)mIndices.size();i+=3) {
 767     const int t1 = mIndices[i];
 768     const int t2 = mIndices[i+1];
 769                 const int t3 = mIndices[i+2];
 770                 const ntlVec3Gfx p1 = mPoints[t1].v;
 771                 const ntlVec3Gfx p2 = mPoints[t2].v;
 772                 const ntlVec3Gfx p3 = mPoints[t3].v;
 773                 const ntlVec3Gfx n1=p1-p2;
 774                 const ntlVec3Gfx n2=p2-p3;
 775                 const ntlVec3Gfx n3=p3-p1;
 776                 const gfxReal len1 = normNoSqrt(n1);
 777                 const gfxReal len2 = normNoSqrt(n2);
 778                 const gfxReal len3 = normNoSqrt(n3);
 779                 const ntlVec3Gfx norm = cross(n1,n2);
 780                 mPoints[t1].n += norm * (1./(len1*len3));
 781                 mPoints[t2].n += norm * (1./(len1*len2));
 782                 mPoints[t3].n += norm * (1./(len2*len3));
 783         }
 784
 785   for(int i=0;i<(int)mPoints.size();i++) {
 786                 normalize(mPoints[i].n);
 787         }
 788 }
 789
 790 // Diffuse a vector field at 1 vertex, weighted by
 791 // a gaussian of width 1/sqrt(invsigma2)
 792 bool IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int src, float invsigma2, ntlVec3Gfx &target)
 793 {
 794         if((neighbors[src].size()<1) || (pointareas[src]<=0.0)) return 0;
 795         const ntlVec3Gfx srcp = mPoints[src].v;
 796         const ntlVec3Gfx srcn = mPoints[src].n;
 797         if(mSCrad1>0.0 && mSCrad2>0.0) {
 798                 ntlVec3Gfx dp = mSCcenter-srcp; dp[2] = 0.0; // only xy-plane
 799                 float rd = normNoSqrt(dp);
 800                 if(rd > mSCrad2) {
 801                         return 0;
 802                 } else if(rd > mSCrad1) {
 803                         // optimize?
 804                         float org = 1.0/sqrt(invsigma2);
 805                         org *= (1.0- (rd-mSCrad1) / (mSCrad2-mSCrad1));
 806                         invsigma2 = 1.0/(org*org);
 807                         //errMsg("TRi","p"<<srcp<<" rd:"<<rd<<" r1:"<<mSCrad1<<" r2:"<<mSCrad2<<" org:"<<org<<" is:"<<invsigma2);
 808                 } else {
 809                 }
 810         }
 811         target = ntlVec3Gfx(0.0);
 812         target += *(field+pointerScale*src) *pointareas[src];
 813         float smstrSum = pointareas[src];
 814
 815         int flag = mFlagCnt;
 816         mFlagCnt++;
 817         flags[src] = flag;
 818         mDboundary = neighbors[src];
 819         while (!mDboundary.empty()) {
 820                 const int bbn = mDboundary.back();
 821                 mDboundary.pop_back();
 822                 if(flags[bbn]==flag) continue;
 823                 flags[bbn] = flag;
 824
 825                 // normal check
 826                 const float nvdot = dot(srcn, mPoints[bbn].n); // faster than before d2 calc?
 827                 if(nvdot <= 0.0f) continue;
 828
 829                 // gaussian weight of width 1/sqrt(invsigma2)
 830                 const float d2 = invsigma2 * normNoSqrt(mPoints[bbn].v - srcp);
 831                 if(d2 >= 9.0f) continue;
 832
 833                 // aggressive smoothing factor
 834                 float smstr = nvdot * pointareas[bbn];
 835                 // Accumulate weight times field at neighbor
 836                 target += *(field+pointerScale*bbn)*smstr;
 837                 smstrSum += smstr;
 838
 839                 for(int i = 0; i < (int)neighbors[bbn].size(); i++) {
 840                         const int nn = neighbors[bbn][i];
 841                         if (flags[nn] == flag) continue;
 842                         mDboundary.push_back(nn);
 843                 }
 844         }
 845         target /= smstrSum;
 846         return 1;
 847 }
 848
 849
 850 // perform smoothing of the surface (and possible normals)
 851 void IsoSurface::smoothSurface(float sigma, bool normSmooth)
 852 {
 853         int nv = mPoints.size();
 854         if ((int)flags.size() != nv) flags.resize(nv);
 855         int nf = mIndices.size()/3;
 856
 857         { // need neighbors
 858                 vector<int> numneighbors(mPoints.size());
 859                 int i;
 860                 for (i = 0; i < (int)mIndices.size()/3; i++) {
 861                         numneighbors[mIndices[i*3+0]]++;
 862                         numneighbors[mIndices[i*3+1]]++;
 863                         numneighbors[mIndices[i*3+2]]++;
 864                 }
 865
 866                 neighbors.clear();
 867                 neighbors.resize(mPoints.size());
 868                 for (i = 0; i < (int)mPoints.size(); i++) {
 869                         neighbors[i].clear();
 870                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
 871                 }
 872
 873                 for (i = 0; i < (int)mIndices.size()/3; i++) {
 874                         for (int j = 0; j < 3; j++) {
 875                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
 876                                 int n1 =  mIndices[i*3+((j+1)%3)];
 877                                 int n2 =  mIndices[i*3+((j+2)%3)];
 878                                 if (std::find(me.begin(), me.end(), n1) == me.end())
 879                                         me.push_back(n1);
 880                                 if (std::find(me.begin(), me.end(), n2) == me.end())
 881                                         me.push_back(n2);
 882                         }
 883                 }
 884         } // need neighbor
 885
 886         { // need pointarea
 887                 pointareas.clear();
 888                 pointareas.resize(nv);
 889                 cornerareas.clear();
 890                 cornerareas.resize(nf);
 891
 892                 for (int i = 0; i < nf; i++) {
 893                         // Edges
 894                         ntlVec3Gfx e[3] = {
 895                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
 896                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
 897                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
 898
 899                         // Compute corner weights
 900                         float area = 0.5f * norm( cross(e[0], e[1]));
 901                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
 902                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
 903                                         l2[1] * (l2[2] + l2[0] - l2[1]),
 904                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
 905                         if (ew[0] <= 0.0f) {
 906                                 cornerareas[i][1] = -0.25f * l2[2] * area /
 907                                                                 dot(e[0] , e[2]);
 908                                 cornerareas[i][2] = -0.25f * l2[1] * area /
 909                                                                 dot(e[0] , e[1]);
 910                                 cornerareas[i][0] = area - cornerareas[i][1] -
 911                                                                 cornerareas[i][2];
 912                         } else if (ew[1] <= 0.0f) {
 913                                 cornerareas[i][2] = -0.25f * l2[0] * area /
 914                                                                 dot(e[1] , e[0]);
 915                                 cornerareas[i][0] = -0.25f * l2[2] * area /
 916                                                                 dot(e[1] , e[2]);
 917                                 cornerareas[i][1] = area - cornerareas[i][2] -
 918                                                                 cornerareas[i][0];
 919                         } else if (ew[2] <= 0.0f) {
 920                                 cornerareas[i][0] = -0.25f * l2[1] * area /
 921                                                                 dot(e[2] , e[1]);
 922                                 cornerareas[i][1] = -0.25f * l2[0] * area /
 923                                                                 dot(e[2] , e[0]);
 924                                 cornerareas[i][2] = area - cornerareas[i][0] -
 925                                                                 cornerareas[i][1];
 926                         } else {
 927                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
 928                                 for (int j = 0; j < 3; j++)
 929                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
 930                                                                                          ew[(j+2)%3]);
 931                         }
 932
 933                         // NT important, check this...
 934 #ifndef WIN32
 935                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
 936                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
 937                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
 938 #else // WIN32
 939                         // FIXME check as well...
 940                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
 941                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
 942                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
 943 #endif // WIN32
 944
 945                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
 946                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
 947                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
 948                 }
 949
 950         } // need pointarea
 951         // */
 952
 953         float invsigma2 = 1.0f / (sigma*sigma);
 954
 955         vector<ntlVec3Gfx> dflt(nv);
 956         for (int i = 0; i < nv; i++) {
 957                 if(diffuseVertexField( &mPoints[0].v, 2,
 958                                    i, invsigma2, dflt[i])) {
 959                         // Just keep the displacement
 960                         dflt[i] -= mPoints[i].v;
 961                 } else { dflt[i] = 0.0; } //?mPoints[i].v; }
 962         }
 963
 964         // Slightly better small-neighborhood approximation
 965         for (int i = 0; i < nf; i++) {
 966                 ntlVec3Gfx c = mPoints[mIndices[i*3+0]].v +
 967                           mPoints[mIndices[i*3+1]].v +
 968                           mPoints[mIndices[i*3+2]].v;
 969                 c /= 3.0f;
 970                 for (int j = 0; j < 3; j++) {
 971                         int v = mIndices[i*3+j];
 972                         ntlVec3Gfx d =(c - mPoints[v].v) * 0.5f;
 973                         dflt[v] += d * (cornerareas[i][j] /
 974                                    pointareas[mIndices[i*3+j]] *
 975                                    exp(-0.5f * invsigma2 * normNoSqrt(d)) );
 976                 }
 977         }
 978
 979         // Filter displacement field
 980         vector<ntlVec3Gfx> dflt2(nv);
 981         for (int i = 0; i < nv; i++) {
 982                 if(diffuseVertexField( &dflt[0], 1,
 983                                    i, invsigma2, dflt2[i])) { }
 984                 else { /*mPoints[i].v=0.0;*/ dflt2[i] = 0.0; }//dflt2[i]; }
 985         }
 986
 987         // Update vertex positions
 988         for (int i = 0; i < nv; i++) {
 989                 mPoints[i].v += dflt[i] - dflt2[i]; // second Laplacian
 990         }
 991
 992         // when normals smoothing off, this cleans up quite well
 993         // costs ca. 50% additional time though
 994         float nsFac = 1.5f;
 995         if(normSmooth) { float ninvsigma2 = 1.0f / (nsFac*nsFac*sigma*sigma);
 996                 for (int i = 0; i < nv; i++) {
 997                         if( diffuseVertexField( &mPoints[0].n, 2, i, ninvsigma2, dflt[i]) ) {
 998                                 normalize(dflt[i]);
 999                         } else {
1000                                 dflt[i] = mPoints[i].n;
1001                         }
1002                 }
1003                 for (int i = 0; i < nv; i++) {
1004                         mPoints[i].n = dflt[i];
1005                 }
1006         } // smoothNormals copy */
1007
1008         //errMsg("SMSURF","done v:"<<sigma); // DEBUG
1009 }
1010
1011 // only smoothen the normals
1012 void IsoSurface::smoothNormals(float sigma) {
1013         // reuse from smoothSurface
1014         if(neighbors.size() != mPoints.size()) {
1015                 // need neighbor
1016                 vector<int> numneighbors(mPoints.size());
1017                 int i;
1018                 for (i = 0; i < (int)mIndices.size()/3; i++) {
1019                         numneighbors[mIndices[i*3+0]]++;
1020                         numneighbors[mIndices[i*3+1]]++;
1021                         numneighbors[mIndices[i*3+2]]++;
1022                 }
1023
1024                 neighbors.clear();
1025                 neighbors.resize(mPoints.size());
1026                 for (i = 0; i < (int)mPoints.size(); i++) {
1027                         neighbors[i].clear();
1028                         neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries
1029                 }
1030
1031                 for (i = 0; i < (int)mIndices.size()/3; i++) {
1032                         for (int j = 0; j < 3; j++) {
1033                                 vector<int> &me = neighbors[ mIndices[i*3+j]];
1034                                 int n1 =  mIndices[i*3+((j+1)%3)];
1035                                 int n2 =  mIndices[i*3+((j+2)%3)];
1036                                 if (std::find(me.begin(), me.end(), n1) == me.end())
1037                                         me.push_back(n1);
1038                                 if (std::find(me.begin(), me.end(), n2) == me.end())
1039                                         me.push_back(n2);
1040                         }
1041                 }
1042         } // need neighbor
1043
1044         { // need pointarea
1045                 int nf = mIndices.size()/3, nv = mPoints.size();
1046                 pointareas.clear();
1047                 pointareas.resize(nv);
1048                 cornerareas.clear();
1049                 cornerareas.resize(nf);
1050
1051                 for (int i = 0; i < nf; i++) {
1052                         // Edges
1053                         ntlVec3Gfx e[3] = {
1054                                 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v,
1055                                 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v,
1056                                 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v };
1057
1058                         // Compute corner weights
1059                         float area = 0.5f * norm( cross(e[0], e[1]));
1060                         float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) };
1061                         float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]),
1062                                         l2[1] * (l2[2] + l2[0] - l2[1]),
1063                                         l2[2] * (l2[0] + l2[1] - l2[2]) };
1064                         if (ew[0] <= 0.0f) {
1065                                 cornerareas[i][1] = -0.25f * l2[2] * area /
1066                                                                 dot(e[0] , e[2]);
1067                                 cornerareas[i][2] = -0.25f * l2[1] * area /
1068                                                                 dot(e[0] , e[1]);
1069                                 cornerareas[i][0] = area - cornerareas[i][1] -
1070                                                                 cornerareas[i][2];
1071                         } else if (ew[1] <= 0.0f) {
1072                                 cornerareas[i][2] = -0.25f * l2[0] * area /
1073                                                                 dot(e[1] , e[0]);
1074                                 cornerareas[i][0] = -0.25f * l2[2] * area /
1075                                                                 dot(e[1] , e[2]);
1076                                 cornerareas[i][1] = area - cornerareas[i][2] -
1077                                                                 cornerareas[i][0];
1078                         } else if (ew[2] <= 0.0f) {
1079                                 cornerareas[i][0] = -0.25f * l2[1] * area /
1080                                                                 dot(e[2] , e[1]);
1081                                 cornerareas[i][1] = -0.25f * l2[0] * area /
1082                                                                 dot(e[2] , e[0]);
1083                                 cornerareas[i][2] = area - cornerareas[i][0] -
1084                                                                 cornerareas[i][1];
1085                         } else {
1086                                 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
1087                                 for (int j = 0; j < 3; j++)
1088                                         cornerareas[i][j] = ewscale * (ew[(j+1)%3] +
1089                                                                                          ew[(j+2)%3]);
1090                         }
1091
1092                         // NT important, check this...
1093 #ifndef WIN32
1094                         if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6;
1095                         if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6;
1096                         if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6;
1097 #else // WIN32
1098                         // FIXME check as well...
1099                         if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6;
1100                         if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6;
1101                         if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6;
1102 #endif // WIN32
1103
1104                         pointareas[mIndices[i*3+0]] += cornerareas[i][0];
1105                         pointareas[mIndices[i*3+1]] += cornerareas[i][1];
1106                         pointareas[mIndices[i*3+2]] += cornerareas[i][2];
1107                 }
1108
1109         } // need pointarea
1110
1111         int nv = mPoints.size();
1112         if ((int)flags.size() != nv) flags.resize(nv);
1113         float invsigma2 = 1.0f / (sigma*sigma);
1114
1115         vector<ntlVec3Gfx> nflt(nv);
1116         for (int i = 0; i < nv; i++) {
1117                 if(diffuseVertexField( &mPoints[0].n, 2, i, invsigma2, nflt[i])) {
1118                         normalize(nflt[i]);
1119                 } else { nflt[i]=mPoints[i].n; }
1120         }
1121
1122         // copy results
1123         for (int i = 0; i < nv; i++) { mPoints[i].n = nflt[i]; }
1124 }
1125
1126