5

我正在尝试使用 scipy.spatial.Delaunay 函数复制由 Python 中的 Matlab delaunayn 函数执行的 N 维 Delaunay 三角剖分。然而，虽然 Matlab 函数给了我想要和期望的结果，但 scipy 给了我一些不同的东西。考虑到两者都是 QHull 库的包装器，我觉得这很奇怪。我假设 Matlab 在其调用中隐式设置了不同的参数。我试图在他们两者之间复制的情况可以在Matlab 的文档中找到。

设置是有一个立方体，中心有一个点，如下所示。我提供的蓝线有助于可视化形状，但它们对这个问题没有任何目的或意义。

一个中心有一个点的立方体

我期望的三角剖分产生 12 个单纯形（在 Matlab 示例中列出），如下所示。

Matlab的三角剖分

然而，这个 python 等价物会产生“额外的”单纯形。

x = np.array([[-1,-1,-1],[-1,-1,1],[-1,1,-1],[1,-1,-1],[1,1,1],[1,1,-1],[1,-1,1],[-1,1,1],[0,0,0]])
simp = scipy.spatial.Delaunay(x).simplices

返回的变量simp应该是一个 M x N 数组，其中 M 是找到的单纯形数（在我的情况下应该是 12），N 是单纯形中的点数。在这种情况下，每个单纯形应该是一个四面体，即 N 为 4。

但我发现 M 实际上是 18 并且额外的 6 个单纯形不是四面体，而是立方体的 6 个面。

这里发生了什么？如何将返回的单纯形限制为仅是四面体？我用这个简单的案例来演示这个问题，所以我想要一个不是针对这个问题量身定制的解决方案。

编辑

感谢 Amro 的回答，我能够弄清楚这一点，并且我可以在 Matlab 和 Scipy 之间得到一个简单的匹配。有两个因素在起作用。首先，正如所指出的，Matlab 和 Scipy 使用不同的 QHull 选项。其次，QHull 返回体积为零的单纯形。Matlab 删除了这些，Scipy 没有。这在上面的示例中很明显，因为所有 6 个额外的单纯形都是立方体的零体积共面面。可以使用以下代码在 N 维中删除这些。

N = 3 # The dimensions of our points
options = 'Qt Qbb Qc' if N <= 3 else 'Qt Qbb Qc Qx' # Set the QHull options
tri = scipy.spatial.Delaunay(points, qhull_options = options).simplices
keep = np.ones(len(tri), dtype = bool)
for i, t in enumerate(tri):
    if abs(np.linalg.det(np.hstack((points[t], np.ones([1,N+1]).T)))) < 1E-15:
        keep[i] = False # Point is coplanar, we don't want to keep it
tri = tri[keep]

我想应该解决其他条件，但我保证我的点已经不包含重复项，并且方向条件似乎对我可以辨别的输出没有影响。

4

1 回答 1

3

比较 MATLAB 和 SciPy 函数的一些注释：

  • 根据 MA​​TLAB 文档，默认情况下它使用 Qt Qbb QcQhull 选项进行 3 维输入，而 SciPy使用 Qt Qbb Qc Qz.

  • 不确定是否重要，但您的 NumPy 数组与ndgrid在 MATLAB 中创建的点的顺序不同。

事实上，如果您查看 MATLAB 代码edit delaunayn.m，您可以看到执行了三个额外的步骤：

  • 首先它合并重复点mergeDuplicatePoints（在你的情况下这不是问题）
  • 然后它强制执行点的方向约定（参见代码）

  • 最后从 Qhull 得到结果（实现为 MEX-function qhullmx）后，在几行代码上方有以下注释：

    去除可能由简并性的存在产生的零体积单纯形。

由于该文件受版权保护，因此我不会在此处发布代码，但您可以自行检查。

于 2016-04-13T17:39:36.373 回答