An improved patch-based regularization method for PET image reconstruction

Introduction

Different methods of iterative image reconstruction have become increasingly widely used in positron emission tomography (PET) imaging, producing promising results (1-5). Assuming the PET sampling data being investigated have a Poisson distribution, iterative statistical reconstruction methods can model the physical detection process of the scan (6). Iterative reconstruction can additionally use noise statistics, accurate system modeling, and prior knowledge of images to improve the accuracy of reconstruction (7). However, the maximum-likelihood expectation-maximization (MLEM) algorithm (8,9) used in the iterative reconstruction process increases noise when the iteration reaches a certain point. This produces an estimation of the image via projection by finding the maximum-likelihood (ML) solution. To reduce this increased noise and improve the quality of the image, statistical reconstruction methods based on the penalized maximum likelihood (PML) principle (10-14) have been used increasingly in PET imaging. The penalty function controls spatial smoothness in PML reconstruction. While the most commonly used quadratic penalty often over-smooths sharp edges in reconstructed images, non-quadratic penalties preserve edges at the expense of blocky artifacts. Images with blocky artifacts therefore limit the use of non-quadratic penalties in clinical settings.

Many potential improvements for these image reconstruction methods have been explored. To suppress noise and improve image quality, Muller et al. (15) improved the MLEM Total Variation (MLEM-TV) algorithm by using the inverse scale space method with Bregman distances (15,16) that restores loss of contrast in the images. Yang et al. (17) designed a shift-variant quadratic penalty function for PML image reconstruction to improve lesion detectability. Karaoglanis et al. (18) explored the effect of low-count statistics on ordered subset expectation maximization regularized with median root prior (OS-MRP-OSL) reconstructed images. Tang et al. (19) proposed the use of dictionary learning (DL)-based sparse signal representation in the formation of the prior for maximum a posteriori (MAP) PET image reconstruction; this improved bias and contrast with comparable noise. Deidda et al. (20) proposed a list-mode-hybrid kernel expectation-maximization (LM-HKEM) reconstruction algorithm to maintain the benefits of anatomically-guided methods and overcome their limitations by iteratively incorporating synergistic information. Finally, Wang and Qi (21) have proposed patch-based edge-preserving penalties (In the text it will be called “patch-based regularization algorithm”) that use neighborhood patches instead of individual pixels in computing non-quadratic penalties. This preserves the sharp edges of images and reduces blocky artifacts.

The MLEM image update in the patch-based regularization algorithm restores the fine structural features of an image but often produces noise. The image smoothing step of the patch-based regularization algorithm suppresses noise in the image domain but may result in the loss of fine structural features. The pixel-by-pixel image fusion step of the patch-based regularization algorithm, which fuses the MLEM updated image and the smoothed image, involves a trade-off between preserving fine structural features and suppressing noise. The patch-based regularization algorithm is required to iterate more times to suppress noise and retain fine structure features. This method of reconstruction is unable to sufficiently preserve the fine structural features of images reconstructed from low-count data.

However, low-count data and fast reconstruction algorithms are required to obtain clear images in clinical diagnoses. In a recent study, we found that the MLEM-TV-FR method suppressed noise and preserved fine structural features of images from low-count data (22). In order to accelerate the algorithm convergence, reduce the number of iterations, and reduce reconstruction time, in this study, we proposed to improve on patch-based regularization algorithm by adding TV regularization (23-27) following the MLEM image update. We then extracted the lost fine structural features from the residual image between the TV regularized image and the fused image based on patch regularization using a feature refinement (FR) step (28,29). The lost fine structural features were then added back to the fused image. In this way, each iterated image obtained more structural information. This improved reconstruction method accelerates the reconstructed image convergence and while retaining fine structural features. A brain phantom simulation experiment and data obtained from a mouse study were used to assess the performance of this proposed improved image reconstruction method.

The remaining sections of this paper are organized as follows. Method section describes penalized likelihood (PL) reconstruction, improved patch-based regularization reconstruction, and evaluation methods. The FR step is introduced in detail in this section. In Results section, details of the brain phantom simulation experiment and mouse study are reported and the reconstruction results and quantitative evaluation are also presented. The conclusions and discussion follow this.

Methods

This study was approved by the University of California, Davis, Institutional Animal Care, and Use Committee.

Penalized likelihood reconstruction

Measured emission sinogram data y is a vector where each element y_i contains the number of counts accumulated in bin i during a PET scan. The measurement y can be assumed to be a collection of independent Poisson random variables due to the annihilation events (8,30). Poisson distribution can model ${\overline{y}}_i$ mean as follows:

$$\text{P}\left(y_i\text{|}{\overline{y}}_i\right)=\frac{e^{-{\overline{y}}_i}{\overline{y}}_i}{y_i ！}$$[1]

where ${\overline{y}}_i$, the mean of the measurement data is related to unknown activity distribution map x.

$$\overline{y}\approx Gx+r$$[2]

where $\text{G}=\left\{g_{i,j}\right\}\in {ℝ}^{I\times J}$ is the system matrix where I and J are the number of PET measurements and image pixels, respectively. r is the model of random and scattered coincidences. PL reconstruction, therefore, estimates activity distribution map x by maximizing the PL function, which is written as follows (8,31):

$$\begin{array}{l}\widehat{x}=\arg ma{x}_{x\ge 0}\text{Φ}\left(x\right),\\ \text{Φ}\left(x\right)=L\left(x\text{|}y\right)-\beta U\left(x\right)\end{array}$$[3]

where β is a regularization parameter to balance data fidelity with spatial smoothness. When β = 0, the image reconstruction method represents the ML estimate. L(x|y) is the log of the likelihood function, which is written as follows:

$$\text{L}\left(x|y\right)={\displaystyle \sum }\left(y_{i}ln\left({\displaystyle \sum }{g}_{i,j}{x}_j\right)-{\displaystyle \sum }{g}_{i,j}{x}_j\right)$$[4]

U(x) is a penalty term with a different method of calculation that suppresses image noise.

Pixel-based regularization

U(x) is the computed pixel -based as follows:

$$U\left(x\right)=\frac{1}{4}{\displaystyle \sum }{\displaystyle \sum }{\omega }_{jk}\psi \left(x_j-x_k\right)$$[5]

where ω_jk is a weighting factor related to the distance between pixel j and pixel k in the neighborhood N_j, and ψ(t) is a penalty function with continuous second-order derivatives written as follows where t is an independent variable:

$$\psi \left(t\right)=\delta \left(\frac{\left|t\right|}{\delta }-\log \left(1+\frac{\left|t\right|}{\delta }\right)\right)$$[6]

where δ is a hyperparameter that controls the shape of the penalty function. This pixel-based regularization reconstruction algorithm is highly sensitive to δ (21).

Patch-based regularization

In a past study, Wang and Qi (21) proposed a patch-based regularization algorithm insensitive to hyperparameter δ. The penalty term U(x) of this algorithm is written as follows:

$$U\left(x\right)=\frac{1}{4}{\displaystyle \sum }{\displaystyle \sum }\psi \left({\Vert E_j\left(x\right)-E_k\left(x\right)\Vert }_a\right)$$[7]

where E_j(x) is a feature vector comprising the intensity values of all pixels in the patch centered on pixel j. ‖E_j(x)–E_k(x)‖_α denotes the patch-based distance between pixel j and pixel k and is defined as follows:

$$\Vert E_j\left(x\right)-E_k{\left(x\right)}_a\Vert =\sqrt{{\displaystyle \sum }{a}_m{\left(x_{j_m}-x_{k_m}\right)}^2}$$[8]

where M is the total number of pixels in a patch, j_m is the m-th pixel in the pixel j patch, and k_m is the m-th pixel in the pixel k patch. Both pixel j and pixel k have the same geometric relationship with respect to their central pixels. α_m is a positive weighting factor equal to the normalized inverse spatial distance between the pixel j_m and pixel j where ${\displaystyle \sum }{a}_m=1$.

In the patch-based regularization algorithm, the MLEM image update restores fine structural features, but the image often contains noise. The smoothing step is used to suppress noise in the image domain, but may result in loss of fine structural features. The pixel-by-pixel image fusion step, which fuses the MLEM updated image and the smoothed image, involves a trade-off between preserving the fine structural features of the image and suppressing noise.

Improved patch-based reconstruction

Patch-based regularization reconstruction requires many iterations to suppress noise and retain the fine structural features of the original image. To accelerate the algorithm convergence, reduce the number of iterations, and reduce reconstruction time, we proposed to improve this patch-based regularization reconstruction method. We did this by adding TV regularization following the MLEM image update. An FR step was then used to extract the lost fine structural features from the residual image between the TV regularized image and the fused image based on patch regularization. These features were then added back to the fused image. With these added steps, each iteration of the image should gain more structural information. Using the sparsity of image gradient magnitude to calculate TV is one of the most commonly used methods.

$$\text{TV}\left(x\right)={\displaystyle \sum }\sqrt{{\left(x_{s,t}-x_{s-1,t}\right)}^2+{\left(x_{s,t}-x_{s,t-1}\right)}^2+\alpha }$$[9]

Here, s and t are indexes of the desired tracer distribution map location, and α is a small constant used to maintain differentiability concerning image intensity. We assigned α a value of 10⁻⁸ in this study.

The FR step mentioned above is described by the following equation:

x_FR = x + f ⊗ ν[10]

where x_FR is the feature-refined image, x is the pixel-by-pixel fused image, and ν is the residual image between the TV-based MLEM-updated image x_TV and the pixel-by-pixel fused image x. The symbol ⊗ denotes pointwise multiplication. f is a feature descriptor, which is defined as follows (28):

$$f=1-\left|\frac{2{\sigma }_{qp}+C}{{\sigma }_p^2+{\sigma }_q^2+C}\right|$$[11]

where the constant C is introduced for numerical stability (C =1.25 × 10^–6 in this study). Local statistics σ_p, σ_q, and σ_pq at pixel j are defined as ${\sigma }_p={\left(\frac{1}{N-1}{\displaystyle \sum }{\left(x\left(j\right)-P\left(j\right)\right)}^2\right)}^{1/2}$, ${\sigma }_q={\left(\frac{1}{N-1}{\displaystyle \sum }{\left(x_d\left(j\right)-Q\left(j\right)\right)}^2\right)}^{1/2}$ and ${\sigma }_{qp}=\frac{1}{N-1}{\displaystyle \sum }\left(x\left(j\right)-P\left(j\right)\right)\left(x_d\left(j\right)-Q\left(j\right)\right)$, where $\text{P}=\frac{1}{N}{\displaystyle \sum }x\left(j\right)$, $\text{Q}=\frac{1}{N}{\displaystyle \sum }{x}_d\left(j\right)$, and p_j and q_j denote two local image patches with a size of $\sqrt{N}\times \sqrt{N}$ centered at pixel j. These images patches were extracted from x and degraded image x_d obtained by applying a 2D Gaussian filter to x, respectively. The parameters of the Gaussian filter function were a filter size of 5×5 and a standard deviation of the Gaussian function of 10. The nature of the proposed model structure descriptor involves a contrast variation component and a structure correlation component. The former calculates the reduction of contrast variation caused by the degrading operation, and the latter is the structural correction between the original image and the degraded image. The value of each element of the feature descriptor image f falls within the interval [0, 1]; a larger value is correlated with a greater likelihood of belonging to part of the structure.

The feature descriptor, designed to distinguish structures from noise and artifacts, plays a vital role in our improved algorithm. As such, several scalar parameters in this algorithm should be carefully tuned. For example, an image patch size of 7 × 7 is a good choice to balance structure-detection capacity and computational efficiency in the FR step. Additionally, parameter C is included to avoid instability when ${\sigma }_p^2+{\sigma }_q^2$ is close to zero. C should therefore be assigned a small constant value (C =1.25 × 10^–6 in our study). The added TV regularization and FR steps do not increase the calculation amount of the algorithm, but accelerate the algorithm convergence, reduce the number of iterations, and reduce reconstruction time. This improved approach is described in Supplementary file 1.

Evaluation

Signal-to-noise ratio (SNR) was used to evaluate the noise reduction capability of different reconstruction algorithms and is defined as follows:

$$\text{SNR}=10lo{g}_{10}\frac{{\displaystyle \sum }{x}_{tru{e}_i}}{{{\displaystyle \sum }}_i^n{\left(x_i-x_{tru{e}_i}\right)}^2}$$[12]

where x is the reconstructed image, x_true is the true image, i is the i-th pixel value of the image, and n is the total number of pixels in the image.

Covariance (COV) was used to measure the correlation between the reconstructed image and the true image. This is defined as follows (30):

$$\text{COV}=\frac{{{\displaystyle \sum }}_{i=1}^n\left(x_i-\overline{x}\right)\left(x_{tru{e}_i}-{\overline{x}}_{true}\right)}{n-1}$$[13]

where $\overline{x}=\frac{1}{n}{\displaystyle \sum }{x}_i$ and ${\overline{x}}_{true}=\frac{1}{n}{\displaystyle \sum }{x}_{tru{e}_i}$ are the average value of reconstructed image x and true image x_true, respectively.

We quantitatively compared the mean of the tumor contrast recovery coefficient (CRC) and the COV of the tumor as a function of the iterative reconstruction process. The CRC of the reconstructed image is defined as follows (21):

$$\text{CRC}=\frac{\left|\overline{T}-\overline{B}\right|}{\overline{B}}/R_0$$[14]

where R₀=8 is the true contrast, is the mean activity of the tumor region, $\overline{B}$ and $\overline{T}$ is the mean activity of background regions.

Regional relative bias and relative variance were used to evaluate the reconstructed images and are defined as follows (32):

$$\text{bias}=\frac{1}{n}{\displaystyle \sum }\frac{\left|x_i-x_{tru{e}_i}\right|}{x_{tru{e}_i}}$$[15]

$$\text{variance}=\frac{1}{n-1}{\displaystyle \sum }{\left(\frac{x_i-x_{tru{e}_i}}{x_{tru{e}_i}}\right)}^2$$[16]

where x_i is the i-th pixel value of the reconstructed image, x_truei is the i-th pixel value of the true image, and n is the total number of pixels in the evaluated region.

Results

We designed a brain phantom simulation and a mouse study to evaluate the performance of our proposed improved method and compare it with the patch-based regularization reconstruction algorithm. These methods involved a patch size and neighborhood size for the patch-based regularization step of 3×3 pixels (21) and an image where all pixel values = 1 was used as the first iteration image. SNR, COV, CRC, regional relative bias, and relative variance were used to evaluate the performance of our improved method.

All codes in this study were run on an HP z840 Workstation with a 64-bit Ubuntu 16.04 operating system, an Intel^® Xeon (R) CPU E5-2687W v4 @ 3.00 GHz × 48.

Brain phantom simulation

This study simulated a GE Discovery ST PET scanner in a 2D mode where each ring has 420 crystals with a cross-section of 6.3 × 6.3 mm². In this experiment, a PET emission image was simulated using a 2D brain phantom, shown in Figure 1A. The computed tomography (CT) image shown in Figure 1B was used to generate the attenuation factors. The brain phantom emission image is shown in Figure 1A was forward projected to obtain a noise-free sinogram. To simulate mean random events and scatter events in 2D mode, a uniform background of 20% total true coincidences was added to the sinogram (21). Independent Poisson noise was then introduced. The total event coincidence number was 500 K.

Figure 1 (A) Simulated PET emission image and (B) the attenuation map generated from a CT image. The blue boxes in (A) represent the background regions. These brain images are shown in the same window [0, 100].

Figure 2 shows images that were reconstructed using the pixel-based regularization method, the patch-based regularization method, and our proposed improved method for three different hyperparameter values: δ =0.1, 0.01, and 0.001. These images were reconstructed using 50 iterations. Figure 2 shows that the pixel-based regularization method (Figure 2A) is sensitive to hyperparameter δ, which controls the shape of the penalty function, while the patch-based regularization method (Figure 2B) and our improved method (Figure 2C) are relatively insensitive to changes in δ. However, the structural features that appear in Figures 2C, are clearer than that of Figure 2B, as indicated by the green arrows.

Figure 2 Image reconstruction results with 50 iterations (structural features indicated by arrows). (A1,A2,A3) Pixel-based reconstructions; (B1,B2,B3) patch-based reconstructions; and (C1,C2,C3) improved reconstructions. From left to right, the images were reconstructed using δ=0.1, 0.01, and 0.001, respectively.

Figure 3 shows brain images that were reconstructed using the hyperparameter δ =10^–9. Figure 3A shows the true image, Figure 3B shows the patch-based regularization image reconstructed with 50 iterations, Figure 3C shows the patch-based regularization image reconstructed with 100 iterations, and Figure 3D shows the reconstruction result obtained using our improved method with 50 iterations. While 50 iterations were used for the images in both Figure 3B and 3D, the fine structural features of Figure 3D are much clearer than that of Figure 3B, as indicated by the green arrows. The sharpness of Figure 3C image is second only to Figure 3D, but the image in Figure 3C required more time reconstruction time.

Figure 3 Image reconstruction results using the parameter δ=10^–9 (structural features indicated by arrows). (A) True image; (B) patch-based reconstruction with 50 iterations; (C) patch-based reconstruction with 100 iterations; and (D) improved method reconstruction with 50 iterations.

As indicated in Table 1, the reconstruction of Figure 3C image took 83.637 seconds, while the reconstruction of Figure 3D image took only 41.851 seconds. The green dotted-line boxes in Figure 3A are marked regions of interest (ROIs). Table 1 indicates that our improved method yields better bias and variance performance for each ROI, as well as for the whole reconstructed image, than the patch-based regularization method. Our improved method also took less time than the patch-based regularization method to reconstruct images to the same quality.

Table 1 Statistical analysis of brain image reconstruction
Full table

A quantitative evaluation of the reconstruction results in Figure 3 is shown in Figure 4. Compared with the patch-based regularization algorithm, our improved method performs convergence more quickly during the image update process; this is shown in Figure 4A. Figure 4B indicates that our improved method can quickly achieve higher contrast recovery and covariance, and therefore better performance, than the patch-based regularization algorithm.

Figure 4 Evaluation of brain image reconstruction shown in Figure 3. (A). SNR as a function of the total number of iterations; and (B) tumor contrast recovery plotted against COV curves for the patch-based regularization method and the improved method. These curves show the results obtained with 20 iterations.

Mouse study

A mouse scan was performed with the approval of the University of California, Davis, Institutional Animal Care, and Use Committee. Regular scanner quality assurance (QA) was performed to ensure the scanner worked appropriately. A mouse weighing 19.4 g was injected with 21 MBq of ¹⁸F-fluorodeoxyglucose and positioned in an Inveon D-PET scanner for 15 minutes starting 30 minutes after administration the injection (33). Four hundred ninety-eight million counts were acquired; counts are defined as the total number of events in a sinogram. The MAP method was implemented in the Siemens software package to reconstruct the mouse image. A 2D coronal image extracted as a slice in the Y direction was taken as the true mouse image, shown in Figure 5. A 2D 159 × 256 sinograms were extracted from this image, taking the same slice in the Y direction. The sinogram count was altered by changing the number of frames scanned. In our study, we examined sinograms with counts of 5,000 K, 250 K, 50 K, 40 K, and 25 K.

Figure 5 This coronal image was taken as the true image in the mouse study. All mouse images are shown in the same window [0, 100].

A sinogram with a count of 50 K was reconstructed using the pixel-based regularization method, the patch-based regularization method, and our improved method using different hyperparameter values: δ =0.1, 0.01, 0.001, and 10^–9. The images in Figure 6 were reconstructed using 50 iterations. Figure 6 shows that the improved method was more stable with changes to hyperparameter δ than the patch-based regularization method. Figure 6 also shows that the fine structural features of the images reconstructed using the improved method were clearer than that reconstructed by patch-based regularization, as indicated by the green arrows.

Figure 6 Images of the scanned mouse are reconstructed with 50 iterations. The sinograms have a count of 50 K. (A1,A2,A3,A4) Pixel-based reconstruction; (B1,B2,B3,B4) patch-based reconstruction; and (C1,C2,C3,C4) improved reconstruction. From left to right, these images are reconstructed using hyperparameter values of δ=0.1, 0.01, 0.001, and 10^–9.

We also examined the effects that sinograms with different counts had on these reconstruction algorithms when using the same hyperparameter value of δ =10^–9, as shown in Figure 7. From left to right, the images were reconstructed from sinograms with counts of 5,000 K, 250 K, 40 K, and 25 K, respectively. From top to bottom, these images were reconstructed using the MLEM method with 50 iterations, the patch-based regularization method with 50 iterations, the patch-based regularization method with 100 iterations, and our improved method with 50 iterations. Noise is present in the images reconstructed using MLEM even when the count is 5,000 K, as shown in the first row and first column of Figure 7. Additionally, the noise increased as the count decreased (see the first row of Figure 7). The fine structural features seen in the images reconstructed using the patch-based method and our improved method, shown in the first column of Figure 7, are very similar. This indicates that both of these methods can preserve the fine structural features of an image when the count was 5,000 K. When the count was 250 K, some structural features in the images reconstructed using the patch-based regularization method are somewhat blurred, while the fine structural features in the images reconstructed using our improved method remain clear, as indicated by the green arrows in the second column of Figure 7. When the count is 40 K or 25 K, the images reconstructed by patch-based regularization are very blurred, mainly when the count is 25 K. However, the fine structural features of the image reconstructed using our improved method remain visible when the count is 25 K, as indicated by the green arrows.

Figure 7 Images of the scanned mouse reconstructed using a hyperparameter value of δ=10^–9. From left to right, the sinograms used to reconstruct these images had counts of 5000 K, 250 K, 40 K, and 25 K. From top to bottom, these images were reconstructed by the MLEM method with 50 iterations, the patch-based regularization method with 50 iterations, the patch-based regularization method with 100 iterations, and the improved method with 50 iterations.

For quantitative comparison, we plotted intensity profiles along the horizontal green dotted line indicated in Figure 7 (see Figure 8). The intensity profiles in Figure 8 show that the MLEM method produced noisy results and, as the count decreased, the intensity of the image became more random, deviating further from the true value. However, the patch-based regularization method with 50 iterations, the patch-based regularization method with 100 iterations, and our improved method with 50 iterations produced results that were close to the true value with a count of 5000 K. When a count of 250 K was used, our improved method performed better than the patch-based regularization method, and particularly better than the patch-based regularization method with 50 iterations, as shown in Figure 8B. When the count was 40 K or 25 K, the results of the patch-based regularization method with 100 iterations were better than the results obtained using the patch-based regularization method with 50 iterations. In comparison, the results of our improved method with 50 iterations were much better than the results obtained using the patch-based regularization method with 100 iterations. This is shown in Figure 8C,D, respectively.

Figure 8 Intensity profiles along the horizontal green dotted line (located on the 86^th row of the image) shown in Figure 7.

Theoretically, the dotted lines in green (the patch-based regularization method with 50 iterations) and blue (the patch-based regularization method with 100 iterations) near the high intensity spot (approximately pixel index 125) in Figure 8 should decrease as count decreases. However, these are higher in the 25 K panel than the 40 K panel. This is due to the blurring of large hotspots in images reconstructed using the patch-based method. This further illustrates that the performance of the patch-based method is relatively poor with low-count data. We added another blue line in Figure 7 to further determine whether image intensity changed with count. The quantitative results are shown in Figure 9, which indicates that, as count decreased, the intensity of the reconstructed images gradually decreased. However, at low counts, the images reconstructed using our improved method were better than that reconstructed by the patch-based method. The improved method with 50 iterations further reduced reconstruction time by half than the patch-based regularization method with 100 iterations, as shown in Table 2.

Figure 9 Intensity profiles along the horizontal blue line (located on the 102nd row of the image) indicated in Figure 7.

Table 2 Time statistics for mouse image reconstruction
Full table

Table 3 lists the green dotted-line boxes indicate the bias and variance values of the ROIs of images reconstructed from sinograms with counts of 40 K. The ROIs in Figure 7. As Table 3 shows, our improved method yielded better bias and variance performance for each ROI and the whole reconstructed image, than the patch-based regularization method.

Table 3 Statistical analysis of mouse image reconstruction with a count of 40 K
Full table

Conclusions

This study shows that our proposed method improves on the patch-based regularization method for PET image reconstruction. The original patch-based regularization algorithm comprises three steps: an MLEM image update, an image smoothing step, and a pixel-by-pixel image fusion step. Two further key steps were introduced in our improved approach. One was TV regularization, which was executed following the MLEM update. This was done to suppress the noise generated by the MLEM update to facilitate the FR step. The FR step involved extracting the fine structural features of the image. This was performed following the MLEM update and TV regularization to restore fine structural features to the residual image between the TV-regularized MLEM-updated image and the pixel-by-pixel fused image. We tested our improved method by designing a brain phantom simulation experiment and a mouse study, and these results demonstrate that our improved method performs better and takes less time to reconstruct images than the patch-based regularization method. The advantage in using our improved method becomes more obvious when reconstructing PET images from low-count data. In conclusion, our improved method can reduce image reconstruction time and improve the quality of PET images. Our method also produces better reconstructions of PET images from low-count data.

Discussion

While the improved method tested in this study has been verified by 2D reconstruction, it has not been tested using 3D reconstruction. Reconstruction of 3D images takes a large amount of data, which greatly increases the calculation burden. Additionally, 3D reconstruction requires consideration of attenuation and scatter correction, which is relatively complicated; we intend to pursue further research in this area in the future. Another reason we only used 2D reconstruction to verify our reconstruction method was that, in most cases, 3D reconstruction requires graphics processing unit (GPU) technology, which was not available to us.

Parameter δ was used to control the shape of the penalty function. This parameter does not depend on 2D reconstruction, 3D reconstruction, or what regions are being scanned. As the pixel-based method is sensitive to δ, its value must be selected carefully when using this method. However, the patch-based method and our improved method were not sensitive to δ (see Figures 2 and 6). Using the patch-based method and the improved method to reconstruct image, relatively small δ values are generally selected, being approximately one-thousandth below the value of the image mean, such as 10^–3, 10^–6, or 10^–9. Parameter α is a small constant used to maintain differentiability concerning image intensities that can be set as any small value and it is not dependent on the object being imaged. This parameter can take the same value no matter what object, for example, digital brain phantoms, mice, or human structures, is being imaged. Additionally, the value of α can remain constant for both whole body (WB) PET imaging and the imaging of individual structures.

FR is a feature extraction operation performed in the image domain and is not dependent on the object being imaged. FR operates in the same manner for both WB PET imaging and the imaging of individual structures. Parameter C is a constant introduced for numerical stability and is also not dependent on the object being imaged. Whether WB PET imaging or individual body structure imaging is being performed, the value of C does not need to change.

Acknowledgments

The authors would like to thank our editor and anonymous reviewers for their constructive comments and suggestions. The authors thank Dr. Jinyi Qi and Guobao Wang of the Department of Biomedical Engineering at the University of California, Davis, for providing us with the patch-based regularization algorithm source code and for their assistance in implementing the algorithm.

Funding: This work was supported by the Guangdong Special Support Program of China (2017TQ04R395), the National Natural Science Foundation of China (81871441, 32022042), the Shenzhen International Cooperation Research Project of China (GJHZ20180928115824168), the Guangdong International Science and Technology Cooperation Project of China (2018A050506064), and the Natural Science Foundation of Guangdong Province in China (2020A1515010733).

Footnote

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at http://dx.doi.org/10.21037/qims-20-19). Dr. DL serves as an unpaid editorial board member of Quantitative Imaging in Medicine and Surgery. The other authors have no conflicts of interest to declare.

Ethical Statement: The mouse study described in this paper was approved by the University of California, Davis, Institutional Animal Care and Use Committee, in compliance with University of California–Davis Institutional Animal Care and Use Committee guidelines for the care and use of animals.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Shooli H, Dadgar H, Wang YXJ, Vafaee MS, Kashuk SR, Nemati R, Jafari E, Nabipour I, Gholamrezanezhad A, Assadi M, Larvie M. An update on PET-based molecular imaging in neuro-oncology: challenges and implementation for a precision medicine approach in cancer care. Quant Imaging Med Surg 2019;9:1597-610. [Crossref] [PubMed]
2. Blake GM, Puri T, Siddique M, Frost ML, Moore AEB, Fogelman I. Site specific measurements of bone formation using [F-18] sodium fluoride PET/CT. Quant Imaging Med Surg 2018;8:47-59. [Crossref] [PubMed]
3. Bhatt G, Jain A, Bhatt A, Civelek AC. Intramedullary spinal cord metastases and whole body F-18-FDG PET-CT-A case report. Quant Imaging Med Surg 2019;9:530-4. [Crossref] [PubMed]
4. Gu B, Xia L, Ge H, Liu S. Preoperative PET/CT score can predict complete resection in advanced epithelial ovarian cancer: a prospective study. Quant Imaging Med Surg 2020;10:743-53. [Crossref] [PubMed]
5. Rivlin M, Navon G. Molecular imaging of tumors by chemical exchange saturation transfer MRI of glucose analogs. Quant Imaging Med Surg 2019;9:1731-46. [Crossref] [PubMed]
6. Valiollahzadeh S, Clark JW, Mawlawi O. Dictionary learning for data recovery in positron emission tomography. Phys Med Biol 2015;60:5853-71. [Crossref] [PubMed]
7. Qi J, Leahy RM. Iterative reconstruction techniques in emission computed tomography. Phys Med Biol 2006;51:R541-78. [Crossref] [PubMed]
8. Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans Med Imaging 1982;1:113-22. [Crossref] [PubMed]
9. Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomogr 1984;8:306-16. [PubMed]
10. Yang L, Ferrero A, Hagge RJ, Badawi RD, Qi J. Evaluation of penalty design in penalized maximum-likelihood image reconstruction for lesion detection. J Med Imaging (Bellingham) 2014;1:035501. [Crossref] [PubMed]
11. Wang G, Qi J. Edge-Preserving PET Image Reconstruction Using Trust Optimization Transfer. IEEE Trans Med Imaging 2015;34:930-9. [Crossref] [PubMed]
12. Tahaei MS, Reader AJ. Patch-based image reconstruction for PET using prior-image derived dictionaries. Phys Med Biol 2016;61:6833-55. [Crossref] [PubMed]
13. Fessler JA. Penalized weighted least-squares image reconstruction for positron emission tomography. IEEE Trans Med Imaging 1994;13:290-300. [Crossref] [PubMed]
14. Yang L, Zhou J, Asma E, Wang WL. Evaluation of Penalized Maximum-Likelihood PET Image Reconstruction for ROI Quantification. 2016 IEEE Nuclear Science Symposium, Medical Imaging Conference And Room-Temperature Semiconductor Detector Workshop (Nss/Mic/Rtsd) 2016.
15. Muller J, Brune C, Sawatzky A, Kosters T, Schafers KP, Burger M. Reconstruction of Short Time PET Scans Using Bregman Iterations. 2011 IEEE Nuclear Science Symposium and Medical Imaging Conference (Nss/Mic) 2011:2383-5.
16. Brune C, Sawatzky A, Burger M. Primal and Dual Bregman Methods with Application to Optical Nanoscopy. Inter J Comput Vis 2011;92:211-29. [Crossref]
17. Yang L, Ferrero A, Hagge RJ, Badawi RD, Qi JY. Evaluation of penalty design in penalized maximum-likelihood image reconstruction for lesion detection. Medical Imaging 2014: Image Perception, Observer Performance, and Technology Assessment 2014;9037.
18. Karaoglanis K, Polycarpou I, Efthimiou N, Tsoumpas C. Appropriately regularized OSEM can improve the reconstructed PET images of data with low count statistics. Hell J Nucl Med 2015;18:140-5. [PubMed]
19. Tang J, Yang B, Wang YH, Ying L. Sparsity-constrained PET image reconstruction with learned dictionaries. Phys Med Biol 2016;61:6347-68. [Crossref] [PubMed]
20. Deidda D, Karakatsanis NA, Robson PM, Tsai YJ, Efthimiou N, Thielemans K, Fayad ZA, Aykroyd RG, Tsoumpas C. 2017 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC), Atlanta, GA, 2017:1-7.
21. Wang G, Qi J. Penalized likelihood PET image reconstruction using patch-based edge-preserving regularization. IEEE Trans Med Imaging 2012;31:2194-204. [Crossref] [PubMed]
22. Gao J, Zhang QY, Liu QG, Zhang XZ, Zhang MX, Yang YF, Liang D, Liu X, Zheng HR, Hu ZL. Positron emission tomography image reconstruction using feature extraction. J Xray Sci Technol 2019;27:949-63. [Crossref] [PubMed]
23. Ahn S, Kim SM, Son J, Lee DS, Lee JS. Gap compensation during PET image reconstruction by constrained, total variation minimization. Medical Physics 2012;39:589-602. [Crossref] [PubMed]
24. Liu Y, Ma JH, Fan Y, Liang ZR. Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction. Phys Med Biol 2012;57:7923-56. [Crossref] [PubMed]
25. Chen Z, Jin X, Li L, Wang G. A limited-angle CT reconstruction method based on anisotropic TV minimization. Phys Med Biol 2013;58:2119-41. [Crossref] [PubMed]
26. Zhu L, Niu T., Iterative CT. Reconstruction Via Minimizing Adaptively Reweighted Total Variation. Med Phys 2013;40:542. [Crossref]
27. Rose S, Andersen MS, Sidky EY, Pan XC. Noise properties of CT images reconstructed by use of constrained total-variation, data-discrepancy minimization. Medical Physics 2015;42:2690-8. [Crossref] [PubMed]
28. Liu QG, Luo JH, Zhu YM. Adaptive Image Decomposition by Improved Bilateral Filter. Int J Comput Appl 2011;23:16-22.
29. Hu Z, Zhang YW, Liu JB, Ma JH, Zheng HR, Liang D. A feature refinement approach for statistical interior CT reconstruction. Phys Med Biol 2016;61:5311-34. [Crossref] [PubMed]
30. Tang J, Yang B, Wang Y, Ying L. Sparsity-constrained PET image reconstruction with learned dictionaries. Phys Med Biol 2016;61:6347-68. [Crossref] [PubMed]
31. Wang G, Qi J. PET Image Reconstruction Using Kernel Method. Ieee Transactions on Medical Imaging 2015;34:61-71. [Crossref] [PubMed]
32. Chen S, Liu H, Shi P, Chen Y. Sparse representation and dictionary learning penalized image reconstruction for positron emission tomography. Phys Med Biol 2015;60:807-23. [Crossref] [PubMed]
33. Yang Y, Bec J, Zhou J, Zhang M, Judenhofer MS, Bai X, Di K, Wu Y, Rodriguez M, Dokhale P, Shah KS, Farrell R, Qi J, Cherry SR. A Prototype High-Resolution Small-Animal PET Scanner Dedicated to Mouse Brain Imaging. J Nucl Med 2016;57:1130-5. [Crossref] [PubMed]