一種基于CTA 圖像的冠狀動脈自動分割算法*

羅宇杰，祝 磊*，馬 駿，2，薛凌云，張子恒，徐 平，劉亦安，嚴 明

(1.杭州電子科技大學自動化學院，浙江 杭州 310018；2.浙江省人民醫院，杭州醫學院附屬人民醫院，信息中心，浙江 杭州 310014)

根據世界衛生組織2005—2015 年的統計數據[1]，心血管疾病是全球的頭號死因，僅2016 年，全球估計有超過900 萬人死于心血管疾病。隨著成像分辨率的不斷提升以及技術本身無創性和便捷性的優點，CTA 逐漸成為臨床上早期篩查冠狀動脈血管疾病(Coronary Artery Disease，CAD)的有效方式[2]。為了輔助醫生提高診斷效率，基于CTA 的冠狀動脈自動分割成為了目前研究的一大熱點，它是分析血管異常(例如動脈瘤，狹窄和斑塊)的重要基礎步驟。

目前主要的血管分割算法包括:活動輪廓模型法[3]和區域生長法[4]。Guo 等人[3]提出一種結合血管形狀約束項的活動輪廓模型法，通過在基于梯度矢量信息構造的邊界能量函數中加入基于Hessian 矩陣的血管形狀約束項，充分融合圖像中的梯度場信息和灰度信息，克服了傳統活動輪廓模型法容易過分割的問題，但是該方法對于細小分叉的分割能力仍然不足，并且需要手動初始化左右冠脈分割輪廓。Wang Y 等人[4]提出了一種非參數形狀約束的活動輪廓模型法，通過在傳統的Chan-Vese活動輪廓模型法的能量函數中加入一個形狀約束項，提高了弱邊緣的分割能力，但是該方法比傳統的Chan-Vese 活動輪廓模型法更加耗時且對細小血管分割能力不足。Wang C 等人[6]結合隱式3D 血管模型和水平集法不斷調整血管中心線和血管估計直徑來實現更加準確的血管分割。然而，這些方法都存在邊界泄漏的問題，在處理血管分叉方面未能取得理想的結果[7]。Sammer 等人[8]通過對升主動脈區域的幾何特征變化分析，實現了對冠狀動脈起點的自動初始化，并且通過雙向分割有效減少了部分冠狀動脈分支的遺漏。Oksuz 等人[9]結合固定的單閾值預處理，Frangi 血管增強函數和傳統3D 區域生長法實現對冠脈的分割，但是此方法中的閾值預處理仍然保留了許多無用的區域并且傳統的3D 區域生長法需要人機交互且缺乏準確性。Tian 等人[9]采用三維多尺度血管濾波、閾值處理以及三維腐蝕等預處理，得到一個種子點集，然后使用基于統計的3D 區域生長法，通過邊緣控制系數改變區域生長法的相似閾值范圍，逐步得到完整的冠狀動脈樹，該方法中三維多尺度血管濾波非常費時且容易遺漏一部分冠狀動脈上的種子點導致血管分割不完整。Eiho等人[10]提出一種帶分叉檢測的區域生長法，并且在每一條分叉的分割中采用自適應生長準則以解決各條冠狀動脈分支灰度值不均的現象。吳明珠等人[11]提出了一種結合八元數矢量積的三維區域生長算法，對斷裂血管進行了有效的連接。程明等人[12]通過對已生長區域各分支根結點的回溯生長，提高了細小血管的分割能力。本文首先采用灰狼優化算法自動選取最優的多級閾值，去除CTA 圖像中大部分的非冠脈組織，無需再手動選取閾值，其次采用光流法識別冠狀動脈起始層，進一步提高了算法的自動化程度。本文提出的結合端點檢測的自適應區域生長法克服了傳統區域生長在分割血管時容易斷裂的情況，提高了冠脈的分割精度。

1 算法原理

針對上述方法仍然存在的局限性，提出了一種改進的面向CTA 圖像的冠狀動脈自動分割算法。本文算法框架如圖1 所示，主要包含圖像預處理、升主動脈分割以及冠狀動脈分割三大部分。第一部分圖像預處理主要包括窗寬窗位調整和基于灰狼優化算法的多級閾值處理，可有效抑制非冠脈組織對后續分割產生的干擾，防止過分割；第二部分升主動脈分割包括兩個階段:首先采用二維迭代區域生長法實現升主動脈的分割，然后利用光流法識別包含冠狀動脈血管的升主動脈層，識別左右冠狀動脈的起始層；第三部分采用結合端點檢測的自適應區域生長法實現左右冠狀動脈的分割和提取，相對于傳統區域生長法，本文算法提取的冠狀動脈分割結果更加完整且準確性更高。

圖1 算法框架圖

1.1 圖像預處理

在心臟的CTA 圖像中，不僅包含了冠狀動脈和升主動脈，還包括了心臟的其他組織如心房、心室以及其他非心臟組織如肺動脈、骨骼等。這些不同物理密度的組織在CTA 圖像上會具有不同的HU(Hounsfield Unit)值。對原始CTA 圖像進行窗寬(Window Width，WW)和窗位(Window Level，WL)調整的預處理可有效抑制肺動脈組織并提升冠狀動脈組織和背景的對比度。本文采用Dicom 頭文件中醫生所設置的WW 和WL，截取HU 值在范圍內的CTA圖像并歸一化到0～255 的灰度級，HU 值小于的像素點設置為0，HU 值大于的像素點設置為255。映射過程如式(1)所示:

式中，y表示歸一化后每個像素點的灰度值，x表示原始CTA 圖像中每個像素點的HU 值。

由于每個病人的個體差異以及醫生所設置的窗寬窗位不同等因素都會導致每個病人CTA 圖像序列中冠狀動脈的灰度值范圍不一，無法采用一個固定的閾值區分冠狀動脈和周圍灰度值相近的組織，所以本文采用基于最大Kapur 熵的多級閾值法[12]對CTA 圖像進行粗分割。在這種方法中，最優閾值選擇的標準是能否最大化基于灰度直方圖的Kapur熵函數之和，如式(2)所示:

式中，dim 是閾值級數，Hm是每一個閾值區間內的熵函數，如式(3)所示

式中，H0，H1，…，Hm為各級閾值的Kapur 熵，t1，t2，…，tm為閾值，p0，p1，…，pi為每一個灰度級的概率，ω0，ω1，…，ωm為每一個閾值區間的概率，L為圖像灰度級數。

隨著閾值級數的增加，多級閾值的計算復雜度呈指數增長，為此本文以式(2)Kapur 熵函數為目標函數J，采用灰狼優化算法[13]提高多級閾值法的全局收斂性和魯棒性。

1.2 升主動脈的分割

由于左右冠狀動脈的根節點位于升主動脈根部，本文算法首先采用二維區域生長法對升主動脈進行分割，然后采用光流法對冠狀動脈起始層進行識別。

1.2.1 迭代區域生長法分割升主動脈

CTA 圖像經過4 級閾值分割后，由于顯影劑的作用，升主動脈和冠狀動脈灰度值往往大于第4 級閾值，而心室、心房等非冠脈組織的灰度值分布往往低于第4 級閾值，所以為了減小非冠脈組織對后續分割的干擾，本文算法將4 級閾值處理后，灰度值低于第4 級閾值像素點都設為0，例如4 級閾值分別為24，92，151，204 時，則將灰度值小于204 的像素點設為0，將灰度值在204 到255 之間的像素點設置為1。然后采用霍夫變換檢測類圓形的升主動脈區域，以被檢測到的圓形中心作為區域生長種子點，利用二維區域生長法分割每一層升主動脈，最后得到完整的升主動脈。

1.2.2 光流法識別冠狀動脈起點

升主動脈分割完成后，結合光流法和心臟的醫學解剖結構先驗知識對左右冠狀動脈的起始位置進行識別。在CTA 圖像中，大部分的升主動脈區域都呈現類圓形，但是升主動脈根部與左右冠狀動脈相連的幾層CTA 圖像上，類圓形升主動脈區域會多一小段凸起的冠狀動脈區域。本文算法采用Lucas-Kanade 光流法[14]檢測相鄰兩層升主動脈間的光流位移距離，當某一升主動脈層包含冠狀動脈分支時，該層相對于前一層在冠狀動脈根結點處的光流位移距離會存在巨大的變化，通過對比相鄰兩層間光流位移距離的最大值即可判斷冠狀動脈的起始層。假設連續的兩層CTA 圖像強度為I(x，y，t)和I(x＋Δx，y＋Δy，t＋Δt)，根據光流法的第一個基本假設，連續的兩層圖像之間光流的強度不變，如式(4)所示:

將式(4)使用泰勒級數展開后可得:

式中，u和v就是代求的光流，但是只有式(6)一個約束方程無法求解2 個未知變量，所以Lucas-Kanade 光流法假設目標像素的周圍像素也存在相同的光流矢量，即在目標像素周圍m×m的窗口內，每個像素點都具有相似的光流矢量u和v。簡單的來說，就是通過增加約束方程來求解光流矢量u和v，如式(7)所示。

通過式(7)求解所得的u和v可以計算出每個像素點在兩層升主動脈間的位移距離Distances，如式(8)所示:

根據每一層Distances 最大的點，篩選出最大的前10 層做為候選冠脈起始層，然后結合三個心臟醫學解剖結構先驗知識進一步篩選出左右冠狀動脈的起始層，第一，升主動脈根部只存在左右兩個冠狀動脈開口，在10 層冠脈候選層中，位移距離最大的是左冠狀動脈起始層；第二，左右冠脈起始層間相隔10 層以上且左冠狀動脈起始層位于右冠狀動脈起始層的Z軸上方；第三，如圖2 所示左冠狀動脈起始點OL位于該層質心C的-90°～45°之間，右冠狀動脈起始點OR位于該層質心C的45°～180°之間。最后在滿足上述三個條件的剩余候選冠脈起始層中選擇層號最小的層做為右側冠狀動脈起始層。

圖2 左右冠狀動脈起點位置示意圖

1.3 冠狀動脈分割

1.3.1 血管多尺度特征提取

冠狀動脈的直徑通常介于2 mm～7 mm 之間，從升主動脈根部出發冠狀動脈會變得越來越細，在單尺度條件下只能對某一特定直徑的冠脈血管進行特征提取，所以本文引入多尺度分析。如式(9)所示，通過將預處理后的圖像Ithresh(x，y)與不同方差的高斯核Gσ(x，y)卷積得到一組不同尺度的圖像，每一幅圖像都代表不同維度的原始圖像，只是比例不同。通過選取不同的σ，可以增強不同大小的血管輪廓。

然后計算每個像素點在不同尺度空間下的Hessian 矩陣Hσ(x，y)，如式(10)所示:

Hessian 矩陣是實對稱矩陣，存在兩個實特征值λ1和λ2(｜λ1｜＜｜λ2｜)。特征值λ1和λ2的符號以及絕對值的大小對應了圖像中不同形狀的結構。在清晰的血管邊緣上，灰度值存在躍變，此時Hessian矩陣的特征值滿足｜λ1｜≈0 且｜λ1｜≤｜λ2｜。最后如式(11)所示，Frangi 血管相似性函數V(σ，x，y)利用不同尺度σ的Hessian 矩陣特征值λ1和λ2對血管狀區域進行增強。V(σ，x，y)的取值范圍在0 到1之間，越接近于1 表示該點處的組織與血管結構的相似性最大。

式中，β、γ為常數

每一個像素的V(x，y)為不同尺度響應中的最大值，如式(12)所示:

式中，σmin取0.5，σmax取2.5，不同σ間的步長為1。

1.3.2 冠狀動脈分割

本文算法從左右冠狀動脈起點開始，在上一節Frangi 血管增強后的CTA 圖像上，使用結合端點檢測的自適應區域生長法(Branch based Adaptive Region Growing，BARG)分割每一層CTA 圖像中的血管區域。冠狀動脈在三維空間中的連通性使得上下兩層血管在二維平面上存在重疊區域，本文算法在前一層血管區域內以一定步長選取15 個種子點作為當前層區域生長的種子點。由于三維空間中血管的拓撲結構復雜，上下層血管之間只有部分重疊區域，所以從前一層血管區域中選取的15 個種子點只有一部分會落在當前層的血管區域內。如果種子點落在血管區域內，區域生長法從該種子點開始，不斷將區域邊緣像素8 鄰域中和種子點灰度值相差在1%以內的像素點歸并到已分割的區域中，直到不能生長為止。為了減少種子點的重復分割，每一個種子點生長完成后，都會剔除那些已經包含在分割完成區域內的種子點；如果種子點落在血管區域外，區域生長法會出現嚴重的過分割現象，所以當超過2 000 個像素點被分割成血管區域時，算法會強制終止生長，并且將該種子點判定為錯誤種子點，所分割的區域也會被丟棄。最后將15 個種子點分割得到的區域之和做為該層的最終分割結果，繼續執行下一層分割。

在完成第一輪從左右冠狀動脈起點開始，沿著Z軸自頂向下的迭代生長后，可以獲得大部分的冠狀動脈樹，但是由于冠狀動脈的三維拓撲結構復雜，會遺漏部分主支末端在三維空間中上下起伏延伸的細小分支。為此提出BARG 算法可利用血管骨架線上體素的連通性特點，尋找漏分割冠脈支的起點，從各起點處繼續分割剩余冠狀動脈。首先使用細化方法[14]提取血管的骨架線，然后利用血管骨架線的26 鄰域單連通性，將骨架線上的體素分為如圖3 所示的三類:①端點P1，在P1的26 鄰域中，少于2 個鄰域點位于骨架線上；②普通連接點P2，在P2的26 鄰域中，有且僅有2 個鄰域點位于骨架線上；③分叉點P3，在該P3的26 鄰域中，有多于2 個鄰域點位于骨架線。根據以上骨架線的連通特性分析，從各個端點P1和分叉點P3開始，使用區域生長法嘗試雙向生長，將圖3中點P1和P3之間的血管分割完整。在每一輪迭代完成后，本文算法都會重新檢測端點P1和分叉點P3，直到沒有新增的端點可以繼續生長而結束分割。

圖3 骨架線體素分類

2 實驗結果與分析

本文采用20 例CTA 數據進行算法測試，數據來源于山東大學齊魯醫院，均帶有醫生標注，由14 例男性和6 例女性病例組成，年齡在51 歲～85 歲之間，平均年齡(65±9.2)歲，每例數據包括大約150 層～250層包含冠脈血管的軸向切片，每層軸向切片的圖像分辨率為512 像素×512 像素，像素間距dx，dy，dz分別為0.35 mm，0.35 mm，0.80 mm。測試平臺采用64 位MacOS 10.14.5 系統，Intel Core i9-9980HK CPU，16G內 存，MATLAB R2018a 軟件。本文使用Dice、Jaccard、MSD (Mean Squared Surface Distance)、MAXSD(Maximum Surface Distance)作為分割評價指標[16]。Dice 和Jaccard 相似性系數用于衡量算法分割區域和專家手工標注區域之間的重合度，其定義如式(13)和式(14)所示:

式中，GT 是專家人工標注的分割區域，Seg 是算法自動分割區域。Dice 和Jaccard 取值都在0 到1 之間，取值越接近1，表明算法自動分割和專家手工分割區域之間的重合度越高，分割精度越好。MSD 和MAXSD 是衡量算法分割輪廓與專家手工標注輪廓間距離差異的評價指標，其定義如式(15)和式(16)所示:

式中，Bseg和BGT分別代表算法分割輪廓和人工標注輪廓點集，d(x，BGT)代表對于Bseg上的任意一點x到BGT的最短距離，d(y，Bseg)代表對于BGT上的任意一點y到Bseg的最短距離。MSD 和MAXSD 越小說明算法自動分割輪廓越接近于手工標注輪廓形狀，表明分割效果越好。

2.1 預處理結果

本文的預處理主要包括窗寬窗位的調整和基于灰狼優化的最大熵多級閾值處理。圖4(a)為原始的16 位Dicom 格式的CTA 數據，其HU 范圍在-1 024到3 071 之間，冠狀動脈和周圍背景組織在CTA 原始數據上對比度較差。窗寬窗位調整后的結果如圖4(b)所示，升主動脈和冠狀動脈被明顯增強的同時，心臟外部的背景組織如肺部血管得到了明顯的抑制。為了進一步區分冠狀動脈組織和周圍其他的心臟組織，本文采用基于灰狼優化的最大熵4 級閾值處理，將圖像分為5 個不同的灰度區間，結果如圖4(c)所示。最后本文算法只保留最高一級區間內的像素點，以排除非冠狀動脈組織對后續升主動脈分割和冠狀動脈分割的干擾。預處理最后結果如圖4(d)所示，該例CTA 數據對應的4 級閾值分別為24，92，151，204 時，本文算法只保留[204-255]區間內的像素點用于下一步升主動脈分割。

圖4 預處理結果

2.2 升主動脈分割

由于左右冠狀動脈的起點位于升主動脈的根部，為此準確的升主動脈分割是識別冠狀動脈起點的重要步驟。如圖5(a)所示，本文算法首先利用霍夫變換檢測類圓形的單層升主動脈區域，并以圓心做為該層區域生長法的種子點，分割得到每一層的升主動脈。當遇到如圖5(b)和5(c)包含左右冠狀動脈起始分叉的不規則升主動脈層時，霍夫變換無法準確檢測到圓心，可選擇上一層的種子點做為該層的起始種子點。圖5(d)所示為所有升主動脈層三維重構圖。

圖5 升主動脈分割結果

升主動脈分割完成后，本文算法結合光流法與心臟醫學解剖結構先驗知識檢測左右冠狀動脈起始層。首先通過光流法計算每一層所有像素相對于前一層的光流位移距離，圖6 所示為每一層中最大的光流位移距離。在該例CTA 數據的90 層升主動脈中，圖5(a)所示的普通類圓形升主動脈層，不會出現位移激增的現象，對應于圖6 曲線中較平緩的區域。但是當遇到圖5(b)和5(c)這些包含冠狀動脈分叉的不規則升主動脈層時，會出現位移距離激增的現象。一般而言，選擇位移距離最大的層(圖6第58 層)做為左冠狀動脈起始層。由于左右冠狀動脈相距10 層以上即在Z軸方向上大于0.8 cm，所以圖6 第62 層和第64 層并不能做為右冠狀動脈的起始層，轉而選擇位移距離不如第62 層大的第81層做為右側冠狀動脈的起始層。最后需要驗證左冠狀動脈起始層中位移最大的像素點是否處于該層質心的-90°～45°之間，右冠狀動脈起始層的最大位移像素點是否處于該層質心的45°～180°之間，在一些極端的病例中，該先驗知識可以有效降低識別的錯誤率。本文在20 例CTA 數據上對該算法進行了測試，結果如表1 所示，左側冠狀動脈起始層檢測成功19 例，準確率達到95%，右側冠狀動脈起始層檢測成功17 例，準確率達到85%。

圖6 冠狀動脈起始層識別

表1 光流法識別冠狀動脈起始層

2.3 冠狀動脈分割

為了增強不同管徑的冠狀動脈，本文算法采用了多尺度融合的Frangi 管狀結構增強濾波。圖7(a)、7(b)、7(c)分別代表δ＝0.5、δ＝1.5、δ＝2.5的管狀結構增強結果，不同的δ可增強不同直徑的血管，圖7(d)為多尺度融合后的結果，冠狀動脈得到了有效增強。

圖7 多尺度Frangi 血管增強

圖8 分別為傳統區域生長法、不帶端點檢測的區域生長法、本文算法和人工標注的三維分割結果。

圖8 分割結果3D 重構圖

如圖8(a)所示，如果不對原始CTA 數據進行預處理，傳統的區域生長法非常容易產生過分割現象。圖8(b)和(c)所示為本文算法使用端點檢測進行迭代生長前后的分割對比結果，箭頭①和②所指的分支末端具有較為復雜的三維拓撲結構，在Z軸方向上存在上下起伏，圖8(c)的分割結果驗證了本文算法對拓撲結構復雜的分叉有更好的分割能力。對于圖8 中箭頭③和④所指的分支，其末端由于顯影劑灌注等原因，灰度值明顯低于上層血管，本文算法能夠自適應選擇多級閾值中最合適的灰度區間進行分割，例如2.1 節預處理后4 級閾值分別為24，92，151，204，如果[204，255]灰度區間無法繼續分割，算法自動選擇[151，255]灰度區間繼續嘗試分割，直到分割完整。

20 例CTA 數據的分割結果如表2 和圖9 所示，采用本文提出的BARG 算法得到的分割結果，Dice和Jaccard 系數總體分布高于不帶端點檢測的區域生長法，分別平均提高了4%和8%，達到了0.70 和0.57，MSD 和MAXSD 總體分布低于不帶端點檢測的區域生長法，分別平均降低了2%和4%，達到了0.40 和2.66。綜上實驗對比結果，本文提出的BARG 算法對灰度不均且拓撲結構復雜的細小分叉有更好的分割能力以及更高的分割精度。

表2 不同方法結果比較

圖9 指標分布盒型對比圖

3 總結與討論

本文提出了一種面向CTA 數據的冠狀動脈自動分割算法。首先，通過圖像預處理有效抑制了非冠脈組織，提升了冠狀動脈和背景的對比度；其次，結合光流法與心臟解剖結構先驗知識檢測帶冠狀動脈分叉的不規則升主動脈層，避免了手動初始化左右冠狀動脈的起始點；最后，相比于傳統區域生長法，本文提出的結合端點檢測的自適應區域生長法對于灰度不均且拓撲結構復雜的細小分支有更好的分割能力和準確性，但是和人工標注的冠狀動脈仍有差距。同時，本文算法仍存在諸多可以改進之處:CTA圖像中細小分支的血管和狹窄部位的血管由于顯影劑計量低，血管邊界處灰度不均，采用區域生長法分割血管時容易產生過分割和欠分割。所以，在后續的研究過程中擬將血管分割與深度學習方法相結合，神經網絡可以根據人工標注過的CTA 圖像，自動學習血管的多種特征，如紋理、邊界、灰度值等，更好地判斷血管的邊界，使得灰度值較弱、邊界模糊的細小血管也能夠被很好的分割。