專利名稱：正電子發(fā)射計(jì)算機(jī)斷層顯像的全變分加權(quán)成像方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種全變分加權(quán)最小二乘法的圖像成建方法，尤其涉及一種正電子發(fā)射計(jì)算機(jī)斷層顯像的全變分加權(quán)成像方法。
背景技術(shù)：
正電子發(fā)射計(jì)算機(jī)斷層顯像(Positron emission tomography，PET)是當(dāng)今最高層次的核醫(yī)學(xué)影像技術(shù)，PET掃描儀是醫(yī)學(xué)界公認(rèn)的最先進(jìn)的大型醫(yī)療診斷成像設(shè)備。PET是一種利用正電子發(fā)射體放射性核素及其標(biāo)記化合物進(jìn)行人體局部或全身成像的現(xiàn)代醫(yī)學(xué)影像技術(shù)。PET成像過程是通過注射或吸入放射性藥物，經(jīng)過一段時(shí)間延時(shí)，當(dāng)放射性藥物傳送到待檢查器官后，開始進(jìn)行掃描，放射性同位素衰減時(shí)，它發(fā)射一個(gè)正電子，正電子經(jīng)過短距離移動(dòng)后，與一個(gè)電子相遇湮滅，并產(chǎn)生兩個(gè)幾乎相反方向傳播的高能光子，如果兩個(gè)光子在一個(gè)很短的時(shí)間被檢測(cè)到，則記錄一個(gè)事件，由探測(cè)器探測(cè)到的這些事件組成一個(gè)斷層上放射性同位素濃度分布的投影數(shù)據(jù)，使用這些投影數(shù)據(jù)，根據(jù)成像方法可以得到該斷層上放射性同位素分布的二維圖像。引起PET成像誤差的原因有很多，如正電子類藥物強(qiáng)度的快速衰減、高計(jì)數(shù)率造成系統(tǒng)死時(shí)間損失、隨機(jī)符合、散射和人體吸收衰減的影響、統(tǒng)計(jì)噪聲等也嚴(yán)重地影響了PET成像質(zhì)量，因此PET圖像高精度成像對(duì)PET的商用化和普及起著決定性地作用。由于PET投影數(shù)據(jù)的不完備性以及重建方法的不適定性，造成重建后的圖像邊緣不規(guī)則，噪化現(xiàn)象明顯。

發(fā)明內(nèi)容
本發(fā)明提供一種能夠提高成像后圖像的質(zhì)量，減弱甚至消除噪聲對(duì)成像影響的正電子發(fā)射計(jì)算機(jī)斷層顯像的全變分加權(quán)成像方法。
本發(fā)明采用如下技術(shù)方案一種最小二乘法的正電子發(fā)射計(jì)算機(jī)斷層顯像的全變分加權(quán)成像方法1)獲取投影數(shù)據(jù)，根據(jù)待重建圖像的尺寸要求，確定初始圖像范圍，給定初始的灰
2)度值大于1，并將2維圖像變換成1維向量進(jìn)行計(jì)算，3)根據(jù)投影數(shù)據(jù)規(guī)模和要求成像圖像x的大小，計(jì)算系統(tǒng)概率矩陣P，3)將系統(tǒng)概率矩陣和初始圖像x相乘，得到前向投影a，4)取自正電子發(fā)射計(jì)算機(jī)斷層顯像掃描儀的每一個(gè)投影數(shù)據(jù)y進(jìn)行平方計(jì)算，得到y(tǒng)2，再將這個(gè)數(shù)除以前向投影a的平方，得到投影數(shù)據(jù)的修正值c，c＝y(tǒng)2/a25)將系統(tǒng)概率矩陣P乘以投影數(shù)據(jù)的修正值c，得到圖像成像迭代過程中圖像的修正值d，6)將離散化的初始圖像進(jìn)行全變分，再將這個(gè)離散化的全變分初始圖像對(duì)每一個(gè)象素點(diǎn)求導(dǎo)，得到目標(biāo)函數(shù)的校正值 7)將圖像修正值d除以1與β倍的目標(biāo)函數(shù)的校正值 得到用于重建圖像的目標(biāo)函數(shù)，再將這個(gè)用于重建圖像的目標(biāo)函數(shù)值乘以初始圖像，得到迭代更新后的圖像，再將這一個(gè)圖像作為初始圖像，返回到第3步，重復(fù)這個(gè)過程直到重建后的圖像收斂。
與現(xiàn)有技術(shù)相比，本發(fā)明具有如下優(yōu)點(diǎn)本發(fā)明用全變分作為正則項(xiàng)并將其融合到加權(quán)最小二乘PET成像方法中，提高成像后圖像的質(zhì)量，消除噪聲對(duì)成像的影響。


圖1為用來測(cè)試成像方法的腹腔模板圖像。
圖2為腹腔模板的二值圖像。
圖3為用現(xiàn)有加權(quán)最小二乘法的PET成像方法成像后的結(jié)果。
圖4為對(duì)測(cè)試模板的投影數(shù)據(jù)加入噪聲后用現(xiàn)有加權(quán)最小二乘法的PET成像方法成像后的結(jié)果。
圖5為圖4二值化后的結(jié)果。
圖6為用本發(fā)明成像后的結(jié)果，其中β為0.002。
圖7為用本發(fā)明成像后的結(jié)果，其中β為0.005。
圖8為用本發(fā)明成像后的結(jié)果，其中β為0.008。
圖9為用本發(fā)明成像后的結(jié)果，其中β為和0.01。
圖10為對(duì)測(cè)試模板加入噪聲后，用本發(fā)明成像后的結(jié)果，其中β為0.002。
圖11為對(duì)測(cè)試模板加入噪聲后，用本發(fā)明成像后的結(jié)果，其中β為0.005。
圖12為對(duì)測(cè)試模板加入噪聲后，用本發(fā)明成像后的結(jié)果，其中β為0.008。
圖13為對(duì)測(cè)試模板加入噪聲后，用本發(fā)明成像后的結(jié)果，其中β為0.01。
圖14為圖10二值化后的結(jié)果。
圖15為圖11二值化后的結(jié)果。
圖16為圖12二值化后的結(jié)果。
圖17為圖13二值化后的結(jié)果。
圖18為投影數(shù)據(jù)不含有噪聲的均方根誤差分析結(jié)果。
圖19為投影數(shù)據(jù)含有噪聲的均方根誤差分析結(jié)果。
圖20為投影數(shù)據(jù)不含有噪聲的偏差分析結(jié)果。
圖21為投影數(shù)據(jù)含有噪聲的偏差分析結(jié)果。
圖22為圖3第48列線圖。
圖23為圖9第48列線圖。
具體實(shí)施例方式
實(shí)施例1一種最小二乘法的正電子發(fā)射計(jì)算機(jī)斷層顯像的全變分加權(quán)成像方法1)獲取投影數(shù)據(jù)，根據(jù)待重建圖像的尺寸要求，確定初始圖像范圍，給定初始的灰度值大于1，并將2維圖像變換成1維向量進(jìn)行計(jì)算，2)根據(jù)投影數(shù)據(jù)規(guī)模和要求成像圖像x的大小，計(jì)算系統(tǒng)概率矩陣P，3)將系統(tǒng)概率矩陣和初始圖像x相乘，得到前向投影a，4)取自正電子發(fā)射計(jì)算機(jī)斷層顯像掃描儀的每一個(gè)投影數(shù)據(jù)y進(jìn)行平方計(jì)算，得到y(tǒng)2，再將這個(gè)數(shù)除以前向投影a的平方，得到投影數(shù)據(jù)的修正值c，c＝y(tǒng)2/a25)將系統(tǒng)概率矩陣P乘以投影數(shù)據(jù)的修正值c，得到圖像成像迭代過程中圖像的修正值d，6)將離散化的初始圖像進(jìn)行全變分，再將這個(gè)離散化的全變分初始圖像對(duì)每一個(gè)象素點(diǎn)求導(dǎo)，得到目標(biāo)函數(shù)的校正值 7)將圖像修正值d除以1與β倍的目標(biāo)函數(shù)的校正值 得到用于重建圖像的目標(biāo)函數(shù)，再將這個(gè)用于重建圖像的目標(biāo)函數(shù)值乘以初始圖像，得到迭代更新后的圖像，再將這一個(gè)圖像作為初始圖像，返回到第3步，重復(fù)這個(gè)過程直到重建后的圖像收斂，上述投影數(shù)據(jù)的獲取是從正電子發(fā)射計(jì)算機(jī)斷層顯像掃描儀上獲取，或者是從仿真模板圖像進(jìn)行雷當(dāng)(Radon)變換，得到的投影數(shù)據(jù)。
實(shí)施例2本發(fā)明是通過對(duì)已有PET成像方法進(jìn)行改進(jìn)后得到的，具體實(shí)施方案的內(nèi)容如下1.現(xiàn)有的加權(quán)最小二乘估計(jì)方法現(xiàn)在商用PET機(jī)上所用的成像模型主要是假定PET掃描儀所探測(cè)到的光子發(fā)射過程是服從泊松分布基礎(chǔ)上的，即yi~Poisson{Σjpijxj}---(1)]]>其中yi表示第i個(gè)探測(cè)器所探測(cè)到的光子數(shù)，0≤i≤m，m為探測(cè)器總數(shù)；xj表示第j個(gè)象素處發(fā)出的光子數(shù)，xj≥0，0≤j≤n，n為象素?cái)?shù)；pij表示第j個(gè)象素處發(fā)出的光子能被第i個(gè)探測(cè)器檢測(cè)到的概率。根據(jù)這個(gè)假設(shè)，我們建立一個(gè)加權(quán)最小二乘估計(jì)的PET重建模型。該模型根據(jù)數(shù)據(jù)的方差來決定具體的權(quán)值。這是因?yàn)榉讲疃糠从沉藰颖敬砜傮w期望的可信度，方差越大數(shù)據(jù)可信度越低，故合理的做法顯然是給予方差較小的數(shù)據(jù)以較大的權(quán)值，現(xiàn)在余下的就是要確定權(quán)值與數(shù)據(jù)方差的定量關(guān)系以使估計(jì)值的方差最小或精度最高，由統(tǒng)計(jì)學(xué)知識(shí)，要做到這點(diǎn)應(yīng)使權(quán)值反比于方差。對(duì)于Poisson統(tǒng)計(jì)誤差，我們知道數(shù)據(jù)的方差等于期望，故現(xiàn)在可以描述加權(quán)最小二乘估計(jì)準(zhǔn)則下的建模問題，也就是說我們可以以下面優(yōu)化問題的解作為最終估計(jì)值。它可以表示為Φ:argmin{(Px-y)TW-1(Px-y)}s.t.:x≥0---(2)]]>Φ(x)=Σi=1m((Px)i-yi)2(Px)i---(3)]]>這里W是一個(gè)m×m的權(quán)對(duì)角矩陣，其第i個(gè)元素為(Px)iwij＝diag((Px)1，(Px)2，....，(Px)m) (4)令Φ(x)的一階偏導(dǎo)為零，根據(jù)Kuhn-Tucker條件，我們有∂∂xj(Φ(x))=Σi=1m(-yi2pij(Px)i2+pij)=0,xj>0---(5)]]>∂∂xj(Φ(x))=Σi=1m(-yi2pij(Px)i2+pij)≥0,xj=0---(6)]]>
因此我們得出一個(gè)固定點(diǎn)的加權(quán)最小二乘法的PET成像算法xj(k+1)=xj(k)Σi=1myi2pij(Σj=1npijxj(k))2,j=1,2,Λ,n---(7)]]>這種PET成像方法得到的圖像對(duì)比度比較高，重建圖像質(zhì)量也比較好，但噪化現(xiàn)象比較嚴(yán)重，邊緣上的偽影難以消除。為了驗(yàn)證這種方法的重建效果，我們用一個(gè)計(jì)算機(jī)仿真的PET腹腔幻影模板來驗(yàn)證。圖1顯示了這個(gè)腹腔模板，模板圖像大小為96×96象素矩陣，數(shù)據(jù)規(guī)模為139×180，即180個(gè)投影角度，每一個(gè)角度上有139條平行投影線，我們使平行線的間距與圖像像素的邊長相等，以便系統(tǒng)概率矩陣P的確定。圖2是這個(gè)幻影模板的二值圖像。圖3是用公式(7)成像后的結(jié)果。從圖3中，我們可以看到邊緣上有明顯的偽影，且圖像較模糊，由此可以看出，投影不加噪聲用現(xiàn)有的方法很難重建出較理想的圖像。圖4是對(duì)這個(gè)模板加噪聲后用公式(7)成像后的結(jié)果。
圖5是圖4二值化后的結(jié)果，它進(jìn)一步說明用現(xiàn)有方法成像后，圖像質(zhì)量很不理想。
2.本發(fā)明的全變分正則化加權(quán)最小二乘PET成像方法為了提高圖像質(zhì)量，降低噪聲并保持邊緣，我們發(fā)明了一種用全變分作為正則項(xiàng)的加權(quán)最小二乘PET成像方法。全變分的使用主要在于它能有效的去噪聲的同時(shí)能保持邊緣盡可能不被破壞掉。全變分的表達(dá)式為TV(f)=∫Ω|▿f|dxdy=∫Ωfx2+fy2dxdy---(8)]]>這里fx=∂∂xf,]]>fy=∂∂yf.]]>上式關(guān)于i，j的差分表達(dá)式為UTV=Σi,j(fi+1,j-fi,j)2+(fi,j+1-fi,j)2+ϵ2---(9)]]>我們發(fā)現(xiàn)參數(shù)ε應(yīng)小于等于1％的f最大值。ε值過大會(huì)平滑掉邊緣。公式(9)的偏導(dǎo)數(shù)為∂UTV∂fi,j=fi,j-fi-1,j(fi,j-fi-1,j)2+(fi-1,j+1-fi-1,j)2+ϵ2]]>+fi,j-fi,j-1(fi+1,j-1-fi,j-1)2+(fi,j-fi,j-1)2+ϵ2]]>-fi+1,j+fi,j+1-2fi,j(fi+1,j-fi,j)2+(fi,j+1-fi,j)2+ϵ2---(10)]]>本發(fā)明的基于全變分正則項(xiàng)的加權(quán)最小目標(biāo)函數(shù)Jβ替換掉公式(3)所示的現(xiàn)有加權(quán)最小二乘目標(biāo)函數(shù)φ(x)，重建后的圖像 由使新的目標(biāo)函數(shù)Jβ(x)最小給出
x)=argminx(Jβ(x))---(11)]]>本發(fā)明這里的目標(biāo)函數(shù)由兩部分組成普通的加權(quán)最小二乘項(xiàng)和全變分正則項(xiàng)，本發(fā)明的目標(biāo)函數(shù)Jβ(x)為Jβ(x)＝φ(y，Px)+βU(12)這里β為權(quán)因子，它將影響全變分正則項(xiàng)在成像過程中的作用程度，對(duì)公式(12)求一階偏導(dǎo)。并將公式(10)帶入到公式(12)，那么針對(duì)每一個(gè)象素xj求Jβ的一階偏導(dǎo)為∂Jβ(x)∂xj=Σi(-yi2pij(Px)2i+pij)+β∂UTV∂xj---(13)]]>因?yàn)?amp;Sigma;i=1mpij=1]]>根據(jù)Kuhn-Tucher條件，解決這個(gè)問題的固定點(diǎn)成像迭代式為xj(k+1)=xj(k)(1+β∂UTV∂xj)Σi=1mpijyi2(Σj=1npijxj(k))2---(14)]]>由于本發(fā)明將全變分正則項(xiàng)加入到普通的加權(quán)最小二乘PET成像方法中，使成像后的PET圖像精度得到較大提高，有效地除去了噪聲，全變分的作用在投影有噪聲的成像中尤為明顯。在公式(14)中參數(shù)β的作用是用來調(diào)節(jié)全變分正則項(xiàng)對(duì)算法的影響程度，隨著β的增加，正則項(xiàng)的功能加強(qiáng)，圖像進(jìn)一步被平滑，當(dāng)β為零時(shí)，公式(14)變成普通的加權(quán)最小二乘成像方法。
圖6和圖7是用本發(fā)明方法(公式(14))成像后的結(jié)果，其中β分別為0.002和0.005。圖8和圖9也是用我們發(fā)明的新算法(公式(14))成像后的結(jié)果，其中β分別為0.008和0.01。由此試驗(yàn)，發(fā)現(xiàn)在投影數(shù)據(jù)沒有噪聲的情況下，參數(shù)β對(duì)重建圖像的影響不大，但重建后的圖像質(zhì)量遠(yuǎn)比普通的加權(quán)最小二乘重建法(圖3)要好。
圖10和圖11是用我們發(fā)明的新方法(公式(14))對(duì)投影數(shù)據(jù)加噪聲后的成像結(jié)果，其中β分別為0.002和0.005。圖12和圖13也是用本發(fā)明(公式(14))對(duì)投影數(shù)據(jù)加噪聲后的成像結(jié)果，其中β分別為0.008和0.01。由這個(gè)試驗(yàn)可以發(fā)現(xiàn)隨著β的增加，成像后圖像的平滑作用增強(qiáng)，邊緣也得到了有效的保持。成像后的圖像遠(yuǎn)比在投影具有同樣噪聲條件的普通加權(quán)最小二乘成像方法(圖4)要好。
圖14是圖10二值化后的結(jié)果，圖15是圖11二值化后的結(jié)果，圖16是圖12二值化后的結(jié)果，圖17是圖13二值化后的結(jié)果，這些二值化后的圖像進(jìn)一步說明了本發(fā)明所述方法能重建出高精度的圖像，這是由于全變分正則項(xiàng)的加入極大地改進(jìn)了原有方法的成像精度。
本發(fā)明對(duì)PET圖像成像特別有效，這是由于PET成像的原理導(dǎo)致了PET成像圖像的分辨率較低，圖像信噪比較低，噪化現(xiàn)象嚴(yán)重，而本發(fā)明的這種方法對(duì)信噪比低的PET投影重建尤為有效，由圖13和圖17可知隨著β的增加，全變分正則化作用增強(qiáng)，去噪效果明顯，圖像二值化后的結(jié)果(圖17)接近原始模板二值化后的結(jié)果(圖2)，且遠(yuǎn)比圖4二值化的結(jié)果要好。圖4的二值化結(jié)果見圖5。
3.本發(fā)明的測(cè)試結(jié)果本發(fā)明的全變分加權(quán)最小二乘法的PET成像方法的測(cè)試是在一臺(tái)Pentium 4 CPU，2.4GHz，1.00GB上進(jìn)行的。為了測(cè)試本發(fā)明方法的有效性，并對(duì)全變分加權(quán)最小二乘法的PET成像方法有一個(gè)比較清楚的認(rèn)識(shí)，用某種準(zhǔn)則來測(cè)試以某種成像方法生成的圖像與原始標(biāo)準(zhǔn)圖像的接近程度來衡量這種成像方法的優(yōu)劣，我們采用以下兩種準(zhǔn)則均方根誤差，即用測(cè)試模板圖像和成像后圖像之間的均方根誤差(root mean square error(RMS))來評(píng)估成像后圖像的質(zhì)量。均方根誤差被定義為RMS=[1nΣj=1n(xjrec-xjorg)2]12---(15)]]>這里xjorg和xjrec分別代表一個(gè)仿真的模板圖像和重建后的圖像在象素j這個(gè)位置的值。好的成像方法將具有較小的均方根誤差。圖18顯示了投影數(shù)據(jù)不含有噪聲的均方根誤差分析結(jié)果，圖19顯示了投影數(shù)據(jù)含有噪聲的均方根誤差分析結(jié)果，這些結(jié)果表明了新發(fā)明的成像方法重建后的圖像遠(yuǎn)比現(xiàn)有的加權(quán)最小二乘方法更接近測(cè)試模板數(shù)據(jù)，因?yàn)楸景l(fā)明成像后的圖像和測(cè)試模板圖像之間的均方根誤差要小。
偏差，即Variance=1n-1Σj=1n(xjrec-x‾org)2---(16)]]>這里xorg代表用來測(cè)試模板圖像的平均灰度值。成像方法越好，偏差越小。圖20顯示了投影數(shù)據(jù)不含有噪聲的偏差分析結(jié)果，圖21顯示了投影數(shù)據(jù)含有噪聲的偏差分析結(jié)果。
圖22顯示了投影數(shù)據(jù)不含有噪聲情況下用現(xiàn)有方法成像后的圖像的第48列輪廓線與原始模板圖像第48列輪廓線接近的程度。
圖23顯示了投影數(shù)據(jù)不含有噪聲情況下用新發(fā)明方法對(duì)投影數(shù)據(jù)成像后圖像的第48列輪廓線與原始模板圖像第48列輪廓線接近的程度。由此可見，本發(fā)明重建的結(jié)果更接近原始測(cè)試模板圖像。
權(quán)利要求
1.一種最小二乘法的正電子發(fā)射計(jì)算機(jī)斷層顯像的全變分加權(quán)成像方法，其特征在于采用下列步驟1)獲取投影數(shù)據(jù)，根據(jù)待重建圖像的尺寸要求，確定初始圖像范圍，給定初始的灰度值大于1，并將2維圖像變換成1維向量進(jìn)行計(jì)算，2)根據(jù)投影數(shù)據(jù)規(guī)模和要求成像圖像x的大小，計(jì)算系統(tǒng)概率矩陣P，3)將系統(tǒng)概率矩陣和初始圖像x相乘，得到前向投影a，4)取自正電子發(fā)射計(jì)算機(jī)斷層顯像掃描儀的每一個(gè)投影數(shù)據(jù)y進(jìn)行平方計(jì)算，得到y(tǒng)2，再將這個(gè)數(shù)除以前向投影a的平方，得到投影數(shù)據(jù)的修正值c，c＝y(tǒng)2/a25)將系統(tǒng)概率矩陣P乘以投影數(shù)據(jù)的修正值c，得到圖像成像迭代過程中圖像的修正值d，6)將離散化的初始圖像進(jìn)行全變分，再將這個(gè)離散化的全變分初始圖像對(duì)每一個(gè)象素點(diǎn)求導(dǎo)，得到目標(biāo)函數(shù)的校正值 7)將圖像修正值d除以1與β倍的目標(biāo)函數(shù)的校正值 得到用于重建圖像的目標(biāo)函數(shù)，再將這個(gè)用于重建圖像的目標(biāo)函數(shù)值乘以初始圖像，得到迭代更新后的圖像，再將這一個(gè)圖像作為初始圖像，返回到第3步，重復(fù)這個(gè)過程直到重建后的圖像收斂。
2.根據(jù)權(quán)利要求1所述的基于隱含活動(dòng)輪廓先驗(yàn)的貝葉斯圖像重建方法，其特征在于投影數(shù)據(jù)的獲取是從正電子發(fā)射計(jì)算機(jī)斷層顯像掃描儀上獲取的。
3.根據(jù)權(quán)利要求1所述的基于隱含活動(dòng)輪廓先驗(yàn)的貝葉斯圖像重建方法，其特征在于投影數(shù)據(jù)的獲取是從仿真模板圖像進(jìn)行雷當(dāng)變換，得到的投影數(shù)據(jù)。
全文摘要
本發(fā)明公開了一種最小二乘法的正電子發(fā)射計(jì)算機(jī)斷層顯像的全變分加權(quán)成像方法，首先獲取投影數(shù)據(jù)，確定初始圖像范圍，計(jì)算系統(tǒng)概率矩陣，再將系統(tǒng)概率矩陣和初始圖像相乘，得到前向投影，將投影數(shù)據(jù)的修正值與系統(tǒng)概率矩陣相乘，得到圖像成像迭代過程中圖像的修正值，然后將離散化的初始圖像進(jìn)行全變分，對(duì)每一個(gè)象素點(diǎn)求導(dǎo)，得到目標(biāo)函數(shù)的校正值，最后經(jīng)運(yùn)算，得到重建圖像的目標(biāo)函數(shù)，將此目標(biāo)函數(shù)值乘以初始圖像，得到迭代更新后的圖像，再將這一個(gè)圖像作為初始圖像，重復(fù)這個(gè)過程直到重建后的圖像收斂，本發(fā)明具有提高成像后圖像的質(zhì)量，消除噪聲對(duì)成像的影響等優(yōu)點(diǎn)。
文檔編號(hào)G06T5/00GK1641700SQ20051003762
公開日2005年7月20日 申請(qǐng)日期2005年1月6日 優(yōu)先權(quán)日2005年1月6日
發(fā)明者朱宏擎, 舒華忠, 羅立民, 周鍵, 李松毅 申請(qǐng)人:東南大學(xué)