CN102682199A - 基于二维索引计算投影系数构造***矩阵的简易方法 - Google Patents

Abstract

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

Description

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

技术领域

本发明涉及单光子发射计算机断层成像术（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所示。

Claims (2)

1.一种基于二维索引计算投影系数构造***矩阵的简易方法，其特征在于包括下列步骤：

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

2.根据权利要求1所示的基于二维索引计算投影系数构造***矩阵的简易方法，其特征在于具体步骤如下：

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

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

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

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

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

，行编号有两种情况，斜率大于零时：；斜率小于零时：

，当斜率不存在时只计算列编号：

；当斜率等于零时只计算行编号：

；

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

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

Images