席雅睿，喬志偉，溫 靜，張艷嬌，楊雯晶，閆慧文

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

（?通信作者電子郵箱zqiao@sxu.edu.cn）

0 引言

計(jì)算機(jī)斷層成像（Computer Tomography，CT）作為最為先進(jìn)的無損檢測(cè)手段，應(yīng)用于醫(yī)學(xué)界和工業(yè)界。在醫(yī)學(xué)領(lǐng)域，其要求在保持重建圖像精度的同時(shí)盡量減少對(duì)于病患的輻射劑量。而目前，基于濾波反投影（Filtered BackProjection，F(xiàn)BP）為經(jīng)典的解析法依然廣泛應(yīng)用于醫(yī)學(xué)和商用CT 中。在傳統(tǒng)奈奎斯特采樣定理的要求下，應(yīng)用解析法稀疏重建高精度的圖像是十分困難的事情，會(huì)產(chǎn)生條狀偽影［1-2］。隨著壓縮感知（Compressed Sensing，CS）理論的提出，全變分（Total Variation，TV）最小算法成為稀疏重建的研究熱點(diǎn)［3］。Sidky等［4］提出了TV 最小重建模型，在仿真狀態(tài)下實(shí)現(xiàn)了高精度稀疏重建。之后，Sidky 等［5］完善了該算法，詳細(xì)給出了ASDPOCS（Adaptive Steepest Descent-Projection Onto Convex Sets）算法的偽代碼，并將其應(yīng)用到三維錐束CT 中。隨后，有大量的TV 變體以及求解的優(yōu)化算法被提出，如：Chambolle 等［6］提出了原始對(duì)偶優(yōu)化算法。Sidky等［7］利用CP（Chambolle-Pock）優(yōu)化算法去求解重建模型，并推導(dǎo)出一系列的算法實(shí)例。Tang 等［8］設(shè)計(jì)了一種基于CP 框架的對(duì)偶鄰近點(diǎn)算法應(yīng)用于全變分圖像重建。Sidky 等［9］又提出了一種約束最小的TpV（Totalp-Variation）模型，并使用CP 算法對(duì)模型進(jìn)行求解。宋文琪等［10］對(duì)于基于CP 算法的TV 重建算法進(jìn)行研究，并且與傳統(tǒng)的FBP 算法進(jìn)行了詳細(xì)的比較。喬志偉［11］設(shè)計(jì)了一種TV 約束的數(shù)據(jù)分離最小TV 模型，推導(dǎo)了CP 算法的實(shí)例，并對(duì)該模型在重建精度、算法收斂性等方面進(jìn)行詳細(xì)研究。這些TV 變體的提出以及優(yōu)化算法的改進(jìn)極大地促進(jìn)了稀疏重建的發(fā)展。

然而，在實(shí)際的重建過程中，人們發(fā)現(xiàn)TV 算法有時(shí)會(huì)引入塊狀偽影，尤其在圖像去噪領(lǐng)域中，此現(xiàn)象尤其明顯，這引起了很多學(xué)者的注意并力求壓制此效應(yīng)，如：Chan 等［12］提出了基于高階TV（High-Order TV，HOTV）的圖像復(fù)原模型；Lysaker 等［13］提出了應(yīng)用高階TV 去除圖像噪聲，該模型有效去除了塊狀偽影，同時(shí)保持了圖像平坦區(qū)域的光滑度；Papafitsoros 等［14］針對(duì)圖像的修復(fù)問題詳細(xì)比較了傳統(tǒng)TV 與HOTV 的區(qū)別，并采用Split Bregman 框架進(jìn)行求解；胡悅等［15］提出了基于增廣拉格朗日乘子的快速高階TV模型，該模型能夠高效去噪并且一定程度上提高了算法的速度。HOTV 算法已經(jīng)在圖像處理領(lǐng)域中得到大量的使用，并且其不僅可以減小梯度效應(yīng)、提高圖像的質(zhì)量而且對(duì)于圖像去噪和圖像復(fù)原都有很好的效果；但是HOTV 算法并沒有在圖像重建領(lǐng)域中展開深刻的研究。

ASD-POCS 算法是一種有效約束的TV 最小圖像重建模型的求解算法，它使用ART（Algebra Reconstruction Technology）算法讓數(shù)據(jù)殘差逐漸下降，并使用自適應(yīng)梯度下降算法讓圖像的TV值逐漸下降，最終找到數(shù)據(jù)保真約束下的TV 范數(shù)最小的最優(yōu)解；但是ASD-POCS 算法往往需要憑借經(jīng)驗(yàn)去調(diào)節(jié)參數(shù)，并且其ART 部分由于需要逐根射線處理，無法并行改造，需要耗費(fèi)大量時(shí)間。然而，CP 算法其本身是嚴(yán)格按照數(shù)學(xué)公式推導(dǎo)而來，不需要根據(jù)經(jīng)驗(yàn)去調(diào)節(jié)參數(shù)。

為壓制TV 算法引起的塊狀偽影，解決ASD-POCS 調(diào)節(jié)參數(shù)問題，本文提出了一種基于CP 框架的HOTV 算法，并詳細(xì)地推導(dǎo)了該算法，進(jìn)一步揭示了算法的性能。本文以二階梯度構(gòu)建二階TV 范數(shù)，并采用CP 算法對(duì)模型進(jìn)行求解。為進(jìn)一步比較其性能，采用基于波浪背景的Shepp-Logan 模型、灰度漸變模體以及真實(shí)CT圖像模體共三種模體作為實(shí)驗(yàn)?zāi)ｓw，分別在理想投影數(shù)據(jù)與含噪投影數(shù)據(jù)下詳細(xì)地對(duì)CP-HOTV算法與CP-TV算法開展了定性和定量的評(píng)估。

1 相關(guān)知識(shí)

1.1 成像模型

本文以二維平行束CT重建為研究對(duì)象。已知離散-離散（Discrete to Discrete，D2D）的成像模型是將投影和圖像均看成離散函數(shù)，則D2D 模型本身就是一個(gè)線性方程組。那么D2D成像模型可以被表示為：

其中：g是大小為I的列向量，表示投影數(shù)據(jù)；u是大小為J的列向量，表示圖像像素；M是大小為I*J的矩陣，表示系統(tǒng)矩陣。Mi，j即為第j個(gè)像素對(duì)第i個(gè)投影測(cè)量的貢獻(xiàn)。本文中用Siddon射線驅(qū)動(dòng)法求取系統(tǒng)矩陣［16］。

D2D成像模型本身是一個(gè)線性方程組。在實(shí)際情況中這個(gè)方程組往往是巨大的、欠定的。如圖像自身大小為256×256，投影數(shù)據(jù)也為256×256，則由式（1）可知其系統(tǒng)矩陣大小為4 GB。這是一個(gè)巨大的且直接求逆不可行的線性方程組。而在稀疏重建中，成像模型就是一個(gè)有無數(shù)解的欠定線性系統(tǒng)。成像系統(tǒng)的噪聲更會(huì)將此線性系統(tǒng)變成一個(gè)病態(tài)的線性系統(tǒng)，鑒于上述原因，為獲得更高的重建精度，需進(jìn)一步將此模型建模為一個(gè)最優(yōu)化問題。

1.2 TV最優(yōu)化模型

根據(jù)上述成像模型，對(duì)于離散圖像u約束的TV 最優(yōu)化模型可表示為：

其受下述不等式約束：

若設(shè)圖像大小為n*n，s∈[1，n]，t∈[1，n]，us，t表示第s行、第t列的像素，則可定義為：

式（5）表示圖像每個(gè)像素的梯度。

1.3 約束的CP算法求解TV模型

CP 算法是一種基本-對(duì)偶的形式，其最小化和對(duì)偶最大化可定義為：

式中：x、y分別為線性空間X、Y中向量；K為X空間到Y(jié)空間的線性變換；G和F為非平滑凸函數(shù)；“?”表示凸共軛。則可得基本-對(duì)偶問題的一般鞍點(diǎn)最優(yōu)化問題可表示為：

則式（3）～（4）的TV最小化模型可表示為：

則利用CP 算法求解TV 模型的整個(gè)目標(biāo)函數(shù)可表示為F(Kx)形式，有：

TV最小化問題的F1函數(shù)、F2函數(shù)可表示為：

那么，式（3）最接近映射proxσ為：

其中，σ為CP 算法的非負(fù)參數(shù)。上述詳細(xì)推導(dǎo)過程可參考Sidky等發(fā)表的相關(guān)研究［7］。

2 本文方法

2.1 HOTV最優(yōu)化模型設(shè)計(jì)

為了進(jìn)一步求解1.1 節(jié)中D2D 巨大的線性方程組，結(jié)合HOTV 算法在圖像處理領(lǐng)域的良好效果，本文提出了一種約束的HOTV最小模型，可表示為：

現(xiàn)定義圖像的HOTV 范數(shù)是圖像二階梯度大小變換的l1范數(shù)：

其中：D1和D2均為N*N的矩陣，分別為沿著x和y軸方向的離散梯度變換。

設(shè)圖像大小為n*n，s∈[1，n]，t∈[1，n]，us，t表示第s行、第t列的像素，則D1變換可表示為：

則D2變換為：

那么圖像的離散梯度變換可表示為：

則圖像的HOTV范數(shù)可表示為：

以128×128 的Shepp-Logan 圖像為例，其二階梯度大小圖像如圖1（c）所示。

圖1 Shepp-Logan模體及其梯度大小圖像Fig.1 Shepp-Logan phantom and its gradient magnitude images

比較圖1（b）和圖1（c）可看出，一階和二階梯度大小圖像均可表示邊緣信息，但是兩者提取的邊緣特征有較大的區(qū)別，這就決定了TV算法和HOTV算法會(huì)有不同的稀疏重建性能。

2.2 HOTV重建模型的Chambolle-Pock求解算法

研究表明，總變差最小化在凸變分成像方法中起著重要作用，其優(yōu)勢(shì)在于允許在求解過程中出現(xiàn)明顯的不連續(xù)，這對(duì)于很多成像問題至關(guān)重要，因?yàn)檫吘壨碇匾奶卣鳌６疚脑O(shè)計(jì)的HOTV 重建模型有更好的保邊性能。病態(tài)的反向成像問題一般可以分為兩類：凸問題和非凸問題。而凸問題相較于非凸問題的優(yōu)點(diǎn)在于其可以計(jì)算出全局最優(yōu)，在很多情況下，具有良好的精度和合理的時(shí)間，且獨(dú)立于初始化，因此凸問題解決方案的質(zhì)量將完全取決于模型的準(zhǔn)確性。本文中利用約束的CP 算法對(duì)HOTV 重建模型進(jìn)行求解。CP 算法是一種原對(duì)偶的優(yōu)化算法，該算法具有收斂性保證，提供了非啟發(fā)式的收斂檢查-對(duì)偶間隙。則根據(jù)上述理論，本文設(shè)計(jì)CP-HOTV算法如算法1所示。

算法1 HOTV-CP算法。

其中：L為系統(tǒng)矩陣M和2.1節(jié)所求D的二范數(shù)；τ和σ為算法的非負(fù)參數(shù)；θ為參數(shù)，θ∈[0，1]；n為迭代次數(shù)；u0，p0，q0和初始值均為0，8）～11）行對(duì)其進(jìn)行更新。算法中所有參

數(shù)詳細(xì)說明請(qǐng)參考文獻(xiàn)［11］。

2.3 圖像質(zhì)量的評(píng)價(jià)標(biāo)準(zhǔn)

本文中將以均方根誤差Root Mean Square Error，RMSE）和結(jié)構(gòu)相似性（Structural SIMilarity，SSIM）作為重建質(zhì)量評(píng)估判據(jù)，計(jì)算式如下：

其中：u為重建圖像；r為真值圖像；μr和μu分別為r、u的平均值；和分別為r、u的方差；c1和c2是為了防止分母為零的常數(shù)。

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

為了進(jìn)一步驗(yàn)證HOTV 算法的性能，本文采用了基于波動(dòng)背景的Shepp-Logan 模體、灰度漸變模體以及真實(shí)CT 圖像模體這三種模體，分別在理想數(shù)據(jù)投影和含噪數(shù)據(jù)投影兩種情形下進(jìn)行稀疏重建。

3.1 理想數(shù)據(jù)稀疏重建

基于波動(dòng)背景的Shepp-Logan 模型是在傳統(tǒng)的Shepp-Logan 模體上加入波動(dòng)背景；灰度漸變模體的外部是線性漸變，內(nèi)部是由小正方形中心到其四條邊中心點(diǎn)的線性漸變；而真實(shí)CT 圖像模體是將現(xiàn)實(shí)中的CT 圖像作為模體。三種模體的大小均為128×128。本文實(shí)驗(yàn)室采用的都是平行束的掃描方式，其中成像的幾何參數(shù)為：在360°的采樣角度范圍內(nèi)，均勻地采集了30個(gè)角度，每個(gè)角度的投影信號(hào)為128個(gè)測(cè)量值；探測(cè)器探源大小歸一化為1；圖像大小也歸一化為1。使用TV 算法和HOTV 算法進(jìn)行稀疏重建，迭代次數(shù)設(shè)為5 000 次，其重建結(jié)果如圖2 所示，表1 為迭代結(jié)束后重建結(jié)果的RMSE值和SSIM值。

圖2 理想數(shù)據(jù)投影下三種模體重建結(jié)果Fig.2 Reconstruction results of three phantoms under ideal data projection

圖2 給出了傳統(tǒng)的TV 算法和HOTV 算法的重建結(jié)果。由圖2 可以看出，對(duì)于基于波動(dòng)背景的Shepp-Logan 模體，其HOTV 重建結(jié)果與TV 的重建結(jié)果視覺質(zhì)量相當(dāng)；但是對(duì)于灰度漸變模體，可以明顯發(fā)現(xiàn)，TV 重建圖像中有很多條形偽影，而HOTV 重建圖像很平滑；對(duì)于真實(shí)CT 圖像模體，由于模體本身復(fù)雜度遠(yuǎn)高于其他兩組，其圖像質(zhì)量均不算太高。比較其矩形框部分，HOTV 重建圖像結(jié)構(gòu)清晰且保留了更多細(xì)節(jié)信息，而TV重建圖像中有多處塊狀偽影。

表1 理想數(shù)據(jù)投影下迭代結(jié)束后不同模體的RMSE值和SSIM值Tab.1 RMSE and SSIM values of different phantoms after iteration under ideal data projection

表1中，當(dāng)?shù)Y(jié)束后，會(huì)發(fā)現(xiàn)對(duì)于三種模體，HOTV重建圖像的RMSE 值均低于TV 重建圖像的RMSE 值；當(dāng)?shù)Y(jié)束，對(duì)于三種模型，相較于TV 重建圖像，HOTV 重建圖像的SSIM 值更接近于1，表明HOTV 重建圖像的圖像結(jié)構(gòu)與真值圖像更相似。

通過定性定量比較可知，對(duì)于基于波動(dòng)背景的Shepp-Logan 模體以及灰度漸變模體這一類含有灰度漸變過渡區(qū)域的模體，HOTV 算法能提高其重建精度；而對(duì)于復(fù)雜結(jié)構(gòu)的醫(yī)學(xué)圖像，相較于傳統(tǒng)的TV 算法，HOTV 算法能有效減少塊狀偽影的同時(shí)取得更高的重建精度。

3.2 含噪數(shù)據(jù)稀疏重建

為進(jìn)一步評(píng)估在特定噪聲強(qiáng)度情形下HOTV 算法的重建性能，本文將對(duì)投影數(shù)據(jù)分別加入方差為0.01、0.02、0.03、0.04 和0.05 的高斯噪聲，通過比較TV 算法和HOTV 算法的重建性能，分析二者的抗噪能力。

為了更深刻地比較兩種算法的抗噪性，本文還將對(duì)投影數(shù)據(jù)分別加入λ值為1、3、5、7 和9 的泊松噪聲，定性、定量地對(duì)HOTV 算法和TV 算法的重建質(zhì)量進(jìn)行比較。為簡(jiǎn)便起見，此部分只給出灰度漸變模體的重建結(jié)果。圖像重建的成像條件及參數(shù)設(shè)定與2.1 節(jié)部分相同，方差為0.01 的高斯噪聲和λ值為1 的泊松噪聲情形下的重建圖像如圖3 所示，表2 為迭代結(jié)束后重建圖像的RMSE值和SSIM值。

圖3 含噪數(shù)據(jù)投影下灰度漸變模體重建結(jié)果Fig.3 Reconstruction results of grayscale gradual changing phantom under noisy data projection

從圖3 和表2 可知，相較于TV 算法，HOTV 算法的重建精度更高。加入其他強(qiáng)度高斯噪聲情形下的圖像重建效果與上述實(shí)驗(yàn)效果類似，故不再展示相應(yīng)重建結(jié)果。表3 給出了各種強(qiáng)度高斯噪聲情形下的HOTV及TV算法重建圖像的RMSE值以及SSIM 值。表4 給出了各種強(qiáng)度泊松噪聲情形下的HOTV 及TV 算法重建圖像的RMSE 值以及SSIM 值。無論是高斯噪聲還是泊松噪聲，均可以看出，HOTV 算法可以更高精度地稀疏重建圖像，這表明HOTV 算法具有更強(qiáng)的抗噪能力。

表2 含噪數(shù)據(jù)投影下重建圖像的RMSE值和SSIM值Tab.2 RMSE and SSIM values of reconstructed images under noisy data projection

表3 不同強(qiáng)度高斯噪聲情形下兩種重建算法的重建精度比較Tab3 Reconstruction accuracy comparison between two algorithms under different intensities of Gaussian noise

表4 不同強(qiáng)度泊松噪聲情形下兩種重建算法的重建精度比較Tab.4 Reconstruction accuracy comparison between two algorithms under different intensities of Poisson noise

3.3 CP-HOTV算法的收斂性

為了描述CP-HOTV 算法的收斂性，通過觀察重建過程中RMSE 值的迭代走勢(shì)，揭示CP-HOTV 算法的迭代規(guī)律。以基于波動(dòng)背景的Shepp-Logan 模型為例，圖4 為CP-TV 和CPHOTV重建圖像的迭代走勢(shì)。

圖4 重建圖像的RMSE隨迭代次數(shù)的變化Fig.4 RMSE of reconstructed image varying with number of iterations

從圖4 可以看出：CP-TV 算法與CP-HOTV 算法的RMSE迭代過程中均在不斷下降；迭代過程中，很明顯CP-HOTV 算法的RMSE 值幾乎一直小于CP-TV 算法的RMSE 值；迭代結(jié)束時(shí)，CP-HOTV 算法的RMSE 值為0.008 7，而CP-TV 算法的RMSE 值為0.009 0。綜上可以得出，CP-HOTV 算法在整個(gè)重建迭代過程中收斂精度更高。

3.4 算法運(yùn)行時(shí)間的比較

在實(shí)際的醫(yī)學(xué)CT 和工業(yè)CT 中，算法的運(yùn)行時(shí)間也是衡量算法的重要指標(biāo)。為了更全面地對(duì)CP-TV 算法與CPHOTV算法進(jìn)行比較，表5中以三種模體在理想數(shù)據(jù)投影下的算法運(yùn)行時(shí)間為例。

表5 不同算法針對(duì)不同模體理在想數(shù)據(jù)投影下的運(yùn)行時(shí)間對(duì)比單位：sTab.5 Running time comparison of different algorithm for different phantoms under ideal data projectionunit：s

由表5 可以得出，在理想數(shù)據(jù)投影下，CP-TV 算法和CPHOTV算法在5 000次迭代結(jié)束后，CP-HOTV算法的運(yùn)行時(shí)間略長(zhǎng)，且隨著模型難度的加深，其算法的運(yùn)行時(shí)間也越來越長(zhǎng)。相較于CP-TV 算法，CP-HOTV 算法中運(yùn)用的是HOTV 范數(shù)，其求解過程較為復(fù)雜，但卻有更高的精確度。

4 結(jié)語

本文提出了一種數(shù)據(jù)保真約束的HOTV 最小圖像重建模型，進(jìn)而設(shè)計(jì)了相應(yīng)的CP求解算法。該算法在壓制塊狀偽影的同時(shí)在一定程度上提高了重建精度。復(fù)雜醫(yī)學(xué)圖像中往往存在很多灰度線性漸變或波浪式漸變的區(qū)域，而HOTV 算法是從數(shù)據(jù)保真項(xiàng)所定義的凸集中選擇HOTV 最小的圖像，所以其更適合重建灰度波動(dòng)特征明顯的圖像。本文研究的目的是探索HOTV 算法的性能及可能的重建優(yōu)勢(shì)，故設(shè)計(jì)的三種模體均含有灰度波動(dòng)特征。本文采用的CP算法，雖然解決了ASD-POCS 算法憑借經(jīng)驗(yàn)調(diào)節(jié)參數(shù)的問題，但是占用大量?jī)?nèi)存。若使用圖形處理（Graphics Processing Unit，GPU）加速算法，將會(huì)使原有算法占用更小的內(nèi)存且有更快的運(yùn)行速度，這將是迭代法日后發(fā)展的趨勢(shì)，也是我們進(jìn)一步的研究方向。