griefith/test2
1
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases Wiki Activity

Code got unreadble due to copy-paste (hint for me: revert point)

Browse Source
This commit is contained in:
Daniel Genrich
2007-11-26 14:59:58 +00:00
parent f50c70bc7d
commit db57ad68cd
1 changed files with 541 additions and 529 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1070
intern/elbeem/intern/isosurface.cpp
View File

@@ -18,28 +18,32 @@
#define round(x) (x)
#endif
#if PARALLEL==1
#include <omp.h>
#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;i<nodes;i++) { mpData[i] = 0.0; }
mpData = new float[nodes];
for(int i=0;i<nodes;i++) { mpData[i] = 0.0; }
// allocate edge arrays (last slices are never used...)
int initsize = -1;
@@ -93,7 +97,7 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e
mpEdgeVerticesZ = new int[sliceNodes];
initsize = 3*sliceNodes;
}
for(int i=0;i<mEdgeArSize;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
for(int i=0;i<mEdgeArSize;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
// WARNING - make sure this is consistent with calculateMemreqEstimate
// marching cubes are ready
@@ -106,7 +110,7 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e
/*! Reset all values */
void IsoSurface::resetAll(gfxReal val) {
int nodes = mSizez*mSizey*mSizex;
for(int i=0;i<nodes;i++) { mpData[i] = val; }
for(int i=0;i<nodes;i++) { mpData[i] = val; }
}
@@ -130,9 +134,9 @@ IsoSurface::~IsoSurface( void )
*****************************************************************************/
void IsoSurface::triangulate( void )
{
double gsx,gsy,gsz; // grid spacing in x,y,z direction
double px,py,pz; // current position in grid in x,y,z direction
// IsoLevelCube cubie; // struct for a small subcube
double gsx,gsy,gsz; // grid spacing in x,y,z direction
double px,py,pz; // current position in grid in x,y,z direction
IsoLevelCube cubie; // struct for a small subcube
myTime_t tritimestart = getTime();
if(!mpData) {
@@ -141,16 +145,16 @@ void IsoSurface::triangulate( void )
}
// get grid spacing (-2 to have same spacing as sim)
gsx = (mEnd[0]-mStart[0])/(double)(mSizex-2.0);
gsy = (mEnd[1]-mStart[1])/(double)(mSizey-2.0);
gsz = (mEnd[2]-mStart[2])/(double)(mSizez-2.0);
gsx = (mEnd[0]-mStart[0])/(double)(mSizex-2.0);
gsy = (mEnd[1]-mStart[1])/(double)(mSizey-2.0);
gsz = (mEnd[2]-mStart[2])/(double)(mSizez-2.0);
// clean up previous frame
mIndices.clear();
mPoints.clear();
// reset edge vertices
for(int i=0;i<mEdgeArSize;i++) {
for(int i=0;i<mEdgeArSize;i++) {
mpEdgeVerticesX[i] = -1;
mpEdgeVerticesY[i] = -1;
mpEdgeVerticesZ[i] = -1;
@@ -159,418 +163,421 @@ void IsoSurface::triangulate( void )
// edges between which points?
const int mcEdges[24] = {
0,1, 1,2, 2,3, 3,0,
4,5, 5,6, 6,7, 7,4,
0,4, 1,5, 2,6, 3,7 };
4,5, 5,6, 6,7, 7,4,
0,4, 1,5, 2,6, 3,7 };
const int cubieOffsetX[8] = {
0,1,1,0, 0,1,1,0 };
const int cubieOffsetY[8] = {
0,0,1,1, 0,0,1,1 };
const int cubieOffsetZ[8] = {
0,0,0,0, 1,1,1,1 };
const int cubieOffsetX[8] = {
0,1,1,0, 0,1,1,0 };
const int cubieOffsetY[8] = {
0,0,1,1, 0,0,1,1 };
const int cubieOffsetZ[8] = {
0,0,0,0, 1,1,1,1 };
const int coAdd=2;
const int coAdd=2;
// let the cubes march
if(mSubdivs<=1) {
if(mSubdivs<=1) {
pz = mStart[2]-gsz*0.5;
for(int k=1;k<(mSizez-2);k++) {
pz += gsz;
py = mStart[1]-gsy*0.5;
for(int j=1;j<(mSizey-2);j++) {
py += gsy;
px = mStart[0]-gsx*0.5;
for(int i=1;i<(mSizex-2);i++) {
px += gsx;
int cubeIndex; // index entry of the cube
float value[8];
int triIndices[12]; // vertex indices
int *eVert[12];
IsoLevelVertex ilv;
pz = mStart[2]-gsz*0.5;
for(int k=1;k<(mSizez-2);k++) {
pz += gsz;
py = mStart[1]-gsy*0.5;
for(int j=1;j<(mSizey-2);j++) {
py += gsy;
px = mStart[0]-gsx*0.5;
for(int i=1;i<(mSizex-2);i++) {
px += gsx;
int cubeIndex; // index entry of the cube
float value[8];
int triIndices[12]; // vertex indices
int *eVert[12];
IsoLevelVertex ilv;
value[0] = *getData(i ,j ,k );
value[1] = *getData(i+1,j ,k );
value[2] = *getData(i+1,j+1,k );
value[3] = *getData(i ,j+1,k );
value[4] = *getData(i ,j ,k+1);
value[5] = *getData(i+1,j ,k+1);
value[6] = *getData(i+1,j+1,k+1);
value[7] = *getData(i ,j+1,k+1);
value[0] = *getData(i ,j ,k );
value[1] = *getData(i+1,j ,k );
value[2] = *getData(i+1,j+1,k );
value[3] = *getData(i ,j+1,k );
value[4] = *getData(i ,j ,k+1);
value[5] = *getData(i+1,j ,k+1);
value[6] = *getData(i+1,j+1,k+1);
value[7] = *getData(i ,j+1,k+1);
// check intersections of isosurface with edges, and calculate cubie index
cubeIndex = 0;
if (value[0] < mIsoValue) cubeIndex |= 1;
if (value[1] < mIsoValue) cubeIndex |= 2;
if (value[2] < mIsoValue) cubeIndex |= 4;
if (value[3] < mIsoValue) cubeIndex |= 8;
if (value[4] < mIsoValue) cubeIndex |= 16;
if (value[5] < mIsoValue) cubeIndex |= 32;
if (value[6] < mIsoValue) cubeIndex |= 64;
if (value[7] < mIsoValue) cubeIndex |= 128;
cubeIndex = 0;
if (value[0] < mIsoValue) cubeIndex |= 1;
if (value[1] < mIsoValue) cubeIndex |= 2;
if (value[2] < mIsoValue) cubeIndex |= 4;
if (value[3] < mIsoValue) cubeIndex |= 8;
if (value[4] < mIsoValue) cubeIndex |= 16;
if (value[5] < mIsoValue) cubeIndex |= 32;
if (value[6] < mIsoValue) cubeIndex |= 64;
if (value[7] < mIsoValue) cubeIndex |= 128;
// No triangles to generate?
if (mcEdgeTable[cubeIndex] == 0) {
continue;
}
if (mcEdgeTable[cubeIndex] == 0) {
continue;
}
// where to look up if this point already exists
int edgek = 0;
if(mUseFullEdgeArrays) edgek=k;
const int baseIn = ISOLEVEL_INDEX( i+0, j+0, edgek+0);
eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
int edgek = 0;
if(mUseFullEdgeArrays) edgek=k;
const int baseIn = ISOLEVEL_INDEX( i+0, j+0, edgek+0);
eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, edgek+1) ];
eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+1) ];
eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, edgek+1) ];
eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+1) ];
eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, edgek+0) ];
eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, edgek+0) ];
eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, edgek+0) ];
eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, edgek+0) ];
eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
// grid positions
ntlVec3Gfx pos[8];
pos[0] = ntlVec3Gfx(px ,py ,pz);
pos[1] = ntlVec3Gfx(px+gsx,py ,pz);
pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
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);
pos[7] = ntlVec3Gfx(px ,py+gsy,pz+gsz);
ntlVec3Gfx pos[8];
pos[0] = ntlVec3Gfx(px ,py ,pz);
pos[1] = ntlVec3Gfx(px+gsx,py ,pz);
pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
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);
pos[7] = ntlVec3Gfx(px ,py+gsy,pz+gsz);
// check all edges
for(int e=0;e<12;e++) {
if (mcEdgeTable[cubeIndex] & (1<<e)) {
for(int e=0;e<12;e++) {
if (mcEdgeTable[cubeIndex] & (1<<e)) {
// is the vertex already calculated?
if(*eVert[ e ] < 0) {
if(*eVert[ e ] < 0) {
// interpolate edge
const int e1 = mcEdges[e*2 ];
const int e2 = mcEdges[e*2+1];
const ntlVec3Gfx p1 = pos[ e1 ]; // scalar field pos 1
const ntlVec3Gfx p2 = pos[ e2 ]; // scalar field pos 2
const float valp1 = value[ e1 ]; // scalar field val 1
const float valp2 = value[ e2 ]; // scalar field val 2
const float mu = (mIsoValue - valp1) / (valp2 - valp1);
const int e1 = mcEdges[e*2 ];
const int e2 = mcEdges[e*2+1];
const ntlVec3Gfx p1 = pos[ e1 ]; // scalar field pos 1
const ntlVec3Gfx p2 = pos[ e2 ]; // scalar field pos 2
const float valp1 = value[ e1 ]; // scalar field val 1
const float valp2 = value[ e2 ]; // scalar field val 2
const float mu = (mIsoValue - valp1) / (valp2 - valp1);
// init isolevel vertex
ilv.v = p1 + (p2-p1)*mu;
ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) +
getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * ( mu) ;
mPoints.push_back( ilv );
ilv.v = p1 + (p2-p1)*mu;
ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) +
getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * ( mu) ;
mPoints.push_back( ilv );
triIndices[e] = (mPoints.size()-1);
triIndices[e] = (mPoints.size()-1);
// store vertex
*eVert[ e ] = triIndices[e];
} else {
*eVert[ e ] = triIndices[e];
} else {
// retrieve from vert array
triIndices[e] = *eVert[ e ];
}
} // along all edges
}
triIndices[e] = *eVert[ e ];
}
} // along all edges
}
if( (i<coAdd+mCutoff) || (j<coAdd+mCutoff) ||
((mCutoff>0) && (k<coAdd)) ||// bottom layer
(i>mSizex-2-coAdd-mCutoff) ||
(j>mSizey-2-coAdd-mCutoff) ) {
if(mCutArray) {
if(k < mCutArray[j*this->mSizex+i]) continue;
} else { continue; }
}
if( (i<coAdd+mCutoff) || (j<coAdd+mCutoff) ||
((mCutoff>0) && (k<coAdd)) ||// bottom layer
(i>mSizex-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<ParticleObject>::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<ParticleObject>::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<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
pit!= mpIsoParts->getParticlesEnd(); pit++) {
if( (*pit).getActive()==false ) continue;
if( (*pit).getType()!=PART_DROP) continue;
for(vector<ParticleObject>::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:"<<pInUse<<", Subdivs:"<<mSubdivs, 9);
pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz;
for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) {
pz += gsz;
const int k = ok/mSubdivs;
if(k<=0) continue; // skip zero plane
debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<<pInUse<<", Subdivs:"<<mSubdivs, 9);
pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz;
for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) {
pz += gsz;
const int k = ok/mSubdivs;
if(k<=0) continue; // skip zero plane
#if PARALLEL==1
#pragma omp parallel for
for(int j=1;j<(mSizey-2);j++) {
for(int i=1;i<(mSizex-2);i++) {
float value[8];
ntlVec3Gfx pos[8];
int cubeIndex; // index entry of the cube
int triIndices[12]; // vertex indices
int *eVert[12];
IsoLevelVertex ilv;
gfxReal orgval[8];
gfxReal subdAr[2][11][11]; // max 10 subdivs!
#endif
for(int j=1;j<(mSizey-2);j++) {
for(int i=1;i<(mSizex-2);i++) {
float value[8];
ntlVec3Gfx pos[8];
int cubeIndex; // index entry of the cube
int triIndices[12]; // vertex indices
int *eVert[12];
IsoLevelVertex ilv;
gfxReal orgval[8];
gfxReal subdAr[2][11][11]; // max 10 subdivs!
orgval[0] = *getData(i ,j ,k );
orgval[1] = *getData(i+1,j ,k );
orgval[2] = *getData(i+1,j+1,k ); // with subdivs
orgval[3] = *getData(i ,j+1,k );
orgval[4] = *getData(i ,j ,k+1);
orgval[5] = *getData(i+1,j ,k+1);
orgval[6] = *getData(i+1,j+1,k+1); // with subdivs
orgval[7] = *getData(i ,j+1,k+1);
orgval[0] = *getData(i ,j ,k );
orgval[1] = *getData(i+1,j ,k );
orgval[2] = *getData(i+1,j+1,k ); // with subdivs
orgval[3] = *getData(i ,j+1,k );
orgval[4] = *getData(i ,j ,k+1);
orgval[5] = *getData(i+1,j ,k+1);
orgval[6] = *getData(i+1,j+1,k+1); // with subdivs
orgval[7] = *getData(i ,j+1,k+1);
// prebuild subsampled array slice
const int sdkOffset = ok-k*mSubdivs;
const int sdkOffset = ok-k*mSubdivs;
for(int sdk=0; sdk<2; sdk++)
for(int sdj=0; sdj<mSubdivs+1; sdj++)
for(int sdi=0; sdi<mSubdivs+1; sdi++) {
subdAr[sdk][sdj][sdi] = ISOTRILININT(sdi*subdfac, sdj*subdfac, (sdkOffset+sdk)*subdfac);
}
for(int sdk=0; sdk<2; sdk++)
for(int sdj=0; sdj<mSubdivs+1; sdj++)
for(int sdi=0; sdi<mSubdivs+1; sdi++) {
subdAr[sdk][sdj][sdi] = ISOTRILININT(sdi*subdfac, sdj*subdfac, (sdkOffset+sdk)*subdfac);
}
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) { // */
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<ParticleObject>::iterator pit= mpIsoParts->getParticlesBegin();
pit!= mpIsoParts->getParticlesEnd(); pit++) { { { {
for(vector<ParticleObject>::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(pfLen<minPfLen) pfLen = minPfLen;
gfxReal pfLen = p->getSize()*1.5*mPartSize; // test, was 1.1
const gfxReal minPfLen = subdfac*0.8;
if(pfLen<minPfLen) pfLen = minPfLen;
//errMsg("ISOPPP"," at "<<PRINT_IJK<<" pp"<<ppos<<" sp"<<PRINT_VEC(spi,spj,spk)<<" pflen"<<pfLen );
//errMsg("ISOPPP"," subdfac="<<subdfac<<" size"<<p->getSize()<<" ps"<<mPartSize );
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(
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<pfLen) {
isoadd = baseIsoVal*1.;
} else {
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<pfLen) {
isoadd = baseIsoVal*1.;
} else {
// falloff linear with pfLen (kernel size=2pfLen
isoadd = baseIsoVal*(1. - (len-pfLen)/(pfLen));
}
if(isoadd<0.) { continue; }
isoadd = baseIsoVal*(1. - (len-pfLen)/(pfLen));
}
if(isoadd<0.) { continue; }
//errMsg("ISOPPP"," at "<<PRINT_IJK<<" sp"<<PRINT_VEC(spi+swi,spj+swj,spk+swk)<<" cellp"<<cellp<<" pp"<<ppos << " l"<< len<< " add"<< isoadd);
const gfxReal arval = subdAr[spk+swk][spj+swj][spi+swi];
if(arval>1.) { 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<mSubdivs;sj++) {
py += gsy;
px = mStart[0]+(((double)i-0.5)*orgGsx)-gsx;
for(int si=0;si<mSubdivs;si++) {
px += gsx;
value[0] = subdAr[0+0][sj+0][si+0];
value[1] = subdAr[0+0][sj+0][si+1];
value[2] = subdAr[0+0][sj+1][si+1];
value[3] = subdAr[0+0][sj+1][si+0];
value[4] = subdAr[0+1][sj+0][si+0];
value[5] = subdAr[0+1][sj+0][si+1];
value[6] = subdAr[0+1][sj+1][si+1];
value[7] = subdAr[0+1][sj+1][si+0];
py = mStart[1]+(((double)j-0.5)*orgGsy)-gsy;
for(int sj=0;sj<mSubdivs;sj++) {
py += gsy;
px = mStart[0]+(((double)i-0.5)*orgGsx)-gsx;
for(int si=0;si<mSubdivs;si++) {
px += gsx;
value[0] = subdAr[0+0][sj+0][si+0];
value[1] = subdAr[0+0][sj+0][si+1];
value[2] = subdAr[0+0][sj+1][si+1];
value[3] = subdAr[0+0][sj+1][si+0];
value[4] = subdAr[0+1][sj+0][si+0];
value[5] = subdAr[0+1][sj+0][si+1];
value[6] = subdAr[0+1][sj+1][si+1];
value[7] = subdAr[0+1][sj+1][si+0];
// check intersections of isosurface with edges, and calculate cubie index
cubeIndex = 0;
if (value[0] < mIsoValue) cubeIndex |= 1;
if (value[1] < mIsoValue) cubeIndex |= 2; // with subdivs
if (value[2] < mIsoValue) cubeIndex |= 4;
if (value[3] < mIsoValue) cubeIndex |= 8;
if (value[4] < mIsoValue) cubeIndex |= 16;
if (value[5] < mIsoValue) cubeIndex |= 32; // with subdivs
if (value[6] < mIsoValue) cubeIndex |= 64;
if (value[7] < mIsoValue) cubeIndex |= 128;
cubeIndex = 0;
if (value[0] < mIsoValue) cubeIndex |= 1;
if (value[1] < mIsoValue) cubeIndex |= 2; // with subdivs
if (value[2] < mIsoValue) cubeIndex |= 4;
if (value[3] < mIsoValue) cubeIndex |= 8;
if (value[4] < mIsoValue) cubeIndex |= 16;
if (value[5] < mIsoValue) cubeIndex |= 32; // with subdivs
if (value[6] < mIsoValue) cubeIndex |= 64;
if (value[7] < mIsoValue) cubeIndex |= 128;
if (mcEdgeTable[cubeIndex] > 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<<e)) {
for(int e=0;e<12;e++) {
if (mcEdgeTable[cubeIndex] & (1<<e)) {
// is the vertex already calculated?
if(*eVert[ e ] < 0) {
if(*eVert[ e ] < 0) {
// interpolate edge
const int e1 = mcEdges[e*2 ];
const int e2 = mcEdges[e*2+1];
@@ -582,79 +589,84 @@ void IsoSurface::triangulate( void )
// init isolevel vertex
ilv.v = p1 + (p2-p1)*mu; // with subdivs
#pragma omp critical
#if PARALLEL==1
#pragma omp critical
#endif
{
mPoints.push_back( ilv );
triIndices[e] = (mPoints.size()-1);
}
// store vertex
*eVert[ e ] = triIndices[e];
} else {
} else {
// retrieve from vert array
triIndices[e] = *eVert[ e ];
}
} // along all edges
}
}
} // along all edges
}
// removed cutoff treatment...
#pragma omp critical
{
// 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] ] ); // with subdivs
mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
//errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() );
}
}
} // triangles in edge table?
}//si
}// sj
}//i
}// j
// Create the triangles...
#if PARALLEL==1
#pragma omp critical
#endif
{
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] ] ); // with subdivs
mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
//errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() );
} // triangles in edge table?
}
}
}//si
}// sj
}//i
}// j
// copy edge arrays
for(int j=0;j<(mSizey-0)*mSubdivs;j++)
for(int i=0;i<(mSizex-0)*mSubdivs;i++) {
for(int j=0;j<(mSizey-0)*mSubdivs;j++)
for(int i=0;i<(mSizex-0)*mSubdivs;i++) {
//int edgek = 0;
const int dst = EDGEAR_INDEX( 0, 0, 0, i,j);
const int src = EDGEAR_INDEX( 0, 0, 1, i,j);
mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ];
mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; // with subdivs
mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ];
mpEdgeVerticesX[ src ]=-1;
mpEdgeVerticesY[ src ]=-1; // with subdivs
mpEdgeVerticesZ[ src ]=-1;
}
// */
const int dst = EDGEAR_INDEX( 0, 0, 0, i,j);
const int src = EDGEAR_INDEX( 0, 0, 1, i,j);
mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ];
mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; // with subdivs
mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ];
mpEdgeVerticesX[ src ]=-1;
mpEdgeVerticesY[ src ]=-1; // with subdivs
mpEdgeVerticesZ[ src ]=-1;
}
// */
} // ok, k subdiv loop
} // ok, k subdiv loop
//delete [] subdAr;
delete [] arppnt;
computeNormals();
} // with subdivs
delete [] arppnt;
computeNormals();
} // with subdivs
// perform smoothing
float smoSubdfac = 1.;
if(mSubdivs>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("<<mSmoothSurface<<","<<mSmoothNormals<<"),"<<
" verts:"<<mPoints.size()<<" tris:"<<(mIndices.size()/3)<<" subdivs:"<<mSubdivs
, 10 );
if(mpIsoParts) debMsgStd("IsoSurface::triangulate",DM_MSG,"parts:"<<mpIsoParts->getNumParticles(), 10);
myTime_t tritimeend = getTime();
debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< getTimeString(tritimeend-tritimestart)<<", S("<<mSmoothSurface<<","<<mSmoothNormals<<"),"<<
" verts:"<<mPoints.size()<<" tris:"<<(mIndices.size()/3)<<" subdivs:"<<mSubdivs
, 10 );
if(mpIsoParts) debMsgStd("IsoSurface::triangulate",DM_MSG,"parts:"<<mpIsoParts->getNumParticles(), 10);
}
@@ -665,8 +677,8 @@ void IsoSurface::triangulate( void )
* Get triangles for rendering
*****************************************************************************/
void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles,
vector<ntlVec3Gfx> *vertices,
vector<ntlVec3Gfx> *normals, int objectId )
vector<ntlVec3Gfx> *vertices,
vector<ntlVec3Gfx> *normals, int objectId )
{
if(!mInitDone) {
debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10);
@@ -675,8 +687,8 @@ void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles,
t = 0.;
//return; // DEBUG
/* triangulate field */
triangulate();
/* triangulate field */
triangulate();
//errMsg("TRIS"," "<<mIndices.size() );
// new output with vertice reuse
@@ -689,20 +701,20 @@ void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles,
//errMsg("NM"," ivi"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
//errMsg("NM"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
for(int i=0;i<(int)mPoints.size();i++) {
for(int i=0;i<(int)mPoints.size();i++) {
vertices->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"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->size() );
//errMsg("N2"," ovs"<<mVertices.size()<<" ons"<<mVertNormals.size()<<" ots"<<mIndices.size() );
for(int i=0;i<(int)mIndices.size();i+=3) {
for(int i=0;i<(int)mIndices.size();i+=3) {
const int smooth = 1;
int t1 = mIndices[i];
int t2 = mIndices[i+1];
int t1 = mIndices[i];
int t2 = mIndices[i+1];
int t3 = mIndices[i+2];
//errMsg("NM"," tri"<<t1<<" "<<t2<<" "<<t3 );
@@ -718,20 +730,20 @@ void IsoSurface::getTriangles(double t, vector<ntlTriangle> *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"<<iniVertIndex<<" ini"<<iniNormIndex<<" vs"<<vertices->size()<<" ns"<<normals->size()<<" ts"<<triangles->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<ntlVec3Gfx> 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<ntlVec3Gfx> 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