Received 29 January 2023, accepted 7 February 2023, date of publication 10 February 2023, date of current version 27 February 2023.

Digital Object Identifier 10.1109/ACCESS.2023.3243957

A Priori Information and Off-Axis Measurements

in Terahertz Tomography
KARL H. MAY 1,2 , ANDREAS KEIL1,3 , (Member, IEEE), GEORG VON FREYMANN 1,2 ,

AND FABIAN FRIEDERICH 1,2

1 Fraunhofer Institute for Industrial Mathematics ITWM, 67663 Kaiserslautern, Germany
2 Research Center OPTIMAS, Department of Physics, Technische Universität Kaiserslautern, 67663 Kaiserslautern, Germany
3 Becker Photonik GmbH, 32457 Porta Westfalica, Germany
Corresponding author: Karl H. May ([email protected])

ABSTRACT Terahertz tomography is a non-contact inspection technique to image objects from multiple
angles and reconstruct their 3D volume from intensity and time-of-flight transmission data, without the
need for radiation protection measures. Unlike X-rays, terahertz radiation is subject to strong diffraction and
refraction when propagating through dielectric materials, which often deteriorate the image reconstruction
quality. Our solution to this problem applies ray tracing, considering the beam shape and an a priori model
of the sample under investigation to predict the beam paths of the terahertz radiation. We present two
reconstruction methods based on the resulting beam path predictions yielding higher image quality. Method
1 filters out beams deviating strongly, thus removing induced artifacts and errors from the reconstruction
image. Method 2 employs off-axis measurements that acquire data along the full detection plane and in
this way detect even strongly deflected beams. Considering these beams and the information they carry
in the reconstruction enhances the image quality. Applying these methods to terahertz tomography, even
complicated structures can be imaged. We display the significant enhancements achieved with the two
methods by comparing the reconstruction results of different polymeric samples.

INDEX TERMS Terahertz radiation, computed tomography, non-destructive testing, a priori informa-
tion, off-axis measurement, reconstruction algorithms, imaging, time-of-flight, conjugate gradient least
square (CGLS).

I. INTRODUCTION Many imaging geometries create a volumetric scan of the

Terahertz radiation (0.1THz − 10THz) enables a wide range sample under test by measuring in a reflection setup with
of applications in the field of non-invasive testing and access to only one side of the sample. These approaches
imaging [1], [2], [3], [4]. Contactless measurements and, sometimes suffer from image distortion and artifacts due to
in contrast to X-Ray inspection, the absence of radiation optical effects occurring often at complicated sample struc-
protection requirements make terahertz scanning techniques tures. Moreover, potentially concealed areas behind other
an important addition to current inspection practices [5]. sample features, or generally a lack of radiation power in
Typical applications in an industrial context are non- certain geometric configurations can lead to errors and arti-
destructive layer thickness determination and defect detection facts when imaging with one-sided approaches. Especially
using volumetric imaging [6], [7], [8], [9], [10]. Materials for additively manufactured [13] or extruded objects [14],
of interest range from ceramics, glass-fiber reinforced com- such as window profiles, performing a full terahertz transmis-
posites [2] and plastics [6] over paints and coatings [9] to sion tomography on the sample from many different angles
paper [11] and wood [12]. can be beneficial. For an overview of tomographic techniques
applying terahertz radiation, refer to [15] and [16].
In classical X-Ray transmission tomography, projec-
The associate editor coordinating the review of this manuscript and tions of an object under test are acquired from different
approving it for publication was Yi Zhang . angles and a volume reconstruction is performed which is

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/
VOLUME 11, 2023 18311
 K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

FIGURE 1. Ray tracing simulation of terahertz beam paths (a) assuming infinitely small ray diameters and (b) considering the Gaussian beam shape (one
ray iconsists of many thin rays). The density of the plotted beams is reduced in comparison to the simulations to improve visibility.

predicated on the Fourier-slice theorem [17]. A practical and

fast algorithm for this reconstruction is the filtered back-
projection (FBP), assuming that the radiation travels through
the sample in straight lines. For tomography using X-rays
applied in a medical or NDT context, this is usually a good
assumption. In contrast to X-ray radiation, terahertz radiation
is subject to significant diffraction and refraction, due to its
larger wavelength. Additionally, the wavelength of terahertz
radiation (in the mm to sub-mm range) is of the same order
of magnitude as the typical features of the samples, often
leading to a deflection from its straight path (see Fig. 1).
In this work we introduce a measurement setup, which detects
radiation passing through the imaging scene on a straight line
and also acquires data off-axis to take beams into account that
are deflected by the sample. The information carried by these
beams would otherwise be lost. FIGURE 2. Measurement setup with rotational stage (purple) and a linear
To utilize the full data set, we resort to a versatile, flex- stage (orange) to turn and move the sample (yellow) in z-direction. Two
independent translational stages for receiver (blue) and transmitter
ible problem formulation, and an iterative reconstruction (green) allow the acquisition of deflected and non-deflected radiation.
algorithm to solve it. We apply the conjugate gradient least The radiation is collimated and focused by the lens pairs in front of the
transmitter/receiver units. The dashed red line indicates a possible beam
square (CGLS) algorithm, a robust and efficient method to path.
solve tomography problems iteratively. It has been applied
to terahertz tomography for the first time in [18]. It solves
large linear systems quickly and offers the possibility to apply
constraints to the reconstruction space to improve the recon-
struction. The CGLS algorithm solves the inverse problem of on enhancing the imaging capabilities of different tomo-
the simple matrix multiplication graphic reconstruction techniques based on a priori infor-
mation can be found for microwave tomography in [19],
Ā⃗x = b⃗ (1)
[20], [21], [22], and [23], for X-Ray tomography in [24], and
where the vector b⃗ contains the measurement data and x⃗ in the context of terahertz tomography in [25], [26], and [27].
represents the unknown reconstruction. The matrix Ā con- Ding et al. [19] distinguish between two categories of a priori
tains information about how the radiation travels through the information: The first category represents information about
imaging scene. Since Ā is derived from a physical model, this the target or sample, such as internal or external boundaries
formulation allows incorporating a priori information about or regions and the respective refractive indices [20], [21],
the measurement system and the sample into the reconstruc- [22], [24]. The second category considers the measurement
tion process and thus, improves reconstruction quality and system, i.e. the number and position of antennas or infor-
resolution. Modeling the sample a priori is a common tool mation about the background [23]. In this work, we sug-
in the field of NDT, to assure the quality of objects made gest a method for the integration of a priori information
out of a well-known material by either extraction or addi- from both categories into the reconstruction formula-
tive manufacturing based on a geometric model. Prior works tion. In section IV-C. the data stemming from the off-axis

18312 VOLUME 11, 2023

K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

measurements are used to improve the reconstruction images Lambert-Beer’s law:

even further.  
I X
τR = = exp − αi |x⃗i | (4)
II. MEASUREMENT SETUP I0 R i∈IR
The measurement system is a single-pixel setup with one
′′
transmitter and one receiver as shown in Fig. 2. The transmit- where αi = 4ni πf /c0 represents the absorption coefficient
ter, consisting of a voltage-controlled oscillator and several and f is the frequency of the radiation. In a logarithmic
active frequency multipliers, continuously emits a frequency- representation of the intensity, measured at a given rotation
modulated continuous-wave (FMCW) signal in the range angle ϑ at a distance d from the center of rotation, (4) can be
from 230 GHz to 320GHz from a Pickett-Potter horn antenna. simplified to
The beam is collimated by a f = 50 mm lens and focused
⃗ i · Iiτ
X X
to the center of the setup by a f = 200 mm lens. The τln (ϑ, d) = ln (τR ) = − αi |x⃗i | = A (5)
receiver has the same lens configuration as the transmitter, i∈IR i∈I
focusing the radiation into the detector antenna. The lenses
have an aperture of 2′′ , which was found to be a good trade-off The last expression represents a product of the image
between resolution and illumination for the samples we used. vector I⃗τ consisting of the entries αi for every pixel of the
To image larger samples, it is possible to use lenses with a imaging scene with the Matrix Ā comprised of the column
larger aperture, broadening the beam. The receiver consists vectors A⃗i (ϑ, d), defining which pixels contribute to which
of a subharmonic mixer, mixing the local oscillator signal ray to which extent. The αi represent the material absorption
with the transmitted radiation, acquiring the signal amplitude coefficient αmat at the pixel position. This model does not take
and phase coherently. This way, the signal attenuation and the losses due to reflection at surfaces into consideration.
time delay caused by traversing the sample are measured. Similarly to (5), (3) can be expressed as:
The mathematical model of the signal after passage through X
⃗ i · IiT
TR (ϑ, d) = A (6)
the sample is shown in the next section. The transmitter and
i∈I
receiver are each positioned on separate translation stages
allowing an independent movement on the Y-axis perpendic- Here, the same matrix Ā is utilized to calculate the distribution
ular to the signal propagation direction. This way, not only of n′ (x) in the imaging area by solving Ā · ⃗I T = T⃗ (ϑ, d) for
radiation traveling straight through the imaging scene but also ⃗I T and Ā · ⃗I τ = τ⃗(ϑ, d) for ⃗I τ .
deflected beams can be detected performing off-axis mea- This model assumes the medium to be non-dispersive, i.e.
surements. Furthermore, the setup features a rotational stage n′ (x) is independent of frequency. T⃗ (ϑ, d) and τ⃗(ϑ, d) rep-
to turn the sample under test and another translation stage resent the vectors containing the measurement data, namely
to lift or lower the sample. By combining cross-sectional the ToF and intensity sinograms, respectively. These formu-
images (B-scans) acquired at different heights, a full 3D las were derived and discussed in detail in [18], [28], [29],
reconstruction of a sample can be created. and [30].
The forward problem of terahertz tomography can thus be
III. A PRIORI INFORMATION described by a simple matrix multiplication
A. MATHEMATICAL MODEL OF TRANSMISSION
TOMOGRAPHY Ā⃗x = b⃗ (7)
Tomographic imaging of dielectric objects reconstructs the
with the vector b⃗ denoting the sinogram of the measurement
complex refractive index
data, the unknown vector x⃗ representing the image to be
reconstructed and the matrix Ā mapping the image pixels to
n(xi ) = n′ (xi ) + in′′ (xi ) (2)
the corresponding sinogram entries. The latter contains the
at each pixel xi ∈ I in the discretized image area I . trajectories of the rays describing how the terahertz rays travel
The time-of-flight (ToF) delay TR of a terahertz signal through the measurement area. To calculate these trajectories
R traversing through an imaging scene is proportional to and thereby model the tomographic process as accurately
(n′i (xi ) − nair ): as possible, a priori information about the measurement
system and sample under test can be incorporated into the
1 X design of Ā.
n′i − nair |x⃗i |


TR = (3)
c0 i∈IR
B. OPTICAL EFFECTS OF THE SAMPLE BOUNDARIES
where c0 is the vacuum speed of light, and x⃗i represents The first category of a priori information we consider takes
the path through the pixels xi in the subset IR ⊂ I the the geometry of the sample under test and its real refractive
Ray R traverses. In the following, the refractive index of index n′ into consideration. The information about the surface
air is approximated with nair ≈ 1. For a given beam the boundaries of the sample under test is used for ray trac-
relative intensity loss τR due to absorption is according to ing simulations including optical effects like refraction and

VOLUME 11, 2023 18313

 K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

Algorithm 1 Ray Tracing Algorithm Returning the Beam

Paths x⃗li,θ
1: for θ = 1 . . . Nθ
2: for i = 1 . . . N
3: determine d⃗li,θ
4: for l = 0 . . . lmax
5: t=0
6: determine thit > 0 values for Bscene and BSUT
7: if x⃗li,θ (min(thit )) = BSUT
i,θ
8: determine x⃗l+1,0 = x⃗li,θ (thit )
i,θ
9: determine d⃗l+1 from equation (8)
10: else if x⃗li,θ (min(thit )) = Bscene (−xbound )
i,θ
11: discard x⃗0...l ; break
FIGURE 3. Intensity distribution of the beam due to the lens system. The
X-axis represents the axial direction, the Y-axis is the radial direction.
12: else if x⃗li,θ (min(thit )) = Bscene (±ybound ) and d⃗li · ⃗ex < 0
i,θ
(Compare Fig. 1 (b)). The blue lines indicate the theoretical boundaries of 13: discard x⃗0...l ; break
a Gaussian beam (defined by an intensity decay to 1/e2 of the value on 14: else
the beam axis) for a wavelength of λ = 1 mm and a focus beam width of 15: break
w0 = 5.2 mm.
16: end
17: end
18: end
reflection (depicted in Fig. 1). If a simulated ray path inter- 19: return all x⃗li,θ
sects with the sample, the diversion of the beam trajectory at
the sample boundary is calculated according to Snell’s law:
n′1 sin (ε1 ) = n′2 sin (ε2 ) (8)
determines the radius of the beam at any distance z from the
where ε1 and ε2 are the incident angle and the angle of πw2
focus. w0 is the center beam waist and zR = λ 0 is the
the refracted beam, respectively, both measured towards the Rayleigh length.
surface normal. n′1 and n′2 are the real refractive indices of The two blue lines in Fig. 3 indicate the 1/e2 intensity
the two materials, that form the boundary. For an incident expected from the Gaussian beam model with a center wave-
angle of length of λ = 1 mm. The focus beam width of 5.2 mm
 ′
n was found by fitting the theoretical curve to the data. This
ε ≥ εc = arcsin 2′ (9) corresponds to a Rayleigh range of 85mm.
n1
The measured intensity agrees well with the theoretical
the critical angle, total internal reflection causes the beam to
distribution, indicating that the Gaussian beam shape is a
be reflected at the surface. Whether it is refracted or reflected,
good approximation for the given beam intensity distribu-
the new propagation direction of the beam is calculated.
tion. The beam shape enters the construction of the matrix
After simulating the necessary beam paths for all transmitter
Ā by including multiple infinitely small elementary beams
positions, the resulting trajectories are incorporated into the
that lie within the propagation of the full beam as indicated
design of the matrix Ā (see section IV-A).
in Fig. 1 (b).
C. BEAM SHAPE
The shape of the beam influences the measurement data and IV. RECONSTRUCTION TECHNIQUES
can be considered when reconstructing an image [25]. The A. RAY TRACING AND CONSTRUCTION OF Ā
lens system of the measurement setup (see Fig. 2) focuses By utilizing a priori information, we improve the matrix Ā
the terahertz radiation into the center of the setup. To measure in comparison to the standard design, which only assumes
the beam shape, the receiver is placed onto an X/Y-translation straight and infinitely thin rays. To do so, just as for conven-
stage, allowing it to move in the axial direction and perpen- tional algorithms the matrix entry aij represents the length of
dicular to it. The beam pattern is shown in Fig. 3. The beam the path the i-th ray travels in the j-th pixel Pj of the imaging
profile can be modeled as a Gaussian beam, whose intensity area. Since the beam paths are influenced by the material
distribution depends on the distance from the beam axis r and parameters, the rays do not propagate straight through the
the distance from the focus point z as follows: sample and thus, the full matrix cannot be calculated in
one step. Due to the finite diameter of the Gaussian beam,
w0 2 2r 2
   
I (r, z) = I0 exp − 2 (10) we build up one wide beam out of several infinitely thin
w (z) w (z) elementary beams. The ray tracing algorithm to determine
where I0 is the maximum intensity at the focal point, and the trajectories of these elementary beams is described in the
s
 2 follwing.
z The pseudo-code for the ray tracing can be found in
w (z) = w0 1 + (11)
zR Algorithm 1. Every part l of an elementary ray i is defined

18314 VOLUME 11, 2023

K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

FIGURE 4. Method 1 (a priori): Image reconstruction considering a priori information, filtering out deflected
beam paths.

i and its direction d

by its start position x⃗l,0 ⃗ i: In order to account for the Gaussian beam shape of the rays
l
and finally build the matrix Ā, we follow Algorithm 2. For
x⃗li (t) = x⃗l,0
i
+ t · d⃗li (12) every sender and receiver position s = 1 . . . S, we consider
the elementary rays, that lie within the Gaussian beam (see
The propagation of the ray is represented by an increase of Fig. 1). The definitions of the maximum beam width wmax and
the parameter t. By finding the intersections of (12) with the the angle extrema θmin and θmax can be found in the appendix.
given boundaries of the sample under test (SUT) BSUT and the Defining the Gaussian beam as a sum of multiple ele-
imaging scene Bscene and calculating the respective t values, mentary beams does not describe the beam perfectly, but it
one can find the first boundary the ray intersects with – repre- is an approximation that serves our purposes. In particular,
sented by the lowest positive thit . If this boundary corresponds the linear relations between the sinogram value τ⃗/T⃗ and the
to the imaging scene bounding box, the ray propagation is distance x⃗ traveled in the imaging scene in (5) and (6) are
completed and we continue with the next ray. Otherwise, if a preserved by this approach.
sample surface boundary is reached, we calculate the new
i
start position x⃗l+1,0 = x⃗li (thit ) and the new d⃗l+1
i accord- B. METHOD 1:A PRIORI
ing to (8). This procedure is repeated until l reaches lmax . The inclusion of the a priori information presented above
To avoid infinite loops, we limit the maximum amount of leads to an improved model of the forward problem of
reflection/refraction incidents to lmax = 5, after that the ray tomography represented by (7). This approach, here called
is discarded. This was found as a suitable choice considering method 1, shares some similarities with [27]: When con-
the strong absorption such a ray would experience due to the structing the matrix Ā from the trajectories simulated by our
long distance of travel. A ray is also discarded if it coincides Ray tracing algorithm (Algorithm 2), we include the deflected
with the back boundary of the imaging window (xbound = beam paths by connecting the measurement value in the
−60 mm in Fig. 1), or if it is heading towards it when leaving sinogram with the pixels the deflected beam passes. Hereby,
the imaging scene bounding box. In Fig. 1(a) we see two in this first method, we include only rays, whose trajectories
strongly diverted blue beams which will be discarded and end at a position, that is covered by the aperture of the receiver
three yellow beams which will be kept since their trajectories placed at the same position on the Y-axis as the transmitter.
can be extrapolated to finally reach the detector plane. This Therefore, only rays with a maximum deflection of 25mm
procedure is repeated for all rays i= 1 . . .N and for all angles from the receiver and sender position y0 are considered, while
θ = 1 . . . Nθ . all other rays are ignored. In contrast to [27] we only take the

VOLUME 11, 2023 18315

 K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

FIGURE 5. Method 2 (a priori + off-axis): Image reconstruction considering a priori information using deflected beams.

outer surface boundaries of the samples into account, when in the reconstruction, especially when reducing the density
performing the ray tracing simulations. This way we can also of the measured angles to increase the measurement speed.
reconstruct defects and features of the sample, whose exis- Therefore, we introduce the second reconstruction method,
tence was a priori not known. This reconstruction procedure, also considering these strongly deflected beams.
schematically depicted in Fig. 4, allows the reconstruction
of samples from measurement data resulting from a parallel C. METHOD 2:A PRIORI + OFF-AXIS
movement of sender and receiver, placing them in front of The second method is schematically depicted in Fig. 5.
each other at any time. Alternatively, the sinogram b⃗ can be As indicated in section II, we designed the measurement setup
extracted from a full off-axis measurement by creating pro- allowing an independent movement of the transceiver and
jections along the diagonal ytrans = yreic , where the position receiver unit, rendering the detection of strongly deflected
of the transmitter equals the position of the receiver on the radiation possible. For every transmitter position, we move
opposite side of the beam axis. The diagonal is indicated in the receiver off-axis along its full range of motion spanning
Fig. 4 by the red line in the lower left full data block. Follow- 400 mm. For a measurement resolution of 0.5 mm for the
ing this procedure for every angle θ yields the sinogram b. ⃗ receiver unit and 1 mm for the transmitter range of 120 mm,
Since strongly deflected beams do not pass the receiver aper- this leads to a full data block of 800 × 120 values for every
ture, the information they carry is lost. This can induce errors angle for a single cross-sectional image. It is schematically

18316 VOLUME 11, 2023

K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

FIGURE 6. Comparison of sinograms; (a), (b) simulated without a priori information; (c), (d) including optical effects on the sample boundaries;
(e), (f) full a priori information (Gaussian beam shape and opt. effects); (g), (h) measured data of Sample 1(a); intensity and time-of-flight (ToF).

Algorithm 2 Calculating the Reconstruction Matrices Ā and both sinograms and analogously to method 1 apply the
and ĀASM (Additional Scatter Matrix Including a Priori CGLS algorithm to solve the inverse problem and reconstruct
Information). See Also Appendix. the image ⃗I :
1: for s = 1 . . . S:
b⃗
   
2: for i = 1 . . . N : Ā ⃗
·I = ⃗ (13)
3: get trajectory x⃗i,θ from Algorithm 1 ĀASM bASM
i,θ
4: calculate distance d = x⃗0,0 − y0 (s) ymax and bASM contain the information of the strongly
5: if d ≤ w(−xbound ) deflected beams about the object’s internal structure. How-
6: calculate θmin (d) and θmax (d) from A.(15) ever, they might be prone to errors due to possible deviations
7: for θ = θmin . .
. θmax of the model from the real
measurement process. We found
8: if end x⃗i,θ − y0 (s) ≤ w(xbound ) that combining Ā, ĀASM and (b, bASM ) yields better results
9: aij = lenght(⃗x i,θ ∩ Pj ) than reconstructing from either of the pairs solely.
10: end Even though (13) is twice the size of the original problem
if end x⃗i,θ − ymax (s) ≤ w (xbound )


11: (e.g. (7)), its inversion can still be solved within seconds,
12: aASM
ij = lenght(⃗x i,θ ∩ Pj ) thanks to the excellent performance of the CGLS algorithm,
13: end which is explained in the next section.
14: end
15: end D. CONJUGATE GRADIENT LEAST SQUARE (CGLS)
16: end ALGORITHM
17: end In literature, there exists a large variety of methods to tackle
18: return Ā, ĀASM inverse problems, such as (7) and (13). For an overview,
refer to [28] and [33]. While some are based on statisti-
cal approaches like the Bayesian methods, others, the so-
depicted in the lower-left corner of Fig. 5 and contains called ‘‘row-action’’ methods like algebraic reconstruction
the attenuation and time-of-flight (ToF) information of all, technique (ART) or Split-Bregman method [31], consider the
including strongly deflected, beams, provided their incident matrix one row at a time. The latter option performs typically
angle is not too steep to pass the detector optics. Based on slow on large matrices. In this work we utilize the conju-
the ray tracing simulations including the a priori information gate gradient least square algorithm (CGLS), which iterates
about the sample under test and the beam shape, we simulate towards an optimal solution following conjugate gradients.
these data blocks for every measurement angle. From the It performs fast and efficiently, even on the typically large
simulated data, we determine the receiver position ymax (s) and non-square system matrices of the tomography problem.
in the detector plane, where we expect most radiation to The CGLS algorithm has been covered extensively in
be detected. The value intensity and ToF value of ymax are [32], [33], and [34] and in the context of terahertz tomography
denoted in the scatter sinogram b⃗ASM , where ASM stands in [18]. Here, we will give a brief outline of the main idea.
for additional scatter matrix. The beam pixels associated The CGLS algorithm is an iterative algorithm for approx-
with all the beams reaching the aperture of the receiver at imating the solution xi of an inverse problem such as (7),
the position ymax , are denoted in the scatter matrix ĀASM , for a known positive-definite (not necessarily square) matrix
following Algorithm 2. Finally, we combine both matrices Ā and an also known measurement vector b. ⃗ It does so by

VOLUME 11, 2023 18317

 K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

FIGURE 7. Off-axis measurement and corresponding simulation of Sample 1(a) for θ = 50◦ .

updating an initially given solution x0 going in steps, which exists in a version (a) without defects and a version (b) with
are all ĀT Ā-conjugate to each other. This way, CGLS iterates defects. In the case of Sample 1(b) (Fig. 8) and 2(b) (Fig. 11),
towards the optimal solution of the problem, provided that the the defects consist of three concentric holes with diameters
measurement vector is noiseless. When noise is present, the of 3, 4, and 5 mm placed with center distances of 8.5 mm
CGLS algorithm becomes semi-convergent, i.e. after finding and 7.5 mm, respectively. The defect of Sample 3(b) (Fig. 13)
the optimal solution, it continuous iterating often deterio- has the shape of a pill with a width of 10 mm and a length
rating the result. Therefore, one has to apply regularization of 50 mm.
measures. This can be done by limiting the solution space (for Projections of the samples were acquired in the angle range
example restriction to R+ ) or stopping the iteration process [0◦ , 358◦ ] in steps of 2◦ . When taking optical effects into
before deterioration starts. The right number of iterations account, it is important to cover the whole circumference and
to stop at can be found by applying the L-curve criterion not only the half-space of 180◦ , since the resulting sinograms
[18], [35]. The reconstructions obtained by applying the do not necessarily have rotational symmetry. For every angle
CGLS algorithm in combination with the above-mentioned position, the transmitter unit is moved along the sample in
reconstruction methods are displayed in the following steps of 1 mm (for example over a width of 120 mm, depend-
section. ing on the size of the sample). The receiver unit continuously
acquires data moving along its full range of 40cm with a
V. RESULTS AND DISCUSSION step size of 0.5 mm. This procedure results in a data set of
A. EXPERIMENT DESCRIPTION 800 × 120 data points for every angle, the relevant subset of
To test the capabilities of the different reconstruction meth- which is shown exemplarily for Sample 1(a) in Fig. 7.
ods, we designed three different samples displayed in Fig. 8, B. SIMULATION VALIDATION
Fig. 11, and Fig. 13. Samples 1 and 2 are cuboids of Before discussing the reconstruction results, we compare the
30 × 40 × 60 mm3 made from Polyvinyl chloride (PVC) and simulated data with the measurements, in order to verify
Polymethylmethacrylate (PMMA), respectively. Sample 3 that the inclusion of the a priori information improved the
consists of a 25 × 75 × 75 mm3 block made from Polyethy- model of the measurement process. Fig. 6(a)-(f) show the
lene (PE). The optical parameters of the three materials at sinograms resulting from ray tracing of Sample 1(a) (see
300 GHz are indicated in Table 1. section IV-A above). Fig. 6(g),(h) show the corresponding
measured sinograms. The upper row of Fig. 6 displays the
TABLE 1. Optical properties of the used materials [36].
intensity sinograms while the data in the lower one repre-
sents the ToF. The simulated sinograms 6(a) and 6(b) assume
straight ray trajectories and infinitely thin beams – the stan-
dard configuration for (X-ray) tomography. Apart from a
factor, the two sinograms are identical, since all rays travel
on straight trajectories interacting proportionally to the path
Each sample was constructed twice. Holes were drilled into length they travel through the sample. This is not the case for
one specimen serving as artificial defects so that every sample Fig. 6(c) and (d) displaying the simulated sinograms when

18318 VOLUME 11, 2023

K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

FIGURE 8. Reconstructions of the absorption coefficient α (a,b,e,f,I,j) and the real refractive index n (c,d,g,h,k,l) of Sample 1(a) and 1(b) made out of PVC.
X- and Y-axis in mm. The reconstructions where performed without a priori information (a)-(d), with method 1 considering beam shape and optical
effects (e)-(h), with method 2 additionally considering scattered radiation from off-axis measurements (i)-(l). Iteration numbers of the CGLS
reconstruction algorithm are given in the lower right corner.

refraction and reflection are considered. The simulations take

the surface boundaries and the refractive index of the sample
as a priori information, assuming thin rays. Following the
simulation shown in Fig. 1(a), for specific angles ‘blind areas’
appear, where no radiation reaches the detector aperture.
They are visible in all following sinograms as dark areas
around the angles θi∈N = 45◦ + i · 90◦ , where the radiation
is mainly deflected by the sample. In Fig. 6(c) and (d) said
areas are too large and too sharp in comparison with Fig. 6(g) FIGURE 9. Photography of Sample 1(a) and 1(b) made from
Polyvinylchlorid (PVC): 30 mm × 40 mm × 60 mm (b) hole diameters:
and (h) since in the experiment the beam has a Gaussian 3, 4, and5 mm.
shape and therefore a finite radius (see section III-C). This
fact is considered in Fig. 6(e) and (f), where both the optical
effects as well as the beam shape are regarded as a priori Similarly, we compare the results from the off-axis
information. As a result, the ‘blind areas’ become smaller and measurements and the corresponding simulations. Fig. 7.
less sharp, since due to the larger beam radius some radiation shows the data resulting from an off-axis measurement of
can still circumvent sharp corners of the sample and reach the Sample 1(a) at an angle of 50◦ . The transmitter was moved
detector aperture. Overall, by considering both sets of a priori from −60 mm to 60 mm and the receiver along its full
information the sinograms Fig. 6(e) and (f) are in very good range of motion from −200 mm to 200 mm. To increase
agreement with the measured data as seen in Fig. 6(g) and (h). visibility, the y-axis was adjusted here to display only the
This promises good reconstruction results because it shows relevant data range. The two graphs on the right side dis-
that the most prominent effects influencing the measurement play the corresponding simulated data from the ray tracing
process have been identified and considered. procedure.

VOLUME 11, 2023 18319

 K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

FIGURE 10. Reconstructions of the absorption coefficient α (a,b,e,f,I,j) and the real refractive index n (c,d,g,h,k,l) of Sample 2(a) and 2(b) made out
of PMMA. X- and Y-axis in mm. The reconstructions where performed without a priori information (a)-(d), with method 1 considering beam shape
and optical effects (e)-(h), with method 2 additionally considering scattered radiation from off-axis measurements (i)-(l). Iteration numbers CGLS
reconstruction algorithm are given in the lower right corner.

The bright areas in the x = y diagonal correspond to

the radiation traveling from the transmitter directly into the
receiver unit. The ToF, or more specifically the time delay
induced by the target, is 0 since the beams do not interfere
with the sample under test. The intensity reaches values
around 0 dB here, indicating maximum transmittance. In the
central part, however, where the sample is placed, the radi-
ation is split into two parts, very much in accordance with
FIGURE 11. Photography of Sample 2(a) and 2(b) made from
Fig. 1(b). Here the actual interaction between the radiation Polymethylmethacrylat (PMMA): 30 mm × 40 mm × 60 mm (b) hole
and the sample becomes evident – in the form of a reduction diameters: 3, 4, and5 mm.
of intensity and an increase of the ToF, respectively. Around
the edges of the sample there are dark areas indicating that,
due to scattering effects, less radiation enters the receiver unit. displayed in Fig. 9, Fig. 10, and Fig. 12. The respective
Fig. 7 shows a good agreement between measurement and number of iterations of the CGLS algorithm is given in the
simulation as well, justifying the model for the forward prob- lower right corner of every reconstruction.
lem and the assumptions and simplifications we make. This
is crucial for the success of reconstruction method 2 since it 1) SAMPLE 1
uses the off-axis simulations to find the position of maximum Fig. 9 depicts the reconstructions of Samples 1(a) and (b)
radiation intensity ymax . made from PVC. The first row Fig. 9(a)-(d) shows the recon-
structions with the CGLS algorithm utilizing the standard
C. IMAGE RECONSTRUCTIONS matrix without regarding any a priori information. The recon-
The reconstructions of the six samples with the two recon- structions from the intensity data as well as from the ToF
struction methods described in sections IV-B. and IV-C. are data indicate how the refraction and scattering effects are

18320 VOLUME 11, 2023

K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

FIGURE 12. Reconstructions of the absorption coefficient α (a,b,e,f,I,j) and the real refractive index n (c,d,g,h,k,l) of Sample 3(a) and 3(b) made out
of PE. X- and Y-axis in mm. The reconstructions where performed without a priori information (a)-(d), with method 1 considering beam shape and
optical effects (e)-(h), with method 2 additionally including scattered radiation from off-axis measurements (i)-(l). Iteration numbers CGLS
reconstruction algorithm are given in the lower right corner.

erroneously interpreted as areas with strong absorption or a

high refractive index n, respectively. For the reconstructions
from the ToF data, n reaches a value of up to 3 leading to
very dark edges (for the expected values refer to Table 1).
Additionally, the sharp corners of the samples are not recon-
structed, so that the shape of the reconstructed object appears
rounded. The defects in Sample 1(b) are visible, but, because
their boundaries strongly scatter the radiation, they appear to FIGURE 13. Photography of Sample 3 (a) and 3(b) made from Polyethylen
have a higher absorption value or refractive index than the (PE): 20 mm × 75 mm × 75 mm (b) drilled-out ’’defect’’ ∼ 50 mm × 10 mm.

surrounding material.
We compensate for these effects by applying reconstruc-
tion method 1 introducing a priori information to the matrix. For the n reconstruction from the TOF of Sample 1(b)
The visibility of the object’s shape in the reconstruction is (Fig. 9(h)), however, the three defects are perfectly visible
much better, as depicted in Fig. 9(e)-(h). The corners are dis- as holes in an otherwise uniform block. Although it would
tinctively visible in all four reconstructions and the edges of be difficult to define the exact outlines of the defects from
the sample appear less pronounced. In the reconstruction of α the reconstruction and thus compare them to their real size,
from the intensity data of Sample 1(b) (Fig. 9(f)), the defects it is obvious that the defects’ dimensions appear exaggerated.
are still visible as strong absorbers. This is due to the fact, that A possible reason for this could be the fact that the rays
the outlines of the defects or even their existence were not part that interfere with the defects are the ones that are refracted
of the a priori information, so they were not considered in the when traveling through the sample. Therefore, in many cases,
matrix. This is advantageous for NDT applications because they do not reach the detector aperture and their information
ultimately the goal of NDT is to find defects, whose existence is not included in the reconstruction. Consequently, it is
is unknown a priori. necessary to consider the deflected radiation as well. We do

VOLUME 11, 2023 18321

 K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

this by applying reconstruction method 2 (Fig. 9(i)-(l)). For

the sample without defects (Fig. 9(i) and (k)), the uniformity
of the reconstruction is improved in comparison to method 1,
since the deflected radiation traveling diagonally through the
samples is considered in the reconstruction. For Sample 1(b),
carrying defects, the reconstruction improves as well. The
defects become sharper and more pronounced than before,
allowing a clear distinction between a sample with and a
sample without defects, especially when one looks at the
reconstruction of the refractive index n.

2) SAMPLE 2
FIGURE 14. 3D reconstruction of Sample 1(b).
The reconstructed images of Samples 2(a) and (b) are dis-
played in Fig. 10. For the classic reconstruction process with-
out a priori information, the deflection of the rays resulting
from the optical effects at the sample surface boundaries in Fig. 12. Neither the consideration of a priori information
leads to even worse reconstructions than for Sample 1. The nor the utilization of off-axis data could cope with this effect,
outer surface boundaries are much too pronounced and the completely covering up the relevant information about the
corners disappear so that the blocks seem to be round. Again, sample interior. Solely the outer shape of the sample is visible
scattering effects make the defects appear as dark areas in in Fig. 12(e) and (i). The TOF measurement data on the other
the image, however, for the reconstruction from the ToF hand, since it is not sensitive to Fresnel losses, leads to very
data, the low reconstruction quality does not allow distin- good reconstruction results.
guishing the holes. The reconstructions from the plain CGLS algorithm show
Using a priori information improves the reconstruction the same problems as before, namely round corners and
significantly, similar to Sample 1. Especially in the recon- overly pronounced edges. The ‘‘defect’’ in Fig. 12(d) is vis-
structions of Sample 2(a) (Fig. 10(e) and 10(g)) the values ible, but its shape is distorted, as it seems to be a rectan-
of α and n are varying less throughout the uniform parts gle rather than a pill with round corners. This changes for
of the object. Two of the three defects of Sample 2(b) are the reconstructions with method 1 (Fig. 12(g) and (h)). The
visible in Fig. 10(f) as dark areas. Scattering is likely to be outer corners of the objects are much sharper and the correct
the reason for the cross-shaped artifacts around the defects. value of n ≈ nmat = 1.52 is determined. In addition, the
In the reconstruction of n (Fig. 10(h)), the defects are vis- shape of the defect is imaged correctly. On either side of the
ible but very blurry. Their geometric dimensions are too defect, there are areas erroneously reconstructed with a lower
large and their boundaries are not clearly distinguishable. absorption coefficient. These are artifacts, which can be sup-
This improves when the scattered radiation is considered. pressed by considering the off-axis data in the reconstruction.
In Fig. 10(l) the defects appear smaller and sharper even Applying method 2 (Fig. 12(k) and (l)), the solid parts of the
though their outlines are still diffuse. In Fig. 10(j), the recon- samples are reconstructed more precisely and appear more
struction of α, the cross-shaped artifacts are reduced, so that uniform, because radiation traveling diagonally in the sample
the reconstruction is closer to being uniform in the parts of the is considered. The defect is reconstructed by this method with
sample without defects. The same holds true for Fig. 10(i) and its correct shape, too.
Fig. 10(k), the reconstructions of Sample 2(a) considering As already stated in [18], when inspecting PE with tera-
a priori information and the off-axis measured data. Overall, hertz tomography one should mainly rely on the ToF data
by combining the results of intensity and ToF reconstructions, analysis for sample quality evaluation, rather than on the
the unknown defects can be detected and a better sample intensity data. The reconstructions of n are more reliable in
representation in the reconstructed image can be achieved. the case of Sample 3. The good results were improved further
However, the defect size and arrangement bring the mea- by applying methods 1 and 2 displaying the defect within the
surement setup to its resolution limit, so that one of the sample very accurately.
defects is only barely visible regardless of the reconstruction
method. VI. CONCLUSION AND OUTLOOK
We have shown that a priori information can have a strong
3) SAMPLE 3 positive impact on the reconstruction quality of terahertz
Regarding its very low absorption coefficient αmat Sample 3 tomography. By including information about the beam shape
differs drastically from the other two samples. As a result, the and optical effects occurring at the sample boundaries,
intensity decrease of the radiation when traversing the sample we were able to get largely accurate reconstructions of the
is much lower. The Fresnel losses occurring at the sample samples. We have shown that defects down to a size of
surface boundary therefore dominate the intensity measure- 2 mm can be resolved. Even defects whose existence was
ment, which leads to the low reconstruction quality visible not included in the set of a priori information were imaged

18322 VOLUME 11, 2023

K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

more accurately and reliably than it is possible without the

application of this new technique. We are able to improve the
good results even further by including the analysis of off-axis
measurements in the reconstruction. This way deflected radi-
ation is utilized for the reconstruction process, thus sharpen-
ing the images of samples and defects. This is also illustrated
in Fig. 14, which shows a 3D reconstruction of Sample 1(b)
created by combining multiple B-scans (slice images) at dif-
ferent heights of the sample. Since the shapes of the samples
under investigation are similar, the differences between the
images resulting from the respective reconstruction schemes
of samples 1 to 3 mainly result from the different optical
properties. This is most evident in the reconstructions from
the intensity data, since the samples’ absorption coefficents
differ substantially. On the other hand, the reconstructions of
FIGURE 15. Definition of the angles to parametrize the Gaussian beam
the real refractive index n (all in the range between 1.5 and from several elementary beams.
1.7) show similar performance for all three samples.
One drawback of terahertz tomography with off-axis mea-
surements is the long acquisition time. Applying this method
highly parallelizable calculations. These procedures have to
the transmitter and the receiver are moving independently
be performed only once for every sample structure or outline
from each other along their full range of motion. The acqui-
and do not have to be repeated for every measurement e.g.
sition time of one of the cross-section reconstructions seen
in an industrial scenario. They can be stored and utilized for
above is around six hours. This time could be decreased
the reconstruction of every sample with the respective outline
significantly since only a small subset of the acquired data
and material. The calculation of the final reconstruction by
is actually used for the reconstruction, even when applying
CGLS takes about 1s for method 1 and 1.5s for method 2,
method 2. Building the sinograms b⃗ and b⃗ASM in (13) requires
depending on the number of iterations. The stated values hold
only the data acquired by the receiver at the opposite side
for around 10 iterations. In comparison with the measurement
of the transmitter position and on the position ymax , i.e. the
time, the evaluation time is therefore negligible.
receiver position at which the ray tracing simulation predicts
the most radiation to be received. Since ymax is determined
APPENDIX
a priori, only the data from these two positions has to be
In IV-A we introduce the concept of simulating a wide
acquired by the receiver.
Gaussian beam as a combination of multiple infinitesimal
Acquiring the full range of receiver data, however, gives
thin elementary beams. We do so by defining a maximum
rise to the opportunity of utilizing even more off-axis data
width at the start position wmax = w(−xbound ), which can
from different positions and thus potentially improving the
be calculated following (11), where 2·x bound is the width of
reconstruction quality further. For future work, we will
the imaging scene. Elementary rays x⃗i,θ are considered for
speed up the acquisition without reducing the field of view. i,θ
We intend to apply detector arrays [37] covering the full range the reconstruction, if the distance d = x⃗0,0 − y0 (s) ≤
of the receiver movement simultaneously. wmax . y0 (s) represents the sth sender (and receiver) center
An additional limitation of our current setup that can position on the y-axis (see Fig. 15). Another condition for
potentially be tackled with the application of detector arrays, an elementary ray to lie within the Gaussian beam is that its
is the limitation of the incident angle. For strongly diverted starting angle θ fulfills
beams, the receiver lenses fail to guide the radiation to
the detector antenna and therefore, method 2 does improve max (θmin1 , θmin2 ) ≤ θ ≤ min(θmax1 , θmax2 ) (14)
the image quality significantly. Applying detector arrays
with a very high acceptance angle per pixel can solve this where:
problem. 
wmax + d

Before an image can be reconstructed from the measure- θmin1 = tan −1
−
2xbound
ment data, the ray-tracing procedure has to be performed.  
w0 + d
This takes about 2 hours for 180 angles and 100 sampling θmin2 = tan−1 −
points on the projection axis to simulate the beam paths. One xbound
 
additional hour has to be invested to build the matrices Ā and w max − d
θmax1 = tan−1
ĀASM . These rough numbers were found on a standard 4-core 2xbound
 
processor, and could be reduced further by using general- w 0−d
θmax2 = tan −1
(15)
purpose graphics processing units (GPGPU) to perform the xbound

VOLUME 11, 2023 18323

 K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

The definitions of the angles are depicted in Fig. 15. The [17] L. A. Shepp and B. F. Logan, ‘‘The Fourier reconstruction of a head
combination of multiple elementary beams to form a Gaus- section,’’ IEEE Trans. Nucl. Sci., vol. NC-21, no. 3, pp. 21–43, Jun. 1974,
doi: 10.1109/TNS.1974.6499235.
sian beam supports our mathematical model of the interac- [18] K. H. May, A. Keil, G. Von Freymann, and F. Friederich, ‘‘The con-
tion of the radiation and the sample under test. The good jugate gradient least square algorithm in terahertz tomography,’’ IEEE
agreement between simulation and measurement, shown in Access, vol. 9, pp. 142168–142178, 2021, doi: 10.1109/ACCESS.2021.
3116801.
Section V-B. lets us conclude that it approximates the exper- [19] L. Ding, K. Xiao, P. Liu, J. Liu, and S. Chai, ‘‘Overview and classifi-
iment well enough. cation of a priori information used in microwave tomography,’’ in IEEE
MTT-S Int. Microw. Symp. Dig., Beijing, China, Jul. 2016, pp. 1–2, doi:
10.1109/NEMO.2016.7561616.
REFERENCES [20] K. C. Tam and L. J. Thomas, ‘‘Digital tomography incorporating a pri-
[1] E. Cristofani, F. Friederich, S. Wohnsiedler, C. Matheis, J. Jonuscheit, ori information,’’ in Review of Progress in Quantitative Nondestructive
M. Vandewal, and R. Beigang, ‘‘Nondestructive testing potential eval- Evaluation, D. O. Thompson and D. E. Chimenti, Eds. Boston, MA, USA:
uation of a terahertz frequency-modulated continuous-wave imager for Springer, 1990, pp. 407–414, doi: 10.1007/978-1-4684-5772-8_50.
composite materials inspection,’’ Opt. Eng., vol. 53, no. 3, Mar. 2014, [21] A. Zakaria, A. Baran, and J. LoVetri, ‘‘Estimation and use of prior
Art. no. 031211, doi: 10.1117/1.OE.53.3.031211. information in FEM-CSI for biomedical microwave tomography,’’ IEEE
[2] I. Amenabar, F. Lopez, and A. Mendikute, ‘‘In introductory review to THz Antennas Wireless Propag. Lett., vol. 11, pp. 1606–1609, 2012, doi:
non-destructive testing of composite mater,’’ J. Infr., Millim., Terahertz 10.1109/LAWP.2012.2237537.
Waves, vol. 34, no. 2, pp. 152–169, Feb. 2013, doi: 10.1007/s10762-012- [22] A. Fhager and M. Persson, ‘‘Using a priori data to improve the recon-
9949-z. struction of small objects in microwave tomography,’’ IEEE Trans.
[3] A. D’Arco, M. D. Di Fabrizio, V. Dolci, M. Petrarca, and S. Lupi, ‘‘THz Microw. Theory Techn., vol. 55, no. 11, pp. 2454–2462, Nov. 2007, doi:
pulsed imaging in biomedical applications,’’ Condens. Matter, vol. 5, no. 2, 10.1109/TMTT.2007.908670.
p. 25, Apr. 2020, doi: 10.3390/condmat5020025. [23] P. M. V. D. Berg and J. T. Fokkema, ‘‘Removal of undesired wave-
[4] Y. Peng, C. Shi, X. Wu, Y. Zhu, and S. Zhuang, ‘‘Terahertz imaging and fields related to the casing of a microwave scanner,’’ IEEE Trans.
spectroscopy in cancer diagnostics: A technical review,’’ BME Frontiers, Microw. Theory Techn., vol. 51, no. 1, pp. 187–192, Jan. 2003, doi:
vol. 2020, pp. 1–11, Jan. 2020, doi: 10.34133/2020/2547609. 10.1109/TMTT.2002.806900.
[5] H. Hintzsche, C. Jastrow, T. Kleine-Ostmann, U. Kärst, T. Schrader, and [24] G. R. Myers, A. M. Kingston, T. K. Varslot, M. L. Turner, and
H. Stopper, ‘‘Terahertz electromagnetic fields (0.106 THz) do not induce A. P. Sheppard, ‘‘Dynamic tomography with a priori information,’’
manifest genomic damage in vitro,’’ PLoS ONE, vol. 7, no. 9, Sep. 2012, Appl. Opt., vol. 50, no. 20, p. 3685, Jul. 2011, doi: 10.1364/AO.50.
Art. no. e46397, doi: 10.1371/journal.pone.0046397. 003685.
[6] M. Naftaly, N. Vieweg, and A. Deninger, ‘‘Industrial applications of tera- [25] B. Recur, J. P. Guillet, I. Manek-Hönninger, J. C. Delagnes, W. Benhar-
hertz sensing: State of play,’’ Sensors, vol. 19, no. 19, p. 4203, Sep. 2019, bone, P. Desbarats, J. P. Domenger, L. Canioni, and P. Mounaix, ‘‘Propa-
doi: 10.3390/s19194203. gation beam consideration for 3D THz computed tomography,’’ Opt. Exp.,
[7] M. Bauer, R. Hussung, C. Matheis, H. Reichert, P. Weichenberger, J. Beck, vol. 20, no. 6, p. 5817, Mar. 2012, doi: 10.1364/OE.20.005817.
U. Matuschczyk, J. Jonuscheit, and F. Friederich, ‘‘Fast FMCW terahertz [26] J. Tepe, T. Schuster, and B. Littau, ‘‘A modified algebraic reconstruction
imaging for in-process defect detection in press sleeves for the paper technique taking refraction into account with an application in terahertz
industry and image evaluation with a machine learning approach,’’ Sensors, tomography,’’ Inverse Problems Sci. Eng., vol. 25, no. 10, pp. 1448–1473,
vol. 21, no. 19, p. 6569, Sep. 2021, doi: 10.3390/s21196569. Oct. 2017, doi: 10.1080/17415977.2016.1267168.
[8] F. Friederich, K. May, B. Baccouche, C. Matheis, M. Bauer, J. Jonuscheit, [27] P. Fosodeder, S. Van Frank, and C. Rankl, ‘‘Highly accurate THz-CT
M. Moor, D. Denman, J. Bramble, and N. Savage, ‘‘Terahertz radome including refraction effects,’’ Opt. Exp., vol. 30, no. 3, p. 3684, Jan. 2022,
inspection,’’ Photonics, vol. 5, no. 1, p. 1, Jan. 2018, doi: 10.3390/pho- doi: 10.1364/OE.444151.
tonics5010001. [28] J. Kaipio and E. Somersalo, Statistical and Computational Inverse Prob-
[9] F. Ellrich, M. Bauer, N. Schreiner, A. Keil, T. Pfeiffer, J. Klier, S. Weber, lems. New York, NY, USA: Springer, 2005.
J. Jonuscheit, F. Friederich, and D. Molter, ‘‘Terahertz quality inspection [29] F. Natterer, The Mathematics of Computerized Tomography. Philadelphia,
for automotive and aviation industries,’’ J. Infr., Millim., Terahertz Waves, PA, USA: Society for Industrial and Applied Mathematics, 2001. [Online].
vol. 41, no. 4, pp. 470–489, Apr. 2020, doi: 10.1007/s10762-019-00639-4. Available: https://books.google.de/books?id=gjSO1hLbcD0C
[10] S. Mohammadzadeh, A. Keil, S. Leuchs, C. Krebs, D. NuBler, [30] A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imag-
J. Seewig, and F. Friederich, ‘‘Hand-guided mobile terahertz 3D ing. Philadelphia, PA, USA: Society for Industrial and Applied Mathemat-
imaging platform with aspherical telecentric f -θ optics,’’ in Proc. ics, 2001, doi: 10.1137/1.9780898719277.
18th Eur. Radar Conf. (EuRAD), Apr. 2022, pp. 377–380, doi: [31] W. Wu, D. Hu, K. An, S. Wang, and F. Luo, ‘‘A high-quality photon-
10.23919/EuRAD50154.2022.9784526. counting CT technique based on weight adaptive total-variation and
[11] T. Wang, K. Wang, K. Zou, S. Shen, Y. Yang, M. Zhang, Z. Yang, image-spectral tensor factorization for small animals imaging,’’ IEEE
and J. Liu, ‘‘Virtual unrolling technology based on terahertz computed Trans. Instrum. Meas., vol. 70, pp. 1–14, 2021, doi: 10.1109/TIM.2020.
tomography,’’ Opt. Lasers Eng., vol. 151, Apr. 2022, Art. no. 106924, doi: 3026804.
10.1016/j.optlaseng.2021.106924. [32] A. Bjorck, Numerical Methods for Least Squares Problems. Philadel-
[12] K. Krugener, E.-M. Stubling, R. Jachim, B. Kietz, M. Koch, and phia, PA, USA: SIAM Society for Industrial and Applied Mathematics,
W. Viöl, ‘‘THz tomography for detecting damages on wood caused 1996.
by insects,’’ Appl. Opt., vol. 58, no. 22, p. 6063, Aug. 2019, doi: [33] P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems: Numer-
10.1364/AO.58.006063. ical Aspects of Linear Inversion. Philadelphia, PA, USA: SIAM,
[13] J. B. Perraud, A. F. Obaton, J. Bou-Sleiman, B. Recur, H. Balacey, 1998.
F. Darracq, J. P. Guillet, and P. Mounaix, ‘‘Terahertz imaging and tomog- [34] M. R. Hestenes and E. Stiefel, ‘‘Methods of conjugate gradients for solving
raphy as efficient instruments for testing polymer additive manufac- linear systems,’’ J. Res. Nat. Bureau Standards, vol. 49, no. 6, p. 409, 1952,
turing objects,’’ Appl. Opt., vol. 55, no. 13, p. 3462, May 2016, doi: doi: 10.6028/jres.049.044.
10.1364/AO.55.003462. [35] P. C. Hansen, ‘‘The L-curve and its use in the numerical treatment of
[14] P. Fosodeder, S. Hubmer, A. Ploier, R. Ramlau, S. van Frank, and C. Rankl, inverse problems,’’ in Computational Inverse Problems in Electrocardiol-
‘‘Terahertz computer tomography for plastics extrusion (TACTICS),’’ pre- ogy. Cambridge, MA, USA: MIT Press, 2000.
sented at the ATTRACT Final Conf., 2020. [36] M. Naftaly, ‘‘Terahertz optical materials,’’ in Proc. Int. Traveling Summer
[15] S. Wang and X.-C. Zhang, ‘‘Pulsed terahertz tomography,’’ J. Phys. D, School, Frankfurt, Germany, 2019, pp. 25–49.
Appl. Phys., vol. 37, no. 4, pp. 1–36, Feb. 2004, doi: 10.1088/0022- [37] R. Hussung, A. Keil, and F. Friederich, ‘‘Handheld millimeter wave
3727/37/4/R01. imaging system based on a two-dimensional multistatic sparse array,’’
[16] J. P. Guillet, ‘‘Review of terahertz tomography techniques,’’ J. Infr., in Proc. 45th Int. Conf. Infr., Millim., Terahertz Waves (IRMMW-
Millim., Terahertz Waves, vol. 35, no. 4, pp. 382–411, 2014, doi: THz), Buffalo, NY, USA, Nov. 2020, pp. 1–2, doi: 10.1109/IRMMW-
10.1007/s10762-014-0057-0. THz46771.2020.9370732.

18324 VOLUME 11, 2023

K. H. May et al.: Priori Information and Off-Axis Measurements in Terahertz Tomography

KARL H. MAY was born in Mannheim, Germany, GEORG VON FREYMANN graduated in physics
in 1992. He received the Dipl.-Phys. degree from Universität Karlsruhe (TH), Karlsruhe,
in physics from Technische Universität Kaiser- Germany, in 1998. He received the Dr.rer.-nat.
slautern, Kaiserslautern, Germany, in 2017, work- degree in physics from TH, in 2001.
ing on sensor fusion of electronic terahertz FMCW In 2002, he was a Postdoctoral Researcher at the
systems, where he is currently pursuing the Ph.D. Institute of Nanotechnology, Forschungszentrum
degree in development of terahertz tomography Karlsruhe, Germany. From 2003 to 2004, he was
systems and software. a Postdoctoral Researcher at the University of
He joined the Fraunhofer Institute for Industrial Toronto, Canada. From 2005 to 2010, he headed
Mathematics ITWM, Kaiserslautern, in 2013, as a an independent DFG Emmy-Noether Research
Research Assistant at the Optimization Department. In 2016, he switched Group, Institute of Nanotechnology, Karlsruhe Institute of Technology
to the Materials Characterization and Testing Department to pursue his (KIT). Since 2010, he has been a Full Professor of experimental physics
Diploma degree. His research interests include (electronic) terahertz trans- with the Technische Universität Kaiserslautern, Germany. Since 2013, he has
mitters and detectors, broadband antennas, terahertz measurement and appli- also been heads the Department of Materials Characterization and Testing,
cations, data fusion, and tomography reconstruction algorithms. Fraunhofer Institute for Industrial Mathematics ITWM. He is the author of
more than 110 articles and more than 15 inventions. He is the co-founder of
two companies. His research interests include 3-D laser-lithography, nano-
photonics, optical quantum simulators, terahertz technology, and spin-wave
optics.

FABIAN FRIEDERICH received the degree in

optoelectronics from Aalen University Sciences,
Germany, in 2007, and the Ph.D. degree in
physics from Goethe University Frankfurt am
Main, Germany, for his work in terahertz imag-
ing with the Ultrafast Spectroscopy and Terahertz
Physics Group, in 2012.
During his studies, he was also associated with
the Centre for Micro-Photonics, Swinburne Uni-
versity of Technology, Melbourne, Australia, and
the Laser Zentrum Hannover, Germany. In 2011, he joined the Institute of
Technical Physics, German Aerospace Center (DLR), Stuttgart, Germany,
ANDREAS KEIL (Member, IEEE) was born in to establish laser-based concepts for monitoring space debris. In 2013, he was
Kiel, Germany, in 1975. He received the degree in granted a Fraunhofer Attract Funding to form a new research group in
physics from Christian Albrechts University, Kiel, millimeter-wave and terahertz measurement techniques at the Fraunhofer
in 2003, and the Ph.D. degree in quantum infor- Institute for Physical Measurement Techniques IPM. Since 2017, his group
mation from the National University of Singapore, has been with the Fraunhofer Institute for Industrial Mathematics ITWM,
Singapore, in 2010. Kaiserslautern, pursuing its millimeter-wave and terahertz activities with an
Since 2009, he has been involved in the research even stronger focus on signal and image processing in nondestructive testing.
and development of terahertz imaging systems; He is an Active Member of the Microwave and Terahertz Methods Technical
where he is currently with Fraunhofer ITWM. His Committee of the German Society for Non-Destructive Testing (DGZfP) and
main research interests include terahertz imag- the VDI/VDE-GMA Technical Committee 8.17 Terahertz Systems of the
ing, image processing algorithms, and high-performance computation on Association of German Engineers (VDI).
general-purpose graphics processing units.

VOLUME 11, 2023 18325