# Beam-Shape Effects and Noise Removal from
THz Time-Domain Images in Reflection Geometry in the 0.25 - 6 THz Range

Marina Ljubenović, Alessia Artesani, Stefano Bonetti, and Arianna Traviglia
This work has been submitted to the IEEE for possible publication. Copyright
may be transferred without notice, after which this version may no longer be
accessible.This project has received funding from the European Union’s Horizon
2020 research and innovation programme under grant agreement No.
101026453.Marina Ljubenović, Alessia Artesani, and Arianna Traviglia are with
the Center for Cultural Heritage Technology, Istituto Italiano di Tecnologia,
Venice, Italy.Stefano Bonetti is with the Department of Molecular Sciences and
Nanosystems, Ca’Foscari University of Venice, Venice, Italy and the Department
of Physics, Stockholm University, Stockholm, Sweden.

###### Abstract

The increasing need of restoring high-resolution Hyper-Spectral (HS) images is
determining a growing reliance on Computer Vision-based processing to enhance
the clarity of the image content. HS images can, in fact, suffer from
degradation effects or artifacts caused by instrument limitations. This paper
focuses on a procedure aimed at reducing the degradation effects, frequency-
dependent blur and noise, in Terahertz Time-Domain Spectroscopy (THz-TDS)
images in reflection geometry. It describes the application of a joint
deblurring and denoising approach that had been previously proved to be
effective for the restoration of THz-TDS images in transmission geometry, but
that had never been tested in reflection modality. This mode is often the only
one that can be effectively used in most cases, for example when analyzing
objects that are either opaque in the THz range, or that cannot be displaced
from their location (e.g., museums), such as those of cultural interest.
Compared to transmission mode, reflection geometry introduces, however,
further distortion to THz data, neglected in existing literature. In this
work, we successfully implement image deblurring and denoising of both
uniform-shape samples (a contemporary 1 Euro cent coin and an inlaid pendant)
and samples with the uneven reliefs and corrosion products on the surface
which make the analysis of the object particularly complex (an ancient Roman
silver coin). The study demonstrates the ability of image processing to
restore data in the 0.25 - 6 THz range, spanning over more than four octaves,
and providing the foundation for future analytical approaches of cultural
heritage using the far-infrared spectrum still not sufficiently investigated
in literature.

###### Index Terms:

THz-TDS, Reflection Geometry, Image restoration, Deblurring, Denoising,
Cultural Heritage

##  I Introduction

Spectroscopy and imaging in the terahertz (THz) spectral range (typically 0.1
to 10 THz) is gaining attention in many disciplines such as biomedicine,
agriculture, security, and communication services [1, 2]. One of the major
applications of THz spectroscopy and imaging is in noninvasive detection of
material composition and retrieval of complex surface topography of objects
[3, 4, 5]. A THz signal exhibits many unique properties: i) its low-photon
energy facilitates nondestructive and non-ionizing testing [6], ii) it has the
ability to penetrate nondielectic materials (e.g., paper, plastic, etc.); iii)
it can be used in reflection geometry for studying objects that are too thick
for signal transmission, or highly reflecting (e.g., objects made of metal).

Terahertz Time-Domain Spectroscopy (THz-TDS) employs short pulses of THz
radiation to probe dielectric response of materials in the far-infrared region
and reconstruct objects’ inner structure, giving as result a 3-dimensional
hyperspectral (HS) data-cube, providing information on both surface and inner
structures. Despite the fact that THz-TDS imaging offers a unique means to
probe materials’ properties, it poses also major challenges in terms of data
processing because of the multiple sources of image degradation that make
difficult the analysis of the THz signals [7]. THz-TDS images are in fact
corrupted by frequency-dependent beam-shape induced blur and noise with
additional signal distortions introduced by scattering losses (more
information on degradation effects will be provided in Section IV). All these
image degradation effects remarkably reduce the frequency region of the
resolved images and limit the recovering of sharp THz images.

Some conventional image restoration algorithms can singularly tackle blur
effects, noise, or low resolution, thus overcoming only partially physical
limitations determined by the THz-TDS imaging technique. The joint removal of
several degradation effects in THz HS images represents still a challenge due
to several reasons: a frequency dependency of blur and noise; the presence of
blocking artefacts introduced by the step size of a scanning system; and
existence of additional effects such as reflection and refraction [7]. At
present, there is a limited number of efficient computational strategies to
address some of the THz imaging problems and the majority of the proposed
approaches have dealt with single-frequency measurements, such as THz images
in continuous wave mode or volumetric reconstruction with THz computed
tomography.

THz-TDS imaging can be carried in transmission and reflection geometries. In
transmission geometry, the temporal profile of the electric field transmitted
through the investigated sample is measured and compared with a reference
spectrum obtained from the free path. However, many samples possess
significantly high absorption properties in the THz range or are highly
reflective, situations that limit or completely impede the use of THz imaging
in transmission geometry. Furthermore, for samples with a complex structure,
absorption of the radiation in the high frequency range limits the bandwidth
to 2-2.5 THz in most cases. All these issues are minimized or completely
removed performing the measurements in reflection geometry, which allows to
measure virtually any sample, and with a much larger bandwidth. Furthermore,
it has been demonstrated that reflection imaging is more sensitive to sample
adhesion (e.g., deviation in thickness of paint and coating) that can cause
variation in the phase of the THz signal [8]. However, this greater
sensitivity is also a potential source of new degradation artefacts in the
resulting image, such as discontinuous spots. Other degradation effects that
could be present in reflection geometry, and are not present in the
transmission one, includes distortion and signal losses due to refraction,
diffraction, edge effects, specular reflection and Fabry-Perot effects, which
might heavily corrupt the THz images. It follows that in order to fully
exploit the great potential of THz imaging in reflection geometry, specific
corrections of these degradation effects need to be implemented.

This work describes the application to THz-TDS images of an algorithm that has
been recently developed and successfully applied for joint THz image
deblurring and denoising [7]. The algorithm was originally developed for THz
imaging in the transmission mode and validated on objects having a simple
shape (e.g., leaf, ring, plate). Here, it is adapted for the first time for
THz reflection images and tested on samples with more complex material
structures and textures. Section II starts by introducing existing image
processing methods that are developed or applied to THz data and their
limitations. The experimental setup is introduced in Section III, followed by
the introduction of degradation effects in Section IV. The efficiency of the
joint deblurring and denoising approach is tested on highly reflective objects
(Section V), and then applied to the study of an ancient coin (Section VI).
The joint deblurring and denoising approach, tailored to frequency- dependent
blur degradation, utilizes the high correlation of spectral bands to increase
their availability for the further analysis.

##  II Image Processing Methods for THz Data

Initial attempts of THz image restoration have employed conventional image
restoration methods originally developed for other imaging modalities (e.g.,
X-ray, RGB images) [8]. These include but are not limited to iterative
backprojection, Richardson-Lucy method [9],[10], and 2D wavelet decomposition
reconstruction [11]. For instance, Li et al. in [12] applied a well-known
Richardson-Lucy deblurring algorithm to remove beam shape effects from a
single-frequency THz images improving the boundary sharpness of binary
samples, but introducing artefacts in the presence of fine details. Those
conventional deblurring algorithms have been designed to improve resolution of
images that follow different statistics: they are demonstrated to be
inadequate for THz data mostly due to the huge difference in the assumed prior
knowledge imposed on a sharp image and the degradation effects (e.g., blur and
noise).

At first, the development of procedures aiming at improving THz data was
mainly focused on specific applications such as THz-based computed tomography
(THz-CT) or subwavelength microscopy. THz-CT data often require sophisticated
algorithms to manage volume reconstruction from a projection domain. One of
the first methods tailored to THz-CT was introduced by Recur et al. in [13].
The method combines the CT reconstruction approaches (e.g., back-projection of
filtered projections, simultaneous algebraic reconstruction technique, and
ordered subsets expectation maximisation) with a convolution filter to mimic
the THz beam effects. The influence of these effects can be further removed by
employing the Wiener deconvolution approach [14, 15]. Some attempts were made
for increasing the resolution of images with a multilayered structure by
applying deconvolution in the time domain [16]. These methods can tackle
noise-free single-frequency THz images or signals, but do not take into
account the differences imposed by the dependence of the noise and blur
effects on the frequency value. These methods are not suitable for a broadband
THz imaging systems and the complex structure of its beam as they are
considering images corresponding to a single THz frequency as independent. To
beat the diffraction limit in subwavelength microscopy, a super-resolution
technique for THz images was proposed in [17]. This technique utilizes
specially designed artificial fluorophores in the form of optimal mask sets
and stochastic fluctuations in the intensity of the fluorophores limiting its
broad usage.

The complex structure of a broadband THz beam and its effects was studied in
[18]. The authors modelled a 3D THz point spread function (PSF) and THz
imaging equation, further utilized to enhance single-frequency THz images by
applying a simple deconvolution approach (i.e., Matlab in-built function,
deconv). Utilisation of a more advanced deconvolution procedure for tackling
blur degradation in frequency modulated continuous wave THz images was
introduced by Wong et al. in [19]. This paper introduced two steps leading to
THz image enhancement: high-precision depth estimation followed by 2D
deconvolution with off-the-shelf methods [20]. Although these approaches
worked reasonably well for noise-free images and simple-shape samples, due to
the ill-conditioned nature of the convolution, they introduced severe error in
a presence of even small amount of noise.

Deep neural networks (DNNs) have also been used to ameliorate THz data using
robust restoration approaches, but rarely considering the frequency-dependence
of the THz-TDS images during the training process. They were mostly used for
2D THz image deblurring [21] and super-resolution [22, 23] or for specific
applications, such as biological product analysis [24]. The main limitation of
DNN-based methods is the need for a large number of THz images (hundreds or
thousands) for accurately training the network: unavailability of a sufficient
number of THz images has, as a consequence, limited their application.

One of the domains in which the removal of blur and noise degradation of HS
images have been significantly developed is that of airborne Remote Sensing
imagery: here, the clarity of the HS images is key for their interpretation
and the quality of the spectral datum pivotal to further image processing.
This has led to the development of several methods to improve the quality of
the image, such as denoising (HySime [25] and FastHyDe [26]) and deblurring
methods [27] that leverage some characteristics of the HS images like low-rank
and self-similarity of image patches. In particular, FastHyDe demonstrated to
be especially successful in denoising HS images affected by types of noise
such as Gaussian independent and identically distributed (i.i.d.), Gaussian
non-i.i.d., and Poisson noise. Notwithstanding their potential transferability
in other domains, these image processing approaches have not found so far
application in close-range THz imaging.

An approach that considers bands of THz HS images jointly has been introduced
in [7]. Similar to FastHyDe, the method is based on the low-rank and self-
similarity properties of HS images and it includes two additional expansions
that make it suitable for THz data: i) a specialized THz beam modelling
process and ii) a deblurring step based on well-known, robust deblurring
approaches. The joint deblurring and denoising approach represents an
effective measure to counter frequency-dependent blur and three different
noise types present in time-domain THz images acquired in the transmission
mode, but, to the authors’ knowledge, it has never been applied in the
reflection mode and on objects with irregular geometries. Our work focuses on
demonstrating the application of the HS algorithms for THz imagery on uneven
object samples, as an archaeological coin. The analysis on these kind of
samples is complicated by the irregular surface, the asymmetry of their
thickness, the presence of irregular reliefs on the samples, the variability
of their material compositions (often unknown and/or multiform) and
contribution of both reflecting and absorbing materials.
Our research has entailed the application of a joint denoising and deblurring
approach to improve the quality of THz images in reflection geometry: the
paper critically discusses the acquisition procedure, the image restoration
requirements, and the limitations of the proposed framework. The degradation
effects induced by the THz beam size are evaluated in order to determine the
point-spread function (PSF) of the THz-TDS system necessary to remove image
degradation. Band-by-band deblurring, HS image denoising, and a method that
remove blur and noise jointly are initially tested on objects characterized by
uniform material composition and relatively unpatterned shape (a 1 Euro cent
coin –€0.01– and an inlaid silver pendant). After demonstrating that only the
joint deblurring and denoising method achieves good results through all bands,
and after evaluating its limitations and advantages, the method is directly
applied to an ancient Roman silver coin and the differences between the two
types of samples (i.e., contemporary and archaeological objects) are
highlighted.

##  III Experimental Setup and Acquisition Procedure

Figure 1: Scheme of the THz-TDS system in the reflection mode. The instrument
is based on an ultrafast laser source, emitter and receiver antennas and
motorized raster-scanning stage, where the sample is positioned for
measurements in reflection geometry.

This work has been completed using a commercial THz-TDS system (TOPTICA
TeraFlash Pro, Fig. 1). The system is based on a mode-locking Erbium-doped
fiber laser ($\lambda$ = 1560 $\pm$ 10 nm), with pulse width of 50-60 fs and
repetition rate of 100 MHz (TOPTICA FemtoFErb THz FD6.5). The ultrafast laser
pulse is delivered via a single-mode fiber onto an InGaAs/InAlAs-based
photoconductive antenna that in turn produces an electromagnetic pulse with
the frequency spectrum in the THz range. The system is capable of generating a
spectrum width that spans 6 THz, has the highest power being emitted at
frequencies around 1 THz. The precise reconstruction of the time axis and the
reduced jitter noise of the system enables to achieve a good dynamic range
below 3.5 THz (in nitrogen atmosphere), with a peak dynamic range of more than
90 dB. The measured THz electric field is Fourier transformed to yield the
phase and the amplitude of the THz pulse.

The system relies on point-like measurements. The analysis of a surface is
obtained in the reflection mode by raster-scanning the sample placed in the
focal plane of the two parabolic mirrors, using a motorized XY stage that is
synchronized with the THz excitation. The image is obtained by continuous line
scanning along the X direction (from left to right) with a step size of
$\Delta$x, and by step line along the Y direction of $\Delta$y. For the
TOPTICA TeraFlash Pro, the step sizes in both directions can be fixed equal to
0.1, 0.2, 0.5, 1 or 2 mm. Whereas the frequency resolution of the THz‐TDS is
given by the time step and the number of data points in the time‐domain
measurement, the spatial resolution is limited by the point spread function
(PSF) of the THz beam (as discussed in Section IV), and also by the
acquisition step size along the X and Y axis.

Figure 2: Example of the amplitude signal $A(\omega)$ of a raster-scanned
sample taken with three different $\Delta$x and $\Delta$y step sizes, equal to
(A) 0.5 mm, (B) 0.2 mm, (C) 0.1 mm. The frequency is fixed at 2.72 THz.

In order to show the difference between the step size choice, a portion of a 1
Euro cent coin has been raster-scanned with step size set equal to 0.5 (A),
0.2 (B), and 0.1 (C) mm, respectively (Fig.2). With the biggest step size, the
detector loses information leading to severe blocking degradation in the
resulting image. Due to this degradation, smoothness of the object edges is
lost. With the smallest step size, the THz images reaches a higher level of
details, but the acquisition procedure lasts significantly longer (e.g., 3 min
(B) vs. 32 min (C) for scanning $5\times 5$ mm) and a data size increases
($\sim 12$ MB (B) vs. $\sim$ 120 MB (C)), while missing in improving
significantly the spatial resolution due to the beam shape effects. For the
considered sample, the optimum between the choice of the step size and the
acquisition time is the option (B), whereas the value of the step size should
be carefully considered taking into account aspects such as dimension of the
samples and its features, and the desired application results.

##  IV Degradation Effects

Figure 3: Typical degradation effects in the representation of the amplitude
$A(\omega)$ images. The images are corrupted by noise that increases with the
frequency, while at lower frequency value, the blur effects are more
dominant.The frequency region with relatively low noise and blur degradation
effects stays between 0.5 and 3.5 THz.

The THz system design, with focusing optics and beam forming process, leads to
frequency-dependent degradation effects that corrupt the amplitude images. The
main causes for the degradation effects that can alter the spatial information
of the THz images are mainly: i) the THz beam-shape; ii) the system-dependent
noise; iii) losses due to the propagation of the THz wave through the sample
(caused by reflection and/or refraction); iv) the scattering and absorption
properties of the sample. Whereas the latter two degradation effects are
intrinsically dependent on the sample and hard to model, the consequences of
the THz beam-shape and the noise of the system can be estimated by introducing
simple approximations. Thus, the restoration algorithms are based on the
following assumptions for the blur and noise effects.

Blur effects. The intensity profile of the THz beam is assumed to have a
Gaussian-shape [13]. Under this approximation, the transverse profile of the
THz intensity is described by

$I(r,z)=\frac{P}{\pi
w^{2}(z)/2}\text{exp}\left(-2\frac{r^{2}}{w^{2}(z)}\right),$ (1)

where $P$ is the beam power, $w(z)$ is the beam radius, $z$ is the propagating
direction and $r$ is the distance from the beam axis. The radius of the beam
varies along the propagating direction as

$w(z)=w_{0}\sqrt{1+\left(\frac{z}{z_{R}}\right)^{2}},$ (2)

where $w_{0}$ is minimum beam radius (half of the beam waist) and $z_{R}=\pi
w_{0}^{2}/\lambda$ represents the Rayleigh length at wavelength $\lambda$ that
determines the length over which the beam can propagate without diverging
significantly. In Gaussian beam approximation, the beam waist $2w_{0}$ depends
on the wavelength $\lambda$ accordingly to the [28]:

$2w_{0}=\frac{4}{\pi}\lambda\frac{f_{L}}{D},$ (3)

where $f_{L}$ and $D$ represent a focal length and a diameter of the focusing
THz lens. The size of the THz beam waist is strongly dependent on the
frequency and the focusing optics of the THz system. The first degradation
effect associated to the beam-shape consequently affects the resolution of the
HS image and introduces blurring. In particular, every pixel in the image
spatial domain is blurred following a 2D Gaussian PSF corresponding to an
intersection of the THz beam with an orthogonal plane (represented by the
surface of the sample). As the THz beam waist is wider at lower frequencies,
the blurring effects will be more important in the lower region of the THz
frequency range, as illustrated in Fig. 3. On the contrary, at higher
frequencies, a THz beam waist is smaller and the resulting images are sharper.

Noise effects. The principal degradation effect at high frequency is the noise
determined by the instrument. The signal from a THz-TDS system is corrupted
by, at least, three noise components [29] [30]. The first noise component is
caused by the transmitting antenna, which represents the principal noise
source. The second contribution is Poisson noise in the detector, which
usually originates from the discrete nature of the electric charge. The third
component summons other signal-independent noise in the detector such as
amplification noise, laser noise, thermal noise, etc. All these noise
components are severely decreasing the signal-to-noise ratio with the increase
of frequency, making the images corresponding to higher frequency bands
impractical for further analysis. For the instrument employed, this occurs for
frequencies above 4.5 THz (see Fig. 3).

In order to validate the above-mentioned Gaussian modelling assumption
implemented in the algorithm for working in reflection geometry, synthetic
reference images of a circle blurred by a Gaussian model have been compared
with the real bands of THz time-domain images of a hole on a metallic plate.
The theoretical THz beam waist, calculated considering that the instrument has
two focusing mirrors with 1-inch diameter and 4-inches length, is
$w_{0}\approx 2.547\lambda$. The experiments were performed considering step
size $\Delta x$ and $\Delta y$ equal to 0.2 mm, and beam diameter equals to
the beam waist. The results are presented in Fig. 4. On the lower and medium
frequencies (i.e., 0.97 and 1.94 THz), the difference in intensity between
synthetic and real bands is mostly visible at the borders of the circle, while
on the higher frequency (i.e., 3.11 THz) the small difference is additionally
visible in the center, this being most likely due to the influence of noise.
Thus, the performance of the overall method is mostly reduced near the edges.
The test shows the theoretical PSFs (intersection of an orthogonal plane with
the beam calculated by applying the equations (1) and (2)) are comparable with
the real blurring degradation.

Figure 4: Comparison between simulated reference images and experimental
measurements of a hole on a metallic plate at three THz frequencies. From left
to right: the denoised bands of raw data on three different frequencies (0.97,
1.94, and 3.11 THz); simulated bands corresponding to the same frequencies;
intensity vectors corresponding to the middle cross-section of real (blue) and
simulated (red) bands; and theoretical PSFs.

##  V Image Restoration for THz Images

The goal of this work is restoring the amplitude images, which have the form
of a hypercube with the third dimension corresponding to the THz frequencies.
These amplitude images are modelled as:

$\textbf{Y}=\textbf{X}\textbf{H}+\textbf{N},$ (4)

where $\textbf{Y}\in\mathbb{R}^{b\times n}$ represents the degraded HS
amplitude image with the rows containing $b$ spectral bands. Every band is an
image with $n$ pixels corresponding to the amplitude intensity.
$\textbf{X}\in\mathbb{R}^{b\times n}$ represents an underlying clean HS image
(without noise and blur effects), while $\textbf{H}\in\mathbb{R}^{n\times n}$
is the blurring operator that models the convolution of a single band and the
corresponding PSF. In the defined model, the same blurring operator is assumed
over bands. Finally, $\textbf{N}\in\mathbb{R}^{b\times n}$ represents additive
white Gaussian i.i.d. noise.

Prior to the application of the joint denoising and deblurring approach, the
restoration of the degraded images is undertaken by conventional band-by-band
deblurring (e.g., Richardson-Lucy or Wiener filtering) and denoising
(FastHyDe) algorithms to demonstrate how the two algorithms only partially
improve the image quality. The joint application of the denoising and
deblurring approach is then described, taking into account the variable
effects for each band. Noteworthy, the joint deblurring and denoising
algorithm does not simple combine the two methods (i.e., FastHyDe and
Richardson-Lucy), but follows a number of carefully designed steps, described
in Subsection V-CC.

###  V-A Band-By-Band Deblurring Algorothm

A band-by-band deblurring is performed to remove blurring effects from THz-TDS
images, where the deblurring problem for one band is modeled as

$\textbf{y}_{i}=\textbf{H}\textbf{x}_{i}+\textbf{n}.$ (5)

Here, $\textbf{y}_{i}$, $\textbf{x}_{i}$, and n represent an observed
vectorized image of the $i$-th band ($i$-th row of Y), a sharp image of the
$i$-th band, and noise respectively.

Three conventional approaches are tested on the 1 cent coin and the silver
pendant to perform band-by-band deblurring, specifically: 1) Richardson-Lucy
deblurring [9], [10]; 2) Krishnan et al.’s deblurring with a prior knowledge
depicting heavy-tailed distribution of gradients in natural scenes (hyper-
Laplacian prior) [31]; 3) Deblurring based on Wiener filtering [14]. These
three approaches are selected as the ones most commonly used in the literature
for THz beam shape effects removal.

Figure 5: Results of band-by-band deblurring approaches for restoring
amplitude images on 1 Euro cent coin (top) and a silver pendant (bottom). The
first rows corresponding to each sample show the optical image of the sample
and selected row bands on seven different frequencies (red dashed line). The
tested deconvolution approaches: 1) Richardson-Lucy iterative deconvolution
[9],[10] (blue dashed line); 2) Fast image deconvolution with a hyper-
Laplacian priors by Krishnan et al. [31] (green dashed line); 3) Image
deblurring using Wiener filtering [32] (yellow dashed line).

The results obtained by three conventional deblurring approaches are presented
in Figure 5 by displaying seven selected bands of THz time-domain images
corresponding to lower frequencies (0.28 and 0.38 THz), medium frequencies
(1.1, 2.1, and 3.3 THz), and high frequencies (4.87 and 5.85 THz). The input
parameters (the number of iterations, regularization parameters, etc.) are
selected according to values reported in the original papers and additionally,
when necessary, are adjusted to provide the best visual results.

The tested methods show some improvements for the selected medium frequencies
and poor results when tested on bands corrupted by severe blur (the band
corresponding to the low frequencies) or severe noise (the band corresponding
to the high frequencies). The additional limitation when applying 2D-based
deblurring approaches is hand-tuning of input parameters for every band
separately (e.g., number of iterations for the Richardson-Lucy approach or a
regularization parameter for the approach from Krishnan et al.). This is
because these band-by-band deblurring methods have been tailored to
conventional RGB images and thus, in the original papers, values for input
parameters are governed by the statistics of natural images, which are quite
different than the one of THz images. In addition, these methods usually
require several input parameters that are different for different HS image
bands and noise levels and therefore, in the case of THz HS images,
challenging to properly set.

Figure 6: Results of the FastHyDe method for restoring amplitude images on 1
Euro cent coin (top) and a silver pendant (bottom).

###  V-B Denoising Algorithm

In order to deal with bands severally corrupted by strong noise, a state-of-
the-art HS image denoising method, FastHyDe, proposed by Zhuang et al. in [26]
is adjusted and tested on THz-TDS images. HS denoising problem is formulated
following (4) with $\textbf{H}=\textbf{I}$, where I represents the identity
matrix.

Considering the high correlation between bands of THz-TDS images, same as in
FastHyDe, the following assumption is considered: the spectral vectors
$\textbf{x}_{i}$, for $i=1,...,n$, live in a $p$-dimensional subspace
$\mathcal{S}_{p}$, where $p\ll b$. Following this assumption and defining
$\textbf{E}=[\textbf{e}_{1},...,\textbf{e}_{p}]\in\mathbb{R}^{b\times p}$ as a
basis for $\mathcal{S}_{p}$, X can be written as

$\textbf{X}=\textbf{E}\textbf{A},$ (6)

where $\textbf{A}\in\mathbb{R}^{p\times n}$ corresponds to the representation
coefficients of X in $\mathcal{S}_{p}$. In FastHyDe, the rows of A are called
eigen-images. The above assumption on the high correlation between bands
(i.e., the assumption that HS data follow a low-rank structure) is crucial, as
matrix E may be learned directly from Y by subspace identification methods
such as HySime [25]. Bands of HS images (including THz time-domain images) are
self-similar in the spatial domain, i.e., they contain many similar patches at
different locations and scales. The exploitation of self-similarity as a form
of prior knowledge is the second main assumption of the FastHyDe method
applied through the exploitation of a patch-based denoising step. The main
steps of FastHyDe are shown in Algorithm 1.

Algorithm 1 FastHyDe

1) Learn the subspace E from Y

2) Compute noisy eigen-images $\textbf{E}^{T}\textbf{Y}$

3) Denoise eigen-images by a state-of-the-art denoiser [33]

4) Compute an estimate of the clean HS image

The results of FastHyDe applied to all bands jointly of two real THz time-
domain images are displayed in Figure 6. Note that only selected bands are
showed. The results are obtained by assuming Gaussian i.i.d. noise and the
number of subspaces $p=10$. The tested method is able to deal even with severe
noise that corrupt bands corresponding to higher frequencies (e.g. frequencies
higher than 4 THz). At lower frequencies, the output is comparable to the
input, as the tested method only removes noise and not beam shape effects.

###  V-C Joint denoising and deblurring Algorithm

The fast deblurring and denoising method [7], tailored to THz time-domain
images, three noise types (Gaussian i.i.d., Gaussian non-i.i.d., and Poisson),
and band-dependent blur is further tested on 1 cent and a silver pendant. The
method follows the same low-rank assumption from (6) but instead of performing
only denoising on subspace components (i.e., eigen-images) it additionally
adds a non-blind deblurring step (i.e., deblurring with a known PSF). Note
that deblurring in this case is performed only on subspace components and not
on all bands separately (as in case of methods presented in Subsection V-AA).
Consequently, the joint deblurring and denoising method consists of steps
presented in Algorithm 2.

Algorithm 2 Joint deblurring and denoising of THz time-domain images

Dimensionality reduction

1) Learn the subspace E from Y

2) Compute noisy eigen-images $\textbf{E}^{T}\textbf{Y}$

Blur and noise removal

3) Analyse subspace domain

4) Create synthetic PSFs

5) Deblur eigen-images with an off-the-shelf non-blind deblurring method

6) Denoise eigen-images by a state-of-the-art denoiser from [33]

Reconstruction

7) Compute an estimate of the clean HS image.

Figure 7: Results of the joint deblurring and denoising method for restoring
amplitude images, as described in Algorithm 2. The raw results are compared
with the images obtained after restoration for single THz frequency values.
The method is applied on two metallic samples: 1 Euro cent coin (top) and one
silver pendant (bottom).

The results obtained with one instance of the joint deblurring and denoising
method with Gaussian i.i.d. noise assumption and Richardson-Lucy deblurring
step (Step 5 in Algorithm 2) are shown in Figure 7. With respect to the
transmission mode, here a new step is included covering detailed analysis of
the subspace domain and corresponding eigen-images (Step 3 in Algorithm 2). In
[7], only samples with simple shapes and textures are tested thus, the number
of subspaces can be fixed for all samples. Here, the number of subspaces is
carefully chosen to include all significant features (e.g., main edges and
texture) and exclude features corresponding to reflection effects (see Fig.
8). Similarly, the values of the minimum beam radius utilized for PSF
modelling is carefully set by analysing blur degradation of eigen-images,
i.e., eigen-images corresponding to main edges (e.g., images 1, 2, and 3 in
Fig. 8) are deblurred by using bigger PSFs relative to the eigen-images
corresponding to textural features (e.g., images 5, 6, and 7 in Fig. 8). The
subspace components corresponding to reflection effects (images 10, 11, and
12) are discarded during the analysis and thus, the number of subspace
components for the 1 cent sample is set to $p=9$. The similar analysis of
subspace components is repeated for each sample.

Figure 8: Subspace components of 1 Euro cent ($p=12$): Subspace components
from 1 to 3 cover the main edges; components from 5 to 7 correspond to
textural features; components 4, 8, and 9 account for flat regions; component
from 10 to 12 are showing effects of reflection.

As the tested joint deblurring and denoising method is based on FastHyDe, the
noise from bands corresponding to higher frequencies is removed equally
successful as previously. Additionally, beam shape effects removal is visible
in bands corresponding to lower frequencies (i.e., frequencies less than 1
THz).

The performance of the joint deblurring and denoising method applied to images
acquired in reflection geometry was measured by comparing standard deviation
of an even region (Fig. 9) and feature sharpness (Fig. 10) of raw and restored
bands. The standard deviation was estimated over a $2\times 2$ mm homogeneous
region selected from the 1 euro cent coin. For each pixel, the variance of
reflected reference intensity (raw data) and reconstructed intensity
(estimated data) is computed and compared. The improvement resulting from the
restoration procedure is visible over all bands, increasing with the increase
of frequency. Similarly, a feature sharpness (in mm) is estimated by measuring
a distance between a high and a low peak that represent a border of the
feature. Decrease in distance indicate the steepest feature border and thus, a
sharper feature. The distance measure is highly affected by the presence of
noise, especially visible on the higher frequency bands.

Figure 9: Performance of the joint deblurring and denoising method measured
through standard deviation. (a) A selected flat region from raw and restored
bands for calculated standard deviation; (b) log standard deviation as a
function of different frequencies. Figure 10: Performance of the joint
deblurring and denoising method measured through a sharpness of a feature (the
number one on the coin). (a) Restored band at 3 THz and its cross section; (b)
Feature sharpness estimated as a distance between a high and a low peak.

Both visual and the numerical results (i.e., standard deviation and feature
sharpness) respectively, demonstrate that the joint deblurring and denoising
approach achieves significantly improved THz HS images of homogeneous
materials without introducing additional distortions. The reconstruction
approach was able to restore both lower (e.g., below 1 THz) and higher (e.g.,
above 3.5 THz) spectral components, making them more suitable for further
analysis. Successful restoration of bands corresponding to the broad THz
spectral range drastically increase the amount of information on a sample,
enabling both spatial and spectral domains to be fully utilized in different
applications such as cultural heritage, biomedical imaging, quality
inspection, etc.

##  VI Application of the method to the study of an ancient coin

Amelioration results obtained on THz-TDS images of the two objects show the
potential of the image restoration procedures in reflection geometry. In
particular, they enable to retrieve clean images at very high (above 3.5 THz)
and very low (below 1 THz) frequencies, otherwise not accessible due to heavy
noise and blur effects. This represents a highly desirable result, especially
when dealing with material characterisation. To test the capacity of the
algorithm on more complex objects (in terms of shape and material
composition), the archaeological coin showed in Fig. 11 is utilized.

Figure 11: Stereoscopic images of the archaeological silver coin. On the front
side the figure of a goddess is inlaid (A), while on the back side the one of
a Roman emperor portrait (B).

The ancient coin used in this study has worn-out reliefs, a surface altered by
the ageing process, and is encrusted with soil residue. The opacity to the THz
radiation of this type of object requires the acquisition of THz images in
reflection geometry, but the presence of alteration layers on the surface can
induce multiple reflections and can heavily affect the amplitude signal. The
analysis was undertaken in the same way as for the contemporary objects
illustrated in the previous sections and both sides of the Roman coin were
included.

Figure 12: Results of the joint denoising and deblurring method on the front
and back faces of the archaeological Roman coin. The results are showed as
integrated amplitude images over the frequency range which starts from the
frequencies presented above the images + 0.4 THz. Figure 13: False RGB image
of raw data and after image restoration, constructed from three integrated
amplitude images. The red channel (R) corresponds to 0.4-0.8 THz, which is a
highly blurry frequency region, the green channel (G) is integrated amplitude
between 1.7-2.1 THz, which shows a good quality result and high dynamic range,
and lastly the blue channel (B) corresponds to 4.5-5.5 THz, more affected by
noise. The false RGB image after restoration delineate better the distribution
of the amplitude intensity in the three frequency regions.

The image restoration approach used in this experiment is the one presented in
Section V-C with assumption of Gaussian i.i.d. noise. The number of subspaces
and $w_{0}$ corresponding to PSFs are chosen by analysing subspace components
as explained in Subsection V-C. As for the other samples previously
investigated, the conventional Richardson-Lucy algorithm is used as a
deblurring step and the state-of-the-art patch-based method from [33] as a
denoising step. The integrated amplitude values for six frequency ranges of
raw and restored THz HS images are presented in Figure 12: the high variety
between bands (from blurrier to noisier ones) demonstrate the need for a
robust approach when dealing with cultural heritage objects. Compared to
contemporary samples, that are made by known and homogeneous composition,
Roman coins are made of bronze or silver alloy, which developed metallic oxide
on the surface. The spectral signatures of these compounds can be highly
challenging to determine in the THz domain. Thus, the results in Figure 12 are
represented as an integrated amplitude over a range of frequencies (starting
from the frequencies presented in the figure plus 0.4 THz). The restoration
procedure results in sharpening the bands corresponding to lower frequencies
(e.g., 0.38-0.78 THz) and fully removing severe noise from bands corresponding
to higher frequencies (e.g., 5.85-6.25 THz), thus expanding the usable
frequency range for further analysis.

The highly modular and robust approach for joint deblurring and denoising
utilized in this work for the analysis of objects with irregular shape and
complex surface texture provides highly promising results. The modularity of
the approach, reflected in its step-like design, ensures a better evaluation
of the performance of each separate step and therefore an easier tuning of
input parameters. The patch self-similarity assumption makes the approach more
robust to irregular shapes of the sample, representative of the shapes
normally found in cultural heritage items: the approach is focused to image
patches (small image parts) rather than the whole image, thus separately
capturing wider complexity of shapes and textures.

The final result is represented by a false-color image, where three amplitude
images at very low frequency (0.4-0.8 THz), in the medium region (1.7-2.1 THz)
and at very high frequency (4.5-5.5 THz), are merged, as illustrated in Fig.
13. Outputting a false-color image is quite common in HS image analysis,
however in this case THz-TDS raw results could in principle induce to wrong
conclusions because the presence of blurring and noising effects can alter the
features of the objects. After the image restoration, the resolution of the
images become similar and the distributional map more reliable. This work does
not include comparison of the results to a reference (e.g., a spectra from a
spectral signature library), and thus the features presented could be
topographical.

##  VII Conclusion

The paper demonstrated the positive outcomes achieved using a computational
procedure for THz-TDS imaging to restore amplitude signals in reflection
geometry and to tailor it to deal with irregular surfaces and uneven reliefs,
such as the ones of the archaeological item illustrated in this work.

Utilizing image restoration approaches that are able to handle strong noise is
crucial when applying THz time-domain imaging to objects that have
distinguishable spectral features in the middle and high frequencies (i.e.,
frequencies corrupted by such noise). The importance of applying image
processing strategies to HS images is most evident for applications that
relies on the possibility of extension of the available frequency range using
far-infrared frequencies (spanning from few hundreds of GHz to 10 THz)
enabling a more comprehensive analysis. This study broke the 3 THz limitation,
defined by the previous studies [5], and showed reliable image reconstruction
over more than 4 octaves, from 0.28 THz and reaching up to 5.85 THz.

Further work in this field will have to entail the development of a novel THz
super-resolution approach to overcome the resolution limitation partially
imposed by a step size of the instrument during acquisition and significantly
increase the size and number of samples suitable for analysis.

## Acknowledgment

Images of the ancient Roman coin are courtesy of the Soprintendenza
Archeologia, belle arti e paesaggio del Friuli Venezia Giulia (prot. SABAP
N.0004987-P, 30/12/20), which the authors wish to thank for the kind support
and availability.

## References

  * [1] S. Wietzke, C. Jansen, M. Reuter, T. Jung, D. Kraft, S. Chatterjee, B. Fischer, and M. Koch, “Terahertz spectroscopy on polymers: A review of morphological studies,” _Journal of Molecular Structure_ , vol. 1006, no. 1, pp. 41–51, 2011, sTRUCTURAL APPLICATIONS OF TERAHERTZ SPECTROSCOPY. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S002228601100593X
  * [2] J. Neu and C. A. Schmuttenmaer, “Tutorial: An introduction to terahertz time domain spectroscopy (thz-tds),” _Journal of Applied Physics_ , vol. 124, no. 23, p. 231101, 2018. [Online]. Available: https://doi.org/10.1063/1.5047659
  * [3] J. B. Jackson, J. Bowen, G. Walker, J. Labaune, G. Mourou, M. Menu, and K. Fukunaga, “A survey of terahertz applications in cultural heritage conservation science,” _IEEE Transactions on Terahertz Science and Technology_ , vol. 1, no. 1, pp. 220–231, Sep. 2011.
  * [4] G. Filippidis, M. Massaouti, A. Selimis, E. J. Gualda, J.-M. Manceau, and S. Tzortzakis, “Nonlinear imaging and thz diagnostic tools in the service of cultural heritage,” _Applied Physics A_ , vol. 106, no. 2, pp. 257–263, 2012\.
  * [5] J.-M. Manceau, A. Nevin, C. Fotakis, and S. Tzortzakis, “Terahertz time domain spectroscopy for the analysis of cultural heritage related materials,” _Applied Physics B_ , vol. 90, no. 3, pp. 365–368, 2008.
  * [6] M. M. Nazarov, A. P. Shkurinov, E. A. Kuleshov, and V. V. Tuchin, “Terahertz time-domain spectroscopy of biological tissues,” _Quantum Electronics_ , vol. 38, no. 7, pp. 647–654, jul 2008. [Online]. Available: https://doi.org/10.1070/qe2008v038n07abeh013851
  * [7] M. Ljubenović, L. Zhuang, J. De Beenhouwer, and J. Sijbers, “Joint deblurring and denoising of thz time-domain images,” _IEEE Access_ , vol. 9, pp. 162–176, 2020.
  * [8] L.-M. Xu, W. Fan, and J. Liu, “High-Resolution Reconstruction for Terahertz Imaging,” _Applied Optics_ , vol. 53, 2014.
  * [9] W. H. Richardson, “Bayesian-Based Iterative Method of Image Restoration,” _Journal of the Optical Society of America_ , vol. 62, no. 1, pp. 55–59, 1972\. [Online]. Available: http://www.osapublishing.org/abstract.cfm?URI=josa-62-1-55
  * [10] L. B. Lucy, “An Iterative Technique for the Rectification of Observed Distributions,” _The Astronomical Journal_ , vol. 79, no. 6, pp. 745–754, 1974.
  * [11] S. Mallat, _A Wavelet Tour of Signal Processing: The Sparse Way_ , 3rd ed., 2009\.
  * [12] Y. Li, L. Li, A. Hellicar, and Y. J. Guo, “Super-resolution reconstruction of terahertz images,” in _Terahertz for Military and Security Applications VI_ , J. O. Jensen, H.-L. Cui, D. L. Woolard, and R. J. Hwu, Eds., vol. 6949, International Society for Optics and Photonics. SPIE, 2008, pp. 144 – 152. [Online]. Available: https://doi.org/10.1117/12.777814
  * [13] B. Recur, J. P. Guillet, I. Manek-Hönninger, J. C. Delagnes, W. Benharbone, P. Desbarats, J. P. Domenger, L. Canioni, and P. Mounaix, “Propagation Beam Consideration for 3D THz Computed Tomography,” _Optics Express_ , vol. 20, no. 6, pp. 5817–5829, 2012. [Online]. Available: http://www.opticsexpress.org/abstract.cfm?URI=oe-20-6-5817
  * [14] A. Dhawan, R. Rangayyan, and R. Gordon, “Image Restoration by Wiener Deconvolution in Limited-View Computed Tomography,” _Applied optics_ , vol. 24, p. 4013, 1986.
  * [15] D. C. Popescu and A. D. Hellicar, “Point Spread Function Estimation for a Terahertz Imaging System,” _EURASIP Journal on Advances in Signal Processing_ , vol. 2010, no. 1, p. 575817, 2010. [Online]. Available: https://doi.org/10.1155/2010/575817
  * [16] J. Takayanagi, H. Jinno, S. Ichino, K. Suizu, M. Yamashita, T. Ouchi, S. Kasai, H. Ohtake, H. Uchida, N. Nishizawa, and K. Kawase, “High-resolution time-of-flight terahertz tomography using a femtosecond fiber laser,” _Opt. Express_ , vol. 17, no. 9, pp. 7533–7539, Apr 2009. [Online]. Available: http://www.osapublishing.org/oe/abstract.cfm?URI=oe-17-9-7533
  * [17] H. Guerboukha, Y. Cao, K. Nallappan, and M. Skorobogatiy, “Super-resolution orthogonal deterministic imaging technique for terahertz subwavelength microscopy,” _ACS Photonics_ , vol. 7, no. 7, pp. 1866–1875, 2020. [Online]. Available: https://doi.org/10.1021/acsphotonics.0c00711
  * [18] K. Ahi and M. Anwar, “Developing terahertz imaging equation and enhancement of the resolution of terahertz images using deconvolution,” in _Terahertz Physics, Devices, and Systems X: Advanced Applications in Industry and Defense_ , M. F. Anwar, T. W. Crowe, and T. Manzur, Eds., vol. 9856, International Society for Optics and Photonics. SPIE, 2016, pp. 57 – 74. [Online]. Available: https://doi.org/10.1117/12.2228680
  * [19] T. Wong, M. Kahl, P. Haring Bolívar, and A. Kolb, “Computational image enhancement for frequency modulated continuous wave (fmcw) thz image,” _Journal of infrared, millimeter and terahertz waves_ , vol. 40, p. 775–800, 07 2019.
  * [20] L. Xu, S. Zheng, and J. Jia, “Unnatural l0 sparse representation for natural image deblurring,” in _2013 IEEE Conference on Computer Vision and Pattern Recognition_ , 2013, pp. 1107–1114.
  * [21] M. Ljubenović, S. Bazrafkan, J. D. Beenhouwer, and J. Sijbers, “Cnn-based deblurring of terahertz images,” in _Proceedings of the 15th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications (VISAPP)_ , vol. 4, 2020, pp. 323–330.
  * [22] Z. Li, Z. Cen, and X. Li, “A terahertz image super-resolution reconstruction algorithm based on the deep convolutional neural network,” in _AOPC 2017: Optical Sensing and Imaging Technology and Applications_ , Y. Jiang, H. Gong, W. Chen, and J. Li, Eds., vol. 10462, International Society for Optics and Photonics. SPIE, 2017, pp. 353 – 361. [Online]. Available: https://doi.org/10.1117/12.2283469
  * [23] Z. Long, T. Wang, C. You, Z. Yang, K. Wang, and J. Liu, “Terahertz image super-resolution based on a deep convolutional neural network,” _Applied Optics_ , vol. 58, no. 10, pp. 2731–2735, Apr 2019. [Online]. Available: http://ao.osa.org/abstract.cfm?URI=ao-58-10-2731
  * [24] T. Lei, B. Tobin, Z. Liu, S.-Y. Yang, and D.-W. Sun, “A terahertz time-domain super-resolution imaging method using a local-pixel graph neural network for biological products,” _Analytica Chimica Acta_ , vol. 1181, p. 338898, 2021\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0003267021007248
  * [25] J. M. P. Nascimento and J. M. Bioucas-Dias, “Hyperspectral Signal Subspace Estimation,” in _IEEE International Geoscience and Remote Sensing Symposium_ , July 2007, pp. 3225–3228.
  * [26] L. Zhuang and J. M. Bioucas-Dias, “Fast Hyperspectral Image Denoising and Inpainting Based on Low-Rank and Sparse Representations,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 11, no. 3, pp. 730–742, 2018.
  * [27] W. Liao, B. Goossens, J. Aelterman, H. Luong, A. Pizurica, N. Wouters, W. Saeys, and W. Philips, “Hyperspectral Image Deblurring with PCA and Total Variation,” in _2013 5th Workshop on hyperspectral image and signal processing : evolution in remote sensing (WHISPERS)_. IEEE, 2013, p. 4.
  * [28] P. F. Goldsmith, _Gaussian Beam Propagation_ , 1998, pp. 9–38.
  * [29] M. van Exter and D. R. Grischkowsky, “Characterization of an optoelectronic terahertz beam system,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 38, no. 11, pp. 1684–1691, 1990.
  * [30] L. Duvillaret, F. Garet, and J.-L. Coutaz, “Influence of Noise on the Characterization of Materials by Terahertz Time-Domain Spectroscopy,” _JOSA B_ , vol. 17, pp. 452–461, 03 2000.
  * [31] D. Krishnan and R. Fergus, “Fast Image Deconvolution using Hyper-Laplacian Priors,” in _Advances in Neural Information Processing Systems_ , 2009, pp. 1033–1041. [Online]. Available: http://papers.nips.cc/paper/3707-fast-image-deconvolution-using-hyper-laplacian-priors
  * [32] N. Wiener, _Extrapolation, Interpolation, and Smoothing of Stationary Time Series_. The MIT Press, 1964.
  * [33] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image Denoising by Sparse 3-D Transform-Domain Collaborative Filtering,” _IEEE Transactions on Image Processing_ , vol. 16, no. 8, pp. 2080–2095, 2007.