使用scipy.spatial.Delaunay替代matplotlib.tri.Triangulation内置版本

8
似乎matplotlib.tri.Triangulation使用了有缺陷并且可能是不正确的Delaunay三角剖分实现,这个实现应该被替换成qHull
我试图使用mpl_toolkits.mplot3d.plot_trisurf()绘制一个trisurf图形,但是遇到了很多无用的异常(主要是IndexErrorKeyError),没有任何关于出错原因的指示。
由于scipy.spatial.Delaunay已经使用了qHull,因此我想知道是否有一种方法可以利用scipy的Delaunay三角剖分实现构建一个matplotlib.tri.Triangulation对象,以便在mpl_toolkits.mplot3d.plot_trisurf()中使用。
我已经尝试将delaunay.points直接传递给matplotlib.tri.Triangulate,通过triangles参数来实现,但结果会导致ValueError: triangles min element is out of bounds

顺便提一下,如果您可以尝试使用当前的主分支,这个问题将在1.4中得到解决。 - tacaswell
2个回答

6
您需要将qhull生成的点和三角形传递给Triangulation构造函数:

http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.spatial.Delaunay.html http://matplotlib.org/dev/api/tri_api.html

import numpy as np
import scipy.spatial
import matplotlib
import math

import matplotlib.tri as mtri
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# First create the x and y coordinates of the points.
n_angles = 20
n_radii = 10
min_radius = 0.15
radii = np.linspace(min_radius, 0.95, n_radii)
angles = np.linspace(0, 2*math.pi, n_angles, endpoint=False)
angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1)
angles[:, 1::2] += math.pi/n_angles
x = (radii*np.cos(angles)).flatten()
y = (radii*np.sin(angles)).flatten()

# Create the Delaunay tessalation using scipy.spatial
pts = np.vstack([x, y]).T
tess = scipy.spatial.Delaunay(pts)

# Create the matplotlib Triangulation object
x = tess.points[:, 0]
y = tess.points[:, 1]
tri = tess.vertices # or tess.simplices depending on scipy version
triang = mtri.Triangulation(x=pts[:, 0], y=pts[:, 1], triangles=tri)

# Plotting
z = x*x + y*y
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(triang, z)
plt.show()

输出结果(使用当前主分支的matplotlib): enter image description here


如果我的输入只有一个二维数组呢? - marco

-1

@Marco很好奇如何在2D数组上运行这个程序。希望这对你有用。根据坐标的顶点列表应该被转换为一个数组,并且可以使用mtri.Triangulation进行三角剖分。 以下是示例代码:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri


verts = np.array([[0.6,0.8],[0.2,0.9],[0.1,-0.5],[0.2,-2]])

triang = mtri.Triangulation(verts[:,0], verts[:,1])

plt.triplot(triang, marker="o")

plt.show()`enter code here`

1
欢迎来到 Stack Overflow。代码配合解释会更有帮助。Stack Overflow 的宗旨是学习,而不是提供盲目复制粘贴的片段。当回答已有答案的老问题时(本问题已经存在八年半了),这一点尤为重要。请编辑您的答案,解释它如何回答具体问题,并改进现有答案。请参见[答案]。 - Chris

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接