一种反投影模型下的图像重建迭代算法

赵芳薇，邱 钧，刘 畅

（北京信息科技大学 应用数学研究所，北京100101）

0 引 言

图像重建技术在工业、医学等很多领域有着广泛的应用［1］，其本质是利用获取的投影数据重构被检测目标内部的高分辨率断层图像.常见的图像重建方法有解析算法和迭代算法［2-3］，其中迭代算法将投影数据对应的积分方程离散化，把图像重建问题直接转化为线性方程组求解问题，在处理不完全投影数据及噪声污染较为严重的投影数据时较解析算法更有优势.目前迭代算法有诸多理论上的研究成果，是一种被广泛应用的图像重建算法［4-8］.

这里求解线性方程组是不适定问题，且由于噪声等其他因素的影响，传统求解线性方程组的方法并不适用.基于超平面投影和目标优化理论建立的迭代算法，在由投影重建图像中应用广泛.ART（Algebraic Reconstruction Technique）［9］、SIRT（Simultaneous Iteration Reconstruction Technique）［10］、CAV（Component Averaging）［11］等都是常见的图像重建迭代算法，不同的迭代算法对应着不同的迭代校正格式，也对应着对误差不同的分配校正方式.

现有的迭代校正格式通常是基于对投影模型及其几何物理意义的分析，在迭代过程中需计算投影系数矩阵的投影系数［12］.其中，对于投影系数的计算，常用的一类方法是计算射线穿过剖分格的截格长度［1，13］，另一类是引入重建点模型，计算重建点对射线的贡献因子［14-15］.通常的迭代重建算法可视为用迭代方法给出广义敏感逆矩阵估计形式［16］，也有研究者通过重新刻画图像重建算子，直接寻找广义敏感逆矩阵的近似表达形式，避免了迭代中出现的半收敛现象，并且有利于节省图像重建时间［17］.

本文从与投影模型相对应的反投影模型出发，研究全体扫描射线对重建点的贡献因子，用迭代方法给出广义逆矩阵估计形式，建立了广义敏感逆矩阵在某种意义下的近似，形成了迭代成像算法.本文关于反投影模型下快速图像重建算法的研究，进一步可以推广到广义敏感逆矩阵的较为精确的迭代估计，以期形成基于扫描系统特征的广义逆敏感矩阵形式的计算模版，从而形成有效的直接计算成像方法.

1 重建点离散化模型

依据图像重建离散化方式的不同，线性方程组也会存在差别.常见的图像重建离散化模型有基于网格剖分的离散化模型［1，13］和重建点离散化模型［14-15］，本文引入重建点离散化模型，如图1所示.

图1 重建点离散化模型示意图 Fig.1 Diagram of discrete model based on reconstruction points

重建点离散化模型淡化剖分网格的概念，重建点对应着像素点，重建点周围的一个理想小区域内像素值视为定值，用重建点到射线的距离的函数刻画像素对扫描射线的贡献因子.较之基于网格剖分的离散化模型，该模型更好地刻画了显示模式、扫描模式以及坐标架之间的关系，减少了离散误差，具有很好的适用性.

记待重建图像被离散化为N维有限向量X=（x1，x2，…，xN）T，xj是第j个重建点处的像素值，投影数据向量为b=（b1，b2，…，bM）T.重建点离散化模型下，图像重建问题可归结为线性方程组求解问题

式中：C=（cij）M×N是投影系数矩阵，投影系数cij=f（dij，φ）表示第j个重建点对第i条射线Li衰减的贡献因子，dij是第j个重建点到射线Li的垂直距离，φ标定射线入射方向，i∈｛1，2，…，N｝，j∈｛1，2，…，M｝，如图2所示.

图2 投影系数矩阵示意图 Fig.2 Diagram of the projection coefficient matrix

2 反投影模型

反投影算法的本质是断层平面内某一点的值可看作全体是投影数据的加权求和.

设待重建区域被离散化为N个理想小区域，有M个有效投影数据.反投影模型刻画了依据不同扫描射线对某重建点的贡献因子大小，进行图像重建的过程.记tji是第i根射线对第j个重建点的贡献因子，则第j个重建点像素值可记为式中：bi是第i根射线的投影数据；tji为反投影系数，i∈｛1，2，…，N｝，jš｛1，2，…，M｝.如图3所示.

图3 反投影系数矩阵示意图 Fig.3 Diagram of the back-projection coefficient matrix

式（2）的矩阵表示形式为

是由全体的反投影系数组成的反投影系数矩阵［18］.

在全体投影数据对重建点的贡献因子已知的情形下，图像重建问题归结为由全体投影向量b和反投影系数矩阵T直接计算图像X的问题.

3 反投影系数矩阵T的估计

利用反投影模型对图像进行重建的关键是对反投影系数矩阵给出精确的估计，本文依据反投影模型中扫描射线和重建点间的几何位置关系，及其反投影系数的几何物理意义给出一种合理的估计方法.

反投影系数tji是第i根射线对第j个重建点的贡献因子，在计算全体射线对某固定重建点的贡献因子时，由于距离重建区域中心不同的扫描射线穿过待重建区域的长度不同，需先对全体投影数据进行均衡.记第i条射线Li被待重建区域所截长度是si，探测器检测到的投影数据为bi，则投影数据可被均衡化为bi／si.均衡化后的M维投影向量可记b′=（b1／s1，…，bM／sM）T∈RM.

对于均衡化后的投影数据，不同射线对重建点的贡献因子与射线到重建点的距离有关，且随着射线到重建点的距离逐渐变小而逐渐变大.记f（dij）是第i根射线到第j个重建点的距离dij的函数，在考虑均衡化后的全体射线对第j个重建点的贡献系数时，第i根射线对第j个重建点的贡献因子，即投影系数tji的表达形式为

式中：si是第i根射线被重建区域截得的长度；dij是第j个重建点到第i根射线的距离；f（dij）是dij的函数.si与dij如图4所示.

图4 反投影系数矩阵影响参数示意图 Fig.4 Diagram of the influence parameters in back-projection coefficient matrix

4 反投影模型下的图像重建迭代算法

反投影模型下的图像重建迭代算法依赖于对反投影系数的估计.算法的分量迭代格式为

式中：M为投影射线总数；N为重建点的总数；bi为投影数据为计算投影值；λk是松弛

因子，i=1，2，…，M，j=0，1，2，…，N，k=0，1，2，….

迭代算法步骤为：

1）赋予待重建图像向量初始值X（0）=（0，0，…，0）T；

2）依据初始值X（0），计算投影值与投影数据bi，按照式（5）进行迭代.k=0，1，2，….

5 算法验证的数据实验

本节采用两组模拟数据验证算法，数据实验中f（dij）选取线性插值函数.

模拟数据角度采样数为576，平行线采样数为721，模型和重建结果的图像分辨率均为600×600.

图5 Shepp-Logan原图及第200行位置 Fig.5 Diagram for the Shepp-Logan and the position of line 200

实验引入归一化均方距离d，归一化平均绝对值距离r和最坏情况距离e对重建结果进行评价.

5.1 模拟数据实验1

模拟实验1采用经典的Shepp-Logan模型.迭代5，10，20，30次的重建结果及第200行对应的像素曲线如图7所示，与原图的距离如表1所示.

图6 Shepp-Logan的投影数据位图 Fig.6 The bitmap of the Shepp-Logan's projection data

图7 实验1反投影模型下迭代算法迭代不同轮次的重建结果 Fig.7 Reconstruction results of the new iterative algorithm under different numbers of iterations for experiment one

表1 实验1反投影模型下迭代算法误差分析 Tab.1 Errors analysis of the new iterative algorithm for experiment one

5.2 模拟数据实验2

模拟数据实验2选取鹦鹉模型如图8和图9所示.迭代5，10，20，30次的重建结果及第200行对应的像素曲线如图10所示，与原图的距离如表2所示.

图8 鹦鹉原图及第200行位置 Fig.8 Diagram for the parrot and the position of line 200

图9 鹦鹉的投影数据位图 Fig.9 The bitmap of the parrot's projection data

图10 实验2反投影模型下迭代算法迭代不同轮次的重建结果 Fig.10 Reconstruction results of the new iterative algorithm under different numbers of iterations for experiment two

两组模拟实验的数据结果表明：随着迭代轮次的增加，成像精度逐步提高，细节也愈加明显.从像素曲线可以看出重建结果有较好的密度分辨率，空间分辨率也随着迭代次数的增加有所改善.

表2 实验2反投影模型下迭代算法误差分析 Tab.2 Errors analysis of the new iterative algorithm for experiment two

6 结 论

本文从与投影模型相对应的反投影模型出发，给出了一种反投影模型下的图像重建迭代算法.算法基于反投影系数矩阵的估计，是广义敏感逆矩阵在某种意义下的近似，反投影系数的估算考虑了全体扫描射线对重建点的贡献因子，实现了误差的全局分配，重建结果具有较好的密度分辨率.全体扫描射线对重建点的贡献因子，是投影与反投影模型中的基本刻画因子，反映了投影过程中重建点的物理与几何的基本性质和特征，建立对其的刻画模型是用迭代方法给出广义敏感逆矩阵估计形式的关键.

反投影系数的估计方法直接影响算法的收敛速度与成像精度，如何给出更合理的反投影系数估计方法，进一步有效地提高重建结果的空间分辨率，值得进一步研究.本文提出的反投影模型下的图像重建迭代方法，可以推广到广义敏感逆矩阵的较为精确的迭代估计，以期形成基于扫描系统特征的广义逆敏感矩阵形式的计算模板，从而形成有效的快速实用的直接计算成像方法，进一步丰富研究的内容.

［1］庄天戈.CT原理与算法［M］.上海：上海交通大学出版社，1992.

［2］Herman G T.Fundamentals of computerized tomography：image reconstruction from projections［M］.2ed.Germany：Springer，2009.

［3］Jiang Ming，Wang Ge.Development of iterative algorithms for image reconstruction［J］.Journal of X-ray Science and Technology，2002，10（1／2）：77-86.

［4］Herman G T.Algebraic reconstruction techniques can be made computationally efficient［J］.Medical Imaging（0278-0062），1993，12（3）：600-609.

［5］Jiang Ming，Wang Ge.Convergence of iterative algo-rithms for image reconstruction［C］.Bio-medical Imaging，2002：685-688.

［6］邱钧，徐茂林.由投影重建图像的对称块迭代算法［J］.电子与信息学报，2007，29（10）：2296-2300.Qiu Jun，Xu Maolin.A method of symmetric block-iterative for image reconstruction［J］.Journal of Electronics ffInformation Technology（1009-5896），2007，29（10）：2296-2300.（in Chinese）

［7］毕丹丹，李筠.有关ART算法中松驰因子的研究［J］.信息技术，2014（9）：121-124，128.Bi Dandan，Li Jun.Research on the relaxation parameter of algebraic reconstruction technique［J］.Information Technology，2014（9）：121-124，128.（in Chinese）

［8］党晓军，韦孟伏.CT代数重建算法的快速实现［J］.核电子学与探测技术，2011，31（10）：1104-1106.Dang Xiaojun，Wei Mengfu.A fast implementation of CT statistical reconstruction methods［J］.Nuclear Electronics ffDetection Technology，2011，31（10）：1104-1106.（in Chinese）

［9］Gorden R，Bender R，Herman G T.Algebraic reconstruction techniques（ART）for three-dimensional electron microscopy and X-ray photograph［J］.Journal of Theoretical Biology，1970，29（3）：471-476.

［10］Gilbert P.Iterative methods for the three-dimensional reconstruction of an object from projections［J］.Journal of Theoretical Biology，1972，36（2）：105-117.

［11］Censor Y，Gordon D，Gordon R.Component averaging：an efficient iterative parallel algorithm for large and sparse unstructured problems［J］.Parallel Computing，2001，27：777-808.

［12］陈建林，闫镔，李雷，等.CT重建中投影矩阵模型研究综述［J］.CT理论与应用研究，2014，23（2）：317-328.Chen Jianlin，Yan Bin，Li Lei，et al.Reviews of the model of projection matrix of CT reconstruction algorithm［J］.CT Theory and Applications，2014，23（2）：317-328.（in Chinese）

［13］查国震.基于正六边形像素的扇束等距滤波反投影及平行束插值代数重建算法研究［D］.北京：首都师范大学，2009.

［14］徐键，邱钧.一种基于计算点的图像重建离散化模型［C］.应用数学暑期研讨会论文集，电子工业出版社，2011.

［15］孙守思，邱钧，桂志国，等.重建点模型的EM迭代成像［J］.中北大学学报（自然科学版），2014，35（2）：209-217.Sun Shousi，Qiu Jun，Gui Zhiguo，et al.The EM iterative imaging on the discrete reconstruction points［J］.Journal of North University of China（Natural Science Edition）.2014，35（2）：209-217.（in Chinese）

［16］王超，王化祥.电阻抗断层图像重建算法研究——预迭代算法提出［J］.信号处理，2002，18（6）：547-550.Wang Chao，Wang Huaxiang.Research on the reconstruction algorithm of electrical impedance tomography——the pre-iteartion reconstruction algorithm［J］.Singal Processing，2002，18（6）：547-550.（in Chinese）

［17］张兆田.不完全数据断层成像理论与应用研究［D］.北京：中国科学院，2005.

［18］柯丽，曹冯秋，杜强.MIT中反投影矩阵的计算与数据处理方法［J］.仪器仪表学报，2014，35（10）：2256-2262.Ke Li，Cao Fengqiu，Du Qiang.Back-projection matrix calculation and data processing methods used in magnetic induction tomography［J］.Chinese Journal of Scientific Instrument，2014，35（10）：2256-2262.（in Chinese"）