低秩矩陣在CT圖像重建中的應用*

馬海英，宣士斌，b，c，向順靈

(廣西民族大學 a.信息科學與工程學院；b.廣西混雜計算與集成電路設計分析重建實驗室；c.中國-東盟研究中心，廣西 南寧 530006)

?

低秩矩陣在CT圖像重建中的應用*

馬海英a，宣士斌a，b，c，向順靈a

(廣西民族大學 a.信息科學與工程學院；b.廣西混雜計算與集成電路設計分析重建實驗室；c.中國-東盟研究中心，廣西 南寧 530006)

CT圖像重建是醫學影像學的重要研究課題，但由于噪聲對醫學 CT 圖像的影響比較大，為了在不犧牲圖像精度和空間分辨率的情況下，重建出噪聲含量最低的圖像，就要選擇合適的去噪方法對圖像進行預處理.針對于此，筆者提出一種新的CT圖像重建算法，重建過程分成兩個步驟：首先用低秩矩陣加權核范數最小化(WNNM)進行圖像去噪，再用低秩矩陣分解(LRMD)更新CT圖像.實驗結果表明，提出的方法具有較強的細節保持能力，低秩矩陣的特性簡化計算過程，降低算法復雜度，同時保證了重建圖像的去噪效果.

低秩矩陣；核范數；CT圖像重建

0 引 言

近年來，CT圖像重建的統計學習方法發展迅速，這是由于CT掃描時需要低劑量的X射線輻射的同時要保留高質量的重建圖像.然而，統計學方法計算量大、耗時長的特點限制了它的實際應用，為了加速統計方法，許多優化技術被提出，這些算法包括：迭代閾值法(Iterative shrinkage/thresholding algorithm,IST)[1]、兩步迭代閾值法(Two step iterative shrinkage/thresholding algorithm,TwIST)[2]、快速迭代閾值法(Fast Iterative shrinkage/thresholding algorithm,FISTA)[3]；分裂 Bregman 方法(Split Bregman algorithm)[4]、Bregman 算子分裂方法(Bregmanized operator splitting，BOS)[5]；低秩矩陣恢復(low-rank matrix recovery，LRMR)技術[6]；Tao[7]等在交替最小化方法的基礎上提出了交替方向乘子法(Alternating direction method of multipliers,ADMM)[8]；低秩矩陣分解(low-rank matrix decomposition，LRMD)技術是近幾年迅速發展起來的一種高維數據分析工具，并在協同過濾(collaborative filtering)、控制(control)、遙感(remote sensing)、量子態層析成像(quantum state tomography)、機器學習和計算機視覺等領域得到廣泛應用.近似低秩矩陣，旨在從它的退化視圖中恢復潛在低秩矩陣，它在計算機視覺和機器學習中有較大應用.例如，通過人臉面部圖像形成矩陣的低秩特性允許我們重建損壞的臉部[9].網飛公司客戶數據矩陣就被認為是一種低秩矩陣，因為客戶的選擇大部分受一些常規因素的影響[10].通過靜態相機捕獲的視頻片段有一個清晰的低秩特性，基于背景建模和前景抽取[11]可以被統計出來.在自然圖像中通過非局部相似塊形成矩陣也是低秩特性.由于凸凹優化技術的迅速發展，近年來在近似低秩矩陣中有一系列的研究，同時提出許多重要模型和算法.

目前低秩矩陣技術主要包含矩陣填充(matrix completion,MC)[12]、魯棒主成分分析(robust principle component analysis,RPCA)[13]和低秩表示(low-mnk representation,LRR)[14]三個方面的內容.該技術的理論基礎是矩陣的仿射秩最小化理論，即在給定線性方程組約束下，以矩陣的秩作為測度對目標矩陣進行分析和處理.然而，秩最小化問題在理論上是NP難(Non-deterministic Polynomial Hard,NP Hard)的.類似于壓縮感知(compressive sensing,CS)中用l1范數代替l0范數[15]，在拓展了約束等距性(restricted isometry property，RIP)條件后，核范數(矩陣的所有奇異值的和)被用來代替秩函數作為原優化問題的目標函數[16].事實上，壓縮感知與秩最小化是密切相關的.當矩陣為對角矩陣時，秩最小化問題就是退化為在矩陣的子空間中找一個最稀疏向量的問題.此時，矩陣的奇異值的和就等同于矩陣的對角元的絕對值之和，即求解核范數最小化問題與l1范數最小化問題是等價的.

由于低秩矩陣分解的凸松弛問題，核范數最小化成為近年來研究的重點.標準核范數最小正則化(NNM)[17]每一個奇異值等同于追求目標函數的凸性問題.然而，這也使其在處理許多實際問題中(如，圖像去噪、圖像恢復)很大程度上限制了它的性能和靈活性，因為有些奇異值有明顯的物理意義，應該區別對待.文中，我們研究了加權核范數最小化(WNNM)問題，此處的奇異值被分配了不同的權重，然后利用圖像非局部自相似性將提出的WNNM算法進行圖像去噪，實驗證明提出的WNNM算法在圖像定量測量和視覺感知質量方面明顯高于許多先進的去噪算法(如BM3D).

文中，我們提出了基于低秩矩陣加權核范數最小化的去噪模型，并將該模型應用于CT圖像去噪，同時將基于低秩矩陣分解應用于CT重建，在重建模型中利用前面提出的圖像去噪模型進行圖像去噪，建立CT數據重建數學模型，利用傅里葉變換和低秩矩陣的特性簡化計算過程，降低算法復雜度.實驗表明本文提出的方法具有較好的去噪效果，且為CT重建中的圖像去噪步驟提供了堅實的基礎，同時具有較強的細節保持能力.

1 相關工作

(1)

我們所要求解的去噪后的圖像為px.其中，p代表對動態圖像投影，‖·‖F為F范數，‖·‖1為l1范數.

1.1 加權核范數進行圖像去噪

核范數最小化(NNM)是一個凸性最優問題.由于許多低秩矩陣能通過NNM方法得到很好的恢復并能高效的解決，因此核范數最小化廣泛應用于低秩矩陣最優化問題中，它能通過F范數測量觀測數據矩陣Y和潛在數據矩陣X的區別，通過奇異值的軟閾值法得到一個分析解.由于相同的軟閾值將會應用到所有奇異值中，NNM方法顯然不太合理，因為不同的奇異值可能有不同的價值，因此他們需要區別對待.為達到這個目的，我們使用加權核范數來正則化X.下式為加權核范數最小化式化：

(2)

(3)

顯然，現在的關鍵問題是權重向量w的確定.對于自然圖像，我們有普遍的先驗知識，即pXj的較大奇異值比較小的更重要，因為他們代表Xj主要部分的能量.去噪應用中，奇異值越大，他就應該縮減得越小.因此，權重分配給σi(Xj)，Xj的第i個奇異值應該和σi(Xj)成反比，我們讓：

(4)

c>0是常數，n是Yj中相似塊的數目，ε=10-16是防止除數為0.

假定噪聲能量跨越基底U和V的每個子空間是均勻分布的，然后最初σi(Xj)估計可以寫成如下：

(5)

σi(Yj)表示Yj第i個奇異值.通過將以上的程序應用到每個塊中然后聚集所有的塊，就能重建圖像x.實際操作中，我們可以多次運行以上程序以提高去噪質量.整體的去噪算法在算法1中總結出來：

算法1：基于WNNM的圖像去噪

輸入：噪聲圖像y

2)for k=1:K do

4)for y(k)中的每個塊yjdo

5)找到相似塊組Yj

6)評估權重向量w

7)奇異值分解 [U,∑，V]=SVD(Yj)

9)結束

11)結束

2 低秩矩陣分解更新CT圖像

在重建模型中利用前面提出的圖像去噪模型進行圖像去噪，建立CT數據重建數學模型，然后利用低秩矩陣分解的特性將得到的清晰圖像運用到錐束CT成像(CBCT)[19]圖像重建中.

Cai等人[19]將時間作為一個維度，利用序列CBCT圖像中潛在的周期性或重復性等時間上的相關性建模并求解.首先將應用于所有不同投影時刻的圖像xi以向量的形式按列依次排成一個矩陣X.矩陣的每一列代表一幅待重建的CBCT圖像，矩陣的列數即為投影的次數.該算法的核心思想是矩陣X的秩遠小于投影的次數，因此對其進行矩陣的乘法分解X=LR.X中的圖像性質分別體現在矩陣L的稀疏性和矩陣R的近似周期性上.首先，矩陣L的列是對矩陣R的秩的約束，無形中對CBCT中所有圖像加了一個時間相容性條件.事實上，矩陣L的每一列都是一幅CBCT圖像，因此L是可被用于表示矩陣X中的所有圖像的一組基.其次，矩陣R的行是矩陣X在基L下的系數，具有一定的周期性或重復性.

可以將CBCT重建看成一個如下最優化問題：

s.t.‖p(LR)-Y‖≤σ2

(6)

其中，p代表對動態圖像投影，Y為投影數據，σ為誤差控制項.考慮到L和R分別具有稀疏性和潛在的周期性等先驗信息，分別采用了在小波緊框架[21]下的稀疏算法d和傅里葉變換f.‖·‖為l1范數，λ為平衡參數.

2.1 算法

首先，Cai使用split Bregman[22]方法來解決這個優化問題，首先引入兩個輔助變量C和D，那么等式(6)就等價于下式：

s.t.p(LR)=F,C=dL,D=fR

(7)

增廣拉格朗日式即：

(8)

此處<·,·>表示內積，Z,Z1和Z2表示拉格朗日乘子.合理的固定Z，Z1，Z2，再通過最小化E(C,D,L，R,Z,Z1,Z2)就能找到最優的C，D，L，R.因此，關鍵是確定Z，Z1，Z2.在增廣拉格朗日算法中，我們使用下式進行交替極小化算法：

(9)

這4個子問題通過軟閾值法和線性方程求解器(如共軛梯度法)求解.這個算法總結到如下算法中，Γ是軟閾值運算符，定義為[Γ(A)]ij=sign([A]ij)·max{‖A‖ij-,0}.

1)通過如下最小化迭代.

返回.

3)Z(k+1)=Z(k)+(p(L(k+1)R(k+1))-F)

6)k=k+1，返回第1步.

3 圖像質量客觀評價標準

為了對圖像去噪效果進行評價，采用峰值信噪比(PSNR)進行客觀評價.令大小為M×N的原圖和有噪聲圖像分別為x和y，PSNR值計算公式如下：

(10)

其中x(i,j)和y(i,j)分別表示圖像x和y在位置(i,j)處的幅值.PSNR的值越高表示圖像和原圖越相似，去噪效果越好.

為了對重建圖像效果進行評價，本文采用均方根誤差(RMSD)進行評價.

(11)

RMSD值越小表示圖像和原圖越接近，重建效果越好.

4 實驗結果與分析

程序仿真基于VirtualBox centos 6.6-32bit 系統下的Matlab編程環境，在CPU為AMD Athlon2.99 GHz П X2 B24 處理器下，內存為1.75 GB的PC機上運行.視頻的分辨率均為480×320.

為了驗證文中所提基于低秩矩陣加權核范數最小化的圖像去噪算法的有效性，對CT圖像進行了仿真測試，并采用峰值信噪比PSNR (Peak Signal Noise Ration)評價標準作為評價圖像去噪的標準.

圖1 CT圖像在WNNM下的去噪效果

為了驗證加權核范數最小化(WNNM)去噪效果的優越性，本文對比算法是NNM和BM3D，所用評價標準為PSNR.從表1可以看出，對于CT圖像的去噪結果，本文所提算法PSNR值比NNM和BM3D (精確已知噪聲標準差的情況)高，因此本文所提去噪算法的去噪結果會比NNM、BM3D好.

為了驗證本文提出算法(LRMD-WNNM)的可行性，由于每一個錐束CT(CBCT)圖像的重建都是基于相應的瞬時投影.該算法首先進行去噪預處理，然后有效地利用潛在的周期性或重復性等時間上的相關性建模并求解，如圖2所示，容易發現本文算法能捕獲解剖圖的運動狀態并恢復得其結構，同時能重構出高分辨率的CT圖像.

表1 不同方法的去噪效果(PSNR)

圖2 提出算法(LRMD-WNNM)對NCAT體模的重建過程效果圖

為了驗證本文提出的算法(LRMD-WNNM)相對于Cai提出的簡單的LRMD方法更具優越性，分別將這兩種方法進行CT重建，由圖3可知，提出的算法(LRMD-WNNM)更能有效地去除偽影，具有較好的去噪能力，重建效果更清晰.這主要是由于CT圖像在低秩矩陣分解之前進行去噪預處理，因此重建的圖像更接近原始圖像.

圖3 原始圖像、簡單的低秩矩陣分解(LRMD)和提出的算法(LRMD-WNNM)對NCAT體模的重建效果比較

為了驗證提出算法的優越性，實驗將三種算法在相同條件下進行的CT圖像重建效果比較.圖4將濾波反投影法(FBP)、可分二次迭代(SQS)和提出的算法(LRMD-WNNM)進行CT圖像重建效果的比較，圖5將這三種方法在前10次迭代的RMSD進行了比較.由圖可知LRMD-WNNM算法的重建效果優于FBP、SQS算法，這主要是因為LRMD-WNNM算法相對于SQS算法具有更好的穩定性且預先進行了更好地去噪處理，這就使得提出的算法在更新圖像的過程中降低了圖像偽影，低秩矩陣的特性簡化計算過程，降低算法復雜度，提高了算法的收斂速率，同時也降低了RMSD值，具有更優越的重建精度.

圖4 濾波反投影法(FBP)、可分二次代理(SQS)與提出的算法(LRMD-WNNM)對NCAT體模的重建效果比較

迭代次數

5 結語

本文提出了基于低秩矩陣加權核范數最小化的去噪模型，并將該模型應用于CT圖像去噪，同時將基于低秩矩陣分解應用于CT重建，在重建模型中利用前面提出的圖像去噪模型進行圖像去噪，建立CT數據重建數學模型，利用傅里葉變換和低秩矩陣的特性簡化計算過程，降低算法復雜度.實驗表明本文提出的方法具有較好的去噪效果，且為CT重建中的圖像去噪步驟提供了堅實的基礎，同時具有較強的細節保持能力.

盡管這種可行性實驗取得了成功，但是對于CBCT的臨床應用仍然存在一些實際問題.首先，當把CB幾何模型建模成一個立體CBCT圖像時，由于涉及極大的數據就會帶來一些潛在問題，計算效率也會降低.這些問題能通過一些更有力的計算平臺(如計算GPU)得到一定緩解.降低圖像質量的另一問題是呼吸模型的奇異性.這種方法在CBCT圖像中有效地利用時間相關的周期性，然而它在患者不規則呼吸運動情況下(如咳嗽)時將有所下降，將來可以對那些不規則運動的情形在仿真數據中進行更深遠的研究.

[1] Figueiredo M A T,Nowak R D.An EM algorithm for wavelet-based image restoration[J].IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society,2003,12(8):906-916.

[2] J M Bioucas-Das,M Figueiredo.A new twist:two step iterative shrinkage /thresholding algorithms for image restoration [J].IEEE Transactions on Image Processing,2007,16(12):2992-3004.

[3] A Beck,M Teboulle.Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems [J].IEEE Transactions on Image Processing,2009,18(11):2419-2434.

[4] T Goldstein,S Osher.The split Bregman algorithm for l1 regularized problems [J].SIAM Journal on Imaging Sciences,2009,2(2):323-343.

[5] X Zhang,M Burger,X Bresson,S Osher.Bregmanized nonlocal regularization for deconvo-lution and sparse reconstruction [R].UCLA CAM Report,2009.

[6] Chandrasekaran V,Recht B,Parrilo P A,et al.The Convex Geometry of Linear Inverse Problems[J].Foundations of Computational Mathematics,2012,12(6):805-849.

[7] Y Wang,J Yang,W Yin,Y Zhang.A new alternating minimization algorithm for total varia-tion image reconstruction [J].SIAM Journal on Imaging Science,2008,1(3):248-272.

[8] M Tao,J Yang,B He.Alternating direction algorithms for total variation deconvolution in image reconstruction [R].Available at Optimization-online,2009.

[9] Zheng Y,Liu G,Sugimoto S,et al.Practical low-rank matrix approximation under robust L1-norm[C]// IEEE Conference on Computer Vision and Pattern Recognition，2012:1410-1417.

[10] Salakhutdinov R,Srebro N.Collaborative Filtering in a Non-Uniform World:Learning with the Weighted Trace Norm[J].Advances in Neural Information Processing Systems,2010:2056- 2064.

[11] Mu Y,Dong J,Yuan X,et al.Accelerated low-rank visual recovery by random projection[C]// IEEE Conference on Computer Vision and Pattern Recognition，IEEE Computer Society,2011:2609-2616.

[12]Gamper U,Boesiger P,Kozerke S.Compressed sensing in dynamic MRI[J].Magnetic Resonance in Medicine Official Journal of the Society of Magnetic Resonance in Medicine,2008,59(2):365-373.

[13]Jung H,Sung K,Nayak K S,et al.k-t FOCUSS:a general compressed sensing framework for high resolution dynamic MRI.[J].Magnetic Resonance in Medicine,2009,61(1):103-116.

[14]Ravishankar S,Bresler Y.MR image reconstruction from highly undersampled k-space data by dictionary learning.[J].IEEE Transactions on Medical Imaging,2011,30(5):1028-1041.

[15] Dabov K ,Foi A ,Katkovnik V ,et al.Image Denoising by Sparse 3-D Transform-Domain Collaborative Filtering[J].IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society,2007,16(8):2080-2095.

[16] Ak?akaya M,Basha T A,Goddu B,et al.Low-dimensional-structure self-learning and thresholding:Regularization beyond compressed sensing for MRI Reconstruction[J].Magnetic Resonance in Medicine Official Journal of the Society of Magnetic Resonance in Medicine,2011,66(3):756-767.

[17] Lin Z,Liu R,Su Z.Linearized Alternating Direction Method with Adaptive Penalty for Low-Rank Representation[J].Advances in Neural Information Processing Systems,2011:612-620.

[18] Dong W,Shi G,Li X.Nonlocal image restoration with bilateral variance estimation:a low-rank approach.[J].Image Processing IEEE Transactions on,2013,22(2):700-711.

[19] Cai J F,Jia X.Cine cone beam CT reconstruction using low-rank matrix factorization:algorithm and a proof-of-princple study[J].IEEE Transactions on Medical Imaging,2012,33(8):1581 - 1591.

[20] Ron A,Shen Z.Ane systems in L2(IR d ):the analysis of the analysis operator[J].Journal of Functional Analysis,1997,148(2):408-447.

[21] Gao X.Penalized weighted low-rank approximation for robust recovery of recurrent copy number variations[J].Bmc Bioinformatics,2015,16(1):1-14.

[22] The split Bregman method for L1 regularized problems[C]// SIAM J.Imag.Sci，2009.

[責任編輯 蘇 琴]

[責任校對 黃招揚]

Low-Rank Matrix Technology for CT Image Reconstruction

MA Hai-yinga,XUAN Shi-bina,b,c,XIANG Shun-linga

(a.CollegeofInformationScienceandEngineering;b.GuangxiKeyLaboratoryofHybridComputationandICDesignAnalysis;c.TheChina-ASEANStudyCenterofGuangxiUniversityforNationalities,GuangxiUniversityforNationalities,Nanning530006,China)

Computed tomography (CT)image reconstruction is an important research subject in field of medical imaging.But as the heavily influence of the noise in medical CT image,We must choose appropriate denoising method for image preproce-ssing to get the lowest noise images,while without sacrificing image precision and spatial resolution.To this problem,this paper proposes a new CT image reconstruction algorithm,the reconstruction process has two steps:first,the low rank weighted nuclear matrix norm minimization(WNNM)which is applied to image denoising.Then a low-rank decomposition of matrix which is used to update CT images.Experimental results show that the proposed method has strong ability to keep the details of the CT images,the characteristics of low-rank matrix to simplify the calculation process,reduce the complexity of the algorithm,and the denoising method has good denoising effect.

low-rank matrix; nuclear norm; CT image reconstruction

2016-04-08.

馬海英(1990-)，女，湖南湘潭人，廣西民族大學在讀碩士，研究方向：圖像處理與模式識別；宣士斌(1964-)，男，安徽無為人，博士，廣西民族大學教授，碩士生導師，研究方向：圖像處理與模式識別.

TP391.4

A

1673-8462(2016)03-0086-07