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.