发明内容
本发明要解决的问题是提供一种PET二维图像重建方法和装置,对现有的OSEM算法结构进行优化,使用该方法和装置重建PET图像,能够显著降低运算时间。
为解决上述问题,本发明提供了一种PET二维图像重建方法,包括以下步骤:
1)为需要重建的一层创建一个二维初始图像,所述初始图像位于第一坐标系(X,Y)中的各像素点的像素值为f1(x,y);
2)将所述需要重建的一层所有的响应线划分为多个子集,每一个所述子集包含多个响应线组,每一个响应线组包含多条相互平行的响应线;
3)从多个子集中选择一个子集,对所选择的子集中的每个响应线组进行预设操作,并根据所选择的子集的多个响应线组对应的所述预设操作的结果更新所述初始图像,其中所述预设操作包括:
3.1)根据第一坐标系(X,Y)中的像素点的像素值f1(x,y),计算位于第二坐标系(S,T)中的第一像素组的每个像素点的第一像素值g1(s,t),所述第一像素组的像素点位于所选择的响应线组对应的多条响应线上;
3.2)根据得到的所述第一像素组的每个像素点的第一像素值g1(s,t)进行正投影计算,以获得所述第一像素组的各投影点的理论投影值p1(s,t);
3.3)将所述第一像素组的各投影点的理论投影值p1(s,t)与实际测量投影值p0(s,t)进行对比,以获得校正系数R(s,t),并根据所述校正系数R(s,t)进行反投影计算,以获得所述第一像素组的每个像素点的第二像素值g2(s,t);
3.4)根据位于第二坐标系(S,T)中的第一像素组的每个像素点的第二像素值g2(s,t),计算位于第一坐标系(X,Y)中的第二像素组的每个像素点的像素值f2(x,y),所述第二像素组的每个像素点分布在所选择的响应线组对应的多条响应线上;
3.5)对所述第二像素组的每个像素点的像素值f2(x,y)做归一化处理;
4)以更新后的图像作为初始图像,继续执行步骤3),直至满足停步规则。
进一步地,所述步骤3.1)中,第二坐标系(S,T)的t轴的方向与所选择的响应线组对应的响应线方向平行。
进一步地,所述步骤3.2)中的正投影计算包括:在第一预设范围内将与投影点相关的像素点的像素值进行加权求和;所述步骤3.3)中的反投影计算包括:在第二预设范围内将与像素点相关的投影点的投影值进行加权求和。
进一步地,所述步骤3.2)中,进行正投影计算之前,还包括:将第一像素组中沿t轴方向分布的m个像素点进行像素合并,其中m为大于1的整数;所述步骤3.3)中,实际测量投影值的获取方法还包括:在第二坐标系中将沿t轴方向分布的m个投影点的实际测量投影值进行合并。
进一步地,所述步骤3.1)中,计算所述第一像素组的每个像素点的第一像素值g1(s,t)的步骤包括:
a1)对于每个所述第一像素组的像素点(s1,t1),在第一坐标系中寻找与之对应的像素点的坐标(x1,y1);
a2)在第一坐标系中,通过与所述像素点(x1,y1)相邻的像素点的像素值获得所述像素点(x1,y1)的像素值,所述第一像素组的像素点(s1,t1)的第一像素值g1(s1,t1)与位于第一坐标系中的所述像素点(x1,y1)的像素值相同。
进一步地,所述步骤a1)包括:根据第一坐标系旋转至第二坐标系的角度,计算与每个所述第一像素组的像素点(s1,t1)相对应的位于第一坐标系中的像素点的坐标(x1,y1),计算方法为:
进一步地,所述步骤a2)中,通过与所述像素点(x1,y1)相邻的像素点的像素值获得所述像素点(x1,y1)的像素值的方法包括:通过与所述像素点(x1,y1)相邻的多个像素点进行插值,获得所述像素点(x1,y1)的像素值。
进一步地,所述步骤3.4)中,计算所述第二像素组的每个像素点的像素值f2(x,y)的方法包括:
b1)对于每个所述第二像素组的像素点(x2,y2),在第二坐标系中寻找与之对应的像素点的坐标(s2,t2);
b2)在第二坐标系中,通过与所述像素点(s2,t2)相邻的第一像素组的像素点的第二像素值获得所述像素点(s2,t2)的像素值,所述第二像素组的像素点(x2,y2)的像素值f2(x2,y2)与位于第二坐标系中的所述像素点(s2,t2)的像素值相同。
进一步地,所述步骤b1)包括:根据第二坐标系旋转至第一坐标系的角度,计算与每个所述第二像素组的像素点(x2,y2)相对应的位于第二坐标系中的像素点的坐标(s2,t2),计算方法为:
进一步地,所述步骤b2)中,通过与所述像素点(s2,t2)相邻的第一像素组的像素点的第二像素值获得所述像素点(s2,t2)的像素值的方法包括:通过与所述像素点(s2,t2)相邻的多个像素点进行插值,获得所述像素点(s2,t2)的像素值。
进一步地,所述步骤4)中,停步规则为:所述多个子集都已被选择过或所述初始图像更新后与更新前的像素值变化小于阈值。
为解决上述问题,本发明还提供了一种PET二维图像重建装置,包括:
初始化单元,用于为需要重建的一层创建一个二维初始图像,所述初始图像位于第一坐标系(X,Y)中的各像素点的像素值为f1(x,y);
子集划分单元,用于将所述需要重建的一层所有的响应线划分为多个子集,每一个所述子集包含多个响应线组,每一个响应线组包含多条相互平行的响应线;
图像更新单元,用于从多个子集中选择一个子集,对所选择的子集中的每个响应线组进行预设操作,并根据所选择的子集的多个响应线组对应的所述预设操作的结果更新所述初始图像,所述响应线组处理单元包括:
第一坐标系转换子单元,用于根据第一坐标系(X,Y)中的像素点的像素值f1(x,y),计算位于第二坐标系(S,T)中的第一像素组的每个像素点的第一像素值g1(s,t),所述第一像素组的像素点位于所选择的响应线组对应的多条响应线上;
正投影子单元,用于根据得到的所述第一像素组的每个像素点的第一像素值g1(s,t)进行正投影计算,以获得所述第一像素组的各投影点的理论投影值p1(s,t);
反投影子单元,用于将所述第一像素组的各投影点的理论投影值p1(s,t)与实际测量投影值p0(s,t)进行对比,以获得校正系数R(s,t),并根据所述校正系数R(s,t)进行反投影计算,以获得所述第一像素组的每个像素点的第二像素值g2(s,t);
第二坐标系转换子单元,用于根据位于第二坐标系(S,T)中的第一像素组的每个像素点的第二像素值g2(s,t),计算位于第一坐标系(X,Y)中的第二像素组的每个像素点的像素值f2(x,y),所述第二像素组的每个像素点分布在所选择的响应线组对应的多条响应线上;
归一化处理单元,用于对所述第二像素组的每个像素点的像素值f2(x,y)做归一化处理;
控制单元,用于以更新后的图像作为初始图像,继续执行步骤3),直至满足停步规则。
本发明技术方案通过图像旋转以避免反复计算每条响应线与图像相交的像素坐标和加权系数;以及在利用飞行时间技术信息时,将多个小像素合并为一个大像素进行计算,减少了运算量,实现了OSEM算法运算结构的优化,显著降低了运算时间。
具体实施方式:
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合附图对本发明的具体实施方式做详细的说明。
图1是一种PET结构示意图。请参考图1,所示PET由封闭多环型探测器(置于机架内),电子学系统(前端放大与符合系统)、主控计算机以及检测床构成。其中多环探测器多是由闪烁体(BGO或LYSO)细条、光电倍增管(PMT)和环间隔板(septa)所组成,它的功能是探测在同一环内正电子湮灭时转换成的一对伽马光子所分别命中的环上晶体条的位置,并把这些位置信号转换成电信号,连带伽马光子的能量信号和到达时刻的时间信息一起送到后续的电子学系统中去。电子学系统的功能是确定符合,即是判定一对伽马光子是否由一次湮灭事例所发出的。此后就把经选出的真实的符合事例所命中的两个基本点探测器条的坐标经计算机接口,送到后面的主控计算机去。主控计算机及相应的各个软件包的功能是完成数据采集、系统监控与校正、图像重建和图像处理,并实现临床上各种扫描操作和诊断的要求。
图2是本发明实施例中PET二维图像重建方法的流程图。请参考图2,利用本发明的方法重建PET二维图像,首先执行步骤S110,为需要重建的一层创建一个二维初始图像,所述初始图像位于第一坐标系(X,Y)中的各像素点的像素值为f1(x,y)。
请参考图3,以0点为圆心的圆是PET探测器的截面,所述探测器的轴向方向垂直于所述截面,所述第一坐标系(X,Y)的x轴和y轴确定的平面与所述探测器的截面平行。为需要重建的一层创建一个二维初始图像,所述初始图像对应于位于第一坐标系(X,Y)中的一个二维矩阵,该二维矩阵可以是一个包含N*N个像素点的二维点阵,为了简化说明,图3中只是展示出了该点阵中的部分像素点。
继续执行步骤S120,将所述需要重建的一层所有的响应线划分为多个子集,每一个所述子集包含多个响应线组,每一个响应线组包含多条相互平行的响应线。
采集正电子发射断层扫描飞行时间响应线数据,飞行时间技术记录同一个符合事件中两个伽马光子被探测器接收的时间差,从而对符合事件的定位更准确。在数据采集过程中,探测器会获得大量的响应线数据,图4中仅显示5条响应线作为简单的示意性说明。
我们将所述需要重建的一层所有的响应线划分成不同的响应线子集,其中每个响应线子集包含多个响应线组,每个响应线组包含多条相互平行的响应线。
例如,其中,一个响应线子集包含12个响应线组,所述12个响应线组对应的响应线延伸方向与y轴的夹角分别为0度、30度、60度……300度和330度。
又例如,其中,另一个响应线子集包含9个响应线组,所述9个响应线组对应的响应线延伸方向与y轴的夹角分别为10度、50度、90度……290度和330度。
如前所述,每个响应线组包含多条相互平行的响应线。例如,当一个响应线组对应的响应线方向与y轴的夹角为20度,所述响应线组包含所有与y轴夹角为20度且相互平行的响应线。
例如,在图4中,响应线L1、L2和L3相互平行,那么响应线L1、L2和L3属于同一响应线组,而相互平行的响应线L4和L5上属于另一个响应线组。
继续执行步骤S130,从多个子集中选择一个子集,对所选择的子集中的每个响应线组进行预设操作,并根据所选择的子集的多个响应线组对应的所述预设操作的结果更新所述初始图像。
在本发明的实施例中,每选择一次子集并根据所选择子集进行运算,对投影数据进行校正,所述初始图像就更新一次。
在本发明的实施例中,要获得所选择的子集的图像合并结果,需要依次选择所述选择的子集中多个响应线组分别进行投影校正,根据校正后的投影数据获取多个响应线组对应的校正后的图像信息,并根据多个响应线组对应的校正后的图像信息进行图像信息合并。
下面对所选择的子集中的一个被选中的响应线组对应的预设操作过程进行详细说明。
步骤S131,根据第一坐标系(X,Y)中的像素点的像素值f1(x,y),计算位于第二坐标系(S,T)中的第一像素组的每个像素点的第一像素值g1(s,t),所述第一像素组的像素点位于所选择的响应线组对应的多条响应线上。
图5的左边为第一坐标系(X,Y)中的像素点的示意图,图5的右边为第二坐标系(S,T)中的第一像素组的示意图。
所述二维初始图像中的像素点为图5左边矩形中的小圆点,像素宽度为d,每个像素点对应的像素值为f1(x,y)。
经研究,发明人发现,在图5左边的示意图中,二维初始图像的像素点并不位于所选择的响应线组对应的响应线上,如果通过现有技术计算每条响应线对每个像素的校正系数,将需要消耗较多的时间成本。
为此,发明人提出将所述第一坐标系(X,Y)旋转从而产生第二坐标系(S,T),并根据二维初始图像的像素点的像素值计算出第一像素组的第一像素值,如果响应线L和y轴的夹角为,那么需要将所述第一坐标系旋转的角度为。
如图5所示,所述旋转后的第二坐标系(S,T)的t轴的方向与所选择的响应线组对应的响应线方向平行。所述第一像素组的像素点分布在与所选择响应线组对应的多条平行的响应线L上。所述第一像素组中的像素点的第一像素值为g1(s,t),所述第一像素值g1(s,t)可以根据所述初始图像位于第一坐标系(X,Y)中的像素点的像素值f1(x,y)来获得。
具体地,在本发明的实施例中,可以先获得所述第一像素组的像素点在第一坐标系(X,Y)中与之对应的坐标,然后根据位于所获得的坐标的像素点的像素值来确定所述第一像素组的像素点的第一像素值。
下面以图5中的第一像素组中的像素点A为例进行说明。在第二坐标系(S,T)中,所述像素点A中的坐标为(s1,t1),由第一坐标系(X,Y)旋转至第二坐标系(S,T)的角度就可以计算出所述像素点A在第一坐标系(X,Y)中的坐标(x1,y1),计算方法为:
根据所述坐标(x1,y1)就可以找出第一坐标系(X,Y)中与坐标(x1,y1)相邻的像素点,然后就可以根据所述相邻像素点的像素值确定出位于坐标(x1,y1)的像素点的像素值,像素点A的第一像素值g1(s1,t1)与像素点(x1,y1)的像素值相等。
在本发明的实施例中,可以根据第一坐标系(X,Y)中与坐标(x1,y1)相邻的多个像素点的像素值进行插值确定出所述像素点A在第二坐标系中的第一像素值g1(s1,t1),具体地,可以根据第一坐标系(X,Y)中与坐标(x1,y1)相邻的四个像素点的像素值进行双线性插值确定出所述像素点A在第二坐标系中的第一像素值g1(s1,t1)。
获得了第二坐标系中第一像素组的像素点的第一像素值g1(s,t)后,可以根据所述第一像素值g1(s,t)进行正投影计算。
在本发明的实施例中,为了进一步减少计算量和提高二维图像重建的效率,在所述正投影计算之前,还可以对所述第一像素组中的像素点进行像素合并。具体地,将所选择响应线组对应的所述第一像素组在第二坐标系中沿t轴方向的m个像素点的第一像素值合并,也就是说,将m个像素点合为一个像素点,新合成的像素点的像素值为所述m个像素点的像素值之和。其中,△t为飞行时间信息处理系统中记录飞行时间的时间间隔,d为像素宽度。
经过像素合并后图像大小不变,像素点在s轴方向上的距离不变,在t轴方向上的距离变为md。显然,在通过上述合并操作后,需要进行计算的像素点的数量减少为原来的倍,将较大地提高图像重建效率。
步骤S132,根据得到的所述第一像素组的每个像素点的第一像素值g1(s,t)进行正投影计算,以获得所述第一像素组的各投影点的理论投影值p1(s,t)。
具体地,在本发明的实施例中,所述正投影计算包括:在第一预设范围内将与投影点相关的像素点的像素值进行加权求和,具体为:
其中,g1(s,tj)是第一像素组在坐标(s,tj)的像素点的第一像素值,p1(s,t)是计算得到的第一像素组在坐标(s,t)的理论投影值,w(tj)为时间飞行加权系数,集合J由时间加权系数w(tj)决定,即对于坐标(s,tj)的加权系数w(tj)大于某一阈值时,如0.1,该坐标属于集合J。
对于坐标(s,t)计算与其相邻的坐标(s,t0)的时间加权系数可由归一化高斯函数计算得到,计算式为:
其中,FWHM为高斯函数的半高宽,它的值是系统的时间分辨率与光速的乘积,d为像素宽度,△t为飞行时间信息处理系统中记录飞行时间的时间间隔,d为像素宽度。
步骤S133,将所述第一像素组的各投影点的理论投影值p1(s,t)与实际测量投影值p0(s,t)进行对比,以获得校正系数R(s,t),并根据所述校正系数R(s,t)进行反投影计算,以获得所述第一像素组的每个像素点的第二像素值g2(s,t)。
如果进行正投影计算之前,将第一像素组中沿t轴方向分布的m个像素点进行像素合并,那么相应地可以在第二坐标系中,将沿t轴方向分布的m个投影点的实际测量投影值进行合并。具体地,将m个投影点合并为一个投影点,合并后的投影点的实际测量投影值为所述m个投影点的实际测量投影值的总和,合并后投影点之间的距离增大了,为之前的m倍。
在本发明的实施例中,所述第一像素组的像素点分布在所选择响应线组对应的多条平行的响应线上,所述第一像素组在所选择响应线组对应的第二坐标系中的各投影点的投影值为所述校正系数值R(s,t)。
需要说明的是,所述反投影计算为所述正投影的逆向过程,所述正投影为根据像素值计算投影值,而所述反投影计算为根据投影值计算像素值。具体地,在本发明的实施例中,所述反投影计算可以包括:在第二预设范围内与像素点相关的投影点的投影值进行加权求和。
所述第二预设范围的确定和所述投影点的权重系数的计算方法和上述第一预设范围的确定以及像素点对应的权重系数的计算方法类似,请参考关于上述正投影计算的描述。
步骤S134,根据位于第二坐标系(S,T)中的第一像素组的每个像素点的第二像素值g2(s,t),计算位于第一坐标系(X,Y)中的第二像素组的每个像素点的像素值f2(x,y),所述第二像素组的每个像素点分布在所选择的响应线组对应的多条响应线上。
在本发明的实施例中,所述第二像素组的像素点分布在第一坐标系中的一个二维矩阵中。
需要说明的是,所述步骤S134是步骤S131的一个逆向过程,在步骤S131中是根据所述第一坐标系中的像素点的像素值,来获取第二坐标系中第一像素组的像素点的第一像素值;所述步骤S134是根据第二坐标系中的第一像素组的像素点的第二像素值,来获得第一坐标系中的第二像素组的像素点的像素值。
在本发明的实施例中,可以先获得所述第二像素组的像素点在所述第二坐标系中的坐标,然后根据与所述第二像素组的像素点在所述第二坐标系中的坐标相邻的第一像素组的像素点的第二像素值确定所述第二像素组的像素点的像素值。
关于所述第一坐标系中的第二像素组的像素点的像素值的计算方法的进一步说明,可参考上述关于步骤S131的描述。
步骤S135,对所述第二像素组的每个像素点的像素值f2(x,y)做归一化处理,图像f2(x,y)与图像f1(x,y)像素尺寸不同,需要对f2(x,y)做归一化处理,具体为:。
通过步骤S131至步骤S135将所选择的子集的多个响应线组对应的第二像素组的每个像素点的像素值f2(x,y)都计算出来后,可以将所选择的子集的多个响应线组对应的第二像素组的每个像素点的像素值f2(x,y)进行合并,以获得所选择的响应线子集对应的图像信息。
具体地,在本发明实施例中,还可以根据所述叠加合并后的结果更新所述初始图像的像素值包括:将所述叠加合并后的结果除以归一化的结果的值乘以所述初始图像像素值的结果替换所述初始图像的像素值。其中,所述归一化的结果为将所选择的响应线组对应的所述校正系数R(s,t)置为1后进行所述反投影计算和叠加合并计算后获得的结果。
继续执行步骤S140,以更新后的图像作为初始图像,继续执行步骤S130,直至满足停步规则。
在本发明的实施例中所述停步规则可以为所述多个响应线子集都已被选择过。在本发明的其他实施例中,所述停步规则也可以为所述初始图像更新后与更新前的像素值变化小于阈值,具体为所述叠加合并后的结果除以所述归一化的结果的值小于预设阈值。
当满足停步规则时,输出更新后的所述初始图像的像素值;当未满足停步规则时,返回步骤S130,选择下一响应线子集,继续执行。
本发明的实施例还提供了一种PET二维图像重建装置,请参考图6,所述PET二维图像重建装置600包括:
初始化单元610,用于为需要重建的一层创建一个二维初始图像,所述初始图像位于第一坐标系(X,Y)中的各像素点的像素值为f1(x,y);
子集划分单元620,用于将所述需要重建的一层所有的响应线划分为多个子集,每一个所述子集包含多个响应线组,每一个响应线组包含多条相互平行的响应线;
图像更新单元630,用于从多个子集中选择一个子集,对所选择的子集中的每个响应线组进行预设操作,并根据所选择的子集的多个响应线组对应的所述预设操作的结果更新所述初始图像,所述响应线组处理单元包括:
第一坐标系转换子单元,用于根据第一坐标系(X,Y)中的像素点的像素值f1(x,y),计算位于第二坐标系(S,T)中的第一像素组的每个像素点的第一像素值g1(s,t),所述第一像素组的像素点位于所选择的响应线组对应的多条响应线上;
正投影子单元,用于根据得到的所述第一像素组的每个像素点的第一像素值g1(s,t)进行正投影计算,以获得所述第一像素组的各投影点的理论投影值p1(s,t);
反投影子单元,用于将所述第一像素组的各投影点的理论投影值p1(s,t)与实际测量投影值p0(s,t)进行对比,以获得校正系数R(s,t),并根据所述校正系数R(s,t)进行反投影计算,以获得所述第一像素组的每个像素点的第二像素值g2(s,t);
第二坐标系转换子单元,用于根据位于第二坐标系(S,T)中的第一像素组的每个像素点的第二像素值g2(s,t),计算位于第一坐标系(X,Y)中的第二像素组的每个像素点的像素值f2(x,y),所述第二像素组的每个像素点分布在所选择的响应线组对应的多条响应线上;
归一化处理单元,用于对所述第二像素组的每个像素点的像素值f2(x,y)做归一化处理;
控制单元640,用于以更新后的图像作为初始图像,继续执行步骤3),直至满足停步规则。
本领域普通技术人员可以理解上述实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成,该程序可以存储于一计算机可读存储介质中,存储介质可以包括:ROM、RAM、磁盘或光盘等。
综上所述,本发明提供的PET二维图像重建方法和装置,通过图像旋转以避免反复计算每条响应线与图像相交的像素坐标和加权系数;以及在利用飞行时间技术信息时,将多个小像素合并为一个大像素进行计算,减少了运算量,实现了OSEM算法运算结构的优化,显著降低了运算时间。
虽然本发明披露如上,但本发明并非限定于此。任何本领域技术人员,在不脱离本发明的精神和范围内,均可作各种更动与修改,因此本发明的保护范围应当以权利要求所限定的范围为准。