如何使用Python计算带插值密度的Delaunay凸包质心

核心思路 #

质心计算的本质是三个凸包域上的二重积分：

• 总质量 $M = \iint_D \rho(x,y) d\sigma$
• 质心x坐标 $C_x = \frac{1}{M}\iint_D x\cdot\rho(x,y) d\sigma$
• 质心y坐标 $C_y = \frac{1}{M}\iint_D y\cdot\rho(x,y) d\sigma$

你提到的griddata插值底层本身就依赖Delaunay三角剖分，完全不需要额外把插值结果映射到规则网格再做积分——直接在Delaunay生成的三角形单元上逐块计算积分再累加，是精度最高、速度最快、完全贴合凸包边界的方案，不会出现网格积分常见的边缘截断误差。

线性插值是这类散点密度插值最常用的选择，此时三角形单元上的积分用3点高斯求积就可以得到精确结果，不需要复杂的数值迭代。

实现步骤 #

• 导入依赖库，用散点生成Delaunay三角剖分，剖分结果天然覆盖所有散点的凸包范围，不需要额外裁剪
• 基于三角剖分构建线性插值器，效果和griddata(method='linear')完全一致，计算效率更高
• 遍历所有三角形单元，跳过面积接近0的退化三角形，对每个单元用高斯求积计算质量、x项积分、y项积分的贡献
• 所有单元贡献累加后除以总质量，得到最终质心坐标

可直接运行的代码示例 #

import numpy as np
from scipy.spatial import Delaunay
from scipy.interpolate import LinearNDInterpolator

# --------------------------
# 替换成你自己的数据即可
# points格式：N行2列的数组，每一行是一个散点的[x,y]坐标
# weights格式：长度为N的数组，对应每个散点的权重w=f(x,y)
np.random.seed(42)
points = np.random.rand(200, 2)  # 示例：生成200个单位正方形内的随机点
weights = points[:,0] + 2*points[:,1] + 1  # 示例：自定义权重，替换为你的实际w值
# --------------------------

# 生成Delaunay三角剖分
tri = Delaunay(points)
# 构建密度插值器，和griddata线性插值结果完全一致
rho_interp = LinearNDInterpolator(tri, weights)

total_mass = 0.0
cx_sum = 0.0
cy_sum = 0.0

for simplex in tri.simplices:
    # 取出当前三角形的三个顶点
    x1, y1 = points[simplex[0]]
    x2, y2 = points[simplex[1]]
    x3, y3 = points[simplex[2]]
    # 计算三角形面积
    area = 0.5 * abs((x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1))
    if area < 1e-12:  # 跳过退化的无效三角形
        continue
    # 三角形3点高斯积分点（三边中点，对二次及以下被积函数精确成立）
    int_pts = np.array([
        [(x1+x2)/2, (y1+y2)/2],
        [(x2+x3)/2, (y2+y3)/2],
        [(x3+x1)/2, (y3+y1)/2]
    ])
    # 计算积分点上的密度值
    rho = rho_interp(int_pts)
    # 累加单元贡献
    total_mass += area * rho.mean()
    cx_sum += area * (int_pts[:,0] * rho).mean()
    cy_sum += area * (int_pts[:,1] * rho).mean()

# 计算最终质心
centroid_x = cx_sum / total_mass
centroid_y = cy_sum / total_mass

print(f"凸包质心坐标：({centroid_x:.4f}, {centroid_y:.4f})")

注意事项 #

• 如果需要更光滑的密度插值，可以把LinearNDInterpolator替换为CloughTocher2DInterpolator（对应griddata的cubic插值模式），此时需要把每个三角形细分或者用更高阶的高斯积分点保证积分精度，一般工程场景下线性插值足够使用。
• 不推荐规则网格积分方案：该方案需要手动mask掉凸包外的网格点，网格分辨率低时积分误差大，分辨率高时计算量陡增，边缘区域的插值误差很难消除。
• 后续如果需要扩展到三维实体场景，逻辑完全一致：用三维Delaunay生成四面体剖分，替换为四面体单元的高斯求积规则，逐单元累加贡献即可。

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