Due to diffraction of light fields, optical imaging systems act as a low-pass filter, which blocks the high-frequency information of the object [1]. A computational super-resolution imaging technique aims to recover the high-frequency information of the object outside the cutoff frequency of the optical imaging system based on the theory of the analytic continuation and non-linear operator [2]. In a conventional imaging system, by assuming the spatial shift-invariant property of the point-spread function (PSF) of the imaging system, deconvolution techniques, such as Wiener filtering [3], Richardson-Lucy algorithm [4], and Tikhonov regularization algorithm [5], have been widely exploited as super-resolution techniques for their practical feasibility [6–8].

Different from conventional imaging, ghost imaging (GI) can nonlocally achieve the image of object by conducting a high-order correlation between the detected signals without spatial resolution (referred as the object arm) and a reference light field that does not pass through the object (referred as the reference arm) [9–12]. GI was originally developed by utilizing the quantum entangled property of two quantum entangled photons [13,14] and was later demonstrated to be realized with thermal light fields by exploiting the classical correlation property of the photons [15,16]. The key advantages of GI lie in that it decouples the detection and imaging processes and encodes the object with a fluctuation of light fields. The former makes GI more flexible than conventional imaging systems, and the latter can naturally combine compressive sensing (CS) [17,18] for image acquisition with high efficiency, and allow objects to be directly imaged in high-dimensional light-field space [19]. Owing to these advantages, GI has been widely applied in GI Lidar [20], neutron GI [21], Fourier-transform GI for X-ray [22], holographic GI [23], phase-contrast GI [24], and high-speed GI nanoscopy [25].

The diffraction-limited imaging resolution of GI is determined by the normalized second-order correlation function g(2) under an ensemble average [22–24,26–28], which is a counterpart to the PSF in conventional imaging systems [29]. Thus, the idea to enhance the imaging resolution of GI by using the deconvolution technique is reported [30]. However, the spatial shift-invariant property of g(2) in GI only holds under the ensemble average, which is not measurable, and it will be distorted within the limited number of samples [31]. Therefore, the classical deconvolution technique fails to improve the imaging resolution of GI, especially for a small number of samples.

Inspired by the prevalent CS theory [17,18], GI via sparsity constraints (GISC) can significantly reduce the number of samples to get the desired image [12]. Specifically, the image of the object in GISC is recovered by solving an ill-posed problem: minx‖x‖0,subjecttoy=Ax,(1)where y∈Rm is the sampling signals measured by a detector at the object arm, x∈Rn is the vectorized image, the sampling matrix A∈Rm×n (m<n) is made up of the reference arm light-field intensity, and ‖·‖0 denotes the ℓ0-norm. In CS, the mutual coherence of sampling matrix is vital to affect the image recovery [32], and smaller mutual coherence generally leads to better image quality. However, the mutual coherence of the sampling matrix is directly related to the g(2) in GI [19], and the mutual coherence is inevitably high when the spatial structures of the object beyond the resolution limited by g(2) need to be resolvable. Thus, for the high-resolution imaging application of GI, its sampling matrix will be a random high-coherence matrix within a limited number of samples, which is bad for CS to recover the image.

In the last decades, preconditioned methods have been proposed to reduce the mutual coherence of the random sampling matrix and improve the reconstruction quality for sparse signals in CS [33–37]. To be specific, a matrix P in the preconditioned method is multiplied to both sides of Eq. (1): minx‖x‖0,subjecttoPy=PAx.(2)By optimizing the matrix P, referred to as the preconditioner in the number linear algebra [38], the mutual coherence of preconditioned sampling matrix PA is smaller than that of sampling matrix A, which contributes to better image recovery through Eq. (2). However, those existing preconditioned methods are somewhat sensitive to the unavoidable detection noise, especially for high-coherence matrices. To mitigate the influence of the detection noise, a regularized method for preconditioning has been proposed [39,40]. Although the regularized preconditioned methods have been successfully applied to practical fluorescence diffuse optical tomography, their preconditioners are designed to be square matrices for the undersampling case (i.e., m<n), which makes their reconstructed performances highly rely on recovery algorithms to solve Eq. (2) and be sensitive to the unknown regularization parameter as well.

In this paper, aiming at the noisy environment in practical applications, we propose a preconditioned deconvolution method to improve the imaging resolution of GI. The theoretical analysis shows that the preconditioned deconvolution method is an extension of the deconvolution technique to GI and reduces to the classical deconvolution technique for conventional imaging systems. The imaging resolution of GI after preconditioning is restricted to the detection noise, which is consistent with the classical deconvolution technique. The effectiveness of the preconditioned deconvolution method is demonstrated on the GISC camera, which is a typical paradigm of the GI technique. Specifically, both simulations and experiments on the GISC camera system display that the imaging resolution and image quality are obviously enhanced by using the preconditioned deconvolution method. Moreover, the experimental results exhibit that 1.4-fold resolution enhancement over Rayleigh criterion of the conventional imaging system is achieved via the preconditioned deconvolution method, which could promote GI in the application of high spatial-resolution imaging scenarios. We would like to mention that the preconditioned deconvolution method utilizes prior imaging systems to overcome the diffraction limit, which is different from those super-resolution imaging methods that mainly focus on the use of sparsity priors of imaging objects in GI [25,41,42], and it can also combine prior imaging objects naturally to improve the imaging resolution furthermore.

To alleviate the influence of detection noise, the objective function of the proposed preconditioned method to solve the preconditioner P is as follows: minP‖PA−I‖F2+γ‖Py−PAx‖22,(3)where ‖·‖F denotes the Frobenius norm, and γ is a weight factor. Different from those existing preconditioned methods, whose objective function aims to make the Gram matrix of matrix PA approach an identity matrix, the first term of Eq. (3) is to let the matrix PA approach an identity matrix to refine the mutual coherence of preconditioned sampling matrix PA, and the second term of Eq. (3) aims to restrict the amplification of the detection noise after preconditioning.

Let the derivative of Eq. (3) on matrix P be 0, and one can find a closed-form solution of Eq. (3): P≔AT(AAT+γσ2I)−1.(4)Specifically, the derivatives of two terms in Eq. (3) are, respectively, ∂∂P‖PA−I‖F2=2(PA−I)AT,(5)and ∂∂P‖Py−PAx‖22=2σ2P,(6)where the detection noise v (namely, y−Ax) is assumed as an additive white Gaussian noise (AWGN), and each element i.i.d. obeys the uncorrelated Gaussian distribution N(0,σ2).

Herein we introduce two important lemmas, which will be useful for our theoretical analysis.

Lemma 1 (Woodbury Formula [43]) The following equality holds if A−1 and B−1 exist: (A+CBCT)−1=A−1−A−1C(B−1+CTA−1C)−1CTA−1.(7)Lemma 2 (Approximations [44]) The following equality holds if A is large and symmetric: (I+A)−1=I−A+A2−A3+⋯.(8)Using Eq. (4), the preconditioned sampling matrix PA can be given as follows: $$\begin{gathered} PA=AT(AAT+γσ2I)−1A \\ =Eq.(7)[ATAγσ2−(1γσ2)2ATA(I+ATAγσ2)−1ATA] \\ =Eq.(8)ATA(ATA+γσ2I)−1· \end{gathered}$$(9)In GI, the imaging resolution is characterized by the normalized second-order correlation function g(2), which plays the role of the PSF in the conventional imaging system [29–31]. Thus, the matrix ATA approaches a block circulant matrix for the large number of samples m, where the matrix ATA is the matrix representation of the g(2). Thus, ATA can be decomposed as [45] ATA=FTΛ1F,(10)where the matrices F and FT are associated with the Fourier and inverse Fourier transforms (satisfying FFT=FTF=I), and Λ1=diag[α1,…,αn] is a diagonal matrix, whose diagonal elements are the Fourier coefficients of the first column of the block circulant matrix ATA.

Substituting Eq. (10) into Eq. (9), the preconditioned sampling matrix PA can be given by $$\begin{gathered} PA=FTΛ1/(Λ1+γσ2I)F \\ =FTdiag[α1α1+γσ2,…,αiαi+γσ2,⋯]F \\ =FTΛ2F, \end{gathered}$$(11)where the diagonal matrix Λ2=diag[α1α1+γσ2,…,αiαi+γσ2,⋯].

Equation (11) shows that the preconditioned sampling matrix PA is also a block circulant matrix, and the frequency coefficient αi/(αi+γσ2) of the matrix PA approaches one as the detection noise reaches 0, which implies that the improvement of spatial resolution of matrix PA is limited to the detection noise. One can interpret from Eq. (11) that directly multiplying the preconditioner P to the sampling signals y formally is the same as the deconvolution technique in the classical convolution model [46,47]. The difference is that the sampling matrix A in GI is a random matrix, whereas that in the classical convolution model is a block circulant matrix. Specifically, for the spatial shift-invariant imaging system, which can be described as a classical convolution model, its sampling matrix A is a block circulant matrix. Thus, A can be similarly decomposed as follows: A=FTΛ1′F,(12)where Λ1′=diag[α1′,…,αn′] is a diagonal matrix whose diagonal elements are the Fourier coefficients of the PSF of the conventional imaging system. Correspondingly, Eq. (9) becomes $$\begin{gathered} PA=FTΛ1′2/(Λ1′2+γσ2I)F \\ =FTdiag[α1′2α1′2+γσ2,…,αi′2αi′2+γσ2,⋯]F, \end{gathered}$$(13)which regresses into the classical deconvolution technique [47].

We would like to mention that the preconditioned method for the noiseless case, i.e., σ2=0, has the same form as the pseudo-inverse GI method [48,49]. Unfortunately, the performances of the pseudo-inverse GI method are sensitive to the detection noise, which hampers their practical applications.

As a typical paradigm of the GI technique, the GISC camera can simultaneously achieve an object’s high-dimensional (e.g., spatial, spectral, polarized) information with one snapshot by encoding the high-dimensional information of light fields irradiating from the object into speckle patterns on a detectable two-dimensional (2D) plane [27]. Figure 1 shows the setup of the GISC camera, which applies a spatial random phase modulator (SRPM) after the image plane of conventional imaging systems, such as a microscope and telescope, to modulate the light field irradiating from the object into the speckle field that is detected by detector 2. In the GISC camera, each pixel of detector 2 acts like a bucket detector or single-pixel detector and contains information from all pixels in the image, which is the main feature of the object beam in GI. Correspondingly, the image of the GISC camera can be reconstructed by conducting a high-order correlation between the detected signals and the predetermined reference light-field patterns under the spatial ensemble average. For GI calculating the ensemble average in the time domain, a bucket detector or single-pixel detector is adopted to measure the modulated signal irradiating from the objects, where the fluctuation of the illuminating light fields needs to be changed per measurement. Because a large number of measurements are needed to obtain the desired image, the GI ensemble averaged in the time domain generally demands long image acquisition time. By converting the calculation of the ensemble average in the time domain to the spatial domain, the GISC camera achieves the image in one snapshot, which overcomes the shortage of long image acquisition time and makes it possible to observe high-speed dynamic processes. Although the GISC camera needs a certain detection signal-to-noise ratio (SNR) to guarantee image recovery, its image acquisition is fast because only one snapshot is required to achieve the image.

In this paper, we focus on improving the spatial resolution of the GISC camera. By calculating the second-order correlation function between the intensity fluctuations at the predetermined reference arm and the object arm, the spatial image of object can be achieved [27]: $$\begin{gathered} ΔG(2,2)(xr)=⟨ΔIt(xt)ΔIr(xt;xr)⟩xt \\ ∝∫T(xr′)g(2)(xr,xr′)dxr′, \end{gathered}$$(14)where T(xr) represents the spatial image of object, It(xt) is the speckle image, and Ir(xt;xr) denotes the speckle pattern of the GISC camera corresponding to one monochrome point light source at the position xr. xr and xt are 2D vectors on the object plane and the detection plane of the GISC camera, respectively, and ⟨·⟩xt means the spatial ensemble average over the coordinate xt. The normalized second-order correlation function g(2)(xr,xr′) describes the spatial resolution of the GISC camera.

The spatial image of object can be also retrieved by solving an inverse problem: minx‖x‖0subjecttoy−y¯=(A−A¯)x,(15)where the sampling matrix A is made up of speckle patterns Ir(xt;xr), the sampling signal y consists of the speckle image It(xt), x represents the spatial image, and ·¯ denotes an average operator to calculate the mean value of the columns of matrix A and vector y. For notational simplicity, denote Φ=A−A¯. As mentioned in Ref. [19], the mutual coherence of matrix Φ is equal to the maximum value of g(2)(xr,xr′), which represents the correlation of speckle patterns corresponding to the two nearest points on the object plane. Thus, for the high-resolution imaging application, the mutual coherence of matrix Φ is naturally large, which is not good for high-resolution image recovery.

Figure 1 shows the experimental setup of the GISC camera. To be specific, an object is projected on the image plane through the conventional imaging system, which is made up of a filter (FL 532-3, Thorlabs), lens 1 (AF70-300 mm, TAmRon) with focus length 300 mm, lens 2 (AF70-300 mm, TAmRon) with focus length 300 mm, and an iris with diameter 25 mm. The image of the object is modulated into a speckle image by an SRPM (DGUV 10-1500, Thorlabs), and is directly recorded on detector 2 (iKon–M, Andor, pixel size of 13 μm). The SRPM is placed after the image plane of the conventional imaging system with distance 9.7 mm, and is set before detector 2 with distance 40 mm. For comparative purposes, the ground-truth image of the object through the conventional imaging system with a large aperture diameter is recorded on detector 1 (Stingray F–504B, Allied Vision) at another optical path split by a 10/90 (R:T) beam splitter (BS043, Thorlabs). The GISC camera system was built on an indoor optical platform, where the objects were illuminated by a white light source through the arc lamp (Arc lamp 66902, Newport).

Before imaging, we need to calibrate the GISC camera for constructing the sampling matrix A. In the calibration process, a fiber (core diameter 16 μm) coupling with incident light (λ=532nm) acts as a point source, which is mounted on an electric translation stage. A series of speckle images are generated by shifting the point source at a step size of 3 μm in the whole 400μm×400μm field of view (FoV). Then, the sampling matrix A is built in the following way, namely, each column of the sampling matrix consists of speckle images by reshaping a 2D speckle image into a column vector. In the imaging process, we just replace the point source with the imaging object and record the corresponding speckle image on detector 2 with one snapshot, which makes up the sampling signal y in the same processing manner as that in the calibration process.

In order to validate the proposed preconditioned deconvolution method for improving the reconstruction performance of the GISC camera, we carry out numerical simulations. In the simulations, we adopt the real speckle patterns obtained in the calibration process, and the sampling signal y is generated via y=Ax+v. The images of object x are constructed as shown in Fig. 2. Herein the detection noise v is assumed as an AWGN, and the noise level is measured by the detection SNR, which is defined as SNR=10log10(y¯/σ2) (dB), where y¯ is the mean values of the sampling signal y, and σ2 is the variance of the distribution of AWGN.

The first simulation aims to show how effective the proposed preconditioned deconvolution method is in reducing the mutual coherence of sampling matrices. To this end, we compare the mutual coherence of the sampling matrix before and after the proposed preconditioned deconvolution method, as well as the normalized mutual correlation function. As observed in Figs. 2(a) and 2(b), the mutual coherence of preconditioned matrix PΦ is uniformly smaller than that of original matrix Φ among the tested sampling rates η=m/n, where m is the number of samples and n is the number of pixels of the vectorized image x, and the width of the normalized mutual correlation curve of preconditioned matrix PΦ is sharper than that of original matrix Φ, which implies that the spatial resolution could be enhanced after preconditioning. One can see that the mutual coherence of Φ approaches a lower bound as the sampling rate increases, which is fundamentally restricted to the resolution limit by the g(2) in GI. Overall, the effectiveness of the proposed preconditioned deconvolution method to refine the sampling matrix is clearly demonstrated in the simulation.

Next, we empirically compare the image reconstruction of the preconditioned deconvolution method with other approaches. For comparative purposes, the representative approaches including differential GI (DGI) [50], two-step iterative shrinkage/thresholding (TwIST) for sparse signal recovery [51], the proposed preconditioned GI (PreGI, namely, the proposed preconditioner P directly multiplies the sampling signal y), and the proposed preconditioned deconvolution method plus TwIST [PreTwIST, namely, solving Eq. (2) with the TwIST from the preconditioned sampling signals Py] algorithm are adopted. Figure 2(c) represents the imaging results of the resolution target and cell image by DGI, PreGI, TwIST, and PreTwIST algorithms under the detection SNR 25 dB at a sampling rate of 0.4, and Figs. 2(d) and 2(e) exhibit the comparison of image quality under different sampling rates by adopting the quantitative evaluation metric: peak SNR (PSNR) [52]. Clearly, the image quality of the GISC camera is greatly improved after the proposed preconditioned operator.

To investigate the improvement of image quality by the preconditioned deconvolution method in the noisy environment, we perform simulations and compare the recovery performance by the abovementioned algorithms under different noisy levels. In the simulation, the detection SNRs are set from 15 dB to 45 dB with increasing step 5 dB, and the sampling rate is η=0.4. Figure 3(a) presents the normalized mutual correlation function of the preconditioned matrix PΦ at different detection SNRs. One can observe that the width of the mutual correlation curve gets small when the detection SNR increases, which matches our theoretical analysis results well. Moreover, both Figs. 3(b)–3(d) show that the image quality of the GISC camera is indeed improved by the proposed preconditioned deconvolution method among the tested detection SNRs, and the improvement gets big as the detection SNR increases. Overall, the simulation results demonstrate the effectiveness of the preconditioned deconvolution method in refining the sampling matrix and improving the image quality of the GISC camera.

In this section, we apply the proposed preconditioned deconvolution method to dealing with the experimental data of the GISC camera. In the experiment, we image five objects, including the transmission objects and reflective objects. Figure 4 shows the imaging results of different objects by the DGI, PreGI, TwIST, and PreTwIST algorithms. The left column in Fig. 4 represents the ground-truth images of objects, which are achieved through a conventional imaging system as shown in Fig. 1 with an aperture diameter 25 mm of the iris. It can be observed that the PreTwIST algorithm achieves the best results compared with the other approaches, especially in clearing the background noise of the image. Moreover, as shown in Fig. 4(d), the preconditioned deconvolution method can obviously improve the spatial resolution of the image compared with the DGI method.

We also test the recovery performance by using the preconditioned deconvolution method under different sampling rates and different detection SNRs. In the experiment, the detection noise is mainly limited to the shot noise, and the detection SNR is measured by SNR=10log10N (dB), where N denotes the average photon number which can be estimated from the number of electronics recorded on detector 2. Figure 5 shows the PSNR of the reconstructed image of the transmission object resolution target under different sampling rates and different detection SNRs. It can be observed that the reconstruction quality is uniformly improved after the preconditioned deconvolution method among the tested sampling rates and detection SNRs.

In this section, we extend the preconditioned deconvolution method to super-resolution imaging. Following the requirement of Ref. [19], where the spatial resolution of the modulation module is much higher than that of the conventional imaging module in the GISC camera, we adjust the parameters of the GISC camera to make the spatial resolution of the GISC camera approach the Rayleigh limit of the conventional imaging system. To be specific, the diameter of aperture is set as 10 mm, the distance between the first image plane and the SRPM is changed into 2.7 mm, and the distance between the SRPM and detector 2 is 15 mm. Figure 6 reveals that the normalized mutual correlation function of the GISC camera approaches the theoretical Rayleigh limit of the conventional imaging system at an incident wavelength of 532 nm. The Rayleigh limit of the conventional imaging system in its image plane before the SRPM of the GISC camera is about 19.5 μm.

Figure 7 shows reconstructed images of the resolution target by the aforementioned four algorithms. Compared with the imaging results of DGI and TwIST algorithms, the proposed PreGI algorithm and PreTwIST algorithm achieve much higher imaging resolution. Particularly, the resolution below 13.8 μm can be obviously recognized by PreGI and PreTwIST algorithms. Considering the Rayleigh criterion 19.5 μm, both the imaging results of PreGI and PreTwIST algorithms achieve a 1.4-fold resolution enhancement over the Rayleigh criterion in this experiment. Moreover, as shown in Figs. 7(d), 7(f), and 7(g), the imaging result of the PreTwIST algorithm has a much clearer background and higher resolution than that of PreGI. Also, a much complicated object butterfly is tested in the experiment. Figure 8 presents images of a butterfly by the aforementioned algorithms. It shows that the PreGI and PreTwIST algorithms achieve higher imaging resolution than DGI and TwIST algorithms. Overall, the proposed preconditioned deconvolution method exhibits competitive advantages in super-resolution imaging.

For conventional imaging systems, by exploiting the sparsity-imposing priors of imaging objects and the prior information of imaging systems with the assumption of the spatial shift-invariant property of the PSF, super-resolution imaging methods via sparsity constraints, which are based on the sparsity priors of imaging objects, and the deconvolution techniques that use the prior of the PSF of imaging systems have been widely studied [7,8,53–55]. However, existing GI techniques to improve the imaging resolution mainly focus on the use of sparsity priors of imaging objects in reconstruction [25,41,42], and rarely take prior imaging systems into account, because the spatial shift-invariant property of g(2) is distorted within the limited number of samples. In this work, with the proposed preconditioned deconvolution method based on prior GI systems, the classical deconvolution method is extended to GI systems for improving the imaging resolution. The simulation and experimental results demonstrate that both the methods using the sparsity prior of imaging objects and the proposed preconditioned deconvolution method applying the prior of imaging systems can effectively improve the imaging resolution of GI. Moreover, the combination of the sparsity prior of imaging objects and the prior of imaging systems achieves the best performance in the improvement of GI’s imaging resolution.

For linear imaging systems, where the imaging model of those systems can be expressed as y=Ax and the sampling matrix A can be structured matrices or random matrices, such as the fluorescent microscopy, telescope, aperture coding imaging system, computational GI, and single-pixel imaging, the preconditioned deconvolution method offers new possibilities to improve the imaging performance of all those systems. Although the computational cost of the inverse operator for a large-scale matrix in the proposed preconditioned deconvolution method is inevitably high, the preconditioner P can be obtained in advance through highly effective computation servers. Moreover, the computational cost for the inverse operator is expected to decrease by dividing a large-scale matrix into several small-scale matrices and further solving a partitioned matrix [56].

Existing preconditioned methods achieving the preconditioner are to make the Gram matrix of the preconditioned sampling matrix PA approach an identity matrix [33–37,39,40]. Those Gram-matrix methods are mathematically strict because the mutual coherence of a preconditioned sampling matrix equals the largest non-diagonal element of its Gram matrix. However, the objective function of those Gram-matrix methods is a non-convex function, and the optimal preconditioner P is obtained in an iterative fashion. Different from those Gram-matrix preconditioned methods, the proposed preconditioned method is to let the preconditioned sampling matrix PA directly approach an identity matrix. Although it is an indirect way in mathematics to make its Gram matrix reach the identity matrix, thereby reducing the mutual coherence, the objective function of the proposed method is a convex function, which exists within the optimal closed-form solution. Moreover, the idea to make the preconditioned sampling matrix PA close to the identity matrix can reduce the interference between the adjacent pixels of image x that pass through the imaginary system corresponding to the preconditioned sampling matrix PA. Therefore, combined with the noise constraint, the objective function of the proposed preconditioned method has clear physical interpretations.

In conclusion, based on prior imaging systems, we presented a preconditioned deconvolution method to improve the spatial resolution of GI by refining the mutual coherence of the sampling matrix corresponding to the imaging system. Our theoretical analysis shows that the spatial resolution of GI after preconditioning is still restricted to the detection noise, and the preconditioning method actually is an extension of the classical deconvolution technique of GI. As a typical paradigm of the GI technique, the GISC camera was implemented to validate the effectiveness of the proposed preconditioned deconvolution method. Both simulation and experimental results demonstrate that the spatial resolution and image quality of the GISC camera are greatly enhanced by using the preconditioned deconvolution method. Owing to the high efficiency of information acquisition, the GISC camera has been successfully applied to fluorescent microscopy [25] and remote sensing [27,57]. Therefore, in conjunction with the GISC camera, the preconditioned deconvolution method is expected to promote the applications of GI for high spatial and temporal resolution imaging scenarios, such as live cells imaging and remote sensing.