1. Introduction

Single-pixel imaging is an imaging technology that has gained much attention in recent years. In 1995, Y.H. Shi et al. used entangled photon pairs for the first time to achieve ghost imaging [1]. In 2002, classical source was used to implement ghost imaging experiments, removing the limitations on light sources [2]. In 2008, Sharpio proposed a computational ghost imaging system (CGI) [3]. CGI uses a spatial light modulator to replace the reference light path in the ghost imaging system, which significantly reduces the complexity of the system. CGI is also called single-pixel imaging (SPI).

In the field of special waveband imaging, Padgett et al. have applied SPI to the detection of methane leaks [4]. Furthermore, researchers have developed the SPI system in the THz band [5] and recently a real-time Thz SPI setup has been demonstrated [6]. SPI has also been widely used in three-dimensional (3D) imaging. For the targets at a distance, SPI is combined with the technique of time-of-flight (ToF) to form a non-scanning single-pixel lidar [7,8]. For the daily Lambeau targets, SPI reconstructs the shading images to form the 3D shape of the target based on photometric stereo vision [9,10].For the case where the target size is small and the depth accuracy is high, the defocus characteristics of SPI can be used to estimate the target depth [11]. Additionally, SPI is also applied in X-ray imaging [12], infrared imaging [13], object tracking [14], ultrafast imaging [15], and Non-line-of-sight imaging [16].

SPI usually uses random matrix, or Hadamard matrix modulation as the measurement matrix. For the traditional imaging, we must obtain the image first, and then extract its features. However, in SPI the features (such as edge images [17], wavelet coefficients [18]) can be directly obtained by changing the measurement matrix, without the need to reconstruct the entire image. Obviously, the frequency spectrum is also an important feature for an image. In 2015, Zhang et al. proposed the technique of Fourier single-pixel imaging (FSI) [19]. FSI projects a series of patterns of Fourier basis onto the target, obtains the corresponding Fourier coefficients through the N-step phase shift method, and reconstructs the target image. FSI can directly obtain the Fourier coefficients of interest and realize image reconstruction with low sampling ratio. Additionally, compared with the SPI based on Hadamard matrix, FSI has higher imaging efficiency [20]. Since its introduction, FSI has expanded to multi-mode imaging [21], microscopic imaging [22] and other fields [23,24].

In 2019, Meng et al. proposed a method called sparse Fourier single-pixel imaging (S-FSI) [25], where the sampling density is inversely proportional to the distance between the sampling point and the center of the Fourier domain. Clearly, it is a variable density sampling method. Compared with FSI, the formed image using S-FSI has a significant improvement in the signal-to-noise ratio (SNR) and the ability to retain details of objects, when the algorithm of compressed sensing(CS) is used for reconstruction. However, it takes a long time to reconstruct an image by using the CS algorithm, especially for high-resolution images. If the inverse Fourier transform (IFT), which is much faster than CS, is applied, then the formed image of S-FSI will suffer significant noise.

To address this problem, in this paper we propose adaptive Fourier single pixel imaging (A-FSI) to solve the trade-off between the quality of the formed image and the reconstruction time. We first study the radial correlation between low-frequency and high-frequency components in the Fourier domain statistically. Secondly, we estimate the possible locations of significant high-frequency components based on the Fourier coefficients of low-frequency. Finally, we measure the coefficients of the estimated Fourier components and form the images. Simulations and experiments show that compared with FSI and S-FSI, the formed image of A-FSI has higher quality by using the IFT.

2. Methodology

In this section, we first introduce the principle of FSI and S-FSI. We then analyze the radial correlation between the significant components in the low-frequency and high-frequency regions statistically. Finally we introduce A-FSI and show the evaluation methods. Note that for FSI, S-FSI, and A-FSI, the IFT is used for image reconstruction.

2.1 Principles

2.1.1 FSI

FSI obtains the Fourier coefficients of an object by projecting a series of Fourier basis images, and forms the image through the IFT. The patterns are shown in Eq. (1) where $\left ( u, v\right )$ denote the spatial frequency, $\left ( x, y\right )$ denote the 2D Cartesian coordinates, $\phi$ is the phase of the projected pattern, $a$ is the DC term governing the average intensity of the projected pattern and $b$ is the contrast.

The corresponding detected signal at the bucket detector is given by

 

Fig. 1. Sampling maps and the reconstructed images of FSI and S-FSI (image size: $256\times 256$, sampling ratio:$~11\%$) (a) Sampling map of FSI; (b) Reconstructed image of FSI; (c) Sampling map of S-FSI; (d) Reconstructed image of S-FSI. Note that both FSI and S-FSI are using the IFT for reconstruction.

Download Full Size | PDF

2.1.2 S-FSI

S-FSI is a random sampling method in the Fourier domain, with a variable density which is governed by the sampling probability $\rho$ as shown in Eq. (5)

It can be seen that the noise of the formed image (Fig. 1(b)) using FSI is relatively low compared to the one formed by S-FSI (Fig. 1(d)), but due to the lack of high-frequency components, significant oscillation is observed. The formed image using S-FSI retains detailed information, but it is noisy. It has been proved that the noise of S-FSI can be suppressed by the compressed sensing algorithm [25], but at the cost of computation time and complexity especially for high-resolution images .

Therefore, the research motivation of this article is to look for a sampling method that only uses the IFT for image reconstruction, which can not only retain key high-frequency information, but also suppress noise in reconstruction when the number of samples is limited.

2.2 Radial correlation of significant Fourier coefficients

Before introduce our proposed A-FSI, we present the radial correlation between low-frequency and high-frequency components, which is the basis of this method.

The main idea of A-FSI is to use the locations of the significant low-frequency components to predict where significant high-frequency components may appear, and then form a sampling map through an iterative process. In order to quantitatively describe the radial correlation between the locations of high-frequency and low-frequency components, we conducted the following verification, using 500 images from the ILSVRC2012 database [26], with each image uniformly scaled to 256x256.

We take a Fourier transformation of each image and uniformly divide the Fourier domain into several regions in polar coordinates, as shown in Fig. 2. Due to the conjugate symmetry of the Fourier spectrum, we only consider half of the area of the image (see Fig. 2). In this area, we construct five concentric circles with radii $\rho _{0}$, $\rho _{1}$, $\rho _{2}$, $\rho _{3}$ and $\rho _{4}$ in sequence. The annular areas formed between concentric circles are labeled $Z^{\left ( 1 \right ) }$, $Z^{\left ( 2 \right ) }$, $Z^{\left ( 3 \right ) }$, $Z^{\left ( 4 \right ) }$, and the central area is denoted $Z^{\left ( 0 \right ) }$. Regions $Z^{\left ( 0 \right ) }$ to $Z^{\left ( 4 \right ) }$ are marked as orange, red, yellow, green and blue respectively.

 

Fig. 2. Schematic diagram of regional division rules

Download Full Size | PDF

Each region is divided into 6 fan-shaped areas. The apex angle of each fan-shaped area is 30 degrees. Each fan-shaped sub-area is marked by $Z_{q}^{\left ( i \right )}$ where the superscript $i$ represents the number of the concentric ring ( $i=0\cdots 4$), and the subscript $q$ indicates the $q$-th ($q=1\cdots 6$) direction in this concentric ring. The dashed box in the Fig. 2 is given as the example to illustrate this naming method.

The method of calculating the radial correlation between significant Fourier coefficients of the high-frequency and low-frequency is as follows.

Table 1. The example of calculating ${P_{Z^{(0)} Z^{(\imath )}}}$ and ${P_{Z^{(i)} Z^{(i+1)}}}$

View Table | View all tables in this article

In the verification, we use the ILSVRC2012 database [26] as the data set. Two ways are considered to select the radius :1. $\rho _{i+1}-\rho _{i}=c_{1}$ ; 2. $\rho _{i+1} / \rho _{i}=c_{2}$;

The test result is shown in Fig. 3. It can be seen that the radial correlation governed by $\overline {P_{Z^{(0)} Z^{(\imath )}}}$ is obvious, and even the region far away from $Z^{(0)}$ has a high correlation (for example: $\overline {P_{Z^{(0)} Z^{(4)}}} > 60\%$ ). At the same time, this correlation gets higher and higher as the distance gets closer.

 

Fig. 3. Test the values of $\overline {P_{Z^{(0)} Z^{(\imath )}}}$ and $\overline {P_{Z^{(i)} Z^{(i+1)}}}$ in the database

Download Full Size | PDF

On the other hand, $\overline {P_{Z^{(i)} Z^{(i+1)}}}$ shows a stronger correlation: the values of $\overline {P_{Z^{(i)} Z^{(i+1)}}}$ are all higher than $70\%$ . More importantly, $\overline {P_{Z^{(i)} Z^{(i+1)}}}$ increases as $i$ grows, which indicates that this correlation is enhanced in the region of high frequency.

According to the behaviors of $\overline {P_{Z^{(0)} Z^{(i)}}}$ and $\overline {P_{Z^{(i)} Z^{(i+1)}}}$ shown in Fig. 3, we can see that significant high-frequency components are also likely to appear in the radial directions where the significant low-frequency components are located, and this correlation is more obvious in the adjacent regions. Based on this phenomenon, the technique of A-FSI is proposed.

2.3 Adaptive sampling method in the frequency domain

The method proposed in this paper contains 4 steps in total. We have defined some variables, as shown in Table 2. Due to the symmetry of the Fourier coefficients, we only consider half of the spectrum points during the sampling process, and the other half is filled by symmetry.

Table 2. The meanings of the variables used in algorithm

View Table | View all tables in this article

The sampling map is obtained iteratively. For the $p$-th iteration, the input is the spectra $F_{\Delta }^{(p-1)}$ obtained from the $p-1$th iteration , and the corresponding sampling map $S_{\Delta }^{(p-1)}$. The algorithm flow chart is shown in Fig. 4, where $p$ is equal to 2.

 

Fig. 4. Flow chart of the proposed method

Download Full Size | PDF

Step 1: Select the seed points. Sort the moduli of the Fourier coefficients in $F_{\Delta }^{(p-1)}$ (Fig. 4(a)) in descending order , and take the largest $m$ Fourier coefficients as seed points (Fig. 4(b)). The set of coordinates of these chosen Fourier coefficients in $F_{\Delta }^{(p-1)}$ is expressed as $\left \{\left (u_{\Delta, n}^{(p-1)}, v_{\Delta, n}^{(p-1)}\right )\right \}$, where $n$ denotes the $n$-th point after sorting ($n=1 \cdots m$), and the superscript $(p-1)$ denotes the $(p-1)$ th iteration. Here, we define the $p$-th growth coefficient $k^{(p)}$ which is given by Eq. (7) where $N_{\Delta }^{(p-1)}$ is the number of sampling point in $S_{\Delta }^{\left ( p-1 \right ) }$

Step 2: Coordinates expansion. The coordinates of the center point are $\left (u_{c}, v_{c}\right )$, the coefficient of expansion is $\alpha$ and the coordinates of the expanded seed points are $\left (u_{e, n}^{(p-1)}, v_{e, n}^{(p-1)}\right )$, as shown in Eq. (8).(See Fig. 4(c)).

Step 3: Growth of sampling points

Step 4: Form the image, and obtain the final $S_{\Delta }^{(p)}$ for the next iteration

Note that the original sampling area of A-FSI is half a circle, where the center of the circle is the center of the image, and its radius is $R_{0}$. The first iteration is based on this original full sampling area.

2.4 Evaluating the efficiency of the sampling method

In order to evaluate the efficiency of the different sampling strategies, two evaluation methods are used. First, we use the root-mean-square error, as shown in Eq. (9)

Secondly, we propose the coverage ratio of key spectrum as an evaluation function. In order to define this measure, we again make use of the area division method shown in Fig. 2 with $\rho_{i+1} / \rho_{i}=c_{2}$ ($c_{2} = 1.414$ for the image size $128 \times 128$ and $c_{2} = 1.515$ for the image size $256 \times 256.\;\rho_{0} = 16$).

Suppose that after $K$ iterations, the final sampling map is $S_{Final}$, and the number of sampling points is $N_{Final}$. The part of $S_{Final}$ in each $Z^{(i)}$ is marked by $S_{Final}^{(i)}$. Then, we sort the moduli of all the Fourier coefficients of the original picture in descending order. After sorting, we record the coordinates of the largest $N_{Final}$ components, and denote the corresponding sampling map $S_{ori}$. This represent the positions of the key components. The part of $S_{ori}$ in $Z^{(i)}$ is labeled $S_{ori}^{(i)}$ , and the number of sampling points in $S_{ori}^{(i)}$ is denoted by $N_{ori}^{(i)}$. We define $S_{cov}$ as the overlap area between $S_{Final}$ and $S_{ori}$, and the number of sampling points in $S_{cov}$ is $N_{cov}$. Similarly, $S_{cov}^{(i)}$ is the overlap area between $S_{Final}^{(i)}$ and $S_{ori}^{(i)}$ , and the number of sampling points in $S_{cov}^{(i)}$ is $N_{cov}^{(i)}$. These definitions are summarized in Table 3.

Table 3. the meanings of variable used in the evaluation

View Table | View all tables in this article

We define the total key spectrum coverage ratio $\eta _{T}$ (as shown in Eq. (10)), and the local coverage ratio $\eta _{Z}^{(i)}$ for each $Z^{\left (i \right )}$ (as shown in Eq. (11))

Obviously, the larger the $\eta _{T}$ and $\eta _{Z}^{(i)}$, the higher the sampling efficiency. Eq. (10 ) and Eq. (11) are also applicable to other sampling methods. In the following section, these ratios will be used to compare FSI, S-FSI, and A-FSI.

3. Simulation and experiment

In our proposed sampling method, the sampling map will grow several times to encompass the high frequency components, and a series of corresponding growth coefficients $k^{(p)}$ (see Eq. (7)) are required. The coefficients can be represented as a coefficient vector $\widehat {\boldsymbol {k}}_{\mathbf {M}}=\left [k^{(1)} \cdots k^{(p)} \cdots k^{(M)}\right ]$ where $M$ denotes the total number of iterations. The values of $\widehat {\boldsymbol {k}}_{\mathbf {M}}$ determine the shape of the sampling map and the quality of the reconstructed image. Therefore, in this section, we first compare three different strategies of selecting the values of $\widehat {\boldsymbol {k}}_{\mathbf {M}}$. Secondly, we then analyze the proposed A-FSI, and compare it with FSI and S-FSI.

3.1 Simulation

In the simulation, the target images are in the size of $256\times 256$. We obtain the sampling maps and form the images with 5 iterations. The corresponding coefficient vector is labeled $\widehat {\boldsymbol {k}}_{\mathbf {5}}$ and the corresponding $\alpha$ equals to 1.8 (see Eq. (8)).

3.1.1 Selecting the values of $\widehat {\boldsymbol {k}}_{\mathbf {M}}$

Three selection strategies for the values of $\widehat {\boldsymbol {k}}_{\mathbf {M}}$ are considered: 1. one where they gradually decrease; 2. one where they gradually increase ;3. one where they remain constant.

The results are shown in Fig. 5, where the three strategies are shown for two example images. In Fig. 5. the first column represents the spectrum sampling map $S_{Final}$ obtained by the proposed A-FSI, with different $\widehat {\boldsymbol {k}}_{\mathbf {5}}$. $S_{ori}$ are shown in the second column. The overlap between $S_{Final}$ and $S_{ori}$ is displayed in the third column, where orange and cyan points indicate those which only appear in $S_{ori}$ and $S_{Final}$ respectively, while red points indicate points which appear in both $S_{ori}$ and $S_{Final}$. Obviously the number of red points is equal to $N_{cov}$. The fourth column show the reconstructed images with different $\widehat {\boldsymbol {k}}_{\mathbf {5}}$. To ensure a fair comparison, we have made the sampling ratio (SR) of the three methods close.

 

Fig. 5. Examining the different strategies for selecting the values of $\widehat {\boldsymbol {k}}_{\mathbf {M}}$.

Download Full Size | PDF

It can be seen that the strategy of gradually decreasing (the first and fourth rows) has obvious advantages, and leads to an $\eta _{T}$ higher than the other two strategies, which in turn gives the lowest RMSE for the formed image. It is easy to see that a higher $\eta _{T}$ will significantly enhance the sampling efficiency. In other words, a higher $\eta _{T}$ means that more key frequency components are measured with a certain number of samples.

3.1.2 Comparison of different sampling methods

We now compare the proposed A-FSI with FSI [19] and S-FSI [25]. For a fair comparison, two constraints are considered. First, the sampling number. We choose sampling number of the A-FSI to be equal to that of the FSI, and lower than that of the S-FSI. Second, the full sampling radius. The initial full sampling radii of A-FSI and S-FSI are the same as 16 pixels ($R_{0} = 16$, equivalent to $R = 0.0884$ in S-FSI)

Figure 6 and 7 show the formed images at different SR by using three methods. The resolution of the target image is in the size of $256\times 256$, and the sampling ratio gradually decreases from the first column to third column. Rows one to three are formed using A-FSI, FSI and S-FSI respectively. The fourth row is the local coverage rate ($\eta _{Z}^{(i)}$) corresponding to the three methods (see Eq. (11)).

 

Fig. 6. Comparing the different sampling methods. Example 1

Download Full Size | PDF

 

Fig. 7. Comparing the different sampling methods. Example 2

Download Full Size | PDF

It can be seen that as the sampling number increases, the RMSE of the three methods drops significantly. Both A-FSI and FSI have the ability to suppress noise. However, when the SR is low, A-FSI can substantially eliminate the ringing effect in the image while preserving the significant high-frequency components (see the regions of the images highlighted by red dashed boxes in Fig. 6 and Fig. 7).

S-FSI can also obtain high-frequency information, which allows the method to retain more details in reconstructed images. However, it is obvious that the formed image is noisy when the IFT is used for reconstruction, and the corresponding RMSE is higher than that of the other two methods.

The above phenomenon can be explained by $\eta _{Z}^{(i)}$. It can be seen that the values of $\eta _{Z}^{(i)}$ of FSI are maintained at a high level in the low-frequency regions, better than the ones of A-FSI. However, A-FSI shows its own advantages in the high-frequency region where the corresponding values of $\eta _{Z}^{(i)}$ are significantly higher than the other two methods. This is why the method in this paper can reconstruct images with better details with using fewer samples. It is worthwhile to compare A-FSI and S-FSI. S-FSI has a lower coverage rate than A-FSI in most areas, so when using the IFT as a reconstruction algorithm, the obtained image is obviously affected by noise.

3.2 Experiments

The experimental scheme is shown in Fig. 8. A series of Fourier base patterns are generated and projected onto the scene by the projector (DLP4500, Texas Instruments), and the reflected light from the target is received by the detector (DET100A2, Thorlabs). The signal collected by the detector is transmitted to the computer and the image is formed. The ground truth image is achieved by FSI with a $100\%$ sampling ratio.

In single-pixel imaging, a commercial digital projector or a DMD is usually used as the spatial light modulator. The advantages of a digital projector are that it is easy to use, low in cost, and suitable for the verification of methods, but it takes a long time to display patterns, resulting in low imaging efficiency. Compared with the commercial projectors, the modulation speed of a DMD is much faster (up to 22000 Hz), which can greatly improve imaging efficiency, but its cost is higher, and the control method is relatively complicated.

As seen in Fig. 8, we use a digital projector as a spatial light modulator, and use a four-step phase shift to obtain a Fourier coefficient. In the given system, the total time for displaying the four phase-shifted patterns and collecting the corresponding intensities of the bucket detector is $\sim 1.3$ s. The formed image is in the size of 128$\times$128. The measurement time in the case of full sampling is about $\sim 3$ hours, the time to calculate the Fourier coefficients of whole image is $\sim 1$ ms , and the reconstruction time is about $\sim 5.4$ ms.

 

Fig. 8. Experimental scheme

Download Full Size | PDF

Fortunately, Zhang et al. [27] successfully simulated the Fourier basis patterns by using the binary patterns based on the Floyd-Steinberg error diffusion dithering method, thus making full use of the high-speed modulation capability of the DMD. With a resolution of 256*256, they achieved 10 frames per second at a sampling ratio of 2$\%$, which strongly promoted the development of Fourier single-pixel imaging.

In the following sections, we first verify the selection strategy of the growth coefficients $\widehat {\boldsymbol {k}}_{\mathbf {M}}$, and then compare the performance of the different sampling methods. The resolution of the target scene is 128$\times$128, the number of iterations $M=4$ and the corresponding $\alpha$ equals to 1.5.

3.2.1 Comparing the selection strategies of the growth coefficients through experiments

In this section, we consider the same three selection strategies for selecting the values of $\widehat {\boldsymbol {k}}_{\mathbf {4}}$ as we did in the simulation: gradually decreasing, gradually increasing, and remaining unchanged. These experimental results are shown in rows one to three of Fig. 9 respectively.

 

Fig. 9. Verifying the coefficient growth strategy through experiments.

Download Full Size | PDF

In Fig. 9, each row contains the sampling maps obtained by different growth coefficients and the corresponding reconstruction results. It can be seen that $\eta _{T}$ of the gradually decreasing growth coefficient (first row) is significantly higher than the other two strategies, which agrees with the results shown in the simulations. It should be noted that even though the SR of the gradual decrease strategy is slightly lower than that of the other two methods, the corresponding reconstructed image can still have the lowest RMSE.

3.2.2 Comparison of different sampling methods through experiments

Finally, we compare the different sampling methods experimentally, with the results shown in Fig. 10. As we did in simulation, the initial full sampling radii of A-FSI and S-FSI are the same as 16 pixels ($R_{0} = 16$, equivalent to $R = 0.1768$ in S-FSI)

 

Fig. 10. Comparing A-FSI, FSI and S-FSI experimentally.

Download Full Size | PDF

Rows one to three show the results of A-FSI, FSI and S-FSI respectively. From the first to the third column, the SR decreases gradually. The fourth row shows curves of local coverage of different methods. It can be seen that A-FSI consistently maintains a high $\eta _{T}$ . When the SR is low (about 10$\%$), the details of the image can still be maintained , and at the same time reconstruction noise and ringing effects can be suppressed, in agreement with the simulations (see the regions highlighted by red dashed boxes in Fig. 10). Similarly, comparing the first and third rows, we can see the SR of A-FSI is slightly lower than that of S-FSI , while the corresponding reconstructed image of A-FSI gives the lowest RMSE. This verifies the advantages of the method proposed in this article. In addition, the behavior of the local coverage is also consistent with the simulations.

4. Conclusion

This paper proposes an adaptive spectrum sampling method, which can obtain key frequency components with a limited sampling number. Furthermore, the proposed method utilizes the IFT to reconstruct images with low noise and rich details. Both image quality and reconstruction efficiency are taken into account.

However, there are some obvious limitations in the proposed method to be studied in future. First of all, the fundamental reasons for the radial correlation in Fourier domain should be studied theoretically in depth.

Secondly, for images whose high-frequency components are relatively independent of their low-frequency components, the performance of our proposed A-FSI will decrease. Therefore, estimating the locations of these high-frequency components in such cases is an important issue for A-FSI.

Furthermore, A-FSI requires several iterations to achieve the sampling map. In this article, we only discussed the influence of the approximate trend of the values of the growth coefficients $\widehat {\boldsymbol {k}}_{\mathbf {M}}$ on the results. However, it would be worthwhile to study how to achieve optimal values of each growth coefficient, to further improve the coverage of the key spectrum components.

The core idea of this method is to extend the idea of adaptive imaging into the Fourier domain. This method helps to promote the practical application of single-pixel imaging technology. On the other hand, it is also helpful to find the deeper internal correlation of images (the correlation between low frequency and high frequency Fourier components), which has potential application in image compression, transmission and other fields.