基于QHull的SciPy凸包体积

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

不管您传入哪些参数，似乎没有任何直接获得所需结果的明显方法。如果您使用的是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

Answer 2

虽然这个问题已经庆祝了它的第二个生日，我想指出，现在，枕木包装自动报告的体积(和面积)计算的Qhull。