一種改進的變步長OMP圖像重建算法

吳小龍 伍松

摘? ? 要：為了快速、高精度的重建圖像，解決濾波反投影（FBP）算法重建圖像精度不高，正交匹配追蹤（OMP）算法運行時間較長的問題，基于改變步長，提出一種步長變換正交匹配追蹤（SCOMP）算法.當殘差不小于閾值時，增大步長進行運算，當殘差小于閾值時，恢復原步長進行運算.研究結(jié)果表明：SCOMP算法重建圖像精度高于OMP算法，且運行時間快于FBP算法.SCOMP算法采用大步長快速添加原子，小步長有效去除原子的方法，使得重建圖像的精度較高且運行時間也較短.

關鍵詞：重建;反投影;正交匹配;步長變換

中圖分類號：TN911.73? ? ? ? ? DOI：10.16375/j.cnki.cn45-1395/t.2019.04.011

0? ? 引言

圖像重建廣泛出現(xiàn)在醫(yī)學掃描，食品檢測，加工裝配等諸多領域，所以對于圖像的重建在生產(chǎn)生活中非常重要.

濾波反投影（Filtered Back Projection，F(xiàn)BP）成像屬于工業(yè)計算機斷層成像（Industrial Computed Tomography，CT）技術.CT技術是一種由外到內(nèi)的檢測技術[1].FBP算法具有運算速度快的優(yōu)點.國外關于濾波反投影的研究包括：Katsevich等[2]研究了不完全投影的濾波反投影重建算法.Pelt等[3]設計了一種關于數(shù)據(jù)的濾波器來減少投影產(chǎn)生的誤差.國內(nèi)西安交通大學，上海交通大學等也開始對該方面內(nèi)容進行了研究.

正交匹配追蹤（Orthogonal Matching Pursuit，OMP）算法屬于壓縮感知（Compressed Sensing，CS），壓縮感知理論首先由Candes等[4]提出.在國外，麻省理工學院、萊斯大學等一些知名大學已經(jīng)成立了專門的課題組對此進行研究.國內(nèi)很多研究單位的學者也已經(jīng)開始對壓縮感知進行研究，如清華大學、西安電子科技大學也都專門成立了課題組.

為快速、高精度的重建出圖像，本文提出一種基于改變步長的步長變換正交匹配追蹤（Step Change Orthogonal Matching Pursuit，SCOMP）算法，該算法運行時間較短，重建圖像的精度較高.

1? ? 壓縮感知理論與OMP算法

1.1? ?壓縮感知的基本原理

Donoho[5]提出并擴展了壓縮傳感理論，他用大量的實驗證明了壓縮傳感在信號處理方面具有廣闊的前景.根據(jù)壓縮感知理論可得知，運用適當?shù)闹亟ㄋ惴軌驈倪@些數(shù)據(jù)中高概率的恢復原始信號.壓縮感知理論核心的思想是：將壓縮與采樣過程合并在一起完成，從而可以顯著的減少采樣點數(shù)，節(jié)省儲存空間[6].

1.2? ?信號的稀疏表示

由文獻[7]可知任意信號X可表示成下式：

[X=ΨΘ]? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（1）

式中[Θ]是投影系數(shù)，其維數(shù)為N×1的列向量，Ψ為變換基.

1.3? ?測量矩陣

由文獻[7]可知，測量信號[Y] 可表示成下式：

[Y=ФΘ=ФΨTX=ACSX]? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （2）

式中[Ф]是測量矩陣，[ACS]是傳感矩陣.

1.4? ?優(yōu)化重構(gòu)

首先，若信號X[∈]RN在某一個變換基[Ψ]上是能夠稀疏表示的，那么變換系數(shù)可以寫作[Θ=ΨTX];然后，需要設計一個穩(wěn)定的、與變換基[Ψ]不相關的測量矩陣[Ф]，對[Θ]進行映射，映射結(jié)果表示為[Y=ACS]（其中[ACS=ФΨT]）;最后，求解一個范數(shù)的優(yōu)化問題[8]，從而解得原始信號X的精確解或者是近似解.重建過程如下式：

[minΘ0s.t.? ACSX=ФΨTX=Y]? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（3）

優(yōu)化求得原信號X在變換基上的最稀疏表示[Θ]，然后作逆變換就可求得原始信號X.

1.5? ?OMP算法

正交匹配追蹤算法是壓縮感知中最為常用的一種算法[9].其核心思想是，在迭代過程中，要從傳感矩陣[ACS]選出與觀測信號Y相關度（內(nèi)積）最大的那一列，然后從[ACS]中去掉該列并加入到擴充矩陣T中，接著求得殘差r_n最小的一個估值aug_y，最后重復運行，達到迭代次數(shù)s為止.

OMP算法的基本核心步驟如下：

輸入：觀測信號Y，傳感矩陣[ACS];

輸出：原信號的稀疏逼近值theta;

初始化：殘差r_n=Y（：，t0），儲存集T為空集，Y的行數(shù)t0=0，迭代次數(shù)s=0;

Step1 找到傳感矩陣[ACS]與殘差r_n最相關的列[index_ACS]（s）=max_index;

Step2 更新索引集T=[T，[ACS]（：，[index_ACS] （s））];

Step3 使殘差最小aug_y=pinv（T）* Y（：，t0），更新殘差r_n=Y（：，t0）-T*aug_y;

Step4 直至迭代次數(shù)結(jié)束.

2? ? FBP算法

2.1? ?傅里葉中心切片定理

假設f（x，y）為待重構(gòu)物體的密度函數(shù)，[p?]（xr）為? f（x，y）在角度?=?0時的平行束投影.該定理的表達式為：

[F1[p?（xr）]=F（ρ，?）|?-?0]? ? ? ?（4）

其中：F1[ ]是一維變換，F(xiàn)（ρ，?）是二維極坐標表示.

傅里葉中心切片定理[10-11]提供了頻域上更簡單的數(shù)學關系，即對物體的投影進行一維變換，如圖1所示[12].

2.2? ?濾波反投影算法原理

濾波反投影重建算法是最常用的CT重建算法[13-14].濾波反投影重建公式如下：

[ftr， θ=fx， y=02πp（xr， ?）*h（xr）d?=02πg（xr， ?）d?]? ? ? ? ? ? ? ?（5）

式中，[ft（r，θ）]為[f（x，y）]沿直線的線積分，[h（xr）=F-11[ρ]]，[p（xr ， ?）=F-11[p（ρ， ?）]].

FBP算法基本核心步驟如下：

Step1 設計合適的濾波器[h];

Step2 把投影數(shù)據(jù)與濾波器進行卷積運算，得到投影數(shù)據(jù)[g（xr ， ?i）];

Step3 對于每一個角度[?i]，把濾波后的投影[g（xr ， ?i）]反投影于滿足[xr]，[rcos（θ-?i）]的射線上的所有點[（r， θ）];

Step4 反投影數(shù)值進行累加，得到重建后的圖像數(shù)據(jù).

3? ? SCOMP算法

3.1? ?連續(xù)小波變換

函數(shù)f （t）連續(xù)小波變換表達式如下：

[ Wfa，b=ft ，? ψa，bt=1aRftψ*（t-ba）dt]? ? ? ? ? ? ? ? ? ? （6）

小波變換的逆變換為：

[ft=1Cψ-∞∞-∞∞1a2Wfa，? bψ（t-ba）dadb]? ? ? ? ? ? ? ? ? ? ? ?（7）

其中a與小波大小有關，b與位置有關，小波函數(shù)[ψ]（t）的傅里葉變換為：

[ψω=-∞∞ψ（t）e-jωtdt]? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （8）

當[Cψ=-∞∞ψ（ω）2ωdω<∞]時滿足小波完全重構(gòu)的條件，此時[ψ]（t）經(jīng)過伸縮和平移得到波基函數(shù)：

[ψa，bt=1aψa，bt-ba ，? a，? b∈R， a≠0]? ? ? ? ? ? ? ? ? ? ? （9）

3.2? ? DWT離散小波變換

在特征值提取時要采用連續(xù)的小波變換在每個可能的尺寸去計算小波系數(shù)，這樣會產(chǎn)生大量的冗余數(shù)據(jù)，因此在實際的應用中連續(xù)小波必須加以離散化，這一離散化是針對尺度參數(shù)a和平移參數(shù)b的，而不是針對時間變量t的.

對a和b離散化，公式為：[a= aj0]，[b= aj0b0]，（j，k∈Z），a0為擴展步長，且a0≠1是固定值，對應的離散小波函數(shù)[ψj，kt]可寫成：

[ψj，kt=a-j20ψt-kaj0b0aj0=a-j20ψ（a-j0t-kb0）]? ? ? ? ? ? ? ? ? ? （10）

[Cj，k=-∞∞ψ*j，k（t）dt=f， ψj，k]? ? ? ? ? ? ? ? ? ? ? ? ? ? （11）

式（11）為離散化小波變換的系數(shù)的求解公式.

3.3? ?SCOMP算法原理

OMP算法的步長是固定值1，這樣就會導致計算的時間變長，從而影響效率.因此提出一種步長變換正交匹配追蹤重建算法，當殘差的長度大于閾值時，增大步長.即殘差的長度若不小于閾值r2（norm（Res）≥r2），則步長會變大，以縮短運行時間，提高重建精度.

SCOMP算法基本核心步驟如下：

輸入：M×N維測量矩陣Phi，觀測信號y;

定義：前向步長α，后向步長β，閥值r1=eps*norm（y），r2=5×r1;

初始化：K=5α，eps=10-8，索引集T=[ ]，迭代次數(shù)n=0，殘差Res=y，α=10，β=2;

Step1 計算殘差Res和標準化測量矩陣Phinorm每一列的內(nèi)積C=（Phinorm’*Res）;

Step2 將更新的α個原子放入索引集中T=（T;ind（1∶ α））并去除β個最小的原子T（ind（1∶ β））=[ ];

Step3 更新殘差，若此時norm（Res）≥r2，則α變?yōu)?0，β變?yōu)?，否則還是按照α=10，β=2運行;

Step4 如果T的長度大于N或者殘差的長度小于r1，則結(jié)束運行，輸出X，否則返回Step2.

對于前向步長與后向步長的取值，表1給出了實驗數(shù)據(jù)：

可以看出，初始前后步長的差值與改變后前后步長的差值若是較小，會使運算時間變長且重構(gòu)圖像的精度變低;初始前后步長的差值與改變后前后步長的差值若是較大，此時運算時間較短但重構(gòu)圖像的精度會變低.

由此可得，要使運行時間較低且重建圖像精度較高，則前后步長差值要適當，過長會導致重建過程中不能去除較差的原子使重建圖像精度低，過短會使算法運行時間變長.且前后步長具體取值過大會導致迭代迅速結(jié)束，圖像精度不高，過短會使算法運行時間變長.經(jīng)試驗，取初始前向步長α=10，后向步長? ?β=2，改變后的前向步長α=20，后向步長β=4，此時可以保證算法運行時間較短，重建圖像精度較高.

對于給定值eps的取值，表2給出了實驗數(shù)據(jù).

eps的取值直接關系r1的取值，eps取值較小，r1較小，對迭代的精度要求高，則算法迭代的次數(shù)會比較多，可能會影響運算時間，eps取值較大，r1較大，則算法迭代的次數(shù)比較少，不能保證重構(gòu)圖像的精度.經(jīng)試驗，取eps為10-8，此時可以保證算法運行時間較短，圖像重建精度較高.

對于r2的取值，表3給出了實驗數(shù)據(jù).

r2取值較大，則儲存集保存的較大原子就多，但會增加運行時間，取值較小，儲存集保存的較大原子就小.經(jīng)試驗，選取r2的值為5倍的r1，此時可以保證算法運行時間較短，圖像重建精度較高.

綜上所述，該算法中，對運行時間與重建精度影響較大的是前后步長的選取，因為前步長過大會直接影響運行時間，前步長過小則添加的原子數(shù)不夠，遺失了較大的原子，影響重建精度;后步長過大會直接去除掉較大的原子，即增加運行時間，又影響精度，后步長過小則可能保留較小的原子，索引集的長度會變大，后面較大的原子可能不會被選入索引集中，影響重建精度.

4? ? 實驗結(jié)果

選取大小為256×256的Lena圖像進行重建，軟件為MatlabR2014b，3種算法重建的圖像如圖3—圖5所示.其中OMP算法和SCOMP算法均用DWT離散小波進行信號稀疏，圖像壓縮比均為0.5，所用測量矩陣均為高斯隨機矩陣，F(xiàn)BP算法測量角度為180°.各算法運算時間和峰值信噪比如表4所示.

圖2為大小256×256的Lena原圖.

圖3為FBP算法重建圖像，圖像整體較為模糊，峰值信噪比較低.

圖4為OMP算法重建圖像，圖像較為清晰，峰值信噪比較高，但從表中可以看出該算法運行時間比較長.

圖5為SCOMP算法重建圖像，取初始前步長α=10，后步長β=2，改變后的前步長α=20，后步長β=4，取定值eps=10-8，取閾值r2=5×r1，此時圖像清晰度稍高于OMP算法重建圖像的清晰度，且該算法的運行時間大大低于OMP算法的運行時間.

SCOMP算法采用大步長快速添加原子，小步長有效去除原子的方法，經(jīng)過試驗，可以使重建圖像精度較高且運行時間較短，解決了FBP算法重建圖像精度不高，OMP算法運行時間較長的問題.

5? ? 結(jié)果分析

在對圖像進行重建時，要平衡好運算時間與重建精度的關系，當通過某種方法使運算時間變短時要充分考慮到此改變對于精度的影響，并通過大量實驗去驗證.反之，當通過某種方法使精度提高時，要考慮該種方法是否會使運算變的冗長從而增加了運算的時間.在某一個因素改變時要對它改變的范圍進行確定，充分了解該因素改變會對運算過程中哪些值產(chǎn)生影響，并思考如何去抑制不良的影響，也要考慮改變的值相互之間是否有影響.有時一個因素的改變會對結(jié)果產(chǎn)生較差的影響，此時可以對兩個或者多個因素進行改變.此外，還要注意對主要因素的控制，因為該因素會對結(jié)果產(chǎn)生最大的影響，要將該因素的值控制在一個合理的范圍，再對其他因素進行改變尋找最合適的取值.SCOMP算法較好的平衡了運算時間與重建精度的關系.

6? ? 結(jié)束語

未來可從以下幾方面進行研究：一是將步長變換正交匹配追蹤算法與濾波反投影算法相結(jié)合進行圖像重建;二是提高步長變換正交匹配追蹤算法重建圖像的精度;三是可以通過除殘差長度外的其他因素來判斷是否要進行步長改變;四是可以進行三維的重建.

參考文獻

[1]? ? 張朝宗，郭志平，張朋，等.工業(yè)CT技術和原理[M].北京：科學出版社，2009.

[2]? ? KATSEVICH A，RAMM A. Filtered back projection method for inversion of incomplete tomographic data[J].Applied Mathematics Letters，1992，5：77-80.

[3]? ? PELT D M，BATENBURG K J. Improving filtered back-projection reconstruction by data-dependent filtering[J].IEEE Transactions on Image Processing，2014，23（11）：4750-4762.

[4]? ? CANDES E，TAO T.Decoding by linear programming[J].IEEE Transactions on Information Theory，2005，51（12）：4203-4215.

[5]? ? DONOHO D.Compressed sensing[J].IEEE Transations on Information Theory，2006，52（4）：1289-1306.

[6]? ? 陶佳偉，李春貴.基于壓縮感知的肝臟CT/PET醫(yī)學圖像融合研究[J].廣西科技大學學報，2015，26（3）：32-33.

[7]? ? 杜寶.基于壓縮感知的平面近場聲全息理論與實驗研究[D].昆明：昆明理工大學，2017.

[8]? ? 汪霜霜，李春貴.一種lp正則化改進的車輛軌跡學習算法[J].廣西科技大學學報，2019，30（2）：58-59.

[9]? ? SHEN Y，LI S. Sparse signals recovery from noisy measurements by orthogonal matching pursuit[J].Inverse Problems & Imaging，2015，9（1）：231-238.

[10]? RONALD N BRACEWELL. The fourier transform and its applications[M].西安：西安交通大學出版社，2005.

[11]? REN N G. Fourier slice photography[J].ACM Transaction on Graphics，2005，24（3）：735-744.

[12]? 莊天戈.CT理論與算法[M].上海：上海交通大學出版社，1992.

[13]? 范慧赟.CT圖像濾波反投影重建算法的研究[D].西安：西北工業(yè)大學，2007.

[14]? 余曉鍔，龔劍，馬建華.CT原理與技術[M].北京：科學出版社，2013.

An improved OMP image reconstruction algorithm based on

changing step size

WU Xiaolong1，2， WU Song*1，2

（1.School of Mechanical and Traffic Engineering， Guangxi University of Science and Technology， Liuzhou 545006， China; 2.Guangxi Key Laboratory of Automotive Components and Vehicle Technology（Guangxi

University of Science and Technology）， Liuzhou 545006， China）

Abstract： A step change orthogonal matching pursuit（SCOMP） algorithm is proposed to reconstruct? image with high speed and precision as the resolution of the filtered back projection（FBP） algorithm is not high and the orthogonal matching pursuit（OMP） algorithm has a long running time. When the? ? ? ?residual is not less than the threshold， the step size is increased. When the residual is less than the threshold， the original step is restored. The results show that the accuracy of SCOMP reconstruction