梁 婕

(山西師范大學(xué)臨汾學(xué)院，山西 臨汾 041000)

0 引言

計(jì)算機(jī)斷層成像(Computed Tomography，CT)技術(shù)，是一門(mén)利用X射線穿透物體的衰減信息進(jìn)行重建來(lái)獲取物體的短促圖像信息的技術(shù)，它被公認(rèn)為20世紀(jì)后期最偉大的科技成果之一.CT技術(shù)的核心是CT圖像重建，其特點(diǎn)是可以在不對(duì)物體造成任何的物理性損傷的基礎(chǔ)上，獲取被檢測(cè)物體內(nèi)部結(jié)構(gòu)的圖像.目前已廣泛應(yīng)用于工業(yè)無(wú)損檢測(cè)、醫(yī)學(xué)診斷、農(nóng)林業(yè)、地球物理等領(lǐng)域[1].

CT圖像重建算法大致可以分為以Radon變換為理論基礎(chǔ)的解析類重建算法和以解方程為主要思想的迭代類重建算法兩大類.解析法是以濾波反投影算法(Filtered Back Projection，F(xiàn)BP)為主，優(yōu)點(diǎn)是運(yùn)算速度快，但其抗噪聲性能差，對(duì)投影數(shù)據(jù)的完備性要求較高.聯(lián)立迭代重建算法(Simultaneous Iterative Reconstruction Technique，SIRT)是迭代重建中典型的算法之一，其優(yōu)勢(shì)在于對(duì)數(shù)據(jù)的完備性要求不高，缺點(diǎn)是計(jì)算量大，重建速度慢，難以用硬件實(shí)現(xiàn)[2].近年來(lái)，計(jì)算機(jī)技術(shù)的發(fā)展十分快速，這使得SIRT算法的缺點(diǎn)已成為次要矛盾，其優(yōu)勢(shì)更加突出，人們更為關(guān)注的是如何提高該算法的重建質(zhì)量.

在實(shí)際應(yīng)用中，由于投影數(shù)據(jù)的不完備或者存在著噪聲，使得圖像重建的效果并不明顯，繼而給一些后繼的圖像處理，如特征提取、圖像分割、邊緣檢測(cè)等帶來(lái)許多困難.因此，需要通過(guò)對(duì)圖像進(jìn)行去噪處理來(lái)提高圖像重建的質(zhì)量.近年來(lái)，偏微分方程(PDE)方法已廣泛應(yīng)用于圖像去噪，在噪聲消除和與邊緣保留方面取得了較好的效果，其中較為常用的去噪模型是PM模型[3].1990年，Perona 和Malik 提出了經(jīng)典的非線性擴(kuò)散模型.基于SIRT重建算法，本文利用PM模型的去噪原理，分析了該模型的優(yōu)缺點(diǎn)并對(duì)其擴(kuò)散系數(shù)進(jìn)行了改進(jìn)，給出了一種新的迭代重建算法，記作PM-SIRT算法.仿真實(shí)驗(yàn)結(jié)果驗(yàn)證了該算法的有效性.

1 背景知識(shí)

1.1 迭代重建算法原理

與解析法不同，迭代重建算法是在離散域中進(jìn)行的.把圖像重建區(qū)域劃分為一個(gè)n×n的正方形網(wǎng)格，且每個(gè)網(wǎng)格(像素)內(nèi)f為常數(shù)，這樣就可以建立一個(gè)線性方程組來(lái)描述像素與投影數(shù)據(jù)之間的關(guān)系，從而把圖像重建問(wèn)題轉(zhuǎn)化成為解該方程組[4].

(1)

式中:每一個(gè)方程代表一個(gè)超平面，pi表示第i條射線的投影值，fj表示第j個(gè)網(wǎng)格內(nèi)的常量，wij代表權(quán)因子，反映了第j個(gè)像素對(duì)第i條射線投影值的貢獻(xiàn)，N(N=n×n)代表像素總數(shù).在實(shí)際應(yīng)用中，為提高重建速度，減少計(jì)算量，對(duì)wij簡(jiǎn)化，假設(shè)像素寬度等于射線寬度，像素值集中在像素的中心，當(dāng)射線穿過(guò)像素中心時(shí)，權(quán)因子定義為1，即

(2)

式(1)可以用矩陣形式表示為：

WF=P

(3)

式中:W為M×N維投影系數(shù)矩陣；F為N維圖像向量；P為M維投影數(shù)據(jù)向量.如何利用W和P求出F是圖像重建的任務(wù)，這樣就將圖像重建問(wèn)題轉(zhuǎn)化為求解一個(gè)大型線性方程組的問(wèn)題.

1.2 SIRT算法

在實(shí)際應(yīng)用中，因?yàn)橥队熬仃嘩非常龐大，求取圖像F時(shí)，一般不直接用M×N維矩陣的逆矩陣運(yùn)算.隨后，由WF=P推導(dǎo)出了ART算法，且Gilbert于1972 年提出了改進(jìn)技術(shù)聯(lián)立迭代重建算法(SIRT).在SIRT算法中，一個(gè)像素的平均修正值是通過(guò)利用這一像素內(nèi)通過(guò)的所有射線修正值來(lái)確定，這樣可以平滑噪聲，相比于ART算法，SIRT能夠使重建圖像更加平滑，更好地壓制帶狀偽影.SIRT算法的迭代公式為：

(4)

式中:λ為松弛因子(0<λ<2)，且λ的選取根據(jù)投影數(shù)的多少和噪聲的大小而有所不同.在迭代過(guò)程中，需要尋求適當(dāng)?shù)摩藖?lái)提高收斂速度和重建質(zhì)量[5].

在方程組(1)中，每一個(gè)方程表示一個(gè)超平面.每按照式(4)進(jìn)行一次迭代，表示要利用所有的射線對(duì)某些像素值進(jìn)行修正，每修正一次即對(duì)重建圖像更新一次，不斷地對(duì)圖像進(jìn)行校正，直到滿足所需要求，然后結(jié)束迭代過(guò)程.

2 非線性擴(kuò)散模型

2.1 PM模型

由Perona和Malik提出的非線性擴(kuò)散模型—PM模型[4]，是一種常用的去噪方法，表達(dá)式如下：

(5)

式中：div為散度算子；?為梯度算子；?u=[ux，uy]T；|?u|表示圖像梯度的模；c(·)為擴(kuò)散系數(shù)；u0為原始圖像；t為時(shí)間變量.

c1(|?u|)=1/1+(|?u|/k)2

(6)

c2(|?u|)=exp[-(|?u|/k)2]

(7)

式中:k為擴(kuò)散門(mén)限，且該參數(shù)的大小由具體的應(yīng)用情況確定，k越大，平滑的力度就越大.

根據(jù)擴(kuò)散系數(shù)c(·)的性質(zhì)可知，在|?u|較小的均勻區(qū)域，為了去除噪聲，對(duì)圖像的平滑作用應(yīng)增強(qiáng)；而在|?u|較大的邊緣附近，對(duì)圖像的平滑作用應(yīng)“停止”，以保護(hù)邊緣細(xì)節(jié)信息.

2.2 PM模型的優(yōu)缺點(diǎn)分析

PM模型是一種基于PDE的非線性擴(kuò)散模型，在對(duì)圖像去噪的過(guò)程中，有以下2個(gè)優(yōu)點(diǎn)[6-7]：

① 保留邊緣細(xì)節(jié).在邊緣附近所對(duì)應(yīng)的c(·)趨于0，擴(kuò)散速度較低，故在邊緣附近，平滑之前像素的灰度值對(duì)平滑之后的影響很小.因此，利用PM模型去噪可以保護(hù)圖像的邊緣.

② 平滑圖像噪聲.在均勻區(qū)域內(nèi)所對(duì)應(yīng)的c(·)趨于1，擴(kuò)散速度較高，故該區(qū)域內(nèi)的像素值做一次平滑后，對(duì)圖像的影響不大.因此，對(duì)于含有噪聲的像素點(diǎn)，利用PM模型可以進(jìn)行有效地平滑.

PM模型有以下3個(gè)缺點(diǎn)[6-7]：

① 有零散斑點(diǎn)的出現(xiàn).因?yàn)镻M模型在去除均勻區(qū)域內(nèi)的孤立噪聲點(diǎn)時(shí)，會(huì)把它們誤認(rèn)為圖像的邊緣而予以保留，這就使得去噪后的圖像中有零散的斑點(diǎn)出現(xiàn).

② 平滑效果一般.若像素p是一個(gè)孤立的噪聲點(diǎn)，如果該點(diǎn)與鄰域中每個(gè)像素q相比，都有較大的灰度值差距，那么對(duì)應(yīng)的梯度一定很大，使得擴(kuò)散速度降低，不利于去除該點(diǎn)的孤立噪聲.

③ 擴(kuò)散系數(shù)是病態(tài)的.因?yàn)閏(·)是趨向于0，但達(dá)不到0，所以給出的擴(kuò)散系數(shù)式(6)和式(7)是病態(tài)的，也就是說(shuō)在圖像的邊緣附近，即使|?u|很大，但由于c(·)≠0，使得這些像素依然會(huì)受到鄰域像素的影響.盡管影響很小，但是隨著迭代次數(shù)的增加會(huì)被放大，這就會(huì)模糊甚至是破壞圖像的邊緣細(xì)節(jié).

對(duì)于PM模型所存在的問(wèn)題，為了更好地降噪和保持邊緣細(xì)節(jié)能力，需對(duì)PM模型進(jìn)行改進(jìn).

2.3 PM模型的改進(jìn)

本文中，針對(duì)PM模型去噪時(shí)會(huì)有零散斑點(diǎn)出現(xiàn)、平滑效果一般的缺點(diǎn)，結(jié)合Catté_PM模型的特點(diǎn)，通過(guò)引進(jìn)高斯算子Gσ來(lái)實(shí)現(xiàn).也就是說(shuō)利用Gσ對(duì)噪聲圖像進(jìn)行高斯平滑，即Gσ*u，這就大大降低了孤立噪聲點(diǎn)處的鄰域像素對(duì)其梯度的影響，且增大擴(kuò)散系數(shù)，有利于對(duì)孤立噪聲點(diǎn)的去除.

其次，針對(duì)擴(kuò)散系數(shù)c(·)是病態(tài)的問(wèn)題，本文通過(guò)對(duì)c(·)進(jìn)行改進(jìn)，構(gòu)造出一個(gè)新的擴(kuò)散系數(shù)：

(9)

綜上所述，改進(jìn)后的PM模型可表示為：

(10)

3 PM-SIRT重建算法

由于傳統(tǒng)的SIRT算法并不能去除投影數(shù)據(jù)及重建圖像中的噪聲[8]，故本文提出了一種改進(jìn)算法，主要是將圖像去噪與圖像重建相結(jié)合，在SIRT算法每次迭代后，利用PM模型對(duì)重建圖像去噪，使得SIRT算法與PM模型可以循環(huán)交替進(jìn)行.該算法的具體步驟如下：

①對(duì)投影數(shù)據(jù)P進(jìn)行加噪處理；

②初始化各參數(shù)，令F0=0，設(shè)定初始迭代值j=1和最大迭代次數(shù)maxiter；

③利用初始值F0及投影數(shù)據(jù)P進(jìn)行SIRT迭代，得到重建圖像Fj；

⑥j=j+1：若j≤maxiter，返回步驟③繼續(xù)進(jìn)行迭代，直到滿足停止迭代條件.

4 仿真實(shí)驗(yàn)結(jié)果與分析

為了驗(yàn)證PM-SIRT算法的有效性，仿真實(shí)驗(yàn)中使用180*180的Shepp-Logan模型進(jìn)行圖像重建，如圖1所示.本文中，針對(duì)投影數(shù)據(jù)中加入高斯白噪聲(均值為0，標(biāo)準(zhǔn)差分別為0.5和1)的情況進(jìn)行仿真.其中，設(shè)定最大迭代次數(shù)maxiter=10，σ=5，通過(guò)反復(fù)實(shí)驗(yàn)、對(duì)比，選取參數(shù)k=20.

圖2給出了SIRT、SIRT+PM和PM-SIRT算法所重建的圖像，其中SIRT+PM算法的具體流程與PM-SIRT算法一致，但其擴(kuò)散系數(shù)選用式(7).

圖2 三種算法的重建圖像(從左到右依次為SIRT算法，SIRT+PM算 法和PM-SIRT算法

從圖2可以看出，在加入標(biāo)準(zhǔn)差為0.5和1的高斯白噪聲的情況下，相比于SIRT+PM算法、SIRT算法，PM-SIRT算法能夠有效地減少了噪聲和偽影且平滑性較好，重建出的圖像與原始圖像更加接近.

三種算法所重建出的圖像在第90列的密度曲線如圖3所示，其中圖3(a)為高斯白噪聲的標(biāo)準(zhǔn)差為0.5時(shí)的密度曲線，圖3(b)為高斯噪聲的標(biāo)準(zhǔn)差為1時(shí)的密度曲線.由此可見(jiàn)，本文算法的重建質(zhì)量?jī)?yōu)于SIRT+PM算法和SIRT算法.

為了進(jìn)一步比較三種算法的重建效果，本文中利用均方誤差(MSE)、信噪比(SNR)和峰值信噪比(PSNR)來(lái)評(píng)估重建圖像的質(zhì)量，衡量重建圖像和原始圖像的逼近程度，計(jì)算結(jié)果如表1所示.

圖3 3種算法的密度曲線圖

表1 三種量化評(píng)估

續(xù)表

由表1可以看出，PM-SIRT算法的MSE最小，而SNR和PSNR最大.表明本文算法的重建效果較好，更加逼近Shepp-Logan模型圖像.

5 結(jié)束語(yǔ)

基于PM模型的優(yōu)缺點(diǎn)分析，本文提出了改進(jìn)方法，構(gòu)造一個(gè)新的擴(kuò)散系數(shù)，將SIRT算法與PM模型相結(jié)合，使之循環(huán)交替進(jìn)行，從而給出了基于PM模型的改進(jìn)SIRT圖像重建算法.實(shí)驗(yàn)結(jié)果表明，本文算法的信噪比和峰值信噪比更高，誤差更小，進(jìn)一步提高了重建圖像的質(zhì)量，有效地減少了噪聲對(duì)圖像的影響，使重建圖像更加逼近于原始圖像，獲得較好的視覺(jué)效果圖.