基于EMD和模極大值的造影圖像血管提取

邵欣,尹清波,魯明羽

(大連海事大學(xué) 信息科學(xué)技術(shù)學(xué)院，遼寧 大連 116026)

?

基于EMD和模極大值的造影圖像血管提取

邵欣,尹清波,魯明羽

(大連海事大學(xué) 信息科學(xué)技術(shù)學(xué)院，遼寧 大連 116026)

摘要：針對冠脈造影圖像噪音大且血管邊緣模糊問題，提出一種基于經(jīng)驗?zāi)B(tài)分解的冠脈造影圖像血管提取方法以較完整有效地提取出造影圖像中的血管信息。該方法首先對圖像預(yù)處理后的造影圖像進行經(jīng)驗?zāi)B(tài)分解，將其分解成若干圖像層，保留包含血管信息的圖像層并丟棄含非血管信息圖像層；其次對保留的圖像層進行血管特征提取，依據(jù)連通區(qū)域的一系列屬性來提取圖像層中血管信息以獲取清晰明顯的血管段；最后通過層間信息相互參照，重構(gòu)血管脈絡(luò)明顯的心血管圖像。實驗證明了方法有很高的精確度和實際的醫(yī)用價值。

關(guān)鍵詞：血管提取；經(jīng)驗?zāi)B(tài)分解；冠脈造影圖像；特征提取；連通區(qū)域；層間參照；血管重構(gòu)；圖像層

中文引用格式：邵欣,尹清波,魯明羽. 基于EMD和模極大值的造影圖像血管提取[J]. 智能系統(tǒng)學(xué)報， 2015, 10(6): 851-857.

冠脈造影圖像是臨床診斷冠心病的“金標(biāo)準(zhǔn)”。其成像過程中，受心臟跳動、成像設(shè)備等因素影響，圖像不可避免地會引入一些噪聲，使得圖像中的血管脈絡(luò)一定程度上模糊不清,尤其是血管邊緣與細節(jié)的模糊，影響醫(yī)生對病情的診斷。

近年來，不少學(xué)者致力于研究冠脈造影圖像中血管脈絡(luò)特征的提取。常見的提取冠脈造影圖像中的血管信息方法主要有數(shù)學(xué)形態(tài)學(xué)方法[1-3]、匹配濾波方法[4-6]、基于中心線的方法[7]等。其中，數(shù)學(xué)形態(tài)學(xué)方法是通過膨脹和腐蝕2種最主要的形態(tài)學(xué)操作實現(xiàn)圖像中目標(biāo)的擴張或收縮，以先腐蝕后膨脹的開操作和先膨脹后腐蝕的閉操作作為基礎(chǔ)，實現(xiàn)圖像中的邊緣檢測、圖像增強和特征提取等相關(guān)操作；匹配濾波方法是將多個線性濾波器與卷積運算相結(jié)合的方式提取圖像中的目標(biāo)信息，通常還需要其他相關(guān)圖像處理操作才能得到圖像中的血管信息；基于中心線方法的基本思想是提取出圖像中整個血管樹的中心線，然后通過測量血管在某位置的半徑重構(gòu)圖像中的血管脈絡(luò)。但由于血管的形狀復(fù)雜，血管的直徑和造影劑的密度也會發(fā)生變化等，使得現(xiàn)有的方法并不能有效且完整地將冠脈造影圖像中的血管信息提取出來。

經(jīng)驗?zāi)B(tài)分解是一種新的非平穩(wěn)信號分析處理方法，可用于冠狀動脈造影圖像的圖像處理過程[8-10]。本文提出一種基于經(jīng)驗?zāi)B(tài)分解的血管特征提取方法：采用經(jīng)驗?zāi)B(tài)分解將造影圖像分解成若干圖像層，形成圖像層組，該組中有的圖像層包含了大量有效的血管信息而有的則包含的是肋骨、肺或其他組織等形成的非血管信息，丟棄含非血管信息的圖像層，對含有血管信息的圖像層進行血管特征提取，依據(jù)連通區(qū)域一系列屬性提取出圖像層中清晰明顯的血管段，將各圖像層提取出的血管信息相互參照并重構(gòu)血管脈絡(luò)明顯的心血管圖像。此外，本方法對冠脈造影圖像經(jīng)驗?zāi)B(tài)分解之前須進行圖像預(yù)處理操作，這是因為經(jīng)驗?zāi)B(tài)分解方法的基本原理是根據(jù)目標(biāo)的變化頻率而將圖像中信息逐層分離出來，但受現(xiàn)有投影設(shè)備的限制，在冠脈造影圖像的四周存在從有像素區(qū)到黑色無像素區(qū)的過渡區(qū)域，該區(qū)域的頻率變化最大，如果不進行圖像預(yù)處理進行圖像像素填充使該區(qū)域?qū)崿F(xiàn)平滑過渡，會在經(jīng)驗?zāi)B(tài)分解結(jié)果中出現(xiàn)嚴(yán)重邊緣效應(yīng)。

1冠脈造影圖像預(yù)處理

在對冠脈造影圖像進行特征提取之前，須對圖像中無像素黑色區(qū)域進行像素填充，使該區(qū)域像素點的灰度值實現(xiàn)平滑過渡，以抑制此后經(jīng)驗?zāi)B(tài)分解過程中邊緣效應(yīng)的出現(xiàn)。

受投影設(shè)備的限制，冠脈造影圖像的四周會出現(xiàn)黑色區(qū)域，該區(qū)域不包含任何血管信息且頻率變化最大，由于經(jīng)驗?zāi)B(tài)分解是根據(jù)目標(biāo)的變化頻率而將圖像中信息逐層分離出來，如果不對此區(qū)域進行像素填充則會在分解結(jié)果中出現(xiàn)嚴(yán)重邊緣效應(yīng)。圖1所示的是臨床常見的冠脈造影圖像。觀察圖像可知在黑色區(qū)域與圓形亮色區(qū)域的相接處，像素點的灰度值變化最為劇烈，在圖像經(jīng)驗?zāi)B(tài)分解過程中此處會出現(xiàn)偽血管信息。

為此，本文對圖像中無像素的黑色區(qū)域進行像素填充，填充的過程如下：

1)將圖1的冠脈造影圖像劃分為如圖2所示的4個區(qū)域，箭頭表示該區(qū)的像素填充方向。每個區(qū)域的填充方向各不相同，依次為左上、右上、左下和右下。

圖1　冠脈造影圖像Fig.1　Coronary angiography images

圖2　冠脈造影圖像分區(qū)Fig.2　Partition of angiography images

2)定義像素填充模板。本文是在待處理像素點的鄰近位置選取大小合適的矩形填充模板，由此模板計算得出待填充點的像素值。例如，對填充方向為左上方的1區(qū)來說，采取的操作是在待填充點的右下方取一個大小為m×n的矩陣。根據(jù)式(1)計算該矩陣的平均灰度值，將該值作為待填充點的灰度值。如圖3所示，P是待填充點，本文中取m=4,n=4。

圖3　4×4大小的矩陣模塊Fig.3　Partition’s 4×4 matrix module

(1)

式中:P(i-k1,j-k2)表示點(i-k1,j-k2)的灰度值。據(jù)此對圖像1區(qū)由下至上、從右到左掃描，以填充黑色區(qū)域中的所有待填充點。對其他待填充區(qū)域進行類似處理，唯一不同的是矩陣模塊與待填充點的鄰接方向發(fā)生了變化。

3)設(shè)定閾值并通過閾值T將待填充區(qū)域內(nèi)的所有點分類成待填充點和保留點。如果一個像素點的灰度級小于閾值T則使用填充模板進行灰度級填充，如果大于閾值T則保留該像素點的原灰度級，不做灰度級填充處理。閾值T的設(shè)定方法是：對于選定的圖像填充區(qū)域，首先統(tǒng)計該區(qū)域內(nèi)每個灰度級像素的出現(xiàn)概率，即得到該區(qū)域的灰度直方圖，然后概率小于某一數(shù)值的所有灰度級中的最大值作為閾值。本文實驗中概率取1.16×10-2，約為0.011 6。

4)依據(jù)圖像填充模對圖像①區(qū)由下至上、從右到左掃描，以填充黑色區(qū)域中的所有待填充點。對其他待填充區(qū)域進行類似處理，唯一不同的是矩陣模塊與待填充點的鄰接方向發(fā)生了變化。

5)圖像二次填充以去除毛刺。經(jīng)過4)對冠脈造影圖像進行初次掃描及像素填充之后，部分像素點的過渡處仍存在小的毛刺，為了更好地實現(xiàn)圖像像素平滑過渡的效果，重復(fù)4)對圖像進行二次掃描填充。冠脈造影圖像經(jīng)像素填充預(yù)處理后的圖像如圖4所示。下文中提到的冠脈造影圖像均指預(yù)處理后的圖像。

圖4　預(yù)處理后圖像Fig.4　The result of image preprocessing

2造影圖像經(jīng)驗?zāi)B(tài)分解

本文采用經(jīng)驗?zāi)B(tài)分解方法將造影圖像分解成若干圖像層，利用了血管以下特性：在三維空間中為管狀結(jié)構(gòu)，經(jīng)X光照射后形成冠脈造影圖像中的血管段，其在血流方向上灰度級變化緩慢，而在垂直于血管血流方向上灰度級變化劇烈。如圖5所示，圖像灰度級在λ1方向上變化較平緩，在λ2方向上變化劇烈。

以圖4為實驗對象，圖像通過經(jīng)驗?zāi)B(tài)分解后的結(jié)果是原圖像被分解成不同的圖像層，每一層中包含圖像中不同變化頻率血管特征信息，圖像分解結(jié)果如圖6所示。通過對照原冠脈造影圖像與圖像分解后的圖像層組發(fā)現(xiàn)，圖像中的血管信息主要分布在imf1、imf2和imf3中。因此，在對各分解層信號的后續(xù)處理過程中，本文只針對前3層進行處理。

圖5　心血管的空間結(jié)構(gòu)Fig.5　Spatial structure of the cardiovascular

(a) imf1

(b)imf2

(c)imf3

(d)imf4

(e)imf5

(f)imf6圖6　圖像EMD分解Fig.6　Image EMD decomposition

3圖像層中血管特征提取

針對冠脈造影圖像強背景噪音的特點，經(jīng)驗?zāi)B(tài)分解雖然能很好地將造影圖像中的血管信息和背景信息分離出來，并能將圖像中的血管根據(jù)變化頻率的不同進行分級，生成含有不用尺度級血管的圖像層組，但是經(jīng)驗?zāi)B(tài)分解算法在分解圖像產(chǎn)生若干個圖像層的過程中某些情況下可能發(fā)生模態(tài)混疊現(xiàn)象。此現(xiàn)象會造成圖像層中的血管特征信息失真，使得圖像層中的信息不能真實反映冠脈造影圖像中血管脈絡(luò)信息，所以采用小波變換模極大值方法實現(xiàn)圖像層間的信息參照以消除模態(tài)混疊造成的信息失真現(xiàn)象。

對冠脈造影圖像分解后的圖像層進行層間信息參照，主要是imf1、imf2和imf3之間信息的參照。以圖6中圖像為例，得到如圖7所示的極大值序列。

(a) imf1

(b)imf1極大值

(c)imf2

(d)imf2極大值

(e)imf3

(f)imf3極大值圖7　各圖像層中的極大值Fig.7　The maxima of each image layer

針對imf3中的極大值點，如果在圖像層imf2相同圖像位置上，該圖像位置周圍一定的區(qū)域范圍內(nèi)檢測到極大值點，則imf3中的此類點屬于血管特征信息，將該類點標(biāo)記并生成標(biāo)記點集合Mark1。同樣的參照方法對imf2圖像層中的極大值進行標(biāo)記，以imf1圖像層中的極大值點為參照信息，對Mark=Mark1∪Mark2中的極大值點，如果在圖像層imf1相同圖像位置上，周圍一定的區(qū)域范圍內(nèi)檢測到極大值點，則imf2中的此類點為標(biāo)記點，生成的標(biāo)記點結(jié)果標(biāo)記點集合Mark2中。由此可知，造影圖像血管信息層間的信息參照結(jié)果是一個標(biāo)記點集合imf1，其中Mark=Mark1∪Mark2。在圖像重構(gòu)過程中標(biāo)記點集合即為結(jié)果圖像中血管樹生成過程的參照信息，包含標(biāo)記點的連通區(qū)域即為血管脈絡(luò)樹中的特征信息，否則該連通區(qū)域?qū)儆诒尘靶畔ⅲ鰜G棄處理。采用類似小波變換的模極大值思想，對imf圖像層進行層間極大值參照處理，結(jié)果得到一個標(biāo)記點集合Mark{mi(x,y)|i=1,2,…,k,1≤x≤Mand1≤y≤N}。依據(jù)標(biāo)記點集合進行圖像血管脈絡(luò)重構(gòu)。

4圖像重構(gòu)

依據(jù)小波變換模極大值的思想對各圖像層間的血管信息進行層間參照，得到標(biāo)記點集合。在對各層圖像進行信息重構(gòu)以實現(xiàn)血管脈絡(luò)樹重構(gòu)時，參照標(biāo)記點集合，實現(xiàn)各層信息的信息綜合，從而得到最終的血管樹，獲取冠脈造影圖像中的完整血管信息。算法1詳細描述了本文中冠脈造影圖像中的完整血管信息的提取過程。

算法1基于EMD和模極大值的血管提取算法輸入：插值預(yù)處理后圖像；

輸出：含有血管脈絡(luò)樹的結(jié)果圖像。

1)插值預(yù)處理后圖像進行條件初始化: r0(t)=x(t),i=1；

2)whileri(t)≥2 // ri(t)是x(t)的“殘余量”

3)do { //提取第i個本征模態(tài)函數(shù)；

4)條件進一步初始化: h0(t)=ri(t),k=1；

5 )do{

6)計算hk-1(t)的局部極大值與局部極小值；

7)利用三次樣條對局部極大值與局部極小值序列進行插值；

8)分別形成hk-1(t)的上包絡(luò)線和下包絡(luò)線；

9)計算hk-1(t)的上、下包絡(luò)線的均值，求得平均包絡(luò)線mk-1(t)；

10 )定義hk(t)=hk-1(t)-mk-1(t)；

11) k=k+1；}

12 ) while hk(t)不滿足本征模態(tài)函數(shù)IMF的標(biāo)準(zhǔn)；

13 )IMFi(t)=hk(t); //生成本征模態(tài)函數(shù)

14)ri(t)=ri-1(t)-IMFi(t)；

15 ) i=i+1;}

16) 分解完成，生成imf圖像層組合殘差圖像；

17)圖像層組分成包含血管信息的有效層和包含背景信息的無效層；

18)有效圖像層進行極大值檢測，層間信息相互參照，生成標(biāo)記點集合；

19)依據(jù)Mark集合重構(gòu)血管脈絡(luò)樹；

20)生成結(jié)果圖像，血管信息提取過程完成。

依據(jù)標(biāo)記點集合對冠脈造影圖像進行重構(gòu)，將各個有效圖像層中的血管信息進行綜合，生成結(jié)果圖像。圖8 給出圖像重構(gòu)過程中不參照Mark集合和參照Mark集合的結(jié)果圖像效果對比。在圖8(a)中，結(jié)果圖像是直接將經(jīng)驗?zāi)B(tài)分解后的有效圖像層中提取出的血管信息直接相加得到，圖8(b)中的結(jié)果圖像則在圖像重構(gòu)過程中，參照標(biāo)記點集合，將包含標(biāo)記點的連通區(qū)域保留在結(jié)果圖像中，不包含標(biāo)記點的連通區(qū)域看做噪音處理。對比發(fā)現(xiàn)，通過模極大值進行圖層間的信息相互參照，生成的標(biāo)記點集合能很好地實現(xiàn)圖像血管脈絡(luò)樹的重構(gòu)并有效抑制噪音。

圖8　結(jié)果對比Fig.8　Result contrast

5實驗結(jié)果分析

為了更好地說明本文提出的算法能有效提取冠脈造影圖像中的血管信息，將本文算法的血管提取結(jié)果與現(xiàn)有常用的Hessian矩陣、Gabor濾波器、直方圖均衡化(HE)、線性反銳化掩模(LUM)等算法結(jié)果作比較。

圖9給出了Hessian矩陣、Gabor濾波器、直方圖均衡化(HE)、線性反銳化掩模(LUM)和本文提出的基于EMD和模極大值的血管提取算法的實驗結(jié)果對比。

圖9　算法對比Fig.9　Algorithm effect contrast

以冠脈造影圖像為例進行血管信息提取結(jié)果對比，圖9(b)給出了Hessian矩陣對造影圖像的血管提取結(jié)果，該方法能將血管與背景對比度非常強的區(qū)域中的血管段提取出來，有很好的血管連續(xù)性，但是不足的是對于血管與背景對比度稍微較弱的區(qū)域則無能為力，能提取出的血管信息非常有限，提取效果不太理想。圖9(c)中給出了Gabor濾波器對造影圖像中的血管信息進行過濾提取，能夠?qū)⒀苊}絡(luò)的輪廓提取出來，但是得到的結(jié)果圖像中血管特征并沒有很好地凸顯出來，血管樹的連續(xù)性也較差。圖9(d)給出了造影圖像直方圖均衡化后得到的結(jié)果圖像，該圖像在一定程度上提高了冠脈造影圖像的圖像效果，但是對于背景噪音強、圖像對比度低的圖像區(qū)域還是存在背景噪音過強的缺點，整體來說還是能夠?qū)⒀苊}絡(luò)結(jié)構(gòu)較好地表達出來。圖9(e)給出了線性反銳化掩模對造影圖像的處理結(jié)果，該方法能夠?qū)崿F(xiàn)造影圖像中血管的細節(jié)和邊緣突出表達，但是分段函數(shù)產(chǎn)生的量化誤差較大時會將圖像中的邊緣區(qū)域?qū)Ρ榷冗^度增強，產(chǎn)生偽影現(xiàn)象。圖9(f)給出了本文算法對造影圖像的血管提取結(jié)果。

6結(jié)束語

針對冠脈造影圖像噪音大和血管邊緣模糊等問題，提出的基于經(jīng)驗?zāi)B(tài)分解的血管特征提取方法，是將造影圖像分解成若干圖像層，形成圖像層組，對圖像層進行血管提取，將各層中有效的血管信息匯總在結(jié)果圖像中，得到血管脈絡(luò)完整的心血管圖像。

對比其他算法的實驗結(jié)果，文中算法得到的結(jié)果圖像有極好的視覺效果，去除了造影圖像中強背景噪音，使得圖像更便于血管信息的直觀分析和觀察，同時提取出的血管樹結(jié)構(gòu)有很好的連續(xù)性和完整性，造影圖像中血管特征信息，如血管方向、狹窄程序和血管分叉點都有很高的精確度，算法具有很好的實際醫(yī)用價值。

參考文獻:

[1]ZHOU Shoujun, YANG Jian, WANG Yongtian, et al. Automatic segmentation of coronary angiograms based on fuzzy inferring and probabilistic tracking[J]. Biomedical Engineering Online, 2010, 9: 40.

[2]YANG Guanyu, KITSLAAR P, FRENAY M, et al. Automatic centerline extraction of coronary arteries in coronary computed tomographic angiography[J]. The International Journal of Cardiovascular Imaging, 2011, 28(4): 921-933.

[3]MARIN D, AQUINO A, GEGUNDEZ-ARIAS M E, et al. A new supervised method for blood vessel segmentation in retinal images by using gray-level and moment invariants-based features[J]. IEEE Transactions on Medical Imaging, 2011, 30(1): 146-158.

[4]Van Der ZWET P M J, NETTESHEIM M, GERBRANDS J J, et al. Derivation of optimal filters for the detection of coronary arteries[J]. IEEE Transactions on Medical Imaging, 1998, 17(1): 108-120.

[5]WANG Li, BHALERAO A, WILSON R. Analysis of retinal vasculature using a multiresolution Hermite model[J]. IEEE Transactions on Medical Imaging, 2007, 26(2): 137-152.

[6]RANGAYYAN R M, AYRES F J, OLOUMI F, et al. Detection of blood vessels in the retina with multiscale Gabor filters[J]. Journal of Electronic Imaging, 2008, 17(2): 023018.

[7]GREENSPAN H, LAIFENFELD M, EINAV S, et al. Evaluation of center-line extraction algorithms in quantitative coronary angiography[J]. IEEE Transactions on Medical Imaging, 2001, 20(9): 928-941.

[8]康文煒. 冠狀動脈造影圖像的分割方法研究[D]. 長春: 吉林大學(xué), 2010.

KANG Wenwei. The segmentation methods of coronary angiogram images[D]. Changchun, China: Jilin University, 2010.

[9]周壽軍, 楊俊, 陳武凡, 等. 一種從X射線血管造影圖像中自動分割冠狀動脈的新方法[J]. 中國科學(xué) E輯: 信息科學(xué), 2007, 37(12): 1531-1543.

ZHOU Shoujun, YANG Jun, CHEN Wufan, et al. New approach to the automatic segmentation of coronary artery in X-ray angiograms[J]. Science in China Series F: Information Sciences, 2008, 51(1): 25-39.

[10]SHIH M Y, TSENG D C. A wavelet-based multiresolution edge detection and tracking[J]. Image and Vision Computing, 2005, 23(4): 441-451.

[11]賈桂鋒, 袁勝發(fā), 唐成文. 基于高階導(dǎo)數(shù)改進的經(jīng)驗?zāi)B(tài)分解算法[J]. 數(shù)據(jù)采集與處理, 2012, 27(1): 117-120.

JIA Guifeng, YUAN Shengfa, TANG Chengwen. Modified empirical mode decomposition algorithm based on higher derivative[J]. Journal of Data Acquisition & Processing, 2012, 27(1): 117-120.

[12]鐘佑明, 趙強, 周建庭. 實時經(jīng)驗?zāi)B(tài)分解的實現(xiàn)方法[J]. 振動、測試與診斷, 2012, 32(1): 68-72.

ZHONG Youming, ZHAO Qiang, ZHOU Jianting. Application of residual-based cumulative sum control charts for detecting faults in autocorrelated processes[J]. Journal of Vibration, Measurement & Diagnosis, 2012, 32(1): 68-72.

[13]林玉榮, 王強. 基于一維經(jīng)驗?zāi)B(tài)分解的圖像細節(jié)提取方法[J]. 吉林大學(xué)學(xué)報：工學(xué)版, 2011, 41(6): 1766-1770.

LIN Yurong, WANG Qiang. Extracting details from images based on 1-DEMD[J]. Journal of Jilin University： Engineering and Technology Edition, 2011, 41(6): 1766-1770.

[14]HARIHARAN H, KOSCHAN A, ABIDI B, et al. Fusion of visible and infrared images using empirical mode decomposition to improve face recognition[C]//Proceedings of IEEE International Conference on Image Processing. Atlanta, GA: IEEE, 2006: 2049-2052. --------------------------------------------------------------------------------------------

英文引用格式：SHAO Xin, YIN Qingbo, LU Mingyu. Vessels extraction in coronary angiogram based on empirical mode decomposition[J]. CAAI Transactions on Intelligent Systems, 2015, 10(6): 851-857.

Vessels extraction in coronary angiogram

based on empirical mode decomposition

SHAO Xin, YIN Qingbo, LU Mingyu

(College of Information Science and Technology, Dalian Maritime University, Dalian 116026, China)

Abstract：In this paper, we present a new vascular characteristics extraction method based on empirical mode decomposition. The proposed method overcomes the problems of noise and the blurred vascular edge generally encountered in the image processing procedure. This method uses empirical mode decomposition to resolve the coronary angiogram into several image layers, retaining the layers that contain blood vessel information and discarding layers that contain only background information. Clear and distinct vessel segments are then extracted from the reserved image layers, based on the connected region's properties. Finally, the coronary angiogram is reconstructed according to the information extracted from the reserved layers. The experimental results prove that this method has high accuracy and practical medical value.

Keywords：vessels extraction; empirical mode composition; coronary angiogram; characteristics extraction; connected region; interlayer reference; vascular remodeling; image layer

通信作者：邵欣.E-mail: shaoxin2123@163.com.

基金項目：國家自然科學(xué)基金資助項目(61175053, 61272369)；中央高校基本科研業(yè)務(wù)費基金主項目(2011QN 126).

收稿日期：2013-03-20.

中圖分類號：TP39

文獻標(biāo)志碼：A

文章編號：1673-4785(2015)06-0851-07

DOI:10.11992/tis.201303032