次梯度投影算子求解图像重建问题的并行策略

王彩芳 殷绪导 张晓媛 徐兆亮

摘要：

将次梯度投影迭代算法应用到数字图像重建问题。将图像重建问题转化为求一个加权最小二乘问题，导出次梯度投影算子在该问题下的具体迭代形式，并采用并行计算策略重建算法。通过三维数值实验对比次梯度投影迭代算法与常用的SART算法，验证算法的可行性和效率。

关键词：

图像重建； 加权最小二乘； 次梯度投影算子； 并行计算

中图分类号： O244；O29

文献标志码： 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

基金项目：

国家自然科学基金（11401372）

作者简介：

王彩芳（1981—），女，江苏苏州人，副教授，博士，研究方向为数字图像重建算法，（E-mail）cfwang@shmtu.edu.cn

0 引 言

数字图像重建问题是根据待测物体边界测量数据获取待测物体内部结构的一类反问题。[1]近年来，重建问题的研究主要集中在如何找到一种合适的扫描方式和快速的成像算法。以三维CT为例，基于FDK公式[2]的近似重建算法和基于Katsevich公式[3]的精确重建算法已经在成像领域得到广泛应用。这类分析重建算法虽然成像速度快，但基本在射线变换和三维Radon变换的框架下推演，对数据的完整性要求高。当部分投影数据缺失或获取存在困难时，分析算法设计极为复杂，因此代数算法渐渐进入研究者的视野。随着计算机技术的高度发展，代数重建算法耗时长的缺点也有所改善。

为利用代数重建算法，很多成像问题都用线性方程组进行描述。

Ax=

b（1）

式中：

A∈Rm×n為系数矩阵，称为成像矩阵；未知数x∈Rn表示待测物体对入射光的吸收系数或其他光学系数；右端项b∈Rm为边界测量数据或修正后的边界测量数据。图像重建问题即为求解光学系数x的反问题。在图像重建问题中，系数矩阵A规模很大，且Ax=b是病态问题。EM迭代算法和Landweber 迭代算法是图像重建问题中应用广泛的两类算法。EM算法主要用于求解观测数据中包含的隐藏参数的极大似然估计或极大先验估计。[4]在成像问题中，其本质为寻找使得观测数据与计算数据之间 Kullback-Leibler（KL）距离最小的解。[5]Landweber 迭代算法包括一系列著名算法：SART算法[6]、CIMMINO方法[7]、CAV算法[8]、DWE算法[9]和DROP方法[10]等。其本质是寻找上述线性系统的加权最小二乘解[5]。Landweber算法的收敛性受迭代过程中松弛因子的影响。次梯度投影迭代算法是求解凸可行性问题的一类重要算法[11]，因而被广泛应用到实际问题的求解中[12-13]。与Landweber算法不同，次梯度投影算法在迭代过程中无需预设松弛因子，迭代步长完全由目标函数和当前估计值所确定。因此，考虑将次梯度迭代算法应用到线性图像重建问题中，研究该算法的快速重建方法。

本文回顾次梯度投影算法，并与重建问题相结合，导出具体的迭代格式，估算迭代步长，研究迭代算法的并行策略，得到具体算法实现过程。利用三维数值实验验证算法的可行性，并将重建结果与SART算法进行比较。

1 重建算法导出

假设f：Rm→R为凸且可微的函数，将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）为函数f（x）在点x处的次梯度；‖·‖为Rm上的2-范数。根据次梯度投影算子的性质，其不动点集合满足Fix（Gf）={x∈Rm|f（x）≤δ}。

针对式（1）病态问题，应考虑其对应的加权最小二乘问题，并导出次梯度投影算子的具体表达形式，计算原理如下。

1.1 考虑加权最小二乘问题

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

式中：W∈Rm×m为对称正定矩阵。该问题导出的次梯度投影算子可描述为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为对称正定矩阵；δ为常数，且满足 {x∈Rn|LW（x）≤δ}≠。

将W和Q写成分解形式W=（W1/2）T（W1/2），Q=（Q1/2）T（Q1/2）。设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}≠。当x=Q1/2xlevδ L时，式（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，因此对应的次梯度投影迭代算子可以表示成式（4）。

式（4）的加权最小二乘问题的次梯度投影算子，对应的迭代算法为

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的加权奇异值可分解为A=UVT，其中UTWU=Im，VTQ-1V=In，σmax=σ1≥σ2≥…≥σr=σmin，为A的r个非零奇异值，则式（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的值域和零空间上的斜投影。由于常数δ满足{x∈Rn|LW（x）≤δ}≠，因此根据

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

A）b‖2W+

12‖PN（

AT）b‖2W-δ （12）

可得12‖PN（

AT）b‖2W-δ<0。

由式（9）迭代步长的定义可知λ（x（k））≥0。另一方面，当LW（x（k））>δ时，

λ（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往往是一个大型稀疏矩阵。若将

A完全调入内存后执行后续计算，则内存开销太大，导致计算无法进行。此外，将A存储在外部设备（如硬盘）上，会导致每次读取耗时太长。文献[15]给出计算第i根直线与三维固定网格相交的交线长度和索引的算法，这里直接采用这一算法，给出A的第i行的非零元素和列标号的索引。

仿照SART算法选取实对称矩阵W和Q。[16]W和Q为对角阵且对角元满足

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

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

采用并行策略实现重建算法，利用次梯度投影算子求解线性方程组Ax=b的并行算法代码如下，其中：输入变量包括常数向量

bm×1，初始的解向量xn×1，最大迭代步数（MaxStep）和步长（delta）；输出为解向量xn×1。

/* Begin： 对所有 PEs 同时执行如下算法*/

/* （1）初始化*/

/*初始化 MPI， 得到相应的 PE 个数：size， */

/*获取当前 PE 号：rank*/

Ite = 0； flag = true；

/* W， Q 初始化为 m 维、n维零向量*/

/* （2）计算初始的Ax-b和LW（x）-δ，并确定是否进行迭代*/

val = 0.0；

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

c[id]= -b[id]；

/* 获取当前第id条射线与三维固定网格交线的长度与index，分别存储在 [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； //不执行迭代*/

/* （3）开始迭代运算*/

while （flag&&（Ite

/* （3.1）计算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）计算 ‖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）计算新的Ax-b和LW（x）-δ， 并确定是否迭代执行步骤 （2） 中除却生成W， Q部分代码； */

/* （3.5）更新迭代步数*/

Ite += 1；

end while

/* （4）结束 MPI*/

MPI_Finalize（）；

end

3 数值试验

采用三维低对比度的Shepp-Logan模型[7]进行数值试验。该模型由10个椭球组成，椭球半轴长度分别为a、b、c，中心坐标为（x0，y0，z0），长半轴与x轴正方向的夹角为，椭球内部修正的吸收系数为A。具体参数见表1。

表 1 低对比度Shepp-Logan模型参数

采集数据方式模拟螺旋CT，采集协议中的各项参数见表2。

表 2 数据采集参数

重建过程中， 将待测物体离散化为

xRes=yRes=zRes=128（16）

最大迭代步数设置为 20， δ=（20/xRes）（m2/2n）。

数值试验在一台Dell desktop机器上完成。处理器为Intel（R） Core（TM） i7-4790 CPU @ 3.60 GHz， 内存为8 GiB。操作系统为Ununtu 15.04，gcc版本为4.9.2。

为对比重建结果，给出SART算法在不同松弛因子下第20步的重建结果和次梯度投影算法的重建结果，见图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 不同算法和松弛因子的重建结果

由于2个算法迭代停止准则不同，所以仅从重建图像上比较次梯度投影算法和SART算法，而没有显示两者的均方误差等信息。从重建图像看，SART算法的重建结果受松弛因子影响很大。SART算法中目标函数随迭代次数变化情况见图2。

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

S=T1/Tp

E=S/p

（17）

式中：p为处理器个数；T1和Tp分别为单个处理器和p个处理器的计算时间，以此确定并行算法的性

能。 利用次梯度投影迭代算法進行数值试验，算法执行20步迭代的计算时间T、加速比S和并行效率E见表3。数值试验中S和E与处理器数量的关系见图3。

图 2 SART算法中目标函数随迭代次数变化情况

表 3 数值试验的计算时间、加速比和并行效率

图 3 数值试验中S和E与处理器数量的关系

从统计结果看，随着处理器数量的增加，S增加，但是E下降非常快。在单机上运行时，随着处理器数量的增加，内存不足尤为明显，内存与交换空间的交互限制了S和E。此外，当处理器数量增加时，处理器的负载不平衡导致部分处理器计算时间长，而其他处理器等待时间长，更显示出总体运行时间减少缓慢。为提升这方面的性能，需要重新分割光源点。

4 结束语

将线性成像问题

Ax=

b转化为加权最小二乘问题进行求解。 利用加权最小二乘问题的目标函数， 导出次梯度投影迭代算子， 从理论上估算次梯度投影迭代算法迭代步长。 根据成像问题， 结合次梯度投影算法， 设计算法的快速实现策略。根据X-Ray CT成像理论，模拟生成三维物体的投影数据，利用快速重建算法进行计算，并与重建问题中常用的SART算法进行比较，得到可信的重建结果。分析重建结果，次梯度投影算法和SART算法都能得到清晰的重建图像，但SART算法需要调节多组松弛因子以确保得到更佳的重建结果。同时，分析次梯度投影快速算法的加速比和并行效率，确定该快速算法在单机上是高效的。

参考文献：

[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] 張晓媛， 王彩芳， 殷绪导， 等. 直线与三维固定网格的交线长度与索引公式计算[J]. 应用数学进展， 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.