基于二维索引计算投影系数构造系统矩阵的简易方法

摘要

本发明提供一种基于二维索引计算投影系数构造系统矩阵的简易方法，首先构造坐标轴与网格，确定放射射线，射线一般使用直线截距公式表示。然后确定放射射线与网格的位置关系，计算出射线与网格相交的各点的横纵坐标，根据坐标计算每个单独网格中所截射线线段长度。最后确定每条所截线段的网格编号，并且与其线段长度按照一定规则保存，得到最终结果。本发明减少了对放射线位置的要求，并使用二维检索的方式加快其计算速度，通过实验也验证该方法的有效性。在实际运用在有很大的提升空间。

说明书

技术领域

本发明涉及单光子发射计算机断层成像术（Single-Photon Emission Computed Tomography, SPECT）和正电子发射断层成像术（Positron Emission Tomography, PET）等主要使用迭代算法成像的技术，具体涉及基于二维引索计算投影系数构造系统矩阵的简易方法，属于核医学成像技术领域。

背景技术

在SPECT和PET成像中主要运用两种迭代算法：似然期望值最大化算法（Maximum Likelihood Expectation Maximization, MLEM）和有序子集-期望值最大化算法(Ordered Sub-sets Expectation Maximization, OSEM)。这两种算法的关键在于事先求得精确的系统矩阵。系统矩阵反映了所要还原的图像像素对放射射线的贡献，即放射射线在还原图像上的概率分布。

目前有许多方法计算系统矩阵，其中最简单的就是先计算投影系数，然后在把投影系数按一定顺序排列就可以得到系统矩阵。投影系数数值上等于放射射线穿过像素的长度。对于投影系数的计算通常是把一幅图像网格化，每个像素对应网格中的一个小格并进行编号，再计算放射射线穿过小格的长度。编号与长度构成了投影系数。

目前投影系数算法中，网格中的小格都是以一维的方式进行编号，比如网格大小为64×64，其中小格一共有4096(64×64)个，即小格就有4096个编号。并且为了保证计算的精确度，放射线每穿过一个小格都要考虑与小格相交点位置的情况，不同的位置计算方法不同，这就导致算法运算量大，计算效率降低。

发明内容

本发明的目的在于在保持计算的精确度不变的基础下，简化投影系数计算，提高系统矩阵的计算效率，提供一种基于二维索引计算投影系数构造系统矩阵的简易方法，通过下列技术方案实现。

一种基于二维索引计算投影系数构造系统矩阵的简易方法，包括下列步骤：

首先构造坐标轴与网格，确定放射射线，射线一般使用直线截距公式表示；然后确定放射射线与网格的位置关系，计算出射线与网格相交的各点的横纵坐标，根据坐标计算每个单独网格中所截射线线段长度；最后确定每条所截线段的网格编号，并且与其线段长度按照一定规则保存，得到最终结果。

具体步骤如下：

（1）根据原始正弦图的尺寸构造网格和坐标轴；

（2）用直线斜截式方程y=kx+b表示放射射线，并且计算出射线与网格各个相交点的坐标；

（3）根据射线斜率，计算出网格所截射线上线段的长度；

（4）根据斜率情况，按照下列公式计算行和列编号：

N表示网格大小，一般为偶数，x为相交点的横坐标点，y为相交点的纵坐标点，计算列编号使用公式：                                                ，行编号有两种情况，斜率大于零时：；斜率小于零时：，当斜率不存在时只计算列编号：；当斜率等于零时只计算行编号：；

（5）线段长度按照计算所得编号保存；

（6）把步骤（5）得到的矩阵转换为列或行向量，得到最终的系统矩阵。

本发明的原理：在计算线段长度时本方法运用经典的欧拉公式。计算出射线与网格交点坐标后，根据射线的斜率不同情况对坐标点保存。当斜率大于零，横纵坐标都按照从小到大的顺序保存；当斜率小于零，横坐标按照从小到大的顺序保存，而纵坐标按照从大到小的顺序保存。对于斜率的另外两种特殊情况：不存在和等于零，没有必要再保存所有的相交坐标点，此时网格所截线段长度都为网格的大小（一般设为1）。线段长度总数要比坐标总数少一个，如图2所示。二维的网格索引方式，就是用行和列对每个小格编号，类似于矩阵中元素的编号方式。N表示网格大小，一般为偶数，x为横坐标点，y为纵坐标点，计算列编号使用公式：，行编号有两种情况，斜率大于零时：；斜率小于零时：。当斜率不存在时只考虑列编号：；当斜率等于零时这考虑行编号：。

本方法的特点就在于二维，例如网格尺寸为64×64，第一维把网格中的小格编为64个，第二维同样编为64个，就相当于矩阵中的行号和列号，每个小格有两个数值的编号表示，处理64个编号比处理4096个编号更简单，并且只考虑放射线斜率情况，降低的计算的复杂性，提高了计算效率。

本发明克服了现有投影系数计算效率低、算法复杂等缺点。二维索引编号方式不但可以简化系统矩阵的计算，提高计算效率，同时也符合一般图像的二维存储方式，易于确定图像像素点的位置。这种编号方式非常直观，并且便于计算，行编号和列编号可以单独计算，大大提高计算效率。

本发明可在CT、SPECT和PET等断层成像中都得到广泛运用。当今CT、SPECT和PET采集到的数据即原始数据都是以正弦图的形式展现出来，所谓重建图像就是指把正弦图转换为正常图像的过程，也可以称为把数据转换为可视图像的过程。迭代重建法是当今断层图像重建的主流方法，其中最主要的参数就是系统矩阵，系统矩阵直接影响重建图像的质量和速度。精确快速的计算系统矩阵，对整个图像重建过程具有至关重要的作用。

与现有算法相比，本发明减少了对放射线位置的要求，并使用二维检索的方式加快其计算速度，通过实验也验证该方法的有效性。在实际运用在有很大的提升空间。

附图说明

图1为某一原始数据正弦图；

图2为实施例1的射线L1与网格的相交图；

图3为实施例1的网格所截射线上线段的长度储存示意图；

图4为实施例1计算行和列的编号示意图；

图5为实施例1线段长度编号储存示意图；

图6为实施例1得到最终的系统矩阵示意图；

图7为实施例2的射线L2与网格的相交图；

图8为实施例2的网格所截射线上线段的长度储存示意图；

图9为实施例2计算行和列的编号示意图；

图10为实施例2线段长度编号储存示意图；

图11为实施例2得到最终的系统矩阵示意图；

图12为实施例3的射线L3与网格的相交图；

图13为实施例3计算行和列的编号示意图；

图14为实施例3线段长度编号储存示意图；

图15为实施例3得到最终的系统矩阵示意图；

图16为实施例4的射线L4与网格的相交图；

图17为实施例4计算行和列的编号示意图；

图18为实施例4线段长度编号储存示意图；

图19为实施例4得到最终的系统矩阵示意图；

图20为射线L1～L4的运动轨迹图；

图21为经本发明方法计算后的重建图像。

具体实施方式

下面结合实施例和附图对本发明做进一步说明。

实施例1

（1）根据原始正弦图的尺寸构造网格和坐标轴；如图1为原始数据正弦图，正弦图的大小为180×4，其中180表示探测器扫描的角度一般为固定值，4表示所要构造网格的大小，所要构造的网格尺寸为4×4，即N=4；

（2）用直线斜截式方程y=kx+b表示放射射线，并且计算出射线与网格各个相交点的坐标；如图2的射线L₁，已知射线角度为45°，所以斜率k=1；网格中心为坐标原点，L₁穿过原点所以b=0，得到直线方程为y=x；由于网格的大小已知，这里设为1，根据y=x可以得到交点坐标A(-2,-2)，B(-1,-1)，C(0,0)，D(1,1)，E(2,2)；

（3）根据射线斜率，使用欧拉公式计算计算出网格所截射线上线段的长度；射线L₁的斜率k=1（大于0），储存示意图如图3；

（4）根据斜率情况，按照下列公式计算行和列编号：

N表示网格大小，一般为偶数，x为横坐标点，y为纵坐标点，计算列编号使用公式：，行编号使用公式：；如图4所示；

（5）线段长度按照计算所得编号保存；这里的编号类似于矩阵中单个元素的位置，L₁储存示意图如图5，空白位置值为0：

（6）把步骤（5）得到的矩阵转换为列或行向量，得到最终的系统矩阵，如图6。

实施例2

（1）根据原始正弦图的尺寸构造网格和坐标轴；如图1为原始数据正弦图，正弦图的大小为180×4，其中180表示探测器扫描的角度一般为固定值，4表示所要构造网格的大小，所要构造的网格尺寸为4×4，即N=4；

（2）用直线斜截式方程y=kx+b表示放射射线，并且计算出射线与网格各个相交点的坐标；图7的射线L₂，已知射线角度为135°，所以斜率k=-1；网格中心为坐标原点，L₂穿过原点所以b=0，得到直线方程为y=-x；由于网格的大小已知，这里设为1，根据y=-x可以得到交点坐标A(-2,2)，B(-1,1)，C(0,0)，D(1,-1)，E(2,-2)；

（3）根据射线斜率，使用欧拉公式计算计算出网格所截射线上线段的长度；射线L₂的斜率k=-1（小于0），储存示意图如图8所示：

（4）根据斜率情况，按照下列公式计算行和列编号：

N表示网格大小，一般为偶数，x为相交点的横坐标点，y为相交点的纵坐标点，计算列编号使用公式：，行编号使用公式：斜率小于零时：；如图9所示；

（5）线段长度按照计算所得编号保存；这里的编号类似于矩阵中单个元素的位置， L₂储存示意图如图10所示，空白位置值为0；

（6）把步骤（5）得到的矩阵转换为列或行向量，得到最终的系统矩阵，如图11。

实施例3

（1）根据原始正弦图的尺寸构造网格和坐标轴；如图1为原始数据正弦图，正弦图的大小为180×4，其中180表示探测器扫描的角度一般为固定值，4表示所要构造网格的大小，所要构造的网格尺寸为4×4，即N=4；

（2）用直线斜截式方程y=kx+b表示放射射线，并且计算出射线与网格各个相交点的坐标；图12中的射线L₃，已知射线角度为90°，所以斜率不存在；网格中心为坐标原点，L₂穿过原点所以b=0，得到直线方程为x=0.5；由于网格的大小已知，这里设为1；

（3）根据射线斜率，使用欧拉公式计算计算出网格所截射线上线段的长度；射线L₃的斜率不存在；

（4）根据斜率情况，按照下列公式计算行和列编号：

N表示网格大小，一般为偶数，x为相交点的横坐标点，y为相交点的纵坐标点，计算列编号使用公式：，行编号使用公式：；如图13所示；

（5）线段长度按照计算所得编号保存；这里的编号类似于矩阵中单个元素的位置，L₃储存示意图如图14所示，空白位置值为0；

（6）把步骤（5）得到的矩阵转换为列或行向量，得到最终的系统矩阵，如图15所示。

实施例4

（1）根据原始正弦图的尺寸构造网格和坐标轴；如图1为原始数据正弦图，正弦图的大小为180×4，其中180表示探测器扫描的角度一般为固定值，4表示所要构造网格的大小，所要构造的网格尺寸为4×4，即N=4；

（2）用直线斜截式方程y=kx+b表示放射射线，并且计算出射线与网格各个相交点的坐标；图16的射线L₄，已知射线角度为180°，所以斜率k=0；网格中心为坐标原点，L₂穿过原点，得到直线方程为y=-0.5；由于网格的大小已知，这里设为1；

（3）根据射线斜率，使用欧拉公式计算计算出网格所截射线上线段的长度；射线L₄的斜率k=0；

（4）根据斜率情况，按照下列公式计算行和列编号：

N表示网格大小，一般为偶数，x为相交点的横坐标点，y为相交点的纵坐标点，计算列编号使用公式：，行编号使用公式：；如图17所示；

（5）线段长度按照计算所得编号保存；这里的编号类似于矩阵中单个元素的位置，L₃储存示意图如图18所示，空白位置值为0；

（6）把步骤（5）得到的矩阵转换为列或行向量，得到最终的系统矩阵，如图19所示。

射线的数目与网格大小相对应，如图20中网格为4×4，相应的射线有L₁~ L₄ 四条，其起始位置如图2、7、12、16中4条实线的位置。射线的运动轨迹是以网格中心为轴心，逆时针方向旋转179°，图20中的虚线是射线L₁~ L₄旋转90°后所在位置。每条射线每转动一个角度就要按照以上步骤计算，最后把行向量合并为一个大的矩阵就为所要求的系统矩阵。在MLEM、OSEM等迭代重建算法中使用系统矩阵就可以进行图像重建，重建后的图像如图21所示。