班曉征,李志華
(江南大學(xué) 物聯(lián)網(wǎng)工程學(xué)院,江蘇 無(wú)錫 214122)
圖像重建是層析成像技術(shù)的基礎(chǔ)和核心,它是一種利用投影信息來反演重建原始圖像特征的過程[1]。目前,圖像重建技術(shù)已經(jīng)被廣泛應(yīng)用于無(wú)損檢測(cè)、雷達(dá)成像、地質(zhì)勘探等多個(gè)領(lǐng)域,但是由于受多方面客觀因素的制約,投影信息往往是欠采樣的。人們一方面要求降低輻射劑量以減少對(duì)環(huán)境的污染和人體的傷害,另一方面復(fù)雜的應(yīng)用環(huán)境又導(dǎo)致可獲得采樣信息有限,因此欠采樣圖像重建問題的研究成為目前圖像重建領(lǐng)域的熱點(diǎn)[2]。
由于缺乏充分的投影數(shù)據(jù),欠采樣圖像重建通常被看作是一個(gè)病態(tài)反問題。針對(duì)這一問題,目前主要的解決方案有兩種,改進(jìn)的濾波反投影算法[3,4]和基于壓縮感知的迭代圖像重建算法[5]。改進(jìn)的濾波反投影算法考慮在保持重建速度優(yōu)勢(shì)情況下提升重建質(zhì)量,但是在數(shù)據(jù)量不充足時(shí),該算法仍然不能重建出質(zhì)量較高的圖像。
基于壓縮感知的迭代圖像重建算法考慮到圖像作為一種特殊的可稀疏信號(hào),通常使用壓縮感知方法提取圖像特征作為正則項(xiàng),從而優(yōu)化重建圖像的質(zhì)量。Yu等人[6]提出了(Total variation-algebraic reconstruction technique,TV-SART)算法,在傳統(tǒng)的統(tǒng)計(jì)ART迭代算法基礎(chǔ)上,引入了梯度約束,交替計(jì)算重建圖像和圖像梯度特征,從而重建出更多圖像細(xì)節(jié)。為了加速圖像重建的速度,劉石等人[7]將Tikhonov正則化算法應(yīng)用于圖像重建問題,提升了重建的實(shí)時(shí)性。文獻(xiàn)[8]將TV和L1結(jié)合在一起,為了獲取更多的梯度特征計(jì)算了八方向的全變分圖像重建模型,給出了基于改進(jìn)權(quán)重選取策略的TV軟閾值方法對(duì)L1-TV模型進(jìn)行求解,基于此提出了一種L1TV圖像重建算法??紤]到TV容易出現(xiàn)階梯效應(yīng),造成更多的噪聲污染,Ke等人[9]對(duì)重建后圖像使用兩級(jí)小波變換,然后再對(duì)小波系數(shù)進(jìn)行維納濾波自適應(yīng)去噪,并使用醫(yī)療 CT 圖像證明了所提出算法的可用性,但是所提出的重建算法效果易受小波閾值影響。
針對(duì)上述算法在稀疏采樣環(huán)境圖像重建方面存在的抗噪聲干擾能力較差的情況,本文提出了一種基于自適應(yīng)小波和TV相結(jié)合的圖像重建(Adaptive wavelet and TV,AWAT)算法。實(shí)驗(yàn)結(jié)果表明,在稀疏采樣環(huán)境中,所提出算法的圖像重建效果較好。
全變差正則化模型在有效去除噪聲的同時(shí),可有效保留圖像的邊緣信息,各向異性TV正則化圖像重建模型可表示為
(1)
式中:X為待重建圖像,DxX為圖像水平方向梯度,DyX為圖像垂直方向梯度。
然而全變差正則化模型本身存在“階梯效應(yīng)”,并且僅使用TV正則化模型并不能獲得圖像的最優(yōu)稀疏表示,重建圖像在邊緣處容易產(chǎn)生偽影。離散小波變換作為一種圖像稀疏表示和邊緣提取方法,可以有效地抑制重建過程中的條紋噪聲。因此向式(1)添加小波先驗(yàn)項(xiàng),得到
(2)
式中:W表示對(duì)圖像進(jìn)行小波變換,WX表示小波系數(shù)。
為了壓縮計(jì)算時(shí)間同時(shí)提高重建質(zhì)量,引入支集檢測(cè)計(jì)算加權(quán)因子T,對(duì)小波系數(shù)進(jìn)行進(jìn)一步稀疏,獲得一個(gè)自適應(yīng)加權(quán)小波模型
(3)
式中:T表示自適應(yīng)權(quán)重矩陣,可通過迭代支集檢測(cè)方法計(jì)算,并且T中的元素只有0或1兩種取值。如果探測(cè)到小波系數(shù)的邊緣處,那么將T的對(duì)應(yīng)位置設(shè)為0,如果不在小波系數(shù)的邊緣,那么T設(shè)為1。通過加權(quán),不僅可以在盡可能保留圖像邊緣的同時(shí),降低先驗(yàn)計(jì)算數(shù)據(jù)量,并且加權(quán)的L1范數(shù)更加逼近L0范數(shù),提升了計(jì)算的準(zhǔn)確性。
考慮圖像重建過程中噪聲的存在,得到混合正則化圖像重建模型
(4)
式中:第一個(gè)正則項(xiàng)表示自適應(yīng)加權(quán)小波系數(shù)項(xiàng),第二個(gè)正則項(xiàng)是TV正則化約束項(xiàng),第三個(gè)正則項(xiàng)表示圖像重建過程中的噪聲。
由于本文所提出的圖像重建模型(即式(4))中有多個(gè)未知量,因此考慮使用增廣拉格朗日方法將式(4)轉(zhuǎn)變?yōu)闊o(wú)約束問題
(5)
(6)
迭代支撐檢測(cè)(Iterative support detection,ISD)[10]是一種增加一維稀疏信號(hào)稀疏性的方法,本文將面向一維信號(hào)的迭代支撐檢測(cè)理論擴(kuò)充到二維圖像信號(hào),設(shè)計(jì)出基于ISD的自適應(yīng)加權(quán)權(quán)重計(jì)算方法,具體描述為算法1。
由于圖像本身不具備稀疏性,不能滿足迭代支撐檢測(cè)理論的前提,但是圖像通過小波變換得到的小波系數(shù)具有稀疏性,因此通過設(shè)置閾值θmax(|Ek(i,j)|),獲取小波系數(shù)矩陣的邊緣信息,并根據(jù)邊緣信息更新權(quán)重,通過步驟6的權(quán)重更新得到小波系數(shù)矩陣的支撐集。
(1)軟閾值算法更新
通過固定式(5)中其他變量,發(fā)現(xiàn)LA(C)問題是典型的加權(quán)L1范數(shù)約束最小化問題,LA(p)問題是普通的L1范數(shù)約束最小化問題,分別表示為
(7)
(8)
對(duì)式(7)和(8)應(yīng)用軟閾值收縮方法,得到
(9)
(2)硬閾值算法更新
LA(σ)問題是對(duì)L0范數(shù)約束的最小化問題,目標(biāo)函數(shù)為
(10)
硬閾值方法可以高效地計(jì)算出式(10)的解析解,得到
(11)
(3)快速傅里葉變換算法更新
在得到C、σ和p后,LA(X,C,σ,p)僅剩下變量X是未知的,此時(shí)更新求解重建圖像X的目標(biāo)函數(shù),可以寫為
(12)
由于式(12)是典型的帶有罰函數(shù)的最小二乘優(yōu)化式,對(duì)式(12)求解可以通過對(duì)X求偏導(dǎo)的方式得到
(13)
使用快速傅里葉變換方法降低式(13)的求解計(jì)算量,寫作
X=F-1·
(14)
式中:F表示傅里葉變換,F-1表示傅里葉逆變換。
由2.2節(jié)的推導(dǎo),本文提出一種欠采樣圖像重建算法,將其稱為AWAT(Adaptive waveletand TV)算法,算法框架見算法2。
實(shí)驗(yàn)中用于對(duì)比的算法選取了FBP[11]算法、LSQR[12]算法以及近年來提出的BTV[13]算法和L1TV[8]算法。為了更好地模擬實(shí)際CT工作環(huán)境,本文向投影數(shù)據(jù)中添加了10%高斯隨機(jī)噪聲,實(shí)驗(yàn)中的初始化參數(shù)如表1所示。
表1 實(shí)驗(yàn)參數(shù)
實(shí)驗(yàn)使用圖1所示的像素大小為256×256的Shepp-Logan體膜、Forbild-Abdomen體膜以及真實(shí)的肺部CT圖像。為了體現(xiàn)算法的應(yīng)用價(jià)值,實(shí)驗(yàn)分為標(biāo)準(zhǔn)樣本實(shí)驗(yàn)和真實(shí)肺部CT樣本實(shí)驗(yàn)。對(duì)于標(biāo)準(zhǔn)樣本實(shí)驗(yàn)組,僅在稀疏投影角度下進(jìn)行實(shí)驗(yàn);而在真實(shí)CT樣本實(shí)驗(yàn)組中,考慮到實(shí)際應(yīng)用情況的限制,在有限和稀疏投影角度兩種采樣環(huán)境下進(jìn)行實(shí)驗(yàn)。兩組實(shí)驗(yàn)采用峰值信噪比(Peak signal-to-noise,PSNR)、結(jié)構(gòu)相似性(Structural similarity index,SSIM)以及條紋指標(biāo)(Streak index,SI)作為定量評(píng)價(jià)指標(biāo)[14],同時(shí)使用輪廓線方法作為主觀視覺評(píng)價(jià)對(duì)比方法。
稀疏實(shí)驗(yàn)設(shè)置的采樣空間在0°~180°范圍內(nèi),每隔10°進(jìn)行一次采樣。圖3和圖4展示了各算法重建Sheep-Logan體膜和Forbild-Abdomen體膜在稀疏采樣角度下的重建結(jié)果。
由圖2和圖3可以看出,FBP算法重建的圖像存在很多條形偽影,這是由于FBP是在頻域上直接通過反radon變換重建的圖像,在稀疏投影角度下缺少采樣數(shù)據(jù),保留了投影的痕跡。由于沒有正則項(xiàng)約束,LSQR算法重建圖像中也有斷續(xù)的條形偽影。由于TV正則化約束平滑了一部分圖像噪聲,BTV算法重建圖像的細(xì)節(jié)比LSQR更清晰。L1TV算法重建的圖像完全沒有圖像偽影,并且圖像細(xì)節(jié)更加清晰,這是由于L1TV算法基于雙約束正則化模型,TV模型平滑了圖像噪聲,同時(shí)L1正則化約束消除了部分偽影。所提出的AWAT算法重建的圖像更加清晰,基本已經(jīng)完全重建出了原始圖像的細(xì)節(jié),完整保留了圖像的邊緣。圖4和表1反映了各算法重建仿真數(shù)據(jù)時(shí)各主觀、客觀評(píng)價(jià)指標(biāo)的對(duì)比結(jié)果。從圖4的輪廓線橫線剖面圖可以看出,所提出的AWAT算法重建圖像的像素灰度值與原圖像非常接近,擬合度很高。
在客觀指標(biāo)方面,本文AWAT算法的PSNR、SSIM指標(biāo)最高,同時(shí)SI最低,說明AWAT算法重建的圖像與原圖像更接近,并且具有比較少的條紋噪聲。
表2 客觀評(píng)價(jià)參數(shù)對(duì)比
考慮到在實(shí)際醫(yī)學(xué)診斷采樣環(huán)節(jié)經(jīng)常會(huì)出現(xiàn)低劑量情況,因此在稀疏投影角度和有限投影角度下分別進(jìn)行重建實(shí)驗(yàn),結(jié)果如圖5所示。對(duì)圖5觀察可知,FBP算法重建結(jié)果中的白噪聲非常嚴(yán)重,LSQR算法和BTV算法重建結(jié)果好于FBP算法,重建圖像和背景肉眼可區(qū)分開,但是LSQR重建圖像背景中還有條紋噪聲,BTV算法本身不夠純凈。L1TV算法和AWAT算法重建情況類似,并且AWAT算法要優(yōu)于L1TV算法,肺部陰影層次更明顯。
同樣地,得到表3所示的對(duì)客觀評(píng)價(jià)指標(biāo)的對(duì)比,以及圖6 SSIM變化曲線判斷各算法迭代收斂情況。表3中所有算法重建真實(shí)CT數(shù)據(jù)和仿真CT數(shù)據(jù)時(shí)所需要的時(shí)間基本相似,但是SI指標(biāo)整體略高于仿真數(shù)據(jù),說明此時(shí)各算法重建的圖像中具有比較明顯的條紋噪聲。BTV算法雖然SSIM指標(biāo)比較高,但是PSNR和SI指標(biāo)甚至低于LSQR算法。L1TV和AWAT算法的結(jié)構(gòu)相似性指標(biāo)非常高。
圖6中,在開始迭代重建時(shí),LSQR算法的SSIM指標(biāo)隨迭代次數(shù)增長(zhǎng),隨后保持穩(wěn)定并略有下降,這是由于LSQR具有比較明顯的“半收斂”性質(zhì);BTV算法中的TV正則化約束隨迭代次數(shù)增長(zhǎng)不能發(fā)揮穩(wěn)定的作用,呈現(xiàn)與LSQR相似的現(xiàn)象;L1TV算法在整個(gè)迭代過程中,SSIM指標(biāo)曲線始終保持增長(zhǎng)趨勢(shì),并且高于前兩種算法;AWAT算法雖然由于設(shè)計(jì)的自適應(yīng)權(quán)重的原因,SSIM指標(biāo)有些波動(dòng),但最終收斂值最穩(wěn)定且最接近1。
表3 客觀評(píng)價(jià)參數(shù)對(duì)比
針對(duì)欠采樣圖像重建過程中容易出現(xiàn)重建精度不高的缺點(diǎn),本文提出了一種基于自適應(yīng)加權(quán)小波和TV正則化模型的AWAT欠采樣圖像重建算法。在仿真實(shí)驗(yàn)中,與其他經(jīng)典圖像重建算法如FBP、LSQR,以及較新的算法BTV、L1TV進(jìn)行對(duì)比,仿真實(shí)驗(yàn)結(jié)果表明,本文算法可以有效地降低TV算法中的條紋噪聲,同時(shí)重建圖像與原圖像有著較高的結(jié)構(gòu)相似性,迭代時(shí)間無(wú)太大差異,具有比較高的重建質(zhì)量和重建速度。