基于密度峰值搜索的腦纖維快速聚類算法

(浙江工業(yè)大學 信息工程學院，浙江 杭州 310023)

彌散張量成像(Diffusion tensor imaging, DTI)是目前活體顯示神經(jīng)纖維走向的唯一方法，能夠對大腦內(nèi)組織微觀結構進行定量分析，在腦外科手術導航及神經(jīng)類疾病的研究中起重要作用。大腦白質(zhì)纖維可視化一直是醫(yī)學可視化領域的熱門話題[1]。基于DTI的全腦纖維束追蹤結果是個非常龐大的多維數(shù)據(jù)集，這些纖維錯綜復雜且相互遮擋，很難直接通過肉眼觀測腦纖維結構，很大程度上限制了纖維束追蹤在臨床中的應用。為了解決該問題，通常根據(jù)先驗知識人工劃定對應的興趣區(qū)，再對該區(qū)域相關的腦纖維進行可視與分析[2]，然而這種方法過度依賴解剖學結構知識且極易受到使用者偏見的影響。

纖維聚類技術通過對人類大腦中幾何分布結構相似的纖維進行歸類，從而更好地展示纖維束之間的空間關系，為提升對大腦組織結構和整體連通性的感知提供可行性方案。合理的纖維聚類方式可以在很大程度上消除不同類別纖維之間的干擾，更好地描述纖維束的走向，從而幫助理解相互混合纖維的空間聯(lián)系。早期研究認為起點和終點相近的纖維是相似的，以兩條纖維的起點之間的距離以及終點之間的距離作為纖維的相似度。但是并不是所有的纖維都有相同的起點和終點，而且這種方法忽略了纖維的形狀信息。為了解決這一問題，Klein等提出了基于網(wǎng)格的纖維相似度測量方法[3]，用網(wǎng)格將纖維分隔成獨立的單元，把構成纖維的數(shù)據(jù)點以權重的形式分配給每一個單元，通過計算單元對應的權重來衡量纖維的相似度。以纖維A上的所有點到纖維B的最短距離的平均值作為纖維相似度的平均最近點距離[4]和以纖維A上的所有點到纖維B的最短距離中的最大值作為纖維相似度的豪斯多夫距離[5]是最常用的纖維相似度衡量方法。圖集算法[6]首先由醫(yī)學專家在大腦皮層中指定一系列興趣區(qū)，然后在興趣區(qū)中隨機生成種子點獲得纖維樣本，最后將得到的特定纖維樣本作為集群中心進行聚類算法。K-means算法[7]以類內(nèi)距離為依據(jù)，迭代地更新集群中心，是一種典型的劃分聚類方法。主成分分析(PCA)以最少的信息損失量實現(xiàn)多維空間向各點對低維空間的投影將彎曲的纖維軌跡映射成PCA空間中的一個點，然后再對PCA空間中的點進行聚類分析[8-9]。譜聚類算法[10]首先根據(jù)纖維之間的相似度形成關聯(lián)矩陣，通過計算該關聯(lián)矩陣的特征向量和特征值，根據(jù)計算得到的特征向量進行聚類。凝聚聚類將每條纖維視為一個單獨的集群，連續(xù)的合并最近的一對集群直到所有的纖維都在一個集群中[11]。分裂聚類將整個纖維集視為一個集群，然后逐漸細分為越來越小的集群，直到每條纖維自成一個集群或者達到了某個終止條件[12]。DBSCAN會設定一個密度閾值，將密度小于該閾值的纖維作為噪聲消除，把整個連續(xù)區(qū)域分成多個高密度區(qū)域進行聚類[13-14]。OPTICS算法首先隨機的選擇一個對象作為種子點，將種子點的近鄰插入隊列之中，然后遞歸的從隊列中選擇種子點直至隊列為空，此時再隨機的選擇一個對象作為種子點，直至所有對象都被計算[15]。

1 纖維相似度計算方法

在纖維跟蹤成像之后，可以得到由一系列連續(xù)的點組成的纖維，然而由于纖維本身的長度不同以及纖維追蹤算法的差別，組成每條纖維序列中的點數(shù)目很難保持一致。基于這種情況，采用動態(tài)時間規(guī)整算法(Dynamic time warping, DTW)計算纖維之間的相似度。

1.1 動態(tài)時間規(guī)整

DTW主要應用于語音識別，用來衡量兩個語音信號的相似度，其原理是要找到兩條時間序列之間的最短折疊路徑[16]。將時間序列沿著時間軸進行拉伸和折疊等操作使其扭曲在一起，設現(xiàn)有時間序列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中時間節(jié)點的連接情況，即wk(i,j)。折疊路徑W遵循的3 個約束條件為

1) 邊界約束條件：w1=(1,1)和wK=(m,n)。這要求折疊路徑從第一個時間點開始，到最后一個時間點結束。

2) 單調(diào)性約束條件：對于任意的wk=(i,j)和wk+1=(i′,j′)有i′-i≥0和j′-j≥0，這使得折疊路徑W上的節(jié)點在時間軸上單調(diào)增加。

3) 連續(xù)性約束條件：對于任意的wk=(i,j)和wk+1=(i′,j′)，有i′-i≤1和j′-j≤1，這使得折疊路徑中的每一步都在時間序列上相鄰。

然而，實際應用中存在許多滿足這3 個約束條件的折疊路徑，為了能夠更好地衡量兩條時間序列之間的差異，選擇用所有可能的折疊路徑中累計距離最短的折疊路徑作為最佳路徑，最佳路徑如圖1所示。最佳路徑的累積距離定義為

(5)

式中d(·)為距離函數(shù)，定義為

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

(6)

圖1 最佳折疊路徑圖Fig.1 Optimal warping path

最佳折疊路徑可以通過動態(tài)算法[17]得到：首先計算一個維數(shù)為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 個相鄰元素中的最小元素將被加入折疊路徑，并且搜索過程將從該元素繼續(xù)進行直至搜索到D(1,1)。如果折疊路徑經(jīng)過成本矩陣中的元素D(i,j)，則表示時間序列X上的第i個點與時間序列Y上第j個點相連。由于一條時間序列單個點可能對應著另一條時間序列上的多個點，因此，該方法可以衡量不同長度時間序列之間的相似度。

1.2 DTW測量纖維相似度

(8)

通過使用該距離函數(shù)，可以將三維空間的纖維相似性問題轉化為時間序列相似性問題。最佳折疊路徑W(w1,…,wk,…,wK)由動態(tài)規(guī)劃算法計算得出，其中K為路徑步長且d(wk)=d(Xi,Yj)。

為了降低由不同纖維的長度差異對纖維相似度測量的影響，定義纖維相似度距離為

(9)

2 快速密度峰值搜索算法

腦纖維聚類算法不受大腦皮層功能區(qū)域的影響，檢測每束纖維的中心、分配等纖維束空間結構信息。在完成纖維相似度計算后，每條纖維可以被看作度量空間中的一個對象。該方法側重于測量纖維束的空間信息，而不是對腦神經(jīng)網(wǎng)絡進行圖像分析。

快速密度峰值搜索算法是一種基于距離和密度的混合聚類算法，這種方法信息源僅基于數(shù)據(jù)點之間的距離，能對任意形狀的數(shù)據(jù)進行聚類并自動得到正確的集群數(shù)量[18]。使用該聚類算法的數(shù)據(jù)集合需滿足兩個條件：1) 集群中心的局部密度高于周圍數(shù)據(jù)點的局部密度；2) 兩個集群中心之間的距離較大。顯然，腦纖維數(shù)據(jù)滿足這兩個假設條件。

對于每條纖維i，需要計算兩個變量，纖維的局部密度ρi和纖維與更高密度纖維之間的最小距離δi。這兩個變量僅取決于兩條纖維之間的距離dij，將纖維i的局部密度定義為

(10)

Rodriguez等[15]基于點數(shù)據(jù)的聚類經(jīng)驗認為使平均點密度為數(shù)據(jù)總點數(shù)1%～2%的密度半徑dc是合適的。但是經(jīng)過實驗發(fā)現(xiàn)：這種基于點數(shù)據(jù)的密度半徑設置方式并不能很好的應用于腦纖維聚類，因為dc的值不僅與樣本數(shù)量有關，而且還與樣本的分布情況有關。根據(jù)腦纖維的結構特性以及降低計算時間成本的需求，筆者提出了一種新的計算密度半徑的方法，該方法既能快速得到密度半徑，又能夠根據(jù)腦纖維結構特性自動調(diào)整密度半徑，該密度半徑計算方法為

1) 隨機選擇樣本空間中的一部分纖維作為種子點，采樣比例為r1,種子點數(shù)量為t。

2) 設定閾值r2，令i=Sum×r2(Sum為總纖維數(shù)，i向上取整)，分別計算每個種子點與其第i近的纖維距離di。

3) 定義密度半徑dc為所有種子點的di的平均值，其計算式為

(11)

實驗表明：當采樣比例r1=4%，閾值r2=0.8%時可達到較好的可視化效果，筆者方法通過從纖維集中隨機取樣計算密度半徑來提高計算效率。實驗證明：筆者方法對于不同的纖維集都能獲得適當?shù)拿芏劝霃絛c，從而得到較好的聚類結果。

定義最小距離δi為纖維i到局部密度高于自身的其他纖維之間的最小距離，其計算式為

(12)

對于局部ρτ密度為最大的纖維τ，定義δτ=max(dij)。ρi當取到極大值時其對應的δi將遠大于一般的纖維，因此以最小距離值為縱軸、以局部密度為橫軸將纖維投影到密度距離決策圖中來交互地選擇纖維中心。用戶可根據(jù)δi值和局部密度都較大的原則，直觀地對纖維是否可能為纖維中心進行判斷，便于選擇纖維中心進行聚類。密度—距離決策圖如圖2所示，其中灰色點表示用戶選擇的纖維中心，剩余纖維由黑色點表示。在選定纖維中心之后，所有未被選擇的纖維將被分配到局部密度高于自身的最近纖維。

圖2 密度 —距離決策圖Fig.2 The density-distance decision graph

3 連續(xù)聚類框架

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

3.1 初 始 化

初始化是指第一次計算集群的參數(shù)模型的過程。該過程當讀入纖維數(shù)量達到用戶設定的閾值N時被觸發(fā)。在初始化過程進行時，程序將暫停纖維數(shù)據(jù)的讀取,同時對初始化數(shù)據(jù)計算密度半徑dc，并使用快速密度峰值搜索算法進行聚類。由于在連續(xù)聚類框架中無法交互地選擇纖維中心，這里設置判定變量γi為

γi=ρj×δi

(13)

當γi>t×γmax時，將纖維選擇為纖維中心，實驗證明閾值t取0.1時，能達到較好的聚類效果。

3.2 標記新纖維路徑

3.3 模型更新的觸發(fā)條件

隨著聚類過程的進行，纖維中心可能會發(fā)生漂移，這將使得參數(shù)模型不再滿足聚類要求，因此需要在適當?shù)臅r候對參數(shù)模型進行更新操作。觸發(fā)模型更新的條件主要有：1) 緩存容器R的容量限制，即當容器R中的纖維數(shù)量到達某一閾值時執(zhí)行參數(shù)模型更新；2) 基于Page-Hinkley test(PHT)，PHT可以通過分析相關性信息來發(fā)現(xiàn)纖維中心的漂移情況。由于第1 種觸發(fā)條件較為簡單，這里著重說明第2 種觸發(fā)條件。

3.4 相關性信息

(14)

3.5 Page-Hinkley測試

PHT算法可以用分析ρτ序列來找出偏移的纖維路徑。將相關性信息的測量值和平均值之間的累計差值設為變量Ut，則

(15)

(16)

式中δ表示為容忍誤差。為了找到累計差值Ut突變的時刻，首先計算|Ut|的最大值mT，再判斷(mT-Ut)>λ，λ為探測閾值。當該條件成立時表明當前纖維中心發(fā)生了偏移，這將觸發(fā)模型更新程序，此時纖維相關性信息將被清零，時間索引τ將被重置為1。

3.6 模型更新

4 實驗結果

4.1 實驗環(huán)境

實驗采用Matlab作為平臺架構，配置：Intel(R) Core(TM) i7-4770K CPU @ 3.50 GHz，16 GB RAM，1 T SATA 硬盤。分別使用DBSCAN聚類方法和連續(xù)聚類方法對同一纖維數(shù)據(jù)集進行聚類，最終通過實驗結果對比證明連續(xù)聚類方法在得到良好結果的同時能顯著減少整體計算時間。

4.2 實驗數(shù)據(jù)

為體現(xiàn)算法的普遍適用性，實驗分別對從臨床數(shù)據(jù)中通過散布矩陣篩選得到的3 組形態(tài)不同的纖維數(shù)據(jù)集進行聚類分析。纖維數(shù)據(jù)集A由1 544 條纖維組成，其中最長的纖維由198 個數(shù)據(jù)點構成，最短的纖維由110 個數(shù)據(jù)點構成，整個數(shù)據(jù)集共計209 597 個數(shù)據(jù)點。纖維數(shù)據(jù)集B由1 354 條纖維組成，其中最長的纖維由276 個數(shù)據(jù)點構成，最短的纖維由70 個數(shù)據(jù)點構成，整個數(shù)據(jù)集共計183 977 個數(shù)據(jù)點。纖維數(shù)據(jù)集C由1 045 條纖維組成，其中最長的纖維由312 個數(shù)據(jù)點構成，最短的纖維由53 個數(shù)據(jù)點構成，整個數(shù)據(jù)集共計151 986 個數(shù)據(jù)點。

4.3 實驗結果和分析

筆者方法的參數(shù)設置：采樣率r1=4%，距離閾值r2=0.8%，纖維中心閾值t=0.1，容器容量閾值R=50，容忍誤差δ=0.01，探測閾值λ=0.25。圖4(a～c)中：左側是未使用聚類算法的纖維結構，右側是使用連續(xù)框架聚類算法得到的結果，當完成纖維聚類后給不同的纖維集群分配不同的灰度值。通過實驗證明連續(xù)框架聚類算法能有效地分類纖維，結構相似的纖維被聚類為纖維集群并用同一種灰度值顯示，各纖維集群之間通過不同的灰度值清晰的區(qū)分。

圖4 纖維聚類結果Fig.4 The fiber clustering results

實驗證明連續(xù)框架聚類算法所需時間明顯低于一般的DBSCAN聚類方法，如表1所示。連續(xù)框架聚類算法的計算量與纖維數(shù)量以及纖維結構復雜度都有關，纖維結構越復雜觸發(fā)模型更新的次數(shù)就越多額外計算量也就越大，在不計模型更新整體計算量與纖維數(shù)呈線性關系。常用的聚類方法的整體計算量則與纖維數(shù)的平方呈線性關系。設纖維數(shù)量為n，構成纖維的體素點數(shù)量為m，聚類數(shù)量為l，則DBSCAN聚類方法的時間復雜度為O(n2m2)，連續(xù)框架聚類算法的時間復雜度為O(ln2m2)，因而當纖維數(shù)量越多時，連續(xù)框架聚類算法計算用時減少地越明顯。

表1 不同聚類方法計算效率對比

Table 1 The comparison of time efficiency using different clustering methods

纖維集連續(xù)聚類用時/sDBSCAN用時/sA1446.5720568.23B1452.5914023.67C1350.4810136.44

5 結 論

對快速密度峰值聚類算法進行改進，采用動態(tài)時間規(guī)整算法測量纖維之間的相似度，相比于常用的最近點平均距離和豪斯多夫距離兼顧了纖維整體形態(tài)上的相似性，能夠更準確地衡量相似性，提高聚類準確度。現(xiàn)有的聚類方法需要計算每對纖維之間的相似度，這將耗費大量的計算成本。相比于常規(guī)聚類方法，采用連續(xù)聚類框架進行腦纖維聚類在得到滿足可視分析需求結果的同時，顯著減少纖維相似度計算量，進一步降低整體聚類時間。