如何绘制数据点高分辨率最外层轮廓？凸包法精度不足求方案

我有一组以黑色散点呈现的数据点，想要绘制它们的最外层轮廓。尝试用凸包（ConvexHull）计算点集的轮廓，但丢失了过多分辨率，形状细节严重损失。

使用的代码如下：

# load usual stuff
from __future__ import print_function
import sys, os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import colors
from scipy.spatial import ConvexHull

# read input file
cbm_contour = sys.argv[1]


def parse_pdb_coords(file):
    f = open(file, "r")

    coords = X = np.empty(shape=[0, 3])

    while True:
        line = f.readline()
        if not line: break

        if line.split()[0] == "ATOM" or line.split()[0] == "HETATM":
            Xcoord = float(line[30:38])
            Ycoord = float(line[38:46])
            Zcoord = float(line[46:54])

            coords = np.append(coords, [[Xcoord, Ycoord, Zcoord]], axis=0)

    return coords

##########################################################################

plt.figure(figsize=(11, 10))

# parse input file
cbm_coords = parse_pdb_coords(cbm_contour)

# consider only x- and y-axis coordinates
flattened_points = cbm_coords[:, :2]
x = cbm_coords[:,0]
y = cbm_coords[:,1]

# Find the convex hull of the flattened points
hull = ConvexHull(flattened_points)

for simplex in hull.simplices:
    plt.plot(flattened_points[simplex, 0], flattened_points[simplex, 1], color='red', lw=2)

plt.scatter(cbm_coords[:,0], cbm_coords[:,1], s=1, c='black')

plt.xlabel('X-axis coordinate ($\mathrm{\AA} $)', size=16)
plt.ylabel('Y-axis distance ($\mathrm{\AA} $)', size=16)

plt.yticks(np.arange(-20, 24, 4),size=16)
plt.xticks(np.arange(-20, 24, 4),size=16)

plt.savefig("example.png", dpi=300, transparent=False)
plt.show()

注：因数据点复杂度无法制作最小可运行示例，需适用于其他数据集的通用解决方案。

凸包仅能捕捉点集的凸边界，对含凹形或细节的点集会丢失大量信息，以下是几种通用替代方案：

1. Alpha形状（Alpha Shape） #

Alpha形状是凸包的扩展，通过调整alpha参数控制轮廓紧致度，可保留点集的凹部细节，适合需要贴合点集边缘的场景。

实现步骤： #

• 先安装依赖库：pip install alphashape
• 代码示例：

from alphashape import alphashape
import matplotlib.pyplot as plt
import numpy as np

# 沿用原数据读取逻辑
cbm_coords = parse_pdb_coords(cbm_contour)
flattened_points = cbm_coords[:, :2]

# 计算Alpha形状，alpha值越小轮廓越贴合细节
alpha = 0.5
hull = alphashape(flattened_points, alpha)

# 可视化
plt.figure(figsize=(11, 10))
plt.scatter(flattened_points[:,0], flattened_points[:,1], s=1, c='black')
plt.plot(*hull.exterior.xy, color='red', lw=2)

plt.xlabel('X-axis coordinate ($\mathrm{\AA} $)', size=16)
plt.ylabel('Y-axis distance ($\mathrm{\AA} $)', size=16)
plt.yticks(np.arange(-20, 24, 4), size=16)
plt.xticks(np.arange(-20, 24, 4), size=16)
plt.savefig("alpha_shape_example.png", dpi=300)
plt.show()

说明： #

• alpha需根据数据集调整：值越大轮廓越接近凸包；值过小会生成细碎轮廓，建议遍历0.3-1.0区间找到最优值。

2. 基于密度的轮廓（Density Contour） #

通过核密度估计（KDE）计算点集的密度分布，绘制最高密度区域的边界，适合分布密集的点集，能自然捕捉外层轮廓。

代码示例： #

from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt
import numpy as np

cbm_coords = parse_pdb_coords(cbm_contour)
flattened_points = cbm_coords[:, :2]
x, y = flattened_points[:,0], flattened_points[:,1]

# 计算核密度估计
kde = gaussian_kde([x, y])
# 创建网格用于密度计算
xgrid = np.linspace(x.min()-1, x.max()+1, 1000)
ygrid = np.linspace(y.min()-1, y.max()+1, 1000)
Xgrid, Ygrid = np.meshgrid(xgrid, ygrid)
Z = kde.evaluate([Xgrid.ravel(), Ygrid.ravel()]).reshape(Xgrid.shape)

# 可视化
plt.figure(figsize=(11, 10))
plt.scatter(x, y, s=1, c='black')
# 取密度最低的1%区域的边界作为外层轮廓（可调整百分位数）
contour_level = np.percentile(Z, 1)
plt.contour(Xgrid, Ygrid, Z, levels=[contour_level], colors='red', linewidths=2)

plt.xlabel('X-axis coordinate ($\mathrm{\AA} $)', size=16)
plt.ylabel('Y-axis distance ($\mathrm{\AA} $)', size=16)
plt.yticks(np.arange(-20, 24, 4), size=16)
plt.xticks(np.arange(-20, 24, 4), size=16)
plt.savefig("density_contour_example.png", dpi=300)
plt.show()

说明： #

• 调整contour_level的百分位数：数值越小，轮廓越紧凑；数值越大，轮廓越宽松。
• 可通过gaussian_kde(bw_method)参数调整密度平滑度，比如bw_method=0.5让密度更精细。

3. DBSCAN聚类+边界提取 #

先用DBSCAN过滤噪声点，提取核心聚类点后再计算Alpha形状或凸包，适合含离散噪声的数据集，结果更干净。

代码示例： #

from sklearn.cluster import DBSCAN
from alphashape import alphashape
import matplotlib.pyplot as plt
import numpy as np

cbm_coords = parse_pdb_coords(cbm_contour)
flattened_points = cbm_coords[:, :2]

# DBSCAN聚类，调整eps和min_samples适配点集密度
dbscan = DBSCAN(eps=1.0, min_samples=5)
labels = dbscan.fit_predict(flattened_points)

# 提取核心聚类点（排除噪声点，label=-1）
core_points = flattened_points[labels != -1]

# 对核心点计算Alpha形状
alpha = 0.6
hull = alphashape(core_points, alpha)

# 可视化
plt.figure(figsize=(11, 10))
plt.scatter(flattened_points[:,0], flattened_points[:,1], s=1, c='black')
plt.plot(*hull.exterior.xy, color='red', lw=2)

plt.xlabel('X-axis coordinate ($\mathrm{\AA} $)', size=16)
plt.ylabel('Y-axis distance ($\mathrm{\AA} $)', size=16)
plt.yticks(np.arange(-20, 24, 4), size=16)
plt.xticks(np.arange(-20, 24, 4), size=16)
plt.savefig("dbscan_alpha_example.png", dpi=300)
plt.show()

说明： #

• eps是邻域半径，min_samples是邻域内最少点数，需根据点集密度调整，确保能正确聚类目标点集。

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