如何在Python散点图中绘制包裹所有点的紧密填充区域？

问题分析 #

scipy.spatial.ConvexHull只能生成包裹点集的最小凸多边形边界，无法捕捉你数据中存在的非凸曲线轮廓。要生成贴合所有点的紧密非凸边界，推荐使用Alpha Shape算法，或分上下边界拟合的方法。

方法一：Alpha Shape算法（推荐） #

Alpha Shape通过Delaunay三角剖分筛选符合参数的三角形边，能生成贴合非凸点集的边界。

1. 读取数据 #

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

# 读取数据（根据你的文件格式调整，示例为CSV）
df = pd.read_csv("你的数据文件.csv")
yloc = df['Points:1'].values
zloc = df['Points:2'].values
points = np.column_stack((yloc, zloc))

2. 实现Alpha Shape工具函数 #

def alpha_shape(points, alpha, only_outer=True):
    assert points.shape[0] > 3, "点集数量必须大于3"
    tri = Delaunay(points)
    edges = set()

    # 遍历所有三角形，筛选符合alpha条件的边
    for ia, ib, ic in tri.simplices:
        pa, pb, pc = points[ia], points[ib], points[ic]
        # 计算三角形外接圆半径
        a, b, c = np.linalg.norm(pa-pb), np.linalg.norm(pb-pc), np.linalg.norm(pc-pa)
        s = (a + b + c) / 2.0
        area = np.sqrt(s * (s-a) * (s-b) * (s-c)) if s*(s-a)*(s-b)*(s-c) > 0 else 0
        circum_r = a*b*c/(4.0*area) if area > 1e-10 else float('inf')
        
        if circum_r <= 1.0/alpha:
            edges.add(tuple(sorted((ia, ib))))
            edges.add(tuple(sorted((ib, ic))))
            edges.add(tuple(sorted((ic, ia))))

    # 构建连续的边界路径
    edge_list = list(edges)
    if not edge_list:
        return []
    
    current_edge = edge_list.pop(0)
    boundary = [current_edge[0], current_edge[1]]
    
    while edge_list:
        next_idx = None
        for i, e in enumerate(edge_list):
            if e[0] == boundary[-1]:
                next_idx = e[1]
                edge_list.pop(i)
                break
            elif e[1] == boundary[-1]:
                next_idx = e[0]
                edge_list.pop(i)
                break
        if next_idx is None:
            break
        boundary.append(next_idx)
    
    # 闭合边界
    if boundary[0] != boundary[-1]:
        boundary.append(boundary[0])
    
    # 过滤内部小环（可选）
    if only_outer:
        def poly_area(vertices):
            area = 0.0
            n = len(vertices)
            for i in range(n):
                x1, y1 = vertices[i]
                x2, y2 = vertices[(i+1)%n]
                area += (x1*y2) - (x2*y1)
            return abs(area)/2.0
        
        boundary_points = points[boundary]
        if poly_area(boundary_points) < 0.1:  # 阈值根据数据调整
            return []
    
    return boundary

3. 绘制填充边界 #

调整alpha参数控制边界紧密程度（值越小越贴合点集）：

plt.scatter(yloc, zloc, s=5, label='原始点集')

# 计算并绘制Alpha Shape边界
alpha = 0.05  # 可根据实际效果调整，比如0.01~0.1之间尝试
boundary_idx = alpha_shape(points, alpha)
if boundary_idx:
    boundary_points = points[boundary_idx]
    plt.fill(boundary_points[:,0], boundary_points[:,1], color='red', alpha=0.2, label='紧密边界')
    plt.plot(boundary_points[:,0], boundary_points[:,1], 'r-', linewidth=1.5)

plt.xlabel('Y坐标')
plt.ylabel('Z坐标')
plt.legend()
plt.show()

方法二：分上下边界拟合 #

如果点集可按Y坐标排序，可分别提取上下边界并平滑拟合，再填充区域：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import splrep, splev

# 读取数据
df = pd.read_csv("你的数据文件.csv")
yloc = df['Points:1'].values
zloc = df['Points:2'].values

# 按Y坐标排序
sorted_idx = np.argsort(yloc)
sorted_y = yloc[sorted_idx]
sorted_z = zloc[sorted_idx]

# 滚动窗口提取上下边界极值
window_size = 20
upper_z = []
lower_z = []
for i in range(len(sorted_y)):
    start = max(0, i - window_size//2)
    end = min(len(sorted_y), i + window_size//2)
    upper_z.append(np.max(sorted_z[start:end]))
    lower_z.append(np.min(sorted_z[start:end]))

# 平滑曲线
upper_tck = splrep(sorted_y, upper_z, s=100)  # s值控制平滑程度
lower_tck = splrep(sorted_y, lower_z, s=100)
smooth_y = np.linspace(sorted_y.min(), sorted_y.max(), 500)
smooth_upper = splev(smooth_y, upper_tck)
smooth_lower = splev(smooth_y, lower_tck)

# 绘制
plt.scatter(yloc, zloc, s=5, label='原始点集')
plt.fill_between(smooth_y, smooth_lower, smooth_upper, color='red', alpha=0.2, label='紧密边界')
plt.plot(smooth_y, smooth_upper, 'r-', linewidth=1.5)
plt.plot(smooth_y, smooth_lower, 'r-', linewidth=1.5)

plt.xlabel('Y坐标')
plt.ylabel('Z坐标')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Atharva Sunil Sathe