技術(shù)領(lǐng)域：

本發(fā)明涉及醫(yī)學(xué)圖像分割技術(shù)領(lǐng)域，特別涉及一種基于atlas的頭部mri圖像自動分割方法。

背景技術(shù)：
：

醫(yī)學(xué)圖像分割就是將目標(biāo)組織從圖像中提取出來，然后進(jìn)行定量、定性的分析。目前廣泛應(yīng)用于疾病診斷，術(shù)前規(guī)劃和介入手術(shù)術(shù)中導(dǎo)航領(lǐng)域中。傳統(tǒng)的醫(yī)學(xué)圖像分割多為醫(yī)生手動進(jìn)行分割，不僅耗時久而且對醫(yī)生的技術(shù)水平和經(jīng)驗有較大依賴?，F(xiàn)有的分割技術(shù)分為半自動分割，和全自動分割。半自動分割仍是以人的經(jīng)驗為主要因素，目前提出的全自動分割方法大多只適用于特定的組織分割情況，并且對于邊緣不清晰的區(qū)域可能會出現(xiàn)輪廓曲線不連續(xù)的情況。

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

本發(fā)明的目的是提供一種基于atlas的頭部mri圖像自動分割方法，建立人體頭部mri圖像的atlas圖像庫，所述的每個atlas圖像中包含一幅灰度圖像和一幅由醫(yī)生從該灰度圖中手動分割出的腦組織二值輪廓圖，根據(jù)atlas圖像庫中的灰度圖像與目標(biāo)圖像的相似性選出最佳模板，然后將模板中的二值輪廓圖經(jīng)變換后作為分割的初始輪廓曲線，最后采用主動輪廓法分割出腦組織輪廓。

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

一種基于atlas的頭部mri圖像分割方法，所述方法的具體實現(xiàn)過程為：

步驟一、建立人體頭部mri圖像的atlas圖像庫，所述的每個atlas圖像中包含一幅灰度圖像mi(v)和一幅由醫(yī)生從mi(v)中手動分割出的腦組織二值輪廓圖

步驟二、統(tǒng)計圖像庫中所有手動分割出的腦組織的灰度值取范圍；

步驟三、采用基于仿射變換的剛性配準(zhǔn)將圖像庫中的每幅atlas灰度圖mi(v)(i＝1,2…n)分別與目標(biāo)圖像f(v)進(jìn)行配準(zhǔn)，并存儲每個配準(zhǔn)的變換參數(shù)ki,θi,δxi,δyi(i＝1,2…n)；

步驟四、采用dice評價算法計算變換后的灰度圖t(mi(v))與目標(biāo)圖像f(v)的相似性，按所得dice值將atlas圖像進(jìn)行排列，將dice值最大的atlas圖像作為下一步的分割模板；

步驟五、將步驟四選出的atlas圖像中的二值輪廓圖用步驟二所得變換參數(shù)ki,θi,δxi,δyi進(jìn)行仿射變換，將變換后的圖像映射到目標(biāo)圖像中作為腦組織分割的初始輪廓；

步驟六、采用主動輪廓法分割腦組織。

在步驟三中，所述剛性配準(zhǔn)具體過程為：

1)對atlas灰度圖像進(jìn)行仿射變換

t(mi(vj))＝kirmi(vj)+δti(1)

t(mi(vj))＝[xj′,yj′]^t(2)

mi(vj)＝[xj,yj]^t(3)

δti＝[δxi,δyi]^t(5)

式中mi(vj)為灰度圖像在位置vj處的灰度值，t(mi(vj))為變換后的灰度圖像在位置vj處的灰度值，ri和δti為變換矩陣，ki,θi,δxi,δyi為配準(zhǔn)參數(shù)；

2)使用平方差相似性測度(ssd)計算變換后的灰度圖像t(mi(v))與目標(biāo)圖像f(v)的相似性，ssd計算公式如下：

ωr為目標(biāo)圖像f(v)的像素集，|ωr|表示目標(biāo)圖像f(v)中像素的總個數(shù)；

3)優(yōu)化配準(zhǔn)參數(shù)ki,θi,δxi,δyi使ssd值最大，并存儲當(dāng)前配準(zhǔn)參數(shù)。

在步驟四中，所述dice評價算法為：

式中為像素值在步驟二中所統(tǒng)計的腦組織灰度值范圍內(nèi)的像素點集合。

在步驟六中，使用主動輪廓法分割腦組織的具體過程為：

1)建立snake主動輪廓模型：

設(shè)目標(biāo)輪廓曲線為：v(s)＝(x(s),y(s))，建立能量泛函，目標(biāo)輪廓曲線為該泛函的局部極小值；

式中為圖像的內(nèi)能項，由步驟四中所得的初始輪廓定義內(nèi)能項的初值，α、β為權(quán)重參數(shù)，用于控制曲線的連續(xù)性和光滑性；eext(v(s))為目標(biāo)圖像的圖像能項，根據(jù)目標(biāo)圖像的灰度梯度信息來自動識別腦組織的輪廓邊界，并對圖像做高斯模糊處理用于去噪和擴(kuò)大捕捉范圍；econ(v(s))為外加約束能，根據(jù)應(yīng)用情況而定并沒有普遍的形式；

2)離散化初始輪廓曲線：

將封閉的二值輪廓線以等間距h離散為一系列點v0,v1,v2…vn其中(v0＝vn,vj＝(xj,yj))，曲線離散后的能量泛函為：

內(nèi)能項表達(dá)式為：

當(dāng)滿足公式(11)時，求得公式(9)的極小值；

求解公式(11)

令f(x,y)＝eext，結(jié)合公式(11)和公式(12)得到關(guān)于x，y的線性方程組：

ax+fx(x,y)＝0(13)

ay+fy(x,y)＝0(14)

假設(shè)公式(13)、公式(14)等式右邊為x、y對時間的導(dǎo)數(shù)，當(dāng)x、y穩(wěn)定時，導(dǎo)數(shù)為0，則有：

axt+fx(xt,yt)＝-xt′＝-(xt-xt-1)·γ(15)

ayt+fy(xt,yt)＝-yt′＝-(yt-yt-1)·γ(16)

式中1/γ為時間步長，t為迭代次數(shù)；

為了簡化計算，假設(shè)外力在一個時間步長內(nèi)是常值，即：

fx(xt-1,yt-1)＝fx(xt,yt)(17)

fy(xt-1,yt-1)＝fy(xt,yt)(18)

整理公式(15)、公式(16)、公式(17)與公式(18)得到關(guān)于目標(biāo)曲線x、y的迭代公式：

xt＝(a+γe)^-1(γxt-1-fx(xt-1,yt-1))(19)

yt＝(a+γe)^-1(γyt-1-fx(xt-1,yt-1))(20)

3)分割過程

a.給定參數(shù)α、β、γ；

b.輸入離散后的初始輪廓曲線；

c.初始化迭代參數(shù)t＝1；x(0),y(0)；

d.計算目標(biāo)圖像的圖像能

對目標(biāo)圖像f做高斯模糊處理，然后據(jù)目標(biāo)圖像f的灰度梯度值檢測腦組織邊緣，所述目標(biāo)圖像的圖像能計算公式如下：

e.給定誤差參數(shù)ε

用公式|(xt,yt)-(xt-2,yt-2)|≤ε作為迭代終止條件；

f.將上述步驟所得參數(shù)帶入公式(19)、公式(20)中進(jìn)行迭代計算，輸出目標(biāo)輪廓曲線。

本發(fā)明有益效果：

本發(fā)明的目的是提供一種自動的頭部圖像分割方法，只需使用者設(shè)定初始參數(shù)，分割過程完全由計算機實現(xiàn)，不需要依賴醫(yī)生的技術(shù)水平和經(jīng)驗，避免圖像分割結(jié)果受人的主觀因素影響。

本發(fā)明使用主動輪廓法分割腦組織，該方法通過驅(qū)動閉合的初始輪廓曲線向腦組織邊界變形從而得到腦組織邊緣輪廓，因此得到的邊緣輪廓曲線也是閉合的，避免了傳統(tǒng)分割方法可能出現(xiàn)的邊緣不連續(xù)問題。

本發(fā)明通過選擇與目標(biāo)圖像最為相似的atlas二值輪廓圖作為初始輪廓曲線，分割過程只需要驅(qū)動初始輪廓曲線進(jìn)行小的變形，因此有更快分割速度和更精確的分割結(jié)果。

本發(fā)明通過建立atlas圖像庫作為目標(biāo)圖像的分割模板，因此本發(fā)明一種基于atlas的頭部mri圖像分割方法可以通過擴(kuò)充所述atlas圖像庫來實現(xiàn)對更多的腦部組織分割情況，有更加廣泛的適用性。

附圖說明：

圖1是本發(fā)明一種基于atlas的頭部mri圖像分割方法的流程示意圖。

圖2是一個atlas圖像中的灰度圖和分割出的二值輪廓圖。

圖3是采用主動輪廓法分割腦組織流程示意圖。

圖4是采用主動輪廓法分割腦組織的效果圖。

具體實施方式：

如圖1所示，本實施方式所述的一種基于atlas的頭部mri圖像分割方法的具體實現(xiàn)過程為：

步驟一、建立人體頭部mri圖像的atlas圖像庫作為目標(biāo)圖像的分割模板，給定一幅目標(biāo)圖像f(v)，通過計算f(v)與圖像庫中灰度圖像的相似性選擇最佳atlas分割模板，如圖2所示，所述的每個atlas圖像中包含一幅灰度圖像mi(v)和一幅由醫(yī)生從mi(v)中手動分割出的腦組織二值輪廓圖

步驟二、統(tǒng)計圖像庫中所有手動分割出的腦組織的灰度值取范圍。

步驟三、由于目標(biāo)圖像與圖像庫中的灰度圖像可能存在空間位置，角度的差別，所以采用基于仿射變換的剛性配準(zhǔn)將圖像庫中的每幅atlas灰度圖mi(v)(i＝1,2…n)分別與目標(biāo)圖像f(v)進(jìn)行空間匹配，并存儲每個配準(zhǔn)的變換參數(shù)ki,θi,δxi,δyi(i＝1,2…n)，所述剛性配準(zhǔn)具體過程為：

1)對atlas灰度圖像進(jìn)行仿射變換

t(mi(vj))＝kirmi(vj)+δti(1)

t(mi(vj))＝[xj′,yj′]^t(2)

mi(vj)＝[xj,yj]^t(3)

δti＝[δxi,δyi]^t(5)

式中mi(vj)為灰度圖像在位置vj處的灰度值，t(mi(vj))為變換后的灰度圖像在位置vj處的灰度值，ri和δti為變換矩陣，ki,θi,δxi,δyi為配準(zhǔn)參數(shù)。

2)使用平方差相似性測度(ssd)計算變換后的灰度圖像t(mi(v))與目標(biāo)圖像f(v)的相似性，ssd計算公式如下：

ωr為目標(biāo)圖像f(v)的像素集，|ωr|表示目標(biāo)圖像f(v)中像素的總個數(shù)。

3)優(yōu)化配準(zhǔn)參數(shù)ki,θi,δxi,δyi使ssd值最大，并存儲當(dāng)前配準(zhǔn)參數(shù)。

步驟四、根據(jù)圖像庫中變換后的灰度圖t(mi(v))與目標(biāo)圖像f(v)中腦組織的相似性選擇最佳的atlas分割模板，假設(shè)t(mi(v))與f(v)中灰度值在步驟二所求的灰度值范圍內(nèi)的像素點為腦組織，采用dice評價算法計算t(mi(v))與f(v)中腦組織的相似性，并按所得dice值將atlas圖像進(jìn)行排列，將dice值最大的atlas圖像作為下一步的分割模板。

在步驟四中，所述dice評價算法為：

式中為像素值在步驟二中所統(tǒng)計的腦組織灰度值范圍內(nèi)的像素點集合。

步驟五、將步驟四選出的atlas圖像中的二值輪廓圖用步驟二所得變換參數(shù)ki,θi,δxi,δyi進(jìn)行仿射變換，將變換后的圖像映射到目標(biāo)圖像中作為腦組織分割的初始輪廓曲線。

步驟六、采用主動輪廓法分割腦組織，具體過程為：

1)建立snake主動輪廓模型：

設(shè)目標(biāo)輪廓曲線為：v(s)＝(x(s),y(s))，建立能量泛函，目標(biāo)輪廓曲線為該泛函的局部極小值：

式中為圖像的內(nèi)能項，由步驟四中所得的初始輪廓定義內(nèi)能項的初值，α、β為權(quán)重參數(shù)，用于控制曲線的連續(xù)性和光滑性；eext(v(s))為目標(biāo)圖像的圖像能項，根據(jù)目標(biāo)圖像的灰度梯度信息來自動識別腦組織的輪廓邊界，并對圖像做高斯模糊處理用于去噪和擴(kuò)大捕捉范圍；econ(v(s))為外加約束能，根據(jù)應(yīng)用情況而定并沒有普遍的形式。

2)離散化初始輪廓曲線：

將封閉的二值輪廓線以等間距h離散為一系列點v0,v1,v2…vn其中(v0＝vn,vj＝(xj,yj))，曲線離散后的能量泛函為：

內(nèi)能項表達(dá)式為：

當(dāng)滿足公式(11)時，求得公式(9)的極小值；

求解公式(11)

令f(x,y)＝eext，結(jié)合公式(11)和公式(12)得到關(guān)于x，y的線性方程組：

ax+fx(x,y)＝0(13)

ay+fy(x,y)＝0(14)

假設(shè)公式(13)、公式(14)等式右邊為x、y對時間的導(dǎo)數(shù)，當(dāng)x、y穩(wěn)定時，導(dǎo)數(shù)為0，則有：

axt+fx(xt,yt)＝-xt′＝-(xt-xt-1)·γ(15)

ayt+fy(xt,yt)＝-yt′＝-(yt-yt-1)·γ(16)

式中1/γ為時間步長，t為迭代次數(shù)；

為了簡化計算，假設(shè)外力在一個時間步長內(nèi)是常值，即：

fx(xt-1,yt-1)＝fx(xt,yt)(17)

fy(xt-1,yt-1)＝fy(xt,yt)(18)

整理公式(15)、公式(16)、公式(17)與公式(18)得到關(guān)于目標(biāo)曲線x、y的迭代公式：

xt＝(a+γe)^-1(γxt-1-fx(xt-1,yt-1))(19)

yt＝(a+γe)^-1(γyt-1-fx(xt-1,yt-1))(20)

3)如圖3所示，具體分割過程為：

a.給定參數(shù)α、β、γ；

b.輸入離散后的初始輪廓曲線；

c.初始化迭代參數(shù)t＝1；x(0),y(0)；

d.計算目標(biāo)圖像的圖像能

對目標(biāo)圖像f做高斯模糊處理，然后據(jù)目標(biāo)圖像f的灰度梯度值檢測腦組織邊緣，所述目標(biāo)圖像的圖像能計算公式如下：

e.給定誤差參數(shù)ε

用公式|(xt,yt)-(xt-2,yt-2)|≤ε作為迭代終止條件；

f.將上述步驟所得參數(shù)帶入公式(19)、公式(20)中進(jìn)行迭代計算，輸出目標(biāo)輪廓曲線。