# Seeing through a Black Box: Toward High-Quality Terahertz Tomographic
Imaging via Multi-Scale Spatio-Spectral Image Fusion

Weng-tai Su†, Yi-Chun Hung†, Ta-Hsuan Chao†, Po-Jen Yu†, Shang-Hua Yang† and
Chia-Wen Lin†
†National Tsing Hua University, Taiwan

###### Abstract

Terahertz (THz) imaging has recently attracted significant attention thanks to
its non-invasive, non-destructive, non-ionizing, material-classification, and
ultra-fast nature for object exploration and inspection. However, its strong
water absorption nature and low noise tolerance lead to undesired blurs and
distortions of reconstructed THz images. The performances of existing
restoration methods are highly constrained by the diffraction-limited THz
signals. To address the problem, we propose a novel Subspace-and-Attention-
guided Restoration Network (SARNet) that fuses multi-spectral features of a
THz image for effective restoration. To this end, SARNet uses multi-scale
branches to extract spatio-spectral features of amplitude and phase which are
then fused via shared subspace projection and attention guidance. Here, we
experimentally construct ultra-fast THz time-domain spectroscopy system
covering a broad frequency range from 0.1 THz to 4 THz for building up
temporal/spectral/spatial/phase/material THz database of hidden 3D objects.
Complementary to a quantitative evaluation, we demonstrate the effectiveness
of our SARNet model on 3D THz tomographic reconstruction applications.

##  1 Introduction

Ever since the first camera’s invention, imaging under different bands of
electromagnetic (EM) waves, especially X-ray and visible lights, has
revolutionized our daily lives [17, 30, 39]. X-ray imaging plays a crucial
role in medical diagnosis, such as cancer, odontopathy, and COVID-19 symptom
[1, 31, 35], based on X-ray’s high penetration depth to great varieties of
materials; visible-light imaging has not only changed the way of recording
lives but contributes to the development of artificial intelligence (AI)
applications, such as surveillance security and surface defect inspection
[38]. However, X-ray and visible-light imaging still face tough challenges.
X-ray imaging is ionizing, which is harmful to biological objects and thus
severely limits its application scope [6]. On the other hand, although both
non-ionizing and non-destructive, visible-light imaging cannot retrieve most
optically opaque objects’ interior information due to the highly absorptive
and intense scattering behaviors between light and matter in the visible light
band. To visualize the 3D information of objects in a remote but accurate
manner, terahertz (THz) imaging has become among the most promising candidates
among all EM wave-based imaging techniques [3, 4].

Figure 1: Illustration of THz data collection with our in-house THz-TDS
tomographic imaging system. Figure 2: Illustration of THz 3D tomographic
imaging based on SARNet. The red block shows the images retrieved from
different spectral frequencies of the amplitude and phase of the THz
measurements of Deer.

THz radiation, between microwave and infrared, has often been regarded as the
last frontier of EM wave [32], which provides its unique functionalities among
all EM bands. Along with the rapid development of THz technology, THz imaging
has recently attracted significant attention due to its non-invasive, non-
destructive, non-ionizing, material-classification, and ultra-fast nature for
advanced material exploration and engineering. As THz waves can partially
penetrate through varieties of optically opaque materials, it carries hidden
material tomographic information along the traveling path, making this
approach a desired way to see through black boxes without damaging the
exterior [22, 14, 23]. By utilizing light-matter interaction within the THz
band, multifunctional tomographic information of a great variety of materials
can also be retrieved even at a remote distance. In the past decades, THz
time-domain spectroscopy (THz-TDS) has become one of the most representative
THz imaging modalities to achieve non-invasive evaluation because of its
unique capability of extracting geometric and multi-functional information of
objects. Owing to its fruitful information in multi-dimensional domains —
space, time, frequency, and phase, THz-TDS imaging has been already allocated
for numerous emerging fields, including drug detection [19], industrial
inspection, cultural heritage inspection [8] and cancer detection [2].

To retrieve temporal-spatio-spectral information of each object voxel, our THz
imaging experiment setup is based on a THz-TDS system as shown in Fig. 1. Our
measured object (a covered 3D printed deer) is placed on the rotation stage in
the THz path between the THz source and detector of the THz-TDS system. During
the scanning, the THz-TDS system profiles each voxel’s THz temporal signal
with 0.1 ps temporal resolution, whose amplitude corresponds to the strength
of THz electric field. Based on the dependency between the amplitude of a
temporal signal and THz electric field, in conventional THz imaging, the
maximum peak of the signal (Time-max) is extracted as the feature for a voxel.
The reconstructed image based on Time-max features can deliver great signal-
to-noise ratio and a clear object contour. However, the conventional THz
imaging based on Time-max features suffers from several drawbacks, such as the
undesired contour in the boundary region, the hollow in the body region, and
the blurs in high spatial-frequency regions. To break this limitation, we
utilize the spectral information of THz temporal signals to supplement the
conventional method based on Time-max features since the voxel of the material
behaviors are encoded in both the phase and amplitude of different frequency
components, according to the Fresnel equation [7].

Due to the large number of spectral bands with measured THz image data, it is
required to sample a subset of the spectral bands to reduce the training
parameters. The THz beam is significantly attenuated at water absorption
frequencies. Thus, reconstructed THz images based on water absorption lines
offer worse details. Besides, the high-speed THz-TDS system offers more than
20 dB SNR in a frequency range of 0.3 THz–1.3 THz. Considering the water
absorption in THz regime [36, 34] and the superior SNR in the range of 0.3
THz–1.3 THz, we select 12 frequencies at 0.380, 0.448, 0.557, 0.621, 0.916,
0.970, 0.988, 1.097, 1.113, 1.163, 1.208, and 1.229 THz. The spectral
information including both amplitude and phase at the selected frequencies is
extracted and then employed to restore clear 2D images. Specifically, Fig. 2
shows multiple 2D THz images of the same object at the selected frequencies,
showing very different contrasts and spatial resolutions as these
hyperspectral THz image sets have different physical characteristics through
the interaction of THz waves with objects. The lower-frequency phase images
offer relatively accurate depth information due to their higher SNR level,
whereas the higher-frequency phase images offer finer contours and edges
because of the shrinking diffraction-limited wavelength sizes (from left to
right in Fig. 2). The phase also contains, however, a great variety of
information of light-matter interaction that could cause learning difficulty
for the image restoration task. To address this issue, we utilize amplitude
spectrum as complementary information. Although the attenuated amplitude
spectrum cannot reflect comparable depth accuracy levels as phase spectrum,
amplitude spectrum still present superior SNR and more faithful contours such
as the location information of a measured object. Besides, as the
complementary information to phase, the lower-frequency amplitude offers
higher contrast, whereas the higher-frequency amplitude offers a better object
mask.

In summary, the amplitude complements the shortcomings of phase. The
advantages of fusing the two signals from low-frequency to high-frequency are
as follows: Since the low-frequency THz signal provides precise depth (the
thickness of an object) and fine edge/contour information in the phase and
amplitude, respectively, they together better delineate and restore the
object. In contrast, the high-frequency feature maps of amplitude and phase
respectively provide better edges/contours and precise position information,
thereby constituting a better object mask from the complementary features.
With these multi-spectral properties of THz images, we can extract rich
information from a wide spectral range in the frequency domain to
simultaneously restore the 2D THz images without any additional computational
cost or equipment, which is beneficial for further development of THz imaging.

We here propose a multi-scale Subspace-and-Attention guided Restoration
Network (SARNet) that fuses complementary spectral features of the THz
amplitude and phase to supplement the Time-max image for restoring clear 2D
images. To this end, SARNet learns common representations in a common latent
subspace shared between the amplitude and phase, and then incorporates a Self-
Attention mechanism to learn the short/long-range dependency of the spectral
features for guiding the restoration task. Finally, from clear 2D images
restored from corrupted images of an object captured from different angles, we
can reconstruct high-quality 3D tomography via inverse Radon transform. Our
main contributions are summarized as follows:

  * •

We are the first research group to merge THz temporal-spatio-spectral data,
data-driven models, and material properties to the best of our knowledge. The
proposed SARNet has an excellent performance in extracting and fusing features
from the light-matter interaction data in THz spectral regime, which
inherently contains fruitful 3D object information and its material behaviors.
Based on the delicately designed architecture of SARNet on feature fusion,
SARNet delivers state-of-the-art performance on THz image restoration.

  * •

With our newly established THz tomography dataset — the world’s first in its
kind, we provide comprehensive quantitative/qualitative analyses among SARNet
and SOTAs. SARNet significantly outperforms Time-max, baseline U-Net, and
multi-spectral U-Net by 11.17dB, 2.86dB, and 1.51dB in average PSNR.

  * •

This pioneering work shows that computer vision techniques can significantly
contribute to the THz community and further open up a new interdisciplinary
research field to boost practical applications, e.g., non-invasive evaluation,
gas tomography, industrial inspection, material exploration, and biomedical
imaging.

##  2 Related Work

###  2.1 Deep Learning-based Image Restoration

In recent years, deep learning methods were first popularized in high-level
visual tasks, and then gradually penetrated into many tasks such as image
restoration and segmentation. Convolutional neural network (CNNs) have proven
to achieve the state-of-the-art performances in fundamental image restoration
problem [41, 20, 43, 42, 29]. Several network models for image restoration
were proposed, such as U-Net [29], hierarchical residual network [20] and
residual dense network [43]. Notably, DnCNN [41] uses convolutions, BN, and
ReLU to build 17-layer network for image restoration which was not only
utilized for blind image denoising, but was also employed for image super-
resolution and JPEG image deblocking. FFDNet [42] employs noise level maps as
inputs and utilizes a single model to develop variants for solving problems
with multiple noise levels. In [20] a very deep residual encoding-decoding
(RED) architecture was proposed to solve the image restoration problem using
skip connections. [43] proposed a residual dense network (RDN), which
maximizes the reusability of features by using residual learning and dense
connections. NBNet [5] employs subspace projection to transform learnable
feature maps into the projection basis, and leverages non-local image
information to restore local image details. Similarly, the Time-max image
obtained from a THz imaging system can be cast as an image-domain learning
problem which was rarely studied due to the difficulties in THz image data
collection. Research works on image-based THz imaging include [25, 26, 37],
and THz tomographic imaging works include [12, 11].

Figure 3: Network architecture of SARNet consisting of five scale branches,
where the finest-scale scale takes the Time-max image as input, and each of
the second to fifth takes 6 images of spectral frequencies (3 bands in
amplitude and 3 bands in phase) as inputs. The two gray blocks show the
detailed structures of the Conv-block and the Channel Attention Module (CAM).

###  2.2 Tomographic Reconstruction

Computer tomographic (CT) imaging methods started from X-ray imaging, and many
methods of THz imaging are similar to those of X-ray imaging. One of the first
works to treat X-ray CT as an image-domain learning problem was [18], that
adopts CNN to refine tomographic images. In [15], U-Net was used to refine
image restoration with significantly improved performances. [44] further
projects sinograms measured directly from X-ray into higher-dimensional space
and uses domain transfer to reconstruct images. The aforementioned works were
specially designed for X-ray imaging.

Hyperspectral imaging [33, 24, 9] constitutes image modalities other than THz
imaging. Different from THz imaging, Hyperspectral imaging collects continuous
spectral band information of the target sample. Typically, the frequency bands
fall in the visible and infrared spectrum; hence, most hyperspectral imaging
modalities can only observe the surface characteristics of targeted objects.

Figure 4: Block diagram of Subspace-and-Attention guided Fusion Module (SAFM).

##  3 Terahertz Image Restoration

###  3.1 Overview

As different EM bands interact with objects differently, THz waves can
partially penetrate through various optically opaque materials and carry
hidden material tomographic information along the traveling path. This unique
feature provides a new approach to visualize the essence of 3D objects, which
other imaging modalities cannot achieve. Although existing deep neural
networks can learn spatio-spectral information from a considerable amount of
spectral cube data, we found that directly learning from the full spectral
information to restore THz images usually leads to an unsatisfactory
performance. The main reason is that the full spectral bands of THz signals
involve diverse characteristics of materials, noises, and scattered signal,
which causes difficulties in model training. To address this problem, our work
is based on extracting complementary information from both the amplitude and
phase of a THz signal. That is, as illustrated in Fig. 2, in the low-frequency
bands, the amplitude images delineate finer edges and object contours while
the phase images offer relatively precise depths of object surfaces. In
contrast, in the high-frequency bands, the amplitude images offer object mask
information while the phase images delineate finer edges and object contours.
Therefore, the amplitude and phase complement to each other in both the low
and high frequency bands.

Motivated by the above findings, we devise a novel multi-scale Subspace-and-
Attention-guided Restoration Network (SARNet) to capture such complementary
spectral characteristics of materials to restore damaged 2D THz images
effectively. The key idea of SARNet is to fuse spatio-spectral features with
different characteristics on a common ground via deriving the shared latent
subspace and discovering the short/long-range dependencies between the
amplitude and phase to guide the feature fusion. To this end, SARNet is based
on U-Net[29] to perform feature extraction and fusion in a multi-scale manner.

###  3.2 Network Architecture

On top of U-Net [29], the architecture of SARNet is depicted in Fig. 3.
Specifically, SARNet is composed of an encoder with 5 branches of different
scales (from the finest to the coarsest) and a decoder with 5 corresponding
scale branches. Each scale branch of the encoder involves a Subspace-and-
Attention-guided Fusion module (SAFM), a convolution block (Conv-block), and a
down-sampler, except for the finest-scale branch that does not employ SAFM. To
extract and fuse multi-spectral features of both amplitude and phase in a
multi-scale manner, the encoder takes a THz 2D Time-max image as the input of
the finest-scale branch as well as receives to its second to fifth scale
branches 24 images of additional predominant spectral frequencies extracted
from the THz signal, where each branch takes 6 images of different spectral
bands (3 bands of amplitude and 3 bands of phase) to extract learnable
features from these spectral bands. To reduce the number of model parameters,
these 24 amplitude and phase images (from low to high frequencies) are
downsampled to 4 different resolutions and fed into the second to fifth scale
branches in a fine-to-coarse manner as illustrated in Fig. 3. We then fuse the
multi-spectral amplitude and phase feature maps in each scale via the proposed
SAFM that learns a common latent subspace shared between the amplitude and
phase features to facilitate associating the short/long-range amplitude-phase
dependencies. Projected into the shared latent subspace, the spectral features
of amplitude and phase components, along with the down-sampled features of the
upper layer, can then be properly fused together on a common ground in a fine-
to-coarse fashion to obtain the final latent code.

The Conv-block($L$) contains two stacks of $L\times L$ convolution, batch
normalization, and ReLU operations. Because the properties of the spectral
bands of amplitude and phase can be very different, we partly use $L=1$ to
learn the best linear combination of multi-spectral features to avoid noise
confusion and reduce the number of model parameters. The up-sampler and down-
sampler perform $2\times$ and $\frac{1}{2}\times$ scaling, respectively. The
skip connections (SC) directly pass the feature maps of different spatial
scales from individual encoder branches to the Channel Attention Modules
(CAMs) of their corresponding branches of the decoder, as indicated by the
gray arrows in Fig. 3. The details of SAFM and CAM will be elaborated later.

In the decoder path, each scale branch involves a up-sampler, a CAM, and a
Conv-block. The Conv-block has the same functional blocks as that in the
encoder. Each decoding-branch receives a “shallower-layer” feature map from
the corresponding encoding-branch via the skip-connection shortcut and
concatenates the feature map with the upsampled version of the decoded
“deeper-layer” feature map from its coarser-scale branch. Besides, the
concatenated feature map is then processed by CAM to capture the cross-channel
interaction to complement the local region for restoration.

Note, a finer-scale branch of SARNet extracts shallower-layer features which
tend to capture low-level features, such as colors and edges. To complement
the Time-max image for restoration, we feed additional amplitude and phase
images of low to high spectral-bands into the fine- to coarse-scale branches
of SARNet. Since the spectral bands of THz amplitude and phase offer
complementary information, as mentioned in Sec. 3.1, besides the Time-max
image SARNet also extracts multi-scale features from the amplitude and phase
images of 12 selected THz spectral bands, which are then fused by the proposed
SAFM.

###  3.3 Subspace-and-Attention guided Fusion Module

How to properly fuse the spectral features of THz amplitude and phase are,
however, not trivial, as their characteristics are very different. To address
the problem, inspired by [5] and [40], we propose the SAFM as shown in Fig. 4.

We denote ${\mathbf{X}}^{A}_{\text{in}}$, ${\mathbf{X}}^{P}_{\text{in}}$ $\in$
$\mathbb{R}^{H\times W\times 3}$ as the spectral features of the THz amplitude
and phase, respectively. The Conv-block $f_{C}(\cdot)$ extracts two
intermediate feature maps $f_{C}({\mathbf{X}}^{A}_{\text{in}})$,
$f_{C}({\mathbf{X}}^{P}_{\text{in}})$ $\in$ $\mathbb{R}^{H\times W\times
C_{1}}$ from ${\mathbf{X}}^{A}_{\text{in}}$ and
${\mathbf{X}}^{P}_{\text{in}}$, respectively. As a result, we then derive the
$\mathbf{k}$ shared basis vectors
$\textbf{V}=[\textbf{v}_{1},\textbf{v}_{2},...,\textbf{v}_{K}]$ from
$f_{C}({\mathbf{X}}^{A}_{\text{in}})$ and
$f_{C}({\mathbf{X}}^{P}_{\text{in}})$, where $\textbf{V}\in\mathbb{R}^{N\times
K}$, $N=HW$ denotes the dimension of each basis vector, and $K$ is the rank of
the shared subspace. The basis set of the shared common subspace is expressed
as

$\textbf{V}=f_{F}(f_{C}({\mathbf{X}}^{A}_{\text{in}}),f_{C}({\mathbf{X}}^{P}_{\text{in}})),$
(1)

where we first concatenate the two feature maps in the channel dimension and
then feed the concatenated feature into the fusion-block $f_{F}(\cdot)$. The
structure of the fusion-block is the same as that of the Conv-block with $K$
output channels as indicated in the red block in Fig. 4. The weights of the
fusion-block are learned in the end-to-end training stage. The shared latent
subspace learning mainly serves two purposes: 1) learning the common
representation between the THz amplitude and phase, and 2) learning the
subspace projection matrix to project the amplitude and phase features into a
shared subspace such that they can be analyzed on a common ground. These both
help identify short/long-range dependencies of amplitude and phase features
for feature fusion in the next stage.

To identify the short/long-range dependencies between the amplitude and phase
features on a common ground, we utilize the orthogonal projection matrix V in
(1) to estimate the self-attentions in the shared feature subspace as

$\beta_{j,i}=\frac{\exp(s_{ij})}{\sum_{i=1}^{N}\exp(s_{ij})}~{},~{}s_{ij}=\mathbf{v}_{i}^{T}\mathbf{v}_{j}$
(2)

where $\beta_{j,i}$ represents the model attention in the $i$-th location of
the $j$-th region. The orthogonal projection matrix P is derived from the
subspace basis V as follows [21]:

$\textbf{P}=\textbf{V}(\textbf{V}^{T}\textbf{V})^{-1}\textbf{V}^{T}$ (3)

where $(\textbf{V}^{T}\textbf{V})^{-1}$ is the normalization term to make the
basis vectors orthogonal to each other during the basis generation process. As
a result, the output of the self-attention mechanism becomes

$\mathbf{o}_{j}=\left(\sum_{i=1}^{N}\beta_{j,i}\mathbf{s}_{i}\right),~{}~{}\mathbf{s}_{i}=\texttt{Concate}(\textbf{P}{\mathbf{X}}^{A}_{\text{in}},\textbf{P}{\mathbf{X}}^{P}_{\text{in}})$
(4)

where the key of $\mathbf{s}_{i}$ $\in$ $\mathbb{R}^{HW\times 6}$ is obtained
by concatenating the two feature maps $\textbf{PX}^{A}_{\text{in}}$ and
$\textbf{PX}^{P}_{\text{in}}$ projected by orthogonal projection matrix P
$\in$ $\mathbb{R}^{HW\times HW}$, and ${\mathbf{X}}^{A}_{\text{in}}$ and
${\mathbf{X}}^{P}_{\text{in}}$ are reshaped to $HW\times 3$. Since the
operations are purely linear with some proper reshaping, they are
differentiable.

Finally, we further fuse cross-scale features in the self-attention output by
adding the down-sampled feature map $\mathbf{X}_{f}$ from the finer scale as

$\mathbf{X_{out}}=f_{s}(\mathbf{o})+\mathbf{X}_{f}$ (5)

where $f_{s}$ is the 1×1 convolution to keep the channel number consistent
with $\mathbf{X}_{f}$.

###  3.4 Channel Attention Module

To fuse multi-scale features from different spectral bands in the channel
dimension, we incorporate the efficient channel attention mechanism proposed
in [27] in the decoder path of SARNet as shown in Fig. 3. In each decoding-
branch, the original U-Net directly concatenates the up-sampled feature from
the coarser scale with the skip-connection feature from the corresponding
encoding-branch, and then fuse the intermediate features from different layers
by convolutions. This, however, leads to poor image restoration performances
in local regions such as incorrect object thickness or details. To address
this problem, we propose a channel attention module (CAM) which adopts full
channel attention in the dimensionality reduction operation by concatenating
two channel attention groups. CAM first performs the global average pooling to
extract the global spatial information in each channel:

$\displaystyle G_{t}=\frac{1}{H\times
W}\sum_{i=1}^{H}\sum_{j=1}^{W}X_{t}(i,j)$ (6)

where $X_{t}(i,j)$ is the $t$-th channel of $X_{t}$ at position $(i,j)$
obtained by concatenating the up-sampled feature map $\textbf{X}_{c}$ of the
coarser-scale and the skip-connection feature map $\textbf{X}_{s}$. The shape
of $G$ is from $C\times H\times W$ to $C\times 1\times 1$.

We directly feed the result through two $1\times 1$ convolution, sigmoid, and
ReLU activation function as:

$\displaystyle\textbf{w}=\sigma\left(\texttt{Conv}\left(\delta\left(\texttt{Conv}(G)\right)\right)\right),$
(7)

where Conv represents the $1\times 1$ convolution, $\sigma$ is the sigmoid
function, and $\delta$ is the ReLU function. In order to better restore a
local region, we group the weights w of different channels as two groups
$\textbf{w}=[\textbf{w}_{1},\textbf{w}_{2}]$ corresponding to two different
sets of input feature maps. Finally, we element-wise multiply the input
$X_{c}$, $X_{s}$ of the weights w and add these two groups features.

###  3.5 Loss Function for THz Image Restoration

To effectively train SARNet, we employ the following mean squared error (MSE)
loss function to measure the dissimilarity between the restored image
$\mathbf{X}_{\text{rec}}$ and its ground-truth $\mathbf{X}_{\text{GT}}$:

$\displaystyle\mathcal{L}_{\text{MSE}}(\mathbf{X}_{\text{GT}},\mathbf{X}_{\text{rec}})=$
$\displaystyle\frac{1}{HW}\sum_{i=1}^{H}\sum_{j=1}^{W}(\mathbf{X}_{\text{GT}}(i,j)-\mathbf{X}_{\text{rec}}(i,j))^{2},$
(8)

where $H$ and $W$ are the height and width of the image.

The 3D tomography of an object can then be reconstructed from the restored THz
2D images of the object scanned in different angles. To this end, we can
directly apply the inverse Radon transform to obtain the 3D tomography, using
methods like filtered back-projection [16] or the simultaneous algebraic
reconstruction technique [28].

##  4 Experiments

Figure 5: Illustration of authentic measured THz data for the seven 3D-printed
HIPS objects in our experiments. Figure 6: Qualitative comparison of THz image
restoration results for Deer, Box, and Robot from left to right: (a) Time-max,
(b) DnCNN-S [41], (c) RED [20], (d) NBNet [5], (e)
$\texttt{U-Net}_{\texttt{base}}$ [29], (f) $\texttt{U-Net}_{\texttt{MS}}$, (g)
SARNet, and (h) the ground-truth.

We conduct experiments to evaluate the effectiveness of SARNet against
existing state-of-the-art restoration methods. We first present our experiment
settings and then evaluate the performances of SARNet and the competing
methods on THz image restoration.

###  4.1 THz-TDS Image Dataset

As shown in Fig. 1, we prepare the sample objects by a Printech 3D printer,
and use the material of high impact polystyrene (HIPS) for 3D-printing the
objects due to its high penetration of THz waves. We then use our in-house
Asynchronous Optical Sampling THz time-domain spectroscopy (ASOPS THz-TDS)
system [13] to measure the sample objects. Each sample object is placed on a
motorized stage between the source and the receiver. With the help of the
motorized stage, raster scans are performed on each object in multiple view
angles. In the scanning phase, we scan the objects covering a rotational range
of 180 degrees (step-size: 6 degrees), a horizontal range of 72mm (step-size:
0.25mm), and a variable vertical range corresponding to the object height
(step-size: 0.25mm). In this way, we obtain 30 projections of each object,
which are then augmented to 60 projections by horizontal flipping. The ground-
truths of individual projections are obtained by converting the original 3D
printing files into image projections in every view-angle. We use markers to
indicate the center of rotation so that we can align the ground-truths with
the measured THz data. In this paper, totally 7 objects are printed, measured,
and aligned for evaluation.

###  4.2 Data Processing and Augmentation

In our experiments, we train SARNet using the 2D THz images collected from our
THz imaging system shown in Fig. 2. The seven sample objects are illustrated
in Fig. 5, consisting of 60 projections per object and 420 2D THz images in
total. In order to thoroughly evaluate the capacity of SARNet, we adopt the
leave-one-out strategy: using the data of 7 objects as the training set, and
that of the remaining object as the testing set. Due to the limited space, we
only present part of the results in this section, and the complete results in
the supplementary material. We will release our code and the THz image dataset
after the work is accepted.

We also perform typical data augmentations to enrich the training set,
involving random transformations on brightness, contrast, and saturation,image
scaling with a scaling factor ranging in $[0.5,1.5]$, and geometric
transformations including random flipping horizontally and vertically.
Finally, the images are randomly cropped to $128\times 128$ patches.

####  4.2.1 Experiment Settings

We initialize SARNet following the initialization method in [10], and train it
using the Adam optimizer with $\beta_{1}=0.9$ and $\beta_{2}=0.999$. We set
the initial learning rate to $10^{-4}$ and then decay the learning rate by
$0.1$ every $300$ epochs. SARNet converges after $1,000$ epochs. For a fair
comparison with the competing methods, we adopt their publicly released codes.

###  4.3 Quantitative Evaluations

To the best of our knowledge, there is no method specially designed for
restoring THz images besides Time-max. Thus, we compare our method against
several representative CNN-based image restoration models, including DnCNN
[41], RED [20], and NBNet [5]. Moreover, we also compare two variants of U-Net
[29]: baseline U-Net ($\texttt{U-Net}_{\texttt{base}}$) and multi-spectral
U-Net ($\texttt{U-Net}_{\texttt{MS}}$). $\texttt{U-Net}_{\texttt{base}}$
extracts image features in five different scales following the original
setting in U-Net [29], whereas $\texttt{U-Net}_{\texttt{MS}}$ incorporates
multi-spectral features by concatenating the features of Time-max image with
additional 12 THz bands for amplitude as the input (i.e., $12+1$ channels) of
the finest scale of U-Net. For objective quality assessment, we adopt two
widely-used metrics including the Peak Signal-to-Noise Ratio (PSNR) and
Structural SIMilarity (SSIM) to respectively measure the pixel-level and
structure-level similarities between a restored image and its ground-truth.

Table 1 shows that our SARNet significantly outperforms the competing methods
on all sample objects in both metrics. Specifically, SARNet outperforms Time-
max, baseline U-Net ($\texttt{U-Net}_{\texttt{base}}$), and the multi-spectral
U-Net ($\texttt{U-Net}_{\texttt{MS}}$) by 11.17 dB, 2.86 dB, and 1.51 dB in
average PSNR, and 0.744, 0.170, and 0.072 in average SSIM for 7 objects. For
qualitative evaluation, Fig. 6 illustrates a few restored views for Deer, Box,
and Robot, demonstrating that SARNet can restore objects with much finer and
smoother details (e.g., the antler and legs of Deer, the depth and shape of
Box, and the body of Robot), faithful thickness of material (e.g., the body
and legs of Deer and the correct edge thickness of Box), and fewer artifacts
(e.g., holes and broken parts). Both the quantitative and qualitative
evaluations confirm a significant performance leap with SARNet over the
competing methods.

Table 1: Quantitative comparison (PSNR and SSIM) of THz image restoration
performances with different methods on Deer, DNA, Box, Eevee, Polarbear,
Robot, and Skull. ($\uparrow$: higher is better)

Method |  PSNR$\uparrow$
---|---
| Deer | DNA | Box | Eevee | Polarbear | Robot | Skull
Time-max | 12.42 | 12.07 | 11.97 | 11.20 | 11.21 | 11.37 | 10.69
DnCNN-S [41] | 19.94 | 23.95 | 19.13 | 19.69 | 19.44 | 19.72 | 17.33
RED [20] | 19.30 | 24.17 | 20.18 | 19.97 | 19.17 | 19.76 | 16.28
NBNet [5] | 20.24 | 25.10 | 20.21 | 19.84 | 20.12 | 20.01 | 19.69
$\texttt{U-Net}_{\texttt{base}}$ [29] | 19.84 | 24.15 | 19.77 | 19.95 | 19.09 | 18.80 | 17.49
$\texttt{U-Net}_{\texttt{MS}}$ | 22.46 | 25.05 | 20.81 | 20.34 | 19.86 | 20.64 | 19.43
SARNet (Ours) | 22.98 | 26.05 | 22.67 | 20.87 | 21.42 | 22.66 | 22.48
Method |  SSIM$\uparrow$
| Deer | DNA | Box | Eevee | Polarbear | Robot | Skull
Time-max | 0.05 | 0.05 | 0.14 | 0.14 | 0.12 | 0.08 | 0.09
DnCNN-S [41] | 0.73 | 0.77 | 0.73 | 0.72 | 0.63 | 0.77 | 0.36
RED [20] | 0.81 | 0.83 | 0.74 | 0.77 | 0.75 | 0.80 | 0.74
NBNet [5] | 0.81 | 0.85 | 0.75 | 0.77 | 0.43 | 0.80 | 0.78
$\texttt{U-Net}_{\texttt{base}}$ [29] | 0.55 | 0.78 | 0.77 | 0.76 | 0.56 | 0.76 | 0.51
$\texttt{U-Net}_{\texttt{MS}}$ | 0.76 | 0.73 | 0.78 | 0.76 | 0.78 | 0.79 | 0.78
SARNet (Ours) | 0.84 | 0.90 | 0.83 | 0.82 | 0.82 | 0.83 | 0.84

Table 2: Quantitative comparison (PSNR and SSIM) of THz image restoration
performances on Deer, DNA and Box under different settings. ($\uparrow$:
higher is better)

Method |  PSNR$\uparrow$ |  SSIM$\uparrow$
---|---|---
| Deer | DNA | Box | Deer | DNA | Box
$\texttt{U-Net}_{\texttt{base}}$ [29] | 19.84 | 25.63 | 19.77 | 0.55 | 0.78 | 0.77
Amp-Unet w/o attention | 22.05 | 25.84 | 20.32 | 0.80 | 0.83 | 0.77
Phase-Unet w/o attention | 21.14 | 24.98 | 20.42 | 0.82 | 0.72 | 0.78
Mix-Unet w/o attention | 21.44 | 25.78 | 20.00 | 0.81 | 0.81 | 0.78
Amp-Unet w/ attention | 20.97 | 26.00 | 21.83 | 0.84 | 0.90 | 0.78
Phase-Unet w/ attention | 22.66 | 25.52 | 21.65 | 0.83 | 0.86 | 0.79
SARNet (Ours) | 22.98 | 26.05 | 22.67 | 0.84 | 0.90 | 0.83

###  4.4 Ablation Studies

To verify the effectiveness of multi-spectral feature fusion, we evaluate the
restoration performances with our SARNet under different settings in Table 2.
The compared methods include (1) $\texttt{U-Net}_{\texttt{base}}$ using a
single channel of data (Time-max) without using features of multi-spectral
bands; (2) Amp-Unet w/o attention employing multi-band amplitude feature
(without the attention mechanism) in each of the four spatial-scale branches,
except for the finest scale (that accepts the Time-max image as the input),
where 12 spectral bands of amplitude (3 bands/scale) are fed into the four
spatial-scale branches with the assignment of the highest-frequency band to
the coarsest scale, and vice versa; (3) Phase-Unet w/o attention employing
multi-spectral phase features with the same spectral arrangements as (2), and
without the attention mechanism; (4) Mix-Unet w/o attention concatenating
multi-spectral amplitude and phase features (without the attention mechanism)
in each of the four spatial-scale branches, except for the finest scale (that
accepts the Time-max image as the input), where totally 24 additional spectral
bands of amplitude and phase (3 amplitude plus 3 phase bands for each scale)
are fed into the four branches; (5) Amp-Unet with attention utilizing
attention-guided multi-spectral amplitude features with the same spectral
arrangements as specified in (2); and (6) Phase-Unet with attention utilizing
attention-guided multi-spectral phase features with the same spectral
arrangements as in (2).

The results clearly demonstrate that the proposed SAFM can benefit fusing the
spectral features of both amplitude and phase with different characteristics
for THz image restoration. Specifically, employing additional multi-spectral
features of either amplitude or phase as the input of the multi-scale branches
in the network (i.e., Amp-Unet or Phase-Unet w/o attention) can achieve
performance improvement over $\texttt{U-Net}_{\texttt{base}}$. Combining both
the amplitude and phase features without the proposed subspace-and-attention
guided fusion (i.e., Mix-Unet w/o attention) does not outperform Amp-Unet w/o
attention and usually leads to worse performances. The main reason is that the
characteristics of the amplitude and phase features are too different to be
fused to extract useful features with direct fusion methods. This motivates
our subspace-and-attention guided fusion scheme, that learns to effectively
identify and fuse important and complementary features on a common ground.

Figure 7: Illustration of 3D tomographic reconstruction results on Deer, from
left to right: (a) Time-max, (b) $\texttt{U-Net}_{\texttt{base}}$ [29], (e)
$\texttt{U-Net}_{\texttt{MS}}$, (d) SARNet, and (e) ground-truth.

###  4.5 3D Tomography Reconstruction

Our goal is to reconstruct clear and faithful 3D object shapes through our THz
tomographic imaging system. In our system, the tomography of an object is
reconstructed from 60 views of 2D THz images of the object, each being
restored by SARNet, via the inverse Radon transform. Fig. 7 illustrates the 3D
reconstructions of Deer, showing that Time-max,
$\texttt{U-Net}_{\texttt{base}}$, and $\texttt{U-Net}_{\texttt{MS}}$ tend to
lose important object details such as holes in the deer’s body with Time-max
and the severely distorted antlers and legs with the three methods. In
contrast, our method reconstructs much clearer and more faithful 3D images
with finer details, achieving by far the best 3D THz tomography reconstruction
quality in the literature. Complete 3D reconstruction results are provided in
the supplementary material.

##  5 Conclusion

We proposed a 3D THz imaging system based on fusing multi-spectral features
which is the first to merge THz spatio-spectral data, data-driven models, and
material properties. Based on the physical characteristics of THz waves
passing through different materials, our novel SARNet efficiently fuses
spatio-spectral features with different characteristics on a common ground via
deriving a shared latent subspace and discovering the short/long-range
dependencies between the amplitude and phase to guide the feature fusion for
boosting restoration performance. Our results have confirmed a performance
leap from the relevant state-of-the-art techniques in the area. We believe our
findings in this work will stimulate further applicable research for THz
imaging with advanced computer vision techniques.

## References

  * [1] A. Abbas, M. Abdelsamea, and M.M. Gaber.  Classification of covid-19 in chest x-ray images using detrac deep convolutional neural network.  Appl. Intell., 51(2):854–864, 2021.
  * [2] T. Bowman, T. Chavez, K. Khan, J. Wu, A. Chakraborty, N. Rajaram, K. Bailey, and M. El-Shenawee.  Pulsed terahertz imaging of breast cancer in freshly excised murine tumors.  J. Biomedical optics, 23(2):026004, 2018.
  * [3] E. braham, A. Younus, T.-C. Delagnes, and P. Mounaix.  Non-invasive investigation of art paintings by terahertz imaging.  Applied Physics A, 100(3):585–590, 2010.
  * [4] Y. Calvin, F. Shuting, S. Yiwen, and P.-M. Emma.  The potential of terahertz imaging for cancer diagnosis: A review of investigations to date.  Quantitative imaging in medicine and surgery, 2(1):33, 2012.
  * [5] S. Cheng, Y. Wang, H. Huang, D. Liu, H. Fan, and S. Liu.  NBNet: Noise basis learning for image denoising with subspace projection.  In Proc. IEEE/CVF Int. Conf. Comput. Vis. Pattern Recognit., pages 4896–4906, 2021.
  * [6] A. B. de Gonzalez and S. Darby.  Risk of cancer from diagnostic x-rays: estimates for the uk and 14 other countries.  The Lancet, 363(9406):345–351, 2004.
  * [7] Timothy D. Dorney, Richard G Baraniuk, and Daniel M. Mittleman.  Material parameter estimation with terahertz time-domain spectroscopy.  JOSA A, 18(7):1562–1571, 2001.
  * [8] K. Fukunaga.  Thz technology applied to cultural heritage in practice.  Springer, 2016.
  * [9] P. Geladi, J. Burger, and T. Lestander.  Hyperspectral imaging: calibration problems and solutions.  Chemometrics and intelligent laboratory systems, 72(2):209–217, 2004.
  * [10] K. He, X. Zhang, S. Ren, and J. Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proc. IEEE/CVF Int. Conf. Comput. Vis., pages 1026–1034, 2015\.
  * [11] Y.-C. Hung and S.-H. Yang.  Kernel size characterization for deep learning terahertz tomography.  pages 1–2, 2019.
  * [12] Y.-C. Hung and S.-H. Yang.  Terahertz deep learning computed tomography.  In Proc. Int. Infrad. Milli. THz. Wav., pages 1–2. IEEE, 2019.
  * [13] Christof Janke, Michael Först, Michael Nagel, Heinrich Kurz, and Albrecht Bartels.  Asynchronous optical sampling for high-speed characterization of integrated resonant terahertz sensors.  Optics Lett., 30(11):1405–1407, 2005.
  * [14] C. Jansen, S. Wietzke, O. Peters, M. Scheller, N. Vieweg, M. Salhi, N. Krumbholz, C. Jördens, T. Hochrein, and M. Koch.  Terahertz imaging: applications and perspectives.  Appl. Optics, 49(19):E48–E57, 2010.
  * [15] K. H. Jin, M. T. McCann, E. Froustey, and M. Unser.  Deep convolutional neural network for inverse problems in imaging.  IEEE Trans. Image Process., 26(9):4509–4522, 2017.
  * [16] A. C. Kak.  Algorithms for reconstruction with nondiffracting sources.  Principles of computerized tomographic imaging, pages 49–112, 2001\.
  * [17] M. Kamruzzaman, G. ElMasry, D.-W. Sun, and P. Allen.  Application of nir hyperspectral imaging for discrimination of lamb muscles.  J. Food Engineer., 104(3):332–340, 2011.
  * [18] E. Kang, J. Min, and J. C. Ye.  A deep convolutional neural network using directional wavelets for low-dose x-ray ct reconstruction.  J. Medical physics, 44(10):e360–e375, 2017.
  * [19] K. Kawase, Y. Ogawa, Y. Watanabe, and H. Inoue.  Non-destructive terahertz imaging of illicit drugs using spectral fingerprints.  Optics Express, 11(20):2549–2554, 2003.
  * [20] X. Mao, C. Shen, and Y.-B. Yang.  Image restoration using very deep convolutional encoder-decoder networks with symmetric skip connections.  In Proc. Adv. Neural Inf. Process. Syst., page 2802–2810, 2016\.
  * [21] C. D. Meyer.  Matrix Analysis and Applied Linear Algebra.  SIAM, 2000.
  * [22] D. Mittleman, M. Gupta, R. Neelamani, R. Baraniuk, J. Rudd, and M. Koch.  Recent advances in terahertz imaging.  Applied Physics B, 68(6):1085–1094, 1999.
  * [23] D. M. Mittleman.  Twenty years of terahertz imaging.  Optics Express, 26(8):9417–9431, 2018.
  * [24] A. Ozdemir and K. Polat.  Deep learning applications for hyperspectral imaging: a systematic review.  Journal of the Institute of Electronics and Computer, 2(1):39–56, 2020.
  * [25] D. C. Popescu and A. D. Ellicar.  Point spread function estimation for a terahertz imaging system.  EURASIP J. Adv. Signal Process., 2010(1):575817, 2010.
  * [26] D. C. Popescu, A. Hellicar, and Y. Li.  Phantom-based point spread function estimation for terahertz imaging system.  pages 629–639, 2009.
  * [27] X. Qin, X. Wang, Y. Bai, X. Xie, and H. Jia.  Ffa-net: Feature fusion attention network for single image dehazing.  In Proc. of the AAAI Conf. on Artificial Intelligence, volume 34, pages 11908–11915, 2020.
  * [28] B. Recur, A. Younus, S. Salort, P. Mounaix, B. Chassagne, P. Desbarats, J. Caumes, and E. Abraham.  Investigation on reconstruction methods applied to 3d terahertz computed tomography.  Optics Express, 19(6):5105–5117, 2011.
  * [29] O. Ronneberger, P. Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Proc. Int. Conf. Medical Image Comput. Comput.-Assisted Intervention, pages 234–241, 2015.
  * [30] H. H. Rotermund, W. Engel, S. Jakubith, A. Von Oertzen, and G. Ertl.  Methods and application of uv photoelectron microscopy in heterogenous catalysis.  Ultramicroscopy, 36(1-3):164–172, 1991.
  * [31] A. R. Round, S. J. Wilkinson, C. J. Hall, K. D. Rogers, O. Glatter, T. Wess, and I. O. Ellis.  A preliminary study of breast cancer diagnosis using laboratory based small angle x-ray scattering.  Physics in Medicine & Biology, 50(17):4159, 2005.
  * [32] D. Saeedkia.  Handbook of terahertz technology for imaging, sensing and communications.  Elsevier, 2013.
  * [33] R. Schultz, T. Nielsen, R. Jeff Zavaleta, R. Wyatt, and H. Garner.  Hyperspectral imaging: a novel approach for microscopic analysis.  Cytometry, 43(4):239–247, 2001.
  * [34] D. M. Slocum, E. J. Slingerland, R .H. Giles, and T. M. Goyette.  Atmospheric absorption of terahertz radiation and water vapor continuum effects.  J. Quantitative Spectroscopy and Radiative Transfer, 127:49–63, 2013.
  * [35] T. M. Tuan, H. Fujita, N. Dey, A. S. Ashour, T. N. Ngoc, D.-T. Chu, et al.  Dental diagnosis from x-ray images: an expert system based on fuzzy computing.  Biomed. Signal Process. Control, 39:64–73, 2018.
  * [36] M. Van Exter, C. Fattinger, and D. Grischkowsky.  Terahertz time-domain spectroscopy of water vapor.  Optics Lett., 14(20):1128–1130, 1989.
  * [37] T. M. Wong, M. Kahl, P. H. Bolívar, and A. Kolb.  Computational image enhancement for frequency modulated continuous wave (fmcw) thz image.  J. Infrared, Millimeter, and Terahertz Waves, 40(7):775–800, 2019\.
  * [38] Xianghua Xie.  A review of recent advances in surface defect detection using texture analysis techniques.  ELCVIA: Electron. Lett. Comput. Vis. Iimage Ana., pages 1–22, 2008\.
  * [39] L. Yujiri, M. Shoucri, and P. Moffa.  Passive millimeter wave imaging.  IEEE Microwave Mag., 4(3):39–50, 2003.
  * [40] H. Zhang, I. Goodfellow, D. Metaxas, and A. Odena.  Self-attention generative adversarial networks.  In Proc. Int. Conf. Mach. learn., pages 7354–7363, 2019.
  * [41] K. Zhang, W. Zuo ana Y. Chen, D. Meng, and L. Zhang.  Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising.  IEEE Trans. Image Process., 26(7):3142–3155, 2017.
  * [42] K. Zhang, W. M. Zuo, and L. Zhang.  FFDNet: Toward a fast and flexible solution for cnn-based image denoising.  IEEE Trans. Image Process., 27(9):4608–4622, 2018.
  * [43] Y. Zhang, Y. Tian, Y. Kong, B. Zhong, and Y. Fu.  Residual dense network for image restoration.  IEEE Trans. Pattern Anal. Mach. Intell., 2020.
  * [44] B. Zhu, J. Z. Liu, S. F. Cauley, R. Bruce Rosen, and M. S.Rosen.  Image reconstruction by domain-transform manifold learning.  Nature, 555(7697):487–492, 2018.


11institutetext: Department of Mechanical and Mechatronics Engineering,
University of Waterloo, Waterloo, Canada
200 University Ave. W., Waterloo, Ontario, N2L 3G1

# Uncertainty quantification of tabulated supercritical thermodynamics for
compressible Navier-Stokes solvers

Sai Praneeth<EMAIL_ADDRESS>
###### Abstract

Non-ideal state equations are needed to compute a growing number of
engineering-relevant problems. The additional computational overhead from the
complex thermodynamics accounts for a significant portion of the total
computation, especially the near-critical or transcritical thermodynamic
regimes. A compromise between computational speed and the accuracy of the
thermodynamic property evaluations results in a propagation of the error from
the thermodynamics to the hydrodynamic computations. This work proposes a
systematic error quantification and computational cost estimate of the various
approaches for equation of state computation for use in compressible Navier-
Stokes solvers in the supercritical regime. We develop a parallelized, high-
order, finite volume solver with highly-modular thermodynamic implementation
to compute the compressible equations in conservative form. Three tabular
approaches are investigated: homogeneous tabulation, block structured adaptive
mesh refinement tabulation, and a n-dimensional Bezier surface patch on an
adaptive structured mesh. We define a set of standardized error metrics and
evaluate the thermodynamic error, table size and computational expense for
each approach. We also present an uncertainty quantification methodology for
tabular equation of state.

##  1 Introduction

The simulation of thermodynamically complex systems present new challenges to
the computational fluid dynamics community. Recent progress has been made in
predictive modelling of trans- and supercritical flows to address design
challenges in supercritical water-cooled reactors (SCWR) [2], understanding
thermoacoustic oscillation in transcritical heat exchangers, efficient
supercritical extraction [1], cryogenic mixing and combustion in liquid rocket
engines [6], and technologies relating to carbon capturing [12]. These flows
are characterized by strong thermodynamic non-linearities that arise,
primarily, in the near-critical region of the fluid and require special
considerations for efficient computational modelling. In its essence, the
Navier-Stokes equations are a set of partial differential equations that
govern the transport of the conservative properties of dynamic fluid system.
Typically, the governing equations are written in terms of density, momentum
and total energy. These governing equations must be closed using an equation
of state (EOS) that is needed to compute the contribution due to local
pressure and thermal gradients in the system. The EOS is also needed to obtain
an estimate of the transport properties (e.g. viscosity, thermal
conductivity), which depend on the local thermodynamic state.

For a single phase, single species fluid, two independent thermodynamic
variables are needed to fully define the local thermodynamic state. The
compressible Navier-Stokes equation fully define the thermodynamic state of
the system as both the density ($\rho$) and internal energy ($e$) can be
computed directly from the conservative variable matrix:
$\left[\rho,\rho\vec{u},\rho E\right]$, where $\rho E$ typically, represents
the total energy (internal plus kinetic). Most CFD solvers require the state
equation to transform the variables from conservative to their primitive
state: $\left[\rho,\vec{u},p,T\right]$. Under a perfect gas assumption, we
trivially find $p=\rho(\gamma-1)e$ which gives us, using the EOS, the solution
of $T=p/(\rho R)$; under an ideal gas assumption the complexity rises slightly
with the non-constant definition of the specific heat. For real fluids, there
is no explicit relation between $(T,p)$ and $(\rho,e)$ from classical
analytical state equations. This typically requires an iterative algorithm to
find the correct $(T,p)$ tuple for the known $(\rho,e)$ from the conservation
equations. Given the complexity of the EOS in the supercritical regime—its
simplest form being a highly non-linear cubic equation—this thermodynamic
computation demands non-negligible computational resources. Furthermore, due
to the inherent thermodynamic non-linearities about pseudo-boiling point,
small errors in the thermodynamic property estimation can lead to large
variation in the computed pressure and temperature and concomitantly directly
affect the thermodynamic gradient calculations in the flow. This, in turn, may
negatively affect in the computed integral properties of the simulation.

The added computational expense required for complex state equations quickly
becomes a limiting factor in large scale computations. Typical real-fluid
simulations with cubic state equations will show slow downs of 50% to 100%
compared to ideal gas simulations depending on the thermodynamic state of the
simulation. Masquelet [9] notes that over 20% of the computational expenses in
the real fluid calculations comes from expensive numerical operations such as
square-roots, logarithms and iterative algorithms. The overall computational
expense rises drastically with more complex analytical state equation
evaluations (e.g. Modified Benedict-Webb-Rubin). For this reason, cubic class
of state equations (e.g. Peng-Robinson-Stryjek-Vera, Soave-Redlich-Kwong etc.)
have become the standard when using analytical thermodynamic definitions for
CFD, in spite of their sub-optimal representation of the actual
thermodynamics.

The predictive capability of fluid solvers is heavily dependent on the
accuracy of the EOS in the transcritical regime. Given the limitations of a
finite-degree polynomial-type function, typical cubic-state equations show a
very significant errors compared to experimental values; errors of over 30%
are expected for the Peng-Robinson equation in the some parts of the state
space [11]. Furthermore, within the subcritical, two-phase, the continuous
state equation has an inversion, meaning that the isobaric lines do not behave
monotonically which can lead to non-physical behavior and numerical
instabilities.

Therefore, a compromise needs to be struck between thermodynamic accuracy and
computational efficiency. To more closely control the thermodynamic error and,
in some cases, accelerate the computational cost related to the EOS
computations, a growing number of codes use pre-tabulated state equations
based on well-defined REFPROP quantities. The tabulation accuracy is generally
proportional to the size of the table. On the downside, the size of the table
may hinder efficient caching and prefetching for memory which has a
significant cost on the computational speed of a simulation. It is also
particularly problematic for large scale, memory-limited computation on high-
performance computers. There is surprising lack of quantitative information
relating the thermodynamic accuracy and computational expense specifically
related to CFD solvers. Most of the comparisons discuss the errors related to
state equations, and do not tackle the specificity of the thermodynamic
calculations of CFD solvers which require solutions for known density and
internal energy. There are some works that address novel ways to tabulate the
state equations for CFD solvers [4, 7, 8] and other works that evaluate the
uncertainty propagation of thermodynamic error [3]. The present work
quantifies different state equation implementations and proposes quantifiable
metrics to evaluate the accuracy and computational expenses of each approach.
The focus of this work lies in the trans- and supercritical regime in which
there is a first-order continuity in the state properties. Finally, this work
will present the resulting implications of thermodynamic error on simulation
results of well defined CFD cases at transcritical thermodynamic conditions.
The following section describes the main details of our numerical framework.
Afterward, the various tabulated approaches used in the current work are
presented followed by a set of metrics to evaluate the thermodynamic accuracy
and computational cost for all tabulated approaches. Finally, statistical
study is undertaken to quantify the propagation of uncertainties from the
thermo- to the hydronamic quantities for well defined numerical simulations
(this last point will be completed in the final paper).

##  2 Numerical Framework

A modular finite-volume code, named ThermoCode, was developed for the study of
complex thermodynamic fluid systems. The code, which solves the compressible
Navier-Stokes equations in fully conservative form, rests on a high-order flux
reconstruction at the faces which requires a structured grid arrangement. The
code is formally of fourth-order accuracy. The flux at each face is computed
as the sum of a left- and a right biased stencil at the face both of second-
order accuracy. The code has been parallelized with efficient flux passing
between the processors which allows for a limited number of ghost cells and
communication overhead. Explicit time advancement is done with low-storage,
fourth-order Runge-Kutta scheme. The purpose of a high fidelity code is to
remove, as much as possible, numerical dissipation from the thermodynamic
error evaluation. The coupling between the Navier-Stokes equations and the
thermodynamics and transport property modules are designed to be extremely
modular. The code contains a number of different analytical cubic state
equations (Peng-Robinson, Suave-Redlich-Kwong) as well as a number of tabular
approaches discussed in this paper.

##  3 Tabulated State Equations

State equation tabulation involves storing the thermodynamic and transport
properties of the species in a library (e.g. binary file, .xml file), by doing
so we move the evaluation of the EOS to a preprocessing step. Thus decoupling
the EOS computation complexity from runtime performance. This approach allows
for the use of accurate high-quality EOS or REFPROP reference data without
major runtime penalty. Tabulated approaches also provide an effective control
on the thermodynamic error of the state equation. For single species fluid,
two independent thermodynamic variables are needed to fully define the local
thermodynamic state in a single phase. Classically, the choice of the
thermodynamic tuple has been pressure and temperature given the convenience of
computing the isothermal and isobaric conditions. For CFD simulations, this
approach is sub-optimal since the known transported quantities are $\rho$ and
$e$; thus making these the ideal basis for any thermodynamic table. Once the
corresponding position in the table is found, all thermodynamic and transport
properties may be interpolated or computed.

Figure 1: A typical cell in a table, the red dot indicating the unknown point

In this work, we consider three different tabulation approaches:

  * •

Homogeneous of logarithmic tabulation

  * •

Block structured Adaptive Mesh Refinement (AMR)

  * •

Adaptive Mesh Refinement with Bezier patch fitting

Details of each tabulation approach is presented in the following:

###  3.1 Homogeneous or logarithmic tabulation

The simplest data structure for tabular reconstruction is one that spans the
region of thermodynamic interest with equally spaced grids, or alternatively
under a known mapping function such as a logarithmic distribution. Equal
spacing provides very efficient table look up because the reconstruction
interval can be computed directly without searching, but this results in very
large table sizes as the grid size is not adapted to the magnitude of the
local error. The level of accuracy is set by the region in which the function
varies most rapidly so that slowly varying regions must be over-refined to
accommodate high-gradient regions. Tabulating in density/energy space, as
shown in figure 2, allows for a simple lookup from the transport equations but
simultaneously the corner regions of the table contain pressure and
temperature values that are far outside the region of interest for most CFD
problems. Therefore, bounds typical isobaric conditions of the simulation are
mapped to the energy/density space.

Figure 2: Uniformly spaced lookup table for oxygen

###  3.2 Block Structured Adaptive Mesh Refinement (AMR)

The adaptive mesh refinement table is based on Cartesian adaptive mesh
structure. We begin by defining a rectangular region in the x-y plane (P-T
plane for this example) that includes the physical domain of interest. The
required properties (internal energy, density, enthalpy, entropy, etc.) are
evaluated at each of the four corners of the rectangle and appropriate
reconstruction is used to evaluate these properties at pre-selected points in
the rectangle. For simplicity, we chose bivariate interpolation to evaluate
the properties at those preselected 25 equidistant points. The reconstructed
(interpolated) properties are compared with the ’exact’ values obtained from
the NIST database to assess the reconstruction (interpolation) error. If the
error is more than the user-specified tolerance, it is subdivided into four
equal rectangles with property evaluations again being computed and compared
with the exact NIST database values. As the interpolation grid is refined,
relatively flatter(linear) portions of the surface will fall within the error
tolerance and will not require further refinement while rectangles in the
regions of larger gradients will be refined further. Keeping track of these
adaptively refined regions in a tree-based structure enables rapid search
procedures that can be comparable to the speed of equally of spaced tables.

Figure 3: Block structured AMR lookup table for oxygen with 0.01% accuracy

####  3.2.1 Lookup Approach

The Structured block based adaptive mesh refinement has been widely used in
problems with highly oscillatory phenomenon, as it provides uniform
interpolation accuracy over the entire thermodynamic region of interest. The
patch (in this case, the rectangular cell) must be located before bilinear
interpolation is used to reconstruct the thermodynamic property. One easy
straightforward method is to search the patches one by one until we zero-in
onto the last child node (patch) that contains the given unknown coordinates
(Input thermodynamic parameters T,P or $\rho$, E). This sequential searching
algorithm’s efficiency strongly depends on the table size thereby inhibits
speed if the thermodynamic domain of interest is large and the depth of
refinement is high.

Figure 4: Sample quadtree subdivision

As an alternative approach, we employ a fast lookup algorithm that requires
additional memory to store the indices in tree data-structure, and the
efficiency doesn’t depend much on the depth of the tree. It is constructed by
the following steps,

A uniform index array based on quadtree subdivision covering the entire
computational domain is generated as per the maximum refinement level of the
generated adaptive table. For example, at a given maximum refinement level of
3, the uniform index table is shown in the figure 4. If the maximum refinement
level is n+1 (the refinement level of the single rectangular cell that cover
the entire thermodynamic region of interest is 0) then there will be $4^{n}$
number of individual patches (in this case rectangular cells), the integer
number stored in the cell indicates the index of that patch.

An overlapping mapped index array is also generated which contains the child
node index in the place of all further sub divisions pertaining to that child
if there was a refinement as illustrated in figure 5. The cell within which a
point $(T_{x},P_{x})$ is located can be determined at an algorithmic
complexity $O(1)$, directly without searching the tree (which is of algortimic
complexity $O(n)$ where n being the refinement level) , by using the mapping
index array. Furthermore, the storage required for the uniform index table is
small for this generation computers, for example, if the maximum refinement
level of the adaptive table is 10, and an integer requires 4 B of memory, the
uniform index table requires $(2^{9}x2^{9}x4)/(1024^{2})MB=1MB$.

84 81 72 69 36 33 24 21 83 82 71 70 35 34 23 22 80 77 76 73 32 29 28 25 79 78
75 74 31 30 27 26 68 65 56 53 52 49 40 37 67 66 55 54 51 50 39 38 64 61 60 57
48 45 44 41 63 62 59 58 47 46 43 42 Figure 5: Uniform index for a maximum
refinement level of 3 Figure 6: Sample Overlapping index for the uniform index
at the maximum refinement level = 3

###  3.3 Bezier patch with AMR

In the block structured AMR table, the underlying criteria for refinement is
based on the bivariate approximation. Rather than tabulating solely on the
basis of bilinear reconstruction error we chose to implement Bezier patches
that are approximated surfaces. The approximated function of two independent
variables is given by [4]

$F(u,v)={{\sum}}_{i=0}^{3}{{\sum}}_{j=0}^{3}B_{i}^{3}(u)B_{j}^{3}(v)b_{ij}$
(1)

Figure 7: Bezier patch based Tabulated EoS for Density($\rho$) adaptively
refined using Quadtree

Here, $B_{i}^{3}(u)$ are the third-degree (fourth-order) Bernstein
polynomials, $b_{ij}$ are the control points.

$B_{i}^{3}(u)=\frac{3!}{i!(3-i)!}(1-u)^{3-i}u^{i}$ (2)

Using the linear combination of these $b_{ij}$ control points, the shape of
the surface is determined. The surface is generated by recursively subdividing
the domain until the resulting Bezier patches are capable of reconstructing
the original value to a specified error threshold or if the maximum refinement
level is reached.

##  4 Uncertainty quantification of Tabular equation of State

An open problem in tabular equation of states is uncertainty visualization and
quantification, how uncertainty and interpolation interact. There are many
established interpolation schemes available, but how reliable are the
interpolated values ? We want to address this problem by analyzing the
interpolation in the tabular equation of state and its impact on gaussian
distribution data.

An important feature of the tabular representation of Equation of state is the
feasibility to represent the underlying uncertainty associated within the
thermodynamic properties. This helps in better understanding of the
uncertainty propagation through the CFD code.

We followed this framework for parametric uncertainty quantification in
analytical equation of state.

  * •

Identify and characterize all sources of input uncertainty, aleatory and
epistemic.

  * •

Propagate input uncertainties through the computation.

  * •

Estimate uncertainty due to numerical approximation.

  * •

Determine total uncertainty in the predicted quality of uncertainty.

Inorder to quantify the uncertainty, we should begin with the underlying
analytical form of the equation of state with which the points are tabulated.
analytical equation as in a mathematical form of equation of state is an
approximate representation of an element’s complex thermodynamic system. Such
mathematical forms (models) use a large number of model parameters ,whose
exact values can’t be determined analytically most of the times, resulting in
the fitting of model predictions with calibration data obtained from
experiments or much more highly defined models. The author has chosen Peng-
Robinson equation of state to demonstrate the UQ, as the highly accurate
Benedict-Webb-Rubin(BWR) equation has more than 11 parametric constraints it
would be intensive.

$\displaystyle P$ $\displaystyle=$
$\displaystyle\frac{RT}{V_{m}-b}-\frac{a\alpha}{{V_{m}}^{2}+2bV_{m}-b^{2}}$
(3) $\displaystyle a$ $\displaystyle=$
$\displaystyle\frac{0.45724R^{2}{T_{c}}^{2}}{p_{c}}$ $\displaystyle b$
$\displaystyle=$ $\displaystyle\frac{0.07780RT_{c}}{p_{c}}$
$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle(1+\kappa(1-{\frac{T}{T_{c}}}^{0.5}))^{2}$ $\displaystyle\kappa$
$\displaystyle=$ $\displaystyle 0.37464+1.54226\omega-0.26992\omega^{2}$

For the Peng Robinson equation of state the uncertain model parameters are
obtained from [10]

Table 1: Parametric Uncertainty in Peng Robinson Equation of state

Parameter | Uncertainty(%)
---|---
Critical Temperature* [T] |  $\scriptstyle\pm$3
Critical Pressure* [p] |  $\scriptstyle\pm$5
Accentric factor [$\omega$] |  $\scriptstyle\pm$5
$k^{\\#}$ |  $\scriptstyle\pm$5

  * •

*Experimental, #-Assumed

Figure 8: The contour plots show the two-dimensional distributions of all
pairs of parameters, where the color scale represents the likelihood value
normalized to the maximum likelihood. The levels corresponding to $1-\sigma$
and $2-\sigma$ regions are indicated in the color bar. Figure 9: Density
computed with UQ enabled Peng Robinson Equation of state for Oxygen

In this paper, one of the objective is to propagate this model parameter
uncertainty through the flux calculations for the solver to understand the
statistical behaviour of output properties.

###  4.1 UQ in Tabular equation of state

We analyzed the effect of bilinear interpolation of a data with four
uncorrelated scalar Gaussian distributions $X_{i}\sim
N(\mu_{i},{\sigma_{i}}^{2})$ at the positions $s_{i}$, where i = 1,2,3,4.

Figure 10: Color coded variance as the result of Bilinear interpolation of
uncorrelated Gaussian distributed variables given on the grid points

By our understanding, bilinear interpolation is in general not well suited to
estimate uncertain values. A key element of interpolation is that the
interpolation result equals the value at the data points. In the case that the
value is uncertain, this property may be unwanted because of that data point
is disregarded. And the result of the variance might be that the best
estimator for that data point is actually not the give value at that data
point.

We point out that the traditional bilinear interpolation, which is widely used
in tabular equation of state may yield very counter-intuitive results when
dealing with uncertain data. If the covariance between neighbouring points is
weaker than the variance at these points, then the resulting interpolated data
is most certain between the grid points (where we have no data) and the
variance has its maxima at the grid points (where the data is given). If one
encounters such data, it should be obvious by now that one needs other
methods.

###  4.2 Thermodynamic consistency

Thermodynamic consistency in property evaluations has been discussed by
Swesty[13], using the Helmholtz free energy as a basis function. The problem
of thermodynamic consistency is very important in many cases where reactive
flows need to be simulated over very long timescales. Thermodynamic
inconsistency may become manifest via the unphysical buildup of entropy, or
temperature.

Thermodynamic consistency requires that the first law of thermodynamics be
satisfied, so by referring to how well the EOS satisfies the constraints posed
by Maxwell relations we can quantify the thermodynamic consistency[5].

If we have internal energy E expressed as a function of specific volume V and
entropy S. The thermodynamic definitions of pressure P and temperature T are

$P=-\frac{\partial E}{\partial V}\Bigg{|}_{S},T=-\frac{\partial E}{\partial
S}\Bigg{|}_{V}$ (4)

With the bilinear interpolation, both the density and its first order
derivatives vary linearly over each local rectangle. Further, enthalpy and
density have no interdependence on each other. Nevertheless, as the
thermodynamic properties stored at the mesh nodes are themselves taken from a
consistent thermodynamic database (in this case NIST), the internal
inconsistencies will be small on a fine grid.

###  4.3 Error estimates and Verification

The T-P tabular grid is sampled into 1000x1000 log-linear grid. In each patch
100 points are sampled and then the maximum relative error in density is
computed. As seen from the relative error contour plot, the maximum error is
along the phase change line or where the density changes rapidly. Except for
that region the uniform accuracy is maintained as a result of adaptive
subdivision. The highest error points lies close to the critical point. The
number of patches required to capture the thermodynamics is different in three
approaches (Homogeneous, Adaptive, Bezier patch) The highest being the
homogeneous due to its uniform spacing, it will take enormous number of
patches to get the error close to 0.01%, althought Bezier patch has less
number of patches it can only accommodate one property at a time, which means
it requires individual table for each property. The Block based AMR has
moderate number of patches that are not oversized but can tabulate all
thermodynamic properties.

Figure 11: Relative error in density for the Block based adaptively refined table Table 2: Error estimates Grid type |  | Block based AMR
---
(0.0001 tol.)
| Bezier Patch
---
AMR
(0.0001 tol.)
| Total
---
patches
12648 | 1120
| % patches
---
with rho
error >0.001
0.101 | -
| Maximum
---
rho error $\%$
1.2 | -
Figure 12: Error plot for speed of sound for Oxygen computed with different
equation of states

The comparison of various tabulation approaches demands a quantitative
thermodynamic error metric. As the exactness thermodynamic properties remain a
subject of scientific relevance, the current study evaluates the error
relative to the NIST database. These data points represent the highest
fidelity data set for the computation of the thermodynamic and thermophysical
properties and constructed from an ensemble of experimental data supplemented
with Modified Benedict-Webb-Rubin curve fitting.

To evaluate different types of equation of states, statistical parameters are
used, The Average Percent Relative Error (ARE%), which is defined below is a
measure of the bias of the correlation; a value of zero indicates a random of
the measured values around the correlation.

$\displaystyle
ARE\%=\frac{100}{N_{d}}{{\sum}}_{i=1}^{N_{d}}\frac{(\rho_{i}^{exp.}-\rho_{i}^{calc.})}{\rho_{i}^{exp.}}$
(5)

The AARE%, is the arithmetic average of the absolute values of the relative
errors; is an indication of the accuracy of the correlation.

$\displaystyle
AARE\%=\frac{100}{N_{d}}{{\sum}}_{i=1}^{N_{d}}\frac{|\rho_{i}^{exp.}-\rho_{i}^{calc.}|}{\rho_{i}^{exp.}}$
(6)

The $R^{2}$, is the correlation coefficient; is a measure of the precision of
fit of the data. If data are perfectly correlated, then $R^{2}=1$. A small
value of AARE% and $R^{2}$ close value to one (simultaneously) denote a good
correlation based on good data.

Another parameter, Sum of Absolute of Residual (SAR) shows the reliability of
correlation for higher order data points.

$\displaystyle SAR\%={{\sum}}_{i=1}^{N_{d}}|\rho_{i}^{exp.}-\rho_{i}^{calc.}|$
(7)

The statistical parameters of the well known EoSs are listed in table 1.

Table 3: Error metrics comparison with NIST | Equation
---
of State
AARE | ARE |  | SAR
---
[Kg/m^3]
R^2
| Peng-
---
Robinson
3.015 | 0.993 | 70.600 | 0.991
| Redlich-
---
Kwong
4.409 | -0.026 | 89.816 | 0.987
| Adap-
---
Table
1.022 | 0.0021 | 91.232 | 0.995
| Soave-
---
Redlich-
Kwong
6.388 | 6.381 | 131.009 | 0.990

##  5 Computational cost estimates

The increase in speed of evaluating properties compared to cubic EOS is one of
the primary motivations. Here we show an example of the speedup of
computations. When pressure is evaluated from randomly generated $\rho$, $e$
combination in the supercritical regime. For each EOS, 100000 loops/calls are
executed and the slowest run time is tabulated.

Searching algorithm and interpolation scheme influence the efficiency of the
table tabular equation of state, the adaptive mesh is comparably denser in the
near critical point regime, because of high non linearity which resulted in
comparatively higher lookup search times for those regions as shown in the
figure 13

Figure 13: Block based tabulated state equation lookup Computational Cost in
$\mu$ sec  Figure 14: Travelling Acoustic wave, one flow-through time period
simulation computational cost comparison with different equation of states

##  6 One dimensional test case:

A one dimensional test case is identified to evaluate the performance of the
proposed tabulated equation of state.

## Harmonic Acoustic wave

: In one-dimensional acoustic wave propagation in supercritical fluid case,
Using periodic boundary conditions at both sides a harmonic wave is
initialized. The computational domain is $x\in{[0,10]}$ m and 100 grid points
are used, giving a uniform grid spacing $\Delta x=0.1m$. The sound pressure
level $L_{p}=23$ (dB) is calculated according to:

$L_{p}=20log_{10}(\frac{\Delta p}{p_{ref}\sqrt{2}})$ (8)

Where $\Delta p$ is the amplitude of the pressure harmonic wave (Pa).

As shown in the figure, there is a significant effect in the form of density
by using most accurate tabular real fluid thermodynamics (AdapTable) rather
than conventional Peng robinson equation of state (PREOS)

Figure 15: Harmonic Acoustic wave propagation

## Conclusion

We have demonstrated three different types of tabulated equation of states for
an efficient calculation of thermodynamic properties of real and supercritical
fluids. The EOS look-up tables are constructed using Block based AMR, Bezier
patch based AMR, homogeneous tabulation using bilinear interpolation. Their
error metrics, computational cost and uncertainty in such tabular
representation of gaussian distributed data is discussed along with a test
case highlighting the use of a tabular real gas equation of state in a CFD
code.

## References

  * [1] O. Aizpurua-Olaizola, M. Ormazabal, A. Vallejo, M. Olivares, P. Navarro, N. Etxebarria, and A. Usobiaga.  Optimization of supercritical fluid consecutive extractions of fatty acids and polyphenols from vitis vinifera grape wastes.  Journal of Food Science, 80(1):E101–E107, 2015.
  * [2] J. Cai, C. Renault, and J. Gou.  Supercritical water-cooled reactors.  Science and Technology of Nuclear Installations, 2014.
  * [3] P. Cinnella, P. M. Congedo, V. Pediroda, and L. Parussini.  Sensitivity analysis of dense gas flow simulations to thermodynamic uncertainties.  Phys. Fluids, 2011.
  * [4] E. Collins and E. Luke.  Fast evaluation of complex equations of state.  Electronic Journal of Differential Equations, 2013.
  * [5] G. A. Dilts.  Consistent thermodynamic derivative estimates for tabular equations of state.  Phys. Rev. E, 73:066704, Jun 2006.
  * [6] J.-P. Hickey and M. Ihme.  Center for Turbulence Research - Annual Research Briefs 2014, chapter Supercritical mixing and combustion in rocket propulsion, pages 21–36.  Stanford University, 2013.
  * [7] Z. Liu, J. Liang, and Y. Pan.  Construction of thermodynamic properties look-up table with block-structured adaptive mesh refinement method.  J. Thermophysics. Heat Transfer, 2014.
  * [8] Z. Liu, J. Liang, and Y. Pan.  Efficient thermodynamic properties reconstruction method with adaptive triangular mesh.  J Thermophys Heat Transfer, 2015.
  * [9] M. M. Masquelet.  Simulations of a sub-scale liquid rocket engine: Transient heat transfer in a real gas environment.  PhD thesis, Georgia Institute of Technology, 2006.
  * [10] R. L. Rowley, W. V. Wilding, J. L. Oscarson, and Y. Yang.  Database tools for evaluating thermophysical property data.  International Journal of Thermophysics, 28(3):805–823, 2007.
  * [11] R. Stryjek and J. Vera.  Prsv: An improved peng-robinson equation of state for pure compounds and mixtures.  The canadian journal of chemical engineering, 64(2):323–333, 1986\.
  * [12] S. Stucki, F. Vogel, C. Ludwig, A. G. Haiduc, and M. Brandenberger.  Catalytic gasification of algae in supercritical water for biofuel production and carbon capture.  Energy & Environmental Science, 2(5):535–541, 2009.
  * [13] F. D. Swesty.  Thermodynamically consistent interpolation for equation of state tables.  Journal of Computational Physics, 127(1):118–127, 1996.