張子恒，祝 磊，馬 駿，徐 平，劉亦安，薛凌云

(杭州電子科技大學(xué)自動化學(xué)院，浙江 杭州 310018)

冠心病(coronary artery disease, CAD)在中國心血管疾病死因中居第2位[1]。診斷CAD主要依靠CT血管造影術(shù)(computed tomography angiography, CTA)和冠狀動脈造影(coronary artery angiography, CAG)。CTA易于三維可視化，掃描時間短、空間分辨率較高，已越來越多用于診斷CAD。CTA診斷CAD時，需要對圖像的血管結(jié)構(gòu)進(jìn)行包括分割血管區(qū)域，提取中心線、分叉點以及測量血管直徑等精確處理，其中提取血管中心線最為重要[2]。

現(xiàn)有中心線提取方法分為全自動、半自動和人工提取3類。人工提取中心線簡單、可靠，但耗時、費力[3]。全自動方法一般指提取中心線過程中不施加任何人工干預(yù)，較經(jīng)典者如利用圓柱體模型擬合提取冠狀動脈血管后進(jìn)一步提取中心線[4]；半自動方法指在血管提取過程中需人為指定1個或數(shù)個點作為參照點，如分割主動脈和冠狀動脈并行三維重建后利用血管局部灰度值及走向選取迭代起始點[5]。此外，還有利用DBSCAN聚類算法進(jìn)行分割后利用冠狀動脈主曲線提取中心線[6]，利用改進(jìn)型Frangi濾波器消除心腔邊界處無關(guān)階躍響應(yīng)后自動提取中心線[7]，利用快速行進(jìn)算法提取圖像中初始點之間的最短路徑，據(jù)此在原始圖像上計算冠狀動脈中心線[8]，以及同時使用基于邊界的距離變換和基于源點(種子點)的距離變換[9]等。隨著人工智能的興起，可利用神經(jīng)網(wǎng)絡(luò)提取中心線，如手動設(shè)置種子點，利用卷積神經(jīng)網(wǎng)絡(luò)(convolutional neural network, CNN)結(jié)合冠狀動脈半徑沿血管不同方向預(yù)測血管中心線走向而提取中心線[10]。針對傳統(tǒng)骨骼細(xì)化方法在血管連接處易產(chǎn)生細(xì)小分支的問題，本研究結(jié)合Dijkstra算法提出一種基于骨骼細(xì)化的全自動中心線提取算法，旨在使所提取的骨架線更加準(zhǔn)確、平滑。

1 提取冠狀動脈中心線

1.1 圖像預(yù)處理 原圖中分割出冠狀動脈進(jìn)行三維重構(gòu)，之后去除噪聲和增強血管區(qū)域。

各向異性擴散濾波器(anisotropic diffusion filter)[11]系空域濾波器，平滑程度取決于擴散函數(shù)，在圖像各個方向上有所不同。各向異性擴散濾波見公式1：

(1)

(2)

其中，k為與噪聲級和邊界強度相關(guān)的常數(shù)，這里選取k=max|u|。

若三維圖像I中存在一點A，想要判斷點A是否為血管上一點就需要分析圖像I在A點的局部特征，利用Hessian矩陣就可以判斷該點的位置。點A的Hessian矩陣H(A)見公式3：

(3)

式中Ixx(A)，Ixy(A)等分別代表圖像I在點A處的二階偏導(dǎo)，并且H(A)三個特征值符合|λ1|≤|λ2|≤|λ3|的關(guān)系。若點A屬于血管上一點利用基于Hessian矩陣的Frangi血管增強函數(shù)[12]就可將血管區(qū)域增強，該血管增強函數(shù)見公式4：

(4)

1.2 提取中心線 三維細(xì)化算法屬三維形態(tài)學(xué)操作，可用“在草地上放火”加以解釋。將圖像視作一塊草地，在其兩端各位置同時點火，草地上會形成兩面火墻，逐漸向草地中央逼近，來自草地兩端火墻相遇位置的連線就是該草地的骨架線，各點火點即為算法的起始點，被燒掉區(qū)域即為算法中可刪除的點。

三維空間中存在一點p，根據(jù)該點到其相鄰一點的歐氏距離的長度，可將其關(guān)系分為6-鄰接、18-鄰接和26-鄰接，見圖1。用Nj(p)(j=6,18,26)表示點p和j鄰接點的集合即鄰域(neighbor)。圖1中用U、D、W、E、S、N和點p表示N6(p)，6-鄰接中的6個方向即為三維立體網(wǎng)格中的6個主要方向；用N6(p)和12個圓圈標(biāo)識表示N18(p)；用N18(p)和6個五角星標(biāo)識表示N26(p)。三維細(xì)化算法通常在二值化圖片上進(jìn)行操作，而二值化圖像只有0，1兩種值，將圖中值為1的點定義為目標(biāo)點(black point)，值為0的點定義為背景點(white point)。如果集合N26(p)/{p}恰好包含一個目標(biāo)點，則目標(biāo)點p是圖片中的曲線端點(curve end-point)；若N6(p)包含至少1對相對背景點，則目標(biāo)點p是圖片中的表面端點(surface end-point)。如果圖像中任意一目標(biāo)點至少與1個背景點6-鄰接，則稱為邊界點(border point)。簡化點(simple point)則是可刪除(將該點的值由1置0)而不改變圖形拓?fù)浣Y(jié)構(gòu)的目標(biāo)點，應(yīng)只由某些類型邊界點構(gòu)成；曲線或曲面的端點與圖形的拓?fù)浣Y(jié)構(gòu)相關(guān)，因而不能被刪除。任何細(xì)化算法中，可刪除的點只有簡化點，而傳統(tǒng)三維細(xì)化算法大多是從體素的6個主要方向并行操作，可能會將一些曲線或曲面的端點當(dāng)作簡化點刪除，導(dǎo)致原圖像的拓?fù)浣Y(jié)構(gòu)發(fā)生改變。為保證細(xì)化算法不會誤刪某些端點，以保證原圖形的拓?fù)浣Y(jié)構(gòu)，本研究改進(jìn)傳統(tǒng)算法，使每個主迭代步驟由多個并行子迭代過程組成，其中每個子迭代過程只能刪除某種特定的邊界點。

圖1 6-鄰接、18-鄰接和26-鄰接示意圖

采用12方向?qū)θS重構(gòu)后的血管進(jìn)行細(xì)化處理，每個迭代步驟根據(jù)圖2所示12個方向進(jìn)行并行細(xì)化。選擇刪除方向的順序為(US，NE，WD)、(ES，UW，ND)、(SW，UN，ED)、(NW，UE，SD)，據(jù)此又可將這12個方向分為4組，每組3個方向，每個小組包含6個主要方向，即在第1個子迭代過程US方向上可刪除某些U和S方向上的邊界點，其余各迭代過程以此類推。需要被細(xì)化的冠狀動脈會在每個方向上均勻收縮，即從外層開始，向內(nèi)被層層腐蝕。當(dāng)12個子迭代過程中不再有可被刪除的目標(biāo)點時，主迭代過程結(jié)束，即中心線提取完成。

圖2 12方向示意圖

迭代過程中，利用匹配模板可判斷某點是否為可被刪除的簡單點。若圖像上一點p的3×3×3鄰域與一組匹配模板中的1個模板相符，則認(rèn)為點p是可在細(xì)化過程中刪除的簡單點。以圖2中US方向為例，展示該方向上的14個基本模板(圖3)，可用這些模板經(jīng)過一系列旋轉(zhuǎn)或反轉(zhuǎn)操作獲得其他方向的模板。

1.3 剔除骨架線分支 血管為不規(guī)則管狀結(jié)構(gòu)，進(jìn)行細(xì)化算法時會產(chǎn)生細(xì)小分支，影響后續(xù)測量冠狀動脈直徑及判定是否產(chǎn)生狹窄的準(zhǔn)確性。為消除干擾，提出基于Dijkstra算法的分支消除算法。

Dijkstra算法是一種解決有向加權(quán)圖中單源最短路問題的算法，特點是迭代從起點向外層層擴張，直至到達(dá)終點[13]。首先將上一步提取出的三維骨架線轉(zhuǎn)換為由節(jié)點和有向邊組成的有向圖，圖中每個節(jié)點代表被提取出的體素(voxel)，不同體素間歐氏的距離由有向邊表示。根據(jù)Dijkstra算法，將所有體素設(shè)為全集U；用點集G表示已找到最短路徑的點，初始狀態(tài)下僅有起點S∈G；設(shè)點集O表示還未找到最短路徑的點的集合，則可以推出O=U-G；設(shè)Lk為當(dāng)前情況下，起點經(jīng)過G中若干點到點k(k∈U)的最短距離，初始狀態(tài)下LS=0，其他均為+∞。當(dāng)算法開始時，從起點S開始，沿某條有向邊(設(shè)權(quán)值為arcs)找到起點的一個鄰居n；此時令Ln=min{Ln,LS+arcs}；按此方式更新起點所有鄰居；在集合O中找到Lk最小的點v，則Lv即起點到v的最短路徑長度；將點v從O中取出加入G，對點v重復(fù)上述所有操作；如此重復(fù)，直到G=U，即O=φ時，算法結(jié)束，Lk即為從起點到各點的最短路徑長度。冠狀動脈中心線的有向圖模型見圖4。最后再將有向圖轉(zhuǎn)換為三維骨架線，即可得到去除細(xì)小分支的冠狀動脈中心線。

2 結(jié)果與評估

2.1 運行環(huán)境與實驗數(shù)據(jù) 算法在Matlab R2019a平臺下編譯、運行。電腦配置為64位Windows 10操作系統(tǒng)，CPU主頻2.20 GHz，內(nèi)存16 GB。所用數(shù)據(jù)均來自鹿特丹冠狀動脈算法評價框架(Rotterdam Coronary Artery Algorithm Evaluation Framework, RCAAEF)提供的公共數(shù)據(jù)集，包含32例，于2005年6月—2006年6月采自荷蘭鹿特丹伊拉斯謨醫(yī)學(xué)中心(Erasmus MC)，其中20例以Siemens 64排螺旋CT、12例以Siemens雙源CT采集。

2.2 評價標(biāo)準(zhǔn) 用2個指標(biāo)評價中心線提取的準(zhǔn)確率[14]。第1個指標(biāo)為重疊率(rate of overlap, OV)，計算見公式5：

圖3 US方向上的14個基本模板 點p為圖中各模板中心點，圖中符號“○”表示背景點，“●”表示目標(biāo)點，“·”表示該位置上可為目標(biāo)點也可為背景點；字母X表示在這些位置上至少有1個為目標(biāo)點，V和W表示這些位置上至少有1個為背景點；如2個不同位置標(biāo)注字母Z，則表示這2個位置上必須為不同的點(如1個為目標(biāo)點，另1個為背景點)，但其位置不固定

圖4 冠狀動脈有向圖模型 虛線部分為提取的骨架線模型，點S為中心線的起點，點E為終點，A為骨架細(xì)化后產(chǎn)生的細(xì)小分支的端點。虛線SE表示實際的中心線走向，和表示三點的有向邊其值等于兩點之間的歐氏距離，所以圖4中的最短路徑即為間的體素點得以保留，A點則被刪除

(5)

式中‖·‖表示為這一類點的集合，TPMov表示本文算法所求出的中心線上的一點，該點到實際中心線的位置小于參考標(biāo)準(zhǔn)中給出的誤差半徑，與之對應(yīng)的實際中心線上的點記為TPRov；若求出的點位于誤差半徑之外，則將其記為FPov，與之對應(yīng)的實際中心線上的點記為FNov，這4種點的位置關(guān)系見圖5。

第2個指標(biāo)為平均距離(average distance， AD)，即所有由本文算法提取出的中心線的位置與實際中心線之間的距離，用distance(i)表示第i個點到該點的實際位置的距離，則AD見公式6：

(6)

圖5 評價中所用術(shù)語的位置說明 血管分為A、B、C、D四段，代表中心線實際位置與算法提取的中心線之間的不同結(jié)果，位于圖片頂端的字符表示該段中算法提取的中心線類型，中心線下方的字符代表實際中心線類型。TPMov、TPRov、FNov均只存在一種情形，而FPov有可能存在于未出現(xiàn)冠狀動脈的區(qū)域(D段)

圖6 圖片預(yù)處理過程 A.原始CTA圖像； B.各向異性擴散濾波處理； C.血管區(qū)域； D.三維重構(gòu)血管

冠狀動脈為直徑逐漸縮小的管狀結(jié)構(gòu)，近心段冠狀動脈直徑較大，而其分支血管遠(yuǎn)端直徑較小，臨床意義較小，故可以血管中首次出錯前的重疊率(overlap until first error, OF)來評估程序算法的穩(wěn)定性。OF見公式7：

(7)

如圖5所示，TPRof存在于A段中，C段之前因為已經(jīng)發(fā)生錯誤，故該段無TPRof。

2.3 實驗結(jié)果 本研究冠狀動脈CTA圖像為512×512，每例包含近300幅圖像，根據(jù)不同長度，每條血管會在40～130幅圖像中連續(xù)出現(xiàn)，圖6A為原始圖像。利用各向異性擴散濾波器對原始圖像進(jìn)行去噪處理(圖6B)，可明顯抑制噪聲；根據(jù)數(shù)據(jù)集中提供的每條冠狀動脈的中心點位置，將無關(guān)區(qū)域刪除(像素值置0，不改變圖像大小)，保留原圖信息的區(qū)域包含該血管在此層切面的完整圖像；再利用Frangi血管增強函數(shù)，增強被保留圖像中的血管區(qū)域，使其與背景形成較明顯的差異，為區(qū)域生長法提供更易區(qū)分的灰度值；最后在該區(qū)域中選取血管區(qū)域的像素作為種子點進(jìn)行區(qū)域生長，獲得完整血管區(qū)域(圖6C)。對每例所有存在該條血管的層切面進(jìn)行上述處理，可在二維圖像中分割出整條血管；再根據(jù)原始圖像中像素與世界系坐標(biāo)的對應(yīng)關(guān)系進(jìn)行三維重建，即得到該血管的三維數(shù)據(jù)(圖6D)。

完成圖片預(yù)處理后，利用改良骨骼細(xì)化算法提取血管中心線，刪除細(xì)小分支(圖7)。

為保證算法具有較強的魯棒性，分別選取6個不同病例的血管進(jìn)行評價(表1、2)。

表1 未改進(jìn)算法提取中心線部分結(jié)果

圖7 提取中心線結(jié)果 A.右冠狀動脈； B.左前降支淺 (灰色部分為重構(gòu)的冠狀動脈，中間藍(lán)色曲線為提取中心線)

表2 改進(jìn)算法提取冠狀動脈中心線部分結(jié)果

相比未移除細(xì)小分支的骨骼細(xì)化算法，本算法將OV提高了2%，AD減少了38.2%，OF提高了16%；中心線提取的平均時間為0.48 s。計算數(shù)據(jù)集中所有血管的OV、AD、OF，平均重疊率、平均距離和平均第1次出錯前重疊率分別為0.985、1.974和0.832。

3 討論

實驗結(jié)果表明本算法用于提取冠狀動脈中心線有效，相比文獻(xiàn)[4-5]提出的重疊率分別為0.847和0.670的2種方法，本算法精度較高。提取冠狀動脈中心線在計算機輔助診斷CAD中具有重要作用，提取速度直接影響運行效率。本算法提取中心線用時較短，可滿足臨床需求。

本算法仍存在需改進(jìn)之處：未對提取完整冠狀動脈血管樹中心線性能進(jìn)行測試，后續(xù)將擴大數(shù)據(jù)樣本進(jìn)一步完善；以區(qū)域生長法分割血管區(qū)域，要求較高，且較依賴圖像閾值，即圖像中因加入對比劑而產(chǎn)生的高亮區(qū)域，若遠(yuǎn)端血管無強化，則分割將會產(chǎn)生較大誤差，從而影響提取冠狀動脈中心線。后續(xù)研究中擬與人工智能相結(jié)合，利用訓(xùn)練完成的神經(jīng)網(wǎng)絡(luò)自動分割血管區(qū)域。