Volume of convex hull with QHull from SciPy

Question

I'm trying to get the volume of the convex hull of a set of points using the SciPy wrapper for QHull.

According to the documentation of QHull, I should be passing the "FA" option to get the total surface area and volume.

Here is what I get.. What am I doing wrong?

> pts
     [(494.0, 95.0, 0.0), (494.0, 95.0, 1.0) ... (494.0, 100.0, 4.0), (494.0, 100.0, 5.0)]


> hull = spatial.ConvexHull(pts, qhull_options="FA")

> dir(hull)

     ['__class__', '__del__', '__delattr__', '__dict__', '__doc__', '__format__', '__getattribute__', '__hash__', '__init__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_qhull', '_update', 'add_points', 'close', 'coplanar', 'equations', 'max_bound', 'min_bound', 'ndim', 'neighbors', 'npoints', 'nsimplex', 'points', 'simplices']

 > dir(hull._qhull)
     ['__class__', '__delattr__', '__doc__', '__format__', '__getattribute__', '__hash__', '__init__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__']

Try to update your question with a real question (not "here's what I get").. It took me a while to figure out that nowhere the total area and volume can be found, despite the fact that you supplied the proper option. — Hugues Fontenelle, Jul 14 '14 at 09:46
My wild guess is that SciPy does not wrap that particular option flag. — Hugues Fontenelle, Jul 14 '14 at 09:48
The hard way is to implement it: http://wiki.scipy.org/Cookbook/Finding_Convex_Hull — Hugues Fontenelle, Jul 14 '14 at 09:54
One thing that would help would be complete `pts`. That way we could try it ourselves. — DrV, Jul 14 '14 at 10:30
It's not implemented in Scipy Qhull wrappers. It could be easily added, if there's a need. — pv., Jul 14 '14 at 12:11

Jaime · Accepted Answer · 2014-07-15T20:46:56.540

15

There does not seem to be any obvious way of directly getting the results you are after, regardless of what parameters you pass in. It shouldn't be too hard to compute yourself if, instead of ConvexHull, you use Delaunay (which also provides most of the convex hull related info).

def tetrahedron_volume(a, b, c, d):
    return np.abs(np.einsum('ij,ij->i', a-d, np.cross(b-d, c-d))) / 6

from scipy.spatial import Delaunay

pts = np.random.rand(10, 3)
dt = Delaunay(pts)
tets = dt.points[dt.simplices]
vol = np.sum(tetrahedron_volume(tets[:, 0], tets[:, 1], 
                                tets[:, 2], tets[:, 3]))

EDIT As per the comments, the following are faster ways of obtaining the convex hull volume:

def convex_hull_volume(pts):
    ch = ConvexHull(pts)
    dt = Delaunay(pts[ch.vertices])
    tets = dt.points[dt.simplices]
    return np.sum(tetrahedron_volume(tets[:, 0], tets[:, 1],
                                     tets[:, 2], tets[:, 3]))

def convex_hull_volume_bis(pts):
    ch = ConvexHull(pts)

    simplices = np.column_stack((np.repeat(ch.vertices[0], ch.nsimplex),
                                 ch.simplices))
    tets = ch.points[simplices]
    return np.sum(tetrahedron_volume(tets[:, 0], tets[:, 1],
                                     tets[:, 2], tets[:, 3]))

With some made up data, the second method seems to be about 2x faster, and numerical accuracy seems very good (15 decimal places!!!) although there has to be some much more pathological cases:

pts = np.random.rand(1000, 3)

In [26]: convex_hull_volume(pts)
Out[26]: 0.93522518081853867

In [27]: convex_hull_volume_bis(pts)
Out[27]: 0.93522518081853845

In [28]: %timeit convex_hull_volume(pts)
1000 loops, best of 3: 2.08 ms per loop

In [29]: %timeit convex_hull_volume_bis(pts)
1000 loops, best of 3: 1.08 ms per loop

edited Jul 15 '14 at 20:46

answered Jul 14 '14 at 10:34

Jaime

65,696
17
124
159

This is great. I was looking for the volume information from qhull to save computational time, but I don't think it can be done in Python any faster than your version. So just to confirm, Delaunay does the tessellation of the convex hull, correct? Also, neat little thing this `einsum` is... Did not know that. :) – Diana Jul 15 '14 at 04:03
1

Delaunay turned out to be very slow, so instead I used ConvexHull, got the extremity points (the ones on the hull), and ran Delaunay on those points only. It is 1-2 orders of magnitude faster this way on my data. – Diana Jul 15 '14 at 08:09
That's smart indeed... It may not improve much further, but you may want to try skipping the call to `Delaunay` altogether, and build a triangulation of your convex hull by choosing a point on the hull, then computing the volume of all tetrahedra that contain that point and the points on each of the convex hull's simplicial facets (i.e. the `.simplices` attribute of `ConvexHull` objects). It will produce more elongated tetrahedra, so it may be less numerically stable. But it has to be faster. – Jaime Jul 15 '14 at 11:55
I tried that too, but the difference in results was too big. Maybe I did something wrong. Numerical instability can't explain two orders of magnitude difference in the volume value, given that it's already convex hull, right? I'll probably take a look at this again. – Diana Jul 15 '14 at 19:35
1

Take a look at the latest edit, for random data I get a 2x speed-up and accuracy to 15 decimal places... – Jaime Jul 15 '14 at 20:47
I just tried this in my code too, and realized what I did wrong last time (confused vertex values and vertex indices). My ConvexHull does not have a `vertices` attribute but I just used the first point in my point list instead. It probably also helps with the numerical stability if it turns out to be an interior point. I think this is as fast as it gets. Thank you! – Diana Jul 15 '14 at 21:13
Though very good and instructive, this is not the latest answer to the question – gota Jul 22 '19 at 10:49

score 13 · Answer 2 · answered Aug 24 '16 at 08:16

13

Although this question has celebrated its second birthday, I would like to point out that now, the scipy wrapper automatically reports the volume (and area) computed by Qhull.

answered Aug 24 '16 at 08:16

Urukann

475
4
10

1

Note, that you need at least the scipy version 0.17.0 for this to work. – awakenting Oct 26 '17 at 17:05

Volume of convex hull with QHull from SciPy

2 Answers2

Linked