IMPROVEMENT OF TERAHERTZ IMAGES BY ADAPTIVE DISCRETE COSINE TRANSFORM (DCT)-BASED DENOISING

1. Introduction

Electromagnetic waves of terahertz (THz) frequencies exhibit an ability to penetrate through dielectric materials, such as plastics, paper, cardboard, fabric, etc. This attractive feature of THz illumination suggests THz imaging as a convenient and noninvasive technique well suitable for implementation in such fields as security  [1, 2], pharmacy and medical diagnostics [3, 4], food industry [5], inspection of ancient art pieces [6] and other versatile areas [7]. Direct implementation and reliability in real operational conditions require compact design solutions, convenience in use and reasonable costs of the imaging system.

In the simplest THz imaging approach – a direct THz imaging – an object is placed in the THz beam coming from the  source and the  signal is registered in transmission and/or reflection geometries. In this case, only the power of THz radiation is detected, which restricts possible applications of the system and decreases the informational value of the obtained THz images, since lots of tiny details, particularly in the case of low-absorbing objects, are not recorded. More promising approaches to THz imaging are interferometric homodyne [8] or heterodyne [9] techniques, which allow one to increase the sensitivity of the system by orders of magnitude and to extract phase information as well. However, these approaches require more elaborate experimental equipment resulting in a higher sensitivity to optical alignment, which is not always possible to achieve, as well as significantly increased economic costs. Thus, there is a need to find rational trade-off between the cost and overall complexity of the THz imaging setup, on the one hand, and the quality of resulting images, on the other one [7–9].

Despite the  rapid development of THz detectors, the  sensitive, room temperature operated detector arrays ensuring reasonable spatial resolution are still lacking. Therefore, it is common to use a single detector scheme where the THz beam is focused sharply on the object using the large numerical aperture elements and the  object is raster scanned in the plane normal to the optical axis of the  system using mechanically moving elements, while the  beam modulation by the  object is registered by a pixel detector. In such a way, a full 2D image of the object under investigation is revealed. The single detector-based imaging is a valid solution for the lack of effective THz sensing arrays; however, it does not offer real-time imaging since the image is acquired pixel-by-pixel. It is worth mentioning that the acquisition time scales approximately quadratically with a change in lateral imaging resolution, resulting thus in even greater time expenses when the  high resolution is needed. Moreover, low frequency commercial THz sources are more compact, user-friendly and have higher optical output power, but suffer from the lack of sufficient spatial resolution. As an example, the Gunn diode can operate at 0.1 THz with output power of tens of milliwatts at room temperature, and still be relatively compact. For higher frequencies Schottky diode-based frequency multipliers or molecular optically pumped lasers are usually used. These devices are rather bulky and require peripheral equipment such as a signal synthesizer for a frequency multiplier; vacuum pumps, cooling system and molecular gas lines for the laser. So, there is a strong need for a compromise between the compactness and convenience of the equipment and the  informational value of the  output images, since the lower the frequency the larger the resolved minimal size of the features in the image.

One more peculiarity of THz imaging is the fact that image quality depends critically on the  illumination conditions, the  material of the  object under study, and the parameters of the recording equipment [7]. Due to this, the images obtained in the THz range may contain noise and specific distortions, which make it impossible to extract valuable information from them even when the image resolution is high enough.

Thus, there are two ways that can help one to enhance the quality of output THz images. The first way is to improve the equipment and try to create optimal conditions for shooting, which is often associated with certain technical and economic difficulties. The other way is to increase the quality of the images by applying modern digital processing methods. Currently, there are many algorithms for image upsampling [10–12], deblurring [13–15] and denoising  [16–18] that have proven themselves well on grayscale and colour optical images, as well as on remote sensing radar images and some types of medical images. However, application of such methods to THz images is quite scarce. The processing of THz images is usually confined to automatic edge and feature detection [19], which, however, is only one of a huge variety of processing tasks that can potentially be performed. It is worth noting that with advances in terahertz imaging equipment the interest in processing of such images grows as well. For example, attempts have been made to use neural networks to increase the resolution of terahertz images [20, 21]. But in general, processing of terahertz images is at the very beginning of its development [22].

The purpose of this work is to study the possibility of improving the quality (and informational value) of THz images obtained using a fairly simple, compact and inexpensive setup operating at a frequency of 0.1 THz through the use of modern digital processing methods.

2. Samples and the experiment setup

The research object (sample) was mounted on a metallic holder with a 50 × 50 mm opening in the middle as shown in Fig. 1. Also, a metal bushing serving as the reference was placed on the holder. The holder was attached to the mechanical xyz translation stage in order to raster scan the sample. The  experiment setup was aligned in reflection geometry as depicted in Fig.  1, where the  Gunn diode equipped with a  horn antenna and radiating ~30  mW power at 0.1  THz to the  free space was employed as an emitter. The output radiation was collimated using a  high numerical aperture off-axis parabolic mirror (OAP2) and split in two beams by a  silicon beamsplitter (BS). The  beam which emanated straight was collected by another off-axis parabolic mirror (OAP1) and focused on the  sample at the  right angle. The  reflected back beam was collimated by the same OAP1 mirror and reached BS which split it once again and directed it to a plasmonic field effect transistor (FET) THz sensor through the OAP3 mirror. The sample was raster scanned and the signal from the detector was read out using the conventional lock-in technique. The  THz source was modulated electronically at a frequency of 1 kHz.

In this experiment, a piece of a car tire was used as a sample in order to evaluate the possibility of non-invasive THz based defect detection critical to industry and road safety in terms of preventing abrupt tire degradation. Moreover, it is a  dielectric sample with rich surface topology allowing for a more accurate assessment of experimental setup and performance of mathematical processing algorithms. A set of THz images of the tire was taken at distances from the OAP1 mirror (z) varying from 42 to 58 mm. It is worth noting that the focus distance of the OAP1 mirror employed in the experiment is 50 mm. The dataset obtained in this experiment contained 16 images in total. Some of them are shown in Fig. 2. The images were transformed into an 8-bit bitmap format. The size of these images is 135 × 185 pixels.

3. Analysis of noise characteristics in the obtained images

The first observation that can be made from the pictures shown in Fig. 2 is that quite many features of the sample can be distinguished, and these features differ for images obtained at different distances. The  second observation is that these images are blurry and contain visible noise, while the ratio of blur and noise differs from image to image. For some images (Fig. 2(b, c)) blur seems to be the main distorting factor, while for other images it is vice versa (Fig. 2(e, f)). Obviously, in order to be able to enhance the quality of these images it is needed to get more information on noise and distortions present.

Let us start with the  analysis of correlation characteristics of the noise. For this, the method described in Ref.  [23] was used. This method is based on estimating the  mode of local kurtosis estimates (M_K) obtained in the  discrete cosine transform (DCT) domain in non-overlapping image blocks of size 8 × 8 pixels. The advantage of this method is its independence of the noise type, thus, it can be applied to images corrupted with signal-independent, signal-dependent or mixed noise, while knowing the noise type in advance is not needed.

The value of M_K characterizes the  width of the  main lobe of the  autocorrelation function of noise and allows one to determine not only the  presence of spatial correlation, but its degree as well. If the M_K value is less than 3.75, the noise can be considered spatially uncorrelated; if M_K is larger than 3.75 and less than 5.25, the noise has a medium correlation level; and, finally, if M_K value exceeds 5.25, the noise in an analyzed image is characterized by a high level of spatial correlation.

The results of M_K estimation for different images are given in Fig. 3. Along the horizontal axis arbitrary distances are given. It is important to stress that here and further the focus point is at a value of 50 mm. Blue and red dashed lines denote medium and high correlation thresholds, respectively.

As seen, there are noticeable fluctuations in the M_K values for distances below 50 mm and the noise in these images is mostly highly correlated (for some images the level of spatial correlation is very high). Starting with z = 50 mm, or the focus point, the plot starts showing a clear tendency with M_K decreasing as distance grows; the values of M_K for this plot area are significantly lower than for the images obtained at z < 50 mm and for most images the level of spatial correlation is medium. However, for all the images the noise is correlated to some extent, which significantly complicates the task of evaluating its characteristics.

The most advanced modern methods for the blind evaluation of noise parameters either operate in the spectral domain [24–26] or are based on the usage of neural networks [27]; however, none of these methods can be directly applied for the evaluation of the characteristics of spatially correlated noise and it is needed to preliminarily eliminate the spatial correlation with additional tools. As shown in Ref. [28], one of the most effective tools allowing one to do this is image downsampling. According to Ref. [28], downsampling by a factor of 2 is usually enough to transform noise with the medium level of spatial correlation into the uncorrelated noise, while downsampling by 3 times is needed if the noise is highly correlated. Yet, an initially low resolution of the images from the dataset under study does not allow using this approach, because even after downsampling by 2 times images become too small to find enough quasi-homogeneous regions larger than the block size used. Since in this task the operability and stability of the evaluation method in the presence of spatially correlated noise is the  primary requirement, it is reasonable to choose from methods operating in the spatial domain [29–31], because among the existing alternatives their performance is the one least affected by correlation characteristics of noise.

Considering the peculiarities of the image forming process, we assume that the noise in the obtained images would contain a signal-independent component as well as a signal-dependent one, obeying the quasi-Poisson law. Thus, to describe the dependence of noise variance on the image intensity the following model can be used:

Here I_ij is the intensity value of the ijth image pixel $(i=\overline{1,N},j=\overline{1,M},$ where N and M are vertical and horizontal image sizes, respectively), ${\widehat{\sigma }}_a^2$ is the additive noise variance estimate, and k̂ is the quasi-Poisson noise parameter.

To evaluate the parameters of this model, it is logical to apply a  method based on the  scatter-plot approach. One of such methods operating in the  spatial domain is the  technique described in Ref. [31]. The main idea of the considered method is to get a scatter-plot of local variance ${\widehat{\sigma }}_{{\text{loc}}_m}^2$ and local mean Î_locm estimates obtained for quasi-homogeneous image blocks, then to estimate the centres of scatter-plot clusters (groups of estimates obtained for close local mean levels) and use these reference points to fit a polynomial curve into them. The parameters of the fitted curve are taken as the estimates of noise characteristics. In Ref. [31], it is recommended to use the double-weighted least mean square algorithm with constraints (DWLMSC), because it possesses a high enough flexibility and robustness.

The examples of quasi-homogeneous regions maps obtained using the  method described in Ref.  [31] for the images with z = 43 mm and z = 54 mm presented previously in Fig. 2 are shown in Fig. 4(a) and Fig. 4(b), respectively. The blocks belonging to detected quasi-homogeneous regions are marked with grey colour. Figures  4(c) and 4(d) show the  corresponding scatter-plots of local estimates with marked cluster centres and fitted polynomial curves.

Although during the visual inspection of the images presented in Fig. 2 it may seem that there are quite a lot of homogeneous regions, in fact, most of these areas turn out to be gradients or contain low contrast objects and thus are dropped by the detector. For most images, the amount of quasi-homogeneous regions is somewhere around 30–40%, which means that the number of blocks used in the process of evaluation of noise parameters is hardly enough for this estimation to be considered reliable. Another interesting feature of the images under study, that follows from the scatter-plots in Fig. 4, is their narrow dynamic range: for most images there are no clusters at local means above 150, which makes it difficult to determine the exact character of signal dependence of noise. To quantitatively measure the quality of curve fitting, we have calculated the values of several criteria that are traditionally used for this purpose: R-square and adjusted R-square. For all the images the values of both criteria turned out to be close to 0.99, thus model (1) describes the dependence of noise on the signal well enough and can be further used when choosing a denoising method and its settings.

The dependences of noise parameters evaluated under the assumption of model (1) on the distance from the source, at which the images were obtained, are shown in Fig. 5. As seen, both signal-independent and signal-dependent noise components show the  tendency of increasing as the  distance from the  source grows, although there are noticeable fluctuations in the values, especially for the quasi-Poisson component. These fluctuations appear due to an insufficient number of suitable blocks that could be found for some images, thus the estimates of noise characteristics for these images turned out to be biased.

Similarly to Fig. 3, plots in Fig. 5 can be divided into two areas. The  first area includes the  images obtained at z < 50 mm. In this area, the values of both noise parameters are relatively low and they slowly grow as z increases. However, starting from z = 50 mm, this growth becomes much faster and the overall level of noise in the images is significantly higher. Although this tendency is more noticeable for the additive noise component, it also exists for the signal-dependent one. However, it is just somewhat masked by the  outliers which appear due to the reasons explained above.

One more interesting observation from Fig. 5 is that the  higher z, the  greater the  contribution of the  additive noise component to the  total noise variance is, which is probably due to the slight misalignment of the  system when going further from the output mirror.

If we analyze the data in Figs 3 and in 5 jointly, we can notice that the noise in the images becomes less correlated as the distance z increases, while its variance, on the contrary, grows. Thus, we can note the existence of a certain value of z, which is quasi-optimal in terms of the visual quality of the obtained images. As follows from the  obtained results, this value is in the vicinity of the focus point. However, the  images obtained at the  quasi-optimal distance will not necessarily be the most informative, since it has been mentioned before that images obtained at different distances reveal different characteristic details and hidden features of the research object.

4. Description and performance analysis of the denoising method used

In the  previous sections, it has been shown that the images from the dataset under study are corrupted with the highly correlated mixed noise. Therefore, the denoising technique described in Ref. [17] seems to be the most adequate solution in this case. It is based on the use of discrete cosine transform (DCT) filtering with a hard spatial spectrum adapted (SSA) threshold. This method is also known as SSA DCT and it contains four main steps as described below.

Firstly, direct 2D DCT is performed for each 8 × 8 pixel block of the image. After this, one has an 8 × 8 array containing 64 coefficients, where the coefficient (0, 0) corresponds to the block mean and the rest of coefficients characterize different spatial frequencies. The obtained coefficients, except (0, 0), are thresholded (in the hard-thresholding mode all the values under threshold are replaced with zeros). The threshold T(m, k, l) is calculated according to the formula

where β denotes the filter parameter commonly set equal to 2.7, I̅_m is the local mean in the mth block, f(I̅_m) is the  function that describes the  dependence of the local variance on the local mean, and W(k, l) is the normalized power DCT spectrum obtained in 8 × 8 pixel blocks, k = 0…7, l = 0…7. After thresholding, the inverse 2D DCT is applied to each block. Since filtering is usually performed in fully overlapping blocks (positions of neighbouring blocks differ by only one pixel in the horizontal or the vertical direction), after step 3 there are usually 64 filtered values for a  single pixel. Thus, at the last step these values are averaged to obtain the final value for each pixel of the filtered image.

As mentioned in Section 3, for the dataset under study the dependence of noise variance on the image intensity can be described by model (1). Thus, function f(I̅_m) can be calculated as

According to Ref.  [32], the  normalized spectrum W(k, l) can be estimated in a blind manner as follows:

Here $D_{\text{H}}(k,l)=\sum _{b=1}^{B_{\text{HB}}}{D}_b(k,l)/B_{\text{HB}},D_b(k,l)\text{ are DCT }$ coefficients for an image block of size 8 × 8 pixels and B_HB is the number of overlapping 8 × 8 pixels blocks belonging to quasi-homogeneous image regions; here the maps of homogeneous regions obtained at the stage of evaluation of noise characteristics are used (see Fig. 4).

Examples of the normalized spectra W(k, l) for two images (with z = 43 mm and z = 54 mm) are presented in Fig. 6.

According to the results presented in Fig. 3, the noise in all the images is correlated. This is also confirmed by analysis of spectra in Fig. 6: for both images the spectra values are obviously non-equal and their energy is maximal for the coefficients corresponding to low spatial frequencies (small indices) and is minimal for the high ones. For the image obtained at z = 43 mm, most of the energy is concentrated in low spatial frequencies, which is typical for the images with highly correlated noise, while for the image obtained at z = 54 mm, the greater part of energy is distributed between the low and medium spatial frequencies, which is a typical spectrum shape for the noise with medium correlation.

Figure 7 shows the  examples of some images from the dataset under study filtered using an SSA DCT filter with the threshold calculated as described above. Comparing these pictures to the ones from Fig. 2, the following observations can be made.

For most images, the noise is suppressed effectively, while edges and features are also well preserved. In general, images after filtering are more pleasant visually and it is easier for a human operator to study small details in them, since smoother images are less irritating for eyes.

Images after filtering appear to be blurrier than their originals. There are several reasons behind this. Firstly, since it is very difficult to distinguish the  noise from the  signal, especially, if the  latter is spatially correlated, filtering inevitably affects the image content, which is usually more noticeable at high spatial frequencies (edges and textures). Thus, chosen filter settings are always a  compromise between a good noise suppression and preserving informational details. In this case, however, images look blurry mostly not because of being overfiltered but because they were very blurred initially, while the noise created a deceptive sharpening effect that got lost after the noise was removed. A possible solution to this problem may be the application of deblurring methods; however, this direction needs to be explored carefully further.

Some images, mostly the ones obtained at high z values, thus more corrupted with noise, look a bit underfiltered: some artifacts can be found near high contrast edges (see Fig. 6(e, f)). Such a situation is quite typical of the DCT filter, since it originates from the core idea that its operation is based on. At a high noise level, its separation from the signal in the upper spatial frequencies turns out to be practically impossible, therefore, when setting filtering parameters, one has to choose between signal over-smoothing and residual artifacts, while artifacts are usually more tolerable visually than excessive blurring.

To quantitatively measure the  quality of original and filtered images, we have applied the metric ILNIQE [33] that is considered one of the best modern no-reference image quality assessment metrics. The values of ILNIQE change from 0 to 150, while higher values indicate a worse visual quality of images. The average value of ILNIQE is about 124 for the  original images and about 114 for the  filtered ones, which confirms our previous conclusions about a  significant improvement of image quality achieved by means of filtering.

In general, the proposed method of automatic filtering of THz images seems to have a good potential in terms of their visual quality enhancement, yet there are several directions for its improvement.

5. Conclusions

Significant advances have been made in the field of terahertz imaging in the last decades, but there are still several problems that affect the quality of the resulting images negatively and are often difficult to solve at the hardware level.

The possibility of improving the quality of terahertz images obtained using a basic imaging setup operating at a frequency of 0.1 THz using modern image processing methods is considered.

A preliminary analysis of the statistical and correlation characteristics of noise in terahertz images has been carried out. It is shown that the noise in these images is spatially correlated, while the level of its spatial correlation decreases as the distance from the setup output mirror increases. The noise in terahertz images has a complex structure and is characterized by the presence of additive and quasi-Poisson components, the  magnitudes of which vary significantly for images obtained at different distances from the output mirror.

A denoising method based on the  use of SSA DCT is proposed. It implies the  automatic determination of image blocks belonging to quasi-homogeneous regions and using these blocks in the process of evaluation of statistical parameters of the  noise and its normalized DCT spectrum. Subsequently, the obtained information is used to calculate the thresholds for the SSA DCT filtre.

It is shown that the proposed denoising method provides a significant noise reduction with a good preservation of edges and details for most images, which makes them more suitable for visual analysis by a human operator. The average improvement of image quality achieved by means of filtering is about 7% according to the metric ILNIQE. Meanwhile, there are several directions in which further improvement of these images is desirable, in particular, increasing their clarity and sharpness.

Acknowledgements

This research was funded by a grant (No. P-LU-PAR-22-54) from the  Research Council of Lith-uania and supported by the  Marius Jakulis Jason Foundation within the special program ‘For Ukraine’s Science’.

References

TERAHERCINIŲ VAIZDŲ PAGERINIMAS NAUDOJANT DISKRETINE KOSINUSINE TRANSFORMACIJA PAREMTĄ TRIUKŠMO SUMAŽINIMO METODĄ

^a  Nacionalinio aerokosmoso universiteto Informacijos ir komunikacijos technologijų departamentas, Charkivas, Ukraina

^b  Fizinių ir technologijos mokslų centro Optoelektronikos skyrius, Vilnius, Lietuva

Santrauka

Žemą terahercinių (THz) vaizdų kokybę dažniausiai lemia aparatūriniai apribojimai, trukdantys naudingos informacijos vaizduose išskyrimui. Šio darbo tikslas buvo ištirti galimybes apeiti iškylančius apribojimus, pasitelkiant skaitmeninio vaizdų apdorojimo metodiką. Tyrimas atliktas naudojant vaizdų, užrašytų bandiniui esant skirtingu atstumu nuo 0,1 THz dažnio terahercinės spinduliuotės šaltinio, rinkinį. Darbe parodyta, kad gautuose THz vaizduose vyrauja kelių skirtingų prigimčių stipriai erdviškai koreliuotas triukšmas. Vaizdų kokybei pagerinti siūloma naudoti automatinį triukšmo sumažinimo metodą, paremtą diskretine kosinusine transformacija su erdviškai pritaikytu spektru. Parodoma, kad tokia metodika leidžia padidinti THz vaizdų kokybę efektyviai sumažinant pradinį triukšmą bei išsaugant ryškius vaizde esančių objektų kontūrus net ir esant žemai erdvinei skyrai ir dideliam pradiniam triukšmo lygiui, taip paruošiant vaizdus tolesnei vizualinei analizei.