© The Author(s), published by EDP Sciences, 2024

1 Introduction

Terahertz (THz) radiation has attracted notable interest in the last twenty years due to its remarkable properties. Active terahertz imaging, for instance, has become an excellent choice for imaging in different applications, notably for security checking in government buildings and airports, thanks to its great penetrability through fabrics and clothing [1, 2]. Additionally, the non-invasive nature of THz radiation for human health renders it more suitable than traditional X-ray systems [3, 4]. This has led to considerable efforts directed to the development of novel terahertz screening systems which are mainly based on active sub-THz or millimeter wave imaging for achieving better penetration performance [5, 6]. These systems are usually equipped with an illumination module, a beam guiding and shaping module and a detection module associated with an automatic detection and identification feature for different possible threats to assist security agents in achieving high efficiency while protecting privacy [7, 8].

Active sub-THz imaging systems for security applications have reached certain maturity with a variety of systems developed and commercialized [9–11]. Synview group realized two active systems operating at room temperature at 645 GHz and 300 GHz. They exploited active illumination combined with the frequency modulated continuous wave (FMCW) to implement unique features such as phase-sensitive detection, suppression of spurious reflections and high-resolution ranging [9]. Pacific Northwest National Laboratory (PNNL), designed an active 350 mW standoff imaging system based on a wide band heterodyne, frequency-multiplier-based transceiver system, coupled with a quasi-optical focusing system [12]. Another system developed by NASA Jet Propulsion Laboratory (JPL) operated at 600 GHz and used a quasi-optical system to achieve the mechanical scanning of the scene [13]. The team of TeraSense have designed a sub-THz security body scanner for standoff screening [10] that is based on the measurement of reflected radiation from the target using a focal-plane array of detectors. The developed system shows comparable performance as the radar approach with a better image frame rate. The same team also proposed a scanning system for mail which our work attempts to improve its performance [11]. Golenkov et al. [14], proposed a similar scanner based on linear detector array of 160 silicon metal-oxide-semiconductor field-effect transistors (Si-MOSFETs) working at a frequency of 140 GHz, with which they demonstrated a resolution better than 5 mm.

Despite recent progress in sub-THz active imaging systems, they still have room for development and performance enhancement. One area of interest is related to the quality of the probing beam, which is a major factor in the detection ability of the system. In fact, the size of the beam and its irradiance distribution must match the used detector module and ensure uniform characteristics through the probed object. Otherwise, this can be the limiting factor for the system resolution instead of usual factors such as: the diffraction limit, the detector size or the low signal to noise ratio [15]. In the case of beam shaping at the THz domain, we can count several works including Kleindienst et al. [16] who proposed a single-element refractive beam shaper made of polypropylene to convert a Gaussian beam at 625 GHz generated by a diagonal horn antenna into a flat top at a distance of 3 m. In another work [17], a double-curved beam shaper based on a single element and on the principles of geometrical optics has been designed for applications in mm-wave standoff active imaging systems. Such a lens was used to shape a divergent quasi-Gaussian beam at 102 GHz radiated from a conical horn antenna into a flat-top beam with a 40 cm radius at a distance of 3 m. Ye et al. [18] on the other hand, used diffractive structures where they designed two thin diffractive phase plates based on the phase retrieve algorithm to obtain a collimated flat-top beam at 0.3 THz. The two-phase plates with a complex surface structure have been fabricated by 3D printing technologies. The resulting sub-THz beam showed a root mean square (RMS) uniformity and a fraction of the total beam power of the output top-hat beam over 80%.

Another area of enhancement for active imaging in the sub-THz band, is based on image processing technics for enhancing the quality of the produced image. In this area, one major challenge is found in multi-element detectors where each element response to the same radiation is different. This is a common problem for most spectral regions and has been thoroughly investigated in the infrared band [19]. For instance, a well-known noise resulting from this problem is stripe noise which is manifested as parallel lines imposed on the image raising from the architecture of the imaging system in a same manner as in our proposed setup. State-of-the-art methods dealing with this type of noise can be roughly divided into three different categories, namely: statistical-based methods, filtering-based methods, and model-based methods. Methods of the first category are based on some assumptions on the statistical properties of the image, which need to be identical for all pixels [20, 21]. Filtering-based approaches attempt to separate the noise from the image in an adapted filtering domaine: block-matching and 3-D (BM3D) filter [22], wavelet domain [23], and spatial 1-D guided filtering [24]. Finally, the model-based methods use an optimization scheme where they encode priors about the noise and the clear image in its regularization term using derived mathematical model of these priors. The mostly used are information about the sparsity and the smoothness of the image [25] and the directionality of the noise [26].

In this paper, we propose an experimental setup for an active sub-THz scanning system dedicated to screening baggage and parcels. The system is based on a source and a detector commercially available and a customized beam shaping module. The latter insures a uniformly distributed intensity over the detection area to maintain the detection performance across it. Furthermore, to ensure high quality images a pre-processing denoising step was adopted to correct the response of each pixel and maintain a uniform representation of the same irradiance across the detector elements. It employs a model-based optimization step and a filtering step to better deal with the noise while preserving the details in the image.

The remainder of this paper is organized as follow: in the next section, we state the problem with the concept of an active sub-THz scanner and the design of the beam shaper associated with it. In Section 3, we describe the proposed system built in our experimental work, in addition to all the implemented features enhancing its performance. Finally, Section 4 is dedicated to the presentation and discussion of experimental results followed by a conclusion in Section 5.

2 Problem statement

2.1 Imaging setup

In order to detect and image objects through common packaging materials, the use of sub-THz radiation is a straightforward solution considering the transparent property of these materials in this electromagnetic range. However, due to the reduced emissivity variation of objects in this range, active imaging is the most effective approach to obtain better detection performance. Thus, to image an object, a sub-THz source illuminates it and a detection module captures its image. In this case, two main setups can be considered: a transmission one where the object is placed between the source and the detector which will receive the transmitted radiation through the object, or a reflective scheme where the detector is placed on the path of reflected radiation from the object. The latter case can be more challenging due to the difficulty to collect radiation and project it on the detector sensitive area, and also due to the possible effect of specular reflection in the case of some objects especially metals.

Hence, we opted for the transmission mode with the addition of a mechanical scanning system to scroll the object through the detector field of view, this latter can be chosen as a linear setup of element detectors for a compact and time efficient system. A schematic view of such system is presented in Figure 1.

Fig. 1 Schematic of the experimental setup of the terahertz package scanner.

The detection module is a sub-THz linear camera placed across the scanning direction close to the object to minimize diffraction effects. The sub-THz source is an IMPATT Diode emitting at a fixed frequency of 100 GHz, with an output beam of a quasi-Gaussian intensity distribution and quasi-circular cross-section. This motivates the use of a beam shaping system to uniformly project the output beam on the active area of the detection module. In other words, convert the circular Gaussian beam into linear flat-top one to ensure an equal distribution of the intensity along the detector elements. Finally, the used belt conveyor has a fixed speed of 20 cm/s with a highly transmissive belt in the sub-THz region.

2.2 Beam shaping

The beam shape usually defines the intensity distribution and phase of a beam along its propagation path. Beam shaping refers to the operation of changing a beam intensity distribution to a desired profile [27]. Common beam shapes include Gaussian, with maximum intensity at the center and decreasing with increased distance from it, and uniform or flat top, which refers to a constant intensity distribution over the majority of the beam cross-section. Beam shapers have been extensively studied especially in the optical range [28]; however, several works have also introduced them to the THz domain. In our setup, the source output beam has a Gaussian profile but the required distribution has to be flat top along the detector array in order to maintain a uniform response for each element of the array.

In a previous work, we demonstrated the design and fabrication of a beam shaper based on Powell lens [29] that can provide flat top shape using the same source. Our approach is based first on a numerical study using COMSOL Multiphysics to validate the design parameters, followed by an experimental validation after 3D printing the lens. The proposed lens is considered aspherical and its surface profile is given by:$$Z\left(y\right)=\frac{C{y}^2}{1+\sqrt{1+(1+Q)C^2{y}^2}}$$(1)where Q is the conic constant, C is the curvature parameter. This surface introduces aberrations that causes the redistribution of the intensity profile from the center to the beam edges, thus providing a linear divergent beam with flat top intensity distribution. Using COMSOL Multiphysics, we can visualize this redistribution by representing the electric field norm emerging from the lens and its profile across the beam at a distance of 260 mm as illustrated in Figure 2.

Fig. 2 COMSOL simulation of the beam shaper (a) 2D electric field normalized distribution, (b) electric field profile across the beam [29].

The used parameters for this simulation are: R = 14, Q = −3.2 and the incident collimated beam diameter is set to 28 mm, where R represents the curvature radius 1/C.

3 Experimental setup

3.1 Illumination module

The implemented sub-terahertz scanner is a quasi-optical imaging system based on a sub-terahertz IMPATT-diode and a linear THz camera.

The IMPATT-diode is associated with a high-gain horn antenna, emitting a continuous wave (CW) at a frequency of 100 GHz with an output power of 90 mW. The experimental measurements indicated that the output beam is quasi-Gaussian with divergence around 14°. The camera is a 256 pixels linear camera with rectangular sensors (1.5 mm vertical along the detector array and 3 mm horizontal) based on plasmonic detection [11], which allows high imaging rates reaching 5000 lps (lines per second) with a minimum detectable power at this rate of 100 nW. It can be interfaced with a USB cable to an operating software that ensures image acquisition and processing.

The provided beam shaper by the constructor is a PTFE (n = 1.45) cylindrical lens ensuring the projection of the sub-THz beam on the linear active area of the camera. However, this solution provides only spatial shaping of the beam without a uniform distribution of the irradiance on the detector array. Additionally, we will compare our proposed beam shaper to this lens in the next section, to highlight the claimed enhancement.

As it can be seen in Figure 2, the Powell lens converts a collimated Gaussian beam to a divergent beam with quasi-uniform intensity. However, as stated before, the sub-THz source beam has a divergent quasi-Gaussian profile which mandate the use of collimating optic before the Powell lens. Hence, the proposed beam shaping module is constructed using a spherical lens associated with the Powell lens as shown in Figure 3.

Fig. 3 Inside view of the experimental setup (illumination module).

The spherical lens has a 50 mm diameter with a focal length of 65 mm ensuring a collimated beam with a 16 mm diameter. Both lenses were constructed by 3D printing using the ZMorph VX 3D printer with the same material which is a Polylactic acid (PLA) filament, due to its characteristics in the THz region [30]. The lenses were then experimentally tested to verify their performance and results are shown in Figure 4 along with results from the cylindrical PTFE lens. Figure 4a shows the collimated beam after the spherical lens captured by a 2D 64 × 64 camera from TeraSense Inc (Tera-4096). Its intensity distribution was then fitted to a Gaussian profile to estimate the beam radius which was found 15 mm ± 1.5 mm. The Powell lens was then placed between the spherical lens and the 2D camera and the resulting beam was captured in Figure 4c. Additionally, we also imaged the resulting beam from the cylindrical lens in Figure 4b to demonstrate the enhancements made by the proposed beam shaper.

Fig. 4 Experimental comparison between the Powell lens and a classical cylindrical lens. (a) intensity distribution of the collimated terahertz beam in 3D (left) and 1D (right) planes respectively, (b) effect of a cylindrical lens on the THz beam, (c) uniform beam shaping after passing across the Powell lens [29].

Figures 4b and 4c depict the resulting beams from the cylindrical lens and Powell lens, respectively, captured at a distance of 120 mm. The cylindrical lens 3D beam profile is a 22 mm wide by 69.7 mm long ellipse, while the intensity distribution resembles a Gaussian distribution with a gradual decrease from its maximum at the beam center to zero. In contrast, the Powell lens beam 3D profile is changed into a rectangular beam with a uniform intensity distribution that is 25 mm wide and 69 mm long.

3.2 Non-uniformity correction

Despite a uniform distribution of the beam intensity the resulting image tend to be degraded with a particular type of noise which manifests under the form of a fixed pattern superimposed on the image. As seen in Figure 5a, horizontal lines appear on the produced image from the scanner due to a different response from line to line in the image. This noise, commonly known as the non-uniformity noise or the fixed-pattern noise, rises essentially from the difference in responses between individual detectors in the array when exposed to the same irradiance.

Fig. 5 Non-uniformity noise (a) noisy image (b) horizontal gradient of the noisy image (c) vertical gradient of the noisy image.

A straightforward approach for dealing with this type of noise is a calibration procedure, where the camera is observing a uniform source and all detectors responses are adjusted to a mean response to provide uniformity. This procedure can be executed using the provided operating software of the camera, however, we observed that the noise rapidly reappears on the image and degrades its quality.

Therefore, we proposed an adaptative approach aimed at tackling this problem continuously and in real-time. To do so, we opted for a software correction of the uniformity using image processing techniques. First, as observed in Figures 5b and 5c, the noise has a directional property depicted using the gradient components of the image where the noise is mainly present in the gradient across the scan direction (the vertical gradient in our case). Based on that, a model that describes the noise can be derived as follows:$$f_{i,j}=u_{i,j}+b_i$$(2)where f _i,j represents the (i, j)th pixel in the observed image, u _i,j is the (i, j)th pixel in the clean image and b _i represents a random bias added to the ith line of the image.

To estimate the clean image from the noisy one, an optimization model is used which is considered as an inverse problem where the estimation is conducted by minimizing an energy functional. This latter is mainly composed of a “fidelity term” ensuring the solution is close to the observed image, and a “regularization term” encoding prior knowledge on the image. The regularization also ensures a unique solution and can directly affect the correction performance. Hence, a proper choice of this term will lead to better denoising results. Traditionally, total variation regularization, for instance the one proposed by Rudin et al. [31] commonly known as the ROF model represents an appropriate choice for image denoising. In the ROF model, the regularization term was chosen to better preserve edge information by dealing with functions that can be discontinuous. This model can be written as:$$\underset{u}{\min }{\Vert u-f\Vert }_2^2+\lambda {\Vert \nabla u\Vert }_1$$(3)where ∇ represents the gradient operator, ‖.‖₁ and ‖.‖₂ are the ℓ₁ and ℓ₂ norms respectively and λ is a regularization parameter.

Inspired from this model and taking into account the noise directional property, we propose to use the 1D total variation to deal with lines and columns separately as shown in the following optimization model:$$\underset{u}{\min }\frac{1}{2}{\Vert {\nabla }_x\left(u-f\right)\Vert }_2^2+\lambda {\Vert {\nabla }_y\left(u\right)\Vert }_1$$(4)where ∇ _x , ∇ _y are the gradient operators along horizontal and vertical direction respectively. The discrete formulation of these operators, considering a 2D N × M image and using finite differences, is approximated by:$${\nabla }_{x}u\left(i,j\right)= \{\begin{array}{cc}u\left(i,j+1\right)-u\left(i,j\right)& \mathrm{if} j<N\\ 0& \mathrm{if} j=N\end{array}$$ $${\nabla }_{y}u\left(i,j\right)=\{\begin{array}{cc}u\left(i+1,j\right)-u(i,j)& \mathrm{if} i<M\\ 0& \mathrm{if} i=M\end{array}$$(5)

To solve the problem in equation (4), we opted for the method in [26] for its computational efficiency and flexibility. Thus, the variational functional is approximated as follows:$$\underset{u}{\min }\frac{1}{2}{\Vert W_1^{1/2}\left({\nabla }_x\left(u-f\right)\right)\Vert }_2^2+\frac{\lambda }{2}{\Vert W_2^{1/2}{\nabla }_y\left(u\right)\Vert }_1$$(6)where W ₁ = diag(2ϕ(∇ _x (u − f))), W ₂ = diag(2Ψ(∇ _y (u))), and diag(v) is the diagonal matrix of the vector v, ϕ(v) and Ψ(v) are given as:$$\varphi \left(v\right)=\{\begin{array}{c}\begin{array}{cc}{\left|v\right|}^{-1}& \left|v\right|>{ϵ}_1,\end{array}\\ \begin{array}{cc}{ϵ}_1^{-1}& \hspace{1em}\left|v\right|\le {ϵ}_1,\end{array}\end{array}\hspace{1em}\mathrm{\Psi }\left(v\right)=\{\begin{array}{c}\begin{array}{cc}{\left|v\right|}^{-1}& \left|v\right|>{ϵ}_2\end{array}\\ \begin{array}{cc}{ϵ}_2^{-1}& \hspace{1em}\hspace{1em}\left|v\right|\le {ϵ}_2\end{array}\end{array}$$(7)

ϵ ₁ and ϵ ₂ are small positive numbers chosen to avoid division by zero-valued components. The solution is then estimated using a gradient descent scheme as follows:$$u_{n+1}=u_n-\alpha \left({\left({\nabla }_x\right)}^T{W}_1\left({\nabla }_x\left(u-f\right)\right)+\lambda {\left({\nabla }_y\right)}^T{W}_2{\nabla }_y\left(u\right)\right)$$(8)where α is the convergence step. The stopping criteria for this iterative process is based on the difference between two successive estimations of u being less than a fixed threshold.

3.3 Guided filtering

The proposed correction scheme as it is, has a major shortcoming due to its smoothing effect on the solution. In fact, the regularization term which helps reduce the noise, will also remove strong edges along the direction of the noise. To avoid this, we propose to add a filtering step where the noise component estimated from the TV minimization is filtered to remove edges wrongly attributed to the noise. Based on previous works [19, 24], this type of structured noise is proven to be linearly related to the clear image. Thus, the best filtering approach to capture this relation is the one based on the guided filter using the produced solution from the TV minimization as a guide. In fact, the guided filter smooths an image I under the guidance of another image g by assuming that the resulting image f is linearly related to g as follows:$$f_i=a_k{g}_i+b_k, \hspace{1em}\forall i\in {\omega }_k$$(9)

where ω _k is a window centred around pixel k. (a _k , b _k ) are constant coefficients in the window ω _k that are estimated by minimizing the following functional:$$\sum _{i\in {\omega }_k}\left({\left(a_k{g}_i+b_k-I_i\right)}^2+ϵa_k^2\right)$$(10)

where ϵ is the regularization term penalizing large a _k . The solution to this least square minimization is:$$\begin{array}{c}{a}_k=\frac{\frac{1}{\left|\omega \right|}\sum _{i\in {\omega }_k}{g}_i{I}_i-{\mu }_k{\overline{I}}_k}{\frac{1}{\left|\omega \right|}\sum _{i\in {\omega }_k}{I}_i^2-{\mu }_k^2+ϵ}\\ b_k={\overline{I}}_k-a_k{\mu }_k\end{array}$$(11)

where μ _k and ${\sigma }_k^2$ are the mean and variance of the guidance image g in window ω _k , ${\overline{I}}_k$ is the mean of input image I in ω _k , and |ω| is the number of pixels in ω _k .

Adding to this the directional propriety of the noise demonstrated in Figure 5, we propose to use a 1D guided filter along the noise direction (the horizontal direction in our case) by using 1D row windows ω _k in the previous equations.

Finally, the proposed correction scheme to deal with the non-uniformity between detector elements is shown in Figure 6.

Fig. 6 The overall scheme of the proposed correction method.

4 Results and discussion

To verify the performance of the proposed sub-THz scanner, a series of experiments were conducted to measure the shape of the probing beam, the spatial resolving power of the system and its detection ability.

4.1 Beam shape study

To highlight the change in the illumination beam we used the previous two beam shapers, the first one is based on the cylindrical lens and the second is the proposed one based on the Powell lens. The intensity profiles along the detector array for both cases are represented in Figure 7 with a distance of 45 cm from the cylindrical lens corresponding to its focal length and a distance of 50 cm for the Powell lens. As it can be seen, the intensity decreases from the center to the edges of the detector module for the cylindrical lens case while it maintains a constant value for the Powell lens case across the majority of the detector array. It should be noted that both lenses were positioned in slightly farther positions from the detector module to capture their characteristic intensity profiles. In an ideal position, the proposed beam shaper can produce a beam that covers all the detector elements with the same intensity. However, for the cylindrical lens, the intensity will always suffer a decrease around the edges regardless of the chosen position.

Fig. 7 Intensity profile of the sub-THz beam: (a) after the cylindrical lens, (b) after the Powell lens.

The decrease in intensity can affect the performance of the detection which we investigated in another experiment that studies the resolution of the system. To do so, a target was printed on PCB substrate containing alternating transmissive (substrate) and opaque (copper) 3 mm lines equally spaced by 3 mm with a total distance corresponding to the length of the detector sensitive area. This target is usually used to test the imaging system response to the chosen spatial frequency, while the 3 mm value is set to be equal to the vertical resolution of the system which we will demonstrate in the next experimental study. The aim is to see the effect of the beam intensity on the resolution of the system by working at the minimum value. Results are shown in Figure 8 for the cylindrical lens-based setup and Figure 9 for the proposed setup. Each figure shows the intensity variation across the array corresponding to the contrast between high and low values from transmissive and opaque regions respectively, in addition to a trace of its envelope. In the case of cylindrical lens, the contrast is more pronounced in the center of the array where the maximum intensity is located, however, when we approach the edges, it tends to decrease sharply which indicates a limited resolution at these places. For the Powell lens, the contrast is relatively constant across the array indicating a good uniformity of the illumination and ensuring a constant resolution across the system.

Fig. Contrast variation in the target image using the cylindrical lens-based system (continuous blue line: contrast variation, dashed red line: its envelope).

Fig. Contrast variation in the target image using the Powell lens-based system (continuous blue line: contrast variation, dashed red line: its envelope).

4.2 Resolution study

To measure the resolving power of the system, a target was printed on PCB substrate designed using the USAF 1951 model adapted to the sub-THz wavelength. The target has different groups of 3 cycles of opaque and transmissive stripes equally spaced with different sizes oriented horizontally and vertically as shown in Figure 10. The resulting image along with pixels values profile along the highlighted column and line are shown in Figure 10, where we compare the ideal profile to the measured one to determine the maximum spatial frequency detected by the system. The vertical resolution is estimated using results from zone 1, 2, 3 and 4, while the horizontal resolution is estimated using results from zone 5. As seen in Figure 10c, the three-cycle target profile is observed on the image for zones 1, 2 and 3 and is lost for zone 4 which suggest that spatial frequency in zone 3 is the limited vertical resolution with a 3 mm step. In the case of Figure 10d, the last three-cycle target image failed to reproduce these cycles which also limit the horizontal resolution to the spatial step of 5 mm of the previous target. Hence, the resolution of the system can be estimated to 3 mm in the vertical direction along the detector array and 5 mm on the horizontal direction perpendicular to the array. This is a clear improvement from the reported resolution in [11], which is 5 mm by 5.6 mm for vertical and horizontal directions respectively.

Fig. 10 Study of the resolution: (a) USAF target (b) produced image (c) horizontal intensity profiles for zones (a), (b), (c) and (d), (d) vertical intensity profile for line (e) (red line for ideal profile).

4.3 Non-uniformity correction performance

Finally, the system was tested using different objects to highlight its detection capacity and the advantages of correcting the fixed pattern noise using the proposed method compared to three state of the art methods, namely: GIF [24], UTV [26], WD [23]. In Figure 11b, the GIF based method have some residual noise on the corrected image from the wide horizontal stripes while in Figure 11c, the UTV method have less residual noise but the objects in the image are slightly smoothed. For the WD method in Figure 11d, the noise is efficiently reduced however the details of the image appear to be also slightly smoothed. Finally, the proposed method achieves comparable denoising performance to the WD method while maintaining the sharpness of the objects present in the image as seen in Figure 11e. These results can be related to the high ability of the TV-based methods to achieve noise reduction compared to the guided filter which explains the performance of the GIF method. On the other hand, TV-based method will produce a smoothed image unless an additional step is added to avoid this problem as it has been done in the proposed method.

Fig. 11 Correction results of the noisy image in (a) using: (b) GIF method, (c) UTV method, (d) WD method, (e) proposed method.

Furthermore, we calculated the roughness index, which measures the high-pass content of an image and is defined as follows:$$\rho =\frac{{\Vert h_1*I\Vert }_1+{\Vert h_2*I\Vert }_1}{{\Vert I\Vert }_1}$$(12)where h ₁ = [1, −1] and h ₂ = [1, −1] ^T are horizontal and vertical masks for image I. We used three sets of noisy images, each set contains several images of the same objects at different positions and orientations. Image set 2 has images containing large smooth objects and image set 3 includes objects with small details, while image set 1 has a combination of both. Roughness scores for each set are given in Table 1.

It should be noted that the lower the roughness of an image the smoother it appears and thus the less noise it contains. Therefore, we can see that besides the UTV method which has the lowest values due to its smoothing effect, our method achieves the best performance compared to other methods. Additionally, as expected, the roughness values for images with small details (image set 3) are greater than those from images of smooth objects (image set 2).

4.4 Runtime comparison

To investigate the complexity of our proposed model, we conducted simulation experiments using various images and averaged the runtime for each method in Table 2. The simulations were performed on Python 3.9 platform with an AMD Ryzen 5 5500U processor operating at 2.10 GHz and 8 GB memory.

The findings showed that the GIF method achieved the shortest running time since it is a single iteration method with simple spatial filters. While the UTV method had the longest time due to its iterative nature. The WD approach showed similar performance.

In conclusion, the runtime is a good metric for the method efficiency and computational complexity, however, denoising performance should be considered to conduct a fair assessment.

5 Conclusion

In this work, an improved sub-terahertz active imaging scanner of parcels and mail is proposed based on a commercially available solution. Such system suffered initially from the change in the distribution of the illuminating module beam intensity across the inspected area, and the non-uniform response of the detectors elements when observing the same irradiance. This can cause a reduction of detection performance especially on the detection module edges where the intensity is decreased. The main contribution of our work, was the design of a quasi-optical beam shaping module to uniformly distribute the beam intensity over the inspected area, in addition to implementing a novel correction method to reduce the non-uniform response across the detectors array. The improved beam intensity had a uniform distribution compared to the initial one and the produced image quality has improved. This has led to better performance and improved detection resolution of up to 3 mm.

Further work can be done to improve the signal to noise ratio by improving the designed Powell lens and using more transmissive materials for the manufacturing of the lenses in the beam shaper. In addition to that, an implemented feature for automatic detection of objects is needed to further study the detection ability of the proposed scanner.

Funding

This work is supported as part of the NATO-SPS Programme Multi-Year Project: Implementation of a terahertz systems dedicated to identification of products and substances, N# G5973.

Conflicts of interest

No potential conflict of interest was reported by the author(s).

Data availability statement

Data sets generated during the current study are available from the corresponding author on reasonable request.

Author contribution statement

B. Bezziou and M. Lazoul designed and fabricated the beam shaping module. A. Zebellah and A. Boutemedjet conducted the experiments. A. Boutemedjet implemented the correction process, analysed the results and wrote the manuscript with input from all authors.

References

All Tables

All Figures

	Fig. 1 Schematic of the experimental setup of the terahertz package scanner.
In the text

	Fig. 2 COMSOL simulation of the beam shaper (a) 2D electric field normalized distribution, (b) electric field profile across the beam [29].
In the text

	Fig. 3 Inside view of the experimental setup (illumination module).
In the text

	Fig. 4 Experimental comparison between the Powell lens and a classical cylindrical lens. (a) intensity distribution of the collimated terahertz beam in 3D (left) and 1D (right) planes respectively, (b) effect of a cylindrical lens on the THz beam, (c) uniform beam shaping after passing across the Powell lens [29].
In the text

	Fig. 5 Non-uniformity noise (a) noisy image (b) horizontal gradient of the noisy image (c) vertical gradient of the noisy image.
In the text

	Fig. 6 The overall scheme of the proposed correction method.
In the text

	Fig. 7 Intensity profile of the sub-THz beam: (a) after the cylindrical lens, (b) after the Powell lens.
In the text

	Fig. Contrast variation in the target image using the cylindrical lens-based system (continuous blue line: contrast variation, dashed red line: its envelope).
In the text

	Fig. Contrast variation in the target image using the Powell lens-based system (continuous blue line: contrast variation, dashed red line: its envelope).
In the text

	Fig. 10 Study of the resolution: (a) USAF target (b) produced image (c) horizontal intensity profiles for zones (a), (b), (c) and (d), (d) vertical intensity profile for line (e) (red line for ideal profile).
In the text

	Fig. 11 Correction results of the noisy image in (a) using: (b) GIF method, (c) UTV method, (d) WD method, (e) proposed method.
In the text