DSA圖像冠狀動脈血管夾角的測量

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

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

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

該研究基于臨床實(shí)踐需要，在Windows7環(huán)境下，以Matlab為平臺，首先對DSA圖像進(jìn)行平滑處理，然后提取冠狀動脈血管的中心線及分叉點(diǎn)，從而實(shí)現(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)平均的過程，每一個(gè)像素點(diǎ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)平均灰度值去替代模板中心像素點(diǎ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)。圖像平滑后的像素點(diǎn)(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）是冠狀動脈血管上的點(diǎn)，λ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 計(jì)算冠狀動脈血管分叉處夾角

2.1 獲取冠狀動脈血管的分叉點(diǎn)

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

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

2.2 尋找分支血管中心線上的兩個(gè)點(diǎn)

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

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

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

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

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

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

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

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

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

3 討論

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

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

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

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

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