The THz image resolution enhancement algorithm based on PSF

Bin XIE

doi:10.11972/j.issn.1001-9014.2023.02.007

Abstract

In recent years, terahertz (THz) imaging technology has received increased attention owing to the promising imaging tools now available for nondestructive testing. However, theoretical models of THz imaging systems are yet to be developed. In this study, we proposed a comprehensive mathematical modeling and simulation theory for the THz imaging system. Gaussian beam distribution was used in the mathematical modeling of the point spread function (PSF). The target function was convolved using the PSF, and the analog transmission of THz images was realized. The PSF of the imaging system was calculated using the optical imaging method. Next, the image was deconvolution-enhanced using the PSF. The image restoration results revealed that the image resolution improved. The restored image contains more details.

Keywords

deconvolution image resolution point spread function terahetrz（THz）imaging transmission imaging

Introduction

Terahertz（THz）imaging technology has received increased attention in recent years. The ability of THz beams to penetrate nonmetallic materials enables the formation of transmitted images similar to X-rays. The THz beam is nonionized，and the system can be safely operated. Therefore，these systems are promising imaging tools for nondestructive testing in fields such as medical imaging ^［1］，encapsulated integrated circuit certification ^［2-3］，nondestructive testing of cultural relics ^［4］，and thickness measurement ^［5-6］. Despite its wide range of applications，this type of imaging system still constitutes an emerging field in imaging science. The generation and detection of THz electromagnetic radiation were not commercially available until a few decades ago. THz imagery was developed less than 20 years ago ^［7-8］. This type of imaging system is still in the nascent stage，and the perfect imaging equation has not been established yet. Their theoretical models are yet to be developed ^［9］. Perfect imaging theory and model are indispensable for advanced research on THz optics and image reconstruction.

1 Point spread function mathematical model

The resolution of an object in an optical imaging system decreases with an increase in the wavelength. The THz wave has a large wavelength（λ = 0.3 mm at f = 1 THz）. The diameter of the point spread function（PSF）and the diffraction-limited resolution are directly related to the wavelength. THz systems with higher signal-to-noise ratios（SNRs）have been developed to achieve better resolution and a larger depth of penetration. However，a trade-off between the resolution and wavelength of the light source is unavoidable. Point diffusion refers to the scattering of light energy encountered in a point light source. If an object is assumed to be a point，the image of the object is at least the width of a reflected point. However，the image of an object is usually made up of a speck of several reflection points. The central part of the speck is the brightest，and the brightness decreases as we move away from the central region. The more concentrated the blob，the better the image resolution. THz imaging equations are formulated using the PSF.

Due to the interaction of the light waves with the optical system，the incident light exhibits diffraction or scattering phenomenon，that is，the ideal point on the object plane is no longer the ideal point on the image plane. Mathematically，a two-dimensional（2D）δ function（impulse function）is used to represent the point light source. The intensity distribution of the output image is known as the impulse response. In incoherent-illumination imaging systems（e.g.，fluorescence microscopy），the imaging process follows spatial translation invariance and satisfies the superposition principle of linear systems. The imaging process for a point can be considered as the convolution of a real point by using the PSF. For common optical systems（wherein all the lenses are ideal lenses and possess circular symmetry，i.e.，all the lenses have a circular aperture），the ideal PSF corresponding to the point source on the axis on the object plane is often identified as an airy disk on the image plane and is represented by the diffraction pattern of concentric rings around the central bright spot，as shown in Fig. 1. However，in practice，spatial translational invariance is an ideal condition that no practical imaging system can satisfy，but this condition is reasonable for high-quality research instruments.

Figure 1.Convolution imaging of a point and the PSF（airy disk）（a）schematic of airy disk formation，（b）convolution imaging diagram of a point and the PSF.

An imaging system can be regarded as a linear shift-invariant system that transforms the ideal image into the image that one observes. The transformation process involves the convolution of the real object image and point diffusion image. The object transmittance function of the object transmittance plane can be mathematically expressed as follows：

O (u_{o}, v_{o}) = \iint O (x, y) δ (u_{o} - x, v_{o} - y) d x d y

，（1）

where O（u， v）is the object transmittance function，O（x， y）represents the object plane，and $δ (u_{o} - x, v_{o} - y)$ denotes the 2D δ function. Figuratively speaking，the 2D δ function is a cuboid in space，and the volume of the cuboid is 1.

In transmission imaging，the image plane domain is calculated as the superposition of each impulse function image，that is，as weighted PSF of the image plane；the superposition of the same weight functions is used in the object plane. Each point on the object surface（O（u，v））is abstracted as a spatial impulse function. The image of this point on the image plane is the convolution of the spatial impulse function with the PSF of the imaging system. Because the imaging system has different PSF for different frequency bands，the PSF of all frequency bands must be integrated. The image can be mathematically expressed as follows：

S (u, v) = \int_{f_{1}}^{f_{2}} O (u, v) \otimes P S F (u, v, f) d f

，（2）

where S（u， v） denotes the image of the point in the imaging system，and PSF（u， v，f ）denotes the image of the 2D δ function. f₁ and f₂ are，respectively the start and end frequencies of the light source.

In the imaging system proposed in this paper，the frequency of the light source is set as 2.52 THz（λ = 0.118 mm），and the output of the system is a Gaussian beam. After passing through two lenses，the beam becomes a collimating beam with a diameter of 80 mm. This collimating beam then passes through the object and shines into the THz camera. The Microxcam384 THz uncooled bolometric camera from INO is used in the proposed system. The primary parameters of this camera are as follows：a 384 × 288-pixel sensor with a pixel pitch of 35 µm，wavelength range：70-3 198 µm（4.25-0.094 THz），F-number：0.7，focal length：44 mm，object distance is 600 mm to infinity ^［10-11］. The proposed THz imaging setup is illustrated in Fig. 2.

Figure 2.Schematic of the THz imaging setup

In Fig. 2，lens 1 and lens 2 uniformly expand the beam to produce a smooth illumination area without interference fringes.

In the imaging system illustrated in Fig. 2，the magnification of the THz camera lens is set as 1：1.56. At this magnification，each pixel of the object’s surface is approximately 54.6 µm（35 µm × 1.56）. Under the 2.52-THz irradiation，the aperture of 54.6 µm is close to the diffraction limit. According to the Huygens–Fresnel principle，an aperture of this size can be treated as a point light source. This point light source produces a spherical wave that is captured by the lens of the THz camera and is refocused on the image plane to form a blurred point. This causes image blurring.

We simulated the point by using a Gaussian beam. The spot diameter of a Gaussian beam is the diameter at which the beam intensity drops to 1/e² of the peak value. The waist radius at a distance between the image point and the point source Z can be expressed as follows ^［12］：

h (z, f) = h (0, f) \sqrt[]{1 + {(\frac{λ z}{π h^{2} (0, f)})}^{2}}

，（3）

where h（0，f）is the spot radius at the center of the beam，and f is the frequency of the beam. The intensity distribution of the THz beam can be expressed using Gaussian distribution as follows：

I (u, v, z, f) = I (0,0, z, f) e x p (- 2 (u^{2} + v^{2}) / w {(z, f)}^{2})

.（4）

The full width at half maxima（FWHM）for the aforementioned Gaussian distribution（Eq. 4）can be obtained as follows

F W H M (z, f) = \sqrt[]{2 l n 2} h (z, f)

.（5）

The FWHM of diffraction-limited focused spot can be expressed as follows

F W H M (0, f) = 0.565 \frac{k}{N A} λ

，（6）

where k depends on the truncation ratio and irradiance level，NA is the numerical aperture，and λ is the wavelength. Substituting Eqs. 3，5，and 6 into Eq. 4，the relationship between PSF（z， f）and the physical parameters of the system can be expressed as follows：

P S F (z, f) = I_{0} e x p (- 2 (u^{2} + v^{2}) / {(\frac{0.565}{\sqrt[]{2 l n 2}} \frac{k}{N A} \frac{c}{f} \sqrt[]{1 + {(\frac{2 l n 2}{c π} {(\frac{N A}{0.565 k})}^{2} f z)}^{2}})}^{2})

,（7）

where $I_{0} = I (0,0, z, f)$ denotes the intensity at the center of the beam.

In the proposed imaging system，the light source is a single-frequency laser（f = 2.52 THz）；thus，the integral of frequency parameter f can be removed in Eq. 2. Accordingly，the transfer function of the proposed THz imaging system can be expressed as follows：

S (ρ) = O (ρ) \otimes I_{0} e x p (- 2 ρ^{2} / {(\frac{0.565}{\sqrt[]{2 l n 2}} \frac{k}{N A} \frac{c}{f} \sqrt[]{1 + {(\frac{2 l n 2}{c π} {(\frac{N A}{0.565 k})}^{2} f z)}^{2}})}^{2})

，（8）

where $ρ = \sqrt[]{u^{2} + v^{2}}$ is the radial distance between a point in the z-plane and the center of the beam. The output of the convolution operation is the simulated THz image.

2 THz image restoration by using PSF

The parameters of the THz imaging system are as follows：NA = 0.022，f = 2.52 THz，pupil diameter = 70 mm，focal length = 44 mm，k = 0.76，and object distance = 736 mm. First，by using Eq. 8，we compared the simulated image and the measured image of the spot with an aperture of 0.3 mm. For test point imaging，we employed the integral average method to reduce random noise and set the number of integrals as 50. The calculation results and test results of point imaging are shown in Fig. 3.

Figure 3.Simulation and test of 0.3-mm point imaging（All data are normalized and displayed in false color for ease of understanding）（a）the point imaging 3D map obtained using Eq. 8，（b）the point imaging 2D map obtained using Eq. 8，（c）the point imaging 3D map，（d）the point imaging 2D map

	MES	SNR	PSNR
Original THz image	160.042 9	60.392 2	16.382 5
Reconstructed image	39.069 5	86.567 7	27.891 1

Table 1. MES, SNR, and PSNR of the original THz image and reconstructed image

View all Tables

As can be seen from Fig. 3，the calculation results were similar to the test results and within the allowable error range. This highlights that the PSF equation established in this paper can be used. Many factors affect the imaging quality in transmission imaging，which cannot be discussed in this paper because of space constraints. However，the comparison of the computed and tested images of points revealed that PSF has a considerable effect. This confirms the previously mentioned point that the transmission imaging result of a target can be obtained as the convolution of the transmission function of the target by using the PSF. Consequently，the deconvolution of the PSF and the THz image yields a higher-resolution representation of the object. The algorithm can be mathematically expressed as follows ^［13-15］：

O (u, v) = S (u, v) \otimes^{- 1} P S F

.（9）

By using Eqs. 8 and 9，we restored the image of a simple transmission resolution plate. The result is shown in Fig. 4.

Figure 4.Image restoration by using the PSF（a）schematic of a simple resolution plate size，where R denotes the radius of the largest circle，w₁ is the width of the light transmittance gap，w₂ is the width of the narrowest part of the mask，and dimension error is ±0.01 mm，（b）the original image outputted by the camera，（c）the image restored using the PSF

By comparing Figs. 4（b）and 4（c），it can be seen that from the perspective of visual effect，the proposed image restoration algorithm based on PSF increased the image resolution to a certain extent. The restored image contained more details. The edges became sharper，and finer shades were found. The thin steel sheet connecting the central ring was not clearly visible in the original THz image. The mean square error（MES），SNR，and peak signal to noise ratio（PSNR）of the original THz image and reconstructed image were calculated and are presented in Table 1. The reconstructed images possessed improved MES，SNR，and PSNR.

3 Conclusions

In this study，we presented a method to improve the THz imaging resolution. This method is based on the fact that the transmission imaging result of the target can be obtained as the convolution of the transmission function of the target by using the PSF. First，the PSF of the imaging system is calculated using the optical imaging method. Next，the image is deconvolution-enhanced using the PSF. Furthermore，the derivation of a clear THz imaging equation and the image reduction results demonstrated that the image resolution improved. In addition，the image possessed location details and structural similarity. The results confirmed that the proposed method could be used to improve the THz image resolution. The reconstructed images possessed improved MES，SNR，and PSNR.

References

[1] S J Oh, J Choi, I Maeng et al. Molecular imaging with terahertz waves. Optics Express, 19, 4009-16(2011).

[2] K Ahi, S Shahbazmohamadi, N Asadizanjani. Quality control and authentication of packaged integrated circuits using enhanced-spatial-resolution terahertz time-domain spectroscopy and imaging. Optics & Lasers in Engineering, 104, 74-84(2018).

[3] Z Ali, F Bonnefoy, R Siragusa et al. Potential of chipless authentication based on randomness inherent in fabrication process for RF and THz(2017).

[4] H Zhang, Y Zhao, C Li et al. Characterization of the lacquerware from the Palace Museum by THz reflectometric imaging(2020).

[5] J Dong, A Locquet, D Citrin. Terahertz deconvolution based on autoregressive spectral extrapolation(2017).

[6] K Su, Y C Shen, J A Zeitler. Terahertz sensor for non-contact thickness and quality measurement of automobile paints of varying complexity. IEEE Transactions on Terahertz Science & Technology, 4, 432-9(2014).

[7] B B Hu, M C Nuss. Imaging with terahertz waves. Optics Letters, 20, 1716-8(1995).

[8] N Sunaguchi, T Yuasa. Imaging with terahertz waves. Optics Express, 17, 9558-70(2009).

[9] K Ahi. Review of GaN-based devices for terahertz operation. Optical Engineering, 56, 090901(2017).

[10] T Pope, M Doucet, F Dupont et al. Uncooled detector, optics, and camera development for THz imaging(2009).

[11] D Dufour, L Marchese, M Terroux et al. Review of terahertz technology development at INO. Journal of Infrared, Millimeter, and Terahertz Waves, 36, 922-46(2015).

[12] S Fathololoumi, D Ban, H Luo et al. Beam pattern investigation of terahertz quantum cascade lasers. PIERS Online, 4, 267-70(2008).

[13] A Liebendorfer. New deconvoluting algorithm offers new approach for mass spectrometry data compression. AIP Scilight(2020).

[14] A A Quadeer, J P Barton, A K Chakraborty et al. Deconvolving mutational patterns of poliovirus outbreaks reveals its intrinsic fitness landscape. Nature Communications, 11, 1-13(2020).

[15] M D Zeiler, G W Taylor, R Fergus. Adaptive deconvolutional networks for mid and high level feature learning(2011).