在 matplotlib 中使用三角剖分时如何处理在我的几何图形边缘之间形成的（不需要的）三角形
How to deal with the (undesired) triangles that form between the edges of my geometry when using Triangulation in matplotlib
我有一个用 space 中的 (x,y) 点列表定义的几何体。我想用这些数据创建一个三角形网格，所以我为此尝试了 Triangulation function in matplotlib。但是，由于我的几何图形有一些曲线，算法在我的零件的边缘之间生成了不需要的三角形：
红色曲线是我的几何图形的边缘。
有什么办法可以解决这个问题吗？可能我需要的不是Triangulation功能，那你有什么推荐的吗？
以下代码来自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)
如果几何形状定义明确，例如曲线，则可以检查每个三角形是否位于该形状内。然后可以定义一个掩码并屏蔽不需要的三角形。我找到了一个使用 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三角剖分原始结果，蓝色为masked三角剖分。
编辑:
我刚刚注意到，过滤后右侧图片中显然没有包含一个轮廓点。这一定是由于数值不准确造成的。解决这个问题的一种方法是使用 buffer() 命令稍微增加轮廓的大小。像这样的事情似乎适合手头的问题：
outline = sPolygon(zip(x,y)).buffer(.01)
但实际缓冲量可能需要调整。
如果您有内部形状的轮廓来绘制三角剖分，您可以应用@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
我有一个用 space 中的 (x,y) 点列表定义的几何体。我想用这些数据创建一个三角形网格，所以我为此尝试了 Triangulation function in matplotlib。但是，由于我的几何图形有一些曲线，算法在我的零件的边缘之间生成了不需要的三角形：
红色曲线是我的几何图形的边缘。
有什么办法可以解决这个问题吗？可能我需要的不是Triangulation功能，那你有什么推荐的吗？
以下代码来自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)
如果几何形状定义明确，例如曲线，则可以检查每个三角形是否位于该形状内。然后可以定义一个掩码并屏蔽不需要的三角形。我找到了一个使用 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三角剖分原始结果，蓝色为masked三角剖分。
编辑:
我刚刚注意到，过滤后右侧图片中显然没有包含一个轮廓点。这一定是由于数值不准确造成的。解决这个问题的一种方法是使用 buffer() 命令稍微增加轮廓的大小。像这样的事情似乎适合手头的问题：
outline = sPolygon(zip(x,y)).buffer(.01)
但实际缓冲量可能需要调整。
如果您有内部形状的轮廓来绘制三角剖分，您可以应用@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