基于密度峰值搜索的脑纤维快速聚类算法

(浙江工业大学 信息工程学院，浙江 杭州 310023)

弥散张量成像(Diffusion tensor imaging, DTI)是目前活体显示神经纤维走向的唯一方法，能够对大脑内组织微观结构进行定量分析，在脑外科手术导航及神经类疾病的研究中起重要作用。大脑白质纤维可视化一直是医学可视化领域的热门话题[1]。基于DTI的全脑纤维束追踪结果是个非常庞大的多维数据集，这些纤维错综复杂且相互遮挡，很难直接通过肉眼观测脑纤维结构，很大程度上限制了纤维束追踪在临床中的应用。为了解决该问题，通常根据先验知识人工划定对应的兴趣区，再对该区域相关的脑纤维进行可视与分析[2]，然而这种方法过度依赖解剖学结构知识且极易受到使用者偏见的影响。

纤维聚类技术通过对人类大脑中几何分布结构相似的纤维进行归类，从而更好地展示纤维束之间的空间关系，为提升对大脑组织结构和整体连通性的感知提供可行性方案。合理的纤维聚类方式可以在很大程度上消除不同类别纤维之间的干扰，更好地描述纤维束的走向，从而帮助理解相互混合纤维的空间联系。早期研究认为起点和终点相近的纤维是相似的，以两条纤维的起点之间的距离以及终点之间的距离作为纤维的相似度。但是并不是所有的纤维都有相同的起点和终点，而且这种方法忽略了纤维的形状信息。为了解决这一问题，Klein等提出了基于网格的纤维相似度测量方法[3]，用网格将纤维分隔成独立的单元，把构成纤维的数据点以权重的形式分配给每一个单元，通过计算单元对应的权重来衡量纤维的相似度。以纤维A上的所有点到纤维B的最短距离的平均值作为纤维相似度的平均最近点距离[4]和以纤维A上的所有点到纤维B的最短距离中的最大值作为纤维相似度的豪斯多夫距离[5]是最常用的纤维相似度衡量方法。图集算法[6]首先由医学专家在大脑皮层中指定一系列兴趣区，然后在兴趣区中随机生成种子点获得纤维样本，最后将得到的特定纤维样本作为集群中心进行聚类算法。K-means算法[7]以类内距离为依据，迭代地更新集群中心，是一种典型的划分聚类方法。主成分分析(PCA)以最少的信息损失量实现多维空间向各点对低维空间的投影将弯曲的纤维轨迹映射成PCA空间中的一个点，然后再对PCA空间中的点进行聚类分析[8-9]。谱聚类算法[10]首先根据纤维之间的相似度形成关联矩阵，通过计算该关联矩阵的特征向量和特征值，根据计算得到的特征向量进行聚类。凝聚聚类将每条纤维视为一个单独的集群，连续的合并最近的一对集群直到所有的纤维都在一个集群中[11]。分裂聚类将整个纤维集视为一个集群，然后逐渐细分为越来越小的集群，直到每条纤维自成一个集群或者达到了某个终止条件[12]。DBSCAN会设定一个密度阈值，将密度小于该阈值的纤维作为噪声消除，把整个连续区域分成多个高密度区域进行聚类[13-14]。OPTICS算法首先随机的选择一个对象作为种子点，将种子点的近邻插入队列之中，然后递归的从队列中选择种子点直至队列为空，此时再随机的选择一个对象作为种子点，直至所有对象都被计算[15]。

1 纤维相似度计算方法

在纤维跟踪成像之后，可以得到由一系列连续的点组成的纤维，然而由于纤维本身的长度不同以及纤维追踪算法的差别，组成每条纤维序列中的点数目很难保持一致。基于这种情况，采用动态时间规整算法(Dynamic time warping, DTW)计算纤维之间的相似度。

1.1 动态时间规整

DTW主要应用于语音识别，用来衡量两个语音信号的相似度，其原理是要找到两条时间序列之间的最短折叠路径[16]。将时间序列沿着时间轴进行拉伸和折叠等操作使其扭曲在一起，设现有时间序列X和Y，它们的长度分别为m和n，即

X=(x1,x2,…,xi,…,xn)

(1)

Y=(y1,y2,…,yi,…,ym)

(2)

设K为折叠路径W的长度，则折叠路径W为

W=w1,w2,…,wk,…,wK

(3)

那么

max(m,n)

(4)

式中：W的元素为时间序列X，Y中时间节点的连接情况，即wk(i,j)。折叠路径W遵循的3 个约束条件为

1) 边界约束条件：w1=(1,1)和wK=(m,n)。这要求折叠路径从第一个时间点开始，到最后一个时间点结束。

2) 单调性约束条件：对于任意的wk=(i,j)和wk+1=(i′,j′)有i′-i≥0和j′-j≥0，这使得折叠路径W上的节点在时间轴上单调增加。

3) 连续性约束条件：对于任意的wk=(i,j)和wk+1=(i′,j′)，有i′-i≤1和j′-j≤1，这使得折叠路径中的每一步都在时间序列上相邻。

然而，实际应用中存在许多满足这3 个约束条件的折叠路径，为了能够更好地衡量两条时间序列之间的差异，选择用所有可能的折叠路径中累计距离最短的折叠路径作为最佳路径，最佳路径如图1所示。最佳路径的累积距离定义为

(5)

式中d(·)为距离函数，定义为

d(wk)=d(i,j)=|xi-yj|

(6)

图1 最佳折叠路径图Fig.1 Optimal warping path

最佳折叠路径可以通过动态算法[17]得到：首先计算一个维数为m×n的成本矩阵D，成本矩阵中的每个元素D(i,j)由距离d(i,j)与相邻的元素的和递归地赋值，其计算式为

D(i,j)=d(i,j)+min{D(i-1,j-1),
D(i,j-1),D(i-1,j)}

(7)

当建立成本矩阵后，可以从终点D(m,n)反向地寻找一条累积距离最短的路径。为了达到这个目的，需要对成本矩阵沿着左、下和对角线3 个方向进行贪婪搜索。这3 个相邻元素中的最小元素将被加入折叠路径，并且搜索过程将从该元素继续进行直至搜索到D(1,1)。如果折叠路径经过成本矩阵中的元素D(i,j)，则表示时间序列X上的第i个点与时间序列Y上第j个点相连。由于一条时间序列单个点可能对应着另一条时间序列上的多个点，因此，该方法可以衡量不同长度时间序列之间的相似度。

1.2 DTW测量纤维相似度

(8)

通过使用该距离函数，可以将三维空间的纤维相似性问题转化为时间序列相似性问题。最佳折叠路径W(w1,…,wk,…,wK)由动态规划算法计算得出，其中K为路径步长且d(wk)=d(Xi,Yj)。

为了降低由不同纤维的长度差异对纤维相似度测量的影响，定义纤维相似度距离为

(9)

2 快速密度峰值搜索算法

脑纤维聚类算法不受大脑皮层功能区域的影响，检测每束纤维的中心、分配等纤维束空间结构信息。在完成纤维相似度计算后，每条纤维可以被看作度量空间中的一个对象。该方法侧重于测量纤维束的空间信息，而不是对脑神经网络进行图像分析。

快速密度峰值搜索算法是一种基于距离和密度的混合聚类算法，这种方法信息源仅基于数据点之间的距离，能对任意形状的数据进行聚类并自动得到正确的集群数量[18]。使用该聚类算法的数据集合需满足两个条件：1) 集群中心的局部密度高于周围数据点的局部密度；2) 两个集群中心之间的距离较大。显然，脑纤维数据满足这两个假设条件。

对于每条纤维i，需要计算两个变量，纤维的局部密度ρi和纤维与更高密度纤维之间的最小距离δi。这两个变量仅取决于两条纤维之间的距离dij，将纤维i的局部密度定义为

(10)

Rodriguez等[15]基于点数据的聚类经验认为使平均点密度为数据总点数1%～2%的密度半径dc是合适的。但是经过实验发现：这种基于点数据的密度半径设置方式并不能很好的应用于脑纤维聚类，因为dc的值不仅与样本数量有关，而且还与样本的分布情况有关。根据脑纤维的结构特性以及降低计算时间成本的需求，笔者提出了一种新的计算密度半径的方法，该方法既能快速得到密度半径，又能够根据脑纤维结构特性自动调整密度半径，该密度半径计算方法为

1) 随机选择样本空间中的一部分纤维作为种子点，采样比例为r1,种子点数量为t。

2) 设定阈值r2，令i=Sum×r2(Sum为总纤维数，i向上取整)，分别计算每个种子点与其第i近的纤维距离di。

3) 定义密度半径dc为所有种子点的di的平均值，其计算式为

(11)

实验表明：当采样比例r1=4%，阈值r2=0.8%时可达到较好的可视化效果，笔者方法通过从纤维集中随机取样计算密度半径来提高计算效率。实验证明：笔者方法对于不同的纤维集都能获得适当的密度半径dc，从而得到较好的聚类结果。

定义最小距离δi为纤维i到局部密度高于自身的其他纤维之间的最小距离，其计算式为

(12)

对于局部ρτ密度为最大的纤维τ，定义δτ=max(dij)。ρi当取到极大值时其对应的δi将远大于一般的纤维，因此以最小距离值为纵轴、以局部密度为横轴将纤维投影到密度距离决策图中来交互地选择纤维中心。用户可根据δi值和局部密度都较大的原则，直观地对纤维是否可能为纤维中心进行判断，便于选择纤维中心进行聚类。密度—距离决策图如图2所示，其中灰色点表示用户选择的纤维中心，剩余纤维由黑色点表示。在选定纤维中心之后，所有未被选择的纤维将被分配到局部密度高于自身的最近纤维。

图2 密度 —距离决策图Fig.2 The density-distance decision graph

3 连续聚类框架

图3 算法框架图Fig.3 The algorithm framework

3.1 初 始 化

初始化是指第一次计算集群的参数模型的过程。该过程当读入纤维数量达到用户设定的阈值N时被触发。在初始化过程进行时，程序将暂停纤维数据的读取,同时对初始化数据计算密度半径dc，并使用快速密度峰值搜索算法进行聚类。由于在连续聚类框架中无法交互地选择纤维中心，这里设置判定变量γi为

γi=ρj×δi

(13)

当γi>t×γmax时，将纤维选择为纤维中心，实验证明阈值t取0.1时，能达到较好的聚类效果。

3.2 标记新纤维路径

3.3 模型更新的触发条件

随着聚类过程的进行，纤维中心可能会发生漂移，这将使得参数模型不再满足聚类要求，因此需要在适当的时候对参数模型进行更新操作。触发模型更新的条件主要有：1) 缓存容器R的容量限制，即当容器R中的纤维数量到达某一阈值时执行参数模型更新；2) 基于Page-Hinkley test(PHT)，PHT可以通过分析相关性信息来发现纤维中心的漂移情况。由于第1 种触发条件较为简单，这里着重说明第2 种触发条件。

3.4 相关性信息

(14)

3.5 Page-Hinkley测试

PHT算法可以用分析ρτ序列来找出偏移的纤维路径。将相关性信息的测量值和平均值之间的累计差值设为变量Ut，则

(15)

(16)

式中δ表示为容忍误差。为了找到累计差值Ut突变的时刻，首先计算|Ut|的最大值mT，再判断(mT-Ut)>λ，λ为探测阈值。当该条件成立时表明当前纤维中心发生了偏移，这将触发模型更新程序，此时纤维相关性信息将被清零，时间索引τ将被重置为1。

3.6 模型更新

4 实验结果

4.1 实验环境

实验采用Matlab作为平台架构，配置：Intel(R) Core(TM) i7-4770K CPU @ 3.50 GHz，16 GB RAM，1 T SATA 硬盘。分别使用DBSCAN聚类方法和连续聚类方法对同一纤维数据集进行聚类，最终通过实验结果对比证明连续聚类方法在得到良好结果的同时能显著减少整体计算时间。

4.2 实验数据

为体现算法的普遍适用性，实验分别对从临床数据中通过散布矩阵筛选得到的3 组形态不同的纤维数据集进行聚类分析。纤维数据集A由1 544 条纤维组成，其中最长的纤维由198 个数据点构成，最短的纤维由110 个数据点构成，整个数据集共计209 597 个数据点。纤维数据集B由1 354 条纤维组成，其中最长的纤维由276 个数据点构成，最短的纤维由70 个数据点构成，整个数据集共计183 977 个数据点。纤维数据集C由1 045 条纤维组成，其中最长的纤维由312 个数据点构成，最短的纤维由53 个数据点构成，整个数据集共计151 986 个数据点。

4.3 实验结果和分析

笔者方法的参数设置：采样率r1=4%，距离阈值r2=0.8%，纤维中心阈值t=0.1，容器容量阈值R=50，容忍误差δ=0.01，探测阈值λ=0.25。图4(a～c)中：左侧是未使用聚类算法的纤维结构，右侧是使用连续框架聚类算法得到的结果，当完成纤维聚类后给不同的纤维集群分配不同的灰度值。通过实验证明连续框架聚类算法能有效地分类纤维，结构相似的纤维被聚类为纤维集群并用同一种灰度值显示，各纤维集群之间通过不同的灰度值清晰的区分。

图4 纤维聚类结果Fig.4 The fiber clustering results

实验证明连续框架聚类算法所需时间明显低于一般的DBSCAN聚类方法，如表1所示。连续框架聚类算法的计算量与纤维数量以及纤维结构复杂度都有关，纤维结构越复杂触发模型更新的次数就越多额外计算量也就越大，在不计模型更新整体计算量与纤维数呈线性关系。常用的聚类方法的整体计算量则与纤维数的平方呈线性关系。设纤维数量为n，构成纤维的体素点数量为m，聚类数量为l，则DBSCAN聚类方法的时间复杂度为O(n2m2)，连续框架聚类算法的时间复杂度为O(ln2m2)，因而当纤维数量越多时，连续框架聚类算法计算用时减少地越明显。

表1 不同聚类方法计算效率对比

Table 1 The comparison of time efficiency using different clustering methods

纤维集连续聚类用时/sDBSCAN用时/sA1446.5720568.23B1452.5914023.67C1350.4810136.44

5 结 论

对快速密度峰值聚类算法进行改进，采用动态时间规整算法测量纤维之间的相似度，相比于常用的最近点平均距离和豪斯多夫距离兼顾了纤维整体形态上的相似性，能够更准确地衡量相似性，提高聚类准确度。现有的聚类方法需要计算每对纤维之间的相似度，这将耗费大量的计算成本。相比于常规聚类方法，采用连续聚类框架进行脑纤维聚类在得到满足可视分析需求结果的同时，显著减少纤维相似度计算量，进一步降低整体聚类时间。