WO2021114306A1 - 一种基于gpu加速的高精度pet重建方法及装置 - Google Patents

一种基于gpu加速的高精度pet重建方法及装置 Download PDF

Info

Publication number
WO2021114306A1
WO2021114306A1 PCT/CN2019/125416 CN2019125416W WO2021114306A1 WO 2021114306 A1 WO2021114306 A1 WO 2021114306A1 CN 2019125416 W CN2019125416 W CN 2019125416W WO 2021114306 A1 WO2021114306 A1 WO 2021114306A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
projection
estimated
iteration
gpu
Prior art date
Application number
PCT/CN2019/125416
Other languages
English (en)
French (fr)
Inventor
胡战利
曾天翼
杨永峰
高娟
高东芳
梁栋
刘新
郑海荣
Original Assignee
深圳先进技术研究院
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 深圳先进技术研究院 filed Critical 深圳先进技术研究院
Priority to PCT/CN2019/125416 priority Critical patent/WO2021114306A1/zh
Publication of WO2021114306A1 publication Critical patent/WO2021114306A1/zh

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis

Definitions

  • the present invention relates to the field of GPU imaging, in particular to a method and device for high-precision PET reconstruction based on GPU acceleration.
  • Small animal PET is a miniaturized imaging machine that displays physiological structures in the clinic, enabling researchers to quantify the biochemical activity in a variety of organs. It has been used in preclinical imaging research, biological and medical basic research, and new radiopharmaceutical development. widely used.
  • Iterative reconstruction algorithm is one of the mainstream methods in PET reconstruction. Compared with analytical image reconstruction, it can improve the quality of reconstructed images by performing more accurate physical and statistical modeling of the high-energy photon generation and detection process.
  • GPU is a highly parallelized multi-threaded multi-core processor designed to solve image processing tasks, and its floating-point computing power is significantly higher than that of the CPU. More and more researchers use GPU for PET reconstruction calculation.
  • Zhou Jian et al. use GPU to calculate the forward projection and back projection functions, which proves that it can greatly save the calculation time.
  • the traditional iterative reconstruction method uses the CPU to calculate the system matrix and performs forward and back projection of the estimated image. It is close to the true distribution image in iteration.
  • the reconstruction speed depends on the dimension of the input sinogram and the dimension of the generated image, as well as the number of iterations.
  • the high-precision PET with DOI function can output the depth information of photons in the crystal, and the resolution is higher than that of general small animal PET. Therefore, the dimensionality of the sinogram is very large, and the calculation time of the traditional algorithm takes several days.
  • the embodiments of the present invention provide a high-precision PET reconstruction method and device based on GPU acceleration, so as to at least solve the existing technical problem of long calculation time for PET reconstruction using GPU.
  • a high-precision PET reconstruction method based on GPU acceleration which includes the following steps:
  • S100 Input an initial estimated image, and perform an orthographic projection operation on the estimated image to obtain an estimated projection;
  • the method further includes the steps:
  • step S500 Return the input estimated image of the next iteration to step S100 for iterative processing, and the estimated image obtained in the last iteration is the result of the final PET reconstruction.
  • the method further includes the steps:
  • S001 Preprocess the initial estimated image, perform statistics and merge the initial estimated image, and generate a four-dimensional sinogram that meets the preset requirements.
  • Each element of the four-dimensional sinogram represents a type of response line Line Of Response, four-dimensional They are the in-plane distance r from the center of the LOR, the in-plane angle phi, the number of the detector ring that receives the first photon, ring1, and the number of the detector ring that receives the second photon, ring2.
  • the intensity value x of each voxel in the image is estimated, and a Poisson distribution A with an expected intensity value x is constructed.
  • the mathematical model is:
  • k is the number of occurrences of the internal source of each voxel in the image
  • e is the natural constant
  • y is the actual measured projection value
  • the matrix p is related to the geometric structure of the detector and indicates the possibility of being detected by the detector at a certain position
  • D is the distance from the center in the plane of the LOR;
  • N is the normalization correction matrix
  • i is the number of iterations
  • T is the matrix transposition
  • d is the random correction matrix
  • x i is the estimated image
  • x i+1 is the result of this iteration and also the input of the next iteration.
  • the GPU is used to calculate the p T N, and the normalization correction coefficient is back-projected to obtain the sensitivity map of the sensitivity image.
  • the image matrix is bound by the three-dimensional texture memory; in the back projection process, the atomic operation is used to ensure the correctness of the calculation: when the atomic operation is executed, the address is locked, and other parallel threads are not allowed Reading and writing operations to the element corresponding to this address can only wait for the end of the current calculation before starting the next access.
  • the calculation process of the method steps are all processed in parallel in the GPU, and the kernel function is calculated in parallel in each thread at the same time; when the LOR is calculated online through the voxel depth, the calculation results are repeated more according to the symmetry of the system geometry. Times used.
  • the estimated image is first convolved before the projection and after the corrected image is obtained by back-projection. Specifically, the voxel value of each estimated image is compared with the equivalent system function near the voxel value. After multiplying and adding together, the new voxel position image value is obtained.
  • the equivalent system response function H(x,y,z) is a three-dimensional Gauss-like function:
  • r, t, and a represent the radial, tangential and axial directions of the coordinate axis, respectively.
  • the propagation parameter sigma is estimated in the fitting. After obtaining the values of the three parameters ⁇ r , ⁇ t , and ⁇ a, these three The parameters represent the model of the system response function, ⁇ r , ⁇ t , and ⁇ a represent the propagation parameters in the radial, tangential and axial directions, respectively.
  • a high-precision PET reconstruction device based on GPU acceleration including:
  • the orthographic projection unit is used to input the initial estimated image, and perform orthographic operation on the estimated image to obtain an estimated projection;
  • the comparison unit is used to compare the estimated projection obtained with the actual measured projection to obtain a comparison result
  • the back-projection unit is used to perform a back-projection operation on the comparison result to obtain a corrected image
  • the iterative unit is used to multiply the corrected image and the corresponding voxel of the original image, and then divide the sensitivity image to complete one iteration, and the resulting image is the input estimated image of the next iteration.
  • the initial estimated image is input, and the estimated image is subjected to orthographic projection operation to obtain the estimated projection; the obtained estimated projection is compared with the actual measured projection to obtain the comparison Result: Perform back-projection operation on the comparison result to obtain a corrected image; multiply the corrected image with the corresponding voxel of the original image, and then divide the sensitivity image to complete one iteration, and the resulting image is the input estimated image for the next iteration.
  • the iterative reconstruction algorithm of PET is improved by GPU acceleration strategy, and fast calculation is performed to obtain high-resolution, high-precision three-dimensional PET images.
  • Fig. 1 is a flowchart of a high-precision PET reconstruction method based on GPU acceleration according to the present invention
  • FIG. 3 is a graph showing the axial translation characteristics of LOR in the high-precision PET reconstruction method based on GPU acceleration of the present invention
  • FIG. 5 is an overall flowchart of the high-precision PET reconstruction method based on GPU acceleration according to the present invention
  • Fig. 6 is a block diagram of a high-precision PET reconstruction device based on GPU acceleration according to the present invention.
  • Fig. 7 is a preferred module diagram of a high-precision PET reconstruction device based on GPU acceleration of the present invention.
  • a high-precision PET reconstruction method based on GPU acceleration is provided.
  • the method includes the following steps:
  • S100 Input an initial estimated image, and perform an orthographic projection operation on the estimated image to obtain an estimated projection;
  • an initial estimated image is input, and the estimated image is subjected to an orthographic projection operation to obtain an estimated projection; the obtained estimated projection is compared with the actual measured projection to obtain a comparison result; Perform a back-projection operation on the comparison result to obtain a corrected image; multiply the corrected image with the corresponding voxel of the original image, and then divide the sensitivity image to complete one iteration, and the resulting image is the input estimated image for the next iteration.
  • the iterative reconstruction algorithm of PET is improved by GPU acceleration strategy, and fast calculation is performed to obtain high-resolution, high-precision three-dimensional PET images.
  • the method further includes the steps:
  • step S500 Return the input estimated image of the next iteration to step S100 for iterative processing, and the estimated image obtained in the last iteration is the result of the final PET reconstruction.
  • the method further includes the steps:
  • S001 Perform preprocessing on the initial estimated image, perform statistics and merge on the initial estimated image, and generate a four-dimensional chord diagram that meets the preset requirements.
  • the high-precision PET reconstruction method based on GPU acceleration proposed in the present invention is divided into preprocessing and iterative reconstruction.
  • the preprocessing process is to count and merge all events obtained by the detector to generate a four-dimensional sinogram (Sinogram) that meets the requirements.
  • Each element of the chord diagram represents a type of response line (Line Of Response).
  • the four dimensions are the distance from the center of the LOR in the plane (r), the angle in the plane (phi), and the number of the detector ring that receives the first photon (ring1). ), the number of the detector ring that receives the second photon (ring2).
  • the principle of iterative reconstruction is that the PET data y conforms to the variables of the Poisson distribution and is composed of all the LORs collected by the system. Each LOR represents the sum of the signals generated by its source.
  • the present invention needs to estimate the endogenous intensity x of each voxel in the image, and construct a Poisson distribution A whose expectation is x, and its mathematical model is:
  • k is the number of occurrences of the internal source of each voxel in the image
  • e is the natural constant
  • y is the actual measured projection value
  • the matrix p is related to the geometric structure of the detector and indicates the possibility of being detected by the detector at a certain position
  • D is the distance from the center in the plane of the LOR.
  • the present invention needs to find a distribution of x’, x’ is a variable, and its dimension is the same as x, so that P(y
  • the maximum expectation method can be used to find the extreme points of formula (4). It is usually used to estimate the parameters of a probability model containing hidden variables or missing data. It is an optimization algorithm for maximum likelihood estimation through iteration. The convergence of the algorithm can be Ensure that the iteration at least approximates the local maximum. Without considering attenuation and scattering:
  • N is the normalization correction matrix
  • i is the number of iterations
  • T is the matrix transposition
  • d is the random correction matrix
  • x i is the estimated image
  • x i+1 is the result of this iteration and also the input of the next iteration.
  • the GPU When operating in the computer, input the initial estimated image. Generally, the voxel value of each image is 1. Before the iteration, the GPU is used to calculate the p T N, that is, the normalization correction coefficient is back-projected to obtain the sensitivity image (sensitivity image). map). According to the formula, perform the forward projection operation on the estimated image x i . The estimated projection obtained is compared with the actual measured projection, and the result obtained is subjected to the back projection operation. Then the obtained image is multiplied by the corresponding voxel of the original image, and then divided by the sensitivity map , After completing one iteration, the image obtained is the input estimated image of the next iteration. All processes are calculated in the GPU.
  • the principle of GPU parallelism is to calculate the kernel function in parallel in each thread at the same time.
  • the estimated image obtained in the last iteration is the final reconstruction result.
  • the fast floating-point calculations of the GPU can be fully utilized.
  • the present invention uses a three-dimensional texture memory to bind the image matrix.
  • Texture memory is a special form of CUDA global memory. Its read and write operations are performed through a special texture cache, which allows multiple threads to access it at the same time, and the access speed is very fast.
  • the present invention uses atomic operations to ensure the correctness of the calculation.
  • the atomic operation is executed, the address is locked, and other parallel threads are not allowed to read and write the element corresponding to the address, and can only wait for the end of the current calculation to start the next access, ensuring safety.
  • the GPU’s simultaneous parallel stream processors are about 1000 to 3000, which means that in fact the GPU can only calculate so many threads at the same time. Therefore, for high-resolution PET systems containing DOI information, it is impossible to calculate all the allocated threads in parallel at the same time. In fact, calculations are performed batch by batch, and this process is serial.
  • the method of the present invention utilizes the symmetry of the geometric structure of the system to reduce the total number of LORs that need to be calculated, thereby reducing the number of parallel threads and speeding up the calculation time.
  • the consideration of symmetry is that in the process of projection and back-projection, that is, when calculating the depth of LOR through voxels online, the calculation results can be used repeatedly according to the geometric structure characteristics of the system.
  • the PET reconstruction method proposed in the present invention not only corrects the uniformity and random events, but also adds the correction of the physical effects of the system's positron free path and photon nonlinearity.
  • the method is to measure the point source in the system to reconstruct the image in the system. Function, using these physical effects to the point source image to be equivalent.
  • the equivalent system response function H(x,y,z) is a three-dimensional Gauss-like function:
  • r, t, and ⁇ represent the radial, tangential, and axial directions of the coordinate axis.
  • the propagation parameter sigma is estimated in the fitting. After obtaining the values of the three parameters ⁇ r , ⁇ t , and ⁇ a, these three The parameters can represent the model of the system response function. ⁇ r , ⁇ t , and ⁇ a represent the propagation parameters in the radial, tangential and axial directions, respectively. In the iterative reconstruction process, the estimated image must be convolved before projection.
  • the voxel value of each estimated image is multiplied by the equivalent system response function near the voxel value and then added to obtain a new The image value of the voxel position.
  • the present invention collects multiple point sources in different positions of the system, and selects representative images among them for fitting to obtain the system response function.
  • the overall flow chart of the reconstruction method is shown in Figure 5.
  • a high-precision PET reconstruction device based on GPU acceleration See FIG. 6, including:
  • the orthographic projection unit 201 is used to input an initial estimated image, and perform an orthographic projection operation on the estimated image to obtain an estimated projection;
  • the comparing unit 202 is configured to compare the obtained estimated projection with the actual measured projection to obtain a comparison result
  • the back-projection unit 203 is configured to perform a back-projection operation on the comparison result to obtain a corrected image
  • the iteration unit 204 is configured to multiply the corrected image and the corresponding voxel of the original image, and then divide the sensitivity image to complete one iteration, and the obtained image is the input estimated image of the next iteration.
  • an initial estimated image is input, and the estimated image is subjected to an orthographic projection operation to obtain an estimated projection; the obtained estimated projection is compared with the actual measured projection to obtain a comparison result; Perform a back-projection operation on the comparison result to obtain a corrected image; multiply the corrected image with the corresponding voxel of the original image, and then divide the sensitivity image to complete one iteration, and the resulting image is the input estimated image for the next iteration.
  • the iterative reconstruction algorithm of PET is improved by GPU acceleration strategy, and fast calculation is performed to obtain high-resolution, high-precision three-dimensional PET images.
  • the device further includes:
  • the returning iterative unit 205 is used to return the input estimated image of the next iteration to the orthographic projection unit 201 for iterative processing, and the estimated image obtained in the last iteration is the result of the final PET reconstruction.
  • the device further includes:
  • the preprocessing unit 200 is configured to preprocess the initial estimated image, perform statistics and merge the initial estimated image, and generate a four-dimensional chord diagram that meets the preset requirements.
  • the GPU-accelerated high-precision PET reconstruction device proposed in the present invention is divided into preprocessing and iterative reconstruction.
  • the preprocessing unit 200 the preprocessing process is to count and merge all events obtained by the detector to generate a four-dimensional chord diagram that meets the requirements ( Sinogram), each element of the four-dimensional sinogram represents a type of line of response (Line Of Response).
  • the four dimensions are the distance from the center in the LOR plane (r), the in-plane angle (phi), and the detection of the first photon received.
  • the principle of iterative reconstruction is that the PET data y conforms to the variables of the Poisson distribution and is composed of all the LORs collected by the system. Each LOR represents the sum of the signals generated by its source.
  • the present invention needs to estimate the endogenous intensity x of each voxel in the image, and construct a Poisson distribution A whose expectation is x, and its mathematical model is:
  • k is the number of occurrences of the internal source of each voxel in the image
  • e is the natural constant
  • y is the actual measured projection value
  • the matrix p is related to the geometric structure of the detector and indicates the possibility of being detected by the detector at a certain position
  • D is the distance from the center in the plane of the LOR.
  • the present invention needs to find a distribution of x’ to make P(y
  • x’ is a variable whose dimension is the same as x
  • the likelihood function is:
  • the maximum expectation method can be used to find the extreme points of formula (4). It is usually used to estimate the parameters of a probability model containing hidden variables or missing data. It is an optimization algorithm for maximum likelihood estimation through iteration. The convergence of the algorithm can be Ensure that the iteration at least approximates the local maximum. Without considering attenuation and scattering:
  • N is the normalization correction matrix
  • i is the number of iterations
  • T is the matrix transposition
  • d is the random correction matrix
  • x i is the estimated image
  • x i+1 is the result of this iteration and also the input of the next iteration.
  • the GPU When operating in the computer, input the initial estimated image. Generally, the voxel value of each image is 1. Before the iteration, the GPU is used to calculate the p T N, that is, the normalization correction coefficient is back-projected to obtain the sensitivity image (sensitivity image). map). According to the formula, perform the forward projection operation on the estimated image x i . The estimated projection obtained is compared with the actual measured projection, and the result obtained is subjected to the back projection operation. Then the obtained image is multiplied by the corresponding voxel of the original image, and then divided by the sensitivity map , After completing one iteration, the image obtained is the input estimated image of the next iteration. All processes are calculated in the GPU.
  • the principle of GPU parallelism is to calculate the kernel function in parallel in each thread at the same time.
  • the estimated image obtained in the last iteration is the final reconstruction result.
  • the fast floating-point calculations of the GPU can be fully utilized.
  • the present invention uses a three-dimensional texture memory to bind the image matrix.
  • Texture memory is a special form of CUDA global memory. Its read and write operations are performed through a special texture cache, which allows multiple threads to access it at the same time, and the access speed is very fast.
  • the present invention uses atomic operations to ensure the correctness of the calculation.
  • the atomic operation is executed, the address is locked, and other parallel threads are not allowed to read and write the element corresponding to the address, and can only wait for the end of the current calculation to start the next access, ensuring safety.
  • the GPU’s simultaneous parallel stream processors are about 1000 to 3000, which means that in fact the GPU can only calculate so many threads at the same time. Therefore, for high-resolution PET systems containing DOI information, it is impossible to calculate all the allocated threads in parallel at the same time. In fact, calculations are performed batch by batch, and this process is serial.
  • the method of the present invention utilizes the symmetry of the geometric structure of the system to reduce the total number of LORs that need to be calculated, thereby reducing the number of parallel threads and speeding up the calculation time.
  • the consideration of symmetry is that in the process of projection and back-projection, that is, when calculating the depth of LOR through voxels online, the calculation results can be used repeatedly according to the geometric structure characteristics of the system.
  • the PET reconstruction method proposed in the present invention not only corrects the uniformity and random events, but also adds the correction of the physical effects of the system's positron free path and photon nonlinearity.
  • the method is to measure the point source in the system to reconstruct the image in the system. Function, using these physical effects to the point source image to be equivalent.
  • the equivalent system response function H(x,y,z) is a three-dimensional Gauss-like function:
  • r, t, and a represent the radial, tangential and axial directions of the coordinate axis, respectively.
  • the propagation parameter sigma is estimated in the fitting. After obtaining the values of the three parameters ⁇ r , ⁇ t , and ⁇ a, these three The parameters can represent the model of the system response function. ⁇ r , ⁇ t , and ⁇ a represent the propagation parameters in the radial, tangential and axial directions, respectively. In the iterative reconstruction process, the estimated image must be convolved before projection.
  • the voxel value of each estimated image is multiplied by the equivalent system response function near the voxel value and then added to obtain a new The image value of the voxel position.
  • the present invention collects multiple point sources in different positions of the system, and selects representative images among them for fitting to obtain the system response function.
  • the overall flow chart of the reconstruction method is shown in Figure 5.
  • the present invention uses the iterative reconstruction algorithm in GPU calculations, and considers the geometric symmetry of the system and the CUDA C texture memory acceleration method, and proposes a fast and high-precision high-precision PET reconstruction method based on GPU acceleration.
  • Parallel acceleration and the geometric characteristics of the detector parallelize the reconstruction of ultra-large-scale data to perform high-precision and rapid reconstruction, that is, through the GPU acceleration strategy to improve the iterative reconstruction algorithm of PET and perform rapid calculations to obtain high-resolution, high-precision 3D PET image.
  • the reconstructed image can be obtained by simply inputting the data collected by the PET system. Improve reconstruction speed and image resolution.
  • the protection points of the method and device are: a GPU-based method for accelerating high-definition iterative reconstruction using geometric symmetry; and a system response function calculation method based on actual collected data.
  • the disclosed technical content can be implemented in other ways.
  • the system embodiment described above is only illustrative.
  • the division of units may be a logical function division, and there may be other divisions in actual implementation.
  • multiple units or components may be combined or integrated into Another system, or some features can be ignored, or not implemented.
  • the displayed or discussed mutual coupling or direct coupling or communication connection may be indirect coupling or communication connection through some interfaces, units or modules, and may be in electrical or other forms.
  • the units described as separate components may or may not be physically separate, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple units. Some or all of the units may be selected according to actual needs to achieve the objectives of the solutions of the embodiments.
  • the functional units in the various embodiments of the present invention may be integrated into one processing unit, or each unit may exist alone physically, or two or more units may be integrated into one unit.
  • the above-mentioned integrated unit can be implemented in the form of hardware or software functional unit.
  • the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer readable storage medium.
  • the technical solution of the present invention essentially or the part that contributes to the existing technology or all or part of the technical solution can be embodied in the form of a software product, and the computer software product is stored in a storage medium.
  • a computer device which can be a personal computer, a server, or a network device, etc.
  • the aforementioned storage media include: U disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), mobile hard disk, magnetic disk or optical disk and other media that can store program codes. .

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Nuclear Medicine (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

一种基于GPU加速的高精度PET重建方法及装置,涉及GPU成像领域。该方法包括:输入初始的估计图像,对估计图像进行正投影操作,获得估计投影(S100);将获得的估计投影与实际测量投影作比较,获得比较结果(S200);对比较结果进行反投影操作,获得校正图像(S300);将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像(S400)。该方法通过GPU加速策略将PET迭代重建算法改进后进行快速计算,获得高分辨率、高精度的三维PET图像。

Description

一种基于GPU加速的高精度PET重建方法及装置 技术领域
本发明涉及GPU成像领域,具体而言,涉及一种基于GPU加速的高精度PET重建方法及装置。
背景技术
小动物PET是小型化在临床上显示生理结构的成像机器,使研究人员能够量化多种器官中的生化活性,目前已经在临床前成像研究、生物学和医学基础研究、新的放射性药物开发上广泛应用。迭代重建算法是PET重建中的主流方法之一,相较于解析图像重建,它可以通过对高能光子产生和探测过程进行更精确的物理和统计建模,提升重建图像质量。GPU是为了解决图像处理任务而被设计出来的高度并行化的多线程多核心处理器,浮点计算能力显著高于CPU。越来越多的研究者使用GPU进行PET重建计算,Zhou Jian等用GPU计算正投影和反投影函数,证明能极大节省计算时间。
传统的迭代重建方法使用CPU计算系统矩阵,对估计图像进行正反投影,在迭代中接近真实分布图像,其重建速度依赖于输入sinogram的维度大小和生成的图像维度大小,以及迭代的次数。高精度带有DOI功能的PET能输出光子在晶体中的深度信息,分辨率比一般小动物PET高,因此sinogram的维度非常大,传统算法计算时间需要数天。
发明内容
本发明实施例提供了一种基于GPU加速的高精度PET重建方法及装置,以至少解决现有使用GPU进行PET重建计算时间长的技术问题。
根据本发明的一实施例,提供了一种基于GPU加速的高精度PET重建方法,包括以下步骤:
S100:输入初始的估计图像,对估计图像进行正投影操作,获得估计投影;
S200:将获得的估计投影与实际测量投影作比较,获得比较结果;
S300:对比较结果进行反投影操作,获得校正图像;
S400:将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像。
进一步地,方法还包括步骤:
S500:将下一次迭代的输入估计图像返回步骤S100进行迭代处理,最后一次迭代获得的估计图像即是最终PET重建的结果。
进一步地,方法还包括步骤:
S001:对初始的估计图像进行预处理,将初始的估计图像进行统计和合并,生成符合预设要求的四维弦图Sinogram,该四维弦图的每一个元素代表一类响应线Line Of Response,四维分别是LOR的平面内距中心距离r、平面内角度phi、接收第一个光子的探测器环编号ring1、接收第二个光子的探测器环编号ring2。
进一步地,在迭代过程中估计出图像中每个体素内源的强度值x,构建一个期望是强度值x的泊松分布A,其数学模型为:
Figure PCTCN2019125416-appb-000001
x和y的期望的关系是:
Figure PCTCN2019125416-appb-000002
其中,k是图像中每个体素内源的出现次数,e是自然常数,y是实际测量的投影值,矩阵p是与探测器几何结构相关的表示在某一位置被探测器探测到的可能性,d是LOR的平面内距中心距离;
在PET重建中找到一个x’的分布,使得P(y|x’)最大,x’是一个变量,其 维度与x相同,似然函数为:
Figure PCTCN2019125416-appb-000003
利用对数似然函数极值点不变的特性,将公式(2)和公式(3)结合成l(y|x)=logL(y|x),其中:
Figure PCTCN2019125416-appb-000004
使用最大期望方法求公式(4)的极值点,在不考虑衰减和散射的情况下:
Figure PCTCN2019125416-appb-000005
其中N是均一化校正矩阵,i是迭代次数,T表示矩阵转置,d是随机校正矩阵,x i为估计图像,x i+1是本次迭代的结果,也是下一次迭代的输入。
进一步地,进行迭代前先利用GPU计算p TN,把均一化校正系数做反投影,得到敏感度图像sensitivity map。
进一步地,在正投影过程中,使用三维纹理内存对图像矩阵进行绑定;在反投影过程中,用原子操作来确保计算的正确性:原子操作执行时,锁定该地址,其他并行线程不允许对该地址对应的元素进行读取和写入操作,只能等待当前计算结束才能开始下一次访问。
进一步地,方法步骤的计算过程均在GPU中并行处理,将核函数在每个线程中同时并行计算;在线计算LOR穿过体素深度时,根据系统几何结构的对称性,将计算结果重复多次使用。
进一步地,迭代重建过程中,在进行投影之前及反投影得到校正图像后先对估计图像进行卷积操作,具体是将每个估计图像体素值与该体素值附近的等效系统函数相乘后相加,得到新的该体素位置图像值。
进一步地,测量点源在系统中重建图像上的系统函数,用这些物理效应对点源图像的影响来等效;等效系统响应函数H(x,y,z)是一个类三维高斯函数:
Figure PCTCN2019125416-appb-000006
其中,r,t,a分别代表坐标轴的径向、切向和轴向,传播参数sigma在拟合中被估计,得到三个参数σ r、σ t、σ a的值后,这三个参数代表系统响应函数的 模型,σ r、σ t、σ a分别表示径向、切向和轴向的传播参数。
根据本发明的另一实施例,提供了一种基于GPU加速的高精度PET重建装置,包括:
正投影单元,用于输入初始的估计图像,对估计图像进行正投影操作,获得估计投影;
比较单元,用于将获得的估计投影与实际测量投影作比较,获得比较结果;
反投影单元,用于对比较结果进行反投影操作,获得校正图像;
迭代单元,用于将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像。
本发明实施例中的基于GPU加速的高精度PET重建方法及装置,输入初始的估计图像,对估计图像进行正投影操作,获得估计投影;将获得的估计投影与实际测量投影作比较,获得比较结果;对比较结果进行反投影操作,获得校正图像;将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像。通过GPU加速策略将PET迭代重建算法改进后进行快速计算,获得高分辨率、高精度的三维PET图像。
附图说明
此处所说明的附图用来提供对本发明的进一步理解,构成本申请的一部分,本发明的示意性实施例及其说明用于解释本发明,并不构成对本发明的不当限定。在附图中:
图1为本发明基于GPU加速的高精度PET重建方法的流程图;
图2为本发明基于GPU加速的高精度PET重建方法的优选流程图;
图3为本发明基于GPU加速的高精度PET重建方法中LOR的轴向平移特性图;
图4为本发明基于GPU加速的高精度PET重建方法中LOR的轴向对称特性图;
图5为本发明基于GPU加速的高精度PET重建方法的整体流程图;
图6为本发明基于GPU加速的高精度PET重建装置的模块图;
图7为本发明基于GPU加速的高精度PET重建装置的优选模块图。
具体实施方式
为了使本技术领域的人员更好地理解本发明方案,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分的实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都应当属于本发明保护的范围。
需要说明的是,本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象,而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换,以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外,术语“包括”和“具有”以及他们的任何变形,意图在于覆盖不排他的包含,例如,包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元,而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
实施例1
根据本发明一实施例,提供了一种基于GPU加速的高精度PET重建方法,参见图1,包括以下步骤:
S100:输入初始的估计图像,对估计图像进行正投影操作,获得估计投影;
S200:将获得的估计投影与实际测量投影作比较,获得比较结果;
S300:对比较结果进行反投影操作,获得校正图像;
S400:将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像。
本发明实施例中的基于GPU加速的高精度PET重建方法,输入初始的估计图像,对估计图像进行正投影操作,获得估计投影;将获得的估计投影与实际测量投影作比较,获得比较结果;对比较结果进行反投影操作,获得校正图像;将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像。通过GPU加速策略将PET迭代重建算法改进后进行快速计算,获得高分辨率、高精度的三维PET图像。
具体地,参见图2,方法还包括步骤:
S500:将下一次迭代的输入估计图像返回步骤S100进行迭代处理,最后一次迭代获得的估计图像即是最终PET重建的结果。
具体地,参见图2,方法还包括步骤:
S001:对初始的估计图像进行预处理,将初始的估计图像进行统计和合并,生成符合预设要求的四维弦图。
下面以具体的实施例,对本发明基于GPU加速的高精度PET重建方法进行详细描述。
本发明提出的基于GPU加速的高精度PET重建方法分为预处理和迭代重建,预处理过程为把探测器得到的所有事件进行统计和合并,生成符合要求的四维弦图(Sinogram),该四维弦图的每一个元素代表一类响应线(Line Of Response),四维分别是LOR的平面内距中心距离(r)、平面内角度(phi)、接收第一个光子的探测器环编号(ring1)、接收第二个光子的探测器环编号(ring2)。
迭代重建的原理是PET数据y符合泊松分布的变量,是由系统收集到的所有LOR组成,每一条LOR代表着经过它的源产生的信号的总和。本发明需要估计出图像中每个体素内源的强度x,构建一个期望是x的泊松分布A,其数学模型为:
Figure PCTCN2019125416-appb-000007
x和y的期望的关系是:
Figure PCTCN2019125416-appb-000008
其中,k是图像中每个体素内源的出现次数,e是自然常数,y是实际测量的投影值,矩阵p是与探测器几何结构相关的表示在某一位置被探测器探测到的可能性,d是LOR的平面内距中心距离。在PET重建中本发明需要找到一个x’的分布,x’是一个变量,其维度与x相同,使得P(y|x’)最大,似然函数为:
Figure PCTCN2019125416-appb-000009
利用对数似然函数极值点不变的特性,可以将公式(2)和公式(3)结合成l(y|x)=logL(y|x),其中:
Figure PCTCN2019125416-appb-000010
使用最大期望方法可以求公式(4)的极值点,通常用于对包含隐变量或缺失数据的概率模型进行参数估计,是通过迭代进行极大似然估计的优化算法,算法的收敛性可以确保迭代至少逼近局部极大值。在不考虑衰减和散射的情况下:
Figure PCTCN2019125416-appb-000011
其中N是均一化校正矩阵,i是迭代次数,T表示矩阵转置,d是随机校正矩阵,x i为估计图像,x i+1是本次迭代的结果,也是下一次迭代的输入。
计算机中操作的时候,输入初始估计图像,一般每个图像体素值都为1,在进行迭代前先利用GPU计算p TN,即把均一化校正系数做反投影,得到敏感度图像(sensitivity map)。按公式先对估计图像x i进行正投影操作,得到的估计投影与实际测量投影作比,得到的结果进行反投影操作,然后把得到的图像与原始图像对应体素相乘,然后除sensitivity map,完成一次迭代,得到的图像就是下一次迭代的输入估计图像,所有过程均在GPU中计算,GPU并行的原理是将核函数在每个线程中同时并行计算。最后一次迭代获得的估计图像就是最终重建的结果。对于重建中所用的函数,由于不需要线程间进行过多的通信,可以充分发挥GPU浮点计算速度快的特点。
在正投影过程中,如果每种LOR分配一个线程将有大于3*109个线程,如此多个线程进行并行,必然会有同时访问一个图像体素内存情况发生,造成数据竞争的问题,降低运行速度,对此本发明使用三维纹理内存对图像矩阵进行绑定。纹理内存是CUDA全局内存的特殊形式,对其读写操作通过专门的纹理缓存进行,可以让多个线程同时进行访问,且访问速度很快。
在反投影过程中,由于多个线程会同时修改同一个图像体素值,因此本发明使用原子操作来确保计算的正确性。原子操作执行时,锁定该地址,其他并行线程不允许对该地址对应的元素进行读取和写入操作,只能等待当前计算结束才能开始下一次访问,确保了安全。
GPU同时并行的流处理器为1000到3000左右,意味着实际上GPU只能同时计算这么多的线程,因此对于含有DOI信息的高分辨率PET系统,不可能同时并行运算所有分配的线程,线程实际上是一批批的进行计算,而这个过程是串行的。本发明的方法是利用系统的几何结构的对称性,可以减小需要计算的总LOR数量,从而减少并行线程数量,加快计算时间。对称性的考虑是在投影和反投影过程中,即在线计算LOR穿过体素深度时,根据系统的几何结构特性,可以将计算结果重复多次使用。有三类对称可以使用,第一类对称称为LOR的轴向平移特性,是指一条LOR可以在轴向平移,在相应的体素中深度不变,如图3所示,将LOR轴向平移,经过对应的体素深度不变,d1=d2=d3,因此可以将这些LOR归为一类进行计算;第二类对称称为LOR的轴向对称特性,是指LOR在轴向上中心对称,因此两条对称的LOR经过对称的像素深度完全相同,如图4所示,LOR轴向中心对称的另一条LOR,经过对应的体素深度不变,d1=d2,因此可以将这2条LOR归为一类进行计算;第三类对称称为LOR的径向对称特性,因为LOR在径向上根据r和phi中心对称,因此可以将4条LOR归为一类,该对称可以类比第二类对称。
本发明提出的PET重建方法除了对均一化和随机事件进行校正以外,还加入了对系统正电子自由程和光子非线性的物理效应的校正,方法是测量点源在系统中重建图像上的系统函数,用这些物理效应对点源图像的影响来等效。等效系统响应函数H(x,y,z)是一个类三维高斯函数:
Figure PCTCN2019125416-appb-000012
其中,r,t,σ分别代表坐标轴的径向、切向和轴向,传播参数sigma在拟合中被估计,得到三个参数σ r、σ t、σ a的值后,这三个参数就能代表系统响应函数的模型,σ r、σ t、σ a分别表示径向、切向和轴向的传播参数。迭代重建过程中,在进行投影之前先要对估计图像进行卷积操作,具体是将每个估计图像体素值与该体素值附近的等效系统响应函数相乘后相加,得到新的该体素位置图像值。在反投影得到校正图像后,也要对其类似操作。本发明采集了系统多个不同位置的点源,选择了其中有代表性的图像进行拟合获得了系统响应函数。重建方法整体流程图如图5所示。
实施例2
根据本发明的另一实施例,提供了一种基于GPU加速的高精度PET重建装置,参见图6,包括:
正投影单元201,用于输入初始的估计图像,对估计图像进行正投影操作,获得估计投影;
比较单元202,用于将获得的估计投影与实际测量投影作比较,获得比较结果;
反投影单元203,用于对比较结果进行反投影操作,获得校正图像;
迭代单元204,用于将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像。
本发明实施例中的基于GPU加速的高精度PET重建装置,输入初始的估计图像,对估计图像进行正投影操作,获得估计投影;将获得的估计投影与实际测量投影作比较,获得比较结果;对比较结果进行反投影操作,获得校正图像;将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像。通过GPU加速策略将PET迭代重建算法改进后进行快速计算,获得高分辨率、高精度的三维PET图像。
具体地,参见图6,装置还包括:
返回迭代单元205,用于将下一次迭代的输入估计图像返回正投影单元201进行迭代处理,最后一次迭代获得的估计图像即是最终PET重建的结果。
具体地,参见图6,装置还包括:
预处理单元200,用于对初始的估计图像进行预处理,将初始的估计图像进行统计和合并,生成符合预设要求的四维弦图。
下面以具体的实施例,对本发明基于GPU加速的高精度PET重建装置进行详细描述。
本发明提出的基于GPU加速的高精度PET重建装置分为预处理和迭代重建,预处理单元200:预处理过程为把探测器得到的所有事件进行统计和合并,生成符合要求的四维弦图(Sinogram),该四维弦图的每一个元素代表一类响应线(Line Of Response),四维分别是LOR的平面内距中心距离(r)、平面内角度(phi)、接收第一个光子的探测器环编号(ring1)、接收第二个光子的探测器环编号(ring2)。
迭代重建的原理是PET数据y符合泊松分布的变量,是由系统收集到的所有LOR组成,每一条LOR代表着经过它的源产生的信号的总和。本发明需要估计出图像中每个体素内源的强度x,构建一个期望是x的泊松分布A,其数学模型为:
Figure PCTCN2019125416-appb-000013
x和y的期望的关系是:
Figure PCTCN2019125416-appb-000014
其中,k是图像中每个体素内源的出现次数,e是自然常数,y是实际测量的投影值,矩阵p是与探测器几何结构相关的表示在某一位置被探测器探测到的可能性,d是LOR的平面内距中心距离。在PET重建中本发明需要找到一个x’的分布,使得P(y|x’)最大,x’是一个变量,其维度与x相同,似然函数为:
Figure PCTCN2019125416-appb-000015
利用对数似然函数极值点不变的特性,可以将公式(2)和公式(3)结合成 l(y|x)=logL(y|x),其中:
Figure PCTCN2019125416-appb-000016
使用最大期望方法可以求公式(4)的极值点,通常用于对包含隐变量或缺失数据的概率模型进行参数估计,是通过迭代进行极大似然估计的优化算法,算法的收敛性可以确保迭代至少逼近局部极大值。在不考虑衰减和散射的情况下:
Figure PCTCN2019125416-appb-000017
其中N是均一化校正矩阵,i是迭代次数,T表示矩阵转置,d是随机校正矩阵,x i为估计图像,x i+1是本次迭代的结果,也是下一次迭代的输入。
计算机中操作的时候,输入初始估计图像,一般每个图像体素值都为1,在进行迭代前先利用GPU计算p TN,即把均一化校正系数做反投影,得到敏感度图像(sensitivity map)。按公式先对估计图像x i进行正投影操作,得到的估计投影与实际测量投影作比,得到的结果进行反投影操作,然后把得到的图像与原始图像对应体素相乘,然后除sensitivity map,完成一次迭代,得到的图像就是下一次迭代的输入估计图像,所有过程均在GPU中计算,GPU并行的原理是将核函数在每个线程中同时并行计算。最后一次迭代获得的估计图像就是最终重建的结果。对于重建中所用的函数,由于不需要线程间进行过多的通信,可以充分发挥GPU浮点计算速度快的特点。
在正投影过程中,如果每种LOR分配一个线程将有大于3*109个线程,如此多个线程进行并行,必然会有同时访问一个图像体素内存情况发生,造成数据竞争的问题,降低运行速度,对此本发明使用三维纹理内存对图像矩阵进行绑定。纹理内存是CUDA全局内存的特殊形式,对其读写操作通过专门的纹理缓存进行,可以让多个线程同时进行访问,且访问速度很快。
在反投影过程中,由于多个线程会同时修改同一个图像体素值,因此本发明使用原子操作来确保计算的正确性。原子操作执行时,锁定该地址,其他并行线程不允许对该地址对应的元素进行读取和写入操作,只能等待当前计算结束才能开始下一次访问,确保了安全。
GPU同时并行的流处理器为1000到3000左右,意味着实际上GPU只能同时计算这么多的线程,因此对于含有DOI信息的高分辨率PET系统,不可能同时并行运算所有分配的线程,线程实际上是一批批的进行计算,而这个过程是串行的。本发明的方法是利用系统的几何结构的对称性,可以减小需要计算的总LOR数量,从而减少并行线程数量,加快计算时间。对称性的考虑是在投影和反投影过程中,即在线计算LOR穿过体素深度时,根据系统的几何结构特性,可以将计算结果重复多次使用。有三类对称可以使用,第一类对称称为LOR的轴向平移特性,是指一条LOR可以在轴向平移,在相应的体素中深度不变,如图3所示,将LOR轴向平移,经过对应的体素深度不变,d1=d2=d3,因此可以将这些LOR归为一类进行计算;第二类对称称为LOR的轴向对称特性,是指LOR在轴向上中心对称,因此两条对称的LOR经过对称的像素深度完全相同,如图4所示,LOR轴向中心对称的另一条LOR,经过对应的体素深度不变,d1=d2,因此可以将这2条LOR归为一类进行计算;第三类对称称为LOR的径向对称特性,因为LOR在径向上根据r和phi中心对称,因此可以将4条LOR归为一类,该对称可以类比第二类对称。
本发明提出的PET重建方法除了对均一化和随机事件进行校正以外,还加入了对系统正电子自由程和光子非线性的物理效应的校正,方法是测量点源在系统中重建图像上的系统函数,用这些物理效应对点源图像的影响来等效。等效系统响应函数H(x,y,z)是一个类三维高斯函数:
Figure PCTCN2019125416-appb-000018
其中,r,t,a分别代表坐标轴的径向、切向和轴向,传播参数sigma在拟合中被估计,得到三个参数σ r、σ t、σ a的值后,这三个参数就能代表系统响应函数的模型,σ r、σ t、σ a分别表示径向、切向和轴向的传播参数。迭代重建过程中,在进行投影之前先要对估计图像进行卷积操作,具体是将每个估计图像体素值与该体素值附近的等效系统响应函数相乘后相加,得到新的该体素位置图像值。在反投影得到校正图像后,也要对其类似操作。本发明采集了系统多个不同位置的点源,选择了其中有代表性的图像进行拟合获得了系统响应函数。重建方法整体流程图如图5所示。
本发明将迭代重建算法使用在GPU计算中,并考虑了系统的几何对称性和CUDA C纹理内存加速方法,提出一种速度快、精度高的基于GPU加速的高精度PET重建方法,通过CUDA C并行加速和探测器几何特性将超大规模数据量的重建工作并行处理,进行高精度快速重建,即通过GPU加速策略将PET迭代重建算法改进后进行快速计算,获得高分辨率、高精度的三维PET图像。
在实际应用中,只需将PET系统采集到的数据输入即可获得重建图像。提高重建速度,提高图像分辨率。该方法及装置的保护点在于:基于GPU的利用几何对称性加速高清迭代重建方法;基于实际采集数据的系统响应函数计算方法。
上述本发明实施例序号仅仅为了描述,不代表实施例的优劣。
在本发明的上述实施例中,对各个实施例的描述都各有侧重,某个实施例中没有详述的部分,可以参见其他实施例的相关描述。
在本申请所提供的几个实施例中,应该理解到,所揭露的技术内容,可通过其它的方式实现。其中,以上所描述的系统实施例仅仅是示意性的,例如单元的划分,可以为一种逻辑功能划分,实际实现时可以有另外的划分方式,例如多个单元或组件可以结合或者可以集成到另一个系统,或一些特征可以忽略,或不执行。另一点,所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口,单元或模块的间接耦合或通信连接,可以是电性或其它的形式。
作为分离部件说明的单元可以是或者也可以不是物理上分开的,作为单元显示的部件可以是或者也可以不是物理单元,即可以位于一个地方,或者也可以分布到多个单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外,在本发明各个实施例中的各功能单元可以集成在一个处理单元中,也可以是各个单元单独物理存在,也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现,也可以采用软件功能单元的 形式实现。
集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时,可以存储在一个计算机可读取存储介质中。基于这样的理解,本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来,该计算机软件产品存储在一个存储介质中,包括若干指令用以使得一台计算机设备(可为个人计算机、服务器或者网络设备等)执行本发明各个实施例方法的全部或部分步骤。而前述的存储介质包括:U盘、只读存储器(ROM,Read-Only Memory)、随机存取存储器(RAM,Random Access Memory)、移动硬盘、磁碟或者光盘等各种可以存储程序代码的介质。
以上所述仅是本发明的优选实施方式,应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。

Claims (10)

  1. 一种基于GPU加速的高精度PET重建方法,其特征在于,包括以下步骤:
    S100:输入初始的估计图像,对估计图像进行正投影操作,获得估计投影;
    S200:将获得的估计投影与实际测量投影作比较,获得比较结果;
    S300:对比较结果进行反投影操作,获得校正图像;
    S400:将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像。
  2. 根据权利要求1所述的基于GPU加速的高精度PET重建方法,其特征在于,所述方法还包括步骤:
    S500:将下一次迭代的输入估计图像返回步骤S100进行迭代处理,最后一次迭代获得的估计图像即是最终PET重建的结果。
  3. 根据权利要求2所述的基于GPU加速的高精度PET重建方法,其特征在于,所述方法还包括步骤:
    S001:对初始的估计图像进行预处理,将初始的估计图像进行统计和合并,生成符合预设要求的四维弦图Sinogram,该四维弦图的每一个元素代表一类响应线Line Of Response,四维分别是LOR的平面内距中心距离r、平面内角度phi、接收第一个光子的探测器环编号ring1、接收第二个光子的探测器环编号ring2。
  4. 根据权利要求3所述的基于GPU加速的高精度PET重建方法,其特征在于,在迭代过程中估计出图像中每个体素内源的强度值x,构建一个期望是强度值x的泊松分布A,其数学模型为:
    Figure PCTCN2019125416-appb-100001
    x和y的期望的关系是:
    Figure PCTCN2019125416-appb-100002
    其中,k是图像中每个体素内源的出现次数,e是自然常数,y是实际测量的投影值,矩阵p是与探测器几何结构相关的表示在某一位置被探测器探测到的 可能性,d是LOR的平面内距中心距离;
    在PET重建中找到一个x’的分布,使得P(y|x’)最大,x’是一个变量,其维度与x相同,似然函数为:
    Figure PCTCN2019125416-appb-100003
    利用对数似然函数极值点不变的特性,将公式(2)和公式(3)结合成l(y|x)=logL(y|x),其中:
    Figure PCTCN2019125416-appb-100004
    使用最大期望方法求公式(4)的极值点,在不考虑衰减和散射的情况下:
    Figure PCTCN2019125416-appb-100005
    其中N是均一化校正矩阵,i是迭代次数,T表示矩阵转置,d是随机校正矩阵,x i为估计图像,x i+1是本次迭代的结果,也是下一次迭代的输入。
  5. 根据权利要求4所述的基于GPU加速的高精度PET重建方法,其特征在于,进行迭代前先利用GPU计算p TN,把均一化校正系数做反投影,得到敏感度图像sensitivity map。
  6. 根据权利要求5所述的基于GPU加速的高精度PET重建方法,其特征在于,在正投影过程中,使用三维纹理内存对图像矩阵进行绑定;在反投影过程中,用原子操作来确保计算的正确性:原子操作执行时,锁定该地址,其他并行线程不允许对该地址对应的元素进行读取和写入操作,只能等待当前计算结束才能开始下一次访问。
  7. 根据权利要求6所述的基于GPU加速的高精度PET重建方法,其特征在于,所述方法步骤的计算过程均在GPU中并行处理,将核函数在每个线程中同时并行计算;在线计算LOR穿过体素深度时,根据系统几何结构的对称性,将计算结果重复多次使用。
  8. 根据权利要求7所述的基于GPU加速的高精度PET重建方法,其特征在于,迭代重建过程中,在进行投影之前及反投影得到校正图像后先对估计图像进行卷积操作,具体是将每个估计图像体素值与该体素值附近的等效系统函数相乘后相加,得到新的该体素位置图像值。
  9. 根据权利要求8所述的基于GPU加速的高精度PET重建方法,其特征在于,测量点源在系统中重建图像上的系统函数,用这些物理效应对点源图像的影响来等效;等效系统响应函数H(x,y,z)是一个类三维高斯函数:
    Figure PCTCN2019125416-appb-100006
    其中,r,t,a分别代表坐标轴的径向、切向和轴向,传播参数sigma在拟合中被估计,得到三个参数σ r、σ t、σ a的值后,这三个参数代表系统响应函数的模型,σ r、σ t、σ a分别表示径向、切向和轴向的传播参数。
  10. 一种基于GPU加速的高精度PET重建装置,其特征在于,包括:
    正投影单元,用于输入初始的估计图像,对估计图像进行正投影操作,获得估计投影;
    比较单元,用于将获得的估计投影与实际测量投影作比较,获得比较结果;
    反投影单元,用于对比较结果进行反投影操作,获得校正图像;
    迭代单元,用于将校正图像与原始图像对应体素相乘,然后除敏感度图像,完成一次迭代,得到的图像就是下一次迭代的输入估计图像。
PCT/CN2019/125416 2019-12-14 2019-12-14 一种基于gpu加速的高精度pet重建方法及装置 WO2021114306A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/CN2019/125416 WO2021114306A1 (zh) 2019-12-14 2019-12-14 一种基于gpu加速的高精度pet重建方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2019/125416 WO2021114306A1 (zh) 2019-12-14 2019-12-14 一种基于gpu加速的高精度pet重建方法及装置

Publications (1)

Publication Number Publication Date
WO2021114306A1 true WO2021114306A1 (zh) 2021-06-17

Family

ID=76328818

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2019/125416 WO2021114306A1 (zh) 2019-12-14 2019-12-14 一种基于gpu加速的高精度pet重建方法及装置

Country Status (1)

Country Link
WO (1) WO2021114306A1 (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113744358A (zh) * 2021-08-30 2021-12-03 上海交通大学 一种基于gpu加速的三维粒子场及速度场重建方法
CN115330900A (zh) * 2022-10-13 2022-11-11 中加健康工程研究院(合肥)有限公司 一种用于pet图像重构的系统矩阵快速迭代计算方法
CN117495990A (zh) * 2024-01-02 2024-02-02 中国科学技术大学 基于单片fpga的pet正弦图数据压缩存储方法、系统及设备

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102831627A (zh) * 2012-06-27 2012-12-19 浙江大学 一种基于gpu多核并行处理的pet图像重建方法
US20130101194A1 (en) * 2011-10-21 2013-04-25 University Of Utah Research Foundation Filtered Backprojection Image Reconstruction with Characteristics of an Iterative Map Algorithm
CN105046744A (zh) * 2015-07-09 2015-11-11 中国科学院高能物理研究所 基于gpu加速的pet图像重建方法
CN107123095A (zh) * 2017-04-01 2017-09-01 上海联影医疗科技有限公司 一种pet图像重建方法、成像系统
CN107220924A (zh) * 2017-04-11 2017-09-29 西安电子科技大学 一种基于gpu加速pet图像重建的方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130101194A1 (en) * 2011-10-21 2013-04-25 University Of Utah Research Foundation Filtered Backprojection Image Reconstruction with Characteristics of an Iterative Map Algorithm
CN102831627A (zh) * 2012-06-27 2012-12-19 浙江大学 一种基于gpu多核并行处理的pet图像重建方法
CN105046744A (zh) * 2015-07-09 2015-11-11 中国科学院高能物理研究所 基于gpu加速的pet图像重建方法
CN107123095A (zh) * 2017-04-01 2017-09-01 上海联影医疗科技有限公司 一种pet图像重建方法、成像系统
CN107220924A (zh) * 2017-04-11 2017-09-29 西安电子科技大学 一种基于gpu加速pet图像重建的方法

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113744358A (zh) * 2021-08-30 2021-12-03 上海交通大学 一种基于gpu加速的三维粒子场及速度场重建方法
CN113744358B (zh) * 2021-08-30 2024-02-20 上海交通大学 一种基于gpu加速的三维粒子场及速度场重建方法
CN115330900A (zh) * 2022-10-13 2022-11-11 中加健康工程研究院(合肥)有限公司 一种用于pet图像重构的系统矩阵快速迭代计算方法
CN117495990A (zh) * 2024-01-02 2024-02-02 中国科学技术大学 基于单片fpga的pet正弦图数据压缩存储方法、系统及设备
CN117495990B (zh) * 2024-01-02 2024-05-17 中国科学技术大学 基于单片fpga的pet正弦图数据压缩存储方法、系统及设备

Similar Documents

Publication Publication Date Title
Merlin et al. CASToR: a generic data organization and processing code framework for multi-modal and multi-dimensional tomographic reconstruction
Shackleford et al. On developing B-spline registration algorithms for multi-core processors
AU2019238323B2 (en) Deep encoder-decoder models for reconstructing biomedical images
WO2021114306A1 (zh) 一种基于gpu加速的高精度pet重建方法及装置
US8620054B2 (en) Image reconstruction based on accelerated method using polar symmetries
CN110811667B (zh) 一种基于gpu加速的高精度pet重建方法及装置
US8464026B2 (en) Method and apparatus for computing massive spatio-temporal correlations using a hybrid CPU-GPU approach
Cui et al. Distributed MLEM: An iterative tomographic image reconstruction algorithm for distributed memory architectures
Pedemonte et al. GPU accelerated rotation-based emission tomography reconstruction
Kösters et al. EMRECON: An expectation maximization based image reconstruction framework for emission tomography data
Keck et al. GPU-accelerated SART reconstruction using the CUDA programming environment
WO2014016626A1 (en) Method, computer readable medium and system for tomographic reconstruction
Chen et al. A hybrid architecture for compressive sensing 3-D CT reconstruction
Nguyen et al. GPU-accelerated 3D Bayesian image reconstruction from Compton scattered data
Liu et al. GPU-based acceleration for interior tomography
Schellmann et al. Parallel medical image reconstruction: from graphics processing units (GPU) to grids
WO2024109762A1 (zh) Pet参数确定方法、装置、设备和存储介质
Xie et al. An effective CUDA parallelization of projection in iterative tomography reconstruction
Weinlich et al. Comparison of high-speed ray casting on GPU using CUDA and OpenGL
Caucci et al. GPU programming for biomedical imaging
Sampson et al. Investigating multi-threaded SIMD for helical CT reconstruction on a CPU
EP3493163B1 (en) Compressive sensing of light transport matrix
Chou et al. Accelerating image reconstruction in dual-head PET system by GPU and symmetry properties
Pratx et al. 3-D tomographic image reconstruction from randomly ordered lines with CUDA
Valencia Pérez et al. Parallel approach to tomographic reconstruction algorithm using a Nvidia GPU

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19955971

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19955971

Country of ref document: EP

Kind code of ref document: A1

32PN Ep: public notification in the ep bulletin as address of the adressee cannot be established

Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205A DATED 20.01.2023)

122 Ep: pct application non-entry in european phase

Ref document number: 19955971

Country of ref document: EP

Kind code of ref document: A1