Gamma-ray imaging systems1–3 are widely used for radiographic diagnosis, such as in digital radiography,4–6 medical imaging,7,8 and inertial confinement fusion (ICF) implosions.9–11 Recent results indicate that gamma-ray imaging systems are also valuable for studies of plasma physics: for example, the photon polarization characteristics of gamma rays could be utilized to reveal novel plasma mechanisms.12,13 The fundamental principle of gamma-ray imaging involves the use of a gamma-ray radiation source to irradiate the sample, followed by recording the gamma-ray radiation image attenuated by the sample. This recorded gamma-ray radiation image provides structural information about the sample. The spatial resolution of the imaging system plays a crucial role in accurately diagnosing the gamma-ray radiation image of the sample and serves as a core metric for evaluating the performance of the imaging system. However, during the process of imaging, the recorded image suffers degradation due to physical limitations inherent in the imaging equipment, resulting in the presence of noise, blurring, and downsampling.14 The noise mainly encompasses statistical fluctuations, camera readout noise, and the impact of scattered gamma rays. The blurring occurs owing to the dispersion of gamma rays during the imaging process. Moreover, limitations such as field of view (FOV), detection efficiency, and cost impose constraints on the pixel numbers of the cameras used in gamma-ray imaging systems, as a consequence of which the operation of the camera is inherently a discrete downsampling process. These factors collectively contribute to a reduction in the spatial resolution of the imaging system, posing challenges in meeting high-precision spatial diagnostic requirements.

Numerous image reconstruction algorithms have been proposed to enhance the spatial resolution of gamma-ray imaging systems. The point spread function (PSF), which reflects the response of an imaging system to a point source, serves as the foundation for image deblurring algorithms. Typically, the PSF of a gamma-ray imaging system can be considered spatially invariant, and the convolution operation is employed to characterize the blur of the imaging system. Consequently, deblurred images can be directly obtained using inverse filtering15 and Wiener filtering16 algorithms. However, these algorithms are highly susceptible to noise, leading to significant artifacts in the reconstructed images. The Richardson–Lucy algorithm,17 which adopts the Poisson noise distribution model, demonstrates some robustness to Poisson noise. Nevertheless, when the noise distribution deviates from a Poisson distribution or as the number of iterations increases, this algorithm will still amplify the noise. To address the noise more effectively, many algorithms based on maximum a posteriori (MAP) estimation18 have been proposed. From a Bayesian perspective, the core of these algorithms is that, in addition to the data fidelity term composed of the imaging model and the degraded image, an additional regularization term is added to describe the prior distribution of the latent clean image. Commonly utilized regularization terms in image reconstruction include Tikhonov regularization19 and total variation (TV) regularization,20 among others. However, these simple regularization terms will excessively smooth the reconstructed image while suppressing the noise. Several other algorithms have been proposed, such as the sparsity prior,21 the genetic algorithm,22 and simulated annealing.23 It is worth noting that although these algorithms can mitigate noise and blurring, they fail to resolve the high-frequency aliasing caused by the downsampling by the camera.

Single-image super-resolution (SISR)24–27 has been a prominent focus within the field of computer vision over the past two decades. The primary objective of SISR is to reconstruct an original high-resolution (HR) image from a noisy, blurred, and downsampled low-resolution (LR) image, and it seems to provide a comprehensive solution to the problem of decreasing spatial resolution in gamma-ray imaging systems. Given that the reconstructed HR image contains more pixels than the obtained LR image, SISR becomes a significantly underdetermined problem, necessitating the incorporation of advanced regularization terms (or image priors). Many sophisticated image priors have been proposed for SISR, such as nonlocally centralized sparse representation (NCSR),28 group-based sparse representation (GSR),29 and nonlocal similarity (NLS).30 With manually designed image priors, however, it is difficult to capture the intricate underlying information of natural images, resulting in poor quality of the reconstructed HR images. With the rapid development of convolutional neural networks (CNNs), their powerful nonlinear expressive ability has been exploited to enable direct modeling of the mapping relationship between LR images and HR images. SRCNN31 is a pioneering CNN-based approach for SISR. Subsequently, CNN-based SISR algorithms such as EDSR,32 SRFlow,33 and BayeSR34 have been proposed. CNN-based image reconstruction methods have also been applied to gamma-ray imaging.35,36 The combination of prior information and CNN-based models further improves the quality of reconstructed images.37 An improved CNN-based model called the fast Fourier transform neural network (FFTNN)38 has been proposed, which can achieve effective noise suppression and structure preservation. Although CNN-based algorithms have demonstrated significant performance improvement compared with traditional algorithms based on MAP estimation, they do entail an inherent limitation. Specifically, pretrained SISR algorithms based on CNNs are restricted to a single imaging model, necessitating distinct network models for different imaging models and representing a lack of flexibility. Therefore, it is worth exploring how to simultaneously retain the flexibility of MAP-based methods and the high performance of CNN-based methods.

Fortunately, the plug-and-play (PnP) framework39 presents an alternative solution. Based on MAP estimation, the PnP framework utilizes variable splitting techniques, such as the alternating direction method of multipliers (ADMM)40 and half-quadratic splitting (HQS),41 to decouple the data fidelity term and the regularization term. Consequently, an image prior trained using CNNs can be flexibly applied to replace the regularization term without affecting the data fidelity term. It has been shown that any image denoiser can serve as an implicit image prior.42 Inspired by the outstanding performance of CNN-based image denoisers, we aim to utilize a pretrained deep image denoiser to replace the traditional image denoiser with what is termed the deep denoiser prior. A variety of deep image denoisers suitable for the PnP framework have already been proposed, including IRCNN,43 FFDNet,44 and DRUNet.45 The combination of the PnP framework and the deep denoiser prior has found application in numerous image reconstruction fields, such as infrared imaging,46 hyperspectral image restoration,47 and snapshot compressive imaging.48

In this paper, we make the first attempt to apply the deep denoiser prior based on the PnP framework for SISR of a gamma-ray imaging system. Through the PnP framework, we maintain the flexibility of the MAP estimation method while introducing the more powerful data-driven image prior. The purpose of establishing the SISR model is to simultaneously address the noise, blurring, and downsampling of the gamma-ray imaging system. In addition, we focus on two issues: one is the boundary artifacts caused by the choice of the image boundary condition, and the other is the impact of the choice of downsampling factor (also super-resolution factor) on the reconstruction results. The rest of this paper is organized as follows. Section II introduces the gamma-ray imaging system and its SISR model, the proposed algorithm with estimated boundaries, and the implementation details. In Sec. III, the results of simulation and real experiments using the proposed algorithm are presented and discussed. Finally, a brief summary is given in Sec. IV.

A schematic of the gamma-ray imaging system is shown in Fig. 1(a). Typically, the gamma-ray source is positioned a long distance from the scintillator, and the gamma rays are almost vertically incident on the entrance surface of the scintillator after passing through the collimator. The sample is positioned in front of the entrance surface of the scintillator, and the gamma-ray radiation image, after attenuation by the sample, conveys two-dimensional structural information about the sample. Since the direct recording of the gamma-ray radiation image is difficult, the scintillator is employed to convert the gamma-ray image into a fluorescence image. For gamma-ray detection, organic crystals are commonly used as scintillators. A lutetium yttrium orthosilicate (LYSO) scintillator is selected as the gamma-ray image conversion screen in our imaging system owing to its high light yield (>30 000 photons/MeV) and high physical density (7.1 g/cm³).49 To prevent direct irradiation of the CMOS camera by gamma rays, the fluorescence photons are deflected 90° by the mirror and focused onto the incident plane of the CMOS camera by the lens; this is referred to as the lens-coupled scintillator.2 Lead shielding is implemented to enclose and protect the imaging system from the effects of scattered gamma rays. The CMOS camera converts the fluorescence photons into electrons and records them as a digital image, which is subsequently transmitted to a computer via an optical fiber for further image processing.

The spatial resolution of the gamma-ray imaging system comprises two components: the optical resolution, which is determined by the blur of the lens-coupled scintillator, and the image resolution, which is determined by the equivalent pixel size of the CMOS camera. To ensure that the obtained image has a sufficient signal-to-noise ratio (SNR), the LYSO scintillator needs to have an appropriate thickness to improve the detection efficiency of the gamma rays. The interaction process between the gamma rays and the LYSO scintillator is shown in Fig. 1(b). Gamma rays have a probability of interacting with the LYSO scintillator, primarily through processes such as Compton scattering, pair production, and the photoelectric effect, in which secondary electrons are produced. These secondary electrons disperse and deposit their energy within the LYSO scintillator, exciting it to emit fluorescence photons. These fluorescence photons continue to propagate within the LYSO scintillator, and a portion of them escaping from the exit surface of the LYSO scintillator can be received by the lens and transmitted to the CMOS camera. Typically, the spatial resolution of the lens-coupled LYSO scintillator can be evaluated by the full width at half maximum (FWHM) of the PSF.

Owing to the typically larger FOV of the gamma-ray imaging system in comparison with the small chip size of the CMOS camera, it is necessary to use a lens to reduce the fluorescence image. Considering the reduction factor of the lens, although the physical pixel size of the CMOS camera is small, the equivalent pixel size of the CMOS camera at the plane of the sample is large. In practice, the intensity distribution of the fluorescence image can be considered continuous, and the CMOS camera essentially performs a downsampling process. To further improve the image resolution of the CMOS camera, it is assumed that the image recorded by the CMOS camera is an LR image, while the image to be reconstructed is a HR image, as shown in Fig. 1(c). When setting the downsampling factor s, the HR image undergoes downsampling by s × s to obtain the LR image. SISR is the reverse of the downsampling process and involves reconstructing s × s pixels from a single downsampled pixel.

Under the assumption that the PSF of the lens-coupled scintillator is spatially invariant, the forward mathematical model of the gamma-ray imaging system can be expressed as$$y=(x\otimes h){\downarrow }_s+\epsilon ,$$(1)where $h\in {\mathbb{R}}^{k\times k}$ represents the PSF (or blur kernel) of the lens-coupled scintillator, k is the support size of the PSF, ⊗ represents the two-dimensional convolution operation, ↓_s represents the downsampling operation of the CMOS camera with a downsampling factor s, $y\in {\mathbb{R}}^{m\times m}$ is the LR image degraded by the gamma-ray imaging system, $x\in {\mathbb{R}}^{sm\times sm}$ represents the HR image of the gamma rays attenuated by the sample, and $\epsilon \in {\mathbb{R}}^{m\times m}$ is the additive noise.

More generally, by representing the blurring operation of the lens-coupled scintillator and the downsampling operation of the CMOS camera in matrix form, the forward mathematical model of the gamma-ray imaging system described by Eq. (1) can be further expressed as a linear equation$$y=SHx+n,$$(2)where $x\in {\mathbb{R}}^{s^2{m}^2\times 1}$, $y\in {\mathbb{R}}^{m^2\times 1}$, and $n\in {\mathbb{R}}^{m^2\times 1}$ are lexicographically stacked representations of the HR image, the LR image, and the additive noise, respectively. $H\in {\mathbb{R}}^{s^2{m}^2\times s^2{m}^2}$ and $S\in {\mathbb{R}}^{m^2\times s^2{m}^2}$ are representations in matrix form of the blurring operation and the downsampling operation, respectively.

SISR can be formulated as a linear inverse problem. Given the degraded LR image y, the blurring operation H, and the downsampling operation S, the objective of SISR is to find the underlying HR image x that best satisfies Eq. (2). Since the SISR is an underdetermined problem, there are theoretically an infinite number of solutions x that satisfy Eq. (2). Generally, image priors can be incorporated through MAP estimation to constrain the solution space of x, resulting in the following optimization problem:$$x=\underset{x}{\mathrm{arg\; min}}\frac{1}{2}\Vert y-SHx{\Vert }_2^2+\lambda \mathrm{\Phi }(x),$$(3)where $\frac{1}{2}\Vert y-SHx{\Vert }_2^2$ represents the data fidelity term, Φ(·) denotes the image prior adopted to constrain the solution x, and λ serves as a hyperparameter to balance the data fidelity term and the image prior.

To make the solution x derived from Eq. (3) as close as possible to the distribution of the original HR image, it is necessary to design a sophisticated image prior. However, this approach faces two challenges: first, it is difficult to design an image prior that can capture rich intrinsic image prior information, and, second, sophisticated image priors are often nonconvex or nondifferentiable, making it challenging to directly optimize Eq. (3) using methods such as gradient descent.50 Therefore, the core idea of the PnP framework is to decouple the data fidelity term and the image prior through variable splitting methods, thereby treating them as two independent subproblems that can be optimized separately.

In this paper, we use the HQS method41 to split the variable x in Eq. (3) into two variables x and z, which can be described as follows:$$x=\underset{x}{\mathrm{arg\; min}}\frac{1}{2}\Vert y-SHx{\Vert }_2^2+\lambda \mathrm{\Phi }(z)\mathrm{s}.\mathrm{t}.x=z,$$(4)where z represents the auxiliary variable corresponding to x. Then, Eq. (4) can be expressed as a bivariate optimization problem:$$L_{\rho }(x,z)=\frac{1}{2}\Vert y-SHx{\Vert }_2^2+\lambda \mathrm{\Phi }(z)+\frac{\rho }{2}\Vert x-z{\Vert }_2^2,$$(5)where ρ is a penalty parameter. Equation (5) can be addressed by iteratively solving the x and z subproblems:$$x_{k+1}=\underset{x}{\mathrm{arg\; min}}\frac{1}{2}\Vert y-SHx{\Vert }_2^2+\frac{\rho }{2}\Vert x-z_k{\Vert }_2^2,$$(6a)$$z_{k+1}=\underset{z}{\mathrm{arg\; min}}\lambda \mathrm{\Phi }(z)+\frac{\rho }{2}\Vert z-x_{k+1}{\Vert }_2^2.$$(6b)

Thus, the PnP framework achieves decoupling of the data fidelity term ${\frac{1}{2}\Vert y-SHx\Vert }_2^2$ and the image prior Φ(·). Regardless of the choice of Φ(·), the x subproblem described in Eq. (6a) remains a least-square problem only associated with the mathematical model of the imaging system.

By reformulating the z subproblem described by Eq. (6b), it can be rewritten as follows:$$z_{k+1}=\underset{z}{\mathrm{arg\; min}}\mathrm{\Phi }(z)+\frac{1}{2{\left(\sqrt{\lambda /\rho }\right)}^2}\Vert z-x_{k+1}{\Vert }_2^2.$$(7)

From a Bayesian perspective, Eq. (7) is precisely a Gaussian image denoising problem on x_k+1 with a Gaussian noise level $\sigma =\sqrt{\lambda /\rho }$. Therefore, we can use any off-the-shelf Gaussian image denoiser as an implicit image prior instead of explicitly designing an image prior.

In our previous work,51 we assessed the performance of three traditional nontrained Gaussian image denoisers, namely, TV,20 NLM,52 and BM3D,53 within the PnP framework for a neutron imaging system. The results demonstrated that stronger image denoisers yield higher-quality reconstructed images. Building upon this insight, to improve the quality of the reconstructed images, it is essential to select the most advanced Gaussian image denoiser whenever feasible. Motivated by the outstanding performance of CNN-based Gaussian image denoisers, we propose utilizing a dataset-trained Gaussian image denoiser to replace traditional nontrained Gaussian image denoisers.

Given that the hyperparameters λ and ρ are both adjustable, the image denoiser employed in each iteration of Eq. (7) needs to operate at a specified Gaussian noise level. This implies that the Gaussian image denoiser employed within the PnP framework should be nonblind and capable of handling a wide range of Gaussian noise levels. Although various CNN-based Gaussian image denoisers have been proposed recently, most models are designed for a single, fixed Gaussian noise level. By contrast, the authors of DRUNet45 incorporated the Gaussian noise level as an additional input when designing the CNN model. On the basis of the noise level, a pure Gaussian noise image with the same dimensions as the noisy image is first generated, concatenated with the noisy image, and then input into the CNN model. In this manner, DRUNet can handle images with arbitrary Gaussian noise levels. Furthermore, DRUNet combines the advantages of UNet54 and ResNet,55 further improving its denoising performance. It is well known that CNNs benefit from the availability of large-scale training data. A large natural image dataset, including BSD,56 DIV2K,57 and Flick2K,58 is used to train DRUNet. To handle a wide range of noise levels, the noise level σ is randomly sampled from the range of (0, 50) during training. The network parameters are optimized by minimizing the loss between noisy images and clean images with the Adam algorithm.59

On the basis of the above analysis, we choose DRUNet as an implicit image prior in this paper, referred to as a deep denoiser prior. Consequently, Eq. (7) can be rewritten as follows:$$z_{k+1}=DRUNet\left(x_{k+1},\sqrt{\lambda /\rho }\right).$$(8)

It is worth noting that the x_k+1 obtained in each iteration of Eq. (6a) is not an image completely corrupted by Gaussian noise, and the role of DRUNet employed in Eq. (8) is not intended to remove Gaussian noise from x_k+1. Instead, the purpose of using DRUNet here is to leverage its learned implicit image prior to make x_k+1 more consistent with the distribution of natural images.

Section II D discussed the use of DRUNet as a deep denoiser prior to solve the subproblem (6b). This subsection will focus on solving the subproblem (6a). Although the subproblem (6a) is ostensibly a simple least-square problem, the boundaries of the image x must be taken into account, owing to the convolution operation H. When performing convolution, the boundaries of an image play a significant role in determining the output. As shown in Fig. 2(a), given the blur kernel $h\in {\mathbb{R}}^{k\times k}$ and the original image $\widehat{x}\in {\mathbb{R}}^{(sm+k-1)\times (sm+k-1)}$, the blurred image is $x\in {\mathbb{R}}^{sm\times sm}$. This means that for a blurred image x, the pixels near the boundaries are influenced by unknown pixels outside the image itself. Therefore, when solving the subproblem (6a), certain boundary conditions usually assumed for x, such as periodic boundaries, zero boundaries, and reflexive boundaries.

The most commonly used assumption for the boundary conditions is that of periodic boundaries, as shown in Fig. 2(b). In this case, the convolution operation H is circulant, and it can be diagonalized by the discrete Fourier transform matrix. Therefore, Eq. (6a) can be efficiently solved using the Fourier transform and can be implemented in the following closed form:60$$x_{k+1}={\rho }^{-1}b-{\rho }^{-1}{(SH)}^T\left[{\mathcal{F}}^{-1}\left(\frac{\mathcal{F}(SHb)}{|\mathcal{F}({\tilde{h}}_0){|}^2+\rho }\right)\right],$$(9)where $\mathcal{F}(\cdot )$ is the Fourier transform operator, b = (SH)y + ρz_k, and ${\tilde{h}}_0=\left[{\mathcal{F}}^{-1}\left(\mathcal{F}(h)\overline{\mathcal{F}(h)}\right)\right]{\downarrow }_s$.

In the simulation experiments, the forward imaging process can be accurately simulated on the basis of the periodic boundary condition. Thus, the use of Eq. (9) to solve the subproblem (6a) provides good performance and high speed. However, for real imaging systems, the periodic boundary condition is obviously an incorrect assumption that will introduce boundary artifacts in the reconstructed image.61,62 To address these boundary artifacts, we do not assume any boundary conditions, but instead directly estimate the boundaries during the reconstruction process. Specifically, we introduce a mask matrix M, as shown in Fig. 2(c), between the blurring operation H and the downsampling operation D. Consequently, Eq. (6a) is rewritten as$${\widehat{x}}_{k+1}=\underset{\widehat{x}}{\mathrm{arg\; min}}\frac{1}{2}\Vert y-SM\widehat{H}\widehat{x}{\Vert }_2^2+\frac{\rho }{2}\Vert \widehat{x}-z_k{\Vert }_2^2,$$(10)where $M\in {\mathbb{R}}^{s^2{m}^2\times {(sm+k-1)}^2}$ is a lexicographically stacked representation of a mask. The entries in the mask can be either 1 or 0, indicating whether to preserve or delete the corresponding pixels. In this paper, the size of the border is used to denote the size of the cropped boundaries. Owing to the presence of the mask M, the reconstructed image $\widehat{x}\in {\mathbb{R}}^{{(sm+k-1)}^2\times 1}$ is larger than the reconstructed image $x\in {\mathbb{R}}^{s^2{m}^2\times 1}$ obtained by Eq. (9), and the corresponding convolution operation $\widehat{H}\in {\mathbb{R}}^{{(sm+k-1)}^2\times {(sm+k-1)}^2}$ is also larger than the convolution operation $H\in {\mathbb{R}}^{s^2{m}^2\times s^2{m}^2}$ in Eq. (9). In other words, Eq. (10) simultaneously reconstructs x while also estimating the missing boundaries of x, which together constitute $\widehat{x}$.

With the introduction of the mask M, Eq. (10) no longer has an exact solution similar to Eq. (9). By reformulating the $\widehat{x}$ subproblem described by Eq. (10), it can be rewritten as$${\widehat{x}}_{k+1}=\underset{\widehat{x}}{\mathrm{arg\; min}}\frac{1}{2}{‖\left[\begin{array}{c}\hfill SM\widehat{H}\hfill \\ \hfill \sqrt{\rho }I\hfill \end{array}\right]\widehat{x}-\left[\begin{array}{c}\hfill y\hfill \\ \hfill \sqrt{\rho }{z}_k\hfill \end{array}\right]‖}_2^2,$$(11)where $I\in {\mathbb{R}}^{s^2{m}^2\times s^2{m}^2}$ is the sparse identity matrix.

Since the S, M, $\widehat{x}$, and I are all large sparse matrices, the least-square QR-factorization (LSQR) method63 can be adopted to solve Eq. (11) quickly. Letting $G=\left[\genfrac{}{}{0.0pt}{}{SM\widehat{H}}{\sqrt{\rho }I}\right]$ and $r=\left[\genfrac{}{}{0.0pt}{}{y}{\sqrt{\rho }{z}_k}\right]$, Eq. (11) can be solved by$${\widehat{x}}_{k+1}=LSQR(G,r,1e^{-3}),$$(12)where LSQR is a wrapped function in Python and the tolerance is 1 × 10⁻³.

In summary, we use Eq. (12) to solve the $\widehat{x}$ subproblem and Eq. (8) to solve the z subproblem.

In this subsection, the hyperparameters, convergence, and stopping criteria of the proposed algorithm are discussed. The performance of the proposed algorithm in this paper relies on two crucial hyperparameters: λ and ρ. Experimental results demonstrate that decreasing the value of λ leads to gradual emergence of speckle noise in the reconstructed images, whereas increasing λ results in overly smooth images.51 Therefore, in all our experiments, we choose relatively small values of λ to maximize the potential of the algorithm while guaranteeing the absence of speckle noise in the reconstructed images. In the PnP framework, the hyperparameter ρ significantly affects the convergence of the algorithm. Although the global convergence of the algorithm within the PnP framework remains unproven, previous research has indicated that the algorithm exhibits fixed-point convergence as long as ρ is nondecreasing.60 We employ an adaptive criterion for updating ρ by considering the residual between two adjacent iterations:$${\mathrm{\Delta }}_k=\frac{1}{\sqrt{N}}\left(\Vert {\widehat{x}}_k-{\widehat{x}}_{k-1}{\Vert }_2^2+\Vert z_k-z_{k-1}{\Vert }_2^2\right)$$(13)where N is the pixel number of ${\widehat{x}}_k$ and z_k. By defining constants ρ₀ = 1, γ = 1.1, and τ = 0.9, if Δ_k > τΔ_k−1, then ρ_k = γρ_k−1, and if Δ_k ≤ τΔ_k−1, then ρ_k = ρ_k−1.

In addition, we define an iteration stopping criterion based on the residual Δ_k. If Δ_k ≤ tol, then the algorithm is considered to reach a local optimal solution and the iteration is stopped. In our tests, it has been found that setting tol ≈ 10⁻⁴ ∼ 10⁻³ is sufficient.

On the basis of the above discussion, we have summarized the proposed algorithm of SISR using DRUNet based on PnP framework, namely, PnP-DRUNet, in Algorithm 1.ALGORITHM 1.

Single image super-resolution using deep denoiser prior based on plug-and play framework (PnP-DRUNet).

We use three metrics to quantitatively evaluate the quality of super-resolution reconstructed images, including the peak-to-noise ratio (PSNR), the structural similarity (SSIM) index, and the contrast of line-pairs.

PSNR is used to evaluate the pixel error between the reconstructed image and the original image. The formula for calculating PSNR is$$\mathrm{P}\mathrm{S}\mathrm{N}\mathrm{R}=10{\log }_{10}\left[\frac{{\mathrm{M}\mathrm{A}\mathrm{X}}_I^2}{\frac{1}{mn}\sum _{i=1}^m\sum _{j=1}^n(x_{ij}-I_{ij})}\right],$$(14)where MAX_I denotes the maximum pixel value in I, and x and I are the reconstructed and original images, respectively. The unit of PSNR is dB. The higher the value of PSNR, the smaller the pixel error between the reconstructed image and the original image.

SSIM is used to evaluate the structural error between the reconstructed image and the original image. SSIM is expressed as$$\mathrm{S}\mathrm{S}\mathrm{I}\mathrm{M}=\frac{(2{\mu }_x{\mu }_I+c_1)(2{\sigma }_{xI}+c_2)}{({\mu }_x^2+{\mu }_I^2+c_1)({\sigma }_x^2+{\sigma }_I^2+c_2)},$$(15)where μ_x and σ_x are the mean and variance of the reconstructed image, μ_I and σ_I are those the original image, $c_1={(K_1{\mathrm{M}\mathrm{A}\mathrm{X}}_I)}^2$, and $c_2={(K_2{\mathrm{M}\mathrm{A}\mathrm{X}}_I)}^2$, with K₁ = 0.01 and K₂ = 0.03.

Both PSNR and SSIM require knowledge of the original image, which is typically unavailable in real experiments. For line-pair images, the contrast of line-pairs can be used to assess the spatial resolution, eliminating the need for the original image. Given the nonuniformity of the reconstructed line-pair images, we define the average contrast in this paper as follows:$$\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{s}\mathrm{t}=\frac{\frac{1}{m}\sum _{i=1}^m{x}_{max,i}-\frac{1}{n}\sum _{j=1}^n{x}_{min,i}}{x_{max}-x_{min}}\times 100\%,$$(16)where x_max,i and x_min,i are the local maxima and local minima, respectively, of the intensity curve of the line-pairs, as shown in Fig. 3m and n are the numbers of local maxima and local minima, respectively, and x_max and x_min are the maximum and minimum, respectively, of the line-pairs.

In all simulation experiments, the original HR gamma-ray images are discretized using a 50 μm grid. For convenience, a Gaussian blur kernel is employed to represent the PSF of the lens-coupled scintillator of the gamma-ray imaging system, with a σ = 200 μm. The FWHM of the PSF of the Gaussian blur kernel is $2\sqrt{2\ln 2}\sigma \approx 470\mu \text{m}$, indicating that the optical resolution of this gamma-ray imaging system is also about 470 μm. The image resolution of the LR images recorded by the gamma-ray imaging system is determined by the downsampling factor s. For example, when s = 4, the image resolution of the LR images is 200 μm.

To demonstrate the impact of boundaries on the reconstructed images, we utilize an original HR image with dimensions of 288 × 288 pixels, as shown in Fig. 4(a1). The pixel intensities have been normalized to a range of 0–1, and pseudocolor coding has been applied to provide additional information to the observer. During the convolution operation, owing to the size of the blur kernel (31 × 31), pixels at the boundaries are lost (border = 15 pixels), as illustrated in Fig. 4(a2). The downsampling factor s of the CMOS camera is set to 3%, and 1% Gaussian noise is added, resulting in a recorded LR image with dimensions of 88 × 88 pixels and a pixel size of 150 μm, as shown in Fig. 4(a3). The information near the boundaries of the recorded LR image still originates from the boundaries lost in the original HR image. Therefore, it is necessary to assume or estimate the boundaries during image reconstruction. To evaluate the performance of the proposed algorithm PnP-DRUNet, the PnP-BM3D51 in our previous work is adopted as a comparison algorithm, which has been proven to outperform the traditional Richardson–Lucy algorithm.17 The second row of Fig. 4 presents the reconstruction results using Eq. (9) under the assumption of periodic boundaries. Figures 4(b2) and 4(c2) are zoomed images near the boundaries of the reconstructed images obtained by PnP-BM3D and PnP-DRUNet, respectively. It can be observed that there are significant boundary artifacts in the reconstructed images. Additionally, ringing artifacts are visible in the central area of Fig. 4(b1). By contrast, the reconstructed images using Eq. (12), which simultaneously estimates the boundaries, do not exhibit boundary artifacts, and the visual quality of the reconstructed images is improved. Both PnP-BM3D and PnP-DRUNet can approximately correctly estimate the missing boundaries, as indicated by the pixels outside the black dashed lines in Figs. 4(d1) and 4(e1).

To quantitatively evaluate the performance of different algorithms, the results of PSNR and SSIM are shown in Table I. When calculating PSNR and SSIM, the pixels at the border are cropped. Regardless of the type of boundaries used, the proposed PnP-DRUNet always outperforms PnP-BM3D. More importantly, the method using estimated boundaries is significantly better than the assumption of periodic boundaries, and the resulting performance improvement is much greater than the improvement brought by the image denoiser itself. This indicates that for images where the pixels at the boundaries are nonzero, such as images cropped from a larger image within a region of interest (ROI), the proposed boundary estimation method in this paper is crucial.

Intuitively, as the downsampling factor s increases, the solution space for SISR also expands, posing a greater challenge to the reconstruction process. Meanwhile, for SISR using an image denoiser prior based on the PnP framework, improving the image denoiser can also enhance reconstruction performance. Nevertheless, the quantitative relationship between the performance improvement brought by advanced image denoisers and the performance degradation caused by an increased downsampling factor s remains unclear. To investigate this, we conduct a quantitative analysis of the relationship between image denoisers and the downsampling factor s. We design an original HR image containing rich shape information with dimensions of 288 × 288 pixels, as shown in Fig. 5(a). To ensure that the blurred image is an integer multiple of the downsampling factor s, we maintain the blur kernel size of 31 × 31 with a border of 15 pixels for downsampling factors s = 2 and 3. For a downsampling factor s = 4, the blur kernel size is chosen as 33 × 33, with a border of 16 pixels. Figures 5(b1)–5(b3) present the recorded LR images with 1% Gaussian noise at s = 2, 3, and 4, respectively. It can be observed that as the downsampling factor s increases, the image resolution of the recorded LR images gradually decreases, accompanied by the emergence of aliasing artifacts. The third and fourth columns of Fig. 5 show the reconstructed images using PnP-BM3D and PnP-DRUNet, respectively. From the reconstruction results with PnP-BM3D, it can be seen that as the downsampling factor s increases, the artifacts in the reconstructed images gradually increase. In comparison with PnP-BM3D, the proposed PnP-DRUNet yields reconstructed images with almost no artifacts, indicating the superior ability of the DRUNet image denoiser to accurately describe the distribution of natural images. Figures 5(e) and 5(f) present the PSNR and SSIM values, respectively, of the reconstructed images. The results indicate that regardless of the choice of image denoiser, the quality of reconstructed images decreases as the downsampling factor s increases. An interesting observation is that the reconstruction performance of PnP-DRUNet at s = 4 is even better than that of PnP-BM3D at s = 2. This indicates that the performance of the image denoiser itself has a greater impact on the reconstruction results than the effect of the downsampling factor s, emphasizing the importance of selecting a more potent image denoiser whenever feasible.

Next, we quantitatively evaluate the spatial resolution of the reconstructed images using the contrast of line-pairs. In the previous simulation experiment, the dimensions of the original HR image were fixed, and different recorded LR images were obtained by adjusting the downsampling factor s. However, in practice, only the LR image to be reconstructed is available, and the corresponding HR image is missing. This is equivalent to fixing the dimensions of the LR image and obtaining different HR images through different downsampling factors s. An original HR image containing four different widths (100, 150, 200, and 250 μm) of line-pairs is designed on a 50 μm grid, as shown in Fig. 6(a1). In this simulation experiment, we fix the pixel size of the recorded LR image at 200 μm, equivalent to a 4 × 4 downsampling of the blurred HR image. Varying the starting positions for downsampling will result in subpixel offsets, as shown in Fig. 6(a2). To demonstrate the impact of subpixel offset on SISR, four LR images are obtained with offsets of 0, 1/4, 2/4, and 3/4, and the addition of 1% Gaussian noise, as shown in Figs. 6(b1), 6(c1), 6(d1), and 6(e1), respectively. For the sake of brevity, only PnP-DRUNet is adopted in this simulation experiment, setting s = 2, 3, 4 to obtain the reconstructed HR images, and the estimated boundaries are cropped. The dimensions of the HR images obtained with different downsampling factors s are different, and PSNR and SSIM are unsuitable for quantitative comparison of the reconstruction results. Hence, we use the contrast of line-pairs to quantitatively evaluate image quality and resolution of the reconstructed images. Figures 6(b5)–6(b8), 6(c5)–6(c8), 6(d5)–6(d8), and 6(e5)–6(e8) respectively show the average intensity curves of the corresponding reconstructed images. It can be seen that all reconstructed images can almost perfectly reconstruct line-pairs with widths of 200 and 250 μm. Notably, except for the case of s = 2 with a subpixel offset of 2/4, all reconstructed images can achieve a resolution of 150 μm, breaking through the Nyquist sampling limit determined by the camera pixel size of 200 μm and achieving subpixel resolution.

It can be observed that different subpixel offsets and downsampling factors s significantly impact the resolution of the reconstructed images. Figure 6(f) shows the contrast of the 150 μm-wide line-pairs in the reconstruction results. With a subpixel offset of 0, the larger the downsampling factor s, the lower is the resolution, which is consistent with the phenomena shown in Figs. 5(e) and 5(f). However, with the subpixel offset, a noticeable decrease in resolution occurs when s = 2. This is because the pixel size of the reconstructed image with s = 2 is only 100 μm, which is relatively low compared with the 150 μm linewidth, making it highly susceptible to subpixel offset. The degree of contrast reduction of the 150 μm-wide line-pairs in the reconstructed images is somewhat alleviated in the cases of s = 3 and 4. This is because the increase in the downsampling factor s decreases the pixel size of the reconstructed image, which in turn reduces the impact of subpixel offset. However, for s = 4, there is no significant reduction in the pixel size of 50 μm in the reconstructed image compared with s = 3 with a 66.6 μm pixel size, while increasing s will increase the difficulty of reconstruction, as shown in Fig. 6(f). Therefore, for a fixed LR image, a smaller downsampling factor s is not always advantageous, which differs from the conclusion from Fig. 5. On the contrary, to fully leverage the subpixel resolution capability of the proposed algorithm, an appropriate downsampling factor s should be selected according to the target resolution.

We carried out real experiments of the gamma-ray imaging system at Northwest Institute of Nuclear Technology’s cobalt-60 radiation source platform, as shown in Fig. 7(a). The imaging system consists primarily of a cobalt-60 gamma-ray source, collimator, sample, LYSO scintillator, mirror, lens, CMOS camera, and shielding. Cobalt-60 emits gamma rays with an average energy of 1.2 MeV during its decay. Since the cobalt-60 is positioned several meters away from the sample and is small in size, the emitted gamma rays can be considered parallel to the optical axis. To verify the subpixel resolution capability of the proposed algorithm, we fabricated a sample made of tungsten, containing line-pairs with various widths, as shown in Fig. 7(b). The sample dimensions were 50 mm in height, 40 mm in width, and 30 mm in thickness. Gamma rays were incident perpendicular to the plane formed by the width and height of the sample. For 1.25 MeV gamma rays, a thickness of 30 mm tungsten attenuates the radiation intensity to 5%. We specifically focused on the line-pairs in the sample with widths of 250 μm (five pairs), 278 μm (six pairs), and 417 μm (six pairs). By adjusting the magnification of the lens, the CMOS camera was calibrated to achieve an equivalent pixel size of 300 μm at the sample plane, setting the image resolution of the gamma-ray imaging system to 300 μm. Consequently, regardless of the optical resolution of the imaging system, it is impossible to distinguish between line-pairs in the sample with widths of 250 and 278 μm solely on the basis of the image resolution. Therefore, employing the proposed SISR method could potentially achieve subpixel resolution, surpassing the limits imposed by the Nyquist sampling limit determined by the equivalent pixel size of the CMOS camera.

Figures 7(c)–7(e) show a raw image, a dark-field image, and a flat-field image, respectively. The raw image was obtained by exposing the sample to cobalt-60 radiation source for 100 s. The dark-field image refers to the image containing background noise of the gamma-ray imaging system, which was obtained with the Co-60 radiation source turned off and was also exposed for 100 s. Three dark-field images were acquired and averaged before being subtracted from the raw image. The flat-field image was acquired by exposing the imaging system to the cobalt-60 radiation source evenly after removing the sample, also for 100 s. Similarly, three flat-field images were captured and averaged during the experiment. From Fig. 7(e), it can be observed that the response of the gamma-ray imaging system exhibits spatial nonuniformity. Therefore, the image with dark field subtracted was divided by the flat-field image. Meanwhile, the corrected image was normalized to a range of 0–1, as shown in Fig. 7(f). An ROI image was selected containing line-pairs with widths of 250, 278, and 417 μm, and the enlarged pseudocolor ROI image is shown in Fig. 7(g). Figure 7(h) illustrates the contrast of the three types of line-pairs within the ROI image. It is observed that the 250 μm-wide and 278 μm-wide line-pairs exhibit high-frequency aliasing and are indistinguishable. The contrast of the 417 μm-wide line-pairs is noticeably reduced. The PSF of the imaging system is the foundation of SISR reconstruction. The knife-edge method64 was utilized to measure the PSF of the gamma-ray imaging system, and the result is shown in Fig. 7(i). The FWHM of the PSF is ∼360 μm, corresponding to an optical resolution of about 360 μm.

The ROI image was reconstructed using SISR with factor s = 2, 3, and 4. The first and second rows of Fig. 8 display the reconstructed images using PnP-BM3D and PnP-DRUNet, respectively, with periodic boundaries, as well as their corresponding average intensity curves of line-pairs. Since the periodic boundary assumption is incorrect for this ROI image, boundary artifacts appear in all reconstructed images. From the average intensity curves of the line-pairs in the reconstructed images, it is evident that the reconstruction results from both PnP-BM3D and PnP-DRUNet are worst when s = 2, best when s = 3, and intermediate when s = 4. This phenomenon can be explained by the conclusions obtained from Fig. 6. The third and fourth rows of Fig. 8 show the reconstructed images using the proposed boundary estimation method with PnP-BM3D and PnP-DRUNet, respectively, along with their corresponding average intensity curves of line-pairs. Different borders are used for different factors s, as marked in Figs. 8(c1)–8(d3). Compared with the results using periodic boundaries, those using estimated boundaries do not exhibit boundary artifacts, and the distinction of the line-pairs is significantly improved. This indicates that artifacts at the boundaries not only affect the boundary regions, but also impact the entire reconstructed image. Additionally, the estimated boundaries are shown outside the white dashed boxes, indicating that the boundaries can be approximately correctly estimated in all cases. Although there are some errors in the estimated boundaries, they do not affect the information within the white dashed boxes.

Figures 9(a)–9(c) show the average contrast of line-pairs with widths of 250, 278, and 417 μm, respectively, in the reconstructed images using different algorithms. From Fig. 9(a), it can be seen that for the 250 μm-wide line-pairs, the proposed PnP-DRUNet algorithm with estimated boundaries outperforms other algorithms across all factors s, particularly when s = 3. When using the same boundary conditions, for the 278 μm-wide and 417 μm-wide line-pairs, both PnP-BM3D and PnP-DRUNet can achieve high resolutions, as shown in Figs. 9(b) and 9(c). Furthermore, the method utilizing the estimated boundaries is significantly better than that using the periodic boundary assumption, indicating that the proposed boundary estimation method is crucial for images with nonzero boundaries. All algorithms achieve optimal results when s = 3, which represents a balance between sampling rate and reconstruction difficulty in SISR.

In this paper, we have proposed a novel SISR algorithm for gamma-ray imaging systems using a deep denoiser prior DRUNet based on the PnP framework, referred to as PnP-DRUNet. To address boundary artifacts resulting from incorrect assumptions about image boundaries, we have introduced a method to simultaneously estimate image boundaries during the reconstruction process. Both simulations and real experiments indicate that image boundaries have a significant impact on the reconstruction results. Employing the proposed PnP-DRUNet algorithm with estimated boundaries can effectively suppress boundary artifacts. Furthermore, the impact of the downsampling factor s on the reconstruction results has been analyzed. The results show that an appropriate downsampling factor s should be selected according to the target resolution. Notably, the proposed algorithm exhibits subpixel resolution capability, surpassing the Nyquist sampling limit determined by the pixel size of recorded images. Specifically, in simulation experiments, a resolution of 150 μm can be achieved for recorded images with a pixel size of 200 μm. In real experiments, a resolution of 250 μm is attained for recorded images with an equivalent pixel size of 300 μm of the CMOS camera.

Acknowledgment. This work was supported by the National Natural Science Foundation of China (Grant No. 12175183). The authors thank the staff of the cobalt-60 radiation source platform of the Northwest Institute of Nuclear Technology (Xi’an) for their support during the experiments.