Calculation of centroid & volume of a polyhedron when the vertices are given

Question

Given the location of vertices of a convex polyhedron (3D), I need to calculate the centroid and volume of the polyhedron. Following code is available at Mathworks site.

function C = centroid(P)
k=convhulln(P);
if length(unique(k(:)))<size(P,1)
    error('Polyhedron is not convex.');
end
T = delaunayn(P);
n = size(T,1);
W = zeros(n,1);
C=0;
for m = 1:n
    sp = P(T(m,:),:);
    [null,W(m)]=convhulln(sp);
    C = C + W(m) * mean(sp);
end
C=C./sum(W);
return
end

The code is elegant but is terribly slow. I need to calculate the volume and centroid of thousands of polyhedrons hundreds of times. Using this code in its current state is not feasible. Does anyone know a better approach or can this code be made faster? There are some minor changes I can think of such as, replacing mean with expression for mean.

Answer 1

There is a much simpler approach to calculate the volume with minimal effort. The first flavour uses 3 local topological information sets of the polyhedron, the tangent unit vector of the edges, the unit vectors of the in-plane normal on this tangent and the unit vector of the facet itself (which are very simple to extract from the vertices). Please refer to Volume of a Polyhedron for further details.

The second flavour uses the face areas, the normal vectors and the face barycenters to calculate the polyhedron volume according to this Wikipedia Article. Both algorithms are rather simple and very easy to implement and through the simple summing structure easy to vectorize too. I suppose that both approaches will be much faster than doing a fully fledged tesselation of the polyhedron.

The centroid of the polyhedron can then be calculated by applying the divergence theorem transferring the integration over the full polyhedron volume into an integration over the polyhedron surface. A detailed description can be found in Calculating the volume and centroid of a polyhedron in 3d. I did not check if the tesselation of the polyhedron into triangles is really necessary or one could work with the more complex polygon surfaces of the polyhedron too, but in any case the area tessellation of the faces is much simpler than the volume tesselation. In total such a combined approach should be much faster than the volume approach.

Answer 2

Thinking your only option if quickhull isn't good enough is cudahull if you want exact solutions. Although, even then you're only going to get about a 40x increase max (it seems).

I'm assuming that your the convex hulls you make each have at least 10 vertices (if it's much less than that, there isn't much you can do). If you don't mind "close enough" solutions. You could create a version of quickhull that limits the number the of vertices per polygon. The number of vertices you limit the calculation to will also allow for calculation of maximum error if needed.

The thing is that as the number of vertices on the convex hull approach infinity, you end up with a sphere. This means due to the way quick hull works, each additional vertex you add to the convex hull has less of an effect* than the ones before it.

*Depending on how quickhull is coded, this may only be true in a general sense. Making this true in practice would require modifying quickhull's recursion algorithm, so while the "next vertex" is always calculated (except for after the last vertex is added, or no points remain for that section), vertices are actually added to the convex hull in the order that maximizes the increase to the polyhedrons volume (likely by order of most distant to least distant). You'll incur some performance cost for keeping track of the order to add vertex's but as long as the ratio of pending convex hull points vs pending points is high enough, it should be worth it. As for error, the best option is probably to stop the algorithm either when the actual convex hull is reached, or the max increase to volume gets smaller than a certain fraction of the current total volume. If performance is more important, then simply limit the number of convex hull points per polygon.

You could also look at the various approximate convex hull algorithms, but the method I outlined above should work well for volume/centroid approximation with ability to determine error.

Answer 3

How much you can speed up your code depends a lot on how you want to calculate your centroid. See this answer about centroid calculation for your options. It turns out that if you need the centroid of the solid polyhedron, you're basically out of luck. If, however, only the vertices of the polyhedron have weight, then you could simply write

[k,volume] = convhulln(P);
centroid = mean(P(k,:));