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

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

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

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

1 提取冠狀動(dòng)脈中心線

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

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

(1)

(2)

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

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

(3)

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

(4)

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

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

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

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

圖2 12方向示意圖

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

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

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

2 結(jié)果與評(píng)估

2.1 運(yùn)行環(huán)境與實(shí)驗(yàn)數(shù)據(jù) 算法在Matlab R2019a平臺(tái)下編譯、運(yùn)行。電腦配置為64位Windows 10操作系統(tǒng)，CPU主頻2.20 GHz，內(nèi)存16 GB。所用數(shù)據(jù)均來自鹿特丹冠狀動(dòng)脈算法評(píng)價(jià)框架(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 評(píng)價(jià)標(biāo)準(zhǔn) 用2個(gè)指標(biāo)評(píng)價(jià)中心線提取的準(zhǔn)確率[14]。第1個(gè)指標(biāo)為重疊率(rate of overlap, OV)，計(jì)算見公式5：

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

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

(5)

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

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

(6)

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

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

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

(7)

如圖5所示，TPRof存在于A段中，C段之前因?yàn)橐呀?jīng)發(fā)生錯(cuò)誤，故該段無TPRof。

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

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

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

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

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

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

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

3 討論

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

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