Super-resolution PET imaging using convolutional neural networks

TL;DR

This paper introduces a CNN-based super-resolution technique for PET imaging that leverages MRI data and spatial location information to significantly improve resolution and quantitative accuracy over classical methods.

Contribution

The study presents a novel CNN architecture incorporating spatial location and MRI data for PET super-resolution, validated through simulation and clinical neuroimaging datasets.

Findings

CNNs outperform classical deconvolution methods in resolution recovery.

Deeper CNNs generally yield better image quality.

Incorporating MRI data enhances super-resolution performance.

Abstract

Positron emission tomography (PET) suffers from severe resolution limitations which limit its quantitative accuracy. In this paper, we present a super-resolution (SR) imaging technique for PET based on convolutional neural networks (CNNs). To facilitate the resolution recovery process, we incorporate high-resolution (HR) anatomical information based on magnetic resonance (MR) imaging. We introduce the spatial location information of the input image patches as additional CNN inputs to accommodate the spatially-variant nature of the blur kernels in PET. We compared the performance of shallow (3-layer) and very deep (20-layer) CNNs with various combinations of the following inputs: low-resolution (LR) PET, radial locations, axial locations, and HR MR. To validate the CNN architectures, we performed both realistic simulation studies using the BrainWeb digital phantom and clinical neuroimaging data analysis. For both simulation and clinical studies, the LR PET images were based on the Siemens HR+ scanner. Two different scenarios were examined in simulation: one where the target HR image is the ground-truth phantom image and another where the target HR image is based on the Siemens HRRT scanner --- a high-resolution dedicated brain PET scanner. The latter scenario was also examined using clinical neuroimaging datasets. A number of factors affected relative performance of the different CNN designs examined, including network depth, target image quality, and the resemblance between the target and anatomical images. In general, however, all CNNs outperformed classical penalized deconvolution techniques by large margins both qualitatively (e.g., edge and contrast recovery) and quantitatively (as indicated by two metrics: peak signal-to-noise-ratio and structural similarity index).

Tables

Table 1. Table 1: CNN architectures

Network	S1	V1	S2	V2	S3	V3	S4	V4
properties
Number	3	20	3	20	3	20	3	20
of layers
Input types	LR PET	LR PET	LR PET	LR PET	LR PET	LR PET	LR PET	LR PET
			HR MR	HR MR	Radial	Radial	HR MR	HR MR
					locations	locations
					Axial	Axial	Radial	Radial
					locations	locations	locations	locations
							Axial	Axial
							locations	locations

Table 2. Table 2: Simulation and experimental studies

Study index	Study type	LR image (input)	HR image (target)
1	Simulation	HR+ PET	True PET
2	Simulation	HR+ PET	HRRT PET
3	Clinical	HR+ PET	HRRT PET

Table 3. Table 3: LR and HR image sources

Image type	Scanner	Spatial resolution	Bore diameter	Axial length
LR	HR+	4.3 - 8.3 mm	562 mm	155 mm
HR	HRRT	2.3 - 3.4 mm	312 mm	250 mm

Table 4. Table 4: Study 1: Performance comparison

Metric	Reference	LR	TV	JE	S1	V1	S2	V2	S3	V3	S4	V4
PSNR	True	21.47	22.35	22.37	26.86	27.71	37.17	37.61	27.17	27.98	37.27	37.93
SSIM	True	0.74	0.85	0.86	0.74	0.86	0.97	0.97	0.83	0.87	0.97	0.98

Table 5. Table 5: Study 2: Performance comparison

Metric	Reference	LR	TV	JE	S1	V1	S2	V2	S3	V3	S4	V4
PSNR	HRRT	27.83	26.68	27.48	35.11	35.69	38.37	38.48	35.48	35.92	38.44	38.69
PSNR	True	21.47	22.35	22.37	22.69	23.06	22.46	23.92	23.07	23.48	23.78	24.29
SSIM	HRRT	0.76	0.83	0.82	0.82	0.80	0.88	0.88	0.85	0.82	0.89	0.90
SSIM	True	0.74	0.85	0.86	0.77	0.75	0.86	0.86	0.79	0.77	0.86	0.87

Table 6. Table 6: Study 3: Performance comparison

Metric	Reference	LR	TV	JE	S1	V1	S2	V2	S3	V3	S4	V4
PSNR	HRRT	29.93	34.98	33.32	35.75	36.93	36.00	37.27	36.29	37.42	36.46	38.02
SSIM	HRRT	0.86	0.93	0.88	0.94	0.95	0.95	0.95	0.95	0.95	0.95	0.96

Equations

$$z_{ij}^{lk}={{\mathbf{w}}^{lk}}^T{\mathbf{x}}_{ij}^l+b^{lk},$$

$$\mathbf{z}=\max (\mathbf{x},0),$$

$${\Phi }_{\text{JE}}(\mathbf{x}|\mathbf{y})=-\sum _{i=1}^M\sum _{j=1}^M\delta u\delta vp(u_i,v_j)\log p(u_i,v_j).$$

$${\Phi }_{\text{TV}}(\mathbf{x})=(\Vert {\Delta }_1\mathbf{x}{\Vert }_1+\Vert {\Delta }_2\mathbf{x}{\Vert }_1+\Vert {\Delta }_3\mathbf{x}{\Vert }_1),$$

$$\text{PSNR}(\bm{\hat{x}},\bm{x})=20\log_{10}\Big{(}\frac{\max{(\bm{\hat{x}})}}{\text{RMSE}(\bm{\hat{x}},\bm{x})}\Big{)},$$

$$\text{RMSE}(\hat{\mathbf{x}},\mathbf{x})=\sqrt{\frac{1}{N}\sum _k({\hat{x}}_k-x_k{)}^2}.$$

$$\text{SSIM}(\hat{\mathbf{x}},\mathbf{x})=\frac{(2{\mu }_x{\mu }_{\hat{x}}+c_1)(2{\sigma }_{x\hat{x}}+c_2)}{({\mu }_x^2+{\mu }_{\hat{x}}^2+c_1)({\sigma }_x^2+{\sigma }_{\hat{x}}^2+c_2)}.$$

Full text

K****eywords Super-resolution  $\cdot$ CNN  $\cdot$ Deep learning  $\cdot$ PET  $\cdot$ MRI  $\cdot$ Multimodality imaging  $\cdot$ Partial volume correction

1 Introduction

Positron emission tomography (PET) is a 3D medical imaging modality that allows in vivo quantitation of molecular targets. While oncology [1] and neurology [2] are perhaps the fields where PET is of the greatest relevance, its applications are expanding to many other clinical domains [3, 4]. The quantitative capabilities of PET are confounded by a number of degrading factors, the most prominent of which are low signal-to-noise ratio (SNR) and intrinsically limited spatial resolution. While the former is largely driven by tracer dose and detector sensitivity, the latter is the driven by a number of factors, including both physical and hardware-limited constraints and software issues. The resolution limitations pose an even greater challenge when the target regions-of-interest (ROIs) are smaller. Physical and hardware-related factors limiting the spatial resolution include the non-collinearity of the emitted photon pairs, intercrystal scatter, crystal penetration, and the non-zero positron range of PET radionuclides [5, 6, 7]. On the software front, resolution reductions are largely a product of smoothing regularizers and filters commonly used within or post-reconstruction for lowering the noise levels in the final images [8]. Together image blurring and tissue fractioning (due to spatial sampling for image digitization) lead to the so-called partial volume effect that is embodied by spillover of estimated activity across different ROIs [9].

Broadly, the efforts to address the resolution challenge encompass both within-reconstruction and post-reconstruction corrections. The former family includes methods that incorporate image-domain or sinogram-domain point spread functions (PSFs) in the PET image reconstruction framework [10, 11, 12] and/or smoothing penalties that preserve edges by incorporating anatomical information [13, 14, 15, 16, 17, 18, 19] or other transform-domain information [20, 21, 22]. The latter family of post-reconstruction filtering techniques includes both non-iterative corrections [23, 24, 9, 25, 26, 27] and techniques that rely on an iterative deconvolution backbone [28, 29, 30] which is stabilized by different edge-guided or anatomically-guided penalty or prior functions [31, 32].

Unlike partial volume correction, strategies for which are often modality-specific, super-resolution (SR) imaging is a more general problem in image processing and computer vision. SR imaging refers to the task of converting a low-resolution (LR) image to a high-resolution (HR) one. The problem is inherently ill-posed as there multiple HR images that may correspond to any given LR image. Classical approaches for SR imaging involve collating multiple LR images with subpixel shifts and applying motion-estimation techniques for combining them into an HR image frame [33, 34]. Less computationally-intensive modern approaches to SR encompass the family of so-called “example-based" techniques [35, 36, 37, 38, 39, 40, 41], which exploit self-similarities or recurring redundancies within the same image by searching for similar “patches" or sub-images within a given image. These methods are particularly effective for super-resolving natural images with fine and repetitive textures and less meaningful for most medical image types. With the proliferation of deep learning techniques, many deep SR models have been proposed and have demonstrated state-of-the-art performance at SR tasks. Deep SR models commenced with a paper that proposed a 3-layer CNN architecture (commonly referred to in literature as the SRCNN) [42]. Subsequently, a very deep SR (VDSR) CNN architecture [43] that had 20 layers and used residual learning [44] was demonstrated to result in much-improved performance over the shallower SRCNN approach. More recently SR performance has been further boosted by leveraging generative adversarial networks (GANs) [45], although GAN training remains notoriously difficult. Our previous work on SR PET spans a wide gamut, including both (classical) penalized deconvolution based on joint entropy [46, 32] and deep learning approach based on the VDSR CNN [47]. This work is an extension of the latter effort.

In this paper, we design, implement, and validate several CNN architectures for SR PET imaging, including both shallow and very deep varieties. As a key adaptation of these models to PET imaging, we supplement the LR PET input image with its HR anatomical counterpart, e.g. a T1-weighted MR image. Unlike uniformly-blurred natural images, PET images are blurred in a spatially-variant manner [11, 48]. We, therefore, further adjust the network to accommodate spatial location details as inputs to assist the SR process. “Ground-truth" HR images are required for the training phase of all supervised learning methods. Though simulation studies are not constrained by this demand, it is a challenge for clinical studies where it is usually infeasible to obtain the “ground-truth" HR counterparts for LR PET scans of human subjects. To ensure clinical utility of this method, we extend it to training based on imperfect target images derived from a higher-resolution clinical scanner, thereby exploiting the SR framework to establish a mapping from an LR scanner’s image domain to an HR scanner’s image domain. In section 2 of this paper, we present the underlying network architecture. In section 3, we describe the simulation data generation steps and the network training and validation procedures. In section 4, we present simulation results comparing the performance of CNN-based SR with two well-studied reference approaches. A summary of our work and a discussion on the application of CNNs for SR PET imaging are presented in section 6.

2 Theory

2.1 Network Design

2.1.1 CNN Basics

CNNs typically contain three types of layers: convolutional layers, nonlinear activation layers, and pooling layers. Pooling layers, which ensure shift invariance, are useful for classification and object recognition tasks. For SR imaging, which is an estimation/regression task, it suffices for the CNN to contain only convolutional and activation layers. A convolutional layer contains an array of convolutional kernels which extracts linear features from an input based on local connectivity. The pixel intensity at location $(i,j)$ in the $k$th feature map for the $l$th layer, $z^{lk}_{ij}$, can be mathematically represented as [49]:

[TABLE]

where $\bm{x}^{l}_{ij}$ is the input patch at the $(i,j)$th pixel location and $\bm{w}^{lk}$ and $b^{lk}$ are the kernel weights and bias respectively. An activation layer introduces nonlinearities that enable the extraction of nonlinear features by the CNN. In deep neural networks, the rectified linear unit (ReLU) is a popular choice for the activation function as, unlike alternatives such as the sigmoid function, it does not exhibit the vanishing gradient problem [50]. The ReLU activation function is defined as:

[TABLE]

where $\bm{x}$ is the input vector and $\bm{z}$ is the output vector. In addition, a ReLU activation function can accelerate computation of a network model because ReLUs set input values that are negative to be zero.

2.1.2 Residual learning

The first CNN applied to SR imaging (SRCNN) [42] sought to directly estimate the SR image as the network output. In contrast to such methods which compute the latent clean image, we adopt the residual learning strategy to remove the latent clean image from the blurry observations. This strategy was originally described in the ResNet [44] architecture for image recognition and has proven particularly beneficial for very deep networks, for which training accuracy degrades with increasing network depth. A generalization of this idea based on residual blocks demonstrated even higher efficacy at solving the SR problem [51].

2.2 Network Inputs

One key contribution of this paper is the tailoring of the CNN inputs to address the needs of SR PET imaging. All CNNs implemented here have the LR PET as the main input. This is similar to the two papers that inspired this work [42, 43], both of which had the LR single- or multi-channel images as their inputs. To further assist the resolution recovery process, additional inputs are incorporated as described below:

2.2.1 Anatomical Inputs

In generating the SR PET images, we seek to exploit the similarities between the PET and its high-resolution anatomical counterpart. Most clinical and preclinical PET scanners come equipped with anatomical imaging capabilities, in the form of computed tomography (CT) or MR imaging, to complement the functional information in PET with structural information. As illustrated in the schematic in Fig. 1, we employ CNNs with multi-channel inputs, that include LR PET and HR MR input channels.

2.2.2 Spatial Inputs

We provide location details to the network via additional input channels in the form of patches representing radial and axial coordinate locations. In light of the cylindrical symmetry of PET scanners, we deem these sufficient for learning the spatially-variant structure of the blurring operator directly from the training data.

2.2.3 Fusion

For natural RGB images, where the three input channels typically have a high degree of structural similarity, usually the same set of kernels are effective for feature extraction. In contrast, our CNN architecture exhibits a higher degree of input heterogeneity. Our initial experiments indicated the need for greater network width to accommodate a diverse set of features based on the very different input channels. We found an efficient solution to this by using separate kernels at the lower levels and fusing this information at higher levels as demonstrated in Fig. 1.

2.3 Network Depth

A key finding of the VDSR paper was that the deeper the network, the better its performance. The paper showed that VDSR outperforms SRCNN by a great margin in terms of PSNR [43]. Here, we implement networks with different depths:

2.3.1 Shallow CNN

We designed and implemented a set of “shallow" networks loosely inspired by the 3-layer SRCNN [42]. The modifications used here are 1) residual learning and 2) modified inputs as described in section 2.2. It only has three convolutional layers, each followed by a ReLU, except for the last output layer. While, as a rule of thumb, the qualifier “deep" is applied to networks with three or more layers, we use the word “shallow" in this paper in a relative way.

2.3.2 Very Deep CNN

We designed and implemented a series of “very deep" CNNs based on the VDSR architecture in [43], but with a different input design described in section 2.2. The very deep SR networks all have 20 convolutional layers, each followed by a ReLU, except for the last output layer.

2.4 Network Types

We implemented, validated, and compared both shallow (3-layer) and very deep (20-layer) CNN architectures with varying numbers of inputs. For the rest of the paper, we refer to these configurations as S1, S2, S3, S4, V1, V2, V3, and V4 as summarized in Table 1.

3 Methods

3.1 Overview

In the following, we describe two simulation studies using the BrainWeb digital phantom and the ${}^{18}\mathrm{F}$-FDG radiotracer and a clinical patient study also based on ${}^{18}\mathrm{F}$-FDG. The inputs and target for the studies are summarized in Table 2. All LR PET images were based on the Siemens ECAT EXACT HR+ scanner. For the first simulation study, the HR PET images were the “ground-truth" images generated from segmented anatomical templates. For the second simulation study and the clinical study, the HR PET images were based on the Siemens HRRT scanner. Bicubic interpolation was used to resample all input and target images to the same voxel size of 1 mm$\times$1 mm$\times$1 mm with a 256$\times$256$\times$207 grid size. The LR and HR scanner properties are summarized in Table 3.

3.2 Simulation Setup

3.2.1 HR+ PSF Measurement

An experimental measure of the true PSF was made by placing $0.5$ mm diameter sources filled with ${}^{18}\mathrm{F}$-FDG inside the HR+ scanner bore $56.2$ cm in diameter and $15.5$ cm in length. The PSF images were reconstructed using ordered subsets expectation maximization (OSEM) with post-smoothing using a Gaussian filter. The PSFs were fitted with Gaussian kernels. We assumed radial and axial symmetry and calculated the PSFs at all other in-between locations as linear combinations of the PSFs measured at the nearest measurement locations. Interpolation weights for the experimental datasets were determined by means of bilinear interpolation over an irregular grid consisting of the quadrilaterals formed by the nearest radial and axial PSF sampling locations from a given point.

3.2.2 Input Image Generation for Studies 1 and 2

Realistic simulations were performed using the 3D BrainWeb digital phantom (http://brainweb.bic.mni.mcgill.ca/brainweb/). 20 distinct atlases with 1 mm isotropic resolution were used to generate a set of “ground-truth" PET images. The atlases contained the following region labels: gray matter, white matter, blood pool, and cerebrospinal fluid. Static PET images were generated based on a  1 hour long ${}^{18}\mathrm{F}$-FDG scan as described in our earlier paper [32]. This “ground-truth" static PET is referred to as “true PET" for the rest of the paper. The geometric model of the HR+ scanner was used to generate sinogram data. Noisy data was generated using Poisson deviates of the projected sinograms, a noise model widely accepted in the PET imaging community [52]. The Poisson deviates were generated with a mean of $10^{8}$ counts for the full scan duration of 3640 s. The data were then reconstructed using the OSEM algorithm (6 iterations, 16 subsets). The images were subsequently blurred using the measured, spatially-variant PSF to generate the LR PET images. In order to match HR PET image grid size, the LR PET images were interpolated into $256\times 256\times 207$ from the HR+ output size of $128\times 128\times 64$ using bicubic interpolation. T1-weighted MR images with 1 mm isotropic resolution derived directly from the BrainWeb database were used as HR MR inputs.

3.2.3 Target Image Generation for Study 1

In Study 1, our purpose was to train the networks to map LR PET images to “ground-truth" image domain. The target HR PET images are, therefore, the true PET images as explained previously in section 3.1.

3.2.4 Target Image Generation for Study 2

In Study 2, we trained the networks using simulated HRRT PET images as our target HR images. The geometric model of the HRRT scanner was used to generate sinogram data. Poisson noise realizations were generated for the projected sinograms with a mean of $10^{8}$ counts for a scan duration of 3640 s. The images were then reconstructed using the OSEM algorithm (6 iterations, 16 subsets). OSEM reconstruction results typically appear grainy due to noise. We, therefore, perform post-filtering with a 2.4 mm full width at half maximum (FWHM) 3D Gaussian filter. Since the intrinsic resolution of the HRRT scanner is in the 2.3-3.4 mm range, this step improves image quality without any appreciable reduction is resolution.

3.3 Experimental Setup

Clinical neuroimaging datasets for this paper were obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI, http://adni.loni.usc.edu/) database, a public repository containing images and clinical data from 2000+ human datasets. We selected 20 HRRT PET scans and the anatomical T1-weighted MPRAGE MR scans for clinical validation of our method. $10$ of the 20 subjects were from the cognitively normal category. The remaining 10 subjects had mild cognitive impairment. The full scan duration was $30$ minutes ($6\times 5$-minute frames). The OSEM algorithm (6 iterations, 16 subsets) was used for reconstruction.

3.3.1 Target Image Generation for Study 3

As with Study 2 (section 3.2.3), the OSEM-reconstructed HRRT PET images, which were grainy, were post-filtered using a 2.4 mm FWHM 3D Gaussian filter to suppress some of the noise without substantially reducing the image resolution. This led to the target HR PET images for the clinical study.

3.3.2 Input Image Generation for Study 3

The LR counterparts of the HRRT images were generated by applying the measured spatially-variant PSF of the HR+ scanner described in section 3.2.1 to the OSEM-reconstructed HRRT images. While, not directly derived from the HR+ scanner, the use of a measured image-domain PSF ensured parity in terms of spatial resolution with true HR+ images. Rigidly co-registered T1-weighted MR images with $1\times 1\times 1$ voxels were used as HR MR inputs. Cross-modality registration was performed using FSL (https://fsl.fmrib.ox.ac.uk)[53, 54].

3.4 Network Implementation, Training, and Validation

All networks were implemented on the PyTorch platform. Training was performed using GPU-based acceleration achieved by using an NVIDIA GTX 1080Ti graphic card. An $L_{1}$ loss function was used for network training. Training was performed using Adam, an algorithm for optimizing stochastic objective functions via adaptive estimates of lower order moments [55].

The cohort size (total $20$ subjects) was the same for all three studies. For all studies, training was performed using data from $15$ out of $20$ available subjects to predict a residual image that is an estimated difference between the input LR PET and the ground truth HR PET. The stride and padding for the convolution kernels were both set to $1$. The kernel size was set to $3\times 3$. All the convolutional layers had $64$ filters except for the last layer, which had only $1$ filter. The batch size was $10$. The learning rate was set to $3\times 10^{-4}$. The network was trained for $400$ epochs. We validated our results by employing the data from the remaining $5$ subjects.

3.5 Reference Approaches and Evaluation Metrics

For the SR PET image, denoted by a vector $\bm{x}\in\mathbb{R}^{N}$, and the HR MR image, denoted by a vector $\bm{y}\in\mathbb{R}^{N}$, where $N$ is the number of voxels, the JE penalty is defined as:

[TABLE]

Here $\bm{u}\in\mathbb{R}^{M}$ and $\bm{v}\in\mathbb{R}^{M}$ are intensity histogram vectors based on the PET and MR images respectively and $M$ is the number of intensity bins. The TV penalty is defined as:

[TABLE]

where $\Delta_{k}$ ($k=1,2,$ or $3$) are finite difference operators along the three Cartesian coordinate directions.

The evaluation metrics used here are defined below. The true and estimated images are denoted $\bm{x}$ and $\bm{\hat{x}}$ respectively. We use the notation $\mu_{x}$ and $\sigma_{x}$ respectively for the mean and standard deviation of $\bm{x}$.

3.5.1 Peak Signal-to-Noise Ratio (PSNR)

The PSNR is the ratio of the maximum signal power to noise power and is defined as:

[TABLE]

where the root-mean-square error (RMSE) is defined as:

[TABLE]

3.5.2 Structural Similarity Index (SSIM)

The SSIM is a well-accepted measure of perceived image quality and is defined as:

[TABLE]

Here $c_{1}$ and $c_{2}$ are parameters stabilizing the division.

4 Results

4.1 Simulation Results: Study 1

Fig. 2a showcases results from Study 1 — transverse slices from the HR MR, HR PET (same as the true PET in this case), LR PET and the SR imaging results from the two deconvolution methods and eight CNNs described in Table 1. The images are for a given subject from the validation dataset. Magnified subimages in Fig. 2b highlight artifacts/inaccuracies indicated by purple arrows that are observed in all the techniques that lack anatomical guidance, namely TV, S1, V1, S3, and V3. Comparison of the subimage pairs (S1, S3) and (V1, V3) illustrates that, in the absence a anatomical information, spatial information greatly enhanced image quality. Comparison of the subimage pairs (S1, V1) and (S3, V3) also shows that the addition of more convolutional layers (increased network depth) is also very effective in the absence of anatomical information. Since JE incorporates anatomical information, it shows better edge recovery than TV. The CNN results show better gray-to-white contrast than both JE and TV.

The PSNR and SSIM of the different methods are tabulated in Table 4. In terms of PSNR, the supervised CNN networks outperform the classical approaches by a wide margin. The networks that have anatomical guidance, S2, S4, V2, and V4 have better performance on both PSNR and SSIM than the networks S1, V1, S3, and V3, which lack anatomical guidance. The PSNR and SSIM figures for S1 vs. S3 shows that these metrics increase noticeably with spatial information for the shallow case.

4.2 Simulation Results: Study 2

Fig. 3a showcases results from Study 2 — transverse slices from the HR MR, true PET, HR PET, LR PET and the SR imaging results from the two deconvolution methods and eight CNNs described in Table 1. The images are for a given subject from the validation dataset. As with Study 1, the CNNs with MR-based anatomical inputs (S2, V2, S4, and V4) still produce the best results. However, since the target image is now a corrupt image with diminished structural similarity with the MR image, the margin of gain from using anatomical information is now reduced. Magnified subimages in Fig. 3b highlight artifacts/inaccuracies indicated by purple arrows that are more preponderant in this study compared to Study 1. Interestingly, for this more challenging problem, the deeper networks (V2 and V4) have reduced background noise variations than their shallower counterparts (S2 and S4), as indicated by purple arrows in the latter. The JE and TV images, which are unsupervised, are the same as those showcased in Fig. 2.

The PSNR and SSIM of the different methods are tabulated in Table 5. An additional goal for this study was to understand the variability in results that could be anticipated when imperfect HR images were used for training. We, therefore, computed two sets of PSNR and SSIM measures: one with respect to the target HR PET and another with respect to the true PET. Our results show that there is an overall reduction in performance when the true PET is used as the reference. This is expected with the HR PET used for training deviates substantially from the true PET. But a key observation here is that the CNNs exhibit consistent relative levels of accuracy for the two reference images. As with Study 1, anatomically-guided networks showed better performance than the non-anatomically guided networks and the two classical methods.

4.3 Experimental Results: Study 3

Fig. 4a showcases results from Study 3 — the HR MR, HR PET, LR PET and the SR imaging results from the two deconvolution methods and eight CNNs described in Table 1. For this study, the deeper networks (V1, V2, V3, and V4) produced visually sharper images than their shallower counterparts (S1, S2, S3, and S4). This is clearly evident from the magnified subimages in Fig. 4b. This is consistent with our previous observation that, in the absence of a strong contribution of the anatomical inputs, the extra layers lead to a stronger margin of improvement. That said, V4, which is deeper and uses anatomical information, led to the highest levels of gray matter contrast as highlighted by red arrows.

The PSNR, SSIM of the different methods are tabulated in Table 6. As displayed in the table, V4 continues to exhibit the best performance in terms of PSNR and SSIM. As in Studies 1 and 2, all CNN based methods outperformed TV and JE.

5 Discussion

Overall, our results indicate that CNNs vastly outperform penalized deconvolution at the SR imaging task. Among different CNN architectures, relative performance depends on the problem at hand. Our simulation and clinical studies all agree that deep CNNs outperform shallow CNNs and that the additional channels contribute to improving overall performance. The relative importance of anatomical and spatial input channel depends on the underlying structural similarity between the HR MR and the true PET.

It should be noted that the utilization of anatomical images in conjunction with functional images required co-registration. To ensure robustness, a well-tested standardized registration tool was used for the task. Since the registration is intra-subject, rigid registration based on mutual information suffices for this application.

A limitation of the CNNs designed and validated in this paper is that they are specific to inputs with noise and blur levels similar to those used in training. In other words, they lack portability. Even for a given LR-HR scanner pair, such as the HR+ and HRRT, the SR performance is expected to drop when the LR inputs are based on a much lower or higher tracer dose than the datasets used for training. Our future work will characterize the performance of these networks when input noise levels are varied. One remedial approach to address noise sensitivity is the use of transfer learning for easy retraining of the networks with a much smaller training dataset containing LR inputs with altered SNR. Another, more sophisticated strategy is to use the CNN output as a prior in reconstruction. We have previously used this approach and produced promising results for denoising [56] and anticipate that it will work for deblurring and SR by extension.

Another limitation of CNN-based SR, and perhaps the most significant one, is that it relies on supervised learning and, therefore, requires paired LR PET and HR PET images for training. This requirement is easy to address while training using simulated datasets. But paired LR and HR clinical scans are rare. To address this limitation, we are currently exploring self-supervised learning strategies based on adversarial training of generative adversarial networks, that circumvent the need for paired training inputs.

6 Conclusion

We have designed, implemented, and validated a family of CNN-based SR PET imaging techniques. To facilitate the resolution recovery process, we incorporated both anatomical and spatial information. In the studies presented here, the anatomical information was provided as an HR MR image. So as to easily provide spatial information, we supplied patches containing the radial and axial coordinates of each voxel as additional input channels. This strategy is well-consistent with standard CNN multi-channel input formats and, therefore, convenient to implement. Both simulation and clinical studies showed that the CNNs greatly outperform penalized deconvolution both qualitatively (e.g., edge and contrast recovery) and quantitatively (as indicated by PSNR and SSIM).

As future work, we will develop new networks based on self-supervised learning that will enable us to circumvent the need for paired training datasets. We will also explore avenues to ensure applicability of this method to super-resolve inputs with noise levels different from those in the training data. From an applications perspective, we are interested in using this technique to super-resolve PET images of tau tangles, a neuropathological hallmark of Alzheimer’s disease, with the goal of developing sensitive image-based biomarkers for tau.