魏志晴,鄭文康,白艷萍,譚秀輝,程 蓉,胡紅萍,王 鵬

(中北大學(xué) 數(shù)學(xué)學(xué)院,山西 太原 030051)

0 引 言

傳統(tǒng)計算機斷層成像算法可分為解析法和迭代法。解析法由于投影數(shù)據(jù)等角度間隔缺失,重建圖像中會引入條紋狀偽影[1],且圖像細節(jié)嚴重模糊,重建圖像的噪聲較大,無法應(yīng)用于實際診斷。迭代法可以將先驗信息融合到重建圖像過程當(dāng)中,具有較好的重建效果。斷層成像中應(yīng)用最廣泛的迭代算法是代數(shù)重建技術(shù)(Algebraic Reconstruction Technique,ART)、聯(lián)合代數(shù)重建技術(shù)(Simultaneous Algebraic Reconstruction Technique,SART)和期望最大化(Expectation Maximization,EM)。稀疏角度背景下,一般選取迭代重建法進行圖像重建。近年來隨著壓縮感知理論的提出,基于稀疏先驗的正則化方法得到廣泛關(guān)注,目前已經(jīng)成為主流[2]。其核心思想是利用重建圖像的稀疏先驗信息來提高CT圖像重建質(zhì)量,在解決稀疏角度問題中獲得了較好的效果。主要包括基于局部平滑的方法、基于字典的正則化方法和基于非局部相似的方法。

基于局部平滑的方法,主要是利用全變分(Total Variation,TV)模型對圖像進行局部平滑。Candes等人[3]最先提出了TV正則化模型,此模型在圖像處理等領(lǐng)域可以有效抑制噪聲并保留圖像的邊緣、紋理等。Sidky等人[4]在Candes等人的研究成果上,設(shè)計了基于約束TV最小化的CT圖像重建模型,在稀疏角度采樣背景下獲得了較好的重建效果。但是TV正則化方法對于具有更多圖像紋理細節(jié)和復(fù)雜邊緣結(jié)構(gòu)的圖像,會使得重建的圖像過度平滑,導(dǎo)致部分細節(jié)丟失。

基于字典的正則化方法,如離散余弦變換等,其單一基函數(shù)無法滿足圖像信號的自適應(yīng)稀疏表示,因此在此基礎(chǔ)上出現(xiàn)了許多改進方法。Chen等人[2]提出了一種基于TV與字典學(xué)習(xí)方法的磁共振成像算法,在保持精細結(jié)構(gòu)方面具有比較好的效果,進一步提高了重建圖像質(zhì)量。Fan[5]提出了基于雙字典的稀疏角度MRI圖像重建方法,不僅提升了圖像質(zhì)量,而且收斂速度更快。但字典正則項的選取對重建圖像的質(zhì)量起到?jīng)Q定性作用,若選取不當(dāng),易造成某些特征被過度懲罰,不能保留微小的特征,使圖像細節(jié)模糊。

基于非局部相似方法主要是挖掘圖像在紋理等結(jié)構(gòu)信息具有重復(fù)性的特性。Buades等人[6]提出了著名的非局部均值(Nonlocal Mean,NLM)濾波算法。Huang[7]將NLM算法應(yīng)用在CT圖像ART算法的重建上,得到了較好的效果。崔樹輝[8]將BM3D (Block Matching and 3D collaborative filtering)加入到ART算法中,在稀疏角度情況下依舊取得了較高分辨率。Zhang[9]提出了群稀疏表示正則化(Group-Sparse Regularization,GSR)算法,具有良好的捕捉圖像奇異性的優(yōu)良特性,且單幅圖像中的稀疏性和非局部相似性是同時存在的,可以很好地達到對稀疏角度下SART重建圖像去除偽影的目的。

為了在平滑圖像的同時保留不同層次的結(jié)構(gòu),Zhang等人[10]提出一種用于保持邊緣的濾波器,即滾動引導(dǎo)濾波(Rolling Guidance Filter,RGF)。裴佩佩等人[11]將RGF與卷積稀疏表示相結(jié)合,應(yīng)用在圖像處理中,克服了邊緣模糊,保留了對比度和邊緣紋理信息。

受上述研究的啟發(fā),該文提出了一種稀疏角度下基于滾動引導(dǎo)濾波和群稀疏表示正則化的SART重建算法,稱為SART-GSR-RGF。該算法首先對稀疏角度采樣下的圖像進行FBP重建,作為SART算法的初始圖像;然后采用SART算法進行重建,再經(jīng)過群稀疏正則化處理;最后采用滾動引導(dǎo)濾波提升圖像的對比度,將所得圖像再次作為SART重建算法的輸入,繼續(xù)迭代重建,直到得到分辨率較高的圖像。這樣既利用了GSR的保護邊緣和平滑圖像的特點,又利用了RGF保留不同層次結(jié)構(gòu)同時傳遞圖像特征的特點。經(jīng)過實驗驗證,所提算法對于圖像分辨率的提升具有很大作用。

1 理論依據(jù)

1.1 FBP重建和SART重建

FBP(Filtered Back-Projection)重建算法的依據(jù)是中心切片定理[12]。即空域進行投影得到的一維傅里葉變換,是空域進行二維傅里葉變換的一個切片。通過取遍所有投影值,再進行傅里葉反變換得到空間域分布的圖像,就可以得到二維函數(shù)。FBP公式表示如下[13]：

(1)

上述公式表明,重建圖像f(x,y)在某一位置的值,是通過該點的所有濾波投影采樣的疊加。流程可以簡化為圖1。

圖1 FBP重建示意圖

SART算法[14]是CT圖像重建中的一類迭代重建算法,在ART和SIRT算法的基礎(chǔ)上,進行一次更新只用到了一個投影采樣角度下的投影數(shù)據(jù),也即該投影角度下所對應(yīng)的系統(tǒng)方程。在修正某個像素值之前,需要計算一下經(jīng)過像素的所有射線的誤差,以對像素進行修正,并進行加權(quán)和歸一化處理,再把上述的結(jié)果更新到該像素上去。重復(fù)此過程,直至滿足收斂條件。算法公式為：

(2)

圖2 SART算法示意圖

其中灰色區(qū)域部分即為SART、SIRT和ART算法步驟中變化的部分。就圖像重建處理時間而言,解析法比迭代算法具有更好的性能。就降噪和有限數(shù)據(jù)而言,迭代法優(yōu)于解析法,因此該文先用FBP對平行透射數(shù)據(jù)進行成像,之后作為SART迭代重建圖像的初值。

1.2 基于群的SR建模

圖像處理的稀疏表示(Sparse Representations,SR)就是將圖像信號變換到另一個變換域中,絕大部分變換系數(shù)的絕對值很小,所得到的變換向量是稀疏的或者近似稀疏的,可以將其作為圖像的簡潔表達。SR模型可以表示為：

(3)

其中,x表示觀測到的信號向量,D是字典,α表示信號的系數(shù),λ是正則化參數(shù),‖·‖0表示l0范數(shù)。SR的目標(biāo)是為訓(xùn)練好的D尋找一個稀疏向量α來表示x。為了更好地用α表示x,必須選擇一個有效的字典D。已經(jīng)提出了一些近似算法來交替優(yōu)化D和α,如MOD,k-SVD和在線學(xué)習(xí)。

1.3 滾動引導(dǎo)濾波

滾動引導(dǎo)濾波器具有尺度感知和邊緣恢復(fù)的優(yōu)良特性[15],能夠保留圖像邊緣,有效避免模糊現(xiàn)象。RGF首先對物體的尺度進行了篩選,小尺度圖像信息主要指的是細節(jié)、紋理、小物體和一些噪聲。大尺度圖像信息主要指的是邊界、顏色過渡和變化不明顯的區(qū)域。小尺度的邊緣被平滑掉,僅僅保留了大尺度物體的邊緣。大多數(shù)邊緣保持濾波器沒有考慮到在濾波器設(shè)計中融入尺度,不能很好地實現(xiàn)小結(jié)構(gòu)與大結(jié)構(gòu)之間的分離。因此,RGF對于雙邊濾波和非線性擴散濾波具有很大優(yōu)勢,可以用來解決偽影噪聲問題。

2 SART-GSR-RGF算法

受群稀疏正則化圖像重建[16]研究的啟發(fā),該文提出一種基于群稀疏性正則化與滾動引導(dǎo)濾波算法相結(jié)合的SART-GSR-RGF重建方法。

圖3 相似群構(gòu)造示意圖

GSR算法對每個相似群fGk尋找合適的稀疏字典DG,在確定稀疏字典后,

(4)

其中,αG稱為結(jié)構(gòu)GSR,在GSR算法中把αG作為正則項,則求解GSR的無約束優(yōu)化模型為：

(5)

此正則化模型利用文獻[9]來求解。

最后,利用具有尺度感知和邊緣恢復(fù)的優(yōu)良特性的RGF進行處理,先要進行小結(jié)構(gòu)去除。

當(dāng)用小尺度的高斯濾波去除小結(jié)構(gòu)的邊緣時,會使得圖像邊緣被嚴重模糊。相反,當(dāng)經(jīng)過相同方式用同樣尺度的高斯濾波去除大結(jié)構(gòu)的邊緣時,只是達到了模糊的效果,仍能復(fù)原出為視覺所能接受的圖像。利用加權(quán)平均的形式來表示高斯濾波器算子,可得到：

(6)

其次進行邊緣恢復(fù)。在雙邊濾波器基礎(chǔ)上將Ft作為引導(dǎo)圖像,迭代過程如式7所示：

(7)

算法1：SART-GSR-RGF

重復(fù)

SART重建：

forp=1,2,…,P

ifp=1

up=SART(u0,ω)

else

up=SART(up-1,ω)

end if

end for

ifup<0

up=0,u=up

GSR正則化：

fort=0,1,…,T

updateft+1,et+1=ft+1-bt

for eachfGk

end for

end for

u0=RGF(DG*αG)

滿足條件時停止

3 實驗結(jié)果

本節(jié)將對提出算法的重建性能進行驗證,測試平臺如下：64位Intel Core i7-10700處理器,主機頻率2.90 GHz,8G運行內(nèi)存的Windows10操作系統(tǒng),數(shù)學(xué)軟件為Matlab R2021a。實驗所用模型為Shepp-Logan。在比較不同算法的重建性能之前,首先利用所提算法進行消融實驗,即對SART算法、SART-GSR算法SART-RGF算法以及所提算法進行對比分析。同時驗證了在不同稀疏角度下各算法的性能,觀察重建效果。圖4和圖5展示了在256個換能器背景下32角度和64角度時不同算法的重建效果;圖6和圖7展示了在512換能器背景下32角度和64角度時不同算法的重建效果,相對應(yīng)地也展示了同一放大部分的重建結(jié)果。

圖4 256換能器32角度下圖像重建效果比較

圖5 256換能器64角度下圖像重建效果比較

圖6 512換能器32角度下圖像重建效果比較

圖7 512換能器64角度下圖像重建效果比較

通過圖4至圖7可以看出,所提算法具有較好的重建效果,但僅通過視覺效果進行評估算法的重建效果時,評估結(jié)果具有一定的主觀性。為了更客觀地評價所提算法的有效性,使用具體的評價指標(biāo)進行評價。選擇PSNR[17],MSE以及FSIM進行客觀評價,其中圖8、圖10和圖12展示了256換能器設(shè)置下32/64稀疏角度的20次迭代過程中的PSNR值、MSE值和FSIM值的比較結(jié)果及趨勢;圖9、圖11和圖13展示了512換能器設(shè)置下32/64稀疏角度的20次迭代過程中的PSNR值、MSE值和FSIM值的比較結(jié)果及趨勢。

圖8 256換能器圖像重建效果PSNR比較

圖9 512換能器圖像重建效果PSNR比較

圖10 256換能器圖像重建效果MSE比較

圖11 512換能器圖像重建效果MSE比較

圖12 256換能器圖像重建效果FSIM比較

圖13 512換能器圖像重建效果FSIM比較

根據(jù)圖8至圖13可以很清晰地看到,所提算法SART-GSR-RGF具有最高的PSNR值和FSIM值,以及最低的MSE值。其中PSNR值越大,意味著圖像重建后的保真效果越好,圖像越清晰;FSIM值越高,意味著當(dāng)前所測試的圖像越接近參考圖像;MSE值越小,意味著與原始重建圖像相差越低,具有較好的重建性能。并且從圖中可以看到,所提算法在迭代開始前期就已經(jīng)達到了比較好的性能,說明所提算法具有很好的收斂性能,即算法收斂速度快。

為比較各類方法的優(yōu)劣性,驗證所提算法的有效性,表1給出了SART算法、TV-POCS算法、SART-NLM算法以及SART-GSR-RGF算法的PSNR,RMSE以及FSIM的數(shù)據(jù)。結(jié)果表明,所提算法在重建結(jié)果上均優(yōu)于其他4個對比算法。

表1 不同算法重建性能比較

4 結(jié)束語

該文提出了一種稀疏角度下基于滾動引導(dǎo)濾波和群稀疏表示正則化的SART重建算法,稱為SART-GSR-RGF。由于群稀疏正則化在去掉偽影的同時,可能將邊緣或細節(jié)過度平滑,使得對比度降低,無法獲得符合人類視覺效果的高分辨率圖像,因此利用滾動引導(dǎo)濾波進行一定的對比度提升,再次作為SART的輸入進行迭代。經(jīng)過算法的分析以及實驗的驗證表明,所提算法具有良好的重建效果。接下來,將對此算法進行進一步研究,以期獲得三維頭骨模型的較高質(zhì)量的重建。