一种基于图像分割的动态PET图像重建及示踪动力学参数估计方法与流程

本发明属于pet成像技术领域，具体涉及一种基于图像分割的动态pet图像重建及示踪动力学参数估计方法。

背景技术：

正电子发射断层成像(positronemissiontomography，简称pet)是核医学成像的一种，在生物医学研究和临床诊疗中正逐渐显示出其重要的作用。pet成像技术通过对注入生物体内示踪剂的放射性分布进行显像，能够从分子水平上发现细胞的代谢异常，为疾病的早期诊断和预防提供了有效依据。相对于只能提供一定时间窗内生物体组织中放射性强度的平均分布的静态pet成像技术而言，动态pet成像基于动力学模型能够获得描述生物体生命活动的定量功能参数，在科学研究和临床应用中具有重要应用价值。

在动态pet的重建问题中，通常的做法是将图像分割、图像重建与动力学参数估计作为三个单独的问题分开求解。随着动态pet成像技术对时间分辨率的要求越来越高，每一帧测量数据中有限的光子计数不断地挑战着动态pet成像质量。此时，由于单帧数据的光子计数不足无法较好地反映出统计特性，利用基于统计特性的重建算法对每一帧数据进行单独重建的方法往往无法得到十分准确的重建结果。此时，通过结合示踪剂动力学模型引入来自时间纵轴的放射性强度分布的先验信息，能够在提高pet重建图像质量的基础上同时实现了对动力学参数的估计。此外，图像分割也是pet成像技术中的一个传统难题。除了对pet图像进行分割以外，还有一类做法是对动态pet的时间放射性曲线(time-activitycurve，tac)来进行聚类的操作，从而区分出生物组织里不同的功能区域。但无论是从pet的图像平面上还是从动力学参数的角度来看待这个问题，基于重建结果的分割算法准确度始终依赖与前一步重建结果的精确度，无法解决分割算法对测量噪声敏感的问题。

实际上，图像分割、图像重建与动力学参数估计这三个问题之间的关系是十分紧密的。通过将三者耦合进入同一个目标方程，通过选择合适的噪声模型能够更好地匹配pet测量数据，此时分割结果对噪声的敏感度将会降低；而分割结果的引入又能够使已有的图像估计结果更加准确，从而总体上获得一个更加符合真实情况的重建结果。

技术实现要素：

鉴于上述，本发明提供了一种基于图像分割的动态pet图像重建及示踪动力学参数估计方法，能够同时实现功能区域的分割以及动态pet联合重建。

一种基于图像分割的动态pet图像重建及示踪动力学参数估计方法，包括如下步骤：

(1)利用探测器对注入有放射性药剂的生物组织进行探测，动态采集得到对应各个时刻的符合计数向量，并建立符合计数矩阵y；

(2)使动态的pet图像序列组合成pet浓度分布矩阵x，根据pet成像原理，建立pet测量方程；

(3)通过对pet测量方程引入全变分(totalvariation，tv)约束，得到基于tv的pet图像重建模型l(x)；

(4)利用房室模型匹配估计示踪动力学参数，建立关于x和ф的同步重建模型s(x，ф)；

(5)基于预处理得到的区域分割结果，对示踪动力学参数矩阵ф进行聚类，得到聚类分割模型c(ф)；

(6)将上述三个模型l(x)、s(x，ф)和c(ф)相结合得到同步重建的目标函数lsc(x，ф)如下：

其中：γ和∈均为权重系数；

(7)对目标函数lsc(x，ф)进行最优化求解后即得到pet浓度分布矩阵x和示踪动力学参数矩阵ф。

所述符合计数矩阵y由各符合计数向量按时序排列组成，所述pet浓度分布矩阵x由各时刻对应的pet浓度分布向量(即一帧pet图像)按时序排列组成。

所述pet测量方程的表达式如下：

y＝gx+r+s

其中：g为系统矩阵，r和s分别为反映随机事件和散射事件的测量噪声矩阵。

所述pet图像重建模型l(x)的表达式如下：

其中：α为权重系数，tv(x)为关于x的全变分正则项，gij为系统矩阵g中第i行第j列元素值(其表示的是第j个像素点处出射的光子被第i个探测器接收到的概率)，yim为符合计数矩阵y中第i行第m列元素值，xjm为pet浓度分布矩阵x中第j行第m列元素值，rim为反映随机事件的测量噪声矩阵r中第i行第m列元素值，sim为反映散射事件的测量噪声矩阵s中第i行第m列元素值，i、j和m均为自然数且1≤i≤n，1≤j≤k，1≤m≤m，n为符合计数向量的维度，k为pet浓度分布矩阵x的行数即pet图像的像素点个数，m为pet浓度分布矩阵x的列数即采样时间长度。

所述全变分正则项tv(x)的表达式如下：

其中：d(xjm)为关于xjm的二维差分向量，该向量第一行元素值为xjm-xj，m+1，第二行元素值为xjm-xj+1，m，xj，m+1为pet浓度分布矩阵x中第j行第m+1列元素值，xj+1，m为pet浓度分布矩阵x中第j+1行第m列元素值，||||2表示2范数。

所述同步重建模型s(x，ф)的表达式如下：

其中：μ为权重系数，ψ为字典矩阵，^t表示转置，||||2表示2范数，||||1表示1范数。

所述聚类分割模型c(ф)的表达式如下：

其中：ch为示踪动力学参数矩阵ф中属于第h类的参数向量集合，φh为参数向量集合ch中的任一参数向量，nh为参数向量集合ch中的参数向量个数，h为自然数且1≤h≤h，h为聚类的类数，^t表示转置。

所述步骤(7)中采用admm(alternatingdirectionmethodofmultipliers，交替方向乘子算法)结合软阈值(soft-thresholding)迭代优化算法对目标函数lsc(x，ф)进行最优化求解；其中，admm针对pet浓度分布矩阵x进行迭代优化求解，软阈值针对示踪动力学参数矩阵ф迭代优化求解。

所述软阈值迭代优化算法基于以下方程进行线性处理以对示踪动力学参数矩阵ф进行迭代更新：

nh＝ih-eh/nh

其中：eh为一个全部由1组成的nh×nh维矩阵，ih为一个与矩阵eh维度相同的单位矩阵，μ为权重系数，^t表示转置，ψ为字典矩阵，ch为示踪动力学参数矩阵ф中属于第h类的参数向量集合，φh为参数向量集合ch中的任一参数向量，nh为参数向量集合ch中的参数向量个数，xh为pet浓度分布矩阵x中与参数向量φh对应的一条tac，||||2表示2范数，||||1表示1范数，h为自然数且1≤h≤h，h为聚类的类数。

所述字典矩阵ψ的表达式如下：

其中：ci(t)和ci(τ)分别为t时刻和τ时刻放射性药剂在血浆中的浓度值，和分别为关于第m组符合计数向量采集的开始时间和结束时间，θc对应为第c个房室组织指数函数的系数，m和c均为自然数且1≤m≤m，1≤c≤z，m为pet浓度分布矩阵x的列数即采样时间长度，z为大于1的自然数；θ1～θn的取值为在区间[θmin，θmax]内按指数级间隔进行选取，θmin和θmax分别为系数的上下限阈值，t和τ均表示时间。

本发明通过将统计学重建模型与有生理意义的模型相结合能够得到更加准确的联合重建结果。当重建和分割耦合到一个联合或同时的求解框架中时，对分割任务而言，可以获得噪声模型对原始投影数据进行建模的信息，而基于分割的结果也能够增强每个区域内的均匀性，从而实现一个更符合真实情况的重建结果。与其他单独重建动态pet图像或估计动力学参数的算法相比，本发明也都能获得较好的重建结果。结合本发明在模拟数据和真实数据实验中的表现，与其他单独重建动态pet图像、估计动力学参数以及图像分割的算法相比，本发明也都能获得较好的重建结果。

附图说明

图1为蒙特卡洛模拟zubal胸腔数据的模板图像。

图2(a)为zubal胸腔数据第4帧的真实图像。

图2(b)为数据计数率为5*10^4下采用ml-em方法对蒙特卡洛模拟zubal胸腔数据重建的第4帧图像结果。

图2(c)为数据计数率为1*10^5下采用ml-em方法对蒙特卡洛模拟zubal胸腔数据重建的第4帧图像结果。

图2(d)为数据计数率为5*10^4下采用本发明方法对蒙特卡洛模拟zubal胸腔数据重建的第4帧图像结果。

图2(e)为数据计数率为1*10^5下采用本发明方法对蒙特卡洛模拟zubal胸腔数据重建的第4帧图像结果。

图3(a)为zubal胸腔数据第6帧的真实图像。

图3(b)为数据计数率为5*10^4下采用ml-em方法对蒙特卡洛模拟zubal胸腔数据重建的第6帧图像结果。

图3(c)为数据计数率为1*10^5下采用ml-em方法对蒙特卡洛模拟zubal胸腔数据重建的第6帧图像结果。

图3(d)为数据计数率为5*10^4下采用本发明方法对蒙特卡洛模拟zubal胸腔数据重建的第6帧图像结果。

图3(e)为数据计数率为1*10^5下采用本发明方法对蒙特卡洛模拟zubal胸腔数据重建的第6帧图像结果。

图4(a)为蒙特卡洛模拟zubal胸腔数据在采用本发明方法下的聚类结果示意图。

图4(b)为蒙特卡洛模拟zubal胸腔数据在采用k均值分类方法下的聚类结果示意图。

图4(c)为蒙特卡洛模拟zubal胸腔数据在采用基于动力学模型谱分类方法下的聚类结果示意图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

本发明基于区域分割的动态pet及动力学参数联合重建方法，包括如下步骤：

(1)根据动态pet扫描的模式建立测量数据矩阵y和系统矩阵g。

将动态pet的扫描过程按照需要划分出一定数量的时间帧，每个时间帧内探测器采集得到的符合计数向量的可以按照时间的顺序构建出一个动态pet的测量数据矩阵y；并统计每一像素点处出射的光子被各探测器接收到的概率，从而得到系统矩阵g。

(2)根据测量数据的标记示踪化合物的正电子核素种类以及时间帧间隔的分布，基于双房室模型理论建立时间基函数的字典矩阵ψ。

若是模拟数据，则根据模拟的示踪同位核素的半衰期设置指数参数θ和探测器扫描时间间隔建立相应的字典矩阵；若是真实数据，则按照实际操作中使用的示踪同位核素和实验设置建立相应的字典矩阵。

根据房室模型可以解出组织中的放射性强度随时间的变化，即时间-放射性函数是输入组织的tac和一个δ(t)函数以及一组指数函数的卷积，即是所有组织房室的tac之和。因此，一个pet示踪房室模型可以用一组一阶线性方程来表示，每个方程表示一个不同的房室中的示踪剂放射性强度，即：

s.t.ψ0(t)＝ci(t)

其中：ct(t)表示t时刻组织中的放射性强度。n表示不同房室的总数，ψc(t)表示第c种房室对应的放射性强度随时间变化的函数，φc是ψc(t)相对应的系数。θc是对应房室的指数函数中的参数。ci(t)表示t时刻输入的放射性强度。

重建的每一帧的动态pet图像可看作在这一帧与上一帧的时间间隔内对房室组织时间-放射性曲线的积分，即：

其中，xjm是pet图像序列中第m帧的像素j所对应的放射性强度，和分别是第m帧的开始时间和结束时间，是方式组织的像素j处对应的tac。基于此关系，我们将一组tac对应的方程组扩展为一个字典矩阵如下：

其中：

(3)初始化，设置权重系数α、γ和μ，软阈值算法中的步长σ，迭代次数k置0，设置最大迭代次数kmax。

权重系数α在[2^7.5,2^8.5]左右，这个是根据原来已有的算法给出的参考范围进行调整得到的。权重系数γ一般在[1,5]左右的范围内，这个是基于实验的结果获得的。权重系数μ一般可在[0.01,0.15]的范围内进行选择得到合适的值(根据参考文献positronemissiontomographycompartmentalmodels&colon；abasispursuitstrategyforkineticmodeling中的实验结果)，当噪声增加的时候，算法不能准确地从字典中确定出正确的时间放射性曲线，此时需要适当增大μ的取值，增强稀疏约束。

(4)进入算法迭代更新过程，对于第k次迭代,固定动力学参数φ^k得到pet图像的更新结果x^k+1。

4.1提取出目标参数中仅包括pet图像x的部分获得以下子优化问题：

其中：tv(x)为全变分正则项，α为权重系数，i＝1，…，n为探测器的编号，m＝1，…，m代表时间帧，j＝1，…，k代表成像平面上像素点的编号；因此yim是第m帧第i个探测器的测量数据，xjm是第m帧第j个像素点的浓度值，而gij则是系统矩阵的一个元素，代表第j个像素点处出射的光子被第i个探测器接收到的概率，rim和sim分别表示相对应的时间帧和探测器测量数据中的随机事件和散射事件。

4.2基于全变分正则项的计算公式定义新变量ωjm并将4.1中的目标方程转化为增广拉格朗日方程如下：

其中：其中向量υjm是拉格朗日乘子，βjm则是惩罚系数，xm是第m帧pet图像像素点组成的向量，则是xm在第j个像素点处的离散偏微算子。

4.3固定x，利用二维收缩公式更新变量ωjm；

4.4固定ωjm，将la(ωjm，x)对x求导后按照如下公式更新pet图像x：

x^k+1＝x^k-ρ^kd^k(x)x^k

其中：ρ^k是通过逆向非单调线搜索确定的最快下降步长，d^k(x)是增广拉格朗日方程对x的导数；

4.5固定x和ωjm，更新拉格朗日乘子；

4.6判断是否满足迭代停止条件：x的变化＜10^-3或达到最大子迭代次数，不满足则返回4.4；若满足则迭代停止，得到pet图像x的更新结果x^k+1。

(5)对于第k次迭代,固定pet图像x^k+1和动力学参数φ^k，使用狄利克雷聚类过程更新聚类结果

(6)对于第k次迭代,固定pet图像x^k+1得到动力学参数的估计结果φ^k+1：

6.1提取出目标参数中仅包括动力学参数φ的部分获得以下子优化问题：

其中：代表着系数矩阵φ中属于h类的像素点的子集，nh是该集合中像素点的个数，ρh由第h类对应的像素点的平均值组成的矩阵，∈为控制子模型的权重系数。

6.2通过定义一个新的矩阵nh＝ih-eh/nh，假设eh是一个全部由1组成的nh×nh维的矩阵，ih是一个与eh维度相同的单位矩阵并定义新矩阵nh＝ih-eh/nh。引入化简上述方程：

6.3对上式进行线性化处理后利用软阈值迭代算子可以得到更新的动力学参数φ^k+1。

(7)判断是否满足迭代停止条件：x和φ的变化＜10^-3或达到最大迭代次数kmax，不满足则返回步骤(4)，若满足则迭代停止，得到pet图像x、示踪动力学参数φ以及聚类结果。

以下我们通过对蒙特卡洛模拟的zubal胸腔模板数据进行实验从而验证本发明系统重建和分割结果的准确性，图1为实验所用的zubal胸腔数据的模板示意图，将不同的区域分为三个感兴趣的区域(regionofinterest，roi)。实验运行环境为：8g内存，3.40ghz，64位操作系统，cpu为inteli7-3770；所模拟的pet扫描仪型号为hamamatsushr-22000，设定的放射性核素及药物为¹⁸f-fdg，设置sinogram为128个投影角度在每个角度下128条射束采集到的数据结果，系统矩阵g的大小为16384×16384。在本次试验中，对5×10⁴、1×10⁵两种不同的计数率下的投影数据进行实验。

本发明重建框架中对pet图像的结果和传统的ml-em(最大似然-期望最大法)的重建结果做比较，二者使用相同的测量数据矩阵y和系统矩阵g以控制结果的可比性。图2(a)和图3(a)均为分布真值，图2(b)～图2(e)分别是在ml-em在计数率为5×10⁴和1×10⁵的情况下以及本发明重建框架在计数率为5×10⁴和1×10⁵的情况下的动态pet图像序列中第4帧的重建结果；图3(b)～图3(e)分别是在ml-em在计数率为5×10⁴和1×10⁵的情况下以及本发明重建框架在计数率为5×10⁴和1×10⁵的情况下的动态pet图像序列中第6帧的重建结果。以明显看出联合重接获得的结果在区域内噪声明显小于ml-em获得的图像，在保证边缘对比的情况下功能区域内更加平滑，表1是其进一步的量化分析结果。

表1

图4(a)～图4(c)是聚类的比较结果，本发明重建框架得到的功能区域分布和已有的k均值分类以及基于动力学模型的谱分类方法得到的结果进行比较。可以看出联合重建的聚类结果中的噪声明显较大，这主要是由于狄利克雷过程在进行分类预处理时是以像素点为单位进行处理的，因此对噪声十分敏感。

上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。