2D笛卡尔坐标系点云孔洞区域提取及聚类数确定技术问询

嗨，我看你现在用固定n_clusters=5的KMeans来分离孔洞区域，但没法提前确定实际的孔洞数量，这确实是KMeans的典型局限——它必须预先指定聚类数。不过有几个实用的办法能解决这个问题，还能帮你顺利计算每个孔洞的最大直径，我给你一步步说：

一、先解决聚类数自动确定的问题 #

针对你的2D孔洞点集，有两种思路：要么给KMeans自动找最优K，要么换用不需要指定聚类数的密度聚类算法（比如DBSCAN）。

方法1：用肘部法则自动找最优K #

KMeans的SSE（误差平方和）会随着K的增大而递减，但当K达到真实聚类数时，SSE的下降速度会突然变缓（形成“肘部”），我们可以利用这个特征自动选K。

代码示例：

from sklearn.cluster import KMeans
import numpy as np
import matplotlib.pyplot as plt

# 你的2D点数据
xycoors = dat[:, 0:2]

# 尝试K从1到10，计算每个K的SSE
sse = []
k_range = range(1, 11)
for k in k_range:
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(xycoors)
    sse.append(kmeans.inertia_)

# 绘制肘部图
plt.plot(k_range, sse)
plt.xlabel('聚类数K')
plt.ylabel('SSE（误差平方和）')
plt.title('肘部法则找最优K')
plt.show()

# 找到肘部对应的K（你可以手动从图里选，或者写简单逻辑自动判断）
best_k = 4  # 根据你的肘部图调整
fit = KMeans(n_clusters=best_k, random_state=42).fit(xycoors)
clus_datas = {i: xycoors[np.where(fit.labels_ == i)] for i in range(fit.n_clusters)}

方法2：用轮廓系数选最优K #

轮廓系数能衡量每个样本在聚类中的“贴合度”，取值范围[-1,1]，越接近1说明聚类效果越好。我们可以计算不同K下的平均轮廓系数，选系数最高的K。

代码示例：

from sklearn.metrics import silhouette_score

sil_scores = []
for k in range(2, 11):  # 轮廓系数至少需要2个聚类
    kmeans = KMeans(n_clusters=k, random_state=42)
    labels = kmeans.fit_predict(xycoors)
    sil_score = silhouette_score(xycoors, labels)
    sil_scores.append(sil_score)

# 找最高轮廓系数对应的K
best_k = np.argmax(sil_scores) + 2  # 因为k从2开始
fit = KMeans(n_clusters=best_k, random_state=42).fit(xycoors)
clus_datas = {i: xycoors[np.where(fit.labels_ == i)] for i in range(fit.n_clusters)}

方法3：换用DBSCAN（更适合孔洞这类密集点群） #

如果你的孔洞是独立的密集点区域，DBSCAN会更合适——它不需要预先指定聚类数，而是基于点的密度自动划分簇，还能过滤掉噪声点。

代码示例：

from sklearn.cluster import DBSCAN

# 调整eps（邻域半径）和min_samples（邻域内最少点数）参数，适配你的数据
dbscan = DBSCAN(eps=0.5, min_samples=10)  # 参数需要根据你的点密度调整
labels = dbscan.fit_predict(xycoors)

# 整理聚类结果（去掉噪声点，标签为-1的是噪声）
clus_datas = {}
unique_labels = set(labels) - {-1}
for label in unique_labels:
    clus_datas[label] = xycoors[np.where(labels == label)]

二、计算每个孔洞的最大直径 #

拿到聚类后的点集后，计算最大直径最高效的方式是先求点集的凸包，然后找凸包上所有点对的最大距离（凸包的直径就是点集的最大直径）：

from scipy.spatial import ConvexHull
import itertools

def calculate_max_diameter(points):
    if len(points) < 2:
        return 0
    # 计算凸包
    hull = ConvexHull(points)
    hull_points = points[hull.vertices]
    # 计算凸包上所有点对的距离，取最大值
    max_dist = 0
    for pair in itertools.combinations(hull_points, 2):
        dist = np.linalg.norm(pair[0] - pair[1])
        if dist > max_dist:
            max_dist = dist
    return max_dist

# 遍历每个聚类，计算最大直径
for cluster_id, points in clus_datas.items():
    diameter = calculate_max_diameter(points)
    print(f"孔洞{cluster_id+1}的最大直径：{diameter:.2f}")

最后调整你的可视化代码 #

把原来固定的聚类替换成自动得到的clus_datas，可以用循环来绘制，避免重复代码：

import ipyvolume.pylab as p

fig = p.figure(width=1000)
fig.xlabel='x'
fig.ylabel='z'
fig.zlabel='y'

# 定义颜色列表
colors = ["black", "red", "green", "blue", "orange", "purple"]
for idx, (cluster_id, points) in enumerate(clus_datas.items()):
    color = colors[idx % len(colors)]
    p.scatter(points[:,1], points[:,1]*0, points[:,0], color=color, size=.1)

p.squarelim()
p.show()

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