王力，孟慶民，孫超，郭永新，姚吉進(jìn)，宗寶良，焦青

1 山東第一醫(yī)科大學(xué)（山東省醫(yī)學(xué)科學(xué)院）醫(yī)學(xué)信息工程學(xué)院，泰安市，271016

2 泰安市中心醫(yī)院 介入放射科，泰安市，271000

3 山東第一醫(yī)科大學(xué)（山東省醫(yī)學(xué)科學(xué)院）放射學(xué)院，泰安市，271016

0 引言

據(jù)《中國心血管病報告2017（概要）》顯示，我國心血管病患病率及死亡率仍處于上升階段，心血管病現(xiàn)患人數(shù)高達(dá)2.9億。心血管病死亡占城鄉(xiāng)居民總死亡原因的首位，農(nóng)村為45.01%，城市為42.61%[1]。心血管病嚴(yán)重影響了人類健康和生活質(zhì)量，給患者及其家庭乃至整個社會都造成沉重的負(fù)擔(dān)。

數(shù)字減影成像（Digital Subtraction Angiography,DSA）是心血管疾病診斷的金標(biāo)準(zhǔn)[2]，是目前臨床廣泛采用的診斷和治療心血管疾病的手段之一[3]。根據(jù)DSA圖像可以判斷冠狀動脈血管狹窄的程度，從而進(jìn)行支架置入，將支架放置前后冠狀動脈的狹窄程度進(jìn)行對比，可評價治療效果。DSA檢查還可清晰地顯示冠狀動脈的側(cè)支循環(huán)情況，對于冠狀動脈搭橋手術(shù)的效果評價具有重要價值[4]。

由于受到流速、壁面剪切力和血管壁壓力的局部血流動力學(xué)因素的影響，冠狀動脈病變位置通常會出現(xiàn)在血管分叉、彎曲以及血管的開口處[5-6]。冠狀動脈分叉病變在臨床工作中十分常見，約占目前所有進(jìn)行冠狀動脈介入治療量的15%～20%，此位置手術(shù)難度大，極具技術(shù)挑戰(zhàn)性[7-8]。有些患者的冠狀動脈血管會產(chǎn)生不同程度的變異，如冠狀動脈血管分支較多，在常規(guī)體位下獲取的冠狀動脈血管圖像的分支，會產(chǎn)生相互掩蓋的現(xiàn)象，此時醫(yī)生難以判斷此處血管是否發(fā)生了病變。如果能得知分叉血管之間的角度，則可向DSA操作人員提供體位設(shè)置參考，因此確定冠狀動脈分叉血管之間的角度對醫(yī)生的決策具有重要指導(dǎo)意義。

該研究基于臨床實踐需要，在Windows7環(huán)境下，以Matlab為平臺，首先對DSA圖像進(jìn)行平滑處理，然后提取冠狀動脈血管的中心線及分叉點，從而實現(xiàn)DSA圖像冠狀動脈血管角度的測量。

1 提取冠狀動脈血管的中心線

1.1 DSA圖像平滑

DSA圖像在采集、量化與傳輸過程中會由于一些干擾，致使圖像包含很多離散隨機(jī)的噪聲，這種情況下通常采用平滑的方法對圖像進(jìn)行處理以消除噪聲。在Matlab平臺中讀入DSA圖像（圖1(a)），使用rgb2gray函數(shù)將圖像轉(zhuǎn)化為8位灰度BMP圖像，然后使用fspecial函數(shù)與imfilter函數(shù)進(jìn)行高斯濾波。將fspecial函數(shù)中的HSIZE參數(shù)設(shè)置成較大數(shù)值（1 000），可獲得圖像的背景（圖1(b)），將HSIZE參數(shù)設(shè)置成較小數(shù)值（10），可得到消除圖像的噪聲（圖1(c)）。高斯濾波是對整幅圖像進(jìn)行加權(quán)平均的過程，每一個像素點的值，都由其本身和鄰域內(nèi)的其他像素值經(jīng)過加權(quán)平均后得到[9]。對于DSA二維圖像，用二維零均值離散高斯函數(shù)做平滑濾波器，表達(dá)式如下：

式（1）中的x、y是DSA圖像對應(yīng)的矩陣坐標(biāo)，G是模板系數(shù)，表示用鄰域內(nèi)像素的加權(quán)平均灰度值去替代模板中心像素點的值。σ是DSA圖像像素的標(biāo)準(zhǔn)差。σ越大，高斯濾波器的頻帶越寬，平滑程度越好，通過調(diào)節(jié)σ的值，可在過平滑與欠平滑間取得折衷。經(jīng)過處理后的圖像既適合人眼的視覺、消除圖像的干擾又能保護(hù)圖像的邊緣。

圖1 圖像處理Fig.1 Processed image

1.2 提取冠狀動脈血管中心線

獲取DSA圖像中冠狀動脈血管中心線，是測量血管夾角的基礎(chǔ)。將冠狀動脈血管的中心線看作線條邊緣，冠狀動脈血管邊界看作階躍邊緣[10]。線條邊緣一階導(dǎo)數(shù)值為零，二階導(dǎo)數(shù)幅度最大，階躍邊緣一階導(dǎo)數(shù)幅度最大，二階導(dǎo)數(shù)值為零。根據(jù)此微分幾何特征，可將血管中心線提取出來。

對于經(jīng)平滑過程去除噪聲后的DSA灰度圖像（圖1(c)），通過Hessian矩陣的增強(qiáng)濾波器對血管進(jìn)行增強(qiáng)。圖像平滑后的像素點(x,y)的 Hessian 矩陣[11]定義如下：

式（2）中的H是曲面的曲率，Ixx,Ixy,Iyx,Iyy分別為平滑后的灰度圖像（圖1(c)）的二階方向?qū)?shù)。

Hessian矩陣的特征向量和特征值表示了圖像的本質(zhì)特征。最大特征值和對應(yīng)的特征向量表示三維曲面最大曲率的強(qiáng)度和方向，最小特征值對應(yīng)的特征方向與最強(qiáng)曲率對應(yīng)的方向垂直。Hessian矩陣的特征值和特征向量分別用λi(x,y)和vi(x,y)（i=1,2）表示，簡記為λi和vi。設(shè)特征值|λ1|≤λ2|，在冠狀動脈DSA圖像上，假設(shè)（x,y）是冠狀動脈血管上的點，λ1數(shù)值非常小，大的曲率垂直于血管，因此v2與局部血管的方向垂直，v1與局部血管軸的走向一致[12]。

通過高斯濾波去除噪聲的干擾，在Matlab中使用struct函數(shù)構(gòu)建多尺度立體數(shù)組，再使用sigma函數(shù)得到當(dāng)前圖像空間的值，在此基礎(chǔ)上對Hessian矩陣使用atan2函數(shù)，獲取不同方向的三維曲面的曲率，進(jìn)而獲得血管的走向。然后使用imerode函數(shù)、filter函數(shù)以及bwareaopen函數(shù)對圖像進(jìn)行腐蝕，消除冠狀動脈DSA圖像中垂直血管走向的干擾和較小的連通分支，最后使用bwmorph函數(shù)得到冠狀動脈血管的中心線（圖1(d)）。

2 計算冠狀動脈血管分叉處夾角

2.1 獲取冠狀動脈血管的分叉點

將DSA血管中心線圖像（圖2(a)）經(jīng)二值化處理，矩陣大小為512×512，冠狀動脈血管的中心線像素值為“1”，圖像背景像素值為“0”。以3×3矩陣（圖2(b)）為單位遍歷整個矩陣，找到滿足中心單元為“1”，周圍8個單位中有3處為“1”的3×3矩陣。經(jīng)過大量數(shù)據(jù)試驗，滿足此條件的是中心線的分叉點，中心線的縱截面最大是2個單位，避免了尋找到的點不是血管的分叉點，因此選擇3×3矩陣。使用Matlab中的size函數(shù)獲取DSA圖像的矩陣大小，然后使用四個for循環(huán)遍歷整個矩陣，以3×3矩陣為單位統(tǒng)計“1”的個數(shù)和位置，獲取冠狀動脈血管的分叉點。

圖2 血管中心線圖像矩陣Fig.2 Matrix of vessel central line image

2.2 尋找分支血管中心線上的兩個點

對于DSA圖像中冠狀動脈血管的每個分叉結(jié)點（圖2(a)中的a點），會存在兩個位于分支血管上且距分叉點一定距離的點（圖2(a)中的b、c點）。將血管交點為出發(fā)點，分上下兩個方向，以2×2的矩陣為單位（圖2中的矩形區(qū)域），尋找2×2矩陣中“1”的位置，如果有“1”，則將初始點移到此點，然后以該點為初始點，以包含該點的相鄰2×2矩陣為單位繼續(xù)尋找下一個為“1”的點，總共找尋8次、12次、16次（圖2中箭頭指示尋找過程），把最后為“1”的點（圖2(a)中的 b、c點）的坐標(biāo)輸出。在Matlab中，使用兩次for循環(huán)分別遍歷矩陣中的分叉點和以分叉點為出發(fā)點的2×2矩陣，再使用if函數(shù)尋找滿足條件的下一個為“1”的坐標(biāo)，最后使用fprintf函數(shù)將目標(biāo)點輸出。

2.3 淘汰較小的冠狀動脈血管

如圖3所示，圖3(a)為原始DSA血管圖像，可見冠狀動脈血管存在著較小的分支（圖3(a)中的a、b處），處理后的圖像會存在較短的中心線（圖3(b)中的a、b、c、d處）。通常這種冠狀動脈血管的細(xì)小分支臨床意義不大，因此將其淘汰掉。

在DSA圖像的中心線矩陣（512×512）中，以血管交點為中心，截取上下60行，左右24列的矩陣，在交點右側(cè)距離12個單位（較小分叉點的長度小于12個單位）做平行于y軸的直線，在這個范圍內(nèi)，如果直線（圖3(b)中的直線1）與血管中心線有兩個交點（圖3(b)中的e、f點），則該點（圖3(b)中g(shù)點）保留；如果直線（圖3(b)中的直線2）與血管中心線有一個或零個交點（圖3(b)中的h點），則說明此處的中心線較短，淘汰掉這樣的分叉點（圖3(b)中的i點）。最后保留下來的點，均為較長血管的交點。

圖3 血管中心線圖像Fig.3 Central line image

在Matlab中，首先使用兩個for循環(huán)遍歷整個矩陣，得到冠狀動脈血管的分叉點，然后以此分叉點為中心，橫坐標(biāo)增加12個單位，統(tǒng)計此長度中“1”的數(shù)量，最后將數(shù)量不是2的分叉點淘汰掉，保留較大的分叉點。

2.4 獲取冠狀動脈血管分叉處的角度

圖4所示計算結(jié)果，其中第1至第14列為經(jīng)過循環(huán)步驟計算后DSA圖像像素點的坐標(biāo)值，第15列為血管分叉處的角度，并在第16列說明是否淘汰了該交叉點。圖4中的1區(qū)域所指為交叉點坐標(biāo)，在這些像素中，2區(qū)域的點所指即較小的冠狀動脈血管，被淘汰掉。只保留冠狀動脈血管的分叉點及與分叉點有一定距離且在血管中心線上的兩個點的坐標(biāo)，在進(jìn)行8次、12次、16次迭代之后，獲得了的這三個點坐標(biāo)（圖4中的3區(qū)域）。在Matlab 中，使用sqrt和acos函數(shù)計算冠狀動脈血管分叉處的夾角，最后使用mean函數(shù)獲得冠狀動脈血管分叉處（圖4A中a、b、c所示）角度的平均值。

圖4 血管分叉角度測量Fig.4 Vessel angle measurement

3 討論

該文在Matlab平臺下，通過對DSA圖像進(jìn)行平滑，提取冠狀動脈血管中心線及分叉點，從而實現(xiàn)了DSA圖像冠狀動脈血管分叉角度的測量。本研究具有圖像處理精確度高、算法簡單而高效以及算法自動化程度高、易于擴(kuò)展等特點。

該文采用高斯濾波和中值濾波的方法對DSA圖像進(jìn)行處理，使得處理后的心血管圖像適合人眼的視覺，消除了圖像的干擾且保護(hù)了血管邊緣。針對Hessian矩陣的特點，消除了垂直血管走向的干擾分支與較短的連通分支，使得提取的冠狀動脈血管中心線精確度高，實用價值大。針對每一個夾角，算法可自動獲取多組坐標(biāo)組合，并求取角度的平均值，進(jìn)一步提高了角度測量的精確度。

根據(jù)血管分叉點處的矩陣特點，本文在Matlab平臺上，采用3×3矩陣獲得冠狀動脈血管分叉處的坐標(biāo)，使用2×2矩陣獲取冠狀動脈血管中心線上的兩個點的坐標(biāo)，用三個點坐標(biāo)求角度，算法簡單而高效。

算法自動化程度高，用戶只需讀入DSA圖像，平臺就會將處理后的圖像、冠狀動脈血管夾角處像素的坐標(biāo)，以及冠狀動脈血管分叉處的角度予以自動顯示。算法易于擴(kuò)展，將算法代碼在Matlab上生成文件包，并與Java web相融合，可簡單高效地運行于醫(yī)院網(wǎng)絡(luò)圖像系統(tǒng)。

該研究可對冠狀動脈血管DSA圖像進(jìn)行自動處理，提取血管分叉點及中心線，并最終給出血管分叉角度，這些信息均豐富了操作人員的經(jīng)驗，可為操作人員設(shè)定體位提供參考。另外該研究可供一線臨床科研人員分析冠狀動脈DSA臨床數(shù)據(jù)，也可以作為教學(xué)工具為學(xué)生展示DSA相關(guān)信息，界面友好，操作方便，具有一定的實用價值。