基于Hessian矩陣的DSA圖像冠狀動(dòng)脈直徑的測(cè)量

王力，孟慶民，郭永新，焦青

1.山東第一醫(yī)科大學(xué)（山東省醫(yī)學(xué)科學(xué)院）醫(yī)學(xué)信息工程學(xué)院，山東泰安271016；2.泰安市中心醫(yī)院介入放射科，山東泰安271000；3.山東第一醫(yī)科大學(xué)（山東省醫(yī)學(xué)科學(xué)院）放射學(xué)院，山東泰安271016

前言

心血管疾病具有高患病率、高致殘率和高死亡率等特點(diǎn)［1］。據(jù)《2017～2022年中國(guó)醫(yī)療服務(wù)市場(chǎng)行情動(dòng)態(tài)及發(fā)展前景預(yù)測(cè)報(bào)告》顯示，我國(guó)心血管疾病患病人數(shù)接近3 億，心血管疾病死亡率居首位，占居民疾病死亡構(gòu)成的40%以上，而且我國(guó)心血管疾病患病率及死亡率處于上升階段，心血管疾病的負(fù)擔(dān)日漸加重，已成為重大的公共衛(wèi)生問(wèn)題，防治心血管疾病刻不容緩［2］。

數(shù)字減影血管造影技術(shù)（Digital Subtraction Angiography,DSA）是一種使X射線序列圖片中的血管可視化的強(qiáng)大技術(shù)，是血管疾病無(wú)創(chuàng)診斷與介入治療手術(shù)導(dǎo)航的重要依據(jù)，是心血管疾病診斷的金標(biāo)準(zhǔn)，廣泛應(yīng)用于X射線序列成像中血管的可視化系統(tǒng)中［3］。冠狀動(dòng)脈血管壁堆積較多脂質(zhì)、糖類物質(zhì)或動(dòng)脈粥樣硬化逐漸明顯，易引起部分血管腔堵塞或狹窄。醫(yī)生根據(jù)狹窄處血管直徑的大小選擇是否置入支架，及置入支架的尺寸［4］。置入過(guò)大的支架直徑，冠狀動(dòng)脈血管可能會(huì)出現(xiàn)新破口或內(nèi)漏，損傷冠狀動(dòng)脈血管管腔內(nèi)壁，導(dǎo)致支架兩端主動(dòng)脈新的病變；置入過(guò)小的支架直徑，會(huì)使支架的擴(kuò)張應(yīng)力不足，支架固定不穩(wěn)固，易出現(xiàn)支架移位［5］。因此精確地判斷冠狀動(dòng)脈血管直徑，是確保支架置入手術(shù)的必要前提，這對(duì)于醫(yī)生的臨床決策具有重要的指導(dǎo)意義，對(duì)于冠狀動(dòng)脈搭橋手術(shù)的效果評(píng)價(jià)具有重要價(jià)值［6］。

該研究基于臨床實(shí)踐需要，在Windows7 環(huán)境下，以MATLAB為平臺(tái)，通過(guò)對(duì)DSA圖像進(jìn)行去噪、均衡化等處理，并通過(guò)Hessian矩陣，提取冠狀動(dòng)脈血管的中心線和骨架，從而實(shí)現(xiàn)DSA 圖像冠狀動(dòng)脈血管直徑的測(cè)量。

1 DSA圖像平滑

由于人為或儀器的各種原因，DSA 圖像在獲取和信號(hào)傳輸過(guò)程中會(huì)受到噪聲的干擾，從而降低DSA 圖像的質(zhì)量，這對(duì)DSA 圖像的處理與視覺(jué)效果均產(chǎn)生不利的影響。通常采用平滑技術(shù)對(duì)DSA圖像進(jìn)行處理，以消除圖像的噪聲［7］。

高斯濾波是一種平滑濾波器，取濾波器窗口內(nèi)像素的均值作為輸出，使用高斯濾波器對(duì)DSA 圖像進(jìn)行濾波，其效果是降低DSA 圖像灰度的“尖銳”變化，能夠有效抑制噪聲，平滑圖像［8］。對(duì)于DSA二維圖像，用高斯函數(shù)做平滑濾波器，表達(dá)式如下：

式中，x、y是整個(gè)DSA 圖像對(duì)應(yīng)的矩陣坐標(biāo)，H是模板的系數(shù)，表示每一個(gè)像素點(diǎn)的值，都由鄰域內(nèi)像素的加權(quán)平均灰度值代替。σ是DSA圖像像素的標(biāo)準(zhǔn)差，σ取值較小時(shí)，高斯濾波器的頻帶較窄，信號(hào)邊緣的定位精度較高，但平滑作用較小，對(duì)噪聲的抑制能力較差（欠平滑）；反之，σ取值較大時(shí)，對(duì)信號(hào)的平滑作用較大，去噪效果較好，但會(huì)使信號(hào)邊緣變得模糊，且邊緣移位現(xiàn)象嚴(yán)重（過(guò)平滑）。通過(guò)調(diào)節(jié)σ值，可在過(guò)平滑與欠平滑間取得折衷［9］。

在MATLAB 平臺(tái)中讀入DSA 圖像（圖1），使用rgb2gray函數(shù)將圖像轉(zhuǎn)化為灰度圖像，在此基礎(chǔ)上使用fspecial函數(shù)與imfilter函數(shù)進(jìn)行高斯濾波，選取較大的fspecial函數(shù)參數(shù)值（HSIZE參數(shù)設(shè)為1 000），得到DSA 圖像的背景（圖2），選取較小的fspecial 函數(shù)參數(shù)值（HSIZE 參數(shù)設(shè)為10），得到消除噪聲后的DSA圖像（圖3）。

圖1 原始DSA圖像Fig.1 Original digital subtraction angiography image

圖2 背景圖像Fig.2 Background image

圖3 消除噪聲后的圖像Fig.3 Image after noise removal

2 提取冠狀動(dòng)脈血管的中心線

Hessian矩陣是由圖像在各個(gè)像素點(diǎn)的二階偏導(dǎo)構(gòu)成［10］，可表示為原圖與高斯模板函數(shù)二次偏導(dǎo)數(shù)的卷積，反映各個(gè)像素點(diǎn)處的局部灰度幾何信息［11］：

式中，x、y是平滑后灰度DSA 圖像對(duì)應(yīng)的矩陣坐標(biāo)，H是曲面的曲率，I是平滑后的灰度圖像，是高斯函數(shù)的二階偏導(dǎo)數(shù)，Ixx、Ixy、Iyy分別為與I進(jìn)行卷積，即灰度圖像（圖3）的二階方向?qū)?shù)。

通過(guò)求取Hessian 矩陣的兩個(gè)特征值和特征向量，可對(duì)圖像中不同特征區(qū)域所具有的結(jié)構(gòu)特性進(jìn)行描述［12］。將Hessian 矩陣的特征值分解，得到兩個(gè)正交的方向向量，分別用ν1、ν2 表示，對(duì)應(yīng)的特征值用λ1、λ2 表示。特征值λ1 與其對(duì)應(yīng)的特征向量ν1分別表示曲面曲率小的強(qiáng)度和方向，即ν1 和局部血管軸平行；特征值λ2 與特征向量ν2 分別表示曲面曲率大的強(qiáng)度和方向，即ν2 和局部血管軸垂直［13］。

在MATLAB 平臺(tái)中，對(duì)于圖3圖像，使用struct函數(shù)構(gòu)建多尺度立體數(shù)組，并采用sigma函數(shù)得到當(dāng)前空間的值，在此基礎(chǔ)上對(duì)Hessian矩陣使用atan2函數(shù)，獲取不同方向三維曲面的曲率，進(jìn)而獲得血管的走向。然后采用imerode 函數(shù)、filter 函數(shù)和bwareaopen 函數(shù)對(duì)圖像進(jìn)行腐蝕，消除DSA 圖像中較小的冠狀動(dòng)脈血管，最后使用bwmorph 函數(shù)得到冠狀動(dòng)脈血管的中心線（圖4）。

圖4 血管中心線圖像Fig.4 Vessel centerline image

3 冠狀動(dòng)脈血管直徑的測(cè)量

3.1 DSA圖像增強(qiáng)

對(duì)經(jīng)過(guò)高斯濾波的DSA圖像進(jìn)行直方圖均衡化處理，該過(guò)程是采用灰度變換函數(shù)修正圖像的直方圖，使其成為灰度均勻分布的均衡直方圖［14］，然后用此均衡直方圖修正灰度圖像，即對(duì)圖像進(jìn)行非線性拉伸，重新分配圖像象元值，使一定灰度范圍內(nèi)象元值的數(shù)量大致相等，此時(shí)圖像的信息熵最大，圖像包含的信息量最大，圖像對(duì)比度增強(qiáng)，層次分明，圖像更加清晰［15］。對(duì)灰度圖像進(jìn)行直方圖均衡化［16］的函數(shù)表達(dá)式為：

式中，nk為圖像中出現(xiàn)該灰度的像素?cái)?shù)（k=0,1,…,n-1），n是圖像中的像素總數(shù)，是概率論的頻數(shù)，si是灰度級(jí)的概率分布。

功能實(shí)現(xiàn)時(shí)，分別采用motion算子與unsharp算子來(lái)消除血管周圍的失真晶格，及增大血管邊緣的對(duì)比度。在MATLAB中分別采用histeq函數(shù)、fspecial函數(shù)及imfilter函數(shù)，對(duì)灰度圖像（圖5）進(jìn)行均衡化處理，并增大冠狀動(dòng)脈血管邊緣對(duì)比度（圖6）。

圖5 灰度DSA圖像Fig.5 Grayscale digital subtraction angiography image

圖6 圖像均衡化Fig.6 Image equalization

3.2 提取冠狀動(dòng)脈血管的骨架

由于高斯濾波模板尺寸選取較大，因此濾波后圖像接近于原圖像的背景［17］（圖2），經(jīng)平滑與直方圖均衡化處理后的圖像，高斯濾波模板尺寸選取較小，起到去噪的作用［18］（圖3），將這兩步獲得的后者矩陣減去前者矩陣，即可得到層次分明的血管圖像（圖7）。

圖7 血管圖像Fig.7 Vessel image

冠狀動(dòng)脈血管圖像在物理意義上是有相同特點(diǎn)的區(qū)域集合，這些區(qū)域的形狀或灰度或顏色大致相同、紋理屬性相近等，通常采用邊緣檢測(cè)算子提取血管骨架。邊緣檢測(cè)算子包括Prewitt、Sobel、Roberts、Canny 和拉普拉斯算子等。Prewitt 算子邊緣的定位準(zhǔn)確度不夠高、檢測(cè)出的邊緣較寬、間斷點(diǎn)略多；Roberts算子檢測(cè)出的邊緣不是很平滑、邊緣比較寬；Canny算子運(yùn)算很復(fù)雜；拉普拉斯算子對(duì)噪聲過(guò)于敏感［19］。而Sobel算子算法可對(duì)圖像進(jìn)行平滑去噪，邊緣檢測(cè)效果良好，且算法簡(jiǎn)單易實(shí)現(xiàn)［20］。該文使用Sobel 算子進(jìn)行圖像邊緣檢測(cè)及提取血管骨架。在MATLAB 中使用edge 函數(shù)對(duì)血管圖像（圖7）處理，進(jìn)行邊緣提取，即可得到冠狀動(dòng)脈血管的骨架（圖8）。

圖8 血管骨架圖像Fig.8 Vessel skeleton image

3.3 測(cè)量冠狀動(dòng)脈血管的直徑

由于冠狀動(dòng)脈血管中心線圖像是以“0”、“1”數(shù)字為單位的矩陣，“0”代表圖像背景，“1”代表冠狀動(dòng)脈血管的中心線。為了區(qū)分冠狀動(dòng)脈血管的骨架與中心線，遍歷整個(gè)矩陣，把冠狀動(dòng)脈血管中心線矩陣中“1”的位置，全部用“2”替換。替換后的冠狀動(dòng)脈血管中心線圖像矩陣與冠狀動(dòng)脈血管骨架圖像矩陣相加，即可得到冠狀動(dòng)脈中心線血管骨架圖像（圖9）。采用Hessian 矩陣提取冠狀動(dòng)脈血管中心線時(shí)，得到Hessian 矩陣的兩個(gè)特征向量（圖10 中的v1、v2），v1 與局部血管軸的走向一致，v2 與局部血管軸垂直［21］。以v2 表示的方向作為斜率，過(guò)冠狀動(dòng)脈血管中心線上的點(diǎn)（圖10中的a點(diǎn)），作一條直線（圖10中的直線1），直線恰好與冠狀動(dòng)脈血管骨架有兩個(gè)交點(diǎn)（圖10中的b、c點(diǎn)），該交點(diǎn)間的距離即冠狀動(dòng)脈血管的直徑。

圖9 骨架中心線圖像Fig.9 Skeleton centerline image

圖10 血管直徑測(cè)量圖示Fig.10 Vessel diameter measurement

在MATLAB中，使用ginput函數(shù)獲取光標(biāo)位置，以此點(diǎn)為中心，使用兩個(gè)for 循環(huán)，遍歷矩陣，找到位于冠狀動(dòng)脈血管中心線上（矩陣中的值為“2”）且距離光標(biāo)位置最近的點(diǎn)(x1,y1)。將此點(diǎn)與特征向量代入直線方程中，得到直線與冠狀動(dòng)脈血管骨架的兩個(gè)交點(diǎn)坐標(biāo)，使用dist函數(shù)求取這兩個(gè)交點(diǎn)的距離，即冠狀動(dòng)脈血管的直徑。

4 討論

本研究通過(guò)對(duì)DSA 圖像進(jìn)行平滑、均衡化及銳化等處理，提取冠狀動(dòng)脈血管的骨架和中心線，從而實(shí)現(xiàn)DSA圖像冠狀動(dòng)脈血管直徑的測(cè)量。本研究具有以下特點(diǎn)：

圖像處理效果好，適用范圍廣泛。在眾多的圖像平滑處理方法中，根據(jù)實(shí)際需求，本研究選取高斯濾波方法對(duì)圖像進(jìn)行處理，使得處理后的心血管圖像適合人眼的視覺(jué)、消除圖像的干擾并且保護(hù)血管邊緣。為了突出血管邊緣并進(jìn)一步削弱圖像中的雜質(zhì)，采用直方圖均衡化的方法，使得圖像對(duì)比度增強(qiáng)，層次分明。為了使冠狀動(dòng)脈血管的邊緣更為清晰，又使用濾鏡消除血管周圍失真的晶格。將處理后的圖像減去高斯濾波后的背景圖像，進(jìn)一步弱化背景的干擾，使得提取的血管骨架更加明顯。本研究采取的思路與方法，極大地弱化或消除圖像中的雜質(zhì)，突出血管邊緣，適用于絕大多數(shù)醫(yī)學(xué)圖像的處理。

直徑測(cè)量精確度高。針對(duì)Hessian 矩陣的特點(diǎn)，將冠狀動(dòng)脈血管的中心線看作線條邊緣，冠狀動(dòng)脈血管邊界看作階躍邊緣。根據(jù)線條邊緣一階導(dǎo)數(shù)值為零，二階導(dǎo)數(shù)幅度最大，階躍邊緣一階導(dǎo)數(shù)幅度最大，二階導(dǎo)數(shù)值為零的微分幾何特征，消除掉較小的血管，使得冠狀動(dòng)脈血管中心線精確度高。冠狀動(dòng)脈血管中心線上每個(gè)點(diǎn)的血管走向可由Hessian矩陣的特征向量得到，使用直線與冠狀動(dòng)脈血管骨架的交點(diǎn)測(cè)量血管直徑，結(jié)果準(zhǔn)確可靠，精度達(dá)到亞像素級(jí)，從而為操作人員提供精準(zhǔn)的冠狀動(dòng)脈血管信息。

本研究的功能實(shí)現(xiàn)操作簡(jiǎn)便且易于擴(kuò)展。用戶只需讀入DSA 圖像，平臺(tái)就會(huì)將處理后的DSA 圖像以及冠狀動(dòng)脈血管中心線圖像予以顯示。醫(yī)生只需用鼠標(biāo)點(diǎn)擊選擇感興趣的血管區(qū)域，平臺(tái)即可自動(dòng)計(jì)算并顯示該位置的血管直徑，極大地提高操作人員的工作效率。將算法代碼在MATLAB上生成文件包，并與Java web 相融合，即可簡(jiǎn)單高效地運(yùn)行在醫(yī)生工作站。也可將冠狀動(dòng)脈血管的骨架、中心線圖像以及直徑等數(shù)據(jù)存放到數(shù)據(jù)庫(kù)中，供科研人員進(jìn)行分析處理。

本研究可供臨床醫(yī)師及科研人員對(duì)冠狀動(dòng)脈血管的DSA 圖像進(jìn)行回顧性分析與處理，也可為學(xué)生展示冠狀動(dòng)脈血管的DSA圖像，易于操作，具有一定的實(shí)用價(jià)值。