# Robot Localization and Mapping
Final Report
Sequential Adversarial Learning for
Self-Supervised Deep Visual Odometry

Akankshya Kar, Sajal Maheshwari, Shamit Lal, Vinay Sameer Raja Kadi
{akankshk,sajalm,shamitl vkadi} @andrew.cmu.edu

##  I INTRODUCTION

††Our Pytorch implementation: BitBucket repo

Visual odometry (VO) and SLAM have been using multi-view geometry via local
structure from motion for decades. These methods have a slight disadvantage in
challenging scenarios such as low-texture images, dynamic scenarios, etc.
Meanwhile, use of deep neural networks to extract high level features is
ubiquitous in computer vision. For VO, we can use these deep networks to
extract depth and pose estimates using these high level features. The visual
odometry task then can be modeled as an image generation task where the pose
estimation is the by-product. This can also be achieved in a self-supervised
manner, thereby eliminating the data (supervised) intensive nature of training
deep neural networks. Although some works tried the similar approach [1], the
depth and pose estimation in the previous works are vague sometimes resulting
in accumulation of error (drift) along the trajectory.

The goal of this work is to tackle these limitations of past approaches and to
develop a method that can provide better depths and pose estimates. To address
this, a couple of approaches are explored:

  1. 1.

Modeling: Using optical flow and recurrent neural networks (RNN) in order to
exploit spatio-temporal correlations which can provide more information to
estimate depth.

  2. 2.

Loss function: Generative adversarial network (GAN) [2] is deployed to improve
the depth estimation (and thereby pose too), as shown in Figure 1. This
additional loss term improves the realism in generated images and reduces
artifacts.

The goal of generator is to generate an image by view synthesis which looks
exactly like the target image and can be used to fool the discriminator. On
the other hand, discriminator in GAN evaluates the generated warped image and
tries to distinguish it from the real target image. The eventual goal is to
train the generator such that it is hard for the discriminator to distinguish.
During inference, only the depth and pose estimators are used to compute the
trajectory for VO task. This project would be an implementation of [3]. It is
built on the Pytorch implementation of SFMlearner [4].

However, using GAN to generate the images does not take into account the
geometric constraints of the scene. Most recent works based on learning-based
techniques have not provided a comparison of their results with those
generated using the traditional geometric analysis. We propose to explore this
aspect of GAN-based synthesis models and analyze their possible
advantages/disadvantages over the traditional methods. We have released our
code here: code repo.

Figure 1: Entire pipeline of the proposed method [3]

To isolate the performance of each component of our model and validate its
contribution to the final goal of high quality depth and pose predictions, we
also perform a number of ablation experiments. Unlike [3], we explore
different GAN models and algorithms like Wasserstein generative adversarial
networks [5] and PatchGANs [6].

The structure of this report is as follows: In Section II, we will cover the
existing literature related to our work. Section III delineates the
methodology, different components and the loss functions used. This will be
followed by the experiments in Section IV. We discuss the challenges faced in
Section V. We finally wrap up the report with Conclusion and Future Work in
Section VI.

##  II Literature Review

###  II-A Pose and Depth estimation

Pose estimation and visual odometry have traditionally been viewed as a multi-
view geometry problems. Initial approaches to solve this problem used graph
based optimizations [7, 8]. These methods typically relied on handcrafted
features for finding corresponding points between images and then using these
points for subsequent pose estimation [9].

Recently, deep learning based methods have taken over the pose and depth
estimation problem as they tend to perform significantly better than the
metric-based approaches. Some of the recent works using deep learning for pose
estimation are [10, 1, 11, 3, 12]. The model proposed by Ummenhofer et al.in
[12] tries to jointly estimate the pose along with the depth in an end-to-end
manner. DeepTAM [10], on the other hand, uses two different networks for pose
estimation and depth estimation. However, these methods require supervision
for training. The models, therefore, can only be trained on synthetic data as
supervised data generation for pose estimation and depth estimation in the
wild is both expensive and heavily prone to inaccuracies.

In order to overcome this problem, there have been numerous self-supervised
methods that have recently been proposed. These methods include SFMLearner by
Zhou et al. [1]. Since the current work that we present in the report [3] is
built over this method, we explain the SFMLearner model in detail.

The SFMLearner model utilizes the correlation between depth and pose to learn
both entities jointly. The model makes use of a photometric and geometric loss
to model the constraints between the source and the target pixel in an image
as a function of its depth, pose change and the camera intrinsic. Self-
supervision is used to train these models by using novel view synthesis. View
synthesis is predicting what a scene will look like from a different view.
This is done by unprojecting the image in source view using camera intrinsics
and depth image, warping it to the pose of the target view using the source to
target relative pose transformation, and finally projecting it to the pixel
space. Differential modules are used for view synthesis task in order to take
advantage of gradient based training [13, 14, 15, 16, 17, 18]. Authors show
the performance of their method on KITTI [19] and Make3D datasets [20]. The
results show that the single view depth prediction gives results equivalent to
the supervised learning approaches on KITTI dataset. Ablation studies showed
that explainability network only offered a small boost in performance,
probably because most scenes in KITTI are static. The results on Make3D
dataset shows the generalization ability of the depth prediction model. For
pose estimation network, authors use the KITTI odometry split, and show that
their network performs comparable to ORB-SLAM.

This work was further extended in [3], by incorporating sequential learning as
part of the VO framework, estimating depth conditioned on multiple previous
frames, in contrast to the single view depth prediction. This work is
explained in more detail in Section III.

###  II-B Generative Adversarial learning

A majority of the self-supervised approaches methods mentioned in Section II-A
need to generate photo-realistic images as part of the target image. This is
done using Generative Adversarial Networks or GANs as they are commonly
referred to. These were first proposed by Goodfellow et al. in [2] to try to
generate realistic images from a latent space by trying to model the loss as a
mini-max function between the two components (Generator and Discriminator) of
the architecture. GANs have been a topic of extensive study for the past few
years and many improvements and variants over the original GAN architecture
have been suggested [5, 21, 22]. GANs have also been used with great success
in a number of image editing applications [23, 24, 25, 26, 27, 28, 29, 30,
31]. In this work, we have used two variants as part of the discriminator in
the adversarial training. These are the PatchGAN architecture proposed in [21]
and the WGAN method in [5]. The details for the same are provided in Section
III.

##  III Methodology

Our network used for estimation of depth and pose in an self-supervised manner
consists of the following 4 components:

###  III-A Sequence representation

Since there is lot of redundant information in the images for estimating pose
and depth, the optical flow between consecutively sampled frames is extracted
and computed into a 128 dimensional vector, using and encoder. Let us denote
it by $c_{t}$. These $c_{t}$ vectors are passed to LSTM to leverage the
correlations in sequence. The output of LSTM at each time stamp, denoted by
${c_{t}}^{\prime}$, is then passed into next module. This ${c_{t}}^{\prime}$
is expected to contain both sequence context and optical flow information.
This way we can encode the spatio-temporal aspect in a video sequence.

###  III-B Depth estimation

Although ${c_{t}}^{\prime}$ contains much information in the form of
correlations, compared to just a single image, we pass the current image also
as an input to depth prediction network at pixel level. In earlier works, they
use a encoder-decoder architecture for estimating depth using a single image.
This is an ill posed problem. But here, instead of just estimating depth using
a single image, we pass the encoded LSTM output at the current time to
estimate the depth using an up-sampling decoder. DispNet [32], an architecture
which used encoder-decoder design with skip connections and multi-scale side
predictions, is used for our Depth prediction.

###  III-C Pose estimation

Once we have the depth estimation, we then use it along with the respective
images to regress for the pose prediction. Pose network predicts the
transformation parameters (Rotation (R) and translation (T)) between a source
frame and a target frame. Considering only one pair of source and target
frames, these parameters are used to calculate the corresponding locations in
source for each location in target frame, i.e take a grid/matrix with each
entry containing the pixel location in target frame and multiply with
intrinsics, depth estimation and transformation parameters to get the
corresponding pixel location in source frame from which we have to sample.
Now, we use a bi-linear sampling kernel to create a mask which is multiplied
with the source image (pixel values) to get the predicted view. This is the
same procedure followed in [33]. Thus the output of PoseNet is used to warp
the source image to get the target image. Now, whatever loss we compute on the
predicted view (such as reconstruction loss) can be back-propagated to train
the PoseNet parameters as the sampling procedure described above is nothing
but a series of matrix multiplications. For robustness, we calculate
explainability mask to apply loss only on non-occluded pixels.

###  III-D Adversarial training

As photo-metric loss can get stuck at a local minima due to occlusions,
dynamic objects and texture-less regions, using an adversarial training step
can lead to a better depth and pose prediction (better learning of
parameters). Since, we have access to the original frame (the current time
step frame) that we are trying to predict (using view synthesis) using depth
and pose inputs, we can use a discriminator to evaluate the quality of
prediction and back-propagate this loss to learn depth and pose network
parameters.

###  III-E Loss Function

We now delineate the losses we use to train our model.

##### Appearance Loss

Appearance loss is given by:

$\mathcal{L}_{ap}=\mathcal{L}_{reg}(\hat{M})+(1-\alpha)\mathcal{L}_{pho}+\frac{1}{N}\mathcal{L}_{ssim}$
(1)

where N is the minibatch size and $\alpha$ is set to 0.85.
$\hat{M}$ is the per-pixel explainability mask which the model predicts and it
is a belief on whether a pixel is occluded in target frame or not.
$\mathcal{L}_{reg}$ therefore helps in preventing the trivial solution where
$\hat{M}$ will just become 0, thereby making $\mathcal{L}_{pho}=0$ too. It
does so by penalizing the network for predicting values smaller than 0. It,
therefore, works as a regularizer.
$\mathcal{L}_{pho}$ is the loss over how well the predicted image matches the
target image. It is given by:

$\mathcal{L}_{pho}=\sum_{<I_{1},...,I_{N}>}\sum_{p}\hat{M}_{t}(p)||\hat{I}_{t}(p)-I_{t}(p)||_{1}$
(2)

where $<I_{1},...,I_{N}>$ are the sequence frames, p is the number of pixels
in the image, $\hat{I}$ represents the predicted image and $I$ is the target
image.
Finally, $\mathcal{L}_{ssim}$ [34] is the loss over perceived quality of the
image and is given by:

$\mathcal{L}_{ssim}=\sum_{x,y}\alpha\frac{1-SSIM(\hat{I}(x,y),I(x,y))}{2}$ (3)

We use a filter size of $10\times 10$ for SSIM.

##### Trajectory Loss

The trajectory loss imposes a geometry constraint over the poses predicted by
the network. Suppose A, B, C, D are 4 frames of a sequence. Let $A\rightarrow
B$ be the 6-DOF relative transformation that transforms from state A to state
B. Given this, the trajectory loss penalizes the network when the relative
transformation $A\rightarrow D$ predicted by the network deviates from the
$A\rightarrow D$ transformation computed by combining the individual relative
transforms $A\rightarrow B$, $B\rightarrow C$, $C\rightarrow D$. The
trajectory loss is given as:

$\mathcal{L}_{tc}=\frac{1}{N}\sum_{i=1}^{N}\sum_{t\in[2,4,8]}||\hat{p}_{d_{i}}^{i+t}-\hat{p}_{r_{i}}^{i+t}||_{1}$
(4)

where $\hat{p}_{d_{i}}^{i+t}$ is the 6-DOF pose estimated directly for
$A\rightarrow D$, and $\hat{p}_{r_{i}}^{i+t}$ is the pose computed by
concatenating individual relative poses. We apply trajectory loss for
intervals of size 2, 4, and 8.

##### Depth Regularization:

Depth regularization enforces local smoothness and discontinuity in the
predicted depth. Depth discontinuity usually happens where there are strong
image gradients. This loss is given by:

$\displaystyle\begin{split}\mathcal{L}_{smo}&=\frac{1}{N}\sum_{x,y}||\nabla_{x}\hat{D}(x,y)||e^{-||\nabla_{x}I(x,y)||}+\\\
&||\nabla_{y}\hat{D}(x,y)||e^{-||\nabla_{y}I(x,y)||}\end{split}$ (5)

##### GAN Loss:

GAN loss helps ensure that the image generated by the model will look real and
will not have weird artifacts that are not present in real images. It is given
by:

$\displaystyle\begin{split}\mathcal{L}_{GAN}&=\min_{G}\max_{D}V(G,D)\\\
&=E_{I_{t}\sim p_{real}}[log(D(I_{t}|I_{t}))]\\\
&+E_{c_{t-1}^{{}^{\prime}},c_{t}^{{}^{\prime}}\sim
p_{code}}[1-log(D(\hat{I}_{t}|I_{t}))]\end{split}$ (6)

The final loss is the weighted sum of the losses mentioned above and is given
by:

$\mathcal{L}_{final}=\lambda_{a}\mathcal{L}_{ap}+\lambda_{s}\mathcal{L}_{smo}+\lambda_{t}\mathcal{L}_{TC}+\lambda_{g}\mathcal{L}_{GAN}$
(7)

##  IV Experiments

###  IV-A Dataset

The model is evaluated on the KITTI dataset [19] in our experiments. The
dataset contains videos of outdoor sequences. The depth dataset contains pairs
of RGB images along with the depth maps for over 90,000 images. The visual
odometry dataset in the benchmark consists of 22 sequences with half the
sequences used for training with ground truth trajectories available.

###  IV-B Implementation Details

####  IV-B1 Optical Flow

To obtain the optical flow, we evaluated two methods, namely, Farneback method
[35] and Coarse2Fine optical flow method [36]. Among these, we noticed that
although the latter has better quality compared to the former, it is much
slower compared to the former. So, we finalized on Farneback method based on
its much faster speed (10X) and its output (along with image input) seemed
good enough to help depth estimation even though it is not as perfect as
Coarse2Fine method. The preliminary results of our experiment with optical
flow can be seen in Figure 2

Figure 2: Comparison of Farneback and Ce Liu’s Coarse2Fine optical flow
methods

####  IV-B2 Encoder

The optical flow images obtained from above are passed to an Encoder network
to obtain the optical flow encoding. The Encoder network consists of 6 down-
sampling convolution layers followed by an adaptive 2D average pooling layer.
Each down-sampling convolution layer in turn consists of two convolution
layers with stride 2 for the first convolution layer. Batch normalization is
used on the output of each of these convolution layers. The number of channels
in the 6 down-sampling layers are $\\{16,32,64,128,128,128\\}$ respectively
and the kernel sizes are $\\{7,5,3,3,3,3\\}$ respectively (kernel size is
(i,i) for i in the mentioned list). The weights are initialized using
xavier_uniform initialization whereas the biases are initialized to 0.

####  IV-B3 LSTM

The optical flow encodings are further passed to LSTM to leverage the
sequential information present in videos. We used a 1 layer LSTM network with
a hidden dimension of 128 and took the h state at every time as the improved
optical flow encodings.

####  IV-B4 Disparity Network

The optical flow encodings and images are fed to a disparity network which
predicts the disparity(=1/depth). Similar to the previous works [3], [1], we
predict depth at 4 different scales. First the image is passed through a
series of 7 downsample convolution layers, with each downsampling conv layer
same as mentioned earlier. The optical flow encoding is added to the lowest
level (output of $7^{th}$ layer) and the resulting tensor is passed to an
unsampling convolution layer to upsample using transposed convolutions. The
upsampled tensor is concatenated with downsampled tensor at the higher level
and passed to an upsampling convolution layer. This is repeated for all the 7
levels and resulting tensors at first 4 levels form our disparity predictions
after doing a bilinear interpolation followed by another conv layer where we
use Sigmoid activation function and a scaling($\alpha$) and shift($\beta$)
hyperparameters to constrain the disparity values to a certain range. In our
experiments, $\alpha$ and $\beta$ are set to 10 and 0.01 respectively. The
weights are initialized using xavier_uniform initialization whereas the biases
are initialized to 0.

####  IV-B5 Pose Network

The output depth maps are fed to the a pose estimation network along with the
sequence image at the current and the previous time step. The pose estimation
network consists of seven convolution layers with the number of activation
maps at each convolution output equal to $\\{16,32,64,128,256,256,256\\}$
respectively. The kernel sizes for the first two layers is equal to 7 and 5
respectively. The rest of the convolution filters have a kernel size of 3. We
use ReLU as the non-linear activation function. The output of the convolution
layer is fed to two sub-modules. These sub-modules consist of a decoder
architecture similar to the Flownet architecture [37] which outputs an
explainability mask to discard the occluding/moving points which might be
producing erroneous correspondences for pose estimation. Since the range of
the explainability mask is [0,1], we use a Sigmoid activation at the at the
output of the decoder layers to generate the explainability mask(at all
scales). Similar to depth network, we predict the explainability mask at 4
different scales. The second sub-module consists of a single convolution layer
with 6 output channels corresponding to the 6 DOF for pose estimations. This
output represents the relative pose change between the source and the target
image, and is used for the camera trajectory estimation. The weights are
initialized using xavier_uniform initialization whereas the biases are
initialized to 0.

####  IV-B6 GAN

We experimented with multiple discriminator architectures for the adversarial
training. We used WGAN and PatchGAN to study the effects of using a global and
local discriminator.
The authors in [3] have suggested WGAN[5] with the same architecture as the
encoding part of the DispNet. SGD optimiser is used and WGAN. The loss
function here is (Average critic loss of fake images) - (Average critic loss
of real images). The PatchGAN architecture used here is the architecture
proposed in [21] and since then has been widely used for image-image
translation problems. The discriminator tries to classify if a $N\times N$
patch is real or fake. The final output of the discriminator is the average
output of the discriminator results over all the input patches. The
implementation has been borrowed from [38]
The weight for each of the loss terms in the total loss (for all ablation
experiments) can be seen in README section in the code link provided. The
default values of the weights are 0.75, 0.1, 0.14 and 0.01 for appearance,
smoothness, trajectory and GAN loss respectively. We use Adam optimizer in our
experiments.

###  IV-C Evaluation metrics

For evaluating the predicted depth images, we use 6 different metrics as
suggested in [3]. The Absolute difference metric calculates the mean of
absolute value of pixel depth value differences between ground-truth and
predicted depth images. The Absolute relative metric differs from the above by
normalizing the depth value differences using the original(ground-truth) depth
values. The Square relative metric differs from Absolute relative metric by
using square instead of absolute values while calculating the differences.

For the rest of the metrics, we first introduce the term $\delta$ given by
$\delta=max\left(\frac{\text{predicted depth}}{\text{ground truth
depth}},\frac{\text{ground truth depth}}{\text{predicted depth}}\right)$. We
then take three thresholds 1.25, $1.25^{2}$, $1.25^{3}$ and for each of these
values, we calculate how many percentage of the pixels satisfy the condition
$\delta<threshold$. Other than these 6 metrics, we use the standard
ATE(Absolute Trajectory Error) metric for evaluating the trajectory and
thereby poses.

###  IV-D Results

In this section, we present the experimental results of the current pipeline
benefiting from the additional information in form of the sequential inputs
over the erstwhile state-of-the-art SFMLearner pipeline [1]. We also try to
study the effects of different components of the pipeline used by performing
extensive ablation studies. The depth estimations are evaluated using the
metrics described above. The pose estimations are measured using the Absolute
Trajectory Error (ATE) for both translation and rotation.

####  IV-D1 Ablation Studies

In this section, we show the results of ablation studies of our model and
visualize the predicted depths and loss curves for each ablation experiment.
All the experiments use a sequence of 15 images of which we calculate the
optical flow. Now these optical flow may contain redundant information and we
encode these to get get a code. This code is the input to all our ablations.

##### SFMLearner + OpticalFlowCode + LSTM

In this ablation experiment, we augment the SFMLearner model with LSTM and
optical flow computation. The significance of this experiment is we want to
leverage the sequential information in a video to predict better pose and
depth. This will result in better warped image. Simple photo-metric loss is
used to check view synthesis. The predicted depths and warped RGBs are shown
in Figure 3.

Figure 3: Qualitative results and loss curves for SFMLearner + OpticalFlowCode + LSTM ablation. | Abs.Diff | Abs.rel | Sq.rel | $\delta<1.25$ | $\delta<{1.25}^{2}$ | $\delta<{1.25}^{3}$
---|---|---|---|---|---|---
SfmLearner | 5.12 | 0.443 | 4.24 | 0.517 | 0.765 | 0.811
SFMLearner + OpticalFlowCode + LSTM | 4.74 | 0.243 | 2.27 | 0.518 | 0.761 | 0.857
SFMLearner + OpticalFlowCode + LSTM + WGAN | 5.08 | 0.259 | 2.7 | 0.519 | 0.736 | 0.83
SFMLearner + OpticalFlowCode + LSTM + PatchGAN | 4.50 | 0.237 | 2.25 | 0.56 | 0.775 | 0.861
SFMLearner + OpticalFlowCode + WGAN | 5.361 | 0.266 | 3.0 | 0.507 | 0.717 | 0.815
SFMLearner + OpticalFlowCode + LSTM + WGAN + Trajectory Loss | 4.089 | 0.295 | 2.18 | 0.461 | 0.706 | 0.821
SFMLearner + OpticalFlowCode + LSTM + PatchGAN + Trajectory Loss | 4.874 | 0.248 | 2.55 | 0.539 | 0.745 | 0.837
TABLE I: Comparison of different approaches on KITTI dataset for depth
estimation model

##### SFMLearner + OpticalFlowCode + LSTM + WGAN

In this experiment, we add the GAN loss in the model. This is introduced by
adding a discriminator network. The task of discriminator is to find the fake
ones. The warped image is given as input. We use WGAN for this ablation. This
loss is the same as used by authors in [3]. This loss helps in predicting the
real or fake warped image without collapsing. Results are shown in Figure 4.

Figure 4: Qualitative results and loss curves for SFMLearner + OpticalFlowCode
+ LSTM + WGAN ablation.

##### SFMLearner + OpticalFlowCode + LSTM + PatchGAN

Here we perform the ablation by replacing WGAN in previous experiment with
PatchGAN. PatchGan is one of the most popular GAN networks. here we use the
discriminator from this architecture and use MSE loss. The idea is that this
will be able to judiciously distinguish real from fake warps and take care of
visual artifacts. The results are shown in Figure 5

Figure 5: Qualitative results and loss curves for SFMLearner + OpticalFlowCode
+ LSTM + PatchGAN ablation.

##### SFMLearner + OpticalFlowCode + PatchGAN

Now, we remove LSTM from our previous experiment model. The significance of
this experiment is that we want to see if the recurrence relation is actually
captured by LSTM, or Optical flow does a good enough job to encode that
information. The results for this ablation are shown in Figure 6.

Figure 6: Qualitative results and loss curves for SFMLearner + OpticalFlowCode
+ PatchGAN ablation.

##### SFMLearner + OpticalFlowCode + LSTM + PatchGAN + Trajectory Loss

In this final ablation, we include all the components in our model, including
the trajectory constraint. The trajectory loss adds geometric consistency
constraint between relative poses predicted for different frames. This is
especially needed as the other losses do not care about this. As all our
objects in the dataset are rigid bodies, we can use this as our final
experiment. The results are in Figure 7.

Figure 7: Qualitative results and loss curves for SFMLearner + OpticalFlowCode
+ LSTM + PatchGAN + Trajectory Loss ablation.

####  IV-D2 Depth estimation

The quantitative results for the depth estimation module can be seen in Table
I. For first 3 lower values are better and last 3 higher values are better. It
can be seen our model outperforms the model proposed in [1]. We also present
the results of our model along with various versions of our implementation. We
observe that the PatchGAN architecture for the discriminator produces better
results than the WGAN model proposed in [3]. A possible reason for this might
be that since the PatchGAN model tries to discriminate real and fake images at
a local scale, it might lead to images that are more consistent.

####  IV-D3 Pose estimation

The quantitative results for the pose estimation module can be seen in Table
II. The results demonstrate that this architecture outperforms the SFMLearner
architecture by a significant margin. It can be seen here that although the
effect of trajectory loss was not quite seen in the depth estimation results,
it is apparent that trajectory loss is helpful for correct pose estimation. We
see that removal of LSTM has the largest impact on the results, and PatchGAN
performs better than WGAN, which is in line with the observations above.

| Seq. 9 | Seq. 10
---|---|---
SFMLearner | 0.006 | 0.0054
B + Optical Flow + LSTM | 0.0048 | 0.0046
B + Optical FLow + LSTM + WGAN | 0.0045 | 0.0044
B + Optical Flow + LSTM + PatchGAN | 0.0041 | 0.0037
B + Optical flow + WGAN | 0.0052 | 0.0050
B + Optical Flow + LSTM + WGAN + TC | 0.0039 | 0.0037
B + Optical Flow + LSTM + PatchGAN + TC | 0.0034 | 0.0031
TABLE II: Comparison of different approaches on KITTI visual odometry dataset
for pose estimation model. Here, ’+’ means that the model consists of the
components present in the operands on which ’+’ operator is applied. ’B’
refers to the baseline SFMLearner and ’TC’ refers to the trajectory
constraint.

##  V Challenges

Some of the challenges we faced during this project are:

  1. 1.

Optical flow computation: Optical flow can either be computed at runtime or
can be pre-computed for the entire dataset. Pre-computing optical flow
consumes a lot of disk space for KITTI dataset but saves time during training.
Computing optical flow at runtime, even though cheap on storage, increases the
training time, based on the algorithm being used.

  2. 2.

Generative adversarial training: Training GANs come with all the problems
related to it. Some of these issues are mode collapse, non-convergence,
diminished gradient, etc. In WGAN we noticed the values fluctuating a lot and
not having proper convergence. Perhaps that is why the results from WGAN were
not the best. In PatchGAN the values were diminishing but at a slower rate.
The training time was almost 1.5 times the WGAN training time.

  3. 3.

LSTMs: Training LSTMs is also challenging due to problems like
vanishing/exploding gradients, unstable training, etc. We did not keep track
of the LSTM values, but experiments without LSTM gave poorer performance.

  4. 4.

Hyper-parameters: There were lot hyper-parameters to be set in this as there
were many loss functions. The values given by the authors in the [3] did not
work for our implementation, especially for depth regularisation and
trajectory loss. The weight values given by the authors for trajectory loss
would overpower all other losses and make PoseNet converge to the trivial
solution of predicting no change in pose. After that the network would never
recover. We tried lot of values in the limited compute capability we had. The
hyper-parameter selection especially the weights of different loss functions
could be explored further and may give better results.

  5. 5.

Code base: There were many sub-modules in the code and they had to be modular
enough to plug and play for different ablation experiments. Maintaining that
and conceptualising the entire pipeline was difficult.

  6. 6.

GPU availability: We all had limited compute available to train all the
networks for multiple iterations.

  7. 7.

Instability in training: We had Encoder, LSTM, Dispnet, PoseNet and
Discriminator network being trained simultaneously. There was a lot of
instability and they may have gone through collaborative self-supervision,
where in they would have corrected each other’s mistakes to some extent.

  8. 8.

SSIM loss: SSIM ranges from $[-1,1]$ with 1 being equal images. The equation
in [3] suggested the SSIM value as a part of the appearance loss. However,
intuitively, it made more sense to use (1 - SSIM) as a part of the loss
function. We referred multiple papers to confirm the correctness of our
change. The logic is that we want to penalise dissimilar images and increase
loss when the warped image is not similar to GT image.

##  VI Conclusion and Future Work

In this project we have tried to implement [3] paper. We try to leverage
temporal consistency in video sequence to solve Visual odometry as a View
synthesis problem. We take up sequence of input RGB frames, find optical flow
between consecutive frames then pass them through an encoder to get the code.
The code is passed to an LSTM which tries to incorporate long term and short
term information in the code. This information is passed to the Generator
module which consists of a DisparityNet and PoseNet. These 2 are used to
generate a warped target image. This target image is then judged for being
real or fake by our Discriminator. In order to add some novelty we
experimented with PatchGAN as well as WGAN, and found that PatchGAN actually
performed better than WGAN for our experiments, the reason for which we have
already provided above. This project had multiple loss functions and outputs
which harmoniously blend together to solve Visual Odometry task. It took us a
long time to get each module working independently and then together with the
other modules. In the future we would like to do a more exhaustive
hyperparameter search. We would also like to explore more losses that can add
better constraints on the overall system. We have thought of another loss,
i.e. cycle consistency loss which can be used in this. Instead of just image
generation we can do an inverse warp to get the original image back and they
should be same. This might make the system more robust. We also did not
experiment with various architectures for PoseNet and DispNet. We mainly tried
to use the architectures borrowed from SFMLearner. We would also like to
explore some other newer Depth predicting architectures as that might increase
the performance further. We would also want to train on Cityscapes dataset
[39] and compare the result for that as well. As a concluding remark, it was a
great project to work on and it gave us exposure in multiple areas in Vision
and SLAM.

## References

  * [1] Tinghui Zhou, Matthew Brown, Noah Snavely, and David G. Lowe,  “Unsupervised learning of depth and ego-motion from video,”  in CVPR, 2017.
  * [2] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio,  “Generative adversarial nets,”  in Advances in Neural Information Processing Systems 27, Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, Eds., pp. 2672–2680. Curran Associates, Inc., 2014.
  * [3] Shunkai Li, Fei Xue, Xin Wang, Zike Yan, and Hongbin Zha,  “Sequential adversarial learning for self-supervised deep visual odometry,”  in Proceedings of the IEEE International Conference on Computer Vision, 2019, pp. 2851–2860.
  * [4] “Pytorch sfmlearner implementation,” https://github.com/ClementPinard/SfmLearner-Pytorch.
  * [5] Martin Arjovsky, Soumith Chintala, and Léon Bottou,  “Wasserstein gan,”  arXiv preprint arXiv:1701.07875, 2017.
  * [6] Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros,  “Image-to-image translation with conditional adversarial networks,”  arxiv, 2016.
  * [7] Raul Mur-Artal, Jose Maria Martinez Montiel, and Juan D Tardos,  “Orb-slam: a versatile and accurate monocular slam system,”  IEEE transactions on robotics, vol. 31, no. 5, pp. 1147–1163, 2015.
  * [8] Long Quan and Zhongdan Lan,  “Linear n-point camera pose determination,”  IEEE Transactions on pattern analysis and machine intelligence, vol. 21, no. 8, pp. 774–780, 1999.
  * [9] Shaharyar Ahmed Khan Tareen and Zahra Saleem,  “A comparative analysis of sift, surf, kaze, akaze, orb, and brisk,”  in 2018 International conference on computing, mathematics and engineering technologies (iCoMET). IEEE, 2018, pp. 1–10.
  * [10] Huizhong Zhou, Benjamin Ummenhofer, and Thomas Brox,  “Deeptam: Deep tracking and mapping,”  in Proceedings of the European conference on computer vision (ECCV), 2018, pp. 822–838.
  * [11] Sen Wang, Ronald Clark, Hongkai Wen, and Niki Trigoni,  “Deepvo: Towards end-to-end visual odometry with deep recurrent convolutional neural networks,”  in 2017 IEEE International Conference on Robotics and Automation (ICRA). IEEE, 2017, pp. 2043–2050.
  * [12] Benjamin Ummenhofer, Huizhong Zhou, Jonas Uhrig, Nikolaus Mayer, Eddy Ilg, Alexey Dosovitskiy, and Thomas Brox,  “Demon: Depth and motion network for learning monocular stereo,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 5038–5047.
  * [13] Ricson Cheng, Ziyan Wang, and Katerina Fragkiadaki,  “Geometry-aware recurrent neural networks for active visual recognition,”  Advances in Neural Information Processing Systems, vol. 31, 2018.
  * [14] Hsiao-Yu Fish Tung, Zhou Xian, Mihir Prabhudesai, Shamit Lal, and Katerina Fragkiadaki,  “3d-oes: Viewpoint-invariant object-factorized environment simulators,”  arXiv preprint arXiv:2011.06464, 2020.
  * [15] Shamit Lal, Mihir Prabhudesai, Ishita Mediratta, Adam W Harley, and Katerina Fragkiadaki,  “Coconets: Continuous contrastive 3d scene representations,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 12487–12496.
  * [16] Zhou Xian, Shamit Lal, Hsiao-Yu Tung, Emmanouil Antonios Platanios, and Katerina Fragkiadaki,  “Hyperdynamics: Meta-learning object and agent dynamics with hypernetworks,”  arXiv preprint arXiv:2103.09439, 2021.
  * [17] Mihir Prabhudesai, Shamit Lal, Darshan Patil, Hsiao-Yu Tung, Adam W Harley, and Katerina Fragkiadaki,  “Disentangling 3d prototypical networks for few-shot concept learning,”  arXiv preprint arXiv:2011.03367, 2020.
  * [18] Mihir Prabhudesai, Shamit Lal, Hsiao-Yu Fish Tung, Adam W Harley, Shubhankar Potdar, and Katerina Fragkiadaki,  “3d object recognition by corresponding and quantizing neural 3d scene representations,”  arXiv preprint arXiv:2010.16279, 2020.
  * [19] Andreas Geiger, Philip Lenz, and Raquel Urtasun,  “Are we ready for autonomous driving? the kitti vision benchmark suite,”  in Conference on Computer Vision and Pattern Recognition (CVPR), 2012.
  * [20] Ashutosh Saxena, Min Sun, and Andrew Y Ng,  “Make3d: Learning 3d scene structure from a single still image,”  IEEE transactions on pattern analysis and machine intelligence, vol. 31, no. 5, pp. 824–840, 2008.
  * [21] Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros,  “Image-to-image translation with conditional adversarial networks,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 1125–1134.
  * [22] Alec Radford, Luke Metz, and Soumith Chintala,  “Unsupervised representation learning with deep convolutional generative adversarial networks,”  arXiv preprint arXiv:1511.06434, 2015.
  * [23] Richard Zhang, Phillip Isola, and Alexei A Efros,  “Colorful image colorization,”  in Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part III 14. Springer, 2016, pp. 649–666.
  * [24] Yongcheng Jing, Yezhou Yang, Zunlei Feng, Jingwen Ye, Yizhou Yu, and Mingli Song,  “Neural style transfer: A review,”  IEEE transactions on visualization and computer graphics, vol. 26, no. 11, pp. 3365–3385, 2019.
  * [25] Leon A Gatys, Alexander S Ecker, and Matthias Bethge,  “Image style transfer using convolutional neural networks,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 2414–2423.
  * [26] Maayan Frid-Adar, Idit Diamant, Eyal Klang, Michal Amitai, Jacob Goldberger, and Hayit Greenspan,  “Gan-based synthetic medical image augmentation for increased cnn performance in liver lesion classification,”  Neurocomputing, vol. 321, pp. 321–331, 2018.
  * [27] Jiahui Yu, Zhe Lin, Jimei Yang, Xiaohui Shen, Xin Lu, and Thomas S Huang,  “Free-form image inpainting with gated convolution,”  in Proceedings of the IEEE/CVF international conference on computer vision, 2019, pp. 4471–4480.
  * [28] Raymond A Yeh, Chen Chen, Teck Yian Lim, Alexander G Schwing, Mark Hasegawa-Johnson, and Minh N Do,  “Semantic image inpainting with deep generative models,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 5485–5493.
  * [29] Shamit Lal, Vineet Garg, and Om Prakash Verma,  “Automatic image colorization using adversarial training,”  in Proceedings of the 9th International Conference on Signal Processing Systems, 2017, pp. 84–88.
  * [30] Christian Ledig, Lucas Theis, Ferenc Huszár, Jose Caballero, Andrew Cunningham, Alejandro Acosta, Andrew Aitken, Alykhan Tejani, Johannes Totz, Zehan Wang, et al.,  “Photo-realistic single image super-resolution using a generative adversarial network,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 4681–4690.
  * [31] Tero Karras, Timo Aila, Samuli Laine, and Jaakko Lehtinen,  “Progressive growing of gans for improved quality, stability, and variation,”  arXiv preprint arXiv:1710.10196, 2017.
  * [32] Nikolaus Mayer, Eddy Ilg, Philip Hausser, Philipp Fischer, Daniel Cremers, Alexey Dosovitskiy, and Thomas Brox,  “A large dataset to train convolutional networks for disparity, optical flow, and scene flow estimation,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 4040–4048.
  * [33] Max Jaderberg, Karen Simonyan, Andrew Zisserman, et al.,  “Spatial transformer networks,”  in Advances in neural information processing systems, 2015, pp. 2017–2025.
  * [34] Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli,  “Image quality assessment: from error visibility to structural similarity,”  IEEE transactions on image processing, vol. 13, no. 4, pp. 600–612, 2004.
  * [35] “Opencv’s farneback method,” https://docs.opencv.org/2.4/modules/video/doc/motion_analysis_and_object_tracking.html#calcopticalflowfarneback.
  * [36] “Ce liu’s c++ implementation of coarse2fine optical flow,” https://people.csail.mit.edu/celiu/OpticalFlow/.
  * [37] Alexey Dosovitskiy, Philipp Fischer, Eddy Ilg, Philip Hausser, Caner Hazirbas, Vladimir Golkov, Patrick Van Der Smagt, Daniel Cremers, and Thomas Brox,  “Flownet: Learning optical flow with convolutional networks,”  in Proceedings of the IEEE international conference on computer vision, 2015, pp. 2758–2766.
  * [38] “Pytorch gan implementation,” https://github.com/eriklindernoren/PyTorch-GAN.
  * [39] Marius Cordts, Mohamed Omran, Sebastian Ramos, Timo Rehfeld, Markus Enzweiler, Rodrigo Benenson, Uwe Franke, Stefan Roth, and Bernt Schiele,  “The cityscapes dataset for semantic urban scene understanding,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 3213–3223.


# Structural Gaussian Priors for Bayesian CT reconstruction of Subsea Pipes

Silja L. Christensen, Nicolai A. B. Riis, Felipe Uribe, and Jakob S. Jørgensen
CONTACT Silja L. Christensen. Email<EMAIL_ADDRESS>Department of Applied
Mathematics and Computer Science, Technical University of Denmark, Kongens
Lyngby, Denmark

###### Abstract

A non-destructive testing (NDT) application of X-ray computed tomography (CT)
is inspection of subsea pipes in operation via 2D cross-sectional scans. Data
acquisition is time-consuming and costly due to the challenging subsea
environment. While reducing the number of projections in a CT scan yields time
and cost savings, this compromises the reconstruction quality when
conventional methods are used. To address this issue we take a Bayesian
approach to CT reconstruction and focus on designing an effective prior that
introduces additional information to the problem. We propose a new class of
_structural Gaussian priors_ to enforce expected material properties in
different regions of the reconstructed image based on available a-priori
information about the pipe’s layered structure. The prior is composed from
Gaussian distributions such that we can define an efficient computational
strategy to sample the resulting posterior distribution. This is essential in
practical NDT applications for fast processing of the large-scale data we see
in CT. Numerical experiments with synthetic and real data show that the
proposed structural Gaussian prior reduces artifacts and enhances contrast in
the reconstruction, compared to using only a conventional Gaussian Markov
Random Field (GMRF) prior or no prior information at all.

###### keywords:

X-ray computed tomography; Uncertainty Quantification; Bayesian Inversion;
Non-destructive testing; Pipeline inspection

##  1 Introduction

Subsea pipelines are used to transport oil and gas around the world. It is
crucial that the subsea pipes are in good condition to avoid environmentally
damaging leaks. Therefore, efficient methods for non-destructive inspection
are in demand. One suitable method is X-ray computed tomography (CT) that is
used to obtain 2D cross-sectional images of the pipes using data from a finite
number of projections [1]. To produce high-quality images with limited noise
and artifacts, standard CT reconstruction methods require lengthy measurements
at many projection (view) angles. To reduce costs, it is desirable to limit
measurement time, and this can be achieved by reducing the number of view
angles. However, if standard iterative reconstruction methods are used, this
compromises the reconstruction quality significantly [2]. Moreover, since the
data is typically scarce and noisy, there is substantial uncertainty affecting
the reconstructions. Therefore, we follow a Bayesian approach, where the idea
is to update a _prior_ probability distribution of the image by including
information about the CT model and projection data to obtain a so-called
_posterior_ distribution describing the reconstruction and its uncertainty.
Our idea is to use prior knowledge of the pipes’ internal structure to improve
the image reconstruction quality, even for cases with measurements taken at
only few view angles.

Previous work has taken a deterministic approach to CT scanning of subsea
pipes [3]. The authors use microlocal analysis to design a favorable offset
scan geometry and propose a reconstruction method based on shearlets with a
weighted sparsity penalty. However, this method does not account for and
quantify the uncertainty in the reconstructions, and thus a Bayesian approach
is needed for a reliable solution of the CT problem.

In Bayesian CT reconstruction [4, 5, 6], it is important to choose an
appropriate prior. Typical choices in this setting are Markov Random Field
(MRF) priors based on multivariate Gaussian and Laplacian distributions. These
represent a statistical parallel to the well-known Tikhonov and total
variation (TV) regularization methods, which are commonly used in inverse
problems and in particular CT [1, 7]. An overview of several MRF type priors
applied to 2D inverse problems (including CT) is given in [8] and more
detailed descriptions can be found in earlier papers [9, 10]. A recent
approach is the Cauchy differences prior (also a MRF-type prior), which is
edge-preserving and discretization invariant [11, 12]. All MRF priors have
regularizing effects, and depending on the scanned object, a smoothness-
promoting or edge-preserving prior can be chosen. However, a common
characteristic for these distributions is that they do not model specific
features, such as the layer structure in subsea pipes, when used as priors in
Bayesian reconstruction. Furthermore, non-Gaussian MRF-type priors require
more advanced sampling techniques, which in turn are also computationally
expensive and can therefore be prohibitive for large-scale problems.

If a particular structure of the target image is known in advance, a stronger
prior can be defined. A series of studies use this idea by forming a texture-
preserving prior for low-dose CT reconstruction using a previous normal-dose
scan [13, 14, 15, 16]. Another study exploits the fact that the human body
consists of a finite number of tissue types with well-known attenuation
coefficients and the study designs the prior accordingly including an inherent
classification problem [17].

The driving idea in this work is to utilize as much prior information as
possible. We know the overall layout of the pipes because they are
manufactured according to specific design standards. This leads to our
proposed structural Gaussian prior that introduces information about the
pipes’ layered structure as well as their constitutive materials. Numerical
experiments show that this prior improves the reconstruction quality compared
to standard priors, especially if only few view angles are available.
Moreover, the Bayesian approach allows probabilistic modelling of noise and
prior information, and through sampling of the resulting posterior
distribution, we can obtain reconstructions with uncertainty estimates.
Another advantage of the structural Gaussian prior is that it can be sampled
efficiently in large-scale problems.

Our paper is organized as follows. In Section 2 we introduce the Bayesian
approach to CT reconstruction. Section 3 describes our proposed structural
Gaussian prior. In Section 4 we derive the posterior and how to sample it
efficiently. In Section 5 we perform numerical experiments with synthetic and
real data, and in Section 6 we discuss and conclude our findings.

##  2 Bayesian CT problem

Deep sea oil pipes are large and consist of dense materials including steel
and concrete. In order to penetrate the pipe with X-rays, it is necessary to
use high power and a narrow beam source. The narrow beam cannot illuminate the
full pipe, but an exterior tomography setup can be used. In this case, the
X-ray source is offset from the pipe centre, which is also the rotation axis.
It can be shown that this setup is ideal, if we are interested in recovering
the pipe layer structure and defects therein rather than the center of the
pipe [3]. We use real data acquired from a prototype set-up of an oil pipe CT
scanner. A linear detector panel with 510 sensors of sizes 0.8 mm each was
used, and 360 equi-angular projections were recorded in a full 360 degree
rotation. The pipe, its specifications, and the offset fan-beam acquisition
geometry are illustrated in Figure 1.

###  2.1 CT model

If we discretize the reconstruction domain in pixels, the inverse problem
associated with CT can be expressed as a linear system of equations,

$\textbf{d}=\textbf{Ax}+\textbf{e},\qquad\textbf{e}\sim\mathcal{N}\left(\textbf{0},\frac{1}{\lambda}\textbf{I}_{n}\right),$
(1)

where $\textbf{d}\in\mathbb{R}^{m}$ is the observed X-ray absorption (known as
a sinogram), $\textbf{x}\in\mathbb{R}^{n}$ is the unknown vector of linear
attenuation coefficients in the reconstruction domain with $N\times N=n$
pixels, $N$ is the domain discretization size, and
$\textbf{A}\in\mathbb{R}^{m\times n}$ contains entries $a_{ij}$ that describe
the intersection of ray $i$ with pixel $j$ (see, e.g., [1] for details). We
model the data noise as additive Gaussian, suitable for the high-intensity
X-rays used in the data collection process. That is,
$\textbf{e}\in\mathbb{R}^{m}$ is independent and identically distributed
Gaussian noise with mean zero and precision parameter $\lambda$, and
$\textbf{I}_{n}\in\mathbb{R}^{n\times n}$ is an identity matrix.

Figure 1: Left: Photograph of subsea oil pipe. Right: Pipe and offset scan
geometry specifications; dimensions in mm. The CT rotation axis is marked with
a star.

###  2.2 Bayesian inversion

We take a Bayesian approach to the CT problem. Using Bayes’ Theorem we can
characterize the solution in terms of the posterior distribution given by the
probability density function (PDF),

$\pi_{\text{post}}(\textbf{x}\,|\,\textbf{d})\propto\pi_{\text{pr}}(\textbf{x})\pi_{\text{li}}(\textbf{d}\,|\,\textbf{x}).$
(2)

Here, $\pi_{\text{pr}}(\textbf{x})$ denotes the prior PDF describing some
known properties of the attenuation coefficients x. The prior is modelled as a
Gaussian distribution with PDF given by

$\pi_{\text{pr}}(\textbf{x})=\frac{1}{\left(2\pi\right)^{n/2}}\left(\left|\textbf{R}^{T}\textbf{R}\right|\right)^{1/2}\exp\left(-\frac{1}{2}\lVert\textbf{R}(\textbf{x}-\bm{\mu})\rVert_{2}^{2}\right),$
(3)

where $\bm{\mu}\in\mathbb{R}^{n}$ is the mean and
$(\textbf{R}^{T}\textbf{R})\in\mathbb{R}^{n\times n}$ the precision matrix
with full rank. We refer to R as the square-root precision matrix, and
$|\cdot|$ denotes the matrix determinant. Below, in Section 3, the specifics
of how we build the pipe structure into the Gaussian prior are given.

From (1), we obtain the distribution of the data d given some known vector of
linear attenuation coefficients x,

$\pi_{\text{li}}(\textbf{d}\,|\,\textbf{x})=\left(\frac{\lambda}{2\pi}\right)^{m/2}\exp\left(-\frac{\lambda}{2}\left\|\textbf{A}\textbf{x}-\textbf{d}\right\|_{2}^{2}\right),$
(4)

which can be seen as a likelihood function of x given some observed data d.

Because the prior and likelihood are Gaussian, we can obtain a closed-form
expression for the posterior, which is also Gaussian (see Section 4.1). In
general, it is numerically impractical to compute the posterior precision
matrix, although the posterior mean can be simple to obtain. This is because
the distribution is Gaussian, and thus the mean is equivalent to the maximum
a-posteriori probability (MAP) estimate, which is often used as a point
estimator to describe the posterior, as it can be easily computed through
optimization. Nevertheless, the MAP estimator itself does not provide any
uncertainty information, and therefore does not reach the full potential of
the Bayesian approach. Indeed, one of the strong points of Bayesian inversion
is the ability to give uncertainty estimates. Therefore, our approach entails
a method for sampling the posterior distribution as described in Section 4.2,
instead of focusing on a single point estimate.

##  3 Structural Gaussian priors

We must specify a prior in order to perform Bayesian CT reconstruction. We
work with Gaussian priors because they have nice mathematical properties
making them easy to combine, and because they can be sampled efficiently
compared to more general distribution models. First, we outline the prior
structural information that we have. Then we describe our proposed method of
incorporating the known pipe geometry into what we call structural Gaussian
priors (SGPs). The SGP combines (stacks) priors for each of the materials in
the pipe, thereby enforcing a different prior in each region. In addition to
the pipe material and geometry information, the SGP also contains a prior on
the pixel differences of the whole image. The purpose of this is to add a
global smoothness assumption on the image, which also ensures full rank of the
precision matrix of the combined Gaussian prior in (3).

###  3.1 Prior structural information

We have prior information about the pipe’s layered structure and constitutive
materials, which we will incorporate in the prior probability model. This is
the core idea of our proposed SGP. To build the SGP we need: 1) masks for each
pipe layer and for the background, and 2) estimates of the linear attenuation
coefficients of the layer materials. These are described in the following
subsections.

Figure 2: Left: Real pipe with four layers of different materials and a
background of air. Right: SGP masks for $i=1,\ldots,5$ that each represent a
material with attenuation coefficient $\alpha_{i}$ given in (5). The white
regions do not have any prior promoting specific attenuation coefficients. The
index $i=0$ is used to denote the prior on the whole image domain (background
and pipe).

####  3.1.1 Pipe layer masks

Given $p$ non-overlapping regions of different materials, we represent each
region by a mask, i.e., a binary vector $\textbf{m}_{i}\in\mathbb{R}^{n}$ for
$i=1,\ldots,p$, in which the value of element $j$ is 1, if pixel $j$ is in
region $i$ and $0$ if not, for $j=1,\ldots,n$. Figure 2 shows a reconstruction
of the real pipe and the corresponding masks (with $p=5$). We use masks that
are slightly smaller than the actual pipe layers to account for modelling
error in the pipe structure and scan geometry. The masks are constructed based
on a high-quality standard CT reconstruction of the real pipe. This may also
be obtained from lower-quality data or from available CAD drawings. We further
define mask matrices $\textbf{M}_{i}\in\mathbb{R}^{l_{i}\times n}$,
$i=1,\ldots,p$, that pick out the $l_{i}$ pixels of region $i$ of the
$n$-dimensional image. Specifically, $\textbf{M}_{i}$ is formed by keeping
only the $l_{i}$ rows of the identity matrix $\textbf{I}_{n}$ corresponding to
where there are 1-values in the mask $\textbf{m}_{i}$.

####  3.1.2 Linear attenuation coefficients of pipe materials

When the X-ray energy is known, we can estimate the linear attenuation
coefficients $\alpha_{i}$ (ignoring scattering) for the materials as [18]:

$\alpha_{i}=1/\ell_{i}=\kappa_{i}(E)\rho_{i},\qquad i=1,\ldots,p,$ (5)

where $\rho_{i}$ is the material density and $\kappa_{i}(E)$ is the mass
attenuation coefficient at X-ray energy $E$ that can be found, e.g., in the
NIST database [19]. Equation (5) also indicates that the linear attenuation
coefficient is sometimes given in terms of the mean free path of the material
$\ell_{i}$.

For thick high-density materials, the scattering of X-rays cannot be ignored
and it becomes necessary to account for a physical effect referred to as
buildup in the material. Consider the X-ray intensity at some point in the
material. Some of the intensity will be due to photons that have reached the
point undisturbed (primary) and some of the intensity will be due to scattered
photons (secondary). The buildup is defined at a given point as the ratio of
the total (primary and secondary) intensity to the primary intensity. With
buildup, both primary and secondary intensity reach the detector, which then
gives a higher intensity reading compared to the case where buildup is not
present, and only primary intensity reaches the detector. With a so-called
buildup factor we can account for the effect.

Assuming a homogeneous material, the X-ray attenuation including buildup is
described as [20]:

$I=I_{0}\,B(E,w_{i}/\ell_{i})\exp\left(-w_{i}/\ell_{i}\right),$ (6)

where $I_{0}$ is the X-ray intensity entering the homogeneous material, $I$ is
the X-ray intensity leaving the homogeneous material, $w_{i}$ is the width of
the homogeneous material $i$, and $B(E,w_{i}/\ell_{i})$ is the buildup factor
which can be found from e.g. [21]. However, in most CT implementations, the
forward operator A is constructed ignoring buildup and assumes that the X-ray
attenuation follows Lambert–Beer’s law:

$I=I_{0}\exp\left(-\alpha_{i}w_{i}\right).$ (7)

Using Lambert–Beer’s law for CT reconstruction of dense materials is
essentially a model error, and it will lead to linear attenuation coefficients
that are lower than expected according to (5) because the secondary intensity
is measured but not included in the model. Instead of correcting the model
error, we adjust our prior expectations about the linear attenuation
coefficient to include the buildup factor. By equating the right-hand sides of
(6) and (7), and inserting (5) instead of $1/\ell_{i}$, we estimate the linear
attenuation coefficient including the buildup factor for a dense material as:

$\alpha_{i}=\kappa_{i}(E)\rho_{i}-\frac{\log\left(B(E,w_{i}\kappa_{i}(E)\rho_{i})\right)}{w_{i}}.$
(8)

In table 1 we will use (5) for prior information about concrete, polyurethane,
and polyethylene, and we will use (8) for prior information about the much
denser material, steel.

###  3.2 Structural Gaussian priors

To form the SGP, we use Gaussian distributions defined in terms of means
$\bm{\mu}_{i}\in\mathbb{R}^{n}$ and a symmetric positive semi-definite
precision matrices
$\left(\textbf{R}_{i}^{T}\textbf{R}_{i}\right)\in\mathbb{R}^{n\times n}$. In
(3) we gave an expression for Gaussians with precision matrices of full rank.
When forming the SGP, we will need intrinsic Gaussian distributions that can
be defined for rank-deficient precision matrices of rank $k$, where $n\geq
k>0$ . The intrinsic Gaussian prior density is defined as [22, Ch. 3]:

$\pi_{i}(\textbf{x})=\frac{1}{\left(2\pi\right)^{k/2}}\left(\left|\textbf{R}_{i}^{T}\textbf{R}_{i}\right|^{*}\right)^{1/2}\exp\left(-\frac{1}{2}\lVert\textbf{R}_{i}(\textbf{x}-\bm{\mu}_{i})\rVert_{2}^{2}\right),$
(9)

where $|\cdot|^{*}$ is a generalized determinant given by the product of the
$k$ non-zero eigenvalues of the precision matrix. Note this density is
multivariate Gaussian if and only if $k=n$. In the SGP, we use two different
Gaussian priors that are specified below:

IID:

The independent and identically distributed (IID) Gaussian prior is a noise-
dampening prior that favors posterior pixel values around the prior mean. We
use it to express prior information of expected linear attenuation
coefficients $\alpha_{i}$ in the $p$ different materials. The prior is defined
as an intrinsic Gaussian parameterized by mean and square-root precision:

$\displaystyle\bm{\mu}_{i}=\alpha_{i}\textbf{1}\qquad\text{and}\qquad\textbf{R}_{i}=\sqrt{\delta_{i}}\textbf{M}_{i},$
(10)

where $\textbf{1}\in\mathbb{R}^{n}$ is a vector of ones and $\delta_{i}$ is a
prior precision parameter.

GMRF:

The Gaussian Markov Random Field (GMRF) prior is a smoothing prior that
promotes neighbouring pixels to take the same value. The pipe materials can be
considered roughly homogeneous and therefore the GMRF prior is enforced on the
entire image domain. We index this prior with $i=0$ corresponding to the whole
image, and in this case the mean and square-root precision are:

$\displaystyle\bm{\mu}_{0}=\textbf{0},\qquad\text{and}\qquad\textbf{R}_{0}=\sqrt{\delta_{0}}\begin{bmatrix}\textbf{D}_{1}\\\
\textbf{D}_{2}\end{bmatrix},$ (11)

where $\textbf{0}\in\mathbb{R}^{n}$ is a vector of zeros,
$\textbf{D}_{1}=\textbf{I}_{N}\otimes\textbf{D}$,
$\textbf{D}_{2}=\textbf{D}\otimes\textbf{I}_{N}$,
$\textbf{D}\in\mathbb{R}^{(N+1)\times N}$ is a finite difference matrix with
backward differences and zero Dirichlet boundary conditions, and $\otimes$
denotes the Kronecker product. Note that $\textbf{R}_{0}$ has full column
rank. Therefore the precision matrix has full rank, and the density is thus
multivariate Gaussian.

To construct the SGP prior, we enforce IID priors on the masked image regions
and a GMRF prior on the whole image domain as described above. Under the
assumption that the priors $\pi_{i}(\textbf{x})$, $i=0,\ldots,p$ are
independent of each other, the SGP prior is given as a product of the
individual priors:

$\pi_{\text{pr}}(\textbf{x})=\prod_{i=0}^{p}\pi_{i}(\textbf{x}).$ (12)

Then, the SGP prior is also Gaussian with mean and square-root precision:

$\textbf{R}_{\text{pr}}=\begin{bmatrix}\textbf{R}_{0}\\\ \vdots\\\
\textbf{R}_{p}\end{bmatrix}\qquad\text{and}\qquad\bm{\mu}_{\text{pr}}=\left(\textbf{R}_{\text{pr}}^{T}\textbf{R}_{\text{pr}}\right)^{-1}\sum_{i=0}^{p}\textbf{R}_{i}^{T}\textbf{R}_{i}\bm{\mu}_{i}.$
(13)

Note, that the full column rank of $\textbf{R}_{0}$ leads to
$\textbf{R}_{\text{pr}}$ also having full column rank, and therefore the prior
precision matrix $\textbf{R}_{\text{pr}}^{T}\textbf{R}_{\text{pr}}$ has full
rank and is invertible.

####  3.2.1 SGP configurations

We use the pipe shown in Figure 2 as a case study for numerical experiments in
Section 5. The pipe consists of four different materials, steel in region 2,
polyurethane (PU) foam in region 3, polyethylene (PE) rubber in region 4, and
concrete in region 5. The pipe is scanned with air as the background and this
will be denoted region 1. Thus we have $p=5$ masked regions where we can
impose IID priors promoting the expected attenuation coefficients. We study
two configurations of the SGP:

SGP-BG:

Assuming we have information about the pipe’s dimensions but not internal
structure, we use a GMRF on the whole domain and add an IID prior on the air
background (‘BG’) such that
$\pi_{\text{pr}}(\textbf{x})=\prod_{i=0}^{1}\pi_{i}(\textbf{x})$.

SGP-F:

Assuming we also have information about the internal pipe structure, we again
use a GMRF on the full (‘F’) domain and add IID priors on the background and
on the pipe layers such that
$\pi_{\text{pr}}(\textbf{x})=\prod_{i=0}^{5}\pi_{i}(\textbf{x})$.

##  4 Gaussian posterior and efficient sampling

In Sections 2 and 3, we assumed Gaussian distributions for the likelihood and
priors respectively, which in turn yields a Gaussian posterior PDF [4]. This
is mathematically convenient because we can derive a closed-form expression
for this distribution and sample it efficiently. Below, we define the
posterior PDF and then show how to sample it by solving a system of linear
equations.

###  4.1 Posterior PDF

Using the likelihood in (4) and the prior in (12)-(13), the posterior PDF is
given by Bayes’ Theorem (2) as:

$\displaystyle\pi_{\text{post}}(\textbf{x}\,|\,\textbf{d})\propto{}$
$\displaystyle\left(\frac{\lambda}{2\pi}\right)^{\frac{m}{2}}\exp\left[-\frac{\lambda}{2}\left\|\textbf{A}\textbf{x}-\textbf{d}\right\|_{2}^{2}\right]\times\prod_{i=0}^{p}\left(\left|\frac{\textbf{R}_{i}^{T}\textbf{R}_{i}}{2\pi}\right|^{*}\right)^{\frac{1}{2}}\exp\left[-\frac{1}{2}\lVert\textbf{R}_{i}(\textbf{x}-\bm{\mu}_{i})\rVert_{2}^{2}\right]$
$\displaystyle\propto{}$
$\displaystyle\exp\left[-\frac{\lambda}{2}\left\|\textbf{A}\textbf{x}-\textbf{d}\right\|_{2}^{2}-\frac{1}{2}\sum_{i=0}^{p}\lVert\textbf{R}_{i}(\textbf{x}-\bm{\mu}_{i})\rVert_{2}^{2}\right].$
(14)

The above expression defines the Gaussian random vector

$\textbf{x}\,|\,\textbf{d}\sim\mathcal{N}(\bm{\mu}_{\mathrm{post}},(\textbf{R}^{T}_{\mathrm{post}}\textbf{R}_{\mathrm{post}})^{-1}),$
(15)

where the posterior square-root precision is given by

$\textbf{R}_{\text{post}}=\begin{bmatrix}\sqrt{\lambda}\textbf{A}\\\
\textbf{R}_{0}\\\ \vdots\\\ \textbf{R}_{p}\end{bmatrix}\in\mathbb{R}^{q\times
n}\qquad\text{with}\qquad
q=\left(m+2N\left(N+1\right)+\sum_{i=1}^{p}l_{i}\right),$ (16)

and the posterior mean is given implicitly by the expression:

$\displaystyle\textbf{R}_{\text{post}}^{T}\textbf{R}_{\text{post}}\bm{\mu}_{\text{post}}=\lambda\textbf{A}^{T}\textbf{d}+\sum_{i=0}^{p}\textbf{R}_{i}^{T}\textbf{R}_{i}\bm{\mu}_{i}.$
(17)

Since $\textbf{R}_{0}$ and $\textbf{R}_{\text{post}}$ have full column rank,
the posterior precision is invertible, and the posterior distribution is
Gaussian with well-defined covariance. That is, the closed-form expression for
the posterior distribution is (15), which can be seen by factoring out (14)
and reordering the equation into the form (3).

The large-scale nature of CT makes it infeasible to form the matrix A
explicitly, which can be avoided by matrix-free implementations of the
application of A and its adjoint. Then, the posterior mean can be obtained by
solving (17). However, to compute the posterior precision, we must form
$\textbf{A}^{T}\textbf{A}$, which is numerically impractical. Instead we
sample the posterior as described below, and statistics of the samples are
used to obtain an uncertainty estimate.

###  4.2 Posterior Sampling

Our goal is to efficiently sample the posterior (15) without explicitly
forming the matrix A. Inspired by [6], we generate posterior realizations by
solving the linear system of equations:

$\textbf{R}_{\text{post}}\textbf{x}^{*}=\begin{bmatrix}\sqrt{\lambda}\textbf{d}\\\
\textbf{R}_{0}\bm{\mu}_{0}\\\ \vdots\\\
\textbf{R}_{p}\bm{\mu}_{p}\end{bmatrix}+\bm{\xi},\quad\bm{\xi}\sim\mathcal{N}(\textbf{0},\textbf{I}_{q}).$
(18)

Note that the solution method must support matrix-free implementation of A.

We can show that $\textbf{x}^{*}$ is distributed according to the posterior by
first forming the normal equations of the linear system above, i.e.

$\displaystyle\left(\lambda\textbf{A}^{T}\textbf{A}+\sum_{i=0}^{p}\textbf{R}_{i}^{T}\textbf{R}_{i}\right)\textbf{x}^{*}$
$\displaystyle=\lambda\textbf{A}^{T}\textbf{d}+\sum_{i=0}^{p}\textbf{R}_{i}^{T}\textbf{R}_{i}\bm{\mu}_{i}+\textbf{R}_{\text{post}}^{T}\bm{\xi}.$
(19)

Inserting (16) and solving for $\textbf{x}^{*}$ we find

$\displaystyle\textbf{x}^{*}$
$\displaystyle=\left(\textbf{R}_{\text{post}}^{T}\textbf{R}_{\text{post}}\right)^{-1}\left(\lambda\textbf{A}^{T}\textbf{d}+\sum_{i=0}^{p}\textbf{R}_{i}^{T}\textbf{R}_{i}\bm{\mu}_{i}\right)+\left(\textbf{R}_{\text{post}}^{T}\textbf{R}_{\text{post}}\right)^{-1}\textbf{R}_{\text{post}}^{T}\bm{\xi}$
$\displaystyle=\bm{\mu}_{\text{post}}+\textbf{R}_{\text{post}}^{\dagger}\bm{\xi},$
(20)

where
$\textbf{R}_{\text{post}}^{\dagger}=(\textbf{R}_{\text{post}}^{T}\textbf{R}_{\text{post}})^{-1}\textbf{R}_{\text{post}}^{T}$
is the Moore–Penrose pseudoinverse. Note that $\textbf{x}^{*}$ is a linear
transformation of the Gaussian random variable $\bm{\xi}$. As $\textbf{R}_{0}$
has full column rank, $\textbf{R}_{\text{post}}^{\dagger}$ also has full row
rank. It follows that $\textbf{x}^{*}$ is Gaussian distributed with mean and
covariance [23]:

$\displaystyle\mathrm{E}(\textbf{x}^{*})$
$\displaystyle=\bm{\mu}_{\text{post}}+\textbf{R}_{\text{post}}^{\dagger}\textbf{0}=\bm{\mu}_{\text{post}}$
(21a) $\displaystyle\mathrm{Cov}(\textbf{x}^{*})$
$\displaystyle=\textbf{R}_{\text{post}}^{\dagger}\textbf{I}(\textbf{R}_{\text{post}}^{\dagger})^{T}=\left(\textbf{R}_{\text{post}}^{T}\textbf{R}_{\text{post}}\right)^{-1}.$
(21b)

This confirms that $\textbf{x}^{*}$ is distributed according to the target
posterior distribution. That is, by solving (18) we obtain samples of CT
reconstruction images distributed according to (15).

##  5 Numerical Experiments

FORCE Technology111forcetechnology.com has provided CT data from a deep-sea
oil pipe acquired in a laboratory experiment using the offset fan-beam scan
geometry illustrated in Figure 1. We conduct synthetic and real data CT
experiments to assess the performance of the proposed SGPs for obtaining
useful CT reconstructions and relevant uncertainty estimates. The results are
compared to a conventional deterministic reconstruction method where (1) is
solved with the conjugate gradient least-squares (CGLS) algorithm [24, Ch. 7]
and stopped at semi-convergence. We interpret the posterior mean as the CT
reconstruction and the posterior $95\%$ credible intervals for each pixel as
associated uncertainty estimates.

###  5.1 Data description

####  5.1.1 Phantom

Figure 3: Left: Pipe phantom modelling the real pipe seen in Figure 1. Right:
SGP masks for $i=1,\ldots,5$ that each represent a material with attenuation
coefficient $\alpha_{i}$ given in Table 1. The white regions do not have any
prior promoting specific attenuation coefficients. The index $i=0$ is used to
denote the prior on the whole image domain (background and pipe). Table 1:
Best estimates for the material constants at beam intensity 2 MeV.

$i$ | Material | $\kappa[\text{cm}^{2}/\text{g}]$ | $\rho[\text{g}/\text{cm}^{3}]$ | $B$ | $w[\text{cm}]$ | $\alpha[\text{cm}^{-1}]$
---|---|---|---|---|---|---
1 | Air | 0.044 | 0.0012 | - | - | $\sim 0$
2 | Steela | 0.042 | 7.9 | 2.013 | 4 | 0.16
3 | PU foam | 0.051b | 0.15 | - | - | 0.0077
4 | PE rubber | 0.051 | 0.94 | - | - | 0.048
5 | Concrete | 0.046 | 2.3 | - | - | 0.11

aSubject to build-up. bEstimated by PE rubber’s $\kappa$ due to missing value
in database.

Based on the known pipe structure and materials given in Figure 2, we
construct a phantom and corresponding SGP mask to perform synthetic
experiments (see Figure 3). To avoid inverse crime, the phantom is defined on
a $1024\times 1024$ pixel image and the reconstruction is performed on a
$512\times 512$ pixel image that represents a $55\times 55$ cm domain leading
to pixels of size $\sim 1\times 1$ mm. The linear attenuation coefficients in
each pipe layer are computed as described in Section 3.1.2, with material
constants from the NIST database [19] and buildup factor from [21], assuming a
monochromatic beam of energy 2 MeV222The beam energy is in fact a spectrum
with 2 MeV mean energy.. This beam intensity is an estimate, but around these
energies the mass attenuation coefficients do not vary significantly.
Therefore, we expect the computed linear attenuation coefficients are good
estimates and they are listed in Table 1.

Note, that we add steel inclusions in the phantom’s concrete layer to simulate
steel reinforcement bars in the real pipe. The inclusions are in the radial or
tangential directions and their widths are varied from 2–7 mm. These
inclusions will not be accounted for in the SGP prior; this in order to test
the prior’s ability to reconstruct unexpected features. Furthermore, with the
offset scan geometry, no X-rays will be parallel to near-radial edges, and
therefore we expect that the radial inclusions are more difficult to
reconstruct than the tangential (see, e.g., [3] for more details).

Figure 4: Left: Real subsea pipe sinogram. Right: Synthetic sinogram with 2%
added noise.

####  5.1.2 Sinograms

Using the pipe phantom, offset acquisition geometry, and 2% noise level (i.e.
$\lVert\textbf{e}\rVert_{2}/\lVert\textbf{Ax}\rVert_{2}\sim 2\%$), we simulate
a sinogram with 360 view angles distributed equally over a full rotation.
Assuming the real sinogram has an equivalent noise distribution, we compare
the synthetic and real data (Figure 4). Five bands in the sinograms mark the
centre region, steel, PE, PU, and concrete layers in that order from the
bottom of the sinograms. The bands in the real data sinogram follow sinusoidal
curves because the pipe center did not coincide with the rotation center, and
the phases of the sinusoidal bands vary because the pipe layers are not
completely concentric. The overall similarity between the real data and
synthetic data sinograms indicates that our phantom is realistic and the scan
geometry close to the real setup. For sparse view angle experiments, we
simulate data with only 50%, 20% or 10% of the total view angles in the
sinograms. Note that, since we reconstruct on a $512\times 512$ pixel image,
the CT problem is under-determined in all our experiments.

###  5.2 Numerical implementation

####  5.2.1 Bayesian parameter tuning

The posterior in (14) requires choosing several parameters. The data noise of
2% corresponds to a precision parameter $\lambda\approx 400$. The SGP includes
several IID Gaussians for which we must determine mean and precision
parameters. The means are defined as the estimated linear attenuation
coefficients $\alpha_{i}$ for each pipe material reported in Table 1. The IID
Gaussian precision parameters are chosen heuristically such that the priors
cover reasonable intervals around the estimated linear attenuation
coefficients, representing some uncertainty in the estimates. We expect steel
reinforcement bars in the concrete layer, which are not modeled in the prior,
so the precision parameter is reduced here to allow a large range of possible
values. Hence, the chosen precision parameters are:
$\delta_{1}=\delta_{2}=\delta_{3}=\delta_{4}=1000$ and $\delta_{5}=500$. The
final parameter choice is the GMRF precision $\delta_{0}$. If we fix all other
precision parameters, we can perform a parameter sweep for $\delta_{0}$ values
and choose the value that minimizes the root-mean-square-error (RMSE) with
respect to the synthetic ground truth. An example of $\delta_{0}$ tuning can
be seen in Figure 5 (left). We transfer the GMRF precision choices from the
synthetic to the real data case, and check qualitatively that the
corresponding reconstructions are acceptable.

####  5.2.2 Sampling algorithm

We use the CGLS algorithm to solve (17) for MAP estimates of the posterior,
and to solve (18) for realizations of the posterior. With this approach we
avoid forming and solving the normal equations. Furthermore, this algorithm
supports matrix-free implementation of A as required. We use the ASTRA toolbox
to define the CT system and perform matrix-free, efficient, and GPU
accelerated forward- and backprojections [25, 26, 27].

If the CGLS algorithm is run until convergence when solving (18), we obtain
independent realizations of the posterior (15). However, in order to reduce
computation time, we apply a warm-start that uses the previous posterior
realization as a starting guess for the next one. We fix the number of
iterations to 10 per realization. This approach makes a burn-in phase
necessary, i.e. some initial samples are discarded. Furthermore, the
independence of the samples might be compromised, since we do not check for
convergence. Therefore, we monitor the integrated autocorrelation time (IACT)
of the samples [8, p. 99] to check for potential correlations. The idea is
that, if $\text{IACT}\approx 1$, the computed samples are approximately
independent.

###  5.3 Synthetic experiments

####  5.3.1 Example of Bayesian inversion with SGP prior

We consider a sparse-angle case with only 72 equi-angular projections out of
the original 360 (i.e. 20%). We employ the SGP-F prior in the Bayesian inverse
problem. First, we tune the GMRF precision parameter as described above in
Section 5.2.1. Figure 5 (left) illustrates the parameter sweep from which we
choose $\delta_{0}=1000$ since the RMSE is near a minimum. Then we compute
3000 samples from the posterior as described in Section 5.2.2 and discard the
first 1000 as burn-in. Since the reconstructed parameter x is high-
dimensional, we compute the IACT at 100 randomly selected pixel chains and
consider this as representative for the full parameter. Figure 5 (right) shows
that IACT is close to or below one in all cases, indicating that the samples
are quasi-independent.

Figure 6 shows the posterior mean, 95% interquantile range (defined as the
difference between the $2.5\%$ and $97.5\%$ quantiles), and 1D slices through
the mean with 95% credibility interval. We observe that the pipe layers are
well reconstructed. However, the contrast of the steel inclusions in the
concrete layer are dampened due to the choice of prior. The prior promotes the
expected linear attenuation coefficient of the concrete, which is lower than
that of steel. We also note higher noise in the concrete layer than in the
remaining parts of the reconstruction due to $\delta_{5}<\delta_{1,\ldots,4}$.
The 95% interquantile range is a direct reflection of the SGP prior action. We
see that imposing IID priors in the pipe layers lowers the uncertainty
estimate of the posterior compared to areas without IID priors on them.
Furthermore, priors with high precision values lead to lower uncertainty in
the posterior. As expected, we also see that the radial steel inclusions in
the concrete layer are more difficult to reconstruct than the tangential
inclusions. For instance, the 3 mm radial inclusion is barely distinguishable,
while the 3 mm tangential inclusion is fairly clear.

Figure 5: Left: Example of GMRF precision parameter sweep from the synthetic
experiment. From this, we choose $\delta_{0}=1000$ near the minimum. Right:
Examples of integrated autocorrelation time (IACT) for 100 randomly selected
pixel chains.

Figure 6: Synthetic test case with 20% view angles and SGP-F prior. Left, center: 2D images of the posterior mean and 95% interquantile range. Right: 1D slices of posterior mean (blue line) and 95% credibility interval (blue shade), compared to ground truth (orange). | Deterministic | GMRF | SGP-BG | SGP-F
---|---|---|---|---
360 view angles |  |  |  |
| (RMSE = 19.5$\times 10^{-3}$) | (RMSE = 14.8$\times 10^{-3}$) | (RMSE = 9.22$\times 10^{-3}$) | (RMSE = 9.14$\times 10^{-3}$)
50% view angles |  |  |  |
| (RMSE = 20.6$\times 10^{-3}$) | (RMSE = 16.7$\times 10^{-3}$) | (RMSE = 10.1$\times 10^{-3}$) | (RMSE = 9.91$\times 10^{-3}$)
20% view angles |  |  |  |
| (RMSE = 25.1$\times 10^{-3}$) | (RMSE = 16.5$\times 10^{-3}$) | (RMSE = 12.2$\times 10^{-3}$) | (RMSE = 11.6$\times 10^{-3}$)
10% view angles |  |  |  |
| (RMSE = 30.6$\times 10^{-3}$) | (RMSE = 18.9$\times 10^{-3}$) | (RMSE = 15.0$\times 10^{-3}$) | (RMSE = 12.5$\times 10^{-3}$)
Figure 7: Posterior means and RMSE with increasingly informative SGPs (in columns) and decreasing number of view angles (in rows). The color range is $[-0.05,0.20]$ for all images. | GMRF | SGP-BG | SGP-F
---|---|---|---
360 view angles |  |  |
| ($\delta_{0}=4000$) | ($\delta_{0}=4000$) | ($\delta_{0}=4000$)
50% view angles |  |  |
| ($\delta_{0}=3000$) | ($\delta_{0}=3000$) | ($\delta_{0}=3000$)
20% view angles |  |  |
| ($\delta_{0}=10000$) | ($\delta_{0}=3000$) | ($\delta_{0}=1000$)
10% view angles |  |  |
| ($\delta_{0}=10000$) | ($\delta_{0}=3000$) | ($\delta_{0}=1000$)
Figure 8: Posterior 95% interquantile range with increasingly informative SGPs
(in columns) and decreasing number of view angles (in rows).

####  5.3.2 Comparative study

We investigate the effect of the SGP when reducing the number of view angles.
Figure 7 shows posterior means on a grid, where the number of view angles
decreases by rows, and the prior information increases by columns. We observe
that, as we add more structural information, the reconstructions improve in
terms of 1) reduced streak artifacts, 2) increased layer contrast, 3) reduced
RMSE, and 4) detection of the steel inclusions that are hidden in artifacts.
The action of the prior is especially clear for 10% and 20% view angles,
indicating that the structural priors become more effective in sparse view
angle cases.

The deterministic reconstruction with 360 view angles produces an image where
all the main features of the pipe are visible. However, we also see artifacts:
1) along the radial steel inclusions, because these edges have no tangential
rays, 2) in the centre of the image, because there is no data, and 3) a dark
ring just outside of the pipe, which is a typical artifact that appears at the
fan-beam edge. As we decrease the number of view angles, we get more
dominating artifacts from the fan-beam edges. When we include the GMRF prior,
the artifacts are smoothed out, but they are still very dominating in the
image. Thus, this standard prior is not sufficient to produce good
reconstructions, especially for sparse view angle settings. Alternatively,
when using the SGP-BG prior, we manage to reduce the artifacts in the
background substantially. Furthermore, we see increased contrast between pipe
materials for the 10% and 20% view angle cases. This is due to the reduced
precision parameter of the smoothing GMRF prior, made possible by the added
information from the SGP-BG prior. When we use the SGP-F prior, we can
suppress artifacts in the background and pipe layers. In the concrete layer,
however, we note that unexpected inclusions (here steel reinforcements) that
are not modeled in the prior, are also removed. Recall, that we did not impose
any IID prior in boundary areas, and therefore the artifacts are still present
there. This we see as a ‘ragged’ boundary, especially for sparse view angle
cases.

We see that, with the SGP-F prior, we can resolve tangential steel inclusions
down to widths of 2 mm using only 20% of the view angles. The radial
inclusions are well reconstructed down to 4 mm. The smaller 2-3 mm radial
inclusions are still visible but faint and blurred. Compared to the
deterministic reconstruction with no prior, or to the reconstruction with a
GMRF prior alone, the results still show an improvement. Note that, in those
cases, the narrow radial steel inclusions disappear among the artifacts.
Moreover, for the 10% view angles case, it becomes difficult to reconstruct
any of the steel inclusions using any of the priors. Only the widest
tangential inclusions appear faintly in the reconstructions.

In Figure 8, we show the 95% interquantile ranges of the posteriors which
reflect the uncertainty in the reconstructions. Note that the colorbar for
each image is different. The GMRF prior yields a posterior uncertainty pattern
that reflects the acquisition geometry and that comes from the uncertainty of
the likelihood. For the SGP priors, the posterior uncertainty pattern reflects
the structure of the prior. For SGP-F, we observe lowest uncertainty in the
background and in the steel, PU, and PE layers (regions 1-4), where we find
the highest prior precision parameters. The concrete layer (region 5) has a
lower prior precision, and thus it reflects higher uncertainty. The highest
uncertainty is found at the boundary regions, where no particular attenuation
coefficient is promoted. We note that the boundary regions seem to have
relatively higher uncertainty for 10% and 20% view angles, than for 50% and
100%. This is caused by the lower GMRF precision $\delta_{0}$ used in the two
former cases. Finally, the likelihood structure is visible for the SGP-F prior
cases with 50% and 100% view angles, but not with 10% and 20%. We suspect this
is due to the prior having a relatively larger influence than the likelihood
when data is sparse.

| Deterministic | GMRF | SGP-BG | SGP-F
---|---|---|---|---
360 view angles |  |  |  |
50% view angles |  |  |  |
20% view angles |  |  |  |
10% view angles |  |  |  |
Figure 9: Posterior means with increasingly informative SGPs (in columns) and decreasing number of view angles (in rows). The color range is $[-0.05,0.20]$ for all images. | GMRF | SGP-BG | SGP-F
---|---|---|---
360 view angles |  |  |
| ($\delta_{0}=4000$) | ($\delta_{0}=4000$) | ($\delta_{0}=4000$)
50% view angles |  |  |
| ($\delta_{0}=3000$) | ($\delta_{0}=3000$) | ($\delta_{0}=3000$)
20% view angles |  |  |
| ($\delta_{0}=10000$) | ($\delta_{0}=3000$) | ($\delta_{0}=1000$)
10% view angles |  |  |
| ($\delta_{0}=10000$) | ($\delta_{0}=3000$) | ($\delta_{0}=1000$)
Figure 10: Real data: Posterior 95% interquantile range with increasingly
informative SGPs (in columns) and decreasing number of view angles (in rows).

###  5.4 Real data experiments

We now apply the SGPs to real data and obtain reconstructions with image
qualities that appear comparable to the synthetic experiments. As for the
synthetic data, we will study the effect of the SGP priors in combination with
sparse view angles data. Figure 9 shows the posterior means on a grid, where
the priors and number of view angles are varied. Our observations in the
synthetic data experiments are in good agreement with the real data case.
Again, we see dominating artifacts when decreasing the number of view angles
in the deterministic reconstruction with no prior information. We also see
that they are reduced by introducing structural information. In the 20% view
angles case, we observe that the SGPs not only reduce artifacts, but also help
increase contrast between materials. In the most extreme case, with only 10%
view angles, we lose a lot of details in the reconstruction, such as the
irregularities between the PU and PE layers. We do, however, still see traces
of the large steel inclusion in the concrete layer.

In Figure 10, we see 95% interquantile ranges associated with the posteriors.
Again, we observe that the acquisition geometry is reflected when using the
GMRF prior, and that the SGP structures are reflected (layers) when using
these particular priors. Note that the steel rods used to fixate the pipe
during measurements do not have any IID prior promoting a particular
attenuation coefficient.

##  6 Discussion and Conclusions

We proposed a method of introducing prior structural information in Bayesian
CT reconstruction. Structural Gaussian priors are constructed based on well-
known geometry and materials to promote expected attenuation coefficients in
CT reconstructions. We applied the SGP in a case study of CT imaging of a
subsea pipe with an offset fan-beam acquisition geometry. Numerical
experiments with real and synthetic data showed that SGPs improve
reconstructions in terms of suppressed artifacts and enhanced contrast between
materials, compared to traditional methods with no prior information or with a
standard smoothing GMRF prior. The SGPs were especially powerful in setups
with sparse view angles. Information about all layers of the pipe is required
to construct the SGP-F. Using this prior yielded the best results in the
numerical experiments. In some cases, we might only have access to the pipe’s
exterior dimensions. Then, the SGP-BG prior is applicable, and also in such
cases, the numerical results showed substantial improvements of the
reconstructions compared to using a GMRF prior or no prior.

Since the posterior is Gaussian and defined in terms of square-root precision
matrices, it can be efficiently sampled by solving a system of linear
equations, yielding approximately independent samples. This is a useful
property compared to conventional Markov chain Monte Carlo (MCMC) methods such
as Metropolis Hastings [28] or more modern Hamiltonian MCMC methods [29],
where independent samples cannot be guaranteed. To overcome this issue, MCMC
methods often require a ‘thinning’ process where many of the samples are
discarded. Therefore, MCMC methods require many more simulations, and thus,
much longer computation times to obtain a sufficiently large quasi-independent
sample size.

If uncertainty estimates are not particularly required, our method based on
SGP priors can also achieve reconstructions using MAP estimates that can be
computed extremely efficiently. However, since one of our main goals is
uncertainty quantification, we chose to sample the posterior.

Interpretation of the posterior uncertainty estimates is difficult and
requires more study. We observed that the GMRF prior gives rise to uncertainty
estimates that reflect the acquisition geometry, and this is in agreement with
the results in [30]. When imposing the SGP, we noticed, that the uncertainty
estimates directly reflected the prior. The prior’s effect on the uncertainty
estimates has also been seen in other studies. For instance, edge-preserving
priors give rise to high uncertainty on discontinuities in CT reconstructions
[5, 6]. In future work, we aim to further improve uncertainty estimates of the
pipe and defect structure, for example to help distinguish artifacts from real
features in the reconstruction.

## Acknowledgements

We would like to thank FORCE Technology for providing real data and
information regarding subsea pipes.

## Funding

This work was supported by The Villum Foundation (grant no. 25893).

## References

  * [1] Hansen PC, Jørgensen JS, Lionheart WRB, editors. Computed tomography: Algorithms, insight, and just enough theory. Philadelphia: SIAM; 2021.
  * [2] Hansen PC. Discrete inverse problems: Insight and algorithms. Society for Industrial and Applied Mathematics (SIAM); 2010.
  * [3] Riis NAB, Frøsig J, Dong Y, et al. Limited-data x-ray CT for underwater pipeline inspection. Inverse Problems. 2018 jan;34(3):034002.
  * [4] Kaipio J, Somersalo E. Statistical and Computational Inverse Problems. Dordrecht: Springer; 2005.
  * [5] Suuronen J, Emzir M, Lasanen S, et al. Enhancing industrial x-ray tomography by data-centric statistical methods. Data-Centric Engineering. 2020;1(3):e10.
  * [6] Uribe F, Bardsley J, Dong Y, et al. A hybrid Gibbs sampler for edge-preserving tomographic reconstruction with uncertain view angles; 2021. ArXiv:2104.06919.
  * [7] Sidky EY, Pan X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Physics in Medicine & Biology. 2008;53(17):4777.
  * [8] Bardsley JM. Uncertainty quantification for inverse problems. Philadelphia: SIAM; 2018.
  * [9] Bardsley JM. Gaussian Markov random field priors for inverse problems. Inverse Problems and Imaging. 2013;7(2):397–416.
  * [10] Bardsley JM. Laplace-distributed increments, the laplace prior, and edge-preserving regularization. Journal of Inverse and Ill-Posed Problems. 2012;20(3):271–285.
  * [11] Markkanen M, Roininen L, Huttunen JMJ, et al. Cauchy difference priors for edge-preserving Bayesian inversion. Journal of Inverse and Ill-posed Problems. 2019;27(2):225–240.
  * [12] Suuronen J, Chada NK, Roininen L. Cauchy Markov random field priors for Bayesian inversion. Statistics and Computing. 2022;32(33):1573–1375.
  * [13] Zhang H, Han H, Ma J, et al. Deriving adaptive MRF coefficients from previous normal-dose CT scan for low-dose image reconstruction via penalized weighted least-squares minimization. Medical physics. 2014 04;41:041916.
  * [14] Zhang H, Han H, Liang Z, et al. Extracting information from previous full-dose CT scan for knowledge-based Bayesian reconstruction of current low-dose CT images. IEEE Transactions on Medical Imaging. 2016 03;35:860–870.
  * [15] Han H, Zhang H, Wei X, et al. Texture-preserved penalized weighted least-squares reconstruction of low-dose CT image via image segmentation and high-order MRF modeling. In: Medical Imaging 2016: Physics of Medical Imaging; 03; 2016. p. 97834F.
  * [16] Gao Y, Liang Z, Moore W, et al. A feasibility study of extracting tissue textures from a previous full-dose CT database as prior knowledge for Bayesian reconstruction of current low-dose CT images. IEEE Transactions on Medical Imaging. 2019 01;38(8):1981–1992.
  * [17] Fukuda W, Maeda Si, Kanemura A, et al. Bayesian x-ray computed tomography using material class knowledge. In: 2010 IEEE International Conference on Acoustics, Speech and Signal Processing; 04; 2010. p. 2126–2129.
  * [18] Buzug T. Computed Tomography: From photon statistics to modern cone-beam CT. Berlin: Springer; 2008.
  * [19] Hubbell JH, Seltzer SM. Tables of x-ray mass attenuation coefficients and mass energy-absorption coefficients (version 1.4) ; 2004. [Online] Available: https://dx.doi.org/10.18434/T4D01F. Accessed: 2021-01-19.
  * [20] Rafiei MM, Tavakoli-Anbaran H. Calculation of the exposure buildup factors for X-ray photons with continuous energy spectrum using Monte Carlo code. Journal of Radiological Protection. 2018;38(1):207–217.
  * [21] Trubey DK. New gamma-ray buildup factor data for point kernel calculations: Ans-6. 4. 3 standard reference data. United States: International Nuclear Information System; 1988. Report ORNL/RSIC–49.
  * [22] Rue H, Held L. Gaussian Markov random fields theory and applications. New York: Chapman & Hall/CRC; 2005.
  * [23] Bickel P, Doksum K. Mathematical statistics: Basic ideas and selected topics. 2nd ed. Vol. 1. New Jersey: Prentice Hall; 2001.
  * [24] Björck A. Numerical methods for least squares problems. Philadelphia: SIAM; 1996\.
  * [25] van Aarle W, Palenstijn WJ, Cant J, et al. Fast and flexible x-ray tomography using the ASTRA toolbox. Opt Express. 2016 Oct;24(22):25129–25147.
  * [26] van Aarle W, Palenstijn WJ, De Beenhouwer J, et al. The astra toolbox: A platform for advanced algorithm development in electron tomography. Ultramicroscopy. 2015;157:35–47. Available from: https://www.sciencedirect.com/science/article/pii/S0304399115001060.
  * [27] Palenstijn WJ, Batenburg KJ, Sijbers J. Performance improvements for iterative electron tomography reconstruction using graphics processing units (GPUs). Journal of Structural Biology. 2011;176(2):250–253.
  * [28] Hastings WK. Monte carlo sampling methods using Markov chains and their applications. Biometrika. 1970;57:97–109.
  * [29] Duane S, Kennedy A, Pendleton BJ, et al. Hybrid monte carlo. Physics Letters B. 1987;195(2):216–222.
  * [30] Bardsley JM. MCMC-based image reconstruction with uncertainty quantification. SIAM Journal on Scientific Computing. 2012;34(3):A1316–A1332.