From: boyan@xafs.phys.nd.edu (Boyan I. Boyanov) Newsgroups: sci.math Subject: Volume of a tetrahedron Date: 27 Mar 1995 00:30:42 GMT Hi, I have a stupid question that I cannot solve myself and cannot find the answer to in any of the math handbooks I've looked at. Is there a general formula for the volume of an ARBITRARY polyhedron? I know everything there is to know about the polyhedron - the coordinates of the vertices, the equations of the planes that form the walls, the number of walls, the number of sides, the number of vertices, etc. About the only restriction that is put on the polyherdon is that it is "convex", i.e. all vertices lie in a single "half-space" relative to any wall. I was trying to do it by breaking up the polyhedron into a collection of tetrahedra and summing up the volumes of the tetrahedra, but I can't seem to come up with an algorithm that allows me to determine which vertices form the tetrahedra that make up the polyhedron. Any ideas? Is there a better nad easier way to do this? Please Cc: me if you post a followup, if at all possible. Thanks! Boyan -- Boyan I. Boyanov boyan@tmnxt1.iit.edu http://tmnxt1.iit.edu/~boyan/ ============================================================================== From: rusin@washington.math.niu.edu (Dave Rusin) Newsgroups: sci.math Subject: Re: Volume of a tetrahedron Date: 27 Mar 1995 22:52:56 GMT In article <3l50vi$1i0@news.nd.edu>, Boyan I. Boyanov wrote: >Is there a general formula for the volume of an ARBITRARY polyhedron? Well, I suppose I ought to answer this since I just advertised Green's theorem for a similar purpose. In this question, the author wants a space integral integral_P 1 dxdydz where P is the polyhedron. The corresponding technique is to use Stokes' Theorem, whose general form _ _ _/ w = _/ dw boundary(S) S (for w any differential form on S) specializes for 2-forms to _ _ _/ P dxdy + Q dydz + R dzdx = _/ (dP/dz + dQ/dx + dR/dy) dxdydz so we can get the space integral of 1, for example, by integrating the 2-form z dxdy over the surface of the polyhedron. We can compute the surface integrals face-by-face. The quick answer here is that the individual surface integrals, when divided by the area of the face, would compute the average value of the z coordinate on the face. Therefore, we can compute the surface integrals easily by averaging the z coordinates of the vertices of the face, and multiplying by an area: the surface integral int( zdxdy ) over a face spanned by (x0 y0 z0), (x1 y1 z1), and (x2 y2 z2) is (z0+z1+z2)/3 . | (1/2)( y1.(x0-x2) + y2.(x1-x0) + y0.(x2-x1) ) | This surface integral is added to the total if the three points are in counterclockwise order as seen from the outside of the polyhedron, subtracted otherwise. A proof is appended as a postscript to this post; it may perhaps reinforce an understanding of surface integrals. It perhaps bears mentioning that this procedure does not assume the polyhedron is convex, nor indeed homeomorphic to a sphere. It does assume that the polyhedron is closed (no boundary) and a triangularization of a manifold (which will then be compact and orientable, conditions which I would otherwise have to add). It is not necessary that all the faces be triangles, really, although any other polygons may be so divided. The method may also be extended to other surface integrals besides volume, and so for example may be used to compute the center of mass. One can remove the absolute values and simply add in the polynomial expression if an orientation is chosen consistent with the "right-hand-rule" for outward pointing normals. As in the previous post I will concede that there is some redundant calculation here, so that if it is important that the computation be done quickly it pays to watch for repeated line integrals and so on. dave __Derivation of formula for surface integral over one face__ Assuming that the face is not vertical, we can write z = Ax+By+C for some constants A, B and C (easily computed from any three non-collinear vertices on the face). Then that face contributes to the surface integral the sum A integral(x) + B integral(y) + C integral(1) where the integral may be taken over the projection of the face to the xy-plane (translation: ignore the z coordinates on the face). There is an orientation problem: the integral over the projection to the xy-plane will be off by a sign if the polyhedron faces "in" rather than "out", so you will need to know whether the interior of the polyhedron contains those points whose z coordinates are just larger than Ax+By+C or just smaller. In a previous post, I discussed integrating these very three surface integrals around a polygon. (These were done with Green's theorem.) The upshot is that these integrals may be computed as simple sums around the perimeter of the polygon. The segment from (x0, y0) to (x1, y1) contributes (x0 + x1)/2 . (y1-y0) to the integral of 1 over the polygon, (x0^2 + x0 x1 + x1^2)/6 . (y1-y0) to the integral of x, and ( (y0 x1 + x0 y1) + 2(x1 y1 + x0 y0) )/6 . (y1 - y0) to the integral of y. So, for example, if one face of your polyhedron is a triangle with vertices (x0,y0,z0), (x1,y1,z1), and (x2,y2,z2), then the equation for the plane is z = Ax + By + C where A, B, and C are determined by ( A ) ( x0 y0 1 )^{-1} ( z0 ) ( B ) = ( x1 y1 1 ) ( z1 ) ( C ) ( x2 y2 1 ) ( z2 ) We compute the line integrals around the triangle in the xy-plane: the integral of 1 is I1 = (1/2)( y1.(x0-x2) + y2.(x1-x0) + y0.(x2-x1) ) as in a followup I made to the previous post; the integral of x is this times (x0+x1+x2)/3, and the integral of y is the same factor I1 times (y0+y1+y2)/3. (All these should have their sign adjusted in case the points (x0,y0), (x1,y1), (x2,y2) do not wrap counterclockwise around the triangle. The correction is easy: since integral(1)=Area > 0, the sign correction to multiply by is just the sign S1 of I1. Thus, the surface integral over this face is the linear combination [ A(x0+x1+x2)/3 + B(y0+y1+y2)/3 + C ] . (S1 I1) As it happens, I1 is precisely the determinant of the matrix to be inverted above, so multiplyng thru we see that A B and C may be replaced by the cofactors, making division unnecessary. Apart from the sign correction, the previous surface integral may thus be written (z0.(y1-y2)+z1.(y2-y0)+z2.(y0-y1)) . (x0+x1+x2)/3 + (z0.(x2-x1)+z1.(x0-x2)+z2.(x1-x0)) . (y0+y1+y2)/3 + (z0.(x1y2-x2y1)+z1.(x2y0-x0y2)+z2.(x0y1-x1y0)) This, in turn, simplifies to just (z0+z1+z2)/3 . (I1) (again, times the sign correction so that S1 I1 > 0 ) This is the formula given above.