磁共振腦影像血管結構提取方法研究

薛方，許朝萍

（1.西安醫學院衛生管理學院，陜西西安 710021；2.西安交通大學機械工程學院，陜西 西安 710049）

腦血管疾病具有發病率高、致殘率高的特點，是全球關注的公共衛生問題[1]。腦卒中是60 歲及以上人群的第二位死亡原因[2]，是腦血管病中最常見者。腦動脈的側支循環與腦梗的發生、發展和預后密切相關。研究表明，腦動脈環又稱威廉環，是腦動脈最主要的側支循環，其結構異常與否可以作為缺血性腦中風的一個風險因素[3-6]。

臨床用于腦疾病檢查的技術包括磁共振成像（Magnetic Resonance Imaging,MRI）、斷層掃描（Computed Tomography,CT）等。磁共振成像相較于斷層掃描可以使患者避免X-射線的危害。同時，因為磁共振成像是通過被試細胞中氫元素在磁場中共振現象而成像，特別是顱腦及腎臟的血管成像通過流動增強效應的原理展開，因此不需要造影劑的注入，可有效避免患者出現過敏及不良反應[7]，且可多方位任意層成像，因此在腦血管疾病篩查和鑒別診斷中具有不可取代的優勢[8-9]。雖然，目前在腦血管病診療中磁共振成像不可或缺，已建立了一定的影像診療策略，但在側支循環血管病精準評估及預后預測等方面還需要進一步完善[10]。

文中從影像解剖的形態學出發，探討各個閾值下分割圖像質量的優劣，提出最佳閾值。同時，通過Hessian 矩陣增強圖像從而更好地保留管狀區域，繼而完成血管輪廓及骨架提取。血管骨架與原圖像疊加之后基本吻合，表明血管骨架的提取效果良好，比較符合實際的血管走向，且提取的輪廓清晰，消除了背景干擾。

1 方法與資料

如圖1 所示，在圖像預處理中，通過圖像分割、消除背景區域等處理去除干擾因素，強化血管結構。圖像分割過程采用基于優化參數的灰度多閾值圖像分割算法，圖像增強過程運用多尺度Hessian 矩陣，通過不斷調整迭代尺度因子尋找最優增強結果。然后通過腐蝕、消除較小連通區域等一系列形態學處理提取血管骨架和輪廓。

圖1 血管結構提取步驟

文中原始數據采集于奧泰1.5T，8 陣列線圈超導磁共振掃描儀。磁共振腦血管圖像數據格式為標準DICOM 格式。數據采集參數：使用時間飛躍（Timeof-Flight,TOF）序列，重復時間（Repetition Time,TR）為21.2 ms，回波時間（Time-to-Echo，TE）為5.0 ms，翻轉角（Flip Angle）為25°，平均次數（Average，AV）為1，相位編碼方向視野為170 mm，頻率編碼方向視野為180 mm，總層數為120。

2 圖像分割

圖像分割的過程要求將目標圖像分割成若干個區域，并提取這些區域的特性。常用的圖像分割技術包括基于閾值的分割法、基于小波變換的分割法、基于區域的分割法以及基于神經網絡的分割法[11]。文中采用基于灰度的多閾值圖像分割算法。

經分析，圖像灰度值范圍為[0,255]，假設灰度值為k的像素點有nk個，總的像素數為：

則灰度值為k的像素點出現的頻率為：

設閾值的個數為m，則閾值集合以Td={td1,td2,…,tdm}表示，修正的閾值集合為Td={t0=-1,td1,td2,…,tdm,tdm+1=255}，當i＜j時，滿足tdi＜tdj。閾值的選擇使類間方差達到最大。

閾值選取完畢后,對像素值G(x,y)進行以下處理：

由式（3）可得，m個閾值可以將圖像分割為m+1個灰度等級。

如圖2 所示，將閾值個數不同的分割結果進行了對比。當m=2 時，圖像中左右的像素取值只有兩種情況，即0 和255。由于腦血管旁邊的組織灰度比較接近血管，故圖像被分割成了只有背景和組織的結果，不顯示血管。當m=3 時，圖像中的像素取值有三種情況，基本將除了腦血管之外的組織處理在同一個灰度水平，以方便濾除，但是血管的輪廓并不理想。當m=4 時，腦血管的輪廓信息比較全面，細小血管與較大血管的連接也比較完整，對組成動脈環腦血管輪廓的提取有優勢，更符合實際需要。當m＞5時，盡管腦血管輪廓信息被保留得很完整，但是在血管附近出現了很多細小的片狀信號，圖像灰度分布復雜，可能影響圖像的后續處理。因此，文中選取閾值個數m=4，批量處理MRI 圖像之后，繼續進行后續的數據圖像處理。

圖2 不同閾值分割結果

3 圖像增強

分割后的圖像還存在許多斑塊狀或片狀的區域，需要通過圖像增強使斑塊狀和片狀區域被濾除或得到減弱，從而更好地將管狀結構保留下來。多尺度濾波可以通過融合目標多個尺度細節特征的方式實現對目標區域的圖像增強[12-13]。

定義一幅二維圖像的Hessian 矩陣為：

式中，Gxx為二維圖像G(x,y)在X軸方向的二階偏導，Gyy為二維圖像G(x,y)在y軸方向的二階偏導，Gxy和Gyx是二維圖像G(x,y)在X、Y軸方向的混合偏導。其中，X軸方向的二階偏導為：

Y軸方向的二階偏導為：

X、Y軸方向的混合偏導為：

由式（7）可知，Hessian 矩陣是實對稱矩陣，具有兩個特征值和兩個特征向量。

在除圖像邊緣之外的每個像素點上，都對應一個Hessian 矩陣，它的特征值和特征向量可以被計算出來。Hessian 矩陣的特征值可以被用來判斷像素點是否為血管輪廓邊界上的點。模最大的特征向量與血管方向是垂直的，因此與這個特征向量垂直的方向就是血管的走向。血管的方向和輪廓邊界可以用Hessian 矩陣的特征向量和特征值來表示。

多尺度濾波器結合高斯函數進行構造，高斯函數的標準偏差σ就是迭代尺度因子。采用多尺度Hessian 矩陣對血管增強，通過不斷調整迭代尺度因子σ得到基于不同尺度的血管增強的最大輸出響應，迭代尺度因子的變化范圍為[σmin,σmax]，迭代步長增量為step。迭代尺度因子越小，對細小血管直徑的增強效果越好；迭代尺度因子越大，對較大直徑血管的增強效果越好。通過固定且較小的步長可以比對不同迭代尺度下的圖像增強效果。

圖3 所示為不同尺度因子范圍同一步長下的腦血管增強對比圖，每次迭代步長為1。其中，圖3（a）中組成動脈環的血管信號并未增強；圖3（b）中組成動脈環的血管信號雖然增強，但是中間有中空的地方；圖3（c）中組成動脈環的血管信號保存完好；圖3（d）中組成動脈環的血管信號雖然保存完好，但是血管也變粗，適合血管骨架提取，但不適合血管輪廓提取；圖3（e）的結果與圖3（d）類似，但是對于細小血管，信息丟失嚴重；圖3（f）血管的形貌已嚴重失真。故迭代尺度因子范圍σ∈[1,3]效果最好。

圖3 不同尺度因子的增強效果對比圖

4 骨架及輪廓提取

腦血管骨架提取的實現主要分為腐蝕、消除面積較小的連通區域、提取骨架、消除血管垂直方向的干擾幾個步驟。腐蝕是對圖像的形態學處理，在數學形態學運算中的作用是消除物體的邊界點，腐蝕半徑的選擇對腐蝕的效果至關重要[14-16]。圖4 所示為不同腐蝕半徑處理結果對比圖，圖4（a）的腐蝕半徑R=3，處理后圖中的腦血管斷點較多；圖4（b）的腐蝕半徑R=2，處理后圖中的腦血管仍存在少數斷點；圖4（c）的腐蝕半徑R=1，效果最優。

圖4 不同腐蝕半徑處理結果對比圖

消除面積較小的連通區域就是對于一塊連通的取值相同的像素區域，如果像素個數值小于設定值，則該區域被去除，以處理圖像中散在的斑塊，處理之后的圖像邊緣清晰，適合用于輪廓提取。圖5 所示為不同設定參數值消除面積較小連通區域的結果對比圖。圖5（a）為腐蝕后未進行該操作的圖像，有些地方未連接。由圖5 可知，低于S=6 700 的參數值設置對該實驗中消除面積較小連通區域的操作影響不大，故實驗選取較小的設定值參數128。

圖5 消除面積較小連通區域結果對比圖

骨架提取即細化處理，在保留原圖像幾何形狀的前提下，盡量減少圖像所包含的信息量。在研究中，腦血管骨架提取實際上是將圖形由管狀血管到線狀進行的變換，骨架就是血管中心線，將提取的腦血管骨架與原圖進行疊加，以比較走向，檢測二者是否吻合，評估血管骨架提取效果。圖6 所示為腦血管骨架提取圖像，其中，圖6（a）為骨架提取圖像；圖6（b）為圖6（a）與提取前圖像疊加對比圖，可以看出，提取出的骨架處于血管中心線的位置，提取效果良好。

圖6 骨架提取結果圖

經過消除面積較小的連通區域處理后得到的腦血管圖像的腦血管輪廓清晰，并且沒有多余的散在斑塊。圖7 所示為腦血管輪廓提取圖像及其與骨架提取圖像疊加對比圖，其中，圖7（a）為最優參數腦動脈環輪廓提取圖像；圖7（b）為圖7（a）中輪廓提取圖像與骨架提取圖像疊加對比圖。可以看出，血管骨架基本處于輪廓的中心位置。

圖7 輪廓提取結果圖

5 實驗結果與分析

文中以某醫院影像科采集的MRI 圖像作為實驗樣本，采用Matlab R2014a 平臺編程實現文中算法，并將其與傳統方法進行比較，結果如圖8 所示。

圖8 實驗對比圖

其中，圖8（b）為傳統方法進行圖像分割、增強并提取骨架的結果，圖8（c）為文中方法提取骨架的結果，圖8（d）為傳統方法進行圖像分割、增強并提取血管輪廓的結果，圖8（e）為文中方法提取血管輪廓的結果。由圖像對比可知，文中方法提取的腦動脈環骨架及輪廓結構與原始圖像更吻合，對于信號較弱的連通區域也能正確分析并提取結構，故處理所得圖像質量更高，更符合臨床需求。

6 結束語

文中在磁共振腦影像分割、增強等預處理的基礎上，通過圖像腐蝕、消除較小連通區域等形態學處理，提取血管骨架及血管輪廓，并將提取的結果與原圖像疊加對比。經比較分析可得，該方法能夠較為精確有效地提取腦血管骨架和邊緣輪廓結構。研究結果可以為血管狹窄的定位、定量和分類提供基礎數據，為下一步基于貝葉斯學習的腦血管異診斷常提供基礎，實現自動、精準的腦血管疾病診斷。