在matplotlib中使用三角测量时如何处理在我的几何边缘之间形成的（不需要的）三角形

Question

我有一个用空间中 (x,y) 点列表定义的几何。我想用这些数据创建一个三角形网格，所以我为此尝试了Triangulation function in matplotlib。但是，由于我的几何图形有一些曲线，因此该算法会在我的零件边缘之间生成不需要的三角形：

红色曲线是我的几何图形的边缘。

有什么办法可以解决这个问题吗？可能三角函数不是我需要的，在这种情况下，你有什么建议吗？

以下代码来自this example。在示例中，他们通过显式命名三个点而不是我想通过调用函数使用的 Delaunay 三角剖分来定义三角形：triang = tri.Triangulation(x, y)，这将使我得到与原始图片相同的行为。

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

xy = np.asarray([
    [-0.101, 0.872], [-0.080, 0.883], [-0.069, 0.888], [-0.054, 0.890],
    [-0.045, 0.897], [-0.057, 0.895], [-0.073, 0.900], [-0.087, 0.898],
    [-0.090, 0.904], [-0.069, 0.907], [-0.069, 0.921], [-0.080, 0.919],
    [-0.073, 0.928], [-0.052, 0.930], [-0.048, 0.942], [-0.062, 0.949],
    [-0.054, 0.958], [-0.069, 0.954], [-0.087, 0.952], [-0.087, 0.959],
    [-0.080, 0.966], [-0.085, 0.973], [-0.087, 0.965], [-0.097, 0.965],
    [-0.097, 0.975], [-0.092, 0.984], [-0.101, 0.980], [-0.108, 0.980],
    [-0.104, 0.987], [-0.102, 0.993], [-0.115, 1.001], [-0.099, 0.996],
    [-0.101, 1.007], [-0.090, 1.010], [-0.087, 1.021], [-0.069, 1.021],
    [-0.052, 1.022], [-0.052, 1.017], [-0.069, 1.010], [-0.064, 1.005],
    [-0.048, 1.005], [-0.031, 1.005], [-0.031, 0.996], [-0.040, 0.987],
    [-0.045, 0.980], [-0.052, 0.975], [-0.040, 0.973], [-0.026, 0.968],
    [-0.020, 0.954], [-0.006, 0.947], [ 0.003, 0.935], [ 0.006, 0.926],
    [ 0.005, 0.921], [ 0.022, 0.923], [ 0.033, 0.912], [ 0.029, 0.905],
    [ 0.017, 0.900], [ 0.012, 0.895], [ 0.027, 0.893], [ 0.019, 0.886],
    [ 0.001, 0.883], [-0.012, 0.884], [-0.029, 0.883], [-0.038, 0.879],
    [-0.057, 0.881], [-0.062, 0.876], [-0.078, 0.876], [-0.087, 0.872],
    [-0.030, 0.907], [-0.007, 0.905], [-0.057, 0.916], [-0.025, 0.933],
    [-0.077, 0.990], [-0.059, 0.993]])
x = np.degrees(xy[:, 0])
y = np.degrees(xy[:, 1])

triang = tri.Triangulation(x, y)
fig1, ax1 = plt.subplots()
ax1.set_aspect('equal')
ax1.triplot(triang, 'bo-', lw=1)

Answer 1

如果几何形状定义明确，比如曲线，可以检查每个三角形是否在形状内。然后可以定义一个蒙版并屏蔽掉不需要的三角形。我找到了一个使用shapely 的解决方案，其中为原始形状（outline）定义一个多边形，并为Triangulation() 的每个生成三角形定义一个多边形，然后我检查它是否位于outline 内：

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

import shapely
from shapely.geometry import Polygon as sPolygon


fig, (ax1,ax2) = plt.subplots(ncols=2)
ax1.set_aspect('equal')
ax2.set_aspect('equal')

##setting up basic shape
phi = np.linspace(0,2*np.pi,20)
r = 1 + 2*np.sin(phi)**2
x = np.cos(phi)*r
y = np.sin(phi)*r
ax1.plot(x,y,'ro-', lw=3, ms=6, zorder= 1, label='edge')
ax2.plot(x,y,'ro-', lw=3, ms=6, zorder= 1)


##original triangulation
triang1 = tri.Triangulation(x, y)
ax1.triplot(triang1, 'ko--', lw=1, ms=4, zorder=2, label='all')

##masking
outline = sPolygon(zip(x,y))
mask = [
    not outline.contains(sPolygon(zip(x[tri], y[tri])))
    for tri in triang1.get_masked_triangles()
]
triang1.set_mask(mask)
ax1.triplot(triang1, 'b-', lw=1, zorder=3, label='inner')

##adding more points
x_extra = np.random.rand(30)*(x.max()-x.min())+x.min()
y_extra = np.random.rand(30)*(y.max()-y.min())+y.min()

x = np.concatenate([x,x_extra])
y = np.concatenate([y,y_extra])

triang2 = tri.Triangulation(x,y)
ax2.triplot(triang2, 'ko--', lw=1, ms=4,  zorder=2)

##masking
mask = [
    not outline.contains(sPolygon(zip(x[tri], y[tri])))
    for tri in triang2.get_masked_triangles()
]
triang2.set_mask(mask)
ax2.triplot(triang2, 'b-', lw=1, zorder=3)

fig.legend()
plt.show()

代码的结果如下所示：

我不太确定 OP 想要什么，所以在左侧我只使用边缘的点，而在右侧我添加了一些随机的额外点用于三角测量。图中，形状的轮廓用红色绘制，Delaunay三角剖分的原始结果用黑色虚线绘制，蒙版三角剖分用蓝色绘制。

编辑：

我刚刚注意到，过滤后的右侧图片中显然没有包含其中一个轮廓点。这一定是由于数字不准确。解决此问题的一种方法是使用buffer() 命令稍微增加轮廓的大小。像这样的东西似乎适合手头的问题：

outline = sPolygon(zip(x,y)).buffer(.01)

但可能需要调整实际的缓冲量。

Answer 2

如果您有内部形状的轮廓来绘制三角剖分，您可以应用@ThomasKühn 的答案。

否则，您可能有一个不应该考虑三角形的点之间的最大距离。在这种情况下，您可以掩盖这些三角形。

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

xy = np.asarray([
    [-0.101, 0.872], [-0.080, 0.883], [-0.069, 0.888], [-0.054, 0.890],
    [-0.045, 0.897], [-0.057, 0.895], [-0.073, 0.900], [-0.087, 0.898],
    [-0.090, 0.904], [-0.069, 0.907], [-0.069, 0.921], [-0.080, 0.919],
    [-0.073, 0.928], [-0.052, 0.930], [-0.048, 0.942], [-0.062, 0.949],
    [-0.054, 0.958], [-0.069, 0.954], [-0.087, 0.952], [-0.087, 0.959],
    [-0.080, 0.966], [-0.085, 0.973], [-0.087, 0.965], [-0.097, 0.965],
    [-0.097, 0.975], [-0.092, 0.984], [-0.101, 0.980], [-0.108, 0.980],
    [-0.104, 0.987], [-0.102, 0.993], [-0.115, 1.001], [-0.099, 0.996],
    [-0.101, 1.007], [-0.090, 1.010], [-0.087, 1.021], [-0.069, 1.021],
    [-0.052, 1.022], [-0.052, 1.017], [-0.069, 1.010], [-0.064, 1.005],
    [-0.048, 1.005], [-0.031, 1.005], [-0.031, 0.996], [-0.040, 0.987],
    [-0.045, 0.980], [-0.052, 0.975], [-0.040, 0.973], [-0.026, 0.968],
    [-0.020, 0.954], [-0.006, 0.947], [ 0.003, 0.935], [ 0.006, 0.926],
    [ 0.005, 0.921], [ 0.022, 0.923], [ 0.033, 0.912], [ 0.029, 0.905],
    [ 0.017, 0.900], [ 0.012, 0.895], [ 0.027, 0.893], [ 0.019, 0.886],
    [ 0.001, 0.883], [-0.012, 0.884], [-0.029, 0.883], [-0.038, 0.879],
    [-0.057, 0.881], [-0.062, 0.876], [-0.078, 0.876], [-0.087, 0.872],
    [-0.030, 0.907], [-0.007, 0.905], [-0.057, 0.916], [-0.025, 0.933],
    [-0.077, 0.990], [-0.059, 0.993]])
x = np.degrees(xy[:, 0])
y = np.degrees(xy[:, 1])

triang = tri.Triangulation(x, y)

fig1, ax1 = plt.subplots()
ax1.set_aspect('equal')

# plot all triangles
ax1.triplot(triang, 'bo-', lw=0.2)

# plot only triangles with sidelength smaller some max_radius
max_radius = 2
triangles = triang.triangles

# Mask off unwanted triangles.
xtri = x[triangles] - np.roll(x[triangles], 1, axis=1)
ytri = y[triangles] - np.roll(y[triangles], 1, axis=1)
maxi = np.max(np.sqrt(xtri**2 + ytri**2), axis=1)
triang.set_mask(maxi > max_radius)

ax1.triplot(triang, color="indigo", lw=2.6)


plt.show()

细线显示所有三角形（点的凸包），粗线仅显示边长不大于某些最大值的三角形（在本例中选择为2）。

此线程可能同样相关：matplotlib contour/contourf of **concave** non-gridded data