Li Wang

Ali Mohammad-Djafari

Nicolas Gac

Mircea Dumitru

3D X-ray Computed Tomography reconstruction using sparsity enforcing Hierarchical Model based on Haar Transformation

Keywords: Computed Tomography (CT), Bayesian Approach, Hierarchical Model, Generalized Student-t distribution, Joint Maximum A Posterior (JMAP)

In this paper, we consider the 3D X-ray CT reconstruction problem by using the Bayesian approach with a hierarchical prior model. A generalized Student-t distributed prior model is used to enforce the sparse structure of the multilevel Haar Transformation of the image. Comparisons with some state of the art methods are presented, showing that the proposed method gives more accurate reconstruction results and a faster convergence. Simulation results are also provided to show the effectiveness of the proposed hierarchical model for a reconstruction with more limited projections.

I. INTRODUCTION

An introduction of the development of Tomography and several analytic and algebraic state of the art reconstruction methods are presented in [START_REF] Clackdoyle | Tomographic reconstruction in the 21st century[END_REF]. In X-ray CT, the intensity of X radiation is attenuated when passing through the object, and the parameters to be reconstructed are the linear attenuation coefficients inside the object under the test. The Radon Transform (RT), presented in detail in [START_REF] Clackdoyle | Tomographic reconstruction in the 21st century[END_REF], is one of the most commonly used forward modeling when treating the 3D X-ray CT projection problem, with the expression:

g(r, φ) = Rf (x, y) = L r,φ f (x, y) dl (1)
where f (x, y) represents the attenuation coefficient, r is the perpendicular length from center point of coordinate and φ is the considered X ray angle. L r,φ is the length of ray (r, φ) passing through the object.

There have been different analytic methods to solve the reconstruction problem, for example the Back-Projection (BP) [START_REF] Stanley R Deans | The Radon transform and some of its applications[END_REF], the Filtered Back-Projection (FBP) [START_REF] Stanley R Deans | The Radon transform and some of its applications[END_REF], [START_REF] Chetih | Tomographic image reconstruction using Filtered Back Projection (FBP) and Algebraic Reconstruction Technique (ART)[END_REF], etc. We may mention here the FBP method, which has a very good performance when dealing with sufficient detected data. It can be summarized as f = BF -1 |Ω| F g where F and F -1 represent the direct and the inverse Fourier Transform (FT), |Ω| a modulated filter and B the Back Projection.

A. Deterministic Regularization based methods

Comparing with the conventional FBP method, the regularization and statistical methods give more satisfactory reconstruction results [START_REF] Chetih | Tomographic image reconstruction using Filtered Back Projection (FBP) and Algebraic Reconstruction Technique (ART)[END_REF]- [START_REF] Jang | Regularization designs for uniform spatial resolution and noise properties in statistical image reconstruction for 3-D X-ray CT[END_REF]. The statistical methods can also give more robust and precise results comparing with FBP reconstruction method, especially in the low dose and time limit situations where insufficient data are given.

In many algebraic and statistical methods, one consider the linear forward model:

g = Hf + (2)
where g ∈ R M ×1 represents the projection data, f ∈ R N ×1 the object and ∈ R M ×1 the additive noise. Matrix H ∈ R M ×N corresponds to the linear projection system. In the regularization methods, the result of reconstruction is obtained by optimizing the criterion which combines a data-model term and the regularization terms.

A general criterion of the regularization method is J (f ) = g -Hf 2 + λR(f ) where R(f ) is the regularization term and λ is called the regularization parameter. We list out some conventional regularization reconstruction methods: (a) the Least Square (LS) method with R(f ) = 0 where the estimating result minimize the disparity but is not guaranteed to be precise because of the ill-posedness; (b) the Quadratic Regularization (QR) method [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF] with R(f ) = Df 2 which enforces the global smoothness of the estimated result and (c) the Total Variation (TV) method [START_REF] Goldstein | The Split Bregman method for L1regularized problems[END_REF], where the difference of neighbour pixels are enforced to be sparse with R(f ) = Df 1 . By using l 1 norm, the sparsity of the penalty term is enforced. The appearance of the non-differentiable l 1 term leads to difficulties for the computation of the gradient of constraints. Many methods have been studied in order to solve this l 1 norm optimization problem, for example the Newton's method [START_REF] Tony F Chan | A nonlinear primal-dual method for Total Variation-based image restoration[END_REF] and the Split Bregman method [START_REF] Goldstein | The Split Bregman method for L1regularized problems[END_REF].

Note that in all the above mentioned regularization methods, there is a parameter λ controlling the trade-off of the disparity and the regularization penalty, which need to be fixed. There have been different methods to choose a suitable value for λ, for example the Cross Validation (CV) and L-curve methods, with details presented in [START_REF] Ramani | Regularization parameter selection for nonlinear iterative image restoration and MRI reconstruction using GCV and SURE-based methods[END_REF], [START_REF] Nikolas | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF]. However, the computation for choosing this parameter should be done for each simulation and different dataset, which is not practical. Bayesian methods, therefore, are often used to estimate the parameters and variables simultaneously.

B. Bayesian Method

As mentioned above, the Bayesian methods solves the reconstruction problem as well as the estimation of the parameters at the same time thanks to appropriate prior modelling. In short, by using the Bayes formula, the complete solution of the inverse problem is provided by the posterior distribution:

p(f , θ|g) = p(g|f , θ)p(f |θ)p(θ) p(g) , (3)
Which is then used to infer both on f and on θ.

II. THE SPARSITY ENFORCING HIERARCHICAL HAAR TRANSFORM BASED METHOD (HHBM)

The HHBM method presented in this paper is an extension of the work presented in [START_REF] Wang | Computed tomography reconstruction based on a hierarchical model and variational Bayesian method[END_REF]. Instead of estimating only the hidden variable z and obtaining the final reconstruction result f by a post-processing, here a more general model is proposed, and the method is tested on 3D object reconstruction simulations.

Typically the sparse property of the Haar transformation coefficient is enforced by using three kinds of distributions: the Generalized Gaussian distributions, the Gaussian Mixture distributions and the heavy tailed distributions.

The standard Student-t distribution is heavy tailed, but from the definition of its variance we easily figure out that there is a lower limit of the variance value of Student-t distribution:

Var[f ] = ν ν-2 > 1, (ν > 2)
. This limit implies that this heavy-tailed distribution can't have a small variance, therefore the sparsity couldn't be intensively enforced. In this paper we use a generalization of Student-t distribution (Stg) which can be obtained by using the marginalization of the Normal-Inverse Gamma bi-variate distribution:

St g (f |α, β) = N (f |0, v)IG(v|α, β) dv. (4)
This generalization of Student-t distribution adds a supplementary parameter compared to the standard one, and hence more capable to control the level of sparsity of the prior distribution.

A. Hierarchical Bayesian Model of HHBM

The additive noise is defined belonging to an i.i.d. Gaussian law with mean equal to zero and unknown variance vector v :

p( |v ) = N ( |0, V ), where V = diag [v ]
is a diagonal matrix. According to the linear forward model shown in Eq.( 2), we obtain an expression for the likelihood:

p(g|f , v ) = N (g|Hf , V ).
In the considered applications, generally the object to be reconstructed is piece-wise homogeneous, considering that it consists of several different materials. With this property, an information which can be considered as a prior knowledge is the sparseness of the contours.

The sparse representation of the prior piece-wise continuous image used in this paper is the multilevel Haar Transformation. A vector z is used to present the l-level Haar transformation coefficients of f . Prior distribution of f depends on z:

p(f |z, v ξ ) = N (f |Dz, V ξ )
where D represents the inverse Multilevel Haar transformation operator, and the additive noise ξ is considered to be i.i.d. Gaussian distributed. Definition of diagonal matrix:

V ξ = diag [v ξ ]. The vector z = [z 1 , z 2 , • • • , z N ],
therefore, is sparse. As mentioned above, the generalized Student-t distribution is used in order to enforce the sparsity. By using the Normal-Inverse Gamma property given in Eq.( 4), the prior distribution for z is:

p(z|v z ) = N (z|0, V z ) where V z = diag [v z ] , p(v z |α z0 , β z0 ) = N j IG(v zj |α z0 , β z0 ), (5)
where the elements of vector v z are supposed to be i.i.d.

On the other hand, when considering the variance of two noises: v and v ξ , by knowing that the variance is positive, and the fact that the majority of the values are small, we choose the Inverse Gamma distribution to model them:

p(v |α 0 , β 0 ) = M i IG(v i |α 0 , β 0 ) and p(v ξ |α ξ0 , β ξ0 ) = N j IG(v ξj |α ξ0 , β ξ0 ).
With all the proposed prior distributions, the model containing all the variables, parameters and hyper-parameters is:

                         p(g|f , v ) ∝ |V | -1 2 exp -1 2 (g -Hf ) T V -1 (g -Hf ) , p(f |z, v ξ ) ∝ |V ξ | -1 2 exp -1 2 (f -Dz) T V -1 ξ (f -Dz) , p(z|vz) ∝ |V z | -1 2 exp -1 2 z T V -1 z z , p(vz|αz 0 , βz 0 ) ∝ N j v -(αz 0 +1) z j exp -βz 0 v -1 z j , p(v |α 0 , β 0 ) ∝ M i v -(α 0 +1) i exp -β 0 v -1 i , p(v ξ |α ξ 0 , β ξ 0 ) ∝ N j v -(α ξ 0 +1) ξ j exp -β ξ 0 v -1 ξ j . (6)
The corresponding directed acyclic graph (DAG) of the proposed model is shown in Fig. [START_REF] Clackdoyle | Tomographic reconstruction in the 21st century[END_REF].

α z0 , β z0 c v z p(vz ) = IG(vz j |αz 0 , βz 0 ) c z p(z|vz ) = N (z|0, V z ) c D α 0 , β 0 c v p(v ) = IG(v i |α 0 , β 0 ) © p(g|f , v ) = N (g|Hf , V ) α ξ0 , β ξ0 c v ξ p(v ξ ) = IG(v ξ j |α ξ 0 , β ξ 0 ) d d d p(f |z, v ξ ) = N (f |Dz, V ξ ) f c

H g

Fig. 1: DAG of proposed model.

B. Bayesian Inference

Via the Bayesian inference, the posterior law is obtained from the likelihood and priors:

p(f , z, v , v ξ , v z |g) ∝ p(g|f , v )p(f |z, v ξ )p(z|v z ) • p(v )p(v ξ )p(v z ) (7)
From the posterior distribution obtained, different estimation methods can be used. Mainly there are two options: Posterior Mean (PM) and Joint Maximum A Posterior (JMAP) [START_REF] Qi | Resolution and noise properties of MAP reconstruction for fully 3-D PET[END_REF]. The first one can be computed either by MCMC methods [START_REF] Walter R Gilks | Markov chain monte carlo[END_REF], Variational Bayesian Approximation (VBA) [START_REF] Dimitris G Tzikas | The variational approximation for Bayesian inference[END_REF] or any other approximation methods. However, the computational costs of these methods are too high for 3D applications. This is the reason why we use the JMAP method to estimate all the variables iteratively in this paper.

C. Estimation

The JMAP computation aims at iteratively and alternately estimate the variables and parameters by maximizing the posterior: Because of the huge size of data, the descend gradient algorithm is used since the computation of big size matrix inversion is too expensive. The iterative updating rule is:

( f , z, v , v ξ , v z ) = arg max f ,z,v ,v ξ ,vz {p(f , z, v , v ξ , v z |g)} (8)
iter : f (k+1) = f (k) -γ (k) f ∇J ( f (k)
);

iter :

z (k+1) = z (k) -γ (k) z ∇J ( z (k) )
;

vz j = βz 0 + 1 2 z 2 j / (αz 0 + 3/2) , ∀j ∈ [1, N ]; v i = β 0 + 1 2 gi -H f i 2 / (α 0 + 3/2) , ∀i ∈ [1, M ]; v ξ j = β ξ 0 + 1 2 fj -[D z] j 2 / (α ξ 0 + 3/2) , ∀j ∈ [1, N ], (9)
where

J (f ) = 1 2 V -1 2 (g -Hf ) 2 + 1 2 V -1 2 ξ (f -Dz) 2 ; J (z) = 1 2 V -1 2 ξ (f -Dz) 2 + 1 2 V -1 2 z z 2 ; γ (k) f = ∇J ( f (k) ) 2 / Y H∇J ( f (k) ) 2 + Y ξ ∇J ( f (k) ) 2 ; γ (k) z = ∇J ( z (k) ) 2 / Y ξ D∇J ( z (k) ) 2 + Y z ∇J ( z (k) ) 2 , (10)
where ∇J (•) is the gradient of J (•). γ f and γ z are obtained by using optimized step length strategy, see [START_REF] Boyd | Convex optimization[END_REF]. The initialization of the hyper parameters αs and βs are:

α z0 = α 0 = α ξ0 = 2 + ε 1 , β z0 = β 0 = β ξ0 = ε 2
where ε 1 and ε 2 are small values close to zero.

III. SIMULATION RESULTS

In the simulations, the Shepp Logan phantom with size 64 3 is used as the original object. It consists of several homogeneous zones, each of which corresponds to a different material. Projections are applied in angles uniformly distributed from 0 • to 180 • . The reconstruction performance is measured in terms of the Normalize Mean Square Error (NMSE), or the relative error δ f , defined as

δ f = NMSE = f -f 2 / f
2 , where f is the reconstructed result.

The comparison is between a) the proposed method, b) the Quadratic Regularization (QR) method and c) the Total Variation (TV) method. The simulation of TV method is realized by using the Split Bregman method [START_REF] Goldstein | The Split Bregman method for L1regularized problems[END_REF]. The descend gradient optimization algorithm is used in HHBM method as mentioned in Eq. [START_REF] Nikolas | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF].

In order to deal with 3D data, considering the big data size constraint, a GPU processor is used. Several toolbox are accessible for usage of GPU in MATLAB. Here we used the ASTRA toolbox [START_REF] Jan Palenstijn | The ASTRA tomography toolbox[END_REF], which facilitates the coding with GPU in MATLAB.

The middle slice of the reconstructed phantom is shown in Fig. [START_REF] Stanley R Deans | The Radon transform and some of its applications[END_REF]. The left figures show the reconstructed middle slices of the 3D Shepp-Logan phantom, and on the right there are the details of the corresponding middle slice. 3) shows the δ f of reconstructed Shepp Logan phantom of size 64 3 by using 64 projections, when SNR=40dB and 20dB. It shows that when the noise is at a low level, the HHBM method converges much faster and gives a more precise reconstructed result, when with a more important noise, the TV method keeps robust, but the HHBM method still converges faster than the others. Fig.( 4) compares the δ f of the reconstructed phantom of size 64 3 with 32 projections. It shows that when there are insufficient projection data, the HHBM method remains robust, and the convergence is efficient. This property of the HHBM method is important in the applications where the reconstruction time is important and the results need to be obtained within a few iterations.

More results are given in Tab.(I), where err f and err g are respectively the NMSE of reconstructed phantom and reconstructed projections after 30 iterations. The value T c shows the computation time for each iteration. From the table and the figures we can see that, even though the computation time for one iteration of HHBM is longer than TV method (about 2 times), the number of iterations needed for reaching convergence by using HHBM is much smaller than TV method.

IV. CONCLUSION

We proposed to enforce sparsity of the multi-level Haar Transformation of the piece-wise continuous images by using a Generalized Student-t distribution, and use the Bayesian method to estimated the image and the transformation coefficients simultaneously. We compared the proposed HHBM method with the Quadratic Regularization and the Total Variation method. We conclude from the simulation results that the HHBM method is more robust for cases with less number of projections and higher noise level. The convergence speed is faster than the state of the art methods. As future work, we search for better sparse representation bases, and a more general Bayesian estimation going to Posterior Mean (PM) computation via Variational Bayesian Approximation (VBA).

Fig. 2 :

 2 Fig. 2: The middle slices and zones of the reconstructed 3D Shepp-Logan phantom by using different methods with 1/4 of total projections with noise of SNR=40dB. Top: QR, middle: TV, bottom: HHBM. The NMSE by using different methods are compared. Fig.(3) shows the δ f of reconstructed Shepp Logan phantom of size 64 3 by using 64 projections, when SNR=40dB and 20dB. It shows that when the noise is at a low level, the HHBM method converges much faster and gives a more precise reconstructed result, when with a more important noise, the TV method keeps robust, but the HHBM method still converges faster than the others. Fig.(4) compares the δ f of the reconstructed phantom of size 64 3 with 32 projections. It shows that when there are insufficient projection data, the HHBM method remains robust, and the convergence is efficient. This property of the HHBM method is important in

Fig. 3 :

 3 Fig. 3: Comparison of NMSE of reconstructed Shepp Logan phantom of size 64 3 by using different methods with 64 projections and noise of (a) 40dB and (b) 20dB.

Fig. 4 :

 4 Fig. 4: Comparison of NMSE of reconstructed Shepp Logan phantom of size 64 3 by using different methods with 32 projections and noise of (a) 40dB and (b) 20dB.

TABLE I :

 I Comparison of NMSE of reconstructed phantom with 30 iterations and the computation time of each iteration by using different methods.


Mircea Dumitru

Nicolas Gac

Li Wang

Ali Mohammad-Djafari

Unsupervised sparsity enforcing iterative algorithms for 3D image reconstruction in X-ray Computed Tomography

Keywords: sparsity, Normal variance mixtures, Bayesian inference, Computed Tomography

Unsupervised iterative reconstruction algorithms based on a Bayesian approach for piecewise constant images are presented in this paper. Such images can be expressed via a sparse representation and the reconstruction problem can be addressed using sparsity enforcing priors. We focus on sparsity enforcing priors expressed as Normal variance mixture, considering three mixing distributions: Inverse Gamma distribution, corresponding to Student-t prior, general inverse Gaussian distribution with the real parameter fixed, corresponding to Normal-inverse Gaussian prior and Gamma distribution corresponding to Variance-Gamma prior. We present and discuss the corresponding iterative algorithms considering the Joint Maximum A Posteriori estimation showing simulations results for 3D X-ray Computed Tomography.

I. INTRODUCTION

A widely used discretized linear forward model in the inverse problem approach for image reconstruction is given by g = Hf + ǫ,

where g represents the N × 1 observed data, H represents the N × M linear projection operator, f represents the M × 1 image to be reconstructed and ǫ accounts for measurement errors and model uncertainties. Piecewise continuous images f can be expressed as a transformation applied on a sparse structure z accounting for the uncertainties,

f = Dz + ξ. (2)
Equations ( 1) and [START_REF] Dumitru | Precise periodic components estimation for chronobiological signals through bayesian inference with sparsity enforcing prior[END_REF] represent the linear forward model considered during this paper. For Equation (2) different sparse representations can be considered. In the sparsity context, two main approaches have been successfully used for estimating the unknowns from Equations ( 1) and ( 2): deterministic methods [START_REF] Wipf | Iterative reweighted l 1 and l 2 methods for finding sparse solutions[END_REF], [START_REF] Hansen | The use of the l-curve in the regularization of discrete ill-posed problems[END_REF] and the Bayesian approach [START_REF] Dobigeon | Hierarchical bayesian sparse image reconstruction with application to mrfm[END_REF]. For the deterministic methods, an important issue is the choice of the regularization parameter.

We consider an inversion based on a Bayesian approach, building an hierarchical model, Subsection (II-A) accounting for the sparse structure of z, using sparsity enforcing priors. Sparsity enforcing priors expressed via conjugate distributions laws lead to analytical expressions for the unknowns of the hierarchical model therefore we focus on heavy tailed Normal (N ) variance mixtures: Student-t (St) prior, corresponding to Inverse Gamma (IG) mixing distribution, Normal-inverse Gaussian (N IG) distribution, corresponding to generalized inverse Gaussian (GIG) mixing distribution with real parameter fixed, p = -1/2 and Variance-Gamma (VG) distribution corresponding to Gamma mixing distribution (with particular case Laplace (L) prior corresponding to Exponential (E) mixing distribution), Subsection (II-B). The likelihood is obtained via the considered linear forward model and the distributions chosen to model the uncertainties ǫ and ξ. The hyperparameters (i.e. the parameters of the assigned variance distributions, θ ǫ , θ ξ and θ z ) are estimated along with the unknowns of the forward model f and z and the corresponding variances v ǫ , v ξ and v z from the posterior distribution obtained via the proportionality given by the Bayes rule

p(f , z, v ǫ , v ξ , v z |g) ∝p(g|f , v ǫ )p(f |z, v ξ )p(z|v z )p(v ǫ |θ ǫ ) p(v ξ |θ ξ )p(v z |θ z ). (3)
The corresponding JAMP iterative algorithms are derived and compared in Section (III). Details of the computations can be found in [START_REF] Dumitru | Precise periodic components estimation for chronobiological signals through bayesian inference with sparsity enforcing prior[END_REF] Simulation results for X-ray Computed Tomography (CT) reconstruction are presented in Section (IV) and conclusions are drawn in Section (V).

II. HIERARCHICAL MODEL AND NORMAL VARIANCE

MIXTURES

First, we briefly introduce the construction of the hierarchical model corresponding to the linear forward model Equations ( 1) and (2), in the context of Normal variance mixtures priors, Subsection (II-A).Then we introduce the sparsity enforcing priors discussed, both corresponding to Normal variance mixtures, Subsection (II-B).

A. Hierarchical model

The construction of the hierarchical model, Figure [START_REF] Dobigeon | Hierarchical bayesian sparse image reconstruction with application to mrfm[END_REF], corresponding to the linear forward model, Equations (1) and ( 2) is done using the assigned distributions. As mentioned, we consider two Normal variance mixtures: therefore, in both cases f , z are modelled as zero mean multivariate N distribution, deriving the likelihoods via the linear forward model. For the variances corresponding to the uncertainties we consider the IG and the GIG distributions. considers a zero-mean N distribution for z j |v zj and an IG distribution for the variance v zj |α z , β z , with the corresponding shape and scale parameters, α z and β z :

✒✑ ✓✏ g ✒✑ ✓✏ f ✒✑ ✓✏ v ǫ p(g|f , v ǫ ) = N (g|Hf , v ǫ ) ❄ ✠ H p(f |z, v ξ ) = N (f |Dz, v ξ ) p(v ǫ |θ ǫ ) ✒ ❅ | IG(v ǫ |θ ǫ ) GIG(v ǫ |θ ǫ ) ✒✑ ✓✏ v ξ ✒✑ ✓✏ z θ ǫ ❅ ❅ ❅ |❄ ❄ D p(z|v z ) = N (z|0, v z ) p(v ξ |θ ξ ) ❅ ■ ✠ IG(v ξ |θ ξ ) GIG(v ξ |θ ξ ) θ ξ ✒✑ ✓✏ v z ❄ ❄ p(v z |θ z ) ✒ ❅ | IG(v z |θ z ) GIG(v z |θ z ) θ z ❄ g = Hf + ǫ f = Dz + ξ
StPM: p(z j |0, v zj ) = N (z j |0, v zj ) p(v zj |α z , β z ) = IG(v zj |α z , β z ) (4)
The marginal of the joint probability distribution p(z j |α z , β z ) is a two parameters Student-t distribution with the probability density:

p(z j |α z , β z ) = Γ(α z + 1 2 ) √ 2πβ z Γ(α z ) 1 + z j 2 2β z -(αz + 1 2 ) . (5
NIGPM: p(z j |0, v zj ) = N (z j |0, v zj ) p(v zj |γ 2 z , δ 2 z ) = GIG(v zj |γ 2 z , δ 2 z , p z = -1 2 ) (6)
The marginal of the joint probability distribution p(z j , |γ z , δ z ) is a N IG distribution with zero location and asymmetry β parameters:

p(z j |γ z , δ z ) = γ z δ z K 1 γ z δ 2 z + z j 2 π δ 2 z + z j 2 exp {γ z δ z } , (7)
where K p denotes the modified Bessel function of the second kind.

Observation: A great interest in the use of the Normal-inverse Gaussian in the sparsity context is the possibility to control the tail heaviness via the parameter γ z . c) Variance Gamma: Considering the Gamma distribution as the mixing distribution for a Normal variance mixture, leads to the Variance-Gamma distribution as prior. The comparison between the standard Normal distribution N (x|0, 1) and the Variance-Gamma distribution VG (x|0, 1) is presented in Figure [START_REF] Wang | Computed tomography reconstruction based on a hierarchical model and variational Bayesian method[END_REF]. The Variance-Gamma Prior Model (VGPM), Equation (8), considers zero-mean Normal distributions for z j |v zj , and Gamma distributions with the corresponding shape and scale parameters α z and β z for the variances v zj |α z , β z :

VGPM: p(z j |0, v zj ) = N (z j |0, v zj ) p(v zj |α z , β z ) = G(v zj |α z , β z ) (8)
The marginal of the joint probability distribution p(z j |α z , β z ) is a VG distribution with the zero location and asymmetry parameters:

VG(x|α, β) = β 2α |x| α-1 2 K α-1 2 (α|x|) √ πΓ (α) (2β) α-1 2 . ( 9
)
Observation: Evidently, GIG(x|2β, δ ց 0, α) = G(x|α, β) so the VG prior can also be viewed as a Normal variance mixture with GIG mixing distribution. Also, we mention that the a particular case of the VG obtained via G mixing distribution is L distribution, obtained via E mixing distribution.

III. ITERATIVE ALGORITHMS In order the favour a sparse solution for z, Equation ( 2), we consider the following sparsity enforcing priors, expressed as Normal variance mixtures:

1) St via StPM, Equation ( 4), IG mixing distribution 2) N IG via NIGPM, Equation ( 6), GIG mixing distribution (with fixed real parameter, p z = 1/2) 3) VG via VGPM, Equation (8), G mixing distribution For JMAP estimation, the posterior distribution, Equation (3), is used. Because of the use of the sparsity enforcing priors expressed as Normal variance mixtures, for all prior models considered the first part of the posterior distribution is the same (i.e. p(g|f , v ǫ ), p(f |z, v ξ ), p(z|v z )). The differences between the resulting algorithms corresponds to the choice of the mixing distribution: IG, GIG or G. Consequently, the analytical expressions corresponding to f and z are the same. Numerically, the differences are induced during the iterations because of different expressions corresponding to variance estimates. Furthermore, different behaviours are associated with each prior model in terms of sensibility to the hyperparameters or rate of convergence. For JMAP, the unknowns are estimated by maximizing the posterior distribution, or minimizing the criterion L:

f , z, v ξ , v ǫ , v z = arg min (f , z, v ξ , vǫ, vz) L (f , z, v ξ , v ǫ , v z ) , (10)
where the criterion L is defined as:

L (f , z, v ξ , v ǫ , v z ) = -ln p (f , z, v ξ , v ǫ , v z |g) (11)
The optimisation algorithm considered is an alternate optimization of the criterion L (f , z, v ξ , v ǫ , v z ) with respect to the each unknown, leading to an interative algorithm,

f = H T V ǫ -1 H + V ξ -1 -1 H T V ǫ -1 g + V ξ -1 D z z = D T V ξ -1 D + V z -1 -1 D T V ξ -1 f Update variances v ξ j , v ǫ i , v z j St (12), N IG (13) or VG (14) prior form V ǫ , V ξ , V z
Fig. 5: Joint MAP iterative algorithm Figure [START_REF] Wipf | Iterative reweighted l 1 and l 2 methods for finding sparse solutions[END_REF]. Depending on the mixing distribution, the update of the variances is done as it follows: 1) For IG mixing distribution (St prior)

               v ξj = β ξ + 1 2 ( f j -Dj z) 2 α ξ + 3 2 v ǫi = βǫ+ 1 2 gi-H i f 2 αǫ+ 3 2 v zj = βz+ 1 2 zj 2 αz + 3 2 (12) 2) For GIG mixing distribution (N IG prior)                  v ξj = 2+ 4+γ 2 ξ δ 2 ξ ( f j -Dj z) 2 γ 2 ξ v ǫi = 2+ 4+γ 2 ǫ δ 2 ǫ gi-H i f 2 γ 2 ǫ v zj = 2+ √ 4+γ 2 z δ 2 z zj 2 γ 2 z (13) 3) For G mixing distribution (VG prior)                  v ξj = (αǫ-3 2 )+ (αǫ-3 2 ) 2 +2βǫ( f j -Dj z) 2 βǫ v ǫi = (αξ-3 2 )+ (αξ-3 2 ) 2 +2β ξ gi-H i f 2 β ξ v zj = (αz-3 2 )+ (αz-3 2 ) 2 +2βzzj 2 βz ( 14
)
f j denotes the element i of vector f , g i denotes the element j of vector g, H i denotes the line i of matrix H, D j denotes the line j of matrix D. For all three prior models, the variance matrices are given by:

v ξ = . . . v ξj . . . ; v ǫ = [. . . v ǫi . . .] ; v z = . . . v zj . . . V ξ = diag [ v ξ ] ; V ǫ = diag [ v ǫ ] ; V z = diag [ v z ] (15)
IV. SIMULATION RESULTS

For simulations the reconstruction of X-ray CT piecewise constant images problem is considered. The Shepp Logan phantom (256 3 ) is used as the original image. 128 projections are simulated uniformly between 0 • and 180 • . For the sparse representation of the image the multilevel Haar transform is used as in [START_REF] Wang | Computed tomography reconstruction based on a hierarchical model and variational Bayesian method[END_REF]. V. CONCLUSION Sparsity enforcing iterative algorithms via a Bayesian approach were considered in this paper. Heavy tailed distributions expressed as Normal variance mixtures were considered in order to obtain analytical expressions for the unknowns of the model. The prior models and their corresponding JMAP iterative algorithms were developed and compared: the reconstruction accuracy is similar for the three prior models considered, but the rate of convergence is different and the sensibility to the prior hyperparameters is different. This results are encouraging. However we are now investigating methods to compute the Posterior mean via Variational Bayesian Approximation (VBA). The structure of the algorithms is the same but the computations need the diagonal elements of the covariance matrices which are too costly for 3D applications.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Hierarchical model corresponding to the linear forward model, Eq. (1) and Eq. (2)

)Fig. 3 :

 3 Fig. 3: Comparison between N and N IG distribution. Model (NIGPM), Equation (6), considers zero-mean N distributions for z j |v zj , and generalized inverse Gaussian distributions for the variances v zj |γ 2 z , δ 2 z , with the corresponding parameters γ 2 z , δ 2 z and p z = -1/2:

Fig. 4 :

 4 Fig. 4: Comparison between the N and VG distribution.

Fig. 7 :

 7 Fig. 6: Slice comparison between the original volume (top left) and JMAP reconstructed volumes St (top right), N IG (bottom left), VG (bottom right). 64 3 phantom size, 128 projections, SNR=30dB

Fig. 8 :

 8 Fig. 8: Influence of the prior hyperparameters: NMSE vs. iterations for St (top) N IG (middle) VG (bottom) priors

M. Dumitru, N. Gac, L. Wang and A. Mohammad-Djafari are with Laboratoire des signaux et systèmes, CentraleSupélec, Gif-sur-Yvette, 91192, France, e-mail: surname.name@lss.supelec.fr


Camille Chapdelaine

Ali Mohammad-Djafari

Nicolas Gac

Estelle Parra-Denis

Estelle Parra

Joint Reconstruction and Segmentation of Real 3D Data in Computed Tomography thanks to a Gauss-Markov-Potts Prior Model

Keywords: 3D Computed Tomography, non-destructive testing, iterative reconstruction algorithm, joint reconstruction and segmentation, Gauss-Markov-Potts

published or not. The documents may come

Introduction

Computed tomography (CT) is widely used in industry for nondestructive testing (NDT). In many applications, X-rays are projected on a plan of detectors which measure the decrease of their intensity implied by their passing through the industrial part to be controlled. By rotating the object, a set of images is acquired during the acquisition process. Then, these images feed a reconstruction algorithm to retrieve the interior of the object, which make the user able to see if there is any default.

In industry, analytical reconstruction methods, based on Radon transform [START_REF] Smith | Mathematical foundations of computed tomography[END_REF][START_REF] Rodet | Algorithmes rapides de reconstruction en tomographie par compression des calculs: application à la tomofluoroscopie 3D[END_REF], are widely applied because of they are fast, the most popular of which is FDK algorithm [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF]. But these methods suffer from many drawbacks, like artefacts, which require post-reconstruction corrections [START_REF] Krumm | Reducing non-linear artifacts of multi-material objects in industrial 3d computed tomography[END_REF]. Moreover, they have a lack of flexibility when one wants to introduce prior knowledge in order to improve the quality of the reconstruction [START_REF] Mohammad-Djafari | Maximum entropy image reconstruction in X-ray and diffraction tomography[END_REF]. On the opposite, algebraic reconstruction techniques (ART) give the possibility to easily take into account this prior knowledge by adding a regularization term to the data matching term. By this way, while the first developed ART, like SIRT [START_REF] Gilbert | Iterative methods for the three-dimensional reconstruction of an object from projections[END_REF] or SART [START_REF] Andersen | Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm[END_REF], simply try to match the data, many recent advanced ART try to find a tradeoff between matching the data and minimization of, for instance, total variation [START_REF] Sidky | Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle & Pock algorithm[END_REF] or l 0 -norm [START_REF] Storath | Joint image reconstruction and segmentation using the Potts model[END_REF] of the volume to reconstruct. In order to work well and to find the best tradeoff, the most performing methods require to fit parameters, which can be very tedious because the influence of each of these parameters on the result is not always fully understandable. To overcome this difficulty, Bayesian framework is a good way to simultaneously introduce prior information in the reconstruction algorithm and define parameters of which the fitting is simple thanks to their very clear role.

In this framework, we have chosen to develop a Bayesian reconstruction method for NDT in industry using a Gauss-Markov-Potts (GMP) prior model, which has been successfully applied in microwave imaging [START_REF] Féron | Microwave imaging of inhomogeneous objects made of a finite number of dielectric and conductive materials from experimental data[END_REF] and image restoration [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational bayesian computation[END_REF]. Thanks to this model, our algorithm is able to reconstruct an industrial part containing several compact regions corresponding to different materials. The innovation of our method is that the segmentation jointly retrieved by our method with the reconstruction can be useful to see defects in the controlled object. In this paper, after having briefly presented the models, we present the results of our method on real 3D phantom, and we compare it to a standard FDK reconstruction and Total Variation (TV) regularization. The goal of this paper is to show the actual efficiency of using Gauss-Markov-Potts prior model on real 3D data : the reconstruction obtained by our method has more constrasted regions than FDK and TV. In [START_REF] Chapdelaine | A 3D Bayesian Computed Tomography Reconstruction Algorithm with Gauss-Markov-Potts Prior Model and its Application to Real Data[END_REF], the details of the algorithm and some discussion about the choice of the estimators can be found.

Prior models 2.1 Forward model

We denote by g the projections, and by f the object to be reconstructed. The linearized discretized forward model for the projections accounting for the noise is :

g = Hf + (1) with i ∼ N (0, v i ), ∀i ∈ {1, . . . , M } . (2)
The variances v i are unknowns and are modeled as following an Inverse-Gamma distribution

v i ∼ IG(α 0 , β 0 ), ∀i ∈ {1, . . . , M } (3)
where α 0 and β 0 are fixed parameters.

Gauss-Markov-Potts prior model

We now present the model for the controlled object f itself. This object is composed of several materials mapped in several compact and more or less homogeneous regions. We label each voxel j depending on the material to which the voxel belongs. This leads to introduce a hidden discrete field define by : z j = k if voxel j composes material k, k ∈ {1, . . . , K}, with K the total number of materials. The material, or label or class, z j of a voxel depends on the ones of its neighbours i ∈ V(j) : that is why we assign markovian model for z, which is a Potts model :

p(z|γ0, α) ∝ exp   N j=1   K k=1 α k δ(zj -k) + γ0 i∈V(j) δ(zj -zi)     (4)
where α and γ 0 are fixed parameters. Potts coefficient γ 0 fits the compacity of the regions. There exits a critical value γ c above which regions are compact, as illustrated in figure 1. This critical value has been exactly computed in 2D for Ising field (K = 2) [START_REF] Giovannelli | Estimation of the Ising field parameter thanks to the exact partition function[END_REF]. Given its label z j = k, voxel j has a gray level f j distributed around a mean m k with a variance v k

f j ∼ N (m k , v k ) if z j = k. (5)
The means m and variances v of the classes are unknown and have the following priors

m k ∼ N (m 0 , v 0 ), ∀k ∈ {1, . . . , K} (6)
and

v k ∼ IG(α 0 , β 0 ), ∀k ∈ {1, . . . , K} (7)
where m 0 , v 0 , α 0 and β 0 are fixed parameters.

Joint Maximization A Posteriori (JMAP) algorithm

After having built our prior models, the reconstruction algorithm consists in estimating each unknown f , z, v , m and v by maximizing their joint posterior probability

( f , ẑ, v , m, v) = arg max (f ,z,v ,m,v) {p(f , z, v , m, v|g; M)} (8)
thanks to the approximate maximization iterative algorithm 1) , m (t-1) , v (t-1) |g; M)

                     f (t) ← arg max f p(f , z (t-1) , v (t-
z (t) ← arg max z p(f (t) , z, v (t-1) , m (t-1) , v (t-1) |g; M) v (t) ← arg max v p(f (t) , z (t) , v , m (t-1) , v (t-1) |g; M) m (t) ← arg max m p(f (t) , z (t) , v (t) , m, v (t-1) |g; M) v (t) ← arg max v p(f (t) , z (t) , v (t) , m (t) , v|g; M) (9)
where index t denotes the iteration t. This algorithm is a type of variational bayesian approximation [START_REF] Mohammad-Djafari | Variational Bayes and mean field approximations for Markov field unsupervised estimation[END_REF] which is a local optimization. By applying Bayes rule

p(f , z, v , m, v|g; M) ∝ p(g|f , v ) p(f |z, m, v) p(v |α 0 , β 0 ) p(z|α; γ 0 ) p(m|m 0 , v 0 ) p(v|α 0 , β 0 ), (10
) each step of the algorithm can be performed, as we have detailed in [START_REF] Chapdelaine | A 3D Bayesian Computed Tomography Reconstruction Algorithm with Gauss-Markov-Potts Prior Model and its Application to Real Data[END_REF]. Thanks to this method, we see that an estimation of the object f and the labels z, i.e. a segmentation of the object, are retrieved. This is a great advantage, because there is only one source of uncertainties in the control process, while, when applying many other methods, the reconstruction algorithm and the post-reconstruction segmentation algorithm are two different sources of uncertainties.

Results on real 3D IQI data

The considered volume is a real 3D IQI phantom [START_REF] Gay | Fantôme destiné à être utilisé pour le contrôle de la qualité d'images tomographiques[END_REF], discretized in 512 × 512 × 256 voxels. The data are 300 projections of size 512 × 512 pixels. The number of materials in the volume to reconstruct is K = 4. We apply our method and compare it with FDK [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF] and Total Variation (TV) minimization reconstruction methods. Concerning TV, it is implemented by following the algorithm proposed in [START_REF] Goldstein | The split Bregman method for L1-regularized problems[END_REF], based on Bregman iteration [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF].

Initialization

To apply JMAP algorithm described in section 3, an initial volume and an initial segmentation are required. As initial volume, because we have many projections, we use FDK reconstruction, shown in figure 2. Next, we need a segmentation of this initial volume as initial segmentation. This segmentation is obtained by a non-uniform thresholding of FDK reconstruction.

To determine the thresholds, we estimate the histogram of the initial volume. Then, we apply the peak-picking algorithm explained in [START_REF] Kurugollu | Color image segmentation using histogram multithresholding and fusion[END_REF] and originally developped for non-parametric clustering [START_REF] Koontz | A graph-theoretic approach to nonparametric cluster analysis[END_REF]. For each bin, we seek the nearest peak in its neighbourhood : this peak is called the parent-bin of the current bin, which is so called a child-bin of this peak. Each parentbin and its children are seen as a class, i.e. the voxels in the children-bins of parent-bin k are assigned the initial class k. Because we know the number of classes (K = 4), each parentbin is assigned as the child-bin of a higher parent-bin in its neighbourhood until there are only K = 4 parents-bins. This gives the segmentation shown in figure 3, used as initial segmentation by JMAP algorithm.

Results and comparison with FDK

Table 1 summarizes the values of most of the parameters we fixed to perform the joint reconstruction and segmentation of IQI data. The controlled volume is composed of 3 materials and air : that is why we fix K = 4. The other parameters are fixed according to a strategy we propose in [START_REF] Chapdelaine | A 3D Bayesian Computed Tomography Reconstruction Algorithm with Gauss-Markov-Potts Prior Model and its Application to Real Data[END_REF]. Thanks to the clear Bayesian prior model, tuning the parameters of JMAP is straightforward, while it is more tedious for TV. Next, we run JMAP algorithm and we obtain the reconstruction and the segmentation shown in figures 6 and 7 respectively. When comparing FDK, TV and JMAP reconstructions in figures 2, 4 and 6 respectively, we can see the external circle surrounding the object in FDK reconstruction is completely erased by JMAP algorithm. Moreover, the JMAP reconstruction is more constrasted than the FDK and TV reconstructions, especially for the white region in the top of the object. This better contrast is still emphasized when in figure 8, we plot the profiles of the little holes in the bottom of the object : the profile of JMAP reconstruction is sharper than the one of FDK reconstruction.

Parameters α 0 K γ 0 v 0 α 0 β 0 Fixed values 2.1 4 3 1 5 0.01
To deepen the comparison, we have defined three reconstruction quality indicators. The first focuses on the compactness of the classes and is defined as :

C omp = 1 K K k=1 1 N k j∈R k 1 N V i∈V(j) δ(k -z i ) (11)
where R k = {j|z j = k}. The second measures the distinguishability of the classes on the contours :

d ist = 1 - 1 K K k=1 1 N k j∈R k ( dist ) j (12)
with, if z j = k :

( dist ) j =      i∈V(j) exp(-(fj -fi) 2 )(1-δ(k-zi)) i∈V(j) (1-δ(k-zi))
if j is on a contour 0 otherwise (13) The third indicator which quantifies how homogeneous the classes are :

h omo = 1 K K k=1 1 N k j∈R k ( ĥomo ) j (14)
with, if z j = k :

( ĥomo ) j = i∈V(j) exp -(f j -f i ) 2 δ(k -z i ) i∈V(j) δ(k -z i ) (15)
if i∈V(j) δ(k -z i ) = 0. Otherwise, ( ĥomo ) j = 0.

In this way, we have three reconstruction quality indicators which respectively measure compactness, distinguishability and homogeneity of the classes in the reconstruction. These indicators are designed such that the higher they are, the better the reconstruction is. One shall notice they need a segmentation of the reconstruction to be applied. In table 2, we compute them on the joint recontruction and segmentation retrieved by JMAP, on the initial reconstruction by FDK in figure 2 and the initial segmentation of this reconstruction in figure 3, and on TV reconstruction and its posterior segmentation in figure 5, obtained thanks to a Gauss-Markov-Potts prior model proposed in [START_REF] Rosu | Potts model parameter estimation in Bayesian segmentation of piecewise constant images[END_REF]. Because the labels of each of the three reconstruction are very similar, C omp is almost the same for the three methods. Concerning d ist and h omo achieves the highest values, which confirms that our method outperform FDK and TV.

Indicator

Conclusion and perspectives

In this paper, we have presented results on real 3D phantom of a reconstruction method which has achieved high compactness, distinguishability and homogeneity of the classes obtained by the segmentation performed jointly with the reconstruction. Thanks to clear prior models, the parameters of the algorithm have been proved easy to fit. The method has shown better results than both FDK and TV reconstructions. An important future work would be to propose a massively parallelized version of the algorithm on GPU.

Figure 1 :

 1 Figure 1: Potts field for different values of γ 0

Figure 2 :Figure 3 :Figure 4 :Figure 6 :Figure 7 :Figure 8 :

 234678 Figure 2: FDK reconstruction (bottom and top)

Table 1 :

 1 Reconstruction parameters for 3D IQI data

Table 2 :

 2 Comparison of FDK, TV and JMAP using reconstruction quality indicators

		C omp	d ist	h omo
	FDK	94.4 % 78.4 % 78.3 %
	TV	94.5 % 77.5 % 77.5 %
	JMAP	94.3 % 79.0 % 78.8 %


Maxime Martelli

Mickael Seznec

Nicolas Heemeryck

Nicolas Gac

An OpenCL design for tomographic reconstruction on FPGA. Innovate Europe Design Contest Winners

3D Computed Tomography Projection

where pupϕ 0 , cq, vpϕ 0 , cqq are the values on the sinogram of the beam passing through c for ϕ ϕ 0 , and wpϕ 0 , cq is the distance weight. [START_REF] Lu | A 3D distance-weighted Wiener filter for Poisson noise reduction in sinogram space for SPECT imaging[END_REF] OpenCL Work-group mechanism There are two OpenCL kernel (device functions) programming models : NDRange (ND), similar to GPU programming and single work-item (SWI), similar to CPU programming. OpenCL allows work-group and work-item partitioning, as shown in Fig. 2, and the challenge is to effectively use the inherent mechanism of each memory level to accelerate our design. This repartition allows effective internal communication and data synchronization between work-items of the same work-group.

Memory Benchmark

For benchmarking purposes, we implemented a custom routine program to measure the mean latency of each memory type on an Altera Cyclone V FPGA.

FPGA obtained speedup and GPU comparison

• The Altera Offline Compiler is effective in guaranteeing no kernel stall.

• Optimizing a SWI kernel lies with optimizing memory handling and data streaming effectiveness in order to increase kernel frequency. Even with the extrapolated results on an Aria 10 FPGA, and with the obtained improvements over the Naive versions, the FPGA has a much longer execution time than the GPUs, mostly due to the back-projection algorithm being adapted to SIMD architectures.

Conclusion

1. With an overall 8.74 speedup between naive and optimized kernels, there is room for more improvements closely related to memory handling. 2. On this application, FPGAs still fall short compared to GPUs, mostly due to the algorithm inadequacy for pipeline programming.

Figure 1 :s

 1 Figure 1: Figure caption

Figure 2 :

 2 Figure 2: Work-group queueing mechanism and OpenCL Memory Architecture

•

  SWI Naive : CPU like naive version. • SWI+SRP+LM : Shift Register Pattern (FIFO queue) with local memory optimization to increase bandwith. • ND+Naive : GPU like naive version (shared local memory). • ND+2CU : Kernel replication of the ND+Naive version. • ND+LM+MP : Last implemented optimization. Focused on a prefetching mechanism bound to the algorithm specificity. The goal is to access the pattern needed by a group of voxel in a single call.

Figure 3 :

 3 Figure 3: Sinogram memory fetching pattern optimization.

  An OpenCL design for tomographic reconstruction on FPGA. Innovate Europe Design Contest Winners Maxime MARTELLI, Mickael SEZNEC, Nicolas HEEMERYCK TSA & L2S & SATIE, Telecom ParisTech, ECE Industrial Context 1. With the end of Moore's Law, the semi-conductor industry seeks a reliable way to pursue the performance improvements of the last decades, and architecture-algorithm adequation is a solution for this new landscape. 2. Manufacturers like Intel and Xilinx are pushing for an FPGA resurgence, offering software suites and FPGAs card focused on a software-like FPGA programming model[1].

Table 1 :

 1 Memory

	Memory structure Mean latency (cycles)
	Global	240
	Constant	10
	Local	15
	Private	2

structure latency on an Altera Cyclone V.

•

  Optimizing a NDRange kernel coincide with reducing the logical footprint to allow more kernel replication within the same chip.

	Device	Power (W) Execution time (ms) Energy (mWh)
	Titan X Pascal	250	12	0.83
	Jetson TK2	15	94	0.39
	Intel Arria 10	2.27	991	0.63

Acknowledgments

We would like to thank Terasic, CNFM, and Intel for providing hardware and software used in this project, and Nicolas Gac for his supervision.


Camille Chapdelaine

Ali Mohammad-Djafari

Nicolas Gac

Estelle Parra

A JOINT SEGMENTATION AND RECONSTRUCTION ALGORITHM FOR 3D BAYESIAN COMPUTED TOMOGRAPHY USING GAUSS-MARKOV-POTTS PRIOR MODEL

Keywords: Gauss-Markov-Potts, 3D Computed Tomography, joint reconstruction and segmentation

Gauss-Markov-Potts models for images and its use in many image restoration, super-resolution and Computed Tomography (CT) have shown their effective use for Non Destructive Testing (NDT) applications. In this paper, we propose a 3D Gauss-Markov-Potts model for 3D CT for NDT applications. Thanks to this model, we are able to perform a joint reconstruction and segmentation of the object to control, which is very useful in industrial applications. First, we describe our prior models for each unknown of the problem. Then, we present results on simulated data and compare them to those of Total Variation (TV) minimization algorithm. Two quality indicators exploiting the segmentation are also proposed.

INTRODUCTION

Computed tomography (CT) is a powerful imaging tool to see the interior of a three dimensional object and has a wide field of applications. In particular, in industry, the most used one is the so-called FDK algorithm [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF] for cone-beam CT, which performs a filtered back-projection (FBP) algorithm. This algorithm is part of analytical reconstruction methods, based on the use of Radon transform [START_REF] Smith | Mathematical foundations of computed tomography[END_REF]. There exist other analytical reconstruction methods [START_REF] Jackson | Selection of a convolution function for Fourier inversion using gridding [computerized tomography application[END_REF] which use Fourier slice theorem but need to apply a filter in order to deal with a tricky interpolation in Fourier domain.

These analytical reconstruction methods suffer from artifacts due to approximations and give poor results with limitedangle projections. From this standpoint, algebraic reconstruction techniques (ART) and iterative methods have been provided for the last decades and consider a linearized discretization of Radon transform [START_REF] Smith | Mathematical foundations of computed tomography[END_REF]. In order to account for the errors and to obtain a better reconstruction, most of these methods apply a regularization which enforces some priors on the object to reconstruct [START_REF] Mohammad-Djafari | Maximum entropy image reconstruction in X-ray and diffraction tomography[END_REF].

In our application, we aim at reconstructing an industrial part which is composed of several materials filling one or several compact and quite homogeneous regions. Gauss-Markov-Potts model has been successfully applied in microwave imaging [START_REF] Féron | Microwave imaging of inhomogeneous objects made of a finite number of dielectric and conductive materials from experimental data[END_REF] and image restoration [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational bayesian computation[END_REF]. In this paper, we present a 3D Gauss-Markov-Potts model for CT for non-destructive testing (NDT) in industry, based on the works done in [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational bayesian computation[END_REF]. A joint maximization a posteriori (JMAP) is also proposed, which jointly retrieves an estimation of the object to reconstruct and a segmentation of the reconstructed object in compact and well-distinguishable regions.

MODELS

The forward model is based on the discretization of Radon transform. Denoting by g the data, by f the object to reconstruct and by H the projection operator, we write the forward projection model :

g = Hf + (1)
Errors i are modeled as zero-mean Gaussian with unknown variances (v i ) i , ∀i, which are modeled as following an Inverse Gamma distribution :

p(v i |α 0 , β 0 ) = IG(v i |α 0 , β 0 ), ∀i (2)
where α 0 et β 0 are fixed hyperparameters.

For the object, as said in introduction, we use a Gauss-Markov-Potts model. We define a hidden field z that labels each voxel by the material it represents : for instance, z j = k if voxel j is part of material k. We assume that the number K of materials in the object is known. Then, homogeneity of each material is enforced by modeling that the gray values of the voxels of a same material k are distributed around a mean value m k , with variance

v k p(f j |z j = k, m k , v k ) = N (f j |m k , v k ) if z j = k. (3)
The priors for the means and the variances of the classes are

p(m k |m 0 , v 0 ) = N (m k |m 0 , v 0 ) (4)
and

p(v k |α 0 , β 0 ) = IG(v k |α 0 , β 0 ) (5)
where m 0 , v 0 , α 0 and β 0 are fixed hyperparameters. To translate the compacity of each material, we assign a Markov-Potts prior model for the hidden field z:

p(z|α, γ 0 ) ∝ exp   j K k=1 α k δ(z j -k)+ γ 0 i∈V(j) δ(z j -z i )     (6)
Parameter γ 0 allows to control the compactness of the regions, as shown in figure 1. Here, γ 0 is tuned sufficiently large so the regions are compact. These prior models are very similar to the ones proposed in [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational bayesian computation[END_REF] for image restoration. The main difficulty here is that we deal with very huge 3D objects : this makes matrix H unstorable in memory.

JOINT MAXIMIZATION AND RESULTS

Our hierachical model is summarized in figure 2. We estimate the unknowns f , z, v , m and v by maximizing the joint posterior distribution of the unknowns (JMAP) :

p(f , z, v , m, v|g; M) ∝ p(g|f , v ) p(f |z, m, v) p(v |α 0 , β 0 ) p(z|α; γ 0 ) p(m|m 0 , v 0 ) p(v|α 0 , β 0 ) (7)
Equation ( 7) highlights the importance of our priors. Thanks to this expression, at iteration t of the algorithm, in order to estimate and update each unknown, we can successively maximize its posterior distribution given the other unknowns : this is an approximate maximization of the posterior distribution which makes possible to achieve joint reconstruction and segmentation of the controlled object.

The method is tested on 256 3 size Shepp-Logan 3D phantom, of which the middle slice is shown in figure (3a). The corresponding segmentation is shown in figure (3b): the number of materials is K = 5. 64 projections with 256 2 pixels are obtained from this phantom and are noisy with signal-tonoise ratio (SNR) equal to 20 db. Then, the algorithm is applied with γ 0 = 3 and the reconstruction and the segmentation, of which middle slices are shown in figures (3c) and (3d), are obtained. The core of the algorithm speedup is the parallelization on GPU of the projector H and backprojector H T , which makes gradient descent in the object estimation step fast. For the reconstruction and segmentation in figures (3c) and (3d), we have achieved a total computation time about 10 minutes. The method is compared to Total Variation (TV) minimization algorithm, to which a posterior Potts segmentation is applied. The results for TV are shown in figures (3e) and (3f). TV reconstruction retrieves the bone better than our algorithm, but joint segmentation in figure (3d) retrieves the detail on the front while this is not the case in figure (3f), although it is retrieved in TV reconstruction. This emphasizes the great interest of joint reconstruction and segmentation, that is not cumulating errors from reconstruction and segmentation algorithms.

Moreover, we have used the segmentation to define two quality indicators designed so that the higher they are, the better the reconstruction is. The compactness indicator gives the average rate of voxels which are only surrounded by neighbours in the same class as them : it is 89.2% for our method reconstruction and 88.9% for TV. The distinguishability indicator measures how two voxels in different classes on the contours of the regions are distinguishable : it is 73.9% for our method reconstruction, and 72.8% for TV.

CONCLUSION AND PERSPECTIVES

In this paper, we have presented a comprehensive prior model to perform joint reconstruction and segmentation in 3D CT. This method has been shown to retrieve results comparable to TV. The two herein proposed quality indicators have made us able to see that joint segmentation is better than posterior segmentation. Future works will focus on full implementation on GPU and some variations in the prior model.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Potts field z for different values of γ 0 , with K = 5 classes

Fig. 3 :

 3 Fig. 3: Shepp-Logan 3D phantom (a) and its segmentation (b), obtained joint reconstruction (c) and segmentation (d) by our method, and obtained reconstruction by TV (e) and its posterior Potts segmentation (f) (middle slices)


Li Wang

Nicolas Gac

Ali Mohammad-Djafari

Bayesian 3D X-ray Computed Tomography image reconstruction with a Scaled Gaussian Mixture prior model

Keywords: Computed Tomography, Limited projections, Non Destructive Testing (NDT), Hierarchical Model, Bayesian JMAP, Variational Bayesian Approximation (VBA), Gaussian, Mixture of Gaussians

.

. Sometimes for a discrete image, we can do the segmentation and reconstruction at the same time, then the reconstruction can be done with less projections.

INTRODUCTION OF COMPUTED TOMOGRAPHY Computed Tomography

X-ray computed tomography (X-ray CT) is a technology that uses computerprocessed X-ray to produce tomographic images of specific areas of a scanned object, allowing the users to see what is inside it without cutting it. Digital geometry processing is used to generate a three-dimensional image of the inside of an object from a large series of two-dimensional radiographic images taken around a single axis of rotation. X-ray tomographic image reconstruction consist of determining an object function from its projections. The main forward problem of the tomography is the Radon transform.

When considering a practical problem, a discretization is necessary. If we discretize f (x, y) into pixels and put all the pixels in a vector f f f and put all the data g(φ φ φ , r r r) of different angles φ φ φ in a vector g g g, we then obtain:

g g g = H H H f f f + ε ε ε (1)
where ε ε ε represents the error, and H H H is the projection operator in which the element H H H i j represents the length of the ray i in the pixel j.

BAYESIAN APPROACH

From this point, the main objective is to infer on f f f given the data g g g assuming the forward model g g g = H H H f f f + ε ε ε. By being Bayesian, we mean to use the Bayes rule:

p ( f f f | g g g) = p(g g g | f f f )p( f f f ) p(g g g) ∝ p(g g g | f f f )p( f f f ) (2)
to obtain what is called the posterior law p( f f f |g g g) from the likelihood p(g g g| f f f ) and the prior p( f f f ). To be able to use the Bayesian approach, first we need to assign p(g g g| f f f ) and p( f f f ).

Then, we can obtain the expression of the posterior law. Finally, we can infer on f f f using this posterior law.

Markov model. In the Markov model, the value of f j has a relationship with the values of its neighbours. For example, in the case of 1D, f j = F( f j1 , f j2 , • • • ), and in an image, the value of a pixel depends on the values of its neighbours pixels.

When the value of f j depends only on the values of the neighbours of distance 1, the model is called to have order 1, and for a Gaussian model we have:

p( f j | f j-1 , σ 2 ) = N f j | f j-1 , σ 2 ∝ exp - 1 2
( f j -f j-1 ) 2 σ 2 f ( 3
)
From this we can write:

p( f f f ) ∝ exp -γ [D D D f f f ] j β
where D D D f f f represents a suitable defined gradient of the image, γ is a scale parameter and β a shape parameter. For β = 2 a Gauss-Markov model is obtained.

For the noise term ε ε ε we choose a Gaussia prior law. This leads to

p(g g g | f f f ) ∝ exp - 1 2σ 2 ε g g g -H H H f f f 2 and so p( f f f | g g g, σ 2 f , σ 2 ε ) ∝ exp - 1 2σ 2 ε J J J( f f f ) with J J J( f f f ) = g g g -H H H f f f 2 + λ D D D f f f β β
, where λ is the parameter of the regularization term. When β = 2 (Gauss-Markov case), we got the analytical solution:

f f f = H H H t t t H H H + λ D D D t t t D D D -1 H H H t t t g g g
We have considered different prior laws with Markov model: the Gaussian law, the Generalized Gaussian law, the Cauchy law and the Huber law.

UNSUPERVISED BAYESIAN

In previous section, there are parameters σ 2 ε and σ 2 f as well as β which have to be assigned. These are called hyper parameters. For a practical applications we need to estimate them also. In the Bayesian approach this can be done via the joint posterior:

p( f f f , θ θ θ | g g g, θ 0 ) = p(g g g | f f f , θ 1 )p( f f f | θ 2 )p(θ θ θ | θ 0 ) p(g g g | θ 0 ) (4)
where θ θ θ = [θ 1 , θ 2 ] and θ 0 is a parameter of θ 1 and θ 2 .

In this paper we set β = 2 but we want to estimate

θ 1 = v ε = σ 2 ε and θ 2 = v f = σ 2 f .
As θ 1 and θ 2 are variances, we use a conjugate prior for them. Here we choose the inverse Gamma law as the conjugate prior:

p(v ε | α ε , β ε ) = I G (v ε | α ε , β ε ) p(v f | α f , β f ) = I G (v f | α f , β f ) (5)
One way to estimate the unknowns of our model is to compute the Joint Maximum A Posteriori (JMAP) [START_REF] Mohammad-Djafari | Gauss-Markov-Potts Priors for Images in Computer Tomography Resulting to Joint Optimal Reconstruction and segmentation[END_REF].

f f f , θ θ θ = arg max ( f f f ,θ θ θ ) p ( f f f , θ θ θ | g g g, θ 0 ) (6)
In the case where the only hyper parameters are v ε and v f , we can apply the following iterative algorithm:

           f f f = arg max f p f f f , v ε , v f | g g g, θ 0 = H H H t t t H H H + v ε v f D D D t t t D D D -1 H H H t t t g g g vε = 1 2 g g g-H H H f f f 2 +β ε α ε + M 2 +1 v f = β f + 1 2 D D D f f f 2 α f + N 2 +1 (7)

JMAP AND VBA WITH STUDENT-T PRIOR

In this section, we will consider the hierarchical model which uses a Student-t distribution for modeling the distribution of sparse signals or images. For an image in which most parts are homogeneous, the gradient of the image is sparse. To enforce sparsity, we propose to use a heavy tailed prior law, for example, the Generalized Gaussian law and the Student-t law. Here we propose to use the Student-t distribution. For the Student-t prior law, we have the property:

S S St t t f j | ν, τ = ∞ 0 N f j | 0, 1 z j G z j | ν 2 , ντ 2 2 dz j
where the z j is a hidden variable which represents the inverse variance of f j [START_REF] Mohammad-Djafari | Hierarchical Markovian models for 3D computed tomography in non destructive testing applications[END_REF][START_REF] Mohammad-Djafari | Bayesian inference with hierarchical prior models for inverse problems in imaging systems[END_REF]. This property of Infinite Gaussian Mixture gives the possibility to propose the following hierarchical model:

p( f j | z j ) = N f j | f j-1 , v f z j p(z j | α z j , β z j ) = G z j | α z j , β z j JMAP.
With the hierarchical model, we can obtain the expression of the joint a posterior:

p ( f f f , Z Z Z, θ θ θ | g g g) ∝ p (g g g | f f f , v ε ) p f f f | Z Z Z, v f p (Z Z Z | α Z Z Z , β Z Z Z ) p (v ε ) p v f where p f f f | Z Z Z, v f = ∏ j p f j | z j , v f and p (Z Z Z) = ∏ j
p z j and for p (v ε ) and p v f we use Inverse Gamma law which is given in [START_REF] Mohammad-Djafari | Gauss-Markov-Potts Priors for Images in Computer Tomography Resulting to Joint Optimal Reconstruction and segmentation[END_REF].

An alternating optimisation of this JMAP criterion results to the following algorithm:

                               f f f (k+1) = H H H t t t H H H + v (k) ε v (k) f Z Z Z (k) -1
H H H t t t g g g ẑ(k+1)

j = α z j -1 2 β z j + f (k+1) j 2 2v (k) f v(k+1) ε = 1 2 g g g-H H H f f f (k+1) 2 +β ε α ε + M 2 +1 v(k+1) f = 1 2 ∑ N j=1 z (k+1) j f (k+1) j 2 +β f α f + N 2 +1 (8)
VBA. As we will see, the main inconvenience of the JMAP approach is that we are summarizing the joint posterior law p( f f f , θ θ θ |g g g) by its mode only. Also, for obtaining this mode, in general, an iterative alternating optimization is used, where at each iteration only the values of the estimates at previous iterations are used without accounting for their corresponding uncertainties.

In the VBA approach, the main idea here is to approximate p( f f f , θ θ θ |g g g) by a separable law q( f f f , θ θ θ ) = q 1 ( f f f )q 2 (θ θ θ ) = ∏ j q 1 j ( f j )q 2 (θ θ θ ) which can then be used for inferring on f f f or θ θ θ . The main criterion used is the Kullback-Leibler divergence [START_REF] Bishop | Variational relevance vector machines[END_REF][START_REF] Miyamoto | Phase diagrams of a variational Bayesian approach with ARD prior in NIRS-DOT[END_REF][START_REF] Molina | Variational posterior distribution approximation in bayesian emission tomography reconstruction using a gamma mixture prior[END_REF]. As we will see, f f f depends on q 2 (θ θ θ ) and q 2 (θ θ θ ) depends on q 1 ( f f f ), thus accounting for uncertainties in both steps.

We assume that

f j | µ j , v f j , z j ∼ N ( f j | µ j , σ f j z j ), z j | α z j , β z j ∼ G (z j | α z j , β z j ), v ε | α ε , β ε ∼ I G (v ε | α ε , β ε ) and v f | α f , β f ∼ I G (v f | α f , β f ).
And we are now considering all the hyper parameters, parameters and unknowns. The alternate optimization of VBA is then given by:

                                                                 λ (k+1) j = α(k) f β (k) f β (k) ε α(k) ε α(k) z j β (k) z j 1 A μ j (k+1) = - B A 1+ λ (k+1) j σ f j (k+1) = 1 αε (k) βε (k) A 1+ λ (k+1) j α(k+1) z j = α (k) z j + 1 2 β (k+1) z j = β (k) z j + 1 2 µ (k+1) j 2 α(k) f β (k) f z(k+1) j = α(k+1) z j β (k+1) z j α(k+1) ε = α (k) ε + M 2 β (k+1) ε = β (k) ε + 1 2 g g g -H H H μ µ µ (k+1) 2 + ∑ N j=1 H H H t t t H H H j j σ (k+1) f j α(k+1) f = α (k) f + N 2 β (k+1) f = β (k) f + 1 2 ∑ N j=1 α(k+1) z j β (k+1) z j ( μ(k+1) j ) 2 + σ (k+1) f j (9)
where N is the length of the vector f f f , M is the length of the vector g g g, and the definition of A and B are:

A = H H H t t t H H H j j B = -H H H t t

IMPLEMENTATION

In the forward model g g g = H H H f f f + ε ε ε , the vector f f f correspond to the object, H H H is projection operator and g g g represents the projections. Here in our simulation, we use the synthetic volume "Shepp and Logan" with size 256 × 256 × 256 voxels. In this report, the projection has been done in 256 angles, and under each angle the screen of detectors will receive an image of size 256 × 256 pixels. The matrix H H H corresponds to the matrix of projection, of which the most component are zero. In the problem 3D, we don't know the matrix H H H, but the projection H H H f f f and back-projection H H H g g g can be accessed. The term ε ε ε represents the noise which follows a Gaussian law. In the algorithm, the most costly parts of the computation are operations on H H H f f f (projection) and the H H H t t t g g g (back-projection). These two operations are implemented using a GPU.

NUMERICAL RESULTS

In the numerical implementation part, we have compared the JMAP method with the results of differents reconstruction methods: simple back-projection, filtered backprojection, lest square method and Bayes method (MAP) with different regularizations (Simple Gaussian, Generalized Gaussian, Cauchy and Huber). The middle slice after reconstruction of different methods (after 200 iterations) is shown in F F FI I IG G GU U UR R RE E E 1 1 1, and we can compare the relatif error of different method, where

δ f = f f f -f f f 2 f f f 2 and δ g = ĝ g g-g g g 2 ĝ g g 2 .
With the theoretical projection H H H and back-projection H H H t t t , the object obtained has a large error, and we can not distinguish the details in the image. With the Feldkamp back-projection, we can obtain an image which is already clear enough to distinguish the different materials, but the boundaries are blurred. When we apply the optimisation with different prior models, we can see that the boundaries are easier to distinguish and different materials are clearly separated. In the JMAP method, the boundaries of different materials are more clear than the others methods. A more accurate compare is showed in F F FI I IG G GU U UR R RE E E 2 2 2, the error of the JMAP method is smaller than the other non-hirarchical methods.

We can see that the method JMAP has the lowest error among the methods that we considered. With the JMAP method, the borders of two different materials are more distinct.

The VBA method for the big data 3D problem is more complicated than the JMAP method. The calculate of the diagonals elements of the matrix H H H t t t H H H for the case 3D is very difficult because of the unknown huge size matrix H H H and the time needed for the projection and back-projection. Thus we have compared the VBA method and other methods for a case 1D where the f f f is of size N = 29 , g g g is of size M = 96 and the matrix H H H is known. The relative error of different methods in this case is showed in f f f i i ig g gu u ur r re e e 3 3 3.

From the compare of different method we can see that the VBA method works better than the others.

CONCLUSIONS

In this paper, JMAP method and VBA method are proposed for doing Bayesian computations for inverse problems where a hierarchical prior modeling is used for the unknowns. A Student-t prior model, which can be written via hidden variables, is considered, and it gives the model a hierarchical structure. En comparing the different reconstruction models, we can say that the method JMAP, which considers a hierarchical problem, has a better property than the other methods. The VBA method also considers the hierarchical problem. The difference is that the VBA method will take into account not only the unknown parameters, but also the uncertainties of the unknowns. The main problem of VBA is to calculate the diagonal elements of the matrix H H H t t t H H H. For a cube of 512 × 512 × 512 voxels, it takes more than 10 days to calculate the diagonals elements. The next step is to optimise the coding of projection and back-projection to reduce the time of calculation and compare these two methods.

4 midle

 4 276.4415 δ g = 794.3637 δ f = 0.2536 δ g = 121 × 10 -δ f = 0.0599 δ g = 4.7671 × 10 -4 δ f = 0.0790 δ g = 7.5320 × 10 -4 δ f = 0.0496 δ g = 3.5128 × 10 -4

FIGURE 1 .

 1 FIGURE 1. Middle slice obtained with methods.

FIGURE 2 .

 2 FIGURE 2. The criterion of different method of reconstruction.

FIGURE 3 .

 3 FIGURE 3. The criterion of different method of reconstruction.