来自 SciPy 的带有 QHull 的凸包体积

Posted 2023-02-14

【中文标题】来自 SciPy 的带有 QHull 的凸包体积

【英文标题】：Volume of convex hull with QHull from SciPy

【发布时间】：2014-09-04 04:07:31

【问题描述】：

我正在尝试使用SciPy wrapper for QHull 获取一组点的凸包体积。

根据documentation of 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__']

【问题讨论】：

尝试用一个真实的问题来更新你的问题（不是“这就是我得到的”）。我花了一段时间才弄清楚总面积和体积在哪里都找不到，尽管事实上您提供了正确的选项。

我的猜测是 SciPy 没有包装那个特定的选项标志。

难的是实现它：wiki.scipy.org/Cookbook/Finding_Convex_Hull

完成pts 会有所帮助。这样我们就可以自己尝试了。

它没有在 Scipy Qhull 包装器中实现。如果需要，可以轻松添加。

【参考方案1】：

无论您传入什么参数，似乎都没有任何明显的方法可以直接获得您所追求的结果。如果您使用@而不是ConvexHull，那么您自己计算应该不会太难987654321@（也提供了大部分凸包相关信息）。

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

---

编辑根据 cmets，以下是获取凸包体积的更快方法：

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

【讨论】：

这很棒。我一直在寻找来自 qhull 的卷信息以节省计算时间，但我认为在 Python 中它不能比你的版本更快地完成。所以只是为了确认一下，Delaunay 做了凸包的镶嵌，对吗？另外，这个einsum 的小东西是……不知道。 :)

Delaunay 结果非常慢，所以我改用 ConvexHull，获取末端点（船体上的点），然后仅在这些点上运行 Delaunay。在我的数据上，这种方式要快 1-2 个数量级。

这确实很聪明......它可能不会进一步改善，但您可能想尝试完全跳过对Delaunay的调用，并通过在凸包上选择一个点来构建凸包的三角剖分hull，然后计算包含该点的所有四面体的体积以及每个凸包的单纯面上的点（即ConvexHull 对象的.simplices 属性）。它将产生更细长的四面体，因此它可能在数值上不太稳定。但它必须更快。

我也试过了，但是结果差别太大了。也许我做错了什么。鉴于它已经是凸包，数值不稳定性无法解释体积值的两个数量级差异，对吧？我可能会再看看这个。

看看最新的编辑，对于随机数据，我得到了 2 倍的加速和精确到小数点后 15 位...

【参考方案2】：

虽然这个问题已经庆祝了它的第二个生日，但我想指出的是，现在 scipy 包装器会自动报告 Qhull 计算的体积（和面积）。

【讨论】：

请注意，您至少需要 scipy 版本 0.17.0 才能正常工作。