喬志偉

(山西大學(xué)計(jì)算機(jī)與信息技術(shù)學(xué)院,太原 030006)

(2018年4月27日收到；2018年7月2日收到修改稿)

1 引 言

在計(jì)算機(jī)斷層成像(computed tomography,CT)圖像重建中,以濾波反投影(filtered backprojection,FBP)法為代表的解析重建方法在商用CT中占據(jù)主流[1].其缺點(diǎn)是需要完備的投影數(shù)據(jù).如稀疏視角下的FBP重建算法會(huì)引起條狀偽影[2].隨著壓縮感知等稀疏優(yōu)化理論的提出,基于優(yōu)化的迭代算法成為了十余年來的研究熱點(diǎn)[3,4].其中,總變差最小(total variation minimization,TV)算法是最成功的優(yōu)化算法之一,其已經(jīng)在平板CT[5]、偏置掃描平板CT[6]、短掃描平板CT[7],C-arm CT[8]、低劑量平板CT[9]及機(jī)載成像儀(on-board imager,OBI)平板CT等[10]三維CT中得到了成功的應(yīng)用,展示了其高精度的稀疏重建能力.

TV算法可以高精度地稀疏重建的原因可以從多個(gè)方面來解釋.TV算法體現(xiàn)了壓縮感知的思想.壓縮感知認(rèn)為,如果物體在某個(gè)變換域是稀疏的,那么可以僅僅使用稀疏采集到的原始數(shù)據(jù)來高精度地重建物體.其實(shí)現(xiàn)策略是在數(shù)據(jù)保真的約束下,使得物體的稀疏變換的?1范數(shù)最小[11].TV算法使用的稀疏變換是梯度大小變換(gradient magnitude transform,GMT)[12].TV算法體現(xiàn)了先驗(yàn)約束的思想.在稀疏重建時(shí),迭代法構(gòu)建的線性方程組是欠定的,使得該線性方程組有無窮多個(gè)解.TV算法可以被認(rèn)為加入了一種先驗(yàn)知識,即物體的GMT是稀疏的.這樣,TV算法是從無窮多個(gè)解中選取TV最小的解.TV算法也體現(xiàn)了低通濾波的作用.TV算法最初被當(dāng)作去噪方法應(yīng)用到圖像處理中,其可以在去噪的同時(shí)保持圖像的邊緣信息[13].FBP的稀疏重建引入的條狀偽影是一種高頻噪聲.TV正是通過對這種噪聲的低通濾波,消除了條狀偽影,實(shí)現(xiàn)了高精度稀疏重建.

TV最優(yōu)化模型有約束形式和非約束形式兩種類型.非約束形式的模型參數(shù)是一個(gè)平衡因子,用來平衡數(shù)據(jù)保真力度和正則力度,其沒有明確的物理意義,難以選擇.約束形式的TV模型,其模型參數(shù)有物理意義,更容易選擇[9].常用的約束的TV模型,以數(shù)據(jù)保真項(xiàng)為約束,TV正則項(xiàng)為目標(biāo)函數(shù),其缺點(diǎn)是數(shù)據(jù)保真項(xiàng)對應(yīng)的數(shù)據(jù)容差限這一參數(shù)對重建質(zhì)量敏感,選取方法魯棒性較差[12].芝加哥大學(xué)Pan研究組從2014年陸續(xù)提出了TV約束的、數(shù)據(jù)分離最小(TV constrained,data divergence minimization,TVcDM)新型TV模型,并將之應(yīng)用到了扇形束ROI重建[14],C-arm CT[8]、短掃描平板CT[7]及正電子發(fā)射斷層掃描(positron emission computed tomography,PET)[15]中.該模型的模型參數(shù)是TV限,其在較大范圍的變化對圖像質(zhì)量影響較小,故易于選取,選取方法魯棒性較好.

TVcDM模型是一種約束的最優(yōu)化模型,傳統(tǒng)的非約束優(yōu)化模型求解算法,如交替方向乘子法(alternating direction method of multipliers,ADMM)和分裂布萊格曼法(split Bregman,SB)[16],不適宜用來求解該模型.Chambolle-Pock(CP)算法框架[17?20]是一種優(yōu)秀的最優(yōu)化算法框架,不但適用于非約束的優(yōu)化模型,而且適用于約束的優(yōu)化模型.它有諸多優(yōu)點(diǎn):其一,算法參數(shù)均可以計(jì)算得到而不需人為設(shè)定；其二,可以解所有凸最優(yōu)化模型,無論模型平滑還是非平滑；其三,算法的每個(gè)子最優(yōu)化問題一般均有閉合解.CP算法框架可以用來推導(dǎo)求解TVcDM模型的CP算法實(shí)例.

Pan研究組[20]于2012年推導(dǎo)了一系列CP算法實(shí)例.這些實(shí)例分別用來求解最小二乘模型、非約束的?2-TV模型、非約束的?1-TV模型、非約束的KL-TV模型以及數(shù)據(jù)分離約束的、TV最小(data divergence constrained,TV minimization,DDcTV)模型,但是當(dāng)時(shí)尚未提出TVcDM模型,所以并沒有給出TVcDM模型的CP算法實(shí)例的推導(dǎo).2014—2017年,該組的一系列使用TVcDM模型的CT和PET圖像重建的文章,主要聚焦在完全重建、有限角度或者投影被截?cái)嗲闆r下的實(shí)際重建問題.2018年,本文作者與Pan研究組及芝加哥大學(xué)EPRI中心研究人員將這一模型應(yīng)用于電子順磁共振成像(electron paramagnetic resonance imaging,EPRI)的稀疏重建中[21].這些研究聚焦在模型和算法在實(shí)際數(shù)據(jù)中的評估,而沒有使用仿真數(shù)據(jù),即存在金標(biāo)準(zhǔn)的情況下,從定量的角度展開收斂規(guī)律、稀疏重建能力評估和參數(shù)選擇等理論和方法層面的研究.

鑒于如上原因,本文在總結(jié)出CP算法實(shí)例推導(dǎo)方法的基礎(chǔ)上,詳細(xì)推導(dǎo)了TVcDM-CP算法實(shí)例,并給出算法的詳細(xì)解釋；以Shepp-Logan和FORBILD仿真模體[22]的CT重建為例,實(shí)現(xiàn)該算法,驗(yàn)證模型和算法的正確性；通過分析多個(gè)度量標(biāo)準(zhǔn)的迭代行為,揭示其收斂規(guī)律；通過稀疏重建,評估該模型的高精度稀疏重建能力；通過對含噪投影數(shù)據(jù)的重建,研究了模型參數(shù)對重建質(zhì)量的影響.最后研究了算法參數(shù)對收斂速率的影響.

2 方 法

連續(xù)-連續(xù)(continuous-to-continuous,CC)成像模型把投影和圖像均看成連續(xù)函數(shù).而離散-離散(discrete-to-discrete,DD)成像模型把投影和圖像均看成離散函數(shù).這樣,DD模型就是一個(gè)線性方程組,基于DD模型的重建就是求方程組的解.一般情況下,該線性方程組是一個(gè)大規(guī)模的、病態(tài)的且往往欠定的(如稀疏重建時(shí))方程組,使得求解成為需要重點(diǎn)攻克的問題.為了得到所需要的有用解,往往需要將此成像模型進(jìn)一步建模為一個(gè)最優(yōu)化問題,使其可以將稀疏先驗(yàn)等先驗(yàn)知識引入模型.其后,需要設(shè)計(jì)求解算法,以解最優(yōu)化模型.

基于如上脈絡(luò),以CT重建為背景,以TVcDM模型為研究對象,在該部分將按照成像系統(tǒng)建模-最優(yōu)化模型-CP算法-參數(shù)選擇這一研究鏈展開敘述.因?yàn)楸疚闹荚谔剿餍路f的TV模型及其求解算法,不失一般性,以平行束CT作為研究對象.別的成像模態(tài)或者CT成像的其他掃描模式遵循同樣的規(guī)律,其研究方法只需要做相應(yīng)的類推或擴(kuò)展.

2.1 成像模型

平行束CT的DD成像模型是一個(gè)線性方程組,如(1)式所示:

這里,g是一個(gè)大小為M的列向量,表示離散的投影數(shù)據(jù)；u是一個(gè)大小為N的列向量,表示離散圖像；A是一個(gè)大小為M×N的矩陣,這里代表2D Radon變換；Am,n表示第n個(gè)像素對第m個(gè)投影測量值的貢獻(xiàn).對2D Radon變換來說,Am,n是第n個(gè)像素(正方形)與第m條射線(直線)的相交長度.

設(shè)圖像的大小為nx×ny,將其一列一列串行化,就可以得到圖像的向量表示,其大小為N=nx×ny.設(shè)投影集大小為np×na,即共有na個(gè)角度下的投影,每個(gè)投影上有np個(gè)點(diǎn),將其一列一列串行化,就可以得到投影數(shù)據(jù)的向量表示,其大小為M=np×na.

實(shí)際上,求取Am,n的方法就是所謂的投影方法,因?yàn)橥队笆窍到y(tǒng)矩陣與圖像向量的乘積.投影方法包括像素驅(qū)動(dòng)法[23],Siddon射線驅(qū)動(dòng)法[24],Joseph射線驅(qū)動(dòng)法[25]和距離驅(qū)動(dòng)法[26]等.傳統(tǒng)的像素驅(qū)動(dòng)法會(huì)引入高頻震蕩偽影,已經(jīng)設(shè)計(jì)了一種精確的像素驅(qū)動(dòng)法[23]；本文使用了精確的像素驅(qū)動(dòng)法來生成投影,即求取系統(tǒng)矩陣.

2.2 最優(yōu)化模型

由于(1)式所示的線性方程組的病態(tài)、大規(guī)模及欠定等因素,需要進(jìn)一步將之建模為一個(gè)最優(yōu)化模型.根據(jù)重建的需要,可以將壓縮感知、低秩矩陣等理論及各種先驗(yàn)知識耦合,以設(shè)計(jì)最優(yōu)化模型.本文研究一種新穎的基于壓縮感知的TV模型——TVcDM.

TVcDM,即TV約束的、數(shù)據(jù)分離最小模型,可以表示為

這里,u?是最優(yōu)化模型的解,即被重建的圖像；∥.∥表示?2范數(shù)的平方；圖像的TV范數(shù)∥u∥TV是圖像梯度大小變換|Du|mag的?1范數(shù):

這里D是一個(gè)大小為2N ×N的矩陣

D1和D2均為N×N的矩陣,分別表示沿著x和y軸方向的離散梯度變換.

令ux,y,x∈[1,nx],y∈[1,ny]表示二維圖像的每個(gè)像素,則D1變換可以表示為

D2變換可以表示為

令|.|mag表示對一個(gè)二維向量求模,可以是?1范數(shù)模,也可以是?2范數(shù)模.在本文中,采用?2范數(shù)模,即

這樣,(3)式的各向同性形式可以表示為

為了提高收斂速度,引入λ和ν兩個(gè)算法參數(shù)以調(diào)整TV范數(shù)和數(shù)據(jù)保真項(xiàng)這兩個(gè)凸函數(shù)的相對大小,如此得到的新的TVcDM模型為[7,8,14]

這兩個(gè)參數(shù)雖然出現(xiàn)在最優(yōu)化模型中,但是這里稱其為算法參數(shù),因?yàn)樗鼈儾荒軟Q定最優(yōu)化模型的解,但是可以影響解的收斂路徑和速度.

2.3 Chambolle-Pock算法及其解釋

表1給出了用于求解最優(yōu)化模型(9)的CP算法偽代碼.

在表1中,∥.∥sv是矩陣的最大奇異值,其求法見文獻(xiàn)[20]中的算法8；向量u,uˉ及s的大小為N；向量p的大小為M；向量q和a的大小為2N；6.3中的向量1I是一個(gè)大小為N 的“1”向量；6.2中的投影操作符ProjectOnto?1Ballνt1的求法見表2的偽代碼；uconv是達(dá)到設(shè)定的實(shí)用收斂條件后的收斂解.

在表2中,x表示一個(gè)長度為N的向量；a是?1球的半徑；m是一個(gè)長度為N的向量,其每個(gè)元素對應(yīng)x中相應(yīng)元素的絕對值；sign(x)是符號函數(shù):正數(shù)的函數(shù)值是1,0的函數(shù)值是0,負(fù)數(shù)的函數(shù)值是?1.

表1 用于求解最優(yōu)化模型(9)的CP算法偽代碼[7,8,14]Table 1.Pseudocode of the CP algorithm for solving the optimization model shown in(9)[7,8,14].

表2 ProjectOnto?1Ball的求解算法偽代碼[14]Table 2. Pseudocode of the solving algorithm for ProjectOnto?1Ball[14].

2.4 CP算法的推導(dǎo)步驟

本節(jié)首先從CP算法框架總結(jié)出算法實(shí)例的推導(dǎo)方法,然后推導(dǎo)TVcDM模型的CP算法實(shí)例.

2.4.1 CP算法框架及其推導(dǎo)方法

CP算法框架可以解形如(10)式所示的最優(yōu)化模型.

在此模型中,x和y是分別屬于向量空間X和Y的多維向量；K是一個(gè)矩陣,表示一個(gè)線性變換,定義了一個(gè)從X到Y(jié)的線性映射；G和F分別是關(guān)于x和y的多元函數(shù),分別定義了從X和Y到實(shí)數(shù)的一個(gè)映射.CP框架要求這兩個(gè)函數(shù)是凸函數(shù),但不要求平滑.

CP算法框架見表3.

表3 CP算法框架(N步迭代)[17]Table 3.The CP algorithm framework(N steps of iterations)[17].

表3中,∥K∥SV是指矩陣K的最大奇異值；F?是指凸函數(shù)F 的共軛凸函數(shù)；proxσ[F?]和proxτ[G]是最鄰近映射操作；T表示轉(zhuǎn)置.

任意一個(gè)凸函數(shù)H(z)的共軛凸函數(shù)定義為

這里,需要注意,max操作符不是求取使得目標(biāo)函數(shù)最大的自變量,而是求取當(dāng)c為最大化值時(shí)目標(biāo)函數(shù)的值.

prox最鄰近映射的表達(dá)式為

該最鄰近映射本質(zhì)上是一個(gè)最優(yōu)化問題,但是因?yàn)閦是一個(gè)變量,所以其結(jié)果是一個(gè)關(guān)于z的函數(shù).

現(xiàn)在,可以總結(jié)出使用CP算法框架推導(dǎo)CP算法實(shí)例的步驟為:

1)構(gòu)造如(10)式所示的最優(yōu)化問題；

2)求取F?；

3)求 取proxσ[F?]和proxτ[G]兩 個(gè) 最 鄰 近映射；

4)代入表3的第4行和第5行形成算法實(shí)例.

在如上步驟中,關(guān)鍵問題是基于(11)式的凸共軛函數(shù)的求取和基于(12)式的最鄰近映射的求取.CP算法的優(yōu)勢之一是每個(gè)子最優(yōu)化問題都有閉合解.這就要求在(12)式中,函數(shù)H(Z)要足夠簡單,以使得可以求得閉合形式的解.現(xiàn)有研究表明,在圖像重建中使用的大部分凸函數(shù)所對應(yīng)的最鄰近映射均有閉合解[20].

2.4.2 TVcDM-CP算法實(shí)例的推導(dǎo)

(9)式所對應(yīng)的最優(yōu)化問題,可以通過如下變換形成(10)式的形式.

其中,δ?1Ball是一個(gè)指示函數(shù),表示如果一個(gè)n維向量在特定半徑的?1球內(nèi)部,則函數(shù)值為0,否則,值為∞,其表達(dá)式為:該式表示,如果向量x在半徑為a的?1球內(nèi)部,則函數(shù)值為0,否則值為∞.

現(xiàn)在,根據(jù)如下所示的一系列變換,將(12)與(10)式對應(yīng)起來.

這樣,

根據(jù)CP算法實(shí)例的推導(dǎo)步驟,現(xiàn)在的任務(wù)是求取兩個(gè)凸函數(shù)F1(p)和F2(q)的共軛凸函數(shù),以及求取FF和G三個(gè)函數(shù)的prox最鄰近映射.

根據(jù)(11)式,可得

根據(jù)(12)式,可得

將(15),(19)—(21)代入表3,可以得到表1所示的TVcDM模型的CP算法實(shí)例.

2.5 重建參數(shù)

最優(yōu)化模型及其算法構(gòu)成了一套完整的重建方案,其中所有的參數(shù)構(gòu)成了重建參數(shù).與模型相關(guān)的參數(shù)稱為模型參數(shù)；與算法相關(guān)的參數(shù)稱為算法參數(shù).

模型參數(shù)包括TV限t1、像素大小、投影方法等.在本文中,像素大小和探元大小都?xì)w一化為1,投影方法采用精確的像素驅(qū)動(dòng)法,故將在3.3節(jié)重點(diǎn)研究TV限對重建的影響.

算法參數(shù)包括σ,τ,λ,ν和θ,因?yàn)棣?τ和θ參數(shù)均可以計(jì)算得到,而不需要人為設(shè)定,所以不需要關(guān)注.λ和ν對算法的收斂速度和路徑有影響,故將在3.4節(jié)研究其對收斂速率的影響.

3 結(jié) 果

在本節(jié)中,設(shè)計(jì)4個(gè)研究點(diǎn):1)重建方法的驗(yàn)證及CP算法的收斂行為分析；2)重建方法的稀疏重建能力評估；3)TV限對重建精度的影響；4)λ和ν的選擇對收斂速率的影響.

3.1 逆犯罪研究——重建方法的驗(yàn)證及CP算法的收斂行為

從數(shù)學(xué)的角度講,圖像重建是一個(gè)反問題.對于平行束CT重建,解析法,如濾波反投影算法,是反Radon變換問題；迭代法是線性方程組的求解問題.逆犯罪是一種評估反問題是否正確的驗(yàn)證方法.當(dāng)仿真數(shù)據(jù)理想且充分的情況下,如果求得的解的誤差在計(jì)算機(jī)精度范圍內(nèi)可以充分小,則認(rèn)為逆犯罪成功,表明設(shè)計(jì)的整套重建方案及其實(shí)現(xiàn)——從成像系統(tǒng)建模、最優(yōu)化問題建模、算法設(shè)計(jì)到計(jì)算機(jī)實(shí)現(xiàn),都是正確的.

在本部分中,采用Shepp-Logan和FORBILD模型為仿真模型；仿真生成理想而充分的投影數(shù)據(jù)；使用第二部分的重建方法進(jìn)行重建；通過逆犯罪研究驗(yàn)證整套重建方法的正確性.

3.1.1 仿真模體及重建條件

分別采用Shepp-Logan模體和FORBILD模體為仿真模型,圖像大小為256×256,每個(gè)像素被看成單位1大小的像素.旋轉(zhuǎn)中心定位于128×128.探測器長度為256,每個(gè)探元大小為單位1.均勻采集[0,π]范圍360個(gè)角度下的投影,所以投影數(shù)據(jù)的大小為256×360.采用精確的像素驅(qū)動(dòng)法定義系統(tǒng)矩陣A,并生成投影數(shù)據(jù).如此構(gòu)建的線性方程組,有256×256個(gè)未知數(shù),有256×360個(gè)方程,屬于過定方程.但是因?yàn)榉匠探M是相容的,所以方程組有惟一的精確解.

在CP算法中,設(shè)定λ=1,b=1即ν=νA,其余算法參數(shù)按照表1中計(jì)算得到.TV限設(shè)定為仿真模型的TV值.

CP迭代算法的收斂條件定義為RMSE(un,utruth)6 10?4.RMSE的定義見(22)和(23)式.

3.1.2 重建結(jié)果的定性及定量分析

重建結(jié)果采用定性觀察和定量分析的方法.一方面視覺檢查圖像質(zhì)量,一方面比較圖像中線位置的波形,以進(jìn)行定性觀察.定量分析,則采用均方根誤差(root-mean-square-error,RMSE)作為度量標(biāo)準(zhǔn),如下式所示:

式中x是待評估向量,s是向量真值,n是向量的長度.

對Shepp-Logan重建來說,CP算法在迭代到2273次時(shí),達(dá)到了設(shè)定的收斂條件.重建結(jié)果見圖1. 圖1(a)是Shepp-Logan仿真模型圖像；圖1(b)是收斂態(tài)的重建圖像.比較這兩個(gè)圖像可以發(fā)現(xiàn)它們幾乎完全一樣,用肉眼難以分辨彼此.圖1(c)是重建圖像和仿真模型圖像的中心線波形的比較,從圖中可以看出,兩條波形幾乎完全重合.圖像和波形的主觀評價(jià)表明取得了高精度重建.

對FORBILD重建來說,CP算法在迭代到1712次時(shí),達(dá)到了設(shè)定的收斂條件.重建結(jié)果見圖2.圖2(a)是FORBILD仿真模型圖像；圖2(b)是收斂態(tài)的重建圖像.比較這兩個(gè)圖像可以發(fā)現(xiàn),它們幾乎完全一樣,用肉眼難以分辨彼此；圖2(c)是重建圖像和仿真模型圖像的中心線波形比較圖.從圖中可以看出,兩條波形幾乎完全重合.圖像和波形的主觀評價(jià)表明取得了高精度重建.

兩種模體的重建圖像的RMSE均達(dá)到了10?4,一方面表示算法的收斂條件達(dá)到了,另一方面表明逆犯罪成功,因?yàn)榇朔磫栴}被高精度地求解了出來.

圖1 Shepp-Logan仿真模體的TVcDM重建結(jié)果 (a)Shepp-Logan仿真模體圖像；(b)重建圖像；(c)模型圖像和重建圖像左右中心線波形的比較:紅色線是真實(shí)波形,藍(lán)色線是重建波形,二者幾乎完全重合；圖像的顯示窗口是[0,1]Fig.1.The TVcDM reconstruction images of the Shepp-Logan simulation phantom:(a)The Shepp-Logan phantom image；(b)the reconstructed image；(c)comparison of the profiles on the vertical center-line of the phantom image and the reconstructed image.In(c),it may be seen that the red truth profile and the blue reconstructed profile are almost completely overlapped.

圖2 FORBILD仿真模體的TVcDM重建結(jié)果 (a)FORBILD仿真模體圖像；(b)重建圖像；(c)模體圖像和重建圖像左右中心線波形的比較圖:紅色線是真實(shí)波形,藍(lán)色線是重建波形,二者幾乎完全重合；圖像的顯示窗口是[0,1]Fig.2.The TVcDM reconstruction images of the fORBILD simulation phantom:(a)the fORBILD phantom image；(b)the reconstructed image；(c)comparison of the profiles on the vertical center-line of the phantom image and the reconstructed image.In(c),it may be seen that the red truth profile and the blue reconstructed profile are almost completely overlapped.

3.1.3 CP算法的收斂行為

為了描述收斂行為,引入3個(gè)度量標(biāo)準(zhǔn),通過觀察3個(gè)度量標(biāo)準(zhǔn)的迭代走勢,揭示CP算法的迭代規(guī)律.這些度量標(biāo)準(zhǔn)如下式所示:

(23)式描述的是重建的圖像與仿真模型圖像的RMSE距離；(24)式描述的是重建圖像生成的投影與原始投影之間的?2范數(shù)距離；(25)式描述重建圖像的TV值與TV真值的相對誤差.在此逆犯罪研究中,因?yàn)槟Ｐ偷耐队昂蛨D像是完全一致的、數(shù)據(jù)是充分且精確的,所以這些度量標(biāo)準(zhǔn)在收斂態(tài)均可以充分小,表示取得了高精度重建.

圖3和圖4分別展示了Shepp-Logan和FORBILD模體的CP算法的收斂行為,圖3(a)—(c)分別展示了(23)—(25)式所示的度量標(biāo)準(zhǔn)的迭代走勢.

從圖3和圖4的(a)可以看出,圖像誤差不斷下降,達(dá)到了設(shè)定的收斂條件,而且可以看到其有繼續(xù)下降的趨勢.如果繼續(xù)推進(jìn)迭代,圖像誤差甚至可以達(dá)到10?6.灰度圖像的灰度只有256個(gè)等級,即從0到255,理論上說,當(dāng)圖像誤差達(dá)到1/256≈3.9‰時(shí),顯示器已經(jīng)無法分辨存在此誤差水平的兩幅圖像.所以,我們設(shè)定逆犯罪標(biāo)志,即收斂條件為圖像誤差小于等于10?4.

圖3和圖4的(b)顯示的是數(shù)據(jù)誤差,即數(shù)據(jù)殘差的走勢.數(shù)據(jù)殘差是TVcDM模型的目標(biāo)函數(shù).可以看出,數(shù)據(jù)殘差不斷下降,而且在當(dāng)前的停止點(diǎn)有繼續(xù)向下的趨勢.

圖3 TVcDM-CP算法重建Shepp-Logan模體時(shí)的收斂行為 (a)重建誤差M1(n)的迭代走勢；(b)數(shù)據(jù)誤差M2(n)的迭代走勢；(c)重建圖像的TV相對誤差M3(n)的迭代走勢.Fig.3.The convergence behavior of the TVcDM-CP algorithm for reconstructing the Shepp-Logan phantom:(a)The iteration trend of the reconstruction error M1(n)；(b)the iteration trend of the data error M2(n)；(c)the iteration trend of the relative TV error M3(n).

圖4 TVcDM-CP算法重建FORBILD模體時(shí)的收斂行為 (a)重建誤差M1(n)的迭代走勢；(b)數(shù)據(jù)誤差M2(n)的迭代走勢；(c)重建圖像的TV相對誤差M3(n)的迭代走勢Fig.4.The convergence behavior of the TVcDM-CP algorithm for reconstructing the fORBILD phantom:(a)The iteration trend of the reconstruction error M1(n)；(b)the iteration trend of the data error M2(n)；(c)the iteration trend of the relative TV error M3(n).

圖3和圖4的(c)顯示了TV值相對于TV限的相對誤差,其整體趨勢向下,但是在兩個(gè)重建案例中,當(dāng)?shù)胶笃?出現(xiàn)了較大抖動(dòng),相對誤差變大.這種現(xiàn)象是該算法迭代行為的一個(gè)特點(diǎn),在我們各種各樣的重建實(shí)踐中經(jīng)常出現(xiàn),這并不影響算法的收斂性.雖然TV值相對誤差出現(xiàn)了向上抖動(dòng),但是圖像誤差和數(shù)據(jù)誤差是整體不斷下降的,表明最優(yōu)化模型的解仍然在向收斂態(tài)逼近.其實(shí),可以觀察到TV值下降到了一個(gè)很小的值,比如10?7之后產(chǎn)生了向上的抖動(dòng),其達(dá)到一定值之后,又會(huì)繼續(xù)向下.此后,即使再有抖動(dòng),其上限也會(huì)保持在一個(gè)很小的值附近,比如10?4.而TV值是指圖像梯度變換的?1范數(shù),TV值在此極小范圍的波動(dòng),分配到圖像每個(gè)像素上(在本節(jié)案例中,圖像包含256×256=65536個(gè)像素),圖像基本沒有變化.

從圖3和圖4的所有子圖可以看出,所有的迭代曲線均存在抖動(dòng)現(xiàn)象,這是TVcDM-CP算法迭代行為的一個(gè)特點(diǎn).圖像誤差和數(shù)據(jù)殘差(目標(biāo)函數(shù))整體會(huì)不斷下降,直到計(jì)算機(jī)浮點(diǎn)數(shù)精度(本文因?yàn)橥队芭c反投影操作使用GPU加速,所以均采用了單精度數(shù))所允許的收斂態(tài).

3.2 重建方法的稀疏重建能力評估

為了評估TVcDM模型的稀疏重建能力,在仿真實(shí)驗(yàn)中,分別從30,60,90和120個(gè)角度重建圖像,分析投影個(gè)數(shù)對重建的影響.重建條件與3.1節(jié)相同,除了Shepp-Logan模體重建的收斂條件變?yōu)榈?000次結(jié)束,而FORBILD模體重建的收斂條件為迭代至2000次結(jié)束.

圖5是Shepp-Logan模體的TVcDM的重建結(jié)果與FBP算法重建結(jié)果的比較.由圖可以看出隨著投影個(gè)數(shù)從120個(gè)減少到30個(gè),FBP重建的條狀偽影越來越嚴(yán)重,30個(gè)角度下的FBP重建引入了明顯的偽影；同時(shí)看到,各個(gè)角度下的TVcDM重建圖像幾乎完全相同,看不出明顯的因?yàn)橥队皞€(gè)數(shù)減少而產(chǎn)生的質(zhì)量退化現(xiàn)象.30個(gè)角度下的投影已經(jīng)可以實(shí)現(xiàn)TVcDM的高精度重建,表明TVcDM模型有高精度稀疏重建能力.

圖5 Shepp-Logan模體重建中TVcDM與FBP算法的重建結(jié)果比較 圖像上面的數(shù)字表示投影個(gè)數(shù),左邊的文字表示所使用的算法；顯示窗口為[0,1]Fig.5.Comparison of the reconstructed Shepp-Logan images by use of TVcDM and FBP algorithm.The numbers above the image indicate projection number and the text at the left of the image indicate the algorithm used.The display window is[0,1].

圖6是FORBILD模體的TVcDM的重建結(jié)果與FBP算法重建結(jié)果的比較.由圖可以看出,隨著投影個(gè)數(shù)從120個(gè)減少到30個(gè),FBP重建圖像的條狀偽影越來越嚴(yán)重,60個(gè)角度下的FBP重建引入了明顯的偽影而30個(gè)角度下的偽影已經(jīng)異常嚴(yán)重；同時(shí)看到,各個(gè)角度下的TVcDM重建圖像幾乎完全相同,看不出明顯的因?yàn)橥队皞€(gè)數(shù)減少而產(chǎn)生的質(zhì)量退化現(xiàn)象.30個(gè)角度下的投影已經(jīng)可以實(shí)現(xiàn)TVcDM的高精度重建,表明TVcDM模型有高精度稀疏重建能力.

為了定量比較TVcDM與FBP算法在不同投影個(gè)數(shù)情形下的重建精度,圖7繪制了重建圖像的RMSE相對投影個(gè)數(shù)的走勢圖.圖7(a)對應(yīng)Shepp-Logan模體重建；圖7(b)對應(yīng)FORBILD模體重建.兩個(gè)模體重建對應(yīng)的RMSE走勢規(guī)律是相同的.隨著投影個(gè)數(shù)從120減小到30,FBP重建結(jié)果的RMSE逐漸增大,表示FBP算法對投影個(gè)數(shù)的變化是敏感的,其沒有高精度稀疏重建能力,當(dāng)投影個(gè)數(shù)為30個(gè)時(shí),其RMSE過大,圖像精度過低,導(dǎo)致重建圖像失去有效性.而TVcDM算法對投影個(gè)數(shù)變化不敏感,投影個(gè)數(shù)為60,90和120時(shí),RMSE非常接近.投影個(gè)數(shù)為30時(shí),因?yàn)橥队耙暯沁^于稀疏,重建誤差增大.但即使TVcDM算法只使用30個(gè)角度下的投影,其重建圖像的精度也遠(yuǎn)遠(yuǎn)高于FBP算法從120個(gè)角度下的投影重建圖像的精度.

圖6 FORBILD模體重建中TVcDM與FBP算法的重建結(jié)果比較 圖像上面的數(shù)字表示投影個(gè)數(shù),左邊的文字表示所使用的算法；顯示窗口為[0,1]Fig.6.Comparison of the reconstructed FORBILD images by use of TVcDM and FBP algorithm.The numbers above the image indicate projection number and the text at the left of the image indicate the algorithm used.The display window is[0,1].

圖7 重建圖像RMSE隨投影個(gè)數(shù)變化的走勢圖 (a)對應(yīng)Shepp-Logan模體重建；(b)對應(yīng)FORBILD模體重建Fig.7.Plot of RMSE of the reconstructed image as function of projection number with(a)corresponding to Shepp-Logan phantom reconstruction and(b)FORBILD phantom reconstruction.

圖5和圖6的定性主觀視覺檢查和圖7的定量誤差分析均表明:TVcDM-CP算法有從稀疏投影中高精度重建圖像的能力.圖7(b)中30個(gè)角度對應(yīng)的重建RMSE比圖7(a)中30個(gè)角度對應(yīng)的重建RMSE大.此現(xiàn)象表明:不同的重建對象,因?yàn)閺?fù)雜程度不同,所需的最少投影個(gè)數(shù)是不同的.FORBILD圖像比Shepp-Logan圖像復(fù)雜,所以要達(dá)到特定的重建精度,需要的投影個(gè)數(shù)更多.

3.3 TV限對重建精度的影響

在仿真實(shí)驗(yàn)中,給投影加45.0 dB的高斯白噪聲,進(jìn)行重建.投影個(gè)數(shù)為30個(gè)；TV限分別設(shè)定為TV真值的0.6,0.8,1.0,1.2和1.4倍；其余重建條件相同.

圖8和圖9是Shepp-Logan模體重建結(jié)果.圖8是不同TV限約束下的對含噪投影的重建圖像.圖9是不同TV限約束下的對含噪投影的重建圖像與仿真模體(真值)的差分絕對值圖像.從圖8(c)和圖9(c)可以看出,當(dāng)TV限選取正確,即選取TV真值時(shí),可以獲得高精度的重建結(jié)果,TVcDM模型不但壓制了因?yàn)橄∈璨杉?只有30個(gè)投影)可能引起的條狀偽影,而且壓制了因?yàn)楦咚拱自肼曃廴径a(chǎn)生的圖像噪聲.從圖8和圖9的(a)和(b)可以看出,當(dāng)TV限選取過小時(shí),TV范數(shù)的平滑作用過大,使得重建圖像產(chǎn)生了低頻過平滑偽影.而從圖8和圖9的(d)和(e)可以看出,當(dāng)TV限選取過大時(shí),TV范數(shù)的平滑作用減小,使得重建圖像產(chǎn)生了高頻噪聲.可見,TV限的選擇對重建有重要影響,其決定了最優(yōu)化模型的收斂解.只有在平滑和細(xì)節(jié)保持之間,尋找到平衡點(diǎn),選取出最適合的TV限,才可以高精度地重建圖像.

圖10(a)是Shepp-Logan模體重建中RMSE隨TV限變化的走勢圖.可以看到TV限過大或者過小,都導(dǎo)致了更大的誤差,選取合適的TV限會(huì)取得較小誤差,實(shí)現(xiàn)高精度重建.

圖11和圖12是FORBILD模體重建結(jié)果.結(jié)合圖10(b),可以看出,過小或者過大的TV限均導(dǎo)致了更大的誤差,合適的TV限可以取得較小誤差,取得高精度重建.

圖8 TV限的不同選擇對重建精度的影響 (a)—(e)分別是TV限為TV真值的0.6,0.8,1.0,1.2和1.4倍時(shí)對應(yīng)的重建圖像；為了突出顯示,采用了偽彩色顯示方式；顯示窗口是[0,1]Fig.8.The impact of the selection of TV tolerance on reconstruction accuracy:(a)–(e)are the reconstructed images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,1].

圖9 TV限的不同選擇對重建精度的影響 (a)—(e)分別是TV限為TV真值的0.6,0.8,1.0,1.2和1.4倍時(shí)對應(yīng)的重建圖像與仿真模體的差分絕對值圖像；為了突出顯示,采用了偽彩色顯示方式；顯示窗口是[0,0.2]Fig.9.The impact of the selection of TV tolerance on reconstruction accuracy:(a)–(e)are the reconstructed absolute-difference images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,0.2].

兩種模體復(fù)雜程度不同,圖像特征不同,但是在此研究中,TV限對重建精度的影響規(guī)律一致.實(shí)驗(yàn)結(jié)果和定性分析一致表明,TV限對重建精度有重要影響,只有選擇最佳TV限才可以在圖像的平滑度和細(xì)節(jié)保持之間找到最佳平衡點(diǎn),實(shí)現(xiàn)高精度重建.

圖10 RMSE隨TV限變化的走勢圖 (a)和(b)分別對應(yīng)Shepp-Logan重建和FORBILD重建Fig.10.Plot of RMSE as function of TV tolerance with(a)corresponding to Shepp-Logan reconstruction and(b)FORBILD reconstruction.

圖11 TV限的不同選擇對重建精度的影響 (a)—(e)分別是TV限為TV真值的0.6,0.8,1.0,1.2和1.4倍時(shí)對應(yīng)的重建圖像；為了突出顯示,采用了偽彩色顯示方式；顯示窗口是[0,1]Fig.11.The impact of the selection of TV tolerance on reconstruction accuracy:(a)–(e)are the reconstructed images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,1].

圖12 TV限的不同選擇對重建精度的影響 (a)—(e)分別是TV限為TV真值的0.6,0.8,1.0,1.2和1.4倍時(shí)對應(yīng)的重建圖像與仿真模體的差分絕對值圖像；為了突出顯示,采用了偽彩色顯示方式；顯示窗口是[0,0.2]Fig.12.The impact of the selection of TV tolerance on reconstruction accuracy.(a)–(e)are the reconstructed absolute-difference images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,0.2].

3.4 λ和ν的選擇對收斂速率的影響

仍然使用3.1.1中設(shè)定的重建條件,只是將投影個(gè)數(shù)變?yōu)?0.采用Shepp-Logan和FORBILD兩個(gè)模體分別進(jìn)行實(shí)驗(yàn).為了分析λ和ν的不同選擇對重建的影響,設(shè)計(jì)5種不同的組合形式,使算法均迭代至1000次,觀察RMSE的下降態(tài)勢,比較其對重建速率的影響.

圖13是收斂速率比較圖,其中圖13(a)對應(yīng)Shepp-Logan重建；圖13(b)對應(yīng)FORBILD重建.

從圖13可以看出,不同的λ和ν對應(yīng)不同的收斂速率.對于兩個(gè)模體重建來說,λ=1,ν=0.1νA均為收斂速率最快的一個(gè)組合.但是如果按照從慢到快排序,對兩個(gè)模體重建而言,順序不完全相同.λ和ν的最佳選擇是依賴于特定重建對象的.即使對于同一個(gè)重建對象,當(dāng)投影個(gè)數(shù)發(fā)生變化時(shí),λ和ν的最快組合也可能發(fā)生變化.對于一個(gè)特定的重建任務(wù),要想找到最快的組合,需要使用不同的組合重建,通過觀察迭代收斂行為,選擇最快的一個(gè)組合.

圖13 不同的λ和ν的組合下RMSE在1000次迭代范圍內(nèi)的收斂速率比較 (a)和(b)分別對應(yīng)Shepp-Logan重建和FORBILD重建Fig.13. Comparison of the convergence rates of RMSE of reconstructions with different λ and ν in the iteration range of 1 to 1000 with(a)and(b)corresponding to Shepp-Logan and FORBILD reconstruction,respectively.

4 討 論

4.1 模型設(shè)計(jì)相關(guān)問題

為了提高TV模型的精確重建能力,人們提出了各種TV范數(shù)的變種:廣義的各向異性TV范數(shù),保邊TV范數(shù),適用于有限角度重建的各向異性TV范數(shù),高階TV范數(shù)；TpV(p∈[0,1])范數(shù)等.類似的最優(yōu)化模型還包括小波變換最小、基于字典學(xué)習(xí)的稀疏模型以及非局部均值等.所有這些模型均可以用一個(gè)統(tǒng)一優(yōu)化模型來描述.在統(tǒng)一優(yōu)化模型具體化時(shí),其一要考慮選取何種正則模型,其二要考慮選取何種數(shù)據(jù)保真模型,其三則是考慮選擇約束形式還是非約束形式.

最優(yōu)化模型是數(shù)學(xué)模型,它是對物理模型的一種近似,因此各種各樣的最優(yōu)化模型沒有優(yōu)劣之分,只有在特定重建場合下的適用性問題.對于圖像重建的一個(gè)特定任務(wù),應(yīng)用驅(qū)動(dòng)地設(shè)計(jì)最適合的最優(yōu)化模型是基于優(yōu)化的重建方案設(shè)計(jì)的最重要的工作.

最優(yōu)化模型決定了最終的解,所以模型參數(shù)對重建有重要影響.當(dāng)前,盡管人們已經(jīng)系統(tǒng)研究了模型參數(shù)對重建的影響,但是模型參數(shù)的自動(dòng)選擇問題應(yīng)該是今后的研究重點(diǎn).

4.2 關(guān)于重建算法

本文所用模型的求解是非常具有挑戰(zhàn)性的:1)規(guī)模巨大,系統(tǒng)矩陣很大,設(shè)計(jì)的算法中不能出現(xiàn)矩陣求逆操作；2)往往欠定,方程個(gè)數(shù)遠(yuǎn)小于未知量個(gè)數(shù)；3)病態(tài),因?yàn)樵肼暫蛿?shù)據(jù)的不一致性,造成方程組是不相容的；4)TV范數(shù)是一個(gè)不平滑的凸函數(shù).

人們已經(jīng)為各種各樣的最優(yōu)化模型設(shè)計(jì)了求解算法,如ASD-POCS算法[12],ADMM算法以及本文使用的CP算法等.需要注意的是,模型決定解,求解算法僅僅決定解的收斂路徑和速度.所以,各種各樣的求解算法的重建精度是沒有可比性的,惟一可以比較的是其收斂速度.本文采用CP算法的原因是該算法保證收斂、算法參數(shù)不需要調(diào)節(jié)且每個(gè)子最優(yōu)化問題均有閉合解.

5 結(jié) 論

本文設(shè)計(jì)了TVcDM-CP重建方法.將圖像重建問題建模為一個(gè)TVcDM最優(yōu)化模型.在該約束的最優(yōu)化模型中,TV正則項(xiàng)是約束項(xiàng),數(shù)據(jù)保真項(xiàng)是目標(biāo)函數(shù).基于CP算法框架,詳細(xì)推導(dǎo)了TVcDM模型的CP算法實(shí)例；經(jīng)仿真模型重建,驗(yàn)證了模型及算法的正確性并分析了算法的收斂行為；評估了模型的稀疏重建能力；分析了模型參數(shù)的選擇對重建的影響.

TVcDM模型相較于DDcTV模型,一個(gè)可能的優(yōu)勢是模型參數(shù)比較容易選擇.DDcTV模型的模型參數(shù)是數(shù)據(jù)容差限[12],其大小決定了重建圖像的投影和實(shí)際的測量投影的距離.在實(shí)際重建中,投影數(shù)據(jù)包含噪聲和不一致性,只有選擇合適的數(shù)據(jù)容差限才可以有效抑制噪聲和不一致性.選取過小,則噪聲被引入；選取過大則重建圖像過于平滑,細(xì)節(jié)丟失.TVcDM的模型參數(shù)是TV限,合理的TV限估計(jì)會(huì)實(shí)現(xiàn)高精度重建,而選取過小則圖像過于平滑,選取過大則圖像被引入噪聲.TV限在小范圍內(nèi)的變化對重建影響不大(參考文獻(xiàn)[8]的圖7),所以該參數(shù)的選擇比較魯棒.

本文研究表明,TVcDM-CP重建方法可以高精度地從稀疏角度下的投影中重建圖像；在CP算法的收斂過程中,會(huì)出現(xiàn)抖動(dòng)現(xiàn)象,但是整體呈現(xiàn)下降趨勢；TV限對重建質(zhì)量有重要影響,合理選取該參數(shù)可以取得高精度的重建質(zhì)量；λ和ν的不同取值組合會(huì)導(dǎo)致不同的收斂速率.

本文所設(shè)計(jì)的TVcDM-CP算法是一種通用的圖像重建模型,可以被應(yīng)用到平行束CT、扇形束CT、錐形束CT中,也可以被應(yīng)用到PET,MRI及EPRI中.依據(jù)本文的研究鏈,即成像系統(tǒng)建模-最優(yōu)化問題建模-CP算法設(shè)計(jì)-算法實(shí)現(xiàn)及驗(yàn)證-實(shí)際重建,各種重建模態(tài)均可以實(shí)現(xiàn)高精度稀疏重建.在此研究鏈中,需要做的改變僅僅是將本文的成像模型中的系統(tǒng)矩陣變?yōu)橄鄳?yīng)的成像模態(tài)對應(yīng)的系統(tǒng)矩陣,然后任務(wù)驅(qū)動(dòng)地選擇模型參數(shù)和算法參數(shù).