本發(fā)明涉及一種圖像分割技術(shù)，具體為一種基于swi圖像的腦靜脈血管分割方法。

背景技術(shù)：

圖像分割就是把圖像分割成若干特定的，具有獨(dú)特性質(zhì)的區(qū)域并提出感興趣目標(biāo)的技術(shù)和過(guò)程，swi((susceptibilityweightedimaging磁敏感加權(quán)成像)是利用相位圖的磁敏感信息對(duì)幅值圖進(jìn)行加權(quán)運(yùn)算來(lái)獲得組織和器官的解剖結(jié)構(gòu)、功能結(jié)構(gòu)和病變狀況。swi腦靜脈血管分割就是利用腦部的swi圖像來(lái)分割腦部靜脈血管。

hessian矩陣的血管增強(qiáng)多尺度濾波器是frangi等人提出的用來(lái)增強(qiáng)線性結(jié)構(gòu)的技術(shù)，由于血管拓?fù)浣Y(jié)構(gòu)呈線性的樹(shù)狀結(jié)構(gòu)，所以其特性通常用線性來(lái)表示。hessian矩陣的基本原理就是求取圖像每一個(gè)像素點(diǎn)的二階偏導(dǎo)矩陣，并通過(guò)hessian矩陣求取特征值和特征向量，根據(jù)特征值大小和檢測(cè)結(jié)構(gòu)之間的關(guān)系來(lái)判斷像素所屬的結(jié)構(gòu)。

matlab是專門(mén)用于算法開(kāi)發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計(jì)算的高級(jí)技術(shù)計(jì)算語(yǔ)言及交互式環(huán)境的開(kāi)發(fā)軟件，能夠高效快速的進(jìn)行圖像相關(guān)的處理計(jì)算。

由于圖像中腦靜脈血管復(fù)雜稠密形狀狹小、灰度大小不均勻且變化較大，圖像對(duì)比度低，并且由于偽影等噪聲影響造成實(shí)際中腦靜脈血管分割困難，目前分割方法主要是基于統(tǒng)計(jì)學(xué)、形狀模型以及尺度空間變換的分割方法，這些方法在管狀結(jié)構(gòu)分割方面效果顯著，但實(shí)際應(yīng)用中不能很好地適應(yīng)腦靜脈血管這種復(fù)雜結(jié)構(gòu)的圖像，不能很好地考慮血管圖像整體特征以及圖像的特異性，難以處理偽影等噪聲影響。

技術(shù)實(shí)現(xiàn)要素：

針對(duì)現(xiàn)有技術(shù)中hessian矩陣血管增強(qiáng)濾波器對(duì)于低對(duì)比度的swi腦靜脈血管圖像分割效果較差、偽影等噪聲影響較大的缺陷,本發(fā)明要解決的問(wèn)題是提供一種可在差異性較大的不同swi腦圖像中以及噪聲和干擾情況下都能夠獲得穩(wěn)的結(jié)果的基于swi圖像的腦靜脈血管分割方法。

為解決上述技術(shù)問(wèn)題，本發(fā)明采用的技術(shù)方案是：

本發(fā)明一種基于swi圖像的腦靜脈血管分割方法，包括以下步驟：

1)通過(guò)磁共振設(shè)備讀取每一幅二維swi腦靜脈血管圖像，調(diào)整圖像分辨率，并通過(guò)算法去除掉圖像中覆蓋在腦四周的頭骨和皮膚。

2)將上述處理后的圖像進(jìn)行各向異性濾波增強(qiáng)處理；

3)將增強(qiáng)處理后的圖像進(jìn)行改進(jìn)的2dhessian矩陣濾波增強(qiáng)處理，進(jìn)一步分割血管結(jié)構(gòu)；

4)通過(guò)面積法去掉因?yàn)槟X偽影造成的大面積噪聲，分割出大部分的細(xì)小靜脈血管，并去除腦圖像的邊界輪廓；

5)通過(guò)isdata閾值法保留因?yàn)槊娣e法導(dǎo)致錯(cuò)誤去掉的細(xì)小靜脈血管；

6)計(jì)算分割后的腦靜脈血管圖像的dsc、ppv、sensitivity以及kappa值，對(duì)分割效果差的區(qū)域用區(qū)域生長(zhǎng)法進(jìn)行優(yōu)化。

步驟1)中，通過(guò)算法去除掉圖像中覆蓋在腦四周的頭骨和皮膚包括以下步驟：

1.1)將swi圖像二維數(shù)據(jù)讀到內(nèi)存中，并調(diào)整分辨率；

1.2)通過(guò)ostu動(dòng)態(tài)閾值法計(jì)算swi圖像的閾值，對(duì)swi圖像進(jìn)行二值化處理，得到二值圖像；

1.3)對(duì)上述二值圖像的各連通區(qū)域進(jìn)行標(biāo)記，并求取各連通區(qū)域大小；

1.4)求最大連通區(qū)域的索引，并獲取最大連通區(qū)域的圖像。

步驟2)中，通過(guò)各向異性濾波增強(qiáng)算法進(jìn)行圖像增強(qiáng)包括以下步驟：

2.1)求當(dāng)前圖像像素的散度，對(duì)二維圖像中x、-x、y以及-y四個(gè)方向分別求偏導(dǎo)數(shù)，得到局部不同方向上的變化量；

2.2)求x、-x、y以及-y四個(gè)方向的導(dǎo)熱系數(shù)，變化越多方向值越小，將變化大的方向的點(diǎn)作為邊界保留，從而加強(qiáng)管狀結(jié)構(gòu)；

2.3)迭代計(jì)算，對(duì)圖像進(jìn)行迭代增強(qiáng)。

步驟2.1)中，對(duì)x、-x、y以及-y分別求偏導(dǎo)數(shù)為：

ni＝img(p-1,q)-img(p,q)；

si＝img(p+1,q)-img(p,q)；

ei＝img(p,q-1)-img(p,q)；

wi＝img(p,q+1)-img(p,q)；

其中，ni為像素點(diǎn)在y方向的變化量，si為像素點(diǎn)在-y方向的變化量,ei為像素點(diǎn)在x方向的變化量，wi為像素點(diǎn)在-x方向的變化量，img代表該圖像，p、q分別為圖像像素點(diǎn)的橫縱坐標(biāo)。

步驟2.2)中，求四個(gè)方向的導(dǎo)熱系數(shù)通過(guò)以下公式計(jì)算：

nd＝exp(-ni^2/(k*k))；

sd＝exp(-si^2/(k*k))；

ed＝exp(-ei^2/(k*k))；

wd＝exp(-wi^2/(k*k))；

其中，nd為像素點(diǎn)y方向的導(dǎo)熱系數(shù)，ni為像素點(diǎn)在y方向的變化量，sd為像素點(diǎn)-y方向的導(dǎo)熱系數(shù)，si為像素點(diǎn)在-y方向的變化量,ed為像素點(diǎn)x方向的導(dǎo)熱系數(shù)，ei為像素點(diǎn)在x方向的變化量，wd為像素點(diǎn)-x方向的導(dǎo)熱系數(shù)，wi為像素點(diǎn)在-x方向的變化量，k為導(dǎo)熱系數(shù)常量；

做為邊界的權(quán)重，四個(gè)方向的變化量越大，求得的值越小，從而達(dá)到保留邊界的目的。

步驟3)中，將增強(qiáng)后的圖像通過(guò)改進(jìn)的2dhessian矩陣濾波增強(qiáng)處理包括以下步驟：

3.1)修改2dhessian矩陣濾波器的血管增強(qiáng)響應(yīng)函數(shù)公式，使響應(yīng)函數(shù)n階可導(dǎo)，從而保證圖像平滑；

3.2)計(jì)算不同血管尺度下，血管增強(qiáng)響應(yīng)函數(shù)的最大值。

步驟3.1)中修改2dhessian矩陣濾波器的血管增強(qiáng)響應(yīng)函數(shù)公式，使響應(yīng)函數(shù)n階可導(dǎo)，從而保證圖像平滑包括以下步驟：

3.1.1)像素點(diǎn)r周?chē)徲虻奶├照归_(kāi)計(jì)算公式為：

其中，▽i(a)為梯度，▽²i(a)為點(diǎn)a處的hessian矩陣，a為像素點(diǎn)r鄰域內(nèi)的一點(diǎn)；

3.1.2)hessian矩陣計(jì)算公式為：

hessian矩陣由點(diǎn)的二階偏導(dǎo)數(shù)構(gòu)成，在二維空間體素點(diǎn)r＝(x,y)處，可以得到i(r)的二階偏導(dǎo)矩陣：

3.1.3)hessian矩陣的響應(yīng)函數(shù)計(jì)算公式為：

hessian矩陣的響應(yīng)函數(shù)為：

其中，rb為特征值λ1、λ2的比值，β和c分別為rb和s影響的權(quán)重系數(shù)，λ1、λ1表示hessian矩陣的特征值；s用來(lái)在背景區(qū)域中提取感興趣區(qū)域；rb將線性結(jié)構(gòu)從團(tuán)狀結(jié)構(gòu)中區(qū)分出來(lái)；hσ為特征值的平方和的根；

步驟3.1.4：血管響應(yīng)函數(shù)的改進(jìn)公式：

其中，c1為影響因子。

步驟4)中，通過(guò)面積法去掉因?yàn)槟X偽影造成的大面積噪聲，分割出大部分的細(xì)小靜脈血管，并去除腦圖像的邊界輪廓，其具體步驟如下：

4.1)對(duì)濾波增強(qiáng)后的圖像二值化，對(duì)各連通區(qū)域進(jìn)行標(biāo)記并求得個(gè)各連通區(qū)域的大?。?/p>

4.2)將連通區(qū)域變化量最大的值作為分割條件，分離出大于該值得區(qū)域作為偽影噪聲的備選區(qū)域；

4.3)將去除噪聲的二值區(qū)域的邊界輪廓取反后與2dhessian矩陣濾波增強(qiáng)后的圖像相乘，去掉被錯(cuò)誤識(shí)別成血管的腦圖像的邊界。

步驟6)中，計(jì)算分割后的腦靜脈血管圖像的dsc,ppv,sensitivity,kappa值，對(duì)分割效果差的區(qū)域用區(qū)域生長(zhǎng)法進(jìn)行優(yōu)化，其具體步驟如下：

6.1)計(jì)算分割后圖像的dsc,ppv,sensitivity,kappa值；

6.2)區(qū)域生長(zhǎng)算法進(jìn)行分割優(yōu)化。

步驟6.1)中，計(jì)算分割后的腦靜脈血管圖像的dsc,ppv,sensitivity,kappa值，分別為：

6.1.1)計(jì)算圖像的dsc值：

dsc＝2tp/(fp+2tp+fn)

其中，dsc為相似性系數(shù)，fp為假陽(yáng)性，tp為真陽(yáng)性，fn為假陰性；

6.1.2)計(jì)算圖像的ppv值；

ppv＝tp/(tp+fp)

其中，ppv為陽(yáng)性結(jié)果的比例；

6.1.3)計(jì)算圖像的sensitivity值；

sensitivity＝tp/(tp+fn)

其中sensitivity表示的是敏感度；

6.1.4：計(jì)算圖像的kappa值；

kappa表示的是計(jì)算分類(lèi)精度的方法，其計(jì)算公式為：

kappa＝(p0-pc)/(1-pc)

p0＝s/n；

pc＝(a1*b1+a0*b0)/(n*n)；

其中，n為像素總數(shù)，地面真值圖像中像素為1表示為a1，像素為0表示為a0，分割圖像中像素為1表示為b1，像素為0表示為b0，對(duì)應(yīng)像素相等的值數(shù)量為s；

kappa值計(jì)算結(jié)果為-1到1，值越大一致性越高。

本發(fā)明具有以下有益效果及優(yōu)點(diǎn)：

1.本發(fā)明方法能夠很好的去除圖像中因磁場(chǎng)不均勻造成的高斯噪聲以及腦間隙、生理運(yùn)動(dòng)等造成的偽影噪聲，能夠很好的抑制板狀結(jié)構(gòu)、團(tuán)狀結(jié)構(gòu)，加強(qiáng)血管的管狀結(jié)構(gòu)，能從灰度一致性較差的背景低對(duì)比度區(qū)域中提取到感興趣區(qū)域，正確率高，魯棒性好，在差異性較大的不同swi腦圖像中以及噪聲和干擾情況下都能夠獲得穩(wěn)定的結(jié)果，假陽(yáng)性率低，避免分割出不屬于靜脈的結(jié)構(gòu)，符合醫(yī)學(xué)圖像對(duì)安全性的要求。

2.發(fā)明方法消除了因?yàn)殪o脈血管低對(duì)比度以及腦偽影等噪聲等對(duì)血管分割的影響，并有效分割出了腦靜脈血管結(jié)構(gòu)。

附圖說(shuō)明

圖1為本發(fā)明一種基于swi圖像的腦靜脈血管分割方法的模塊流程圖；

圖2為本發(fā)明一種基于swi圖像的腦靜脈血管分割方法的詳細(xì)流程圖。

具體實(shí)施方式

下面結(jié)合說(shuō)明書(shū)附圖對(duì)本發(fā)明作進(jìn)一步闡述。

本發(fā)明一種基于swi圖像的腦靜脈血管分割方法，運(yùn)行在intel內(nèi)核的windows7系統(tǒng)環(huán)境中，使用matlab軟件進(jìn)行圖像分割處理。如圖1所示是血管分割三部分，包括預(yù)處理部分、分割部分以及優(yōu)化部分。圖像分割使用基于hessian矩陣的血管增強(qiáng)濾波的核心算法進(jìn)行初步的血管分割，通過(guò)各向異性濾波增強(qiáng)技術(shù)進(jìn)行預(yù)處理增強(qiáng)管狀結(jié)構(gòu)，通過(guò)面積法和isdata閾值法進(jìn)行后期優(yōu)化處理得到精確地腦靜脈血管分割圖像。

如圖2所示，本發(fā)明一種基于swi圖像的腦靜脈血管分割方法包括以下步驟：

1)通過(guò)磁共振設(shè)備讀取每一幅二維swi腦靜脈血管圖像，調(diào)整圖像分辨率，并通過(guò)算法去除掉圖像中覆蓋在腦四周的頭骨和皮膚，具體為：

1.1)將二維數(shù)據(jù)讀到內(nèi)存中，并調(diào)整分辨率，啟動(dòng)matlab軟件，調(diào)用dicomread()函數(shù)讀取dicom圖像，保存圖像數(shù)據(jù)矩陣并調(diào)用ngray()函數(shù)調(diào)整圖像分辨率；

1.2)對(duì)圖像進(jìn)行二值化處理，通過(guò)ostu動(dòng)態(tài)閾值法計(jì)算閾值，調(diào)用ostu函數(shù)求取閾值；

1.3)對(duì)上述二值圖像的各連通區(qū)域進(jìn)行標(biāo)記，并求取各連通區(qū)域大小，調(diào)用bwlabel()函數(shù)；

1.4)求最大連通區(qū)域的索引，并獲取最大連通區(qū)域的圖像；

本實(shí)施例中，求各連通區(qū)域的大小，調(diào)用regionprops(imlabel,’area’)，將不同連通區(qū)域大小記錄在stats結(jié)構(gòu)體中；

求最大連通區(qū)域的索引，調(diào)用find函數(shù)找到連通區(qū)域中最大的區(qū)域；

獲取最大連通區(qū)域，調(diào)用ismenber()函數(shù)，獲得最大的連通區(qū)域圖像。

2)將步驟1處理后的圖像通過(guò)各向異性濾波增強(qiáng)算法進(jìn)行圖像增強(qiáng)，包括以下步驟：

2.1)求當(dāng)前像素的散度，對(duì)四個(gè)方向分別求偏導(dǎo)數(shù)，求得局部不同方向上的變化量:

圖像不同方向的偏導(dǎo)數(shù)和變化量的計(jì)算公式為：

ni＝img(p-1,q)-img(p,q)；

si＝img(p+1,q)-img(p,q)；

ei＝img(p,q-1)-img(p,q)；

wi＝img(p,q+1)-img(p,q)；

ni為像素點(diǎn)在y方向的變化量，si為像素點(diǎn)在-y方向的變化量,ei為像素點(diǎn)在x方向的變化量，wi為像素點(diǎn)在-x方向的變化量，img代表該圖像，p、q分別為圖像像素點(diǎn)的橫、縱坐標(biāo)；

2.2)求四個(gè)方向的導(dǎo)熱系數(shù)，變化越多的方向值越小，作為邊界保留，從而加強(qiáng)管狀區(qū)域；

四個(gè)方向的導(dǎo)熱系數(shù)計(jì)算公式為：

nd＝exp(-ni^2/(k*k))；

sd＝exp(-si^2/(k*k))；

ed＝exp(-ei^2/(k*k))；

wd＝exp(-wi^2/(k*k))；

其中，nd為像素點(diǎn)y方向的導(dǎo)熱系數(shù)，ni為像素點(diǎn)在y方向的變化量，sd為像素點(diǎn)-y方向的導(dǎo)熱系數(shù)，si為像素點(diǎn)在-y方向的變化量,ed為像素點(diǎn)x方向的導(dǎo)熱系數(shù)，ei為像素點(diǎn)在x方向的變化量，wd為像素點(diǎn)-x方向的導(dǎo)熱系數(shù)，wi為像素點(diǎn)在-x方向的變化量，k為導(dǎo)熱系數(shù)常量；做為邊界的權(quán)重，四個(gè)方向的變化量越大，求得的值越小，從而達(dá)到保留邊界的目的。

2.3)迭代計(jì)算，對(duì)圖像進(jìn)行迭代增強(qiáng)。

3)將增強(qiáng)后的圖像進(jìn)行改進(jìn)的2dhessian矩陣濾波增強(qiáng)處理，進(jìn)一步分割血管結(jié)構(gòu)，包括：

3.1)將加強(qiáng)后的圖像進(jìn)行2dhessian矩陣濾波增強(qiáng)處理；

3.1.1)像素點(diǎn)r周?chē)徲虻奶├照归_(kāi)計(jì)算公式為：

其中，▽i(a)為點(diǎn)a的梯度，▽²i(a)為點(diǎn)a處的二階偏導(dǎo)數(shù)，a為像素點(diǎn)r鄰域內(nèi)的一點(diǎn)

3.1.2)hessian矩陣計(jì)算公式為：

hessian矩陣由該點(diǎn)的二階偏導(dǎo)數(shù)構(gòu)成，在二維空間體素點(diǎn)r＝(x,y)處，得到i(r)的二階偏導(dǎo)矩陣：

用λ1，λ2表示hessian矩陣的特征值，e1，e2表示相應(yīng)的特征向量，i(r)為圖像i中的像素點(diǎn)r。

表1：hessian矩陣特征值和檢測(cè)出的結(jié)構(gòu)之間的關(guān)系表

表格顯示了frangi的hessian矩陣?yán)碚撝笑舜笮∨袛嘌芙Y(jié)構(gòu)。

3.1.3)hessian矩陣的響應(yīng)函數(shù)計(jì)算公式為：

hessian矩陣的響應(yīng)函數(shù)為：

其中，rb為特征值λ1、λ2的比值β和c分別表示的是λ2和s影響的權(quán)重系數(shù)。rb將線性結(jié)構(gòu)從團(tuán)狀結(jié)構(gòu)中區(qū)分出來(lái)，s用來(lái)在背景區(qū)域中提取感興趣區(qū)域，hσ為特征值的平方和的根；

3.1.4)血管響應(yīng)函數(shù)的改進(jìn)公式：

c1為影響因子，改進(jìn)的公式是為了保證響應(yīng)函數(shù)連續(xù)，n階可導(dǎo)，從而使圖像平滑。

4)通過(guò)面積法去掉因?yàn)槟X偽影造成的大面積噪聲，保留大部分的細(xì)小靜脈血管，并去除血管圖像邊界，具體為：

4.1)通過(guò)ostu閾值法求取閾值，對(duì)濾波后的圖像二值化，通過(guò)bwlabel()函數(shù)對(duì)各連通區(qū)域進(jìn)行標(biāo)記并通過(guò)regionprops()函數(shù)求得個(gè)各連通區(qū)域的大小，并將數(shù)據(jù)保存到stats數(shù)據(jù)結(jié)構(gòu)中；

4.2)根據(jù)連通區(qū)域變化量最大的值作為分割條件，分離出大于該值得區(qū)域作為偽影噪聲的備選區(qū)域；

4.3)對(duì)去除噪聲的區(qū)域通過(guò)imfill()函數(shù)進(jìn)行填洞，并求取bwperim()函數(shù)求取邊界輪廓，將該輪廓取反與2dhessian矩陣濾波增強(qiáng)函數(shù)frangifilter2d()后的圖像相乘，去掉被錯(cuò)誤識(shí)別成血管的圖像的邊界。

5)通過(guò)isdata閾值法保留因?yàn)槊娣e法而錯(cuò)誤去掉的細(xì)小靜脈血管，具體為：

5.1)通過(guò)isdata()閾值法求得2dhessian矩陣濾波增強(qiáng)的腦靜脈血管閾值；

5.2)將面積法過(guò)濾的偽影噪聲圖像通過(guò)閾值分割，分出較大面積的被錯(cuò)誤分割為偽影的血管區(qū)域；

5.3)將上一步分離的血管區(qū)域疊加到步驟4.3所處理后的圖像中，得到最終的血管分割結(jié)果。

6)計(jì)算分割圖像的dsc(dicesimilaritycoefficient),ppv(positivepredictivevalue),sensitivity,kappa值，對(duì)分割效果差的區(qū)域用區(qū)域生長(zhǎng)法進(jìn)行優(yōu)化，具體為：

6.1)計(jì)算分割后圖像的dsc(dicesimilaritycoefficient，相似性系數(shù))、ppv(positivepredictivevalue，敏感度)、sensitivity(敏感度)以及kappa(一致性評(píng)價(jià)系數(shù))值，具體為：

6.1.1)計(jì)算圖像的dsc值，其計(jì)算公式為：

dsc＝2tp/(fp+2tp+fn)

其中，dsc為相似性系數(shù)，fp為假陽(yáng)性，tp為真陽(yáng)性，fn為假陰性。

6.1.2)計(jì)算圖像的ppv值，計(jì)算公式為：

ppv＝tp/(tp+fp)

其中，ppv表示的是陽(yáng)性結(jié)果的比例；

6.1.3)計(jì)算圖像的sensitivity值，具體計(jì)算公式為：

sensitivity＝tp/(tp+fn)

其中sensitivity表示的是敏感度；

6.1.4)計(jì)算圖像的kappa值；

kappa為一致性評(píng)價(jià)系數(shù)，用于計(jì)算分類(lèi)精度的方法中，其計(jì)算公式為：

kappa＝(p0-pc)/(1-pc)

p0＝s/n；

pc＝(a1*b1+a0*b0)/(n*n)；

其中，n為像素總數(shù)，地面真值圖像中像素為1表示a1，像素為0表示為a0，分割圖像中像素為1表示為b1，像素為0表示為b0，對(duì)應(yīng)像素相等的值數(shù)量為s；

kappa值計(jì)算結(jié)果為-1到1，值越大一致性越高。

6.2)區(qū)域生長(zhǎng)算法進(jìn)行分割優(yōu)化，具體為：

6.2.1)通過(guò)ginput()函數(shù)獲取區(qū)域生長(zhǎng)起始種子點(diǎn)；

6.2.2)設(shè)定閾值為0.5，將種子點(diǎn)八鄰域內(nèi)符合閾值范圍內(nèi)像素保存。

本發(fā)明通過(guò)dsc、ppv、sensitivity以及kappa值來(lái)定量分析分割效果，如果分割效果達(dá)到期望效果，則輸出分割圖像，否則對(duì)區(qū)域長(zhǎng)進(jìn)行迭代優(yōu)化，直到達(dá)到期望效果。

通過(guò)應(yīng)用本發(fā)明方法，消除了因?yàn)殪o脈血管低對(duì)比度以及腦偽影等噪聲等對(duì)血管分割的影響，并有效分割出了腦靜脈血管結(jié)構(gòu)。

以上所述僅為本發(fā)明的優(yōu)選實(shí)施例，并不用于限制本發(fā)明，對(duì)于本領(lǐng)域的技術(shù)人員來(lái)說(shuō)，依然可以對(duì)前述實(shí)施例中的技術(shù)方案進(jìn)行修改，或?qū)Σ糠旨夹g(shù)進(jìn)行等效替換，而這些修改或替換都應(yīng)屬于本發(fā)明所附權(quán)利要求的保護(hù)范圍。