如何将单纯形添加到 scipy Delaunay 三角剖分对象答案

【问题标题】：How to add simplices to a scipy Delaunay triangulation object如何将单纯形添加到 scipy Delaunay 三角剖分对象
【发布时间】：2019-02-16 14:11:32
【问题描述】：

我已经有一个由scipy.spatial.Delaunay() 对象三角剖分的矩形。我设法拉伸和弯曲它，使它看起来像一个沿着一条线切割的环。下面是一些代码来制作具有相同拓扑的东西：

from scipy.spatial import Delaunay

NR = 22
NTheta = 36

Rin = 1
Rout = 3
alphaFactor = 33/64
alpha = np.pi/alphaFactor # opening angle of wedge

u=np.linspace(pi/2, pi/2 + alpha, NTheta)
v=np.linspace(Rin, Rout, NR)
u,v=np.meshgrid(u,v)
u=u.flatten()
v=v.flatten()

#evaluate the parameterization at the flattened u and v
x=v*np.cos(u)
y=v*np.sin(u)

#define 2D points, as input data for the Delaunay triangulation of U
points2D=np.vstack([u,v]).T
xy0 = np.vstack([x,y]).T
triLattice = Delaunay(points2D) #triangulate the rectangle U
triSimplices = triLattice.simplices
plt.figure()
plt.triplot(x, y, triSimplices, linewidth=0.5)

从这个拓扑开始，我现在想把两个开放的边连接起来，做一个封闭的环（改变拓扑，也就是说）。如何手动将新三角形添加到现有三角剖分？

【问题讨论】：

之后你想如何使用三角测量？它必须是有效的scipy.spatial.qhull.Delaunay 对象（例如使用find_simplex 方法），还是拥有三角形列表就足够了？
您能分享您的代码或精简版吗？这样做的一种方法是更改节点引用，但如果您的代码不方便，就很难看到存储在哪里
@xdze2 感谢您提出这个问题。它不需要是有效的Delaunay 对象。我只需要修改后的单纯形列表。

标签： python scipy computational-geometry scipy-spatial

【解决方案1】：

解决方案是合并间隙周围的点。这是一种通过跟踪相应点的索引来做到这一点的方法：

import matplotlib.pylab as plt
from scipy.spatial import Delaunay
import numpy as np

NR = 4
NTheta = 16

Rin = 1
Rout = 3
alphaFactor = 33/64  # -- set to .5 to close the gap
alpha = np.pi/alphaFactor  # opening angle of wedge

u = np.linspace(np.pi/2, np.pi/2 + alpha, NTheta)
v = np.linspace(Rin, Rout, NR)
u_grid, v_grid = np.meshgrid(u, v)
u = u_grid.flatten()
v = v_grid.flatten()

# Get the indexes of the points on the first and last columns:
idx_grid_first = (np.arange(u_grid.shape[0]),
                  np.zeros(u_grid.shape[0], dtype=int))

idx_grid_last = (np.arange(u_grid.shape[0]),
                 (u_grid.shape[1]-1)*np.ones(u_grid.shape[0], dtype=int))

# Convert these 2D indexes to 1D indexes, on the flatten array:
idx_flat_first = np.ravel_multi_index(idx_grid_first, u_grid.shape)
idx_flat_last = np.ravel_multi_index(idx_grid_last, u_grid.shape)

# Evaluate the parameterization at the flattened u and v
x = v * np.cos(u)
y = v * np.sin(u)

# Define 2D points, as input data for the Delaunay triangulation of U
points2D = np.vstack([u, v]).T
triLattice = Delaunay(points2D) # triangulate the rectangle U
triSimplices = triLattice.simplices

# Replace the 'last' index by the corresponding 'first':
triSimplices_merged = triSimplices.copy()
for i_first, i_last in zip(idx_flat_first, idx_flat_last):
    triSimplices_merged[triSimplices == i_last] = i_first

# Graph
plt.figure(figsize=(7, 7))
plt.triplot(x, y, triSimplices, linewidth=0.5)
plt.triplot(x, y, triSimplices_merged, linewidth=0.5, color='k')
plt.axis('equal');

plt.plot(x[idx_flat_first], y[idx_flat_first], 'or', label='first')
plt.plot(x[idx_flat_last], y[idx_flat_last], 'ob', label='last')
plt.legend();

给出：

也许您必须调整alphaFactor 的定义，以使间隙具有正确的大小。

【讨论】：