CN110269638B - 图像重建方法、系统、可读存储介质和设备 - Google Patents
图像重建方法、系统、可读存储介质和设备 Download PDFInfo
- Publication number
- CN110269638B CN110269638B CN201910555318.9A CN201910555318A CN110269638B CN 110269638 B CN110269638 B CN 110269638B CN 201910555318 A CN201910555318 A CN 201910555318A CN 110269638 B CN110269638 B CN 110269638B
- Authority
- CN
- China
- Prior art keywords
- image
- kernel function
- matrix
- acquiring
- projection
- 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
- 238000000034 method Methods 0.000 title claims abstract description 45
- 239000011159 matrix material Substances 0.000 claims abstract description 127
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 50
- 230000006870 function Effects 0.000 claims description 165
- 230000004044 response Effects 0.000 claims description 55
- 238000007476 Maximum Likelihood Methods 0.000 claims description 32
- 238000012937 correction Methods 0.000 claims description 13
- 230000008569 process Effects 0.000 claims description 10
- 238000002059 diagnostic imaging Methods 0.000 abstract description 2
- 238000004422 calculation algorithm Methods 0.000 description 5
- 239000013078 crystal Substances 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 4
- 238000003384 imaging method Methods 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 239000000700 radioactive tracer Substances 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000002591 computed tomography Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000009206 nuclear medicine Methods 0.000 description 1
- 238000002600 positron emission tomography Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/037—Emission tomography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/42—Arrangements for detecting radiation specially adapted for radiation diagnosis
- A61B6/4208—Arrangements for detecting radiation specially adapted for radiation diagnosis characterised by using a particular type of detector
- A61B6/4241—Arrangements for detecting radiation specially adapted for radiation diagnosis characterised by using a particular type of detector using energy resolving detectors, e.g. photon counting
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/42—Arrangements for detecting radiation specially adapted for radiation diagnosis
- A61B6/4208—Arrangements for detecting radiation specially adapted for radiation diagnosis characterised by using a particular type of detector
- A61B6/4258—Arrangements for detecting radiation specially adapted for radiation diagnosis characterised by using a particular type of detector for detecting non x-ray radiation, e.g. gamma radiation
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Medical Informatics (AREA)
- Engineering & Computer Science (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Physics & Mathematics (AREA)
- High Energy & Nuclear Physics (AREA)
- Heart & Thoracic Surgery (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Image Processing (AREA)
Abstract
本发明涉及一种图像重建方法、系统、可读存储介质和设备，属于医疗影像技术领域，其方案是将飞行时间核函数进行多级分解，得到多个子函数，利用多个子函数对初始化图像进行正向投影，获得正投影图像，并根据参考图像、正投影图像和多个子函数进一步获取反投影图像；根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像。在此方案中，卷积的运算复杂度与飞行时间核函数的矩阵大小有关，核函数的半高宽直接决定了飞行时间核函数的矩阵大小，用多个子函数对图像进行多层卷积运算，其半高宽相比于直接用飞行时间核函数的半高宽小，因而卷积的运算复杂度也小，可以明显提升运算速度。
Description
技术领域
本发明涉及医疗影像技术领域，特别是涉及一种图像重建方法、系统、可读存储介质和设备。
背景技术
PET(Positron Emission Tomography，正电子发射计算机断层扫描)是核医学领域中较为先进的临床检查影像技术，目前已被广泛应用于医学领域的诊断和研究。
在通过PET系统对生物体进行扫描前，先给生物体注射含有放射性核素的示踪剂，示踪剂的放射性核素在生物体内会发生衰变并产生正电子，接着衰变后产生的正电子在行进十分之几毫米到几毫米后，与生物体内的负电子相遇，发生正负电子对湮灭反应，从而生成一对方向相反、能量相同的伽马光子，这一对伽马光子穿过生物体组织，PET系统的探测器接收这一对伽马光子并分别记录它们到达探测器的时间，当它们的时间差小于给定的时间差(即符合窗)，则认为这对伽马光子是有效的，湮灭事件发生在这一对伽马光子的连线上，即响应线。
在TOF-PET中，探测器能够分辨的伽马光子飞行时间差远远小于符合窗，如目前临床用PET符合窗为4ns，而TOF时间分辨率可以达到600ps。TOF(Time of flight，飞行时间)信息有效提高了对于湮灭事件的定位精度。在图像重建中利用TOF信息可以提升图像信噪比、增强对于数据的鲁棒性以及获得更快的迭代收敛速度。
TOF重建的缺点是速度较慢。探测器在记录湮灭事件的TOF信息时，只能定位中心区域(即Tbin)，而湮灭事件发生的概率分布为高斯函数，其展宽远大于中心区域。以600ps系统时间分辨率为例，按±3σ计算单个Tbin的高斯函数展宽在图像上覆盖范围为46cm。通常600ps时间分辨率对应着10～15个Tbin，由此可计算出总的覆盖范围为460～690cm，远远大于PET扫描视野的覆盖范围(60cm)，所以含TOF信息的重建比无TOF信息的重建更花时间。而随着TOF分辨率的提升，这个问题可能会变的更加严重。
为了解决这个问题，传统技术提出了一种快速TOF重建算法。它的核心思想是将TOF核函数矩阵化，通过图像域的卷积来实现TOF重建中的正反投影运算。该算法避免了反复调用耗时的射线追踪算法，可以有效地提升重建速度。该算法在TOF分辨率高的时候优势非常明显，因为TOF分辨率高意味着TOF核函数矩阵维度较小，它与三维图像的卷积运算速度很快；但TOF分辨率较差时，TOF核函数矩阵维度较大，而卷积的运算量与矩阵维度的立方成正比，这导致运算速度明显降低。
发明内容
基于此，有必要针对传统的TOF重建算法在TOF分辨率较差时，运算速度低的问题，提供一种图像重建方法、系统、可读存储介质和设备。
一种图像重建方法，包括以下步骤：
获取飞行时间核函数，对飞行时间核函数进行多级分解，获得多个子函数；
获取初始化图像，根据多个子函数对初始化图像进行正向投影，获得正投影图像；
获取参考图像，根据参考图像、正投影图像和多个子函数获取反投影图像；
根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像。
根据上述的图像重建方法，将飞行时间核函数进行多级分解，得到多个子函数，利用多个子函数对初始化图像进行正向投影，获得正投影图像，并根据参考图像、正投影图像和多个子函数进一步获取反投影图像；根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像。在此方案中，卷积的运算复杂度与飞行时间核函数的矩阵大小有关，核函数的半高宽直接决定了飞行时间核函数的矩阵大小，用多个子函数对图像进行多层卷积运算，其半高宽相比于直接用飞行时间核函数的半高宽小，因而卷积的运算复杂度也小，可以明显提升运算速度。
在其中一个实施例中，获取飞行时间核函数的步骤包括以下步骤：
获取零角度下的高斯核函数，其中，零角度表示符合响应线的投影角度为零；
根据符合响应线的不同非零投影角度和高斯核函数获取飞行时间核函数。
在其中一个实施例中，对飞行时间核函数进行多级分解，获得多个子函数的步骤包括以下步骤：
根据分解级数和高斯核函数获取零角度下的分解核函数；
根据符合响应线的不同非零投影角度和分解核函数获取多个子函数。
在其中一个实施例中，获取参考图像的步骤包括以下步骤：
获取三维参考图像矩阵，三维参考图像矩阵的前两维图像矩阵与初始化图像的图像矩阵等同，三维参考图像矩阵的第三维表示按投影角度对符合响应线进行归类的总的角度数；
若初始化图像的图像矩阵是二维图像矩阵，根据符合响应线的符合数据获取极大似然点位置，将极大似然点位置添加至三维参考图像矩阵的对应位置，得到参考图像；
或者，
若初始化图像的图像矩阵是三维图像矩阵，获取四维参考图像矩阵，四维参考图像矩阵的前三维图像矩阵与初始化图像的图像矩阵等同，四维参考图像矩阵的第四维表示按投影角度对符合响应线进行归类的总的角度数；
根据符合响应线的符合数据获取极大似然点位置，将极大似然点位置添加至四维参考图像矩阵的对应位置，得到参考图像。
在其中一个实施例中，图像重建方法还包括以下步骤：
获取记录的符合事件，将符合事件的原始记录坐标转换为笛卡尔坐标系下的第一坐标和极坐标系下的第二坐标；
根据符合响应线的符合数据获取极大似然点位置的步骤包括以下步骤：
根据第一坐标和第二坐标获取极大似然点位置。
在其中一个实施例中，根据参考图像、正投影图像和多个子函数获取反投影图像的步骤包括以下步骤：
根据参考图像和正投影图像的矩阵比值获取校正因子，根据校正因子和多个子函数的卷积获取反投影图像。
在其中一个实施例中，根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像的步骤包括以下步骤：
根据反投影图像和飞行时间核函数获取迭代因子，根据迭代因子对初始化图像进行修正，获得待迭代图像；
将初始化图像替换成待迭代图像，根据迭代因子对待迭代图像进行修正，直至达到预设迭代条件，获得重建图像。
一种图像重建系统，包括：
函数分解单元，用于获取飞行时间核函数，对飞行时间核函数进行多级分解，获得多个子函数；
图像正投影单元，用于获取初始化图像，根据多个子函数对初始化图像进行正向投影，获得正投影图像；
图像反投影单元，用于获取参考图像，根据参考图像、正投影图像和多个子函数获取反投影图像；
迭代重建单元，用于根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像。
根据上述的图像重建系统，函数分解单元将飞行时间核函数进行多级分解，得到多个子函数，图像正投影单元利用多个子函数对初始化图像进行正向投影，获得正投影图像，图像反投影单元根据参考图像、正投影图像和多个子函数进一步获取反投影图像；迭代重建单元根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像。在此方案中，卷积的运算复杂度与飞行时间核函数的矩阵大小有关，核函数的半高宽直接决定了飞行时间核函数的矩阵大小，用多个子函数对图像进行多层卷积运算，其半高宽相比于直接用飞行时间核函数的半高宽小，因而卷积的运算复杂度也小，可以明显提升运算速度。
在其中一个实施例中，函数分解单元用于获取零角度下的高斯核函数，根据符合响应线的不同非零投影角度和高斯核函数获取飞行时间核函数；其中，零角度表示符合响应线的投影角度为零。
在其中一个实施例中，函数分解单元用于根据分解级数和高斯核函数获取零角度下的分解核函数；根据符合响应线的不同非零投影角度和分解核函数获取多个子函数。
在其中一个实施例中，图像反投影单元用于在初始化图像的图像矩阵是二维图像矩阵时，获取三维参考图像矩阵，三维参考图像矩阵的前两维图像矩阵与初始化图像的图像矩阵等同，三维参考图像矩阵的第三维表示按投影角度对符合响应线进行归类的总的角度数；根据符合响应线的符合数据获取极大似然点位置，将极大似然点位置添加至三维参考图像矩阵的对应位置，得到参考图像；
或者，
图像反投影单元用于在初始化图像的图像矩阵是三维图像矩阵时，获取四维参考图像矩阵，四维参考图像矩阵的前三维图像矩阵与初始化图像的图像矩阵等同，四维参考图像矩阵的第四维表示按投影角度对符合响应线进行归类的总的角度数；根据符合响应线的符合数据获取极大似然点位置，将极大似然点位置添加至四维参考图像矩阵的对应位置，得到参考图像。
在其中一个实施例中，图像重建系统还包括坐标获取单元，用于获取记录的符合事件，将符合事件的原始记录坐标转换为笛卡尔坐标系下的第一坐标和极坐标系下的第二坐标；
图像反投影单元用于根据第一坐标和第二坐标获取极大似然点位置。
在其中一个实施例中，图像反投影单元用于根据参考图像和正投影图像的矩阵比值获取校正因子，根据校正因子和多个子函数的卷积获取反投影图像。
在其中一个实施例中，迭代重建单元用于根据反投影图像和飞行时间核函数获取迭代因子，根据迭代因子对初始化图像进行修正，获得待迭代图像；将初始化图像替换成待迭代图像，根据迭代因子对待迭代图像进行修正，直至达到预设迭代条件，获得重建图像。
一种可读存储介质，其上存储有可执行程序，可执行程序被处理器执行时实现上述的图像重建方法的步骤。
上述可读存储介质，通过其存储的可执行程序，可以实现将飞行时间核函数进行多级分解，得到多个子函数，卷积的运算复杂度与飞行时间核函数的矩阵大小有关，核函数的半高宽直接决定了飞行时间核函数的矩阵大小，用多个子函数对图像进行多层卷积运算，其半高宽相比于直接用飞行时间核函数的半高宽小，因而卷积的运算复杂度也小，可以明显提升运算速度。
一种图像重建设备，包括存储器和处理器，存储器存储有可执行程序，处理器执行可执行程序时实现上述的图像重建方法的步骤。
上述图像重建设备，通过在处理器上运行可执行程序，可以实现将飞行时间核函数进行多级分解，得到多个子函数，卷积的运算复杂度与飞行时间核函数的矩阵大小有关，核函数的半高宽直接决定了飞行时间核函数的矩阵大小，用多个子函数对图像进行多层卷积运算，其半高宽相比于直接用飞行时间核函数的半高宽小，因而卷积的运算复杂度也小，可以明显提升运算速度。
附图说明
图1为一个实施例中的图像重建方法的流程示意图；
图2为一个实施例中的飞行时间核函数多级分解后，利用多层卷积进行图像的正投影示意图；
图3为一个实施例中的图像重建系统的结构示意图；
图4为另一个实施例中的图像重建系统的结构示意图。
具体实施方式
为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。
需要说明的是，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域技术人员在没有付出创造性劳动前提下所获的所有其他实施例，都属于本发明保护的范围。
在本发明实施例中使用的术语是仅仅出于描述特定实施例的目的，而非旨在限制本发明。在本发明实施例和所附权利要求书中所使用的单数形式的“一种”、“所述”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义。
本申请提供的图像重建方法，可以应用于PET扫描成像的应用场景中。
参见图1所示，为本发明一个实施例的图像重建方法的流程示意图。该实施例中的图像重建方法包括以下步骤：
步骤S110：获取飞行时间核函数，对飞行时间核函数进行多级分解，获得多个子函数；
在本步骤中，飞行时间核函数可以表示湮灭事件发生的概率，飞行时间核函数可以是高斯函数、指数函数等中心不低于边缘的函数，可以通过傅里叶变换对飞行时间核函数进行多级分解；
步骤S120：获取初始化图像，根据多个子函数对初始化图像进行正向投影，获得正投影图像；
在本步骤中，初始图像是进行迭代处理前最初使用的图像，可以是PET扫描系统中默认的图像；在进行正向投影时，可以将多个子函数与初始化图像的像素进行卷积运算，得到正投影图像；
步骤S130：获取参考图像，根据参考图像、正投影图像和多个子函数获取反投影图像；
在本步骤中，PET扫描系统的探测器接收光子以后，通过光电转换可以得到相应的电信号，对电信号进行数据处理后，可以得到PET扫描的符合事件数据，参考图像包括符合事件数据，在进行反向投影时，可以将多个子函数与参考图像、正投影图像进行卷积运算，得到反投影图像；
步骤S140：根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像。
在本步骤中，在飞行时间核函数的基础上对初始化图像进行迭代更新，可以将飞行时间所包含的位置信息准确显示在PET重建图像中。
根据上述的图像重建方法，将飞行时间核函数进行多级分解，得到多个子函数，利用多个子函数对初始化图像进行正向投影，获得正投影图像，并根据参考图像、正投影图像和多个子函数进一步获取反投影图像；根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像。在此方案中，卷积的运算复杂度与飞行时间核函数的矩阵大小有关，核函数的半高宽直接决定了飞行时间核函数的矩阵大小，用多个子函数对图像进行多层卷积运算，其半高宽相比于直接用飞行时间核函数的半高宽小，因而卷积的运算复杂度也小，可以明显提升运算速度。
进一步的，如飞行时间核函数为高斯函数可以将其分解成M个子函数将M个子函数应用于正向投影图像和反向投影图像的获取过程。
在一个实施例中，获取飞行时间核函数的步骤包括以下步骤：
获取零角度下的高斯核函数，其中，零角度表示符合响应线的投影角度为零；
根据符合响应线的不同非零投影角度和高斯核函数获取飞行时间核函数。
在本实施例中，飞行时间核函数可以是高斯函数，其可用零角度下的高斯核函数来获取，零角度表示符合响应线的投影角度为零，不同符合响应线的投影角度会有所不同，这一投影角度是飞行时间核函数的参数之一，利用符合响应线的不同非零角度可以将零角度下的高斯核函数转换为飞行时间核函数，便于飞行时间核函数的卷积计算。
具体的，可以表示角度/>下的飞行时间核函数，/>可以表示零角度下的高斯核函数，符合响应线的投影角度/>可以是符合响应线与系统坐标x轴之间的夹角；可以通过以下公式将零角度下的高斯核函数转换为飞行时间核函数：
在一个实施例中，对飞行时间核函数进行多级分解，获得多个子函数的步骤包括以下步骤：
根据分解级数和高斯核函数获取零角度下的分解核函数；
根据符合响应线的不同非零投影角度和分解核函数获取多个子函数。
在本实施例中，多级分解后的多个子函数依然包括投影角度这一参数，多个子函数的卷积与飞行时间核函数等效，因此可以先利用分解级数对高斯核函数进行分解，得到零角度下的分解核函数，再结合符合响应线的不同非零投影角度得到多个子函数，如此可以简化飞行时间核函数的分解过程，提高运算速度。
具体的，如htof(0)可以表示零角度下的高斯核函数，如分解级数为M，得到零角度下的分解核函数利用符合响应线的不同非零角度将分解核函数/>转换为子函数/>
在一个实施例中，获取参考图像的步骤包括以下步骤：
若初始化图像的图像矩阵是二维图像矩阵，获取三维参考图像矩阵，三维参考图像矩阵的前两维图像矩阵与初始化图像的图像矩阵等同，三维参考图像矩阵的第三维表示按投影角度对符合响应线进行归类的总的角度数；
根据符合响应线的符合数据获取极大似然点位置，将极大似然点位置添加至三维参考图像矩阵的对应位置，得到参考图像；
或者，
若初始化图像的图像矩阵是三维图像矩阵，获取四维参考图像矩阵，四维参考图像矩阵的前三维图像矩阵与初始化图像的图像矩阵等同，四维参考图像矩阵的第四维表示按投影角度对符合响应线进行归类的总的角度数；
根据符合响应线的符合数据获取极大似然点位置，将极大似然点位置添加至四维参考图像矩阵的对应位置，得到参考图像。
在本实施例中，初始化图像可以是二维或三维的，若初始化图像是二维的，三维参考图像矩阵中的前两维图像矩阵与初始化图像的图像矩阵等同，建立三维参考图像矩阵与初始化图像的图像矩阵之间的联系，三维参考图像矩阵中的第三维表示按投影角度对符合响应线进行归类的总的角度数，与符合响应线的投影角度相对应；利用符合响应线的符合数据获取极大似然点位置，将极大似然点位置添加至三维参考图像矩阵的对应位置，可以将符合事件数据添加至参考图像，在后续迭代过程中提供数据支持；若初始化图像可以是三维的，情况与二维的类似，可以将符合事件数据添加至参考图像，以支持三维成像。
在一个实施例中，图像重建方法还包括以下步骤：
获取记录的符合事件，将符合事件的原始记录坐标转换为笛卡尔坐标系下的第一坐标和极坐标系下的第二坐标；
根据符合响应线的符合数据获取极大似然点位置的步骤包括以下步骤：
根据第一坐标和第二坐标获取极大似然点位置。
在本实施例中，PET扫描系统在对符合事件进行记录时，一般是采用列表模式记录符合事件的原始记录坐标，原始记录坐标一般是检测光子的探测器晶体和飞行时间有关的信息，无法直接用于数据计算，因此需要将其转换成笛卡尔坐标系下的第一坐标和极坐标下的第二坐标，通过第一坐标和第二坐标的计算可以得到湮灭事件发生的极大似然点位置，实现符合数据的正确处理。
具体的，原始记录坐标可以是ia,ib,ta,tb等，其中ia,ib表示检测到符合事件的一对探测器晶体编号，ta,tb表示检测到符合事件的一对探测器记录的伽马光子的飞行时间信息，笛卡尔坐标系下的第一坐标可以是xa，ya，xb，yb，ta，tb，其中，xa，ya，xb，yb表示一对探测器晶体的物理位置，极坐标系下的第二坐标可以是s，dt，其中，s表示符合响应线到坐标系中心的距离，/>表示符合响应线与坐标系x轴的夹角，dt表示符合事件的飞行时间信息，dt＝ta-tb。
在一个实施例中，根据参考图像、正投影图像和多个子函数获取反投影图像的步骤包括以下步骤：
根据参考图像和正投影图像的矩阵比值获取校正因子，根据校正因子和多个子函数的卷积获取反投影图像。
在本实施例中，将参考图像和正投影图像的矩阵比值作为校正因子，并与多个子函数进行卷积运算，由于使用多个子函数的卷积的运算复杂度小，因此可以快速获得反投影图像。
在一个实施例中，根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像的步骤包括以下步骤：
根据反投影图像和飞行时间核函数获取迭代因子，根据迭代因子对初始化图像进行修正，获得待迭代图像；
将初始化图像替换成待迭代图像，根据迭代因子对待迭代图像进行修正，直至达到预设迭代条件，获得重建图像。
在本实施例中，利用反投影图像和飞行时间核函数得到迭代因子，对初始化图像进行修正，获得待迭代图像，可以在迭代过程中不断完善符合事件的数据，由于迭代因子中涉及多个子函数的卷积运算，因此，可以加快图像迭代最终获得重建图像的过程。
将飞行时间核函数分解为多个子函数进行卷积运算，这种方式可以减少运算时间的原理如下：
TOF核函数为高斯函数，为
其对应的Fourier变换为：
简单地推导可得：
即一个高斯函数可以写成两个高斯函数卷积的形式，其标准差满足：
高斯函数的半高宽与标准差的关系为：
则
更一般的结论是：一个高斯函数可以写成多个高斯函数卷积的形式，其半高宽满足：
高斯函数的半高宽直接决定了TOF核函数的矩阵大小。
卷积的运算复杂度与TOF核函数的矩阵大小有关，对于一个维度为D(σ)×D(σ)×D(σ)的矩阵，其运算复杂度为：
k·D3(σ)
其中k为常数，表示投影运算中矩阵卷积操作的次数；
将该矩阵分解为M个相同的矩阵卷积后，其运算复杂度为：
从上式可以看出采用M层卷积后，原卷积的运算复杂度下降为原来的
举例来说，假定一个PET系统的TOF时间分辨率为600ps，它可以分解为两个424ps的TOF核函数(即子函数)(如图2所示)，或者3个346ps的TOF核函数，或者4个300ps的TOF核函数，它们的卷积与原核函数等价，而总运算时间分别降低为原来的70％，58％和50％。
在一个实施例中，以图像域上进行二维TOF重建为例，其实现步骤为：
1、确定二维重建图像矩阵fN×N,及像素大小d×d；
2、确定三维参考图像矩阵前两维图像矩阵及像素大小等同于重建图像，第三维表示按角度对响应线进行归类的总的角度数；
3、将列表模式记录的符合事件进行坐标转换，由原始记录的(ia,ib,ta,tb)通过查表转为笛卡尔坐标系下的(xa,ya,xb,yb,ta,tb)和极坐标系下的按角度将符合事件进行归类；
其中：
ia,ib表示一对探测器晶体编号；
ta,tb表示一对探测器记录的伽马光子的飞行时间信息；
xa,ya,xb,yb表示一对探测器晶体的物理位置；
s表示响应线到坐标系中心的距离；
表示响应线与坐标系x轴的夹角；
dt表示符合事件的飞行时间信息，即dt＝ta-tb；
4、在相同角度下，对于每一个符合事件计算其MLP(Maximum likelihood point，极大似然点)位置并累加到参考图像中对应位置；遍历所有数据后得到完整的参考图像；
其中MLP计算公式为：
5、初始化重建图像，将所有像素值赋为1；
6、对重建图像进行正向投影：
7、计算校正因子：
8、反向投影，得到
9、重复步骤6-8直至遍历当前子集sk中所有角度，得到子集是对归类后的数据进行划分得到的；
10、图像迭代更新：
其中矩阵I的大小等同于重建图像，像素值均为1；矩阵I是所有可能的响应线的MLP点的集合。
重复步骤6-10直到达到预设条件。常见的预设条件包括迭代次数，或者图像更新前后的差异小于某个阈值。
上述步骤中表示在角度/>下的TOF核函数，计算公式为:
其中htof(0；x,y)为零角度下的高斯核函数，它是连续高斯函数在数字图像上的离散化形式：
根据本方案的原理分析可以对TOF核函数进行M级分解，得到：
相应的：
从而重建步骤6更改为：
重建步骤8更改为：
根据以上实现方案，本领域技术人员可以直接推出到图像三维TOF重建方案以及含物理校正项的重建方案。在此不加赘述。
根据上述图像重建方法，本发明实施例还提供一种图像重建系统，以下就图像重建系统的实施例进行详细说明。
参见图3所示，为一个实施例的图像重建系统的结构示意图。该实施例中的图像重建系统包括：
函数分解单元210，用于获取飞行时间核函数，对飞行时间核函数进行多级分解，获得多个子函数；
图像正投影单元220，用于获取初始化图像，根据多个子函数对初始化图像进行正向投影，获得正投影图像；
图像反投影单元230，用于获取参考图像，根据参考图像、正投影图像和多个子函数获取反投影图像；
迭代重建单元240，用于根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像。
在本实施例中，函数分解单元210将飞行时间核函数进行多级分解，得到多个子函数，图像正投影单元220利用多个子函数对初始化图像进行正向投影，获得正投影图像，图像反投影单元230根据参考图像、正投影图像和多个子函数进一步获取反投影图像；迭代重建单元240根据反投影图像和飞行时间核函数对初始化图像进行迭代更新，获取重建图像。在此方案中，卷积的运算复杂度与飞行时间核函数的矩阵大小有关，核函数的半高宽直接决定了飞行时间核函数的矩阵大小，用多个子函数对图像进行多层卷积运算，其半高宽相比于直接用飞行时间核函数的半高宽小，因而卷积的运算复杂度也小，可以明显提升运算速度。
在一个实施例中，函数分解单元210用于获取零角度下的高斯核函数，根据符合响应线的不同非零投影角度和高斯核函数获取飞行时间核函数；其中，零角度表示符合响应线的投影角度为零。
在一个实施例中，函数分解单元210用于根据分解级数和高斯核函数获取零角度下的分解核函数；根据符合响应线的不同非零投影角度和分解核函数获取多个子函数。
在一个实施例中，图像反投影单元230用于在初始化图像的图像矩阵是二维图像矩阵时，获取三维参考图像矩阵，三维参考图像矩阵的前两维图像矩阵与初始化图像的图像矩阵等同，三维参考图像矩阵的第三维表示按投影角度对符合响应线进行归类的总的角度数；根据符合响应线的符合数据获取极大似然点位置，将极大似然点位置添加至三维参考图像矩阵的对应位置，得到参考图像；
或者，
图像反投影单元230用于在初始化图像的图像矩阵是三维图像矩阵时，获取四维参考图像矩阵，四维参考图像矩阵的前三维图像矩阵与初始化图像的图像矩阵等同，四维参考图像矩阵的第四维表示按投影角度对符合响应线进行归类的总的角度数；根据符合响应线的符合数据获取极大似然点位置，将极大似然点位置添加至四维参考图像矩阵的对应位置，得到参考图像。
在一个实施例中，如图4所示，图像重建系统还包括坐标获取单元250，用于获取记录的符合事件，将符合事件的原始记录坐标转换为笛卡尔坐标系下的第一坐标和极坐标系下的第二坐标；
图像反投影单元230用于根据第一坐标和第二坐标获取极大似然点位置。
在一个实施例中，图像反投影单元230用于根据参考图像和正投影图像的矩阵比值获取校正因子，根据校正因子和多个子函数的卷积获取反投影图像。
在一个实施例中，迭代重建单元240用于根据反投影图像和飞行时间核函数获取迭代因子，根据迭代因子对初始化图像进行修正，获得待迭代图像；将初始化图像替换成待迭代图像，根据迭代因子对待迭代图像进行修正，直至达到预设迭代条件，获得重建图像。
本发明实施例的图像重建系统与上述图像重建方法一一对应，在上述图像重建方法的实施例阐述的技术特征及其有益效果均适用于图像重建系统的实施例中。
一种可读存储介质，其上存储有可执行程序，可执行程序被处理器执行时实现上述的图像重建方法的步骤。
上述可读存储介质，通过其存储的可执行程序，可以实现将飞行时间核函数进行多级分解，得到多个子函数，卷积的运算复杂度与飞行时间核函数的矩阵大小有关，核函数的半高宽直接决定了飞行时间核函数的矩阵大小，用多个子函数对图像进行多层卷积运算，其半高宽相比于直接用飞行时间核函数的半高宽小，因而卷积的运算复杂度也小，可以明显提升运算速度。
一种图像重建设备，包括存储器和处理器，存储器存储有可执行程序，处理器执行可执行程序时实现上述的图像重建方法的步骤。
上述图像重建设备，通过在处理器上运行可执行程序，可以实现将飞行时间核函数进行多级分解，得到多个子函数，卷积的运算复杂度与飞行时间核函数的矩阵大小有关，核函数的半高宽直接决定了飞行时间核函数的矩阵大小，用多个子函数对图像进行多层卷积运算，其半高宽相比于直接用飞行时间核函数的半高宽小，因而卷积的运算复杂度也小，可以明显提升运算速度。
本领域普通技术人员可以理解实现上述实施例用于医学图像重建的数据处理方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，程序可存储于一非易失性的计算机可读取存储介质中，如实施例中，该程序可存储于计算机系统的存储介质中，并被该计算机系统中的至少一个处理器执行，以实现包括如上述图像重建方法的实施例的流程。其中，存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。
以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成。所述的程序可以存储于可读取存储介质中。该程序在执行时，包括上述方法所述的步骤。所述的存储介质，包括：ROM/RAM、磁碟、光盘等。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。
Claims (10)
1.一种图像重建方法，其特征在于，包括以下步骤：
获取飞行时间核函数，对所述飞行时间核函数进行多级分解，获得多个子函数；所述飞行时间核函数是中心不低于边缘的函数，用于表示湮灭事件发生的概率；
获取初始化图像，根据多个所述子函数对所述初始化图像进行正向投影，获得正投影图像；
获取参考图像，根据所述参考图像、所述正投影图像和多个所述子函数获取反投影图像；
根据所述反投影图像和所述飞行时间核函数对所述初始化图像进行迭代更新，获取重建图像。
2.根据权利要求1所述的图像重建方法，其特征在于，所述获取飞行时间核函数的步骤包括以下步骤：
获取零角度下的高斯核函数，其中，所述零角度表示符合响应线的投影角度为零；
根据所述符合响应线的不同非零投影角度和所述高斯核函数获取所述飞行时间核函数。
3.根据权利要求2所述的图像重建方法，其特征在于，所述对所述飞行时间核函数进行多级分解，获得多个子函数的步骤包括以下步骤：
根据分解级数和所述高斯核函数获取零角度下的分解核函数；
根据所述符合响应线的不同非零投影角度和所述分解核函数获取多个子函数。
4.根据权利要求3所述的图像重建方法，其特征在于，所述获取参考图像的步骤包括以下步骤：
若所述初始化图像的图像矩阵是二维图像矩阵，获取三维参考图像矩阵，所述三维参考图像矩阵的前两维图像矩阵与所述初始化图像的图像矩阵等同，所述三维参考图像矩阵的第三维表示按投影角度对所述符合响应线进行归类的总的角度数；
根据所述符合响应线的符合数据获取极大似然点位置，将所述极大似然点位置添加至所述三维参考图像矩阵的对应位置，得到所述参考图像；
或者，
若所述初始化图像的图像矩阵是三维图像矩阵，获取四维参考图像矩阵，所述四维参考图像矩阵的前三维图像矩阵与所述初始化图像的图像矩阵等同，所述四维参考图像矩阵的第四维表示按投影角度对所述符合响应线进行归类的总的角度数；
根据所述符合响应线的符合数据获取极大似然点位置，将所述极大似然点位置添加至所述四维参考图像矩阵的对应位置，得到所述参考图像。
5.根据权利要求4所述的图像重建方法，其特征在于，还包括以下步骤：
获取记录的符合事件，将所述符合事件的原始记录坐标转换为笛卡尔坐标系下的第一坐标和极坐标系下的第二坐标；
所述根据所述符合响应线的符合数据获取极大似然点位置的步骤包括以下步骤：
根据所述第一坐标和所述第二坐标获取所述极大似然点位置。
6.根据权利要求4所述的图像重建方法，其特征在于，所述根据所述参考图像、所述正投影图像和多个所述子函数获取反投影图像的步骤包括以下步骤：
根据所述参考图像和所述正投影图像的矩阵比值获取校正因子，根据所述校正因子和多个所述子函数的卷积获取反投影图像。
7.根据权利要求4所述的图像重建方法，其特征在于，所述根据所述反投影图像和所述飞行时间核函数对所述初始化图像进行迭代更新，获取重建图像的步骤包括以下步骤：
根据所述反投影图像和所述飞行时间核函数获取迭代因子，根据所述迭代因子对所述初始化图像进行修正，获得待迭代图像；
将所述初始化图像替换成所述待迭代图像，根据所述迭代因子对所述待迭代图像进行修正，直至达到预设迭代条件，获得所述重建图像。
8.一种图像重建系统，其特征在于，包括：
函数分解单元，用于获取飞行时间核函数，对所述飞行时间核函数进行多级分解，获得多个子函数；所述飞行时间核函数是中心不低于边缘的函数，用于表示湮灭事件发生的概率；
图像正投影单元，用于获取初始化图像，根据多个所述子函数对所述初始化图像进行正向投影，获得正投影图像；
图像反投影单元，用于获取参考图像，根据所述参考图像、所述正投影图像和多个所述子函数获取反投影图像；
迭代重建单元，用于根据所述反投影图像和所述飞行时间核函数对所述初始化图像进行迭代更新，获取重建图像。
9.一种可读存储介质，其上存储有可执行程序，其特征在于，所述可执行程序被处理器执行时实现权利要求1至7中任意一项所述的图像重建方法的步骤。
10.一种图像重建设备，包括存储器和处理器，所述存储器存储有可执行程序，其特征在于，所述处理器执行所述可执行程序时实现权利要求1至7中任意一项所述的图像重建方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910555318.9A CN110269638B (zh) | 2019-06-25 | 2019-06-25 | 图像重建方法、系统、可读存储介质和设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910555318.9A CN110269638B (zh) | 2019-06-25 | 2019-06-25 | 图像重建方法、系统、可读存储介质和设备 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110269638A CN110269638A (zh) | 2019-09-24 |
| CN110269638B true CN110269638B (zh) | 2023-07-25 |
Family
ID=67963180
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910555318.9A Active CN110269638B (zh) | 2019-06-25 | 2019-06-25 | 图像重建方法、系统、可读存储介质和设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110269638B (zh) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109564692A (zh) * | 2016-08-03 | 2019-04-02 | 皇家飞利浦有限公司 | 使用局部修改的飞行时间（tof）内核进行tof pet图像重建 |
Family Cites Families (17)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7227149B2 (en) * | 2004-12-30 | 2007-06-05 | General Electric Company | Method and system for positron emission tomography image reconstruction |
| JP5220617B2 (ja) * | 2006-01-09 | 2013-06-26 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | Ｔｏｆｆｏｖを介するランダム低減 |
| US8000513B2 (en) * | 2008-09-22 | 2011-08-16 | Siemens Medical Solutions Usa, Inc. | System and method for 3D time of flight PET forward projection based on an exact axial inverse rebinning relation in fourier space |
| EP2433261B1 (en) * | 2009-05-20 | 2014-03-05 | Koninklijke Philips N.V. | Continuous time-of-flight scatter simulation method |
| US8903152B2 (en) * | 2012-06-29 | 2014-12-02 | General Electric Company | Methods and systems for enhanced tomographic imaging |
| US9189870B2 (en) * | 2012-07-23 | 2015-11-17 | Mediso Orvosi Berendezes Fejleszto es Szerviz Kft. | Method, computer readable medium and system for tomographic reconstruction |
| WO2014066629A1 (en) * | 2012-10-26 | 2014-05-01 | The Regents Of The University Of California | Image-based point-spread-function modelling in time-of-flight positron-emission-tomography iterative list-mode reconstruction |
| KR101912715B1 (ko) * | 2012-11-20 | 2018-10-29 | 삼성전자주식회사 | 방사선이 방출된 위치의 분포를 추정하는 방법 및 장치 |
| CN103908280B (zh) * | 2013-01-08 | 2017-07-28 | 上海联影医疗科技有限公司 | Pet散射校正的方法 |
| US9474495B2 (en) * | 2014-12-22 | 2016-10-25 | General Electric Company | System and method for joint estimation of attenuation and activity information |
| US10663608B2 (en) * | 2015-09-21 | 2020-05-26 | Shanghai United Imaging Healthcare Co., Ltd. | System and method for calibrating a PET scanner |
| US9990741B2 (en) * | 2015-09-28 | 2018-06-05 | Siemens Medical Solutions Usa, Inc. | Motion correction in a projection domain in time of flight positron emission tomography |
| WO2017184681A1 (en) * | 2016-04-19 | 2017-10-26 | The General Hospital Corporation | Systems and methods for data-driven respiratory gating in positron emission tomography |
| CN106353786B (zh) * | 2016-09-30 | 2019-06-28 | 上海联影医疗科技有限公司 | 正电子发射断层成像系统飞行时间性能检测方法及装置 |
| CN106333702A (zh) * | 2016-09-30 | 2017-01-18 | 上海联影医疗科技有限公司 | 利用正电子发射断层成像系统进行有源模体定位方法 |
| CN107978002B (zh) * | 2016-10-25 | 2021-01-05 | 上海东软医疗科技有限公司 | 一种pet图像重建方法和装置 |
| US10438378B2 (en) * | 2017-08-25 | 2019-10-08 | Uih America, Inc. | System and method for determining an activity map and an attenuation map |
-
2019
- 2019-06-25 CN CN201910555318.9A patent/CN110269638B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109564692A (zh) * | 2016-08-03 | 2019-04-02 | 皇家飞利浦有限公司 | 使用局部修改的飞行时间（tof）内核进行tof pet图像重建 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110269638A (zh) | 2019-09-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| EP3067864B1 (en) | Iterative reconstruction with enhanced noise control filtering | |
| US7312455B2 (en) | Method and system for scatter correction in a positron emission tomography system | |
| US8265365B2 (en) | Time of flight scatter distribution estimation in positron emission tomography | |
| US8457380B2 (en) | PET local tomography | |
| US20070075249A1 (en) | Distributed iterative image reconstruction | |
| EP1949136A1 (en) | Pet imaging using anatomic list mode mask | |
| CN111709897B (zh) | 一种基于域变换的正电子发射断层图像的重建方法 | |
| US11234667B2 (en) | Scatter correction using emission image estimate reconstructed from narrow energy window counts in positron emission tomography | |
| CN107705261B (zh) | 一种图像重建方法和装置 | |
| CN110415310B (zh) | 医学扫描成像方法、装置、存储介质及计算机设备 | |
| CN110415311B (zh) | Pet图像重建方法、系统、可读存储介质和设备 | |
| CN109712213B (zh) | Pet图像重建方法、系统、可读存储介质和设备 | |
| EP3685354A1 (en) | A relaxed iterative maximum-likelihood expectation maximization for positron emission tomography random coincidence estimation | |
| US8359345B2 (en) | Iterative algorithms for variance reduction on compressed sinogram random coincidences in PET | |
| US11164344B2 (en) | PET image reconstruction using TOF data and neural network | |
| US10217250B2 (en) | Multi-view tomographic reconstruction | |
| CN110269638B (zh) | 图像重建方法、系统、可读存储介质和设备 | |
| CN114862980A (zh) | 散射校正方法、pet成像方法、装置、设备及存储介质 | |
| CN110264537B (zh) | Pet图像重建方法、系统、可读存储介质和设备 | |
| CN111080737B (zh) | 图像重建方法、装置及pet扫描系统 | |
| CN110070588B (zh) | Pet图像重建方法、系统、可读存储介质和设备 | |
| CN109636873B (zh) | 用于医学图像重建的数据处理方法和医学图像重建方法 | |
| CN112052885A (zh) | 图像处理方法、装置、设备及pet-ct系统 | |
| Boquet-Pujadas et al. | PET Rebinning with Regularized Density Splines | |
| US11908085B2 (en) | Entropy-dependent adaptive image filtering |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB02 | Change of applicant information |
Address after: 201807 Shanghai City, north of the city of Jiading District Road No. 2258 Applicant after: Shanghai Lianying Medical Technology Co.,Ltd. Address before: 201807 Shanghai City, north of the city of Jiading District Road No. 2258 Applicant before: SHANGHAI UNITED IMAGING HEALTHCARE Co.,Ltd. |
|
| CB02 | Change of applicant information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |