 |
Blender V2.61 - r43446
|
isosurface.cpp
Go to the documentation of this file.
00001 00004 /****************************************************************************** 00005 * 00006 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method 00007 * Copyright 2003-2006 Nils Thuerey 00008 * 00009 * Marching Cubes surface mesh generation 00010 * 00011 *****************************************************************************/ 00012 00013 #include "isosurface.h" 00014 #include "mcubes_tables.h" 00015 #include "particletracer.h" 00016 #include <algorithm> 00017 #include <stdio.h> 00018 00019 #ifdef sun 00020 #include "ieeefp.h" 00021 #endif 00022 00023 // just use default rounding for platforms where its not available 00024 #ifndef round 00025 #define round(x) (x) 00026 #endif 00027 00028 /****************************************************************************** 00029 * Constructor 00030 *****************************************************************************/ 00031 IsoSurface::IsoSurface(double iso) : 00032 ntlGeometryObject(), 00033 mSizex(-1), mSizey(-1), mSizez(-1), 00034 mpData(NULL), 00035 mIsoValue( iso ), 00036 mPoints(), 00037 mUseFullEdgeArrays(false), 00038 mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL), 00039 mEdgeArSize(-1), 00040 mIndices(), 00041 00042 mStart(0.0), mEnd(0.0), mDomainExtent(0.0), 00043 mInitDone(false), 00044 mSmoothSurface(0.0), mSmoothNormals(0.0), 00045 mAcrossEdge(), mAdjacentFaces(), 00046 mCutoff(-1), mCutArray(NULL), // off by default 00047 mpIsoParts(NULL), mPartSize(0.), mSubdivs(0), 00048 mFlagCnt(1), 00049 mSCrad1(0.), mSCrad2(0.), mSCcenter(0.) 00050 { 00051 } 00052 00053 00054 /****************************************************************************** 00055 * The real init... 00056 *****************************************************************************/ 00057 void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent) 00058 { 00059 // range 1-10 (max due to subd array in triangulate) 00060 if(mSubdivs<1) mSubdivs=1; 00061 if(mSubdivs>10) mSubdivs=10; 00062 00063 // init solver and size 00064 mSizex = setx; 00065 mSizey = sety; 00066 if(setz == 1) {// 2D, create thin 2D surface 00067 setz = 5; 00068 } 00069 mSizez = setz; 00070 mDomainExtent = extent; 00071 00072 /* check triangulation size (for raytraing) */ 00073 if( ( mStart[0] >= mEnd[0] ) && ( mStart[1] >= mEnd[1] ) && ( mStart[2] >= mEnd[2] ) ){ 00074 // extent was not set, use normalized one from parametrizer 00075 mStart = ntlVec3Gfx(0.0) - extent*0.5; 00076 mEnd = ntlVec3Gfx(0.0) + extent*0.5; 00077 } 00078 00079 // init 00080 mIndices.clear(); 00081 mPoints.clear(); 00082 00083 int nodes = mSizez*mSizey*mSizex; 00084 mpData = new float[nodes]; 00085 for(int i=0;i<nodes;i++) { mpData[i] = 0.0; } 00086 00087 // allocate edge arrays (last slices are never used...) 00088 int initsize = -1; 00089 if(mUseFullEdgeArrays) { 00090 mEdgeArSize = nodes; 00091 mpEdgeVerticesX = new int[nodes]; 00092 mpEdgeVerticesY = new int[nodes]; 00093 mpEdgeVerticesZ = new int[nodes]; 00094 initsize = 3*nodes; 00095 } else { 00096 int sliceNodes = 2*mSizex*mSizey*mSubdivs*mSubdivs; 00097 mEdgeArSize = sliceNodes; 00098 mpEdgeVerticesX = new int[sliceNodes]; 00099 mpEdgeVerticesY = new int[sliceNodes]; 00100 mpEdgeVerticesZ = new int[sliceNodes]; 00101 initsize = 3*sliceNodes; 00102 } 00103 for(int i=0;i<mEdgeArSize;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; } 00104 // WARNING - make sure this is consistent with calculateMemreqEstimate 00105 00106 // marching cubes are ready 00107 mInitDone = true; 00108 debMsgStd("IsoSurface::initializeIsosurface",DM_MSG,"Inited, edgenodes:"<<initsize<<" subdivs:"<<mSubdivs<<" fulledg:"<<mUseFullEdgeArrays , 10); 00109 } 00110 00111 00112 00114 void IsoSurface::resetAll(gfxReal val) { 00115 int nodes = mSizez*mSizey*mSizex; 00116 for(int i=0;i<nodes;i++) { mpData[i] = val; } 00117 } 00118 00119 00120 /****************************************************************************** 00121 * Destructor 00122 *****************************************************************************/ 00123 IsoSurface::~IsoSurface( void ) 00124 { 00125 if(mpData) delete [] mpData; 00126 if(mpEdgeVerticesX) delete [] mpEdgeVerticesX; 00127 if(mpEdgeVerticesY) delete [] mpEdgeVerticesY; 00128 if(mpEdgeVerticesZ) delete [] mpEdgeVerticesZ; 00129 } 00130 00131 00132 00133 00134 00135 /****************************************************************************** 00136 * triangulate the scalar field given by pointer 00137 *****************************************************************************/ 00138 void IsoSurface::triangulate( void ) 00139 { 00140 double gsx,gsy,gsz; // grid spacing in x,y,z direction 00141 double px,py,pz; // current position in grid in x,y,z direction 00142 IsoLevelCube cubie; // struct for a small subcube 00143 myTime_t tritimestart = getTime(); 00144 00145 if(!mpData) { 00146 errFatal("IsoSurface::triangulate","no LBM object, and no scalar field...!",SIMWORLD_INITERROR); 00147 return; 00148 } 00149 00150 // get grid spacing (-2 to have same spacing as sim) 00151 gsx = (mEnd[0]-mStart[0])/(double)(mSizex-2.0); 00152 gsy = (mEnd[1]-mStart[1])/(double)(mSizey-2.0); 00153 gsz = (mEnd[2]-mStart[2])/(double)(mSizez-2.0); 00154 00155 // clean up previous frame 00156 mIndices.clear(); 00157 mPoints.clear(); 00158 00159 // reset edge vertices 00160 for(int i=0;i<mEdgeArSize;i++) { 00161 mpEdgeVerticesX[i] = -1; 00162 mpEdgeVerticesY[i] = -1; 00163 mpEdgeVerticesZ[i] = -1; 00164 } 00165 00166 ntlVec3Gfx pos[8]; 00167 float value[8]; 00168 int cubeIndex; // index entry of the cube 00169 int triIndices[12]; // vertex indices 00170 int *eVert[12]; 00171 IsoLevelVertex ilv; 00172 00173 // edges between which points? 00174 const int mcEdges[24] = { 00175 0,1, 1,2, 2,3, 3,0, 00176 4,5, 5,6, 6,7, 7,4, 00177 0,4, 1,5, 2,6, 3,7 }; 00178 00179 const int cubieOffsetX[8] = { 00180 0,1,1,0, 0,1,1,0 }; 00181 const int cubieOffsetY[8] = { 00182 0,0,1,1, 0,0,1,1 }; 00183 const int cubieOffsetZ[8] = { 00184 0,0,0,0, 1,1,1,1 }; 00185 00186 const int coAdd=2; 00187 // let the cubes march 00188 if(mSubdivs<=1) { 00189 00190 pz = mStart[2]-gsz*0.5; 00191 for(int k=1;k<(mSizez-2);k++) { 00192 pz += gsz; 00193 py = mStart[1]-gsy*0.5; 00194 for(int j=1;j<(mSizey-2);j++) { 00195 py += gsy; 00196 px = mStart[0]-gsx*0.5; 00197 for(int i=1;i<(mSizex-2);i++) { 00198 px += gsx; 00199 00200 value[0] = *getData(i ,j ,k ); 00201 value[1] = *getData(i+1,j ,k ); 00202 value[2] = *getData(i+1,j+1,k ); 00203 value[3] = *getData(i ,j+1,k ); 00204 value[4] = *getData(i ,j ,k+1); 00205 value[5] = *getData(i+1,j ,k+1); 00206 value[6] = *getData(i+1,j+1,k+1); 00207 value[7] = *getData(i ,j+1,k+1); 00208 00209 // check intersections of isosurface with edges, and calculate cubie index 00210 cubeIndex = 0; 00211 if (value[0] < mIsoValue) cubeIndex |= 1; 00212 if (value[1] < mIsoValue) cubeIndex |= 2; 00213 if (value[2] < mIsoValue) cubeIndex |= 4; 00214 if (value[3] < mIsoValue) cubeIndex |= 8; 00215 if (value[4] < mIsoValue) cubeIndex |= 16; 00216 if (value[5] < mIsoValue) cubeIndex |= 32; 00217 if (value[6] < mIsoValue) cubeIndex |= 64; 00218 if (value[7] < mIsoValue) cubeIndex |= 128; 00219 00220 // No triangles to generate? 00221 if (mcEdgeTable[cubeIndex] == 0) { 00222 continue; 00223 } 00224 00225 // where to look up if this point already exists 00226 int edgek = 0; 00227 if(mUseFullEdgeArrays) edgek=k; 00228 const int baseIn = ISOLEVEL_INDEX( i+0, j+0, edgek+0); 00229 eVert[ 0] = &mpEdgeVerticesX[ baseIn ]; 00230 eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ]; 00231 eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ]; 00232 eVert[ 3] = &mpEdgeVerticesY[ baseIn ]; 00233 00234 eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ]; 00235 eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, edgek+1) ]; 00236 eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+1) ]; 00237 eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ]; 00238 00239 eVert[ 8] = &mpEdgeVerticesZ[ baseIn ]; 00240 eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, edgek+0) ]; 00241 eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, edgek+0) ]; 00242 eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ]; 00243 00244 // grid positions 00245 pos[0] = ntlVec3Gfx(px ,py ,pz); 00246 pos[1] = ntlVec3Gfx(px+gsx,py ,pz); 00247 pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz); 00248 pos[3] = ntlVec3Gfx(px ,py+gsy,pz); 00249 pos[4] = ntlVec3Gfx(px ,py ,pz+gsz); 00250 pos[5] = ntlVec3Gfx(px+gsx,py ,pz+gsz); 00251 pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz); 00252 pos[7] = ntlVec3Gfx(px ,py+gsy,pz+gsz); 00253 00254 // check all edges 00255 for(int e=0;e<12;e++) { 00256 if (mcEdgeTable[cubeIndex] & (1<<e)) { 00257 // is the vertex already calculated? 00258 if(*eVert[ e ] < 0) { 00259 // interpolate edge 00260 const int e1 = mcEdges[e*2 ]; 00261 const int e2 = mcEdges[e*2+1]; 00262 const ntlVec3Gfx p1 = pos[ e1 ]; // scalar field pos 1 00263 const ntlVec3Gfx p2 = pos[ e2 ]; // scalar field pos 2 00264 const float valp1 = value[ e1 ]; // scalar field val 1 00265 const float valp2 = value[ e2 ]; // scalar field val 2 00266 const float mu = (mIsoValue - valp1) / (valp2 - valp1); 00267 00268 // init isolevel vertex 00269 ilv.v = p1 + (p2-p1)*mu; 00270 ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) + 00271 getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * ( mu) ; 00272 mPoints.push_back( ilv ); 00273 00274 triIndices[e] = (mPoints.size()-1); 00275 // store vertex 00276 *eVert[ e ] = triIndices[e]; 00277 } else { 00278 // retrieve from vert array 00279 triIndices[e] = *eVert[ e ]; 00280 } 00281 } // along all edges 00282 } 00283 00284 if( (i<coAdd+mCutoff) || (j<coAdd+mCutoff) || 00285 ((mCutoff>0) && (k<coAdd)) ||// bottom layer 00286 (i>mSizex-2-coAdd-mCutoff) || 00287 (j>mSizey-2-coAdd-mCutoff) ) { 00288 if(mCutArray) { 00289 if(k < mCutArray[j*this->mSizex+i]) continue; 00290 } else { continue; } 00291 } 00292 00293 // Create the triangles... 00294 for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) { 00295 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] ); 00296 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); 00297 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] ); 00298 } 00299 00300 }//i 00301 }// j 00302 00303 // copy edge arrays 00304 if(!mUseFullEdgeArrays) { 00305 for(int j=0;j<(mSizey-0);j++) 00306 for(int i=0;i<(mSizex-0);i++) { 00307 //int edgek = 0; 00308 const int dst = ISOLEVEL_INDEX( i+0, j+0, 0); 00309 const int src = ISOLEVEL_INDEX( i+0, j+0, 1); 00310 mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ]; 00311 mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; 00312 mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ]; 00313 mpEdgeVerticesX[ src ]=-1; 00314 mpEdgeVerticesY[ src ]=-1; 00315 mpEdgeVerticesZ[ src ]=-1; 00316 } 00317 } // */ 00318 00319 } // k 00320 00321 // precalculate normals using an approximation of the scalar field gradient 00322 for(int ni=0;ni<(int)mPoints.size();ni++) { normalize( mPoints[ni].n ); } 00323 00324 } else { // subdivs 00325 00326 #define EDGEAR_INDEX(Ai,Aj,Ak, Bi,Bj) ((mSizex*mSizey*mSubdivs*mSubdivs*(Ak))+\ 00327 (mSizex*mSubdivs*((Aj)*mSubdivs+(Bj)))+((Ai)*mSubdivs)+(Bi)) 00328 00329 #define ISOTRILININT(fi,fj,fk) ( \ 00330 (1.-(fi))*(1.-(fj))*(1.-(fk))*orgval[0] + \ 00331 ( (fi))*(1.-(fj))*(1.-(fk))*orgval[1] + \ 00332 ( (fi))*( (fj))*(1.-(fk))*orgval[2] + \ 00333 (1.-(fi))*( (fj))*(1.-(fk))*orgval[3] + \ 00334 (1.-(fi))*(1.-(fj))*( (fk))*orgval[4] + \ 00335 ( (fi))*(1.-(fj))*( (fk))*orgval[5] + \ 00336 ( (fi))*( (fj))*( (fk))*orgval[6] + \ 00337 (1.-(fi))*( (fj))*( (fk))*orgval[7] ) 00338 00339 // use subdivisions 00340 gfxReal subdfac = 1./(gfxReal)(mSubdivs); 00341 gfxReal orgGsx = gsx; 00342 gfxReal orgGsy = gsy; 00343 gfxReal orgGsz = gsz; 00344 gsx *= subdfac; 00345 gsy *= subdfac; 00346 gsz *= subdfac; 00347 if(mUseFullEdgeArrays) { 00348 errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!"); 00349 } 00350 00351 // subdiv local arrays 00352 gfxReal orgval[8]; 00353 gfxReal subdAr[2][11][11]; // max 10 subdivs! 00354 ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex]; 00355 00356 // construct pointers 00357 // part test 00358 int pInUse = 0; 00359 int pUsedTest = 0; 00360 // reset particles 00361 // reset list array 00362 for(int k=0;k<(mSizez);k++) 00363 for(int j=0;j<(mSizey);j++) 00364 for(int i=0;i<(mSizex);i++) { 00365 arppnt[ISOLEVEL_INDEX(i,j,k)] = NULL; 00366 } 00367 if(mpIsoParts) { 00368 for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin(); 00369 pit!= mpIsoParts->getParticlesEnd(); pit++) { 00370 if( (*pit).getActive()==false ) continue; 00371 if( (*pit).getType()!=PART_DROP) continue; 00372 (*pit).setNext(NULL); 00373 } 00374 // build per node lists 00375 for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin(); 00376 pit!= mpIsoParts->getParticlesEnd(); pit++) { 00377 if( (*pit).getActive()==false ) continue; 00378 if( (*pit).getType()!=PART_DROP) continue; 00379 // check lifetime ignored here 00380 ParticleObject *p = &(*pit); 00381 const ntlVec3Gfx ppos = p->getPos(); 00382 const int pi= (int)round(ppos[0])+0; 00383 const int pj= (int)round(ppos[1])+0; 00384 int pk= (int)round(ppos[2])+0;// no offset necessary 00385 // 2d should be handled by solver. if(LBMDIM==2) { pk = 0; } 00386 00387 if(pi<0) continue; 00388 if(pj<0) continue; 00389 if(pk<0) continue; 00390 if(pi>mSizex-1) continue; 00391 if(pj>mSizey-1) continue; 00392 if(pk>mSizez-1) continue; 00393 ParticleObject* &pnt = arppnt[ISOLEVEL_INDEX(pi,pj,pk)]; 00394 if(pnt) { 00395 // append 00396 ParticleObject* listpnt = pnt; 00397 while(listpnt) { 00398 if(!listpnt->getNext()) { 00399 listpnt->setNext(p); listpnt = NULL; 00400 } else { 00401 listpnt = listpnt->getNext(); 00402 } 00403 } 00404 } else { 00405 // start new list 00406 pnt = p; 00407 } 00408 pInUse++; 00409 } 00410 } // mpIsoParts 00411 00412 debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<<pInUse<<", Subdivs:"<<mSubdivs, 9); 00413 pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz; 00414 for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) { 00415 pz += gsz; 00416 const int k = ok/mSubdivs; 00417 if(k<=0) continue; // skip zero plane 00418 for(int j=1;j<(mSizey-2);j++) { 00419 for(int i=1;i<(mSizex-2);i++) { 00420 00421 orgval[0] = *getData(i ,j ,k ); 00422 orgval[1] = *getData(i+1,j ,k ); 00423 orgval[2] = *getData(i+1,j+1,k ); // with subdivs 00424 orgval[3] = *getData(i ,j+1,k ); 00425 orgval[4] = *getData(i ,j ,k+1); 00426 orgval[5] = *getData(i+1,j ,k+1); 00427 orgval[6] = *getData(i+1,j+1,k+1); // with subdivs 00428 orgval[7] = *getData(i ,j+1,k+1); 00429 00430 // prebuild subsampled array slice 00431 const int sdkOffset = ok-k*mSubdivs; 00432 for(int sdk=0; sdk<2; sdk++) 00433 for(int sdj=0; sdj<mSubdivs+1; sdj++) 00434 for(int sdi=0; sdi<mSubdivs+1; sdi++) { 00435 subdAr[sdk][sdj][sdi] = ISOTRILININT(sdi*subdfac, sdj*subdfac, (sdkOffset+sdk)*subdfac); 00436 } 00437 00438 const int poDistOffset=2; 00439 for(int pok=-poDistOffset; pok<1+poDistOffset; pok++) { 00440 if(k+pok<0) continue; 00441 if(k+pok>=mSizez-1) continue; 00442 for(int poj=-poDistOffset; poj<1+poDistOffset; poj++) { 00443 if(j+poj<0) continue; 00444 if(j+poj>=mSizey-1) continue; 00445 for(int poi=-poDistOffset; poi<1+poDistOffset; poi++) { 00446 if(i+poi<0) continue; 00447 if(i+poi>=mSizex-1) continue; 00448 ParticleObject *p; 00449 p = arppnt[ISOLEVEL_INDEX(i+poi,j+poj,k+pok)]; 00450 while(p) { // */ 00451 /* 00452 for(vector<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin(); 00453 pit!= mpIsoParts->getParticlesEnd(); pit++) { { { { 00454 // debug test! , full list slow! 00455 if(( (*pit).getActive()==false ) || ( (*pit).getType()!=PART_DROP)) continue; 00456 ParticleObject *p; 00457 p = &(*pit); // */ 00458 00459 pUsedTest++; 00460 ntlVec3Gfx ppos = p->getPos(); 00461 const int spi= (int)round( (ppos[0]+1.-(gfxReal)i) *(gfxReal)mSubdivs-1.5); 00462 const int spj= (int)round( (ppos[1]+1.-(gfxReal)j) *(gfxReal)mSubdivs-1.5); 00463 const int spk= (int)round( (ppos[2]+1.-(gfxReal)k) *(gfxReal)mSubdivs-1.5)-sdkOffset; // why -2? 00464 // 2d should be handled by solver. if(LBMDIM==2) { spk = 0; } 00465 00466 gfxReal pfLen = p->getSize()*1.5*mPartSize; // test, was 1.1 00467 const gfxReal minPfLen = subdfac*0.8; 00468 if(pfLen<minPfLen) pfLen = minPfLen; 00469 //errMsg("ISOPPP"," at "<<PRINT_IJK<<" pp"<<ppos<<" sp"<<PRINT_VEC(spi,spj,spk)<<" pflen"<<pfLen ); 00470 //errMsg("ISOPPP"," subdfac="<<subdfac<<" size"<<p->getSize()<<" ps"<<mPartSize ); 00471 const int icellpsize = (int)(1.*pfLen*(gfxReal)mSubdivs)+1; 00472 for(int swk=-icellpsize; swk<=icellpsize; swk++) { 00473 if(spk+swk< 0) { continue; } 00474 if(spk+swk> 1) { continue; } // */ 00475 for(int swj=-icellpsize; swj<=icellpsize; swj++) { 00476 if(spj+swj< 0) { continue; } 00477 if(spj+swj>mSubdivs+0) { continue; } // */ 00478 for(int swi=-icellpsize; swi<=icellpsize; swi++) { 00479 if(spi+swi< 0) { continue; } 00480 if(spi+swi>mSubdivs+0) { continue; } // */ 00481 ntlVec3Gfx cellp = ntlVec3Gfx( 00482 (1.5+(gfxReal)(spi+swi)) *subdfac + (gfxReal)(i-1), 00483 (1.5+(gfxReal)(spj+swj)) *subdfac + (gfxReal)(j-1), 00484 (1.5+(gfxReal)(spk+swk)+sdkOffset) *subdfac + (gfxReal)(k-1) 00485 ); 00486 //if(swi==0 && swj==0 && swk==0) subdAr[spk][spj][spi] = 1.; // DEBUG 00487 // clip domain boundaries again 00488 if(cellp[0]<1.) { continue; } 00489 if(cellp[1]<1.) { continue; } 00490 if(cellp[2]<1.) { continue; } 00491 if(cellp[0]>(gfxReal)mSizex-3.) { continue; } 00492 if(cellp[1]>(gfxReal)mSizey-3.) { continue; } 00493 if(cellp[2]>(gfxReal)mSizez-3.) { continue; } 00494 gfxReal len = norm(cellp-ppos); 00495 gfxReal isoadd = 0.; 00496 const gfxReal baseIsoVal = mIsoValue*1.1; 00497 if(len<pfLen) { 00498 isoadd = baseIsoVal*1.; 00499 } else { 00500 // falloff linear with pfLen (kernel size=2pfLen 00501 isoadd = baseIsoVal*(1. - (len-pfLen)/(pfLen)); 00502 } 00503 if(isoadd<0.) { continue; } 00504 //errMsg("ISOPPP"," at "<<PRINT_IJK<<" sp"<<PRINT_VEC(spi+swi,spj+swj,spk+swk)<<" cellp"<<cellp<<" pp"<<ppos << " l"<< len<< " add"<< isoadd); 00505 const gfxReal arval = subdAr[spk+swk][spj+swj][spi+swi]; 00506 if(arval>1.) { continue; } 00507 subdAr[spk+swk][spj+swj][spi+swi] = arval + isoadd; 00508 } } } 00509 00510 p = p->getNext(); 00511 } 00512 } } } // poDist loops */ 00513 00514 py = mStart[1]+(((double)j-0.5)*orgGsy)-gsy; 00515 for(int sj=0;sj<mSubdivs;sj++) { 00516 py += gsy; 00517 px = mStart[0]+(((double)i-0.5)*orgGsx)-gsx; 00518 for(int si=0;si<mSubdivs;si++) { 00519 px += gsx; 00520 value[0] = subdAr[0+0][sj+0][si+0]; 00521 value[1] = subdAr[0+0][sj+0][si+1]; 00522 value[2] = subdAr[0+0][sj+1][si+1]; 00523 value[3] = subdAr[0+0][sj+1][si+0]; 00524 value[4] = subdAr[0+1][sj+0][si+0]; 00525 value[5] = subdAr[0+1][sj+0][si+1]; 00526 value[6] = subdAr[0+1][sj+1][si+1]; 00527 value[7] = subdAr[0+1][sj+1][si+0]; 00528 00529 // check intersections of isosurface with edges, and calculate cubie index 00530 cubeIndex = 0; 00531 if (value[0] < mIsoValue) cubeIndex |= 1; 00532 if (value[1] < mIsoValue) cubeIndex |= 2; // with subdivs 00533 if (value[2] < mIsoValue) cubeIndex |= 4; 00534 if (value[3] < mIsoValue) cubeIndex |= 8; 00535 if (value[4] < mIsoValue) cubeIndex |= 16; 00536 if (value[5] < mIsoValue) cubeIndex |= 32; // with subdivs 00537 if (value[6] < mIsoValue) cubeIndex |= 64; 00538 if (value[7] < mIsoValue) cubeIndex |= 128; 00539 00540 if (mcEdgeTable[cubeIndex] > 0) { 00541 00542 // where to look up if this point already exists 00543 const int edgek = 0; 00544 const int baseIn = EDGEAR_INDEX( i+0, j+0, edgek+0, si,sj); 00545 eVert[ 0] = &mpEdgeVerticesX[ baseIn ]; 00546 eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ]; 00547 eVert[ 2] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ]; 00548 eVert[ 3] = &mpEdgeVerticesY[ baseIn ]; 00549 00550 eVert[ 4] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ]; 00551 eVert[ 5] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+1,sj+0) ]; // with subdivs 00552 eVert[ 6] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+1) ]; 00553 eVert[ 7] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ]; 00554 00555 eVert[ 8] = &mpEdgeVerticesZ[ baseIn ]; 00556 eVert[ 9] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+0) ]; // with subdivs 00557 eVert[10] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+1) ]; 00558 eVert[11] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ]; 00559 00560 // grid positions 00561 pos[0] = ntlVec3Gfx(px ,py ,pz); 00562 pos[1] = ntlVec3Gfx(px+gsx,py ,pz); 00563 pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz); // with subdivs 00564 pos[3] = ntlVec3Gfx(px ,py+gsy,pz); 00565 pos[4] = ntlVec3Gfx(px ,py ,pz+gsz); 00566 pos[5] = ntlVec3Gfx(px+gsx,py ,pz+gsz); 00567 pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz); // with subdivs 00568 pos[7] = ntlVec3Gfx(px ,py+gsy,pz+gsz); 00569 00570 // check all edges 00571 for(int e=0;e<12;e++) { 00572 if (mcEdgeTable[cubeIndex] & (1<<e)) { 00573 // is the vertex already calculated? 00574 if(*eVert[ e ] < 0) { 00575 // interpolate edge 00576 const int e1 = mcEdges[e*2 ]; 00577 const int e2 = mcEdges[e*2+1]; 00578 const ntlVec3Gfx p1 = pos[ e1 ]; // scalar field pos 1 00579 const ntlVec3Gfx p2 = pos[ e2 ]; // scalar field pos 2 00580 const float valp1 = value[ e1 ]; // scalar field val 1 00581 const float valp2 = value[ e2 ]; // scalar field val 2 00582 const float mu = (mIsoValue - valp1) / (valp2 - valp1); 00583 00584 // init isolevel vertex 00585 ilv.v = p1 + (p2-p1)*mu; // with subdivs 00586 mPoints.push_back( ilv ); 00587 triIndices[e] = (mPoints.size()-1); 00588 // store vertex 00589 *eVert[ e ] = triIndices[e]; 00590 } else { 00591 // retrieve from vert array 00592 triIndices[e] = *eVert[ e ]; 00593 } 00594 } // along all edges 00595 } 00596 // removed cutoff treatment... 00597 00598 // Create the triangles... 00599 for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) { 00600 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] ); 00601 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); // with subdivs 00602 mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] ); 00603 //errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() ); 00604 } 00605 00606 } // triangles in edge table? 00607 00608 }//si 00609 }// sj 00610 00611 }//i 00612 }// j 00613 00614 // copy edge arrays 00615 for(int j=0;j<(mSizey-0)*mSubdivs;j++) 00616 for(int i=0;i<(mSizex-0)*mSubdivs;i++) { 00617 //int edgek = 0; 00618 const int dst = EDGEAR_INDEX( 0, 0, 0, i,j); 00619 const int src = EDGEAR_INDEX( 0, 0, 1, i,j); 00620 mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ]; 00621 mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; // with subdivs 00622 mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ]; 00623 mpEdgeVerticesX[ src ]=-1; 00624 mpEdgeVerticesY[ src ]=-1; // with subdivs 00625 mpEdgeVerticesZ[ src ]=-1; 00626 } 00627 // */ 00628 00629 } // ok, k subdiv loop 00630 00631 //delete [] subdAr; 00632 delete [] arppnt; 00633 computeNormals(); 00634 } // with subdivs 00635 00636 // perform smoothing 00637 float smoSubdfac = 1.; 00638 if(mSubdivs>0) { 00639 //smoSubdfac = 1./(float)(mSubdivs); 00640 smoSubdfac = pow(0.55,(double)mSubdivs); // slightly stronger 00641 } 00642 if(mSmoothSurface>0. || mSmoothNormals>0.) debMsgStd("IsoSurface::triangulate",DM_MSG,"Smoothing...",10); 00643 if(mSmoothSurface>0.0) { 00644 smoothSurface(mSmoothSurface*smoSubdfac, (mSmoothNormals<=0.0) ); 00645 } 00646 if(mSmoothNormals>0.0) { 00647 smoothNormals(mSmoothNormals*smoSubdfac); 00648 } 00649 00650 myTime_t tritimeend = getTime(); 00651 debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< getTimeString(tritimeend-tritimestart)<<", S("<<mSmoothSurface<<","<<mSmoothNormals<<"),"<< 00652 " verts:"<<mPoints.size()<<" tris:"<<(mIndices.size()/3)<<" subdivs:"<<mSubdivs 00653 , 10 ); 00654 if(mpIsoParts) debMsgStd("IsoSurface::triangulate",DM_MSG,"parts:"<<mpIsoParts->getNumParticles(), 10); 00655 } 00656 00657 00658 00659 00660 00661 /****************************************************************************** 00662 * Get triangles for rendering 00663 *****************************************************************************/ 00664 void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles, 00665 vector<ntlVec3Gfx> *vertices, 00666 vector<ntlVec3Gfx> *normals, int objectId ) 00667 { 00668 if(!mInitDone) { 00669 debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10); 00670 return; 00671 } 00672 t = 0.; 00673 //return; // DEBUG 00674 00675 /* triangulate field */ 00676 triangulate(); 00677 //errMsg("TRIS"," "<<mIndices.size() ); 00678 00679 // new output with vertice reuse 00680 int iniVertIndex = (*vertices).size(); 00681 int iniNormIndex = (*normals).size(); 00682 if(iniVertIndex != iniNormIndex) { 00683 errFatal("getTriangles Error","For '"<<mName<<"': Vertices and normal array sizes to not match!!!",SIMWORLD_GENERICERROR); 00684 return; 00685 } 00686 //errMsg("NM"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() ); 00687 //errMsg("NM"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() ); 00688 00689 for(int i=0;i<(int)mPoints.size();i++) { 00690 vertices->push_back( mPoints[i].v ); 00691 } 00692 for(int i=0;i<(int)mPoints.size();i++) { 00693 normals->push_back( mPoints[i].n ); 00694 } 00695 00696 //errMsg("N2"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() ); 00697 //errMsg("N2"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() ); 00698 00699 for(int i=0;i<(int)mIndices.size();i+=3) { 00700 const int smooth = 1; 00701 int t1 = mIndices[i]; 00702 int t2 = mIndices[i+1]; 00703 int t3 = mIndices[i+2]; 00704 //errMsg("NM"," tri"<<t1<<" "<<t2<<" "<<t3 ); 00705 00706 ntlTriangle tri; 00707 00708 tri.getPoints()[0] = t1+iniVertIndex; 00709 tri.getPoints()[1] = t2+iniVertIndex; 00710 tri.getPoints()[2] = t3+iniVertIndex; 00711 00712 /* init flags */ 00713 int flag = 0; 00714 if(getVisible()){ flag |= TRI_GEOMETRY; } 00715 if(getCastShadows() ) { 00716 flag |= TRI_CASTSHADOWS; } 00717 00718 /* init geo init id */ 00719 int geoiId = getGeoInitId(); 00720 if(geoiId > 0) { 00721 flag |= (1<< (geoiId+4)); 00722 flag |= mGeoInitType; 00723 } 00724 00725 tri.setFlags( flag ); 00726 00727 /* triangle normal missing */ 00728 tri.setNormal( ntlVec3Gfx(0.0) ); 00729 tri.setSmoothNormals( smooth ); 00730 tri.setObjectId( objectId ); 00731 triangles->push_back( tri ); 00732 } 00733 //errMsg("N3"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() ); 00734 return; 00735 } 00736 00737 00738 00739 inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) { 00740 // WARNING - this requires a security boundary layer... 00741 ntlVec3Gfx ret(0.0); 00742 ret[0] = *getData(i-1,j ,k ) - 00743 *getData(i+1,j ,k ); 00744 ret[1] = *getData(i ,j-1,k ) - 00745 *getData(i ,j+1,k ); 00746 ret[2] = *getData(i ,j ,k-1 ) - 00747 *getData(i ,j ,k+1 ); 00748 return ret; 00749 } 00750 00751 00752 00753 00754 /****************************************************************************** 00755 * 00756 * Surface improvement, inspired by trimesh2 library 00757 * (http://www.cs.princeton.edu/gfx/proj/trimesh2/) 00758 * 00759 *****************************************************************************/ 00760 00761 void IsoSurface::setSmoothRad(float radi1, float radi2, ntlVec3Gfx mscc) { 00762 mSCrad1 = radi1*radi1; 00763 mSCrad2 = radi2*radi2; 00764 mSCcenter = mscc; 00765 } 00766 00767 // compute normals for all generated triangles 00768 void IsoSurface::computeNormals() { 00769 for(int i=0;i<(int)mPoints.size();i++) { 00770 mPoints[i].n = ntlVec3Gfx(0.); 00771 } 00772 00773 for(int i=0;i<(int)mIndices.size();i+=3) { 00774 const int t1 = mIndices[i]; 00775 const int t2 = mIndices[i+1]; 00776 const int t3 = mIndices[i+2]; 00777 const ntlVec3Gfx p1 = mPoints[t1].v; 00778 const ntlVec3Gfx p2 = mPoints[t2].v; 00779 const ntlVec3Gfx p3 = mPoints[t3].v; 00780 const ntlVec3Gfx n1=p1-p2; 00781 const ntlVec3Gfx n2=p2-p3; 00782 const ntlVec3Gfx n3=p3-p1; 00783 const gfxReal len1 = normNoSqrt(n1); 00784 const gfxReal len2 = normNoSqrt(n2); 00785 const gfxReal len3 = normNoSqrt(n3); 00786 const ntlVec3Gfx norm = cross(n1,n2); 00787 mPoints[t1].n += norm * (1./(len1*len3)); 00788 mPoints[t2].n += norm * (1./(len1*len2)); 00789 mPoints[t3].n += norm * (1./(len2*len3)); 00790 } 00791 00792 for(int i=0;i<(int)mPoints.size();i++) { 00793 normalize(mPoints[i].n); 00794 } 00795 } 00796 00797 // Diffuse a vector field at 1 vertex, weighted by 00798 // a gaussian of width 1/sqrt(invsigma2) 00799 bool IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int src, float invsigma2, ntlVec3Gfx &target) 00800 { 00801 if((neighbors[src].size()<1) || (pointareas[src]<=0.0)) return 0; 00802 const ntlVec3Gfx srcp = mPoints[src].v; 00803 const ntlVec3Gfx srcn = mPoints[src].n; 00804 if(mSCrad1>0.0 && mSCrad2>0.0) { 00805 ntlVec3Gfx dp = mSCcenter-srcp; dp[2] = 0.0; // only xy-plane 00806 float rd = normNoSqrt(dp); 00807 if(rd > mSCrad2) { 00808 return 0; 00809 } else if(rd > mSCrad1) { 00810 // optimize? 00811 float org = 1.0/sqrt(invsigma2); 00812 org *= (1.0- (rd-mSCrad1) / (mSCrad2-mSCrad1)); 00813 invsigma2 = 1.0/(org*org); 00814 //errMsg("TRi","p"<<srcp<<" rd:"<<rd<<" r1:"<<mSCrad1<<" r2:"<<mSCrad2<<" org:"<<org<<" is:"<<invsigma2); 00815 } else { 00816 } 00817 } 00818 target = ntlVec3Gfx(0.0); 00819 target += *(field+pointerScale*src) *pointareas[src]; 00820 float smstrSum = pointareas[src]; 00821 00822 int flag = mFlagCnt; 00823 mFlagCnt++; 00824 flags[src] = flag; 00825 mDboundary = neighbors[src]; 00826 while (!mDboundary.empty()) { 00827 const int bbn = mDboundary.back(); 00828 mDboundary.pop_back(); 00829 if(flags[bbn]==flag) continue; 00830 flags[bbn] = flag; 00831 00832 // normal check 00833 const float nvdot = dot(srcn, mPoints[bbn].n); // faster than before d2 calc? 00834 if(nvdot <= 0.0f) continue; 00835 00836 // gaussian weight of width 1/sqrt(invsigma2) 00837 const float d2 = invsigma2 * normNoSqrt(mPoints[bbn].v - srcp); 00838 if(d2 >= 9.0f) continue; 00839 00840 // aggressive smoothing factor 00841 float smstr = nvdot * pointareas[bbn]; 00842 // Accumulate weight times field at neighbor 00843 target += *(field+pointerScale*bbn)*smstr; 00844 smstrSum += smstr; 00845 00846 for(int i = 0; i < (int)neighbors[bbn].size(); i++) { 00847 const int nn = neighbors[bbn][i]; 00848 if (flags[nn] == flag) continue; 00849 mDboundary.push_back(nn); 00850 } 00851 } 00852 target /= smstrSum; 00853 return 1; 00854 } 00855 00856 00857 // perform smoothing of the surface (and possible normals) 00858 void IsoSurface::smoothSurface(float sigma, bool normSmooth) 00859 { 00860 int nv = mPoints.size(); 00861 if ((int)flags.size() != nv) flags.resize(nv); 00862 int nf = mIndices.size()/3; 00863 00864 { // need neighbors 00865 vector<int> numneighbors(mPoints.size()); 00866 int i; 00867 for (i = 0; i < (int)mIndices.size()/3; i++) { 00868 numneighbors[mIndices[i*3+0]]++; 00869 numneighbors[mIndices[i*3+1]]++; 00870 numneighbors[mIndices[i*3+2]]++; 00871 } 00872 00873 neighbors.clear(); 00874 neighbors.resize(mPoints.size()); 00875 for (i = 0; i < (int)mPoints.size(); i++) { 00876 neighbors[i].clear(); 00877 neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries 00878 } 00879 00880 for (i = 0; i < (int)mIndices.size()/3; i++) { 00881 for (int j = 0; j < 3; j++) { 00882 vector<int> &me = neighbors[ mIndices[i*3+j]]; 00883 int n1 = mIndices[i*3+((j+1)%3)]; 00884 int n2 = mIndices[i*3+((j+2)%3)]; 00885 if (std::find(me.begin(), me.end(), n1) == me.end()) 00886 me.push_back(n1); 00887 if (std::find(me.begin(), me.end(), n2) == me.end()) 00888 me.push_back(n2); 00889 } 00890 } 00891 } // need neighbor 00892 00893 { // need pointarea 00894 pointareas.clear(); 00895 pointareas.resize(nv); 00896 cornerareas.clear(); 00897 cornerareas.resize(nf); 00898 00899 for (int i = 0; i < nf; i++) { 00900 // Edges 00901 ntlVec3Gfx e[3] = { 00902 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v, 00903 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v, 00904 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v }; 00905 00906 // Compute corner weights 00907 float area = 0.5f * norm( cross(e[0], e[1])); 00908 float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) }; 00909 float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]), 00910 l2[1] * (l2[2] + l2[0] - l2[1]), 00911 l2[2] * (l2[0] + l2[1] - l2[2]) }; 00912 if (ew[0] <= 0.0f) { 00913 cornerareas[i][1] = -0.25f * l2[2] * area / 00914 dot(e[0] , e[2]); 00915 cornerareas[i][2] = -0.25f * l2[1] * area / 00916 dot(e[0] , e[1]); 00917 cornerareas[i][0] = area - cornerareas[i][1] - 00918 cornerareas[i][2]; 00919 } else if (ew[1] <= 0.0f) { 00920 cornerareas[i][2] = -0.25f * l2[0] * area / 00921 dot(e[1] , e[0]); 00922 cornerareas[i][0] = -0.25f * l2[2] * area / 00923 dot(e[1] , e[2]); 00924 cornerareas[i][1] = area - cornerareas[i][2] - 00925 cornerareas[i][0]; 00926 } else if (ew[2] <= 0.0f) { 00927 cornerareas[i][0] = -0.25f * l2[1] * area / 00928 dot(e[2] , e[1]); 00929 cornerareas[i][1] = -0.25f * l2[0] * area / 00930 dot(e[2] , e[0]); 00931 cornerareas[i][2] = area - cornerareas[i][0] - 00932 cornerareas[i][1]; 00933 } else { 00934 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]); 00935 for (int j = 0; j < 3; j++) 00936 cornerareas[i][j] = ewscale * (ew[(j+1)%3] + 00937 ew[(j+2)%3]); 00938 } 00939 00940 // NT important, check this... 00941 #ifndef WIN32 00942 if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6; 00943 if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6; 00944 if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6; 00945 #else // WIN32 00946 // FIXME check as well... 00947 if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6; 00948 if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6; 00949 if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6; 00950 #endif // WIN32 00951 00952 pointareas[mIndices[i*3+0]] += cornerareas[i][0]; 00953 pointareas[mIndices[i*3+1]] += cornerareas[i][1]; 00954 pointareas[mIndices[i*3+2]] += cornerareas[i][2]; 00955 } 00956 00957 } // need pointarea 00958 // */ 00959 00960 float invsigma2 = 1.0f / (sigma*sigma); 00961 00962 vector<ntlVec3Gfx> dflt(nv); 00963 for (int i = 0; i < nv; i++) { 00964 if(diffuseVertexField( &mPoints[0].v, 2, 00965 i, invsigma2, dflt[i])) { 00966 // Just keep the displacement 00967 dflt[i] -= mPoints[i].v; 00968 } else { dflt[i] = 0.0; } //?mPoints[i].v; } 00969 } 00970 00971 // Slightly better small-neighborhood approximation 00972 for (int i = 0; i < nf; i++) { 00973 ntlVec3Gfx c = mPoints[mIndices[i*3+0]].v + 00974 mPoints[mIndices[i*3+1]].v + 00975 mPoints[mIndices[i*3+2]].v; 00976 c /= 3.0f; 00977 for (int j = 0; j < 3; j++) { 00978 int v = mIndices[i*3+j]; 00979 ntlVec3Gfx d =(c - mPoints[v].v) * 0.5f; 00980 dflt[v] += d * (cornerareas[i][j] / 00981 pointareas[mIndices[i*3+j]] * 00982 exp(-0.5f * invsigma2 * normNoSqrt(d)) ); 00983 } 00984 } 00985 00986 // Filter displacement field 00987 vector<ntlVec3Gfx> dflt2(nv); 00988 for (int i = 0; i < nv; i++) { 00989 if(diffuseVertexField( &dflt[0], 1, 00990 i, invsigma2, dflt2[i])) { } 00991 else { /*mPoints[i].v=0.0;*/ dflt2[i] = 0.0; }//dflt2[i]; } 00992 } 00993 00994 // Update vertex positions 00995 for (int i = 0; i < nv; i++) { 00996 mPoints[i].v += dflt[i] - dflt2[i]; // second Laplacian 00997 } 00998 00999 // when normals smoothing off, this cleans up quite well 01000 // costs ca. 50% additional time though 01001 float nsFac = 1.5f; 01002 if(normSmooth) { float ninvsigma2 = 1.0f / (nsFac*nsFac*sigma*sigma); 01003 for (int i = 0; i < nv; i++) { 01004 if( diffuseVertexField( &mPoints[0].n, 2, i, ninvsigma2, dflt[i]) ) { 01005 normalize(dflt[i]); 01006 } else { 01007 dflt[i] = mPoints[i].n; 01008 } 01009 } 01010 for (int i = 0; i < nv; i++) { 01011 mPoints[i].n = dflt[i]; 01012 } 01013 } // smoothNormals copy */ 01014 01015 //errMsg("SMSURF","done v:"<<sigma); // DEBUG 01016 } 01017 01018 // only smoothen the normals 01019 void IsoSurface::smoothNormals(float sigma) { 01020 // reuse from smoothSurface 01021 if(neighbors.size() != mPoints.size()) { 01022 // need neighbor 01023 vector<int> numneighbors(mPoints.size()); 01024 int i; 01025 for (i = 0; i < (int)mIndices.size()/3; i++) { 01026 numneighbors[mIndices[i*3+0]]++; 01027 numneighbors[mIndices[i*3+1]]++; 01028 numneighbors[mIndices[i*3+2]]++; 01029 } 01030 01031 neighbors.clear(); 01032 neighbors.resize(mPoints.size()); 01033 for (i = 0; i < (int)mPoints.size(); i++) { 01034 neighbors[i].clear(); 01035 neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries 01036 } 01037 01038 for (i = 0; i < (int)mIndices.size()/3; i++) { 01039 for (int j = 0; j < 3; j++) { 01040 vector<int> &me = neighbors[ mIndices[i*3+j]]; 01041 int n1 = mIndices[i*3+((j+1)%3)]; 01042 int n2 = mIndices[i*3+((j+2)%3)]; 01043 if (std::find(me.begin(), me.end(), n1) == me.end()) 01044 me.push_back(n1); 01045 if (std::find(me.begin(), me.end(), n2) == me.end()) 01046 me.push_back(n2); 01047 } 01048 } 01049 } // need neighbor 01050 01051 { // need pointarea 01052 int nf = mIndices.size()/3, nv = mPoints.size(); 01053 pointareas.clear(); 01054 pointareas.resize(nv); 01055 cornerareas.clear(); 01056 cornerareas.resize(nf); 01057 01058 for (int i = 0; i < nf; i++) { 01059 // Edges 01060 ntlVec3Gfx e[3] = { 01061 mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v, 01062 mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v, 01063 mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v }; 01064 01065 // Compute corner weights 01066 float area = 0.5f * norm( cross(e[0], e[1])); 01067 float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) }; 01068 float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]), 01069 l2[1] * (l2[2] + l2[0] - l2[1]), 01070 l2[2] * (l2[0] + l2[1] - l2[2]) }; 01071 if (ew[0] <= 0.0f) { 01072 cornerareas[i][1] = -0.25f * l2[2] * area / 01073 dot(e[0] , e[2]); 01074 cornerareas[i][2] = -0.25f * l2[1] * area / 01075 dot(e[0] , e[1]); 01076 cornerareas[i][0] = area - cornerareas[i][1] - 01077 cornerareas[i][2]; 01078 } else if (ew[1] <= 0.0f) { 01079 cornerareas[i][2] = -0.25f * l2[0] * area / 01080 dot(e[1] , e[0]); 01081 cornerareas[i][0] = -0.25f * l2[2] * area / 01082 dot(e[1] , e[2]); 01083 cornerareas[i][1] = area - cornerareas[i][2] - 01084 cornerareas[i][0]; 01085 } else if (ew[2] <= 0.0f) { 01086 cornerareas[i][0] = -0.25f * l2[1] * area / 01087 dot(e[2] , e[1]); 01088 cornerareas[i][1] = -0.25f * l2[0] * area / 01089 dot(e[2] , e[0]); 01090 cornerareas[i][2] = area - cornerareas[i][0] - 01091 cornerareas[i][1]; 01092 } else { 01093 float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]); 01094 for (int j = 0; j < 3; j++) 01095 cornerareas[i][j] = ewscale * (ew[(j+1)%3] + 01096 ew[(j+2)%3]); 01097 } 01098 01099 // NT important, check this... 01100 #ifndef WIN32 01101 if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6; 01102 if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6; 01103 if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6; 01104 #else // WIN32 01105 // FIXME check as well... 01106 if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6; 01107 if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6; 01108 if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6; 01109 #endif // WIN32 01110 01111 pointareas[mIndices[i*3+0]] += cornerareas[i][0]; 01112 pointareas[mIndices[i*3+1]] += cornerareas[i][1]; 01113 pointareas[mIndices[i*3+2]] += cornerareas[i][2]; 01114 } 01115 01116 } // need pointarea 01117 01118 int nv = mPoints.size(); 01119 if ((int)flags.size() != nv) flags.resize(nv); 01120 float invsigma2 = 1.0f / (sigma*sigma); 01121 01122 vector<ntlVec3Gfx> nflt(nv); 01123 for (int i = 0; i < nv; i++) { 01124 if(diffuseVertexField( &mPoints[0].n, 2, i, invsigma2, nflt[i])) { 01125 normalize(nflt[i]); 01126 } else { nflt[i]=mPoints[i].n; } 01127 } 01128 01129 // copy results 01130 for (int i = 0; i < nv; i++) { mPoints[i].n = nflt[i]; } 01131 } 01132 01133