Defocusing Recovery Technology of Terahertz Image Based on 3-D PSF Simulations

Terahertz (THz) imaging has broad application prospects in nondestructive testing (NDT). However, compared to visible light, it has a low resolution and limited depth of field (DOF) due to its long wavelength. To address these limitations, a 3-D point spread function (PSF) based on a simulated THz time-domain spectroscopy (THz-TDS) system combined with the Lucy–Richardson method is proposed for image defocus restoration. The experimental results show that the contrast of THz images at different defocus positions is 2–5 times higher than that of the original image after restoration using the proposed method, and the resolution is improved by approximately 60%–80%. Moreover, the maximum contrast enhancement at the boundaries of tiny features can reach 1.81 times. For a hole-defect sample, the relative error of the defect area at different depths is < 5%. The proposed method solves the problem of defocus in THz images, enhances the resolution of these images, and facilitates NDT.


I. INTRODUCTION
I N RECENT years, terahertz (THz) waves (0.1-10 THz) have become very attractive due to their broad range of applications in basic science and technology [1], [2], [3], [4], [5], [6], [7], [8], [9], [10]. One of the most prominent advantages of THz waves is that they can penetrate various nonpolar materials, such as plastics [1], composites [2], textiles [3], and wood [4]. In addition, THz photons (1 THz is equivalent to 4 meV) have nonionizing properties [5], [6], [7], [8], [9] and many materials exhibit spectral fingerprints in this frequency range [10]. Thus, THz imaging has been used in various fields, including medical imaging [5], composite defect detection [11], [12], and safety screening [13]. Conversely, THz imaging also faces many challenges, such as limited resolution and depth of field (DOF). THz near-field technology using a subwavelength aperture probe can improve spatial resolution beyond the diffraction limit [14]. However, the Manuscript  object must be located very close to the sensor, which limits the practical applications of the approach. In addition, the temporal resolution is typically low. The use of nondiffracting beams can increase the DOF of a THz imaging system [15], but this leads to a decrease in resolution and the introduction of noise artifacts, such as sidelobes. Therefore, there is an urgent need to develop an algorithm-based approach to recover defocused THz images. Given that THz imaging is a relatively new field, the theory used to evaluate THz imaging systems [16] and the mathematical model of the point spread function (PSF) [17] are not well-established. Ning et al. [18] used the PSF measured at the focus position to perform deconvolution processing of terahertz images, and the resolution was enhanced by a factor of 0.32 times the physical size of the focused beam. The experimentally obtained PSF increases the uncertainty of this parameter due to the influence of the pinhole diameter and noise and cannot attenuate the terahertz beam in the object. Ahi [19] mathematically modeled the PSF of a THz time-domain spectroscopy (THz-TDS) system and used the results to simulate and improve the resolution of a transmission THz imaging system. Using this concept, Zhang et al. [20] used a wavelet-based approach to denoise the data to obtain a THz image of a packaged integrated circuit. The PSF model was then used to restore the denoised image. However, in the aforementioned methods, the integrated circuit pins are not clear, and they are only applied to THz transmission images. Ljubenovic et al. [21] used a Gaussian spot as the PSF and combined the FastHyDe denoising algorithm and the Lucy-Richardson (L-R) method to restore THz reflection images. However, it is unreasonable to directly use a Gaussian spot as the PSF because the excitation and transmission of the THz beam is a complex process. To date, extensive research has been conducted on THz image restoration at the focal plane [20], [21], but only a few reports have been published on defocused THz image restoration. Moreover, there is a need for 3-D PSF simulation techniques that accurately model a THz-TDS system.
To restore the defocused THz image, a recovery method based on a 3-D PSF of a simulated THz-TDS system is proposed. The effects of the focal length of the lens, spectrum, multilayer medium transmission, and absorption loss in the THz-TDS system were introduced into the simulation model, and the 3-D PSF of THz wave transmission in different media was simulated using the angular spectrum method. The PSFs at different positions in the simulated THz-TDS system were This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ compared to the experimentally obtained PSFs to evaluate the accuracy of the simulation model. On this basis, the PSF of the simulation results combined with the L-R method was used to restore THz images at different defocus positions, and the image restoration results were analyzed. In addition, the method was applied to the detection of a coin to evaluate its effectiveness for small feature recognition. The method was also applied to hole-defect sample detection to investigate its applicability in nondestructive testing (NDT).

A. Basic Imaging Concept
The reconstructed THz image by the THz-TDS system operated in the grating scanning mode can be considered as the convolution of the PSF with the function of the measured object. However, the image often contains additive noise. The output of the imaging system is expressed as follows [16]: where g(x, y) is the image function, N is the additive noise, f (x, y) is the function of the measured object, * denotes the convolution operator, and PSF(x, y, z) is the PSF of the position in space.
In the THz-TDS reflective imaging system, it should be considered that the THz pulse is not a single-frequency beam, and the spectrum of the pulse is composed of frequencies ranging from hundreds of GHz to several THz. Therefore, the reconstruction of the THz beam should consider the superposition of different frequency results, and the PSF should consider the influence of distance from the focal plane PSF(x, y; z) = ω n ω=ω 1 PSF(ω, x, y; z) (2) where ω is the THz frequency and z is the distance between the PSF plane and the lens.

B. 3-D PSF Simulation Model
To obtain the PSF at different positions in space, the beam propagation model of the THz-TDS system was established, as shown in Fig. 1. The THz wave is excited by the transmitting antenna, collimated by a small hyper-hemispherical silicon lens directly coupled to a photoconductive antenna (PCA) substrate, and focused by a plano-convex lens and a multilayer sample. Because the excitation THz wavelength (the effective spectral range of the THz-TDS system used in this study is 0.1-1.5 THz) is much larger than the spot size of the pump light irradiated on the PCA (10 µm), it can be simulated as an oscillating small electric dipole. The electric field generated by a small electric dipole in a homogeneous medium is given by [22] where µ 0 is the permeability in a vacuum, c is the speed of light in a vacuum, k = 2π/λ is the wave vector, λ is the wavelength of the incident light, r is the distance from the dipole to the observation plane,n is the unit vector of the observation direction, and p is the dipole moment. Since the refractive index of the coupling lens material is greater than the refractive index of air, there will be a total internal reflection angle. The angle is determined only by the refractive index of the lens. This results in only one circular aperture on the surface of the lens that can pass through radiation. We can calculate the electric field at anywhere W outside the lens by the Fresnel-Kirchhoff diffraction integral [23] θ 0 represents the exit direction relative to the surface normal and θ i represents the incident direction relative to the surface normal. E surf is the electric field outside the lens surface.
The propagation of the THz wave through optical devices can be expressed as where u i (x, y) is the complex amplitude distribution of the input light field before the optical device and φ(x, y) and a(x, y) are the phase difference change and absorption introduced by the device, respectively. Mathematically, a thin lens can be represented by where x and y are the Cartesian coordinates with the center of the lens as the origin, F is the focal length of the lens (if it is a focusing lens, then F > 0), and R is the radius of the lens. According to the angular spectrum theory, the complex amplitude distribution of the light field on the observation plane at a distance z from the aperture plane can be obtained using the following formula: where A( f x , f y ; 0) is the angular spectrum of the incident beam in the aperture plane, H ( f x , f y ) is the transfer function, FFT represents the fast Fourier transform, FFT −1 represents the inverse fast Fourier transform Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
where (x 0 , y 0 ) and (x, y) are the Cartesian coordinates of the aperture plane and the observation plane, f x = x/(λ z), f y = y/(λ z), and ( f x, f y) are the frequency-domain coordinates, and n is the refractive index of the current transmission medium.
When the beam passes through the multilayer sample, each layer of the multilayer sample is regarded as a thin phase plate. Its propagation in the current sample is then calculated using the angular spectrum method, and its phase change function is given as follows: where n sample is the refractive index of the current layer sample in the multilayer sample and d is the thickness of the current layer sample. The broadband THz wave can be decomposed into a linear combination of a series of monochromatic light wave components [24]; therefore, the complex amplitude distribution of the light field at any point in space can be expressed as follows: where is the weight of the monochromatic wave component corresponding to the frequency of the incident light and is the Fourier spectrum coefficient of the system reference signal. The entire observation plane U (x, y; z) is a broadband terahertz wave at this plane PSF(x, y; z).

C. Image Restoration
There are many methods for image restoration, such as the Wiener filtering algorithm [25], regularization filtering algorithm [26], and the L-R method [27], [28]. When the Wiener filtering algorithm or regularization filtering algorithm is applied, the a priori acquisition of noise power plays an important role in the process of image deblurring. The effect of restoration is often unsatisfactory without the accurate acquisition of noise power. When the L-R method is applied, prior knowledge of the noise power is unnecessary. Because THz images are limited by different imaging systems and practical application scenarios, the noise power of each terahertz image cannot be accurately obtained. Therefore, only the L-R method was used for image restoration in the present study. The algorithm was based on the maximization of the likelihood that the resulting image is an instance of the original image based on the use of Poisson statistics. The number of iterations of different images has a certain range. For example, the number of iterations of THz images used in this investigation was 10-50. After the original THz image was obtained, the iteration times were determined by combining the range of the THz image iteration times, the image quality judgment function [29], and the image quality subjective judgment [30]. When the number of iterations k is determined, the iterative equation can be expressed as follows: . (12) The flowchart of the method used in this investigation is shown in Fig. 2. After obtaining the relevant system parameters, the PSF corresponding to the original THz image position was obtained using the 3-D PSF simulation model. If the system was unchanged, only one simulation of the 3-D PSF model was required. The initial value of the iteration number k = k i was then determined based on the range of the iterations and inputted into the L-R method together with the THz original image and PSF to determine whether the image quality was optimum. If not, then k = k i+1 , and the restoration was performed again. Otherwise, the THz restored image was outputted.

III. EXPERIMENTAL RESULTS AND DISCUSSION
The detection system used in this study was a T-Ray 5000 produced by Advanced Photonix, Inc., (API), as shown in Fig. 3(a). It consists of a femtosecond laser, a PCA transmitter, a PCA receiver, a time delay line, and an optical system. The THz spectrum detected using the THz-TDS system had a width of 0.02-5 THz, a spectral resolution of 3.1 GHz, a signal-tonoise ratio >70 dB, a fast scanning range of 160 ps, and a temporal resolution of 0.1 ps. It can be observed from the frequency-domain reference signal of the system in Fig. 3(b) that the effective spectrum of the system is in the range of 0.1-1.5 THz, and a signal greater than 1.5 THz is affected by noise, from which the weight ratio of different frequency monochromatic waves can also be obtained. The THz emitter and receiver were connected using a collinear adapter and fixed to a 2-D guide rail. Combined with the movement of the 2-D guide rail, the detection of the sample reflection can be completed point-by-point using a THz-TDS.

A. Simulation of 3-D PSF
Based on the above theoretical model, the 3-D PSF of the system at different defocus positions under a focusing lens with a focal length equal to 3 in was simulated, as shown in Fig. 4(a). In addition, through the receiving antenna installed on the X -Y -Z displacement platform, the actual PSF at the corresponding position of the THz-TDS system and the simulated PSF were determined as the superpositions of all frequency beams, as shown in Fig. 4(c). The weights of the different frequencies are shown in Fig. 3(b), and the simulation interval and experimental scanning step of the XY plane were set to 0.1 mm. Fig. 4(b) shows the comparison of strength between the center sections of the actual and simulated PSFs at the corresponding positions outcomes can show more intuitively the consistency between the actual and simulated PSF. To avoid additional interference, all the color images in this article are normalized pseudo-color images. It can be observed from the PSF at the focus position and the defocusing position at 5 mm in Fig. 4 that the simulated and actual PSFs are almost the same in the middle position, and the actual PSF energy at the boundary position is higher, which is mainly due to the receiving antenna. The super-hemispherical silicon lens couples into the obliquely incident THz wave; however, when the defocusing amount is large, as shown in Fig. 4, at the defocused 10-mm position, the simulated and real PSF will be significantly different at the boundary, mainly because of additional noise and resolution. In general, the simulated PSF corresponds well to the actual PSF, and the simulated PSF is not affected by noise.
To analyze further the accuracy of the simulated spot, the lateral resolution can be defined according to the full-widthat-half-maximum (FWHM) of the spot [31], and the Gaussian function is used to extract the beam diameter where µ and σ are the expectation and variance of the Gaussian function, respectively, and the beam diameter is  given by [32].
Table I shows that the FWHM increases as the defocusing amount increases, and the FWHM error of the simulated spot is smaller than that of the measured spot. In addition, according to the Rayleigh distance formula, 1 THz was used as the center frequency of the THz-TDS system, and the DOF of the system was approximately equal to ±2.8 mm.

B. Defocus Recovery Based on Simulated PSF
A focusing lens with a focal length of 3 in was used to perform reflective grating scanning on the resolution plate shown in Fig. 5(a) with a scanning step of 0.1 mm at the focal plane position, the defocusing 5-mm position, and the defocusing 10-mm position. The resolutions at different positions on the resolution plate are listed in Table II. In this study, the THz time-domain (TD) minimum imaging method was used to obtain the THz images. The calculation is conducted using (15). The selected time interval was in the range of 3 ps before and after the characteristic peak of the sample's TD signal, as shown in Fig. 5(b) f (x, y) = min(E(x, y, t)), t ∈ [t min , t max ] where E(x, y, t) is the THz TD signal and t min and t max are the minimum and maximum times of the selected time range, respectively.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  Fig. 6(a) shows that the clarity of the original THz image decreases as a function of the amount of defocus. For example, when the focal plane position is located, the 1-1 line pair set in the original image can also be distinguished. In the defocus 5-mm position, this pair set cannot be distinguished. When the defocus is 10 mm, two stripes are obtained. The sharpness of the original image decreased as the defocus increased, which also led to a decrease in the sharpness of the restored image, as shown in Fig. 6(d). When the defocus amount reaches 10 mm, as shown in Fig. 6(d3), due to the poor quality of the original THz image, the restored image has a ringing effect. Many methods have been used to eliminate the ringing effect in the frequency and space domains, which is not the focus of this study; therefore, it is not analyzed for the time being. As the amount of defocus increased, the noise in the original THz image increased. This hindered the restoration effect of the Wiener method, which is sensitive to noise, and the deconvolution method as the effect of noise is not considered in this approach. Finally, using the local detail map at the same position in Fig. 6(a)-(e), it is evident that the fringes of the THz original image are merged, the restoration effect of the proposed method is the best at the same line set position, and the fringes can be clearly distinguished. The most important point is that the remaining steps of the method should be conducted based on the actual PSF, which requires considerable time and effort. Moreover, considering the discreteness of the actual measurement, the measurement results at different locations cannot be guaranteed. The proposed method only requires one simulation to obtain the required 3-D PSF simulation results without changing the system's parameters.
To analyze further the restoration results of different methods, the pixel values of the original and the restored images at different defocus positions were extracted from the red line position in Fig. 5(a), as shown in Fig. 7. As the peak-valley change at each feature position becomes more abrupt, the feature sharpness increases. Fig. 7 shows that compared with the original image at different defocus positions, the sharpness of the restored image at different line pair sets has been greatly improved, especially in the areas where the original resolution reaches the limit. However, there is no distortion or serious deformation, such as 1-1-to-1-3 line pair sets at the focus position and 0-4-to-0-6 line pair sets at the defocus position at 5 mm. When the original image is distorted extensively, it is difficult to recover it by an algorithm, such as 1-1-to-1-6 line pairs at a defocusing position of 10 mm. The proposed method is better than other approaches when the resolution limit is reached.
To quantitatively evaluate the resolution improvement of the restored image, the image contrast is defined as follows: where I max and I min are the maximum and minimum intensities, respectively, at the set of different group lines of the resolution plate in the image. Based on (16), the contrast of different line pair positions of the resolution plate in the original THz image and the restored image at different defocus positions is calculated, as shown in Fig. 8. Fig. 8 indicates that the restoration effect and the stability of the proposed method are better than those of the other methods, especially in the position with large defocus. For the presented approach, the maximum contrast at the focus position increases from 35.4% to 51.63%, and the most obvious position is the 1-2 line pair set (from 8.87% to 48.82%, which corresponds to an increase of approximately 4.5 times). The maximum contrast at the defocused position at 5 mm increased from 33.34% to 46.26%, and the most obvious position is the 0-5 line pair set, which increased from 11.53% to 35.20% (approximately equal to 2.1 times). The maximum contrast at the defocused position at 10 mm increased from 24.78% to 65.15%, and the most obvious position was the 0-2 line pair set, which exhibited an increase from 17.52% to 65.15% (approximately equal to 2.7 times). The most obvious position of the increase is also the original resolution limit described above; however, the area without serious distortion and deformation corresponds well. In addition, the contrast at the other line pairs of the restored image is 2-5 times higher than that of the original image. A contrast higher than 20% was defined as the resolution of the system [35]. As shown in the blue-dotted line in Fig. 8, it can be observed that the system resolution of the focus position is increased from 0.561 to 1 lp/mm, which corresponds to an increase of approximately 78.3%; the system resolution at the defocused 5-mm position is increased from 0.445 to 0.707 lp/mm, which corresponds to an increase of ∼58.9%.
The system resolution at the defocused position at 10 mm increased from 0.354 to 0.561 lp/mm, which corresponds to an increase of ∼58.5%. The actual system resolution was slightly larger than the theoretical resolution because of noise, and   the resolution of the restored image was compared with the theoretical resolution.

C. Applications on Tiny Features and Irregular Shapes
Typically, small features and irregular shapes increase the difficulty of image restoration. Therefore, this study detected the positive and negative sides of a coin at the same detection conditions, as shown in Fig. 9(a) and (b). Here, we only detected the coin at the focal plane position and the defocusing position at 5 mm. The THz original images and restored images obtained at different positions are shown in Fig. 9(c) and (d). Fig. 9(c1) and (c2) shows that even if it is detected at the focal plane position, it is still difficult to distinguish part of the text in the front side of the coin (as shown in the blue box in Fig. 9(a), the text size is 2 × 2 mm) and some of the petal boundaries with weak features on the posterior side of the coin in the original THz image (as shown in the blue box in Fig. 9(b), the local boundary size is 1.5 × 2 mm). After restoration by the method in this study, as shown in Fig. 9(d1) and (d2), the actual structure of the local petal boundary with weaker text and features becomes clearly visible. When detected at a defocusing position of 5 mm, due to the redecline in the quality of the original THz image, as shown in Fig. 9(c3) and (c4), some simple-structured text in front of the coin (as shown in the green box in Fig. 9(a), the text size is 1 × 3 mm) and some features in the back of the coin are relatively strong. The petal boundary (as shown in the green box in Fig. 9(b), the local boundary size was 2 × 2 mm) became blurred. After restoration based on the use of the proposed method, as shown in Fig. 9(d3) and (d4), the simple structure of the text and the characteristics of the relatively strong petal boundary become clear again.
Similarly, the pixel values of the original and the restored images at different defocus positions were extracted from the positions of the red mark lines in Fig. 9(a) and (b), as shown in Fig. 10. The black-dotted line in Fig. 10 represents the pixel value of the marked line position of the original THz image at different defocus positions, shown in Fig. 9, and the red solid line is the pixel value of the marked line position of the THz restored image at different defocus positions, as shown in Fig. 9. It can be shown in Fig. 10 that the restored image has a more obvious feature contrast, which is reflected in the small feature area and feature boundary area in the image. The calculation outcomes based on (16) indicate that In addition, the eigenvalues of the images were obtained by calculating the variance of the entire image [21], and improved results for the overall target were obtained by comparing the eigenvalues before and after restoration, as shown in Table III. This table shows that the overall eigenvalues of the restored image are greatly improved, especially when the defocusing amount increases. Due to the decrease of the overall eigenvalue of the original image, the improvement of the overall eigenvalue of the restored image was more obvious, and the maximum enhancement of the overall eigenvalue was approximately 46.40% (e.g., the front image of the coin with 5-mm defocusing).

D. Application to a Hole-Defect Sample
In actual THz NDT, the sample is often thick, and the defect location is unknown. At this time, because of the influence of the DOF of the detection system, the focus position can only be changed, and the analysis is conducted one-by-one after multiple detections, which consumes considerable time   and energy. During the detection of the internal defects of the sample, the multiple reflection loss and absorption loss caused by the internal defects of the sample or uneven material will inevitably increase the difficulty of detection and image restoration. Therefore, in this study, a multilayer hole-defect sample was designed, as shown in Fig. 11(a). The sample  TABLE III  COMPARISON OF OVERALL FEATURE VALUES OF IMAGES BEFORE AND AFTER RESTORATION   TABLE IV  AREA OF DEFECTS OF TERAHERTZ ORIGINAL IMAGE AND RESTORED IMAGE AT DIFFERENT DEPTHS was fabricated using a Raise3D Pro3 Plus printer with a displacement accuracy of 0.01 mm and a printing accuracy of 0.1 mm. The printing accuracy of the printer is less than the minimum wavelength of the THz-TDS system we used, so it can be approximately considered that the sample made by the 3-D-printed sample is a uniform structure. The biodegradable material polylactic acid was used, and its optical parameters were calculated using the terahertz optical parameter extraction algorithm [36], as shown in Fig. 11(b). The dimensions of the sample were 100 × 73 × 15 mm, and five types of square hole defects with side lengths of 3, 5, 7, 10, and 15 mm were defined as defects A-E, respectively. The defect depths from top to bottom were 3, 6, and 10 mm. The defect depths of the same line are the same and correspond to defocus values of 3, 6, and 10 mm. The sample was placed on the metal plate, and the defects were located below; the system was focused on the bottom part of the sample, and it was detected at the same detection conditions. Due to the half-wave loss, the characteristic waveform of the hole defect was at its peak. Therefore, the maximum imaging method was used to obtain the terahertz images. The calculation formula is shown in (17). The selected time interval was also 3 ps before and after the TD signal characteristic peak of the sample.
The original and restored THz images obtained at different defocus positions are shown in Fig. 12. Fig. 12(a) shows that with an increase in the defocusing amount, the defect position image in the original THz image is blurred, particularly in the position with a small defect area, such as a class A defect. The sharp features become smoother, especially at the angle positions of the two defect boundaries. After the restoration of the method in this study [see Fig. 12(d)], the defect position image and sharp features recover, particularly when the defect area is small and the defocusing amount is large, as shown in Fig. 12(b) and (c).
To quantitatively evaluate the performance of the method in this study on the actual sample, the adaptive threshold segmentation algorithm was used to calculate the rectangular area [37], and the THz original image and restored image area of different depth defects were obtained, as shown in Table IV. This table shows that as the amount of defocus increases, the deviation between the defect area of the original THz image and the theoretical value gradually increases. In addition, as the area becomes smaller, the defect category deviation increases. In particular, when the defocus was 10 mm, the defect area deviation of the original image was approximately 71.11%, which is also consistent with the previous analysis. After the restoration of the method in this study, the relative error of the defect area at different depths was less than 5% except for type A defects. After the method proposed in this study is restored, the maximum relative error can be reduced by 45.11%.

IV. CONCLUSION
In summary, this study proposed a method based on the angular spectrum method to simulate the 3-D PSF of the THz-TDS system in conjunction with the L-R method for image defocusing. The contrast of the restored images at different defocus positions was 2-5 times higher than that of the original image. The resolution of the terahertz image at the focus position was increased from 0.561 to 1 lp/mm, the resolution of the THz image at the defocused 5-mm position was increased from 0.445 to 0.707 lp/mm, and the resolution of the THz image at the defocused 10-mm position was increased from 0.354 to 0.561 lp/mm. At the same time, the maximum contrast enhancement at the feature boundary of the tiny features reached 1.81 times. In the hole-defect sample, the maximum relative error of the defect area was increased by 45.11%, and the relative error of the defect area at different depths was <5%. The THz image restored by the proposed method can meet practical application requirements, solve the problem of defocusing of THz images in practical NDT applications, and expand the DOF of THz imaging systems. In addition, the method proposed in this study can also recover the THz images at the required single frequency and specific frequency band.
In future work, we will investigate different methods for the removal of the ringing effect in terahertz images and the influence of interference fringes on the imaging process, to further improve the resolution of restored terahertz images. In addition, combined with the signal processing starting from the original data, the original THz image contains more information and suppresses the noise in the original THz image. Finally, an image quality evaluation algorithm that is suitable for THz images will be investigated, and the L-R method will be optimized so that the correct number of iterations can be quickly determined to achieve improved restoration.