Spatial Domain Terahertz Image Reconstruction Based on Dual Sparsity Constraints

Terahertz time domain spectroscopy imaging systems suffer from the problems of long image acquisition time and massive data processing. Reducing the sampling rate will lead to the degradation of the imaging reconstruction quality. To solve this issue, a novel terahertz imaging model, named the dual sparsity constraints terahertz image reconstruction model (DSC-THz), is proposed in this paper. DSC-THz fuses the sparsity constraints of the terahertz image in wavelet and gradient domains into the terahertz image reconstruction model. Differing from the conventional wavelet transform, we introduce a non-linear exponentiation transform into the shift invariant wavelet coefficients, which can amplify the significant coefficients and suppress the small ones. Simultaneously, the sparsity of the terahertz image in gradient domain is used to enhance the sparsity of the image, which has the advantage of edge preserving property. The split Bregman iteration scheme is utilized to tackle the optimization problem. By using the idea of separation of variables, the optimization problem is decomposed into subproblems to solve. Compared with the conventional single sparsity constraint terahertz image reconstruction model, the experiments verified that the proposed approach can achieve higher terahertz image reconstruction quality at low sampling rates.


Introduction
Terahertz band refers to the electromagnetic spectrum region with frequency from 100 GHz to 10 THz, which is between millimeter wave and infrared light. Terahertz radiation has unique characteristics and applications in the field of imaging due to its perspective of many optically opaque materials, low-energy lossless and spectral resolution [1][2][3]. Therefore, terahertz imaging technology has great application potential in many fields, such as nondestructive testing, identification of hidden materials and food production quality monitoring [4][5][6][7][8]. In particular, terahertz radiation could recognize biomolecules and its photon energy is too low to cause atoms ionization, so it is attractive for noninvasive biomedical imaging [9][10][11]. In the field of biomedical imaging, the image quality is the most important standard. High-quality images require high resolution and high signal-to-noise ratio, which can be easily obtained by terahertz time domain spectroscopy (THz-TDS) imaging systems.
THz-TDS imaging is one of the simplest terahertz imaging modes. It can acquire high-quality object images with high spatial resolution by raster scanning. However, THz-TDS technology requires space scanning to perform imaging, which has the problems of long image acquisition time and massive data processing. The fast-imaging techniques of THz-TDS are desirable for its practical application. To reduce the long image acquisition time by raster scanning, various techniques for THz-TDS fast imaging have been proposed. Array detectors have been used to overcome the image acquisition speed limitations of sequential data acquisition [12,13]. However, these methods have higher complexity and operational cost. A fast-pulsed terahertz Fourier imaging method based on compressed sensing (CS) was proposed in [14], which allows image reconstruction with a lower sample rate than the traditional method required by utilizing the image sparsity in the frequency domain. However, the image acquisition speed of this method is still limited because the detector is performing a raster scan in the Fourier plane. Subsequently, a single pixel terahertz imaging system based on CS was proposed [15,16]. This system does not require mechanical scanning of the terahertz receiver on the image plane. The spatial profile of the terahertz beam passing the object is modulated by a random pattern, and the resulting beam is focused by the lens onto the single fixed detector. Then, by changing the random patterns, data corresponding to the different patterns are collected, and the image can be reconstructed by CS. Although the number of measurements used for image reconstruction is much smaller than the traditional raster scan system, an additional problem is that this method has slow translation of one random pattern to another. Although other spatial modulation schemes of a terahertz beam driven by optically or electrically were researched [17,18], this method still requires additional hardware devices for the spatial modulation.
In addition, fast spatial domain terahertz imaging using block-based CS was proposed in [19,20]. This method can shorten scan time and speed up the imaging processing of the conventional terahertz imaging systems without any hardware addition or modification. However, this method only uses the sparsity of the terahertz image in the frequency domain for reconstruction, and the reconstruction quality is degraded when the sampling rate is reduced. How to further improve the terahertz image reconstruction quality by exploiting more prior information of the image will be a crucial issue. In order to solve this issue, a novel terahertz imaging method from undersampled data which fuses the dual sparsity constraints of the terahertz image in wavelet and gradient domains is proposed in this paper. To enhance the sparsity of the terahertz image in the wavelet domain, a non-linear exponentiation transform is introduced into the shift invariant wavelet coefficients, which can amplify the significant coefficients and suppress the small ones. Simultaneously, the sparsity of the terahertz image in gradient domain is used to enhance the sparsity of the image for ensuring high-quality image reconstruction, which has the advantage of edge preserving property [21,22]. Therefore, terahertz image can be accurately reconstructed from the undersampled data.
The rest of the paper is organized as follows. Section 2 presents the spatial domain signal model of THz-TDS system. In Section 3, a novel terahertz imaging method that fuses the dual sparsity constraints of the terahertz image in wavelet and gradient domains is described in detail. The performance of the proposed method is investigated in Section 4. Section 5 gives a brief conclusion.

Spatial Domain Signal Model of THz-TDS System
Suppose that the sample image x has m × n pixels. Let N = m × n, and randomly select the M positions from the N pixels for terahertz detection; then, a complete terahertz time domain waveform data will be obtained at each detection position in THz-TDS system. Select the peak value of each time domain waveform as the pixel value of the corresponding detection position; then, the sparse terahertz imaging system model in the spatial domain can be expressed as where y is the M dimensional terahertz measured data. x is the sample image, and its m × n elements are arranged in an N dimensional column vector. R is an M × N measurement matrix with only one element equal to 1 and the others are equal to 0 in each row, and the positions of the elements with the value of 1 are determined by the detection positions.
In conclusion, the objective of the terahertz imaging is to reconstruct the sample image x from the sparse measured data y. As the spectral density of an ordinary terahertz image is usually distributed in a low-frequency band, representing strong sparsity, the terahertz imaging can be transformed into the problem of sparse signal reconstruction, and the CS-based method can be used for terahertz image reconstruction [23]. Among the various frequency domain transform methods, wavelet transform has good spatial and frequency characteristics, and it is often used as the sparse transform for image reconstruction. By utilizing the sparsity of the terahertz image in the wavelet transform domain, the terahertz image can be reconstructed by min where F denotes the discrete wavelet transform.

Proposed DSC-THz Model
In order to further improve the terahertz image reconstruction quality, a novel DSC-THz model based on the dual sparsity constraints of the terahertz image in wavelet and gradient domains is proposed in this paper. The orthogonal wavelet basis is a set of functions obtained by the basic wavelet function through translation and stretching transform. With the increase in scale, the displacement sampling interval increases with the power of 2, which could not match the local structure characteristics of the signal from the multi-scale. Therefore, oscillation and block effect may occur in the region where the signal changes sharply. In order to effectively eliminate this artificial oscillation phenomenon, the shift invariant wavelet transform is used in this paper to decompose the terahertz image, and the sparse representation of the terahertz image is obtained by where W(x) denotes the shift invariant wavelet transform of the image x, Ψ denotes the shift invariant wavelet transform matrix, and r is the shift invariant wavelet coefficient. The performance of CS depends on the sparsity of the image in the sparse domain. To enhance the sparsity, the exponentiation transform is introduced into the wavelet coefficients and has been proven to be more efficient in sparse representation [24]. Inspired by exponentiation transform [24], we transformed the shift invariant wavelet coefficients via the non-linear exponential function to enhance the sparsity of the terahertz images in the wavelet transform domain, which could amplify the significant coefficients and suppress the small ones. The proposed exponential shift invariant wavelet coefficients can be written as: where W e x denotes the exponential shift invariant wavelet transform of the image x, Ψ denotes the shift invariant wavelet transform matrix, and the wavelet coefficients are normalized here. a is a constant greater than 1, which is set to 10 in this paper. Figure 1 gives a simple example to compare the wavelet coefficients obtained by traditional wavelet transform and the proposed method. Figure 1a is a one-dimensional signal taken from a line of the two-dimensional terahertz image. Figure 1b,c are the wavelet coefficients obtained by traditional wavelet transform and the proposed method, respectively. From Figure 1, it is clear that the proposed method could further enhance the sparsity of the wavelet coefficients. Furthermore, in order to improve the reconstruction quality of the terahertz image, we also exploit the sparsity of the terahertz image in gradient domain to enhance the sparsity, which has the advantage of edge preserving property [21,22]. The gradient image of the terahertz image is obtained as where Tx represents gradient operation on the image x. T x and T y denote the gradient operators on the horizontal and vertical directions, respectively. We obtain In conclusion, combined with the sparsity constraints of the terahertz image in wavelet and gradient domains, the proposed DSC-THz model can be expressed as To solve the optimization problem (10), we convert it into an unconstrained optimization problem by adding the penalty function term where µ is regularization parameter. As the gradient operator is nonsmooth in x, it is very difficult to solve the optimization problem (11) involving multiple l 1 -norm terms.

The Proposed Algorithm
Split Bregman iteration is a method that originated in functional analysis for finding extrema of convex functionals, which can split a complex optimization problem into a small number of unconstrained subproblems to solve. Moreover, an advantage of split Bregman iteration is that the value of the regularization parameters could remain constant in the iterations, resulting in fast convergence for the optimization method [25][26][27].
By applying the split Bregman iteration scheme to our imaging method, we first define r e = W e x, d x = T x x and d y = T y x; then, the split Bregman formulation of the optimization problem (11) becomes where λ and γ are regularization parameters, and Using the idea of separation of variables, the optimization problem (12) can be decomposed into four unconstrained optimization subproblems, as follows: Subproblem 1: Solving the x subproblem Fixing r e , d x , d y , b r , b x and b y , the optimization function of the x is derived by splitting (12) As we have decoupled x from the l 1 portion of the optimization problem (12), the subproblem (16) that we must solve for x is now differentiable, and optimality conditions for x are easily derived. By differentiating with respect to x and setting the result equal to zero, we get the update rule Set Then Subproblem 2: Solving the r e subproblem Fixing x, d x , d y , b r , b x and b y , the optimization function of the r e is derived by splitting (12) which can be effectively solved by the shrinkage operator [28][29][30] r i+1 Subproblem 3: Solving the d x subproblem Fixing x, r e , d y , b r , b x and b y , the optimization function of the d x is derived by splitting (12) d i+1 which can be effectively solved by the shrinkage operator Subproblem 4: Solving the d y subproblem Fixing x, r e , d x , b r , b x and b y , the optimization function of the d y is derived by splitting (12) which can be effectively solved by the shrinkage operator The main steps of the proposed algorithm are summarized in the Algorithm 1. Additionally, the whole processing flow of the proposed spatial domain terahertz image reconstruction method is shown in Figure 2.

Algorithm 1. Proposed DSC-THz Imaging Algorithm
Input: measurement y, measurement matrix R, exponential shift invariant wavelet basis W e , horizontal gradient operator T x , vertical gradient operator T y . Initialization:

Convergence Analysis
The convergence of Algorithm 1 is described as Theorem 1, and the proof of Theorem 1 is shown in the Appendix A.

Theorem 1.
Assume that there exists one solution x * to the optimization problem (11). We then have lim

Experiments and Discussion
In this section, real terahertz data were processed to verify the performance of the proposed DSC-THz model in practice.
A standard THz-TDS laboratory setup, using reflection geometry developed by Zomega terahertz company in USA, was used in our experiment. The measurement range of this system is 5 cm × 5 cm. A typical THz-TDS reflection imaging system is shown in Figure 3. The pulsed terahertz beam driven by a femtosecond laser Ti-sapphire has a central wavelength and pulse width of 800 nm and 100 fs, respectively. The beam passes through a half wave plate and then is split into pump and probe beams by a beam splitter. The half wave plate is used to adjust the beam splitter to change the intensity of the two separate beams. The pump beam is irradiated on a photoconductive switch fabricated on a LT-GaAs wafer for generation of the terahertz waves, and the probe beam is focused onto an electro-optic ZnTe crystal for detection of the of terahertz waves [31][32][33]. The terahertz pulse emitted from the generator is focused on the sample by two metal parabolic mirrors. It is then reflected by the sample via two additional parabolic mirrors and guided to the ZnTe crystal, where it is overlapped with the probe beam. The probe beam is modulated by the terahertz field within the ZnTe crystal [34,35]. The modulated probe beam then passes through a quarter wave plate to make the phase difference π/2 between o-light and e-light, and then it is divided into two beams with mutually perpendicular polarization directions by polarization beam splitter to incident on the detector. The sample moves in a raster scanning mode, and the experiment is implemented at room temperature. The proposed sparse terahertz imaging system can be easily implemented from the conventional THz-TDS imaging system by programming the scanner to move according to the sampling positions defined by the measurement matrix. In this experiment, the sparse measurement data are obtained from the raster scan according to the sampling positions of the measurement matrix. Let N denote the number of pixels of the sample image, M denote the number of the pixels scanned, P denote the time taken to scan a pixel, and Q be the total time taken to move the sample in the raster scanning mode. The whole imaging time of the proposed sparse terahertz imaging system can be written as MP + Q, and M = N for the conventional THz-TDS imaging system. Therefore, under the certain system parameters, the imaging time is determined by the sampling rate M/N. Fast imaging can be achieved by reducing M.
In order to evaluate the performance of our algorithm, two samples are used in the experiment. Sample 1 contains two circular solids made from wheat flour under 10 MPa pressure. The thickness of the circular solid is 1.2 mm, and the diameter is 13 mm. The experimental humidity is 22%, and the experimental temperature is 24 • C. Sample 2 is a wheat seed. The experimental humidity is 15%, and the experimental temperature is 25 • C. The samples move in a raster scanning mode, and the full scan terahertz images for the two samples are shown in Figure 4. In addition to characterizing the performance quantitatively, the peak signal to noise ratio (PSNR) is used as the evaluation for the reconstruction quality of the terahertz image. The PSNR is defined as [19] PSNR = 10 log 10 peakval 2 MSE(x,x) where peakval is the maximum value of the image. MSE(x,x) is the mean squared error between the true image x and the estimated imagex. In order to analyze the performance of the proposed DSC-THz, conventional single sparsity constraint terahertz image reconstruction model (SSC-THz) [20] is given for comparison. Figure 5 shows the reconstructed terahertz images of the circular solids made from wheat flour at different sampling rates. Figure 5a,b show the reconstruction results using the proposed DSC-THz at sampling rates of 10% and 30%, respectively. Additionally, Figure 5c,d present the reconstruction results obtained by the conventional SSC-THz. Compared with the reconstructed images obtained by SSC-THz, the proposed method has obvious advantages. As seen in Figure 5, the reconstructed image obtained by the proposed method has the perceptually equivalent quality to that achieved using full scan at the sampling rate of 30%. When the sampling rate decreases, the image quality decays. In particular, the reconstructed image obtained by SSC-THz is severely degraded at the sampling rate of 10%. In order to illustrate the superiority property of the proposed method, Figure 6 shows the PSNR curves of DSC-THz and SSC-THz changing with the different sampling rates for the terahertz image of circular solids. It can be seen from Figure 6 that the proposed DSC-THz provides larger PSNR value compared to SSC-THz, which means the proposed method has better reconstruction performance. Applying the proposed DSC-THz to the terahertz imaging of wheat seed is shown in Figure 7. Figure 7a-c show the reconstruction results using DSC-THz with different sampling rates from 20% to 40%, respectively. Additionally, Figure 7d-f present the reconstruction results obtained by SSC-THz. It can be seen that the reconstructed images obtained by DSC-THz appear similar to that by full scan at the sampling rates of 30% and 40%. However, some degradations are observed in the reconstructed images obtained by SSC-THz at the same sampling rates. Moreover, when the sampling rate decreases to 20%, DSC-THz still has better reconstruction capacity, as can be seen in Figure 7a,d. It also can be seen from Figure 7 that there are some partial losses in the reconstruction images details by using SSC-THz, while the proposed method can better reconstruct the image and preserve more image details. For a close-up comparison, we enlarged the selected regions with the red rectangles in Figure 8 to evaluate the image quality. Figure 8 presents the reconstruction results of the selected regions. From Figure 8, it is obvious that the reconstructed images obtained by SSC-THz lost some detail information of the images.  To further illustrate the performance of DSC-THz, Figure 9 shows the PSNR curves of DSC-THz and SSC-THz changing with the different sampling rates for the terahertz image of wheat seed. As seen in Figure 9, DSC-THz provides higher PSNR compared to SSC-THz. Compared with Figures 6 and 9, it can be seen that the PSNR of the wheat seed image is a little lower than that of circular solids in low sampling rates, that is because the image of wheat has a more complex structure than the circular solids. Overall, DSC-THz is able to reconstruct images with better PSNR than SSC-THz. The above results demonstrated that simultaneously using the sparsity constraints of the terahertz image in wavelet and gradient domains can achieve better reconstruction results.

Conclusions
THz-TDS imaging technology has been used widely in many fields such as nondestructive testing, medical imaging and food production quality monitoring. However, THz-TDS imaging systems suffer from long image acquisition time and massive data processing because of their raster-scanning mechanism. A novel DSC-THz model based on the dual sparsity constraints has been proposed to effectively reconstruct terahertz image from undersampled data in this paper. The sparsity constraints of the terahertz image in wavelet and gradient domains are applied to the proposed model simultaneously, whose advantage is that it could enhance the sparsity and has better edge preserving property. Furthermore, we employ the split Bregman iteration scheme to tackle the optimization problem effectively. By using the idea of separation of variables, the optimization problem of DSC-THz can be decomposed into a series of subproblems to solve. Various experiment results with the real data confirm that the proposed method has the superior performance for terahertz image reconstruction. How to further achieve fast and accurate terahertz image reconstruction under lower sampling rate is the next focus of our research work.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Proof of Theorem 1. The first order optimality condition of Algorithm 1 yields that where J i+1 = ∂||d i+1 y || 1 , and ∂ denotes the subdifferential.
By Theorem 1, x * is the solution of the optimization problem (11), so x * satisfies From (A1) and (A3), we can see that (x * , r * e , d * x , d * y ) is a fixed point of the optimization problem (12) [36].
Subtracting the first equation of (A3) from the first equation of (A1), and then taking inner product operation with respect to x i+1 − x * in both sides, we obtain Next, subtracting the second equation of (A3) from the second equation of (A1), and then taking inner product operation with respect to r i+1 e − r * e in both sides, we obtain r i+1 ee , J i+1 Then, performing the above similar operations for the third equation to the seventh equation in (A3) and (A1), we have Taking the l 2 -norm on both sides of (A8), we obtain Similarly, we can obtain By summing up (A15) from 0 to K, we have From (A16), we can obtain ||T y x i+1 ee − d i ye || 2 2 + ||d K+1 ye || 2 2 (A17) As ||r e || 1 , ||d x || 1 and ||d y || 1 are convex functions, it has been proven that r i+1 ee , J i+1 [30]. Therefore, the above inequality implies that Then, considering the Bregman distance here [25] where P is in the subgradient of the convex function E at v, and D P E (u, v) ≥ 0. Let convex function E(x) = 1 2 ||Rx − y|| 2 2 . By (A22), we have As D P E (x i , x * ) ≥ 0 and D P E (x * , x i ) ≥ 0 [28], Equations (A18) and (A23) imply that Let convex function E(r e ) = ||r e || 1 . By (A22), we have where the second equality in (A29) is obtained by applying (A2). Therefore, the convergence theorem for (27) has been proven.