CN105212957B - 一种基于TV Merge的晶体级PET系统时间修正方法 - Google Patents
一种基于TV Merge的晶体级PET系统时间修正方法 Download PDFInfo
- Publication number
- CN105212957B CN105212957B CN201510526548.4A CN201510526548A CN105212957B CN 105212957 B CN105212957 B CN 105212957B CN 201510526548 A CN201510526548 A CN 201510526548A CN 105212957 B CN105212957 B CN 105212957B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- msup
- time
- msubsup
- 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.)
- Active
Links
Landscapes
- Nuclear Medicine (AREA)
- Measurement Of Radiation (AREA)
Abstract
本发明公开了一种基于TV Merge的晶体级PET系统时间修正方法，包括：(1)扫描生物组织获取single model数据；(2)分别对探测器单元的行方向和列方向进行晶体级分割；(3)对数据进行预处理；(4)将TV约束加到线性方程式中；(5)使用ADM法求得时间修正序列；(6)使用Merge方法将两组时间修正序列合并成一个最终的时间修正序列。本发明方法通过将TV约束引入到时间修正线性过程中，并通过merge法将时间修正水平从探测器级提升到晶体级，来获取更精确的时间修正序列。这种方法有效地加强了修正算法对噪声的过滤，提高了PET时间修正的准确率，提升了PET系统的空间分辨率。
Description
技术领域
本发明属于PET成像技术领域，具体涉及一种基于TVMerge(Total VariationMerge，全变差融合法)的晶体级PET系统时间修正方法。
背景技术
PET全称为Positron emission tomography，也就是通常所说的正电子发射断层成像，是一种基于核物理学和分子生物学的医学影像技术，它能够从分子层面上观察细胞的新陈代谢活动，为早期疾病尤其是肿瘤的检测和预防提供了有效依据。PET本质上是对病人体内药物的浓度分布进行成像，被注射入病人体内的放射性同位核素标记药物通过血液进入循环系统，这些物质在人体内各组织器官中将形成一定的浓度分布。由于放射性同位核素的半衰期较短，且极其不稳定，将很快发生衰变，衰变过程中所释放的正电子与附近的自由电子发生湮灭反应，产生一对方向几乎相反、能量相等，能量大小为511kev的伽玛光子对。这些光子对被PET系统中的探测器环接收，生成记录有光子能量，探测时间，计数率和探测器编号的有效数据(Single model或者sinogram)。之后，这些数据被用于生理图像的重建或生理参数的估计。
近几年PET在实际医学领域的应用日趋广泛，但与此同时，很多医学领域都需要PET能提供更高的空间分辨率，以实现更加精确的医疗诊断。为了获得更高的空间分辨率，一种被称为TOF(Time-of-Flight)-PET的新系统正被广泛的应用在相关的临床医学领域。TOF-PET的基本原理是通过记录探测器探测到光子的精确时间来提升空间分辨率。因此，TOF-PET对PET系统的时间分辨率有着很高的要求，但是在实际情况中，PET系统的时间分辨率常会收到探测晶体的延时、探测器部分的延时和后端电路的延时的影响，使PET系统的时间分辨率变差。所以，PET系统的时间修正对实现高分辨率PET成像是十分必要的。此外，随着闪烁晶体的尺寸不断缩小和PET系统探测器数量的不断增加，准确的大型PET系统时间修正也正在变的越来越困难。
目前，PET时间修正方法大致可分为三种：参考探测器法、特殊散射源法和线性转化法。第一种方法主要是用一个快速光电倍增管作为参考探测器，通过同时记录同一事件在参考探测器和PET系统探测器所记录的探测时间，来估计两者之间的时间差，从而获取PET系统的准确延时并借此对系统进行时间修正；第二种方法则是使用一个经过特殊设计的放射性源来获取系统时间修正序列。在这个特殊设计的放射性源中，每个事件在源内的具体位置都是已知的，由于所有伽马光子的传播速度都是光速，因此我们可以通过计算其飞行时间来获取理论上无延时的探测时间，这个理论值和实际PET系统测量的时间的偏差就是我们所求的系统的时间修正序列；第三种方法则是将PET时间修正问题转化为一个线性过程，通过最小二乘法来求解这一线性过程来获取时间修正序列的估计值。
但是以上这三种方法都有着各自的局限性。第一种方法需要一个额外的参考探测器来获取时间修正序列，而这个额外的探测器会增加整个系统的构建成本，并且为了获取较为准确的估计值，这种方法往往需要较长的采集时间，这都制约了其的普遍应用。而第二种方法虽然不需要额外添加一个参考探测器，但是它却需要一个特别设计的放射性源来计算时间修正序列，这也限制了它的广泛使用。第三种方法虽然对探测器和放射性源都没有特殊的要求，但是由于所采集到的信号存在噪声，所以这种方法的准确度不是很好，并且还受到系统矩阵大小的制约，无法应用于大型的PET系统中。
此外，随着临床医疗对医学图像精确程度要求的不断提高，与之相关的闪烁晶体的尺寸也随之不断变小，晶体数量不断增加，PET的各种修正也不断的向着晶体级和大系统级迈进。但是，以上的方法大多都是探测器级的修正，且无法很好的应用在晶体数巨大的PET系统中，所以急需一种新的晶体级的时间修正方法。
发明内容
针对现有技术所存在的上述技术问题，本发明提供了一种基于TV Merge的晶体级PET系统时间修正方法，能够有效提高PET系统的时间分辨率。
一种基于TV Merge的晶体级PET系统时间修正方法，包括如下步骤：
(1)对PET系统中的每个探测器按行方向进行晶体级分割，每个探测器通过分割对应得到n个晶体单元，n为除去1以外m的任一约数，m为探测器原晶体阵列的维度；
(2)利用晶体级分割后的探测器对注入放射性示踪剂的生物组织进行扫描探测，得到多组LOR数据；
(3)对每一组LOR数据进行预处理，以剔除每组LOR数据中时间信息有极大偏差的Single model数据记录，并确定每组LOR数据的探测时延；
(4)将PET系统时间修正过程转化为线性方程，通过TV(Total Variation，全变差)对该线性方程进行约束得到以下目标函数L；进而根据由各组LOR数据探测时延组成的探测时延序列ΔT，对目标函数L进行最小化求解得到PET系统的时间修正序列x；
其中：A为系统矩阵，|| ||为L2范数，λ和β均为预设的权重系数，θi和ε均为权重系数向量，Di()表示括号内的向量中第i个元素对应的梯度向量，μi为离散梯度向量μ中的第i个元素，i为自然数且1≤i≤N，T表示转置，N为时间修正序列x的维度且为PET系统内所有探测器晶体单元的总个数；
(5)根据步骤(1)～(4)对PET系统中的每个探测器按列方向进行晶体级分割，并计算得到对于列方向分割的时间修正序列y，进而将时间修正序列x和y进行融合得到一个nN维的时间修正序列z；
(6)对PET系统中的每个探测器进行晶体级分割，每个探测器通过分割对应得到一个由n×n个晶体单元组成的晶体单元阵列，利用晶体级分割后的探测器对注入放射性示踪剂的生物组织进行扫描探测，得到多组LOR数据；最后根据所述的时间修正序列z对各组LOR数据Singlemodel数据记录中晶体单元的探测时间进行修正。
进一步地，所述的步骤(2)中的每组LOR数据对应一对探测到同一偶合事件且分属于不同探测器内的晶体单元JA和JB，每组LOR数据包含多条Single model数据记录，其中每条Single model数据记录对应一个探测时间差即晶体单元JA和JB对于同一偶合事件的探测时间之差，每条Single model数据记录包含晶体单元JA和JB的编号、对应探测时间差精度范围内晶体单元JA和JB对于同一偶合事件的探测计数以及对应每次探测计数晶体单元JA和JB对于同一偶合事件各自的探测时间。
进一步地，所述的步骤(3)中对每组LOR数据进行预处理的具体方法为：对于任一组LOR数据，将该组LOR数据中所有Single model数据记录按探测计数做成直方图，并计算所有Single model数据记录探测计数的均值，进而剔除探测计数小于均值的Single model数据记录。
进一步地，所述的步骤(3)中确定每组LOR数据探测时延的具体方法为：对于预处理后的任一组LOR数据，将其中所有Single model数据记录对应探测时间差的最大值作为该组LOR数据的探测时延。
进一步地，所述的步骤(4)中线性方程的表达式为A·x＝ΔT。
进一步地，所述的步骤(4)中通过以下迭代方程对目标函数L进行最小化求解，具体算式如下：
εk+1＝εk-λ(Axk-ΔT)
其中：xk和xk+1分别为第k次迭代和第k+1次迭代的时间修正序列，μk和μk+1分别为第k次迭代和第k+1次迭代的离散梯度向量，和为权重系数向量θi对应第k次迭代和第k+1次迭代的结果，εk和εk+1为权重系数向量ε对应第k次迭代和第k+1次迭代的结果，μi k为离散梯度向量μk中的第i个元素，k为迭代次数。
所述的时间修正序列xk+1通过以下迭代方程求解：
xk+1＝xk-αkg(xk)
其中：αk为第k次迭代的求解步长，其通过Amijo线性搜索获取，g(xk)为关于时间修正序列xk的梯度函数。
所述的离散梯度向量μk+1通过以下迭代方程求解：
所述系统矩阵A的维度为M×N，M＝(N-1)N/2，系统矩阵A的具体表达如下：
其中：A1～AN-1均为系统矩阵A的子矩阵，对于任一子矩阵Aj，其维度为(N-j)×N，j为自然数且1≤j≤N-1；所述的子矩阵Aj中第j列向量的所有元素均为1，前j-1列向量的所有元素均为0，后N-j列向量所组成的方阵为主对角线元素均为-1的对角矩阵。
进一步地，所述的步骤(5)中时间修正序列x和y进行融合的具体方法如下：
5.1对于PET系统中的任一探测器，从时间修正序列x提取该探测器各晶体单元对应的n个时间修正值x1～xn，从时间修正序列y提取该探测器各晶体单元对应的n个时间修正值y1～yn；
5.2计算时间修正值y1～yn的平均值对于时间修正值x1～xn中的任一时间修正值xp，根据以下算式将其扩展成n个时间修正值xp1～xpn；依此遍历时间修正序列x中每个探测器晶体单元对应的时间修正值，从而使时间修正序列x扩展成一个nN维的时间修正序列x*；
其中：p和q均为自然数且1≤p≤n，1≤q≤n；
5.3计算时间修正值x1～xn的平均值对于时间修正值y1～yn中的任一时间修正值yq，根据以下算式将其扩展成n个时间修正值yq1～yqn；依此遍历时间修正序列y中每个探测器晶体单元对应的时间修正值，从而使时间修正序列y扩展成一个nN维的时间修正序列y*；
5.4将时间修正序列x*和y*相加后除以2，即得到时间修正序列z。
本发明PET系统时间修正方法通过利用TV约束与晶体级分割和Merge方法重构，有效地提升了传统线性转换时间修正法对噪声的处理能力，并解决了原有方法不能应用于大型PET系统和系统大小受到计算机内存大小制约的问题，成功将PET系统时间修正提升到了晶体级别，提高了PET系统时间修正序列估计值的准确率，使PET系统能获得更加优秀的时间分辨率，从而能更好的利用TOF信息来获得更好的空间分辨率，使PET系统能为临床诊断和药物研发提出更加精准有用的信息。
附图说明
图1为本发明PET系统时间修正方法的步骤流程示意图。
图2为本发明探测器晶体级分割的示意图。
图3为本发明Merge法的流程示意图。
图4为数据采集和验证的流程示意图。
图5(a)为本发明方法与传统最小二乘法在探测器单元级上时间修正效果对比结果示意图。
图5(b)为本发明基于TV merge法在不同晶体分割级别的时间修正结果和时间分辨率随着分割数的变化趋势示意图。
具体实施方式
为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明PET系统时间修正方法进行详细说明。
如图1所示，本发明基于TV merge法的晶体级PET时间修正方法，包括如下步骤：
(1)利用探测器对注入放射性示踪剂(18F-FDG)的圆柱形体模(Cylinderphantom)进行扫描探测，记录探测器探测到的每一件偶合事件的探测时间和相应探测器编号，生成含有时间，探测器对和计数率信息的Single model数据。
(2)根据计算内存空间和具体的探测器数及每个探测器模块所含的晶体数，对每个探测器进行次晶体级分割；其中，对于每一个探测器模块，根据与其耦合的晶体阵列的数量对其进行次晶体级分割，例如：当晶体阵列为32×32时，可选取的分割数e为32的约数即2,4,8,16,32，其中对于晶体阵列的行与列都进行相同分割数e的分割，其具体过程如图2所示。具体的分割数取决于探测器总数、数据采集总时间和计算使用的电脑内存的大小。之后，按照新的行和列分割获取对应的两组行方向和列方向的时间，分割后探测单元对和计数率的数据。
(3)对数据进行预处理，将每一条LOR上记录的事件数按探测时间做成直方图，求出直方图的均值，将计数率小于均值的数据剔除，把直方图峰值处所对应的时间记作这条LOR所对应的时间。这样变换可以剔除时间信息有极大偏差的噪声数据，获取每一条LOR上真实的时间，探测器对和计数率信息。
(4)将PET系统时间修正过程转化为线性过程，其具体表达式为：
A·x＝ΔT
其中：x为所求的时间修正序列，其维数为系统中最小探测单元的总数；ΔT为探测单元的延时序列，其具体表达式为ΔT＝Tl-Tk，这里Tl和Tk分别代表对于某一特定事件探测单元l和k分别记录的时间信息，其维数为系统中记录的LOR的总数；A为系统矩阵，用于表示系统中每一条LOR对应的探测器单元对信息，其具体形式如下：
其中：un(·)为算法中最小的探测单元，n为系统中最小探测单元的总数，m为PET系统中记录的LOR的总数，理论上m的大小受到n大小的约束，可以表达为m＝(n-1)n/2。
为了获取精确的时间修正序列估计值，将TV约束加在线性方程中，则其表达式变为：
其中：|| ||为L2范数，TV()为TV约束，这里TV约束可以定义为即为序列x中每一个元素的离散梯度之和；λ为权重系数，用于平衡式中第一项和第二项的估计精度。Dix为二维向量且该二维向量中的两个元素值分别为xi,+1x-xi和xi,+1y-xi；xi为x中第i个元素，xi,+1x为xi所在行里xi后一行元素，xi,+1y为xi所在行里xi后一列元素。
(5)使用迭代法求解TV约束后的线性方程，其具体过程描述如下：由于无法直接求解添加TV约束后线性方程，所以我们引入新参数μ，其可以定义为μi＝Dix，然后原方程式可以转换为：
为了求解这一问题，我们将其转换为求解一个拉格朗日方程式最小值的问题，用L(x,μi)表示这一问题，其具体表达式为：
其中：上式中第一项为约束项，其他四项为惩罚项，用于保证最终估计值的准确性和鲁棒性。这里，第二项和第四项为线性惩罚项，以保证估计值满足线性方程的线性特征；第三项和第五项为二次惩罚项，用于使估计值具有一定的鲁棒性。θi、β、ε和λ分别对应四个惩罚项的权重系数，其中，β和λ为选定的常数，θi和ε的值则随着每次迭代改变。为了能有效快速地求解这一问题，我们使用可变方向法(Alternating direction method,ADM)来求解这一问题。在ADM中，首先，我们假定μi为已知来求解x的值，其具体计算过程如下：
其中，k和k+1分别代表第k和k+1次迭代，这个方程的梯度g(xk)可以表示为：
由于我们要求的是方程的最小值，所以使g(xk)＝0，这时，我们就可以得到xk的解，其表达式如下：
xk+1＝xk-αkg(xk)
其中，αk为求解步长，我们通过Amijo非单调线性搜索来获取。
之后，我们假定x为已知来求解μi的值，其具体计算过程如下：
这一方程可以通过二维阈收缩值法求解，其具体表达式为：
最后，我们更新方程中的惩罚项权重θi、β、ε和λ。两个常数权重β和λ的理论取值范围为24到213，而的θi和ε值则是在每次迭代中更新，其具体的表达式如下：
εk+1＝εk-λ(Axk-ΔT)
之后，经过多次迭代之后，获取PET时间修正序列。迭代次数达到300次或满足迭代停止判断是迭代结束，迭代停止判断的表达式如下：
(6)如图3所示，本实施方式通过Merge方法将两组时间修正序列合并成一个最终的时间修正序列具体描述如下：在分别获得行方向和列方向的时间修正序列后，假设两个序列分别表示为和其中，j为每个探测器单元内所对应的晶体阵列的晶体分割数，Di(i＝1～n)为探测器编号，n为探测器总数；让后按探测器编号对行方向和列方向上的时间修正序列进行merge操作，合成最终的时间修正序列。其具体过程如下：以行方向上同一探测器编号下时间修正序列为基底，以列方向上同一探测器编号下时间修正序列为分布，对所有探测器编号，将基底按分布从新投影成一组新的的序列具体投影方式如下：
其中：为分布序列的均值；之后再交换基底和分布，用一样的方法获取一组新的序列求两个序列的平均值，所得的新时间修正序列即为最终的时间修正序列。
以下我们通过实验来验证本实施方式的实用性和可靠性，图4所示了实验的基本流程。
实验中，我们使用的PET系统是日本滨松光电子有限公司的大脑PET系统(HITS-655000)，其包含5个探测器环，每个探测器环分为4的探测深度，每个探测深度由32个探测器单元组成，而每个探测器单元又和一个的32×32的晶体阵列耦合。这一系统总共包括655,360枚晶体，理论上最大LOR数为140,928,614,400。
我们所使用的Cylinder Phantom的尺寸为直径200毫米，长度220毫米。所注射的18F-FDG的放射性浓度为282.3Mbq，整体扫描时间为10小时。这里使用这样长时间的扫描是为了保证每个LOR上的数据量足够充分
为了检验时间修正序列的准确性，我们需要通过对比修正前后的系统时间分辨率。这里我们将放射性点源置于PET系统中心，记录其Single model数据信息。之后使用之前获得的时间修正序列对数据进行修正，列出修正后的数据的时间直方图，来获取修正后的时间分辨率。本实验中所使用的放射性点源为22-Na，其直径为0.25毫米，被置于PET系统中心，总扫描时间为5分钟。
这里为了能充分理解真正晶体级修正的好处，我们对晶体阵列做了多种分割用于比较。由于我们所使用的晶体阵列的尺寸为32×32，所以我们选取的分割数为32的全部约数即2,4,8,16,32。然后我们分别对每一种分割数都进行了计算，获得了其修正后的系统时间分辨率。
为了能够充分的验证本实施方式，我们设计了两组实验来验证我们的方法，具体情况如图5所示；其中，图5(a)为本发明方法与传统最小二乘法在探测器单元级上时间修正效果对比结果，图5(b)为TV merge在不同晶体分割级别的时间分辨率随着分割数变化的趋势。其具体的时间分辨率如表1所示：
表1
通过以上的实验结果我们可以看出，本发明基于TV merge法的晶体级PET时间修正方法的修正结果有效地提高了PET时间修正的准确率，提升了PET系统的时间分辨率，并使其可以真正的做到晶体级的修正，并可以应用于大系统中，不受计算机内存的限制。
Claims (10)
1.一种基于TV Merge的晶体级PET系统时间修正方法，包括如下步骤，TV Merge表示全变差融合法；
(1)对PET系统中的每个探测器按行方向进行晶体级分割，每个探测器通过分割对应得到n个晶体单元，n为除去1以外m的任一约数，m为探测器原晶体阵列的维度；
(2)利用晶体级分割后的探测器对注入放射性示踪剂的生物组织进行扫描探测，得到多组LOR数据；
(3)对每一组LOR数据进行预处理，以剔除每组LOR数据中时间信息有极大偏差的Single model数据记录，并确定每组LOR数据的探测时延；
(4)将PET系统时间修正过程转化为线性方程，通过全变差对该线性方程进行约束得到以下目标函数L；进而根据由各组LOR数据探测时延组成的探测时延序列ΔT，对目标函数L进行最小化求解得到PET系统的时间修正序列x；
<mrow>
<mi>L</mi>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mo>{</mo>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
<mo>|</mo>
<mo>|</mo>
<mo>-</mo>
<msubsup>
<mi>&theta;</mi>
<mi>i</mi>
<mi>T</mi>
</msubsup>
<mrow>
<mo>(</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mi>&beta;</mi>
<mn>2</mn>
</mfrac>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
<mo>|</mo>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
<mo>}</mo>
<mo>-</mo>
<msup>
<mi>&epsiv;</mi>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>A</mi>
<mi>x</mi>
<mo>-</mo>
<mi>&Delta;</mi>
<mi>T</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mi>&lambda;</mi>
<mn>2</mn>
</mfrac>
<mo>|</mo>
<mo>|</mo>
<mi>A</mi>
<mi>x</mi>
<mo>-</mo>
<mi>&Delta;</mi>
<mi>T</mi>
<mo>|</mo>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
</mrow>
其中：A为系统矩阵，|| ||为L2范数，λ和β均为预设的权重系数，θi和ε均为权重系数向量，Di()表示括号内的向量中第i个元素对应的梯度向量，μi为离散梯度向量μ中的第i个元素，i为自然数且1≤i≤N，T表示转置，N为时间修正序列x的维度且为PET系统内所有探测器晶体单元的总个数；
(5)根据步骤(1)～(4)对PET系统中的每个探测器按列方向进行晶体级分割，并计算得到对于列方向分割的时间修正序列y，进而将时间修正序列x和y进行融合得到一个nN维的时间修正序列z；
(6)对PET系统中的每个探测器进行晶体级分割，每个探测器通过分割对应得到一个由n×n个晶体单元组成的晶体单元阵列，利用晶体级分割后的探测器对注入放射性示踪剂的生物组织进行扫描探测，得到多组LOR数据；最后根据所述的时间修正序列z对步骤(6)中获得的各组LOR数据Single model数据记录中晶体单元的探测时间进行修正。
2.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(2)中的每组LOR数据对应一对探测到同一偶合事件且分属于不同探测器内的晶体单元JA和JB，每组LOR数据包含多条Single model数据记录，其中每条Single model数据记录对应一个探测时间差即晶体单元JA和JB对于同一偶合事件的探测时间之差，每条Single model数据记录包含晶体单元JA和JB的编号、对应探测时间差精度范围内晶体单元JA和JB对于同一偶合事件的探测计数以及对应每次探测计数晶体单元JA和JB对于同一偶合事件各自的探测时间。
3.根据权利要求2所述的PET系统时间修正方法，其特征在于：所述的步骤(3)中对每组LOR数据进行预处理的具体方法为：对于任一组LOR数据，将该组LOR数据中所有Singlemodel数据记录按探测计数做成直方图，并计算所有Single model数据记录探测计数的均值，进而剔除探测计数小于均值的Single model数据记录。
4.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(3)中确定每组LOR数据探测时延的具体方法为：对于预处理后的任一组LOR数据，将其中所有Singlemodel数据记录对应探测时间差的最大值作为该组LOR数据的探测时延。
5.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(4)中线性方程的表达式为A·x＝ΔT。
6.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(4)中通过以下迭代方程对目标函数L进行最小化求解，具体算式如下：
<mrow>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>=</mo>
<munder>
<mi>min</mi>
<mi>x</mi>
</munder>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mo>{</mo>
<mo>-</mo>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&theta;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>-</mo>
<msubsup>
<mi>&mu;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mi>&beta;</mi>
<mn>2</mn>
</mfrac>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msubsup>
<mi>&mu;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mo>|</mo>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
<mo>}</mo>
<mo>-</mo>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mi>&epsiv;</mi>
<mi>k</mi>
</msup>
<mo>)</mo>
</mrow>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>A</mi>
<mi>x</mi>
<mo>-</mo>
<mi>&Delta;</mi>
<mi>T</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mi>&lambda;</mi>
<mn>2</mn>
</mfrac>
<mo>|</mo>
<mo>|</mo>
<mi>A</mi>
<mi>x</mi>
<mo>-</mo>
<mi>&Delta;</mi>
<mi>T</mi>
<mo>|</mo>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>=</mo>
<munder>
<mi>min</mi>
<mi>&mu;</mi>
</munder>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mo>{</mo>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
<mo>|</mo>
<mo>|</mo>
<mo>-</mo>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&theta;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>)</mo>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mi>&beta;</mi>
<mn>2</mn>
</mfrac>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
<mo>|</mo>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
<mo>}</mo>
</mrow>
<mrow>
<msubsup>
<mi>&theta;</mi>
<mi>i</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>&theta;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mo>-</mo>
<mi>&beta;</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<msup>
<mi>x</mi>
<mi>k</mi>
</msup>
<mo>)</mo>
<mo>-</mo>
<msubsup>
<mi>&mu;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mo>)</mo>
</mrow>
</mrow>
εk+1＝εk-λ(Axk-ΔT)
其中：xk和xk+1分别为第k次迭代和第k+1次迭代的时间修正序列，μk和μk+1分别为第k次迭代和第k+1次迭代的离散梯度向量，和为权重系数向量θi对应第k次迭代和第k+1次迭代的结果，εk和εk+1为权重系数向量ε对应第k次迭代和第k+1次迭代的结果，为离散梯度向量μk中的第i个元素，k为迭代次数。
7.根据权利要求6所述的PET系统时间修正方法，其特征在于：所述的时间修正序列xk+1通过以下迭代方程求解：
xk+1＝xk-αkg(xk)
<mrow>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mi>k</mi>
</msup>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mo>{</mo>
<msubsup>
<mi>&beta;D</mi>
<mi>i</mi>
<mi>T</mi>
</msubsup>
<mrow>
<mo>(</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<msup>
<mi>x</mi>
<mi>k</mi>
</msup>
<mo>)</mo>
<mo>-</mo>
<msubsup>
<mi>&mu;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msubsup>
<mi>D</mi>
<mi>i</mi>
<mi>T</mi>
</msubsup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&theta;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mo>}</mo>
<mo>+</mo>
<msup>
<mi>&lambda;A</mi>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<msup>
<mi>Ax</mi>
<mi>k</mi>
</msup>
<mo>-</mo>
<mi>&Delta;</mi>
<mi>T</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msup>
<mi>A</mi>
<mi>T</mi>
</msup>
<msup>
<mi>&epsiv;</mi>
<mi>k</mi>
</msup>
</mrow>
其中：αk为第k次迭代的求解步长，其通过Amijo线性搜索获取，g(xk)为关于时间修正序列xk的梯度函数。
8.根据权利要求6所述的PET系统时间修正方法，其特征在于：所述的离散梯度向量μk+1通过以下迭代方程求解：
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>=</mo>
<mi>max</mi>
<mrow>
<mo>(</mo>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>)</mo>
<mo>-</mo>
<mfrac>
<msubsup>
<mi>&theta;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mi>&beta;</mi>
</mfrac>
<mo>|</mo>
<mo>|</mo>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mi>&beta;</mi>
</mfrac>
<mo>,</mo>
<mn>0</mn>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<mfrac>
<mrow>
<mo>(</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>)</mo>
<mo>-</mo>
<mfrac>
<msubsup>
<mi>&theta;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mi>&beta;</mi>
</mfrac>
<mo>)</mo>
</mrow>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>D</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<msubsup>
<mi>&theta;</mi>
<mi>i</mi>
<mi>k</mi>
</msubsup>
<mi>&beta;</mi>
</mfrac>
<mo>|</mo>
<mo>|</mo>
</mrow>
</mfrac>
<mo>.</mo>
</mrow>
9.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述系统矩阵A的维度为M×N，M＝(N-1)N/2，系统矩阵A的具体表达如下：
<mrow>
<mi>A</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mn>2</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mrow>
<mi>N</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中：A1～AN-1均为系统矩阵A的子矩阵，对于任一子矩阵Aj，其维度为(N-j)×N，j为自然数且1≤j≤N-1；所述的子矩阵Aj中第j列向量的所有元素均为1，前j-1列向量的所有元素均为0，后N-j列向量所组成的方阵为主对角线元素均为-1的对角矩阵。
10.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(5)中时间修正序列x和y进行融合的具体方法如下：
5.1对于PET系统中的任一探测器，从时间修正序列x提取该探测器各晶体单元对应的n个时间修正值x1～xn，从时间修正序列y提取该探测器各晶体单元对应的n个时间修正值y1～yn；
5.2计算时间修正值y1～yn的平均值对于时间修正值x1～xn中的任一时间修正值xp，根据以下算式将其扩展成n个时间修正值xp1～xpn；依此遍历时间修正序列x中每个探测器晶体单元对应的时间修正值，从而使时间修正序列x扩展成一个nN维的时间修正序列x*；
<mrow>
<msub>
<mi>x</mi>
<mrow>
<mi>p</mi>
<mi>q</mi>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<msub>
<mi>x</mi>
<mi>p</mi>
</msub>
<msub>
<mi>y</mi>
<mi>q</mi>
</msub>
</mrow>
<mover>
<mi>y</mi>
<mo>&OverBar;</mo>
</mover>
</mfrac>
</mrow>
其中：p和q均为自然数且1≤p≤n，1≤q≤n，yq为时间修正值y1～yn中的任一时间修正值；
5.3计算时间修正值x1～xn的平均值对于时间修正值y1～yn中的任一时间修正值yq，根据以下算式将其扩展成n个时间修正值yq1～yqn；依此遍历时间修正序列y中每个探测器晶体单元对应的时间修正值，从而使时间修正序列y扩展成一个nN维的时间修正序列y*；
<mrow>
<msub>
<mi>y</mi>
<mrow>
<mi>q</mi>
<mi>p</mi>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<msub>
<mi>y</mi>
<mi>q</mi>
</msub>
<msub>
<mi>x</mi>
<mi>p</mi>
</msub>
</mrow>
<mover>
<mi>x</mi>
<mo>&OverBar;</mo>
</mover>
</mfrac>
</mrow>
5.4将时间修正序列x*和y*相加后除以2，即得到时间修正序列z。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510526548.4A CN105212957B (zh) | 2015-08-25 | 2015-08-25 | 一种基于TV Merge的晶体级PET系统时间修正方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510526548.4A CN105212957B (zh) | 2015-08-25 | 2015-08-25 | 一种基于TV Merge的晶体级PET系统时间修正方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105212957A CN105212957A (zh) | 2016-01-06 |
| CN105212957B true CN105212957B (zh) | 2018-01-16 |
Family
ID=54982501
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510526548.4A Active CN105212957B (zh) | 2015-08-25 | 2015-08-25 | 一种基于TV Merge的晶体级PET系统时间修正方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105212957B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108320793B (zh) * | 2018-01-17 | 2020-08-28 | 江苏赛诺格兰医疗科技有限公司 | 确定扫描时长的方法和装置 |
| CN111965691B (zh) * | 2020-09-14 | 2022-12-23 | 明峰医疗系统股份有限公司 | 一种pet中的时间游走修正方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102755172A (zh) * | 2011-04-28 | 2012-10-31 | 株式会社东芝 | 核医学成像方法以及核医学成像装置 |
| CN103065313A (zh) * | 2012-12-27 | 2013-04-24 | 清华大学 | 一种晶体位置表建立方法 |
| CN103810731A (zh) * | 2014-01-20 | 2014-05-21 | 浙江大学 | 一种基于tv范数的pet图像重建方法 |
| CN104851080A (zh) * | 2015-05-08 | 2015-08-19 | 浙江大学 | 一种基于tv的三维pet图像重建方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6852978B2 (en) * | 2002-10-31 | 2005-02-08 | General Electric Company | Crystal-based coincidence timing calibration method |
| US7332721B2 (en) * | 2005-04-13 | 2008-02-19 | Photodetection Systems, Inc. | Separation of geometric system response matrix for three-dimensional image reconstruction |
| EP2024902A4 (en) * | 2006-02-13 | 2012-06-13 | Univ Chicago | IMAGE RECONSTRUCTION FROM INCOMPLETE OR LIMITED NUMBERS |
| JP5360418B2 (ja) * | 2010-01-27 | 2013-12-04 | 株式会社島津製作所 | 放射線断層撮影装置 |
-
2015
- 2015-08-25 CN CN201510526548.4A patent/CN105212957B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102755172A (zh) * | 2011-04-28 | 2012-10-31 | 株式会社东芝 | 核医学成像方法以及核医学成像装置 |
| CN103065313A (zh) * | 2012-12-27 | 2013-04-24 | 清华大学 | 一种晶体位置表建立方法 |
| CN103810731A (zh) * | 2014-01-20 | 2014-05-21 | 浙江大学 | 一种基于tv范数的pet图像重建方法 |
| CN104851080A (zh) * | 2015-05-08 | 2015-08-19 | 浙江大学 | 一种基于tv的三维pet图像重建方法 |
Non-Patent Citations (2)
| Title |
|---|
| Accurate EM-TV Algorithm in PET with Low SNR;Alex Sawatzky,et al;《2008 IEEE Nuclear Science Symposium Conference Record》;20081231;第5133-5137页 * |
| Time-of-Flight Image Reconstruction with TV Minization Constranit for a Dual-Head Small Animal PET System;Hung-Yi Chou,et al;《2012IEEE Nuclear Science Symposium and Medical Imaging Conference Record》;20121231;第3443-3446页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105212957A (zh) | 2016-01-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10765382B2 (en) | Method for mixed tracers dynamic PET concentration image reconstruction based on stacked autoencoder | |
| Green et al. | High resolution PET, SPECT and projection imaging in small animals | |
| US10360699B2 (en) | Correcting count loss | |
| CN106108934A (zh) | 多伽马光子同时发射药物时间符合核医学成像系统及方法 | |
| CN102007430A (zh) | 放射线断层摄影装置 | |
| CN108986916B (zh) | 基于栈式自编码器的动态pet图像示踪剂动力学宏参数估计方法 | |
| CN109993808B (zh) | 一种基于dsn的动态双示踪pet重建方法 | |
| CN109683188A (zh) | 一种契连柯夫事件与伽马事件符合成像装置和方法 | |
| US9851456B2 (en) | TOF-PET tomograph and a method of imaging using a TOF-PET tomograph, based on a probability of production and lifetime of a positronium | |
| CN104700438A (zh) | 图像重建方法及装置 | |
| CN103295207A (zh) | 一种基于h无穷滤波的双示踪剂pet浓度的动态重建方法 | |
| CN102938154B (zh) | 一种基于粒子滤波的动态pet图像重建方法 | |
| CN103810731A (zh) | 一种基于tv范数的pet图像重建方法 | |
| CN105212957B (zh) | 一种基于TV Merge的晶体级PET系统时间修正方法 | |
| Qing et al. | Separation of dual-tracer PET signals using a deep stacking network | |
| CN105212956B (zh) | 一种基于ist的次晶体级pet系统时间修正方法 | |
| CN103400403B (zh) | 一种pet浓度与衰减系数的同时重建方法 | |
| CN103230282A (zh) | 一种pet浓度均值与方差的估计方法及系统 | |
| CN109009198A (zh) | 多模态成像系统、方法和存储介质 | |
| Williams | Anniversary paper: nuclear medicine: fifty years and still counting | |
| CN104598356B (zh) | 一种事件排序方法及装置 | |
| CN104146726B (zh) | Pet系统符合探测响应的生成方法 | |
| CN106859686A (zh) | 成像方法和成像系统 | |
| CN106388845A (zh) | 一种正电子发射切伦科夫-伽玛双辐射的成像方法与装置 | |
| Naunheim et al. | Improving the timing resolution of positron emission tomography detectors using boosted learning—a residual physics approach |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |