Calculating the exact volume of a polygon-based 3d object
I'm not sure where you could need the exact volume of an object, but however here you'll learn how to calculate it quite easily and fast.
Calculating the exact area of a vector-based 2d object
Lets start from the easier thing. Here's an object:
You could triangulate the object:
And then calculate the area of the triangles, but we're not gonna do this, because it is harder to "triangulate" in 3d.
Instead, we divide the object with vertical lines at each end vertex point.
Now we can easily calculate the area A between each of the vertical lines, because the width of the polygon changes linearly as the vertical line moves to right:
A = (a + b) h / 2
So you'll only need to know how long are the both ends of each segment. This is very easy to calculate: just check which edges the vertical line crosses and at what points, and do substraction between adjacent intersection points.
The real deal
The 3d case isn't as intuitive as the 2d one is, as was demonstrated by Jan van Kranendonk when he showed that my previous algorithm was utterly wrong. But this new one should work. And it comes with a proof, as relying on intuition proved to be a bad practice.
The 3d version is quite similar to the 2d one. In 2d case we took a line, kept its orientation and segmented the 2d polygon into pieces along its vertices, where the width of the polygon, as a function of the line's position, changed linearly. In 3d, we take a plane, keep its orientation at all times and segment the 3d object into pieces along its vertices where the area of the object, as a function of the plane's position, changes as a second degree polynom. Now that's a mouthful, so it's best to start with an example, the needed formula, and then give the proof for those interested.
First, lets take a simple shape for demonstration - a slightly tilted pyramid:
Now, just like the 2d object was divided among its vertices, the same can be done here. The splitting is presented below with multiple pictures for clarity.
Level 0
Level 1
Level 2
Level 3
The images show how the plane moves upwards, keeping its orientation. In each of the pictures it is at the same level with some vertex of the 3d object. In the images I've marked the distances between adjacent levels, calculated along the plane's normal and denoted as h1, h2 and h3. I've also drawn the intersection areas between the plane and the object, A1 and A2. These can be calculated by checking the points where the plane hits the object's edges, and using the 2d technique described at the top of this document to calculate the area between these points.
What I said earlier was that as the plane moves from (e.g.) level 1 to level 2, the area changes as a second degree polynom. That is, using the picture as an aid, we have a function A:
| A(x) = a x2 + b x + c, where x
in range [0, h2] A(0) = A1 A(h2) = A2 |
To solve a function with three unknowns, a, b and c, knowing the function's value only for two points, 0 and h2 isn't enough. A third point needs to be introduced, which differentiates the 3d method from the 2d one.
Level 1 and a half
Here the plane has been placed exactly half-way between level 1 and level 2. The area here, Ah, is once again calculated by finding where the plane intersects the object's edges and by using the 2d method for this 2d shape. The equation is finally complete:
| A(x) = a x2 + b x + c, where x
in range [0, h2] A(0) = A1 A(h2 / 2) = Ah A(h2) = A2 |
Once A(x) has been solved, it can be integrated from 0 to h2, giving the result of volume between levels 1 and 2:
| V = (A1 + A2 + 4 Ah) h2 / 6 |
Similarly you can calculate volume between all adjacent plane levels and then sum them up to get the volume of the whole object.
Proof
As the edges of a 3d mesh are line segments, it has to be shown that if a 2d polygon's vertices move linearly, pi = xi + vit for each vertex i, its area will change according to a second degree polynom.
A polygonal shape can be triangulated so that it consists of a set of triangles. Its area can then be calculated as a sum of the areas of the individual triangles. As the vertices in the polygon move, some of the triangles will "warp", i.e. the order of its vertices as counted clockwise will change. At these critical situations, the polygonal shape has to be triangulated again so that warping will not happen. This way we can avoid any special cases. Thus, all we need to show now is that the area of each of these triangles is a second degree polynom as a function to t, which is the parameter by which the vertices move.
A triangle's three points can be represented by position vectors a, b and c (corresponding to some vertices pi). Their linear displacements will be denoted as these first degree polynoms:
| a' = a
+ n1t b' = b + n2t c' = c + n3t |
Then the area of a triangle can be calculated with a 2-dimensional cross product,
| A = 1/2 | e | | f | sin g = 1/2 | e x f | = 1/2 | ex fy - ey fx | |
where e and f are the triangle's sides as vectors, g is the angle between them and ex etc. are the vectors' components in cartesian coordinates. We want to calculate the area of a linearly displaced triangle, so the side vectors are e.g. a' - b' and c' - b'. Taken the definition for these vectors and using the formula for A, it can be seen that 1/2 | (a' - b') x (c' - b') | yields a second order polynom. Without the boring calculation, it suffices to note that we're just substracting first order polynoms, yielding first order polynoms at most and multiplying first order polynoms with each other, yielding second order polynoms at most.
Solving V
General solution for the set of equations
| A(x) = a x2 + b x + c A(0) = A1 A(d / 2) = Ah A(d) = A2 |
is
| a = (2 A1 + 2 A2 - 4 Ah) / d2 b = -(3 A1 + A2 - 4 Ah) / d c = A1 |
Volume V is calculated as
| V | = | Integrate 0 -> d [A(x)dx] |
| = | Integrate 0 -> d [(a x2 + b x + c) dx] | |
| = | 0 -> d [(a x3 / 3 + b x2 / 2 + c x)] | |
| = | (a d3 / 3 + b d2 / 2 + c d) | |
| = |
(A1 + A2 + 4 Ah) d / 6 |
If you have questions or have spotted an error, contact me.