Influence of Beam Distribution on the Quality of Compressed Sensing-Based THz Imaging

This paper investigates the impact of incident beam inhomogeneity on the quality of THz imaging based on the compressed sensing (CS) method. Image sampling and reconstructions under point and Gaussian beams with various geometric parameters are compared with the standard one. The simulation results show that the geometric parameters of beams strongly affect the peak signal-to-noise ratio (PSNR) of the reconstructed images. Especially for the Gaussian one, expanding the beam size at the position of the mask (BZPM) dramatically increases the PSNR. To achieve high sensitivity and resolution, new measurement matrices correlating to the incident beam distribution are proposed and the simulation results are demonstrated. Experiments under VDI THz source reveal that by using the new matrices, the PSNR of CS-based imaging at 100 GHz is evidently improved from 6 dB to 13 dB, informing the new measurement matrices are highly efficient and accurate in removing the beam effect on CS-based THz imaging. Our results may provide a new way for the high-quality CS-based THz imaging for target recognition.


I. INTRODUCTION
The recent developments of Terahertz (THz) imaging technology have shown broad application prospects in various areas such as security, material analysis, and biomedical research [1]. However, it is difficult to obtain high-performance focal-plane-array imaging because of the high manufacture cost and design difficulties of the THz receiver chip [2]. An interesting approach in the use of compressed sensing (CS)-based THz imaging provides a cheaper solution to the problem. In the CS-based THz technique, only a low-cost single-pixel THz detector is required. The object is first compressively sampled with a sub-Nyquist sampling ratio by the single-pixel detector and then recovered by solving an underdetermined linear problem with optimization strategy [3]- [5]. Commonly, CS is used in DT-CMR and MRI clinical applications combining with deep learning techniques to reduce the scanning cost, ease the patient burden, and even yield better image quality [6], [7]. Since the sampling and reconstruction techniques directly determine the performance and quality of the imaging, developing advanced spatial light modulators (SLM) and The associate editor coordinating the review of this manuscript and approving it for publication was Fang Yang . high-performance reconstruction algorithms such as Total Variation Augmented Lagrangian Alternating-direction Algorithm (TVAL3), orthogonal matching pursuit (OMP), compressive sampling matching pursuit algorithm (CoSaMP), etc. have been intensively investigated [8]- [11]. It is worthwhile mentioning that, in real practice, the imaging noise from the sampling systems, including incident beam inhomogeneity, the non-strictly vertical instrumentation arms, and the lens material defect may seriously degrade the CS-based object recognition. However, under these real conditions, much work so far cannot completely consider the influence of system noise, particularly incident beam inhomogeneity on imaging quality. To reduce the system noise of the CS-based imaging, the ghost imaging techniques are highly suggested because the real mask correlated with the mask plane and beam distribution is calculated, showing obvious benefits of the improved signal-to-noise ratio (SNR) and shorter acquisition times [12], [13]. Nevertheless, most of the time, due to special features of the targets, the set-up of the ghost imaging system for sampling is very complicated, especially for the THz frequency range, which limits its applications.
In our work, by setting up the simulation models and conducting the experiments, we have studied the impact of THz beam distribution on the quality of the CS-based THz imaging. The sampling data are decoded using the algorithm TVAL3 to recover the original object. We further propose a practical optimization method to suppress the imaging noise from the beam inhomogeneity. Both simulation and experimental results reveal that the peak signal-to-noise ratio (PSNR) of the image is significantly improved by the optimization.  Fig.1 shows a schematic diagram of the compressive THz imaging system. Provided that the distance between the source and the Lens1 is fixed in all the simulations, the incident beam is collimated by Lens 1, passing through the object and masks, and then focusing on the single-pixel detector by Lens 2. The mask is an amplitude-only spatial modulator with a random pattern. To demonstrate how the incident THz beam distribution affects the quality of the CS-based imaging, we perform a rigorous optical simulation using the optic-design software Zemax. In the simulation set-up, two lenses with TPX material are adapted in the optic path and the focal distances are set as 100 mm. Both of the two lenses have 50.8 mm diameters and 46.13 mm radii of curvature, respectively. The masks with random binary patterns are designed by 3D-CAD software, exhibiting Gaussian distribution, and the object is set as a cross shape. The masks consist of 5 × 5 pixels and the size of each pixel is 6 mm×6 mm. The ratio of the measurement times to the number of pixels is defined as the sampling rate, which is set to 40%. The non-sequential mode (NSC) is used in Zemax to trace the ray and monitor the total power of the detector.

II. SIMULATION UNDER VARIOUS INCIDENT BEAMS
The collected data are then reconstructed using TVAL3 [9] to obtain the object image in accordance with the following equation: In details, where y i represents the measured value corresponding to the i th pattern and k is the number of the random pattern. For an image of the original object with n pixels, the value of each pixel is x j . M ij is the j th pixel of the i th mask. The ratio of k/n is the sampling rate for CS-based imaging. Due to k < n, (1) is actually an under-sampling equation. To obtain the optimal solution, total variation (TV) minimization is usually adopted as follows: Here, D i x is the discrete gradient of x at pixel i. In the TVAL3 algorithm, the Augmented Lagrangian method (ALM) is selected to solve total variation (TV) minimization problem via minimizing the Lagrangian function L(λ, µ) [14]: where λ, a k × 1 zero vector, is the Lagrange multiplier and µ is the penalty parameter. It is very difficult to solve (3) directly due to the non-differentiability and non-linearity of the TV term. To solve the equation, ω i which is a slack variable of D i x is introduced, hence the (3) can be transformed to: Here ν i , 2 × 1 zero vector, is a new Lagrange multiplier and β i is a penalty parameter. The iteration by the alternating direction method is used to obtain the minimization of the Lagrangian function L(ω i , x) in (4). The parameters used in this work is provided by [15], listed in Tabel 1. In which, tol is the stopping tolerance and maxit is the maximum total iterations. A peak signal-to-noise Ratio (PSNR) [5], [13], [16] is chosen as the main evaluation metric, focusing on the SNR of the reconstructed images. The PSNR is defined by (5): where m is the bit-width of the pixel value. MSE is the mean square error between the reconstruction data and the ground VOLUME 8, 2020 truth data of the objects, which is expressed by (6), Here, n is the number of pixels and x j is the value of the original object pixel, and x j is the value of the corresponding pixel in the reconstructed image. The PSNR is the similarity metric and the larger value indicates higher image similarity. Structural Similarity Index metric (SSIM) [17] and Semantic Interpretability Score (SIS) [18] are the other two metrics, which are commonly used to evaluate the image quality. For SSIM, there has a similarity between images and the real objects, which is determined by (7): where µ a and µ b are the mean value of x j and x j respectively. σ a and σ b are the variances for x j and x j respectively. σ ab is the variance between x j and x j , and m is the bit-width of the pixel value. For the default, k 1 = 0.01, and k 2 = 0.03. The range of SSIM is from 0 to 1. The larger SSIM means the more similarity of two images. For the SIS metric, it is usually evaluated as the overlap between the reconstructed image and the original image. To demonstrate how the incident THz beam distribution affects the quality of the CS-based imaging, the point and Gaussian beams with different geometric parameters at the frequency of 100 GHz are used in the simulation. We first simulate the configurations with a point beam defined by the cone angle θ in Zemax software. Fig.2 (a) shows the PSNR of the reconstructed image against the cone angle of point beam and the distance d between the object and the mask. It is found that the PSNR exhibits a cone angle dependence and exists a turning point near 20 • . Since the size of the object is 30 mm×30mm, the diagonal length is 42.43 mm. When the cone angle of the point light source is 20 • , the diameter of the beam size at the position of mask is approximately 46 mm, which just can fully cover the object in the simulation. When θ is larger than 20 • , the PSNR values exhibit a weak dependence of θ , however, the PSNR dramatically decreases with θ when it is below 20 • . This is due to the fact that the narrower light area cannot completely cover the whole object area, leading to a loss of the object information. On the other hand, the distance d between the mask and the object has a large influence on PSNR, as shown by the red line in Fig.2 (a). With the increase of d, PSNR decreases significantly. A small distance d may reduce the effect of diffraction and nonparallel light caused by spherical aberration between the object and mask, resulting in the high PSNR values. Under the minimized d and a cone angle larger than 20 • , a high-quality image could be reconstructed with the PSNR of about 50 dB. Fig.2 (b) and (c) illustrate the reconstructed images under a point source with 20 • and 10 • cone angles. It is seen that the quality of images under point source with θ = 20 • is much higher than that under the one with θ = 10 • . Fig.2 (d) and (e) demonstrate the reconstructed images under mixture beams consisting of two and three-point sources respectively. It is obvious that the more point sources included, the poorer quality of the reconstructed images. Fig.3 (a) shows the PSNR of the reconstructed images under the Gaussian beam against the beam size at the position of the mask (BZPM). The beam size is the beam radius at the 1/e 2 point in intensity. Illustrated from the insert in Fig.3 (a), the beam distribution is inhomogeneous with the beam size where BZPM is 13 mm. It is found that the PSNR starts at a low level of about 10 dB when BZPM = 13 mm. When the BZPM is larger than 40 mm, it dramatically increases with the BZPM and maintains at a high level of about 50 dB. Compared with the results of the standard point beam, the PSNR under the Gaussian beam is quite low when the BZPM is small. Fig.3 (b)-(e) shows the reconstructed images with the BZPM of 13 mm, 20 mm, 47 mm, and 51 mm respectively. The quality of the reconstructed image under the Gaussian beam with BZPM = 13 mm is the worst one while it is dramatically increased with the BZPM. The results are consistent with the data of PSNR in Fig.3 (a), confirming that the beam inhomogeneity significantly destroys the quality of CS-based imaging.
As demonstrated above, it is difficult to achieve a highquality reconstructed image under the real beam environment where the beam inhomogeneity noise is added to the imaging system. At present, there are several methods to remove beam and optical path noise in the imaging. Ghost imaging is one of the main methods and has a good effect on denoising [12], [13]. Such a technique is essentially a multi-photon coincidence sampling and their transverse correlations are applied for the reconstruction process through a multi-correlation sampling matrix. Unfortunately, in general, the ghost imaging system for sampling in the THz frequency range is too complex to be established.
In this work, the real mask at the position of mask plane shown in Fig.1 can be obtained by the product of beam distributions and the mask patterns to eliminate noise from beam inhomogeneity. The i th measurement data y i collected by the detector follows (8) where d j represents the beam distribution intensities of the corresponding pixel at the position of mask plane. We may obtain a normalized matrix whose size is the same as the one of the mask and each matrix element is d j . Through scanning the beam distribution at the position of the mask by a single-pixel detector, the measurements corresponding to the location in each pixel are collected to constitute matrix elements. Normalizing the matrix above, the result represents the beam distribution at the position of the mask plane and each element d j of the normalized matrix can be obtained.
The product of the beam distribution and the mask is used to calculate the matrix of the real mask. The vectors [d j M ij ] in (8) have been reshaped from the matrix of the real mask. They are utilized in the new measurement for imaging reconstruction. Since the mask pattern has a Gaussian distribution with zero mean, the far-field distribution of real beams is near Gaussian or superimposed by several Gaussian distributions, therefore the probability density function (PDF) for the product of the beam and mask also follows Gaussian distribution, which meets the requirements of the CS reconstruction conditions [20]- [22].   4 compares the PSNR and SSIM between the original and optimized images under the Gaussian beams. After optimization, the value of PSNR significantly increases to about 60 dB, which is comparable to the maximum PSNR under the point beam in Fig.2 and the SSIM reaches a high level of about 1.0. Furthermore, both optimized PSNR and SSIM are independent with BZPM of the Gaussian beam, informing that the noise is completely eliminated.
The effect of our proposed method on the reconstructed images under the extreme inhomogeneous beam is further investigated. Fig.5 (a) makes the comparisons of PSNR between the original and the optimized images under the mixed two-point beams with a 10 • cone angle, seen in the inset of Fig.5 (a). The maximum PSNR for point source comes from the simulation in Fig.2 (a). It is found that the original PSNR shows independence with the distance between two beams. The average PSNR following the process of (1) is about 10 dB, which is much lower than that of the single one. The red lines in the plots exhibit the results of the optimized image by (5), where the beam distribution elements d j is measured at the position of mask plane. It is found that after the optimization procedure, the PSNR values dramatically increase from 10 dB to approximately 50 dB. Fig.5 (b) and (c) depict the reconstructed images under mixture beam without/with optimization method respectively. It is clearly seen that the proposed method is very effective to remove the beam noise effect on the CS-based imaging.

III. EXPERIMENT AND RESULTS
On the basis of the above simulation results, the experiments for CS-based THz imaging are performed to certify the effects of our proposed method on imaging denoising. Fig.6 (a) demonstrates the set-up of the measurement system, which is consisted of a terahertz continuous wave source, two flat-convex mirrors, a single-pixel THz detector, and a lock-in amplifier. The beam is provided by VDI Company [23] with 1.4 mW at 100 GHz. The beam is collimated by Lens1, passes through the mask and object, and then focuses on the single-pixel detector by Lens 2. The random patterned mask is made of aluminum tape on a transparent quartz plate, where the opaque area accounts for half of the mask. Each mask contains 5 × 5 pixels and every pixel is 6 mm×6 mm. Placed at a distance of 50 mm from Lens 1, the mask is close to the object and the distance between them is 0 mm. Mask and Lens2 were placed on a mobile platform. For optimization, the mask, the object, and Lens2 were first moved out from the optical path briefly in order to scan the collimated THz beam distribution at the position of the mask by the detector. After that, the three removed components were moved back in the optical path to finish the object detection. The selected TVAL3 parameters are the same as those listed in table.1. Fig.6 (b) shows the measured 2D beam distribution at the position of the mask plane. It is noted that the beam distribution is quite inhomogeneous. Fig.6 (c) illustrates the beam distribution along the red dot line in Fig.5 (b), here the origin of the x-axis is set at the left endpoint of the red dashed in Fig.6 (b). It is seen that the beam at the mask plate is expanded and can be described by the Gaussian with BZPM = 20 mm approximately, as shown in the blue line in the figure. Based on the analysis of Fig.3, the beam will bring large system noise in the imaging. Fig.6 (d) and (e) present the reconstructed image by using the original and optimized measurement matrices respectively. The PSNR and SSIM increase from 6.05 and 0.43 to 13.27 and 0.87 respectively, confirming that the distortion of the reconstruction performance is much improved. In this case of Fig.6 (e), the quality of the reconstructed images is not very satisfactory. We think the limited experimental results are due to the low pixel number of masks used in the experiments. Table 2 compares the PSNR between our results and the previously reported works in [13] and [24]. It can be seen that the PSNR in [13] is similar to our results. However, the used optical path in [13] is much more complicated than ours. Further study is now going on in order to improve our experimental results.

IV. CONCLUSION
Despite Compressed sensing (CS) has been applied for THz wavelength, the current research of the high-quality CS-based THz imaging remains to be explored. In this paper, by modeling and simulation in the optical software Zemax, the influence of beam distribution on CS-based THz imaging has been investigated. The results show that under inhomogeneous illumination, the PSNR of the reconstructed image is significantly reduced. An improved quantitative recovery technique which combines the real beam distribution in the measurement matrices is then proposed. Simulation and experimental results confirm that the distortion of the reconstruction performance is much improved. Compared with the results and methods in [13] and [24], our experimental results are more valid and our method may provide a new way for high-quality CS-based THz imaging to target recognition in real-engineering applications. ZEQI ZHU received the bachelor's degree from the School of Electronic Information, Nantong University, Nantong, China, in 2017, and the master's degree from Nanjing University, Nanjing, China.
FENG YAN, photograph and biography not available at the time of publication.
XIAOLI JI (Member, IEEE) received the B.S. degree from the Department of Physics, Jilin University, and the Ph.D. degree from the University of Tsukuba, Japan. She is currently a Professor with Nanjing University, China. Her research interest includes microelectronics.