吳小龍 伍松

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

關(guān)鍵詞：重建;反投影;正交匹配;步長變換

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

0? ? 引言

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

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

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

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

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

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

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

1.2? ?信號(hào)的稀疏表示

由文獻(xiàn)[7]可知任意信號(hào)X可表示成下式：

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

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

1.3? ?測量矩陣

由文獻(xiàn)[7]可知，測量信號(hào)[Y] 可表示成下式：

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

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

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

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

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

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

1.5? ?OMP算法

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

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

輸入：觀測信號(hào)Y，傳感矩陣[ACS];

輸出：原信號(hào)的稀疏逼近值theta;

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

Step1 找到傳感矩陣[ACS]與殘差r_n最相關(guā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? ?傅里葉中心切片定理

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

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

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

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

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

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

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

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

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

Step1 設(shè)計(jì)合適的濾波器[h];

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

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

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

3? ? SCOMP算法

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

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

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

小波變換的逆變換為：

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

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

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

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

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

3.2? ? DWT離散小波變換

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

對(duì)a和b離散化，公式為：[a= aj0]，[b= aj0b0]，（j，k∈Z），a0為擴(kuò)展步長，且a0≠1是固定值，對(duì)應(yīng)的離散小波函數(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，這樣就會(huì)導(dǎo)致計(jì)算的時(shí)間變長，從而影響效率.因此提出一種步長變換正交匹配追蹤重建算法，當(dāng)殘差的長度大于閾值時(shí)，增大步長.即殘差的長度若不小于閾值r2（norm（Res）≥r2），則步長會(huì)變大，以縮短運(yùn)行時(shí)間，提高重建精度.

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

輸入：M×N維測量矩陣Phi，觀測信號(hào)y;

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

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

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

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

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

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

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

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

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

對(duì)于給定值eps的取值，表2給出了實(shí)驗(yàn)數(shù)據(jù).

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

對(duì)于r2的取值，表3給出了實(shí)驗(yàn)數(shù)據(jù).

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

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

4? ? 實(shí)驗(yàn)結(jié)果

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

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

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

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

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

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

5? ? 結(jié)果分析

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

6? ? 結(jié)束語

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

參考文獻(xiàn)

[1]? ? 張朝宗，郭志平，張朋，等.工業(yè)CT技術(shù)和原理[M].北京：科學(xué)出版社，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ī)學(xué)圖像融合研究[J].廣西科技大學(xué)學(xué)報(bào)，2015，26（3）：32-33.

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

[8]? ? 汪霜霜，李春貴.一種lp正則化改進(jìn)的車輛軌跡學(xué)習(xí)算法[J].廣西科技大學(xué)學(xué)報(bào)，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].西安：西安交通大學(xué)出版社，2005.

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

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

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

[14]? 余曉鍔，龔劍，馬建華.CT原理與技術(shù)[M].北京：科學(xué)出版社，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