楊 琴，姜自波，劉承超

(曲阜師范大學(xué) 信息科學(xué)與工程學(xué)院,山東 日照 276826)

基于頻域快速解卷的血管提取算法

楊 琴，姜自波，劉承超

(曲阜師范大學(xué) 信息科學(xué)與工程學(xué)院,山東 日照 276826)

為了方便識別新生血管，從而有效地對疾病進(jìn)行診斷，所以需要將血管圖像從背景圖像中分割出來。文中提出一種在提高血管分辨率的基礎(chǔ)上進(jìn)行血管分割的算法。首先通過特定公式在梯度圖像上對血管圖像進(jìn)行紋理增強(qiáng)，這是為了在之后的圖像分割中，能識別細(xì)微的血管，使其不被忽略掉；然后對模糊圖像進(jìn)行頻域上的快速解卷去模糊，消除成像儀器在成像過程中對圖像清晰度的影響；由于圖像在傳輸過程中會產(chǎn)生噪聲，因此為了去除噪聲對血管分割的影響，接著對血管圖像進(jìn)行了中值濾波操作；最后使用最大類間方差法來進(jìn)行血管的分割操作，因?yàn)樽畲箢愰g方差法可以有效地將圖像前景和背景分割開。通過實(shí)驗(yàn)對比，直觀上證明了該算法在血管分割中的有效性。

頻率域解卷；紋理增強(qiáng)；中值濾波；OTSU圖像分割

1 概 述

雖然如今用于醫(yī)學(xué)圖像成像的儀器已經(jīng)達(dá)到了很高的精度，例如光顯微鏡，它在橫向空間中限制光的衍射，從而能提供微米甚至亞微米的分辨率[1-4]，但是成像儀器也具有其局限性。比如，由于受成像設(shè)備的影響、局部體效應(yīng)(同一個(gè)體素中包含多種組織)、患者的體位運(yùn)動和檢查床的勻速直線運(yùn)動，使得醫(yī)學(xué)圖像不可避免地盤現(xiàn)噪聲和偽影，邊緣模糊和信噪強(qiáng)度不均勻。此外，在圖形成像和傳輸?shù)倪^程中，圖像的像質(zhì)也會受到一定影響。這些都給醫(yī)生下達(dá)準(zhǔn)確的診斷造成了一定的障礙。

為了提高醫(yī)學(xué)圖像的可讀性，使得醫(yī)生可以對人體的解剖結(jié)構(gòu)以及病變部位進(jìn)行更有效的觀察和診斷，提高診斷的準(zhǔn)確率，醫(yī)學(xué)圖像的后期處理成為一個(gè)必不可少的過程?，F(xiàn)如今許多疾病都伴隨著新生血管的產(chǎn)生，例如視網(wǎng)膜表面血管的變化可能暗示著糖尿病、高血壓、動脈硬化的發(fā)生。但是在疾病發(fā)生初期，血管都處于微米甚至納米層，因此，正確提取和分析血管圖像在預(yù)防和診斷病患中起到了相當(dāng)大的作用。最大類間方差法[5]是較早提出的一種圖像分割技術(shù)，但它對噪音十分敏感，并且僅對類間方差為單峰的圖像有較好的分割效果。目前，進(jìn)行血管提取的方法有很多種。例如，利用Hessian矩陣的特征值來提取血管[6-7]，該方法是一種多分辨率的提取方法。另一種進(jìn)行血管提取的方法，如文獻(xiàn)[8]中的方法，整合了多種圖像增強(qiáng)的方法來進(jìn)行血管提取，雖然得出的效果比較好，但是實(shí)驗(yàn)過程繁雜。

去卷積是一種計(jì)算密集型圖像處理技術(shù)，它試圖從模糊圖像中得到對應(yīng)的清晰圖像。通過該技術(shù)，可以提高顯微圖像的對比度和清晰度[9]。目前，卷積技術(shù)可以分為兩種：盲解卷積和非盲解卷積。盲解卷積需要事先知道模糊因子，而非盲解卷積則不需要知道模糊因子就可以從模糊圖像中恢復(fù)出清晰圖像。卷積技術(shù)的發(fā)展是值得關(guān)注的問題，目前已經(jīng)得到了各研究機(jī)構(gòu)的廣泛關(guān)注[10-12]。

文中提出了基于頻率域的非盲解卷積技術(shù)，用于提高血管的分辨率。其中的關(guān)鍵技術(shù)是在解卷之前使用了紋理增強(qiáng)技術(shù)，強(qiáng)化了血管邊緣。該技術(shù)還能使血管和背景圖像較清晰地分隔開，使得Otsu算法能夠順利進(jìn)行。為了證明文中算法的可行性，分別對腦血管和視網(wǎng)膜血管進(jìn)行了對比實(shí)驗(yàn)，結(jié)果顯示，該算法能將血管分割出來的同時(shí)保留了部分細(xì)小血管，并且對于血管邊緣的分割也比較清晰。

2 算法原理

2.1 圖像解卷

對于一幅模糊圖像，若直接對其進(jìn)行分割操作，有些細(xì)節(jié)則不能很好地分割出來，所以需要在分割之前對其進(jìn)行清晰化操作。而解卷技術(shù)可以提高顯微圖像的對比度和清晰度[1]，在頻率域上進(jìn)行解卷操作，可以大大提高操作的速度[13-14]。

假設(shè)I是原始未退化的清晰灰度圖像，g是模糊或加噪聲后的退化圖像。文中使用卷積核k對圖像I進(jìn)行卷積模糊并對其添加零均值高斯噪聲來獲得退化圖像g，用矩陣方式可表示為：g=I?k+n。假設(shè)g和k是已知的，目的就是構(gòu)建最佳的清晰圖像I。從概率方面考慮，就是為了尋找I的最大后驗(yàn)概率：p(I|g,k)。要得到p(I|g,k)的最大值，相當(dāng)于求-logp(I|g,k)的最小值：

(1)

其中：?為二維卷積操作；矩陣ds,s∈{1,2,…,5}為圖像I的一階和二階導(dǎo)數(shù):{dx,dy,dxx,dxy,dyy}；δs為正權(quán)值，用來控制正規(guī)化的強(qiáng)度；ωs為理想圖像I的過濾操作。

令τ為閾值，對于像素i，當(dāng)|dsg|i<τ時(shí)，ωs=0，否則，ωs=dsg[13]。其中

(2)

(3)

(4)

其中：K=F(k),G=F(g),Ds=F(ds),Ws=F(ωs)；F(?)表示二維離散傅里葉變換；°表示矩陣元素相乘。

令：

(5)

則式(3)可寫為：

(6)

(7)

其中，./表示矩陣元素相除。

2.2 紋理增強(qiáng)

由于直接求圖像梯度，梯度值較小的點(diǎn)在梯度圖像上看不清，所以有必要對圖像紋理做相應(yīng)處理，使其在梯度圖像上能較清楚地顯現(xiàn)出來。

圖像中梯度值較大的地方，光照強(qiáng)度值的改變也較大，需要對梯度值較大的地方進(jìn)行限制，這樣才能保證整幅圖像具有較高的對比度。目的是減小梯度值的大小，并且還能保持它們的方向不變。首先對原始圖像g(x,y)取對數(shù)，再進(jìn)行相應(yīng)的紋理處理。這是由于對數(shù)圖像的亮度和人眼所觀察到的圖像亮度很接近，并且對數(shù)域里的梯度圖像和人眼所觀察到的局部對比度很接近[15]。定義紋理增強(qiáng)后的特征圖為：

F(x,y)=L'(x,y)Θ(x,y)

(8)

其中：L'(x,y)為原始圖像g(x,y)的對數(shù)形式；Θ(x,y)為尺度函數(shù)，用于改變梯度的大小，但不會改變梯度方向。

這里，令

(9)

然而，真實(shí)世界的圖像邊界是多尺度的，使用一個(gè)多分辨率的邊緣檢測機(jī)制可以準(zhǔn)確、有效地檢測到有意義的圖像強(qiáng)度。但是，又不能在每次檢測到圖像強(qiáng)度的時(shí)候都能對其梯度值進(jìn)行相應(yīng)地增強(qiáng)或減弱，因?yàn)檫@會在邊界強(qiáng)度較強(qiáng)的地方出現(xiàn)光環(huán)效應(yīng)[15]。因此，可將一幅圖像分為多層，在每一層中使用一個(gè)梯度減弱系數(shù)。這樣，所有的梯度增強(qiáng)或減弱操作都在單個(gè)層中進(jìn)行，而每一層是針對整幅圖像進(jìn)行的，這樣就不會出現(xiàn)光環(huán)效應(yīng)。

首先，構(gòu)建一個(gè)k層的高斯金字塔，在每一層中，計(jì)算其梯度值：

在第k層中，尺度函數(shù)定義為：

(10)

正如文獻(xiàn)[15]中所說的，參數(shù)α決定哪個(gè)梯度不被改變，參數(shù)β決定哪個(gè)梯度被減小(β<1)，即梯度大的被減小，比α小的梯度被稍微放大。尺度函數(shù)Θ(x,y)使用遞歸的方法，將每層得到的尺度函數(shù)進(jìn)行逐個(gè)相乘并相加：

Θ1(x,y)=?1(x,y)

Θk+1(x,y)=Θk(x,y)+Θk(x,y)?k+1(x,y)

(11)

2.3 去 噪

在進(jìn)行血管分割之前，首先需要進(jìn)行一定程度的降噪，因?yàn)樵肼晻ρ芊指町a(chǎn)生一定程度的影響。目前用于去噪的方法有很多種。例如，均值濾波器，它在濾波過程中增加了圖像的平滑度，但也丟失了一部分細(xì)節(jié)；自適應(yīng)維納濾波，它對保留圖像邊緣和高頻部分有較大的作用，但是計(jì)算量較大；形態(tài)學(xué)濾波，它適用于圖像中的對象尺寸都比較大，且沒有細(xì)小的細(xì)節(jié)的圖像。

文中采用中值濾波器。中值濾波是基于排序統(tǒng)計(jì)理論的一種能有效抑制噪聲的非線性數(shù)字濾波技術(shù)，經(jīng)常用于去除圖像或者其他信號中的噪聲。它可以在濾除噪聲的同時(shí)，保護(hù)信號的邊緣，使之不被模糊。中值濾波的基本原理是把數(shù)字圖像或數(shù)字序列中一點(diǎn)的值用該點(diǎn)的一個(gè)鄰域中各點(diǎn)值的中值代替，讓周圍的像素值接近真實(shí)值，從而消除孤立的噪聲點(diǎn)。其步驟為：

(1)從圖像中的某個(gè)采樣窗口取出奇數(shù)個(gè)數(shù)據(jù)并從小到大排序；

(2)用排序后的中值取代要處理的數(shù)據(jù)。

2.4 血管分割

基于去模糊和去噪的血管圖像，將血管和背景圖像分割開，需要選擇合適的分割算法?，F(xiàn)如今，用于圖像分割的方法有很多種。如閾值分割，這種方法計(jì)算簡單、運(yùn)算效率高、速度快，它的關(guān)鍵技術(shù)在于閾值的選?。换趨^(qū)域的分割，區(qū)域生長技術(shù)需要選擇一組能正確代表同一區(qū)域的種子像素，確定在生長過程中的相似性準(zhǔn)則，制定讓生長停止的條件或準(zhǔn)則。這種技術(shù)對噪聲比較敏感，可能導(dǎo)致區(qū)域內(nèi)有空洞；區(qū)域分裂合并技術(shù)，它是區(qū)域生長的逆過程。這種方法對復(fù)雜圖像的分割效果較好，但算法較復(fù)雜，計(jì)算量大，分裂還可能破壞區(qū)域的邊界；邊緣檢測技術(shù)，它常用微分算子進(jìn)行邊緣檢測，但是這些算子對噪聲敏感，因此這種技術(shù)只適合于噪聲較小、不太復(fù)雜的圖像。

文中選用Otsu法(最大類間方差法)。它使用的是聚類的思想，把圖像的灰度數(shù)按灰度級分成兩個(gè)部分，使得兩個(gè)部分之間的灰度值差異最大，每個(gè)部分內(nèi)部的灰度差異最小。Otsu法通過方差的計(jì)算來尋找一個(gè)合適的灰度級別進(jìn)行分類，并且此算法不受圖像亮度和對比度的影響。前景和背景圖像的方差f由下面的過程得到：

記t為前景與背景的分割閾值，前景點(diǎn)數(shù)占圖像總點(diǎn)數(shù)的比例為w0，平均灰度為u0；背景點(diǎn)數(shù)占圖像總點(diǎn)數(shù)的比例為w1，平均灰度為u1。

圖像的總平均灰度為：

u=w0u0+w1u1

前景和背景圖像的方差為：

f=w0*(u0-u)*(u0-u)+w1*(u1-u)*(u1-u)

從N個(gè)灰度級遍歷t，使得t為某個(gè)值的時(shí)候，前景和背景的方差f(t)最大，則這個(gè)t值便是要求得的閾值。

基于頻率域解卷的血管提取方法描述如下：

數(shù)據(jù)：g為所獲得的退化圖像；k為卷積核。

步驟1：對g求對數(shù)L=log(g+1);

步驟2：求對數(shù)圖像L的梯度，利用式(9)得到梯度圖像;

步驟3：構(gòu)建高斯金字塔，利用式(10)計(jì)算每層的尺度函數(shù)?，利用式(11)計(jì)算最終的尺度函數(shù)Θ；

步驟4：利用式(8)得到最終的紋理增強(qiáng)圖像F;

步驟5：對圖像F求一階和二階導(dǎo)數(shù)；

步驟6：利用式(2)和計(jì)算得到的F的各階導(dǎo)數(shù)計(jì)算w；

步驟7：利用式(5)、(7)得到最終的解卷圖像I；

步驟8：對I進(jìn)行平滑濾波，得到去噪圖像I';

步驟9：對去噪后的圖像進(jìn)行血管分割。

3 實(shí) 驗(yàn)

文中算法流程如圖1所示。

圖1 文中算法流程

圖中，(a)為模糊的原始圖像;(b)為紋理增強(qiáng)的梯度圖像;(c)為基于頻域解卷的圖像;(d)為濾波后的圖像;(e)為血管分割圖像。

首先，對退化圖像進(jìn)行紋理增強(qiáng)，這有助于細(xì)微血管的提??；然后，對模糊圖像進(jìn)行清晰化處理，由于噪音會對血管分割產(chǎn)生一定的影響，所以有必要在提取血管前進(jìn)行濾波操作。

3.1 實(shí)驗(yàn)結(jié)果

將Horacio等提出的圖像恢復(fù)算法直接應(yīng)用到血管提取中，會發(fā)現(xiàn)有些血管被忽略了，如圖2所示。

圖2 腦血管分割

可以發(fā)現(xiàn)，在圖(c)中，只分割出血管主干，某些細(xì)微血管則被認(rèn)為是背景忽略掉。在圖(d)中，某些灰度值較小的血管沒有能分割出來，出現(xiàn)了血管的不連續(xù)性。對于文中算法，它在將血管主干提取出來的同時(shí)，仍然保留了某些細(xì)微血管，這對血管疾病的診斷是非常有必要的。

為了進(jìn)一步說明文中算法的優(yōu)越性，使用了另一幅血管圖像進(jìn)行了對比實(shí)驗(yàn)，如圖3所示。

可以發(fā)現(xiàn)，使用Otsu算法可以將主干血管分割出來，但是細(xì)微血管卻被忽略甚至多根血管被重疊為一根血管。而圖(d)的血管邊緣模糊，這對判斷血管粗細(xì)程度造成了影響。通過對比發(fā)現(xiàn)，文中算法分割出來的血管紋理清晰。

3.2 參數(shù)選擇

利用文獻(xiàn)[8]中所使用的19×19卷積核進(jìn)行實(shí)驗(yàn)。針對腦血管圖像，設(shè)置α等于0.000 05倍的平均梯度，β設(shè)為0.9時(shí)紋理增強(qiáng)效果較好。解卷過程中，對于所有的s，將λ均設(shè)為0.000 05，對于一階導(dǎo)數(shù){dx,dy}，設(shè)置τ(1)=0.065，對于二階導(dǎo)數(shù){dxx,dxy,dyy}，設(shè)置τ(2)=0.0325時(shí)，實(shí)驗(yàn)結(jié)果較好。針對視網(wǎng)膜毛細(xì)血管，將λ設(shè)為0.000 01，其他參數(shù)不變。

圖3 視網(wǎng)膜毛細(xì)血管分割

4 結(jié)束語

文中從造成圖像不清晰的因素考慮，提出了基于頻率域解卷的血管提取方法。同時(shí)考慮到Otsu算法進(jìn)行圖像分割的局限性，所以在使用Otsu算法進(jìn)行血管分割之前，先將血管紋理進(jìn)行了加強(qiáng)，將背景圖像進(jìn)行了一定程度的灰度值變化，便于后面運(yùn)用Otsu算法對血管進(jìn)行分割。紋理增強(qiáng)的關(guān)鍵在于α和β的選取，它們控制背景和前景血管的亮度對比。實(shí)驗(yàn)結(jié)果表明，在進(jìn)行血管分割之前，對圖像進(jìn)行解卷以提高圖像的分辨率是有必要的。

該算法的局限性在于：只能對背景色相對血管較暗的血管圖片進(jìn)行血管提取。因此，適用性更廣、精確度更高的血管提取算法是進(jìn)一步研究的目標(biāo)。

[1]MaR,S?ntgesS,ShohamS,etal.Fastscanningcoaxialoptoacousticmicroscopy[J].BiomedicalOpticsExpress,2012,3(7):1724-1731.

[2]MaslovK,ZhangHF,HuS,etal.Optical-resolutionphotoacousticmicroscopyforinvivoimagingofsinglecapillaries[J].OpticsLetters,2008,33(9):929-931.

[3]SongL,MaslovK,WangLV.Multifocaloptical-resolutionphotoacousticmicroscopyinvivo[J].OpticsLetters,2011,36(7):1236-1238.

[4]HajirezaP,ShiW,ZempR.Label-freeinvivoGRIN-lensopticalresolutionphotoacousticmicro-endoscopy[J].LaserPhysicsLetters,2013,10(5):055603.

[5]OtsuN.Athresholdselectionmethodfromgray-levelhistograms[J].IEEETransonSystems,ManandCybernetics,1979,9(1):62-66.

[6]SatoY,ShiragaN,ShiragaN,etal.Three-dimensionalmulti-scalelinefilterforsegmentationandvisualizationofcurvilinearstructuresinmedicalimages[J].MedicalImageAnalysis,1998,2(2):143-168.

[7]FrangiAF,NiessenWJ,VinckenKL,etal.Multiscalevesselenhancementfiltering[C]//ProcofMICCAI.[s.l.]:[s.n.],1998:130-137.

[8]YangZY.Multi-parametricquantitativemicrovascularimagingwithoptical-resolutionphotoacousticmicroscopyinvivo[J].OpticsExpress,2014,22(2):1500-1511.

[9]ChenJH,LinRQ,WangHN.Blind-deconvolutionoptical-resolutionphotoacousticmicroscopyinvivo[J].OpticsExpress,2013,21(6):7316-7327.

[10]BanhamM,KatsaggelosA.Digitalimagerestoration[J].IEEESignalProcessingMagzine,1997,14(2):24-41.

[11]CampisiP,EgiazarianK.Blindimagedeconvolution:theoryandapplications[M].BocaRaton:CRCPress,2007.

[12]LevinA,WeissY,DurandFW.Efficientmarginallikelihoodoptimizationinblinddeconvolution[C]//ProcofCVPR.[s.l.]:[s.n.],2011:2657-2664.

[13]FortunatoHE,OliveiraMM.Fasthigh-qualitynon-blinddeconvolutionusingsparseadaptivepriors[J].VisualComputer,2014,30(6):661-671.

[14]DilipK,RobF.Fastimagedeconvolutionusinghyper-laplacianpriors[C]//ProcofNIPS.[s.l.]:[s.n.],2009:1033-1041.

[15]FattalR,LischinskiD,WermanM.Gradientdomainhighdynamicrangecompression[J].ACMTransonGraphics,2002,21(3):249-256.

Vessel Extraction Algorithm Based on Fast Frequency-domain Deconvolution

YANG Qin,JIANG Zi-bo,LIU Cheng-chao

(College of Information Science and Engineering,Qufu Normal University,Rizhao 276826,China)

In order to identify new blood vessels,thus effectively diagnosing disease,it is necessary to segment the vascular image from the background image.A segmentation algorithm based on improved resolution of the vessel is proposed.First the texture is enhanced for blood vessel image through specific formula in the gradient image,in order to identify small blood vessels and not ignore it in the image segmentation later.Secondly,the image deconvolution in frequency domain is made to obtain a sharp image,eliminating the influence on image clarity of imaging instrument in the imaging process.Then so as to remove the influence of noise on the vessel segmentation,a filtering operation is carried on vascular image.Finally,the maximum class variance method is used to segment vessels,because it can effectively distinguish the image foreground and background.The effectiveness of the algorithm in vessel segmentation is verified intuitively by experimental comparisons.

frequency domain deconvolution;texture enhancement;median filtering;OTSU image segmentation

2015-06-16

2015-09-18

時(shí)間：2016-02-18

國家自然科學(xué)基金資助項(xiàng)目(61572283)；山東省自然科學(xué)基金(ZR2011FQ026)

楊 琴(1990-)，女，碩士研究生，CCF會員，研究方向?yàn)閳D像處理與模式識別。

http://www.cnki.net/kcms/detail/61.1450.TP.20160218.1634.050.html

TP391

A

1673-629X(2016)03-0113-04

10.3969/j.issn.1673-629X.2016.03.027