總變差約束的數(shù)據(jù)分離最小圖像重建模型及其Chambolle-Pock求解算法?

喬志偉

(山西大學計算機與信息技術學院,太原 030006)

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

1 引 言

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

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

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

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)點:其一,算法參數(shù)均可以計算得到而不需人為設定；其二,可以解所有凸最優(yōu)化模型,無論模型平滑還是非平滑；其三,算法的每個子最優(yōu)化問題一般均有閉合解.CP算法框架可以用來推導求解TVcDM模型的CP算法實例.

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

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

2 方 法

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

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

2.1 成像模型

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

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

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

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

2.2 最優(yōu)化模型

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

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

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

這里D是一個大小為2N ×N的矩陣

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

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

D2變換可以表示為

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

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

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

這兩個參數(shù)雖然出現(xiàn)在最優(yōu)化模型中,但是這里稱其為算法參數(shù),因為它們不能決定最優(yōu)化模型的解,但是可以影響解的收斂路徑和速度.

2.3 Chambolle-Pock算法及其解釋

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

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

在表2中,x表示一個長度為N的向量；a是?1球的半徑；m是一個長度為N的向量,其每個元素對應x中相應元素的絕對值；sign(x)是符號函數(shù):正數(shù)的函數(shù)值是1,0的函數(shù)值是0,負數(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算法的推導步驟

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

2.4.1 CP算法框架及其推導方法

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

在此模型中,x和y是分別屬于向量空間X和Y的多維向量；K是一個矩陣,表示一個線性變換,定義了一個從X到Y的線性映射；G和F分別是關于x和y的多元函數(shù),分別定義了從X和Y到實數(shù)的一個映射.CP框架要求這兩個函數(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表示轉置.

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

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

prox最鄰近映射的表達式為

該最鄰近映射本質上是一個最優(yōu)化問題,但是因為z是一個變量,所以其結果是一個關于z的函數(shù).

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

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

2)求取F?；

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

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

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

2.4.2 TVcDM-CP算法實例的推導

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

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

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

這樣,

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

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

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

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

2.5 重建參數(shù)

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

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

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

3 結 果

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

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

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

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

3.1.1 仿真模體及重建條件

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

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

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

3.1.2 重建結果的定性及定量分析

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

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

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

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

兩種模體的重建圖像的RMSE均達到了10?4,一方面表示算法的收斂條件達到了,另一方面表明逆犯罪成功,因為此反問題被高精度地求解了出來.

圖1 Shepp-Logan仿真模體的TVcDM重建結果 (a)Shepp-Logan仿真模體圖像；(b)重建圖像；(c)模型圖像和重建圖像左右中心線波形的比較:紅色線是真實波形,藍色線是重建波形,二者幾乎完全重合；圖像的顯示窗口是[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重建結果 (a)FORBILD仿真模體圖像；(b)重建圖像；(c)模體圖像和重建圖像左右中心線波形的比較圖:紅色線是真實波形,藍色線是重建波形,二者幾乎完全重合；圖像的顯示窗口是[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個度量標準,通過觀察3個度量標準的迭代走勢,揭示CP算法的迭代規(guī)律.這些度量標準如下式所示:

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

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

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

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

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

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

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

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

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

圖5 Shepp-Logan模體重建中TVcDM與FBP算法的重建結果比較 圖像上面的數(shù)字表示投影個數(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的重建結果與FBP算法重建結果的比較.由圖可以看出,隨著投影個數(shù)從120個減少到30個,FBP重建圖像的條狀偽影越來越嚴重,60個角度下的FBP重建引入了明顯的偽影而30個角度下的偽影已經異常嚴重；同時看到,各個角度下的TVcDM重建圖像幾乎完全相同,看不出明顯的因為投影個數(shù)減少而產生的質量退化現(xiàn)象.30個角度下的投影已經可以實現(xiàn)TVcDM的高精度重建,表明TVcDM模型有高精度稀疏重建能力.

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

圖6 FORBILD模體重建中TVcDM與FBP算法的重建結果比較 圖像上面的數(shù)字表示投影個數(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隨投影個數(shù)變化的走勢圖 (a)對應Shepp-Logan模體重建；(b)對應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個角度對應的重建RMSE比圖7(a)中30個角度對應的重建RMSE大.此現(xiàn)象表明:不同的重建對象,因為復雜程度不同,所需的最少投影個數(shù)是不同的.FORBILD圖像比Shepp-Logan圖像復雜,所以要達到特定的重建精度,需要的投影個數(shù)更多.

3.3 TV限對重建精度的影響

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

圖8和圖9是Shepp-Logan模體重建結果.圖8是不同TV限約束下的對含噪投影的重建圖像.圖9是不同TV限約束下的對含噪投影的重建圖像與仿真模體(真值)的差分絕對值圖像.從圖8(c)和圖9(c)可以看出,當TV限選取正確,即選取TV真值時,可以獲得高精度的重建結果,TVcDM模型不但壓制了因為稀疏采集(只有30個投影)可能引起的條狀偽影,而且壓制了因為高斯白噪聲污染而產生的圖像噪聲.從圖8和圖9的(a)和(b)可以看出,當TV限選取過小時,TV范數(shù)的平滑作用過大,使得重建圖像產生了低頻過平滑偽影.而從圖8和圖9的(d)和(e)可以看出,當TV限選取過大時,TV范數(shù)的平滑作用減小,使得重建圖像產生了高頻噪聲.可見,TV限的選擇對重建有重要影響,其決定了最優(yōu)化模型的收斂解.只有在平滑和細節(jié)保持之間,尋找到平衡點,選取出最適合的TV限,才可以高精度地重建圖像.

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

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

圖8 TV限的不同選擇對重建精度的影響 (a)—(e)分別是TV限為TV真值的0.6,0.8,1.0,1.2和1.4倍時對應的重建圖像；為了突出顯示,采用了偽彩色顯示方式；顯示窗口是[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倍時對應的重建圖像與仿真模體的差分絕對值圖像；為了突出顯示,采用了偽彩色顯示方式；顯示窗口是[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].

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

圖10 RMSE隨TV限變化的走勢圖 (a)和(b)分別對應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倍時對應的重建圖像；為了突出顯示,采用了偽彩色顯示方式；顯示窗口是[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倍時對應的重建圖像與仿真模體的差分絕對值圖像；為了突出顯示,采用了偽彩色顯示方式；顯示窗口是[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ù)變?yōu)?0.采用Shepp-Logan和FORBILD兩個模體分別進行實驗.為了分析λ和ν的不同選擇對重建的影響,設計5種不同的組合形式,使算法均迭代至1000次,觀察RMSE的下降態(tài)勢,比較其對重建速率的影響.

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

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

圖13 不同的λ和ν的組合下RMSE在1000次迭代范圍內的收斂速率比較 (a)和(b)分別對應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 模型設計相關問題

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

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

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

4.2 關于重建算法

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

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

5 結 論

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

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

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

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