我正在尝试使用QHull的SciPy包装器来获取一组点的凸包体积.
根据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__']
Jaime.. 15
似乎没有任何明显的方法可以直接获得你所追求的结果,无论你输入什么参数.如果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
小智.. 11
虽然这个问题庆祝了它的第二个生日,但我想指出现在,scipy包装器会自动报告Qhull计算的音量(和面积).