冠狀動脈造影（CAG）作為一種經典的血管影像模態，在臨床上廣泛應用于心血管疾病的診斷。對于基于CAG的計算機輔助診斷系統而言，準確的血管分割至關重要。但是，對安裝心臟起搏器患者進行CAG檢查時，所獲得的圖像中往往存在起搏器的干擾。在這種情況下，血管分割將變得更為困難。為了去除起搏器的干擾以優化血管分割的結果，本文提出了一種基于偽CAG圖像生成并與待分割CAG圖像配準后移除起搏器干擾的方法。首先，通過CAG視頻中的部分圖像序列生成僅包含起搏器的偽CAG圖像；然后，利用一個圖像局部特征描述子配準偽CAG圖像和待分割CAG圖像中起搏器的相對位置；最后，綜合配準結果和起搏器、主要血管的分割結果，移除重疊起搏器的干擾，優化主要血管的分割。在臨床采集的帶起搏器干擾的11幅CAG圖像上評估所提出的方法，戴斯系數的平均優化率達到12.03%，這表明本文方法可有效移除重疊起搏器，從而達到優化主要血管分割的目的。綜上，本文研究結果或可助力基于CAG的計算機輔助診斷系統，進一步提高其診斷準確性和效率。
引用本文： 黃藝, 楊虹波, 夏夢華, 曲亞男, 郭翌, 周國輝, 張峰, 汪源源. 移除冠脈造影圖像中重疊起搏器以優化主要血管的分割. 生物醫學工程學雜志, 2022, 39(5): 853-861. doi: 10.7507/1001-5515.202104023 復制
引言
近年來，全球心血管疾病患者數持續增加，中國雖然在心血管疾病防治工作上取得了一些成績，但據估計，中國現有心血管疾病患者數仍高達3.3億。其中，由冠狀動脈發生粥樣硬化、血管腔狹窄，進而造成心肌缺血、缺氧乃至壞死的一類心臟病，稱為冠狀動脈粥樣硬化性心臟病（后文簡稱：冠心病）。冠心病作為一種典型的心臟病，其患者數達到1100萬[1-3]。冠狀動脈造影（coronaryangiography，CAG）將導管通過外周動脈血管送至造影位點，然后注射少量造影劑再利用X線對心臟血管進行顯像，其作為冠心病診斷的有效手段，在臨床上得到了廣泛的應用[4]。自動分割CAG圖像中的冠狀動脈對于冠心病的計算機輔助診斷（computer-aided diagnosis，CAD）而言是一個基礎，在此基礎上進行狹窄評估[5]、三維重建[6]和動態評估[7]，可為醫生診斷提供參考。此外，CAG還可以與其他檢查手段相結合，如血管內超聲（intravascularultrasound，IVUS）[8]、血管光學相干斷層掃描成像（opticalcoherence tomography，OCT）[9]，從而對血管狀況進行更全面的評估。
如圖1所示，CAG圖像中前景（血管）和背景區域之間的對比度較低，并且造影劑會隨著血液流動而產生分布不均的現象，使得血管區域的灰度分布也不盡相同，這使得分割CAG圖像中的主要血管存在一定困難。為了有效地分割CAG圖像中的血管，研究者們做了許多的探索，形成較為常用且有效的“兩步走”策略，即先對CAG圖像進行優化處理，主要目的是提高血管與背景區域之間的對比度，使得兩者之間更容易區分，然后再對優化后的CAG圖像進行血管分割。在“兩步走”策略的基礎上，涌現出許多優秀的方法。Cruz-Aceves等[10]提出了一種兩階段方法，首先使用多尺度伽柏（Gabor）濾波器和頂帽（Top-hat）運算優化處理圖像，而后使用基于多目標優化的新閾值方法對優化結果進行二值化，實現對冠狀動脈的分割。這一方法的優點在于多尺度Gabor濾波器能夠有效提升灰度不均勻CAG圖像的對比度，且閾值方法中同時引入了三種加權優化準則，大大提高了分割的魯棒性。Cervantes-Sanchez等[11]同樣使用了Gabor濾波器，并結合差分演化的優化策略選擇濾波器參數以達到最佳的對比度增強效果，最后也使用閾值方法分割血管。Kerkeni等[12]提出了另一種基于多尺度區域生長（multiscale region growing，MSRG）的方法，該方法根據血管特性及血管方向信息設計了一種特殊的區域生長規則，基于此規則進行迭代多尺度搜索最終完成血管分割。通過引入的血管方向信息，MSRG能夠有效克服造影劑在血管內分布不均勻所導致的血管分割不連續的現象。對于CAG圖像優化處理后如何實現血管分割，也有很多研究者將目光投向了深度學習方法。Nasr-Esfahani等[13]提出的方法首先優化CAG圖像的對比度，然后利用一個兩階段的卷積神經網絡（convolutional neural network，CNN）實現血管分割。第一階段將CAG圖像分為若干像素塊作為第一階段CNN的輸入，輸出為與CAG圖像具有相同大小的概率圖，顯示每個像素屬于血管的概率。第一階段CNN輸出的概率圖和坎尼（Canny）邊緣檢測器的結果作為第二階段CNN的輸入，使得第二階段CNN更注意血管邊緣的位置，輸出更具魯棒性的概率圖，最后對概率圖進行閾值化得到最終的血管分割結果。Jun等[14]設計了一個T型網絡（T-Net）來分割CAG圖像中的主要血管，該網絡具有特殊的嵌套式編碼器—解碼器架構。具體來說，T-Net在編解碼過程中通過合理地使用池化和上采樣，能在單個塊中獲得各種大小的特征圖，并在編碼器和解碼器間共享提取到的各級特征。通過這樣的操作，即能夠克服U型網絡（U-Net）[15]在編碼器和解碼器之間只有一組連接層的限制，實現比U-Net更優的性能。
圖1
CAG圖像示例
Figure1.
Examples of CAG images
但是，在眾多需要CAG檢查的心血管患者中，有一些患者（例如心動過緩）安裝有心臟起搏器[16-17]。如圖1 所示，在這些患者的CAG圖像中，起搏器往往會與主要血管重疊在一起。在這種情況下，起搏器對于血管的自動分割會產生負面的影響。如果能夠移除重疊起搏器的干擾，優化血管分割結果，那么基于CAG圖像的CAD將對這些帶有起搏器的患者更加有益。起搏器的灰度分布和形狀與血管較為接近，而且由于心臟運動，起搏器和血管的相對位置也在時刻變化，這也增加移除重疊起搏器干擾的難度。為了克服這些困難，移除重疊起搏器干擾進而優化血管分割結果，本文提出了一種基于偽CAG圖像生成并與待分割CAG圖像配準后移除起搏器干擾的方法。在基于CAG的計算機輔助診斷系統中，該方法有潛力成為一個有益組件，提高其診斷準確性和效率。
1 方法
整個方法的流程圖如圖2所示，首先通過CAG視頻中的部分圖像序列生成一幅偽CAG圖像，偽CAG圖像中無造影劑，故無法觀察血管，但起搏器仍然清晰可見。而后，通過一種圖像局部特征描述子將偽CAG圖像和一幅待分割CAG圖像進行配準以確定起搏器在兩幅圖像中的相對位置。最后，通過“兩步走”策略分割偽CAG圖像中的起搏器和待分割CAG圖像中的血管，并結合配準結果從血管分割結果中移除重疊的起搏器，從而達到優化血管分割的目的。
圖2
基于偽CAG圖像生成并與待分割CAG圖像配準后移除起搏器干擾的方法
Figure2.
The removal of pacemaker interference based on the generation of pseudo CAG image and the registration between the generation result and the CAG image to be segmented
1.1 偽CAG圖像生成
完整的CAG視頻包含注入造影劑前、造影劑顯影到造影劑流出的全過程，其中包含豐富的信息。在本文所提出的方法中，基于注入造影劑前的CAG圖像序列(
)，生成無血管顯影僅有起搏器的偽CAG圖像。
首先，如式(1)所示，
包含n幀連續的CAG圖像序列。在定義降采樣因子p基礎上，自第一幀開始，每隔p幀采樣一幀圖像，得到圖像對應標號(i)為P的稀疏圖像序列，而后對稀疏圖像序列的灰度值求和、取平均得到圖像
。而后，如式(2)所示，對圖像
和稀疏圖像序列中的每一幀圖像逐一進行差分運算Diff(·)并求和，得到差分運動圖像Id。進行差分的主要目的是利用其能有效檢測運動的特性，突出圖像序列中包括起搏器在內的運動部分。最后，如式(3)所示，再將
與
相乘、求和得到偽CAG圖像(
)。這一步的目的在于使If能最大程度地整合Is中的運動部分，即起搏器的部分。
 |
 |
 |
其中，對于降采樣因子p的取值，應當根據CAG圖像序列的采集參數來決定。若CAG圖像序列的幀率較高，即幀間關聯性強，則降采樣因子p可以適當調大，以減少冗余信息，加快運算速度。而當幀率較低時，即幀間關聯性弱，此時降采樣因子p應適當調小，以獲取更多信息。
1.2 偽CAG圖像配準
利用局部特征描述子——尺度不變特征變換（scale invariant feature transform，SIFT）[18]對偽CAG圖像和待分割CAG圖像進行配準。SIFT具有很好的穩定性和不變性，因此廣泛應用于特征檢測和匹配工作中。主要包含四個步驟。第一步，如式(4)所示，利用高斯核函數G(x, y, σ)構造多尺度空間L(x, y, σ)：
 |
其中，I(x, y)為原圖像位于坐標(x, y)的像素灰度值，σ為尺度空間因子，
為卷積運算。第二步，在構造尺度空間中利用梯度算子求特征點。利用如式(5)所示的差分高斯D(x,y, σ)來求特征點并通過如式(6)所示計算對應尺度σ：
 |
 |
其中，σ0（設為1.6）為初始尺度，kS為尺度擴張系數，S（設為3）為尺度層數。
第三步，在求得特征點和尺度σ的基礎上就可以得到特征點所在的尺度圖像L(x, y)。然后按照如式(7)、式(8)所示，對以特征點為中心，以4.5σ為半徑的尺度圖像區域，求各像素點的梯度方向θ(x,y)和梯度幅值m(x, y)。
 |
 |
第四步，建立特征描述符。將特征點周圍區域像素點的梯度方向和梯度幅值進行直方圖統計，并求出直方圖的最大值，該最大值用來表征該區域像素點的梯度方向，該像素點對應的梯度幅值用來表征該區域的梯度幅值。配準兩幅圖像時，計算待配準的兩幅圖像中任意兩個特征點（P1和P2）的梯度方向和幅值之間的差值。在設定閾值T的情況下，當滿足θP1<TθP2和mP1<TmP2時，則判定P1和P2為配對的特征點。
1.3 血管分割優化
本文的分割方法同樣基于“兩步走”策略。首先，對偽CAG圖像和待分割CAG圖像進行對比度優化，使起搏器和血管能更容易從背景區域中區分出來，在此基礎上，再進行起搏器和導管的分割。本研究所用的對比度優化方法是一種基于伽馬校正的優化方法，自適應伽馬校正公式如式(9)~式(11)所示[19]：
 |
 |
 |
其中，l和lmax分別為圖像的灰度值和最大灰度值，cdf(l)為基于概率密度函數pdf的累積分布函數，pdfω為灰度值的權重(ω)分布函數，MN為圖像的像素點個數，nl為圖像中灰度值為l的像素點個數。另外，
為可調參數。這一方法有效克服了直方圖均衡類方法對部分灰度域的對比度欠優化或過度優化的缺點，能在整個灰度域上全面地提升CAG圖像的對比度。之后，再使用如式(12)、式(13)所示的濾波器進一步增強血管結構邊緣，該濾波器常用于血管圖像（如數字剪影血管造影）中，并取得了不錯的效果[20]。接下來，對圖像進行高斯濾波，再對圖像中的某一像素求取其黑塞（Hessian）矩陣（矩陣維度為2 × 2），其對應特征值為λ1和λ2，且λ1<λ2。Hessian矩陣特征值和對應的特征向量代表著該像素沿某一方向上圖形曲率大小和方向，λ1對應的是曲率較小的方向，λ2對應的是曲率較大的方向。在此基礎上，定義變量Rb和S如式(12)所示：
 |
并定義如式(13)所示血管結構響應函數：
 |
其中，β和c為控制濾波器敏感度的閾值，在本研究中設為0.5和15。又定義如式(14)所示的變量A和B：
 |
圖像中各像素可分為三類，即背景、孤立噪聲點和血管結構。其中，背景灰度分布較均勻，任意方向上曲率都較小；孤立噪聲點任意方向上曲率都很大；而圖像中的類圓柱形態的血管結構，沿徑向方向λ2上的曲率較大，沿軸向方向λ1上的曲率較小。因此，對于背景區域，|A|難以確定，|B|較小；對于孤立噪聲點，|A|趨于0，|B|較大；對于血管結構，則|A|趨于1，|B|較大。因此，弗蘭吉(Franji)濾波器能夠有效地區分血管結構、背景和孤立噪聲點，從而進一步實現對血管結構的增強，并抑制背景和孤立噪聲點。
對血管結構增強后的CAG圖像，通過模糊C均值聚類（FuzzyC-mean，FCM）方法[21]進行分割。將n個樣本分為C類時，FCM方法的目標函數如式(15)所示：
 |
其中，
表示第m次迭代時第i個樣本屬于第j類的隸屬度，xi表示第i個樣本，cj表示第j類的聚類中心。FCM方法是一個不斷迭代計算隸屬度并更新聚類中心的過程，其隸屬度和聚類中心的更新公式如式(16)～(17)所示：
 |
 |
對于本文的圖像分割任務而言，即對圖像做二分類任務，聚類數為兩類，xi為圖像中各像素點的灰度值。初始化時聚類中心cj為前景和背景區域的平均灰度值，而后如式(16)計算各像素的隸屬度，再如式(17)更新聚類中心，直到滿足如式(18)所示，即迭代后各個像素之間的隸屬度變化之差小于誤差閾值
，或者達到最大迭代次數，迭代終止。
 |
其中，k表示迭代步數。
至此，已完成起搏器分割和主要血管粗分割。而后，基于偽CAG圖像和待分割CAG圖像間配準結果，如式(19)計算起搏器位置偏移量：
 |
其中，
和
分別表示偽CAG圖像(
)和待分割CAG圖像(
)中配對的第m個SIFT特征點的
坐標，類似地，
和
對應
坐標。mall為配對的SIFT特征點數目，(Mx,My)表示起搏器位置偏移量。心臟搏動會導致冠狀動脈的形狀和位置發生較為復雜的變化，與之相比，起搏器的變化則相對明顯且簡單。因此，根據計算得到的偏移量，移動偽CAG圖像中起搏器的最小包圍矩形區域至新位置（待分割CAG圖像中的位置），就能在緊鄰位置上從血管分割結果中移除重疊起搏器的干擾，從而達到優化血管分割的目的。
2 實驗
2.1 實驗數據
本文實驗數據采集自復旦大學附屬中山醫院心內科，涉及5名患者，1名女性（雙腔永久起搏器置入）、4名男性，年齡在（58.25 ± 4.96）歲。其中，收集11段CAG視頻，均來自安裝起搏器的患者，所用成像系統為InnovaIGS 5（General Electric Inc.，美國）和IGS 6-GE（General Electric Inc.，美國）。從中選出11幅冠狀動脈結構顯示較為完整的CAG圖像進行測試，每幅圖像的大小均為512像素×512像素。此外，也收集了6段不安裝起搏器患者的CAG視頻，以相同標準選出6幅CAG圖像在文中展示，以說明所提出的“兩步走”策略能實現較有效的主要血管分割。本文所涉及的實驗均與患者簽署了知情同意書，研究整體通過復旦大學附屬中山醫院倫理委員會的審查并獲得授權。
2.2 實驗設計
由于冠狀動脈結構復雜，難以有效、準確地進行血管標注，且不同觀察者間往往存在比較大的觀察誤差。因此，為了避免報告具有誤導性的血管分割結果，本研究通過下述方式來設置金標準以評估所提出方法的可行性。首先，對11幅圖像直接用“兩步走”策略進行血管分割。相比于血管，起搏器更容易辨識，因此，將分割結果中重疊起搏器的干擾擦除，將其作為金標準使用。
2.3 評價方法
采用戴斯系數（Dice coefficient）（以符號Dice表示）評估本文提出的方法去除重疊起搏器干擾、優化主要血管分割的性能。戴斯系數計算如式(20)所示：
 |
其中，真陽性（true positive，TP）像素的數量表示為TP，假陰性（false negative，FN）像素的數量表示為FN，假陽性（false positive，FP）像素的數量表示為FP。因為金標準生成方法特殊，且本研究更關注去除重疊起搏器干擾對血管分割結果優化的性能，所以本研究只計算戴斯系數的變化率，即將直接分割的血管分割結果和優化后的血管分割結果分別與金標準間計算評價參數，再計算兩者之間的差異百分比，以此來評價所提出方法的性能。
2.4 實驗環境
所有實驗均在相同的環境下進行。軟件環境為MATLAB 2019a（Mathwork Inc.，美國）；硬件環境為Intel? Xeon? Silver 4215CPU 2.50 GHz（Intel Inc.，美國）。參數設置上，降采樣因子p設置為5；尺度空間因子σ為1.6，S設置為3，閾值T設置為0.2；可調閾值α設置為0.3；聚類數C設置為2，閾值ε設置為1×10?5，最大迭代次數設置為100次。
3 結果與討論
如圖3所示，隨機選取的圖A~圖F中生成的偽CAG圖像中起搏器保留完整，而且與待分割CAG圖像相比，除了由心臟運動而導致的差異外，其余部分仍較為相似，這為后續進行配準和分割優化打下了非常好的基礎。對于后續的偽CAG圖像與待分割CAG圖像配準，SIFT最多定位到29個匹配特征點，最少定位到1個，平均定位到5個。
圖3
生成偽CAG圖像結果與待分割CAG圖像示例
Figure3.
Examples of pseudo CAG images and CAG images for segmentation
如圖4所示，靛藍色實線指示出生成的偽CAG圖像與待分割CAG圖像之間成功配對的特征點，基于這些匹配的特征點，SIFT能夠確定起搏器在兩幅圖像中的相對位置關系。如圖5所示，以隨機選取的圖A~圖F為例，從所提出的“兩步走”分割策略來看，對于偽CAG圖像中起搏器的分割而言，“兩步走”策略能有效地分割出其中的起搏器。另一方面，如圖6所示，以隨機選取的圖G~圖L為例，對于主要血管的分割，分割結果較為準確，有效地保留血管的主要部分，能夠滿足諸如血管中心線提取在內的后續應用要求。總體而言，“兩步走”策略能夠有效實現起搏器和主要血管分割。從整個方法的層面來看，如圖7所示，隨機選取的圖M~圖R為例，為優化前后主要血管分割結果的例子。從中可以發現，重疊起搏器干擾被有效地移除，主要血管結構保留完整，達到了預想的優化血管分割的目的。
圖4
SIFT配準偽CAG圖像和待分割CAG圖像結果示例
Figure4.
Examples of SIFT registration results between pseudo CAG images and CAG images for segmentation
圖5
“兩步走”策略實現起搏器分割結果示例
Figure5.
Examples of pacemaker detection results based on “two-step” strategy
圖6
“兩步走”策略實現主要血管分割結果示例
Figure6.
Examples of main vessel segmentation results based on “two-step” strategy
圖7
CAG圖像粗分割、優化分割結果及金標準示例
Figure7.
Examples of rough and improved vessel segmentation results and corresponding labels
整體上，本文的方法去除重疊起搏器干擾后所獲得的血管分割結果的戴斯系數平均優化率達到了12.03%，最大優化率可以達到41.55%。這證明本文方法能夠通過去除重疊起搏器干擾，有效地提升血管分割的準確度。從執行效率來看，本文方法實現血管分割并優化的平均耗時約為54 s，這是由于SIFT特征需要遍歷整幅圖像來搜索極值點，因此產生了一定的計算耗時。雖然已有的CAG圖像血管分割方法已經實現了較優的性能，但是這些方法所使用的數據集中，往往沒有安裝心臟起搏器患者的圖像數據，針對這類情況的血管分割優化方法也鮮有研究。在本文中，從臨床中遇到的實際情況出發，結合CAG檢查的特點，提出了優化帶有重疊起搏器干擾的CAG圖像的血管分割方法。從實驗結果來看，本文方法能有效從血管分割結果中移除重疊起搏器的干擾，從而達到優化血管分割結果的目的。但是，本文方法也存在一些局限。對于過短的CAG視頻，造影劑注入前的CAG圖像序列也較短，這就導致在偽CAG圖像生成的過程中難以獲取足夠的多幀信息，生成圖像難以接近于待分割CAG圖像的灰度分布。另外，一些CAG視頻中會存在嚴重的噪聲干擾，這導致生成的偽CAG圖像和待分割CAG圖像均不可避免地會存在干擾。這種情況下，起搏器和主要血管的分割將變得更困難。
4 結論
為了更有效地實現對于帶有重疊起搏器干擾的CAG圖像的血管分割，本文提出了一種基于偽CAG圖像生成并與待分割CAG圖像配準后移除起搏器干擾的方法。首先利用CAG圖像序列信息生成帶有起搏器且與待分割圖像類似的偽CAG圖像，而后通過SIFT配準起搏器在兩幅CAG圖像中的相對位置，最后使用“兩步走”策略實現對血管和起搏器分割，并基于配準結果將重疊起搏器從血管分割結果中移除。在臨床獲取的數據上進行了實驗，其結果表明本文方法對于無重疊起搏器干擾的圖像能有效地分割出血管的主要部分，對于有重疊起搏器干擾的圖像則能有效地優化最終的分割結果。
由于SIFT特征配準計算耗時較大，未來的研究中本課題組將考慮以更為高效的配準方法來實現起搏器的相對位置配準，從而進一步提高方法的執行效率。除了相對位置外，也可以嘗試對起搏器的形變進行分析和配準。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：黃藝、楊虹波、夏夢華負責方法設計和論文撰寫；曲亞男負責數據收集、整理和方法測試；郭翌、周國輝、張峰、汪源源指導方法理論和實驗設計，并指導論文撰寫。
倫理聲明：本研究通過了復旦大學附屬中山醫院倫理委員會的審批［批文編號：B2019-251R］。
引言
近年來，全球心血管疾病患者數持續增加，中國雖然在心血管疾病防治工作上取得了一些成績，但據估計，中國現有心血管疾病患者數仍高達3.3億。其中，由冠狀動脈發生粥樣硬化、血管腔狹窄，進而造成心肌缺血、缺氧乃至壞死的一類心臟病，稱為冠狀動脈粥樣硬化性心臟病（后文簡稱：冠心病）。冠心病作為一種典型的心臟病，其患者數達到1100萬[1-3]。冠狀動脈造影（coronaryangiography，CAG）將導管通過外周動脈血管送至造影位點，然后注射少量造影劑再利用X線對心臟血管進行顯像，其作為冠心病診斷的有效手段，在臨床上得到了廣泛的應用[4]。自動分割CAG圖像中的冠狀動脈對于冠心病的計算機輔助診斷（computer-aided diagnosis，CAD）而言是一個基礎，在此基礎上進行狹窄評估[5]、三維重建[6]和動態評估[7]，可為醫生診斷提供參考。此外，CAG還可以與其他檢查手段相結合，如血管內超聲（intravascularultrasound，IVUS）[8]、血管光學相干斷層掃描成像（opticalcoherence tomography，OCT）[9]，從而對血管狀況進行更全面的評估。
如圖1所示，CAG圖像中前景（血管）和背景區域之間的對比度較低，并且造影劑會隨著血液流動而產生分布不均的現象，使得血管區域的灰度分布也不盡相同，這使得分割CAG圖像中的主要血管存在一定困難。為了有效地分割CAG圖像中的血管，研究者們做了許多的探索，形成較為常用且有效的“兩步走”策略，即先對CAG圖像進行優化處理，主要目的是提高血管與背景區域之間的對比度，使得兩者之間更容易區分，然后再對優化后的CAG圖像進行血管分割。在“兩步走”策略的基礎上，涌現出許多優秀的方法。Cruz-Aceves等[10]提出了一種兩階段方法，首先使用多尺度伽柏（Gabor）濾波器和頂帽（Top-hat）運算優化處理圖像，而后使用基于多目標優化的新閾值方法對優化結果進行二值化，實現對冠狀動脈的分割。這一方法的優點在于多尺度Gabor濾波器能夠有效提升灰度不均勻CAG圖像的對比度，且閾值方法中同時引入了三種加權優化準則，大大提高了分割的魯棒性。Cervantes-Sanchez等[11]同樣使用了Gabor濾波器，并結合差分演化的優化策略選擇濾波器參數以達到最佳的對比度增強效果，最后也使用閾值方法分割血管。Kerkeni等[12]提出了另一種基于多尺度區域生長（multiscale region growing，MSRG）的方法，該方法根據血管特性及血管方向信息設計了一種特殊的區域生長規則，基于此規則進行迭代多尺度搜索最終完成血管分割。通過引入的血管方向信息，MSRG能夠有效克服造影劑在血管內分布不均勻所導致的血管分割不連續的現象。對于CAG圖像優化處理后如何實現血管分割，也有很多研究者將目光投向了深度學習方法。Nasr-Esfahani等[13]提出的方法首先優化CAG圖像的對比度，然后利用一個兩階段的卷積神經網絡（convolutional neural network，CNN）實現血管分割。第一階段將CAG圖像分為若干像素塊作為第一階段CNN的輸入，輸出為與CAG圖像具有相同大小的概率圖，顯示每個像素屬于血管的概率。第一階段CNN輸出的概率圖和坎尼（Canny）邊緣檢測器的結果作為第二階段CNN的輸入，使得第二階段CNN更注意血管邊緣的位置，輸出更具魯棒性的概率圖，最后對概率圖進行閾值化得到最終的血管分割結果。Jun等[14]設計了一個T型網絡（T-Net）來分割CAG圖像中的主要血管，該網絡具有特殊的嵌套式編碼器—解碼器架構。具體來說，T-Net在編解碼過程中通過合理地使用池化和上采樣，能在單個塊中獲得各種大小的特征圖，并在編碼器和解碼器間共享提取到的各級特征。通過這樣的操作，即能夠克服U型網絡（U-Net）[15]在編碼器和解碼器之間只有一組連接層的限制，實現比U-Net更優的性能。
圖1
CAG圖像示例
Figure1.
Examples of CAG images
但是，在眾多需要CAG檢查的心血管患者中，有一些患者（例如心動過緩）安裝有心臟起搏器[16-17]。如圖1 所示，在這些患者的CAG圖像中，起搏器往往會與主要血管重疊在一起。在這種情況下，起搏器對于血管的自動分割會產生負面的影響。如果能夠移除重疊起搏器的干擾，優化血管分割結果，那么基于CAG圖像的CAD將對這些帶有起搏器的患者更加有益。起搏器的灰度分布和形狀與血管較為接近，而且由于心臟運動，起搏器和血管的相對位置也在時刻變化，這也增加移除重疊起搏器干擾的難度。為了克服這些困難，移除重疊起搏器干擾進而優化血管分割結果，本文提出了一種基于偽CAG圖像生成并與待分割CAG圖像配準后移除起搏器干擾的方法。在基于CAG的計算機輔助診斷系統中，該方法有潛力成為一個有益組件，提高其診斷準確性和效率。
1 方法
整個方法的流程圖如圖2所示，首先通過CAG視頻中的部分圖像序列生成一幅偽CAG圖像，偽CAG圖像中無造影劑，故無法觀察血管，但起搏器仍然清晰可見。而后，通過一種圖像局部特征描述子將偽CAG圖像和一幅待分割CAG圖像進行配準以確定起搏器在兩幅圖像中的相對位置。最后，通過“兩步走”策略分割偽CAG圖像中的起搏器和待分割CAG圖像中的血管，并結合配準結果從血管分割結果中移除重疊的起搏器，從而達到優化血管分割的目的。
圖2
基于偽CAG圖像生成并與待分割CAG圖像配準后移除起搏器干擾的方法
Figure2.
The removal of pacemaker interference based on the generation of pseudo CAG image and the registration between the generation result and the CAG image to be segmented
1.1 偽CAG圖像生成
完整的CAG視頻包含注入造影劑前、造影劑顯影到造影劑流出的全過程，其中包含豐富的信息。在本文所提出的方法中，基于注入造影劑前的CAG圖像序列()，生成無血管顯影僅有起搏器的偽CAG圖像。
首先，如式(1)所示，包含n幀連續的CAG圖像序列。在定義降采樣因子p基礎上，自第一幀開始，每隔p幀采樣一幀圖像，得到圖像對應標號(i)為P的稀疏圖像序列，而后對稀疏圖像序列的灰度值求和、取平均得到圖像。而后，如式(2)所示，對圖像和稀疏圖像序列中的每一幀圖像逐一進行差分運算Diff(·)并求和，得到差分運動圖像Id。進行差分的主要目的是利用其能有效檢測運動的特性，突出圖像序列中包括起搏器在內的運動部分。最后，如式(3)所示，再將與相乘、求和得到偽CAG圖像()。這一步的目的在于使If能最大程度地整合Is中的運動部分，即起搏器的部分。
|
|
|
其中，對于降采樣因子p的取值，應當根據CAG圖像序列的采集參數來決定。若CAG圖像序列的幀率較高，即幀間關聯性強，則降采樣因子p可以適當調大，以減少冗余信息，加快運算速度。而當幀率較低時，即幀間關聯性弱，此時降采樣因子p應適當調小，以獲取更多信息。
1.2 偽CAG圖像配準
利用局部特征描述子——尺度不變特征變換（scale invariant feature transform，SIFT）[18]對偽CAG圖像和待分割CAG圖像進行配準。SIFT具有很好的穩定性和不變性，因此廣泛應用于特征檢測和匹配工作中。主要包含四個步驟。第一步，如式(4)所示，利用高斯核函數G(x, y, σ)構造多尺度空間L(x, y, σ)：
|
其中，I(x, y)為原圖像位于坐標(x, y)的像素灰度值，σ為尺度空間因子，為卷積運算。第二步，在構造尺度空間中利用梯度算子求特征點。利用如式(5)所示的差分高斯D(x,y, σ)來求特征點并通過如式(6)所示計算對應尺度σ：
|
|
其中，σ0（設為1.6）為初始尺度，kS為尺度擴張系數，S（設為3）為尺度層數。
第三步，在求得特征點和尺度σ的基礎上就可以得到特征點所在的尺度圖像L(x, y)。然后按照如式(7)、式(8)所示，對以特征點為中心，以4.5σ為半徑的尺度圖像區域，求各像素點的梯度方向θ(x,y)和梯度幅值m(x, y)。
|
|
第四步，建立特征描述符。將特征點周圍區域像素點的梯度方向和梯度幅值進行直方圖統計，并求出直方圖的最大值，該最大值用來表征該區域像素點的梯度方向，該像素點對應的梯度幅值用來表征該區域的梯度幅值。配準兩幅圖像時，計算待配準的兩幅圖像中任意兩個特征點（P1和P2）的梯度方向和幅值之間的差值。在設定閾值T的情況下，當滿足θP1<TθP2和mP1<TmP2時，則判定P1和P2為配對的特征點。
1.3 血管分割優化
本文的分割方法同樣基于“兩步走”策略。首先，對偽CAG圖像和待分割CAG圖像進行對比度優化，使起搏器和血管能更容易從背景區域中區分出來，在此基礎上，再進行起搏器和導管的分割。本研究所用的對比度優化方法是一種基于伽馬校正的優化方法，自適應伽馬校正公式如式(9)~式(11)所示[19]：
|
|
|
其中，l和lmax分別為圖像的灰度值和最大灰度值，cdf(l)為基于概率密度函數pdf的累積分布函數，pdfω為灰度值的權重(ω)分布函數，MN為圖像的像素點個數，nl為圖像中灰度值為l的像素點個數。另外，為可調參數。這一方法有效克服了直方圖均衡類方法對部分灰度域的對比度欠優化或過度優化的缺點，能在整個灰度域上全面地提升CAG圖像的對比度。之后，再使用如式(12)、式(13)所示的濾波器進一步增強血管結構邊緣，該濾波器常用于血管圖像（如數字剪影血管造影）中，并取得了不錯的效果[20]。接下來，對圖像進行高斯濾波，再對圖像中的某一像素求取其黑塞（Hessian）矩陣（矩陣維度為2 × 2），其對應特征值為λ1和λ2，且λ1<λ2。Hessian矩陣特征值和對應的特征向量代表著該像素沿某一方向上圖形曲率大小和方向，λ1對應的是曲率較小的方向，λ2對應的是曲率較大的方向。在此基礎上，定義變量Rb和S如式(12)所示：
|
并定義如式(13)所示血管結構響應函數：
|
其中，β和c為控制濾波器敏感度的閾值，在本研究中設為0.5和15。又定義如式(14)所示的變量A和B：
|
圖像中各像素可分為三類，即背景、孤立噪聲點和血管結構。其中，背景灰度分布較均勻，任意方向上曲率都較小；孤立噪聲點任意方向上曲率都很大；而圖像中的類圓柱形態的血管結構，沿徑向方向λ2上的曲率較大，沿軸向方向λ1上的曲率較小。因此，對于背景區域，|A|難以確定，|B|較小；對于孤立噪聲點，|A|趨于0，|B|較大；對于血管結構，則|A|趨于1，|B|較大。因此，弗蘭吉(Franji)濾波器能夠有效地區分血管結構、背景和孤立噪聲點，從而進一步實現對血管結構的增強，并抑制背景和孤立噪聲點。
對血管結構增強后的CAG圖像，通過模糊C均值聚類（FuzzyC-mean，FCM）方法[21]進行分割。將n個樣本分為C類時，FCM方法的目標函數如式(15)所示：
|
其中，表示第m次迭代時第i個樣本屬于第j類的隸屬度，xi表示第i個樣本，cj表示第j類的聚類中心。FCM方法是一個不斷迭代計算隸屬度并更新聚類中心的過程，其隸屬度和聚類中心的更新公式如式(16)～(17)所示：
|
|
對于本文的圖像分割任務而言，即對圖像做二分類任務，聚類數為兩類，xi為圖像中各像素點的灰度值。初始化時聚類中心cj為前景和背景區域的平均灰度值，而后如式(16)計算各像素的隸屬度，再如式(17)更新聚類中心，直到滿足如式(18)所示，即迭代后各個像素之間的隸屬度變化之差小于誤差閾值，或者達到最大迭代次數，迭代終止。
|
其中，k表示迭代步數。
至此，已完成起搏器分割和主要血管粗分割。而后，基于偽CAG圖像和待分割CAG圖像間配準結果，如式(19)計算起搏器位置偏移量：
|
其中，和分別表示偽CAG圖像()和待分割CAG圖像()中配對的第m個SIFT特征點的坐標，類似地，和對應坐標。mall為配對的SIFT特征點數目，(Mx,My)表示起搏器位置偏移量。心臟搏動會導致冠狀動脈的形狀和位置發生較為復雜的變化，與之相比，起搏器的變化則相對明顯且簡單。因此，根據計算得到的偏移量，移動偽CAG圖像中起搏器的最小包圍矩形區域至新位置（待分割CAG圖像中的位置），就能在緊鄰位置上從血管分割結果中移除重疊起搏器的干擾，從而達到優化血管分割的目的。
2 實驗
2.1 實驗數據
本文實驗數據采集自復旦大學附屬中山醫院心內科，涉及5名患者，1名女性（雙腔永久起搏器置入）、4名男性，年齡在（58.25 ± 4.96）歲。其中，收集11段CAG視頻，均來自安裝起搏器的患者，所用成像系統為InnovaIGS 5（General Electric Inc.，美國）和IGS 6-GE（General Electric Inc.，美國）。從中選出11幅冠狀動脈結構顯示較為完整的CAG圖像進行測試，每幅圖像的大小均為512像素×512像素。此外，也收集了6段不安裝起搏器患者的CAG視頻，以相同標準選出6幅CAG圖像在文中展示，以說明所提出的“兩步走”策略能實現較有效的主要血管分割。本文所涉及的實驗均與患者簽署了知情同意書，研究整體通過復旦大學附屬中山醫院倫理委員會的審查并獲得授權。
2.2 實驗設計
由于冠狀動脈結構復雜，難以有效、準確地進行血管標注，且不同觀察者間往往存在比較大的觀察誤差。因此，為了避免報告具有誤導性的血管分割結果，本研究通過下述方式來設置金標準以評估所提出方法的可行性。首先，對11幅圖像直接用“兩步走”策略進行血管分割。相比于血管，起搏器更容易辨識，因此，將分割結果中重疊起搏器的干擾擦除，將其作為金標準使用。
2.3 評價方法
采用戴斯系數（Dice coefficient）（以符號Dice表示）評估本文提出的方法去除重疊起搏器干擾、優化主要血管分割的性能。戴斯系數計算如式(20)所示：
|
其中，真陽性（true positive，TP）像素的數量表示為TP，假陰性（false negative，FN）像素的數量表示為FN，假陽性（false positive，FP）像素的數量表示為FP。因為金標準生成方法特殊，且本研究更關注去除重疊起搏器干擾對血管分割結果優化的性能，所以本研究只計算戴斯系數的變化率，即將直接分割的血管分割結果和優化后的血管分割結果分別與金標準間計算評價參數，再計算兩者之間的差異百分比，以此來評價所提出方法的性能。
2.4 實驗環境
所有實驗均在相同的環境下進行。軟件環境為MATLAB 2019a（Mathwork Inc.，美國）；硬件環境為Intel? Xeon? Silver 4215CPU 2.50 GHz（Intel Inc.，美國）。參數設置上，降采樣因子p設置為5；尺度空間因子σ為1.6，S設置為3，閾值T設置為0.2；可調閾值α設置為0.3；聚類數C設置為2，閾值ε設置為1×10?5，最大迭代次數設置為100次。
3 結果與討論
如圖3所示，隨機選取的圖A~圖F中生成的偽CAG圖像中起搏器保留完整，而且與待分割CAG圖像相比，除了由心臟運動而導致的差異外，其余部分仍較為相似，這為后續進行配準和分割優化打下了非常好的基礎。對于后續的偽CAG圖像與待分割CAG圖像配準，SIFT最多定位到29個匹配特征點，最少定位到1個，平均定位到5個。
圖3
生成偽CAG圖像結果與待分割CAG圖像示例
Figure3.
Examples of pseudo CAG images and CAG images for segmentation
如圖4所示，靛藍色實線指示出生成的偽CAG圖像與待分割CAG圖像之間成功配對的特征點，基于這些匹配的特征點，SIFT能夠確定起搏器在兩幅圖像中的相對位置關系。如圖5所示，以隨機選取的圖A~圖F為例，從所提出的“兩步走”分割策略來看，對于偽CAG圖像中起搏器的分割而言，“兩步走”策略能有效地分割出其中的起搏器。另一方面，如圖6所示，以隨機選取的圖G~圖L為例，對于主要血管的分割，分割結果較為準確，有效地保留血管的主要部分，能夠滿足諸如血管中心線提取在內的后續應用要求。總體而言，“兩步走”策略能夠有效實現起搏器和主要血管分割。從整個方法的層面來看，如圖7所示，隨機選取的圖M~圖R為例，為優化前后主要血管分割結果的例子。從中可以發現，重疊起搏器干擾被有效地移除，主要血管結構保留完整，達到了預想的優化血管分割的目的。
圖4
SIFT配準偽CAG圖像和待分割CAG圖像結果示例
Figure4.
Examples of SIFT registration results between pseudo CAG images and CAG images for segmentation
圖5
“兩步走”策略實現起搏器分割結果示例
Figure5.
Examples of pacemaker detection results based on “two-step” strategy
圖6
“兩步走”策略實現主要血管分割結果示例
Figure6.
Examples of main vessel segmentation results based on “two-step” strategy
圖7
CAG圖像粗分割、優化分割結果及金標準示例
Figure7.
Examples of rough and improved vessel segmentation results and corresponding labels
整體上，本文的方法去除重疊起搏器干擾后所獲得的血管分割結果的戴斯系數平均優化率達到了12.03%，最大優化率可以達到41.55%。這證明本文方法能夠通過去除重疊起搏器干擾，有效地提升血管分割的準確度。從執行效率來看，本文方法實現血管分割并優化的平均耗時約為54 s，這是由于SIFT特征需要遍歷整幅圖像來搜索極值點，因此產生了一定的計算耗時。雖然已有的CAG圖像血管分割方法已經實現了較優的性能，但是這些方法所使用的數據集中，往往沒有安裝心臟起搏器患者的圖像數據，針對這類情況的血管分割優化方法也鮮有研究。在本文中，從臨床中遇到的實際情況出發，結合CAG檢查的特點，提出了優化帶有重疊起搏器干擾的CAG圖像的血管分割方法。從實驗結果來看，本文方法能有效從血管分割結果中移除重疊起搏器的干擾，從而達到優化血管分割結果的目的。但是，本文方法也存在一些局限。對于過短的CAG視頻，造影劑注入前的CAG圖像序列也較短，這就導致在偽CAG圖像生成的過程中難以獲取足夠的多幀信息，生成圖像難以接近于待分割CAG圖像的灰度分布。另外，一些CAG視頻中會存在嚴重的噪聲干擾，這導致生成的偽CAG圖像和待分割CAG圖像均不可避免地會存在干擾。這種情況下，起搏器和主要血管的分割將變得更困難。
4 結論
為了更有效地實現對于帶有重疊起搏器干擾的CAG圖像的血管分割，本文提出了一種基于偽CAG圖像生成并與待分割CAG圖像配準后移除起搏器干擾的方法。首先利用CAG圖像序列信息生成帶有起搏器且與待分割圖像類似的偽CAG圖像，而后通過SIFT配準起搏器在兩幅CAG圖像中的相對位置，最后使用“兩步走”策略實現對血管和起搏器分割，并基于配準結果將重疊起搏器從血管分割結果中移除。在臨床獲取的數據上進行了實驗，其結果表明本文方法對于無重疊起搏器干擾的圖像能有效地分割出血管的主要部分，對于有重疊起搏器干擾的圖像則能有效地優化最終的分割結果。
由于SIFT特征配準計算耗時較大，未來的研究中本課題組將考慮以更為高效的配準方法來實現起搏器的相對位置配準，從而進一步提高方法的執行效率。除了相對位置外，也可以嘗試對起搏器的形變進行分析和配準。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：黃藝、楊虹波、夏夢華負責方法設計和論文撰寫；曲亞男負責數據收集、整理和方法測試；郭翌、周國輝、張峰、汪源源指導方法理論和實驗設計，并指導論文撰寫。
倫理聲明：本研究通過了復旦大學附屬中山醫院倫理委員會的審批［批文編號：B2019-251R］。
