如何从scipy.spatial.Delaunay实例计算欧拉示性数？

前言 #

scipy.spatial.Delaunay可用于计算n维空间中的Delaunay三角剖分。要计算所得结构的欧拉示性数，需统计不同维度下的单纯形数量（顶点、边、面等），再按欧拉公式计算。

欧拉示性数公式 #

对于d维单纯复形，欧拉示性数的计算公式为：

$\chi = \sum_{k=0}^d (-1)^k N_k$
其中$N_k$是去重后的k维单纯形数量（$N_0$为顶点数，$N_1$为边数，$N_2$为面数，以此类推）。

各维度单纯形统计方法 #

顶点数（$N_0$） #

直接通过delaunay.points.shape[0]获取总顶点数。注意：原问题中提到的scipy.spatial.Delaunay().vertices.shape[0]是单纯形（如2维中的三角形）的数量，并非顶点数，需区分开。

k维单纯形数（$k \geq 1$） #

由于Delaunay剖分的单纯形会共享低维面，必须对所有k维面去重后计数：

1. 生成每个d维单纯形的所有k维子集：用itertools.combinations对每个单纯形的顶点索引取$k+1$个元素（k维单纯形由$k+1$个顶点构成）。
2. 去重：将每个顶点组合排序后转为元组，存入集合中，集合的大小即为$N_k$。

代码实现 #

import numpy as np
from scipy.spatial import Delaunay
import itertools

def euler_characteristic(delaunay):
    dim = delaunay.ndim
    chi = 0
    # 统计0维顶点
    n_vertices = delaunay.points.shape[0]
    chi += (-1)**0 * n_vertices
    # 统计1到dim维的单纯形
    for k in range(1, dim + 1):
        faces = set()
        for simplex in delaunay.simplices:
            # 生成当前单纯形的所有k维面
            for face in itertools.combinations(simplex, k + 1):
                # 排序后转元组，确保同一面的不同顺序被视为同一元素
                sorted_face = tuple(sorted(face))
                faces.add(sorted_face)
        n_k = len(faces)
        chi += (-1)**k * n_k
    return chi

# 2维示例
points = np.random.rand(20, 2)
tri = Delaunay(points)
print("欧拉示性数:", euler_characteristic(tri))

# 3维示例
points_3d = np.random.rand(30, 3)
tet = Delaunay(points_3d)
print("3维欧拉示性数:", euler_characteristic(tet))

注意事项 #

• 高维空间下该方法同样适用，只需调整k的遍历范围至空间维度即可。
• 排序步骤不可省略，否则同一面的不同顶点顺序会被误判为不同面。
• 对于大规模点集，该方法可能存在性能瓶颈，可考虑使用更高效的去重方式（如利用数组排序后去重）。

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