針對冠脈造影圖像,本文提出了新的血管分割算法和血管結構識別算子。血管分割中采用了新的血管函數獲取血管特征圖,并在區域增長算法中實施了自動化的種子點選擇、主分支提取、血管細節修復等步驟。血管結構識別中采用了模糊算子,可以判別血管段、分支、交叉、末梢等結構。實驗結果表明:新分割方法能夠提取的血管面積較單純區域增長算法平均增加了5%左右,并能夠從低灰度區域中恢復出一些較細的血管段。實驗對模糊識別算子驗證后得出:對仿真血管結構的識別正確率達100%。對于實際分割出的各級血管結構進行判別,總體識別的正確率平均達90.59%。
引用本文： 梅川, 吳桂良, 楊媛, 謝斕, 何家駒, 李紹林, 周壽軍. 一種基于區域增長和結構識別的心血管X線造影圖像分割方法. 生物醫學工程學雜志, 2014, 31(2): 413-420. doi: 10.7507/1001-5515.20140077 復制
引言
分割并描述X線造影圖像中的冠狀動脈樹，對血管的三維重建、測量和診斷具有重要意義。自20世紀90年代以來，X線在醫學領域受到廣泛關注。X線冠脈造影技術被認為是冠脈疾病診斷的“金標準”。X線造影成像是利用空間X線錐形束穿過目標區域后得到投影圖像，受限于成像條件和環境，造影圖像往往受到呼吸、心臟搏動造成的偽影和各種噪聲(如環境噪聲和人體不同組織對X線的不同衰減作用)干擾。
造影圖像中除了能夠清晰觀察的血管之外，還存在骨骼、脊柱、軟組織等，這些結構形態以及造影劑的非均勻顯像效果使得血管的自動檢測和計算機輔助介入過程比較困難。因此，有必要利用圖像處理技術獲取精確的血管結構并提取血管。血管造影圖像的分割[1-2]涉及血管圖像預處理、動脈樹軸線估計、結構和寬度測量等內容，其中預處理用來增強血管目標、降低背景和噪聲的影響。各種血管分割技術中[1-2]，提取血管特征并利用血管先驗知識構造算法模型能夠有效地分割出造影圖像中的血管樹。然而，由于造影劑的分布不均和成像過程中的環境噪聲以及復雜的血管形態(尤其是血管節點的多樣性和大曲率的血管段)，血管樹的分割具有極大的挑戰性。
近來有很多血管造影圖像分割方面的報道，經典的算法包括直接跟蹤法[3-4]、動態輪廓法[5-6]、中心線優化法[7-8]、Markov標記點過程[9]等。文獻[1-2]總體回顧和評論了上述方法。基于血管跟蹤的方法進行血管分割時，通常采用一個局部算子作用于血管樹的某個初始點并開始跟蹤，執行過程中周而復始地進行測量與決策。區域增長算法用于血管分割的基本思想是從一些血管區域內的種子點和子區域開始，按照某種排除準則，不斷地接受新的血管像素對當前區域進行擴充。通常有兩種分割準則：數值相似性和空間相近[10]。采用基于數值相似性的區域增長算法時，簡單的方法是利用一個單純的閾值法則[11]，即：假使像素具有相近的強度時很可能屬于同一目標，優化的閾值選擇需經過試驗確認。該方法的缺點是易受噪聲和對比劑不均勻性的干擾，常常產生過分割或者孔洞效果。空間一致性條件下的波擴張技術[10-11]是基于空間相似性的典型區域增長方法，也可以被當作一種順序區域增長策略，其中候選像素的排除與否依據了增長波幾何校正的局部策略與擬合準則。空間一致性擴張使得拓撲分析更加容易。例如，擴張界面的連續性可以被用來探測血管分支。此外，在X線造影圖像中識別冠狀動脈樹結構，對于定性和定量描述CAT結構、三維幾何重建、計算機輔助血管介入治療也具有重要作用。本文的血管結構判別算子作用于區域增長過程中順序獲得的骨架線，并通過特征測量進行結構的判別；測量中利用了血管預處理、特征監測、先驗知識等手段進行識別。首先，利用Hessian特征值分析算法[12]和方向濾波器[13-14]技術，經過多尺度血管特征提取和增強后獲得血管圖、矢量場。其次，結合Gabor能量描述子[14]，血管的形狀和方向特征能夠被準確地測量和描述。
本文的主要目的在于：① 建立在多尺度濾波之上的新血管函數，為區域增長提供更好的特征圖；② 提出基于細節修復的血管分割方法；③ 提出基于特征測量的血管結構識別算法。
1 造影圖像預處理
1.1 多尺度血管增強
我們可以將血管看做彎曲的帶狀目標，其直徑沿彎曲的主軸漸進變化。由于血管直徑分布在有限的范圍，我們可以采用尺度σ相當于血管半徑的一系列高斯核函數G(x,y,σ)卷積原始圖像IO(x,y)從而得到不同尺度下的血管響應：
| $\begin{align} & {{I}_{\sigma }}\left( x,y;\sigma \right)={{I}_{O}}\left( x,y \right)\otimes G\left( x,y;\sigma \right), \\ & 其中G\left( x,y;\sigma \right)=\frac{1}{\sqrt{2\pi \sigma }}{{e}^{\frac{-{{x}^{2}}+{{y}^{2}}}{2{{\sigma }^{2}}}}} \\ \end{align}$ |
利用上式對圖像卷積后，寬度相當于σ的血管結構獲得較大增強，反之則被抑制。圖像點p0的局部特性可以由該點處的Hessian矩陣分析得到，Hessian矩陣定義為
| $H({{p}_{0}},\sigma )={{\sigma }^{2}}*\left\{ -{{\partial }^{2}}G\left( x,y;\sigma \right) \right\}\otimes {{I}_{O}},$ |
其中σ2用來歸一化微分結果。當血管相對背景為低暗目標時，對應于血管點處的Hessian矩陣具有較大的負特征值λ2和一個絕對值較小的特征值λ1，且滿足關系。Frangi等[12]曾提出一種包含三個參數的血管測度函數，然而這些參數值不易確定。我們對比了當前各種血管函數[12-14],并提出新的血管函數：
| $\begin{align} & f({{p}_{0}})=\underset{{{\sigma }_{min}}\le \sigma \le {{\sigma }_{max}}}{\mathop{\max }}\,f({{p}_{0}},\sigma )/exp\sqrt{2\sigma },其中f({{p}_{0}},\sigma ) \\ & \left\{ \begin{array}{*{35}{l}} {{\lambda }_{2}} & {{\lambda }_{2}}>{{\sigma }^{2}}/4 \\ 0 & 其它 \\ \end{array} \right. \\ \end{align}$ |
利用公式(3) 的血管函數作用于整幅血管造影圖像后，不妨稱為“血管特征圖”。我們利用尺度參數σ2/4作為閾值。σmin和σmax為最小和最大尺度，其范圍的選擇可依據圖像中血管的最小和最大直徑確定。f(p0)為血管特征圖。圖 1中描述了新血管函數和傳統的Frangi函數增強血管的效果對比。明顯可見，新函數對較粗的和較細的血管具有相似的增強效果。該血管函數利用了特征圖閾值，能夠為后續的區域增算法提供自適應的初始結果。
圖1
不同血管函數對造影圖像增強效果比較
(a)原始造影圖像；(b)Frangi的血管函數增強效果；(c)利用新血管函數得到的血管特征圖
Figure1. Comparison of the angiography image enhancement based on different vessel functions(a) original angiography image; (b) result from Frangi’s vessel function; (c) result from the proposed vessel function
1.2 Gabor濾波器響應
在血管特征測量與結構識別算法設計中，利用Gabor濾波器[14]作用于血管軸線上的各個觀察點，濾波器響應最大的方向便為局部血管方向，從而有效避免了上述問題對結構判別的影響。Gabor濾波器函數為
| $\psi ({{X}_{0}})=\frac{\sum{{{W}_{i,j}}\left( \omega ,\theta ,S \right)\centerdot {{I}_{＜{{X}_{0}}>s}}\left( i,j \right)}}{\sum{{{W}_{i,j}}\left( \omega ,\theta ,S \right)}},$ |
其中Wi,j(ω,θ,S)為Gabor掩膜函數，具有可變的頻率帶寬ω、方向θ和尺度S=M×N，且有-ZX0≤Wi,j(ω,θ,S)≤ZX0。因此，式(4) 滿足：-1＜ψ(X0)＜1。掩膜與圖像進行相關時，每一點對應尺寸的區域中血管的特性如果和掩膜相近，則有最大響應輸出。令，式(4) 可歸一化為0~1之間，最大響應輸出表示為：，其中θ取血管方向場中點的X0對應值V2(X0)時，可極大地減少計算時間。實際處理中，作為預處理手段，血管相似性函數作用于整幅圖像分析血管的亮度特征和方向矢量；而Gabor響應函數作用于血管軸線上的各個觀察點，并測量血管的局部強度和方向響應。
2 區域增長和細節修復
為了分割出血管樹，我們建立了一種全自動的區域增長方法。主要過程為：自動種子點選擇、主分支分割、血管細節修復等。在進行種子點的選擇時，理論上可以自動搜索血管特征圖中的最大強度點，該點自然位于血管區域內。然而對實際造影圖像進行處理時，有時會遇到X線造影圖像存在邊界偽影(見圖 2a)，致使最大特征強度點位于造影圖像中血管區域之外的邊界，產生錯誤結果。因此我們通常利用圖 2(b)所示的掩膜保證自動搜索的正確性。圖 2(c)為特征圖和種子點自動選擇的位置(種子點位于下方的圓圈中心處)。
圖2
種子點的選擇和多尺度血管濾波結果
(a)原始造影圖像；(b)掩膜；(c)血管特征圖；(d)多尺度響應
(a) original angiography image; (b) mask image; (c) vascular feature map; (d) multi-scale response
針對特征圖，區域增長算法從探測到的種子點開始，通過迭代計算和比較，當周圍像素點滿足條件f(p0)≥T時，不斷被加入目標區域的像素列表。期間，我們利用Otsu閾值算法[15]計算灰度直方圖的零階和一階累積增量，自適應地獲取優化的閾值T。以上過程運行的結果是：通常只能獲得精確的血管主分支，而血管的細節和末梢丟失了。作為補償，我們設計采用一階方向導數獲取血管像素點周圍的圖像強度變化，從而完成血管細節修復，方法是：在不同尺度(σ=1,5,...,9)下，梯度｜?Iσ｜的幅度代表了該點強度的斜度。多尺度下的結果計算如下：
| $\left\{ \begin{align} & \partial {{I}_{\sigma }}\left( x,y;\sigma \right)=I\left( x,y \right)\otimes \sigma *\partial G\left( x,y;\sigma \right) \\ & \left| \nabla {{I}_{\sigma }} \right|=\sqrt{{{\left( {{\partial }_{x}}{{I}_{\sigma }} \right)}^{2}}+{{\left( {{\partial }_{y}}{{I}_{\sigma }} \right)}^{2}}} \\ & g({{p}_{0}})=\underset{{{\sigma }_{min}}\le \sigma \le {{\sigma }_{max}}}{\mathop{\max }}\,\left( \left| \nabla {{I}_{\sigma }} \right| \right) \\ \end{align} \right.$ |
血管細節修復過程主要依據了g(p0),圖 2(d)顯示了造影圖像的梯度幅度和g(p0)。類似以前的過程，區域增長算法被用于進一步的分割；來自主分支分割后的所有的血管點被當作血管點。作為血管像素的8鄰域分類可采用條件f(p0)≥T-δv進行，其中δv是f(p0)的標準偏差。由于g(p0)提供了造影圖像中血管的邊緣信息，在區域增長算法作用下，條件f(p0)≥T-δv可以將主分支的區域增長結果擴展至血管邊緣和末梢。
3 特征測量和血管結構的模糊識別
在不同成像角度下，冠狀動脈空間結構的二維投影具有復雜的形態。因此局部血管的測量采用一個圓周探測器沿血管路徑掃描，并進行逐點測量與結構判別，探測器與血管局部骨架線之間的關系如圖 3所示。結構識別時利用多種局部特征數據進行測量，包括：圓周探測器與血管軸線的交點Pm,m=1,…,Mt (也稱為掃描點)，特征圖或交點Pm處的亮度特征、探測器中心Ot到點Pm的矢量VO,m，Gabor濾波器響應，血管方向場Vm。上述測量數據構成了血管的多特征測度[14]：令n=1,2,3,4分別代表血管的末梢、段、分支、交叉四種局部結構模式。區域增長過程中，沿血管骨架線的第t個掃描點處：當掃描點Pm的特征{Mt}、{}、{VO,m}被看作第n類模式的模糊子集時，各特征對應的模糊隸屬度可以寫為
| ${\mu _{M,n}} = \left\{ {\matrix{ {\exp \left[ {{{{{({M_t} - n)}^2}} \over { - \sigma _M^2}}} \right]} \hfill & {{M_t} \ge n} \hfill \cr 0 \hfill & {{M_t} ＜ n} \hfill \cr } } \right.,$ |
| ${{\mu }_{V,n}}\left( m \right)=\exp \left[ \frac{1-\left| \cos \right.{{\left( \angle {{V}_{O,m}}-\left. \angle {{V}_{m}} \right| \right)}^{2}}}{-\sigma _{v}^{2}} \right],$ |
| ${{\mu }_{\tilde{\psi },n}}\left( m \right)=\exp \left[ \frac{{{\left| {\tilde{\psi }} \right|}_{\max }}-{{\left| {{{\tilde{\psi }}}_{m}} \right|}^{2}}}{-\sigma _{\psi }^{2}} \right],$ |
其中σM=2,σv=σψ=1,。采用上述測度進行血管結構識別時，我們提出模糊血管判別方法。血管結構模式的判斷取決于掃描點Pm(m=1，…，Mt)上各類特征的聯合測度，基于多特征的模糊判別函數識別一點的真實結構時，首先可以假設式(3) 的計算結果的可信度受到特征強度的影響--較細和較暗的血管分支在視覺和診斷上意義不大；如果將特征強度作為決策的風險，則可采用可靠性系數有效控制權重，越大時，可靠性越大，風險越小。因此不難定義基于上述多特征模糊測度的血管結構判別公式如下：
圖3
血管的結構探測示意圖和仿真模型
(a) 示意圖，利用圓周探測器測量血管局部結構；(b)右冠狀動脈血管仿真模型
Figure3. Illustration of vascular structure detection and simulation model(a) an illustration of vessel structure inference using circle detector; (b) vascular simulation of right coronary artery
令Bn(n=1,…,4)對應于第n類模式的模糊子集，且其作用于Ω上：Ω={Mt,,,∠VO,m｜m∈[1,Mt],Zm∈[0,1],∈[0,1],∠VO,m∈[0,2π],t∈?T}。那么，探測到的第t個狀態集{Xt,Pm｜m=1,…,Mt}屬于第n類模式的隸屬度為
| $\begin{align} & {{\mu }_{B,m}}={{w}_{1}}\centerdot {{\mu }_{M,n}}+{{w}_{2}}\sum\limits_{m=1:{{M}_{t}}}{{{\alpha }_{m}}\centerdot {{\mu }_{M,n}}}\left( m \right) \\ & +{{w}_{3}}\centerdot \sum\limits_{m=1:{{M}_{t}}}{{{\alpha }_{m}}\centerdot {{\mu }_{\tilde{\psi },n}}\left( m \right)}={{w}_{1}}\centerdot {{\mu }_{M,n}} \\ & +\sum\limits_{m=1:{{M}_{t}}}{{{\alpha }_{m}}\centerdot }\left[ {{w}_{2}}\centerdot {{\mu }_{V,n}}\left( m \right)+{{w}_{3}}\centerdot {{\mu }_{\tilde{\psi },n}}\left( m \right) \right], \\ \end{align}$ |
上述加權和的構造形式考慮到了各特征的差異，加權系數滿足；可靠性系數滿足參數通過血管特征強弱控制隸屬度的計算。沿血管中心線探測血管結構過程中，不斷產生多特征樣本模糊集。依據隸屬度的最大貼近度原理，樣本集的分類識別定義為：在相同論域Ω下，令Bi(i=1,…,4)代表第i類血管結構的模糊子集，如果模糊集合A可具體化為狀態序列的當前模式，當滿足，則模糊集合A最接近模式Bm。
在背景灰度為零的圖像中，仿真模型的構造利用了變尺度和強度的二維高斯函數沿特定的血管骨架線逐點卷積生成，結果如圖 3(b)所示。仿真血管結構識別實驗中，可靠性系數一般選擇αm=1/Mt。仿真模型中標注有感興趣血管點，矩形窗口中的結構判別結果比較在表 1中。表 1統計了Xt對應各點處的實際血管模式，除了空間形態特異性點N4之外，計算結果Bm與真實血管模式一致。
表1
血管仿真模型中感興趣點的偽模式、模糊隸屬度、與結構判別結果的比較
Table1.
Comparison of pseudo patterns,degree of fuzzy memberships,inferred results of the interest sites in the vascular simulation structure
4 實驗結果與討論
實驗采用了10套DSA冠脈造影序列，圖像數據來自解放軍第四五八醫院放射科和北京朝陽醫院心臟中心。實驗分為兩部分，一是圖像預處理、區域增長和細節修復；二是血管結構識別。實驗中，血管分割部分采用Visual C++編程；結構識別部分采用Matlab7.0編程，電腦配置為：PC Intel(R) Core(TM) i5 CPU3.2GHz,8G-DDR-RAM。
4.1 血管圖像分割
實驗基于本文提出的血管函數和基于單純閾值[11]下的區域增長算法，并采用和Frangi[12]的算法對比的方式進行驗證。結果表明我們所提方法的分割結果優于Frangi函數作用下的分割結果，以上血管分割實驗如圖 4所示。分析其主要原因發現：Frangi的血管函數作用于造影圖像后，獲得的血管特征圖在大血管和主血管分支區域具有較高的響應，而在較細的血管和末梢便喪失了血管特征，此外血管邊界處的效應往往較弱，很難完整和真實地判斷血管的結構形態。
圖4
血管增強和區域增長效果比較
(a1)和(b1)為兩例原始圖像; (a2)和(b2)分別為 Frangi血管函數處理后的相應效果;(a3)和(b3)分別為對應于(a2)和(b2)的單純閾值血管區域增長結果； (a4)和(b4)分別為采用本文血管函數的增強效果；(a5)和(b5)分別為新的血管函數作用下的單純閾值區域增長結果
Figure4. Comparison of the effects of vessel enhancement and region growing(a1),(b1) were two different original images; (a2),(b2) were the corresponding results by Frangi’s vessel function,based on which (a3),(b3) were the corresponding results with region growing algorithm,respectively; (a4),(b4) were the results of the proposed vessel function on which (a5),(b5) were the corresponding results of region growing algorithm,respectively
血管細節修復過程如圖 5所示，(a1)、(b1)、(c1)為3例原始血管造影圖像；(a2)、(b2)、(c2)為血管主分支分割結果;(a3)、(b3)、(c3)為修復效果。比較圖中橢圓強調的區域可見，血管修復后效果顯著提高：位置標識1處原始前處理遺留的孔洞被修復填充；且在標識2和3處，小血管得以增長。實驗采用了從不同投影方向下獲取的X線血管造影圖像進行嘗試，都獲得了類似的良好效果。
圖5
三幅造影圖像的血管分割結果比較
(a1)、(b1)、(c1)：3例原始造影圖像；(a2)、(b2)、(c2)：利用區域增長算法得到的相應血管的主分支；(a3)、(b3)、(c3)對應的血管細節修復結果(白色圈為修復后增加的血管細節)
Figure5. Comparisons of the vessel segmentation results for three angiography images(a1),(b1),(c1): three original angiography images; (a2),(b2),(c2): the corresponding segmentation results of vessel main branch using the region-growing algorithm; (a3),(b3),(c3): the corresponding repaired results
以上結果比較明顯地證明了本方法的優越性，即新的血管函數能夠從圖像中對比度高和對比度低的區域自適應地增強血管目標；在此基礎上實施區域增長和細節修復后能夠較大程度地提取血管網絡。統計表明：本方法提取的血管面積較單獨使用區域增長算法后平均增加了5%左右，能夠從低灰度區域中恢復出一些較細的血管段。與文獻[14]的血管中心線提取結果相比，我們的方法提取的血管長度為2 523個血管點，文獻[14]的結果為1 477個血管點，平均增加了1 046個血管點的長度。
4.2 血管結構識別
針對血管結構的判別，實驗采用了1個仿真圖像和10套實際造影圖像，手段與評價方法如下：① 順序標記血管段并進行結構探測：將預處理得到的血管樹鏈碼按由低到高的級別分解為多條血管段，并用鏈碼的形式存儲每一條血管段。② 按照血管樹鏈碼順序，圓周探測器逐點、順序探測各種血管特征測度，并識別血管結構。③ 由心內科醫生逐點評價識別的結果。實驗內容包括：基于模糊識別算法的仿真血管和實際造影圖像識別和結構驗證。結構識別前的預處理是：將實際造影圖像轉變為規定尺寸；調整圖像中血管目標灰度為高亮值，背景灰度為低暗值；計算多尺度血管特征圖、方向場V,并對特征圖二值化和細化，獲取初始骨架CE，上述預處理過程的時耗為30 s。代表性效果顯示如圖 6所示。
圖6
利用區域增長算法分割血管、提取骨架線、血管結構自動識別結果
(a)原始造影圖像；(b)區域增長結果；(c)疊加在原始圖像上的血管中心線
Figure6. Experimental results respectively from algorithms of region growing,extracting skeleton line,and vascular structure recognition(a) the original angiography image; (b) region growing; (c) vascular skeleton lines overlapping in the original image
針對實際造影圖像中的冠狀動脈樹結構分析，表 2統計了10套臨床DSA造影圖像中不同級別血管段的識別效率。對應于各級血管段，在原始噪聲條件下，模糊識別算法對冠脈樹的第1、2級血管段逐點的識別正確率平均達到100%；對第3~5級血管段逐點進行識別的平均正確率均高于89%；總體正確率達到90.59%，相比文獻[14]中所述識別算法的平均率增加了1.5%。本文所提出的方法對于冠脈樹三維重建的價值主要反映在識別算法對特異性結構(如血管分支點)和非特異性結構(如交叉點)的分辨力和判別效率，對此我們統計了10套數據的結構分類識別情況(見表 3)：表中統計出左、右冠脈平均分支數分別為10.2與8.2，數量誤差小于1，并且符合醫生臨床測量的實際情況。
表2
10套造影數據和仿真模型的血管結構識別的效果比較(%)
Table2.
Comparison of effects of vascular structure recognition in 10 sets data and vessel phantom (%)
表3
血管結構識別分類比較
Table3.
Classified comparison of vascular structure recognition
5 結論
本文提出了一種全自動的血管多尺度濾波和區域增長方法,在此基礎上實現了基于模糊多特征測度的血管結構判別算法，以此分割X線造影圖像中的血管并取得良好效果。血管預處理中采用的多尺度函數具有更少的參數以便自動提取血管，在血管修復過程中，采用了特定掩膜和種子點自動探測技術，以此分割出血管邊緣、末梢，及更細的血管網絡。模糊多特征血管結構判別算法采用了血管骨架線和特征圖的Gabor響應算子，在仿真圖像和實際圖像中獲得良好的實驗效果。統計表明：本方法提取的血管面積較單獨使用區域增長算法后平均增加了5%左右；提取的血管長度為2 523，比文獻[14]增加了1 046個血管點的長度。對第2級以上血管段逐點進行識別的總體正確率平均達到90.59%，比文獻[14]所述識別算法平均增加了1.59%。總之本文提出的方法在血管造影圖像分析和血管提取、結構判別方面有較好的研究意義，未來將結合血管三維重建、輔助診斷、介入手術仿真等應用領域不斷研究探索。
引言
分割并描述X線造影圖像中的冠狀動脈樹，對血管的三維重建、測量和診斷具有重要意義。自20世紀90年代以來，X線在醫學領域受到廣泛關注。X線冠脈造影技術被認為是冠脈疾病診斷的“金標準”。X線造影成像是利用空間X線錐形束穿過目標區域后得到投影圖像，受限于成像條件和環境，造影圖像往往受到呼吸、心臟搏動造成的偽影和各種噪聲(如環境噪聲和人體不同組織對X線的不同衰減作用)干擾。
造影圖像中除了能夠清晰觀察的血管之外，還存在骨骼、脊柱、軟組織等，這些結構形態以及造影劑的非均勻顯像效果使得血管的自動檢測和計算機輔助介入過程比較困難。因此，有必要利用圖像處理技術獲取精確的血管結構并提取血管。血管造影圖像的分割[1-2]涉及血管圖像預處理、動脈樹軸線估計、結構和寬度測量等內容，其中預處理用來增強血管目標、降低背景和噪聲的影響。各種血管分割技術中[1-2]，提取血管特征并利用血管先驗知識構造算法模型能夠有效地分割出造影圖像中的血管樹。然而，由于造影劑的分布不均和成像過程中的環境噪聲以及復雜的血管形態(尤其是血管節點的多樣性和大曲率的血管段)，血管樹的分割具有極大的挑戰性。
近來有很多血管造影圖像分割方面的報道，經典的算法包括直接跟蹤法[3-4]、動態輪廓法[5-6]、中心線優化法[7-8]、Markov標記點過程[9]等。文獻[1-2]總體回顧和評論了上述方法。基于血管跟蹤的方法進行血管分割時，通常采用一個局部算子作用于血管樹的某個初始點并開始跟蹤，執行過程中周而復始地進行測量與決策。區域增長算法用于血管分割的基本思想是從一些血管區域內的種子點和子區域開始，按照某種排除準則，不斷地接受新的血管像素對當前區域進行擴充。通常有兩種分割準則：數值相似性和空間相近[10]。采用基于數值相似性的區域增長算法時，簡單的方法是利用一個單純的閾值法則[11]，即：假使像素具有相近的強度時很可能屬于同一目標，優化的閾值選擇需經過試驗確認。該方法的缺點是易受噪聲和對比劑不均勻性的干擾，常常產生過分割或者孔洞效果。空間一致性條件下的波擴張技術[10-11]是基于空間相似性的典型區域增長方法，也可以被當作一種順序區域增長策略，其中候選像素的排除與否依據了增長波幾何校正的局部策略與擬合準則。空間一致性擴張使得拓撲分析更加容易。例如，擴張界面的連續性可以被用來探測血管分支。此外，在X線造影圖像中識別冠狀動脈樹結構，對于定性和定量描述CAT結構、三維幾何重建、計算機輔助血管介入治療也具有重要作用。本文的血管結構判別算子作用于區域增長過程中順序獲得的骨架線，并通過特征測量進行結構的判別；測量中利用了血管預處理、特征監測、先驗知識等手段進行識別。首先，利用Hessian特征值分析算法[12]和方向濾波器[13-14]技術，經過多尺度血管特征提取和增強后獲得血管圖、矢量場。其次，結合Gabor能量描述子[14]，血管的形狀和方向特征能夠被準確地測量和描述。
本文的主要目的在于：① 建立在多尺度濾波之上的新血管函數，為區域增長提供更好的特征圖；② 提出基于細節修復的血管分割方法；③ 提出基于特征測量的血管結構識別算法。
1 造影圖像預處理
1.1 多尺度血管增強
我們可以將血管看做彎曲的帶狀目標，其直徑沿彎曲的主軸漸進變化。由于血管直徑分布在有限的范圍，我們可以采用尺度σ相當于血管半徑的一系列高斯核函數G(x,y,σ)卷積原始圖像IO(x,y)從而得到不同尺度下的血管響應：
| $\begin{align} & {{I}_{\sigma }}\left( x,y;\sigma \right)={{I}_{O}}\left( x,y \right)\otimes G\left( x,y;\sigma \right), \\ & 其中G\left( x,y;\sigma \right)=\frac{1}{\sqrt{2\pi \sigma }}{{e}^{\frac{-{{x}^{2}}+{{y}^{2}}}{2{{\sigma }^{2}}}}} \\ \end{align}$ |
利用上式對圖像卷積后，寬度相當于σ的血管結構獲得較大增強，反之則被抑制。圖像點p0的局部特性可以由該點處的Hessian矩陣分析得到，Hessian矩陣定義為
| $H({{p}_{0}},\sigma )={{\sigma }^{2}}*\left\{ -{{\partial }^{2}}G\left( x,y;\sigma \right) \right\}\otimes {{I}_{O}},$ |
其中σ2用來歸一化微分結果。當血管相對背景為低暗目標時，對應于血管點處的Hessian矩陣具有較大的負特征值λ2和一個絕對值較小的特征值λ1，且滿足關系。Frangi等[12]曾提出一種包含三個參數的血管測度函數，然而這些參數值不易確定。我們對比了當前各種血管函數[12-14],并提出新的血管函數：
| $\begin{align} & f({{p}_{0}})=\underset{{{\sigma }_{min}}\le \sigma \le {{\sigma }_{max}}}{\mathop{\max }}\,f({{p}_{0}},\sigma )/exp\sqrt{2\sigma },其中f({{p}_{0}},\sigma ) \\ & \left\{ \begin{array}{*{35}{l}} {{\lambda }_{2}} & {{\lambda }_{2}}>{{\sigma }^{2}}/4 \\ 0 & 其它 \\ \end{array} \right. \\ \end{align}$ |
利用公式(3) 的血管函數作用于整幅血管造影圖像后，不妨稱為“血管特征圖”。我們利用尺度參數σ2/4作為閾值。σmin和σmax為最小和最大尺度，其范圍的選擇可依據圖像中血管的最小和最大直徑確定。f(p0)為血管特征圖。圖 1中描述了新血管函數和傳統的Frangi函數增強血管的效果對比。明顯可見，新函數對較粗的和較細的血管具有相似的增強效果。該血管函數利用了特征圖閾值，能夠為后續的區域增算法提供自適應的初始結果。
圖1
不同血管函數對造影圖像增強效果比較
(a)原始造影圖像；(b)Frangi的血管函數增強效果；(c)利用新血管函數得到的血管特征圖
Figure1. Comparison of the angiography image enhancement based on different vessel functions(a) original angiography image; (b) result from Frangi’s vessel function; (c) result from the proposed vessel function
1.2 Gabor濾波器響應
在血管特征測量與結構識別算法設計中，利用Gabor濾波器[14]作用于血管軸線上的各個觀察點，濾波器響應最大的方向便為局部血管方向，從而有效避免了上述問題對結構判別的影響。Gabor濾波器函數為
| $\psi ({{X}_{0}})=\frac{\sum{{{W}_{i,j}}\left( \omega ,\theta ,S \right)\centerdot {{I}_{＜{{X}_{0}}>s}}\left( i,j \right)}}{\sum{{{W}_{i,j}}\left( \omega ,\theta ,S \right)}},$ |
其中Wi,j(ω,θ,S)為Gabor掩膜函數，具有可變的頻率帶寬ω、方向θ和尺度S=M×N，且有-ZX0≤Wi,j(ω,θ,S)≤ZX0。因此，式(4) 滿足：-1＜ψ(X0)＜1。掩膜與圖像進行相關時，每一點對應尺寸的區域中血管的特性如果和掩膜相近，則有最大響應輸出。令，式(4) 可歸一化為0~1之間，最大響應輸出表示為：，其中θ取血管方向場中點的X0對應值V2(X0)時，可極大地減少計算時間。實際處理中，作為預處理手段，血管相似性函數作用于整幅圖像分析血管的亮度特征和方向矢量；而Gabor響應函數作用于血管軸線上的各個觀察點，并測量血管的局部強度和方向響應。
2 區域增長和細節修復
為了分割出血管樹，我們建立了一種全自動的區域增長方法。主要過程為：自動種子點選擇、主分支分割、血管細節修復等。在進行種子點的選擇時，理論上可以自動搜索血管特征圖中的最大強度點，該點自然位于血管區域內。然而對實際造影圖像進行處理時，有時會遇到X線造影圖像存在邊界偽影(見圖 2a)，致使最大特征強度點位于造影圖像中血管區域之外的邊界，產生錯誤結果。因此我們通常利用圖 2(b)所示的掩膜保證自動搜索的正確性。圖 2(c)為特征圖和種子點自動選擇的位置(種子點位于下方的圓圈中心處)。
圖2
種子點的選擇和多尺度血管濾波結果
(a)原始造影圖像；(b)掩膜；(c)血管特征圖；(d)多尺度響應
(a) original angiography image; (b) mask image; (c) vascular feature map; (d) multi-scale response
針對特征圖，區域增長算法從探測到的種子點開始，通過迭代計算和比較，當周圍像素點滿足條件f(p0)≥T時，不斷被加入目標區域的像素列表。期間，我們利用Otsu閾值算法[15]計算灰度直方圖的零階和一階累積增量，自適應地獲取優化的閾值T。以上過程運行的結果是：通常只能獲得精確的血管主分支，而血管的細節和末梢丟失了。作為補償，我們設計采用一階方向導數獲取血管像素點周圍的圖像強度變化，從而完成血管細節修復，方法是：在不同尺度(σ=1,5,...,9)下，梯度｜?Iσ｜的幅度代表了該點強度的斜度。多尺度下的結果計算如下：
| $\left\{ \begin{align} & \partial {{I}_{\sigma }}\left( x,y;\sigma \right)=I\left( x,y \right)\otimes \sigma *\partial G\left( x,y;\sigma \right) \\ & \left| \nabla {{I}_{\sigma }} \right|=\sqrt{{{\left( {{\partial }_{x}}{{I}_{\sigma }} \right)}^{2}}+{{\left( {{\partial }_{y}}{{I}_{\sigma }} \right)}^{2}}} \\ & g({{p}_{0}})=\underset{{{\sigma }_{min}}\le \sigma \le {{\sigma }_{max}}}{\mathop{\max }}\,\left( \left| \nabla {{I}_{\sigma }} \right| \right) \\ \end{align} \right.$ |
血管細節修復過程主要依據了g(p0),圖 2(d)顯示了造影圖像的梯度幅度和g(p0)。類似以前的過程，區域增長算法被用于進一步的分割；來自主分支分割后的所有的血管點被當作血管點。作為血管像素的8鄰域分類可采用條件f(p0)≥T-δv進行，其中δv是f(p0)的標準偏差。由于g(p0)提供了造影圖像中血管的邊緣信息，在區域增長算法作用下，條件f(p0)≥T-δv可以將主分支的區域增長結果擴展至血管邊緣和末梢。
3 特征測量和血管結構的模糊識別
在不同成像角度下，冠狀動脈空間結構的二維投影具有復雜的形態。因此局部血管的測量采用一個圓周探測器沿血管路徑掃描，并進行逐點測量與結構判別，探測器與血管局部骨架線之間的關系如圖 3所示。結構識別時利用多種局部特征數據進行測量，包括：圓周探測器與血管軸線的交點Pm,m=1,…,Mt (也稱為掃描點)，特征圖或交點Pm處的亮度特征、探測器中心Ot到點Pm的矢量VO,m，Gabor濾波器響應，血管方向場Vm。上述測量數據構成了血管的多特征測度[14]：令n=1,2,3,4分別代表血管的末梢、段、分支、交叉四種局部結構模式。區域增長過程中，沿血管骨架線的第t個掃描點處：當掃描點Pm的特征{Mt}、{}、{VO,m}被看作第n類模式的模糊子集時，各特征對應的模糊隸屬度可以寫為
| ${\mu _{M,n}} = \left\{ {\matrix{ {\exp \left[ {{{{{({M_t} - n)}^2}} \over { - \sigma _M^2}}} \right]} \hfill & {{M_t} \ge n} \hfill \cr 0 \hfill & {{M_t} ＜ n} \hfill \cr } } \right.,$ |
| ${{\mu }_{V,n}}\left( m \right)=\exp \left[ \frac{1-\left| \cos \right.{{\left( \angle {{V}_{O,m}}-\left. \angle {{V}_{m}} \right| \right)}^{2}}}{-\sigma _{v}^{2}} \right],$ |
| ${{\mu }_{\tilde{\psi },n}}\left( m \right)=\exp \left[ \frac{{{\left| {\tilde{\psi }} \right|}_{\max }}-{{\left| {{{\tilde{\psi }}}_{m}} \right|}^{2}}}{-\sigma _{\psi }^{2}} \right],$ |
其中σM=2,σv=σψ=1,。采用上述測度進行血管結構識別時，我們提出模糊血管判別方法。血管結構模式的判斷取決于掃描點Pm(m=1，…，Mt)上各類特征的聯合測度，基于多特征的模糊判別函數識別一點的真實結構時，首先可以假設式(3) 的計算結果的可信度受到特征強度的影響--較細和較暗的血管分支在視覺和診斷上意義不大；如果將特征強度作為決策的風險，則可采用可靠性系數有效控制權重，越大時，可靠性越大，風險越小。因此不難定義基于上述多特征模糊測度的血管結構判別公式如下：
圖3
血管的結構探測示意圖和仿真模型
(a) 示意圖，利用圓周探測器測量血管局部結構；(b)右冠狀動脈血管仿真模型
Figure3. Illustration of vascular structure detection and simulation model(a) an illustration of vessel structure inference using circle detector; (b) vascular simulation of right coronary artery
令Bn(n=1,…,4)對應于第n類模式的模糊子集，且其作用于Ω上：Ω={Mt,,,∠VO,m｜m∈[1,Mt],Zm∈[0,1],∈[0,1],∠VO,m∈[0,2π],t∈?T}。那么，探測到的第t個狀態集{Xt,Pm｜m=1,…,Mt}屬于第n類模式的隸屬度為
| $\begin{align} & {{\mu }_{B,m}}={{w}_{1}}\centerdot {{\mu }_{M,n}}+{{w}_{2}}\sum\limits_{m=1:{{M}_{t}}}{{{\alpha }_{m}}\centerdot {{\mu }_{M,n}}}\left( m \right) \\ & +{{w}_{3}}\centerdot \sum\limits_{m=1:{{M}_{t}}}{{{\alpha }_{m}}\centerdot {{\mu }_{\tilde{\psi },n}}\left( m \right)}={{w}_{1}}\centerdot {{\mu }_{M,n}} \\ & +\sum\limits_{m=1:{{M}_{t}}}{{{\alpha }_{m}}\centerdot }\left[ {{w}_{2}}\centerdot {{\mu }_{V,n}}\left( m \right)+{{w}_{3}}\centerdot {{\mu }_{\tilde{\psi },n}}\left( m \right) \right], \\ \end{align}$ |
上述加權和的構造形式考慮到了各特征的差異，加權系數滿足；可靠性系數滿足參數通過血管特征強弱控制隸屬度的計算。沿血管中心線探測血管結構過程中，不斷產生多特征樣本模糊集。依據隸屬度的最大貼近度原理，樣本集的分類識別定義為：在相同論域Ω下，令Bi(i=1,…,4)代表第i類血管結構的模糊子集，如果模糊集合A可具體化為狀態序列的當前模式，當滿足，則模糊集合A最接近模式Bm。
在背景灰度為零的圖像中，仿真模型的構造利用了變尺度和強度的二維高斯函數沿特定的血管骨架線逐點卷積生成，結果如圖 3(b)所示。仿真血管結構識別實驗中，可靠性系數一般選擇αm=1/Mt。仿真模型中標注有感興趣血管點，矩形窗口中的結構判別結果比較在表 1中。表 1統計了Xt對應各點處的實際血管模式，除了空間形態特異性點N4之外，計算結果Bm與真實血管模式一致。
表1
血管仿真模型中感興趣點的偽模式、模糊隸屬度、與結構判別結果的比較
Table1.
Comparison of pseudo patterns,degree of fuzzy memberships,inferred results of the interest sites in the vascular simulation structure
4 實驗結果與討論
實驗采用了10套DSA冠脈造影序列，圖像數據來自解放軍第四五八醫院放射科和北京朝陽醫院心臟中心。實驗分為兩部分，一是圖像預處理、區域增長和細節修復；二是血管結構識別。實驗中，血管分割部分采用Visual C++編程；結構識別部分采用Matlab7.0編程，電腦配置為：PC Intel(R) Core(TM) i5 CPU3.2GHz,8G-DDR-RAM。
4.1 血管圖像分割
實驗基于本文提出的血管函數和基于單純閾值[11]下的區域增長算法，并采用和Frangi[12]的算法對比的方式進行驗證。結果表明我們所提方法的分割結果優于Frangi函數作用下的分割結果，以上血管分割實驗如圖 4所示。分析其主要原因發現：Frangi的血管函數作用于造影圖像后，獲得的血管特征圖在大血管和主血管分支區域具有較高的響應，而在較細的血管和末梢便喪失了血管特征，此外血管邊界處的效應往往較弱，很難完整和真實地判斷血管的結構形態。
圖4
血管增強和區域增長效果比較
(a1)和(b1)為兩例原始圖像; (a2)和(b2)分別為 Frangi血管函數處理后的相應效果;(a3)和(b3)分別為對應于(a2)和(b2)的單純閾值血管區域增長結果； (a4)和(b4)分別為采用本文血管函數的增強效果；(a5)和(b5)分別為新的血管函數作用下的單純閾值區域增長結果
Figure4. Comparison of the effects of vessel enhancement and region growing(a1),(b1) were two different original images; (a2),(b2) were the corresponding results by Frangi’s vessel function,based on which (a3),(b3) were the corresponding results with region growing algorithm,respectively; (a4),(b4) were the results of the proposed vessel function on which (a5),(b5) were the corresponding results of region growing algorithm,respectively
血管細節修復過程如圖 5所示，(a1)、(b1)、(c1)為3例原始血管造影圖像；(a2)、(b2)、(c2)為血管主分支分割結果;(a3)、(b3)、(c3)為修復效果。比較圖中橢圓強調的區域可見，血管修復后效果顯著提高：位置標識1處原始前處理遺留的孔洞被修復填充；且在標識2和3處，小血管得以增長。實驗采用了從不同投影方向下獲取的X線血管造影圖像進行嘗試，都獲得了類似的良好效果。
圖5
三幅造影圖像的血管分割結果比較
(a1)、(b1)、(c1)：3例原始造影圖像；(a2)、(b2)、(c2)：利用區域增長算法得到的相應血管的主分支；(a3)、(b3)、(c3)對應的血管細節修復結果(白色圈為修復后增加的血管細節)
Figure5. Comparisons of the vessel segmentation results for three angiography images(a1),(b1),(c1): three original angiography images; (a2),(b2),(c2): the corresponding segmentation results of vessel main branch using the region-growing algorithm; (a3),(b3),(c3): the corresponding repaired results
以上結果比較明顯地證明了本方法的優越性，即新的血管函數能夠從圖像中對比度高和對比度低的區域自適應地增強血管目標；在此基礎上實施區域增長和細節修復后能夠較大程度地提取血管網絡。統計表明：本方法提取的血管面積較單獨使用區域增長算法后平均增加了5%左右，能夠從低灰度區域中恢復出一些較細的血管段。與文獻[14]的血管中心線提取結果相比，我們的方法提取的血管長度為2 523個血管點，文獻[14]的結果為1 477個血管點，平均增加了1 046個血管點的長度。
4.2 血管結構識別
針對血管結構的判別，實驗采用了1個仿真圖像和10套實際造影圖像，手段與評價方法如下：① 順序標記血管段并進行結構探測：將預處理得到的血管樹鏈碼按由低到高的級別分解為多條血管段，并用鏈碼的形式存儲每一條血管段。② 按照血管樹鏈碼順序，圓周探測器逐點、順序探測各種血管特征測度，并識別血管結構。③ 由心內科醫生逐點評價識別的結果。實驗內容包括：基于模糊識別算法的仿真血管和實際造影圖像識別和結構驗證。結構識別前的預處理是：將實際造影圖像轉變為規定尺寸；調整圖像中血管目標灰度為高亮值，背景灰度為低暗值；計算多尺度血管特征圖、方向場V,并對特征圖二值化和細化，獲取初始骨架CE，上述預處理過程的時耗為30 s。代表性效果顯示如圖 6所示。
圖6
利用區域增長算法分割血管、提取骨架線、血管結構自動識別結果
(a)原始造影圖像；(b)區域增長結果；(c)疊加在原始圖像上的血管中心線
Figure6. Experimental results respectively from algorithms of region growing,extracting skeleton line,and vascular structure recognition(a) the original angiography image; (b) region growing; (c) vascular skeleton lines overlapping in the original image
針對實際造影圖像中的冠狀動脈樹結構分析，表 2統計了10套臨床DSA造影圖像中不同級別血管段的識別效率。對應于各級血管段，在原始噪聲條件下，模糊識別算法對冠脈樹的第1、2級血管段逐點的識別正確率平均達到100%；對第3~5級血管段逐點進行識別的平均正確率均高于89%；總體正確率達到90.59%，相比文獻[14]中所述識別算法的平均率增加了1.5%。本文所提出的方法對于冠脈樹三維重建的價值主要反映在識別算法對特異性結構(如血管分支點)和非特異性結構(如交叉點)的分辨力和判別效率，對此我們統計了10套數據的結構分類識別情況(見表 3)：表中統計出左、右冠脈平均分支數分別為10.2與8.2，數量誤差小于1，并且符合醫生臨床測量的實際情況。
表2
10套造影數據和仿真模型的血管結構識別的效果比較(%)
Table2.
Comparison of effects of vascular structure recognition in 10 sets data and vessel phantom (%)
表3
血管結構識別分類比較
Table3.
Classified comparison of vascular structure recognition
5 結論
本文提出了一種全自動的血管多尺度濾波和區域增長方法,在此基礎上實現了基于模糊多特征測度的血管結構判別算法，以此分割X線造影圖像中的血管并取得良好效果。血管預處理中采用的多尺度函數具有更少的參數以便自動提取血管，在血管修復過程中，采用了特定掩膜和種子點自動探測技術，以此分割出血管邊緣、末梢，及更細的血管網絡。模糊多特征血管結構判別算法采用了血管骨架線和特征圖的Gabor響應算子，在仿真圖像和實際圖像中獲得良好的實驗效果。統計表明：本方法提取的血管面積較單獨使用區域增長算法后平均增加了5%左右；提取的血管長度為2 523，比文獻[14]增加了1 046個血管點的長度。對第2級以上血管段逐點進行識別的總體正確率平均達到90.59%，比文獻[14]所述識別算法平均增加了1.59%。總之本文提出的方法在血管造影圖像分析和血管提取、結構判別方面有較好的研究意義，未來將結合血管三維重建、輔助診斷、介入手術仿真等應用領域不斷研究探索。
