稀疏角度CT圖像重建算法研究

林澤田 王 單

1(北京航空航天大學數學與系統科學學院 北京 100191)2(北京航空航天大學經濟管理學院 北京 100191)

0 引 言

計算機斷層成像CT(Computed Tomography)技術是指利用X射線穿透物體的衰減信息進行重建來獲得物體的斷層圖像信息的技術，該技術被公認為20世紀后期最偉大的科技成果之一[1]。在醫學診斷領域，CT技術已在醫學成像領域確立了其不可動搖的地位。不僅如此，在工業無損檢測、樹木年輪測定、地球資源勘探等領域，CT技術也成為了有力的手段。然而，在實際的應用過程中，由于受輻射劑量、對比度等的約束，抑或是受現實條件限制，只能采集到部分數據。而在不完全數據的CT圖像重建中，大多涉及到稀疏角度圖像重建?；诖吮尘?，本文展開研究。

2006年，文獻[2]提出了壓縮感知理論。該理論基于已知信號具有稀疏性或可壓縮性，對信號數據進行采集、編碼以及解碼。當信號具有稀疏性或可壓縮性時，該理論表明通過采集少量的信號投影值便可以實現信號的準確性或者近似重建。該理論突破了傳統采樣定理的局限，實現了壓縮和采樣的統一，在信息與信號處理方面帶來革命性的突破，并已成功應用在人臉識別、地震勘探、通信與信號處理等諸多領域。

信號重構算法是一種由投影矩陣以及投影測量值重建出來原始信號，一般把問題轉化為求解l0范數的優化問題。目前為止，典型的重建算法有貪婪算法[3]、l1最小化算法以及全變差TV(Total Variation)最小化算法[4]。自從2006年，文獻[2]提出壓縮感知理論，發現l1最小化與l0最小化的等價之后，在如安檢、核磁共振等方面，壓縮感知和l1最小化算法已經被廣泛應用。進而，對于圖像重建問題而言，文獻[5]研究表明應用全變差最小化算法能夠更好地保留圖像的邊界，而這對于特征圖像很重要。全變差最小化算法的優點來自一個很重要的性質，那就是它能夠恢復非稀疏的但具有分片連續性質的圖像。換言之，全變差正則化能夠成功地應用于具有梯度稀疏的信號或圖像。近些年，變量分離思想也為解決全變差最小化問題提供了新的有效解決途徑，這強力推動了單像素相機等硬件的發展。

近幾年，諸多學者在研究CT圖像重建的稀疏表示模型的過程中，取得了一些重要研究成果。比較普遍的做法就是首先建立合適的目標函數，之后尋找求解目標函數的快速算法。文獻[6]的目標函數是最小化基于壓縮感知的全變差正則化，基于全變差最速下降法和凸集投影約束，提出二者相結合的ASD-POCS算法。陳慶貴等人又基于總變差最小化原則，提出了PICS算法。文獻[7]改進了目標函數，在全變差正則化目標函數中加入約束條件重建出的先驗圖像，用引用參數平衡兩個目標函數的權重，之后采用迭代重建算法和最速下降法進行模型求解，提出了先驗約束壓縮感知PICCS(Prior Image Constrained Compressed Sensing)算法。文獻[8]又基于總變差最小化原則，提出了PICS算法。這些算法雖然對圖像的整體信息有較好的恢復，但是仍然存在一些偽影。

Bregman迭代正則化最早由文獻[9]引入圖像處理領域，該算法主要分為線性Bregman算法和分離Bregman算法兩種，后來被廣泛應用于基于壓縮感知的核磁共振成像和基于小波變換的圖像去噪中。文獻[10]將分離Bregman算法用于全變差正則化圖像重建，分析了算法可行性，得到了更快的收斂速度。文獻[11]將分離Bregman算法與先驗圖像相結合，也取得了很好的效果。

近年來有學者采用增廣Lagrangian方法進行研究，該方法應用了加入先驗信息等思想。文獻[12]提出了有序線性增廣Lagrangian(OS-LALM)算法，模型的目標函數為全變差正則化。實驗顯示，有序線性增廣Lagrangian算法能夠有效加速CT圖像重建的收斂速度，同時，當子集個數較多時，能有效減少偽影。但該算法需要首先對系統矩陣進行選擇排序，這大大降低了效率。之后又提出了改進過的松弛的有序線性增廣拉格朗日算法。但是，這些算法或者降低了重建效率，或者存在小范圍的偽影。

為了提升CT圖像重建的質量效率，本文提出外點懲罰函數增廣Lagrangian EPFALM算法。首先，提出了新的稀疏角度CT圖像重建的目標函數并求解；然后，通過仿真實驗，將所提出的算法與前人算法作對比，驗證了EPFALM算法用于稀疏CT圖像重建時的質量和效率優勢。

1 CT圖像典型算法綜述

濾波反投影FBP算法盡管已提出多年，目前依然是廣泛應用的圖像算法之一。其重建算法原理為先對各個視角θ的投影數據進行濾波，然后反投影，累加計算f(x,y)。該算法雖然計算效率高，重建速度快，但對數據完整性要求很高。

代數重建ART算法由Kaczmarz提出，由Gordon等人用于圖像重建領域。Hounsfield的第一臺CT機實際上用的是ART算法。該重建算法原理是給定重建區域初值，將所得投影值殘差沿射線方向反投影從而對圖像進行校正，進而滿足要求、結束迭代。該算法較適用于數據不完全的情況，但對較大圖像重建速度較慢。

分離 Bregman算法由文獻[13]提出，并被文獻[14]用于CT圖像重建。分離 Bregman算法是解決l1范數目標函數最小化的有效算法，已被廣泛應用與圖像處理的各個領域。該算法相對于傳統的算法，能夠較好地抑制由于數據不足帶來的條狀偽影，但是收斂速度較慢，重建時間較長。

本文提出一種新的稀疏角度CT圖像重建算法，即EPFALM算法。該重建算法的原理是采用增廣Lagrangian 罰函數方法，先將約束問題非約束化，再通過在交替方向上更新迭代求解目標函數的等價子問題，從而求得全局最優解。通過計算機仿真實驗，可以發現該算法能夠很好地平衡重建速度與圖像質量，在實際應用中有一定價值。

2 稀疏角度CT圖像重建算法

基于壓縮感知理論的CT圖像重建是近年來研究的熱點，其基本思路是先給出目標函數，然后尋找求解目標函數有效算法。

即定義一個稀疏變換為Φ，重建算法可以通過求解以下約束最小問題而實現：

(1)

但是，該問題是數值不穩定的NP問題?；诖?，文獻[15]于2006年提出了壓縮感知理論，證明了l1范數解可以等價于l0范數解。從此，l1范數重建算法被廣泛應用。因此可以考察l1范數最小化問題：

(2)

文獻[15]還說明精確重建與所在頻域的位置無關，這對于稀疏角度重建問題有著至關重要的意義。

本文提出改進的稀疏角度CT圖像重建算法，即EPFALM算法，算法流程圖如圖1所示。

圖中算式編號在下文演算過程中均有對應算式。算法的具體步驟如下：

步驟1：

(1) 初始化參數v0、ζ0、λ0、β0、μ0、α、ρ；

(2) 令Xk=XXp當條件滿足時停止迭代；

(3) 令αk=ραk，當滿足Armijo條件時停止。

步驟2：

(1) 由式(14)，計算Xk+1；

(2) 由式(19)，計算wk+1；

(3) 由式(21)，計算zk+1；

(4) 由式(7)-式(9)，計算新乘子。

2.1 目標函數的建立

在稀疏角度重建問題中，有效方程數量遠少于未知數個數，屬于不確定問題。文獻[7]因此提出了先驗圖像約束壓縮感知PICCS算法，該算法通過在目標函數中加入先驗圖像，從而彌補圖像的不完整，平滑重建偽影。PICCS設定的目標函數為：

(3)

式中：p為先驗圖像，Φ1、Φ2為任意稀疏變換，α為相對權重。Chen等[7]通過兩個獨立的步驟數值計算：第一步，使用ART重建算法得到先驗圖像p；第二步，將式(3)中差圖像的TV與目標的圖像的TV加權求和。兩步計算交替迭代，直到迭代終止。

受Chen等[7]、Xu等[16]和王歡等[17]的啟發，提出以下優化目標函數：

minXfs.t.AX=b

(4)

2.2 外點懲罰函數增廣Lagrangian乘子法求解

文獻[18]的研究說明，用增廣Lagrangian方法求解約束問題最優化是適合的。因此，本文提出外點懲罰函數增廣Lagrangian乘子法(EPFALM算法)。其本質是連續不斷的改變Lagrangian乘子和懲罰因子來求解各異的Lagrangian函數，即使用無約束最小優化方法得到此Lagrangian函數的極小值點。

對于本文提出的優化目標函數式(4)，應用EPFALM，可以將其化為如下問題：

(5)

式中：wi=Di(X-Xp)2,zi=DiX。

相應的Lagrangian函數為：

LA(wi,zi,X)=

(6)

設wi,k、zi,k、Xi,k代表第k次迭代的最優解。通過式(7)-式(9)更新乘子：

vi,k+1=vi,k-βk(DiXk-wi,k)

(7)

ζi,k+1=ζi,k-βk(Di(Xk-Xp)-zi,k)

(8)

λk+1=λk-μk(AXk-b)

(9)

應用EPFALM，可得增廣Lagrangian函數為：

(10)

minXfs.t.AX=b

如果我們引入wi∈R2替換微分變換DiX，并將兩者之間的差作為懲罰項加入到目標函數中；同時引入zi∈R2替換微分變換Di(X-Xp)，并將兩者之間的差作為懲罰項加入到目標函數中，那么有:

LA(wi,zi,X)=

(11)

將式(11)代入式(10)中,易見目標函數與式(6)相同。

在本算法中，每一步迭代最小化十分重要。文獻[10]將變量分離技術首次引入壓縮感知，把交替最小化算法用于圖像卷積去噪。本文將Lagrangian函數LA(wi,zi,X)的最優化問題轉化為三個子問題，通過在交替方向上求解子問題，進而得到最優解Xk+1，計算過程如下。

在X-子問題中，由于

minuLA(wi,k,zi,k,X)=

(12)

式中：wi,k、zi,k、Xp已知。根據式(11)，可求解Xk+1，此問題等價于X-子問題：

minXQk(X)=

(13)

因為Qk(X)是一個二次函數，因此利用最速下降法求解Qk(X)最小化問題:

Xk+1=Xk-γkdk

(14)

式中：dk=dk(Xk)為Qk(X)的梯度。

μkAT(AX-b)-ATλk

(15)

又因為wi,k、zi,k為上一次迭代的最優解，故不能保證全局最優，所以采用一步最速下降法來逼近Qk(u)的最優解。

在迭代過程中，采取Goldstein步長規則，即:

f(xk)+(1-c)αk▽f(xk)Tdk≤f(xk+αkdk)≤

f(xk)+cαk▽f(xk)Tdk

(16)

驗證αk是否滿足Armijo條件，若α滿足:

f(xk+αdk)≤f(xk)+cα▽f(xk)Tdk

(17)

則取ak←a，否則通過αk=ραk(其中0<ρ<1)來減小步長，直到滿足條件為止。

在w-子問題中，根據式(11)，同樣可求解wi,k+1，此問題等價于w-子問題：

(18)

其解為：

(19)

在z-子問題中，同理，此問題等價于z-子問題：

(20)

其解為：

(21)

值得一提的是，參數的選擇也非常重要，本文采用的是文獻[19]的方法。即在初次迭代給予參數初始值之后，以后每次迭代都根據需要調整參數，這樣極大地提高了算法效率。而其他細節，如滿足Armijo條件，如何選取步長等，也需要在迭代過程中根據相應的公式作出調整。這一點在MATLAB仿真實驗中也有說明。

3 仿真實驗與分析

為了驗證算法的有效性，本文進行仿真實驗，將使用算法對醫學圖像進行重建，并與FBP算法、ART算法、Split-Bregman算法與文獻[20]提出的UIAL算法做比較。測試計算機的配置為Intel i7-5700HQ CPU@2.70 GHz，16 GB內存。實驗使用512×512的Shepp-Logan體模模擬生成的圖像。CT圖像掃描采取平行束方式進行。實驗中所有程序均由MATLAB R2014a編程實現。

為了定量分析圖像重建的效果，文獻[3]提出，可利用相對均方誤差(RRMSE)和條紋指標(SI)來量化重建的效果，其中：

(22)

(23)

式中:X為重建的圖像，Xref為原始圖像。顯然，RRMSE和SI的值越小，表示重建效果越好。

算法中，設定參數如下：

v0=0ζ0=0λ0=0β0=1μ0=256

α=0.5ρ=0.6wi 0=0zi,0=0X0=ATb

每個算法迭代次數為1 000次。

其中，部分參數可以通過式(7)-式(9)更新計算。

vi,k+1=vi,k-βk(DiXk-wi,k)

ζi,k+1=ζi,k-βk(Di(Xk-Xp)-zi,k)

λk+1=λk-μk(AXk-b)

另外經計算得知，當α=0.5時，RRMSE與SI的值最小，因此取α=0.5。當α>1時，RRMSE、SI遠高于α<1時得到的值，這也說明了本方法對于CT圖像重建是相當有用的。

實驗采取的肺部、肝部理想Shepp-Logan體模圖像如圖2-圖4所示。

圖2 實驗圖像

兩幅圖像用不同算法重建的實驗結果為：

圖3 肺部重建圖像

圖4 肝部重建圖像

圖3、圖4分別是肺部CT圖像、肝部CT圖像在稀疏角度采樣的情況下重建，五種算法分別為FBP算法、ART算法、Split-Bregman算法、文獻[20]提出的UIAL算法和EPFALM算法。從圖像中定性分析可知，FBP算法重建圖像有嚴重的偽影，Split-Bregman算法重建的圖像細節模糊，而ART算法、UIAL算法和EPFALM算法重建的圖像沒有偽影，細節較清晰。

下面用評價指標來量化比較各算法下的重建結果。

表1 肺部圖像重建效果

表2 肝部圖像重建效果

在重建效果上，EPFALM算法與FBP算法相比，在兩幅圖像中，RRMSE指標值分別降低97.6%、98.6%，SI指標值分別降低97.5%、92.0%；與Split-Bregman算法相比，RRMSE指標值分別降低94.6%、98.2%，SI指標值分別降低96.4%、98.2%；與UIAL算法相比，RRMSE指標值分別降低74.6%、89.1%，SI指標值分別降低84.9%、85.3%。在重建時間上，EPFALM算法與ART算法相比，降低87.1%；與Split-Bregman算法相比，降低76.6%。

4 結 語

醫用CT在有著便捷、精確等優勢的同時，也伴隨著射線輻射對患者的傷害。因此，本文提出了研究稀疏角度CT圖像重建的EPFALM算法，有利于在心臟成像、血管成像、脊柱成像等快速成像應用中獲得更好的成像效果。

基于Chen等[7]的研究，本文提出了新的目標函數，并采取新的算法進行求解。新的目標函數加入了先驗信息，在原全變差正則化的基礎上加入上一步迭代重建的圖像，用參數衡量兩個目標函數的權重。

用MATLAB仿真模擬稀疏角度采樣時，比較了FBP算法、ART算法、Split-Bregman算法、UIAL算法、EPFALM算法重建的質量速度，應用的兩幅體模圖像為肺部圖像和肝部圖像。仿真實驗結果表明，本算法重建的圖像很好地抑制了偽影，細節也較為清晰，有一定的實際應用價值。