Variational Bayesian Approach in Model-Based Iterative Reconstruction for 3D X-Ray Computed Tomography with Gauss-Markov-Potts prior

Variational Bayesian Approach in Model-Based Iterative Reconstruction for 3D X-Ray Computed Tomography with Gauss-Markov-Potts prior Camille Chapdelaine 1* , Ali Mohammad-Djafari 2, Nicolas Gac2 and Estelle Parra 1 1 SAFRAN SA, Safran Tech, Pôle Technologie du Signal et de l’Information, Magny-Les-Hameaux, France; 2 Laboratoire des signaux et systèmes, CNRS, CentraleSupélec-Univ Paris Saclay, Gif-sur-Yvette, France * Correspondence: camille.chapdelaine@safrangroup.com Version July 3, 2019 submitted to Proceedings Abstract: 3D X-ray Computed Tomography (CT) is used in medicine and non-destructive testing 1 (NDT) for industry to visualize the interior of a volume and control its healthiness. Compared to 2 analytical reconstruction methods, model-based iterative reconstruction (MBIR) methods obtain 3 high-quality reconstructions while reducing the dose. Nevertheless, usual Maximum-A-Posteriori 4 (MAP) estimation does not enable to quantify the uncertainties on the reconstruction, which can 5 be useful for the control performed afterwards. Herein, we propose to estimate these uncertainties 6 jointly with the reconstruction by computing Posterior Mean (PM) thanks to Variational Bayesian 7 Approach (VBA). We present our reconstruction algorithm using a Gauss-Markov-Potts prior model 8 on the volume to reconstruct. For PM calculation in VBA, the uncertainties on the reconstruction are 9 given by the variances of the posterior distribution of the volume. To estimate these variances in our 10 algorithm, we need to compute diagonal coefficients of the posterior covariance matrix. Since this 11 matrix is not available in 3D X-ray CT, we propose an efficient solution to tackle this difficulty, based 12 on the use of a matched pair of projector and backprojector. In our simulations using the Separable 13 Footprint (SF) pair, we compare our PM estimation with MAP estimation. Perspectives for this work 14 are applications to real data as improvement of our GPU implementation of SF pair. 15


18
In 3D X-ray CT, MBIR methods enforce a prior model on the volume to image, so the reconstruction   44 We consider a cone-beam acquisition process : X-rays are sent from a source through the object to control and hit a flat detector which measures the decrease of intensity they have undergone inside the volume. Several perspectives of the volume are acquired by rotating the object around its vertical axis. The M collected measurements g are called the projections and are connected to volume f , of size N, by the linear forward model taking uncertainties into account [16]

Variational Bayesian approach
where H is called the projection operator.
Precisions ρ ζ = (ρ ζ i ) i are assigned Gamma conjugate prior [5] : The prior model on the volume is a Gauss-Markov-Potts prior which consists in labelling each voxel j according to its material z j = k ∈ {1, . . . , K}, where K is the number of materials. Then, the distribution of value f j of voxel j depends on its material z j : Means m = (m k ) k and inverses ρ = (ρ k ) k of variances of the classes have to be estimated and are assigned conjugate priors [5] : A Potts model is assigned to labels z in order to favour compact regions in the volume [5] : denoting by V (j) the neighbourhood of voxel j, we have, according to Hammersley-Clifford theorem [17], From our prior model M, the posterior distribution of unknowns ψ = (f , z, ρ ζ , m, ρ) is given by Bayes' rule [5] p(f , z, ρ ζ , m, ρ|g; M) ∝ p(g|f , ρ ζ )p(f |z, m, ρ)p(z|α, γ 0 )p(ρ ζ |α ζ 0 , β ζ 0 )p(m|m 0 , v 0 )p(ρ|α 0 , β 0 ), (7) where α = (α k ) k . Based on this distribution, JMAP can be performed [5] but does not provide uncertainties on the result. MCMC methods for joint computation of the means and the variances of the posterior distribution are too computationally costly for 3D applications [5,18]. For this reason, we apply VBA which consists in approximating the true posterior distribution p by a simpler distribution q on which posterior means and variances can be easily estimated. Approximating distribution q minimizes Kullback-Leibler (KL) divergence KL(q||p) on a chosen set of simple distributions [12]. The choice we make for q is a factorizable approximation, which only preserves a dependence between value f j of voxel j and its label [19] : Minimizing KL divergence with respect to each factor while fixing the others leads to [13,19] The VBA algorithm turns into the iterative updating of the parameters of these distributions with respect to the others. The updating formulae and the order of their applications are given in [13]. In particular, at iteration t, the variances of the approximating distribution for the volume are updated bỹ . Moreover, the updating formula for intensity parameter of the approximating Gamma distribution for To compute approximate posterior variances, formula (10) needs the computation of diagonal Therefore, in order to implement VBA for 3D X-ray CT, we need to find a way to compute diagonal 48 coefficients in formulae (10) and (11) efficiently. We propose a strategy which is detailed in the next 49 section.

Computation of diagonal coefficients
At one iteration of the algorithm, for any voxel j, diagonal coefficient used to compute v jk by (10) is where e (j) has the size of a volume, formula (13) implies to compute N projections, which is very long, even if the projector implemented on GPU is very fast. We calculated that, if we have to reconstruct a volume of size N = 256 3 voxels from 64 projections of size 256 2 pixels, and if one projection takes only 10 milliseconds, computing all dialgonal coefficients d v j , ∀j, for only one iteration of proposed VBA algorithm [13], would require more than 40 hours. Due to this huge computational cost, we prefer to consider the algebraic formula : From this formula, diagonal coefficients d v appear to be similar to a backprojection ofṽ −1 ζ = (ṽ −1 ζ i ) i , except that coefficients H ij are replaced by their squares H 2 ij , ∀i, j. Similarly, diagonal coefficients appear like a projection of volumeṽ, with H 2 ij instead of H ij . Given formulae (14) and (15)

67
The simulated phantom is of size 256 3 voxels and contains K = 5 classes. It is shown in figure 2. 68 We reconstruct this volume from 64 projections of size 256 2 pixels, uniformly distributed over [0, 2π].

69
These projections are noisy with SNR equal to 20 db.

104
In this paper, we have presented an application for 3D X-ray CT of variational Bayesian approach 105 (VBA) with Gauss-Markov-Potts prior model. By computing posterior mean (PM) thanks to VBA, 106 we have been able to jointly perform the reconstruction and the estimation of the posterior variances, 107 which give the uncertainties on the reconstruction. To compute these variances, we have seen that