基于ART算法的投影系数快速计算方法

侯慧玲，王明泉，杨 娟，李世虎

（中北大学 仪器科学与动态测试教育部重点实验室，山西 太原030051）

0 引 言

CT技术是无损检测领域广泛使用的检测手段，通过设计和搭建CT成像系统完成对被检工件的扫描［1］.对采集到的投影进行重建的算法主要有两大类：解析算法和迭代算法.滤波反投影法（Filtered Back Projection，FBP）是最基本的解析算法，运算量小，重建速度快，但对投影数据的完备性要求高.代数重建技术（ART）是经典的代数迭代算法，在投影数据部分缺失的情况下仍然适用，且重建图像质量高，但其计算量大，收敛速度慢，限制了其在实时工业检测系统中的应用［2-3］.

影响ART算法重建速度的因素有很多，如：投影数据的访问方式、图像初值的选择、松弛因子的选择、投影系数的计算等等［4-6］.其中，投影系数的计算在时间上占比最大［7］，因而有不少学者关注投影数据计算方法的研究，不断有新的方法提出［8-12］，其中最为熟知的是Siddon于1985年提出的基于长度加权的基本计算方法，计算射线穿过像素网格的长度，计算耗时比较长［8］；2006年张顺利提出投影系数快速计算方法，通过增量计算来确定射线穿过的网格编号及交线长度，重建速度比Siddon算法快了将近6倍，但引入频繁的分支判断［9］；2012年杨文良提出的快速计算方法，使用简单的加减法及比较运算来完成投影系数计算，重建效率进一步提高［11］.

本文针对ART重建中投影系数计算时间冗长的问题，从研究射线穿过网格的相交规律出发，提出新的计算方法，力求在保证投影重建图像质量的前提下，减少运算量及分支判断，进一步提高计算速度，加快重建效率.

1 ART算法基本原理

式中：i=1，2，3，…，M；j=1，2，3，…，N；pi为射线i的投影值；pij为像素i对像素j的投影贡献；wij为投影系数；xj为像素j的像素值，即该像素对X射线的衰减系数.直接求解该线性方程组极其困难.ART算法是通过迭代方式来求解该方程组的最优解的过程，其迭代公式［2-3］为

式中：k是迭代次数；wij为投影系数；λk为松弛因子，λ∈（0，2）.

计算其中的投影系数wij有三种数学模型：点加权型，面积加权型以及长度加权型.

长度加权模型如图1所示.在模型中，射线被看作是一条直线穿过网格，投影系数wij定义为射线i穿过像素j的交线长度.

就重建精度而言，长度加权模型低于面积加权型，优于点加权型.在保证一定重建精度的情况下，计算速度快，计算复杂度低，所以在实际中运用最多.本文便是基于长度加权模型研究投影系数的快速计算方法.

图1 长度加权模型 Fig.1 Length-weighted model

2 投影系数的快速计算

2.1 二维情况下的投影系数快速计算

投影系数矩阵里非零值的个数较少，可用稀疏矩阵进行描述，为节省存储空间，运算中仅保存非零值和对应的空间坐标（算法执行中以记录网格编号代替）.

图2为射线投影对应的某空间二维坐标系，其重建区域面积为N×N，网格编号规则为从左下角开始，分别命名为0，1，2…N2-1，共有N2个网格，其中每个网格长度为L，文中为简化分析，设L=1（即对射线源坐标和探测器坐标进行归一化处理）.

为了求解射线在重建区域不同位置的投影系数，这里设区域中点O为原点坐标，点M和N分别为入射点和出射点位置，Ni中i为出射点个数，ki为不同射线对应的斜率，PXc，PYc分别表示射线与X轴、Y轴的交点，下标c代表交点的个数.下面以图2给出的投影示意图为例，分析射线穿过重建区域的网格编号以及对应的交线长度.

从图2中可以看出，射线与网格的交点主要有两类，一类是与Y轴的交点，一类是与X轴的交点.图中射线MN1与网格的交点不仅包含与Y轴的交点，也包含与X轴的交点，但MN2为其特例，仅包含与Y轴的交点.为了具有普遍意义，本节主要以MN1为例，求解其与X轴以及Y轴的点坐标.

ART算法是将连续图像离散化为J=N×N的像素单元，待求解图像f（x，y），所有穿过图像的射线条数总计为M，则总的射线投影方程为

图2 投影系数分析原理图 Fig.2 Schematic diagram of projection coefficients

设MN1的斜率0＜k1≤1，与Y轴交点集合用数组PY（PYcx，PYcy）表示，PYcx表示与X轴相交的第c点的横坐标值，代表交线长度；PYcy表示与Y轴相交的第c点的纵坐标值，代表网格编号；与X轴交点集合用数组PX（PXcx，PXcy）表示，PXcx表示与X轴相交的第c点的横坐标值，代表交线长度；PXcy表示与X轴相交的第c点的纵坐标值，代表网格编号.

依据MN1在图2中的位置关系，式（4）给出了数组PY中坐标点的表达式

式中：c∈［1，2n］，表示向下取整，（PY1x，PY1y）表示入射点M的起始坐标.

式（5）给出了数组PX中交线长度的坐标点表达式

式中：（PX1x，PX1y）表示入射点M的起始坐标，这里（PX1x，PX1y）=（PY1x，PY1y）.

式（5）中没有给出数组PX中点的网格编号，依据图2分析发现：A与B两点的网格编号相差N.由于有许多类似B点的点出现，会导致式（4）中数组的交点数增加，这里定义新集合数组为PY′（PY′cx，PY′cy），式（4）的坐标表达式更新为

数组PY′中点的集合包含射线MN1穿过重建区域所有点的交线长度及对应的网格编号，由此即可得到射线在投影区域的所有投影系数.

式（4）～式（6）给出了图2中射线斜率0＜k1≤1时的投影系数计算过程，对于k1的其它几种情况，下面给出了分析说明.

1）当-1≤k1＜0时，从图2可以看出-1≤k1＜0时的射线与0＜k1≤1时的射线关于X轴对称，因此可将式（4）～式（6）中的纵坐标进行对称变换即可得到-1≤k＜0时的投影系数，其结果按小到大的顺序可排列为

2）当k1＞1时，从图2可以看出射线与X轴的交点较多，因此可按X轴方向进行计算，结果即是将式（6）中的横坐标改为纵坐标；

3）当k＜-1时，即是将k1＞1中的网格编号做关于X轴的对称运算即可；

4）当k=0或k→∞时，射线与Y轴或X轴平行，只要确定了M点的坐标，之后所有点的坐标经过简单叠加即可.

由于ART算法主要是被射线驱动，一次迭代仅记录一条射线的信息，实际处理中可分次处理不同射线的信息，这有利于节省内存空间，提高算法重建速度.

2.2 三维情况下的投影系数快速计算

上面分析了二维情况下的投影系数计算方法，实际中被投影物体往往是三维状态，为此需将图2中的二维坐标系扩展为三维坐标系，其中O为物体中心，也是坐标原点，三维空间可将物体划分为N层，每层包含N×N个立方体（1个立方体即可视为物体的一个体素单元），N层共N×N×N个立方体，每个正方体的具体编号可根据坐标轴正方向进行排列.

三维投影系数求解的基本思想为：首先将射线在三维坐标下的投影映射到XOY和Z两个二维平面，其次分别求出射线在XOY平面、Z平面的投影系数，最后将两者获得的结果进行综合分析，即可得到射线的三维投影系数.其中XOY平面的投影系数为层内索引，具体求解可采用2.1节方法；Z平面的投影系数为层间索引，层间索引（Z方向索引）主要是计算射线穿过被测物体的层信息.实际处理中，层信息通过计算射线在XOZ平面或YOZ平面的投影系数获得，当射线在X轴投影较长，选择XOZ平面；反之，当射线在Y轴投影较长，选择YOZ平面.

如图3所示，将射线MN的三维投影系数求解转化为：①求解M′N′与XOY平面相交的投影系数；②求解M″N″与YOZ平面相交的投影系数；③将①和②的结果综合处理.下面对这三个过程分别进行说明.

图3 三维射线投影示意图 Fig.3 3D projection schematic diagram

1）M′N′与XOY投影系数的计算

分析图3可以得出，射线M′N′投影系数的计算与图2中射线MN1的计算完全相同，因此可采用式（4）～式（7）完成.

2）M″N″与YOZ投影系数的计算

分析图3可以得出，射线M″N″与水平网格线l（图3中的点划线）相交的信息即可反映层索引的变化，所以这里仅分析M″N″与l相交的情况.

设M″N″的斜率0＜k′≤1，与Z轴交点集合用数组PZ（PZcy，PZcz）表示，PZcy表示与Y轴平行线相交的第c点的横坐标值，代表交线长度；PZcz表示与Y轴平行线相交的第c点的纵坐标值，代表网格编号.按照式（4）的分析方法可得数组PZ中坐标点的表达式为

式中：c∈［1，2n］，（PZ1y，PX1z）表示入射点M的起始坐标.最后，数组PZ内的元素即为层间索引结果，包含层信息以及对应的长度信息.

3）综合1）和2）结果的分析

射线MN的完整投影系数是将1）结果和2）结果进行综合，其中1）得出了射线M′N′在XOY平面内的网格编号和交线长度；2）得出了射线M″N″在YOZ内的网格编号和交线长度；由于两者都含有相同的横坐标轴Y，因此可利用其将XOY平面内网格编号和YOZ平面内编号相融合，具体融合方法为：

设a，b，c，e，f，h，i为射线M′N′在Y轴上的映射点，其网格编号分别对应为p1，p2，p3，p4，p5，p6，p7，点d，g，j为射线M″N″在Y轴上的映射点，其分别对应的层间网格编号为q1，q2，q3.根据图3中各点在Y轴投影关系可知相邻点间的线段网格编号规则为：

1）当相邻点在Y轴坐标都小于d坐标，则相邻点间线段网格编号为其原有网格编号和d的网格编号相加之和.如：［a，b］间线段网格编号为p1+q1，［b，c］间线段网格编为p2+q1.

2）当相邻点间在Y轴坐标一个小于d坐标，另一个大于d坐标，则表明这两点不在同一层，需要分别进行网格编号.如［c，e］间线段网格编号分为p3+q1和p3+q2.

依据上述规则，可得［e，f］，［f，g］，［g，h］，［h，i］间线段的网格编号分别为：p4+q2，p5+q2，p5+q3，p6+q3.另外，对于这些点的交线长度可由距离公式获得.

考虑到系统设计时，即将射线源、被测物体中心以及探测器中心置于同一平面内的同一直线上，所以探测器每一列对应的射线在XOY平面的层内索引结果是一样的，只是层间索引结果不同，因此，在计算一个投影角度下射线的投影系数时，可以采取列优先的方法，一个列上的投影射线只需计算一次层内索引.而且，利用存在的对称关系［13］，在列优先的基础上，可以只计算Z轴正半轴的射线的投影系数，负半轴方向的射线的投影系数则直接根据对称关系由正半轴得到.由此，又可减少一半的运算量.

3 实验结果与分析

以Sheep＿Logan模型为例进行算法仿真实验，重建图像大小为128×128×128，体素尺寸为0.256 mm，探测器矩阵大小为128×128，像元尺寸0.512 mm.射线源距旋转中心的距离以及探测器中心距旋转中心的距离均设定为780 mm.每绕旋转中心旋转1°采集一幅投影，一共得到360幅投影.采用ART算法进行重建，设定图像初值x（0）=0，松弛因子λ=0.025，顺序访问方式，采用本文提出的投影系数快速计算方法.迭代3次的重建结果切片与原始模型切片做对比，如图4所示.图5为对应第64行重建切片与原始模型切片灰度曲线的对比图.

图4 原始模型与重建结果图 Fig.4 Original model and reconstructed result

图5 切片的灰度对比图 Fig.5 Grayscale contrast of original model and reconstructed result

由图4，图5对比可得，经过3次迭代后，重建图像切片与原始模型切片的灰度曲线相近，说明在ART算法中采用本文提出的投影系数计算方法，重建图像质量能够得到保障.

下面从投影重建时间的角度列表分析，将本文方法与传统的Siddon算法以及张顺利所提的投影系数计算方法进行比较.重建时间对比如表1所示.

表1 不同投影系数据算方法重建时间对比 Tab.1 Reconstruction time contrast of different projection coefficient computation method

从表1可以看出，本文所提方法较传统Siddon算法重建速度提高约13倍；相比张顺利所提的方法，重建速度提高1.26倍.并且由于所提算法中分支判断大为减少，所以后期在进行CUDA算法加速时，有利于并行加速的实现［14］.利用存在的对称关系，进行对称优化后，速度可再提高约3倍，加速效果更加明显.

4 结 论

迭代算法是基础的CT重建算法，本文以ART算法为例，介绍了ART算法的基本原理，其中，投影系数的计算是重建速度和重建质量的最大影响因素.本文分别针对二维情况和三维情况下，以Siddon算法和张顺利所提算法为基础，提出投影系数的快速计算方法，从射线穿过网格的相交规律出发，通过顺序增量计算，快速推导出穿过的网格编号并计算其交线长度，极大地减少了运算量，分支判断也大为减少.仿真实验证明本文所提方法在保证重建图像重建质量的基础上，算法速度大幅提高，相比其他常用算法有明显的速度优势.

随着迭代算法重建速度的提高，鉴于其在不完备投影数据情况下的优异表现，必定会在实时的工业成像检测系统中发挥重要的作用.下一步的工作应重点将放在对算法的硬件加速改进.

［1］Li J，Li LT，Cong P，et al.Rotating polar-coordinate ART applied in industrial CT image reconstruction［J］.NDTffE International，2007，40（4）：333-336.

［2］张朝宗，郭志平，张鹏，等.工业CT技术与原理［M］.北京：科学出版社，2009.

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

［4］冷骏.ART算法中关于松弛因子的研究［J］.光学仪器，2014，36（1）：46-51.Leng Jun.Research on the relaxation factor in algebraic reconstruction technique［J］.Optical Instruments，2014，36（1）：46-51.（in Chinese）

［5］Yang Juan，Hou Huiling，Shi Lang.A method based on vector type to sparsely store and quickly access projection matrix［J］.Journal of Measurement Science and Instrumentation，2015，6（1）：53-56.

［6］Xue Z，Zhang L L.A new algorithm for calculating the radiological path in CT image reconstruction［C］.International conference on Electronic ff Mechanical Engineering and Information Technology，2011：4527-4530.

［7］Miao C.Comparative studies of different system models for iterative CT image reconstruction［D］.Master Dissertation.Winston-Salem：Wake Forest University，2013.

［8］Siddon R L.Fast calculation of the exact radiological path for a three-dimensional CT array［J］.Med.Phys.，1985，12：252-255.

［9］张顺利，张定华，赵歆波.代数重建法中的一种快速投影系数计算方法［J］.计算机应用研究，2007，24（5）：38-40.Zhang Shunli，Zhang Dinghua，Zhao Xinbo.Approach for fast projection coefficient computation in algeb reconstruction technique［J］.Application Research of Computers，2007，24（5）：38-40.（in Chinese）

［10］张顺利，张定华，黄魁东，等.锥束ART算法快速图像重建［J］.仪器仪表学报，2009，30（4）：887-892.Zhang Shunli，Zhang Dinghua，Huang Kuidong，et al.Fast image reconstruction with cone-beam ART algorithm［J］.Chinese Journal of Scientific Instrument，2009，30（4）：887-892.（in Chinese）

［11］杨文良，魏东波.一种改进投影系数计算的快速ART算法［J］.CT理论与应用研究，2012，21（2）：187-195.Yang Wenliang，Wei Dongbo.A fast algebraic reconstruction algorithm based on improved projection coefficient computation［J］.CT Theory and Applications，2012，21（2）：187-195.（in Chinese）

［12］Zhang S L，Zhang D H，Gong H，et al.Fast and accurate computation of system matrix for area integral model-based algebraic reconstruction technique［J］.Optical Engineering.2014，53（11）：1-9.

［13］肖波.基于离散化模型对称结构改进的EM算法研究［D］.北京：北京交通大学，2009.

［14］雷德川，许州，陈浩.基于CUDA的GPU加速代数迭代重建算法［J］.无损检测，2012，34（8）：5-9.Lei Dechuan，Xu Zhou，Chen Hao.Accelerating simultaneous algebraic reconstruction technique based on CUDA-enabled GPU［J］.NDT，2012，34（8）：5-9.（in Chinese）