基于Harr小波的CS-MRI典型重構算法的性能分析

任筱倩，湯敏（南通大學電子信息學院，江蘇 南通，226019）

基于Harr小波的CS-MRI典型重構算法的性能分析

任筱倩，湯敏
（南通大學電子信息學院，江蘇 南通，226019）

目的 壓縮感知理論（Compressed Sensing，CS）與磁共振成像（Magnetic Resonance Imaging，MRI）相結合，縮短MRI圖像數(shù)據(jù)的掃描時間，提高成像質量。方法 以Harr小波進行稀疏表達，分別利用基追蹤（Basis Pursuit，BP）算法、正交匹配追蹤（Orthogonal Matching Pursuit，OMP）算法和分段正交匹配追蹤（Stagewise Orthogonal Matching Pursuit，StOMP）算法實現(xiàn)CS-MRI的二維重構。結果 在采樣率較低（10%-50%）時，以峰值信噪比（Peak Signal to Noise Ratio，PSNR）、平均結構相似度（Mean Structure Similarity，MSSIM）、相對誤差（Relative L2 Norm Error，RLNE）和傳輸邊緣信息（Transferred Edge Information，TEI）四個指標來定性、定量地評價和比較上述三種算法的重構質量，BP算法性能最佳。結論 BP算法能精確重構原始圖像，與完整采樣圖像相比，圖像質量并無明顯下降，同時大大減少MRI采集時間，具有重要的理論意義和臨床應用價值。

壓縮感知；重構算法；小波變換；Harr小波；MRI圖像

0　引言

傳統(tǒng)的信號獲取和處理過程分為信號采樣、變換、壓縮和重構這四個方面，其中信號獲取過程必須滿足香農(nóng)采樣定理，即信號的采樣頻率不得小于模擬信號頻譜中最高頻率的2倍，然后在信號壓縮過程中，進行某一特定變換后少數(shù)絕對值較大的數(shù)據(jù)進行壓縮，丟棄系數(shù)為零或接近零的數(shù)據(jù)。顯然，這種以香農(nóng)采樣定理為標準的高速采樣，然后再對采集獲得的信號進行壓縮的過程，浪費了大量的采樣資源，同時在數(shù)據(jù)傳輸和存儲中也導致巨大的資源浪費。針對上述傳統(tǒng)信號獲取和處理過程中存在的問題，Donoho、Candes、Romberg和Tao等人在大量理論研究的基礎上正式提出壓縮感知理論（Compressed Sensing，CS），可以通過對有限個采樣數(shù)據(jù)的重構，實現(xiàn)近似精確或者精確地恢復原始信號或圖像［1-3］。

磁共振成像（Magnetic Resonance Imaging，MRI），是根據(jù)磁共振現(xiàn)象從人體獲取電磁信號，并且由該信號重建出人體信息的一種斷層成像技術，已成為常用的無創(chuàng)醫(yī)學成像技術之一，在臨床醫(yī)學診斷中具有重要地位。但是，由于MRI掃描儀器所采集到的信號是經(jīng)過全局傅里葉變換后得到的頻域圖像，從而導致MRI成像時間過長，對于老人、小孩等無法長時間保持靜止的人群而言檢測效果不佳，對肺、胃、腸道等運動性器官而言，也可能導致圖像模糊。針對上述MRI的不足之處，Lustig提出CS-MRI成像理論，即從MRI的K空間數(shù)據(jù)的子集中準確重構出原始信號或圖像［4］。CS-MRI改變了醫(yī)學信息的獲取方式，將采樣速度提高到原來的幾十甚至幾百倍，使得掃描時間大為減少，重構圖像具有較高的空間分辨率，在臨床應用中備受關注。

本文研究基于Harr小波的CS-MRI圖像的壓縮感知二維重構，對典型算法進行性能分析和比較。算法流程如下：首先采用Harr小波變換實現(xiàn)MRI圖像的稀疏表達，然后分別利用基追蹤（Basis Pursuit，BP）算法、正交匹配追蹤（Orthogonal Matching Pursuit，OMP）算法和分段正交匹配追蹤（Stagewise Orthogonal Matching Pursuit，StOMP）算法實現(xiàn)CS-MRI圖像的二維重構，并以峰值信噪比（Peak Signal to Noise Ratio，PSNR）、平均結構相似度（Mean Structure Similarity，MSSIM）、相對誤差（Relative L2 Norm Error，RLNE）和傳輸邊緣信息（Transferred Edge Information，TEI）這四項指標來定性、定量地評價和比較上述三種重構算法的性能。

1　壓縮感知理論框架

壓縮感知理論的主要內(nèi)容是指：在某一變換域內(nèi)具有稀疏表示的信號，可以以遠低于奈奎斯特采樣定理標準的方式來采集數(shù)據(jù)，通過與變換基非相干的隨機投影采樣，運用合適的優(yōu)化算法高概率精確重構原始信號。該理論框架主要包含稀疏表示、測量矩陣和重構算法三個方面，詳見圖1［1-3，5］。

由圖1可見，原始信號在某些變換域中具有稀疏性，是壓縮感知精確重構原始信號的先驗條件；測量矩陣的選取與設計，直接影響隨機測量值是否保留了足夠的原始信號信息，能夠精確重構原始信號重構算法則是壓縮感知理論的核心，直接影響信號重構質量和重構速度。目前壓縮感知理論的重構算法主要分成貪婪算法、凸優(yōu)化算法、組合算法和統(tǒng)計優(yōu)化算法這四類。

（1）貪婪算法——基本思想是以當前情況為基礎作出最優(yōu)選擇，通過迭代方法逐步更新重構信號的支撐集，進而求解出最小l0范數(shù)問題的解。貪婪算法主要包括：子空間追蹤算法、匹配追蹤算法、正交匹配追蹤算法、分段正交匹配追蹤算法、正則化正交匹配追蹤算法等。根據(jù)選擇的支撐集元素是否采用回驗機制，還可以將這幾種重構算法劃分為兩類：非回驗算法，即對于已選入支撐集的元素不會被再次選入或者刪除；回驗算法，對已選中的某些元素可將其重新加入支撐集或不將其刪除［6］。

（2）凸優(yōu)化算法——基本思想是將非凸問題轉化成凸優(yōu)化問題來求最優(yōu)解，即將求解l0范數(shù)問題轉化為求解最小l1范數(shù)的線性問題，逐漸逼近原始信號以得到精確的重構結果。常用的凸優(yōu)化算法包括：基追蹤算法、內(nèi)點法、迭代閾值算法、加權迭代算法、梯度投影稀疏重建算法等。

（3）組合算法——首先對原始信號進行高度采樣，再對采樣結果進行群測試以獲得信號支撐。常用的組合算法包括：鏈追蹤算法、稀疏傅里葉描述算法、HHSP追蹤（Heavy Hitters On Steroids Pursuit，HHSP）算法等。

圖1　壓縮感知基本理論框架Fig.1 Theory Framework of Compressed Sensing

（4）統(tǒng)計優(yōu)化算法——利用典型信號的訓練集對信號的主要成分或者獨立成分進行分析，找出信號的最優(yōu)線性投影集合。目前的統(tǒng)計優(yōu)化方法主要分為基于訓練集合學習的統(tǒng)計優(yōu)化算法和基于貝葉斯統(tǒng)計理論的稀疏重構算法。

2　典型重構算法

2.1BP算法

Donoho等學者提出：在一定條件下，求解最小l0范數(shù)問題與求解最小l1范數(shù)問題是等價的。BP算法是對由測量向量組成的不同組合，采用線性優(yōu)化模型進行處理，以求得最優(yōu)解，獲得較佳的重構結果，但其計算復雜度較高。

BP算法的求解過程中，原始信觀測矩陣，基矩陣φ已知，求解s 等同于求解x ，其中s， x∈RN，y∈RM，M＜N。算法求解，滿足φx?=y ，其中。在含噪情況下，求解，其中n 為噪聲，σ≥0為噪聲強度，且，其約束條件可改為，滿足，對其進一步求解就可以重構出原始信號。其中參數(shù)σ表示算法允許誤差的范圍，其值為0時可認為是一般的無噪BP算法。

2.2OMP算法

Gilbert等人提出正交匹配追蹤算法的基本思想是：在信號稀疏度為的條件下，通過貪婪迭代算法選擇重構矩陣中的列，使得當前矩陣余量與每次迭代過程中選擇的列最大程度相關，從而去掉矩陣中的相關部分，然后對矩陣中的余量再次進行迭代直至迭代次數(shù)滿足殘差條件為止［7］。這種強制停止迭代過程的方法使得OMP算法需要較多的線性測量值來保證精確重構。OMP算法的具體步驟如下［7-8］：

輸入：采樣向量s ，傳感矩陣φ，稀疏度K；

輸出：x的K稀疏逼近x?；

初始化：索引集Λ0=φ，殘值r0=s，t=1；

循環(huán)執(zhí)行Step1～Step5：

Step 1：尋找殘值r 和傳感矩陣的列ψj積中最大元素的索引λ；

Step 2：更新索引集Λ0，記錄經(jīng)由傳感矩陣φ進行重構后的原子集合；

Step 3：進行最小二乘法；

Step 4：更新殘值，t=t+1；

Step 5：如Step 4中t＞K，則迭代停止；否則繼續(xù)執(zhí)行Step 1。

2.3StOMP算法

Donoho等學者在正交匹配追蹤算法的基礎上，采用分階段的思想，提出了分段正交匹配追蹤算法，首先從最初的原子集合中篩選出滿足閾值條件且與迭代余量相匹配的原子，然后對現(xiàn)有的支撐集和原子進行更新，最后利用最小二乘法求解得出近似精度，并且對迭代余量的更新完畢［9-10］。StOMP算法在分段處理階段中，首先構建一個殘差矢量序列r1，r2，……，然后建立一個逼近序列X0，X1，……Xs，Is標記Xs中非零元素的坐標值，初始解X0=0，Is為空，具體迭代步驟如下：

Step 1：初始化，設定最大迭代步長、最大迭代誤差e、s=1；

Step 2：將小波變換后獲得的稀疏矩陣記為y ，進行隨機測量；

Step 3：對Θ進行規(guī)范化處理，K為Θs信號長度；

Step 4：對Θ進行閾值處理；

Step 5：最新合并兩次變換所得的坐標索引，并對Is集合進行一致化處理；

Step 6：求解線性方程組min（||x||s 1）s.t.y=φXS ，計算公式如下；

Step 7：計算殘差rs；

Step 8：稀疏解的精度判定，如rS＜e，則X0=Xs；否則s=s+1，返回Step 3。

3　實驗結果與分析

3.1評價指標

醫(yī)學圖像的質量要求較高，需要算法能夠有效重構圖像，本文采用峰值信噪比PSNR、平均結構相似度MSSIM、相對誤差RLNE、傳輸邊緣信息TEI這四個指標共同定性、定量地評價上述算法的重構性能。

PSNR計算公式如下：

其中，f（ x， y）和f′（x， y）分別為完整采樣圖像和重構圖像，M×N為圖像空間分辨率。PSNR值越大，表示重構圖像與完整采樣圖像越接近，重構效果越好，重構精度越高。

結構相似度（Structural Similarity，SSIM），通過對完整采樣圖像和重構圖像進行對比，來判斷兩者之間的相似度，更符合人眼對圖像的視覺判斷，計算公式如下［11-12］：

其中，μx和σx是原始圖像的均值和標準差，μy和σy是重構圖像的均值和標準差，σxy是兩幅圖像的聯(lián)合方差，c1和c2是為了防止分母接近零時產(chǎn)生不穩(wěn)定現(xiàn)象而添加的固定常數(shù)。

平均結構相似度MSSIM是在SSIM基礎上發(fā)展而來的，首先利用大小相同的窗口對完整采樣圖像X 和重構圖像Y 進行分塊，確保各分塊之間不重疊，然后計算各塊的結構相似度，計算公式如下：

RLNE用來衡量兩幅圖像的結構不相符度，其值越小，表示重構圖像與完整采樣圖像之間的誤差越小，兩者更接近，計算公式如下［13-14］：

TEI用于衡量重構圖像與完整采樣圖像在保持圖像邊緣信息方面的差異，計算公式如下［15］：

其中，Qgxy和xy分別表示邊緣保持和方向保持的數(shù)據(jù)；G其值越大表示重構圖像邊緣信息保持越好。

3.2實驗結果分析

本實驗在處理器為Intel Core i3-2330M CPU@ 2.20GHz的機器上，采用MATLAB R2010b編程實現(xiàn)。首先以世界公認的Shepp Logan模型作為測試圖像，比較BP、OMP、StOMP算法在10%-50%欠采樣率下的重構效果，性能指標參見表1。同時，為從人眼視覺系統(tǒng)直觀感受重構圖像質量，圖2-圖4給出10%采樣率下的重構圖像、誤差圖像和局部放大圖像。

表1　Shepp Logan模型在不同欠采樣率下的重構性能指標Table 1 Performance Comparison of BP， OMP， StOMP Algorithms for Shepp Logan Phantom

由表1和圖2—圖4（第20頁），可以得出以下結論：（1）隨著欠采樣率的增大，BP算法和OMP算法重構圖像的PSNR逐步增大，TEI和MSSIM越來越大并接近1，RLNE越來越小，這顯然與壓縮感知理論完全符合；（2）StOMP重構算法重構出的圖像呈現(xiàn)PSNR、TEI和MSSIM越來越小，而RLNE卻越來越大的趨勢，這是因為StOMP算法在迭代過程中不一定能夠找到信號的最佳表示方法，使得重構精度大大降低，且稀疏度選取也不是很合適；（3）在相同欠采樣率條件下，OMP算法重構圖像的PSNR小于BP算法，RLNE指標略高于BP算法，這說明在相同的欠采樣率下，BP算法對于Shepp Logan模型的重構圖像質量明顯優(yōu)于OMP算法；（4）在相同采樣率情況下，BP重構算法運行時間最長，但圖像重構質量最佳，重構誤差最小，重構圖像與原圖更接近。

圖2　采樣率10%時Shepp Logan 模型的重構圖像Fig.2 CS-MRI Recovery Results of Shepp Logan Phantom at 10% Sampling Ratio

圖3　采樣率10%時Shepp Logan模型的重構誤差圖像Fig.3 Recovery Error Images of Shepp Logan Phantom at 10% Sampling Ratio

圖4　采樣率10%時Shepp Logan模型的局部放大圖像Fig.4 Local Magnification Images of Shepp Logan Phantom at 10% Sampling Ratio

由于Shepp Logan圖像是嚴格數(shù)學解析的，因此圖像重構誤差均由算法產(chǎn)生，上述實驗充分說明BP、OMP算法均能很好地實現(xiàn)欠采樣條件下的壓縮感知二維重構。下面選取圖5所示的真實MRI圖像，進一步驗證上述算法的有效性和實用性。表2-表4分別是圖5所示三幅MRI圖像在不同欠采樣率下圖像重構的性能指標。

圖5　真實MRI圖像Fig.5 Truly MRI Images

從表2-表4可以看出，對于真實MRI圖像，在10%-50%欠采樣率下，BP算法總能獲得最佳的性能指標。以MRI-Brain.jpg圖像為例，圖6-圖8是其在25%欠采樣率下的重構圖像、誤差圖像和局部放大圖像。

表2　MRI.jpg圖像不同采樣率下的重構性能指標Table 2 Performance Comparison of BP， OMP， StOMP Algorithms for MRI.jpg

表3　Heart.png圖像不同采樣率下的重構性能指標Table 3 Performance Comparison of BP， OMP， StOMP Algorithms for Heart.png

表4　MRI-Brain.jpg圖像不同采樣率下的重構性能指標Table 4 Performance Comparison of BP， OMP， StOMP Algorithms for MRI-Brain.jpg

圖6　采樣率25%時的MRI-Brain.jpg圖像重構結果Fig.6 CS-MRI Recovery Results of MRI-Brain.jpg at 25% Sampling Ratio

圖7　采樣率25%時的MRI-Brain.jpg圖像重構誤差圖像Fig.7 Recovery Error Images of MRI-Brain.jpg at 25% Sampling Ratio

圖8　采樣率25%的MRI-Brain.jpg局部放大圖像Fig.8 Local Magnification Images of MRI-Brain.jpg at 25% Sampling Ratio

由表2-表4和圖6-圖8可知，對于真實MRI圖像而言，（1）隨著欠采樣率增大，BP算法和OMP算法重構出的圖像PSNR、TEI和MSSIM越來越大，RLNE越來越小，同時，BP算法的各項性能指標均最優(yōu)；（2）欠采樣率≥30%以后，BP算法重構圖像與完整采樣圖像之間的誤差基本＜10%，可應用于臨床真實MRI圖像的壓縮感知重構。

4　討論與結論

伴隨信息技術的蓬勃發(fā)展，生活中所需處理的信息量越來越大，使得信息在采集、存儲及傳輸?shù)确矫嫠媾R的壓力也越來越大，傳統(tǒng)信號的采集與壓縮過程在實際運用過程具有許多不適應的地方，壓縮感知理論則較好地解決了目前面臨的采樣速度慢、存儲及傳輸量巨大的問題。

本文在簡要闡述壓縮感知理論框架和重構算法原理的基礎上，利用BP、OMP和StOMP算法對嚴格解析的Shepp Logan模型和真實MRI圖像進行了不同欠采樣率的壓縮感知二維重構，結論如下：（1）在相同欠采樣率的情況下，BP算法的重構結果最佳，雖運算時間最長，但仍可滿足實時性的要求；StOMP算法重構結果最差，不適用于臨床應用；（2）在采樣率較低的情況下，利用壓縮感知理論可以近似精確地重構出原始MRI圖像，解決了傳統(tǒng)采樣時間長和采集后數(shù)據(jù)大量浪費的問題，縮短采樣時間，提高圖像質量，這對于MRI臨床應用具有重要的理論價值和實用前景。

今后的研究方向包括：（1）在圖形處理器GPU上實現(xiàn)算法的硬件加速，從而在保證BP算法重構誤差最小，重構精度最高的同時，盡量縮短其運行時間；（2）將本文算法應用于更多不同部位、不同疾病的CS-MRI圖像重構中，檢驗算法的普遍適用性。

［1］ David Donoho. Compressed Sensing ［J］. IEEE Transactions on Information Theory， 2006， 52（4）: 1289-1306.

［2］ Emmanuel J. Candes， Justin K Romberg， Terence Tao. Stable Signal Recovery from Incomplete and Inaccurate Measurements ［J］. Communications on Pure and Applied Mathematics， 2006， 59（8）: 1207-1223.

［3］ Emmanuel J. Candes， Justin K Romberg， Terence Tao. Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Fourier Information ［J］. IEEE Transactions on Information Theory， 2006， 52（2）: 489-509.

［4］ Michael Lustig， David Donoho， John Pauly. Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging ［J］. Magnetic Resonance in Medicine， 2007， 58（6）: 1182-1195.

［5］ Yonina Eldar， Gitta Kutyniok. Compressed Sensing: Theory and Applications ［M］. Cambridge University Press， 2012

［6］ 王超. 基于壓縮感知的貪婪迭代重構算法［J］. 數(shù)據(jù)采集與處理，2012，27（S2）：298-303.

［7］ Tropp J， Gilbert A. Signal Recovery from Random Measurements via Orthogonal Matching Pursuit ［J］. IEEE Transactions on Information Theory， 2007， 53（12）: 4655-4666.

［8］ 楊真真，楊震，孫林慧. 信號壓縮重構的正交匹配追蹤類算法綜述［J］. 信號處理，2013，29（4）：486-496.

［9］ Maleki A， Donoho D L. Optimally Tuned Iterative Reconstruction Algorithms for Compressed Sensing ［J］. IEEE Journal of Selected Topics in Signal Processing， 2010， 4（2）: 330-341.

［10］ 解成俊，張鐵山. 基于壓縮感知理論的圖像重構算法研究［J］. 計算機應用與軟件，2012，29（4）: 49-52.

［11］ Zhou Wang， Alan C. Bovik， Hamid R. Sheikh， Eero P. Simoncelli. Image Quality Assessment: From Error Visibility to Structural Similarity ［J］. IEEE Transaction on Image Processing， 2004， 13（4）: 600-612.

［12］ Xiao Chen， Michael Salerno， Yang Yang， Frederick H. Epstein. Motion-Compensated Compressed Sensing for Dynamic Contrast-Enhanced MRI Using Regional Spatiotemporal Sparsity and Region Tracking: Block LOw-rank Sparsity with Motion-guidance （BLOSM）［J］. Magnetic Resonance in Medicine， 2014， 72（4）: 1028-1038.

［13］ Xiaobo Qu， Yingkun Hou， Fan Lam， Di Guo， Jianhui Zhong， Zhong Chen. Magnetic Resonance Image Reconstruction From Undersampled Measurements Using a Patch-Based Nonlocal Operator ［J］. Medical Image Analysis， 2014， 18（6）: 843-856.

［14］ Yunsong Liu， Zhifang Zhan， Jian-Feng Cai， Di Guo， Zhong Chen， Xiaobo Qu. A Simple and Fast Iterative Soft-thresholding Algorithm for Tight Frames in Compressed Sensing Magnetic Resonance Imaging ［J］. arXiv preprint arXiv:1504.07786， 2015.

［15］ C.S. Xydeas， V. Petrovic. Objective Image Fusion Performance Measure ［J］. Electronics Letters， 2000， 36（4）: 308-309.

Performance Analysis of CS-MRI Recovery Algorithms Based on Harr Wavelet

REN Xiao-Qian， TANG Min
（School of Electronics and Information Engineering， Nantong University， Nantong 226007， China）

Objective The CS （Compressed Sensing） theory is introduced into MRI （magnetic resonance imaging） to solve the shortcomings of the long acquisition time of MRI， and to change the way of medical imaging techniques. Methods Three typical CS-MRI recovery algorithms， including basis pursuit （BP）， orthogonal matching pursuit （OMP） and stagewise orthogonal matching pursuit （StOMP） are analyzed and compared based on the Harr wavelet. Results Four performance indices-peak signal to noise ratio （PSNR）， mean structural similarity （MSSIM）， relative l2 norm error （RLNE） and transferred edge information （TEI） are applied to evaluate the recovery image’s quality under different sampling ratios. It is clear that the BP algorithm achieves the best performance. Conclusion The original reference images can be accurately recovered based on BP algorithm under the low sampling ratios （10%-50%）， which has important theoretical significance and clinical application foreground in MRI technique.

Compressed sensing； Reconstruction algorithm； Wavelet transform； Harr wavelet； MRI images

國家自然科學基金資助項目（61401239）；江蘇省自然科學基金資助項目（BK20130393）；江蘇高校品牌專業(yè)建設工程資助項目（PPZY2015B135）

湯敏，女，漢族，江蘇南通人，副教授，碩士生導師，從事醫(yī)學圖像處理與分析等相關研究。