尹寶才 郭曉明 施云惠 丁文鵬

（北京工業(yè)大學(xué)城市交通學(xué)院，多媒體與智能軟件技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京，100124）

引 言

目前研究表明，壓縮感知理論可以在圖像少量觀測值下完成對(duì)該圖像的重建［1-3］。通常，圖像具有塊結(jié)構(gòu)特性并且可以在確定小波基下稀疏表示［4，5］。TV-Wavelet-L1（TVWL1）模型采用全變分正則化約束（Total-variation，TV）和小波L1范式正則化約束的組合稀疏表示方式，使得基于壓縮感知理論的圖像重建算法精確度得到提高［6，7］。

對(duì)應(yīng)于圖像的分解和重建過程，圖像的稀疏表示方法分為分析模型和綜合模型。例如：綜合的稀疏表示方法是一個(gè)圖像x∈Rd可以描述為x＝Dα，這里D∈Rd×n是變換矩陣，α是x的稀疏表示系數(shù)，‖α‖0＝k?n。而分析稀疏表示方法是向量Ωx可以被稀疏化，Ω是一個(gè)線性算子，相對(duì)稀疏性可以定義為向量Ωx的零元素的數(shù)量。

當(dāng)綜合模型被廣泛學(xué)習(xí)研究時(shí)，分析模型的稀疏表示方法卻很少被關(guān)注［8］。但是，最近幾年，分析稀疏表示方法得到了人們的廣泛關(guān)注。基于TVWL1模型，相較于分析方式，綜合稀疏表示方式的TV正則化約束更為復(fù)雜并且求解過程更加困難。所以，通常應(yīng)用基于分析稀疏表示方式的TVWL1模型［9］來進(jìn)行圖像重建工作，該模型如下

對(duì)于簡單正則化約束，例如單獨(dú)的TV或L1范數(shù)約束，迭代收縮閾值算法（Iterative shrinkagethresholding algorithm，ISTA）和快速迭代收縮閾值算法都可以進(jìn)行高效的圖像重建［10］。但是它們不適用于TV和L1范數(shù)組合的正則化約束。并且，TV和L1范數(shù)是不平滑的，它使得求解問題（1）變得更加困難。

目前，一些算法［11，12］已經(jīng)應(yīng)用到求解TVWL1模型。組合分裂算法（Composite splitting algorithm，CSA）和快速組合分裂算法（Fast composite splitting algorithm，F(xiàn)CSA）［9］是兩個(gè)高效的求解問題（1）的算法，它們應(yīng)用于圖像重建，并且有著出色的圖像重建質(zhì)量。這兩個(gè)算法都是把原問題分解為兩個(gè)子問題：L1正則化約束和TV正則化約束。然后，采用已有技術(shù)分別求解這兩個(gè)子問題來完成圖像重建。然而，在少量觀測值的情況下，這兩個(gè)算法不能夠完成對(duì)MR圖像的高質(zhì)量重建。

本文提出了一種新的求解TVWL1模型的圖像重建算法，該算法將問題（1）分解為幾個(gè)子問題，并利用分析稀疏表示特性構(gòu)建子問題的求解算法，再通過組合子問題的解來獲得重建圖像。實(shí)驗(yàn)結(jié)果表明，與現(xiàn)有的一些算法相比，本文提出的新算法可以顯著提高重建圖像的主客觀質(zhì)量。

1 TVWL1模型求解算法

1.1 近似梯度算法

給出一個(gè)連續(xù)的凸函數(shù)g（x）和任意標(biāo)量ρ＞0，凸函數(shù)g的prox算子定義如下［10，13］

無約束優(yōu)化可以分為兩部分

式（3）中f（x）是平滑的凸函數(shù)，g（x）是包含prox算子的凸函數(shù)。近似梯度算法定義如下

▽表示梯度算子，k是迭代步長。

1.2 組合分裂算法

可用于求解問題（1）的CSA算法［9］是一種基于變量和算子分裂技術(shù)的方法。它把組合正則化問題（1）分解為兩個(gè)子問題：（1）把變量x分解為兩個(gè)變量｛xi｝i＝1，2；（2）對(duì)｛xi｝i＝1，2獨(dú)立執(zhí)行算子分裂求解；（3）獲得由｛xi｝i＝1，2線性組合而成的重建結(jié)果x。FCSA是CSA加入加速過程的改進(jìn)算法。算法1給出了用于MR圖像重建的FCSA/CSA。

算法1 快速組合分裂算法

（用rk＋1＝xk替換步驟（6）和（7）為CSA）

為了更加高效地求解問題（1），關(guān)鍵是如何應(yīng)用相對(duì)稀疏度限制去改善圖像的重建質(zhì)量。

2 基于分析稀疏表示方法的圖像重建算法

在學(xué)習(xí)求解問題（1）之前，首先回顧一下分析稀疏表示模型中的貪心分析追蹤算法［14］。

2.1 貪心分析追蹤算法

GAP算法是求解下面分析模型的L1最小化問題

式中x采用分析稀疏表示方式。GAP的核心思想是找到Φx的零元素的位置，并且從Φ中移除與Φx非零元素對(duì)應(yīng)的行，用Φ的子集去執(zhí)行迭代計(jì)算。算法2給出了GAP算法的詳細(xì)描述。

算法2 貪心分析追蹤算法

表達(dá)式（6）解析式為

GAP算法需要一個(gè)大的內(nèi)存空間去存儲(chǔ)大尺度矩陣并且在計(jì)算大尺度矩陣的逆矩陣時(shí)，有著很高運(yùn)算復(fù)雜度，這樣就很難實(shí)現(xiàn)大尺度圖像的重建。

2.2 基于分析稀疏表示算法

本文提出的算法把圖像重建問題分解為兩個(gè)子問題：由近似梯度算法求解的TV正則化約束問題和由貪心分析追蹤算法求解的L1正則化約束問題。這個(gè)算法被命名為基于分析稀疏表示算法（Analysis-sparse-based algorithm，ASBA）。圖像重建算法框架如算法3所示。

算法3 基于分析稀疏表示算法

ASBA是一種迭代算法，它分為近似梯度求解和分析追蹤求解兩個(gè)過程。TV正則化約束問題很容易通過類似于CSA的近似梯度算法求解。而本算法的核心部分是由分析追蹤算法求解L1正則化約束問題，即求解問題（6）的過程。

實(shí)際上，問題（6）可以轉(zhuǎn)換為下面的形式

式中：A＝RTR＋，A是正定對(duì)稱矩陣。所以，式（7）可以通過共軛梯度算法（Conjugate gradient algorithm，CG）高效求解。于是，結(jié)合CG算法和GAP算法得到一個(gè)新的算法——基于共軛梯度的GAP算法（Conjugate-gradient-based GAP，CGGAP），它在算法4中詳細(xì)給出。

算法4 基于共軛梯度GAP算法

算法4中，x0是由近似梯度過程計(jì)算得到初始值是重建圖像；n是共軛梯度算法迭代次數(shù)；K是GAP算法迭代次數(shù)；共軛梯度算法（CG）在算法5中給出；其他參數(shù)在第一部分有詳細(xì)描述。

算法5 共軛梯度算法

在算法4和算法5中，由于部分傅里葉變換矩陣R和小波變換矩陣Φ尺度很大并且非常稠密，與其相關(guān)的計(jì)算需要很大的內(nèi)存才能順利完成。所以，在程序?qū)崿F(xiàn)過程中，采用把大矩陣分塊存儲(chǔ)技術(shù)來降低運(yùn)算對(duì)內(nèi)存的需求。

實(shí)驗(yàn)部分將展示重建圖像的主客觀質(zhì)量。實(shí)驗(yàn)結(jié)果表明，本文提出的新算法的圖像重建質(zhì)量要明顯好于以前提出的一些算法。

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

實(shí)驗(yàn)的軟件環(huán)境是 MATLAB7.6，實(shí)驗(yàn)代碼是基于文獻(xiàn)［12］提供的程序編寫，實(shí)驗(yàn)內(nèi)容是觀測和重建自然圖像和MR圖像。實(shí)驗(yàn)圖像包括自然圖像Lena和Boat，醫(yī)學(xué)圖像 Heart，Brain，Chest和Artery，它們的尺度均為128＊128。觀測方法采用部分傅里葉觀測。采樣率定義為m/n。噪聲模型是由 Matlab產(chǎn)生的高斯噪聲，定義為σ＊randn（m，1），σ＝0.01。正則化參數(shù)α和β分別設(shè)置為0.001和0.035。為了實(shí)驗(yàn)的公平性，在對(duì)比不同算法重建結(jié)果時(shí)，每個(gè)算法的迭代次數(shù)設(shè)置為10。小波變換矩陣Φ設(shè)置為2D小波。實(shí)驗(yàn)將比較每個(gè)算法的主客觀圖像重建質(zhì)量，客觀評(píng)價(jià)標(biāo)準(zhǔn)使用峰值信噪比（Peak-signal-to-noise-ratio，PSNR）。實(shí)驗(yàn)算法是ASBA，CSA和FCSA。

3.1 ASBA圖像重建效果

圖1是在不同采樣率下ASBA重建自然圖像Lena效果圖，圖1（a～d）對(duì)應(yīng)的重建圖像PSNR分別為21.68，28.72，32.34和38.62dB。圖2是在不同采樣率下ASBA重建MR圖像Brain效果圖，圖2（a～d）對(duì)應(yīng)的重建圖像PSNR分別為18.13，24.42，27.85和35.21dB。從圖中可以看出ASBA可以很好地求解TVWL1模型，并且隨著采樣率的提高，ASBA重建圖像的主觀質(zhì)量明顯提高。

圖1 不同采樣率下ASBA重建Lena圖像效果Fig.1 The results of the reconstructed Lena images by ASBA at different sample rates

圖2 不同采樣率下ASBA重建Brain圖像效果Fig.2 The results of the reconstructed Brain images by ASBA at different sample rates

3.2 圖像重建主觀評(píng)價(jià)結(jié)果

通過求解問題（1）來重建圖像，從視覺效果上對(duì)本文提出的新算法ASBA和已有算法CSA＼FCSA進(jìn)行評(píng)估。圖3，4給出了采樣率在25%時(shí)，ASBA，F(xiàn)CSA和CSA的圖像重建質(zhì)量對(duì)比。從圖中可以看出，無論是自然圖像還是醫(yī)學(xué)圖像，本文提出的ASBA重建的圖像主觀質(zhì)量要明顯好于FCSA和CSA。

圖3 采樣率25%不同算法重建Lena圖像Fig.3 Lena image is reconstructed by different algorithm at sampling ratio near 25%

圖4 采樣率25%不同算法重建Brain圖像Fig.4 Brain image is reconstructed by different algorithm at sampling ratio near 25%

圖3（a）是Lena圖像原圖，圖3（b～d）分別是由ASBA，F(xiàn)CSA，CSA重建得到的Lena圖像，它們對(duì)應(yīng)的PSNR分別是32.35，28.06和28.01dB。圖4（a）是Brain圖像原圖，圖4（b～d）分別是由ASBA，F(xiàn)CSA，CSA重建得到的Brain圖像，它們對(duì)應(yīng)的PSNR分別是27.92，24.19和23.53dB。

3.3 圖像重建客觀評(píng)價(jià)結(jié)果

圖5是不同采樣率下不同算法重建圖像的PSNR對(duì)比曲線圖。從圖中可以看出，在不同采樣率下，ASBA重建圖像的PSNR與FCSA和CSA相比，可以提升2～5dB。表1列出的實(shí)驗(yàn)結(jié)果顯示ASBA在所有采樣率下的圖像重建質(zhì)量最好。

圖5中曲線圖上對(duì)應(yīng)的四個(gè)采樣率點(diǎn)是11%，16%，25%和43%。圖5（a，b）圖中的3條曲線分別對(duì)應(yīng)ASBA，F(xiàn)CSA和CSA三種算法重建圖像Lena和Brain的PSNR-Sampling Ratio曲線。

圖5 重建圖像PSNR-Sampling Ratio曲線Fig.5 PSNR-Sampling Ratio curve of the reconstruction images

表1 不同采樣率下不同算法重建圖像實(shí)驗(yàn)結(jié)果Table 1 Experimental results for image reconstruction over different sampling ratios with different methods

4 結(jié)束語

本文提出了一個(gè)基于分析稀疏表示的圖像重建算法。該算法將TVWL1模型分解為TV和L1正則化約束兩個(gè)子問題并分別進(jìn)行求解，從而達(dá)到高效圖像重建的目的。實(shí)驗(yàn)部分，通過與CSA和FCSA重建圖像的比較，可以看出新提出的算法可以高質(zhì)量重建MR圖像和自然圖像。

［1］Candès E J，Tao T.Near-optimal signal recovery from random projections：Universal encoding strate-gies［J］.Information Theory，IEEE Transactions on，2006，52（12）：5406-5425.

［2］Candès E J，Romberg J，Tao T.Robust uncertainty principles：Exact signal reconstruction from highly incomplete frequency information［J］.Information Theory，IEEE Transactions on，2006，52（2）：489-509.

［3］Donoho D L.Compressed sensing［J］.Information Theory，IEEE Transactions on，2006，52（4）：1289-1306.

［4］Candès E，Romberg J.Sparsity and incoherence in compressive sampling［J］.Inverse Problems，2007，23（3）：969.

［5］Romberg J.Imaging via compressive sampling［J］.Signal Processing Magazine，IEEE，2008，25（2）：14-20.

［6］Lustig M，Donoho D L，Santos J M，et al.Compressed sensing MRI［J］.Signal Processing Magazine，IEEE，2008，25（2）：72-82.

［7］Yang J，Zhang Y，Yin W.A fast alternating direction method for TVL1-L2signal reconstruction from partial Fourier data［J］.Selected Topics in Signal Processing，IEEE，2010，4（2）：288-297.

［8］Elad M，Milanfar P，Rubinstein R.Analysis versus synthesis in signal priors［J］.Inverse Problems，2007，23（3）：947-968.

［9］Huang J，Zhang S，Metaxas D.Efficient MR image reconstruction for compressed MR imaging［J］.Med-ical Image Analysis，2011，15（5）：670-679.

［10］Beck A，Teboulle M.A fast iterative shrinkagethresholding algorithm for linear inverse problems［J］.SIAM Journal on Imaging Sciences，2009，2（1）：183-202.

［11］Lustig M，Donoho D，Pauly J M.Sparse MRI：The application of compressed sensing for rapid MR imaging［J］.Magnetic Resonance in Medicine，2007，58（6）：1182-1195.

［12］Ma S，Yin W，Zhang Y，et al.An efficient algorithm for compressed MR imaging using total variation and wavelets［C］∥Computer Vision and Pattern Recognition，CVPR 2008.［S.L.］：IEEE，2008：1-8.

［13］Beck A，Teboulle M.Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems［J］.Image Processing，IEEE Transactions on，2009，18（11）：2419-2434.

［14］Nam S，Davies M E，Elad M，et al.The cosparse analysis model and algorithms［J］.Applied and Computational Harmonic Analysis，2013，34（1）：30-56.