List mode 3D PET image reconstruction method based on learnable gradient descent algorithm

Abstract

The invention discloses a list mode 3D PET image reconstruction method based on a learnable gradient descent algorithm, the list mode is directly used for image reconstruction, and compared with a traditional method for reconstruction through sinogram data, the method occupies less memory and can integrate flight time difference information. Therefore, the data processing efficiency is improved, the requirement for hardware resources is lowered, and the method is more suitable for the resource-limited environment. Meanwhile, the method also has the capability of reconstructing a 3D PET image, and compared with a 2D list mode reconstruction method, the method has the advantages that video memory occupation and model parameter quantity are optimized, and the method can be operated on an existing hardware platform. The method is unique in that the method does not need to depend on labeled training data, greatly reduces the difficulty of medical image data acquisition, and can learn the internal structure and characteristics of the image from non-labeled data, thereby reducing the threshold in practical application.

Description

Technical Field

[0001] The present invention belongs to the field of medical image processing based on deep learning, and particularly relates to a list-mode 3D PET image reconstruction method based on a learnable gradient descent algorithm. Background Art

[0002] Positron Emission Tomography (PET) is a nuclear medicine imaging technology that reflects physiological function information. A tracer labeled with a positron radioactive isotope (such as fluorine-18 ( 18 F), oxygen-15 ( 15 O), nitrogen-13 ( 13 N), and carbon-11 ( 11 C)) will decay in living tissues after being injected into a patient. The released positrons will annihilate with negative electrons in the body, generating a pair of γ photon pairs with an energy of 511 keV and opposite directions. A PET detector ring captures a pair of gamma photon pairs as a coincidence event. The coincidence event indicates that the annihilation reaction occurs on the line of response (LOR), but the specific location cannot be determined. Therefore, PET imaging requires an image reconstruction process.

[0003] Due to its uniqueness in measuring physiological functions and molecular imaging, PET imaging technology plays a wide and irreplaceable role in clinical applications. For example, in oncology [1-2] , PET imaging is crucial for tumor detection, cancer grading, treatment response evaluation, and cancer recurrence detection; in neurology, PET imaging, with its advantage of non-invasively observing brain activities, plays a key role in diagnosing neurological diseases such as Alzheimer's disease [3] , epilepsy [4] . In addition, the clinical application of PET in cardiology [5] is also continuously expanding, and it is widely used in basic research, such as brain function [6] .

[0004] List-mode data is to record various information of each annihilation reaction in the form of a list, such as the angle, energy, photon reception position of the detector, and occurrence time of the LOR. The time-of-flight (TOF) information is also very important for improving the quality of PET image reconstruction. If the distance from the annihilation reaction position to the midpoint of the LOR is denoted as d, according to symmetry:

[0005] Δt = t2 - t1;

[0006]

[0007] Ideally, the position of the coincidence event on the LOR can be determined by Δt. However, due to insufficient time precision, in the Non-TOF PET model, it is necessary to assume that the coincidence events are uniformly distributed along the LOR, resulting in more detected noise. With the improvement of time resolution, the TOF-PET method divides the LOR into multiple TOF-bins, determines the position with the highest probability of coincidence events according to the Δt information, and establishes a Gaussian distribution probability model to improve the image reconstruction quality.

[0008] If we want to introduce TOF information into the sinogram data of PET, the original two dimensions of the sinogram represent the projection angle and the distance from the center point respectively, and cannot accommodate TOF information. Then a new dimension needs to be introduced into the sinogram to represent which TOF-bin the coincidence event is within this LOR. When using the sinogram for PET image reconstruction, the data volume is too large to be stored in memory. [7] In this case, the list-mode data does not have the problem of data empty bins, and there is no need to expand the time dimension information, becoming a promising solution strategy for TOF-PET image reconstruction.

[0009] With the emergence of deep learning, some progress has been made in list-mode PET image reconstruction. FastPET proposed by Whiteley et al. [8] uses an end-to-end network, but uses a histogram, does not directly use list-mode data, and has no physical constraints, so it is not stable enough. LMPD-Net proposed by Li et al. [9] first applies the model-based deep learning method, but this network only focuses on the reconstruction of 2D images, and integrating TOF information occupies a large amount of memory. LM-DIPRecon proposed by Ote et al. [7] directly reconstructs from List-mode data and improves the image quality. This network combines the DRAMA and DIP algorithms, but this is a post-denoising method. The neural network does not incorporate TOF information and the projection process, and the reconstruction time is long.

[0010] References:

[0011] [1] SOTOUDEH H, SHARMA A, FOWLER K J, et al. Clinical application of PET / MRI in oncology. [J]. Journal of magnetic resonance imaging: JMRI, 2016, 44(2): 265 - 276.

[0012] [2]MEIKLE S R,SOSSI V,RONCALI E,et al.Quantitative PET in the 2020s:aroadmap.[J].Physics inmedicine andbiology,2021,66(6):06RM01.

[0013] [3]HERHOLZ K.PET studies in dementia.[J].Annals ofnuclearmedicine,2003,17(2):79-89.

[0014] [4]SARIKAYA I.PET studies in epilepsy.[J].Americanjournal ofnuclearmedicine and molecular imaging,2015,5(5):416-430.

[0015] [5]BENGEL F M,HIGUCHI T,JAVADI M S,et al.Cardiac Positron EmissionTomography[J].Journal oftheAmerican College ofCardiology,2009,54(1):1-15.

[0016] [6]ONISHI Y,ISOBE T,ITO M,et al.Performance evaluation ofdedicatedbrain PET scanner with motioncorrection system.[J].Annals ofnuclearmedicine,2022,36(8):746-755.

[0017] [7]OTE K,HASHIMOTO F,ONISHI Y,et al.List-Mode PET ImageReconstruction Using Deep ImagePrior[J].IEEE Transactions onMedical Imaging,2023,42(6):1822-1834.

[0018] [8] WHITELEY W, PANIN V, ZHOU C, et al. FastPET: Near Real-Time Reconstruction of PET Histo-Image Data Using a Neural Network[J]. IEEE Transactions on Radiation and Plasma Medical Sciences, 2021, 5(1): 65-77.

[0019] [9] LI C, HU R, CUI J, et al. LMPDNet: TOF-PET list-mode image reconstruction using model-based deep learning method[A]. arXiv, 2023.

[0020]

[10] CHEN Y, LIU H, YE X, et al. Learnable Descent Algorithm for Nonsmooth Nonconvex Image Reconstruction[A]. arXiv, 2022.

[0021]

[11] COCOSCO CA, KOLLOKIAN V, KWAN R KS, et al. BrainWeb: Online Interface to a 3D MRI Simulated Brain Database[J]. NeuroImage, 1997.

[0022]

[12] MEHRANIAN A, READER A J. Model-Based Deep Learning PET Image Reconstruction Using Forward–Backward Splitting Expectation–Maximization[J]. IEEE Transactions on Radiation and Plasma Medical Sciences, 2021, 5(1): 54-64.

[0023]

[13] NESTEROV Yu. Smooth minimization of non-smooth functions[J]. Mathematical Programming, 2005, 103(1): 127-152.

[0024]

[14] SCHRAMM G, THIELEMANS K. PARALLELPROJ—an open-source framework for fast calculation of projections into tomography[J]. Frontiers in Nuclear Medicine, 2024, 3: 1324562.

[0025]

[15] POLSON L A, FEDRIGO R, LI C, et al. PyTomography: A Python Library for Medical Image Reconstruction[A]. arXiv, 2024.

[0026]

[16] JOSEPH P M. An Improved Algorithm for Reprojecting Rays through Pixel Images[J]. IEEE Transactions on Medical Imaging, 1982, 1(3): 192-196. Summary of the Invention

[0027] Based on the above, in the present invention, an interpretable variational model is proposed for reconstructing TOF-PET 3D images from list-mode data, a list-mode 3D PET image reconstruction method based on a learnable gradient descent algorithm. Combining the learnable l 2,1 -norm to extract sparse features in the PET image, which makes the network interpretable, robust, and can further improve the image reconstruction quality. The architecture of the proposed multi-stage neural network fully follows the convergent learning gradient descent algorithm (Learned Descent Algorithm, LDA)

[10] , thus inheriting the convergence characteristics of the algorithm.

[0028] The present invention is realized by the following technical solutions:

[0029] The present invention discloses a list-mode 3D PET image reconstruction method based on a learnable gradient descent algorithm, including the following steps:

[0030] Obtain at least hundreds of list-mode PET data;

[0031] Calculate the corresponding normalized ground truth image through a traditional PET image reconstruction algorithm and the obtained list-mode PET data;

[0032] Generate list-mode - image data pairs from list-mode PET data and corresponding normalized ground-truth images;

[0033] Divide the processed list-mode - image data pairs into a training set, a validation set, and a test set;

[0034] Construct an objective function using a fidelity term and a regularization term;

[0035] Define the fidelity term as the list-mode log-likelihood function;

[0036] Parameterize the regularization term as the l 2,1 norm and extract sparse features using a learnable feature extraction operator;

[0037] Introduce a smoothing parameter ε to handle the non-smoothness problem of the objective function caused by the non-differentiability of the l 2,1 norm at the origin;

[0038] Perform residual-type updates on the fidelity term and the regularization term respectively and calculate the convex combination of the two updates as the alternative parameter u k+1 ;

[0039] Calculate the standard gradient descent of the objective function at the reconstructed image as another alternative parameter v k+1 ;

[0040] Select the parameter corresponding to the smaller objective function among u k+1 and v k+1 as the update result and update the smoothing parameter;

[0041] Compare the image reconstructed using sparse features with the corresponding normalized ground-truth image and iteratively update the network weights using the loss between them.

[0042] As a further improvement, the objective function constructed using the fidelity term and the regularization term in the present invention formulates the inverse problem of PET image reconstruction as an optimization problem in a variational framework, specifically:

[0043]

[0044] where L(l|x) is the fidelity term, representing the list-mode log-likelihood function, r(x; θ) is the regularization term, x is the PET activity map to be reconstructed, θ is the learnable parameter, and l is the measured list-mode data expressed as:

[0045] l = {i(t)|t = 1, 2,..., N};

[0046] where i(t) is the index of the t-th event measured on the line where the annihilation reaction occurs, and N is the total number of events.

[0047] As a further improvement, in the present invention, the regularization term is parameterized as an l 2,1 -norm and a learnable feature extraction operator is used to extract sparse features, obtaining a robust and effective regularization method, specifically:

[0048] The r in the objective function is parameterized as an l 2,1 -norm, and a feature extraction operator g(x) learned from the training data is used to extract sparse features, expressed as follows:

[0049] r(x; θ) = ∑ j ‖g j (x; θ)‖₂ = ‖g(x; θ)‖₂;

[0050] where θ is a learnable parameter, and g j (x; θ) is the vector at position j on all channels. A multi-layer convolutional neural network with a non-linear activation function σ(x) is selected as g, and σ(x) is the smoothed ReLU function:

[0051]

[0052] As a further improvement, the present invention introduces a smoothing parameter ε to handle the non-smoothness problem of the objective function caused by the non-differentiability of the l 2,1 -norm at the origin, enabling the direct calculation of the gradient of the regularization term, specifically:

[0053] Due to the non-differentiability of the l 2,1 -norm at the origin, Nesterov's smoothing technique is first used to handle this problem:

[0054]

[0055] where <·,·> represents the inner product operation, is a dual variable, represents a vector space;

[0056] For any ε > 0, consider using the perturbed dual form to smooth r to obtain r ∈ :

[0057]

[0058] Note that the perturbed dual form of the above formula has a corresponding closed-form solution:

[0059]

[0060] Solving the above formula gives the closed-form solutions for each part in :

[0061]

[0062] Substitute into r ε The processed formula of (x; θ) is as follows:

[0063]

[0064] where J0 = {j ∈ [m] | ‖g j (x; θ)‖ ≤ ε}, J1 = [m]\J0; Calculate the gradient as follows:

[0065]

[0066] The parameter ε is a smoothing factor that controls the approximation degree of the original term. By applying a smooth approximation to the l 2,1 norm, the non-smoothness problem is effectively solved, and the gradient is calculated directly.

[0067] As a further improvement, in the present invention, the fidelity term and the regularization term are respectively updated in a residual type, and the convex combination of the two updates is calculated as the alternative parameter u k+1 for more effectively updating the reconstructed image, specifically:

[0068] Rewrite the objective function using the learnable regularization term given above:

[0069]

[0070] In each iteration stage k, first construct the gradient descent of -L at x k denoted as z k+1 :

[0071]

[0072] In this case, separate -L and r εk and let them participate in the residual type update respectively;

[0073] Now consider the gradient descent at z k+1 denoted as p k+1 :

[0074]

[0075] Finally, consider the convex combination of the two functions z k+1 and p k+1 as follows:

[0076]

[0077] where, u k+1 plays a key role in the algorithm, achieving high efficiency through step-by-step use of the residual type updates of -L and r ∈ and avoiding the vanishing gradient when minimizing the loss function.

[0078] As a further improvement, the standard gradient descent of the objective function calculated by the present invention at the reconstructed image is used as another update alternative parameter v for the algorithm k+1 , ensuring the convergence of the algorithm, specifically:

[0079] The existing gradient descent of -L at x k is denoted as z k+1 . Now consider the gradient descent at x k is defined as v k+1 , as an alternative for x k+1 , and the corresponding closed-form solution is obtained as follows:

[0080]

[0081] v k+1 is the standard gradient descent of the objective function φ εk at x k , ensuring the convergence of the algorithm.

[0082] As a further improvement, the present invention selects the parameter corresponding to the smaller objective function among u k+1 and v k+1 as the update result and updates the smoothing parameter to improve the training efficiency, specifically as follows:

[0083] Select the one with the smaller objective function value φ k+1 between u k+1 and v εk as x k+1 :

[0084]

[0085] Meanwhile, the smoothing parameter ε k decays according to the following criterion:

[0086]

[0087] As a further improvement, the present invention compares the image reconstructed using sparse features with the corresponding normalized ground truth image and iteratively updates the network weights using the loss between them to gradually reduce the reconstruction error, specifically as follows:

[0088] The formula of the loss function is as follows:

[0089] L = MSE(xK, x ref ) + μ(1 - SSIM(x K , x ref ));

[0090] Among them, x ref is the reference PET image, x K is the output image, MSE refers to the Mean - Square Error, μ is the weight of the Structural Similarity (SSIM) loss;

[0091] During the training process, the generated reconstructed image is compared with the corresponding normalized ground - truth image, the loss between them is calculated, and the weights of the network are iteratively updated through the back - propagation algorithm to gradually reduce the reconstruction error until the network performance reaches a certain satisfactory level or reaches the preset number of iterations.

[0092] The beneficial effects of the present invention are as follows:

[0093] The method of the present invention directly uses list - mode for image reconstruction. Compared with the traditional method of reconstructing through sinogram data, this method occupies less memory and can integrate time - of - flight difference information. This not only improves the efficiency of data processing but also reduces the requirements for hardware resources, making the method more suitable for resource - constrained environments. At the same time, the method of the present invention also has the ability to reconstruct 3D PET images. Compared with the 2D list - mode reconstruction method, this method has optimized both the video - memory occupancy and the number of model parameters and can run on the existing hardware platform. The regularization term in the objective function mentioned in the present invention combines the learnable l 2,1 norm, which can more effectively extract the sparse features in the PET image and is more robust, significantly improving the resistance of the method to noise and anomalies. Secondly, the training strategy mentioned in the present invention belongs to the category of unsupervised learning. Its uniqueness lies in that it does not rely on labeled training data, greatly reducing the difficulty of obtaining medical imaging data. If through traditional supervised learning, medical imaging data needs to be labeled, which is not only costly but also requires professional knowledge, resulting in limited medical imaging data available for deep learning; through unsupervised learning, the method of the present invention can learn the internal structure and features of the image from unlabeled data, thus reducing the threshold in practical applications. It should be noted that the smoothing technique mentioned in the present invention cleverly solves the problem of non - smoothness of the objective function caused by the non - differentiability of the l 2,1 norm at the origin, greatly reducing the number of network parameters and making the entire reconstruction process more efficient. In addition, the method proposed in the present invention performs well in both quantitative evaluation metrics and visual results. Description of the Drawings

[0094] Figure 1 It is a flow chart of the steps of the learnable gradient descent algorithm;

[0095] Figure 2 It is a schematic diagram of the learnable gradient descent algorithm for list-mode PET image reconstruction;

[0096] Figure 3 They are the normalized true-value images and the reconstructed images of the 1st, 3rd, and 12th phantoms by different methods;

[0097] Figure 4 They are the normalized true-value images and the reconstructed images of the 1st, 3rd, and 12th phantoms obtained by the proposed method at different count levels;

[0098] Figure 5 It is a trade-off curve graph between the contrast and noise of PET image reconstruction by different methods;

[0099] (a) Contrast recovery coefficient - standard deviation curve; (b) Tumor pixel value ratio - signal-to-noise ratio curve. Detailed implementation manners

[0100] To describe the present invention more specifically, the technical solutions of the present invention will be described in detail below in conjunction with the accompanying drawings and specific implementation manners.

[0101] This method provides a list-mode 3D PET image reconstruction method based on the learnable gradient descent algorithm, which can improve the quality of PET image reconstruction. As Figure 1 shown, it includes the following steps:

[0102] (1) Preprocessing of data.

[0103] The entire process of PET image reconstruction starts with data acquisition using a PET scanner, and then through a series of complex processing steps, finally reconstructs the concentration image required by clinicians. In this method, simulation experiments are used to verify the performance of the algorithm, and this verification process is different from the actual process of obtaining real data, specifically manifested in obtaining PET list-mode data by projecting phantoms.

[0104] (1-1) Segmenting the image boundary

[0105] Use 20 3D brain models in BrainWeb

[11] to simulate 3D 18 F-FDG PET images. According to the simulation settings described in previous research

[12] , segment the T1-weighted MR image in BrainWeb into gray matter (GM), white matter (WM), cerebrospinal fluid, skull, and skin.

[0106] (1-2) Assign and generate PET images

[0107] For each image, the FDG PET model is generated as follows: Random uptake values of 96.0 ± 5.0 and 32.0 ± 5.0 are assigned to the GM and WM regions, respectively. Ten spherical hot regions with a radius ranging from 2 mm to 10 mm are embedded in all phantoms. The uptake value of the hot lesion is 144 (1.5 × GM), and the uptake value of the cold lesion is 48.0 (0.5 × WM). The PET image is generated from voxels with a resolution of 2.086 × 2.086 × 2.031 mm 3 and the generated matrix size is 136 × 136 × 127.

[0108] (1-3) Select slices to produce list-mode data

[0109] Two sets of ten consecutive slices are selected from each of the three orthogonal views of each model to generate TOF list-mode data with different photon count levels (5e5, 1e6, and 5e6). The relationship between the PET image and the raw data can be described by the following linear equations and Poisson distribution model:

[0110] y = Poisson(Ax + b);

[0111] That is

[0112]

[0113] where x = (x1, x2,..., x J ) T is the vector of image voxel values, i.e., the PET image obtained by simulation here, y = (y1, y2,..., y I ) T is the vector of actual sampling values of the projection data, b = (b1, b2,..., b I ) T is the vector of error values caused by scattering and random coincidence, etc. during the projection process, A ∈ R I×J is the system matrix, where each element a ij represents the probability that the photons emitted by the j-th voxel are detected by the i-th LOR.

[0114] (1-4) Produce list-mode - image data pairs and partition the dataset

[0115] Finally, every 10 slice data are integrated into a 3D image. The input 3D image has 128×128×10 pixels, and the voxel size is 2mm×2mm×2.1mm. The coincidence time resolution of the detector is set to 200 pS. Each LOR is divided into 21 TOF-bins, and the width of each TOF-bin is set to 12.5 mm. A total of 120 3D list-mode - image pairs are produced. Finite spatial resolution, attenuation, and 20% uniform noise are considered to simulate random and scattered coincidence events. Among them, 17 phantoms (102 list-mode - image pairs) are selected as training data, 2 phantoms (12 list-mode - image pairs) are used as test data, and 1 phantom (6 list-mode - image pairs) is used as validation data.

[0116] (2) Construct the objective function.

[0117] In TOF-PET list-mode reconstruction, the measured list-mode data are represented as:

[0118] l = {i(t)|t = 1, 2,..., N};

[0119] where i(t) is the index of the t-th event measured on the LOR, and N is the total number of events. PET image reconstruction, as a classical inverse problem, aims to invert the concentration distribution map x of the radioactive tracer from the detection data under the influence of various noises, and can be formulated as an optimization problem in the variational framework as follows:

[0120]

[0121] This objective function includes a fidelity term -L(l|x) based on the PET imaging model and a regularization term r(x; θ) for extracting the sparse features of the image.

[0122] (2-1) Fidelity term module

[0123] The fidelity term L(l|x) represents the list-mode log-likelihood function and is expressed as follows:

[0124]

[0125] where x is the PET activity map to be reconstructed. A LM ∈R I×J represents the list-mode TOF forward projection model, also known as the system response matrix, and a ij represents the response of voxel j to the i-th data bin. b is the expectation of the background events (random and scattered events) in the i-th data bin. s = ∑ j s j = ∑ j ∑ j a ijIt is a sensitivity image.

[0126] (2-2) Regular term module

[0127] Parameterize the regular term r as the l 2,1 norm, and use the feature extraction operator g(x) that can be learned from the training data to extract sparse features, expressed as follows:

[0128]

[0129] where θ is the learnable parameter, and g j (x; θ) is the vector at position j on all channels. Select a convolutional neural network with a non-linear activation function σ(x) as g. g contains 4 convolutional layers and 96 channels, and σ(x) is the smoothed ReLU:

[0130]

[0131] (3) Handling the non-smoothness problem of the objective function caused by using the l 2,1 norm.

[0132] Since the non-differentiability of the l 2,1 norm at the origin leads to the non-smoothness of the objective function, adopt Nesterov's smoothing technique

[13] to handle this problem:

[0133]

[0134] where <·,·> represents the inner product operation, is a dual variable, represents a vector space;

[0135] The parameter ε is a smoothing factor that controls the approximation degree to the original term. For any ε > 0, consider using the perturbed dual form to smooth r to get r ε :

[0136]

[0137] Note that the perturbed dual form of the above formula has a corresponding closed-form solution:

[0138]

[0139] Solving the above formula can obtain the closed-form solution of each part in :

[0140]

[0141] Substitute into r ε(x; θ) is processed as follows:

[0142]

[0143] Where J0 = {j ∈ [m] | ‖g j (x; θ)‖ ≤ ε}, J1 = [m] \ J0; Calculate the gradient as follows:

[0144]

[0145] By applying a smooth approximation to the l 2,1 norm, the non-smoothness problem is effectively solved, and the regular term gradient can be directly calculated.

[0146] (4) Construct a learnable gradient descent algorithm for list-mode PET.

[0147] Rewrite the objective function using the learnable regular term given above:

[0148]

[0149] (4-1) Apply the proximal gradient method to calculate the updated alternative parameter u k+1

[0150] As Figure 2 shown, assume the algorithm iterates for a total of K stages. In each iteration stage k, apply the proximal gradient method to solve the above optimization problem. First, construct the gradient descent of -L at x k denoted as z k+1 :

[0151]

[0152] Define the proximal operator for further iterative solution:

[0153]

[0154] To find a closed-form solution for x that satisfies the above conditions within the neighborhood of z k+1 , perform a Taylor expansion of r εk at z k+1 , and denote the expansion as

[0155]

[0156] Substitute into the definition formula of the proximal operator, and combine the high-order terms to obtain the following formula:

[0157]

[0158] To obtain the closed-form solution of the minimum value, take the derivative of \(x\) in the above equation and set its derivative to 0. Denote That is where \(\alpha_0 = 0.01\) and \(\beta_0 = 0.02\), and the following formula is obtained:

[0159]

[0160] Denote the closed-form solution obtained from the above equation as \(u\) k+1 , as follows:

[0161]

[0162] The process of obtaining \(u\) above k+1 can also be regarded as separating \(-L\) and \(r\) εk and each participating in the residual-type update; considering the gradient descent at \(z\) k+1 , denoted as \(p\) k+1 :

[0163]

[0164] Finally, consider the convex combination of the two functions \(z\) k+1 and \(p\) k+1 , as follows:

[0165]

[0166] (4 - 2) Calculate the standard gradient descent of the objective function at \(x\) k as another update alternative parameter \(v\) k+1

[0167] Similarly, we can consider another closed-form solution \(v\) k+1 , as follows:

[0168]

[0169] \(v\) k+1 is the standard gradient descent of the objective function \(\varphi\) εk at \(x\) k , which guarantees the convergence of the algorithm.

[0170] (5) Update \(x\) k+1 and the smoothing parameter \(\varepsilon\) k

[0171] Select the one with the smaller objective function value \(\varphi\) k+1 in \(u\) k+1 and \(v\) εk as the next \(x\) k+1 :

[0172]

[0173] Meanwhile, the smoothing parameter ε k decays according to the following criteria:

[0174]

[0175] (6) Training process

[0176] The model was trained using PyTorch 2.0 on an NVIDIA Quadro RTX 8000 GPU (48GB). The loss function used to train the regularization network is defined as follows:

[0177] L = MSE(x K , x ref ) + μ(1 - SSIM(x K , x ref ));

[0178] where x ref is the reference PET image, x K is the output image, MSE refers to the Mean-Square Error, μ is the weight of the Structural Similarity SSIM loss, which is set to 150 in the experiment. The Adam optimizer was used with a learning rate of 3×10 -6 . The model was trained for 300 epochs with a batch size of 1, unfolded into a total of 4 stages, and the image obtained after 10 iterations using the Maximum Likelihood Expectation Maximization (MLEM) algorithm was used as the initial image x0. The open-source framework named parallelproj [14-15] was used to implement accelerated projection, which is based on the Joseph projection method

[16] to parallelize the calculation of each LOR of the incoming list-mode events.

[0179] (7) Testing process

[0180] The trained model was compared with the List-Mode Ordered Subsets Expectation Maximization (LM-OSEM), List-Mode Stochastic Primal-Dual Hybrid Gradient (LM-SPDHG), and the histogram-based list-mode image reconstruction method Fast-PET. The subset of LM-OSEM was set to 4, and the subset of LM-SPDHG was set to 112. For LM-OSEM, the algorithm converges when the number of iterations is 5; while for LM-SPDHG, the algorithm converges when the number of iterations is 25. For LM-SPDHG, the total variation regularization parameter was set to 5.0. The Fast-PET neural network uses a U-NET architecture consisting of 31 convolutional layers. The histogram and attenuation map were used as inputs to Fast-PET simultaneously. The learning rate of Fast-PET was set to 1×10 -5, the batch size is 1, and the total number of training epochs is 500. The processing effects of using different reconstruction methods compared with the proposed method are as Figure 3 shown, Figure 3 The photon counting level is 1e6, and the TOF resolution is 200 pS. From left to right: normalized ground truth image, LM-OSEM, LM-SPDHG, Fast-PET, proposed method. The fourth row shows the enlarged image of the red square area.

[0181] The results of the learnable gradient descent algorithm have good visual effects at three different photon counting levels, as Figure 4 shown. From left to right: normalized ground truth image, 5e5 photon counting level, 1e6 photon counting level, and 5e6 photon counting level.

[0182] Analyze and compare the Contrast Recovery Coefficient (CRC), Signal-to-Noise Ratio (SNR), Normalized Standard Deviation (NSTD), and Tumor Value Ratio (TR) to measure the model performance. The related definitions are as follows:

[0183]

[0184] K a is the number of Regions of Interest (ROIs) on gray matter and tumors, and a total of 5 regions are selected; K b is the number of ROIs on white matter, and a total of 10 regions are selected; N a,k is the voxel value of the k-th ROI region on gray matter and tumors, and N b,k is the voxel value of the k-th ROI region on white matter.

[0185]

[0186] Among them, x ref is the ground truth image, R is all regions of the brain, and N R is the total number of voxels in the brain.

[0187]

[0188] R tumor is the set of 4 tumor regions, and N tumor is the number of pixels in the 4 tumor regions.

[0189] The related curves are as Figure 5As shown, for LM-OSEM, one point is plotted for each iteration. For LM-SPDHG, one point is plotted for every 5 iterations. For Fast-PET, one point is plotted for every 100 iterations, and for the proposed method (LM-LDA), one point is plotted for every 50 iterations.

[0190] Figure 5 (a) is the curve of contrast recovery coefficient - standard deviation. The smaller the standard deviation and the larger the contrast recovery coefficient, that is, the closer the result in the curve is to the upper left corner, the better the image reconstruction effect of the method; Figure 5 (b) is the tumor pixel value ratio - signal-to-noise ratio image. The higher the signal-to-noise ratio and the higher the tumor pixel value ratio, that is, the closer the result in the curve is to the upper right corner, the better the tumor recovery effect of the method and the higher the image reconstruction quality. As can be seen from Figure 5 it, in the evaluation of the above indicators, the proposed method achieves better results than other methods.

[0191] The trained model is used for PET image reconstruction testing, and quantitative analysis is performed using Peak Signal-to-Noise Ratio (PSNR) and SSIM. The results are shown in Table 1, and the relevant definitions are as follows:

[0192]

[0193] Among them, MAX represents the maximum value of the pixel points on the true image, μ rec and σ rec represent the mean and variance of the reconstructed image; μ ref and σ ref represent the mean and variance of the true image. σ rec,ref represents the covariance between the reconstructed image and the true image. The constants c1 and c2 are defined as c1 = 0.01×max(x rec ) and c2 = 0.03×max(x rec ).

[0194] Table 1 Quantitative analysis of different methods at different counting levels.

[0195]

[0196] The above description of the embodiments is for those of ordinary skill in the art of this technology to understand and apply the present invention. Those familiar with the technology can obviously make various modifications to the above embodiments easily and apply the general principles described herein to other embodiments without creative labor. Therefore, the present invention is not limited to the above embodiments, and the improvements and modifications made by those skilled in the art according to the disclosure of the present invention should be within the protection scope of the present invention.

Claims

1. A list-mode 3D PET image reconstruction method based on a learnable gradient descent algorithm, characterized in that It includes the following steps: Obtain at least hundreds of list-mode PET data; Calculate the corresponding normalized ground-truth image through traditional PET image reconstruction algorithms and the obtained list-mode PET data; Make the list-mode PET data and the corresponding normalized ground-truth image into a list-mode-image data pair; Divide the processed list-mode-image data pairs into a training set, a validation set, and a test set; Construct an objective function using a fidelity term and a regularization term; Define the fidelity term as a list-mode log-likelihood function; Parameterize the regularization term as the l 2,1 norm and extract sparse features using a learnable feature extraction operator; Introduce a smoothing parameter ε to handle l 2,1 The non-smoothness problem of the objective function caused by the non-differentiability of the norm at the origin; perform residual-type updates on the fidelity term and the regularization term respectively and calculate the convex combination of the two updates as the alternative parameter u for algorithm update k+1 ; Calculate the standard gradient descent of the objective function at the reconstructed image as another update alternative parameter v k+1 ; Select u k+1 and v k+1 Among them, the parameter corresponding to the smaller objective function is used as the update result and the smoothing parameter is updated; Compare the image reconstructed using sparse features with the corresponding normalized ground-truth image and iteratively update the network weights using the loss between them.

2. The list-mode 3D PET image reconstruction method based on the learnable gradient descent algorithm according to claim 1, wherein The constructing of the objective function using a fidelity term and a regularization term formulates the inverse problem of PET image reconstruction as an optimization problem in a variational framework, specifically: where L(l|x) is the fidelity term, representing the list-mode log-likelihood function, r(x;θ) is the regularization term, x is the PET activity map to be reconstructed, θ is the learnable parameter, and l is the measured list-mode data expressed as: l = {i(t)|t = 1, 2,..., N}; where i(t) is the index of the t-th event measured on the line where the annihilation reaction occurs, and N is the total number of events.

3. The list-mode 3D PET image reconstruction method based on the learnable gradient descent algorithm according to claim 1, characterized in that Parameterize the regularization term as an l 2,1 -norm and use a learnable feature extraction operator to extract sparse features, resulting in a robust and effective regularization method, specifically: Parameterize the r in the objective function as l 2,1 norm, and use the feature extraction operator g(x) learned from the training data to extract sparse features, which is expressed as follows: r(x; θ) = ∑ j ||g j (x; θ)||² = ||g(x; θ)||²; where θ is a learnable parameter, and g j (x; θ) is the vector at position j over all channels. A multi-layer convolutional neural network with a non-linear activation function σ(x) is selected as g, and σ(x) is the smoothed ReLU function:

4. The list-mode 3D PET image reconstruction method based on the learnable gradient descent algorithm according to claim 1 or 2 or 3, characterized in that, The introduction of the smoothing parameter ε to process l 2,1 The non-smoothness problem of the objective function caused by the non-differentiability of the Due to the non-differentiability of the l 2,1 norm at the origin, Nesterov's smoothing technique is first used to handle this problem: where <·,·> represents the inner product operation, is a dual variable, ∥y j ∥≤1, j ∈ [m]} represents a vector space; For any ε > 0, consider r after smoothing r using the perturbed dual form ∈ : Note that the perturbed dual form of the above equation has a corresponding closed-form solution: Solving the above equation gives the closed-form solution for each part in Substitute into r ε (x; θ), the processed formula is as follows: where \(J_0=\{j\in[m]|\|\mathbf{g}\) j (\mathbf{x};\boldsymbol{\theta})\|\leq\varepsilon\}\), \(J_1 = [m]\setminus J_0\); calculate the gradient as follows: The parameter ε is a smoothing factor that controls the degree of approximation to the original term. By applying a smooth approximation to the l 2,1 norm, the non-smoothness problem is effectively solved, and the gradient is calculated directly.

5. The list-mode 3D PET image reconstruction method based on the learnable gradient descent algorithm according to claim 4, wherein Performing residual type updates on the fidelity term and the regularization term respectively and calculating a convex combination of the two updates as the alternative parameter u for the algorithm update k+1 , which is used to update the reconstructed image more effectively, specifically: Rewrite the objective function using the learnable regularization term given above: In each iteration stage k, first construct the gradient descent of -L at x k and denote it as z k+1 : In this case, -L and r εk are separated and each participates in the residual type update; Now consider r εk Gradient descent at zk +1 is denoted as p k+1 : Finally, consider z k+1 and p k+1 The convex combination of the two functions is as follows: Among them, u k+1 plays a key role in the algorithm, achieving high efficiency through the step-by-step update of the residual type of -L and r ∈ to avoid the vanishing gradient when minimizing the loss function.

6. The list mode 3D PET image reconstruction method based on the learnable gradient descent algorithm according to claim 5, wherein Calculating the standard gradient descent of the objective function at the reconstructed image as another update alternative parameter v k+1 , which guarantees the convergence of the algorithm, specifically: The gradient descent of the existing - L at x k is denoted as z k+1 , and now consider the gradient descent at x k is defined as v k+1 , as an alternative to x k+1 , and find the corresponding closed - form solution as follows: v k+1 is the standard gradient descent of the objective function φ εk at x k which guarantees the convergence of the algorithm.

7. The list-mode 3D PET image reconstruction method based on the learnable gradient descent algorithm according to claim 1 or 2 or 3 or 5 or 6, characterized in that, The selected u k+1 and v k+1 Among them, the parameter with the smaller objective function is used as the updated result to update the smoothing parameter, which is used to improve the training efficiency. Specifically, the following steps are taken: Select u k+1 and v k+1 Among them, the one with the smaller objective function value φ εk is used as x k+1 : Meanwhile, the smoothing parameter ε k decays according to the following criterion:

8. The list-mode 3D PET image reconstruction method based on the learnable gradient descent algorithm according to claim 7, wherein The comparing of the image reconstructed using sparse features with the corresponding normalized ground-truth image and the iteratively updating of the network weights using the loss between them to gradually reduce the reconstruction error is as follows: The formula of the loss function is as follows: L = MSE(x K , x ref ) + μ(1 - SSIM(x K , x ref )); where x ref is the reference PET image, x K is the output image, MSE refers to the Mean-Square Error, and μ is the weight of the Structural Similarity (SSIM) loss; During the training process, compare the generated reconstructed image with the corresponding normalized ground-truth image, calculate the loss between them, and iteratively update the weights of the network through the backpropagation algorithm to gradually reduce the reconstruction error until the network performance reaches a certain satisfactory level or reaches the preset number of iterations.

Images