Computing 3D polygon volume

I wanted to compute the volume of a mesh in Maya. It was surprisingly simple and elegant to do as well! Using the Divergent Theorem (which is unreadable to me when written mathematically) the only constraints are: the mesh must be closed (no holes, borders or tears; Maya’s fill-hole can help), the mesh must be triangulated (using the Maya API you can already query triangles so no need to manually triangulate in this case).

Now imagine to compute the volume of a prism. All you need is the area of the base triangle * the height. To compute the base I use Heron’s formula as described here.

def distance(a, b):
    return sqrt((b[0]-a[0])*(b[0]-a[0])+
      (b[1]-a[1])*(b[1]-a[1]))

def getTriangleArea(pt0, p1, p2):
    a = distance(pt1, pt0)
    b = distance(pt2, pt0)
    c = distance(pt2, pt1)
    s = (a+b+c) * 0.5
    return sqrt(s * (s-a) * (s-b) * (s-c))

Now notice how this only computes the triangle area in the XY plane. This works simply because the 2D projection of the triangle area is all we need. The height is then defined by the triangle’s center Z.

def getTriangleHeight(pt0, pt1, pt2):
    return (pt0[2] + pt1[2] + pt2[2]) * 0.33333333

Consider any triangle, extrude it down to the floor, and see that this works for any prism defined along the Z axis this way.

A rotated triangle’s area in the XY plane is smaller than the actual area, but by using the face-center the volume will remain accurate.

prisms

Now these prisms have the same volume. The trick is to consider every triangle as such a prism, so call getTriangleVolume on each triangle. The last problem is negative space, for this we compute the normal. I use maya’s normals, so the volume is negative if all normals are inversed, but you can compute them all the same.

def getTriangleVolume(pt0, pt1, pt2):
    area = getTriangleArea(pt0, pt1, pt2) * getTriangleHeight(pt0, pt1, pt2)
    # this is an optimized 2D cross product
    sign = (pt1[0]-pt0[0]) * (pt2[1]-pt0[1]) - (pt1[1]-pt0[1]) * (pt2[0]-pt0[0])
    if not sign:
        return 0
    if sign < 0:
        return -area
    return area

prisms2

The selected wireframe shows the prism defined by the bottom triangle, because the normal.z points downwards it will become negative volume. So adding the initial prism volume and this prism volume will give the accurate volume of this cut off prism. Now consider this:

prisms3

To avoid confusion I placed the object above the grid; but below the grid a negative normal * a negative height will still add volumes appropriately.

So that's it.

from math import sqrt
from maya.OpenMaya import MItMeshPolygon, MDagPath, MSelectionList, MPointArray, MIntArray


def distance(a, b):  
    return sqrt((b[0]-a[0])*(b[0]-a[0]) +   
      (b[1]-a[1])*(b[1]-a[1]))  
      
def getTriangleArea(pt0, pt1, pt2):  
    a = distance(pt1, pt0)  
    b = distance(pt2, pt0)  
    c = distance(pt2, pt1)  
    s = (a+b+c) * 0.5  
    return sqrt(s * (s-a) * (s-b) * (s-c))  
    
def getTriangleHeight(pt0, pt1, pt2):  
    return (pt0[2] + pt1[2] + pt2[2]) * 0.33333333  

def getTriangleVolume(pt0, pt1, pt2):  
    area = getTriangleArea(pt0, pt1, pt2) * getTriangleHeight(pt0, pt1, pt2)  
    # this is an optimized 2D cross product  
    sign = (pt1[0]-pt0[0]) * (pt2[1]-pt0[1]) - (pt1[1]-pt0[1]) * (pt2[0]-pt0[0])  
    if not sign:  
        return 0  
    if sign < 0:  
        return -area  
    return area

def getPolygonVolume(shapePathName):
	volume = 0
	li = MSelectionList()
	li.add(shapePathName)
	path = MDagPath()
	li.getDagPath(0, path)
	iter = MItMeshPolygon(path)
	while not iter.isDone():
		points = MPointArray()
		iter.getTriangles(points, MIntArray())
		for i in range(0, points.length(), 3):
			volume += getTriangleVolume(points[i], points[i+1], points[i+2])
		iter.next()
	return volume

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>