改進型骨骼細化算法提取冠狀動脈中心線

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

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

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

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

1 提取冠狀動脈中心線

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

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

(1)

(2)

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

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

(3)

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

(4)

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

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

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

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

圖2 12方向示意圖

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

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

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

2 結果與評估

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

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

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

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

(5)

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

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

(6)

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

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

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

(7)

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

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

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

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

表1 未改進算法提取中心線部分結果

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

表2 改進算法提取冠狀動脈中心線部分結果

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

3 討論

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

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