如何用scipy.spatial.Delaunay生成五边形及以上多边形的完整三角网格？

我已生成正五边形的顶点坐标，代码如下：

import math
import numpy as np

n_angles = 5
r = 1

angles = np.linspace(0, 2 * math.pi, n_angles, endpoint=False)
x = (r * np.cos(angles))
y = (r * np.sin(angles))
corners = np.array(list(zip(x, y)))

尝试使用scipy.spatial.Delaunay生成三角网格，绘图代码如下：

import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

tri = Delaunay(points=corners)

plt.triplot(corners[:,0], corners[:,1], tri.simplices)
plt.plot(corners[:,0], corners[:,1], 'o')
for corner in range(len(corners)):
    plt.annotate(text=f'{corner + 1}', xy=(corners[:,0][corner] + 0.05, corners[:,1][corner]))

plt.axis('equal')    
plt.show()

但生成的网格中缺少(1, 3)、(1, 4)、(2, 4)这些simplices。尝试在Delaunay()中指定incremental=True并使用add_points方法，却出现如下错误：

QhullError: QH6239 Qhull precision error: initial Delaunay input sites are cocircular or cospherical.  Use option 'Qz' for the Delaunay triangulation or Voronoi diagram of cocircular/cospherical points; it adds a point "at infinity".  Alternatively use option 'QJ' to joggle the input

请问：

1. 缺失这些simplices的原因是什么，该如何添加？
2. 如何解决这个QhullError错误？

一、缺失simplices的原因及解决办法 #

原因 #

Delaunay三角剖分有个核心规矩：所有三角形的外接圆内部不能包含其他输入点。正五边形的5个顶点本来就在同一个圆上，那些跨越多条边的连线（比如1-3、1-4）对应的三角形，外接圆肯定会包住其他顶点，不符合Delaunay的规则，所以默认不会生成这些边。

解决办法 #

如果一定要得到这些对角线构成的三角网格，不能直接用Delaunay的默认输出，有两种实用方案：

1. 手动定义三角形
   直接把你需要的三角形列出来，替换掉tri.simplices来绘图。比如正五边形可以拆成3个三角形，示例代码：

   # 注意顶点索引是从0开始的，对应原问题里的1-5
   custom_simplices = np.array([[0,1,2], [0,2,3], [0,3,4]])
   plt.triplot(corners[:,0], corners[:,1], custom_simplices)
   # 后面的标注、坐标轴设置代码不变
   

   要是需要所有对角线构成的三角化，也可以自己组合所有符合要求的三角形。

2. 给正五边形加个中心点
   在正五边形中心(0,0)加一个点，这样Delaunay剖分会自动生成从中心到各顶点的连线，同时也会包含你要的那些跨边三角形。示例代码：

   # 添加中心点到顶点集合里
   center = np.array([[0,0]])
   corners_with_center = np.vstack([corners, center])
   tri = Delaunay(points=corners_with_center)
   # 绘图代码和之前一样就行
   

二、QhullError错误的解决办法 #

错误提示里已经给了两种官方解决方式：

1. 用Qz选项
   初始化Delaunay的时候加上qhull_options="Qz"，这个选项会在无穷远处加一个虚拟点，专门处理共圆的输入点：

   tri = Delaunay(points=corners, qhull_options="Qz")
   

2. 用QJ选项
   加上qhull_options="QJ"，会给输入点加一点点微小的扰动，打破共圆的特性，避免精度报错：

   tri = Delaunay(points=corners, qhull_options="QJ")
   

   要是想用incremental=True模式，初始化的时候就得带上这些选项解决共圆问题，比如：

   tri = Delaunay(points=corners, incremental=True, qhull_options="Qz")
   # 之后就能正常用add_points加新点了
   

内容的提问来源于stack exchange，提问作者HJA24