Rendering Metaballs in 3D

Example of metaballs rendering in 3D

3D Metaballs are a valuable tool for modeling "blobby shapes", or something that should look "organic". They could even be used for fake fluid simulations. In this article, we'll explore how to render Metaballs in 3D, using ray casting tecniques.

The Metaballism 1 video series by Alexander Mitzkus showcases the versatility of metaballs in 3D graphics.

If you're new to metaballs, it's recommended to read the previous article in the series, Rendering Metaballs in 2D.

Info

In this page you'll find some shaders written with Shadertoy, and Desmos graphs. Read how to use the interactive content in this site.

The code in this article is took from the shader below, an implementation of Metaballs I made in ShaderToy:

3D Metaballs

Rendering metaballs in 3D

Rendering Metaballs in 3D is more complex than in 2D, but can be achieved using ray casting tecniques. Colors and shading enhance the perception of three-dimensionality, while animation creates dynamic shapes that continuously glue and split.


Isosurfaces of a 3D scalar field

Similar to 2D Metaballs

3D metaballs are isosurfaces of a three dimensional scalar field.

The field in 3D is defined as the sum of 3D scalar functions, each one of them has its own falloff. The following picture shows an example of a the scalar field that is the sum of three functions: $f1, f2, f3$ (Just imagine that this is a section of the the 3D space :). The value of each on decreases as the points are far from its "center".

The 3D scalar field

An isosurface is made from all the points in space that have equal value. The next figure shows an example (imagine that is a section in 2D :). The isosurfaces of a value of choice, are the metaballs that will be rendered. Chosing small values will make big metaballs, and vice versa.

Isosurface

The following code illustrates the definition of the field function. The position is the center of the function (its greatest value), while the radius is the lenght of the falloff. In this example eight fields are defined.

 1// A field function
 2struct Field
 3{
 4    vec3 position;
 5    float radius;
 6};
 7
 8// Define the number of functions (scalar fields) that generate the total field
 9#define N_FIELDS 8
10Field fields[N_FIELDS] = Field[N_FIELDS]
11(
12     Field(vec3(0,0,0),1.)
13    ,Field(vec3(0,0,0),2.)
14    ,Field(vec3(0,0,0),3.)
15    ,Field(vec3(0,0,0),4.)
16    ,Field(vec3(0,0,0),5.)
17    ,Field(vec3(0,0,0),2.)
18    ,Field(vec3(0,0,0),3.)
19    ,Field(vec3(0,0,0),4.)
20);

The following snippet computes the total field, summing up all the functions. It uses getFallOff() to get the value at a specific distance from the center. The falloff will be defined in the following section.

 1    float f, df;
 2    value = 0.;
 3    
 4    // Compute the total field summing all the scalar field ...
 5    for(int i=0;i<N_FIELDS;i++)
 6    {
 7        // The distance from the center of the field
 8	float r = length(p-fields[i].position);
 9       
10        // The distance from the center of the function is normalized with the radius of the field
11        float d = length((p-fields[i].position)/fields[i].radius);
12        
13        // Get the value of the falloff ...
14        getFalloff(d, f);
15        
16        // ... and sum it to the total
17        value += f;        
18    }

The falloff function

In 2D we used a cubic polynomial. In 3D instead it's better to choose the quintic polynomial $$f(x) = 1 - (6x^5 - 15x^4 + 10x^3)$$. It has the first and second derivative equals to zero, so we can avoids discontinuities in the isosurfaces normals 2 3.

The following graph shows a comparison of the quintic and the cubic falloffs.

These polynomials are defined in the interval [0,1], and have a value of $0$ outside of it. It is possible to modify the lenght of the falloff multipling the argument $x$ for a value. If $value<1$, the falloff become larger, otherwise it shrinks.

 1// Get the value of the falloff function for a value of x in [0,1]
 2void getFalloff(float x, out float f)
 3{
 4    f = 0.;
 5    if(x<0.01||x>1.) return; // Field is 0 outside [0,1]
 6    
 7    // Quintic falloff 1-6x^5 - 15x^4 + 10x^3
 8    f = 1.-(x*x*x*(6.*x*x-15.*x + 10.));
 9}
10
11// Get the falloff derivative
12void getFalloffDerivative(float x, out float df)
13{
14    df = 0.;
15    if(x<0.01||x>1.) return; // Field is 0 outside [0,1]
16   
17    // Quintic fallof derivative 1-(30x^4 - 60x^2 + 30x)
18    df = -(x*x*(30.*x*x - 60.*x + 30.));
19}

Ray tracing the bounding spheres

Bounding spheres are used to optimize ray tracing. These spheres encapsulate the Metaballs and allow for efficient boundary detection. Ray tracing determines the hit point on these bounding spheres. The next picture illustrate this concept:

Ray tracint to the field bounds, then ray march till the isosurface

The next snippet implements this optimization. The function traceSphere is defined in the ShaderToy shader.

 1// Ray trace the bounding sphere and returns the hit point (if any)
 2Hit rayTraceFieldBoudingSpheres(vec3 origin, vec3 dir)
 3{
 4    float minT = INF;
 5    float point;
 6    Hit hit, result;
 7    
 8    // Get nearest point from the spheres
 9    result.t = INF;
10    for(int i=0; i<N_FIELDS; i++)
11    {
12        if(traceSphere(origin, dir, fields[i].position, fields[i].radius, hit) && hit.t<result.t) 
13        {
14            result.t = hit.t;
15            result.point = origin + result.t*dir;
16            result.objId = 1;
17        }
18    }
19    
20    return result;
21}

Ray marching inside the field

Ray marching is employed inside the scalar field to locate the isosurface. By taking small steps along the ray direction, the field value is computed at each step. When the desired isosurface value is reached, the normal and shading are computed.

The following code illustrates the ray marching.

 1
 2// Cast a ray into the scene
 3vec3 castRay(vec3 origin, vec3 dir)
 4{
 5    // Update the position of the fields in the scene
 6    updateField();
 7    
 8    vec3 fieldNormal;
 9    
10    float MAX_DISTANCE = 20.;
11    int MAX_ITERATIONS = 200;
12    
13    // getDistanceToField Ray-Trace the field spheres bounding boxes Get to the field nearest point, ray
14    Hit hit = rayTraceFieldBoudingSpheres(origin, dir);
15    
16    // Ray march inside the scalar field
17    float t = hit.t;
18    if(t==INF) return vec3(1); // No hit -> return white background
19   
20    vec3 p;
21    vec3 color;
22    float value=0.;
23    float threshold = 0.4;
24    for(int i=0; i<MAX_ITERATIONS; i++)
25    {
26        if(t>MAX_DISTANCE) break; // No hit
27		
28		// Take little stesp inside the field
29        t+=0.05;
30        p = origin + t*dir;
31        
32		// Get the field value, normal and color
33        if(getFieldValues(p, threshold, value, fieldNormal,color))
34        {
35            hit.t=t;
36            hit.point=p;
37            hit.normal=fieldNormal;
38            hit.objId=1;            
39            hit.color=color;
40            break;
41            
42        }
43    }
44    
45    // If the ray didn't intersect the isosurface (here I choose the value 0.4)
46    // -> return white background 
47    if(value < 0.4) return vec3(1); 
48    
49    return render(origin,hit);
50}

Computing the normals

The normals to the isosurface are equal to the sum of each function normals. All the normals have the direction from the function center to the application point. The magnitude is equals to the derivative of the falloff function.

 1    
 2    // If the value of the field is greater than threshold we reach the isosurface,
 3    // so compute the normal
 4    if(value>=threshold)
 5    {
 6        float df=0.;
 7        
 8        // The resulting normal is the sum of all the gradient of the fields function summed up.
 9        // The gradient of a single field is the derivative of its falloff in this point (that is a scalar)
10        // multiplied the normal vector to the surface, that is the differece between the current position 
11        // and the center of the field (it's a spherical field).
12        for(int i=0;i<N_FIELDS;i++) 
13        {
14            // Get the value of the falloff derivative in this point ...
15            float d = length((p-fields[i].position)/fields[i].radius);
16            getFalloffDerivative(d, df);
17            
18            // Compute the normal of this specific field in this position and sum to the others
19            normal += df*normalize(fields[i].position-p);
20        }
21        
22        // Normalize the result
23        normal = normalize(normal);
24        
25        return true;
26    }

Rendering the metaballs

Shading is applied to the Metaballs using a simple shading model that considers diffuse, specular, and Fresnel contributions. The shading is implemented in the render() function.


Animating the metaballs

Animating Metaballs involves changing the positions of the functions that generate the field over time. This animation creates dynamic, evolving shapes.

 1// Update the position of the objects in the scene
 2void updateField()
 3{
 4    for(int i=0; i<N_FIELDS; i++)
 5    {
 6        float id = float(i);
 7        float a = 2.;
 8		
 9        // Just change the position of each function using a random combination of sin, cos, *, +  and scalar values
10        fields[i].position = vec3(a*sin(float(iFrame)/200.+id*244.),1.4*abs(sin(id*0.1*iTime+323.3))-1.,a*sin(float(iFrame)/100.+id*1724.)) * 2.;   
11    }
12}

Experimenting with it

With the basics covered, there are endless possibilities for experimentation. For example, different functions can be used to generate the field, leading to unique visual effects. The provided ShaderToy shader includes additional features, such as a field generated from a plane, for further exploration.




That concludes our guide to rendering 3D Metaballs using ray casting techniques. Experiment with the provided shaders and explore the creative possibilities of Metaballs in 3D graphics.

Posts in this series