陳正慧，韋民紅

(淮北師范大學(xué)物理與電子信息學(xué)院，安徽 淮北 235000)

0 引 言

CT(industrial computed tomography)技術(shù)的核心是投影重建圖像的理論，現(xiàn)在圖像重建的方法主要有解析和迭代算法[1]。迭代算法近幾年來由于計算機(jī)的高速發(fā)展，重建速度作為其最大的缺點得到了很大的改善，這個缺點已經(jīng)不再制約著迭代算法的發(fā)展，因此迭代算法再度引起了人們的關(guān)注。

為了進(jìn)一步的發(fā)展，中國工程物理研究首先建成了基于反應(yīng)堆的冷中子和熱中子CT裝置，并取得了研究成果，這個裝置也是以解析算法為重建算法。但是在人們的實際生活中，經(jīng)常有一些奇形怪狀的樣品，這些樣品不能完整的獲取投影數(shù)據(jù)，這就是解析算法的短板了，解析算法對于投影數(shù)據(jù)的完整性要求性高，如果仍用解析算法的話難以得到理想的重建效果了。因此迭代算法在CT中的研究和發(fā)展就勢在必行了。迭代算法最具有代表性的是ART算法，改進(jìn)ART算法也就成為了研究熱點[2-3]。

1 模擬退火算法

模擬退火算法最開始是Metropplis等一批學(xué)者提出來的[4]。這是由物理中固體的降溫過程啟發(fā)而來的。這個算法具有較好的魯棒性和全局收斂性。該算法主要有以下三個過程。

加溫過程。粒子在高溫的時候具有較高的能量，給粒子比較高的溫度，這樣粒子就會重新排序并且會偏離原來的平衡位置。

等溫過程。這個過程中系統(tǒng)封閉，在系統(tǒng)中，粒子會慢慢往自由能小的方向自發(fā)轉(zhuǎn)化，但這個過程中溫度不變而且會與周圍環(huán)境交換熱量,最終粒子會降到自由能最小，而這就是粒子達(dá)到平衡的過程。

圖1

冷卻過程。這個過程就是降溫的過程，當(dāng)粒子完全冷卻時就會形成低能狀態(tài)的晶體。

值得注意的是，這個等溫過程就是算法Metroplis的抽樣過程，而最重要的也是這個過程，這個也是模擬退火能得到全局最優(yōu)解的關(guān)鍵，模擬退火以較低的概率接受惡化解，較高的概率接受優(yōu)解。這樣就可以跳出局部最優(yōu)解而得到全局最優(yōu)解。

假設(shè)是解決一個尋找最小值得優(yōu)化問題，模擬退火尋優(yōu)算法就是將模擬退火的思想應(yīng)用在尋優(yōu)問題上。

如果想要優(yōu)化這樣一個函數(shù)F:x→R，x∈A，x是優(yōu)化問題其中的一個可行解，R={y|y∈A,y>0},定義域用A表示，N(x)?A就是x的一個鄰域集合。

值得注意的是，降溫的過程得足夠慢，不然最后很可能錯過最優(yōu)解，得不到全局最優(yōu)解。

2 CRRT算法原理

首先把目標(biāo)圖像f(x,y)看成是一個包含N個網(wǎng)格的正方形網(wǎng)格，其中N=n×n。

如圖1所示。

這樣劃分的其中一個小網(wǎng)格就可以看成一個像素，小網(wǎng)格的內(nèi)部f(x,y)可以當(dāng)成一個常數(shù)值，射線就以一定的方式射入待測目標(biāo)圖像的區(qū)域網(wǎng)格，對應(yīng)的投影值pi就可以表示為第i條射線。目標(biāo)圖像重建的問題就轉(zhuǎn)化為了一個線性方程組

(1)

這里的M是射線的總線條數(shù)，N表示總的單元格數(shù)，N=n×n,有n行n列。而式中的f'j是對第j個單元格的真實像素fj的估計值。Wij是加權(quán)因子，這時候要解決的關(guān)鍵問題是求取加權(quán)因子Wij，在這里先設(shè)單元格的寬度為ε，而設(shè)射線的寬度是0。把單元格第i條射線和第j個像素相交的長度的比值定義為Wij，Pi表示第i條射線的線積分，這個也就是投影數(shù)據(jù)。

峰值信噪比

(2)

t為圖像的位數(shù)，m和n分別代表圖像的高度和寬度；I，K分別為參考圖像和迭代圖像。P值越大，表示重建效果越好。

歸一化平均距離判據(jù)

(3)

Ii,j表示參考圖像第i行第j列的像素濃度，Ki,j表示重建后第i行，第j列的像素濃度。r的值越小，表示誤差越小，重建的效果越好。

CRRT算法步驟如下：

(1)取圖像的初始估計值fj0=0(一般取0);

(2)根據(jù)模擬退火算法的原理，隨機(jī)選取一個松弛因子λk,這個松弛因子應(yīng)該服從正態(tài)分布(λ0,1/k)。k是迭代次數(shù)，λ0選取的是0.25。

(3)迭代算法如下:

(4)

(4)迭代接受判別依據(jù)

圖2 原始頭模型 圖3 傳統(tǒng)算法 圖4 CRRT算法

如果Pk+1>Pk,說明這一迭代結(jié)果比上一迭代結(jié)果好，選的松弛因子比較合適，可以繼續(xù)下一步迭代。反之，如果Pk+1>Pk，說明迭代朝著惡劣的結(jié)果發(fā)展，就放棄該次迭代結(jié)果。

(5)算法結(jié)束依據(jù)

如果||fjk+1-fjk||<ξ(ξ很小，一般取10-6)，結(jié)束算法，這個fjk+1就是這次的重建結(jié)果。

3 CRRT算法的檢驗與結(jié)果

使用MATLAB實現(xiàn)了該次實驗，實驗選用了Shepp Logan[5]頭模型，Shepp Logan頭模型是一個大的橢圓里面有九個小的橢圓，而且這九個橢圓的位置和大小都不一樣，這個Sheep Logan頭模型就像是一個人體大腦斷層圖像。實驗重建的Shepp Logan頭模型大小是180×180，探測元的個數(shù)是260個，取180個投影角采樣。用文中的算法和傳統(tǒng)的算法進(jìn)行了對比[6-9],實驗結(jié)果如圖2-4所示。

由圖2、圖3、圖4所示，傳統(tǒng)算法相比較于CRRT算法有較明顯的偽影和模糊，這是人們?nèi)庋劭梢杂^察得出來的，但是肉眼觀察會有誤差，也無法用數(shù)據(jù)證明，所以用歸一化平均距離判據(jù)r作為評判標(biāo)準(zhǔn)，并將傳統(tǒng)算法和CRRT算法的歸一化平均距離判據(jù)r作為對比，得出的數(shù)據(jù)如圖5所示。

a 傳統(tǒng)ART算法 b CRRT算法

圖5中a圖是r迭代了十次的變幻趨勢，可知r一直趨于穩(wěn)步降低的趨勢。圖5中b圖是CRRT算法迭代了十次后再優(yōu)化了十次松弛因子，由實驗數(shù)據(jù)可知，在選到更好的松弛因子后，r值明顯降低并趨于平穩(wěn)狀態(tài)。

實驗結(jié)果表明，CRRT算法可以優(yōu)化松弛因子。

4 結(jié) 論

將ART算法和模擬退火相結(jié)合，采用模擬退火來優(yōu)化ART算法的松弛因子。結(jié)果表明，新的算法明顯比傳統(tǒng)算法的偽影要少的多，圖像的質(zhì)量也比較好。值得指出的是，實驗采用的是模擬退火算法優(yōu)化的松弛因子，今后的研究可能還會有更好的優(yōu)化算法來替代模擬退火算法。