霍其潤 李建武 陸 耀 秦 明

計算機斷層成像(Computed tomograpy,CT)技術(shù)作為一種高性能的無損檢測手段在醫(yī)療等多種領(lǐng)域得到了廣泛的應(yīng)用[1?2].受實際成像系統(tǒng)的限制,重建后的CT 圖像中通常會存在一些環(huán)形偽影,環(huán)形偽影的產(chǎn)生原因復(fù)雜,如何處理CT 圖像中的環(huán)形偽影引起了業(yè)界廣泛的關(guān)注.一部分研究[3?4]在數(shù)據(jù)采集時通過校正各探測器的響應(yīng)率或是改變設(shè)備的掃描方式來避免環(huán)形偽影的出現(xiàn),該類方法實現(xiàn)時存在一定的局限性和難度,且重建后生成的CT 圖像中仍有殘留偽影的可能.

近年來,基于圖像處理的環(huán)形偽影校正方法研究取得了較大的進展,這類方法通過軟件實現(xiàn),實現(xiàn)代價較小,同時也可以達到較好的校正效果.該類方法又可劃分為前處理方法和后處理方法兩種.前處理方法是指在系統(tǒng)采集到投影數(shù)據(jù)后,通過對重建前的投影正弦圖進行處理,來避免重建后CT 圖像中出現(xiàn)環(huán)形偽影.Raven[5]提出在投影正弦圖上使用傅里葉變換的方法.Münch 等[6]將小波分解和傅里葉變換下的低通濾波相結(jié)合,實現(xiàn)了將偽影信息與圖像信息進行更為嚴格的區(qū)分.Ashrafuzzaman等[7]提出了自適應(yīng)的可變窗口均值濾波用于正弦圖中偽影的校正.Hasan 等[8]提出了基于形態(tài)學(xué)濾波的偽影檢測和校正方法.Anas 等[9]通過對正弦圖中的偽影進行檢測實現(xiàn)分類校正.Rashid 等[10]將偽影根據(jù)不同的強度和寬度進行不同的處理.Kim等[11]通過計算投影圖像域中各列數(shù)據(jù)間的比率來估計出各探測單元的不一致性并求出校正因子加以平衡,來實現(xiàn)偽影的去除.近些年還出現(xiàn)了一些結(jié)合壓縮感知的方法作用在重建前的投影圖像上[12?13]或圖像重建的過程中[14]實現(xiàn)對偽影的抑制.雖然前處理方法的研究取得了很大的進展,但由于方法本身需要獲取原始的投影數(shù)據(jù),故其應(yīng)用會受到一定的限制.后處理方法就是對重建后的CT 圖像進行直接處理來實現(xiàn)偽影去除的效果.Sijbers 等[15]在CT 圖像中提取出偽影較為明顯的感興趣區(qū)域,然后在極坐標系下利用滑動窗口對偽影進行濾除.Prell等[16]提出了一種結(jié)合均值和中值濾波的高分辨率CT 圖像偽影去除算法,并通過比較說明將該算法作用于極坐標系下的重建圖像校正效果更佳.Wei等[17]將Münch 等所提方法[6]用于重建圖像上,在極坐標系下對圖像進行小波分解和傅里葉低通濾波來實現(xiàn)偽影的校正工作.Yan 等[18]提出了稀疏約束的模型用于極坐標系下的重建圖像來實現(xiàn)對偽影的校正.和作用在投影正弦圖的前處理方法相比較,后處理方法由于是直接作用于重建后的CT 圖像,校正效果更加直觀,且實現(xiàn)時無需使用掃描過程中采集的原始投影數(shù)據(jù),使用更加方便,應(yīng)用場合也不受限制,無論是在成像環(huán)節(jié)還是在后續(xù)的使用環(huán)節(jié),都可以使用基于重建圖像的后處理方法來進行偽影校正,具有更好的應(yīng)用前景.目前后處理的偽影校正研究受到了越來越多的關(guān)注[19?21],本論文工作也是從后處理的角度進行研究的.

目前較為成熟的校正方法中常用的一些頻域或空域濾波方式,在去除偽影的過程中不可避免會對圖像中的邊緣信息產(chǎn)生一定的影響,圖像一些原始的細節(jié)信息也可能會丟失,因此在不損失圖像信息的情況下進行環(huán)形偽影的去除仍然是一個很有挑戰(zhàn)性的任務(wù).近年來基于變分的相關(guān)方法在圖像去噪、復(fù)原等領(lǐng)域取得了很大進展[22?25],該類方法的優(yōu)勢是在去除噪聲或干擾的同時能夠較好地保護邊緣信息,故本文基于變分的思想,把環(huán)形偽影校正問題轉(zhuǎn)化為一個能量優(yōu)化的問題來進行建模和解決,盡可能緩解保持圖像信息和去除偽影之間的矛盾;同時考慮到CT 圖像的環(huán)形偽影并非普通意義上的噪聲,結(jié)合CT 圖像環(huán)形偽影的產(chǎn)生機理和特征表現(xiàn),設(shè)計有針對性的優(yōu)化模型來有效地去除環(huán)形偽影,從而得到更高質(zhì)量的CT 圖像.與大多數(shù)后處理方法類似,為了使處理更加簡單有效,我們構(gòu)造的模型和算法也是在極坐標系下對偽影進行分離的,本文的校正流程如圖1 所示.整個過程中,作為核心步驟的分離偽影操作是在極坐標系下使用本文構(gòu)造的變分模型對圖像進行處理的,為了減少坐標變換操作對整體圖像信息的影響,最后一步是在笛卡爾坐標系下從原始CT 圖像中對得到的環(huán)形偽影成分直接進行減除的,這樣就相當(dāng)于在整個過程中僅偽影信息部分受到了兩次坐標間變換操作的影響,而CT 圖像中其他信息幾乎沒受到影響.

圖1 本文的校正算法流程示意圖Fig.1 The flow chart of the proposed method

本文的主要貢獻如下:

1)分析環(huán)形偽影在極坐標系下的幾何特性,提出了更適合的模型梯度保真形式,盡可能降低處理過程中圖像其他細節(jié)信息的丟失.

2)分析并提出環(huán)形偽影與圖像結(jié)構(gòu)性信息的邊緣特性差異,模型中引入對偽影信息更有針對性抑制作用的相對全變分約束形式,從而增強模型對圖像結(jié)構(gòu)性信息的保護.

本文后續(xù)內(nèi)容安排如下:第1 節(jié)詳細描述了基于變分方法的思想并結(jié)合CT 圖像環(huán)形偽影的結(jié)構(gòu)特性來構(gòu)造偽影校正模型的過程;第2 節(jié)介紹了模型的優(yōu)化求解思路及其偽影校正算法;第3 節(jié)分別在模擬數(shù)據(jù)和真實CT 圖像上進行了實驗和比較,并對實現(xiàn)效果進行了進一步分析;第4 節(jié)做出了最后的總結(jié).

1 變分模型設(shè)計

變分方法的思想是將一幅二維圖像看作是一個能量系統(tǒng),當(dāng)沒有噪聲或其他干擾信息時,認為圖像是平滑的,其對應(yīng)的能量也很小.因而在解決圖像干擾問題時,可以建立一個能量函數(shù)模型,引入某些約束條件,通過求取這個能量函數(shù)的最小值來獲取理想圖像的估計,從而去除圖像中的噪聲和干擾.由Rudin 等首次提出的一種經(jīng)典全變分模型[26],又稱ROF 模型,形式如下:

其中,(x,y)為像素點,I為初始的含噪聲圖像,S為處理后的平滑圖像,? 代表圖像域空間,式中采用了全變分(Total variation,TV)的約束形式:

如式(1)所示,變分模型通常由保真項、正則項和正則化參數(shù)構(gòu)成.模型中的第一項為保真項,其作用是保證處理后圖像與帶噪聲的初始圖像之間不要有太大的差別.第二項為正則項,也叫平滑項或約束項,其作用是去除噪聲,使圖像變得更加平滑,但同時也會使得圖像信息在一定程度上有所損失.λ是用于平衡保真項和正則項的正則化系數(shù),實際使用時通過合理選取正則化系數(shù)λ的值來達到圖像去噪和保持圖像信息之間的平衡.

對于CT 圖像而言,環(huán)形偽影可以被看作是一種噪聲污染,但是一種比較特殊的噪聲形式,直接采用傳統(tǒng)變分模型來處理顯然不太適合,為了能夠更有針對性地抑制環(huán)形偽影,我們結(jié)合環(huán)形偽影在實際CT 圖像中的結(jié)構(gòu)特性,來設(shè)計更為合理的能量函數(shù)模型,從而有效地實現(xiàn)對CT 環(huán)形偽影的校正.

1.1 環(huán)形偽影的幾何特性分析及梯度保真設(shè)計

在笛卡爾坐標系下的CT 重建圖像中,環(huán)形偽影表現(xiàn)為一系列圓心與重建圖像中心相重合的同心圓環(huán),為了降低偽影的提取與處理難度,我們可將待處理的CT 重建圖像從笛卡爾坐標系轉(zhuǎn)至極坐標系下進行處理.在極坐標系下,同心圓環(huán)的環(huán)形偽影就變成了平行的垂直條紋,如圖2 所示.針對條帶模式的噪聲,許多不同領(lǐng)域也有著相關(guān)的研究,并根據(jù)自身特定的情況和應(yīng)用產(chǎn)生了各種各樣的方法[27?29].受這些方法的啟發(fā),我們對極坐標系下的偽影特性進行分析,設(shè)計更加有效的極坐標系下的去偽影模型.

在極坐標系下,包含偽影的CT 圖像可以用數(shù)學(xué)公式表示為:

其中,I表示轉(zhuǎn)至極坐標系下的CT 圖像,即含有偽影的待處理圖像,S表示極坐標系下不含偽影的理想圖像,n為極坐標系下CT 圖像中的偽影信息.如圖2(b)所示,在極坐標系下偽影表現(xiàn)出了明顯的垂直條狀,因此可以被視為是一種有結(jié)構(gòu)的特殊噪聲,其引起的圖像梯度變化主要表現(xiàn)在水平方向,而圖像在垂直方向上受到的梯度干擾幾乎為0,即

對式(3)分別在水平和垂直方向上求偏導(dǎo)可得:

其中,?x和?y分別代表水平和垂直方向的梯度算子,在此由于I、S和n均為極坐標系下的圖像信息,故式中x和y實際分別指示了圖2(b)中的徑向ρ和角度方向θ.結(jié)合式(4),故有:

由式(6)中描述的包含偽影圖像的梯度情況可知,條狀偽影很明顯地影響了圖像水平方向的梯度變化,對圖像垂直方向的梯度變化影響甚微,因此在構(gòu)造偽影校正模型時,我們可以利用這種鮮明的幾何特性.

圖2 極坐標轉(zhuǎn)換Fig.2 Polar coordinate transformation

式(1)所示的傳統(tǒng)ROF 模型中,第一項使用了圖像處理前后像素灰度差的最小二乘形式來作為保真項,根據(jù)CT 圖像的重建原理可知,在誤差相差不大的情況下,靠近圓心的環(huán)形偽影將表現(xiàn)得更為明顯,可見偽影的表現(xiàn)強度不滿足高斯分布,根據(jù)研究者提出在去除非高斯加性噪聲時,保真項采用L1范數(shù)形式比傳統(tǒng)L2范數(shù)形式更有效[30?31],因此保真項可改為L1范數(shù)形式.此外,通過對環(huán)形偽影成因和實際醫(yī)學(xué)CT 圖像的分析可知,僅當(dāng)CT 系統(tǒng)中出現(xiàn)了個別探測器通道損壞時,圖像中會出現(xiàn)少量強度較大的單環(huán)偽影,對于這種情況通常可以通過硬件修復(fù)或圖像檢測等方法實現(xiàn)有效的處理,這部分不屬于本文的研究范圍;而大多數(shù)情況下由于探測器響應(yīng)不一致等原因引起的多環(huán)偽影通常會具有相當(dāng)?shù)臄?shù)量,雖表現(xiàn)強度不會太大,但偽影的去除勢必會產(chǎn)生一定程度的圖像灰度信息變化.根據(jù)上述偽影分析和式(6)中的描述,即理想情況下去除偽影后的結(jié)果圖像對應(yīng)的垂直梯度信息應(yīng)與原圖的基本保持一致,因此在設(shè)計模型時,利用保真項來限制圖像處理前后垂直梯度信息盡可能不變更為合適,在此我們將常用的圖像灰度保真修改為圖像的單向梯度保真.第二項中用到的正則化約束是為了最大限度地去噪,普通的去噪處理要對圖像中水平和垂直方向的梯度均進行約束,針對極坐標下偽影的單一方向性,去偽影時僅考慮一個方向的梯度約束即可,在此只對水平方向進行梯度約束.由此,可以得出以下基于梯度保真的變分模型形式:

這里保真項和正則項中分別選用了垂直和水平兩個不同的方向,第一項保真項用來保持圖像在垂直方向的梯度信息,顯然這些梯度主要來自圖像本身的內(nèi)容,因此更大程度地保留了圖像的原始細節(jié);第二項正則項中只對水平方向的梯度進行了約束,從而更有針對性地去除垂直偽影.

1.2 環(huán)形偽影的邊緣特性分析及正則項設(shè)計

上述梯度保真模型(7)可以較好地保留垂直方向圖像的原始細節(jié)信息,其中的正則項僅對圖像水平方向的梯度進行了約束.但水平方向的梯度顯然不是完全由垂直偽影引起的,其中也包含了圖像本身內(nèi)容的信息,因此在做約束時,若能將兩者區(qū)分開,處理效果將會得到更有效的提升.

從環(huán)形偽影的成因和表現(xiàn)可知,連續(xù)多個探測器出現(xiàn)同樣的響應(yīng)異常是比較少見的,即環(huán)形偽影每個環(huán)的寬度有限,其灰度值與周圍像素存在一定的差異,表現(xiàn)為較細的邊緣形式.對于包含環(huán)形偽影的圖像,其上的邊緣信息有來自于圖像本身的結(jié)構(gòu)性信息,也有來自環(huán)形偽影的信息.按照邊緣像素與其周圍像素的灰度關(guān)系,圖像中的邊緣信息可大致分為兩類:屋脊型邊緣和階躍型邊緣.屋脊型邊緣位于灰度增加與減少的交界處;階躍型邊緣在階躍兩邊的灰度值有明顯的變化.圖3 為兩類邊緣剖面示意圖,分別對應(yīng)了兩類邊緣在垂直于邊緣方向上可能出現(xiàn)的灰度變化趨勢,其中橫坐標表示垂直于邊緣的方向,縱坐標表示像素值.環(huán)形偽影的每個環(huán)較細傾向為屋脊型邊緣;而圖像中的階躍型邊緣主要來自于圖像結(jié)構(gòu)性信息.相應(yīng)兩類邊緣產(chǎn)生的圖像梯度情況也有所不同,如圖4 所示.在邊緣所在的局部區(qū)域,一個屋脊型邊緣產(chǎn)生的圖像梯度既包含正向梯度也包含反向梯度;而一個階躍型邊緣產(chǎn)生的圖像梯度通常方向是一致的,要么全為正,要么全為負.

根據(jù)上述對不同類型邊緣的特性分析可以看出,環(huán)形偽影與圖像結(jié)構(gòu)性信息在產(chǎn)生圖像局部梯度的效果上表現(xiàn)有明顯的差異.在此要去除圖像中的環(huán)形偽影,若能將環(huán)形偽影與結(jié)構(gòu)性信息表現(xiàn)出的差異性引入到對圖像梯度的正則約束項中,可以更有針對性地對偽影信息進行抑制,從而減少處理過程中對其他圖像信息的影響.

圖3 邊緣剖面示意圖Fig.3 Edge profile diagram

圖4 兩類邊緣產(chǎn)生的梯度效果示意圖Fig.4 Gradient diagram of two types of edges

受Xu 等在研究圖像結(jié)構(gòu)–紋理分離方法時提出的相 對全變分(Relative total variation,RTV)概念[32]的啟發(fā),發(fā)現(xiàn)相對全變分的計算結(jié)果能較好地體現(xiàn)出上述分析的差異性.在極坐標系下偽影表現(xiàn)為平行的垂直條紋,基于模型(7)的構(gòu)造,第二項的正則項中只考慮約束圖像水平方向的梯度變化.在水平方向上,圖像中所有邊緣均會產(chǎn)生梯度信息,由于各條垂直偽影寬度有限,為屋脊型邊緣,其相關(guān)像素產(chǎn)生的局部水平梯度會有正有負,而圖像結(jié)構(gòu)性信息更多地傾向為階躍型邊緣,其相關(guān)像素產(chǎn)生的局部水平梯度方向是大體一致的.為突出對垂直偽影的約束,減少對圖像結(jié)構(gòu)信息的影響,我們結(jié)合相對全變分的定義,引入像素局部窗的單向相對全變分形式作為模型的正則項,實現(xiàn)偽影與圖像內(nèi)容更嚴格的分離,得到最終的目標函數(shù)如下所示:

其中,

式(8)中,ε取很小的正實數(shù)以免分母為0,R(p)表示圖像中以像素p為中心的局部區(qū)域窗口,q為該局部窗中的像素點,gp,q為根據(jù)區(qū)域內(nèi)各梯度所在位置定義的權(quán)重函數(shù),表示如下:

其中,σ控制了局部區(qū)域的窗口尺寸.Dx(p)是像素p在x方向的局部窗口全變分,它計算了R(p)中各點水平梯度絕對值的總和.無論是條狀偽影還是圖像的結(jié)構(gòu)性信息邊緣,都對應(yīng)較大的Dx.Lx用來度量局部窗口內(nèi)水平梯度的總體變化,在求和時未對窗口中各梯度單獨取絕對值,故求得的總和大小取決于區(qū)域內(nèi)各梯度是否具有一致的方向性.比較而言,傾向為屋脊型邊緣的條狀偽影更易使區(qū)域中的水平梯度出現(xiàn)有正有負的情況,傾向為階躍型邊緣的圖像結(jié)構(gòu)性信息會產(chǎn)生更多方向一致的水平梯度值,因此僅包含條狀偽影的區(qū)域,其對應(yīng)的Lx值通常小于包含圖像結(jié)構(gòu)性信息的區(qū)域.將Dx和Lx結(jié)合起來構(gòu)成的水平方向相對全變分Dx(p)/(Lx(p)+ε)能夠凸顯出圖像中的條狀偽影信息(條狀偽影信息Dx較大Lx較小;結(jié)構(gòu)性信息Dx較大Lx也較大;平滑區(qū)域Dx和Lx均較小),故將其作為最小化模型中的正則約束項,如式(8)所示,可更好地滿足僅對圖像中條狀偽影進行抑制的需求.

2 模型求解及其校正算法

式(8)所示的目標函數(shù)是非凸的,無法直接求解,由于二次項便于線性優(yōu)化,我們參照文獻[32]中給出的方法將公式中的RTV 約束項分解成一個非線性項和二次項的乘積形式,如下所示:

上述推導(dǎo)中的第三行為避免分母為0 添加了一個小的正實數(shù)εs,故為近似相等.可見約束項被分解為了二次項和非線性部分uxqwxq的乘積,其中

式(13)中求像素的ux時用到了周圍區(qū)域的梯度信息,Gσ表示標準差為σ的高斯核函數(shù),?是一個卷積運算符,式中的除法為點除運算.式(14)中求像素的wx時計算比較簡單,只涉及當(dāng)前像素的梯度信息.

此外,近似可將式(8)中第一項進行如下轉(zhuǎn)換:

基于以上轉(zhuǎn)換,式(8)相應(yīng)可寫為以下矩陣形式:

式中,R為一對角矩陣,對角線元素取值R[i,i]=1/(|?y(Si ?Ii)|+εs),vS和vI分別是圖像S和I的列向量表示形式,Cx和Cy分別是在水平方向和垂直方向做前向差分的梯度算子矩陣,Ux和Wx均為對角矩陣,它們對角線元素的值分別為Ux[i,i]=uxi,Wx[i,i]=wxi.

要對式(16)求最小值,令其對vS求導(dǎo)為0,得到如下線性方程:

故有

結(jié)合上述設(shè)計的模型及其求解過程,基于變分模型的CT 圖像環(huán)形偽影校正算法描述如下:

算法1.基于變分模型的環(huán)形偽影校正算法

輸入.含有環(huán)形偽影的CT 重建圖像F,參數(shù)λ和σ

輸出.去除偽影后的CT 圖像R

1)將圖像F轉(zhuǎn)至極坐標下生成圖像I

2)S(0)←I,t ←0,=1×10?3

3)repeat

基于S(t),利用式(13)和式(14)計算權(quán)值w和u,生成矩陣L(t);

基于L(t),利用式(18)計算,生成圖像S(t+1);

4)條狀偽影K=I ?S

5)將K轉(zhuǎn)至笛卡爾坐標系下生成環(huán)形偽影A

6)無偽影的結(jié)果圖像R=F ?A.

3 實驗結(jié)果與分析

基于本文所提算法,我們分別在模擬數(shù)據(jù)圖像和實際CT 圖像上進行了偽影校正的有效性驗證.實驗中待處理初始圖像均歸一化為[0,1]的灰度圖,在笛卡爾坐標系和極坐標系下的圖像尺寸分別設(shè)定為512×512 和360×360.為了驗證算法的有效性,我們在此選擇了三個有代表性的后處理方法—WF[17]算法、RCP[16]算法和VDM[18]算法作為對比,其中WF 算法雖最初是基于投影正弦圖設(shè)計的[6],由于作用原理相似,后續(xù)也被直接作用于轉(zhuǎn)換至極坐標系下的CT 重建圖像[17].由于這些方法都是在極坐標下對重建圖像進行處理的,公平起見,我們實驗中采用了統(tǒng)一的坐標轉(zhuǎn)換方式.

3.1 模擬數(shù)據(jù)實驗

為了客觀地評價各種算法的性能,我們采用一些常用的定量評估方法來評價去除偽影后的圖像效果,在此引入兩種常用的圖像評價指標:峰值信噪比(Peak signal-to-noise ratio,PSNR)和平均結(jié)構(gòu)相似性(Mean structural similarity,MSSIM).PSNR是最普遍和使用最為廣泛的一種圖像質(zhì)量評價指標,是基于對應(yīng)像素點間的誤差,將處理結(jié)果與參考圖像的像素值進行比較,處理結(jié)果的PSNR 值越大說明與參考圖像的像素值差異越小,即算法性能越好.MSSIM 是一種更接近于人類視覺的圖像質(zhì)量評測方法,指標值體現(xiàn)了處理結(jié)果與參考圖像間的亮度、對比度以及結(jié)構(gòu)的相似性,更多地與圖像視覺特征相關(guān),值越大說明處理結(jié)果視覺上越接近理想的參考圖像.考慮上述評價指標計算時均需用到理想的無偽影圖像做參考,而現(xiàn)實情況下,對于包含偽影的CT 圖像,這些理想的參考圖像往往無法獲得,故在此使用人工生成的一些模擬數(shù)據(jù)進行實驗.

我們首先采用的是Shepp-Logan 模型測試圖,如圖5 所示,其中包含許多均勻平滑的區(qū)域,便于對處理結(jié)果進行有效的觀察.圖5(a)為不包含偽影的標準參考圖像,將人工模擬的環(huán)形偽影疊加在參考圖像上生成了待處理的模擬圖像,如圖5(b).圖5(c)顯示了采用本文算法處理后的整體效果,圖5(d)～5(f)分別對應(yīng)圖5(a)～5(c)中相同部分區(qū)域的放大顯示效果,從中可以看出本文算法的處理效果視覺上非常接近理想的參考圖像.圖5(g)～5(i)分別是WF、RCP 和VDM 三種算法的處理結(jié)果,可以看出WF 算法的處理結(jié)果(圖5(g))中有新增偽影的出現(xiàn),RCP 算法的處理結(jié)果(圖5(h))中大部分偽影均被有效地去除了,但在圖像中心位置處,由于偽影強度較大,校正后仍有一定的殘留偽影存在.VDM 算法的處理結(jié)果(圖5(i))視覺效果較好,僅從視覺上已看不出殘留偽影的存在.除了視覺觀察外,為了更客觀地評價算法性能,我們基于參考圖像,計算了各種算法處理結(jié)果對應(yīng)的峰值信噪比PSNR 和平均結(jié)構(gòu)相似性MSSIM,如表1所示.通過比較可以看出,本文算法較WF 和RCP兩種算法,性能指標值提高明顯,這也與之前的視覺觀察的結(jié)果相一致;此外雖然本文算法與VDM 算法的處理結(jié)果在視覺比較中看不出太大的差異,但通過定量的性能指標計算,兩者還是有一定的差異的,本文算法的處理結(jié)果對應(yīng)了更高的峰值信噪比和平均結(jié)構(gòu)相似性,結(jié)合視覺觀察和定量指標的分析可以看出本文算法在去除偽影和保持圖像信息方面具有更好的性能.

我們的第二組模擬數(shù)據(jù)選取了Lena 圖像,如圖6 所示,雖然該圖像不是CT 圖像,由于圖像中既有平滑區(qū)域又包含豐富的邊緣及紋理細節(jié)信息,在此選用可以更有效地體現(xiàn)處理過程對多種圖像成分的影響效果.理想的無偽影圖像如圖6(a)所示,圖6(b)為添加了環(huán)形偽影的待處理模擬圖像,采用本文算法處理后的整體效果如圖6(c)所示,視覺上看來偽影已被有效去除,同時圖像中的細節(jié)信息也保留得很好.為了便于更清晰的視覺觀察,取圖6(a)～6(c)中相同部分區(qū)域進行放大顯示,依次對應(yīng)圖6(d)～6(f)所示.圖6(g)～6(i)分別是WF、RCP 和VDM 三種算法的處理結(jié)果,總體上來看,各種算法均有效地去除了圖像中的大部分偽影,同時也較好地保留了圖像中的細節(jié)信息.但在某些局部區(qū)域,WF 算法和RCP 算法的處理結(jié)果中仍可看出有少量偽影的殘留,如圖6(g)和圖6(h)中箭頭所指位置.VDM 算法的處理結(jié)果明顯好于WF和RCP 兩種算法,基本看不出有殘留偽影的存在,但和本文算法的處理結(jié)果比較起來,在一些灰度均勻的區(qū)域,如圖6(i)中箭頭所指的臉頰位置,其處理后的效果比較起來略顯得有些粗糙,而本文算法的結(jié)果看起來更加平滑.通過表1 中所列的評價指標值可以看出,和其他三種算法比較起來,本文算法的處理結(jié)果具有最高的PSNR 和SSIM 值.

圖5 Shepp-Logan 圖像處理結(jié)果Fig.5 Experimental results on the Shepp-Logan phantom

3.2 真實數(shù)據(jù)實驗

為了進一步驗證本文算法的實際校正效果,我們利用上述幾種校正算法對兩幅真實CT 圖像進行了去除偽影的處理,如圖7 和圖8 分別顯示了一幅腦部CT 和一幅頸部CT 的各算法實驗情況對比.

表1 各算法結(jié)果的圖像質(zhì)量評價指標值Table 1 Quantitative comparison for the different methods

圖6 Lena 圖像處理結(jié)果Fig.6 Experimental results on the Lena image

在這兩組實驗圖中,圖7(a)和圖8(a)為處理前的原始圖像,在此為了讓環(huán)形偽影看起來更明顯,我們適當(dāng)?shù)卦鰪娏藞D像的對比度.可以看出在接近圖像中心位置有較強的偽影,偏離中心的位置上偽影強度較弱.圖7(b)和圖8(b)分別是圖7(a)和圖8(a)中局部區(qū)域的放大效果,圖7(c)～7(f)和圖8(c)～8(f)依次放大顯示了四種算法處理結(jié)果中的對應(yīng)區(qū)域.圖7(c)為本文算法的處理結(jié)果,環(huán)形偽影沒有明顯的殘留,同時圖像結(jié)構(gòu)及平滑度均得到了較好的保留.圖7(d)和圖8(d)是WF 處理的效果,可以看出靠近圖像中心的環(huán)形偽影去除的較好,但在偏離中心的位置會有一些新的偽影出現(xiàn).圖7(e)和圖8(e)是RCP 處理的結(jié)果,圖中大多數(shù)偽影都被有效的去除,同時細節(jié)信息也基本得到保留.然而由于算法本身需要考慮偽影去除與細節(jié)保留之間的平衡,故一些強度較大的環(huán)形偽影還是會有些殘留的痕跡.在VDM 算法的處理結(jié)果中,如圖7(f)和圖8(f)所示,盡管沒有明顯的環(huán)形偽影存在,但本該平坦的區(qū)域看起來卻不夠平滑.由此可見,對于實際CT 圖像中環(huán)形偽影的校正,本文提出的算法也體現(xiàn)了更優(yōu)的性能,在去除偽影的同時能更好地保持原始圖像的質(zhì)量.

圖7 腦部CT 圖像處理結(jié)果Fig.7 Experiments on a brain CT image

圖8 頸部CT 圖像處理結(jié)果Fig.8 Experiments on a neck CT image

視覺觀察是檢驗偽影去除與圖像細節(jié)保護最直接的途徑,但要更精確地對比處理結(jié)果的差異,仍需客觀的評價指標來進行衡量.但對于實際待處理的CT 圖像,其對應(yīng)的理想無偽影圖像并不存在,故無法利用傳統(tǒng)的峰值信噪比或平均結(jié)構(gòu)相似性等指標來進行評價.在此我們在相對平滑的腦部CT 圖像中選取幾個平坦區(qū)域作為感興趣區(qū)域(Region of interest,ROI),如圖9 所示,計算這些區(qū)域的局部信噪比(Local signal-to-noise ratio,LSNR),即區(qū)域內(nèi)像素值的標準差與平均值之間的比值,公式如下:

對于圖像中的平坦區(qū)域,其標準差很大程度上與圖像中的偽影相關(guān),標準差的降低相應(yīng)地反映了圖像偽影的減少[16].顯然對于處理后的圖像,其平坦區(qū)域的LSNR 值越大,表明圖像偽影的去除效果越好.對應(yīng)圖9 各方框標識的區(qū)域位置,我們依次計算出原始圖像及其各算法處理結(jié)果中相應(yīng)區(qū)域的LSNR 值,如表2 所列.可以看出相對于原始圖像,各種算法處理后相應(yīng)區(qū)域的LSNR 值均明顯變大,其中本文算法較其他算法來說更為有效.

表2 各算法結(jié)果相應(yīng)局部區(qū)域(圖9)的LSNR 值Table 2 LSNRs of the ROIs circled in Fig.9 for different methods

3.3 實現(xiàn)分析

3.3.1 梯度保真與傳統(tǒng)灰度保真的效果對比

為了單獨說明校正模型中的保真項采用梯度保真形式在去除偽影時會產(chǎn)生更好的效果,我們以添加了偽影的Lena 圖像為例,將采用梯度保真形式的模型(7)與如下采用傳統(tǒng)灰度保真形式的模型(20)就示例圖像在極坐標系下的處理效果進行實際比較.

由于Lena 圖像中的邊緣及紋理細節(jié)信息更豐富,在此選用該圖視覺上可以更明顯地看出處理過程對多種圖像成分的影響效果.如圖10 所示,兩種模型的處理均明顯消除了偽影,但比較而言,使用灰度保真的模型(20)處理后的圖像細節(jié)損失相對嚴重,且有階梯現(xiàn)象出現(xiàn);而使用梯度保真的模型(7),處理結(jié)果中細節(jié)信息得到了更好的保留.

圖9 腦部CT 圖像中的ROI 選取Fig.9 A brain CT image with ROIs

為了進一步說明模型保真項的改進在實際CT圖像上的近似效果,同樣在極坐標系下對比了上述兩種模型(20)和(7)在包含細節(jié)信息相對較多的頸部CT 圖像上的處理效果,如圖11 所示.可以看出在消除了垂直偽影的同時,使用灰度保真的模型(20)處理后的圖像細節(jié)損失較多;而使用梯度保真的模型(7)處理后的結(jié)果中細節(jié)信息保留較好(圖中箭頭所指區(qū)域的效果對比更為明顯).

3.3.2 RTV 與TV 的約束效果對比

為了直觀顯示模型中正則項采用RTV 約束形式優(yōu)于傳統(tǒng)TV 形式的效果,我們將模型(7)與模型(8)的處理效果進行了對比,圖12 為Lena 圖像在極坐標系下的效果展示,顯然在單向梯度保真的作用下,兩者在消除垂直條狀偽影的同時均較好地保留了圖像中的細節(jié)信息.但仔細觀察不難發(fā)現(xiàn),在圖像中較為顯著的垂直結(jié)構(gòu)邊緣位置(如圖12(a)中箭頭所指),TV 約束模型(7)處理后其灰度和對比度均受到了一定的影響,而RTV 約束模型(8)的處理結(jié)果中看不出明顯的變化.為方便觀察,我們截取該位置附近的一小段水平線段(如圖12(d)中的

圖10 Lena 圖像偽影去除效果對比Fig.10 Comparison of artifact removal effects on the Lena image

圖11 頸部CT 圖像偽影去除效果對比Fig.11 Comparison of artifact removal effects on the neck CT image

圖12 Lena 圖像上TV 和RTV 的約束效果對比Fig.12 Comparison of constraint effects between TV and RTV on the Lena image

線段標識),將標準原圖、TV 約束模型處理結(jié)果圖和RTV 約束模型處理結(jié)果圖的相應(yīng)位置灰度信息以曲線形式顯示并加以對比,如圖12(e)所示,RTV 約束模型的結(jié)果與標準原圖更為接近.

同樣將上述兩種模型(7)與(8)作用在極坐標系下的頸部CT 圖像,處理效果對比見圖13,可以看出,兩種模型均基于單向梯度保真,故在消除垂直條狀偽影的同時均較好地保留了圖像中的細節(jié)信息.但仔細對比不難發(fā)現(xiàn),使用TV 約束模型(7)處理后,圖像部分區(qū)域(如圖12(b)中箭頭所指)的灰度和對比度出現(xiàn)了一些變化,而使用RTV 約束模型(8)的處理結(jié)果并沒有出現(xiàn)明顯的變化.

3.3.3 RTV 與TV 的約束效果對比

式(8)中的λ作為模型的約束項系數(shù),用來調(diào)整梯度保真程度和偽影去除程度的平衡,其值越大圖像中被去除的原始信息就越多,值越小去除的信息越符合條狀偽影的幾何特性,可較好地保留圖像原始信息.根據(jù)實際CT 圖像中環(huán)形偽影的情況,λ大小通常設(shè)為1×10?4～1×10?5較為合適,可在去除掉偽影信息的同時較好地保留圖像原始信息.

式(13)中尺度參數(shù)σ用于控制相對全變分計算時涉及到的局部窗口大小.σ值的大小取決于要去除偽影的寬度,它實際上是區(qū)分圖像結(jié)構(gòu)性邊緣和偽影邊緣的尺度標準.σ值較小時,僅可去除較細的偽影,比較寬的偽影則會被作為結(jié)構(gòu)性信息予以保留;σ值較大時,則可去除較寬的偽影,但同時也會影響到更多的圖像結(jié)構(gòu)信息.因此要根據(jù)實際偽影情況對σ的取值進行設(shè)定,通常為(0,6]范圍內(nèi),在上述實驗中我們設(shè)定值為3,同時在算法迭代時逐漸減小σ的值可以改善處理后的邊緣輪廓清晰度,每次迭代之后執(zhí)行σ ←max(σ/2,0.5).

此外式(13)中的ε值和式(14)、(15)中的εs值,這兩個參數(shù)設(shè)定為較小的正數(shù)主要是為了避免公式中出現(xiàn)分母為0 的情況.實驗中,我們設(shè)定ε值為1×10?3,εs值為2×10?2,εs設(shè)得稍微大些更有助于保留平穩(wěn)變化的結(jié)構(gòu)性信息.

圖13 頸部CT 圖像上TV 和RTV 的約束效果對比Fig.13 Comparison of constraint effects between TV and RTV on neck CT image

3.3.4 算法收斂性

算法1 中在對公式迭代求解來去除圖像偽影的步驟中,一般經(jīng)過第一次迭代即可明顯消除偽影,為了使偽影消除得更徹底,通常迭代3～5 次足以.我們通過實驗分析了該算法的收斂性,計算每次迭代后結(jié)果的歸一化能量變化(Normalized step difference energy,NSDE),公式如下:

以模擬數(shù)據(jù)Lena 圖像的為例,轉(zhuǎn)至極坐標系下利用式(8)進行處理時,在各次迭代后計算的NSDE 值所構(gòu)成的變化曲線如圖14 所示.從中可以看出隨著迭代次數(shù)的增加,NSDE 的值迅速下降,可見模型求解時能夠快速收斂.

圖14 算法迭代的收斂性趨勢Fig.14 Convergence of iteration algorithm

4 結(jié)論

本文針對CT 圖像常見的環(huán)形偽影問題,基于變分方法的思想,提出了一種作用于重建后CT 圖像的校正算法來有效地去除環(huán)形偽影.構(gòu)造校正模型時,結(jié)合醫(yī)學(xué)CT 圖像環(huán)形偽影的幾何特性及邊緣特性,引入了梯度保真形式和相對全變分約束形式,從而提高了對圖像信息的保真性,增強了對偽影抑制的針對性.同時通過有效的優(yōu)化求解,保證了校正算法的準確性和實用性.實驗驗證本文算法無論在模擬數(shù)據(jù)圖像還是真實的CT 圖像上均展示了較好的校正性能.后續(xù)研究工作中,我們將考慮校正更多的噪聲成分,在現(xiàn)有模型基礎(chǔ)上繼續(xù)完善,不斷提高校正算法的實用性.