绘制数据点的外部轮廓而不失去分辨率

Question

我有一组数据点（黑色的散点），我想画出它们的外部轮廓。我尝试过计算这些点的凸包（下面的代码），但是这样会失去太多细节，形状也变得不那么准确了。

# 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()

需要注意的是，由于数据点的复杂性，这个问题不能简化成一个“最小可工作示例”，但我的数据集可以在这里下载。我的想法是希望能找到一个适用于其他数据集的通用解决方案。

有没有人有什么建议呢？

Answer 1

你可以用一个叫做 Alpha Shape 的东西来代替你的凸包，举个例子，你可以使用 alphashape 这个库。

你需要尝试不同的 alpha 值，找到一个适合你需求的值。在这里，alpha=1 的效果不错。

# pip install alphashape
# pip install descartes

import matplotlib.pyplot as plt
from alphashape import alphashape
from descartes import PolygonPatch

ax = plt.subplot()
ax.scatter(arr[:, 0], arr[:, 1], s=1, color='k')
ax.set_aspect(1)

ashape = alphashape(arr[:, :2].tolist(), alpha=1)

ax.add_patch(PolygonPatch(ashape, fc='none', ec='red'))

输出结果：

假设输入是 arr (X, Y, Z)：

array([[-18.1  ,   1.418,   4.103],
       [-17.712,   2.202,   5.819],
       [-17.58 ,   1.998,   5.527],
       ...,
       [  9.472,  -1.901,  11.157],
       [  9.472,  -1.295,  11.615],
       [  9.472,  -0.895,  12.261]])

---

注意

如果你遇到以下错误：

IndexError: too many indices for array: array is 0-dimensional, but 2 were indexed

那么你的 shapely 版本太新了，不兼容 descartes。你可以按照 这里 的说明来解决这个问题。