diff --git a/intern/elbeem/intern/isosurface.cpp b/intern/elbeem/intern/isosurface.cpp index 77530d413c7..704f3fdacfe 100644 --- a/intern/elbeem/intern/isosurface.cpp +++ b/intern/elbeem/intern/isosurface.cpp @@ -18,28 +18,32 @@ #define round(x) (x) #endif +#if PARALLEL==1 +#include +#endif + /****************************************************************************** * Constructor *****************************************************************************/ IsoSurface::IsoSurface(double iso) : - ntlGeometryObject(), - mSizex(-1), mSizey(-1), mSizez(-1), - mpData(NULL), - mIsoValue( iso ), - mPoints(), - mUseFullEdgeArrays(false), - mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL), - mEdgeArSize(-1), - mIndices(), + ntlGeometryObject(), + mSizex(-1), mSizey(-1), mSizez(-1), + mpData(NULL), + mIsoValue( iso ), + mPoints(), + mUseFullEdgeArrays(false), + mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL), + mEdgeArSize(-1), + mIndices(), - mStart(0.0), mEnd(0.0), mDomainExtent(0.0), - mInitDone(false), - mSmoothSurface(0.0), mSmoothNormals(0.0), - mAcrossEdge(), mAdjacentFaces(), - mCutoff(-1), mCutArray(NULL), // off by default - mpIsoParts(NULL), mPartSize(0.), mSubdivs(0), - mFlagCnt(1), - mSCrad1(0.), mSCrad2(0.), mSCcenter(0.) + mStart(0.0), mEnd(0.0), mDomainExtent(0.0), + mInitDone(false), + mSmoothSurface(0.0), mSmoothNormals(0.0), + mAcrossEdge(), mAdjacentFaces(), + mCutoff(-1), mCutArray(NULL), // off by default + mpIsoParts(NULL), mPartSize(0.), mSubdivs(0), + mFlagCnt(1), + mSCrad1(0.), mSCrad2(0.), mSCcenter(0.) { } @@ -71,11 +75,11 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e // init mIndices.clear(); - mPoints.clear(); + mPoints.clear(); int nodes = mSizez*mSizey*mSizex; - mpData = new float[nodes]; - for(int i=0;i0) && (kmSizex-2-coAdd-mCutoff) || - (j>mSizey-2-coAdd-mCutoff) ) { - if(mCutArray) { - if(k < mCutArray[j*this->mSizex+i]) continue; - } else { continue; } - } + if( (i0) && (kmSizex-2-coAdd-mCutoff) || + (j>mSizey-2-coAdd-mCutoff) ) { + if(mCutArray) { + if(k < mCutArray[j*this->mSizex+i]) continue; + } else { continue; } + } // Create the triangles... - for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) { - mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] ); - mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); - mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] ); - } + for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) { + mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] ); + mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); + mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] ); + } - }//i - }// j + }//i + }// j // copy edge arrays - if(!mUseFullEdgeArrays) { - for(int j=0;j<(mSizey-0);j++) - for(int i=0;i<(mSizex-0);i++) { + if(!mUseFullEdgeArrays) { + for(int j=0;j<(mSizey-0);j++) + for(int i=0;i<(mSizex-0);i++) { //int edgek = 0; - const int dst = ISOLEVEL_INDEX( i+0, j+0, 0); - const int src = ISOLEVEL_INDEX( i+0, j+0, 1); - mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ]; - mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; - mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ]; - mpEdgeVerticesX[ src ]=-1; - mpEdgeVerticesY[ src ]=-1; - mpEdgeVerticesZ[ src ]=-1; - } - } // */ + const int dst = ISOLEVEL_INDEX( i+0, j+0, 0); + const int src = ISOLEVEL_INDEX( i+0, j+0, 1); + mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ]; + mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; + mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ]; + mpEdgeVerticesX[ src ]=-1; + mpEdgeVerticesY[ src ]=-1; + mpEdgeVerticesZ[ src ]=-1; + } + } // */ - } // k + } // k // precalculate normals using an approximation of the scalar field gradient - for(int ni=0;ni<(int)mPoints.size();ni++) { normalize( mPoints[ni].n ); } + for(int ni=0;ni<(int)mPoints.size();ni++) { normalize( mPoints[ni].n ); } - } else { // subdivs + } else { // subdivs #define EDGEAR_INDEX(Ai,Aj,Ak, Bi,Bj) ((mSizex*mSizey*mSubdivs*mSubdivs*(Ak))+\ - (mSizex*mSubdivs*((Aj)*mSubdivs+(Bj)))+((Ai)*mSubdivs)+(Bi)) + (mSizex*mSubdivs*((Aj)*mSubdivs+(Bj)))+((Ai)*mSubdivs)+(Bi)) #define ISOTRILININT(fi,fj,fk) ( \ - (1.-(fi))*(1.-(fj))*(1.-(fk))*orgval[0] + \ - ( (fi))*(1.-(fj))*(1.-(fk))*orgval[1] + \ - ( (fi))*( (fj))*(1.-(fk))*orgval[2] + \ - (1.-(fi))*( (fj))*(1.-(fk))*orgval[3] + \ - (1.-(fi))*(1.-(fj))*( (fk))*orgval[4] + \ - ( (fi))*(1.-(fj))*( (fk))*orgval[5] + \ - ( (fi))*( (fj))*( (fk))*orgval[6] + \ - (1.-(fi))*( (fj))*( (fk))*orgval[7] ) + (1.-(fi))*(1.-(fj))*(1.-(fk))*orgval[0] + \ + ( (fi))*(1.-(fj))*(1.-(fk))*orgval[1] + \ + ( (fi))*( (fj))*(1.-(fk))*orgval[2] + \ + (1.-(fi))*( (fj))*(1.-(fk))*orgval[3] + \ + (1.-(fi))*(1.-(fj))*( (fk))*orgval[4] + \ + ( (fi))*(1.-(fj))*( (fk))*orgval[5] + \ + ( (fi))*( (fj))*( (fk))*orgval[6] + \ + (1.-(fi))*( (fj))*( (fk))*orgval[7] ) // use subdivisions - gfxReal subdfac = 1./(gfxReal)(mSubdivs); - gfxReal orgGsx = gsx; - gfxReal orgGsy = gsy; - gfxReal orgGsz = gsz; - gsx *= subdfac; - gsy *= subdfac; - gsz *= subdfac; - if(mUseFullEdgeArrays) { - errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!"); - } + gfxReal subdfac = 1./(gfxReal)(mSubdivs); + gfxReal orgGsx = gsx; + gfxReal orgGsy = gsy; + gfxReal orgGsz = gsz; + gsx *= subdfac; + gsy *= subdfac; + gsz *= subdfac; + if(mUseFullEdgeArrays) { + errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!"); + } - ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex]; + ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex]; // construct pointers // part test - int pInUse = 0; - int pUsedTest = 0; + int pInUse = 0; + int pUsedTest = 0; // reset particles // reset list array - for(int k=0;k<(mSizez);k++) - for(int j=0;j<(mSizey);j++) - for(int i=0;i<(mSizex);i++) { - arppnt[ISOLEVEL_INDEX(i,j,k)] = NULL; - } - if(mpIsoParts) { - for(vector::iterator pit= mpIsoParts->getParticlesBegin(); - pit!= mpIsoParts->getParticlesEnd(); pit++) { - if( (*pit).getActive()==false ) continue; - if( (*pit).getType()!=PART_DROP) continue; - (*pit).setNext(NULL); - } + for(int k=0;k<(mSizez);k++) + for(int j=0;j<(mSizey);j++) + for(int i=0;i<(mSizex);i++) { + arppnt[ISOLEVEL_INDEX(i,j,k)] = NULL; + } + if(mpIsoParts) { + for(vector::iterator pit= mpIsoParts->getParticlesBegin(); + pit!= mpIsoParts->getParticlesEnd(); pit++) { + if( (*pit).getActive()==false ) continue; + if( (*pit).getType()!=PART_DROP) continue; + (*pit).setNext(NULL); + } // build per node lists - for(vector::iterator pit= mpIsoParts->getParticlesBegin(); - pit!= mpIsoParts->getParticlesEnd(); pit++) { - if( (*pit).getActive()==false ) continue; - if( (*pit).getType()!=PART_DROP) continue; + for(vector::iterator pit= mpIsoParts->getParticlesBegin(); + pit!= mpIsoParts->getParticlesEnd(); pit++) { + if( (*pit).getActive()==false ) continue; + if( (*pit).getType()!=PART_DROP) continue; // check lifetime ignored here - ParticleObject *p = &(*pit); - const ntlVec3Gfx ppos = p->getPos(); - const int pi= (int)round(ppos[0])+0; - const int pj= (int)round(ppos[1])+0; - int pk= (int)round(ppos[2])+0;// no offset necessary + ParticleObject *p = &(*pit); + const ntlVec3Gfx ppos = p->getPos(); + const int pi= (int)round(ppos[0])+0; + const int pj= (int)round(ppos[1])+0; + int pk= (int)round(ppos[2])+0;// no offset necessary // 2d should be handled by solver. if(LBMDIM==2) { pk = 0; } - if(pi<0) continue; - if(pj<0) continue; - if(pk<0) continue; - if(pi>mSizex-1) continue; - if(pj>mSizey-1) continue; - if(pk>mSizez-1) continue; - ParticleObject* &pnt = arppnt[ISOLEVEL_INDEX(pi,pj,pk)]; - if(pnt) { + if(pi<0) continue; + if(pj<0) continue; + if(pk<0) continue; + if(pi>mSizex-1) continue; + if(pj>mSizey-1) continue; + if(pk>mSizez-1) continue; + ParticleObject* &pnt = arppnt[ISOLEVEL_INDEX(pi,pj,pk)]; + if(pnt) { // append - ParticleObject* listpnt = pnt; - while(listpnt) { - if(!listpnt->getNext()) { - listpnt->setNext(p); listpnt = NULL; - } else { - listpnt = listpnt->getNext(); - } - } - } else { + ParticleObject* listpnt = pnt; + while(listpnt) { + if(!listpnt->getNext()) { + listpnt->setNext(p); listpnt = NULL; + } else { + listpnt = listpnt->getNext(); + } + } + } else { // start new list - pnt = p; - } - pInUse++; - } - } // mpIsoParts + pnt = p; + } + pInUse++; + } + } // mpIsoParts - debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<=mSizez-1) continue; - for(int poj=-poDistOffset; poj<1+poDistOffset; poj++) { - if(j+poj<0) continue; - if(j+poj>=mSizey-1) continue; - for(int poi=-poDistOffset; poi<1+poDistOffset; poi++) { - if(i+poi<0) continue; - if(i+poi>=mSizex-1) continue; - ParticleObject *p; - p = arppnt[ISOLEVEL_INDEX(i+poi,j+poj,k+pok)]; - while(p) { // */ + const int poDistOffset=2; + for(int pok=-poDistOffset; pok<1+poDistOffset; pok++) { + if(k+pok<0) continue; + if(k+pok>=mSizez-1) continue; + for(int poj=-poDistOffset; poj<1+poDistOffset; poj++) { + if(j+poj<0) continue; + if(j+poj>=mSizey-1) continue; + for(int poi=-poDistOffset; poi<1+poDistOffset; poi++) { + if(i+poi<0) continue; + if(i+poi>=mSizex-1) continue; + ParticleObject *p; + p = arppnt[ISOLEVEL_INDEX(i+poi,j+poj,k+pok)]; + while(p) { // */ /* - for(vector::iterator pit= mpIsoParts->getParticlesBegin(); - pit!= mpIsoParts->getParticlesEnd(); pit++) { { { { + for(vector::iterator pit= mpIsoParts->getParticlesBegin(); + pit!= mpIsoParts->getParticlesEnd(); pit++) { { { { // debug test! , full list slow! - if(( (*pit).getActive()==false ) || ( (*pit).getType()!=PART_DROP)) continue; - ParticleObject *p; - p = &(*pit); // */ + if(( (*pit).getActive()==false ) || ( (*pit).getType()!=PART_DROP)) continue; + ParticleObject *p; + p = &(*pit); // */ - pUsedTest++; - ntlVec3Gfx ppos = p->getPos(); - const int spi= (int)round( (ppos[0]+1.-(gfxReal)i) *(gfxReal)mSubdivs-1.5); - const int spj= (int)round( (ppos[1]+1.-(gfxReal)j) *(gfxReal)mSubdivs-1.5); - const int spk= (int)round( (ppos[2]+1.-(gfxReal)k) *(gfxReal)mSubdivs-1.5)-sdkOffset; // why -2? + pUsedTest++; + ntlVec3Gfx ppos = p->getPos(); + const int spi= (int)round( (ppos[0]+1.-(gfxReal)i) *(gfxReal)mSubdivs-1.5); + const int spj= (int)round( (ppos[1]+1.-(gfxReal)j) *(gfxReal)mSubdivs-1.5); + const int spk= (int)round( (ppos[2]+1.-(gfxReal)k) *(gfxReal)mSubdivs-1.5)-sdkOffset; // why -2? // 2d should be handled by solver. if(LBMDIM==2) { spk = 0; } - gfxReal pfLen = p->getSize()*1.5*mPartSize; // test, was 1.1 - const gfxReal minPfLen = subdfac*0.8; - if(pfLengetSize()*1.5*mPartSize; // test, was 1.1 + const gfxReal minPfLen = subdfac*0.8; + if(pfLengetSize()<<" ps"< 1) { continue; } // */ - for(int swj=-icellpsize; swj<=icellpsize; swj++) { - if(spj+swj< 0) { continue; } - if(spj+swj>mSubdivs+0) { continue; } // */ - for(int swi=-icellpsize; swi<=icellpsize; swi++) { - if(spi+swi< 0) { continue; } - if(spi+swi>mSubdivs+0) { continue; } // */ - ntlVec3Gfx cellp = ntlVec3Gfx( + const int icellpsize = (int)(1.*pfLen*(gfxReal)mSubdivs)+1; + for(int swk=-icellpsize; swk<=icellpsize; swk++) { + if(spk+swk< 0) { continue; } + if(spk+swk> 1) { continue; } // */ + for(int swj=-icellpsize; swj<=icellpsize; swj++) { + if(spj+swj< 0) { continue; } + if(spj+swj>mSubdivs+0) { continue; } // */ + for(int swi=-icellpsize; swi<=icellpsize; swi++) { + if(spi+swi< 0) { continue; } + if(spi+swi>mSubdivs+0) { continue; } // */ + ntlVec3Gfx cellp = ntlVec3Gfx( (1.5+(gfxReal)(spi+swi)) *subdfac + (gfxReal)(i-1), (1.5+(gfxReal)(spj+swj)) *subdfac + (gfxReal)(j-1), (1.5+(gfxReal)(spk+swk)+sdkOffset) *subdfac + (gfxReal)(k-1) ); //if(swi==0 && swj==0 && swk==0) subdAr[spk][spj][spi] = 1.; // DEBUG // clip domain boundaries again - if(cellp[0]<1.) { continue; } - if(cellp[1]<1.) { continue; } - if(cellp[2]<1.) { continue; } - if(cellp[0]>(gfxReal)mSizex-3.) { continue; } - if(cellp[1]>(gfxReal)mSizey-3.) { continue; } - if(cellp[2]>(gfxReal)mSizez-3.) { continue; } - gfxReal len = norm(cellp-ppos); - gfxReal isoadd = 0.; - const gfxReal baseIsoVal = mIsoValue*1.1; - if(len(gfxReal)mSizex-3.) { continue; } + if(cellp[1]>(gfxReal)mSizey-3.) { continue; } + if(cellp[2]>(gfxReal)mSizez-3.) { continue; } + gfxReal len = norm(cellp-ppos); + gfxReal isoadd = 0.; + const gfxReal baseIsoVal = mIsoValue*1.1; + if(len1.) { continue; } - subdAr[spk+swk][spj+swj][spi+swi] = arval + isoadd; - } } } + const gfxReal arval = subdAr[spk+swk][spj+swj][spi+swi]; + if(arval>1.) { continue; } + subdAr[spk+swk][spj+swj][spi+swi] = arval + isoadd; + } } } - p = p->getNext(); - } - } } } // poDist loops */ + p = p->getNext(); + } + } } } // poDist loops */ - py = mStart[1]+(((double)j-0.5)*orgGsy)-gsy; - for(int sj=0;sj 0) { + if (mcEdgeTable[cubeIndex] > 0) { // where to look up if this point already exists - const int edgek = 0; - const int baseIn = EDGEAR_INDEX( i+0, j+0, edgek+0, si,sj); - eVert[ 0] = &mpEdgeVerticesX[ baseIn ]; - eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ]; - eVert[ 2] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ]; - eVert[ 3] = &mpEdgeVerticesY[ baseIn ]; + const int edgek = 0; + const int baseIn = EDGEAR_INDEX( i+0, j+0, edgek+0, si,sj); + eVert[ 0] = &mpEdgeVerticesX[ baseIn ]; + eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ]; + eVert[ 2] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ]; + eVert[ 3] = &mpEdgeVerticesY[ baseIn ]; - eVert[ 4] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ]; - eVert[ 5] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+1,sj+0) ]; // with subdivs - eVert[ 6] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+1) ]; - eVert[ 7] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ]; + eVert[ 4] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ]; + eVert[ 5] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+1,sj+0) ]; // with subdivs + eVert[ 6] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+1) ]; + eVert[ 7] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ]; - eVert[ 8] = &mpEdgeVerticesZ[ baseIn ]; - eVert[ 9] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+0) ]; // with subdivs - eVert[10] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+1) ]; - eVert[11] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ]; + eVert[ 8] = &mpEdgeVerticesZ[ baseIn ]; + eVert[ 9] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+0) ]; // with subdivs + eVert[10] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+1) ]; + eVert[11] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ]; // grid positions - pos[0] = ntlVec3Gfx(px ,py ,pz); - pos[1] = ntlVec3Gfx(px+gsx,py ,pz); - pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz); // with subdivs - pos[3] = ntlVec3Gfx(px ,py+gsy,pz); - pos[4] = ntlVec3Gfx(px ,py ,pz+gsz); - pos[5] = ntlVec3Gfx(px+gsx,py ,pz+gsz); - pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz); // with subdivs - pos[7] = ntlVec3Gfx(px ,py+gsy,pz+gsz); + pos[0] = ntlVec3Gfx(px ,py ,pz); + pos[1] = ntlVec3Gfx(px+gsx,py ,pz); + pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz); // with subdivs + pos[3] = ntlVec3Gfx(px ,py+gsy,pz); + pos[4] = ntlVec3Gfx(px ,py ,pz+gsz); + pos[5] = ntlVec3Gfx(px+gsx,py ,pz+gsz); + pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz); // with subdivs + pos[7] = ntlVec3Gfx(px ,py+gsy,pz+gsz); // check all edges - for(int e=0;e<12;e++) { - if (mcEdgeTable[cubeIndex] & (1<0) { + float smoSubdfac = 1.; + if(mSubdivs>0) { //smoSubdfac = 1./(float)(mSubdivs); - smoSubdfac = pow(0.55,(double)mSubdivs); // slightly stronger - } - if(mSmoothSurface>0. || mSmoothNormals>0.) debMsgStd("IsoSurface::triangulate",DM_MSG,"Smoothing...",10); - if(mSmoothSurface>0.0) { - smoothSurface(mSmoothSurface*smoSubdfac, (mSmoothNormals<=0.0) ); - } - if(mSmoothNormals>0.0) { - smoothNormals(mSmoothNormals*smoSubdfac); - } + smoSubdfac = pow(0.55,(double)mSubdivs); // slightly stronger + } + if(mSmoothSurface>0. || mSmoothNormals>0.) debMsgStd("IsoSurface::triangulate",DM_MSG,"Smoothing...",10); + if(mSmoothSurface>0.0) { + smoothSurface(mSmoothSurface*smoSubdfac, (mSmoothNormals<=0.0) ); + } + if(mSmoothNormals>0.0) { + smoothNormals(mSmoothNormals*smoSubdfac); + } - myTime_t tritimeend = getTime(); - debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< getTimeString(tritimeend-tritimestart)<<", S("<getNumParticles(), 10); + myTime_t tritimeend = getTime(); + debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< getTimeString(tritimeend-tritimestart)<<", S("<getNumParticles(), 10); } @@ -665,8 +677,8 @@ void IsoSurface::triangulate( void ) * Get triangles for rendering *****************************************************************************/ void IsoSurface::getTriangles(double t, vector *triangles, - vector *vertices, - vector *normals, int objectId ) + vector *vertices, + vector *normals, int objectId ) { if(!mInitDone) { debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10); @@ -675,8 +687,8 @@ void IsoSurface::getTriangles(double t, vector *triangles, t = 0.; //return; // DEBUG - /* triangulate field */ - triangulate(); + /* triangulate field */ + triangulate(); //errMsg("TRIS"," "< *triangles, //errMsg("NM"," ivi"<size()<<" ns"<size()<<" ts"<size() ); //errMsg("NM"," ovs"<push_back( mPoints[i].v ); } - for(int i=0;i<(int)mPoints.size();i++) { + for(int i=0;i<(int)mPoints.size();i++) { normals->push_back( mPoints[i].n ); } //errMsg("N2"," ivi"<size()<<" ns"<size()<<" ts"<size() ); //errMsg("N2"," ovs"< *triangles, if(getCastShadows() ) { flag |= TRI_CASTSHADOWS; } - /* init geo init id */ - int geoiId = getGeoInitId(); - if(geoiId > 0) { - flag |= (1<< (geoiId+4)); - flag |= mGeoInitType; - } + /* init geo init id */ + int geoiId = getGeoInitId(); + if(geoiId > 0) { + flag |= (1<< (geoiId+4)); + flag |= mGeoInitType; + } - tri.setFlags( flag ); + tri.setFlags( flag ); - /* triangle normal missing */ - tri.setNormal( ntlVec3Gfx(0.0) ); - tri.setSmoothNormals( smooth ); - tri.setObjectId( objectId ); - triangles->push_back( tri ); + /* triangle normal missing */ + tri.setNormal( ntlVec3Gfx(0.0) ); + tri.setSmoothNormals( smooth ); + tri.setObjectId( objectId ); + triangles->push_back( tri ); } //errMsg("N3"," ivi"<size()<<" ns"<size()<<" ts"<size() ); return; @@ -743,11 +755,11 @@ inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) { // WARNING - this requires a security boundary layer... ntlVec3Gfx ret(0.0); ret[0] = *getData(i-1,j ,k ) - - *getData(i+1,j ,k ); + *getData(i+1,j ,k ); ret[1] = *getData(i ,j-1,k ) - - *getData(i ,j+1,k ); + *getData(i ,j+1,k ); ret[2] = *getData(i ,j ,k-1 ) - - *getData(i ,j ,k+1 ); + *getData(i ,j ,k+1 ); return ret; } @@ -769,13 +781,13 @@ void IsoSurface::setSmoothRad(float radi1, float radi2, ntlVec3Gfx mscc) { // compute normals for all generated triangles void IsoSurface::computeNormals() { - for(int i=0;i<(int)mPoints.size();i++) { + for(int i=0;i<(int)mPoints.size();i++) { mPoints[i].n = ntlVec3Gfx(0.); } - for(int i=0;i<(int)mIndices.size();i+=3) { - const int t1 = mIndices[i]; - const int t2 = mIndices[i+1]; + for(int i=0;i<(int)mIndices.size();i+=3) { + const int t1 = mIndices[i]; + const int t2 = mIndices[i+1]; const int t3 = mIndices[i+2]; const ntlVec3Gfx p1 = mPoints[t1].v; const ntlVec3Gfx p2 = mPoints[t2].v; @@ -792,7 +804,7 @@ void IsoSurface::computeNormals() { mPoints[t3].n += norm * (1./(len2*len3)); } - for(int i=0;i<(int)mPoints.size();i++) { + for(int i=0;i<(int)mPoints.size();i++) { normalize(mPoints[i].n); } } @@ -903,86 +915,86 @@ void IsoSurface::smoothSurface(float sigma, bool normSmooth) // Edges ntlVec3Gfx e[3] = { mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v, - mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v, - mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v }; + mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v, + mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v }; // Compute corner weights - float area = 0.5f * norm( cross(e[0], e[1])); - float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) }; - float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]), - l2[1] * (l2[2] + l2[0] - l2[1]), - l2[2] * (l2[0] + l2[1] - l2[2]) }; - if (ew[0] <= 0.0f) { - cornerareas[i][1] = -0.25f * l2[2] * area / - dot(e[0] , e[2]); - cornerareas[i][2] = -0.25f * l2[1] * area / - dot(e[0] , e[1]); - cornerareas[i][0] = area - cornerareas[i][1] - - cornerareas[i][2]; - } else if (ew[1] <= 0.0f) { - cornerareas[i][2] = -0.25f * l2[0] * area / - dot(e[1] , e[0]); - cornerareas[i][0] = -0.25f * l2[2] * area / - dot(e[1] , e[2]); - cornerareas[i][1] = area - cornerareas[i][2] - - cornerareas[i][0]; - } else if (ew[2] <= 0.0f) { - cornerareas[i][0] = -0.25f * l2[1] * area / - dot(e[2] , e[1]); - cornerareas[i][1] = -0.25f * l2[0] * area / - dot(e[2] , e[0]); - cornerareas[i][2] = area - cornerareas[i][0] - - cornerareas[i][1]; - } else { - float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]); - for (int j = 0; j < 3; j++) - cornerareas[i][j] = ewscale * (ew[(j+1)%3] + - ew[(j+2)%3]); - } + float area = 0.5f * norm( cross(e[0], e[1])); + float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) }; + float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]), + l2[1] * (l2[2] + l2[0] - l2[1]), + l2[2] * (l2[0] + l2[1] - l2[2]) }; + if (ew[0] <= 0.0f) { + cornerareas[i][1] = -0.25f * l2[2] * area / + dot(e[0] , e[2]); + cornerareas[i][2] = -0.25f * l2[1] * area / + dot(e[0] , e[1]); + cornerareas[i][0] = area - cornerareas[i][1] - + cornerareas[i][2]; + } else if (ew[1] <= 0.0f) { + cornerareas[i][2] = -0.25f * l2[0] * area / + dot(e[1] , e[0]); + cornerareas[i][0] = -0.25f * l2[2] * area / + dot(e[1] , e[2]); + cornerareas[i][1] = area - cornerareas[i][2] - + cornerareas[i][0]; + } else if (ew[2] <= 0.0f) { + cornerareas[i][0] = -0.25f * l2[1] * area / + dot(e[2] , e[1]); + cornerareas[i][1] = -0.25f * l2[0] * area / + dot(e[2] , e[0]); + cornerareas[i][2] = area - cornerareas[i][0] - + cornerareas[i][1]; + } else { + float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]); + for (int j = 0; j < 3; j++) + cornerareas[i][j] = ewscale * (ew[(j+1)%3] + + ew[(j+2)%3]); + } // NT important, check this... #ifndef WIN32 - if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6; - if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6; - if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6; + if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6; + if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6; + if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6; #else // WIN32 // FIXME check as well... - if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6; - if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6; - if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6; + if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6; + if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6; + if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6; #endif // WIN32 - pointareas[mIndices[i*3+0]] += cornerareas[i][0]; - pointareas[mIndices[i*3+1]] += cornerareas[i][1]; - pointareas[mIndices[i*3+2]] += cornerareas[i][2]; + pointareas[mIndices[i*3+0]] += cornerareas[i][0]; + pointareas[mIndices[i*3+1]] += cornerareas[i][1]; + pointareas[mIndices[i*3+2]] += cornerareas[i][2]; } } // need pointarea - // */ + // */ float invsigma2 = 1.0f / (sigma*sigma); vector dflt(nv); for (int i = 0; i < nv; i++) { if(diffuseVertexField( &mPoints[0].v, 2, - i, invsigma2, dflt[i])) { + i, invsigma2, dflt[i])) { // Just keep the displacement - dflt[i] -= mPoints[i].v; - } else { dflt[i] = 0.0; } //?mPoints[i].v; } + dflt[i] -= mPoints[i].v; + } else { dflt[i] = 0.0; } //?mPoints[i].v; } } // Slightly better small-neighborhood approximation for (int i = 0; i < nf; i++) { ntlVec3Gfx c = mPoints[mIndices[i*3+0]].v + - mPoints[mIndices[i*3+1]].v + - mPoints[mIndices[i*3+2]].v; + mPoints[mIndices[i*3+1]].v + + mPoints[mIndices[i*3+2]].v; c /= 3.0f; for (int j = 0; j < 3; j++) { int v = mIndices[i*3+j]; ntlVec3Gfx d =(c - mPoints[v].v) * 0.5f; dflt[v] += d * (cornerareas[i][j] / - pointareas[mIndices[i*3+j]] * - exp(-0.5f * invsigma2 * normNoSqrt(d)) ); + pointareas[mIndices[i*3+j]] * + exp(-0.5f * invsigma2 * normNoSqrt(d)) ); } } @@ -990,7 +1002,7 @@ void IsoSurface::smoothSurface(float sigma, bool normSmooth) vector dflt2(nv); for (int i = 0; i < nv; i++) { if(diffuseVertexField( &dflt[0], 1, - i, invsigma2, dflt2[i])) { } + i, invsigma2, dflt2[i])) { } else { /*mPoints[i].v=0.0;*/ dflt2[i] = 0.0; }//dflt2[i]; } } @@ -1062,58 +1074,58 @@ void IsoSurface::smoothNormals(float sigma) { // Edges ntlVec3Gfx e[3] = { mPoints[mIndices[i*3+2]].v - mPoints[mIndices[i*3+1]].v, - mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v, - mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v }; + mPoints[mIndices[i*3+0]].v - mPoints[mIndices[i*3+2]].v, + mPoints[mIndices[i*3+1]].v - mPoints[mIndices[i*3+0]].v }; // Compute corner weights - float area = 0.5f * norm( cross(e[0], e[1])); - float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) }; - float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]), - l2[1] * (l2[2] + l2[0] - l2[1]), - l2[2] * (l2[0] + l2[1] - l2[2]) }; - if (ew[0] <= 0.0f) { - cornerareas[i][1] = -0.25f * l2[2] * area / - dot(e[0] , e[2]); - cornerareas[i][2] = -0.25f * l2[1] * area / - dot(e[0] , e[1]); - cornerareas[i][0] = area - cornerareas[i][1] - - cornerareas[i][2]; - } else if (ew[1] <= 0.0f) { - cornerareas[i][2] = -0.25f * l2[0] * area / - dot(e[1] , e[0]); - cornerareas[i][0] = -0.25f * l2[2] * area / - dot(e[1] , e[2]); - cornerareas[i][1] = area - cornerareas[i][2] - - cornerareas[i][0]; - } else if (ew[2] <= 0.0f) { - cornerareas[i][0] = -0.25f * l2[1] * area / - dot(e[2] , e[1]); - cornerareas[i][1] = -0.25f * l2[0] * area / - dot(e[2] , e[0]); - cornerareas[i][2] = area - cornerareas[i][0] - - cornerareas[i][1]; - } else { - float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]); - for (int j = 0; j < 3; j++) - cornerareas[i][j] = ewscale * (ew[(j+1)%3] + - ew[(j+2)%3]); - } + float area = 0.5f * norm( cross(e[0], e[1])); + float l2[3] = { normNoSqrt(e[0]), normNoSqrt(e[1]), normNoSqrt(e[2]) }; + float ew[3] = { l2[0] * (l2[1] + l2[2] - l2[0]), + l2[1] * (l2[2] + l2[0] - l2[1]), + l2[2] * (l2[0] + l2[1] - l2[2]) }; + if (ew[0] <= 0.0f) { + cornerareas[i][1] = -0.25f * l2[2] * area / + dot(e[0] , e[2]); + cornerareas[i][2] = -0.25f * l2[1] * area / + dot(e[0] , e[1]); + cornerareas[i][0] = area - cornerareas[i][1] - + cornerareas[i][2]; + } else if (ew[1] <= 0.0f) { + cornerareas[i][2] = -0.25f * l2[0] * area / + dot(e[1] , e[0]); + cornerareas[i][0] = -0.25f * l2[2] * area / + dot(e[1] , e[2]); + cornerareas[i][1] = area - cornerareas[i][2] - + cornerareas[i][0]; + } else if (ew[2] <= 0.0f) { + cornerareas[i][0] = -0.25f * l2[1] * area / + dot(e[2] , e[1]); + cornerareas[i][1] = -0.25f * l2[0] * area / + dot(e[2] , e[0]); + cornerareas[i][2] = area - cornerareas[i][0] - + cornerareas[i][1]; + } else { + float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]); + for (int j = 0; j < 3; j++) + cornerareas[i][j] = ewscale * (ew[(j+1)%3] + + ew[(j+2)%3]); + } // NT important, check this... #ifndef WIN32 - if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6; - if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6; - if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6; + if(! finite(cornerareas[i][0]) ) cornerareas[i][0]=1e-6; + if(! finite(cornerareas[i][1]) ) cornerareas[i][1]=1e-6; + if(! finite(cornerareas[i][2]) ) cornerareas[i][2]=1e-6; #else // WIN32 // FIXME check as well... - if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6; - if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6; - if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6; + if(! (cornerareas[i][0]>=0.0) ) cornerareas[i][0]=1e-6; + if(! (cornerareas[i][1]>=0.0) ) cornerareas[i][1]=1e-6; + if(! (cornerareas[i][2]>=0.0) ) cornerareas[i][2]=1e-6; #endif // WIN32 - pointareas[mIndices[i*3+0]] += cornerareas[i][0]; - pointareas[mIndices[i*3+1]] += cornerareas[i][1]; - pointareas[mIndices[i*3+2]] += cornerareas[i][2]; + pointareas[mIndices[i*3+0]] += cornerareas[i][0]; + pointareas[mIndices[i*3+1]] += cornerareas[i][1]; + pointareas[mIndices[i*3+2]] += cornerareas[i][2]; } } // need pointarea