專利名稱：基于隱含活動輪廓先驗的貝葉斯圖像重建方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種圖像重建方法，尤其涉及一種基于隱含活動輪廓先驗的貝葉斯圖像重建方法。
背景技術(shù)：
正電子發(fā)射計算機斷層顯像(Positron emission tomography，PET)重建方法的研究近年來越來越受到人們的重視?，F(xiàn)有的圖像重建方法主要為解析法和迭代法。解析法主要是以中心切片定理為基礎(chǔ)的濾波反投影方法，其特點是速度快，但當(dāng)測量噪聲較大或采樣不充分時，這類方法的成像效果不理想。迭代法主要包括代數(shù)法，最大似然法等等。迭代法的特點是可以根據(jù)具體成像條件引入與空間幾何有關(guān)的或與測量值大小有關(guān)的約束和條件因子，但迭代法收斂速度慢，運算時間長。正則化技術(shù)或者最大后驗重建(Maximum A Posteriori，MAP)方法的使用是必要的。最大后驗重建方法對于提高圖像的空間分辨率和噪聲特性的效果是非常顯著的，但經(jīng)常會出現(xiàn)過平滑現(xiàn)象，也就是圖像在被平滑的同時，邊緣也被平滑掉了，因為噪聲和邊緣都是高頻分量。而且使用這種方法的前提是要選擇一種合適的先驗，不恰當(dāng)?shù)南闰灧植紩?dǎo)致完全錯誤的重建結(jié)果。由于先驗函數(shù)中存在超參數(shù)估計問題，有文獻表示用其它的高質(zhì)量成像方法如CT、MRI所提供的人體形態(tài)學(xué)結(jié)構(gòu)圖像來構(gòu)造一種分割模板先驗，以保證先驗的可靠性，重建結(jié)果采用基于動態(tài)模擬來計算，這種方法可以很好解決超參數(shù)的選擇，但這種方法也必須事先知道CT、MRI重構(gòu)圖像，在通常情況下要病人同時進行PET，CT或MRI掃描會增加病人的痛苦，因此這種方法是不實用的。

發(fā)明內(nèi)容
本發(fā)明提供一種重建后的圖像即可以保持高的分辨率，又能夠降低噪聲特性且重建后圖像的邊緣清晰的基于隱含活動輪廓先驗的貝葉斯圖像重建方法。
本發(fā)明采用如下技術(shù)方案一種基于隱含活動輪廓先驗的貝葉斯圖像重建方法
1)獲取投影數(shù)據(jù)，根據(jù)待重建圖像的尺寸要求，確定初始圖像范圍，給定初始的灰度值大于1，并將2維圖像變換成1維向量進行計算，2)根據(jù)投影數(shù)據(jù)規(guī)模和要求成像圖像x的大小，計算系統(tǒng)概率矩陣P，3)將系統(tǒng)概率陣和初始圖像x相乘，得到前向投影a，4)將投影數(shù)據(jù)除以前向投影a，得到對投影數(shù)據(jù)的校正值c，5)將系統(tǒng)概率矩陣P乘以投影數(shù)據(jù)的校正值c，得到圖像成像迭代過程中的修正值xd，6)將高斯函數(shù)對初始圖像作卷積計算，然后求這個計算結(jié)果的梯度模，再將這個計算結(jié)果除以0～1之間的任意一個數(shù)β1的平方，得到第一個數(shù)e，將第一個數(shù)e加1的結(jié)果的倒數(shù)乘以圖像的梯度模|x|，得到一個擴散速度g，將圖像的梯度x除以圖像的梯度模|x|，并對這個結(jié)果求散度，這個散度和氣球力F0之和乘以上述擴散速度g，得到能量函數(shù) 用差分方法將這個能量函數(shù)離散化，7)將系統(tǒng)概率矩陣P對它的每一行求和，把這個對行的和的β倍加上能量函數(shù)，得到一個權(quán)值，將初始圖像乘以修正值xd除以這個權(quán)值，得到了一個重建圖像，再將此重建圖像作為下一次迭代的初始圖像，返回到第3步，反復(fù)迭代直到重建后的圖像收斂。
與現(xiàn)有技術(shù)相比，本發(fā)明具有如下優(yōu)點本發(fā)明利用一種隱含活動輪廓模型作為一個新的能量函數(shù)，并通過使這些能量函數(shù)最小這種先驗信息來進行Bayesian重建。由于這種方法利用了圖像的邊緣梯度信息，重建后的圖像即可以保持高的分辨率，又能夠降低噪聲特性，而且重建后圖像的邊緣得到了有效地保持，不需要病人再進行PET，CT或MRI掃描，可以減輕病人的痛苦。


圖1為用來測試成像方法的胸腔模板圖像。
圖2為用來測試重建方法的投影數(shù)據(jù)。
圖3為用來測試重建方法的含泊松噪聲的投影數(shù)據(jù)。
圖4為用現(xiàn)有的貝葉斯(Bayesian)成像方法成像后的結(jié)果，此時投影數(shù)據(jù)不含有泊松噪聲。
圖5為用現(xiàn)有的貝葉斯(Bayesian)成像方法成像后的結(jié)果，此時投影數(shù)據(jù)含有泊松噪聲。
圖6為用本發(fā)明的能量函數(shù)替代現(xiàn)有的貝葉斯(Bayesian)成像方法中的能量函數(shù)成像后的結(jié)果，此時投影數(shù)據(jù)不含有泊松噪聲。
圖7為用本發(fā)明的能量函數(shù)替代現(xiàn)有的貝葉斯(Bayesian)成像方法中的能量函數(shù)成像后的結(jié)果，此時投影數(shù)據(jù)含有泊松噪聲。
具體實施例方式
一種基于隱含活動輪廓先驗的貝葉斯圖像重建方法1)獲取投影數(shù)據(jù)，根據(jù)待重建圖像的尺寸要求，確定初始圖像范圍，給定初始的灰度值大于1，并將2維圖像變換成1維向量進行計算，2)根據(jù)投影數(shù)據(jù)規(guī)模和要求成像圖像x的大小，計算系統(tǒng)概率矩陣P，3)將系統(tǒng)概率陣和初始圖像x相乘，得到前向投影a，4)將投影數(shù)據(jù)除以前向投影a，得到對投影數(shù)據(jù)的校正值c，5)將系統(tǒng)概率矩陣P乘以投影數(shù)據(jù)的校正值c，得到圖像成像迭代過程中的修正值xd，6)高斯函數(shù)對初始圖像作卷積計算，然后求這個計算結(jié)果的梯度模，再將這個計算結(jié)果除以0～1之間的任意一個數(shù)β1的平方，得到第一個數(shù)e，將第一個數(shù)e加1的結(jié)果的倒數(shù)乘以圖像的梯度模|x|，得到一個擴散速度g，將圖像的梯度x除以圖像的梯度模|x|，并對這個結(jié)果求散度，這個散度和氣球力F0之和乘以上述擴散速度g，得到能量函數(shù) 用差分方法將這個能量函數(shù)離散化，7)將系統(tǒng)概率矩陣P對它的每一行求和，把這個對行的和的β倍加上能量函數(shù)，得到一個權(quán)值，將初始圖像乘以修正值xd除以這個權(quán)值，得到了一個重建圖像，再將此重建圖像作為下一次迭代的初始圖像，返回到第3步，反復(fù)迭代直到重建后的圖像收斂，上述投影數(shù)據(jù)的獲取是從正電子發(fā)射計算機斷層顯像掃描儀上獲取，或者是從仿真模板圖像進行雷當(dāng)(Radon)變換，得到的投影數(shù)據(jù)。
實施例2本發(fā)明是通過對已有PET重建方法進行改進后得到的，具體實施方案的內(nèi)容如下1.已有的Bayesian重建方法一個正則化的最大似然估計的優(yōu)化問題為maxx≥0logP(Y=y|x)+φ(x)---(1)]]>
這里φ(x)即為懲罰項，它的意義是指圖像的先驗概率，因此公式(1)的具體形式可寫成maxx≥0logP(Y=y|x)+logP(x)---(2)]]>這里P(x)即表示先驗概率，根據(jù)Bayesian定理有P(x|Y=y)=P(Y=y|x)·P(x)P(Y=y)---(3)]]>我們注意到由于P(Y＝y(tǒng))與待估參數(shù)x無關(guān)，公式(2)等價于最大化后驗概率P(x|Y＝y(tǒng))，這是最大后驗估計的含義，亦稱Bayesian估計Bayesian重建方法首先要解決的問題就是如何選擇一個合適的先驗，先驗?zāi)Ｐ椭饕侵笀D像的光滑模型，圖像的光滑性主要表現(xiàn)為一種局部性質(zhì)，即鄰域內(nèi)的象素間的相互作用，這種作用越強，相鄰近的象素值越均勻，圖像越光滑，因此象素取某一值僅與其鄰域象素有關(guān)，我們可以通過給予鄰域象素值較均勻的構(gòu)型以大的概率來定量反映這種光滑特性，這就構(gòu)成Markov先驗概率模型。一般所采用的先驗都是基于某種鄰域結(jié)構(gòu)的Markov隨機場，但由于Markov隨機場沒有顯式地給出圖像的總體概率分布，因此大多用Gibbs隨機場來替代，這是因為Gibbs可以顯式地給出圖像的概率分布，而且在一定的條件下Markov隨機場等價于Gibbs隨機場。Gibbs隨機場給出的圖像概率分布為P(x)＝Z-1exp(-βV(x)) (4)這里V(x)=Σj∈SΣi∈NU(xi,xj)---(5)]]>稱為系統(tǒng)的能量函數(shù)，Z是歸一化因子，也稱分隔函數(shù)Z=Σxexp(-βV(x))---(6)]]>有了Gibbs先驗概率模型，就可以實現(xiàn)最大后驗估計，結(jié)合Poisson數(shù)據(jù)模型，最大后驗估計可以重新寫成maxx≥0Σi(yilogΣjpijxj-Σjpijxj-logyi!)-βV(x)---(7)]]>這里能量函數(shù)為 常用的結(jié)構(gòu)主要是八鄰域結(jié)構(gòu)，對于八鄰域結(jié)構(gòu)，一般我們還要對最鄰近象素與對角相鄰象素通過引入一定的權(quán)值加以區(qū)別，如果是八鄰域結(jié)構(gòu)，能量函數(shù)表示為 其中wij為權(quán)值，最鄰近點被賦值為1，對角鄰點則被賦值 函數(shù)(x)一般為偶函數(shù)。例如(x)＝log(cosh(x)) (10)利用期望最大方法對公式(7)求得的Bayesian重建方法為xj(k+1)=xj(k)(ΣiyipijΣjpijxj(k))/(Σipij+β∂V(x(k))∂xj),j=1,2,...n---(11)]]>其中yi表示第i個探測器所探到的光子數(shù)，0≤i≤m，m為探測器總數(shù)；xj表示第j個象素處發(fā)出的光子數(shù)。xj≥0，0≤j≤n，n為象素數(shù)；pij表示第j個象素處發(fā)出的光子能被第i個探測器檢測到的概率。
2.基于隱含活動輪廓先驗的Bayesian圖像的重建方法活動輪廓模型是一種有效的圖像分割方法，基于活動輪廓的圖像分割實質(zhì)上就是用活動輪廓逼近物體的邊緣，此過程可以通過能量最小來實現(xiàn)，其形變過程就是活動輪廓在外部能量和內(nèi)能(內(nèi)力)的作用下向物體邊緣靠近的過程，外力推動活動輪廓向著物體邊緣運動，而內(nèi)力則保持活動輪廓的光滑性和拓撲性。當(dāng)能量最小時，活動輪廓收斂到所要檢測的物體邊緣。由于這種方法同時考慮了幾何約束條件、圖像數(shù)據(jù)、輪廓形狀有關(guān)的能量最小等約束條件，所以用這種方法能得到滿意的分割效果。
由公式(3)我們可以得到p(x|y)∝p(x)p(y|x) (12)本發(fā)明定義Bayesian重建方法中的一個能量函數(shù)為V(x)＝V內(nèi)(x)+V外(x)(13)關(guān)于圖像x的先驗概率分布為 這里Z內(nèi)為歸一化因子，也就是內(nèi)力分隔函數(shù)，對于一個待重建的圖像x的觀測數(shù)據(jù)Y的似然為
這里Z外為歸一化因子，也就是外力分隔函數(shù)，因此 在這個新發(fā)明里，我們給出了一種隱含的活動輪廓模型，并利用它來表示能量函數(shù)V(x)。
本發(fā)明首先定義一個內(nèi)部的能量函數(shù)為V(x)=12∫Ω2λ(∂x∂t)2dx---(17)]]>其中λ為系數(shù)，外部的能量函數(shù)假設(shè)為0。此時∂x∂t=g(|▿Gσ*x|)|▿x|(div▿x|▿x|+F0)---(18)]]>用有限差分方法求解(18)式，并將公式(18)代入公式(17)，并將其離散化，得到新的能量函數(shù)，將這個能量函數(shù)代入公式(11)得到本發(fā)明中的Bayesian重建方法。
利用這種模型構(gòu)成的能量函數(shù)作為Bayesian先驗，可使得函數(shù)x的水平曲線沿著垂直x方向以g(|Gσ*x|)速度擴散，由此可見，圖像在邊緣處，也就是梯度最大處的鄰域?qū)嵭腥豕饣?，對邊緣點本身實行較小程度的光滑，而對其它區(qū)域?qū)嵭袕姽饣?。F0為常數(shù)，也稱氣球力，用來平衡內(nèi)外力的作用，它的取值范圍在0～10之間中的任意數(shù)，使重建后的圖像更清晰，其中Gσ為方差為σ的高斯函數(shù)。這里g(x,y)=11+(|▿Gσ(x,y)*I(x,y)|/β1)2---(19)]]>在活動輪廓模型中，g(i，j)為時間獨立的停止因子。這種活動輪廓作為Bayesian重建方法中的先驗函數(shù)的能量函數(shù)，重建后的圖像質(zhì)量較好，而且邊緣定位較準(zhǔn)確。
3.基于隱含活動輪廓先驗的Bayesian圖像重建方法的實驗結(jié)果本發(fā)明用一個計算機仿真的PET胸腔幻影模板來驗證本發(fā)明方法的可靠性。圖1顯示了這個胸腔幻影模板，模板大小為128×128象素矩陣，數(shù)據(jù)規(guī)模為185×180，即180個投影角度，每一個角度上有185條平行線，這里使平行線的間距與圖像象素的邊長相等，以便系統(tǒng)概率陣的確定。圖2和圖3分別表示投影數(shù)據(jù)不包含Poisson噪聲和包含Poisson噪聲時的情況。
圖4和圖5分別表示投影數(shù)據(jù)不包含噪聲和包含噪聲時用現(xiàn)有Bayesian方法對模板進行重建后的結(jié)果。圖6和圖7分別表示投影數(shù)據(jù)不包含噪聲和包含噪聲時用本發(fā)明方法對模板進行重建后的結(jié)果。
權(quán)利要求
1.一種基于隱含活動輪廓先驗的貝葉斯圖像重建方法，其特征在于采用以下步驟1)獲取投影數(shù)據(jù)，根據(jù)待重建圖像的尺寸要求，確定初始圖像范圍，給定初始的灰度值大于1，并將2維圖像變換成1維向量進行計算，2)根據(jù)投影數(shù)據(jù)規(guī)模和要求成像圖像x的大小，計算系統(tǒng)概率矩陣P，3)將系統(tǒng)概率陣和初始圖像x相乘，得到前向投影a，4)將投影數(shù)據(jù)除以前向投影a，得到對投影數(shù)據(jù)的校正值c，5)將系統(tǒng)概率矩陣P乘以投影數(shù)據(jù)的校正值c，得到圖像成像迭代過程中的修正值xd，6)將高斯函數(shù)對初始圖像作卷積計算，然后求這個計算結(jié)果的梯度模，再將這個計算結(jié)果除以0～1之間的任意一個數(shù)β1的平方，得到第一個數(shù)e，將第一個數(shù)e加1的結(jié)果的倒數(shù)乘以圖像的梯度模|x|，得到一個擴散速度g，將圖像的梯度x除以圖像的梯度模|x|，并對這個結(jié)果求散度，這個散度和氣球力F0之和乘以上述擴散速度g，得到能量函數(shù) 用差分方法將這個能量函數(shù)離散化，7)將系統(tǒng)概率矩陣P對它的每一行求和，把這個對行的和的β倍加上能量函數(shù)，得到一個權(quán)值，將初始圖像乘以修正值xd除以這個權(quán)值，得到了一個重建圖像，再將此重建圖像作為下一次迭代的初始圖像，返回到第3步，反復(fù)迭代直到重建后的圖像收斂。
2.根據(jù)權(quán)利要求1所述的基于隱含活動輪廓先驗的貝葉斯圖像重建方法，其特征在于投影數(shù)據(jù)的獲取是從正電子發(fā)射計算機斷層顯像掃描儀上獲取的。
3.根據(jù)權(quán)利要求1所述的基于隱含活動輪廓先驗的貝葉斯圖像重建方法，其特征在于投影數(shù)據(jù)的獲取是從仿真模板圖像進行雷當(dāng)變換，得到的投影數(shù)據(jù)。
全文摘要
本發(fā)明公開了一種基于隱含活動輪廓先驗的貝葉斯圖像重建方法，首先獲取投影數(shù)據(jù)，確定初始圖像范圍，計算系統(tǒng)概率矩陣，得到前向投影，再將投影數(shù)據(jù)除以前向投影，得到對投影數(shù)據(jù)的校正值，乘以系統(tǒng)概率矩陣，得到圖像成像迭代過程中的修正值，然后經(jīng)計算后得到能量函數(shù)，用差分方法將這個能量函數(shù)離散化，最后將系統(tǒng)概率矩陣對它的每一行求和乘以β加上能量函數(shù)，得到權(quán)值，將初始圖像乘以修正值除以權(quán)值，得到重建圖像，作為下一次迭代的初始圖像，反復(fù)迭代直到重建后的圖像收斂，本發(fā)明具有保持高分辨率，降低噪聲，且保持圖像邊緣，減輕病人的痛苦等優(yōu)點。
文檔編號A61B8/13GK1640362SQ20051003762
公開日2005年7月20日 申請日期2005年1月6日 優(yōu)先權(quán)日2005年1月6日
發(fā)明者朱宏擎, 舒華忠, 周鍵, 羅立民, 李松毅 申請人:東南大學(xué)