基于BLS熵与t-SNE的形状聚类：从网络熵到无监督分类的实践

1. 项目概述：一种基于BLS熵的形状聚类新思路

在计算机视觉和模式识别领域，我们经常面对一个看似简单实则棘手的问题：如何让机器像人一样，理解并归类千变万化的形状？无论是从卫星图像中识别不同种类的云层，还是在医学影像中区分病变细胞，亦或是在工业质检中筛选出有缺陷的零件，其核心都离不开对形状的精准分析和分类。传统的解决方案，比如傅里叶描述符或泽尔尼克矩，为我们提供了有力的数学工具，但它们有时就像拿着固定尺寸的筛子去筛分不同大小的沙粒——在面对旋转、缩放或非刚性形变时，其效果往往会大打折扣。

我最近深入研究了一篇来自韩国国家数学科学院的论文，其中提出了一种令人耳目一新的方法。它不再仅仅依赖形状的轮廓点坐标或区域像素，而是转向了一种名为“分支长度相似性熵”的度量。简单来说，BLS熵衡量的是一个形状边界上所有点之间连接关系（想象成一个网络）的均匀程度。这个指标天生就对平移、旋转和缩放免疫，因为它关注的是形状内部点与点之间的相对关系，而非绝对位置。基于此，研究者构建了一个仅由三个核心统计量（熵值、均值、方差）组成的三维特征向量。这个向量的精妙之处在于，它用极低的维度，捕捉了形状最本质的几何“指纹”。

接下来的流程更是将现代数据科学的精髓发挥得淋漓尽致：先用t-SNE这个强大的非线性降维工具，将高维（这里是三维）特征向量投影到二维平面，其目的是在降维过程中最大限度地保持形状之间的“邻里关系”——相似形状在二维空间中依然紧挨着。最后，再对这片二维“星图”施加经典的k-means聚类算法，自然而然地划分出不同的形状族群。

这套组合拳的效果如何？论文在MPEG-7（1400个形状）、瑞典树叶（1125片叶子）和七边形（1936个不规则七边形）三个经典数据集上进行了测试。结果令人振奋，尤其是在树叶分类上，聚类得分达到了0.885，显著优于传统方法。这不仅仅是一个数字的提升，更意味着我们找到了一种更稳健、更通用的形状描述方式。无论你是正在研究目标检测的算法工程师，还是处理生物形态学的科研人员，亦或是任何需要对图像形状进行自动化分析的从业者，理解这套方法的原理与实现，都将为你打开一扇新的大门。接下来，我将带你深入这套方法的每一个技术细节，并分享从原理到代码实现的完整路径。

2. 核心原理深度剖析：BLS熵为何是形状的“本质指纹”？

要理解这个方法的威力，我们必须先抛开对形状的直观像素或轮廓认知，进入一个更抽象的层面：网络与熵。传统方法大多直接处理坐标点，但BLS熵另辟蹊径，它先将一个形状转化为一个特殊的网络，再从这个网络中提取出稳定不变的特征。

2.1 从形状边界到复杂网络

第一步是获取形状的纯净边界。给定一个二值化图像（物体为白色，背景为黑色），我们通过轮廓追踪算法（如OpenCV中的findContours）提取出所有边界像素点的有序序列。假设一个形状有N个边界点，我们将每个边界点视为网络中的一个“节点”。

关键操作：对于边界上的每一个像素点i，我们将其与边界上所有其他N-1个像素点j用“虚拟线段”连接起来。这样，以点i为中心，就形成了一个星状网络，拥有N-1条“分支”。每一条分支的长度，就是点i到点j的欧氏距离。请注意，这个网络是动态的——我们为每一个边界点i都构建这样一个专属的星状网络。这个过程的直观理解是，从形状边界上每一个“观察点”出发，去测量它到形状轮廓上所有其他点的距离分布。这个分布，完整编码了该点在形状中的几何上下文信息。

2.2 BLS熵的计算与物理意义

对于一个以点i为中心的星状网络，我们有了N-1条分支长度{L_i1, L_i2, ..., L_i(N-1)}。BLS熵的计算公式如下：

1. 计算分支长度比例：首先，计算所有分支的总长度L_total = Σ L_ij。然后，计算每条分支的长度占总长度的比例：p_ij = L_ij / L_total。这步归一化确保了熵值不受形状绝对大小的影响。
2. 计算熵值：接着，我们计算这个比例分布的香农熵，并对其进行归一化：S_i = - Σ (p_ij * log(p_ij)) / log(N-1)

这个熵值S_i揭示了什么？它衡量的是从点i出发，到形状轮廓上所有其他点的距离分布的均匀性。

• 如果S_i接近1：意味着从点i到其他所有点的距离几乎相等。想象一个正圆的圆心（虽然圆心不在边界上，但有助于理解），或者一个规则多边形边界上某些对称点，其距离分布会非常均匀。
• 如果S_i接近0：意味着距离分布极度不均，可能从点i到某几个点特别近，而到其他点特别远。这通常发生在形状的尖角、凸起或凹陷等特征显著的位置。

通过遍历形状边界上的每一个像素点，我们就能得到一个熵值序列E = [S_1, S_2, ..., S_N]，这就是该形状的“BLS熵轮廓”。这个轮廓可以看作形状的一维“签名”，它对旋转（起点变化）是循环平移的，对缩放是归一化的，对平移则是完全不变的。

2.3 构建三维特征向量：从轮廓到特征

直接使用整个熵轮廓（一个高维向量）进行聚类计算量大且可能包含冗余。论文的精髓在于，从这个轮廓中提炼出三个核心统计特征，构成一个三维向量[F1, F2, F3]：

1. F1: 整体熵值 (S)：这是整个熵轮廓的均值，即F1 = mean(E)。它代表了形状整体上的“分支均匀性”。一个更接近圆形的、光滑的形状，其F1值会更高；而一个充满尖刺、凹凸不平的形状，其F1值会更低。
2. F2: 相似性均值 (μ)：这里引入了一个中间量γ_ij，它定义为形状i的熵轮廓E_i与形状j的熵轮廓E_j经过循环移位对齐后的最大相关系数。γ_ij衡量了两个形状之间的相似度。对于形状i，我们计算它与数据集中所有形状（包括自身）的γ_ij，然后取均值得到F2 = mean(γ_ij)。这个特征反映了该形状与数据集中其他形状的平均相似程度。
3. F3: 相似性方差 (σ)：同样基于γ_ij，计算其方差F3 = var(γ_ij)。方差大，说明该形状与某些形状很相似，与另一些则差异很大，可能属于特征鲜明的类别；方差小，则说明它与大多数形状的相似度都差不多，可能处于类别边界或是非常常见的形状。

为什么是这三个特征？这是一个非常巧妙的设计。F1捕捉形状的内在均匀性，F2和F3则通过与该形状所在数据集的交互，捕捉其相对的独特性和一致性。它们共同构成了一个既包含绝对几何属性，又包含相对分布信息的紧凑描述符。论文后续实验也表明，增加更多特征（如γ_ij的最大最小值）对性能提升微乎其微，证明了这三个特征的信息充分性。

2.4 t-SNE与k-means的角色：可视化与分组

得到三维特征向量后，我们面临一个典型的无监督学习问题。直接在三/高维空间进行聚类，其几何关系不直观，且距离度量可能失效。

• t-SNE的作用：非线性降维与可视化。t-SNE的核心目标是保持局部结构。它将高维空间中点与点之间的相似性（概率分布）尽可能地在低维（通常是2D或3D）空间中重现。在我们的场景中，它把每个形状的三维特征向量[F1, F2, F3]映射到一个二维点(x, y)。在二维图上，几何结构相似的三维向量点会聚集在一起，而不同的则会分开。这步不仅为了可视化，更重要的是为后续的k-means聚类提供了一个结构更清晰、更易于分离的数据表示。论文中的对比实验（有/无t-SNE）清晰表明，t-SNE能显著提升聚类效果，尤其是当形状类别本身区分度不高时。
• k-means的作用：划定边界。当t-SNE将数据点清晰地呈现在二维平面上后，k-means这个经典的划分式聚类算法就可以大显身手了。它根据我们预设的类别数k，寻找k个中心点，并将每个数据点分配给最近的中心点，迭代优化直至中心点稳定。在t-SNE预处理后的空间中，由于相似形状已经自然聚集，k-means能更准确、更稳定地划分出聚类边界。

3. 完整实现流程与实操要点

理解了原理，我们来看如何动手实现这一整套流程。我将以Python为主要工具，结合常用的科学计算和图像处理库，分步拆解。

3.1 环境准备与数据获取

首先，确保你的Python环境安装了必要的库。我们将主要依赖numpy,opencv-python,scikit-learn,matplotlib和scipy。

pip install numpy opencv-python scikit-learn matplotlib scipy

数据集准备： MPEG-7和Swedish Leaf是学术界广泛使用的基准数据集，通常可以从大学或研究机构的网站下载。Heptagon数据集在论文中是通过算法生成的。为了复现，我们可以直接使用论文中引用的Stack Overflow算法来生成自己的不规则七边形数据集。这里的关键是，无论数据来源如何，我们需要将其统一处理为二值图像（单通道，0为背景，255或1为前景）。

3.2 核心步骤一：计算BLS熵轮廓

这是整个流程中最核心、计算量最大的一步。我们需要为一张二值图像计算其边界上每个点的熵值。

import cv2 import numpy as np from scipy.spatial.distance import cdist from scipy import stats def compute_bls_entropy_profile(binary_image): """ 计算二值图像的BLS熵轮廓。 参数: binary_image: 二维numpy数组，值为0（背景）和255或1（前景）。 返回: entropy_profile: 一维numpy数组，形状的BLS熵轮廓。 contour: 边界点坐标数组。 """ # 1. 提取轮廓 # 使用cv2.RETR_EXTERNAL只取最外层轮廓，cv2.CHAIN_APPROX_NONE存储所有点 contours, _ = cv2.findContours(binary_image.astype(np.uint8), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE) if len(contours) == 0: raise ValueError("未在图像中找到轮廓。") # 取最长的轮廓（假设只有一个主要形状） contour = contours[0].squeeze() # 形状为 (N, 2) N = len(contour) entropy_profile = np.zeros(N) # 2. 为轮廓上每一个点计算熵 for i in range(N): # 点i到轮廓上所有点的欧氏距离 # 使用cdist高效计算所有 pairwise 距离 distances = cdist([contour[i]], contour).squeeze() # 形状 (N,) # 注意：点到自身的距离为0，在计算比例时应排除或处理 # 在BLS熵的原始定义中，是连接到所有其他点，因此我们排除自身 mask = np.ones(N, dtype=bool) mask[i] = False branch_lengths = distances[mask] # 计算分支长度比例 total_length = branch_lengths.sum() if total_length == 0: p = np.ones_like(branch_lengths) / len(branch_lengths) # 避免除零，均匀分布 else: p = branch_lengths / total_length # 计算香农熵并归一化 # 加一个小量防止log(0) epsilon = 1e-10 entropy = -np.sum(p * np.log(p + epsilon)) normalized_entropy = entropy / np.log(len(branch_lengths)) entropy_profile[i] = normalized_entropy return entropy_profile, contour

实操要点与避坑指南：

1. 轮廓提取的稳定性：cv2.findContours对输入图像类型敏感。确保传入uint8类型。参数cv2.CHAIN_APPROX_NONE至关重要，它保留了轮廓上所有的点。如果使用cv2.CHAIN_APPROX_SIMPLE，会压缩水平、垂直、对角方向的线段，丢失大量边界点信息，导致熵轮廓失真。
2. 计算效率：上述代码的双重循环（对每个点i计算到所有点j的距离）复杂度是O(N^2)，对于边界点很多（N很大）的形状会非常慢。优化策略：在循环外一次性计算所有点对之间的距离矩阵dist_matrix = cdist(contour, contour)，然后在循环中直接切片dist_matrix[i, :]。这样将复杂度从O(N^3)降低到O(N^2)。
3. 边界处理与噪声：实际图像可能存在噪声，导致轮廓不平滑或有小毛刺。这会在局部产生异常的距离，影响熵值。建议在计算熵轮廓前，先对轮廓进行适当的平滑处理，例如使用高斯滤波或轮廓近似（cv2.approxPolyDP），但要注意过度平滑会损失形状细节。

3.3 核心步骤二：构建三维特征向量

得到熵轮廓后，我们需要计算形状间的相似度γ_ij，并最终得到特征向量。

def compute_similarity_gamma(profile_i, profile_j): """ 计算两个形状熵轮廓之间的相似度 γ_ij。 通过循环移位找到最大相关系数。 """ n = len(profile_i) max_corr = -1 # 对 profile_j 进行循环移位 for shift in range(n): shifted_j = np.roll(profile_j, shift) # 计算皮尔逊相关系数 corr, _ = stats.pearsonr(profile_i, shifted_j) # 由于浮点误差，corr可能略大于1或小于-1，需要裁剪 corr = np.clip(corr, -1.0, 1.0) if corr > max_corr: max_corr = corr # γ_ij 定义为最大相关系数，论文中似乎直接使用相关系数，可能未取绝对值。 # 根据论文公式(2)，是 max{corr(...)}，我们直接返回最大值。 return max_corr def extract_feature_vector(entropy_profile, all_profiles): """ 为单个形状提取三维特征向量 [F1, F2, F3]。 参数: entropy_profile: 当前形状的熵轮廓。 all_profiles: 列表，包含数据集中所有形状的熵轮廓。 返回: feature_vector: 三维numpy数组 [F1, F2, F3]。 """ # F1: 熵轮廓的均值 F1 = np.mean(entropy_profile) # 计算当前形状与数据集中所有形状（包括自身）的相似度 gamma_values = [] for profile_j in all_profiles: gamma = compute_similarity_gamma(entropy_profile, profile_j) gamma_values.append(gamma) gamma_values = np.array(gamma_values) # F2: γ_ij 的均值 F2 = np.mean(gamma_values) # F3: γ_ij 的方差 F3 = np.var(gamma_values) return np.array([F1, F2, F3])

注意事项：

1. 计算复杂度警告：compute_similarity_gamma函数包含一个内层循环（移位），而extract_feature_vector需要对数据集中的每对形状调用它。对于包含M个形状的数据集，计算所有γ_ij的复杂度是 O(M^2 * N^2)，这是整个流程中最耗时的部分。对于大型数据集，这是不可行的。论文中可能采用了优化策略，或者在小规模数据集上运行。在实际应用中，必须考虑优化，例如：
   • 使用快速傅里叶变换计算循环互相关，可以将复杂度从 O(N^2) 降至 O(N log N)。
   • 对熵轮廓进行下采样，减少N。
   • 使用近似最近邻搜索或基于采样的方法来估算相似度。
2. γ_ij的对称性：理论上γ_ij应该等于γ_ji，因为相关系数运算满足交换律。但在实际数值计算中，由于浮点误差和移位操作，可能会有微小差异。确保你的compute_similarity_gamma函数是确定性的，或者考虑计算一次后存储为上三角矩阵。

3.4 核心步骤三：t-SNE降维与k-means聚类

当为数据集中所有形状都提取了三维特征向量后，我们就可以进行降维和聚类了。

from sklearn.manifold import TSNE from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler import matplotlib.pyplot as plt def cluster_shapes(feature_vectors, n_clusters, perplexity=30, random_state=42): """ 对特征向量进行t-SNE降维和k-means聚类。 参数: feature_vectors: 二维numpy数组，形状为 (M, 3)，M是形状数量。 n_clusters: k-means的聚类数量。 perplexity: t-SNE的困惑度参数。 random_state: 随机种子，保证结果可复现。 返回: labels: 每个形状的聚类标签。 tsne_embeddings: t-SNE降维后的二维坐标。 kmeans_model: 训练好的k-means模型。 """ # 1. 标准化特征 (可选，但通常是个好习惯，尤其是对k-means) scaler = StandardScaler() features_scaled = scaler.fit_transform(feature_vectors) # 2. t-SNE降维 tsne = TSNE(n_components=2, perplexity=perplexity, random_state=random_state, init='random', learning_rate='auto') tsne_embeddings = tsne.fit_transform(features_scaled) # 注意：t-SNE对超参数敏感，perplexity通常建议在5到50之间，需要根据数据调整。 # 3. k-means聚类 kmeans = KMeans(n_clusters=n_clusters, random_state=random_state, n_init='auto') labels = kmeans.fit_predict(tsne_embeddings) # 在降维后的2D空间聚类 # 4. 可视化 plt.figure(figsize=(10, 8)) scatter = plt.scatter(tsne_embeddings[:, 0], tsne_embeddings[:, 1], c=labels, cmap='tab20', s=20, alpha=0.7) plt.colorbar(scatter, label='Cluster ID') plt.title(f't-SNE Visualization with K-means Clustering (k={n_clusters})') plt.xlabel('t-SNE dimension 1') plt.ylabel('t-SNE dimension 2') plt.tight_layout() plt.show() return labels, tsne_embeddings, kmeans

关键参数解析与调优经验：

1. t-SNE的perplexity：这是t-SNE最重要的参数。它大致表示每个点考虑多少邻居。值太小（如5），t-SNE会过度关注局部结构，形成许多分散的小簇；值太大（如50），则会更多关注全局结构，可能使不同类别混在一起。对于形状数量在1000左右的数据集，论文中使用30是一个合理的起点。我的经验是：在最终分析前，用几个不同的perplexity值（例如5, 30, 50）运行t-SNE，观察可视化结果的稳定性。稳定的聚类结构在不同perplexity下应该大致保持。
2. 特征标准化：在应用t-SNE和k-means之前，对三维特征进行标准化（减去均值，除以标准差）是强烈推荐的。因为F1、F2、F3的量纲和数值范围可能不同，标准化可以防止某个特征因其数值大而主导降维和聚类过程。
3. k-means的k值选择：这是无监督聚类的经典难题。论文中使用逻辑曲线拟合来确定最佳k值，这是一个非常聪明的方法。具体做法是：遍历一系列k值（如1到100），对每个k运行聚类，然后计算一个“聚类得分”（例如，类内平均相似度）。随着k增加，得分会上升并逐渐饱和。拟合一个逻辑增长曲线，取达到饱和值99%对应的k作为最佳值。我们也可以结合手肘法（看类内距离平方和下降的拐点）和轮廓系数来综合判断。
4. 聚类是在2D空间进行的：请注意，我们是在t-SNE降维后的二维空间进行k-means聚类，而不是原始的三维特征空间。这是因为t-SNE后的空间局部几何关系更清晰，更适合基于欧氏距离的k-means。

3.5 性能评估与对比实验

为了验证方法的有效性，我们需要像论文中一样，设计对比实验。

def evaluate_clustering_score(feature_vectors, labels): """ 计算聚类得分。这里采用一种简化的类内平均相似度。 更严谨的做法是使用轮廓系数(Silhouette Score)或Calinski-Harabasz指数。 """ from sklearn.metrics import silhouette_score, calinski_harabasz_score # 使用轮廓系数，值在[-1,1]之间，越大越好 score_sil = silhouette_score(feature_vectors, labels) # 使用Calinski-Harabasz指数，值越大表示聚类越紧密、分离越好 score_ch = calinski_harabasz_score(feature_vectors, labels) return score_sil, score_ch def compare_with_fourier_descriptors(images, method='bls'): """ 与傅里叶描述符方法进行对比。 简化版：计算傅里叶描述符，进行k-means，比较轮廓系数。 """ all_features = [] for img in images: if method == 'fourier': # 提取轮廓，计算傅里叶描述符（例如前10个系数） contour = extract_contour(img) fourier_desc = compute_fourier_descriptors(contour, num_descriptors=10) all_features.append(fourier_desc) elif method == 'bls': profile, _ = compute_bls_entropy_profile(img) # 这里需要所有轮廓来计算F2,F3，简化起见，可以先用所有数据计算一次gamma矩阵。 # 为演示，我们假设有一个预计算好的gamma矩阵 `gamma_matrix` # F1 = mean(profile) # F2 = mean(gamma_matrix[i, :]) # F3 = var(gamma_matrix[i, :]) # all_features.append([F1, F2, F3]) pass all_features = np.array(all_features) # 标准化 from sklearn.preprocessing import StandardScaler scaler = StandardScaler() features_scaled = scaler.fit_transform(all_features) # 使用最佳k值或一系列k值进行聚类和评估 # ... (这里可以循环k值，计算评估指标并绘图，如论文中的图8、9、10)

实现对比实验的要点：

1. 傅里叶描述符的实现：需要实现轮廓的傅里叶变换。通常将轮廓点表示为复数序列z(t) = x(t) + i*y(t)，然后进行离散傅里叶变换，取前几个低频系数的幅度值作为描述符。注意需要对轮廓进行重采样，使得点数一致，并对起始点进行对齐（相位归一化）。
2. 泽尔尼克矩的实现：OpenCV提供了cv2.moments和cv2.HuMoments，但泽尔尼克矩需要自己实现或使用第三方库（如mahotas或pymoment）。计算复杂度较高。
3. 评估指标的选择：论文中使用的是自定义的“聚类得分”（类内平均相似度）。在学术界，更通用的内部评估指标是轮廓系数和Calinski-Harabasz指数。外部评估指标（如果有真实标签）则可以使用调整兰德指数或归一化互信息。建议在报告中同时给出多种指标，使对比更全面。

4. 实战中的挑战、调优与问题排查

将论文方法付诸实践时，你会遇到一系列在理想化描述中不会提及的挑战。以下是我在复现和类似项目中总结的经验和解决方案。

4.1 计算效率：从理论到实践的鸿沟

正如前面提到的，BLS熵轮廓和γ_ij的完整计算复杂度是灾难性的。对于一个有1000个形状，每个形状轮廓有500个点的数据集，计算所有γ_ij的理论操作次数约为1000^2 * 500^2 = 2.5e11次距离和相关系数计算，这几乎无法完成。

优化策略实录：

1. 轮廓下采样：这是最直接有效的方法。原始轮廓可能包含数百甚至上千个点，但形状的主要几何信息通常可以由少得多的点来近似。使用cv2.approxPolyDP函数进行轮廓近似，在保持形状整体外观的前提下，将点数减少到50-100个。这可以将计算量降低2个数量级，而对最终特征向量的影响通常很小。
2. 加速相似度计算：计算两个序列的循环互相关最大值，可以利用卷积定理通过FFT加速。scipy.signal库中的correlate函数（设置mode='full'并寻找最大值）比手动循环移位快得多。
3. 并行化：计算γ_ij矩阵是一个典型的令人尴尬的并行问题。每个γ_ij的计算独立于其他。可以使用Python的multiprocessing库或joblib进行多进程并行计算，充分利用多核CPU。
4. 近似与抽样：如果数据集极大，可以考虑不计算全对的γ_ij。例如，只为每个形状计算它与一个随机子集（如10%）的其他形状的相似度，来近似估计F2和F3。或者，使用基于锚点的方法。

4.2 参数敏感性与稳定性

1. t-SNE的随机性：t-SNE的优化过程是随机的，每次运行结果可能略有不同。random_state参数固定可以保证单次实验可复现，但为了评估方法的鲁棒性，必须多次运行（例如10次）并统计聚类指标的平均值和方差。如果结果波动很大，说明降维结果不稳定，可能需要调整perplexity或learning_rate，或者考虑使用UMAP等更稳定的降维方法作为对比。
2. k-means的初始化和局部最优：KMeans的n_init参数（默认已改为'auto'）控制着以不同质心种子运行的次数，最终取结果最好的一个。对于形状复杂的数据，建议显式设置n_init=10或更高，以避免陷入局部最优解。
3. 最佳聚类数k的确定：论文的逻辑曲线拟合法很有启发性，但实现起来需要定义合理的“聚类得分”。我推荐结合以下方法综合判断：
   • 手肘法：绘制不同k值下k-means的惯性（inertia，即样本到其最近聚类中心的距离平方和）。惯性下降的拐点就是“手肘”。
   • 轮廓系数法：计算每个k值对应的平均轮廓系数，取最大值对应的k。
   • 间隔统计法：比较实际数据的惯性与其在参考分布（如均匀分布）下的惯性之差。 将三种方法的结果画在同一张图上，寻找共识区间。

4.3 常见问题与排查清单

4.4 超越论文：方法的扩展与改进思考

论文给出了一个坚实的方法框架，但在实际应用中，我们可以思考如何使其更强健、更高效。

1. 特征向量的增强：论文实验表明增加F4（最大γ_ij）和F5（最小γ_ij）收益不大。但我们可以尝试其他统计量，例如熵轮廓的偏度和峰度，或者将熵轮廓本身通过自编码器压缩成一个低维向量。也可以引入多尺度思想，在不同高斯平滑级别下计算BLS熵，形成多尺度特征。
2. 聚类流程的改进：t-SNE + k-means是经典流程，但并非唯一。可以尝试：
   • 直接在高维空间聚类：使用DBSCAN或HDBSCAN这类密度聚类算法，它们不需要指定k，且能发现任意形状的簇。但前提是特征空间中的距离度量要有意义。
   • 谱聚类：基于特征向量构建相似度矩阵（可以用γ_ij矩阵本身！），然后在其上进行谱聚类，对于非凸形状的簇可能效果更好。
3. 处理更大规模数据：对于海量形状数据集，上述方法需要进一步革新。可以考虑：
   • 学习型特征提取：训练一个卷积自编码器或Siamese网络，直接从形状图像中学习一个紧致的、对变换不变的特征表示，完全绕过手工设计特征和昂贵的相似度计算。
   • 近似最近邻搜索：在需要检索相似形状时，使用Faiss、Annoy等库来快速查找高维空间中的近似最近邻，替代全对计算。

这个基于BLS熵的形状聚类方法，其魅力在于它从一个非常新颖的角度（网络熵）来刻画形状，并且用极简的特征向量取得了不俗的效果。它提醒我们，在深度学习席卷一切的时代，精巧的数学模型和特征工程依然有其不可替代的价值，尤其是在数据量有限、可解释性要求高的场景下。复现它的过程，本身就是一次对形状本质、特征不变性和无监督学习的深度思考。希望这篇详尽的拆解，能帮助你不仅理解这篇论文，更能掌握将一篇学术论文转化为可运行、可改进、可应用的实践项目的能力。