Terahertz image enhancing based on the physical model and multiscale retinex algorithm

To settle the THz image degradation problem, we propose an effective enhancement method based on the physical model and multiscale retinex (MSR) algorithm. The overall enhancing process involves two parts: reconstruction and enhancement. Firstly, the original THz images are reconstructed by a mathematical model, which is built and considered the THz absorption variate and Gaussian distribution of the beam. Then, the original images are processed by the proposed algorithm, which is combined the atmospheric scattering model and optimized MSR algorithm. The proposed algorithm not only recover the image scene radiance and remove haze, but also can make a compromise of the dynamic range of grayscale and edge enhancement of image. Results on a variety of THz images demonstrate our method can effectively improve the quality of THz images and retain sufficient image details.


Introduction
Terahertz imaging technology has been widely researched in the field of biomedical science, security control, and non-destructive detection [1].Non-destructive detection of industrial products has great significance for improving production efficiency [2].The optical inspection technology cannot observe the inside of the integrated circuit because the material of packaging cannot be penetrated.Ultrasonic testing is time-consuming and laborious, and the ultrasonic couplant used in the testing process can contaminate the products [3].Although X-ray technology is widely used in the industrial field, its ionization characteristics can cause damage to the internal circuit structure [4].However, the THz image degradation is still seriously affected by the fluctuation effects from the source, the intrinsic noise from the detector, and the limit of the THz system.The effective THz imaging technology has significant value in the applicant of non-destructive detection.
In recent years, the image enhancing techniques based on image enhancement and physical models have a significant development.In the enhancing approaches based on physical models, there are several methods.Schildknecht et al. used the bind deconvolution and a point-spread function develop the biomedical image of a histopathological sample [5].In order to distinguish counterfeit integrated circuit, Kiarash built a THz mathematical modeling for the point spread function [6].Olivieri et al. improved the resolution of single-pixel imaging by combining nonlinear THz generation with time-resolved field measurements [7].Kaiming He et al. used a dark channel prior removed the haze and recover a high quality image from a single input image [8].These approaches produce impressive results.However, it cannot obtain high quality images and recover sufficient image details.To enhance the THz image contrast and weaken these effects, many enhancing algorithms had been proposed in THz image processing.Trofimov et al. used high contrast correlation function to improve the image quality of the passive THz imaging devices [9].Chen et al. proposed the combined approach of the discrete wavelet transform and retinex algorithm to enhance continuous wave terahertz (CW-THz) scanning images contrast and denoising [10].Wong et al. enhanced the THz image resolution beyond the diffraction limit, which through reconstructing the signal using sinc-envelope and feeding it into the process of a 2D blind deconvolution procedure [11].Burford et al. applied the high-pass error function filters to obtain the balance of image clarity and signal distortion [12].Zhang et al. proposed a THz amplitude polynomial principle component regression algorithm for the inspection of aramid-basalt hybrid composite laminates [13].Muniyappan et al. enhanced images by using contrast limited adaptive histogram equalization method [14].Zhang et al. distinguished different biological tissues through adopting a composite multiscale entropy approach of THz signal complexity analysis [15].The single-scale retinex (SSR) algorithm based on a homomorphic filtering way was proposed by Jobson et al [16].Rahman et al. proposed the multiscale retinex (MSR) algorithm by linearly weighting multiple fixed-scale color channels [17].Although the methods have made the great improvements in the field of the THz image resolution, the enhanced images may amplifier noise details or generate distortion.It are still not sufficient for industrial detection.Based on the above discussion, the physical model and image enhancement for image processing have their disadvantages.We propose an effective enhancement method based on the physical model and MSR algorithm.The overall enhancing process involves two parts: reconstruction and enhancement.Firstly, the original THz images are reconstructed by a mathematical model, which is built and considered the THz absorption variate and Gaussian distribution of the beam.The original images are still hazy.Secondly, the original images are processed by the enhancing algorithm, which is based on the atmospheric scattering model and optimized MSR algorithm.The atmospheric scattering model can recover the image scene radiance and remove haze.The optimized MSR algorithm balances the dynamic range of the grayscale and edge enhancement of the images.Results on a variety of THz images demonstrate our method can effectively improve the quality of THz images and retain sufficient image details.The rest of this paper is organized as follows.In Section 2, the atmospheric scattering model and MSR algorithm are briefly reviewed.Section 3 shows details of the mathematical model and proposed algorithm.Section 4 presents the experimental and results, which include a comparison of quantitative analysis.Finial, Section 5 gives the conclusions of this work.

Theoretical background
A: Atmospheric scattering model In computer vision, it is widely used to describe the formation of blurred images using an optical attenuation model [18].It is expressed as follows: where I() is the observed image, () is the scene transmission map, J() is the scene radiance, and A is the global atmospheric light value.() describes the portion of light which is not scattered.The goal is to recover J(), (), and A from I().The first term on the right-hand side is the direct attenuation and related to the media transmission coefficient.The second term is the airlight.The direct attenuation expresses the scene radiance and its decay in the medium, and the airlight describes scattered light leads to the shift of the scene colors.The degradation of image quality is exacerbated by the brightness of the scene being exhausted by fog and distant objects.

B: MSR algorithm
To better explain the advantages of the proposed method based on the MSR, we briefly present the details of MSR algorithm.According to retinex theory, the brightness perceived by the human eye depends on the illumination of the environment and the reflection of light from the surface of the object.The mathematical expression is where I(x, y) is the observed image, L(x, y) is its illuminance function which depends on the source of illumination and resolve the dynamic range of the image, and R(x, y) is its reflectance function which depend on the target object that carries image detail and is independent of lighting.The R(x, y) represents the reflection component of a target object carrying image detail information.It corresponds to high frequency component, such as the texture information and structure of the image.The L(x, y) is estimated by low-filtering the I(x, y) using   (, ).The goal of retinex algorithm is to obtain the reflection image (x, y) by removing L(x, y) from observed image I(x, y).
The SSR could not simultaneously deal with dynamic range or tonal contrast, and may even cause halo phenomenon.For compensating the deficiency of the SSR algorithm, the MSR algorithm adopts the sum of different weighted scales of SSR.The MSR algorithm is a balance between image color fidelity and dynamic range, which is expressed as follows.
R  (x, y) = ∑   • (ln   (, ) − ln   (, ) * ln   (, )) =1 (3) where R  is the output of the th component image I(x, y) converted using MSR algorithm, and G(, ) is the Gaussian function of N different scales.Generally, N is the number of scales and set to 3 in actual applications.Typically, σ < 50 is the small scale, 50 ≤σ<100 is the medium scale, σ ≥ 100 is the large scale, and ∬ G n (x, y)    = 1.The ∑   = 1  =1 and each value of   is equal.

A: Reconstructed model
The reconstructed model based on the THz absorption variate is built.The Gaussian beam distribution was considered including the numerical aperture of the THZ-TDS system, spot radius at the THz beam waist, wavelength, and distance of the waist from the target into the model.The detailed description of the THz reconstructed model, incorporating the Gaussian distribution [19] is as follows.In a typical THz imaging system, the THz beam travels over the THz broadband.Therefore, the effective frequency of the beam must be included as a variable in the reconstructed function.Thus, I() of the reconstructed model is represented as: where (, , ) is the Gaussian beam shape coefficient, and  0 is the absorption coefficient of the THz beam.
Duvillaret proposed a method in which the absorption coefficient and refractive index were calculated after the THz pulse passed through the sample [20].
where () is the ratio of the amplitude of the THz pulse through the sample and reference THz pulse; is the thickness of the sample;  is the corresponding angular frequency; and () is the difference in the phase between the reference and sample signals.
Jepsen demonstrated that the THz beam, including side lobes, is in the form of a Bessel beam or Airy disk in the THz-TDS systems [21].As reported by Sagan, the profile of the beam becomes a pure Gaussian profile when the ratio of the Gaussian beam diameter to the diameter of the truncating aperture is set to 1 [22].The beam source is a circular aperture lens-coupled antenna, whose output is approximated by the Gaussian illumination distribution.After leaving the circular aperture and cylindrical lens of the imaging system, this beam distribution still shows a Gaussian distribution.The THz beam can be approximately expressed by a Gaussian beam.The mathematical model of the Gaussian beam is expressed as [19]: where  is factor determined by the truncation ratio and irradiance; NA is the numerical aperture of the THZ-TDS system;  is the frequency of the THz beam; z is the distance between the object and waist on the z axis; and  is the speed of light.
Finally, substituting Eq. ( 8) into Eq.( 5) gives the THz image reconstructed function (Eq.( 9)) as: but also can make a compromise of the dynamic range of grayscale and edge enhancement of the image.
The algorithm scheme is described in detail.
Step 1: Data input I(, ) The data input I(x, y) is obtained from the original THz image, which is reconstructed from the (9).
Step 2: Optimized transmission () dark is the minimum value of three channels in each pixel.According to (9), the color of the air in a hazy image is usually very similar to the atmospheric light A. µ is set to 0.3.
Step 4: Optimized MSR According the (4) and ∬ G n (x, y)    = 1 , calculate the G n (x, y) (n=3) for each scale: σ 1 = 40 , σ 2 = 80 , and σ 3 = 150 .The G n (x, y) is regarded as a low-pass filter.Put the I (x, y) and G n (x, y) into (13), compute the R  (x, y).The α and β are respectively added and regarded as the gain factor and adjustment coefficient in optimized MSR.In this work, the weighting coefficients   are 0.333, 0.333, and 0.334.

𝑅(x, y) = β • ln (
Where the α is • √ and β is 0.8.M is the number of pixels of observed image.

Results
The experiment is implemented on the T-Gauge 5000 manufactured by Advanced Photonics.The focal length of the TDS system is 75 mm and the diameter of the lens is 38 mm.The point-by-point raster scanning is adopted in this work.The scanning step size and speed of the mobile station were 0.25 mm and 50 mm s-1, respectively.The signal noise ratio of THz wave above 1 THz is small, and the available frequency is range from 0.3 to 0.8 THz.The original THz images are performed by the different approaches, including the contrast-limited adaptive histogram equalization (CLAHE), grey level transformation (GLT), multiscale retinex (MSR), and our method.Our specimens were STC89C52RC IC, PCB board, copper-coated iron bookmarks of flower and butterfly.
Quantitative evaluation of image enhancement methods is discussed in the following.respectively, with an error margin of 5%.The values of rulers X4 and T4 are 13.66 mm and 14.9 mm, respectively, with an error margin of 8%.Given these results, it can be concluded that the proposed algorithm is better.3, the EN of the image obtained from the GLT method is larger than it achieved from the CLAHE and MSR methods.The SF of the image enhanced by the MSR method is the highest.Therefore, it shows a lot of noise and distortion, which seriously affects the value of PSNR.Comparing the original THz image and enhanced images, it can be concluded that our method works well.4, the EN of the image obtained from our proposed method is larger than others.The SF of the image obtained from the proposed algorithm is slightly smaller than it from the MSR method.But, the PSNR of the image obtained from the proposed algorithm is the best.
From above discussion, the proposed algorithm shows very well able to balance the image contrast, white balance, and signal to noise ratio.

Discussion
This paper presents an effective enhancement method based on the physical model and MSR algorithm.
The overall enhancing process involves two parts: reconstruction and enhancement.Firstly, the original algorithm The original THz images obtained from the reconstructed model are hazy and further enhanced.The proposed algorithm based on the atmospheric scattering model and MSR algorithm is used to enhance the THz images.The proposed algorithm not only recovers the image scene radiance and removes haze,

𝛼 3 + 1
) • ∑   • (ln  ′(, ) − ln   (, ) * ln  ′(, ))  =1 Several indicators are used to quantitatively evaluate the performance of the different image enhancement methods.For example, the standard deviation (SD) is mainly used to measure the overall contrast of the image.The entropy (EN) is mainly used to measure the expected value of a random variable in an image.In other words, EN reflects the degree of chaos in grey-scale map.The spatial frequency (SF) is an indicator of the intensity of gray mutation in an image.For the image, the edge part of the image is the abrupt part.The more rapid the mutation, the sharper the edge.The peak signal-to-noise ratio (PSNR) describes the mean square error (MSE) relationship between the original image and the enhanced image using the enhancement method.The MSE represents the expected error between the original image and enhanced image.A: IC imaging In Figure. 1, we show the enhanced results of THz IC images.The original image is relatively blurry with low image contrast.The high image contrast of the enhanced image is achieved from the proposed algorithm.The proposed algorithm obtains a better quality with high image contrast, clear IC pins, and a rational white balance in the grayscale images.In Figure. 1 (b) and (d), the image using the CLAHE and MSR methods are close and slightly worse in the sight of visual sense.The Figure.1(c) shows good image contrast and therefore the overall enhanced image looks darker.As shown in Table 1, the EN of the image using the GLT method is larger than other methods and therefore the PSNR is large.Correspondingly, the PSNR represented is low.Comparing the original THz image and enhanced images, the image in Figure. 1 (b), (c), and (d) is overall slightly bad white balance and low contrast.The IC X-ray image is shown in Figure. 2 (a).The Figure. 2 (b) shows the enhanced IC image using the proposed algorithm.According to the results, the sizes of the outside pins are highly consistent with those obtained from the X-ray image.For example, the values of rulers X1 and T1 are 49.68 mm and 50.18 mm, respectively, with a 1% margin of error.The values of rulers X2 and T2 are 3.92 mm and 4.35 mm, respectively, with an error margin of 10%.The values of rulers X3 and T3 are 3.3 mm and 3.48 mm,
Figure. 2 (e) are controlled well.As shown in Table2, the EN of the image using the GLT method is larger than images using the CLAHE and MSR methods.The SF of the image enhanced by MSR is the highest and therefore shows a lot of noise and distortion.Correspondingly, its PSNR represented is lowest.Comparing the original THz image and enhanced images, the proposed algorithm shows a very well ability of balance of the image contrast, white balance, and signal to noise ratio.

C:
Figure. 3 (e) are controlled well.As shown in Table3, the EN of the image obtained from the GLT method

D:
Butterfly bookmark imagingThe enhanced THz images of the butterfly bookmark are shown in Figure.4.The size is 6 cm × 6.5 cm × 0.05 cm, including the stamen and texture of the petal.The original image is seriously blurry with low image contrast and saturation.The Figure. 4 (c) shows low image contrast and saturation.In Figure.4 (b), (c), and (d), the image contrast obtained from the CLAHE, MSR, and proposed algorithm are almost close in the sight of visual sense.As shown in Table
THz images are reconstructed by a mathematical model, which is built and considered the THz absorption variate and Gaussian distribution of the beam.The original THz images obtained from the reconstructed model are hazy and further enhanced.Then, the original images are processed by the proposed algorithm, which is combined the atmospheric scattering model and optimized MSR algorithm.The atmospheric scattering model can recover the image scene radiance and remove haze.The optimized MSR algorithm balances the dynamic range of the grayscale and edge enhancement of the images.Besides a subjective evaluation, the performance of the different image enhancement methods is quantitatively evaluated using several indicators, SD, EN, SF, MSE, and PSNR.The experimental results clearly show our enhancing algorithm can effectively improve the quality of THz images.Furthermore, this method can be applied to other fields as well, including solid-state, chemical and biological systems, and body security screening.