基于Chambolle-Pock算法框架的高階TV圖像重建算法

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

（山西大學計算機與信息技術學院，太原 030006）

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

0 引言

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

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

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

為壓制TV 算法引起的塊狀偽影，解決ASD-POCS 調節參數問題，本文提出了一種基于CP 框架的HOTV 算法，并詳細地推導了該算法，進一步揭示了算法的性能。本文以二階梯度構建二階TV 范數，并采用CP 算法對模型進行求解。為進一步比較其性能，采用基于波浪背景的Shepp-Logan 模型、灰度漸變模體以及真實CT圖像模體共三種模體作為實驗模體，分別在理想投影數據與含噪投影數據下詳細地對CP-HOTV算法與CP-TV算法開展了定性和定量的評估。

1 相關知識

1.1 成像模型

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

其中：g是大小為I的列向量，表示投影數據；u是大小為J的列向量，表示圖像像素；M是大小為I*J的矩陣，表示系統矩陣。Mi，j即為第j個像素對第i個投影測量的貢獻。本文中用Siddon射線驅動法求取系統矩陣［16］。

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

1.2 TV最優化模型

根據上述成像模型，對于離散圖像u約束的TV 最優化模型可表示為：

其受下述不等式約束：

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

式（5）表示圖像每個像素的梯度。

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

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

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

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

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

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

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

其中，σ為CP 算法的非負參數。上述詳細推導過程可參考Sidky等發表的相關研究［7］。

2 本文方法

2.1 HOTV最優化模型設計

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

現定義圖像的HOTV 范數是圖像二階梯度大小變換的l1范數：

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

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

則D2變換為：

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

則圖像的HOTV范數可表示為：

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

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

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

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

研究表明，總變差最小化在凸變分成像方法中起著重要作用，其優勢在于允許在求解過程中出現明顯的不連續，這對于很多成像問題至關重要，因為邊緣往往代表著重要的特征。而本文設計的HOTV 重建模型有更好的保邊性能。病態的反向成像問題一般可以分為兩類：凸問題和非凸問題。而凸問題相較于非凸問題的優點在于其可以計算出全局最優，在很多情況下，具有良好的精度和合理的時間，且獨立于初始化，因此凸問題解決方案的質量將完全取決于模型的準確性。本文中利用約束的CP 算法對HOTV 重建模型進行求解。CP 算法是一種原對偶的優化算法，該算法具有收斂性保證，提供了非啟發式的收斂檢查-對偶間隙。則根據上述理論，本文設計CP-HOTV算法如算法1所示。

算法1 HOTV-CP算法。

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

數詳細說明請參考文獻［11］。

2.3 圖像質量的評價標準

本文中將以均方根誤差Root Mean Square Error，RMSE）和結構相似性（Structural SIMilarity，SSIM）作為重建質量評估判據，計算式如下：

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

3 實驗與結果分析

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

3.1 理想數據稀疏重建

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

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

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

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

表1中，當迭代結束后，會發現對于三種模體，HOTV重建圖像的RMSE 值均低于TV 重建圖像的RMSE 值；當迭代結束，對于三種模型，相較于TV 重建圖像，HOTV 重建圖像的SSIM 值更接近于1，表明HOTV 重建圖像的圖像結構與真值圖像更相似。

通過定性定量比較可知，對于基于波動背景的Shepp-Logan 模體以及灰度漸變模體這一類含有灰度漸變過渡區域的模體，HOTV 算法能提高其重建精度；而對于復雜結構的醫學圖像，相較于傳統的TV 算法，HOTV 算法能有效減少塊狀偽影的同時取得更高的重建精度。

3.2 含噪數據稀疏重建

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

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

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

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

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

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

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

3.3 CP-HOTV算法的收斂性

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

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

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

3.4 算法運行時間的比較

在實際的醫學CT 和工業CT 中，算法的運行時間也是衡量算法的重要指標。為了更全面地對CP-TV 算法與CPHOTV算法進行比較，表5中以三種模體在理想數據投影下的算法運行時間為例。

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

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

4 結語

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