霍其潤 李建武 陸 耀 秦 明

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

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

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

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

本文的主要貢獻(xiàn)如下:

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

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

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

1 變分模型設(shè)計(jì)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中,

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

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

2 模型求解及其校正算法

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

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

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

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

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

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

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

故有

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

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

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

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

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

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

3)repeat

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

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

4)條狀偽影K=I ?S

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

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

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

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

3.1 模擬數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

3.2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

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

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

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

對(duì)于圖像中的平坦區(qū)域,其標(biāo)準(zhǔn)差很大程度上與圖像中的偽影相關(guān),標(biāo)準(zhǔn)差的降低相應(yīng)地反映了圖像偽影的減少[16].顯然對(duì)于處理后的圖像,其平坦區(qū)域的LSNR 值越大,表明圖像偽影的去除效果越好.對(duì)應(yīng)圖9 各方框標(biāo)識(shí)的區(qū)域位置,我們依次計(jì)算出原始圖像及其各算法處理結(jié)果中相應(yīng)區(qū)域的LSNR 值,如表2 所列.可以看出相對(duì)于原始圖像,各種算法處理后相應(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 實(shí)現(xiàn)分析

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

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

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

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

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

3.3.2 RTV 與TV 的約束效果對(duì)比

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

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

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

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

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

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

3.3.3 RTV 與TV 的約束效果對(duì)比

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

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

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

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

3.3.4 算法收斂性

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

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

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

4 結(jié)論

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