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

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

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

1 纖維相似度計(jì)算方法

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

1.1 動(dòng)態(tài)時(shí)間規(guī)整

DTW主要應(yīng)用于語音識(shí)別，用來衡量兩個(gè)語音信號(hào)的相似度，其原理是要找到兩條時(shí)間序列之間的最短折疊路徑[16]。將時(shí)間序列沿著時(shí)間軸進(jìn)行拉伸和折疊等操作使其扭曲在一起，設(shè)現(xiàn)有時(shí)間序列X和Y，它們的長度分別為m和n，即

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

(1)

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

(2)

設(shè)K為折疊路徑W的長度，則折疊路徑W為

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

(3)

那么

max(m,n)

(4)

式中：W的元素為時(shí)間序列X，Y中時(shí)間節(jié)點(diǎn)的連接情況，即wk(i,j)。折疊路徑W遵循的3 個(gè)約束條件為

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

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

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

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

(5)

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

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

(6)

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

最佳折疊路徑可以通過動(dòng)態(tài)算法[17]得到：首先計(jì)算一個(gè)維數(shù)為m×n的成本矩陣D，成本矩陣中的每個(gè)元素D(i,j)由距離d(i,j)與相鄰的元素的和遞歸地賦值，其計(jì)算式為

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

(7)

當(dāng)建立成本矩陣后，可以從終點(diǎn)D(m,n)反向地尋找一條累積距離最短的路徑。為了達(dá)到這個(gè)目的，需要對成本矩陣沿著左、下和對角線3 個(gè)方向進(jìn)行貪婪搜索。這3 個(gè)相鄰元素中的最小元素將被加入折疊路徑，并且搜索過程將從該元素繼續(xù)進(jìn)行直至搜索到D(1,1)。如果折疊路徑經(jīng)過成本矩陣中的元素D(i,j)，則表示時(shí)間序列X上的第i個(gè)點(diǎn)與時(shí)間序列Y上第j個(gè)點(diǎn)相連。由于一條時(shí)間序列單個(gè)點(diǎn)可能對應(yīng)著另一條時(shí)間序列上的多個(gè)點(diǎn)，因此，該方法可以衡量不同長度時(shí)間序列之間的相似度。

1.2 DTW測量纖維相似度

(8)

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

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

(9)

2 快速密度峰值搜索算法

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

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

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

(10)

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

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

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

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

(11)

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

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

(12)

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

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

3 連續(xù)聚類框架

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

3.1 初 始 化

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

γi=ρj×δi

(13)

當(dāng)γi>t×γmax時(shí)，將纖維選擇為纖維中心，實(shí)驗(yàn)證明閾值t取0.1時(shí)，能達(dá)到較好的聚類效果。

3.2 標(biāo)記新纖維路徑

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

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

3.4 相關(guān)性信息

(14)

3.5 Page-Hinkley測試

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

(15)

(16)

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

3.6 模型更新

4 實(shí)驗(yàn)結(jié)果

4.1 實(shí)驗(yàn)環(huán)境

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

4.2 實(shí)驗(yàn)數(shù)據(jù)

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

4.3 實(shí)驗(yàn)結(jié)果和分析

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

圖4 纖維聚類結(jié)果Fig.4 The fiber clustering results

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

表1 不同聚類方法計(jì)算效率對比

Table 1 The comparison of time efficiency using different clustering methods

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

5 結(jié) 論

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