Metaballs

As a little exercise to get back to 3D-programming I wrote a small app which can render scalar fields, e.g. metaballs. The app is static, and I would need more optimization to get an animation:







The core “polygonisator” class takes a given scalar field and gives back a trimesh. The used algorithm is “marching tetrahedrons”. Don’t know if my tex coords make sense, but the normals seem to be OK.



package metaballs;

import com.jme.math.FastMath;
import com.jme.math.Triangle;
import com.jme.math.Vector3f;
import com.jme.scene.TriMesh;
import com.jme.scene.batch.TriangleBatch;
import com.jme.util.geom.BufferUtils;
import java.util.ArrayList;
import java.util.List;

/**
 * Uses  Paul Bourke's code from "Polygonising a Scalar Field Using Tetrahedrons"
 * http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/
 */
public class Polygonisator {

    /** do not instantiate */
    private Polygonisator() {
    }

    public static TriMesh getMesh(ScalarField field, Vector3f boxSize, float cubeSize, float iso) {
        List<Triangle> triangles = new ArrayList<Triangle>();

        for (float x = -boxSize.x; x <= boxSize.x; x += cubeSize) {
            for (float y = -boxSize.y; y <= boxSize.y; y += cubeSize) {
                for (float z = -boxSize.z; z <= boxSize.z; z += cubeSize) {
                    GridCell cell = new GridCell(field, x, y, z, cubeSize);
                    if (cell.isActive(iso)) {
                        calculateCell(cell, triangles, iso);
                    }
                }
            }
        }

        if (triangles.isEmpty()) {
            return null;
        }


        TriMesh mesh = new TriMesh("mesh");
        TriangleBatch batch = mesh.getBatch(0);

        int[] indices = new int[triangles.size() * 3];
        for (int i = 0; i < indices.length; i++) {
            indices[i] = i;
        }
        batch.setIndexBuffer(BufferUtils.createIntBuffer(indices));

        batch.setVertexBuffer(BufferUtils.createVector3Buffer(batch.getVertexBuffer(), triangles.size() * 3));
        batch.setNormalBuffer(BufferUtils.createVector3Buffer(triangles.size() * 3));
        batch.setVertexCount(triangles.size() * 3);
        batch.getTextureBuffers().set(0, BufferUtils.createVector2Buffer(triangles.size() * 3));

        for (Triangle triangle : triangles) {
            batch.getVertexBuffer().put(triangle.get(0).x).put(triangle.get(0).y).put(triangle.get(0).z).put(triangle.get(1).x).put(triangle.get(1).y).put(triangle.get(1).z).put(triangle.get(2).x).put(triangle.get(2).y).put(triangle.get(2).z);
            triangle.calculateNormal();
            for (int i = 0; i < 3; i++) {
                batch.getNormalBuffer().put(triangle.getNormal().x).put(triangle.getNormal().y).put(triangle.getNormal().z);

                Vector3f v = triangle.get(i).normalize();
                Vector3f vp = new Vector3f(v.x, 0, v.y).normalizeLocal();
                batch.getTextureBuffer(0).put(Vector3f.UNIT_X.angleBetween(vp) / FastMath.PI).put(v.angleBetween(vp) / FastMath.PI);
            }
        }

        return mesh;

    }

    private static void calculateCell(GridCell cell, List<Triangle> triangles, float iso) {

        calculateTetra(cell, iso, triangles, 0, 4, 7, 6);
        calculateTetra(cell, iso, triangles, 0, 4, 6, 5);
        calculateTetra(cell, iso, triangles, 0, 2, 6, 3);
        calculateTetra(cell, iso, triangles, 0, 1, 6, 2);
        calculateTetra(cell, iso, triangles, 0, 3, 6, 7);
        calculateTetra(cell, iso, triangles, 0, 1, 5, 6);

    }

    private static void calculateTetra(GridCell cell, float iso, List<Triangle> triangles,
            int v0, int v1, int v2, int v3) {
        int triindex =
                (cell.values[v0] < iso ? 1 : 0) + (cell.values[v1] < iso ? 2 : 0)
              + (cell.values[v2] < iso ? 4 : 0) + (cell.values[v3] < iso ? 8 : 0);

        Vector3f[] vec = new Vector3f[4];

        /* Form the vertices of the triangles for each case */
        switch (triindex) {
            case 0x00:
            case 0x0F:
                break;
            case 0x0E:
            case 0x01:
                vec[0] = interpolate(iso, cell, v0, v1);
                vec[1] = interpolate(iso, cell, v0, v2);
                vec[2] = interpolate(iso, cell, v0, v3);
                triangles.add(triindex == 0xE ? new Triangle(vec[0], vec[1], vec[2]) : new Triangle(vec[0], vec[2], vec[1]));
                break;
            case 0x0D:
            case 0x02:
                vec[0] = interpolate(iso, cell, v1, v0);
                vec[1] = interpolate(iso, cell, v1, v3);
                vec[2] = interpolate(iso, cell, v1, v2);
                triangles.add(triindex == 0x0D ? new Triangle(vec[0], vec[1], vec[2]) : new Triangle(vec[0], vec[2], vec[1]));
                break;
            case 0x0C:
            case 0x03:
                vec[0] = interpolate(iso, cell, v0, v3);
                vec[1] = interpolate(iso, cell, v0, v2);
                vec[2] = interpolate(iso, cell, v1, v3);
                vec[3] = interpolate(iso, cell, v1, v2);
                triangles.add(triindex == 0x0C ? new Triangle(vec[0], vec[2], vec[1]) : new Triangle(vec[0], vec[1], vec[2]));
                triangles.add(triindex == 0x0C ? new Triangle(vec[2], vec[3], vec[1]) : new Triangle(vec[2], vec[1], vec[3]));
                break;

            case 0x0B:
            case 0x04:
                vec[0] = interpolate(iso, cell, v2, v0);
                vec[1] = interpolate(iso, cell, v2, v1);
                vec[2] = interpolate(iso, cell, v2, v3);
                triangles.add(triindex == 0x0B ? new Triangle(vec[0], vec[1], vec[2]) : new Triangle(vec[0], vec[2], vec[1]));
                break;
            case 0x0A:
            case 0x05:
                vec[0] = interpolate(iso, cell, v0, v1);
                vec[1] = interpolate(iso, cell, v2, v3);
                vec[2] = interpolate(iso, cell, v0, v3);
                vec[3] = interpolate(iso, cell, v1, v2);
                triangles.add(triindex == 0x0A ? new Triangle(vec[0], vec[1], vec[2]) : new Triangle(vec[0], vec[2], vec[1]));
                triangles.add(triindex == 0x0A ? new Triangle(vec[0], vec[3], vec[1]) : new Triangle(vec[0], vec[1], vec[3]));
                break;
            case 0x09:
            case 0x06:
                vec[0] = interpolate(iso, cell, v0, v1);
                vec[1] = interpolate(iso, cell, v1, v3);
                vec[2] = interpolate(iso, cell, v2, v3);
                vec[3] = interpolate(iso, cell, v0, v2);
                triangles.add(triindex == 0x09 ? new Triangle(vec[0], vec[1], vec[2]) : new Triangle(vec[0], vec[2], vec[1]));
                triangles.add(triindex == 0x09 ? new Triangle(vec[0], vec[2], vec[3]) : new Triangle(vec[0], vec[3], vec[2]));
                break;
            case 0x07:
            case 0x08:
                vec[0] = interpolate(iso, cell, v3, v0);
                vec[1] = interpolate(iso, cell, v3, v2);
                vec[2] = interpolate(iso, cell, v3, v1);
                triangles.add(triindex == 0x07 ? new Triangle(vec[0], vec[1], vec[2]) : new Triangle(vec[0], vec[2], vec[1]));
                break;
        }
    }

    private static Vector3f interpolate(float iso, GridCell cell, int v1, int v2) {
        float ratio = 0.5f;
        if (cell.values[v1] != cell.values[v2]) {
            ratio = (cell.values[v1] - iso) / (cell.values[v1] - cell.values[v2]);
        }
        if (ratio < 0 || ratio > 1) {
            throw new IllegalArgumentException("Ratio: " + ratio);
        }
        return cell.points[v1].mult(1 - ratio).add(cell.points[v2].mult(ratio));
    }

    private static class GridCell {

        Vector3f[] points;
        float[] values;

        public GridCell(ScalarField field, float x, float y, float z, float cubeSize) {
            points = new Vector3f[]{
                new Vector3f(x, y, z),
                new Vector3f(x + cubeSize, y, z),
                new Vector3f(x + cubeSize, y, z + cubeSize),
                new Vector3f(x, y, z + cubeSize),
                new Vector3f(x, y + cubeSize, z),
                new Vector3f(x + cubeSize, y + cubeSize, z),
                new Vector3f(x + cubeSize, y + cubeSize, z + cubeSize),
                new Vector3f(x, y + cubeSize, z + cubeSize)
           
               ,
            };
            values = new float[8];
            for (int i = 0; i < values.length; i++) {
                values[i] = field.calculate(points[i]);
            }
        }

        private boolean isActive(float iso) {
            int v = 0;
            for (float value : values) {
                v += (value > iso) ? 1 : -1;
            }
            return -8 < v && v < 8;
        }
    }
}




package metaballs;

import com.jme.math.Vector3f;

interface ScalarField {
     public float calculate(Vector3f point);
}



The screenshot uses just the sum of weight/distance over some metaballs as scalar field.
1 Like

Cool man…  Now we just need a spaghetti generator and we're set!

Metaballs are awesome…if we could only integrate with jME-Physics it would be completely awesome. :slight_smile:

After some optimization they billow with ~10FPS over my screen.







The grid size has to be real low for animation, but I improved the normals, so I can use now env mapping.



There is still a lot of work to do: I don’t reuse vertices but recalculate them for every single triangle. So I guess I could get around 20FPS with further tweaking…

man, I hope you can optimize further or you have a crappy machine. :wink:

Old machine, crappy graphic card. I’m currently at 13…15 FPS and it starts to look smooth. But if the balls are too small, they look ugly…



Here are the sources:

http://rapidshare.com/files/78527097/metaballs.zip.html

I know how to get it faster:



The current Polygonisator class knows nothing about the scalar field it draws. If we want faster metaballs, I could write a specialized version, using the following assumptions:

  • we know a point inside each separated part (in our case the metaball centers)
  • the "inside" volume (scalar field bigger than iso value) is much smaller than the "outside" volume (scalar field smaller than iso value)



    To meet this conditions it is required that the "force" of a meta ball center is always positive and decreasing with the distance.



    That configuration would allow to start with the mini cubes inside a part and explore its neighbors until finding all mini cubes at the iso surface. All "outside" mini cubes are not calculated.

Does trigger of many a thought for its use

I played aroud with both approaches.



The non specialized "brute force" version is now at 13~16 FPS, quite independent from the size of the surface.



The specialized Metaball-version is not so fast as I expected, very sensitive concerning the size of the surface, 8~30FPS. However I have the feeling that there is still room for improvement…

Added my current code (faster than the download above) to the Wiki:



http://www.jmonkeyengine.com/wiki/doku.php?id=rendering_scalar_fields

Please go ahead!  :smiley:



Maybe the marching tetrahedron class could be used later as a part of a jmex.voxel package or so…



Thanks for your improvements!

Cool, thanks. Added to SVN.



I've added texture coords and color btw - I thought deleting it if it's not needed would be easier than adding for someone who might want to use it.

1 Like

I found it at http://code.google.com/p/jmonkeyengine/source/browse/trunk/src/jmetest/scalarfields/?r=4101

EDIT: I have it updated to JME3, how can I contribute it back?