Ghost imaging (GI) is an imaging technique that may reconstruct the object nonlocally by correlating the light intensity passing through or reflected by the object and the illuminating light field^[1–7]. GI has been already demonstrated with entangled photons and thermal light^[8–10]. It has become an interesting research area in the last 10 years since it has many potential applications such as remote sensing^[11–13], imaging through scattering media^[14], object authentication^[15] and object tracking^[16], imaging with a colored object^[17], and three-dimensional imaging^[18].

One important topic in GI is how to improve sampling efficiency and quality of the reconstructed image. Different methods have been proposed, including homodyne detection GI^[19–21], high-order correlation GI^[22–24], imaging with averaged speckle patterns^[25], and compressive ghost imaging (CGI) which is done by introducing the compressive sense (CS) techniques into the GI^[26]. Under the framework of ghost imaging via sparsity constraints (GISC)^[11], the influence of sparsity property of images on GI with thermal light has been studied^[27]. A speckle light modulator (SLM) working in the amplitude regime is used both to display the object and the intensity distribution that can recover a good-quality image^[28]. Also the use of measurement matrices with multi-correlation scales was proposed and proven to be more efficient than those with a uniform correlation length^[29].

In the theory of compressive sensing, if the numbers of sampling needed for the reconstruction of the signal with overwhelming probability is $m$, then $m\ge \mathrm{C}·{\mu }^2(\mathbf{A},\mathbf{\Psi })·S·\log n$^[30], where C is a constant, $S$ is the sparse level of the signal, and $n$ is the length. Term $\mu$ is the coherence between the measuring matrix $\mathbf{A}$ and dictionary $\mathbf{\Psi }$. It is obvious that if we reduce the coherence, fewer samplings would be needed for the reconstruction. The quality of the reconstructed image is better under the same sampling. In this Letter, the compressed sensing algorithm is used for our reconstruction of the image. Inspired by the CS theory that the quality of the reconstructed result is better if the coherence between the measuring matrix and sparsifying dictionary matrix is smaller^[31–34], we show that the speckle patterns in GI can be optimized with respect to the sparsifying dictionary by minimizing the mutual coherence between them. In the rest of the Letter, we will first describe the experimental setup and the theory of speckles optimization, followed by the simulation and experimental results, and last the conclusion.

Figure 1 shows the experimental setup of GISC using a digital micromirror device (DMD). The light from the halogen lamp is focused by a lens and then modulated by the DMD. The resulting speckle pattern is then projected onto the object plane by lens ${\mathrm{L}}_1$. The light transmitting from the object is focused to the single-pixel detector by lens ${\mathrm{L}}_2$. Each speckle pattern ($m\mathrm{pixels}\times n\mathrm{pixels}$) can be reshaped to one row ($1\times N$, $N=m\times n$) of a measuring matrix denoted by $\mathbf{A}$. At the same time a light intensity $y$ is recorded by the single-pixel detector. Repeating the process $K$ times, a measuring matrix $\mathbf{A}$ of $K\times N$ elements and a light intensity vector denoted by $\mathbf{Y}$ of length $K$ can be obtained. The number of samplings can be far smaller than the pixels of the image. Denoting the unknown object by a column vector $\mathbf{X}$ of length $N$, the sampling process can be described as $$\mathbf{Y}=\mathbf{AX}+\mathbf{\sigma },$$(1)where $\mathbf{\sigma }$ is the measurement noise. According to CS theory, it is desirable to sparsely represent the vector $\mathbf{X}$ to obtain a better reconstruction, which is usually feasible for a natural image^[26,27]. Suppose we have the suitable dictionary $\mathbf{\Psi }$ (e.g., a wavelet dictionary) and $\mathbf{X}$ can be sparsely represented as $\mathbf{X}=\mathbf{\Psi }\mathit{\theta }$, where $\mathit{\theta }$ is the sparse representation vector; with the data $\mathbf{Y}$, the reconstruction of $\mathbf{X}$ can been regarded as the following optimization problem $$\tilde{\mathbf{A}}=\arg \underset{\mathbf{X}}{\min }\frac{1}{2}{\Vert \mathbf{Y}-\mathbf{AX}\Vert }_2^2+\tau {\Vert \mathit{\theta }\Vert }_1,\text{s.t.}\mathbf{X}=\mathbf{\Psi }\mathit{\theta },$$(2)where $\tau$ is a nonnegative parameter and its value should be chosen according to the sparsity of the object. The sparser the object’s representation, the larger $\tau$ should be. Term ${\Vert \mathbf{v}\Vert }_1$ is the taxicab norm and ${\Vert \mathbf{v}\Vert }_2$ the Euclidean norm.

Besides the general convex optimization algorithms^[35], there are many other methods specifically developed to solve the problem described by Eq. (2), such as greedy algorithms^[36], iterative thresholding algorithm^[37], fast iterative shrinkage-thresholding algorithm^[38], and the gradient projection for sparse reconstruction (GPSR) algorithm^[39]. In our work, we adopted the widely used GPSR algorithm.

As pointed out in Ref. [40], the coherence between the Matrix $\mathbf{A}$ and $\mathbf{\Psi }$ has a significant impact on the quality of the reconstruction. The coherence of $\mathbf{A}$ and $\mathbf{\Psi }$ can be characterized by the absolute value of the inner products of the direction vectors of the columns of $\mathbf{D}=\mathbf{A}\mathbf{\Psi }$, i.e. $${\mu }_{ij}=\frac{|{\mathbf{d}}_i^{\mathrm{T}}{\mathbf{d}}_j|}{{\parallel {\mathbf{d}}_i\parallel }_2{\parallel {\mathbf{d}}_j\parallel }_2},1<i,j<N,i\ne j,$$(3)where ${\mathbf{d}}_i$ is the $i$th column of $\mathbf{D}$. The value of ${\mu }_{ij}$ is confined between [0,1]. Approximately speaking, the smaller the overall value of ${\mu }_{ij}$, the better reconstruction one may expect^[31]. Among the different definitions of the overall mutual coherence value, we choose the sum of ${\mu }_{ij}^2$ to be minimized for its ease of computation $$\mu =\sum _{i\ne j}{\mu }_{ij}^2.$$(4)

To find the (approximate) solution of $\mathbf{A}$ that minimizes $\mu$, first one may define $\mathbf{G}={\mathbf{D}}^{\mathrm{T}}\mathbf{D}$, $\mathbf{E}={\Vert \mathbf{G}-\mathbf{I}\Vert }_2^2$ and then solve the following optimization problem $${\mathbf{A}}_{\text{opt}}=\arg \underset{\mathbf{A}}{\min }\mathbf{E},\text{s.t.}\mathbf{A}\ge 0,$$(5)where $\mathbf{A}\ge 0$ means that the elements of $\mathbf{A}$ must be nonnegative since they are the light intensity values. Generally the quantity $\mathbf{E}$ cannot be exactly regarded as $\mu$. However, as $\mathbf{E}$ is decreased, the diagonal elements of $\mathbf{G}$, which are ${\Vert {\mathbf{d}}_j\Vert }_2^2$, are gradually pushed to 1 and then the off-diagonal elements of Matrix $\mathbf{G}$ gradually approach ${\mu }_{ij}$, and thus $\mathbf{E}\approx \mu$. The problem described by Eq. (5) can be solved using the gradient descent algorithm^[29](1)Set the initial matrix ${\mathbf{A}}_{\text{init}}$ to be a Gaussian matrix.(2)Repeat the following iteration until convergence is achieved and the optimized matrix ${\mathbf{A}}_{\text{opt}}$ is obtained •Calculate the gradient $\partial \mathbf{E}/\partial \mathbf{A}$;•Update the measuring matrix $\mathbf{A}$ with $\mathbf{A}\leftarrow \mathbf{A}-\lambda \partial \mathbf{E}/\partial \mathbf{A}$, with $\lambda$ the fixed stepsize;•Update the measuring matrix by imposing the constraints $\mathbf{A}\ge 0$.

Note that to work with the DMD, first we need to rescale the optimized Matrix ${\mathbf{A}}_{\text{opt}}$ to match the dynamic range of the DMD and then round down the elements to get a matrix with integer elements only so that the speckle patterns can be used on the DMD.

The proposed method is tested by simulation first. The Haar wavelet with decomposition level of 6 has been chosen as our dictionary Matrix $\mathbf{\Psi }$. Following the procedure described previously and starting with ${\mathbf{A}}_{\text{init}}$, which is Gaussian-distributed, we obtain ${\mathbf{A}}_{\text{opt}}$ within 50 iterations. Figure 2 shows the distributions of ${\mu }_{ij}$ of the original matrix ${\mathbf{D}}_{\text{init}}={\mathbf{A}}_{\text{init}}\mathbf{\Psi }$ (blue line) and the optimized matrix ${\mathbf{D}}_{\text{opt}}={\mathbf{A}}_{\text{opt}}\mathbf{\Psi }$ (red line). It can be seen that the coherence ${\mu }_{ij}$ of ${\mathbf{D}}_{\text{opt}}$ are concentrated in the region with smaller coherence values if compared with that of ${\mathbf{D}}_{\text{init}}$. We further calculated the averages of the coherence ${\mu }_{ij}$ of ${\mathbf{D}}_{\text{init}}$ and ${\mathbf{D}}_{\text{opt}}$, and the results show that the average is reduced from 0.0126 to 0.0096 by optimization.

With the two Matrices ${\mathbf{A}}_{\text{init}}$ and ${\mathbf{A}}_{\text{opt}}$, we simulated the imaging process using the object (pattern of “tai chi”) with a size of $64\mathrm{pixels}\times 64\mathrm{pixels}$ as shown in Fig. 3(a). The gray level is 8. Figure 3(b)–3(f) are the results with different numbers of samplings (500, 1000, 1500, 2000, and 2500). Shown in the top line are the results obtained with Gaussian Matrix ${\mathbf{A}}_{\text{init}}$ and in the bottom line the optimized Matrix ${\mathbf{A}}_{\text{opt}}$ from which we can see obvious improvements.

To quantitatively evaluate the qualities of reconstructed images, the relative mean square error (RMSE) and peak signal-to-noise ratio (PSNR) is calculated $$\mathrm{RMSE}={\sum }_{ij}{(I_{ij}^{\prime }-I_{ij})}^2/{\sum }_{ij}{I}_{ij}^2,$$(6)where $I_{ij}$ is the gray level of the $(i,j)$th pixel of the original image and $I_{ij}^{\prime }$ is the reconstructed image $$\mathrm{MSE}=\frac{1}{mn}\sum _{i=0}^{m-1}\sum _{j=0}^{n-1}{\Vert I^{\prime }(i,j)-I(i,j)\Vert }^2\mathrm{PSNR}=10·{\log }_{10}\left(\frac{{\mathrm{MAX}}_I^2}{\mathrm{MSE}}\right),$$(7)where ${\mathrm{MAX}}_I=255$ is the gray level of the image.

The RMSE and PSNR values of the reconstructed images with different numbers of samplings by ${\mathbf{A}}_{\text{init}}$ and ${\mathbf{A}}_{\text{opt}}$ are shown in Fig. 4, from which we can see that for the same numbers of samplings, the RMSE value obtained with ${\mathbf{A}}_{\text{opt}}$ is reduced by more than 20% compared with the corresponding ${\mathbf{A}}_{\text{init}}$. The PSNR value obtained with ${\mathbf{A}}_{\text{opt}}$ is also larger than the corresponding ${\mathbf{A}}_{\text{init}}$.

We then test the method with experiments. Figure 5(a) shows the image of the test object. The image is obtained by the same system operating in a pixel-by-pixel scanning mode and thus may serve as the reference image. Figures 5(b)–5(e) show the results reconstructed from different number of samplings (1000, 2000, 3000, and 4000), from which we may see that the result obtained with 2000 samplings using ${\mathbf{A}}_{\text{opt}}$ is of comparable quality to that obtained with 4000 samplings using ${\mathbf{A}}_{\text{init}}$. The RMSE and PSNR values of the results are calculated using Fig. 5(a) as the reference image and are shown in Fig. 6, which confirms the observed improvements in Fig. 5.

In conclusion, we show that the speckle patterns in GISC can be optimized with respect to the sparsifying dictionary by minimizing the mutual coherence between the measuring matrix and the dictionary matrix. The results show that the quality of the reconstructed images can be much improved using the optimized matrix. We expect this method, when combined with, e.g., wavefront shaping techniques^[41], can be used to design GI systems with high performance.