cad - How do I calculate the volume of an object stored in STL files?

Question

I have .stl (STL is a file format native to the stereolithography CAD software created by 3D Systems) files, from which I must calculate the volume. How do I do this? I'm using the calculation below but the volume is not equal to that calculated by other software

float x1,y1,z1;
float x2,y2,z2;
float x3,y3,z3;

The above are the vertices. triangles is just a data structure holding an object with the vertices of each triangle.

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;

Do you also involve the calculation with the normal vector?

Answer 1

The volume of polyhedra is discussed in the polyhedron article in Wikipedia. Any valid STL tessellation must be precisely a polyhedron where each facet is a side, so the formula discussed there holds.

Assuming the vertices are oriented counter-clockwise (looking down the outward normals) the following expansion holds:

Equation 1 is just the divergence formula for polyhedra. The p[i,j] are the jth vertex (as a vector from the origin) of the ith triangle.

Equation 2 expands this by using the fact that the the cross product of the two triangle sides (a parallelogram) is a normal vector of the triangle with magnitude twice the area of the triangle.

By the way, this is why the triangle normal doesn't show up in your expression even though your intuition says it should be there somehow. It's already in the soup!

Equation 3 just expands equation 2 using the distributive properties of the cross product, and that fact that a vector cross itself is the zero vector.

The volume you get from equation 3 is signed according to the orientation of the triangles. I've aligned my formulation to agree with your code, but you may have a negative result. If so, just take the absolute value.

Now writing out the summand in your notation (except that I don't bother putting in the repeated triangle[i]) I get:

(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 - x2 y1 z3 + x1 y2 z3)/6.;

This exactly matches what you've written (except for order)!

So there are a couple of possibilities:

(1) Perhaps your STL file is flawed and the triangles aren't all oriented consistently. You can check this by verifying that each edge is used by a triangle once in the forward direction, and once in the opposite direction. By and "edge" I mean a pair of vertices (p[r], p[s]) that are part of the same triangle. I'm changing notation and using the subscript to indicate the uniquely indexed vertex in the file, not the index of the vertex relative to the face as I did with the two-index notation above).

If a triangle has edge (p[r], p[s]) in its list of oriented edges, then some other triangle must have (p[s], p[r]) (same vertices, reversed order) in it's list of oriented edges.

(2) The tessellation isn't closed. This is the same test as 1, except that it's open if only one triangle uses the edge, regardless of orientation.

(3) Dare I mention it? Are you sure you initialized currentVolume?

(4) You get the right answer, but the volume is negative and you're fooled into thinking the method is wrong. See my comments above.