一种PET图像的重建方法和装置
技术领域
本发明涉及医学图像处理技术领域,尤其涉及一种PET图像的重建方法和装置。
背景技术
在进行PET(Positron Emission Tomograph,正电子发射断层扫描装置)常规检查中,由于向人体注射的放射性元素(如F-18等)对人体将产生一定的伤害,注射剂量越大,扫描时间越长,将明显降低PET临床检查效率且对人体的伤害就越大。但是减少药物剂量或缩短扫描时间将导致重建图像出现大量噪声。
因此,为了能够获取符合临床要求重建图像的前提下尽量减少药物剂量和缩短扫描时间,以提高PET检查效率,需要寻找一种合适PET图像的重建方法。
现有技术中,针对药物剂量低和扫描时间短导致的PET图像噪声大的问题,常用的PET图像的重建方法是在迭代重建过程中对待重建图像X加以限制,例如先验条件和全变差(Total Variation,TV)稀疏变换约束等。迭代重建算法可将待重建图像的先验信息和稀疏变换转化为约束条件或优化准则,使得PET重建图像优化问题转化为求解具有约束条件的目标函数优化问题,从而重建出相对较好的PET图像。现有技术中PET图像重建时采用的目标函数通常包含两项:全变差项||X||TV和最小二乘项此时目标函数可以表示为α1和α2为权重因子,再采用梯度下降法(GradientDescent,GD)进行图像迭代更新。
但不完整投影数据中的噪声影响较之完整投影数据更加明显,此时采用测得的投影数据Y计算最小二乘项进行图像迭代更新会将投影数据中的固有噪声不断引入到重建图像X中,从而导致重建图像一直含有噪声,使图像背景不够均匀,影响重建质量。
发明内容
有鉴于此,本发明的第一方面提供了一种PET图像的重建方法,以提供重建图像精度,提高重建图像质量。
基于本发明的第一方面,本发明的第二方面提供了一种PET图像的重建装置。
为了解决上述技术问题,本发明采用了如下技术方案:
一种PET图像的重建方法,包括:
获取PET投影数据Y;
根据所述PET投影数据Y获取初始估计图像;
根据所述初始估计图像获取目标估计图像;
基于全变差法和最小二乘方法,根据所述目标估计图像建立目标函数;所述目标函数的表达式为:
其中,X为重建图像;α1和α2分别为全变差项||X||TV和最小二乘项的权重因子;Xk为第k步重建出的图像;Xk-1为第k-1步重建出的图像;所述全变差项||X||TV包括所述目标估计图像的每一个像素点的全变差值;
在A·X=Y的约束条件下求解所述目标函数的最小解,所述目标函数的最小解即为PET最终重建图像。
可选地,所述基于全变差法和最小二乘方法,根据所述目标估计图像建立目标函数,具体包括:
基于全变差法根据所述目标估计图像建立所述目标函数的全变差项;
基于最小二乘方法根据所述目标估计图像建立所述目标函数的最小二乘项;
对所述全变差项和所述最小二乘项分别配置第一权重因子α1和第二权重因子α2;
计算所述第一权重因子α1与所述全变差项的乘积与所述第二权重因子α2与所述最小二乘项的乘积之和,得到的结果即为所述目标函数。
可选地,所述基于全变差法根据所述目标估计图像建立所述目标函数的全变差项,具体包括:
对所述目标估计图像进行全变差稀疏变换,得到全变差稀疏变换后的目标估计图像;
计算所述全变差稀疏变换后的目标估计图像的l1范数;
根据所述l1范数分别计算所述全变差稀疏变换后的目标估计图像的每一个像素点的全变差TV(Xs,t)值,所述全变差稀疏变换后的目标估计图像的每一个像素点的全变差TV(Xs,t)值构成所述目标函数的全变差项。
可选地,所述基于最小二乘方法根据所述目标估计图像建立所述目标函数的最小二乘项,具体包括:
计算所述目标估计图像与所述初始估计图像的差值;
根据所述目标估计图像与所述初始估计图像的差值,计算范数;所述范数即为所述目标函数的最小二乘项。
可选地,所述每一个像素点的全变差值TV(Xs,t)值作为所述目标函数最小解求解过程的加权因子;所述求解所述目标函数的最小解具体为:
根据所述目标估计图像通过以下公式(II)计算下一步迭代重建的图像
其中,Xk为第k步重建出的图像;
α1和α2分别为全变差项||X||TV和最小二乘项的权重因子;
和分别为全变差项和最小二乘项在第k步迭代重建时的迭代步长;
判断迭代重建图像是否满足预设条件或迭代重建次数是否达到预设次数,如果是,确定所述迭代重建图像为PET最终重建图像和所述目标函数的最小解,如果否,将迭代重建的图像作为所述初始估计图像,返回执行所述根据所述初始估计图像获取目标估计图像的步骤。
一种PET图像的重建装置,包括:
第一获取单元,用于获取PET投影数据Y;
第二获取单元,用于根据所述PET投影数据Y获取初始估计图像;
第三获取单元,用于根据所述初始估计图像获取目标估计图像;
目标函数建立单元,用于基于全变差法和最小二乘方法,根据所述目标估计图像建立目标函数;所述目标函数的表达式为:
其中,X为重建图像;α1和α2分别为全变差项||X||TV和最小二乘项的权重因子;Xk为第k步重建出的图像;Xk-1为第k-1步重建出的图像;所述全变差项||X||TV包括所述目标估计图像的每一个像素点的全变差值;
计算单元,用于在A·X=Y的约束条件下求解所述目标函数的最小解,所述目标函数的最小解即为PET最终重建图像。
可选地,所述目标函数建立单元具体包括:
全变差项建立子单元,用于基于全变差法根据所述目标估计图像建立所述目标函数的全变差项;
最小二乘项建立子单元,用于基于最小二乘方法根据所述目标估计图像建立所述目标函数的最小二乘项;
权重因子配置子单元,用于对所述全变差项和所述最小二乘项分别配置第一权重因子α1和第二权重因子α2;
计算子单元,用于计算所述第一权重因子α1与所述全变差项的乘积与所述第二权重因子α2与所述最小二乘项的乘积之和,得到的结果即为所述目标函数。
可选地,所述全变差项建立子单元包括:
稀疏变换子单元,用于对所述目标估计图像进行全变差稀疏变换,得到全变差稀疏变换后的目标估计图像;
第一计算子单元,用于计算所述全变差稀疏变换后的目标估计图像的l1范数;
第二计算子单元,用于根据所述l1范数分别计算所述全变差稀疏变换后的目标估计图像的每一个像素点的全变差TV(Xs,t)值,所述全变差稀疏变换后的目标估计图像的每一个像素点的全变差TV(Xs,t)值构成所述目标函数的全变差项的加权因子。
可选地,所述最小二乘项建立子单元包括:
第三计算子单元,用于计算所述目标估计图像与所述初始估计图像的差值;
第四计算子单元,用于根据所述目标估计图像与所述初始估计图像的差值,计算范数;所述范数即为所述目标函数的最小二乘项。
可选地,所述每一个像素点的全变差值TV(Xs,t)值作为所述目标函数最小解求解过程的加权因子;所述计算单元包括:
第五计算子单元,用于根据所述目标估计图像通过以下公式(II)计算下一步迭代重建的图像
其中,Xk为第k步重建出的图像;
α1和α2分别为全变差项||X||TV和最小二乘项的权重因子;
和分别为全变差项和最小二乘项在第k步迭代重建时的迭代步长;
判断子单元,用于判断迭代重建图像是否满足预设条件或迭代重建次数是否达到预设次数,如果是,确定所述迭代重建图像为PET最终重建图像和所述目标函数的最小解,如果否,将迭代重建的图像作为所述初始估计图像,并将迭代重建的图像作为所述初始估计图像的信号传送至所述第三获取单元。
相较于现有技术,本发明具有以下有益效果:
本发明提供的PET图像的重建方法中,其采用的目标函数f(X)中的最小二乘项采用相邻连续两次迭代重建图像进行计算得到,替代了现有技术中采用投影数据Y计算最小二乘项的方式,因而,采用本发明的目标函数不会将投影数据中的固有噪声引入到目标图像中。
另外,由于在迭代过程中,目标估计图像逐渐趋近于目标图像。因此,本发明中的目标函数的最小二乘项的计算结果也逐渐趋近于零,因而,相较于现有技术,利用本发明实施例提供的PET图像的重建方法,能够使得图像背景变得均匀,提高了PET图像的图像质量。
附图说明
为了清楚地理解本发明的技术方案,下面将描述本发明的具体实施方式时用到的附图做一简要说明。显而易见地,这些附图仅是本发明的部分实施例,本领域普通技术人员在不付出创造性劳动的前提下还可以获得其它的附图。
图1是本发明实施例一提供的PET图像的重建方法流程示意图;
图2是本发明实施例一提供的建立目标函数的方法流程示意图;
图3是本发明实施例二提供的PET图像的重建方法流程示意图;
图4是本发明实施例三提供的PET图像的重建装置的结构示意图;
图5是本发明实施例三提供的目标函数建立单元的结构示意图;
图6是本发明实施例三提供的全变差项建立子单元的结构示意图;
图7是本发明实施例三提供的最小二乘项建立子单元的结构示意图;
图8是本发明实施例三提供的计算单元的结构示意图。
具体实施方式
为了清楚地理解本发明的技术方案,下面结合附图对本发明的具体实施方式进行描述。
为了清楚地理解本发明提供的PET图像的重建方法的发明构思,在介绍本发明的PET图像的重建方法之前,首先介绍本发明提供的PET图像的重建方法的理论基础。
压缩感知(Compressive Sensing,CS)理论可利用图像具有稀疏表示的先验知识由少量观测值重建出原始图像。由于PET重建图像经过梯度或小波变换后绝大多数的数据值等于或接近于零,可看作是稀疏的。因此,在PET图像重建中可以利用将压缩感知与普通重建方法相结合的方法实现对不完整投影数据的图像重建。
本发明实施例提供的PET图像的重建方法就是基于压缩感知理论提出的方法。
将压缩感知理论应用于PET重建是由获取到的少量的PET投影数据重建出目标图像,即在满足所测得的PET投影数据的条件下寻求最稀疏解的过程,公式如下:
公式(1)中,零范数||ψ·X||0是指非零元素个数,X指重建出的图像,ψ指稀疏变换,Y指获得的投影数据,A指投影矩阵。
在实际计算中,一般将l1范数代替l0范数,上述公式(1)更换为如下公式(2):
其中,
当采用有限差分进行稀疏变换时,采用全变差(total Variation,TV)来衡量待重建图像的梯度稀疏性,此时,||ψ·X||1即为||X||TV,此时,上述公式(2)变为公式(3):
其中,
除了采用压缩感知中常用的稀疏变换外,通过将连续相邻两次迭代重建图像进行相减,并对其求解l2范数,可以进一步将图像稀疏化,此时目标函数可以表示为:
公式(4)即为本发明实施例提供的PET图像的重建方法中采用的目标函数。
下面介绍本发明实施例提供的PET图像的重建方法的具体实施方式。
实施例一
图1是本发明实施例一提供的PET图像的重建方法的流程示意图。如图1所示,该方法包括以下步骤:
S11、获取PET投影数据Y:
需要说明的是,本发明实施例提供的PET图像的重建方法适用于获取的投影数据Y为不完整投影数据的情形。当在注射药量不足或采集时间变短的情况下进行数据采集得到的数据为不完整投影数据。
S12、根据所述PET投影数据Y获取初始估计图像:
采用本领域常用的解析重建算法如FBP算法对所述PET投影数据Y进行重建或采用全一矩阵的方法,获取初始估计图像。
S13、根据所述初始估计图像获取目标估计图像:
采用本领域常用的迭代重建算法对所述初始估计图像进行迭代重建,得到目标估计图像。
具体地,本步骤可以利用不完整投影数据Y和系统矩阵A通过迭代重建算法如EM和ART等获取目标估计图像。即求解约束条件欠定线性系统AX=Y。
需要说明的是,目标估计图像比初始估计图像更接近最终的目标图像。
S14、基于全变差法和最小二乘方法,根据所述目标估计图像建立目标函数:
如上述所述的公式(4)可知,本发明实施例采用的目标函数包括全变差项和最小二乘项。所以,在建立目标函数时,包括建立目标函数的全变差项和建立目标函数的最小二乘项两个步骤。
具体地,如图2所示,上述基于全变差法和最小二乘方法,根据所述目标估计图像建立目标函数,具体包括:
S141、基于全变差法根据所述目标估计图像建立所述目标函数的全变差项:
具体地,该步骤包括以下分步骤:
S141a、对所述目标估计图像进行全变差稀疏变换,得到全变差稀疏变换后的目标估计图像。
S141b、计算所述全变差稀疏变换后的目标估计图像的l1范数。
S141c、根据所述l1范数分别计算所述全变差稀疏变换后的目标估计图像的每一个像素点的全变差TV(Xs,t)值,所述全变差稀疏变换后的目标估计图像的每一个像素点的全变差TV(Xs,t)值构成所述目标函数的全变差项。
换句话说,所述全变差项||X||TV包括所述目标估计图像的每一个像素点的全变差值。具体地说,所述全变差项||X||TV可以为所述目标估计图像的所有像素点的全变差值的阵列或集合。进一步地,为了提高重建图像的质量,所述全变差项||X||TV可以为带加权的全变差项。全变差项的加权可以在求解目标函数最小化计算中通过对全变差项求偏导后与对应像素点的全变差TV(Xs,t)值相乘来实现。
其中,每个像素点的全变差TV(Xs,t)的计算公式如下:
S142、基于最小二乘方法根据所述目标估计图像建立所述目标函数的最小二乘项:
需要说明的是,本发明采用的目标函数的最小二乘项是根据相邻两次迭代重建图像相减得到的。
本步骤具体包括以下步骤:
S142a、计算所述目标估计图像与所述初始估计图像的差值:
需要说明的是,由于目标估计图像是根据初始估计图像进行迭代重建一次后得到的图像,所以,目标估计图像与初始估计图像为相邻两次迭代重建图像。
S142b、根据所述目标估计图像与所述初始估计图像的差值,计算范数;所述范数即为所述目标函数的最小二乘项:
需要说明的是,在步骤中,Xk为目标估计图像,Xk-1为初始估计图像。
S143、对所述全变差项和所述最小二乘项分别配置第一权重因子α1和第二权重因子α2:
配置的第一权重因子α1决定了全变差项在目标函数中所占的影响大小,同样,配置的第二权重因子α2决定了最小二乘项在目标函数中所占的影响大小。因此,在配置权重因子α1和α2时,需要根据具体的投影数据和目标重建图像进行选择配置。
S144、计算所述第一权重因子α1与所述全变差项的乘积与所述第二权重因子α2与所述最小二乘项的乘积之和,得到的结果即为所述目标函数:
所述目标函数的表达式为上述公式(4)所示,即
其中,X为重建图像;α1和α2分别为全变差项||X||TV和最小二乘项的权重因子;Xk为第k步重建出的图像;Xk-1为第k-1步重建出的图像。
在本发明采用的目标函数中,其最小二乘项为相邻两次迭代重建图像的差值的范数。因而,本发明提供的PET图像的重建方法不会将投影数据Y中的固有噪声引入到目标图像中。
而且,由于迭代重建过程中,第k步重建出的图像比第k-1步重建出的图像更接近目标图像。所以,随着迭代重建的进行,本发明的目标函数的最小二乘项的计算结果逐渐接近零。因此,通过本发明提供的PET图像的重建方法得到的PET图像背景变得均匀,进而,提高了PET图像的图像质量。
S15、在A·X=Y的约束条件下求解所述目标函数的最小解,所述求解所述目标函数的最小解的过程即为初始估计图像的迭代重建过程:
用公式表示步骤S15即为:
由于在上述求解目标函数的最小解的过程中,将待重建图像的先验信息和稀疏变换转化为约束条件或优化准则,所以,上述求解目标函数的最小解的过程即为PET图像的迭代重建过程。所述目标函数的最小解即为PET最终重建图像。
在本发明实施例中,优选梯度下降法对目标估计图像进行迭代重建。其迭代公式为:
其中,hk表示迭代步长,pk=-▽f(Xk)表示搜索方向。因为梯度方向是函数值下降最快的方向,因此,优选沿梯度方向查找。
在进行迭代重建时,选择的迭代步长hk,利用公式(7)对待重建图像逐步进行迭代重建。将利用公式(7)得到的图像作为初始估计图像,返回步骤S13,然后顺序执行步骤S13至步骤S15,直到得到的迭代重建图像满足预设条件或迭代次数满足预设次数。最后一次迭代得到的重建图像即为待重建的目标图像,该目标图像即视为所述目标函数的最小解。
以上为本发明实施例一提供的PET图像的重建方法的具体实施方式。
需要说明的是,在上述目标函数(4)的全变差项||X||TV中,对全变差项求偏导数后,计算结果的数量级均在1左右,但在通常情况下PET图像中不同位置的像素值的数量级与1相差一个或若干个数量级,用固定数量级结果对不同位置不同数量级的像素值进行更新将影响重建图像精度,最终也将直接影响重建图像质量。
为了使得在对目标函数的求解过程中,图像重建计算公式中对应的全变差项的数量级与不同位置像素值的数量级相适应,进而提高图像清晰度,本发明实施例还提供了PET图像的重建方法的另外一种实施方式,具体参见实施例二。
实施例二
图3是本发明实施例二提供的PET图像的重建方法的流程示意图。如图3所示,该方法包括以下步骤:
S31至步骤S34与实施例一中的步骤S11至步骤S14相同,为了简要起见,在此不再详细描述,具体参见实施例一的描述。
该方法还包括:
S35、根据所述目标估计图像以及所述全变差稀疏变换后的目标估计图像的每一个像素点的全变差TV(Xs,t)值,通过以下公式(8)计算下一步迭代重建的图像
所述每一个像素点的全变差值TV(Xs,t)值作为所述目标函数最小解求解过程的加权因子。具体地,该每一个像素点的全变差值TV(Xs,t)值与其对应的全变差项的偏导数的乘积作为图像重建计算公式中的全变差项。具体如公式(8)所示。
其中,Xk为第k步重建出的图像;
α1和α2分别为全变差项||X||TV和最小二乘项的权重因子;
和分别为全变差项和最小二乘项在第k步迭代重建时的迭代步长。
在图像重建计算公式(8)中,采用对每一个像素点的全变差值TV(Xs,t)值与其对应的全变差项的偏导数相乘的方法进行加权,使得相乘后的乘积数量级与重建图像X不同位置的像素值相适应。而且,利用该加权方法能够克服图像边缘模糊的问题,能够提高图像的清晰度。这是因为:位于重建图像X的边缘的像素对应的全变差值比位于重建图像背景均匀区域内的像素对应的全变差值要大,所以,利用不同像素点的全变差值TV值可以使得图像边缘位置全变差项偏导数值进一步变大,且使背景均匀区域内的全变差项偏导数值缩小,将使两个区域全变差项偏导数差值更大,然后再由公式(8)进行迭代更新时可使迭代结果在图像边缘位置与背景区域对比度较之原先方法变好,图像相对更清晰。
S36、判断所述迭代重建图像是否满足预设条件,或者判断所述迭代重建次数是否达到预设次数,如果是,执行步骤S37,如果否,执行步骤S38。
S37、确定所述迭代重建图像为PET图像和目标函数的最小解。
S38、将通过上述公式(8)得到的迭代重建的图像作为所述初始估计图像,返回步骤S33。
基于上述实施例一和实施例二所述的PET图像的重建方法,本发明还提供了PET图像的重建装置。具体参见实施例三。
实施例三
图4是本发明实施例三提供的PET图像的重建装置的结构示意图。如图4所示,该重建装置包括以下单元:
第一获取单元41,用于获取PET投影数据Y;
第二获取单元42,用于根据所述PET投影数据Y获取初始估计图像;
第三获取单元43,用于根据所述初始估计图像获取目标估计图像;
目标函数建立单元44,用于基于全变差法和最小二乘方法,根据所述目标估计图像建立目标函数;所述目标函数的表达式为:
其中,X为重建图像;α1和α2分别为全变差项||X||TV和最小二乘项的权重因子;Xk为第k步重建出的图像;Xk-1为第k-1步重建出的图像;所述全变差项||X||TV包括所述目标估计图像的每一个像素点的全变差值;
计算单元45,用于在A·X=Y的约束条件下求解所述目标函数的最小解,所述目标函数的最小解即为重建的PET图像。
通过以上所述的重建装置,能够避免将投影数据中的固有噪声引入到重建图像中,使得重建出的图像背景变得均匀,提高了PET图像的图像质量。
作为本发明的一个具体实施例,如图5所示,所述目标函数建立单元44具体包括:
全变差项建立子单元441,用于基于全变差法根据所述目标估计图像建立所述目标函数的全变差项;
最小二乘项建立子单元442,用于基于最小二乘方法根据所述目标估计图像建立所述目标函数的最小二乘项;
权重因子配置子单元443,用于对所述全变差项和所述最小二乘项分别配置第一权重因子α1和第二权重因子α2;
计算子单元444,用于计算所述第一权重因子α1与所述全变差项的乘积与所述第二权重因子α2与所述最小二乘项的乘积之和,得到的结果即为所述目标函数。
进一步地,如图6所示,所述全变差项建立子单元441包括:
稀疏变换子单元441a,用于对所述目标估计图像进行全变差稀疏变换,得到全变差稀疏变换后的目标估计图像;
第一计算子单元441b,用于计算所述全变差稀疏变换后的目标估计图像的l1范数;
第二计算子单元441c,用于根据所述l1范数分别计算所述全变差稀疏变换后的目标估计图像的每一个像素点的全变差TV(Xs,t)值,所述全变差稀疏变换后的目标估计图像的每一个像素点的全变差TV(Xs,t)值构成所述目标函数的全变差项。
进一步地,如图7所示,所述最小二乘项建立子单元442可以包括:
第三计算子单元442a,用于计算所述目标估计图像与所述初始估计图像的差值;
第四计算子单元442b,用于根据所述目标估计图像与所述初始估计图像的差值,计算范数;所述范数即为所述目标函数的最小二乘项。
进一步地,为了提高图像的清晰度,所述每一个像素点的全变差值TV(Xs,t)值作为所述目标函数最小解求解过程的加权因子;如图8所示,所述计算单元45包括:
第五计算子单元451,用于根据所述目标估计图像通过以下公式(8)计算下一步迭代重建的图像
其中,Xk为第k步重建出的图像;
α1和α2分别为全变差项||X||TV和最小二乘项的权重因子;
和分别为全变差项和最小二乘项在第k步迭代重建时的迭代步长;
判断子单元452,用于判断迭代重建图像是否满足预设条件或迭代重建次数是否达到预设次数,如果是,确定所述迭代重建图像为PET最终重建图像和所述目标函数的最小解,如果否,将迭代重建的图像作为所述初始估计图像,并将迭代重建的图像作为所述初始估计图像的信号传送至所述第三获取单元。
对于装置实施例而言,由于其基本对应于方法实施例,所以相关之处参见方法实施例的部分说明即可。以上所描述的装置实施例仅仅是示意性的,其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性劳动的情况下,即可以理解并实施。
以上所述仅是本发明的优选实施方式,应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。