A method to calculate the centre of mass from a .stl (stereo lithography) file?

Question

I am trying to calculate the centre of mass (x,y,z) coordinates of an object defined in an STL file (stereo lithography, not to be confused with the standard template library). The STL file contains a closed object (or objects) defined by a boundary made of triangles. The triangles themselves are not necessarily in any order, the file is simply the coordinates 3 vertices of each triangle floating in 3D space plus a normal vector to the triangle (the normal should be disregarded as it is not always done properly). There is nothing that links each triangle to one another, it is assumed that the object is closed.

One simple approach would be to divide a volume (in this case, a box) into millions of elements and determine if each element is inside the object defined in the STL file or not, then sum up the moments and calculate the centre of mass. This would work but its far from elegant and extremely slow.

Another method would be to convert the boundary representation into a number of packed tetrahedron solids. Form that I could calculate the centre of mass of each tetrahedron, its volume, and resulting moment and thus calculate the overall centre of mass from the sum of all tetrahedrons. The problem with this is that I don't know how to convert a surface representation of triangles into a volume representation of tetrahedrons (I'm assuming its a fairly non trivial task).

Dose anyone know of any methods or can think up of any methods that I could try? Or maybe even any reference material that talks about this?

For more information about STL files (only the first 2 sections are important, everything else is useless): http://en.wikipedia.org/wiki/STL_%28file_format%29

Answer 1

After a lot of thinking and experimentation I have the answer!

First we add a 4th point to each triangle in to make them into tetrahedrons with a volume centroid. We calculate the volumes and centres of masses and multiply them by each other to get our moments. We sum the moments and divide by total volume to get our overall centroid.

We calculate volumes using the determinate method shown here (equation 32): http://mathworld.wolfram.com/Tetrahedron.html

The centroids of each of the tetrahedrons is simply the average of the 4 points.

The trick here is that due to the way the STL file is created, the triangles have a normal that point outwards from the part surface, following the right hand rule of the 3 verticies used to create the triangle. we can use this to our advantage by allowing us to have a consistent convention in which to determine if a volume of the tetrahedron should be added or subtracted from our net part (this is because the reference point we chose may not necessarily be inside the part and the overall part is not necessarily convex, it is, however a closed object).

Using the determine method to calculate the volume, the first three coordinate points will represent the three points of our triangle. The fourth point would be our common origin. If the normal created by the triangle (following the right hand rule going from point 1, 2, 3) points towards our common reference point, that volume will be calculated as not part of our overall solid, or negative volume (by pointing towards, i mean the vector created by the triangle's normal is pointing loosely towards the same side as a normal plane created by the vector from our reference point to the centroid of the tetrahedron). If the vector is pointing away from the reference point, it is then positive volume or inside the part. If it is normal, then the volume goes to zero as the triangle is in the same plane as the reference point.

We don't need to worry about actually keeping track of any of this as if we are consistent with our inputs (as in the triangles follow the right hand rule with normal facing outwards from the part) the determine will give us the correct sign.

Anyways, heres the code (its even more simple than the explanation).

class data // 3 vertices of each triangle
{
public:
    float x1,y1,z1;
    float x2,y2,z2;
    float x3,y3,z3;
};

int main ()
{
    int numTriangles; // pull in the STL file and determine number of triangles
    data * triangles = new triangles [numTriangles];
    // fill the triangles array with the data in the STL file

    double totalVolume = 0, currentVolume;
    double xCenter = 0, yCenter = 0, zCenter = 0;

    for (int i = 0; i < numTriangles; i++)
    {
        totalVolume += currentVolume = (triangles[i].x1*triangles[i].y2*triangles[i].z3 - triangles[i].x1*triangles[i].y3*triangles[i].z2 - triangles[i].x2*triangles[i].y1*triangles[i].z3 + triangles[i].x2*triangles[i].y3*triangles[i].z1 + triangles[i].x3*triangles[i].y1*triangles[i].z2 - triangles[i].x3*triangles[i].y2*triangles[i].z1) / 6;
        xCenter += ((triangles[i].x1 + triangles[i].x2 + triangles[i].x3) / 4) * currentVolume;
        yCenter += ((triangles[i].y1 + triangles[i].y2 + triangles[i].y3) / 4) * currentVolume;
        zCenter += ((triangles[i].z1 + triangles[i].z2 + triangles[i].z3) / 4) * currentVolume;
    }

    cout << endl << "Total Volume = " << totalVolume << endl;
    cout << endl << "X center = " << xCenter/totalVolume << endl;
    cout << endl << "Y center = " << yCenter/totalVolume << endl;
    cout << endl << "Z center = " << zCenter/totalVolume << endl;
}

Extremely fast for calculating centres of mass for STL files.

Answer 2

EDIT: Look up "winding number algorithm" or "crossing number algorithm" - what I try to describe below is a 3-d crossing number algorithm.

I've got a feeling that something like this will work, but I don't have the ability to test it out, right now:

Build the filled-in 3-d structure from the triangles in the STL file iteratively. Start by picking a single point to use as a basis for the 3-d structure. Then, begin your structure by creating a triangular pyramid, with base defined by the first triangle in the STL file, and vertex your chosen point. Each such component of your iteratively built volume would also contain an "intersection parity" - initialize it to 0.

For each subsequent triangle in the STL file, create a similar pyramid, and see if it intersects with the 3-d structure that you've built so far. If it does, calculate the intersection, and segment the existing structure and the new pyramid so that no two components overlap. Keep the "intersection parity" of the outermost part of the new polyhedron 0, but toggle it on all inner portions of the intersection -- if it was 0, make it 1, if it was 1, make it 0.

At the end, you'll have the closed polyhedron defined by all the portions of your structure that have intersection parity 0. Calculate the moments of all of these polyhedrons, and average them together to get your center of mass. I think the complexity would be something like O(n^2).