CN109636869B - Dynamic PET Image Reconstruction Method Based on Non-local Total Variation and Low-Rank Constraints - Google Patents

Abstract

The invention discloses a dynamic PET image reconstruction method based on non-local total variation and low-rank constraint, which utilizes the segmentation smooth characteristic and the space-time correlation of a PET image and introduces the low-rank constraint and the non-local total variation constraint at the same time, realizes the reconstruction of the dynamic PET image, can remove the step effect and keep fine details, and is beneficial to improving early focus detection. The method optimizes the space-time correlation of the solution through the constraint of low rank and sparsity, namely, the background components and the detail components of the image are reconstructed through a combined model; meanwhile, in the reconstruction of the dynamic PET image, in order to fully consider the structural smoothness characteristic of image data, the invention introduces non-local total variation constraint, and recovers more image details and removes the step effect by utilizing the image redundancy characteristic. Compared with the prior art, the method can provide more accurate reconstructed images to improve lesion detection, and has better robustness to noise.

Description

Dynamic PET Image Reconstruction Method Based on Non-local Total Variation and Low-Rank Constraints

technical field

The invention belongs to the technical field of PET imaging, and in particular relates to a dynamic PET image reconstruction method based on non-local total variation and low rank constraints.

Background technique

Dynamic positron emission tomography (DPET) enables monitoring of the in vivo spatiotemporal distribution of radiolabeled tracers, which has the potential to improve early detection, cancer characterization, and assessment of treatment response. DPET uses a system of detectors placed around the object in order to obtain many different angular views over a series of time frames. Using these projection data, a map of the radiotracer concentration can be reconstructed, a dynamic PET image reconstruction, which passes through the time series Images can better evaluate the dynamic process of tracer uptake and metabolism, and have important application value in scientific research and clinical applications.

PET image reconstruction is an ill-posed inverse problem, and the standard practice is to use a regular term to constrain the solution of the problem, making the inverse problem well-posed. In addition to the traditional ML-EM and other algorithms, the first type of algorithm is the spatial smoothness constraint, one of which is the maximum a posteriori (MAP) algorithm, and the penalty term designed to preserve the smoothness of the region and the sharp change of the edge has been studied in PET reconstruction. Important point; another approach is to use total variation constraints (TV) to improve the structural smoothness of PET images, however TV-based models assume that each image pixel always has edges and gradient diffusion directions, which may lead to stair-step effects . The second category of algorithms is to use both temporal and spatial information to improve the quality of dynamic PET image reconstruction. Temporal constraints are usually used to improve the additional robustness of spatial solutions, including spline spline models, tracer dynamics, wavelet transforms, etc. However, Tracer dynamics methods assume that all voxels are well modeled by the same set of model dynamics, which may not be the case in practice. In addition, the wavelet transform method still leaves room for choice in choosing the E-spline wavelet parameter vector and the appropriate wavelet coefficients, in terms of the problem of early lesion detection using PET images, although the supplementation obtained from CT and MRI images during reconstruction Information was incorporated into the canonical constraints, but it did not provide information on organ and lesion metabolism, thus lacking sufficient information to guide the reconstruction of dysfunctional regions. Several studies have shown that using anatomical boundaries without accurate lesion contours does not improve lesion detection or quantification tasks.

SUMMARY OF THE INVENTION

In view of the above, the present invention provides a dynamic PET image reconstruction method based on non-local total variation and low-rank constraints, which utilizes the piecewise smoothness and spatiotemporal correlation of PET images, while introducing low-rank constraints and non-local total Variational constraints enable dynamic PET image reconstruction, which can remove staircase effects and maintain fine details, which is beneficial for improving early lesion detection.

A dynamic PET image reconstruction method based on non-local total variation and low-rank constraints, comprising the following steps:

(1) Use a detector to detect the biological tissue injected with radiopharmaceuticals, dynamically collect the coincidence count vectors corresponding to each moment, and combine these coincidence count vectors into a coincidence count matrix Y;

(2) Combining the dynamic PET image sequence into a PET concentration distribution matrix X, and establishing a PET measurement equation according to the PET imaging principle;

(3) By introducing a low-rank constraint to the PET measurement equation, a dynamic PET image reconstruction model M1 based on the low-rank constraint is obtained;

(4) The dynamic PET image reconstruction model M2 based on the non-local total variation constraint is further obtained by performing non-local total variation constraints on the low-rank and sparse parts of each frame of images;

(5) The objective function of dynamic PET reconstruction is obtained by combining the dynamic PET image reconstruction models M1 and M2 as follows:

s.t.L+S=X

where: || ||* represents the nuclear norm, || || ₁ represents the 1-norm, L is the low-rank part that contains the periodically varying radioactive concentration in the image background, S is the sparse part that contains the radioactivity with different metabolic rates The radioactive concentration of heterogeneous tissue, J _NLTV (L) and J _NLTV (S) are the results of the low-rank part L and the sparse part S after non-local total variation, respectively, λ, μ, α _L and α _S are weights coefficient, Ψ(Y|X) is the likelihood function about X and Y;

(6) After the above objective function is optimized and solved, the PET concentration distribution matrix X is obtained, thereby restoring the dynamic PET image sequence.

Further, the expression of the PET measurement equation is as follows:

Y=DX+R+S

Where: D is the system matrix, R and S are the measurement noise matrices reflecting random events and scattering events, respectively.

Further, the expression of the dynamic PET image reconstruction model M1 is as follows:

M1=||L|| _* +λ||S|| ₁ +μΨ(Y|X)

s.t. L+S=X

Further, the expression of the dynamic PET image reconstruction model M2 is as follows:

M2=α _L J _NLTV (L)+α _S J _NLTV (S)+μΨ(Y|X)

s.t. L+S=X

Further, the expression of the likelihood function Ψ(Y|X) is as follows:

Where: d _ij is the element value of the i-th row and the j-th column in the system matrix D, y _im is the element value of the i-th row and the m-th column in the coincidence count matrix Y, and x _j m is the j-th row of the PET concentration distribution matrix X. Column element value, rim is the element value of the i-th row and m-th column in the measurement noise matrix R, sim is the i-th row and m-th column element value of the measurement noise matrix S, i, j, and m are all natural numbers and _1≤i≤ I, 1≤j≤J, 1≤m≤M, I is the dimension conforming to the count vector, J is the row number of the PET concentration distribution matrix X, that is, the number of pixels of the PET image, M is the column of the PET concentration distribution matrix X The number is the sampling time length.

Further, in the step (6), ADMM (Alternating Direction Method of Multipliers, Alternating Direction Multipliers) algorithm is used to optimize the solution of the objective function; wherein, the low-rank constraint problem utilizes the singular value threshold algorithm and the soft shrinking algorithm for iterative optimization. To solve, the likelihood function problem is solved by iterative optimization using the EM (ExpectationMaximizationAlgorithm, expectation maximization) algorithm, and the non-local total variation constraint problem is solved iteratively by the gradient descent method.

Different from the existing methods that use temporal prior and spatial prior as two different constraints, the present invention optimizes the spatial-temporal correlation of the solution through the constraints of low rank and sparseness, that is, reconstructs the image through a joint model At the same time, in dynamic PET image reconstruction, in order to fully consider the structural smoothness of image data, the present invention introduces non-local total variation (NLTV) constraints, and uses image redundancy to restore more images. detail and remove the staircase effect. Compared with TV, NLTV integrates local and non-local correlation of image matrix, the present invention can give more accurate reconstructed images to improve lesion detection, and has better robustness to noise.

The present invention introduces non-local total variation constraints into the low-rank matrix restoration analysis framework to ensure the structural smoothness and clear boundaries of regions of interest (ROIs) in PET images; the low-rank constraints of image sequences are achieved by inherent averaging within the tissue It is also an innovation of the present invention to introduce non-local total variation to improve the spatial resolution of PET images; low-rank and sparse matrix decomposition can provide random noise components for NLTV constraints, with the confirmation of random noise information, NLTV constraints It can provide enhanced clean PET images, and in turn help the decomposition of low-rank and sparse matrices, resulting in a more realistic reconstruction result. Combined with the performance of the present invention in the simulated data experiment, by comparing with the results of the ML-EM algorithm and the LRTV algorithm (based on low rank and total variational constraints), the present invention can obtain more accurate reconstruction results, which is useful for improving early lesions. Detecting and evaluating the dynamic process of tracer uptake and metabolism has important practical application value.

Description of drawings

FIG. 1 is a schematic flowchart of a dynamic PET image reconstruction method according to the present invention.

Figure 2(a) is a template image of Monte Carlo simulation Zubal chest data.

Figure 2(b) is a template image of Hoffman's brain data.

Figure 3(a) is the real image of the eighth frame of Zubal's chest data.

Figure 3(b) is the 8th frame PET image result of reconstruction of Monte Carlo simulated Zubal thoracic data by ML-EM method under the data count rate of 1×10 ⁷ .

Figure 3(c) shows the results of the 8th frame PET image reconstructed by the LRTV method on the Monte Carlo simulated Zubal thoracic data at a data count rate of 1×10 ⁷ .

Fig. 3(d) is the result of the eighth frame of PET image reconstructed from the Monte Carlo simulated Zubal thoracic cavity data using the method of the present invention at a data count rate of 1×10 ⁷ .

Fig. 4(a) is the enlarged image result of the part framed in Fig. 3(a).

Fig. 4(b) is the enlarged image result of the part framed in Fig. 3(b).

Fig. 4(c) is the enlarged image result of the part framed in Fig. 3(c).

Fig. 4(d) is the enlarged image result of the part framed in Fig. 3(d).

Figure 5(a) is a line graph of the deviation of each frame of Zubal thoracic data image results in ROI2 at a data count rate of 1×10 ⁷ .

Figure 5(b) is a line graph of variance of the results of each frame of Zubal thoracic data image in ROI2 at a data count rate of 1×10 ⁷ .

Figure 6(a) is the 14th frame PET image result reconstructed by the ML-EM method on the Monte Carlo simulated Zubal chest data at a data count rate of 1×10 ⁷ .

Figure 6(b) shows the result of the 14th frame PET image reconstructed from the Monte Carlo simulated Zubal thoracic data using the LRTV method at a data count rate of 1×10 ⁷ .

FIG. 6( c ) is the result of the 14th frame PET image reconstructed by the method of the present invention on the Monte Carlo simulated Zubal chest data under the data count rate of 1×10 ⁷ .

Figure 7(a) is the 14th frame PET image result reconstructed from Monte Carlo simulated Zubal thoracic data by ML-EM method at a data count rate of 1×10 ⁶ .

Figure 7(b) is the 14th frame PET image result of reconstruction of Monte Carlo simulated Zubal thoracic data by LRTV method under the data count rate of 1×10 ⁶ .

FIG. 7( c ) is the result of the 14th frame PET image reconstructed from the Monte Carlo simulated Zubal thoracic cavity data using the method of the present invention under the data count rate of 1×10 ⁶ .

Figure 8(a) is a real image of the 11th frame of Hoffman's brain data.

Figure 8(b) is the 11th frame PET image result reconstructed from the Monte Carlo simulated Hoffman brain data using the ML-EM method at a data count rate of 1×10 ⁷ .

Figure 8(c) is the 11th frame PET image result reconstructed from the Monte Carlo simulated Hoffman brain data using the LRTV method at a data count rate of 1×10 ⁷ .

Fig. 8(d) shows the result of the 11th frame of PET image reconstructed from Monte Carlo simulated Hoffman brain data using the method of the present invention at a data count rate of 1×10 ⁷ .

Detailed ways

In order to describe the present invention more specifically, the technical solutions of the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

As shown in FIG. 1, the present invention is based on a dynamic PET image reconstruction method based on non-local total variation and low-rank constraints, including the following steps:

S1. Establish a measurement data matrix Y and a system matrix D according to the dynamic PET scan mode.

The scanning process of dynamic PET is divided into a certain number of time frames according to the needs, and the measurement data matrix Y of dynamic PET can be constructed according to the order of time according to the counting vector collected by the detector in each time frame; The probability that the photons emitted at the pixel point are received by each detector, so as to obtain the system matrix D.

S2. Establish a dynamic PET imaging model.

For dynamic PET imaging, it requires sampling a series of consecutive time frames of the tracer's temporal information. For each independent time frame, the projection data y represents the sum of the coincident events captured by each detector during that time period, ie y={y _i , i=1,2,...,I}, where I is the total number of detectors; the corresponding radioactivity concentration image is recorded as x={ _xj ,j=1,2,...,J}, where J is the total number of pixels; the measured data y is related to the unknown concentration The relationship between graph x is:

Due to the independence of the Poisson assumption, we have the likelihood function for y:

For ease of solution, we take the logarithm of this likelihood function and minimize the negative logarithm of the likelihood function:

y_im代表第m帧被第i个探测器所捕获的投影数据；同样的，

x_jm代表第m帧第j个像素点的放射性浓度；我们对每一帧的负对数似然函数求和可得到：For dynamic PET, we can combine all the time frame information to get a data matrix. For the mth frame of projection data, the vector y _m represents the mth column of the data matrix Y,

y _im represents the projection data captured by the ith detector at the mth frame; similarly,

x _jm represents the radioactive concentration of the jth pixel in the mth frame; we sum the negative log-likelihood functions of each frame to get:

是第m帧

的第i个条目。where: y _im represents the value of the i-th detector of the m-th frame projection data y _m ,

is the mth frame

the ith entry of .

S3. Low rank and sparsity constraints.

A low-rank and sparse decomposition of PET images is performed, where the L component includes periodic changes in the radioactive intensity in the background, and the S component refers to those inhomogeneous tissues with different metabolic rates. Based on this, the following decomposition model is established:

Scale the objective function (5) to a convex optimization problem:

Where: ||L|| _* represents the kernel norm of matrix L, that is, the sum of the singular values of L. ||S|| ₁ represents the l ₁ norm of S.

S4. Nonlocal Total Variation Constraints.

是第m帧图像向量。将每一列

转换至其矩阵形式

其中U×V＝J，被转换的

现在在空间

中。我们依次对每一帧图像矩阵

进行非局部全变分约束：Let P=[p ₁ ,p ₂ ,...,p _m ,...,p _M ] represent an image sequence of M frames, where

is the image vector of the mth frame. put each column

Convert to its matrix form

where U×V=J, the converted

in space now

middle. For each frame of the image matrix in turn, we

Make a nonlocal total variation constraint:

Where: u, v are pixels in Ω space, which are natural numbers, w(u, v) is a non-negative symmetric non-local weight function, G _δ is a Gaussian kernel with standard deviation δ, and h is a filtering parameter.

S5. LRNLTV reconstruction of dynamic PET images.

Applying the non-local total variation constraints to the decomposed L and S respectively, combined with the previous low-rank and sparse constraints, we have the following objective function:

s.t. L+S=X

Where: λ, μ, α _L and α _S are weight coefficients.

S6. An optimization algorithm based on the augmented Lagrange multiplier method.

By introducing auxiliary variables P and Q, write the augmented Lagrangian function of the objective function (8):

Among them: Z, Z _L , Z _S are Lagrange multipliers, β, β _L , β _S are penalty parameters; the problem is decomposed into five sub-problems to solve, and according to the nature of each sub-problem, it is divided into Three categories.

6.1 Solve the L, S subproblems:

For the L sub-problem, let other variables be fixed, only the items related to L are extracted, and the formulas are sorted out, as follows:

Using the singular value threshold method to solve Equation (10), the iterative update formula of L is:

Where: A _ε (Γ)=UB _ε (s)V ^T , UsV ^T is the singular value decomposition of Γ, and the soft shrinkage algorithm B _ε (s) is:

Also for the S subproblem, after simplifying the formula, we have:

Using the soft shrinkage algorithm to solve Equation (12), the iterative update formula of S is:

6.2 Solve the X subproblem:

表示第m帧的符合线i处检测到的来自像素j的发射光子：First write the negative likelihood function Ψ(w|X) for the hidden variable w, where

represents the emitted photon from pixel j detected at coincidence line i for the mth frame:

Next, use the EM algorithm to optimize X.

并将其插入Ψ(w|X)；然后我们有替代函数φ(X；X^k)：Step E: Take the conditional expectation of w,

and plug it into Ψ(w|X); then we have the surrogate function φ(X; X ^k ):

Step M: Calculate the derivative of φ(X; X ^k ) with respect to x _jm and make it equal to 0; after simplification, it is found that the solution of x _jm is the root of the following quadratic equation:

Equation (15) is a convex function, and we update x _jm by taking the larger heel:

st a _jm = β,

where: [ ] _jm represents the jth entry of the matrix.

6.3 Solve P, Q subproblems:

For the P sub-problem, the other variables are also fixed, only the items related to L are extracted, and the formulas are sorted out, as follows:

We use gradient descent to solve equation (17), first we write the Euler-Lagrange of equation (7):

For a fixed u, there are:

转换至其矩阵形式

Z_L的每一列也同样被转换至其矩阵形式，其中U×V＝J，并且L和Z_L的每一帧被依次处理；则对于图像的每一帧，有

的迭代更新式：We put each column of the L matrix

Convert to its matrix form

Each column of Z _L is also converted to its matrix form, where U × V = J, and each frame of L and Z _L is processed in turn; then for each frame of the image, there are

The iterative update formula of :

是梯度下降的步进。在所有帧都被处理之后，将矩阵

转换回至向量形式

即可重新得到L，用同样的方式恢复P和Z_L，则P子问题被解决。in,

is the step of gradient descent. After all frames have been processed, the matrix is

Convert back to vector form

Then L can be obtained again, and P and Z _L can be recovered in the same way, then the P sub-problem is solved.

For the Q subproblem, since it has the same form as the P subproblem, it is solved in the same way; where the Lagrange multipliers are updated as usual.

The following is our experiment on the Zubal thoracic cavity and Hoffman brain template data simulated by Monte Carlo to verify the accuracy of the reconstruction results of the system of the present invention. Figure 2(a) and Figure 2(b) are schematic diagrams of the Zubal thoracic and Hoffman brain data templates used in the experiment, respectively. Different regions are divided into three regions of interest (where ROI4 is the background region). The experimental operating environment is: 8G memory, 3.40GHz, 64-bit operating system, CPU is intel i7-3770; the model of the simulated PET scanner is Hamamatsu SHR-22000, the set radionuclide and drug are ¹⁸ F-FDG, The sinogram is set as the result of data collected by 64 beams at each angle with 64 projection angles, and the size of the system matrix D is 4096×4096. In this experiment, the projection data under two different count rates of 1×10 ⁶ and 1×10 ⁷ are tested.

The PET image results of the reconstruction framework of the present invention are respectively compared with the image results of the two reconstruction methods ML-EM and LRTV, both of which use the same measurement data matrix Y and system matrix D to control the comparability of the results. It can be seen from Figures 3(a) to 3(d) that the image results of the reconstruction framework of the present invention have significantly less noise in the area than the images reconstructed by the other two methods. Under the condition of ensuring edge contrast, the functional area is smoother. Figures 4(a) to 4(d) are the enlarged image results of the parts framed in Figures 3(a) to 3(d) respectively. It can be clearly seen that the reconstructed image results of the present invention have a clearer ROI The boundary information has a significant improvement effect on early lesion detection. Figures 5(a) to 5(b) show the quantization error of each frame of the ROI2 of the Zubal thoracic cavity data, which further illustrates the accuracy of the reconstruction result of the present invention. Figures 6(a) to 6(c) and Figures 7(a) to 7(c) are the results of the 14th frame of Zubal thoracic data with three reconstruction methods at a count rate of 1 × 10 ⁷ and 1 × 10 ⁶ , respectively. The reconstruction result verifies that the reconstruction result of the present invention is robust to noise, and Table 1 shows its further quantitative result analysis. Figures 8(a) to 8(d) are reconstruction experiments performed on Hoffman brain data, which show that the reconstruction results of the present invention are robust to different template data.

Table 1

The above description of the embodiments is for the convenience of those of ordinary skill in the art to understand and apply the present invention. It will be apparent to those skilled in the art that various modifications to the above-described embodiments can be readily made, and the general principles described herein can be applied to other embodiments without inventive effort. Therefore, the present invention is not limited to the above-mentioned embodiments, and improvements and modifications made to the present invention by those skilled in the art according to the disclosure of the present invention should all fall within the protection scope of the present invention.

Claims (6)

1. A dynamic PET image reconstruction method based on non-local total variation and low-rank constraints, comprising the following steps: (1) Use a detector to detect the biological tissue injected with radiopharmaceuticals, dynamically collect the coincidence count vectors corresponding to each moment, and combine these coincidence count vectors into a coincidence count matrix Y; (2) Combining the dynamic PET image sequence into a PET concentration distribution matrix X, and establishing a PET measurement equation according to the PET imaging principle; (3) By introducing a low-rank constraint to the PET measurement equation, a dynamic PET image reconstruction model M1 based on the low-rank constraint is obtained; (4) The dynamic PET image reconstruction model M2 based on the non-local total variation constraint is further obtained by performing non-local total variation constraints on the low-rank and sparse parts of each frame of images; (5) The objective function of dynamic PET reconstruction is obtained by combining the dynamic PET image reconstruction models M1 and M2 as follows:

s.t.L+S=X where: || || _* represents the nuclear norm, || || ₁ represents the 1-norm, L is the low-rank part containing the periodically varying radioactive concentration in the image background, S is the sparse part containing the radioactivity with different metabolic rates The radioactive concentration of heterogeneous tissue, J _NLTV (L) and J _NLTV (S) are the results of the low-rank part L and the sparse part S after non-local total variation, respectively, λ, μ, α _L and α _S are weights coefficient, Ψ(Y|X) is the likelihood function about X and Y; (6) After the above objective function is optimized and solved, the PET concentration distribution matrix X is obtained, thereby restoring the dynamic PET image sequence. 2. dynamic PET image reconstruction method according to claim 1 is characterized in that: the expression of described PET measurement equation is as follows: Y=DX+R+S Where: D is the system matrix, R and S are the measurement noise matrices reflecting random events and scattering events, respectively. 3. dynamic PET image reconstruction method according to claim 1, is characterized in that: the expression of described dynamic PET image reconstruction model M1 is as follows: M1=||L|| _* +λ||S|| ₁ +μΨ(Y|X) s.t.L+S=X. 4. dynamic PET image reconstruction method according to claim 1, is characterized in that: the expression of described dynamic PET image reconstruction model M2 is as follows: M2=α _L J _NLTV (L)+α _S J _NLTV (S)+μΨ(Y|X) s.t.L+S=X. 5. The dynamic PET image reconstruction method according to claim 2, wherein the expression of the likelihood function Ψ(Y|X) is as follows:

where: d _ij is the element value of the i-th row and the j-th column in the system matrix D, y _im is the element value of the i-th row and the m-th column in the coincidence count matrix Y, and x _jm is the j-th row and the m-th column of the PET concentration distribution matrix X Element value, rim is the element value of the ith row and mth column in the measurement noise matrix R, _sim is the element value of the ith row and the mth column in the measurement noise matrix S, i, j and m are all natural numbers and _1≤i≤ I, 1≤j≤J, 1≤m≤M, I is the dimension conforming to the count vector, J is the row number of the PET concentration distribution matrix X, that is, the number of pixels of the PET image, M is the column of the PET concentration distribution matrix X The number is the sampling time length. 6. dynamic PET image reconstruction method according to claim 1, is characterized in that: in described step (6), adopt ADMM algorithm to carry out optimal solution to objective function; Wherein, low rank constraint problem utilizes singular value threshold algorithm and soft The shrinkage algorithm is used for iterative optimization, the likelihood function problem is solved iteratively with the EM algorithm, and the non-local total variational constraint problem is solved iteratively with the gradient descent method.

Images