使用SciPy中的QHull计算凸包体积

Question

我正在尝试使用SciPy对QHull的封装来计算一组点的凸包体积。

根据QHull的文档，我应该传递"FA"这个选项来获取总的表面积和体积。

这是我得到的结果……我哪里做错了？

> 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__']

Answer 1

虽然这个问题已经存在两年了，但我想提一下，现在scipy这个工具会自动报告Qhull计算出来的体积和面积。

Answer 2

看起来没有直接的方法可以得到你想要的结果，无论你传入什么参数。如果你用Delaunay代替ConvexHull，自己计算应该也不会太难，因为Delaunay也能提供大部分和凸包相关的信息。

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]))

---

编辑 根据评论，下面是获取凸包体积的更快的方法：

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]))

用一些虚构的数据，第二种方法似乎快了大约2倍，而且数值精度也很好（有15位小数！！！），不过也可能会有一些更复杂的情况：

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