A Bayesian approach for CT reconstruction with defect detection for subsea pipelines

Subsea pipelines can be inspected via 2D cross-sectional x-ray computed tomography (CT). Traditional reconstruction methods produce an image of the pipe’s interior that can be post-processed for detection of possible defects. In this paper we propose a novel Bayesian CT reconstruction method with built-in defect detection. We decompose the reconstruction into a sum of two images; one containing the overall pipe structure, and one containing defects, and infer the images simultaneously in a Gibbs scheme. Our method requires that prior information about the two images is very distinct, i.e. the first image should contain the large-scale and layered pipe structure, and the second image should contain small, coherent defects. We demonstrate our methodology with numerical experiments using synthetic and real CT data from scans of subsea pipes in cases with full and limited data. Experiments demonstrate the effectiveness of the proposed method in various data settings, with reconstruction quality comparable to existing techniques, while also providing defect detection with uncertainty quantification.


Introduction
Subsea pipelines (figure 1) are critical components of offshore oil and gas infrastructure, and their proper functioning is vital to ensure the safe transport of hydrocarbons [1].These pipelines are exposed to a variety of environmental stresses and can suffer from various types of defects, such as cracks, corrosion, and deformation, which can lead to leaks and other catastrophic failures.Detecting these defects early is essential for pipeline maintenance and the prevention of costly downtime.However, manual inspection of subsea pipelines is challenging, time-consuming, and often requires the removal of outer layers, making it impractical and expensive.
To address these challenges, x-ray computed tomography (CT) has emerged as a promising technology for non-destructive testing of subsea pipelines [2][3][4][5].CT allows for the noninvasive inspection of the internal structure of pipelines by generating cross-sectional images, enabling the detection of defects with high accuracy.However, the detection and quantification of defects in CT images can be challenging due to various factors, such as image noise, artifacts, and limited resolution, which can lead to uncertainty in the detection of critical defects.
Traditionally CT problems are addressed by formulating and solving deterministic linear inverse problems [6,7].In that framework a single reconstruction is obtained with no information about its accuracy or uncertainty, which can arise from e.g.noise, missing data, or model errors.In contrast, by formulating and solving the CT problem in a probabilistic, Bayesian framework, however, uncertainty quantification (UQ) can be achieved [8,9].That is, the Bayesian framework allows us to infer a distribution of solutions to the inverse problem, from which point estimates with associated uncertainty can be computed.UQ in CT has been exemplified through a range of methodologies, from assumption-driven and data-driven Bayesian models using Markov-Chain Monte Carlo (MCMC) sampling [10][11][12][13][14] to iterative maximum A Posteriori (MAP)-based approaches [15][16][17][18], and Bayesian neural networks [19].
In this work, we propose a novel Bayesian CT reconstruction method with built-in defect detection, and apply it to detect defects in subsea pipelines.Our approach is based on decomposing the reconstruction into the sum of two images; one containing the large-scale pipe structure and the other containing small-scale defects.Our method offers two key advantages: First, it eliminates the need for advanced image analysis methods in a post-processing step (see e.g.[20]).Second, it enables us to separately formulate prior information for the large-scale pipe layers and small-scale defects, which in turn allows us to promote sparsity in the defect reconstruction and better capture the internal structure and materials of the pipes.
Our method also enables UQ of detected defects, which can be a critical aspect of the pipeline inspection.UQ helps decision makers determine the accuracy of a reconstructed defect by providing a quantitative measure of the uncertainty associated with the defect detection, ultimately reducing the risk of costly downtime and increasing the safety and reliability of offshore oil and gas infrastructure.The importance of UQ is well-established in the context of Bayesian approaches to inverse problems.The Bayesian framework is equipped with a powerful mathematical decision theory that allows the rigorous use of image data as evidence to inform decisions under uncertainty [21][22][23].
Our approach bears similarities to the widely applied 'sparse and low-rank matrix decomposition' method described in [24], in which an unknown matrix is decomposed into a sum of a sparse and a low-rank matrix.During reconstruction of the unknown matrix the problem is regularized by minimizing penalty terms consisting of the 1-norm of the sparse matrix and the nuclear norm of the low-rank matrix.This methodology has been applied in several areas including dynamic magnetic resonance imaging [25], synthetic aperture radar imaging [26], and spectral CT [27].In our application, we are also interested in one image being sparse (the defect image).The other image, however, should not be low-rank as in [24], rather it must capture the large-scale internal structure of a subsea pipe.Another key difference between our approach and [24] is that we employ a fully Bayesian formulation in which prior information is encoded as probability distributions and not as penalty terms in an optimization problem.This enables us to formulate and quantify the posterior distribution for both pipe and defect images through MCMC sampling, which is the key to quantify uncertainties.To our knowledge, this is the first time a Bayesian decomposition method has been employed in CT enabling reconstruction and quantification of uncertainties directly on defects.
To promote sparsity and pipe structure in the Bayesian setting we must choose appropriate priors for both images.The general internal structure and materials of the pipes are typically well known and this information can be used to formulate a prior for the pipe reconstruction [28].The other important piece of information we have, is that we expect few and small defects in the pipe, i.e. sparsity in the image.This makes a sparsity prior appropriate for the defect reconstruction.The authors in [29] identifies two main categories of sparsity promoting priors: (1) simple heavy-tailed distributions which include e.g.Laplace, student-t, and Cauchy, and (2) mixture models with hidden variables that give rise to hierarchical priors, which include e.g. the class of shrinkage priors [30] and the hidden gamma Markov random field (MRF) [31,32].Of these, the gamma MRF is particularly interesting for our application, because it can be encoded to promote correlation between neighbouring pixels, which means defects will appear as spatially coherent structures in the reconstructed image.

Overview and notation
This paper is organized such that the relevant CT and Bayesian background theory is summarized in section 2. The proposed methodology of Bayesian CT reconstruction with defect detection is presented in section 3.In section 4 we apply the methodology in test cases with synthetic and real data.We demonstrate how the methodology can successfully detect defects in cases with both complete and limited data.In section 5 we give concluding remarks and an outlook to future work.
For convenience, we define notation used in the paper: Random variable vectors are denoted with bold, lower case letters, e.g.b, and with elements b i .Realizations of random vectors are denoted with a descriptive text in a subscript, e.g.b obs .Random vectors rearranged into image arrays are denoted with a tilde, e.g. the random vector s ∈ R n 2 is equivalent to the array s ∈ R n×n .Matrices are denoted with bold, upper case letters, e.g.A, with elements A i,j , and sometimes with descriptive subscripts, e.g.R SGP .Scalars are denoted with roman or greek lower case letters and sometimes with descriptive subscripts, e.g.δ or δ 0 .Distributions are denoted; Gaussian: N (mean, covariance), gamma: G (shape, scale), and inverse gamma: IG (shape, scale).Probability density functions (pdf) are denoted by p(•), and likelihood functions are denoted by L(•).

The inverse problem of computed tomography
X-ray CT is a non-destructive imaging technique that utilizes different material's x-ray attenuation coefficients to create contrast.Measurement data in CT consists of projections and is acquired by letting an x-ray source with known intensity illuminate the object of interest and recording the attenuated x-ray intensity on the other side of the object.Projections are recorded at different angles around the object.By modelling attenuation using the Beer-Lambert law with line integrals [7], we can formulate a discretized linear inverse problem b = Au + e. ( Here, b ∈ R m represents the observed x-ray absorption (also known as a sinogram), e ∈ R m is additive noise, and u ∈ R n contains the unknown attenuation coefficients making up the image to be reconstructed with n = N × N pixels.A ∈ R n×m is the system matrix where the entry A i,j describes the intersection length of the i'th x-ray with the j'th pixel.For more details on the general CT inverse problem see [7], for example, for a motivation of the form of (1), which arises from the Beer-Lambert law of x-ray attenuation in matter combined with a postlogarithm Gaussian approximation to the general Poisson noise model for transmission CT.

Bayesian inversion
In the Bayesian approach to inverse problems, u and b are considered random variables and assigned a probability distribution.Given some observed realization of the data b obs , we perform Bayesian inference on u by forming the posterior distribution of solutions to the inverse problem according to Bayes' theorem: Here we use the notation L(u | b = b obs ) for the likelihood function, which expresses how well a realization of u explains the observed data [33, chapter 1.5.2], and it can be specified by assuming a distribution for the data noise.The prior distribution denoted by p(u) captures prior information about the unknown variable we might have.Generally, quantifying the prior information as a distribution is a challenging task, and ultimately, also a subjective choice as discussed in [8, chapter 3.3] and [34, chapter 3.1].We note that the form of Bayes' Theorem used here omits the normalization factor p(b) and writes the posterior as proportional to the product on the right hand side.This is a commonly used abbreviation, see e.g.[8], that does not affect the quantitative values of inferred parameters (the posterior distribution has unit probability mass by construction, so it is uniquely specified).We remark that in the context of CT, with high-intensity x-rays, it is reasonable to assume additive and independent Gaussian measurement noise, which leads to a Gaussian likelihood function [8, chapter 3.2.1].Also, priors commonly used in CT and other imaging problems include Gaussian and Laplace MRF's, which represent the Bayesian equivalents of Tikhonov and Total-Variation regularization respectively [35,36].
In practice, the posterior distribution is often impossible or impractical to form and evaluate.Instead, the posterior distribution can be approximated by sampling from it using Markov Chain Monte Carlo (MCMC) methods.A sample from the posterior consists of chains of realizations of the random variable u | b = b obs , and it can be used to estimate statistics (such as mean, standard deviation, and credibility intervals) characterizing the posterior distribution.To obtain the best estimates, a high number of independent realizations is necessary, especially for higher-order moments such as variance and thus standard deviation.However, MCMC methods often produce samples where consecutive realizations are correlated to each other.Therefore, we specify the concepts of sample size and effective sample size (ESS) [37, chapter 11.8].The sample size is defined as the number of realizations in the sample, and it carries no information regarding correlation in the chains.On the other hand, the ESS represents the number of independent realizations in the sample, and it is thus a more useful metric, when evaluating the sample's ability to reliably estimate posterior statistics.

Bayesian inversion with multiple unknown parameters
Sometimes, there is more than one unknown parameter in the Bayesian inverse problem.If we denote the unknown parameters u 1 , . . ., u K , we can form a joint posterior distribution where p(u 1 , . . ., u K ) is the joint prior distribution over parameters to be inferred u 1 , . . ., u K that can be of arbitrary finite dimensions given by the problem.Note, if the unknown variables are independent, the joint prior can be decomposed into a product of priors p(u 1 , . . ., u K ) = K k=1 p(u k ), and if the unknown variables are conditionally dependent on each other, the joint prior can be decomposed according to the chain rule of probability [38, chapter 5.1].Bayesian inverse problems that contain conditional dependencies between the unknown random variables are sometimes referred to as hierarchical Bayesian inverse problems (e.g. in [9, chapter 5.2]).
Sampling in the framework of hierarchical Bayesian inverse problems is often carried out using the Gibbs sampling algorithm [39].The general concept is that the joint posterior distribution can be approximated by sampling sequentially from the conditional distributions of the unknown random variables A property of Gibbs sampling is that each unknown variable's marginal posterior statistics can be computed based on the sample from its conditional distribution [40].For instance, if the posterior conditionals (4) have been sampled in a Gibbs scheme, the posterior statistics characterizing the unknown random variable u 1 can be computed purely based on the chain of realizations from p(u Finally, we remark that some combinations of distributions defining the posterior are easier to sample than others.If the product of distributions belong to the same family, the property of conjugacy [41, chapter 3] can be utilized to form closed-form expressions for the conditional distributions.Often this will substantially simplify numerical computations and lead to more efficient samplers in the sense of computation time to reach the same ESS.

Methodology
In this section, we describe our proposed methodology for defect detection in CT pipeline inspection using a Bayesian framework.We specify how we decompose the problem into pipe and defect parts, the Bayesian problem formulation including likelihood and pipe and defect priors, as well as efficient posterior sampling and analysis.

Decomposition into pipe and defect images
In order to include defect detection directly in the reconstruction method, we decompose the reconstructed image as the following sum, where we intend z to contain the expected large-scale structure of the pipe and d to contain any defects in the pipe.Inserting ( 5) into (1) leaves us with the following inverse problem: Solving (6) for z and d simultaneously without any prior information will lead to identifiability issues.This means that even if the sum z + d can be determined uniquely, its parts z and d can not, and hence the estimation problem is not well posed.To remedy this problem we introduce more information about the expected structures and properties of the unknown quantities involved.In the Bayesian framework this is achieved by constructing appropriate and informative priors for each variable.We formulate a prior on z that promotes expected large-scale pipe structure and reduces the probability of z representing small-scale defect structures.On the other hand, the prior on d is formulated such that it promotes few and small structures, and assigns low probability to large structures.This is intended to allow us to reconstruct z and d simultaneously and thus obtain defect detection as part of the reconstruction.Since u = z + d represents linear attenuation coefficients, one natural approach would be to enforce nonnegativity as part of the prior.As discussed in [42] this is not trivial in the Bayesian setting; approaches include using priors natively defined only the nonnegative orthant such as the Gamma distribution, using truncated distributions or so-called projected densities as proposed by [42], however these approaches are either expensive or not directly compatible with the priors considered in the present work.We leave nonnegativity for future work and focus on demonstrating utility with the proposed splitting-based approach.

Bayesian formulation
As described in section 2, we can obtain a posterior distribution of solutions to the CT inverse problem based on a likelihood function and prior information about the unknown variables.Below, we specify the likelihood function and prior distributions for z and d, and explain how the resulting posterior distribution can be sampled in a Gibbs scheme.
3.2.1.Likelihood.We assume that high-intensity x-rays are used for the data collection and therefore we model the additive data noise as being independent and identically distributed (IID) Gaussian with mean zero and variance λ −1 , i.e., e ∼ N 0, λ −1 I , (7) where I ∈ R n×n is an identity matrix.From this follows that the (joint) likelihood function for the inverse problem (6) given an observation of the data b obs is given by 3.2.2.Pipe prior.We assume some structural information about the pipe is known a-priori from design specifications or previous reference scans.Alternatively, the structural information can be obtained via a deterministic reconstruction as described in [28].Specifically, we assume the pipe consists of layers with approximately known layer boundaries, and that the materials in these layers are well known, such that we can obtain good a-priori estimates of their linear attenuation coefficients.The structural Gaussian prior (SGP) was proposed in [28] and it exploits exactly this information.Therefore, it is employed here as the pipe prior.The pipe prior specifies that z follows the Gaussian distribution: with The SGP is made of two main components, where the index k = 0 represents a global smoothing Gaussian MRF prior [35], and k = 1, . . ., p represents regions in the image where IID Gaussian priors with means µ k are used to promote expected attenuation coefficients according to known material properties.Each square-root precision matrix 5 has the form R k = √ δ k M k , where δ k is a scalar precision parameter that controls the strength of the prior, and M k contains the structure of the square-root precision matrix.For more details see [28].To illustrate the strong information encoded in the SGP, we show two realizations from (9) in figure 2, where the prior is specified according to the pipe described in section 4.2.Notice that the prior information is weaker around the layer boundaries, which is indicated by the higher variability between the two prior realizations there.The pipe prior pdf is: where |•| is the determinant.

Defect prior.
We seek a defect prior that promotes small and spatially coherent structures.The desired properties can be achieved with a hidden gamma MRF, that was introduced in [31] for modelling sparsity and local dependence between energies of time-frequency expansion coefficients in audio source modelling.In another study [32], the gamma MRF was applied as a prior for a nonlinear hyperspectral unmixing problem, where non-linearity coefficients were expected to be sparse and correlated.Not only does the gamma MRF have the desired prior properties, but it is given by conjugate distributions, which we will see in section 3.3 is convenient for sampling.
To simplify presentation of the gamma MRF for the defect prior, we rearrange the defect image random vector d into a pixel image represented in the random array d ∈ R N×N .We also introduce two hyperparameters; s ∈ R N×N and w ∈ R (N+1)×(N+1) .Note, these hyperparameters can also be rearranged into random vectors s ∈ R n and w ∈ R (N+1) 2 , which is more convenient when formulating the joint posterior distribution.
The gamma MRF is expressed in a hierarchical manner: where i = 1, . . ., N, j = 1, . . ., N, ω is a scalar parameter, and gi,j (w) = (w i,j + wi+1,j + wi,j+1 + wi+1,j+1 ) /4, We see here, that the defect image prior is IID Gaussian distributed with zero mean and locally varying variance given by the hyperparameter s.This prior pushes the pixel values of the defect image towards zero, but with locally varying 'strength', i.e. in areas where the variance is large, the prior is not very influential on the posterior, and this is where defects can occur.The hyperprior on s is the inverse gamma distribution, which is heavy-tailed for small ω.Thus, s most likely takes small values, but large values can occur.Since s represents the pixel-wise variance of d, this leads to sparsity in the defect image.The role of the hyperparameter w is to introduce positive correlation between neighbour elements in s, which promotes spatial regularity in d.Ultimately, the gamma MRF promotes a defect image that is both sparse and spatially correlated.
As described in [32], the prior's structure is defined from a bipartite conditional independence graph between s and w.The graph's vertex set V = V s ∪ V w can be partitioned into V s and V w , and the edge set E connects elements si,j to four neighbours in w and vice versa.The graph is illustrated in figure 3, where we also see that edges are associated with the coupling parameter ω, which acts as a regularization parameter in the prior controlling spatial smoothness.
The pdf describing the joint distribution of the unknown variables in the defect prior is given [32]: where Z(ω) is a normalization factor.When applying the gamma MRF prior, the unknown variables d, s, and w are all included as a part of the Bayesian inference.That means the only parameters left to specify are ω and boundary conditions for the gamma MRF.The defect prior parameter ω has a double role and affects both the spatial smoothness and the sparsity of the defect image.For small ω, high sparsity and low spatial smoothness is promoted.To compute hi,j (s) for i = 1, j = 1, i = N + 1, and j = N + 1, we implement boundary conditions in the gamma MRF such that s0,• = s•,0 = sN+1,• = s•,N+1 = s BC .Letting s BC be small leads to low variance at the boundaries of the defect image and thus essentially resulting in defect pixels that take values close to zero with very high probability.

Posterior sampling
With the likelihood, pipe prior, and hierarchical defect prior defined above, and assuming the pipe and defect priors are independent, we get the joint posterior distribution: , d, s, w) We sample from this distribution using the Gibbs scheme defined in algorithm 1.This requires, that we can sample from the conditional distributions for each of the four unknown random variables z, d, s, and w, which are given: Due to our choices of priors, the expressions above all have closed form solutions and can be sampled efficiently.Gaussian distributed because their pdf's are defined as products of Gaussians.We use linear RTO ('randomize-thenoptimize') [43] to sample efficiently from these high dimensional Gaussian distributions, as described in the following.A realization z * from z | d, b = b obs is obtained by solving the following system of linear equations [28]: Similarly, we can obtain a realization d * from d | z, s, b = b obs by solving: In both cases we use the conjugate gradient least squares algorithm (CGLS) [44, The conditional w | s is gamma distributed and given directly in the prior (12c).We can obtain realizations w * i directly by sampling from: We use linear RTO as an established sampler that is sufficiently efficient to handle images at the desired resolution of 500-by-500 pixels.Based on the iterative CGLS method, the method scales well to larger images.If required, one may explore further acceleration techniques such as chromatic Gibbs sampling [45] but this is beyond the scope of the present work.

Analysing posterior samples
In this section we address how the posterior samples are analysed and evaluated.We compute statistics of the marginal posterior distributions of the unknown parameters based on the Gibbs samples of the posterior obtained from algorithm 1. Specifically we compute, for each unknown parameter, the pixel-wise mean, which is a point estimate that represents the optimal mean-squared-error reconstruction [34, chapter 2.5.1], and the pixel-wise standard deviation (std), which provides an indication of the uncertainty related to the reconstruction at the scale of the individual pixel (we expect defects to be possibly as small as a pixel).
To ensure that the posterior statistics are estimated well, we compute the ESS of our samples.If the ESS is close to the actual sample size, the sample realizations are close to independent.
We compare our reconstruction method qualitatively to total-variation (TV) reconstruction [46] to evaluate the performance.In the deterministic, variational approach to CT reconstruction, TV is a widely used regularizer for edge-preservation in the image.We use the isotropic TV regularizer, and then the reconstruction is obtained by solving the optimization problem where D is a finite difference operator in the image's horizontal and vertical directions and ∥•∥ 2,1 is the mixed L 2,1 norm.TV reconstruction is carried out with the Core Imaging Library (CIL) [47,48], see these references for details.The optimization problem is solved using the fast iterative shrinkage-thresholding algorithm [49,50] provided in CIL and the regularization constant α is chosen by visual inspection of the image quality.

Numerical experiments
We conduct numerical experiments to demonstrate our methodology applied to defect detection in subsea pipes.The main goals of the experiments are to showcase; (1) how our proposed methodology provides automatic detection of defects in a subsea pipes, and (2) how the Bayesian approach allows uncertainty estimates.This section is organized as follows: In section 4.1 we give details about the numerical implementation of our methodology.In section 4.2 we introduce and describe the data used for the experiments.In section 4.3 we state how the problem's parameters are chosen.In section 4.4 we show results from numerical experiments with synthetic and real data.This includes a case using 'complete' synthetic data, where we give a detailed analysis of the posterior samples.Furthermore, we apply the methodology and show the resulting posterior statistics in more challenging limited data cases using synthetic data, and in a case with real CT data from a scan of a subsea pipe.

Numerical implementation
We use the python library CUQIpy [51][52][53] to specify and sample the Bayesian problem.While linear, our CT inverse problem is large-scale and to accelerate computations we do not form the forward operator matrix explicitly, but instead use the ASTRA library [54,55] for efficient, matrix-free and GPU-accelerated forward and adjoint operations.We observe issues with numerical stability when sampling from the posterior as described above.Realizations of the defect variance s i , are occasionally very close to zero with values down to machine precision.We found that setting a lower bound s LB on s i while sampling improves the stability.This means realizations of s i < s LB are forced to take the value s LB .
The python code used to produce the results in this work is openly available in the GitHub repository [56,57].

Real data description
We present here the practical case that has motivated this work.FORCE Technology has provided a CT dataset from a scan of a real subsea oil pipe and it is openly available from [58].A photograph of the pipe was shown in figure 1.The pipe's diameter is 460 mm and it has four layers of different materials that from outermost to innermost layers are: concrete, polyethylene (PE) rubber, polyurethane (PU) foam, and steel.In the concrete layer we also  see steel reinforcement bars.The CT data was acquired in a laboratory setting and projections were recorded at 360 equidistant angles over 360 • with a 410 mm detector of 510 detector elements.An overview of the acquisition geometry parameters can be found in table 1. Due to the pipe's large size compared to the detector size, a special offset fan-beam geometry was used as shown in figure 4. Compared to the standard centred fan beam geometry, this geometry allows resolving defects in the outer layers of the pipe better [5].The observed data b obs obtained from the CT scan of the pipe is shown in figure 5(a) as a sinogram.We note that the proposed method is not tied to or specially motivated by an offset scan geometry.This case was simply considered here due to availability of the real-data pipe data set, which was acquired in the offset geometry of figure 4. The method itself can be applied equally well in conventional and other scan geometries.

Synthetic data.
A synthetic phantom was constructed for simulation experiments based on the well-known internal structure and materials of the real subsea pipe.See e.g.[6] for a description of how linear attenuation coefficients can be computed for general materials, or see [28] for coefficients specific to our case.Inclusions representing defects are added to the phantom.The inclusions are designed such that we can investigate effects of the shape, direction, width, contrast, and crossing of layer boundaries.The phantom is shown in figure 5(c) and each defect is numbered for easy referencing.Note that the phantom is defined on a 1024 × 1024 pixel image, while reconstructions will be performed on a 500 × 500 pixel image to avoid inverse crime [59].The image domain covers 55 × 55 cm which means that the reconstruction pixel size is 1.1 × 1.1 mm.
We simulate a synthetic dataset by using the pipe phantom u, offset acquisition geometry, and additive Gaussian noise with a relative power of 2% (i.e.∥e∥ 2 /∥Au∥ 2 ≈ 2%).Comparing the synthetic and the real data sinograms in figure 5 we see similar thick vertical bands from the pipe layers.The bands from left to right represent air, concrete, PE rubber, PU foam, steel, and air again.This corresponds to the real pipe structure from outside to inside.We note that the pipe layers follow sinusoidal curves in the real data sinogram, which is due to the fact that the pipe centre and rotation centre did not coincide exactly, and the phase shifts between the curves are due to that fact that the pipe layers are not completely concentric.In both the synthetic and real data sinograms we see defects show up as thin, high amplitude, half sinusoidal curves.The similarity between the sinograms indicate that the synthetic phantom models the pipe well.

Parameter tuning
There are several parameters controlling the performance of the proposed methodology and these must be tuned.Below we describe our parameter choices.However, parameter tuning is not the main focus of this work.This means the parameters choices are good, but not necessarily optimal.Therefore, further work with parameter tuning might lead to improved results.
The noise precision parameter λ is fixed to match the expected data noise level of 2% which leads to λ ≈ 400.The parameters controlling the priors are tuned manually and the specific values used in the numerical experiments can be found in table 2. The pipe prior precision parameters δ 0 , . . ., δ p represent our confidence in the SGP.In the framework of this study, the pipe prior precision parameters must be set to large values such that z cannot represent defects in the pipes.This helps remedy ambiguity between z and d.The defect prior parameter ω acts as a regularization parameter that controls the heaviness of the tails of the defect prior distribution, i.e. it controls the sparsity of the image.The parameter s BC defines the boundary conditions in the gamma MRF and should be set to a low number given that we expect near zero values on the image boundary of the defect reconstruction.Finally, for numerical stability, we set a lower bound s LB on the realizations from s, as described in section 4.1.This value is also given in table 2. In general, we found that a fairly broad range of values for δ 0 , . . ., δ p have a similar effect, whereas the results are more sensitive to the values of ω, s BC and s LB .In all cases we obtain 10 000 realizations from the posterior as described in section 3.3 and we discard the first 2000 as burn-in.Total computation time was 9059 s, 8762 s, 7783 s and 11 475 s for Case 1 through 4, i.e. around 1 s per realization, with the amount of data being the biggest factor affecting computation time.We compare our reconstructions qualitatively to TV reconstructions as described in section 3. 4.
We note that the comparison with TV reconstruction provides a reference point in terms of obtainable quality for the considered data sets.The comparison also highlights the additional information provided by the proposed method: direct reconstruction of a defect image without a subsequent segmentation step as required by e.g.TV-reconstruction, as well as UQ for both pipe and defect images.
4.4.1.Case 1: detailed analysis of the posterior samples using complete, synthetic data.In this test case, we demonstrate our methodology using complete, synthetic data, and we give a detailed analysis of the results.
With our methodology we sample from the joint posterior distribution of four random vari- ables (z, d, s, and w).To visualize the posterior, we report the pixel-wise means and standard deviations of the marginal posterior distributions in figure 6.The sum of the reconstructed pipe and defect images is discussed later and shown in figure 9. Generally, we note that the posterior marginal distributions seem to represent the intended structures introduced by the priors.That is, the mean of z contains layered pipe structure, the mean of d contains small structures that we interpret as defects, and the mean of the hyperparameters s and w mirror the structure of the d standard deviation, that they define.
To confirm the validity of the posterior statistics presented above, we assess the ESS of our sample as discussed in section 3.4.We give pixelwise ESS values in figures 7(a)-(d) for each unknown parameter in the vicinity of defect 14 (reconstruction shown in figure 7(e) for reference).This area is deemed representative of the entire image, as it comprises regions with positive, negative, and no defects.The ESS values for the pipe image are consistently high, approaching the actual sample size (8000) for all pixels.Similarly, the ESS values for pixels without defects in the defect image are also high, at approximately 7000.However, for defect  pixels, the ESS is lower with values near 100.Finally, the ESS values for the hyperparameters are generally larger than 500 and take values up to 5000 on defect pixels.Generally, we consider these sample sizes adequate to obtain reasonable estimates of the posterior statistics.Now, we give a more detailed interpretation of the results from the posterior sampling.First of all, figures 6(a) and (b) indicate that the method is able to decompose the reconstruction such that the posterior means of z and d contain the pipe structure and defects respectively.We observe that when using the complete dataset, we are able to reconstruct almost all defects.Furthermore, most of the detected defects coincide closely with the ground truth (figure 5(c)).
We also observe some limitations.Figures 6(a) and (b) show that defects 10, 11 and 12, that cross layer boundaries, tend to be partly picked up by z instead of d.This is due to the fact that the SGP does not promote any attenuation coefficient near layer boundaries, and it is therefore possible for z to represent defects here.A stronger prior on the pipe image should remedy this issue (see section 4.4 for more details).Furthermore, the defect reconstructions tend to overestimate the pixel values slightly.We expect this is related to the choice of the ω parameter that affect the defect prior.However, for the task of defect detection, this value of the attenuation is not of big concern.Defect 7 and partly defects 11 and 12 can not be reconstructed because the contrast is too low.The signal from these defects is close to the noise level in the data, and thus, extremely difficult to reconstruct.However, this is not specific to our method, and for instance, we also see that it is an issue for the TV reconstruction shown in figure The standard deviation of the posterior images shown in figure 6 indicate the uncertainty associated with the reconstructions.We observe that the pipe standard deviation reflect the masks of the SGP.That is, it is low in regions where a particular absorption coefficient is promoted, and higher in the layer boundary areas.For the defect reconstruction, we observe that the standard deviation is low in the background and increase roughly an order of magnitude on the defects.The standard deviation is also slightly increased compared to the background in a 'halo' around the defects (this is seen more clearly in figure 8, that we discuss next).This might be related to the fact that correlation between the variances of neighbouring d pixels are enforced by the gamma MRF prior.Thus, we see smoothing of the transition between high and low values in the d posterior standard deviation.
We inspect a few interesting defects more closely in figure 8. Defects 4 and 16 are particularly difficult to reconstruct.Defect 4 is only 1 pixel wide and radially oriented and defect 16 is very small.The method is able to detect both defects, with associated uncertainties comparable to other defects.Defects 14-15 are round and with holes in the middle.In fact, defect 14 consists of two defects; a positive valued defect enclosed in a negative valued defect.Defects 14 and 15 are both reconstructed well.Finally, we study the reconstruction of defect 13.We see that it is difficult to accurately detect the concave corners.This is also reflected in the standard deviation, where we see that the shape of the reconstruction is more uncertain compared to the other defects.
Finally, in figure 9, we compare the sum of our pipe and defect reconstructions to a TV reconstruction of the same data.By visual inspection we note that the reconstructions are of similar quality.However, our method has separate defect reconstruction as well as the possibility of uncertainty quantification.This may aid with simple defect detection and quantification as well as risk assessment and decision making regarding the condition of the pipe.Our method is still able to decompose the reconstruction into pipe structure and defects, and most of the defects are detected.Compared to the TV reconstruction, our method is missing a few of the more difficult defects (6,12,16).On the other hand, our method does not smooth the edges of the defects as TV does.Another challenging test case is the limited-angle scan geometry.In this case we have no measurement data in large regions of the image, due to the off-set scan geometry, and standard reconstruction methods give images corrupted by artefacts.In figure 11(f), we for instance see that the TV reconstruction has limitedangle artefacts and staircase artefacts.The results using our method are shown in figures 11(a)-(e).We are again able to decompose the pipe and the defects and we can generate a reconstruction nearly free of artefacts, due to the information from the pipe prior.We see that only defects in a certain region are detected.This result is consistent with the fact that no data was measured outside this region.Finally, we observe that defect 9 and partly defect 13 are not detected even though they are within the scanned area.This is due to the direction of their long edges, which have no tangential x-rays in the data (see i.e. [7, chapter 8.3]).
In general for the synthetic data, we find that our method is able to obtain information about defects in areas with sufficient measurement data while at the same time avoid introducing artefacts into the reconstructed image.
4.4.4.Case 4: real data.We apply our method to real CT data from a scanned pipe with 360 equidistant view angles in a full rotation (see description of the data in section 4.2).The results are shown in figure 12. Like for the synthetic data case, our method is able to decompose the pipe and defect reconstructions.We obtain a detailed defect image that for instance contains steel inclusions, holes/irregularities between the PU foam and PE rubber layers, and the nonhomogeneous structure of the concrete.In figure 13 we show zooms of some of these defects.Defects 1 and 3 are large and small steel inclusions, defect 2 is a hole/irregularity between two pipe layers, and defect 4 shows a positive valued defect embedded in a negative valued defect.Comparing our methodology to TV reconstruction (figures 12(c) and (f)), we find that we are able to eliminate or reduce many artefacts.We completely eliminate the wavy artefacts in the background.There are five steel rods used for holding the pipe in place during data acquisition, and especially the four external rods cast shadows, in terms of streak artefacts, in the TV reconstruction.With our methodology we are able to reduce these compared to TV.The centre of the pipe is empty/air but we see significant streak artefacts here.It is not completely clear what causes these artefacts, but similar streaks occur in both TV and our method.We suspect model errors in the scan geometry or errors in the x-ray absorption model for steel and leave it to future work to study and reduce the centre artefacts.
To evaluate the uncertainty of the defect reconstruction we consider the standard deviation of the posterior shown in figure 12(e).We observe that the uncertainty is generally higher than for the synthetic data experiments, especially in the background.This might be due to a higher noise level (which is unknown) or modelling error.As for the synthetic data, the standard deviation also here indicate that the defect reconstruction is more certain about areas with no defects, than areas with defects.

Conclusions and future work
We proposed a novel Bayesian methodology for detecting defects in CT scans of subsea pipes.By decomposing the reconstruction into two images and using priors that promote specific structures in each image, we were able to successfully reconstruct and detect defects in both real and synthetic CT scans.For our application, this was obtained by using a SGP prior for promoting large-scale pipe structure in one image, and a gamma MRF prior for promoting sparse and coherent defects in the other image.The hierarchical Bayesian problem arising from these choices was efficiently solved using a Gibbs sampling scheme with the Python library CUQIpy, where each conditional was sampled by exploiting conjugacy.
The numerical experiments demonstrated that our methodology can successfully decompose the reconstruction and detect defects in real and synthetic CT scans of subsea pipes in the 'complete' data setting.The reconstruction quality of our method was comparable to TV.Furthermore, we obtained uncertainty estimates that suggested the defect reconstruction was most uncertain at the detected defects.
Our methodology also proved effective in more challenging settings, such as full rotation with few view angles and limited-angle rotation, where artefacts were reduced substantially.However, limitations of our approach still exists, such as the identifiability issues in the boundary regions between pipe layers.In these regions the SGP and the gamma MRF priors were too similar, which meant that defects crossing these areas were reconstructed in the pipe image instead of in the defect image.To address this issue, future work should involve the construction of a new prior for the pipe image that eliminates any possibility of identifiability issues, e.g. by parameterizing the pipe with a basis that can never represent a defect type structure.
Another limitation of the methodology is the manual parameter tuning necessary to obtain good results.We must tune 10 parameters: the noise precision, pipe prior precisions, the defect prior regularization parameter, and the boundary condition value and lower bound on the defect prior hyperparameter s.Future work will benefit from implementing automatic parameter tuning methods, for instance using maximum marginal likelihood estimation of the hyperparameters [60,61].Automatic parameter tuning will be a significant step toward practical application of the proposed method in engineering practice.

Figure 2 .
Figure 2. Two realizations from the SGP prior on z for synthetic data with zooms showing higher variability near layer boundaries.

Figure 3 .
Figure 3. Dependency structure between the hyperparameters s and w.
The conditionals z | d, b = b obs and d | z, s, b = b obs are both

Figure 4 .
Figure 4. Pipe and offset scan geometry specifications; dimensions in mm.The illustration shows the scan geometry at view angle 0 • .The rotation axis is marked with an orange star, and the rotation direction is marked with an orange arrow.

Figure 5 .
Figure 5. (a) Real subsea pipe sinogram.(b) Synthetic sinogram corrupted with additive Gaussian noise with a relative power of 2%.(c) Pipe phantom modelling the real subsea pipe.Defects are annotated with numbers for easy referencing.

Figure 6 .
Figure 6.Statistics of the posterior for the test case using complete synthetic data.The first row shows posterior means, and the second row posterior standard deviations.The columns from left to right give statistics for the inferred parameters z, d, s, w.All images have sizes 500 × 500 pixels.

Figure 7 .
Figure 7. (a)-(d) Pixelwise effective sample size (ESS) of each unknown parameter in area around defect 14.(e) Mean of z + d in same area for reference.Results using complete synthetic data.All images have sizes 26 × 26 pixels.

Figure 8 .
Figure 8. Zoom on posterior mean and standard deviation of the defect sample.All images have sizes 50 × 50 pixels, and are centred at the defects indicated by the column titles.Defect numbers are defined in figure 5(c).Results are obtained using complete synthetic data.

Figure 9 .
Figure 9. Reconstructions using complete synthetic data.(a) Posterior mean of z + d.(b) TV reconstruction of the same data.Both images have sizes 500 × 500 pixels.

Figure 10 .
Figure 10.(a)-(e) Posterior mean and standard deviation images for the test case using few view angles and synthetic data.(f) TV reconstruction of the same data.All images have sizes 500 × 500 pixels.

Figure 11 .
Figure 11.(a)-(e) Posterior mean and standard deviation images for the test case using limited view angles and synthetic data.(f) TV reconstruction of the same data.All images have sizes 500 × 500 pixels.

Figure 12 .
Figure 12. (a)-(e) Posterior mean and standard deviations for the test case using real complete data.(f) TV reconstruction of the same data.Images of sizes 500 × 500 pixels.

Figure 13 .
Figure 13.Zoom on posterior mean and standard deviation of the defect sample.All images are centered at the defects indicated by the column titles.Defect numbers are defined in figure 12(c).Defect 1 images have sizes 80 × 80 pixels and defect 2, 3, and 4 images have sizes 50 × 50 pixels.Results are obtained using complete real data.
is inverse gamma and conjugate to the Gaussian p(d i | s i ) with variance s i .The conjugate relationship leads to an inverse gamma distribution for s i | d i from which we can obtain realizations s chapter 7]to solve the linear system.The conditional s | d, w is inverse gamma distributed since p(s i | w) * i directly:

Table 1 .
Parameters defining the reconstruction geometry and acquisition geometry for a complete scan.The scan geometry is illustrated in figure4.

Table 2 .
Manually optimized parameter values used in numerical experiments.We demonstrate the proposed methodology in several cases below.With synthetic data, we study the results from three common measurement setups.Case 1: Complete dataset with 360 equidistant view angles in a full rotation.Case 2: Few view angles dataset where the sinogram has been subsampled to only contain every fifth view angle, i.e. 72 equidistant view angles in a full rotation.Case 3: Limited view angle dataset where the sinogram has been cropped to only contain the view angles 15 • -105 • , i.e. 90 view angles spread on a 90 • span.Finally, we present the results found using real CT data in case 4.