This tutorial demonstrates how to use the marching cubes algorithm to extract isosurfaces from volume data.
class App_MarchingCubes: public BaseDemo
{
public:
App_MarchingCubes(): mTest(4) {}
void initEvent()
{
openglContext()->setContinuousUpdate(true);
srand((
unsigned int)time(
NULL));
mThreshold = 0.44f;
mFountainSpeed = 2.5f;
mMetaball.resize(mParticleCount);
mMetaballVelocity.resize(mParticleCount);
mMetaballsFrames.resize(mFrameCount);
for(int iframe=0; iframe<mFrameCount; ++iframe)
{
mMetaballsFrames[iframe].resize(mParticleCount);
for(int iball=0; iball<mParticleCount; ++iball)
{
mMetaballsFrames[iframe][iball].x() = (rand()%100 / 100.0f)*0.6f+0.2f;
mMetaballsFrames[iframe][iball].y() = (rand()%100 / 100.0f)*0.6f+0.2f;
mMetaballsFrames[iframe][iball].z() = (rand()%100 / 100.0f)*0.6f+0.2f;
}
}
rendering()->
as<
vl::Rendering>()->transform()->addChild( mTransform.get() );
mText->setFont(
vl::defFontManager()->acquireFont(
"/font/bitstream-vera/VeraMono.ttf", 10) );
mText->setMargin(5);
mText->setColor(vl::white);
mText->setBackgroundColor(
vl::fvec4(0,0,0,.75f));
mText->setBackgroundEnabled(true);
mText->setText("Marching Cubes Demo");
sceneManager()->tree()->addActor( mTextActor.get() );
setupTest();
}
virtual void updateScene()
{
if (mTest == 4)
runTest4();
else
if (mTest == 5)
runTest5();
}
void runTest4()
{
return;
int frame1 = int(t * mFrameCount);
int frame2 = (frame1+1) % mFrameCount;
float Ha = 2*t*t*t - 3*t*t + 1;
float Hb = -2*t*t*t + 3*t*t;
for(unsigned iball=0; iball<mMetaball.size(); ++iball)
mMetaball[iball] = mMetaballsFrames[frame1][iball]*Ha + mMetaballsFrames[frame2][iball]*Hb;
int idx = 0;
float* values = &mMarchingCubes.volumeInfo()->at(0)->volume()->value(0);
for(int z=0; z<mMetaballsResolution; ++z)
{
float pz = (float)z/mMetaballsResolution;
for(int y=0; y<mMetaballsResolution; ++y)
{
float py = (float)y/mMetaballsResolution;
for(int x=0; x<mMetaballsResolution; ++x, ++idx)
{
values[idx] = metaballFunction((float)x/mMetaballsResolution, py, pz);
}
}
}
mMarchingCubes.volumeInfo()->at(0)->volume()->setDataDirty();
mMarchingCubes.run(false);
}
void runTest5()
{
for(unsigned iball=0; iball<mMetaball.size(); ++iball)
{
mMetaballVelocity[iball] -=
vl::fvec3(0,1.5f,0)*(float)mTimer.elapsed()*mFountainSpeed;
mMetaball[iball] = mMetaball[iball] + mMetaballVelocity[iball]*(float)mTimer.elapsed();
if (mMetaball[iball].y() < 0)
{
mMetaball[iball].
y() = 0;
mMetaball[iball].x() = 0.5f;
mMetaball[iball].z() = 0.5f;
mMetaballVelocity[iball].x() = (float)
vl::random(-12,+12);
mMetaballVelocity[iball].z() = (float)
vl::random(-12,+12);
mMetaballVelocity[iball].y() = 100;
mMetaballVelocity[iball].normalize();
mMetaballVelocity[iball] *= (17.0f + (float)
vl::random(0,5))*0.05f*mFountainSpeed;
}
}
mTimer.start();
int idx = 0;
float* values = &mMarchingCubes.volumeInfo()->at(0)->volume()->value(0);
for(int z=0; z<mMetaballsResolution; ++z)
{
float pz = (float)z/mMetaballsResolution;
for(int y=0; y<mMetaballsResolution; ++y)
{
float py = (float)y/mMetaballsResolution;
for(int x=0; x<mMetaballsResolution; ++x, ++idx)
values[idx] = metaballFunction((float)x/mMetaballsResolution, py, pz);
}
}
mMarchingCubes.volumeInfo()->at(0)->volume()->setDataDirty();
mMarchingCubes.run(false);
}
float metaballFunction(float x, float y, float z)
{
float val = 0;
for(unsigned i=0; i<mMetaball.size(); ++i)
{
float rx = x-mMetaball[i].x();
float ry = y-mMetaball[i].y();
float rz = z-mMetaball[i].z();
float r2 = rx*rx+ry*ry+rz*rz;
if (r2 == 0.0f)
return 1.0e+10f;
}
return val;
}
void keyPressEvent(
unsigned short ch,
vl::EKey key)
{
BaseDemo::keyPressEvent(ch,key);
bool update = false;
{
mTest++;
update = true;
}
else
{
mTest--;
update = true;
}
if (update)
{
if (mTest > 6) mTest = 0;
if (mTest < 0) mTest = 6;
setupTest();
updateText();
}
}
{
geom->
drawCalls().push_back(mMarchingCubes.mDrawElements.get());
sceneManager()->tree()->actors()->clear();
sceneManager()->tree()->addActor(act.
get());
sceneManager()->tree()->addActor( mTextActor.get() );
volume->
setup( (
float*)mVolumeImage->pixels(),
false,
true,
vl::fvec3(-5,-5,-5),
vl::fvec3(+5,+5,+5),
vl::ivec3(mVolumeImage->width(), mVolumeImage->height(), mVolumeImage->depth()) );
mMarchingCubes.reset();
if (test1)
{
}
mMarchingCubes.run(true);
stats +=
vl::Say(
"vertices = %n\n") << mMarchingCubes.mVertsArray->size();
}
{
{
return;
}
#if 1
#endif
#if 0
for(
int i=0; i<tex_array->
size(); ++i)
#endif
}
{
for(size_t i=0; i<color_array->size(); ++i)
{
px.
x() *= mColorImage->width();
px.
y() *= mColorImage->height();
px.
z() *= mColorImage->depth();
color_array->at(i) = mColorImage->sampleLinear(px.
x(), px.
y(), px.
z());
}
}
{
sceneManager()->tree()->actors()->clear();
sceneManager()->tree()->addActor( mTextActor.get() );
mMarchingCubes.reset();
mMarchingCubes.volumeInfo()->at(0)->volume()->
setup(
NULL,
false,
false,
vl::fvec3(-10,-10,-10),
vl::fvec3(+10,+10,+10),
vl::ivec3(mMetaballsResolution,mMetaballsResolution,mMetaballsResolution) );
geom->
drawCalls().push_back(mMarchingCubes.mDrawElements.get());
{
}
}
{
sceneManager()->tree()->actors()->clear();
sceneManager()->tree()->addActor( mTextActor.get() );
mMarchingCubes.reset();
mMarchingCubes.volumeInfo()->at(0)->volume()->
setup(
NULL,
false,
false,
vl::fvec3(-10,-10,-10),
vl::fvec3(+10,+10,+10),
vl::ivec3(mMetaballsResolution,mMetaballsResolution,mMetaballsResolution) );
geom->
drawCalls().push_back(mMarchingCubes.mDrawElements.get());
for(unsigned iball=0; iball<mMetaball.size(); ++iball)
{
mMetaball[iball].y() = 0;
mMetaball[iball].x() = 0.5f;
mMetaball[iball].z() = 0.5f;
mMetaballVelocity[iball].x() = (float)
vl::random(-15,+15);
mMetaballVelocity[iball].z() = (float)
vl::random(-15,+15);
mMetaballVelocity[iball].y() = 100;
mMetaballVelocity[iball].normalize();
mMetaballVelocity[iball] *= (5.0f + (float)
vl::random(0,20))*0.05f*mFountainSpeed;
}
mTimer.start();
}
{
public:
virtual float operator()(float x, float y, float z) const
{
return -x/5.0f*
sin(z/5.0f)+
exp(y*y*y/5.0f/5.0f/5.0f);
}
};
void setup3Dplot()
{
sceneManager()->tree()->actors()->clear();
sceneManager()->tree()->addActor( mTextActor.get() );
float range = 5.0f;
}
{
sceneManager()->tree()->actors()->clear();
sceneManager()->tree()->addActor( mTextActor.get() );
volume->
setup( (
float*)mDropImage->pixels(),
true,
false,
vl::fvec3(-5,-5,-5),
vl::fvec3(+5,+5,+5),
vl::ivec3(mDropImage->width(), mDropImage->height(), mDropImage->depth()) );
#if 0
#endif
mMarchingCubes.reset();
mMarchingCubes.volumeInfo()->push_back(
new vl::VolumeInfo(volume.
get(), mThreshold) );
mMarchingCubes.run(false);
mIsosurfGeom->setNormalArray(mMarchingCubes.mNormsArray.get());
mIsosurfGeom->drawCalls().push_back(mMarchingCubes.mDrawElements.get());
#if 0
ps.
simplify( 0.1f, mIsosurfGeom.get() );
mIsosurfGeom->computeNormals();
#endif
fx->shader()->gocLightModel()->setTwoSide(true);
fx->shader()->gocMaterial()->setBackDiffuse(vl::green);
#if defined(VL_OPENGL)
fx->shader(0,1)->gocPolygonOffset()->set(-1.0f, -1.0f);
fx->shader(0,1)->gocColor()->setValue(vl::royalblue);
#endif
sceneManager()->tree()->addActor( actor.
get() );
updateText();
}
void setupTest()
{
sceneManager()->tree()->eraseAllChildren();
if (mTest == 1)
showVolumes(true);
else
if (mTest == 2)
textureVolumeColor( showVolumes(false) );
else
if (mTest == 3)
vertexVolumeColor( showVolumes(false) );
else
if (mTest == 4)
setupMetaballs();
else
if (mTest == 5)
setupFountain();
else
if (mTest == 6)
setup3Dplot();
updateText();
}
void updateText()
{
if(mTest == 0)
str =
vl::Say(
"Marching Cubes Test #0 - Volume Viewer - threshold = %.2n\n(drop a .dat file)") << mThreshold;
else
if(mTest == 1)
str = "Marching Cubes Test #1 - Multiple Volumes & Transparency";
else
if(mTest == 2)
str = "Marching Cubes Test #2 - Texture Colorized Volumes";
else
if(mTest == 3)
str = "Marching Cubes Test #3 - Vertex Colorized Volumes";
else
if(mTest == 4)
str = "Marching Cubes Test #4 - Metaballs";
else
if(mTest == 5)
str = "Marching Cubes Test #5 - Fountain";
else
if(mTest == 6)
str = "Marching Cubes Test #6 - 3D Function Plotting";
str += "\n(press the <- or -> key to change test)";
mText->setText( str );
mText->setDisplayListDirty(true);
}
void mouseWheelEvent(int w)
{
if (mTest != 0)
return;
if (w>0)
mThreshold += 0.010f;
else
mThreshold -= 0.010f;
mThreshold =
vl::clamp(mThreshold, 0.0f, 1.0f);
mMarchingCubes.volumeInfo()->at(0)->setThreshold(mThreshold);
mMarchingCubes.run(false);
if (mIsosurfGeom)
mIsosurfGeom->setBufferObjectDirty(true);
updateText();
openglContext()->update();
}
void fileDroppedEvent(const std::vector<vl::String>& files)
{
mTest = 0;
if(files.size() == 1)
{
if (files[0].endsWith(".dat"))
{
loadVolume(vol_img);
}
}
else
{
std::vector<vl::String> files_sorted = files;
std::sort(files_sorted.begin(), files_sorted.end());
std::vector< vl::ref<vl::Image> > images;
for(unsigned int i=0; i<files_sorted.size(); ++i)
if (vol_img)
loadVolume(vol_img);
}
}
void resizeEvent(int w, int h)
{
BaseDemo::resizeEvent(w,h);
mText->setDisplayListDirty(true);
}
protected:
float mThreshold;
std::vector<vl::fvec3> mMetaball;
std::vector<vl::fvec3> mMetaballVelocity;
std::vector< std::vector<vl::fvec3> > mMetaballsFrames;
static const int mMetaballsResolution = 32;
static const int mParticleCount = 25;
static const int mFrameCount = 20;
float mFountainSpeed;
int mTest;
};