From: andersoc@execpc.com (Chris Anderson)
Newsgroups: sci.image.processing,sci.math
Subject: volume of an irregular polyhedron? (how to find)
Date: Sun, 06 Aug 95 07:33:20 GMT
Summary: how can i find teh volume of a part in an .STL file?
Keywords: polyhedron stl volume polygon


I need to find a method for determination of the volume of a clsoed 
tesilated object in a CAD file. I have the x y and z values of each 
point of each of the triangles that represent the surface of the part, 
and the unit normal for each triangle which tells which side of the 
triangle faces the inside of the object. 

The best I'v come up with is maybe slicing the part in the z axis and
finding the area of the oplygons formed. Then multiplying that by the
slice height to aproximate the volume. Isn't there a more elegant way
to use simply the vertices points and normals?

Could anyone sugjest any good books on manipulating tesalated 3D 
objects? I've found true surface area of the object by adding all the 
areas of the triangles. I'll need to still find volume, and projected 
vertical area as well.

Please email replies if posible, but I'll watch the group of course.

Thanks

	Chris{}Anderson
	andersoc@execpc.com or andersonc@msoe.edu
	MSOE Rapid Prototyping Center
	414-277-7285
==============================================================================

From: rusin@washington.math.niu.edu (Dave Rusin)
Newsgroups: sci.image.processing,sci.math
Subject: Re: volume of an irregular polyhedron? (how to find)
Date: 11 Aug 1995 17:41:57 GMT

In article <401r80$9q4_001@charon.execpc.com>,
Chris Anderson <andersoc@execpc.com> wrote:
>
>I need to find a method for determination of the volume of a clsoed 
>tesilated object in a CAD file. I have the x y and z values of each 
>point of each of the triangles that represent the surface of the part, 
>and the unit normal for each triangle which tells which side of the 
>triangle faces the inside of the object. 

I mailed a version of this to the original poster, but since I haven't
seen complete answers posted I thought I'd do so. This is an edited
version of a post I made some time ago on the same topic. It gives a
method requiring simple computations involving the coordinates of the points.

The poster wants the volume of a region. I will give a geometric
derivation of the answer at the end, but I know some feel uneasy when
I present it, thinking that the diagram is possibly glossing over some
special cases. So I will first view the volume as a space integral   
integral_P  1  dxdydz  where  P  is the polyhedron. This can be
evaluated with 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.

In doing so I will assume that the boundary is presented as a collection
of triangles, oriented in such a way that the boundaries cancel out.
For example, a tetrahedron with vertices  A B C D  has a boundary consisting
of triangles   A B C,  C B D,  D B A,  and  A C D; you can permute the vertices
in any triangle cyclically, but you cannot replace the last, for example,
with the ordering  A D C. For a triangulated surface in  R^3, it is
possible to choose a consistent orientation; we may, for example,
always arrange it so that the points are taken in counterclockwise order
as seen from the outside of the figure. (If the orientations are consistent
but not counterclockwise, the computation below will be off by a sign, so
it's safe to choose any consistent orientation, perform the indicated
computations, and then take absolute values.)

Now let us compute the surface integral on one face. One can parameterize
the face with the triangle tracing from (0,0) to (1,0) to (0,1) by sending
the point  (u,v) inside this triangle to the point
  P0 (1-u-v) + P1 (u) + P2 (v) =
  (x0(1-u-v) + x1 u + x2 v, y0(1-u-v) + y1 u + y2 v, z0(1-u-v) + z1 u + z2 v)
So  z  is the third coordinate,  dx = (x1-x0) du + (x2-x0) dv, and
dy = (y1-y0) du + (y2-y0) dv. To compute  dx dy  you need to remember that
there is an implicit wedge product sign  ( /\ ) in here, as in the 
discussion of antisymmetric products; one finds
  dx dy = [(x1-x0)(y2-y0) - (x2-x0)(y1-y0)]  du dv
so that  z dx dy  is easily integrated over the triangle as a certain
linear combination of the integrals of   1 du dv  (the area, 1/2),
u du dv, and  v du dv. Use Fubini's theorem, or Green's theorem, to
evaluate the last two integrals; you get  1/6  in either case.

So the surface integral of  z dx dy  over this (oriented) triangle is
	[(x1-x0)(y2-y0) - (x2-x0)(y1-y0)]   times
(1/2)(z0) + (1/6)(z1-z0) + (1/6)(z2-z0)  = (1/6)(z0+z1+z2)
(the first factor may also be written a little more symmetrically as
	[ x0(y1-y2) + x1(y2-y0) + x2(y0-y1) ]. )

This then is the conclusion: the volume of the region is found by summing,
over all triangles in the surface, the product
	(1/6)(z0+z1+z2) [ x0(y1-y2) + x1(y2-y0) + x2(y0-y1) ]
where the three vertices of each triangle are taken in counterclockwise
order as seen from the outside of the figure.

It perhaps bears mentioning that this procedure does not assume the
polyhedron is convex, nor indeed homeomorphic to a sphere. The boundary
can even be disconnected (e.g. the shape can be a ball with a hole in
its core) as long as the orientations are done consistently. 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. 

There is some redundancy here: each edge is part of precisely two surface
triangles, taken with opposite orientation in the two cases. So for
example, the term  (1/6) z0 x0 y1  which occurs in one triangle  (P0, P1, P2)
will equal the term  (1/6) z0' x0' y2'  in the adjacent triangle 
(P0'=P0, P1', P2'=P1). These occur with opposite sign and so cancel.
What remains then are only those terms involving one coordinate from
each of the three points of each triangle. Therefore it really
suffices to add only
    (1/6)[ x0 y1 z2 - x0 y2 z1 + x1 y2 z0 - x1 y0 z2 + x2 y0 z1 - x2 y0 z2 ]
for each triangle.

Now, where have we seen this before? Ah, yes: this is just 
             (  x0   x1   x2  )
   (1/6) det (  y0   y1   y2  )
             (  z0   z1   z2  )
One interpretation of this (1/6) of a determinant is that this is the
signed volume of the region spanned by these three vertices and the
origin; the sign is positive or negative according to a sort of right-hand
rule. If we consider the effect of adding several of these together, we
get an alternative proof of the accuracy of the proposed formula. It
may be easier to visualize if you imagine a shape which does not include
the origin. Then for each triangle on the surface, extend linearly to
make a tetrahedron which includes the origin as the fourth vertex.
Add the signed area of all these regions. The triangles on the furthest
side of the figure contribue sums which total the volume of the figure
plus the volume of the region between the figure and the origin; the
triangles on the nearer side of the figure subtract away this excess.
(The easiest figure to see this with is perhaps the tetrahedron with
vertices at (4,0,0), (0,4,0), (0,0,4), and (1,1,1).)

dave (Stokes' Theorem  fan)