CN117830373B - PET image reconstruction method based on acceleration precondition adjacent gradient algorithm - Google Patents
PET image reconstruction method based on acceleration precondition adjacent gradient algorithm Download PDFInfo
- Publication number
- CN117830373B CN117830373B CN202311817363.XA CN202311817363A CN117830373B CN 117830373 B CN117830373 B CN 117830373B CN 202311817363 A CN202311817363 A CN 202311817363A CN 117830373 B CN117830373 B CN 117830373B
- Authority
- CN
- China
- Prior art keywords
- subset
- order
- pet image
- appga
- matrix
- 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
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 54
- 238000000034 method Methods 0.000 title claims abstract description 44
- 230000001133 acceleration Effects 0.000 title claims abstract description 16
- 238000012636 positron electron tomography Methods 0.000 claims description 72
- 239000011159 matrix material Substances 0.000 claims description 47
- 239000000700 radioactive tracer Substances 0.000 claims description 19
- 239000013598 vector Substances 0.000 claims description 15
- 230000008569 process Effects 0.000 claims description 14
- 238000012879 PET imaging Methods 0.000 claims description 8
- 238000005457 optimization Methods 0.000 claims description 7
- 238000009499 grossing Methods 0.000 claims description 4
- 238000007476 Maximum Likelihood Methods 0.000 claims description 3
- 230000007246 mechanism Effects 0.000 claims description 2
- 238000010200 validation analysis Methods 0.000 claims 1
- 238000005516 engineering process Methods 0.000 abstract description 19
- 238000003745 diagnosis Methods 0.000 abstract description 5
- 201000010099 disease Diseases 0.000 abstract description 3
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 abstract description 3
- 230000006870 function Effects 0.000 description 23
- 238000003384 imaging method Methods 0.000 description 8
- 210000004556 brain Anatomy 0.000 description 6
- 230000005855 radiation Effects 0.000 description 6
- 238000002059 diagnostic imaging Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 230000006378 damage Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 239000000243 solution Substances 0.000 description 4
- 238000004088 simulation Methods 0.000 description 3
- 206010028980 Neoplasm Diseases 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 201000011510 cancer Diseases 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000003902 lesion Effects 0.000 description 2
- 230000007918 pathogenicity Effects 0.000 description 2
- 238000002600 positron emission tomography Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 208000027418 Wounds and injury Diseases 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000001627 detrimental effect Effects 0.000 description 1
- 238000009509 drug development Methods 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 238000001727 in vivo Methods 0.000 description 1
- 238000002347 injection Methods 0.000 description 1
- 239000007924 injection Substances 0.000 description 1
- 208000014674 injury Diseases 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 230000004060 metabolic process Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 238000012831 peritoneal equilibrium test Methods 0.000 description 1
- 238000012877 positron emission topography Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000005295 random walk Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 210000001835 viscera Anatomy 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/50—Depth or shape recovery
- G06T7/55—Depth or shape recovery from multiple images
- G06T7/593—Depth or shape recovery from multiple images from stereo images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10104—Positron emission tomography [PET]
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Quality & Reliability (AREA)
- Nuclear Medicine (AREA)
Abstract
The invention discloses a low-dose high-quality PET image reconstruction method based on a smooth high-order isotropy total variation regularization and acceleration precondition adjacent gradient algorithm (APPGA), which combines a plurality of technologies including an adjacent gradient method, a generalized Nesterov momentum technology and a precondition technology and aims at rapidly reconstructing a high-quality PET image under a low-dose condition. Compared with the traditional algorithm, APPGA has the advantages of improving the image quality, accelerating the reconstruction speed, reducing noise and the like, is expected to provide a faster and higher-quality PET image in the field of medical images, and promotes the development of disease diagnosis and treatment.
Description
Technical Field
The invention relates to the technical field of biomedical image analysis, in particular to a low-radiotracer high-quality PET image reconstruction method based on smooth high-order isotropy total variation regularization and acceleration precondition adjacent gradient algorithm.
Background
Positron emission tomography (Positron Emission Tomography, PET) imaging technology is an important medical imaging technology and is widely applied to the fields of clinical diagnosis, disease research, drug development and the like. The PET image generates a functional image of three-dimensional in-vivo tissues and organs by detecting positron annihilation events generated after the injection of the radiotracer into the body. The imaging technology has obvious advantages for early cancer diagnosis, and can detect focus at abnormal stage of physiological characteristics (such as excessive metabolism) before tissue density change occurs at lesion sites, thereby greatly improving the discovery rate of early cancer. However, the imaging technology has certain disadvantages, and the radioactive tracer used by the imaging technology has a risk of radiation damage to human bodies, and especially the radiation damage to the imaging means caused by infants is more sensitive. In addition, the process of reconstructing the 3D image in the human body by using PET is huge in calculation amount, and the diagnosis efficiency of instruments and equipment is greatly affected by the advantages and disadvantages of a reconstruction algorithm.
Therefore, the rapid PET imaging method with low radiotracer dosage and high imaging quality is provided with great significance in clinic, not only can provide data reference for doctor diagnosis, but also can reduce the injury of the radioactivity of the radiotracer to patients and medical staff; the use efficiency of the instrument and equipment can also be improved, thereby reducing the medical cost of patients. In addition, theoretical research on PET imaging technology can also provide certain theoretical support for independently developing instruments and equipment of this type in China and reducing the imported quantity of instruments, thereby further reducing medical cost.
Disclosure of Invention
In view of the balance problem between speed and image quality in the current PET imaging and the limitations of the prior art, the invention provides a low-radiotracer dose high-quality PET image reconstruction method based on smooth high-order isotropy total variation regularization and acceleration preconditioned adjacent gradient algorithm.
In order to solve the technical problems, the technical scheme of the invention is as follows:
A low-radiotracer high-quality PET image reconstruction method based on smooth high-order isotropy total variation regularization and acceleration precondition adjacent gradient algorithm comprises the following specific implementation processes:
(S101) obtaining from the PET imaging mechanism and the instrument parameters Background noisePhoton pair projection data received by detectorThe Poisson model of PET image reconstruction can be expressed as:
g=Possion(Af+γ),
where Possion (x) represents a Poisson's distributed random vector expected to be x; representing a tracer profile to be solved;
(S102) establishing a smooth higher-order isotropic total variation (SHOITV) regularized PET image reconstruction model, i.e., a PET image reconstruction model with SHOITV regularization term, based on step (S101), the model being expressed as:
F is the fidelity term of the objective function, For SHOITV regularization terms, iota is a non-negative constraint indicating function term;
(S103) solving the proposed PET image reconstruction model;
an accelerated preconditioning proximity gradient algorithm APPGA with a high convergence order is introduced, the iterative format of which is shown below:
Due to Therefore, the above formula is written as:
wherein f k+1 refers to a PET reconstructed image generated after the APPGA is iterated in k steps;
representing Fidelity term F in PET image reconstruction model with respect to/> Is used for the gradient of (a),
I.e.
Representing canonical terms/>, in a PET image reconstruction modelConcerningIs used for the gradient of (a),
I.e.
Precondition matrixThe diagonal elements of (a) are positive numbers, and specifically are constructed as follows:
θ k is a momentum parameter, and θ k is defined as:
tk-1=a(k-1)ω+b,ω∈(0,1];
In order to secure the convergence and the convergence order of APPGA, the precondition matrix P and the setting of the momentum parameter θ k are required to satisfy the following conditions (i) and (ii):
(ii) wherein P max is the diagonal element of the precondition matrix Pmax, L is/> Lipschitz constant of (C);
(ii)ω∈(0,1], when ω is E (0, 1), a is equal to (0), ++ infinity a) is provided; when ω=1,/>
(S104) initializing f 0 and f 1, and directly iterating by using APPGA to obtain a reconstructed image f k+1 with k steps of iteration completed.
Preferably, the specific implementation process of the step (S102) is as follows:
The tracer distribution map F can be reconstructed according to g based on the PET image reconstruction model, and the Possion model reconstructed on the PET image is converted into a corresponding Possion maximum likelihood estimation optimization model, namely the following objective function F is minimized:
F(f)==<Af,1m>-<ln(Af+γ),g> (1)
Wherein < ·, · > represents the inner product of the vector, 1 m represents an m-dimensional vector with components 1, ln (·) represents a component-by-component natural logarithm operation;
introducing Gao Jiequan variation HOTV regularization, PET image model reconstruction is noted as follows:
Wherein the definition of the function F is shown as a formula (1); And/> Respectively a first-order total variation TV regularization term and a second-order total variation TV regularization term, wherein lambda 1 and lambda 2 are corresponding regularization parameters; iota is a non-negative constraint indicating function term used for ensuring that all components of the obtained f are greater than or equal to 0, and is specifically defined as:
With respect to the specific definition of the first and second order TV regularization terms, for the sake of symbol simplicity, the definition in the case of two-dimensional images is employed; recording device For the matrix Cronecker product, the first order differential matrixAnd second order differential matrixRespectively defined as:
Wherein/>
The Isotropic (TV) mode is chosen, in particular:
The smooth approximation function of l 2 norms is adopted, and is used for smoothing the regularization of the isotropic TV, and the smooth approximation function phi ε of l 2 norms II 2 is defined as:
The smooth higher-order isotropic total variation regularized (Smoothed HOITV, SHOITV) PET image reconstruction model is written as follows:
Wherein for the following
Simply writing the model (4) into a convex optimization model of three terms; recording device
Model (4) is equivalent to:
The step (S103) is replaced with:
APPGA, ROS-APPGA, which introduces a loosely ordered subset of ROS versions, the iterative format is shown below:
where M is the number of ordered subsets and, correspondingly, Is the gradient corresponding to the first subset, a l,gl,γl is the first subset of a, g, γ, l=1, 2, … M, respectively;
The precondition matrix P k,l is set to:
(iii) Gradient of fidelity term:
the gradients of F l in the iterative format of F and ROS-APPGA in equation (5) are:
(iv) SHOITV gradient of regularized term
First, the gradient of the smooth approximation function φ ε of l 2 norms ε IIII 2 is given:
then combine AndA gradient of SHOITV canonical terms is derived:
Wherein i=1, 2, …, d; r=0, 1; j=0, 1,2,3.
Preferably, the selection manner of the ordered subset data is as follows: the j-th column of the Subset matrix represents a Subset with an index j, and the components in the j-th column represent indexes corresponding to the sub-data selected by the j-th Subset; when the number of subsets is M, setting a system matrix A as M rows and d columns, setting each Subset (sub matrix) of the system matrix as M/M rows and d columns, and setting indexes of a j-th Subset as (j-1) M/m+1 to jm/M, namely, subset (: j) = (j-1) M/m+1:jm/M, j=1, 2 and … M, wherein at the moment, the subsystem matrix A j, the sub projection data g j and the sub background noise gamma j are respectively taken:
Aj=A(Subset(:,j),:);gj=g(Subset(:,j));γj=γ(Subset(:,j))。
Preferably, the order of subset at the inner layer iteration is confirmed as follows: adopting an array OrderedSUB to preset the sequence of subset iteration; if the number of subsets is less than or equal to 3, the i-th element of the array OrderedSUB correspondingly places the subset with index i; otherwise, alternately selecting the data with the subset corresponding to the farther projection angle.
The present invention is an innovative solution based on an accelerated preconditioning proximity gradient algorithm (ACCELERATED PRECONDITIONED PROXIMITY GRADIENT ALGORITHM, APPGA). The method combines the l 2 norm smoothing technology, the proximity gradient method, the generalized Nesterov momentum technology and the precondition technology, improves the quality of PET images under the condition of low radiation dose, and simultaneously accelerates the convergence speed of algorithms.
A smooth higher-order isotropic total variation regularization (SHOITV) PET image reconstruction model is adopted. By converting HOITV regularization terms into a smooth approximation form, smoothness and resolution of images are effectively maintained in the reconstruction process, reconstruction quality is improved, and a more efficient solving algorithm is convenient to design. The algorithm provided by the invention can not only improve the convergence order of the algorithm by effectively combining the generalized Nesterov momentum technology and the preconditioning technology, but also effectively improve the pathogenicity of the PET imaging problem by introducing the preconditioning matrix.
Compared with the traditional PET image method, the invention has the beneficial effects that: the PET image reconstruction method based on SHOITV regularization model and accelerating precondition adjacent gradient algorithm (APPGA) provided by the invention has the following advantages:
first, it can improve the quality of low radiation dose PET images, effectively maintain the smoothness and resolution of the images, and reduce the effects of noise and artifacts.
Second, it can significantly reduce the computational resources and time requirements for PET image reconstruction, thereby enabling a faster reconstruction process. APPGA has wide application prospect, and can be directly applied to solving optimization problems in other fields such as machine learning, image processing and the like.
Third, the method provided by the invention is expected to provide more accurate and high-quality PET images for disease diagnosis and treatment more efficiently in the field of medical imaging, and reduce radiation damage of PET scanning to patients.
Drawings
FIG. 1 is a schematic diagram of an embodiment of the process performed APPGA.
FIG. 2 is a schematic diagram of a specific process for performing ROS-APPGA.
Fig. 3 (a) is a brain phantom: high quality clinical PET brain map; 3 (b) is a uniform phantom: a uniform thermal sphere schematic with a uniform background and six different radii.
FIG. 4 (a) is a graph of Normalized Objective Function Values (NOFV) of three algorithms ROS-PKMA, ROS-PPGA, and ROS-APPGA at low dose, as a function of the number of iterative steps of the algorithm; 4 (b) is a plot of Normalized Root Mean Square Error (NRMSE) of three algorithms as a function of algorithm iteration steps at low dose.
Fig. 5 is a graph of brain phantom reconstructions obtained after 5, 10 and 50 iterations of the three algorithms ROS-PKMA, ROS-PPGA and ROS-APPGA (ω=1), respectively.
FIG. 6 (a) is a graph of Normalized Relative Contrast (NRC) of the maximum thermal sphere reconstructed by the three algorithms as a function of iteration steps for reconstruction of a homogeneous phantom at low dose; 6 (b) a graph of Normalized Relative Contrast (NRC) of the minimum thermal sphere reconstructed by the three algorithms as a function of iteration steps.
FIG. 7 (a) is a centerline profile (CLP) of a reconstructed image obtained after 5 iterations of three algorithms; and 7 (b) performing 10 steps of iteration on the three algorithms to obtain a Central Line Profile (CLP) of the reconstructed image.
FIG. 8 is a graph of homogeneous phantom reconstructions obtained by 5, 10 and 50 iterations of the three algorithms ROS-PKMA, ROS-PPGA and ROS-APPGA, respectively.
Detailed Description
The invention is further described below with reference to the drawings and detailed description.
1. Smooth high-order isotropy total variation regularized PET image reconstruction model
PET imaging is a functional and metabolically active medical imaging technique for observing internal organs and tissues of the human body. The imaging technology firstly injects radioactive tracer into a human body, positron generated after decay of the tracer generates annihilation reaction in the human body to generate a large number of photon pairs, and then a tracer distribution diagram in the human body is reconstructed according to photon pair information received by a detector. In PET imaging, the relationship between the tracer distribution and the photon information received by the detector can be characterized by a system matrix. The system matrix comprises complex physical processes such as positron random walk, photon projection, photon attenuation, detector sensitivity characterization and the like. Recording deviceIs a PET system matrix,For the profile of the tracer (vectorisation)Is background noise,Projection data is pairs of photons received for the detector. The basic model of the PET image can be expressed as:
g=Possion(Af+γ),
where Possion (x) represents a poisson distribution random vector expected to be x. The objective of PET image reconstruction is to reconstruct the tracer distribution profile f from the projection data g, which is a classical inverse problem. For this model, the traditional solution is to transform the model into a corresponding Possion maximum likelihood estimation optimization model, i.e. to minimize the following objective function F:
F(f):=<Af,1m>-<ln(Af+γ),g>, (1)
Where < ·, · > represents the inner product of the vector, 1 m represents an m-dimensional vector with elements 1, and ln (·) represents an element-wise natural logarithm operation. However, the minimization problem of directly solving F is an ill-posed problem that can lead to severe noise in the reconstructed image (especially at low doses) due to overfitting, and to overcome this problem, it is often necessary to introduce a regularization term for adding a priori information and improving the quality and stability of the reconstructed image. Total Variation (TV) regularization has met with great success in image denoising and has been successfully applied to medical imaging problems. However, conventional first-order TV regularization, while effective in removing noise, results in reconstructed images with characteristics of the slice constants, which are not suitable for PET images. Because the PET reconstructed image is a profile of the tracer with a certain smoothness, the first order TV regularization of the slicing constant characteristics can bring serious step artifacts to the reconstructed image.
To solve this problem, the ladder artifact caused by the first-order TV is suppressed by introducing a high-order TV regularization term, while the effects of effectively removing noise and retaining details are played. Specifically, the high-order Quan Bianfen (Higher-Order Total Variation, HOTV) regularized PET image reconstruction model can be written in the form:
Wherein the definition of the function F is shown as a formula (1); And/> Respectively a first-order TV regularization term and a second-order TV regularization term, wherein lambda 1 and lambda 2 are corresponding regularization parameters; iota is a non-negative constraint indicating function term used for ensuring that all components of the obtained f are greater than or equal to 0, and is specifically defined as:
With respect to specific definitions of the first and second order TV regularization terms, for the sake of symbolic simplicity, a definition in the case of a two-dimensional image is given here.
Recording deviceFor the matrix Cronecker product, the first order differential matrixSum second order differential matrixRespectively defined as:
Wherein/>
WhileThere are two options corresponding to anisotropic (Anisotropic) TV and Isotropic (Isotropic) TV, respectively. For anisotropic TV,AndAre all/>, in their respective dimensionsNorms, i.e.
For an isotropic TV,
For PET image models (2) with HOTV regularization terms, whether anisotropic or isotropic TV is employed, the regularization terms are not conductive.
Although many algorithms are currently available for solving such PET image reconstruction models with non-guided regularization terms, it remains a significant challenge to propose a numerical algorithm with a high convergence order to solve the model. To overcome the problem of non-impossibility, the existing work has focused on smoothing the l 1 norms in the anisotropic TV regularization term. However, in PET image models, using anisotropic TV as regularization term results in a round object being reconstructed into a square structure, which is detrimental to the accurate reconstruction of lesions in PET images, whereas isotropic TV does not have this problem.
Therefore, in the PET image reconstruction problem, the selection of the isotropic TV regularization method is more capable of maintaining shape consistency and retaining important detail information. Notably, the present embodiment originally proposesA smooth approximation function of the norm and used to smooth the isotropic TV regularization. Specifically, we defineThe smooth approximation function φ ε for norm II 2 is:
The smooth high-order isotropic total variation regularization (Smoothed HOITV, SHOITV) PET image reconstruction model proposed by the present invention can be written as follows:
Wherein for the following
For ease of solution, the model (4) is further simply written as a convex optimization model of three terms. Recording device
Model (4) is equivalent to:
2. solving the proposed PET image reconstruction model
(1) Acceleration precondition proximity gradient algorithm APPGA
Before introducing the accelerating preconditioning proximity gradient algorithm APPGA, the definition of the proximity operator is reviewed.
For the convex function ψ: it is at/> The proximity operator at a point is defined as:
The conventional adjacent gradient algorithm for solving the model (5) is:
For large-scale and pathological problems such as PET image reconstruction, the traditional adjacent gradient algorithm (6) solves very slow and low-order convergence without adding any preconditioning and acceleration techniques.
Example 1
In order to solve this problem, the present embodiment first proposes an Accelerating Preconditioned Proximity Gradient Algorithm (APPGA) with high-order convergence, which skillfully combines acceleration techniques such as proximity gradient method, generalized Nesterov momentum technique, and preconditioning technique. The generalized Nesterov momentum technology is based on a generalized version (Y.Lin,S.Li and Y.Zhang.Convergence rate analysis of accelerated forward-backward algorithm with generalized Nesterov momentum Scheme.International Journal of Numerical Analysis and Modeling.2023.20(4):518-537), proposed by the traditional Nesterov momentum technology, and has higher degree of freedom and flexibility in selecting momentum parameters. Although this study gives convergence and convergence order results of the adjacent gradient algorithm with generalized Nesterov momentum acceleration, the introduction of a preconditioned acceleration matrix is not considered nor is the algorithm applied to PET medical imaging problems. The introduction of the preconditioned acceleration matrix brings little difficulty to the convergence order analysis of APPGA.
The present invention overcomes these difficulties and theoretically demonstrates the convergence order of APPGA.
The introduction of the generalized Nesterov momentum acceleration technology can enable the algorithm to have higher convergence speed and higher theoretical convergence order in the iterative process; the introduction of the precondition technology can effectively improve the pathogenicity of the original problem, and the gradient is properly adjusted in each iteration, so that the convergence speed and stability of the algorithm are further improved. In this embodiment, the proposed iterative format of the accelerating preconditioning proximity gradient algorithm with high convergence order is as follows:
Wherein the method comprises the steps of For a precondition matrix that is positive for both diagonal elements, the definition of the momentum parameter θ k is:
tk-1=a(k-1)ω+b,ω∈(0,1].
Due to So (7) can be written as:
Wherein the preconditioning matrix P can be set to a classical EM preconditioning matrix:
In order to perform APPGA, two initial vectors f 0 and f 1 need to be given in advance, and f k+1 is calculated.
APPGA high-order convergence
In order to obtain APPGA convergence and convergence order, the parameters in the algorithm are required to satisfy the following conditions:
(i) wherein P max is the diagonal element of the precondition matrix Pmax, L is/> Lipschitz constant of (C);
(ii)ω∈(0,1], When ω is E (0, 1), a is equal to (0), ++ infinity a) is provided; and when ω=1,/>
In the present invention, it has been theoretically demonstrated that if the algorithm parameters in APPGA satisfy the conditions (i) and (ii) above, the objective function can be converged on for any of the initial vectors f 0 and f 1, APPGAAnd the algorithm has the following convergence order: first, according to the convergence order of the objective functionSecond, the convergence order according to the difference II f k+1-fk‖2 between the front and back steps is
That is, it is proved that the smooth high-order isotropic total variation regularized PET image reconstruction model can be solved with a higher convergence order by using APPGA algorithm. This result is of great importance for both theoretical analysis and practical application of the algorithm.
Example 2
To address the problem of PET image reconstruction, a further embodiment proposes to introduce a relaxed ordered subset APPGA.
PET image reconstruction algorithms employed in clinically used PET instruments are typically further accelerated using ordered subsets. The core idea is to divide the system matrix and projection data into a plurality of subsets according to a certain sequence, and then each time the algorithm iteration is carried out, the inner layer iteration and the outer layer iteration are carried out. Each step of inner layer iteration only adopts data of a subset to iterate, and after all subset data completes one round of iteration, the inner layer iteration is equivalent to one outer layer iteration, and the subset data inner layer iteration is continuously circulated in each outer layer iteration by adopting a preset sequence. The ordered subset strategy can greatly accelerate the convergence rate of the algorithm on the problem of PET image reconstruction, because the data corresponding to the adjacent projection angles of PET are relatively close, the similar data are divided into a plurality of subsets, and the step parameters of the iteration process are controlled, so that each subset data iteration can have the effect of approaching the complete data iteration. To further better apply APPGA to clinic, the present invention also proposes a relaxed ordered subset (Relaxed Ordered Subsets, ROS) version of APPGA (ROS-APPGA), whose iterative format is as follows:
where M is the number of ordered subsets. In response to this, the control unit, Is the gradient corresponding to the first subset, a l,gl,γl is the first subset of a, g, γ, l=1, 2, … M, respectively. And precondition matrix P k,l is set to:
I.e., the specific implementation of introducing the relaxed ordered subset APPGA:
(v) Gradient of fidelity term:
the gradient of F in formula (5) and F l in formula (9) is:
(vi) SHOITV gradient of regularized term
First, the gradient of the smooth approximation function φ ε of l 2 norms ε IIII 2 is given:
then combine AndA gradient of SHOITV canonical terms is derived:
Wherein i=1, 2, …, d; r=0, 1; j=0, 1,2,3.
3. Initialization of reconstructed image f
After the PET image reconstruction model is proposed by the solution of example 1 or example 2, a PET reconstructed image f can be obtained.
In PET images, only objects within the circular Field of View (FOV) are reconstructed. In order to better and faster provide the initial vector f initial for the iterative algorithm, the invention calculates the corresponding real average photon count (TMC) of the reconstructed image according to the projection data obtained by the detector, the projection angle number of the PET instrument, photon attenuation information, background noise and other information. The initial vector f initial is then set to a vectorized uniform disk that is the same size as the FOV, its components in the FOV are set to TMC values, and the components outside the FOV are set to 0.
The specific process and operation of implementation APPGA are as follows:
① Input: inputting a system matrix A, calculating A T1m, calculating Λ=A T1m, and Λ (Λ is less than or equal to 0) =1; projection data g; a background noise term γ; regularization parameters λ 1 and λ 2; a threshold epsilon; a step size parameter beta; momentum order parameter ω e (0, 1), let a=1/8, b=1, maximum iteration step number MaxIter.
② Initializing: let t 1=1;f0=f1=finitial.
③ And (3) iteration loop:
for k=1:MaxIter
tk=akω+b
θk=(tk-1-1)/tk
end for
④ And (3) outputting: f k+1.
The specific process and operation for achieving ROS-APPGA are as follows:
① Input: inputting a system matrix A, calculating A T1m, calculating Λ=A T1m, and Λ (Λ is less than or equal to 0) =1; projection data g; a background noise term γ; regularization parameters λ 1 and λ 2; a threshold epsilon; momentum order parameter omega epsilon (0, 1), let a=1/8, b=1, ordered Subset number M, maximum iteration step number MaxIter, M dimension vector OrderedSUB, wherein each element is used to store the index of each Subset, subset matrix Subset.
Subset data selection: the j-th column of the Subset matrix represents a Subset with an index j, and the components in the j-th column represent indexes corresponding to the sub-data selected by the j-th Subset. For example, when the number of subsets is M, since the system matrix a is M rows and d columns, each Subset (sub-matrix) of the system matrix is M/M rows and d columns, the index of the j-th Subset is (j-1) M/m+1 to jm/M, that is, subset (: j) = (j-1) M/m+1:jm/M, j=1, 2, … M. At this time, the sub-system matrix a j, the sub-projection data g j, and the sub-background noise γ j are respectively taken:
Aj=A(Subset(:,j),:);gj=g(Subset(:,j));γj=γ(Subset(:,j)).
② Order of subsets in inner layer iteration: because the subset data of the adjacent projection angles of the PET are relatively close, in order to enable the subset alternating iteration process to have a better acceleration effect, adjacent sub-iterations should select subset data far away for iteration. To achieve this, we use array OrderedSUB to preset the order of subset iterations. If the number of subsets is less than or equal to 3, the i-th element of the array OrderedSUB correspondingly places the subset with index i; otherwise, the data with the subset corresponding to the farther projection angle is selected alternately, for example, after the first subset is selected in sequence, the second subset can select the intermediate value of the number of the subsets farther from the first subset. The method comprises the following steps:
③ Initializing: let t 1,0=a+b;f1,1=f1,0=finitial;
④ And (3) iteration loop:
⑤ And (3) outputting: f k+1,1.
The beneficial effects of the invention are: the method can remarkably reduce the calculation time required by the reconstruction process while improving the quality and accuracy of the low-dose PET image reconstruction, overcomes the problem that the speed and the image quality in the prior PET image field are difficult to reach balance to a certain extent, and breaks the limitation of the prior art.
The following describes the embodiments of the present invention further with reference to the drawings and the examples.
(1) And realizing a 2D simulation experiment of PET image reconstruction based on Matlab. The simulation platform adopted by the simulation experiment is developed by the inventor 2015-2016 during the visit period of the Stonekote cancer center of New York souvenir in U.S. The invention uses ROS-APPGA to reconstruct PET image under low radiation dose condition, and uses two phantoms (phantoms) which are respectively Brain phantoms (Brain phantoms) and Uniform phantoms (uniformity phantoms), and the original diagrams of the two phantoms are respectively shown in (a) and (b) in figure 3.
(2) Experiments on brain phantoms are first shown. In fig. 4, the present example demonstrates ROS-APPGA in comparison to the convergence results of the existing ROS-PKMA and ROS-PPGA algorithms that do not drive the amount acceleration, using the number of iterations as the abscissa, normalized objective function values (Normalized Objective Function Value, NOFV) and normalized root mean square error (Normalized Root Mean Square Error, NRMSE) as the ordinate, and the observation that ROS-APPGA (ω=1) exhibits a faster convergence rate than ROS-PPGA and ROS-PKMA; and as the value of ω increases, the convergence rate of ROS-APPGA increases. Fig. 5 shows a comparison of experimental results obtained by performing 5, 10 and 50 iterations of the three algorithms, respectively, it can be observed that ROS-APPGA can obtain a higher quality reconstructed image with a smaller number of iteration steps, and the 50-iteration reconstructed image exhibits more prominent details, higher accuracy and fidelity.
(3) Secondly, experiments of uniform body molds are carried out. FIG. 6 shows the relationship between the number of iterations of the maximum and minimum thermal spheres in a uniform phantom and normalized relative contrast (Normalized Relative Contrast, NRC); FIG. 7 shows a centerline Profile (CLP) of the reconstructed image after 5 and 10 iterations; in addition, fig. 8 shows a comparison of experimental results obtained by performing 5, 10 and 50 iterations of the three algorithms, respectively. These experimental results clearly demonstrate that ROS-APPGA exhibits superior performance compared to the other two algorithms. Wherein, the reconstructed image of ROS-APPGA (omega=1) obtained by 50 steps of iteration clearly shows 6 uniform hot spheres, which are obviously superior to ROS-PPGA and ROS-PKMA.
Based on the above examples, it can be seen that the image generated by ROS-APPGA designed by the present invention exhibits more prominent detail, higher accuracy and fidelity than ROS-PKMA and ROS-PPGA, exceeding the performance levels of ROS-PKMA and ROS-PPGA, at the same number of iteration steps. Furthermore, as shown in fig. 3 and 6, the reconstruction results of ROS-APPGA exhibit a higher imaging quality, which means a more efficient workflow and more timely diagnostic results for the clinician.
The embodiments of the present invention described above do not limit the scope of the present invention. Any modifications, equivalent substitutions and improvements made within the spirit principles of the present invention should be included in the scope of the claims of the present invention.
Claims (6)
1. A low-radiotracer high-quality PET image reconstruction method based on smooth high-order isotropy total variation regularization and acceleration precondition adjacent gradient algorithm is characterized by comprising the following specific implementation processes:
(S101) obtaining from the PET imaging mechanism and the instrument parameters Background noisePhoton pair projection data received by detectorThe Poisson model of PET image reconstruction can be expressed as:
g=Possion(Af+γ),
where Possion (x) represents a Poisson's distributed random vector expected to be x; representing a tracer profile to be solved;
(S102) establishing a smooth higher-order isotropic total variation (SHOITV) regularized PET image reconstruction model, i.e., a PET image reconstruction model with SHOITV regularization term, based on step (S101), the model being expressed as:
F is the fidelity term of the objective function, For SHOITV regularization terms, iota is a non-negative constraint indicating function term;
(S103) solving the proposed PET image reconstruction model;
an accelerated preconditioning proximity gradient algorithm APPGA with a high convergence order is introduced, the iterative format of which is shown below:
Due to Therefore, the above formula is written as:
wherein f k+1 refers to a PET reconstructed image generated after the APPGA is iterated in k steps;
representing Fidelity term F in PET image reconstruction model with respect to/> Is used for the gradient of (a),
I.e.
Representing canonical terms/>, in a PET image reconstruction modelConcerningIs used for the gradient of (a),
I.e.
Precondition matrixThe diagonal elements of (a) are positive numbers, and specifically are constructed as follows:
θ k is a momentum parameter, and θ k is defined as:
In order to secure the convergence and the convergence order of APPGA, the precondition matrix P and the setting of the momentum parameter θ k are required to satisfy the following conditions (i) and (ii):
(i) wherein P max is the diagonal element of the precondition matrix Pmax, L is/> Lipschitz constant of (C);
(ii)ω∈(0,1], when ω is E (0, 1), a is equal to (0), ++ infinity a) is provided; when ω=1, the number of the cells,
(S104) initializing f 0 and f 1, and directly iterating by using APPGA to obtain a reconstructed image f k+1 with k steps of iteration completed.
2. The method according to claim 1, wherein the specific implementation procedure of the step (S102) is as follows:
The tracer distribution map F can be reconstructed according to g based on the PET image reconstruction model, and the Possion model reconstructed on the PET image is converted into a corresponding Possion maximum likelihood estimation optimization model, namely the following objective function F is minimized:
F(f):=<Af,1m>-<ln(Af+γ),g>, (1)
Wherein < ·, · > represents the inner product of the vector, 1 m represents an m-dimensional vector with components 1, ln (·) represents a component-by-component natural logarithm operation;
introducing Gao Jiequan variation HOTV regularization, PET image model reconstruction is noted as follows:
Wherein the definition of the function F is shown as a formula (1); And/> Respectively a first-order total variation TV regularization term and a second-order total variation TV regularization term, wherein lambda 1 and lambda 2 are corresponding regularization parameters; iota is a non-negative constraint indicating function term used for ensuring that all components of the obtained f are greater than or equal to 0, and is specifically defined as:
With respect to the specific definition of the first and second order TV regularization terms, for the sake of symbol simplicity, the definition in the case of two-dimensional images is employed; recording device For the matrix Cronecker product, the first order differential matrixAnd second order differential matrixRespectively defined as:
Wherein/>
The Isotropic (TV) mode is chosen, in particular:
The smooth approximation function of l 2 norms is adopted, and is used for smoothing the regularization of the isotropic TV, and the smooth approximation function phi ε of l 2 norms II 2 is defined as:
The smooth higher-order isotropic total variation regularized (Smoothed HOITV, SHOITV) PET image reconstruction model is written as follows:
Wherein for the following
Simply writing the model (4) into a convex optimization model of three terms; recording device
Model (4) is equivalent to:
3. the method according to claim 1, wherein the step (S103) is replaced by:
APPGA, ROS-APPGA, which introduces a loosely ordered subset of ROS versions, the iterative format is shown below:
where M is the number of ordered subsets and, correspondingly, Is the gradient corresponding to the first subset, a l,gl,γl is the first subset of a, g, γ, l=1, 2, … M, respectively;
The precondition matrix P k,l is set to:
(i) Gradient of fidelity term:
the gradients of F l in the iterative format of F and ROS-APPGA in equation (5) are:
(ii) SHOITV gradient of regularized term
First, the gradient of the smooth approximation function φ ε of l 2 norms ε IIII 2 is given:
then combine AndA gradient of SHOITV canonical terms is derived:
Wherein i=1, 2, …, d; r=0, 1; j=0, 1,2,3.
4. A method according to claim 3, wherein the ordered subset data is selected in the following manner: the j-th column of the Subset matrix represents a Subset with an index j, and the components in the j-th column represent indexes corresponding to the sub-data selected by the j-th Subset; when the number of subsets is M, setting a system matrix A as M rows and d columns, setting each Subset (sub matrix) of the system matrix as M/M rows and d columns, and setting indexes of a j-th Subset as (j-1) M/m+1 to jm/M, namely, subset (: j) = (j-1) M/m+1:jm/M, j=1, 2 and … M, wherein at the moment, the subsystem matrix A j, the sub projection data g j and the sub background noise gamma j are respectively taken:
Aj=A(Subset(:,j),:);gj=g(Subest(:,j));γj=γ(Subset(:,j))。
5. the method of claim 4, wherein the order validation of subsets at the inner layer iteration is as follows: adopting an array OrderedSUB to preset the sequence of subset iteration; if the number of subsets is less than or equal to 3, the i-th element of the array OrderedSUB correspondingly places the subset with index i; otherwise, alternately selecting the data with the subset corresponding to the farther projection angle.
6. The method according to any one of claims 1-5, wherein f 0=f1=finitial is set up, wherein the reconstructed image f initial is initialized in such a way that: calculating real average photon pair count TMC corresponding to the PET reconstructed image based on projection data, the projection angle number of the PET instrument, photon attenuation information and background noise information; the initial vector f initial is then set to a vectorized uniform disk that is the same size as the circular field of view FOV, its components in the circular field of view FOV are set to TMC values, and the components outside the circular field of view FOV are set to 0.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202311817363.XA CN117830373B (en) | 2023-12-26 | 2023-12-26 | PET image reconstruction method based on acceleration precondition adjacent gradient algorithm |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202311817363.XA CN117830373B (en) | 2023-12-26 | 2023-12-26 | PET image reconstruction method based on acceleration precondition adjacent gradient algorithm |
Publications (2)
Publication Number | Publication Date |
---|---|
CN117830373A CN117830373A (en) | 2024-04-05 |
CN117830373B true CN117830373B (en) | 2024-06-14 |
Family
ID=90522374
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202311817363.XA Active CN117830373B (en) | 2023-12-26 | 2023-12-26 | PET image reconstruction method based on acceleration precondition adjacent gradient algorithm |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN117830373B (en) |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103998202A (en) * | 2011-10-21 | 2014-08-20 | 考泰克公司 | Non-symmetric multiple layer injection molded products and methods |
Family Cites Families (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2810245A4 (en) * | 2012-02-03 | 2015-12-16 | Univ New York State Res Found | Methods and systems for inverse problem reconstruction and application to ect reconstruction |
CN108550172B (en) * | 2018-03-07 | 2020-05-19 | 浙江大学 | PET image reconstruction method based on non-local characteristics and total variation joint constraint |
US11087508B2 (en) * | 2018-11-30 | 2021-08-10 | Canon Medical Systems Corporation | Method and apparatus for acceleration of iterative reconstruction of a computed tomography image |
CN112001981B (en) * | 2020-08-28 | 2023-12-12 | 南京医科大学 | Compressed sampling MR image reconstruction method based on generalized Fresnel-fast-gradient acceleration conjugate gradient algorithm |
CN115018950A (en) * | 2022-07-15 | 2022-09-06 | 重庆师范大学 | Rapid finite angle fan-beam CT image reconstruction method based on preconditioned matrix |
CN115880387A (en) * | 2022-11-14 | 2023-03-31 | 郑州轻工业大学 | Sparse angle sampling CT iterative reconstruction method based on discrete NLHTV regular term |
-
2023
- 2023-12-26 CN CN202311817363.XA patent/CN117830373B/en active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103998202A (en) * | 2011-10-21 | 2014-08-20 | 考泰克公司 | Non-symmetric multiple layer injection molded products and methods |
Non-Patent Citations (1)
Title |
---|
基于PICCS的非局部TGV能谱CT重建算法;雷蕾;孔慧华;;数学的实践与认识;20200308(05);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN117830373A (en) | 2024-04-05 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Cai et al. | Cine cone beam CT reconstruction using low-rank matrix factorization: algorithm and a proof-of-principle study | |
CN111627082B (en) | PET image reconstruction method based on filtering back projection algorithm and neural network | |
CN103514629B (en) | Method and apparatus for iterative reconstruction | |
CN110148215B (en) | Four-dimensional magnetic resonance image reconstruction method based on smooth constraint and local low-rank constraint model | |
Jia et al. | GPU-based fast low-dose cone beam CT reconstruction via total variation | |
CN104700438B (en) | Image rebuilding method and device | |
Shao et al. | A learned reconstruction network for SPECT imaging | |
Zhang et al. | Accurate and robust sparse‐view angle CT image reconstruction using deep learning and prior image constrained compressed sensing (DL‐PICCS) | |
Zhang et al. | PET image reconstruction using a cascading back-projection neural network | |
CN114067013A (en) | System and method for reprojection and backprojection via a homographic resampling transform | |
Shao et al. | SPECTnet: a deep learning neural network for SPECT image reconstruction | |
Li et al. | Sparse CT reconstruction based on multi-direction anisotropic total variation (MDATV) | |
Li et al. | Training end-to-end unrolled iterative neural networks for SPECT image reconstruction | |
CN107146263B (en) | A kind of dynamic PET images method for reconstructing based on the constraint of tensor dictionary | |
Sang et al. | A conditional registration network for continuous 4D respiratory motion synthesis | |
Moriakov et al. | Lire: Learned invertible reconstruction for cone beam ct | |
CN117830373B (en) | PET image reconstruction method based on acceleration precondition adjacent gradient algorithm | |
Jiao et al. | Fast PET reconstruction using multi-scale fully convolutional neural networks | |
Yang et al. | Cycle-consistent learning-based hybrid iterative reconstruction for whole-body PET imaging | |
Cheng et al. | Maximum likelihood activity and attenuation estimation using both emission and transmission data with application to utilization of Lu‐176 background radiation in TOF PET | |
Germer et al. | Limited-angle tomography reconstruction via deep end-to-end learning on synthetic data | |
Chen et al. | Efficient and robust 3D CT image reconstruction based on total generalized variation regularization using the alternating direction method | |
Rao et al. | A novel supervised learning method to generate CT images for attenuation correction in delayed pet scans | |
Zhong et al. | 3D‐2D Deformable Image Registration Using Feature‐Based Nonuniform Meshes | |
KR102329938B1 (en) | Method for processing conebeam computed tomography image using artificial neural network and apparatus therefor |
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 | ||
GR01 | Patent grant | ||
GR01 | Patent grant |