Algorithm for finding the volume of a convex polytope

Question

It's easy to find the area of a convex polygon by division into triangles, but what is the optimal way of finding the volume of higher-dimensional convex bodies? I tried a few methods for dividing them into simplices, but gave it up and went with a Monte Carlo estimation scheme instead. Bonus question: How to find the surface area of those same convex bodies?

EDIT: To answer David's question: the data set is a Voronoi tessellation of an n-dimensional volume (n usually 4) with a periodic boundary (like a torus). So I have the coordinates of the vertices of the convex bodies as well as the connectivity of all the facets, faces, etc. For the Monte Carlo I mentioned, I did convert everything to half-spaces, so I think that was not very difficult.

Answer 1

As a follow-up to Barton's response: for hardness results, see I. Bárány & Z. Füredi, Computing the volume is difficult, Discrete and Computational Geometry, 1987. But there are polynomial time approximation schemes for volume of convex bodies independent of dimension, based on random walks within the body: see e.g. M. Dyer, A. Frieze, & R. Kannan, A random polynomial-time algorithm for approximating the volume of convex bodies, J. ACM 1991 and R. Kannan, L. Lovász, & M. Simonovits, Random walks and an O*(n^5) volume algorithm for convex bodies, Random structures and algorithms, 1997.

Before getting to these answers, though, it's important to ask: how is your input represented? Is it a convex hull of a set of points, and all you know is the points? Is it an intersection of halfspaces? Is it given to you as an entire face lattice? Also, what range of dimensions do you care about The upper and lower bounds above are quite general but the dimension and input representation may still make a difference to the answer.

Answer 2

I think this problem is hard--the known algorithms are both slow and nontrivial to implement. See Exact Volume Computation for Polytopes for a survey. An interesting feature is that there are various algorithms which are well suited for different kinds of polytopes.

As a practical answer, Qhull can compute volumes and surface areas.

Answer 3

Matthias Beck and Dennis Pixton used almost 17 GigaHertz years to compute the value of the 10-dimensional Birkhoff polytope. Read about it, and find the answer, here.

If you can get the same result quicker, I'm sure they'd be delighted to know how you did it.

Answer 4

There is a randomized algorithm with O(n^4) running time based on the simulated annealing using the hit-and-run convex body sampling algorithm along with isotropy approximation by Vempala and Lovasz: here. The running time is the number of oracle calls for membership in the convex body and has an implicit polylogarithmic term.

Answer 5

Here's a fairly straightforward solution for polyhedra (3 dimensions), with running time O(v+ve), where v is the number of vertices and e is the number of edges. I suppose it could be extended to higher dimensions, but it would probably have much worse running time (I fear roughly exponential as in O(vⁿ), where n is the number of dimensions).

Let our polyhedron have n vertices, defined by their x,y,z coordinates: v₁, v₂, ..., v_n and let the lowest point be v₁ and the origin (modify the values for the others accordingly), and let it have e edges, defined by the vertices which they connect. Then, since we have coordinates for the vertices (as that is how we defined them), there must be a "ground level" plane p₀ running through the x and z axes (the y-axis being height, and the ground never having an elevation). Then, let v₂ be the point closest to the ground plane (shortest line perpendicular to the plane), and let v₃ be the next closest, etc, through v_n.

Through each of the points v₂ through v_n, draw a plane perpendicular to the ground, and let them be numbered p_m, where m is the subscript of the vertex through which it was drawn. Then, the volume of our polyhedron is equal to the sum of the volumes of the figures between the planes. We should have something resembling this:

Let the heights between the segments be h₁ through h_n-1, where height h_j is the height between planes p_j and p_j+1.

Now, through each plane, we have a polygon (or more, if the figure is concave), whose vertices' coordinates can be calculated easily as follows:
Let the edge that runs through the plane p_j have endpoints v_a and v_b. Then, the displacement vector is v_b - v_a (assuming the coordinates of v are in vector-form), and the percentage travelled up is $\frac{h_j-h_a}{h_{b-1}-h_a}$. Multiply this by v_b - v_a and add to v_a to calculate the new point of intersection for that edge:
Intersection point = $(v_b-v_a)\frac{h_j-h_a}{h_{b-1}-h_a}+v_a$
The area of these polygons can be determined using triangles, or a simplification of this very process in just 2 dimensions.

PlanetMath says that the volume of a prismatoid (which is the type of figure contained between sequential planes) is $h\frac{B_1 + B_2 + 4M}{6}$, where the Bs are the areas of the parallel polygons and M is the area of the midway polygon, which is exactly halfway between them (and parallel to them). Since we already know the area of each of the end polygons, and we can easily calculate the vertices of the midway polygon (using the previous paragraph's method), we can calculate the volume of the resulting prismatoids. Adding them up yields the total volume of the polyhedron.

I suppose that the only real issue in this case, then, is, via code, determining which edges run through any particular plane, but if we were to actually look at it, we could tell very easily.

A simpler version of this can be used to figure out the area of any polygon; simply draw lines through the vertices parallel to the x-axis and calculate the area of the resulting trapezoids as (b₁+b₂)/2

Answer 6

An implementation of randomized algorithms using hit-and-run and comparison with exact software presented in Bueler, Enge and Fukuda is presented here. It works efficiently in high dimensions (say few hudrend) but it only computes an estimation, not the exact volume. The software is available on github.

Answer 7

There is a result by J. Borwein et al on volumes of convex polytopes constructed by intersecting symmetric hyperplanes. These volumes are obtained from multivariate integrals involving the function sin(x)/x. These integrals have very nice closed forms. Here are two references on this work:

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.40.8186

and

http://algo.inria.fr/seminars/sem01-02/borwein1.pdf

Answer 8

As documented above the paper by Bueler, Enge and Fukuda is a good place to start and give many techniques for computing volumes of polytopes. In my experience using lrs to build the triangulation is a very good method. If your polytope has symmetries, then it is possible to use them to reduce the size of the computation. See here for the relevant programs.

Answer 9

Since you're given so much information about the faces of the polytope, I think that the division into simplices method ought to be straightforward (if tedious) to carry out. The idea is to divide the polytope into orthoschemes, with signed volumes, as described below. This has the advantage that the formula for the volume of an orthoscheme is very simple, and computing the coordinates of the orthoschemes is simple linear algebra. The disadvantage is that there will be one orthoscheme for each flag of the polytope, so the number of them could be quite large. Computing the surface area of the boundary could be done in tandem, since one would just take the area of each orthoschemes intersection with the supporting hyperplanes. I suspect this method won't work for you though, if the polytope is complicated (and maybe this is the approach you already tried).

To compute the vertices of the orthoschemes, first choose a point, then take its orthogonal projection to each supporting hyperplane of the polytope. Inductively do this for each facet. The for each flag, take the corresponding points in each face of the flag, and form an orthoscheme. The sign of the volume of the orthoscheme will be determined in each dimension by whether the vertex lies inside or outside the corresponding hyperplane times the sign of the lower dimensional one it is a cone on. Some of the orthoschemes will lie partly outside of the polytope, but the volumes outside will cancel with this sign convention.

Answer 10

You could use a finite element method. If you stuff the surface with tetrahedra http://www.cs.berkeley.edu/~jrs/papers/stuffing.pdf‎ (pdf link not allowed), you could then sum the volumes of those tetrahedra.

Answer 11

This is a very old topic, but I think there is a very simple way that can be generalized to polytopes easily. The method handles even non-convex cases easily.

I do think this might be what Ian Agol has described. I'm being a bit more explicit here how it may be done.

For 2-polytopes (polygons), the '2-volume' (area) can be calculated as follows: Given n vertices (x_i, y_i) in order against the clock around the polygon, where vertex nr. n is understood to be connected to vertex nr. 1. These define n edges, where edge i is the edge nr. i is the edge from vertex i to vertex i+1 looping back to vertex 1 for the last edge. Define the origin to be the point (0,0). Then each vertex nr i. together with the origin, will span a triangle with area 1/2 * (x_i * y_(i+1) - x_(i+1) * y_i) Where it is understood that x_(n+1) = x_1 and the same for y. (In fact any origin may be chosen which would yield the same result. Area of the polygon is of course translation invariant. The point (0,0) is the most efficient computationally.)

To get the area you sum this quantity for all vertices.

When the origin is inside the polygon, it is clear this yields the correct result. When origin lies outside the polygon, some triangle area's will get opposite sign, in such a manner that they exactly cancel to yield to correct total area.

Generalising to a 3-polytope calculating the 3-volume would be like the following: Given the n faces of the 3-polytope, calculate the oriented 3-volume of the face by adding the origin, to each face and calculating the volume spanned by face and origin. Then sum all oriented 3-volumes of all faces, that would be the total volume. In the case where you have a face that is a triangle (x1,y1,z1),(x2,y2,z2),(x3,y3,z3), the spanned volume would be

1/2 * det |(x1 x2 x3), (y1 y2 y3), (z1 z2 z3)|

If the face has more than 3 vertices, you could go around the face keeping the first two vertices fixed fixed, and summing such a half-determinant for each of the non fixed vertices of the face. Orientation (sign) of the individual determinants will automatically handle non convex cases correctly.

Note that the vertex order of each of the faces has to be oriented the same with respect to the inside of the 3-polytope. So looking from inside the polytope to outside, all face vertices should go e.g. anti-clockwise.

I suspect generalisation to an n-polytope should be straightforward.