Matlab delaunayn 和 Scipy Delaunay 之间的区别
Difference between Matlab delaunayn and Scipy Delaunay
我正在尝试使用 scipy.spatial.Delaunay 函数复制 Python 中由 Matlab delaunayn 函数执行的 N 维 Delaunay 三角剖分。然而,虽然 Matlab 函数给了我想要和期望的结果,但 scipy 给了我一些不同的东西。考虑到两者都是 QHull 库的包装器,我觉得这很奇怪。我假设 Matlab 在其调用中隐式设置了不同的参数。在 Matlab's documentation.
中可以找到我试图在他们两者之间复制的情况
设置是在中心有一个点的立方体,如下所示。我提供的蓝线是为了帮助可视化形状,但它们对这个问题没有任何作用或意义。
我期望的三角剖分结果为 12 个单纯形(在 Matlab 示例中列出),如下所示。
但是这个 python 等价物产生 "extra" 单纯形。
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 returns 单纯形为零体积。 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]
我想应该解决其他条件,但我保证我的点已经不包含重复项,并且方向条件似乎对我可以辨别的输出没有影响。
比较 MATLAB 和 SciPy 函数的一些笔记:
根据 MATLAB 文档,默认情况下它 uses Qt Qbb Qc
Qhull options for 3-dimensional input, while SciPy uses Qt Qbb Qc Qz
.
不确定这是否重要,但您的 NumPy 数组与在 MATLAB 中使用 ndgrid
创建的点的顺序不同。
事实上,如果您查看 edit delaunayn.m
中的 MATLAB 代码,您可以看到执行了三个额外的步骤:
- 首先它合并重复的点
mergeDuplicatePoints
(这不是你的问题)
- 然后它对点强制执行方向约定(参见代码)
最后从Qhull中得到结果(作为MEX函数实现qhullmx
),在几行代码上面有如下注释:
Strip the zero volume simplices that may have been created by the presence of degeneracy.
由于文件受版权保护,我不会post这里的代码,但你可以在你这边查看。
我正在尝试使用 scipy.spatial.Delaunay 函数复制 Python 中由 Matlab delaunayn 函数执行的 N 维 Delaunay 三角剖分。然而,虽然 Matlab 函数给了我想要和期望的结果,但 scipy 给了我一些不同的东西。考虑到两者都是 QHull 库的包装器,我觉得这很奇怪。我假设 Matlab 在其调用中隐式设置了不同的参数。在 Matlab's documentation.
中可以找到我试图在他们两者之间复制的情况设置是在中心有一个点的立方体,如下所示。我提供的蓝线是为了帮助可视化形状,但它们对这个问题没有任何作用或意义。
我期望的三角剖分结果为 12 个单纯形(在 Matlab 示例中列出),如下所示。
但是这个 python 等价物产生 "extra" 单纯形。
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 returns 单纯形为零体积。 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]
我想应该解决其他条件,但我保证我的点已经不包含重复项,并且方向条件似乎对我可以辨别的输出没有影响。
比较 MATLAB 和 SciPy 函数的一些笔记:
根据 MATLAB 文档,默认情况下它 uses
Qt Qbb Qc
Qhull options for 3-dimensional input, while SciPy usesQt Qbb Qc Qz
.不确定这是否重要,但您的 NumPy 数组与在 MATLAB 中使用
ndgrid
创建的点的顺序不同。
事实上,如果您查看 edit delaunayn.m
中的 MATLAB 代码,您可以看到执行了三个额外的步骤:
- 首先它合并重复的点
mergeDuplicatePoints
(这不是你的问题) - 然后它对点强制执行方向约定(参见代码)
最后从Qhull中得到结果(作为MEX函数实现
qhullmx
),在几行代码上面有如下注释:Strip the zero volume simplices that may have been created by the presence of degeneracy.
由于文件受版权保护,我不会post这里的代码,但你可以在你这边查看。