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

Question

看起来 matplotlib.tri.Triangulation 这个工具用的 Delaunay 三角剖分方法有点问题，可能是个bug，而且它的实现方式也不太对，未来可能会换成 qHull。

我想用 mpl_toolkits.mplot3d.plot_trisurf() 来画一个三维表面图，但遇到了一堆错误，主要是 IndexError 和 KeyError，这些错误信息都没告诉我到底哪里出错了。

因为 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 的错误。

Answer 1

@Marco 想知道如何在二维数组上运行这个。我希望这对你有帮助。根据坐标的顶点列表应该变成一个数组，并可以使用 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`

Answer 2

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

所以你需要把点和三角形都从qhull传递给Triangulation构造函数：

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当前的主版本）：