General formula to calculate Polyhedron volume

Question

Given a list of vertices (v), and a list of edges connecting the vertices (e), and a list of surfaces that connect the edges (s), how to calculate the volume of the Polyhedron?

Answer 1

Actually it's not necessary to divide the polyhedron's polygons into triangles. All that is needed is to sum the volumes of the pyramids formed by each polygon with an apex at the origin using the one-third base area times height rule. The base area can be computed using the shoelace method and the height can be accounted for by forming the dot product of that (directed) area with any one vertex of the given polygon. It is apparent that the volume of each pyramid is identical to the sum of the volumes formed from a set of triangles that make up its base polygon, so it must generate the same answer. Potentially this is a much faster method if the polyhedron has faces that are intricate, perhaps reentrant, polygons when dividing them into triangles is not a simple task.

After some thought I worked out there is a simple formula for the volume contribution of an n-sided polygon:

formula for volume contribution of n-sided polygon in terms of the co-ordinates of its vertices xi, yi, zi

Answer 2

1. Take the polygons and break them into triangles.
2. Consider the tetrahedron formed by each triangle and an arbitrary point (the origin).
3. Sum the signed volumes of these tetrahedra.

Notes:

1. This will only work if you can keep a consistent CW or CCW order to the triangles as viewed from the outside.
2. The signed volume of the tetrahedron is equal to 1/6 the determinant of the following matrix:

    [ x1 x2 x3 x4 ]
    [ y1 y2 y3 y4 ]
    [ z1 z2 z3 z4 ]
    [  1  1  1  1 ]

where the columns are the homogeneous coordinates of the verticies (x,y,z,1).

It works even if the shape does not enclose the origin by subracting off that volume as well as adding it in, but that depends on having a consistent ordering.

If you can't preserve the order you can still find some way to break it into tetrahedrons and sum 1/6 absolute value of the determinant of each one.

Edit: I'd like to add that for triangle mesh where one vertex (say V4) of the tetrahedron is (0,0,0) the determinante of the 4x4 matrix can be simplified to the upper left 3x3 (expansion along the 0,0,0,1 column) and that can be simplified to Vol = V1xV2.V3 where "x" is cross product and "." is dot product. So compute that expression for every triangle, sum those volumes and divide by 6.

Answer 3

First, break every face into triangles by drawing in new edges.

Now look at one triangle, and suppose it's on the "upper" surface (some of these details will turn out to be unimportant later). Look at the volume below the triangle, down to some horizontal plane below the polyhedron. If {h1, h2, h3} are the heights of the three points, and A is the area of the base, then the volume of the solid will be A(h1+h2+h3)/3. Now we have to add up the volumes of these solids for the upper faces, and subtract them for the lower faces to get the volume of the polyhedron.

Play with the algebra and you'll see that the height of the polyhedron above the horizontal plane doesn't matter. The plane can be above the polyhedron, or pass through it, and the result will still be correct.

So what we need is (1) a way to calculate the area of the base, and (2) a way to tell an "upper" face from a "lower" one. The first is easy if you have the Cartesian coordinates of the points, the second is easy if the points are ordered, and you can combine them and kill two birds with one stone. Suppose for each face you have a list of its corners, in counter-clockwise order. Then the projection of those points on the x-y plane will be counterclockwise for an upper face and clockwise for a lower one. If you use this method to calculate the area of the base, it will come up positive for an upper face and negative for a lower one, so you can add them all together and have the answer.

So how do you get the ordered lists of corners? Start with one triangle, pick an ordering, and for each edge the neighbor that shares that edge should list those two points in the opposite order. Move from neighbor to neighbor until you have a list for every triangle. If the volume of the polyhedron comes up negative, just multiply by -1 (it means you chose the wrong ordering for that first triangle, and the polyhedron was inside-out).

EDIT: I forgot the best part! If you check the algebra for adding up these volumes, you'll see that a lot of terms cancel out, especially when combining triangles back into the original faces. I haven't worked this out in detail, but it looks as if the final result could be a surprisingly simple function.

Answer 4

Here's a potential implementation for that in Python. Can anyone please check if it's correct? I believe that I am missing permutations of the points because my second test (cube) gives 0.666 and not 1. Ideas anyone?

Cheers EL

class Simplex(object):
    '''
    Simplex
    '''


    def __init__(self,coordinates):
        '''
        Constructor
        '''

        if not len(coordinates) == 4:
            raise RuntimeError('You must provide only 4 coordinates!')

        self.coordinates = coordinates



    def volume(self):
        '''
        volume: Return volume of simplex. Formula from http://de.wikipedia.org/wiki/Tetraeder
        '''
        import numpy

        vA = numpy.array(self.coordinates[1]) - numpy.array(self.coordinates[0])
        vB = numpy.array(self.coordinates[2]) - numpy.array(self.coordinates[0])
        vC = numpy.array(self.coordinates[3]) - numpy.array(self.coordinates[0])

        return numpy.abs(numpy.dot(numpy.cross(vA,vB),vC)) / 6.0  


class Polyeder(object):

    def __init__(self,coordinates):
        '''
        Constructor
        '''

        if len(coordinates) < 4:
            raise RuntimeError('You must provide at least 4 coordinates!')

        self.coordinates = coordinates


    def volume(self):

        pivotCoordinate = self.coordinates[0]

        volumeSum = 0

        for i in xrange(1,len(self.coordinates)-3):

            newCoordinates = [pivotCoordinate]

            for j in xrange(i,i+3):
                newCoordinates.append(self.coordinates[j])

            simplex = Simplex(newCoordinates)
            volumeSum += simplex.volume()

        return volumeSum


coords = []

coords.append([0,0,0])
coords.append([1,0,0])
coords.append([0,1,0])
coords.append([0,0,1])

s = Simplex(coords)
print s.volume()

coords.append([0,1,1])
coords.append([1,0,1])
coords.append([1,1,0])
coords.append([1,1,1])

p = Polyeder(coords)
print p.volume()

Answer 5

I have done this before, but the surface mesh I used always had triangular facets. If your mesh has non triangular facets, you can easily break them up into triangular facets first. Then I fed it to TetGen to obtain a tetrahedralization of the interior. Finally, I added up all the volumes of the tetrahedra. TetGen is reasonably easy to use, and is the only library other than CGAL I know of that can handle complicated meshes. CGAL is pretty easy to use if you don't mind installing a gigantic library and use templates like crazy.

Answer 6

Similarly with a polygon where we can split it into triangles and sum the areas,
you could split a polyhedron into pyramids and sum their volumes. But I'm not sure how hard is to implement an algorithm for that.

(I believe there is a mathematical way/formula, like using vectors and matrices.
I suggest to post your question also on http://mathoverflow.net)