王彩芳 殷緒導(dǎo) 張曉媛 徐兆亮

摘要：

將次梯度投影迭代算法應(yīng)用到數(shù)字圖像重建問題。將圖像重建問題轉(zhuǎn)化為求一個(gè)加權(quán)最小二乘問題，導(dǎo)出次梯度投影算子在該問題下的具體迭代形式，并采用并行計(jì)算策略重建算法。通過三維數(shù)值實(shí)驗(yàn)對比次梯度投影迭代算法與常用的SART算法，驗(yàn)證算法的可行性和效率。

關(guān)鍵詞：

圖像重建； 加權(quán)最小二乘； 次梯度投影算子； 并行計(jì)算

中圖分類號(hào)： O244；O29

文獻(xiàn)標(biāo)志碼： A

Parallel strategy of subgradient projection operator

for image reconstruction

WANG Caifang， YIN Xudao， ZHANG Xiaoyuan， XU Zhaoliang

（College of Arts and Sciences， Shanghai Maritime University， Shanghai 201306， China）

Abstract：

The subgradient projection iterative algorithm is applied to digital image reconstruction. The problem of image reconstruction is transformed into a weighted least-square problem， the specific iterative form of the subgradient projection operator is derived， and parallel computing strategy reconstruction algorithm is adopted. The iterative algorithm of subgradient projection and the common SART algorithm are compared by three-dimensional numerical experiments. The feasibility and efficiency of the algorithm are verified.

Key words：

image reconstruction； weighted least-square； subgradient projection operator； parallel computation

收稿日期： 2018-01-26

修回日期： 2018-02-05

基金項(xiàng)目：

國家自然科學(xué)基金（11401372）

作者簡介：

王彩芳（1981—），女，江蘇蘇州人，副教授，博士，研究方向?yàn)閿?shù)字圖像重建算法，（E-mail）cfwang@shmtu.edu.cn

0 引 言

數(shù)字圖像重建問題是根據(jù)待測物體邊界測量數(shù)據(jù)獲取待測物體內(nèi)部結(jié)構(gòu)的一類反問題。[1]近年來，重建問題的研究主要集中在如何找到一種合適的掃描方式和快速的成像算法。以三維CT為例，基于FDK公式[2]的近似重建算法和基于Katsevich公式[3]的精確重建算法已經(jīng)在成像領(lǐng)域得到廣泛應(yīng)用。這類分析重建算法雖然成像速度快，但基本在射線變換和三維Radon變換的框架下推演，對數(shù)據(jù)的完整性要求高。當(dāng)部分投影數(shù)據(jù)缺失或獲取存在困難時(shí)，分析算法設(shè)計(jì)極為復(fù)雜，因此代數(shù)算法漸漸進(jìn)入研究者的視野。隨著計(jì)算機(jī)技術(shù)的高度發(fā)展，代數(shù)重建算法耗時(shí)長的缺點(diǎn)也有所改善。

為利用代數(shù)重建算法，很多成像問題都用線性方程組進(jìn)行描述。

Ax=

b（1）

式中：

A∈Rm×n為系數(shù)矩陣，稱為成像矩陣；未知數(shù)x∈Rn表示待測物體對入射光的吸收系數(shù)或其他光學(xué)系數(shù)；右端項(xiàng)b∈Rm為邊界測量數(shù)據(jù)或修正后的邊界測量數(shù)據(jù)。圖像重建問題即為求解光學(xué)系數(shù)x的反問題。在圖像重建問題中，系數(shù)矩陣A規(guī)模很大，且Ax=b是病態(tài)問題。EM迭代算法和Landweber 迭代算法是圖像重建問題中應(yīng)用廣泛的兩類算法。EM算法主要用于求解觀測數(shù)據(jù)中包含的隱藏參數(shù)的極大似然估計(jì)或極大先驗(yàn)估計(jì)。[4]在成像問題中，其本質(zhì)為尋找使得觀測數(shù)據(jù)與計(jì)算數(shù)據(jù)之間 Kullback-Leibler（KL）距離最小的解。[5]Landweber 迭代算法包括一系列著名算法：SART算法[6]、CIMMINO方法[7]、CAV算法[8]、DWE算法[9]和DROP方法[10]等。其本質(zhì)是尋找上述線性系統(tǒng)的加權(quán)最小二乘解[5]。Landweber算法的收斂性受迭代過程中松弛因子的影響。次梯度投影迭代算法是求解凸可行性問題的一類重要算法[11]，因而被廣泛應(yīng)用到實(shí)際問題的求解中[12-13]。與Landweber算法不同，次梯度投影算法在迭代過程中無需預(yù)設(shè)松弛因子，迭代步長完全由目標(biāo)函數(shù)和當(dāng)前估計(jì)值所確定。因此，考慮將次梯度迭代算法應(yīng)用到線性圖像重建問題中，研究該算法的快速重建方法。

本文回顧次梯度投影算法，并與重建問題相結(jié)合，導(dǎo)出具體的迭代格式，估算迭代步長，研究迭代算法的并行策略，得到具體算法實(shí)現(xiàn)過程。利用三維數(shù)值實(shí)驗(yàn)驗(yàn)證算法的可行性，并將重建結(jié)果與SART算法進(jìn)行比較。

1 重建算法導(dǎo)出

假設(shè)f：Rm→R為凸且可微的函數(shù)，將f的δ-水平集定義為levδ f={

x∈Rm|f

（x）≤δ}，若滿足levδ f≠條件，則f定義的次梯度投影算子[11-14]可描述為

Gf（x）=

x-f（x）-δ‖s（x）‖2s（x），f（x）>δ

x，f（x）≤δ（2）

式中：s（x）∈f（x）為函數(shù)f（x）在點(diǎn)x處的次梯度；‖·‖為Rm上的2-范數(shù)。根據(jù)次梯度投影算子的性質(zhì)，其不動(dòng)點(diǎn)集合滿足Fix（Gf）={x∈Rm|f（x）≤δ}。

針對式（1）病態(tài)問題，應(yīng)考慮其對應(yīng)的加權(quán)最小二乘問題，并導(dǎo)出次梯度投影算子的具體表達(dá)形式，計(jì)算原理如下。

1.1 考慮加權(quán)最小二乘問題

LW（x）=12‖Ax-b‖2W（3）

式中：W∈Rm×m為對稱正定矩陣。該問題導(dǎo)出的次梯度投影算子可描述為T：Rn→Rn，

Tx=x-LW（x）-δ‖Q-1ATW（Ax-b）‖2Q

Q-1ATW（Ax-b），LW（x）>δ

x，LW（x）≤δ （4）

式中：

Q∈Rn×n為對稱正定矩陣；δ為常數(shù)，且滿足 {x∈Rn|LW（x）≤δ}≠。

將W和Q寫成分解形式W=（W1/2）T（W1/2），Q=（Q1/2）T（Q1/2）。設(shè)Q-1/2為Q1/2的逆矩陣，令

A

=W1/2AQ-1/2，x

=Q1/2x，

b

=W1/2b，則

‖Ax-b‖W=‖Ax-b‖

R‖x‖Q=‖x‖

（5）

記

L（x）=LW（Q-1/2x）=12‖

Ax-b‖2（6）

δ滿足levδ L={x=Q1/2x|L（x）≤δ，x∈Rn}≠。當(dāng)x=Q1/2xlevδ L時(shí)，式（2）中的次梯度投影算子可寫為

GL（x）=x-L（x）-δ‖Δ

L（x）‖2

Δ

L（x）=

x-L（x）-δ‖

AT（

Ax-b）‖2

AT（

Ax-b）=

Q1/2x-LW（x）-δ‖Q-1ATW（Ax-b）‖2QQ-1/2ATW（Ax-b）=

Q1/2x-LW（x）-δ‖Q-1ATW（Ax-b）‖2QQ-1ATW（Ax-b）

（7）

由于x=Q1/2x，因此對應(yīng)的次梯度投影迭代算子可以表示成式（4）。

式（4）的加權(quán)最小二乘問題的次梯度投影算子，對應(yīng)的迭代算法為

x（k+1）=T（x（k））=x（k）-

λ（x（k））Q-1ATW（Ax（k）-b）， k=0，1，2，…

（8）

式中：

λ（x（k））=

LW（x（k））-δ‖Q-1ATW（Ax（k）-b）‖2Q，LW（x（k））>δ

0，LW（x（k））≤δ

（9）

1.2 迭代算法步長的估算

若A的加權(quán)奇異值可分解為A=UVT，其中UTWU=Im，VTQ-1V=In，σmax=σ1≥σ2≥…≥σr=σmin，為A的r個(gè)非零奇異值，則式（8）的迭代步長λ（x（k））∈[0，12σ2min]。

利用奇異值分解，對x∈Rn和b∈Rm，有

x=PR（AT）（x）+P

N（A）（x）=

rj=1〈x，

Q-1

vj〉QQ-1vj+nj=r+1〈x，Q-1vj〉Q

Q

-1vj （10）

和

b=PR（

A）b+PN（

AT）b=

rj=1〈b，

uj〉W

uj+mj=r+1〈b，

uj〉W

uj（11）

式中：PR（

A）和PN（

A）分別為不同向量在A的值域和零空間上的斜投影，PR（

AT）和PN（

AT）分別為不同向量在AT的值域和零空間上的斜投影。由于常數(shù)δ滿足{x∈Rn|LW（x）≤δ}≠，因此根據(jù)

LW（x）-δ=12‖Ax-PR（

A）b‖2W+

12‖PN（

AT）b‖2W-δ （12）

可得12‖PN（

AT）b‖2W-δ<0。

由式（9）迭代步長的定義可知λ（x（k））≥0。另一方面，當(dāng)LW（x（k））>δ時(shí)，

λ（x（k））≤12‖Ax（k）-PR（A）b‖2W‖

Q-1ATW（Ax（k）-PR（

A）b）‖2Q=

12rj=1（〈x（k），Q-1vj〉QAQ-1vj-〈b，uj〉Wuj）2Wrj=1（〈x（k），Q-1vj〉QQ-1ATWAQ-1vj-〈b，uj〉WQ-1ATWuj）2Q=

12rj=1（〈x（k），Q-1vj〉Qσj-〈b，uj〉W）uj2Wrj=1（〈x（k），Q-1vj〉Qσ2j-〈b，uj〉Wσj）Q-1vj2Q=

12rj=1（〈x（k），Q-1vj〉Qσj-〈b，uj〉W）2rj=1（〈x（k），Q-1vj〉Qσj-〈b，uj〉W）2σ2j≤12σ2min （13）

因此λ（x（k））∈[0，12σ2min]。

2 快速迭代重建算法

在圖像重建問題中，A往往是一個(gè)大型稀疏矩陣。若將

A完全調(diào)入內(nèi)存后執(zhí)行后續(xù)計(jì)算，則內(nèi)存開銷太大，導(dǎo)致計(jì)算無法進(jìn)行。此外，將A存儲(chǔ)在外部設(shè)備（如硬盤）上，會(huì)導(dǎo)致每次讀取耗時(shí)太長。文獻(xiàn)[15]給出計(jì)算第i根直線與三維固定網(wǎng)格相交的交線長度和索引的算法，這里直接采用這一算法，給出A的第i行的非零元素和列標(biāo)號(hào)的索引。

仿照SART算法選取實(shí)對稱矩陣W和Q。[16]W和Q為對角陣且對角元滿足

Qj=mi=1|Aij|， j=1，2，…，n （14）

W-1i=nj=1|Aij|， i=1，2，…，m （15）

采用并行策略實(shí)現(xiàn)重建算法，利用次梯度投影算子求解線性方程組Ax=b的并行算法代碼如下，其中：輸入變量包括常數(shù)向量

bm×1，初始的解向量xn×1，最大迭代步數(shù)（MaxStep）和步長（delta）；輸出為解向量xn×1。

/* Begin： 對所有 PEs 同時(shí)執(zhí)行如下算法*/

/* （1）初始化*/

/*初始化 MPI， 得到相應(yīng)的 PE 個(gè)數(shù)：size， */

/*獲取當(dāng)前 PE 號(hào)：rank*/

Ite = 0； flag = true；

/* W， Q 初始化為 m 維、n維零向量*/

/* （2）計(jì)算初始的Ax-b和LW（x）-δ，并確定是否進(jìn)行迭代*/

val = 0.0；

for id = rank*m/size to （rank+1）*m/size-1 do

c[id]= -b[id]；

/* 獲取當(dāng)前第id條射線與三維固定網(wǎng)格交線的長度與index，分別存儲(chǔ)在 [matrix， ind]中*/

[matrix， ind]=

SystemMatrix3d（source[id]，detector[id]）；

len = length（ind）；

for j = 0 to len do

c[id]=c[id]+matrix[j]*x[ind[j]]；

W[id]=W[id]+matrix[j]；

Q[ind[j]]=Q[ind[j]]+matrix[j]；

end for

val += tc[id]*tc[id]/W[id]；

end for

/*用Allreduce操作求出所有PE中c， val， W， Q的值， 并廣播到所有PE中； */

val = val/2.0-delta；

/*若 val<=0， 則 flag = false； //不執(zhí)行迭代*/

/* （3）開始迭代運(yùn)算*/

while （flag&&（Ite

/* （3.1）計(jì)算Ac=AT（Ax-b）*/

/*初始化Ac為n維零向量*/

for id = rank*m/size to id<（rank+1）*m/size-1 do

[matrix， ind] =

SystemMatrix3d（source[id]，detector[id]）；

len = length（ind）；

for j = 0 to len do

Ac[ind[j]] += matrix[j]*c[id]；

end for

end for

/*用Allreduce操作求出所有PE中 Ac 的值，并廣播到所有PE 中； */

/* （3.2）計(jì)算 ‖Q-1AT（Ax-b）‖2Q*/

val2 = 0.0；

for i = rank*n/size to （rank+1）*n/size-1 do

val2 += Ac[i]*Ac[i]/Q[i]；

end for

/*用Allreduce操作求出所有PE中 val2 的值，

并廣播到所有PE中； */

/* （3.3）更新x的值*/

/*初始化tempX為n維零向量*/

for i = rank*n/size to （rank+1）*n/size-1 do

tempX[i]= x[i]- Ac[i]*val/val2/Q[i]；

end for

/*用Allreduce操作求出tempX的值， */

/*廣播到所有PE， 賦值給x； */

/* （3.4）計(jì)算新的Ax-b和LW（x）-δ， 并確定是否迭代執(zhí)行步驟 （2） 中除卻生成W， Q部分代碼； */

/* （3.5）更新迭代步數(shù)*/

Ite += 1；

end while

/* （4）結(jié)束 MPI*/

MPI_Finalize（）；

end

3 數(shù)值試驗(yàn)

采用三維低對比度的Shepp-Logan模型[7]進(jìn)行數(shù)值試驗(yàn)。該模型由10個(gè)橢球組成，橢球半軸長度分別為a、b、c，中心坐標(biāo)為（x0，y0，z0），長半軸與x軸正方向的夾角為，橢球內(nèi)部修正的吸收系數(shù)為A。具體參數(shù)見表1。

表 1 低對比度Shepp-Logan模型參數(shù)

采集數(shù)據(jù)方式模擬螺旋CT，采集協(xié)議中的各項(xiàng)參數(shù)見表2。

表 2 數(shù)據(jù)采集參數(shù)

重建過程中， 將待測物體離散化為

xRes=yRes=zRes=128（16）

最大迭代步數(shù)設(shè)置為 20， δ=（20/xRes）（m2/2n）。

數(shù)值試驗(yàn)在一臺(tái)Dell desktop機(jī)器上完成。處理器為Intel（R） Core（TM） i7-4790 CPU @ 3.60 GHz， 內(nèi)存為8 GiB。操作系統(tǒng)為Ununtu 15.04，gcc版本為4.9.2。

為對比重建結(jié)果，給出SART算法在不同松弛因子下第20步的重建結(jié)果和次梯度投影算法的重建結(jié)果，見圖1，所有圖像均在灰度窗[0.95，1.10]顯示。

（a1） 次梯度投影重建算法

（b1） 次梯度投影重建算法

（c1） 次梯度投影重建算法

（a2） SART，λ=0.5

（b2） SART，λ=0.5

（c2） SART，λ=0.5

（a3） SART，λ=1.0

（b3） SART，λ=1.0

（c3） SART，λ=1.0

（a4） SART，λ=1.5

（b4） SART，λ=1.5

（c4） SART，λ=1.5

（a5） SART，λ=2.0

a） z=-3.750 0 cm

（b5） SART，λ=2.0

b） y=-3.125 0 cm

（c5） SART，λ=2.0

c） x=0.312 5 cm

圖 1 不同算法和松弛因子的重建結(jié)果

由于2個(gè)算法迭代停止準(zhǔn)則不同，所以僅從重建圖像上比較次梯度投影算法和SART算法，而沒有顯示兩者的均方誤差等信息。從重建圖像看，SART算法的重建結(jié)果受松弛因子影響很大。SART算法中目標(biāo)函數(shù)隨迭代次數(shù)變化情況見圖2。

針對次梯度投影迭代算法的并行策略，算法的加速比S和并行效率E表示為

S=T1/Tp

E=S/p

（17）

式中：p為處理器個(gè)數(shù)；T1和Tp分別為單個(gè)處理器和p個(gè)處理器的計(jì)算時(shí)間，以此確定并行算法的性

能。 利用次梯度投影迭代算法進(jìn)行數(shù)值試驗(yàn)，算法執(zhí)行20步迭代的計(jì)算時(shí)間T、加速比S和并行效率E見表3。數(shù)值試驗(yàn)中S和E與處理器數(shù)量的關(guān)系見圖3。

圖 2 SART算法中目標(biāo)函數(shù)隨迭代次數(shù)變化情況

表 3 數(shù)值試驗(yàn)的計(jì)算時(shí)間、加速比和并行效率

圖 3 數(shù)值試驗(yàn)中S和E與處理器數(shù)量的關(guān)系

從統(tǒng)計(jì)結(jié)果看，隨著處理器數(shù)量的增加，S增加，但是E下降非?？?。在單機(jī)上運(yùn)行時(shí)，隨著處理器數(shù)量的增加，內(nèi)存不足尤為明顯，內(nèi)存與交換空間的交互限制了S和E。此外，當(dāng)處理器數(shù)量增加時(shí)，處理器的負(fù)載不平衡導(dǎo)致部分處理器計(jì)算時(shí)間長，而其他處理器等待時(shí)間長，更顯示出總體運(yùn)行時(shí)間減少緩慢。為提升這方面的性能，需要重新分割光源點(diǎn)。

4 結(jié)束語

將線性成像問題

Ax=

b轉(zhuǎn)化為加權(quán)最小二乘問題進(jìn)行求解。 利用加權(quán)最小二乘問題的目標(biāo)函數(shù)， 導(dǎo)出次梯度投影迭代算子， 從理論上估算次梯度投影迭代算法迭代步長。 根據(jù)成像問題， 結(jié)合次梯度投影算法， 設(shè)計(jì)算法的快速實(shí)現(xiàn)策略。根據(jù)X-Ray CT成像理論，模擬生成三維物體的投影數(shù)據(jù)，利用快速重建算法進(jìn)行計(jì)算，并與重建問題中常用的SART算法進(jìn)行比較，得到可信的重建結(jié)果。分析重建結(jié)果，次梯度投影算法和SART算法都能得到清晰的重建圖像，但SART算法需要調(diào)節(jié)多組松弛因子以確保得到更佳的重建結(jié)果。同時(shí)，分析次梯度投影快速算法的加速比和并行效率，確定該快速算法在單機(jī)上是高效的。

參考文獻(xiàn)：

[1] KAK A C， SLANEY M. Principles of computerized tomographic imaging[M]. Philadelphia： Society of Industrial and Applied Mathematics， 2001.

[2] FELDKAMP L A， DAVIS L C， KRESS J W. Practical cone-beam algorithm[J]. Journal of the Optical Society of American： A， 1984， 1（6）： 612-619.

[3] KATSEVICH A. Theoretically exact filtered backprojection-type inversion algorithm for spiral CT[J]. SIAM Journal on Applied Mathematics， 2002， 62（6）： 2012-2026.

[4] DEMPSTER A P， LAIRD N M， RUBIN D B. Maximum likelihood from incomplete data via EM algorithm[J]. Journal of the Royal Statistical Society： Series B（Methodological）， 1977， 39（1）： 1-38.

[5] JIANG M， WANG G. Development of iterative algorithms for image reconstruction[J]. Journal of X-Ray Science and Technology， 2002， 10（1）： 77-86.

[6] ANDERSEN A H， KAK A C. Simultaneous algebraic reconstruction technique（SART）： A superior implementation of ART algorithm[J]. Ultrasonic Imaging， 1984， 6（1）： 81-94.

[7] CIMMINO G. Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari[J]. La Ricerca Scientifica， 1938（9）： 326-333.

[8] CENSOR Y， DAN G， GORDON R. Component averaging： An efficient iterative parallel algorithm for large and sparse unstructured problems[J]. Parallel Computing， 2001， 27（6）： 777-808.

[9] CENSOR Y， ELFVING T. Block-iterative algorithms with diagonally scaled oblique projections for linear feasibility problem[J]. SIAM Journal on Matrix Analysis and Applications， 2002， 24（1）： 40-58.

[10] CENSOR Y， ELFVING T， HERMAN G T， et al. On diagonally relaxed orthogonal projection methods[J]. SIAM Journal on Scientific Computing， 2007， 30（1）： 473-504.

[11] CEGIELSKI A. Iterative Methods for Fixed Point Problems in Hilbert Spaces[M]. Berlin： Springer-Verlag， 2013.

[12] BAUSCHKE H H， BORWEIN J M. On projection algorithms for solving convex feasibility problems[J]. SIAM Review， 1996， 38（3）： 367-426.

[13] COMBETTES P L. Convex feasibility problem in image recovery[J]. Advances in Imaging and Electron Physics， 1996， 95（8）： 155-270. DOI： 10.1016/S1076-5670（08）70157-5.

[14] BAUSCHKE H H， WANG C F， WANG X F， et al. On subgradient projectors[J]. SIAM Journal on Optimization， 2015， 25（2）： 1064-1082.

[15] 張曉媛， 王彩芳， 殷緒導(dǎo)， 等. 直線與三維固定網(wǎng)格的交線長度與索引公式計(jì)算[J]. 應(yīng)用數(shù)學(xué)進(jìn)展， 2017， 6（4）： 496-503. DOI： 10.12677/AAM.2017.64059.

[16] JIANG M， WANG G. Convergence of simultaneous algebraic reconstruction technique（SART）[J]. IEEE Transactions on Image Processing， 2003， 12（8）： 957-961. DOI： 10.1109/TIP.2003.815295.

[17] WANG G， LIN T H， CHENG P， et al. A general cone-beam reconstruction algorithm[J]. IEEE Transactions on Medical Imaging， 1993，12（3）： 486-496. DOI： 10.1109/42.241876.