CN102682199A - 基于二维索引计算投影系数构造***矩阵的简易方法 - Google Patents
基于二维索引计算投影系数构造***矩阵的简易方法 Download PDFInfo
- Publication number
- CN102682199A CN102682199A CN2012101287048A CN201210128704A CN102682199A CN 102682199 A CN102682199 A CN 102682199A CN 2012101287048 A CN2012101287048 A CN 2012101287048A CN 201210128704 A CN201210128704 A CN 201210128704A CN 102682199 A CN102682199 A CN 102682199A
- Authority
- CN
- China
- Prior art keywords
- ray
- grid
- slope
- system matrix
- calculate
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 239000011159 matrix material Substances 0.000 title claims abstract description 41
- 238000000034 method Methods 0.000 title claims abstract description 21
- 230000005855 radiation Effects 0.000 claims abstract description 9
- 238000004364 calculation method Methods 0.000 claims abstract description 4
- 238000005304 joining Methods 0.000 claims description 17
- 230000014509 gene expression Effects 0.000 claims description 15
- 238000004513 sizing Methods 0.000 claims description 9
- 239000013598 vector Substances 0.000 claims description 7
- 238000006243 chemical reaction Methods 0.000 claims description 6
- 238000002474 experimental method Methods 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 20
- 238000002603 single-photon emission computed tomography Methods 0.000 description 6
- 238000003860 storage Methods 0.000 description 6
- 238000012636 positron electron tomography Methods 0.000 description 3
- 230000002285 radioactive effect Effects 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 2
- 238000012831 peritoneal equilibrium test Methods 0.000 description 2
- 238000002600 positron emission tomography Methods 0.000 description 2
- 238000012877 positron emission topography Methods 0.000 description 2
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 238000012879 PET imaging Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000009206 nuclear medicine Methods 0.000 description 1
- 239000010902 straw Substances 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Landscapes
- Apparatus For Radiation Diagnosis (AREA)
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 1,已知射线角度为45°,所以斜率k=1;网格中心为坐标原点,L 1穿过原点所以b=0,得到直线方程为y=x;由于网格的大小已知,这里设为1,根据y=x可以得到交点坐标A(-2,-2),B(-1,-1),C(0,0),D(1,1),E(2,2);
(3)根据射线斜率,使用欧拉公式计算计算出网格所截射线上线段的长度;射线L 1的斜率k=1(大于0),储存示意图如图3;
(4)根据斜率情况,按照下列公式计算行和列编号:
(5)线段长度按照计算所得编号保存;这里的编号类似于矩阵中单个元素的位置,L 1储存示意图如图5,空白位置值为0:
(6)把步骤(5)得到的矩阵转换为列或行向量,得到最终的***矩阵,如图6。
实施例2
(1)根据原始正弦图的尺寸构造网格和坐标轴;如图1为原始数据正弦图,正弦图的大小为180×4,其中180表示探测器扫描的角度一般为固定值,4表示所要构造网格的大小,所要构造的网格尺寸为4×4,即N=4;
(2)用直线斜截式方程y=kx+b表示放射射线,并且计算出射线与网格各个相交点的坐标;图7的射线L 2,已知射线角度为135°,所以斜率k=-1;网格中心为坐标原点,L 2穿过原点所以b=0,得到直线方程为y=-x;由于网格的大小已知,这里设为1,根据y=-x可以得到交点坐标A(-2,2),B(-1,1),C(0,0),D(1,-1),E(2,-2);
(3)根据射线斜率,使用欧拉公式计算计算出网格所截射线上线段的长度;射线L 2的斜率k=-1(小于0),储存示意图如图8所示:
(4)根据斜率情况,按照下列公式计算行和列编号:
N表示网格大小,一般为偶数,x为相交点的横坐标点,y为相交点的纵坐标点,计算列编号使用公式:,行编号使用公式:斜率小于零时:;如图9所示;
(5)线段长度按照计算所得编号保存;这里的编号类似于矩阵中单个元素的位置, L 2储存示意图如图10所示,空白位置值为0;
(6)把步骤(5)得到的矩阵转换为列或行向量,得到最终的***矩阵,如图11。
实施例3
(1)根据原始正弦图的尺寸构造网格和坐标轴;如图1为原始数据正弦图,正弦图的大小为180×4,其中180表示探测器扫描的角度一般为固定值,4表示所要构造网格的大小,所要构造的网格尺寸为4×4,即N=4;
(2)用直线斜截式方程y=kx+b表示放射射线,并且计算出射线与网格各个相交点的坐标;图12中的射线L 3,已知射线角度为90°,所以斜率不存在;网格中心为坐标原点,L 2穿过原点所以b=0,得到直线方程为x=0.5;由于网格的大小已知,这里设为1;
(3)根据射线斜率,使用欧拉公式计算计算出网格所截射线上线段的长度;射线L 3的斜率不存在;
(4)根据斜率情况,按照下列公式计算行和列编号:
(5)线段长度按照计算所得编号保存;这里的编号类似于矩阵中单个元素的位置,L 3储存示意图如图14所示,空白位置值为0;
(6)把步骤(5)得到的矩阵转换为列或行向量,得到最终的***矩阵,如图15所示。
实施例4
(1)根据原始正弦图的尺寸构造网格和坐标轴;如图1为原始数据正弦图,正弦图的大小为180×4,其中180表示探测器扫描的角度一般为固定值,4表示所要构造网格的大小,所要构造的网格尺寸为4×4,即N=4;
(2)用直线斜截式方程y=kx+b表示放射射线,并且计算出射线与网格各个相交点的坐标;图16的射线L 4,已知射线角度为180°,所以斜率k=0;网格中心为坐标原点,L 2穿过原点,得到直线方程为y=-0.5;由于网格的大小已知,这里设为1;
(3)根据射线斜率,使用欧拉公式计算计算出网格所截射线上线段的长度;射线L 4的斜率k=0;
(4)根据斜率情况,按照下列公式计算行和列编号:
(5)线段长度按照计算所得编号保存;这里的编号类似于矩阵中单个元素的位置,L 3储存示意图如图18所示,空白位置值为0;
(6)把步骤(5)得到的矩阵转换为列或行向量,得到最终的***矩阵,如图19所示。
射线的数目与网格大小相对应,如图20中网格为4×4,相应的射线有L 1~ L 4 四条,其起始位置如图2、7、12、16中4条实线的位置。射线的运动轨迹是以网格中心为轴心,逆时针方向旋转179°,图20中的虚线是射线L 1~ L 4旋转90°后所在位置。每条射线每转动一个角度就要按照以上步骤计算,最后把行向量合并为一个大的矩阵就为所要求的***矩阵。在MLEM、OSEM等迭代重建算法中使用***矩阵就可以进行图像重建,重建后的图像如图21所示。
Claims (2)
1.一种基于二维索引计算投影系数构造***矩阵的简易方法,其特征在于包括下列步骤:
首先构造坐标轴与网格,确定放射射线,射线一般使用直线截距公式表示;然后确定放射射线与网格的位置关系,计算出射线与网格相交的各点的横纵坐标,根据坐标计算每个单独网格中所截射线线段长度;最后确定每条所截线段的网格编号,并且与其线段长度按照一定规则保存,得到最终结果。
2.根据权利要求1所示的基于二维索引计算投影系数构造***矩阵的简易方法,其特征在于具体步骤如下:
(1)根据原始正弦图的尺寸构造网格和坐标轴;
(2)用直线斜截式方程y=kx+b表示放射射线,并且计算出射线与网格各个相交点的坐标;
(3)根据射线斜率,计算出网格所截射线上线段的长度;
(4)根据斜率情况,按照下列公式计算行和列编号:
N表示网格大小,一般为偶数,x为相交点的横坐标点,y为相交点的纵坐标点,计算列编号使用公式: ,行编号有两种情况,斜率大于零时:;斜率小于零时:,当斜率不存在时只计算列编号:;当斜率等于零时只计算行编号:;
(5)线段长度按照计算所得编号保存;
(6)把步骤(5)得到的矩阵转换为列或行向量,得到最终的***矩阵。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210128704.8A CN102682199B (zh) | 2012-04-28 | 2012-04-28 | 基于二维索引计算投影系数构造***矩阵的简易方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210128704.8A CN102682199B (zh) | 2012-04-28 | 2012-04-28 | 基于二维索引计算投影系数构造***矩阵的简易方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN102682199A true CN102682199A (zh) | 2012-09-19 |
CN102682199B CN102682199B (zh) | 2015-08-26 |
Family
ID=46814115
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201210128704.8A Expired - Fee Related CN102682199B (zh) | 2012-04-28 | 2012-04-28 | 基于二维索引计算投影系数构造***矩阵的简易方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN102682199B (zh) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107730579A (zh) * | 2016-08-11 | 2018-02-23 | 深圳先进技术研究院 | 一种锥束ct投影矩阵的计算方法及*** |
CN109498048A (zh) * | 2019-01-04 | 2019-03-22 | 南京航空航天大学 | 一种加速正电子图像重建的***矩阵生成与处理方法 |
CN109783771A (zh) * | 2019-01-22 | 2019-05-21 | 清华大学 | 将轨迹序列转换为图像矩阵的处理方法、装置和存储介质 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CA2324832A1 (en) * | 2000-10-31 | 2002-04-30 | Roger Lecomte | Real-time image reconstruction for computed tomography systems |
CN1839758A (zh) * | 2005-02-15 | 2006-10-04 | 西门子公司 | 用圆-线快速扫描算法重构计算机断层分析图像的方法 |
CN101034479A (zh) * | 2006-03-10 | 2007-09-12 | Ge医疗***环球技术有限公司 | 图像重建方法和x射线ct设备 |
-
2012
- 2012-04-28 CN CN201210128704.8A patent/CN102682199B/zh not_active Expired - Fee Related
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CA2324832A1 (en) * | 2000-10-31 | 2002-04-30 | Roger Lecomte | Real-time image reconstruction for computed tomography systems |
CN1839758A (zh) * | 2005-02-15 | 2006-10-04 | 西门子公司 | 用圆-线快速扫描算法重构计算机断层分析图像的方法 |
CN101034479A (zh) * | 2006-03-10 | 2007-09-12 | Ge医疗***环球技术有限公司 | 图像重建方法和x射线ct设备 |
Non-Patent Citations (4)
Title |
---|
ROBERT L. SIDDON: "Fast calculation of the exact radiological path for a threedimensional CT array", 《MEDICAL PHYSICS》, vol. 12, no. 2, 31 December 1985 (1985-12-31), pages 252 - 255 * |
周斌等: "代数重建算法中一种快速投影系数计算方法", 《计算机工程与应用》, vol. 44, no. 25, 31 December 2008 (2008-12-31), pages 46 - 47 * |
王旭等: "联合代数重建算法中基于像素的投影计算方法", 《核电子学与探测技术》, vol. 25, no. 6, 30 November 2005 (2005-11-30), pages 785 - 788 * |
胡小舟等: "一种基于POCS约束的图像代数重建算法", 《模式识别与人工智能》, vol. 22, no. 5, 31 October 2009 (2009-10-31), pages 763 - 767 * |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107730579A (zh) * | 2016-08-11 | 2018-02-23 | 深圳先进技术研究院 | 一种锥束ct投影矩阵的计算方法及*** |
CN107730579B (zh) * | 2016-08-11 | 2021-08-24 | 深圳先进技术研究院 | 一种锥束ct投影矩阵的计算方法及*** |
CN109498048A (zh) * | 2019-01-04 | 2019-03-22 | 南京航空航天大学 | 一种加速正电子图像重建的***矩阵生成与处理方法 |
CN109783771A (zh) * | 2019-01-22 | 2019-05-21 | 清华大学 | 将轨迹序列转换为图像矩阵的处理方法、装置和存储介质 |
Also Published As
Publication number | Publication date |
---|---|
CN102682199B (zh) | 2015-08-26 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
EP1989683B1 (en) | Image reconstruction using data ordering | |
CN110097611A (zh) | 图像重建方法、装置、设备及存储介质 | |
CN103732147A (zh) | X射线计算机断层摄影装置以及图像重建方法 | |
CN101908231A (zh) | 处理含有主平面场景的三维点云重建方法和*** | |
Zhang et al. | Fast and memory‐efficient Monte Carlo‐based image reconstruction for whole‐body PET | |
CN104408756B (zh) | 一种pet图像重建方法及装置 | |
US7332721B2 (en) | Separation of geometric system response matrix for three-dimensional image reconstruction | |
Asharindavida et al. | Study on hexagonal grid in image processing | |
CN102682199B (zh) | 基于二维索引计算投影系数构造***矩阵的简易方法 | |
CN109498048A (zh) | 一种加速正电子图像重建的***矩阵生成与处理方法 | |
US20130181989A1 (en) | Efficiently Reconstructing Three-Dimensional Structure and Camera Parameters from Images | |
CN106859686B (zh) | 成像方法和成像*** | |
Álvarez-Gómez et al. | Fast energy dependent scatter correction for list-mode PET data | |
GB2505998A (en) | Conversion of X-Ray intensity distribution data | |
JP6014738B2 (ja) | 三次元画像の投影方法 | |
KR101128566B1 (ko) | Gpu를 이용한 양전자 방출 단층 촬영 영상에서의 감마선 산란 추정 방법 및 장치 | |
Cecchetti et al. | Accurate and efficient modeling of the detector response in small animal multi-head PET systems | |
Alhassen et al. | Ultrafast multipinhole single photon emission computed tomography iterative reconstruction using CUDA | |
Nassiri et al. | Fast GPU-based computation of spatial multigrid multiframe LMEM for PET | |
Zhang et al. | Characterization of the ringing artifacts in rotator‐based reconstruction with Monte Carlo‐based resolution compensation for PET | |
WO2023228910A1 (ja) | 画像処理装置および画像処理方法 | |
Dimmock et al. | An OpenCL implementation of pinhole image reconstruction | |
CN112381903B (zh) | 一种流水线式tof正电子图像快速重建方法 | |
Wang et al. | A parallel method for texture reconstruction in large-scale 3D automatic modeling based on oblique photography | |
Musmann et al. | List-mode image reconstruction for the high resolution ClearPET/spl trade/Neuro system |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20150826 Termination date: 20210428 |
|
CF01 | Termination of patent right due to non-payment of annual fee |