# Metaballs with Marching cubes

This topic is 1641 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

Hi everyone. I have implemented the marching cubes algorithm and I am trying to render metaballs.

The problem I have is no matter how many metaballs I define, it will always only draw 1, so I guess my question is how to properly integrate the metaballs data with marching cubes.

Each metaball  has a position and radius and then stored in a vector container.

Then in my volume function for marching cubes I have this:

 for(i=0; i<m_volume_depth; ++i)
{
for(j=0; j<m_volume_height; ++j)
{

for(k=0; k<m_volume_width; ++k)
{
float x = (float)i/(m_volume_depth) - 0.5;
float y = (float)k/(m_volume_width) - 0.5;
float z = (float)j/(m_volume_height) - 0.5;

Vec3 p = Vec3(x,y,z);
m_volumeData[i*m_volume_width*m_volume_height + j*m_volume_width + k] = meta(p);
}
}
}



Then here is the meta function:

 float sum=0;

for(int x=0; x<m_metaBalls.size(); ++x)
{

sum +=  m_metaBalls.at(x)->equation_sphere(_pos);
}

return sum;


Here is my equation code:

 return ( (_pos.m_x - m_position.m_x)*(_pos.m_x - m_position.m_x)
+ (_pos.m_y - m_position.m_y)*(_pos.m_y - m_position.m_y)
+ (_pos.m_z - m_position.m_z)*(_pos.m_z - m_position.m_z)


_pos is the passed in position;

m_position is the origin of the metaball;

I don't think there is anything wrong with the metaball class becuase it draws a sphere perfectly but only 1.

How would you integrate the metaball data with marching cubes?

If you require any other parts of the code then please let me know, this has been driving me crazy for a while.

##### Share on other sites

If you are not using marching cubes already how exactly are you rendering the data in m_volumeData at all?  You said it draws 1 sphere so it must be drawing something with some type of algorithm.  There is a good chance you reinvented a marching cubes type algorithm without realizing it.  Obviously there is a bug but from what I've seen so far it's likely an implementation bug not a conceptual one...you appear to have the right idea.

Marching cubes is really simply.  The basic premise is that you chop your volume up into voxles and evaluate your potential function at each corner of every voxel.  By computing where along a voxel edge you transition form being inside the isosurface (also called the level-set or isocountour) to being outside the isosurface geometry can be generated that approximates the boundary.

Here's a good place to start:  http://paulbourke.net/geometry/polygonise/ which also covers a slightly simpler and more aesthetically pleasing variant:  Marching Tetrahedrons.

Edited by nonoptimalrobot

##### Share on other sites

If you are not using marching cubes already how exactly are you rendering the data in m_volumeData at all?  You said it draws 1 sphere so it must be drawing something with some type of algorithm.  There is a good chance you reinvented a marching cubes type algorithm without realizing it.  Obviously there is a bug but from what I've seen so far it's likely an implementation bug not a conceptual one...you appear to have the right idea.

Marching cubes is really simply.  The basic premise is that you chop your volume up into voxles and evaluate your potential function at each corner of every voxel.  By computing where along a voxel edge you transition form being inside the isosurface (also called the level-set or isocountour) to being outside the isosurface geometry can be generated that approximates the boundary.

Here's a good place to start:  http://paulbourke.net/geometry/polygonise/ which also covers a slightly simpler and more aesthetically pleasing variant:  Marching Tetrahedrons.

I have marching cubes already implemented and thats how its drawing 1 sphere from the equation.

I dunno how to explain this but basically how would I pass the metaball data to draw multiple spheres instead of just 1.

The marching cubes algorithm works and my equations work.

For example if I just passed the equation of a circle to my volume data

m_volumeData[i*m_volume_width*m_volume_height + j*m_volume_width + k] = x*x + y*y +z*z - r*r

It draws a sphere perfectly, Im just unsure how would I pass multiple spheres to m_volumeData.

##### Share on other sites

I dunno how to explain this but basically how would I pass the metaball data to draw multiple spheres instead of just 1.

The marching cubes algorithm works and my equations work.

For example if I just passed the equation of a circle to my volume data

m_volumeData[i*m_volume_width*m_volume_height + j*m_volume_width + k] = x*x + y*y +z*z - r*r

You just add more point sources of potential which you appear to be doing with this bit of code:

 float sum=0;

for(int x=0; x<m_metaBalls.size(); ++x)
{

sum +=  m_metaBalls.at(x)->equation_sphere(_pos);
}

return sum;


Adding more balls at different locations should give you a blobby topology, if you have two balls at the same location you should end up with a sphere that's the sum of the radii.  Double check you meta(p) function to verify that it is iterating over all the items in m_metaBalls and returning a sum of all the potential functions.  If that's in order perhaps you can post some more code or explain what happens when you add items to m_metaBalls.

Edited by nonoptimalrobot

##### Share on other sites

Unless I've misunderstood your code, there's something wrong with your potential function. It should go to zero as you go further away from the sphere's center, and reach zero at some distance, not be negative inside the sphere and then tend to infinity, otherwise it will be really hard for a sphere to interact with the potential field of another (since the "most negative" value is negative radius squared, whereas far away from the first sphere your potential becomes arbitrarily large).

The basic "metaball" potential function is 1 / (distance to center squared). There are faster/smoother approximations out there, check out http://en.wikipedia.org/wiki/Metaball.

##### Share on other sites

According to my pen and paper Bacterius is correct.  This is actually very interesting, it never occurred to me that you could pick a simple monotonically increasing potential function that could result in such an intractable field given multiple sources.

This bit of code:

 return ( (_pos.m_x - m_position.m_x)*(_pos.m_x - m_position.m_x)
+ (_pos.m_y - m_position.m_y)*(_pos.m_y - m_position.m_y)
+ (_pos.m_z - m_position.m_z)*(_pos.m_z - m_position.m_z)


Should be something like this:

float delX  = pos.m_x - m_position.m_x;
float delY  = pos.m_y - m_position.m_y;
float delZ  = pos.m_z - m_position.m_z;


?

And you will want to tweak your Marching Cube code such that you tessellate the isocontour at meta(p) = 1 instead of meta(p) = 0.  This is a rather naïve potential and has a singularity but it will get the job done.  The Wikipedia link Bacterius specified will point you toward a number of better choices.

##### Share on other sites

Thanks guys got it working :) Although I have a couple of queries

After changing the equation, I had to reduce the iso value quite alot to 0.008 is it normal to be that low?

Also I can only specify metaball positions between 0 -> 1 anything outside doesn't get drawn and I want to be able to draw them anywhere on screen, Im not sure which part of the algorithm decides the range?

Finally just a general question, eventually I wanted to use this to render fluid. Some people seem to use marching cubes to render fluid, however It took my app like 4seconds just to render those four static metaballs. What kind of technique is good for fluid rendering?

Thanks for the help guys.

##### Share on other sites

After changing the equation, I had to reduce the iso value quite alot to 0.008 is it normal to be that low?

The isovalue really tells you where the limit between solid regions and empty space (the "contour") is to be located. If it is too high or too low, everything will be considered solid or empty, and as you vary it more (or less) of the density field is considered solid, which leads to interesting results. I don't think there is really any "good" value for an isovalue, it depends a lot on your metaball implementation and on the effect you want to achieve, so I wouldn't be worried if you need to tweak it a bit to get good-looking metaballs.

Also I can only specify metaball positions between 0 -> 1 anything outside doesn't get drawn and I want to be able to draw them anywhere on screen, Im not sure which part of the algorithm decides the range?

This shouldn't be a problem with the metaball algorithm itself. I would check your marching cubes code to verify it is working outside the 0..1 unit cube. If it isn't, then obviously everything outside it will never be polygonized, and therefore never rendered. It's not possible to polygonize things arbitrarily far away, unfortunately, there is always a tradeoff between range and resolution, because marching cubes is a discrete algorithm. Other methods like ray marching can sort of scale to arbitrary distances, but have their own set of drawbacks.

To make sure your balls can move around, you can compute the bounding box of your metaballs and use that to define the range for the marching cubes algorithm, but you can't have them go too far away from one another else you will lose in accuracy. You might wonder why, since in that case the marching cubes algorithm would spend most of its time on empty voxels, and, yes, it is possible to do better, but it gets rather nasty as you then need to find an approximate bounding box for every set of connected metaballs and run the marching cubes algorithm on each of them, separately, which sounds great on paper but isn't that efficient in practice. I think most people deal with this by making reasonable tradeoffs between the size of their worlds and how precise the polygonization should be.

Finally just a general question, eventually I wanted to use this to render fluid. Some people seem to use marching cubes to render fluid, however It took my app like 4seconds just to render those four static metaballs. What kind of technique is good for fluid rendering?

Marching cubes can be rather expensive, especially if you want really good resolution. Looking at your screenshot, your mesh is quite smooth, meaning it's probably trying to polygonize up to 128x128x128 voxels, which is *a lot* especially done on the CPU. If you wanted to go interactive - it is possible - you would move the marching cubes code to the GPU, in a shader, and perhaps scale back the resolution a notch. It isn't too hard at all, in fact, and the speed boost is huge since doing the same calculation over a lot of different locations is what the graphics card does best. Then you can add lots of algorithmic optimizations to avoid calculating every single sphere's potential for each voxel - which obviously won't do when you start rendering dozens of metaballs - by using techniques like octrees or spatial hashing. It can get rather efficient, really, but it has its limits. As we all know, the really cool stuff is done by cleverly combining different techniques

I don't know if metaballs would be my first choice for fluid rendering, though. It stills seems like an inefficient method overall, I would just use a grid-based fluid dynamics system if I wanted to do it properly (though a fractal heightmap or even FFT water works well for static water bodies). Marching cubes itself sounds good to display the results, however. Remember to dissociate what you are rendering (metaballs) from how you are rendering it (marching cubes), the two have nothing in common except being often used together.

... great, now I want to write a metaball renderer

##### Share on other sites

After changing the equation, I had to reduce the iso value quite alot to 0.008 is it normal to be that low?

The isovalue really tells you where the limit between solid regions and empty space (the "contour") is to be located. If it is too high or too low, everything will be considered solid or empty, and as you vary it more (or less) of the density field is considered solid, which leads to interesting results. I don't think there is really any "good" value for an isovalue, it depends a lot on your metaball implementation and on the effect you want to achieve, so I wouldn't be worried if you need to tweak it a bit to get good-looking metaballs.

Also I can only specify metaball positions between 0 -> 1 anything outside doesn't get drawn and I want to be able to draw them anywhere on screen, Im not sure which part of the algorithm decides the range?

This shouldn't be a problem with the metaball algorithm itself. I would check your marching cubes code to verify it is working outside the 0..1 unit cube. If it isn't, then obviously everything outside it will never be polygonized, and therefore never rendered. It's not possible to polygonize things arbitrarily far away, unfortunately, there is always a tradeoff between range and resolution, because marching cubes is a discrete algorithm. Other methods like ray marching can sort of scale to arbitrary distances, but have their own set of drawbacks.

To make sure your balls can move around, you can compute the bounding box of your metaballs and use that to define the range for the marching cubes algorithm, but you can't have them go too far away from one another else you will lose in accuracy. You might wonder why, since in that case the marching cubes algorithm would spend most of its time on empty voxels, and, yes, it is possible to do better, but it gets rather nasty as you then need to find an approximate bounding box for every set of connected metaballs and run the marching cubes algorithm on each of them, separately, which sounds great on paper but isn't that efficient in practice. I think most people deal with this by making reasonable tradeoffs between the size of their worlds and how precise the polygonization should be.

Finally just a general question, eventually I wanted to use this to render fluid. Some people seem to use marching cubes to render fluid, however It took my app like 4seconds just to render those four static metaballs. What kind of technique is good for fluid rendering?

Marching cubes can be rather expensive, especially if you want really good resolution. Looking at your screenshot, your mesh is quite smooth, meaning it's probably trying to polygonize up to 128x128x128 voxels, which is *a lot* especially done on the CPU. If you wanted to go interactive - it is possible - you would move the marching cubes code to the GPU, in a shader, and perhaps scale back the resolution a notch. It isn't too hard at all, in fact, and the speed boost is huge since doing the same calculation over a lot of different locations is what the graphics card does best. Then you can add lots of algorithmic optimizations to avoid calculating every single sphere's potential for each voxel - which obviously won't do when you start rendering dozens of metaballs - by using techniques like octrees or spatial hashing. It can get rather efficient, really, but it has its limits. As we all know, the really cool stuff is done by cleverly combining different techniques

I don't know if metaballs would be my first choice for fluid rendering, though. It stills seems like an inefficient method overall, I would just use a grid-based fluid dynamics system if I wanted to do it properly (though a fractal heightmap or even FFT water works well for static water bodies). Marching cubes itself sounds good to display the results, however. Remember to dissociate what you are rendering (metaballs) from how you are rendering it (marching cubes), the two have nothing in common except being often used together.

... great, now I want to write a metaball renderer

Haha you should write one then you can help me make mine better :)

Anyway I have been messing about with it for a while but still cannot draw outside the 0->1 range. I found this source online

https://github.com/kamend/3D-Metaballs-with-Marching-Cubes/tree/master/src (Note its not mine) to try and see how they control range and I see it uses gridX gridY gridZ which is basically the same thing as my m_volume_depth etc...

Can anyone have a look through that source and see which bit alters range so I can draw more than a unit cube.

Thanks

##### Share on other sites

[snip]

Haha you should write one then you can help me make mine better

Anyway I have been messing about with it for a while but still cannot draw outside the 0->1 range. I found this source online

https://github.com/kamend/3D-Metaballs-with-Marching-Cubes/tree/master/src (Note its not mine) to try and see how they control range and I see it uses gridX gridY gridZ which is basically the same thing as my m_volume_depth etc...

Can anyone have a look through that source and see which bit alters range so I can draw more than a unit cube.

Thanks

Well I don't have time to look into it right now (gotta go in five minutes) but reading your original post's code again, you are normalizing your i, j, k counters so that they always fall in (-0.5 .. 0.5), and are always evaluating the metaballs there (so increasing volume_width and so on only increases resolution but not range). A quick fix is to keep the code the same, but multiply your x, y, z values by some factor like "scale". Then scale = 2 would render in (-1.. 1), scale = 4 between (-2..2) and so on. Try it and see if that works.

##### Share on other sites

Well I don't have time to look into it right now (gotta go in five minutes) but reading your original post's code again, you are normalizing your i, j, k counters so that they always fall in (-0.5 .. 0.5), and are always evaluating the metaballs there (so increasing volume_width and so on only increases resolution but not range). A quick fix is to keep the code the same, but multiply your x, y, z values by some factor like "scale". Then scale = 2 would render in (-1.. 1), scale = 4 between (-2..2) and so on. Try it and see if that works.

Wow such a simple fix haha thanks man it worked. Although the resolution automatically decreases when I increase the scale and range. In your opinion what would be the best way to specify those x,y,z values. Should I leave it like it is or is there an alternative method you think would be better.

Again thanks for the help, appreciate it :)

##### Share on other sites

Well I don't have time to look into it right now (gotta go in five minutes) but reading your original post's code again, you are normalizing your i, j, k counters so that they always fall in (-0.5 .. 0.5), and are always evaluating the metaballs there (so increasing volume_width and so on only increases resolution but not range). A quick fix is to keep the code the same, but multiply your x, y, z values by some factor like "scale". Then scale = 2 would render in (-1.. 1), scale = 4 between (-2..2) and so on. Try it and see if that works.

Wow such a simple fix haha thanks man it worked. Although the resolution automatically decreases when I increase the scale and range. In your opinion what would be the best way to specify those x,y,z values. Should I leave it like it is or is there an alternative method you think would be better.

Again thanks for the help, appreciate it

No worries! Yes, these two parameters aren't ideal because they both influence resolution and range. You can however define two parameters which will only influence one of the two variables each, by using the concept of "point per unit length" where a point is where you evaluate the metaballs to create a control point for the marching cubes algorithm. It would go as follows:

Parameters:
-----------

resolution  [points/length]
range       [length]
offset      [length]

(for simplicity each parameter applies to all three dimensions, separate as needed)

Notes:
------

Number of points in each dimension is just resolution * range  [points/length * length = points]

Code:
-----

for (i = 0; i <= resolution * range; ++i)
{
for(j = 0; j <= resolution * range; ++j)
{
for(k = 0; k <= resolution * range; ++k)
{
// corresponds to point (i, j, k), use resolution and offset
float x = (float)i / resolution + offset;
float y = (float)j / resolution + offset;
float z = (float)k / resolution + offset;

// evaluate point at (x, y, z)
}
}
}


Note that now, if you modify the range then the resolution remains unchanged (indeed, you will simply have *more* points to evaluate at the same resolution) and similarly changing the resolution does not modify the effective range (the points will just be more spread out over the original range, trading accuracy for speed). The code above will generate control points between "offset" and "range + offset" inclusive (in every dimension), so tweak as needed. Note you can of course have multiple parameters, one for each dimension if you want to (like offsetY, resolutionX, etc..)

You can have fractional resolution and range by truncating (or rounding? not sure, try it out) the value "resolution * range" inside your for loop. In fact an even more elegant implementation - and perhaps more intuitive - would go as follows:

for (float x = offset; x <= range + offset; x += resolution)
{
for (float y = offset; y <= range + offset; y += resolution)
{
for (float z = offset; z <= range + offset; z += resolution)
{
// evaluate point at (x, y, z)
}
}
}


Which is really saying, start at "offset", and increment the control points in steps of "resolution" until you reach "range + offset". After all, for loops aren't only for integers  (though this may be less optimized and might be subject to floating-point inaccuracies, try it out anyway)

Edited by Bacterius

##### Share on other sites

Well I don't have time to look into it right now (gotta go in five minutes) but reading your original post's code again, you are normalizing your i, j, k counters so that they always fall in (-0.5 .. 0.5), and are always evaluating the metaballs there (so increasing volume_width and so on only increases resolution but not range). A quick fix is to keep the code the same, but multiply your x, y, z values by some factor like "scale". Then scale = 2 would render in (-1.. 1), scale = 4 between (-2..2) and so on. Try it and see if that works.

Wow such a simple fix haha thanks man it worked. Although the resolution automatically decreases when I increase the scale and range. In your opinion what would be the best way to specify those x,y,z values. Should I leave it like it is or is there an alternative method you think would be better.

Again thanks for the help, appreciate it

No worries! Yes, these two parameters aren't ideal because they both influence resolution and range. You can however define two parameters which will only influence one of the two variables each, by using the concept of "point per unit length" where a point is where you evaluate the metaballs to create a control point for the marching cubes algorithm. It would go as follows:

Parameters:
-----------

resolution  [points/length]
range       [length]
offset      [length]

(for simplicity each parameter applies to all three dimensions, separate as needed)

Notes:
------

Number of points in each dimension is just resolution * range  [points/length * length = points]

Code:
-----

for (i = 0; i <= resolution * range; ++i)
{
for(j = 0; j <= resolution * range; ++j)
{
for(k = 0; k <= resolution * range; ++k)
{
// corresponds to point (i, j, k), use resolution and offset
float x = (float)i / resolution + offset;
float y = (float)j / resolution + offset;
float z = (float)k / resolution + offset;

// evaluate point at (x, y, z)
}
}
}


Note that now, if you modify the range then the resolution remains unchanged (indeed, you will simply have *more* points to evaluate at the same resolution) and similarly changing the resolution does not modify the effective range (the points will just be more spread out over the original range, trading accuracy for speed). The code above will generate control points between "offset" and "range + offset" inclusive (in every dimension), so tweak as needed. Note you can of course have multiple parameters, one for each dimension if you want to (like offsetY, resolutionX, etc..)

You can have fractional resolution and range by truncating (or rounding? not sure, try it out) the value "resolution * range" inside your for loop. In fact an even more elegant implementation - and perhaps more intuitive - would go as follows:

for (float x = offset; x <= range + offset; x += resolution)
{
for (float y = offset; y <= range + offset; y += resolution)
{
for (float z = offset; z <= range + offset; z += resolution)
{
// evaluate point at (x, y, z)
}
}
}


Which is really saying, start at "offset", and increment the control points in steps of "resolution" until you reach "range + offset". After all, for loops aren't only for integers  (though this may be less optimized and might be subject to floating-point inaccuracies, try it out anyway)

Sorry for the late reply, I work on this as a side project, anyway I tried your methods and just couldn't get it to work haha :( the balls just don't show at all and Im not really fussed to debug and see why. I just left it like it was previously.

I added an update function and random velocity to see the performance and it is not good. I can only do about 32 metaballs with voxelisation at 32^3. It's most likely my bad coding haha so I need to research further into rendering fluid.

Thanks for all the help.

##### Share on other sites
I doubt you can get interactive rates with a CPU approach, unfortunately. For fluids you probably even have to do the simulation on the GPU.

There's source code for GPU marching cubes available, though: Procedural Terrains using the GPU. Can't find the source code download, you probably need to register (it's free).

There's also this "fake" approach as an alternative: Fluid simulation in Alice Madness Returns.