E Gouillart

F Krzakala

M Mézard

L Zdeborová

Belief Propagation Reconstruction for Discrete Tomography

We consider the reconstruction of a two-dimensional discrete image from a set of tomographic measurements corresponding to the Radon projection. Assuming that the image has a structure where neighbouring pixels have a larger probability to take the same value, we follow a Bayesian approach and introduce a fast messagepassing reconstruction algorithm based on belief propagation. For numerical results, we specialize to the case of binary tomography. We test the algorithm on binary synthetic images with different length scales and compare our results against a more usual convex optimization approach. We investigate the reconstruction error as a function of the number of tomographic measurements, corresponding to the number of projection angles. The belief propagation algorithm turns out to be more efficient than the convex-optimization algorithm, both in terms of recovery bounds for noisefree projections, and in terms of reconstruction quality when moderate Gaussian noise is added to the projections.

Introduction

X-ray computed tomography (CT) [START_REF] Slaney | Principles of computerized tomographic imaging[END_REF][START_REF] Herman | Fundamentals of Computerized Tomography: Image Reconstruction from Projections[END_REF] is a classical 3-D imaging technique in materials science or for medical applications. The X-ray beam transmitted through the sample is recorded on a planar detector, for different angles of incidence of the beam on the sample. Using the Beer-Lambert law for photons absorption, the reconstruction of the absorption x of the sample from the measurements y is a linear problem

y = Fx + w, (1)
with F the tomography projection matrix and w the measure noise. For a parallel X-ray beam, each pixel line of the detector measures the total absorption integrated along a horizontal light-ray through the sample. The geometry of tomographic measurements in shown on the schematic of Fig. 1. For a detector line of L pixels, a square image of N = L 2 pixels has to be reconstructed from M = L × n θ measurements, where n θ is the number of angles, also called number of projections.

Classical reconstruction algorithms such as filtered back-projection (FBP) [START_REF] Shepp | The Fourier reconstruction of a head section[END_REF][START_REF] Gn Ramachandran | Three-dimensional reconstruction from radiographs and electron micrographs: application of convolutions instead of Fourier transforms[END_REF], 1] require the linear system (1) to be sufficiently determined. Therefore, the number of angles n θ needed is generally close to the number of pixels along a line, L. Nevertheless, many applications would benefit from being able to reconstruct from a smaller number of angles. In medical imaging for example, reducing the total X-ray dose absorbed by the patient is highly desirable. For in-situ tomography where the sample is evolving while it is imaged [START_REF] Baruchel | Advances in synchrotron radiation microtomography[END_REF][START_REF] Buffiere | In situ experiments with X-ray tomography: an attractive tool for experimental mechanics[END_REF], reducing the number of angles increases the acquisition speed, hence the time-resolution of the observations. Filtered back-projection fares poorly for a limited number of angles, and algorithms incorporating prior information on the sample are needed to improve the quality of the reconstruction.

For example, a positivity constraint on the absorption incorporated in the Algebraic Reconstruction Technique (ART) improves the result of the reconstruction [START_REF] Gordon | Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography[END_REF]. More specific priors can also be used. In particular, samples with a discrete number of phases and absorption values are frequently encountered in materials science. Using this information can make up for the lack of measurements and result in a satisfying reconstruction. The field of discrete tomography [START_REF] Herman | Discrete Tomography: Foundations, Algorithms, and Applications[END_REF] encompasses the large class of algorithms that have been proposed to solve this problem. Various approaches have been proposed, from heuristic methods [START_REF] Batenburg | DART: A practical reconstruction algorithm for discrete tomography[END_REF] to the approximate resolution of the combinatorial problem [START_REF] Carvalho | Binary tomography for triplane cardiography[END_REF], or of its convex relaxation [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF].

In this article, we propose a new algorithm for discrete tomography, that estimates the marginal probability for every pixel value. We use a Bayesian approach in order to write the a-posteriori probability distribution of discrete images, and a fast approximate message passing scheme relying on belief propagation to find the marginal probabilities. Within this framework, constraints imposed by the measurements are satisfied by iteratively solving a graphical model along the light-rays crossing the sample, known as the Potts model in statistical physics. For the specific case of binary images, this amounts to solving 1-D Ising models for all light-rays and measurements. The Potts and Ising models are specific examples of Markov random fields (MRF). In the case of image segmentation, these models can be solved using graph-cuts methods [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF].

We start this article with a short review of the state of the art on discrete tomography, where we introduce the convex algorithm based on total variation, used for comparison with our algorithm. Then we move on to the presentation of our beliefpropagation tomographic reconstruction algorithm (BP-tomo) and its relation to the field of statistical physics. Numerical results on synthetic data are presented for the case of binary images. We establish empirically the recovery phase diagram by decreasing the number of projections and determining the undersampling rate separating successful and failed reconstruction; it is found that the transition corresponds to an undersampling rate of the order of the density of gradients in the image. This bound on undersampling rates is a significant improvement over performances of traditional compressed sensing. In the realistic case of noisy measurements, reconstructions obtained with BP-tomo are compared with results obtained with a constrained total-variation minimization. For moderate noise, BP-tomo is found to give a much more accurate reconstruction above the critical undersampling rate of the noise-free case.

State of the art in discrete tomography

Discrete tomography algorithms perform at the same time the tomographic reconstruction and the segmentation of the reconstructed image into objects with a finite number of known absorption values. In addition to the constraints on the discrete pixel levels, that strongly restrict the space of possible images, additional constraints on the gradients of the image are often enforced in order to obtain smooth images. A large range of methods has been proposed for discrete tomography, for a review see [START_REF] Herman | Discrete Tomography: Foundations, Algorithms, and Applications[END_REF]. Here, we only cite a few representative examples, in order to discuss briefly their scope of application and performance.

Theoretical studies [START_REF] Gardner | Geometric tomography[END_REF][START_REF] Herman | Discrete Tomography: Foundations, Algorithms, and Applications[END_REF] have addressed the problem of the minimal number of angles needed to reconstruct a binary image, in the noise-free case. An important result of geometric tomography [START_REF] Gardner | Geometric tomography[END_REF] is that four different projections are enough to reconstruct uniquely a uniform convex object [START_REF] Gardner | On Hammer's X-ray problem[END_REF]. For more complicated images, degenerate binary solutions may exist for the same set of projections, but it is possible to bound the distance between two solutions using the projection error of the thresholded filteredback projection solution [START_REF] Batenburg | Bounds on the difference between reconstructions in binary tomography[END_REF]. It is also possible [START_REF] Aharoni | Binary vectors partially determined by linear equation systems[END_REF] to determine the set of pixels uniquely determined (in the absence of noise) by their projections.

As for reconstruction algorithms, a first class of methods [START_REF] Carvalho | Binary tomography for triplane cardiography[END_REF][START_REF] Liao | Automated estimation of the parameters of gibbs priors to be used in binary tomography[END_REF][START_REF] Liao | A coordinate ascent approach to tomographic reconstruction of label images from a few projections[END_REF][START_REF] Liao | A method for reconstructing label images from a few projections, as motivated by electron microscopy[END_REF][START_REF] Mohammad-Djafari | Gauss-markov-potts priors for images in computer tomography resulting to joint optimal reconstruction and segmentation[END_REF] adopts a Bayesian approach in order to make the most of the available information on the sample. The a posteriori probability distribution of a discrete image x given measurements y is

P (x|y) = P (y|x)P 0 (x) P (y) , (2)
with P 0 (x) the a priori probability distribution on the space of images. In the Maximum A Posteriori (MAP) setting, the image maximizing the posterior distribution is searched for:

x MAP = argmax P (x|y).

Finding x MAP is an NP-hard combinatorial problem. An MCMC (Monte-Carlo Markov Chain) algorithm is used in [START_REF] Carvalho | Binary tomography for triplane cardiography[END_REF][START_REF] Liao | Automated estimation of the parameters of gibbs priors to be used in binary tomography[END_REF][START_REF] Liao | A coordinate ascent approach to tomographic reconstruction of label images from a few projections[END_REF][START_REF] Liao | A method for reconstructing label images from a few projections, as motivated by electron microscopy[END_REF] in order to sample the Gibbs distribution and to obtain an approximation of x MAP . Gibbs sampling being very costly in terms of numbers of iterations, such algorithms are limited to small images. Therefore, these methods are more suitable for techniques with lower spatial resolution than computed tomography, such as positron emission tomography (PET) for example. For faster computations, a variational Bayes approach was used in [START_REF] Mohammad-Djafari | Gauss-markov-potts priors for images in computer tomography resulting to joint optimal reconstruction and segmentation[END_REF]. The Bayesian approach is very flexible in terms of prior distribution on the images, allowing for a tuned weighting of local binary patterns in the image [START_REF] Liao | Automated estimation of the parameters of gibbs priors to be used in binary tomography[END_REF]. For binary images, a popular model for smooth images is the Ising model from statistical physics, that we shall use as well in our belief-propagation algorithm.

A second class of algorithms [START_REF] Delaney | Globally convergent edge-preserving regularized reconstruction: an application to limited-angle tomography[END_REF][START_REF] Weber | A linear programming relaxation for binary tomography with smoothness priors[END_REF][START_REF] Schüle | Discrete tomography by convex-concave regularization and DC programming[END_REF][START_REF] Weber | A benchmark evaluation of large-scale optimization approaches to binary tomography[END_REF][START_REF] Weber | Binary tomography with deblurring[END_REF] consider relaxations of the MAP estimation falling into the scope of linear optimization. In this setting, the searched-for image x is found by solving a set of linear constraints, such as the constraint y = Fx, or an interval constraint. For mixtures of convex and concave constraints, methods from DC (Difference of Convex functions) programming are used to converge to a local minimum of the objective function [START_REF] Schüle | Discrete tomography by convex-concave regularization and DC programming[END_REF][START_REF] Weber | Binary tomography by iterating linear programs from noisy projections[END_REF][START_REF] Weber | A benchmark evaluation of large-scale optimization approaches to binary tomography[END_REF][START_REF] Weber | Binary tomography with deblurring[END_REF].

Since the seminal papers by Candes, Romberg, Tao and Donoho [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF][START_REF] Donoho | Compressed sensing. Information Theory[END_REF], the field of compressed sensing has breathed a new lease of life into the use of convex optimization for finding sparse solutions to underdetermined inverse problems. Compressed sensing proposes to find sparse solutions using a regularization with the convex ℓ 1 norm, instead of the ℓ 0 norm measuring sparsity: problems can then be solved using classical techniques of convex optimization [START_REF] Boyd | Convex optimization[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. In general, convex optimization methods can be computationally less intensive than other types of methods, hence more suitable for large images. Candes et al. [START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF] proved under certain hypotheses that perfect reconstruction can be obtained with the ℓ 1 norm in the noiseless case, and that the error is at most proportional to the noise level. These solid mathematical foundations have contributed greatly to the success of compressed sensing. Although tomographic projection does not satisfy the restrictive hypotheses of the compressed sensing setting, compressedsensing techniques have proven to work well for tomographic reconstruction. In fact, the founding paper of compressed sensing [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] demonstrated perfect reconstruction with ℓ 1 regularization using the example of incomplete tomographic projections of the piecewise-constant Shepp-Logan phantom [START_REF] Shepp | The Fourier reconstruction of a head section[END_REF]. For the piecewise-constant images of discrete tomography, the regularization is on the ℓ 1 norm of the gradient of the image, that is its total variation (TV)

TV(x) = N i=0 |∇x i | , (4)
where ∇x i is the discrete gradient of the image evaluated at pixel i. Therefore, several papers [START_REF] Song | Sparseness prior based iterative image reconstruction for retrospectively gated cardiac micro-CT[END_REF][START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF][START_REF] Gt Herman | Image reconstruction from a small number of projections[END_REF][START_REF] Tang | Performance comparison between total variation (tv)-based compressed sensing and statistical iterative reconstruction algorithms[END_REF][START_REF] Sidky | Image reconstruction exploiting object sparsity in boundary-enhanced X-ray phase-contrast tomography[END_REF][START_REF] Jia | Gpu-based fast cone beam ct reconstruction from undersampled and noisy projection data via total variation[END_REF][START_REF] Anthoine | On the efficiency of proximal methods in cbct and pet[END_REF] have used the following convex minimization to reconstruct discrete tomography images:

x TV = argmin x 1 2 ||y -Fx|| 2 + βTV(x) (5)
where ||y -Fx|| 2 is the ℓ 2 data fidelity term (corresponding to a Gaussian noise prior in the MAP setting), and the parameter β controls the amount of regularization. It has been shown [START_REF] Song | Sparseness prior based iterative image reconstruction for retrospectively gated cardiac micro-CT[END_REF][START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF][START_REF] Gt Herman | Image reconstruction from a small number of projections[END_REF][START_REF] Tang | Performance comparison between total variation (tv)-based compressed sensing and statistical iterative reconstruction algorithms[END_REF][START_REF] Sidky | Image reconstruction exploiting object sparsity in boundary-enhanced X-ray phase-contrast tomography[END_REF][START_REF] Jia | Gpu-based fast cone beam ct reconstruction from undersampled and noisy projection data via total variation[END_REF][START_REF] Anthoine | On the efficiency of proximal methods in cbct and pet[END_REF] that the total variation reconstruction greatly improves the quality of the reconstruction of discrete images for a reduced number of projections. There is often a trade-off (controlled by the parameter β) between removing the noise and smoothing out small features [START_REF] Tang | Performance comparison between total variation (tv)-based compressed sensing and statistical iterative reconstruction algorithms[END_REF]. For enforcing spatial regularization, it is also possible to use frame representations [START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF][START_REF] Pustelnik | Parallel algorithm and hybrid regularization for dynamic pet reconstruction[END_REF] -especially for medical images that are often more complex than piecewise-constant -or to combine frames and total variation [START_REF] Pustelnik | Parallel algorithm and hybrid regularization for dynamic pet reconstruction[END_REF].

In this paper, we use a total-variation minimization algorithm to benchmark the quality of the reconstruction of the belief-propagation algorithm. Compared to Eq. ( 5), we also constrain the solution to have values between the minimal discrete value a and the maximal value b [START_REF] Pustelnik | Parallel algorithm and hybrid regularization for dynamic pet reconstruction[END_REF]. We solve

x cTV = argmin x 1 2 ||y -Fx|| 2 + βTV(x) + I [a,b] (x), (6)
where I [a,b] (x) is the multivariate indicator function of the interval [a, b]. In order to solve the minimization problem, we use proximal methods [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] that perform subgradient iterations on non-differentiable terms. Several methods [START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF][START_REF] Raguet | Generalized forward-backward splitting[END_REF] are suited for a sum of non-differentiable functionals (such as I [a,b] and TV), we use here the generalized forward-backward splitting [START_REF] Raguet | Generalized forward-backward splitting[END_REF]. We have verified experimentally that with the additional interval constraint, better reconstruction results are obtained with Eq. ( 6) than with Eq. ( 5). For the total variation proximal operator, we use the second-order FISTA scheme of Beck and Teboulle [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] for the isotropic total variation. Finally, several heuristic algorithms have been proposed for discrete tomography. One of the simplest [START_REF] Myers | The binary dissector: phase contrast tomography of two-and three-material objects from few projections[END_REF] consists in alternating gradient iterations on the data fidelity term and clipping (thresholding) continuous values on the discrete labels. Batenburg [START_REF] Batenburg | A network flow algorithm for reconstructing binary images from discrete X-rays[END_REF] proposed a network-flow algorithm where constraints imposed by pairs of directions are satisfied iteratively for binary pixel values. For binary image, Batenburg also introduced DART [START_REF] Batenburg | DART: A practical reconstruction algorithm for discrete tomography[END_REF], a variation on the gradient descent -clipping method: pixels are progressively labeled from the interior of objects to the more uncertain boundaries. This algorithm takes into account the spatial continuity of the objects as well as their binary nature. Along the same lines, another algorithm was recently proposed by Roux et al. [START_REF] Roux | Tomographic reconstruction of binary fields[END_REF], where the reconstruction is done for the continuous probability of pixels values, reducing the risk of clipping too early the value of a pixel.

Compared to the existing literature on discrete tomography, we shall see that our belief propagation algorithm combines the full power of the Bayesian approach and the efficiency of iterative reconstruction algorithms.

Image reconstruction as a statistical physics problem

The Bayesian approach

Let us first fix the notations. The image that we want to reconstruct sits on a two-dimensional grid, and each of its N pixels i takes a value in a q-ary alphabet:

x i = χ 1 , . . . , χ q . For instance, binary tomography corresponds to the case where q = 2, and χ 1 = -1, χ 2 = 1. We suppose that M tomographic measurements, µ = 1 . . . M, have been made. Each such measurement gives a number y µ which is the sum of variables along a light ray µ:

y µ = i∈µ x i + w µ ( 7
)
where i ∈ µ denotes the set of all variables x i that belong to the measurement line µ, and w µ is the additive noise on measure µ. Notice that, in a slightly more general setting, one would measure y µ = i∈µ F µi x i + w µ , where F µi is a geometrical factor giving the fraction of pixel i that is covered by the ray corresponding to measurement µ [START_REF] Joseph | An improved algorithm for reprojecting rays through pixel images[END_REF]. This more general setting can be handled as well by our method. Furthermore, it does not make much difference when the image has a structure in domains, as soon as the typical domain size is much larger than one pixel. To keep notations simple we thus keep hereafter to the case where F µi = 1, the straightforward extension to more general cases is left to the reader. The problem is now to reconstruct the original image (that is the N-dimensional vector x) from the values of the projection (the M-dimensional vector y).

Here, we shall adopt a Bayesian approach. As a prior information, we use the fact that the image to measure is not a random one, but has a structure in domains, such that if a pixel i has the value χ r , then the neighboring pixel has a large probability to have the same value. We shall enforce this bias in each direction µ in which a measurement has been made. This is a very natural prior for a large class of images. Following the Bayesian approach, our goal is thus to sample from the a posteriori distribution given by

P (x|y) = P (y|x)P 0 (x) P (y) = 1 Z(y) µ δ y µ - i∈µ x i P 0 (x) (8)
where, following the convention in statistical physics, we have rewritten the normalization factor of P (x|y) as Z(y). (In the following, we shall incorporate all normalization factors into Z.) Here the δ function implements the equality y µ = i∈µ x i if the measurement is perfect (w µ = 0). In the case of noisy measurement, this δ function should be smoothed in order to take into account our information on the noise. For instance, when the measurement has Gaussian noise of variance ∆ one could use a function

δ(x) = exp(-x 2 /(2∆))/ √ 2π∆.
Let us denote by p µ the a priori probability for neighboring pixels to have the same value along the line µ. Then

P 0 (x) = µ (ij)∈µ p δx i ,x j µ (1 -p µ ) 1-δx i ,x j = µ (1 -p µ ) Nµ µ e (ij)∈µ log pµ 1-pµ δx i ,x j = cst µ e (ij)∈µ Jµδx i ,x j (9)
where N µ is the number of pixels along line µ, (ij) ∈ µ denote pairs of neighboring pixels along direction µ (for instance, in the case of line µ = 3 in Fig. 2, there would be an interaction between pixels (i) and (ii) and another one between (ii) and (iii)) and cst is a constant that does not depend on the values of the variables, and that we shall thus incorporate into the normalization factor through a redefinition of Z(y). In ( 9) we have introduced the so-called "coupling constant" J µ = log pµ 1-pµ . With these notations, the a posteriori probability can thus be written as a product of M terms, each corresponding to one tomographic projection:

P (x|y) = 1 Z M µ=1 δ y µ - i∈µ x i e Jµ (ij)∈µ δx i ,x j (10)
The value of J µ can be fixed based on previous knowledge of the spatial structure of the sample. It would also be possible to learn its value in the course of the algorithm, by maximizing the probability of J µ given the measures y (this amounts to maximizing the partition function Z in Eq. ( 10)). This can be done for instance using the Expectation Maximization approach [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF].

A message passing approach: Belief Propagation

Sampling from the distribution in Eq. ( 8) is a notoriously difficult problem [START_REF] Carvalho | Binary tomography for triplane cardiography[END_REF][START_REF] Herman | Discrete Tomography: Foundations, Algorithms, and Applications[END_REF]. We shall here use an approach based on a message-passing algorithm called "Belief-Propagation" (BP) [START_REF] Pearl | Reverend bayes on inference engines: A distributed hierarchical approach[END_REF][START_REF] Kschischang | Factor graphs and the sum-product algorithm[END_REF][START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF][START_REF] Mézard | Information, Physics, and Computation[END_REF]. It is not an exact algorithm (at least for the present problem), in the sense that it is not guaranteed to sample correctly from the distribution [START_REF] Herman | Discrete Tomography: Foundations, Algorithms, and Applications[END_REF]. However, we shall see that the approximate sampling obtained from BP performs empirically very well.

Let us first reformulate the problem as a graphical model [START_REF] Pearl | Reverend bayes on inference engines: A distributed hierarchical approach[END_REF][START_REF] Mézard | Information, Physics, and Computation[END_REF]. We have N variables x i and M factor nodes, each one implementing the weight δ(y µ -i∈µ x i )e Jµ (ij)∈µ x i x j . The corresponding factor graph is shown in Fig. 2.

We shall not review here the derivation of the belief-propagation algorithm, and refer instead to classical books and articles on the subject [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF][START_REF] Mézard | Information, Physics, and Computation[END_REF]. There are essentially two equivalent ways to derive the BP recursion: one is to see it as a recursion that would be exact if the factor graph were a tree (i.e. if it had no loop). Another one is to see it as an Ansatz in the usual variational Bayesian inference approach that usually improves on the factorized "mean-field" one, leading to an expression for the log-likelihood whose maximization leads to the belief-propagation equations. The result is a set of message passing equations relating messages that go from variables to nodes and from nodes to 8). Here, for instance, the variables denoted (i),(ii) and (iii) are all linked to the factor node number 3.

! ! EE E EEE EE E EEE
variables factor nodes µ =1,...,M

x 1,...,N

mµ → i ( x i ) m j → γ ( x j ) Figure 3
. Message passing on the graphical model. Each variable j = 1, . . . , N sends a message m j→γ (x i ) to all the factor nodes γ that are connected with it, and each factor node µ = 1, . . . , M sends a message mµ→i (x i ) to all variables i that are connected with it.

variables. The message m i→γ from a variable i to a node γ is the marginal probability of the variable x i in absence of the constraint γ (see Fig. 3). On the other hand, the message mµ→i from the factor µ to the variable i is the marginal probability of the variable x i when only constraint µ is present (see Fig. 3). The message m i→γ is easily written in terms of the messages mµ→i :

m i→γ (x i ) ∝ µ∈i =γ mµ→i (x i ) (11)
where µ ∈ i denotes all the factor nodes (tomographic lines) that contain variable i.

Note that m i→γ (x i ) should be normalized , in the sense that x m i→γ (x) = 1 (we prefer however not to write explicitly the normalisation factor and will therefore use the sign ∝ instead of =). The message from a factor node to a variable, on the other hand, reads:

mµ→i (x i ) ∝ x j ;j∈µ =i δ(y µ - j∈µ x j )e Jµ (jk)∈µ δx j ,x k j∈µ =i m j→µ (x j ) (12)
Equations ( 11) and ( 12) build a closed set of equations that should be solved for the values of the messages. Usually one seeks a solution by iterating these equations starting from some randomly chosen initial condition. Assuming that a fixed point of this iteration is reached, one obtains a set of messages that solve the BP equations. The BP estimate for the marginal distribution of x i is then given by

m i (x i ) ∝ µ∈i mµ→i (x i ) (13)
Let us now take a closer look at the recursion: While the iteration from variables to factors (in eq. ( 11)) is easy to compute (a simple multiplication), the one from nodes to variables (in eq. ( 12)) is unfortunately more complicated, as it involves a sum over all the q possible values of each variable x in the line µ. Moreover, we have to do it for each of the N µ variables in this line, and for each value it can take. At first sight, this seems to give a complexity of order N µ q Nµ for each iteration. However, one can compute the message in [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] in a much more efficient way, namely in O(qN µ ) iterations. Indeed, in physics language, one can recognize that estimating mµ→i (x i ) in Eq. ( 12) is nothing else that estimating the marginal of x i in a one-dimensional graphical model known in physics as the Potts ferromagnet in a random field at fixed total magnetization, which can be done efficiently by transfer matrix techniques, which are equivalent for this problem to using BP. The solution is thus to use BP in order to compute the messages of the initial BP problem.

Computing messages: belief propagation within belief propagation

Let us now focus on the following sub-problem: we want to estimate the marginal of a variable x i coupled in a chain of N µ variables, with a statistical weight given by

P µ→i (x) ∝ δ(y µ - j∈µ x j )e Jµ (jk)∈µ δx j ,x k j∈µ =i m j→µ (x j ) . (14)
m i→µ (x i )

m j→µ (x j ) η R j→j-1 (x j ) η L i→i+1 (x i )
x i x j There are two methods to deal with the delta function which we shall refer to (using the statistical physics language) as the "canonical method" and the "micro-canonical method". The exact solution is to use the microcanonical method where one is indeed summing only over the configurations of the variables that satisfy the delta function in Eq. ( 14). This can be done efficiently using a transfer matrix approach, and we have thus implemented this method. It turns out, however, that the canonical method is faster. It is also a more natural approach when there exist errors or noise in the measurement process. We shall therefore use the latter method in the following.

The canonical method amounts to relaxing the δ(y µ -j∈µ x j ) constraint, and treating it only as an approximate constraint that should hold on average, up to some small fluctuations. More precisely, in the canonical method we replace the delta function using a Lagrange multiplier H (the "field" in physics language) that we shall fix later on in order to enforce the constraint on average. The statistical weight thus becomes, instead of ( 14):

P H µ→i (x) ∝ e -H(yµ-j∈µ x j ) e Jµ (jk)∈µ δx j ,x k j∈µ =i m j→µ (x j ) , . (15)
The value of H is computed using E H ( j∈µ x j ) = y µ , where E H denotes the expectation value with respect to the measure P H . Of course, the value of H depends on both µ, and i. Now, it turns out that one can compute H with an error that becomes negligible (in the limit where the number of variables involved in the measurement µ is large), by considering the complete problem obtained from ( 15) by including all the messages m j→µ . This complete problem is defined by:

P H µ (x) ∝ e -H(yµ-j∈µ x j ) e Jµ (jk)∈µ δx j ,x k j∈µ m j→µ (x j ) . ( 16
)
We now write the nested belief propagation equations for this complete sub-problem. Since the interaction is a two-body one, one can write it using a version that involves only one type of message, from variables to variables (see [START_REF] Zdeborová | Phase transitions in the coloring of random graphs[END_REF][START_REF] Mézard | Information, Physics, and Computation[END_REF] and Fig. 4). It is also equivalent to the so-called transfer matrix approach in statistical mechanics [START_REF] Mézard | Information, Physics, and Computation[END_REF]. The recursion for the BP messages is easily written. In order to keep notations simple, let us assume that the numbering of variables is such that variable i is neighbour to i + 1 along direction µ. We shall denote the messages going from right to left of the line as the right-moving messages η R i→i-1 (x i ), i = 2, . . . , N µ and the ones going from left to right as the left-moving messages η L i→i+1 (x i ), i = 1, . . . , N µ -1 (see Fig. 4). η R i→i-1 (x i ) is the marginal of x i in the absence of the left part of the chain beyond x i . Then one has (for the left-moving messages):

η L i→i+1 (x i ) = m i→µ (x i ) Z e Hx i   η L i-1→i (x i )e Jµ + x i-1 =x i η L i-1→i (x i-1 )   = m i→µ (x i ) Z e Hx i 1 + η L i-1→i (x i )(e Jµ -1) , (17)
where we have used that x i η L i→i+1 (x i ) = 1 and Z is the corresponding normalization factor. The recursion for the right-moving messages is similar, with the i + 1 and i -1 playing reverse roles.

Once all the η L and η R are computed, one can then obtain the desired mµ→i (x i ), taking care of the fact that the message m i→µ (x i ) should not be included (see eq. 12), using mµ→i

(x i ) = e Hx i Z 1 + η L i-1→i (x i )(e Jµ -1) 1 + η R i+1→i (x i )(e Jµ -1) , (18)
where Z is again a normalization.

Finally, the expected value of the variable y µ = i∈µ x i can be computed using the average value of x i in presence of the random field (and this time therefore including the term m i→µ ) given by

E(x i ) = x i x i mµ→i (x i )m i→µ (x i ) x i mµ→i (x i )m i→µ (x i ) , (19)
from which we can compute the expected value of the sum of the variables M(H) along the line µ:

M(H) = E i∈µ x i = i∈µ E(x i ) . (20)
A standard result of statistical physics is that M(H) is a monotonous function of H. This is a consequence of the fact that there are no phase transitions in a one-dimensional Potts model [START_REF] Gabor | Mathematical Physics in One Dimension: ExactlySoluble Models of Interacting Particless[END_REF][START_REF] Landau | Statistical Physics[END_REF]. This means that one can apply safely a dichotomy algorithm to find the value of H such that y µ = M(H).

With the canonical method, one cannot fix exactly i∈µ x i but rather its expectation, and this means that we are allowing fluctuations (typically of order O( √ y µ )) around the strict tomographic constraint. Of course, one might also want to enforce strictly this constraint, in which case one has to resort to the microcanonical method.

Given these equations, the complete marginals for each variable can be computed using Eq. [START_REF] Gardner | Geometric tomography[END_REF]. The marginal m i (x i ) denotes the probability that variable i is assigned the value x i ; and the most likely labeling is thus given, at each step of the algorithm, by

x * i = argmax x i m i (x i ). ( 21
)
After convergence of the algorithm, the values x * i thus represent the segmentation into discrete labels that is looked for.

A special case: binary tomography

We shall now specialize for the rest of this paper to binary tomography, that is, the case q = 2, where the variables take values x i = -1, 1 (after a suitable rescaling of experimental data) ‡. In physics terminology, this amounts to the Ising model, and we shall sometimes use the term "spins" for pixels. In this case, it is convenient to rewrite the expressions in a simpler way. Indeed, for binary variables x one can write a probability using a real variable h (which we will call the field) via:

P (x) = e hx 2 cosh h . (22)
Therefore, we define the fields h i→µ and hµ→i by

m i→µ (x i ) = e h i→µ x i 2 cosh h i→µ , (23) mµ→i
(x i ) = e hµ→i x i 2 cosh hµ→i . (24)
Using this notation, we can rewrite the message from a variable to a node eq. ( 11) as follows:

h i→γ = µ∈i =γ hµ→i . (25)
The message from a node to a variable, on the other hand, reads: mµ→i (x i ) ∝

x j ;j∈µ =i δ(y µ -j∈µ x j )e Jµ (jk)∈µ x j x k j∈µ =i

m j→µ (x j ) (26)

∝

x j ;j∈µ =i δ(y µ -j∈µ x j )e Jµ (jk)∈µ x j x k + j∈µ =i x j h j→µ .

The marginal mµ→i (x i ) thus corresponds to the one of a one-dimensional random field Ising model, where the external fields are given by h j→µ , and constrained to have a total magnetization y µ . The recursion "within a line" (corresponding to Eq. ( 17) in the general case) can also be written in terms of fields, using η L (x) = e u L x 2 cosh u L and η R (x) = e u R x 2 cosh u R , we have

u L i→i+1 = atanh tanh J µ tanh (H + h i→µ + u L i-1→i ) , (28)
u R i→i-1 = atanh tanh J µ tanh (H + h i→µ + u R i+1→i ) . (29)
Again, once these are computed, the marginal can be expressed as mµ→i

(x i ) = e hµ→i x i 2 cosh hµ→i , ( 30
) hµ→i = u L i-1→i + u R i+1→i + H . ( 31
)
‡ Even if the exact values of the discrete levels are not known, a method for estimating the discrete values has been proposed by Batenburg et al. [START_REF] Kj Batenburg | A semi-automatic algorithm for grey level estimation in tomography[END_REF] Finally, the expected value of the variable y µ = i∈µ x i , that is needed to fix the correct value of H, is given by summing the average value of x i (with all random fields), each of them given by

E(x i ) = tanh (h i→µ + hµ→i ) , (32)
which is the analog of eq. ( 19) for the binary case. Finally, the most likely assignment for the values x i at each step of the algorithm is given by

x * i = sign tanh µ∈i hµ→ı = sign µ∈i hµ→ı . ( 33
)

Implementation of the binary algorithm

We shall now summarize the complete algorithm for the binary case. Let us start by a few comments about the practical implementation.

• An efficient way to initialize the fields hµ→i is to compute their equilibrium value for J = 0:

hµ→i = atanh (y µ /N µ ) , (34)
where N µ is the number of variables on the line µ. This initialization is also used in [START_REF] Roux | Tomographic reconstruction of binary fields[END_REF].

• Along a light ray slanted with respect to the pixel grid, the distance between successive spins can be larger than 1. In this case, we use a coupling between the spins that decreases when the distance between neighbors increases. We use J µ = atanh tanh D J where D is the ℓ 1 distance (Manhattan or taxi-cab distance) between successive variables on the line µ. This choice, suggested by the so-called high-temperature approximation [START_REF] Feynman | Statistical Mechanics: A Set of Lectures[END_REF], mimics the effective correlation between two variables at distance D in the two-dimensional Ising model.

• We have observed that the results x * i depend only weakly on the value of J, as long as J = O(1). Hence, we did not optimize the value of J: all the results presented here were obtained with J = 0.2.

• In the case of noise-free measures, we stop the algorithm when all line-sums constraints are satisfied by the x * i . For noisy measures, we stop the algorithm when the number of flipping spins saturates; this criterion is discussed in Sec. 5.

• We have observed empirically that damping the update of the fields is sometimes necessary for a proper convergence of the algorithm. We thus replaced eq.( 31) by

ht+1 µ→i = s ht µ→i + (1 -s) u L i-1→i + u R i+1→i + H . ( 35
)
In our implementation, we have fixed empirically the damping value s to 1-1.6/n θ . This scaling was found to ensure convergence for all the values of n θ and N tested.

• In order to avoid overflow or underflow problems, we clip the value of the different fields inside an interval [-B, B]. We set B to 400, a value large enough for the result to be independent of B.

• In order to find the value of the external field H, a few dichotomy iterations are enough to find H up to a good precision ǫ (we use ǫ = 0.05). In order to speed up the computation of H, we use Newton's method when the total magnetization is close enough to y µ .

BP-TOMO(y µ , criterion, t max )

1 Initialize the messages h as described above. Update the hµ→i according to eq. ( 35). [START_REF] Gardner | Geometric tomography[END_REF] do Compute the new values x * i from eq.( 33) and check for convergence. 14 return signal components x * i .

Note that this algorithm is of the embarrassingly parallel type, since each line µ can be handled separately when computing line 7 in the above pseudo-code. This could be used to speed up the algorithm by solving each line separately on different cores. A Python code (with some C extensions) of the current version is available on https://github.com/eddam/bp-for-tomo, under a free-software BSD license. 4. Numerical test in the zero-noise case

Generation of synthetic data

In order to measure the performance of our algorithm, we generate synthetic images with different levels of sparsity of the gradients, but a similar binary microstructure. We obtain such a set of images with a controlled level of gradient sparsity as follows: p 2 pixels are picked at random inside the image and set to 1 while the background is set to 0, then a Gaussian filter of width proportional to 1/p is applied to the image. The resulting continuous-valued image is finally thresholded at the mean intensity to obtain a binary image. Examples of binary images with different values of p are shown in Fig. 5. A simple measure of the sparsity of our images is given by

ρ(x) = ℓ(x) N , (36)
where ℓ(x) is the length of the internal boundary of one of the two phases in the image. The internal boundary is the set of pixels that have within their 4-connected neighbours at least one pixel belonging to the other phase §. ρ therefore measures the fraction of pixels lying on the boundary. We choose this measure since the original image can be retrieved from the knowledge of the internal boundary (and the value of a single pixel). It would also have been possible to choose the set of pixels for which the finite-differences gradient of the image is non-zero. However, the latter set of pixels is redundant for retrieving the image, and leads to a slight over-estimation of the sparsity of the image. Using the definition of Eq. ( 36), we find that ρ scales linearly with p.

Reconstruction in the absence of noise

Let us first investigate the effect of the number of angles on the reconstruction of a given image, when no noise is added to the projections: y = Fx. We define the undersampling § The internal boundary can be computed by simple mathematical morphology operations [START_REF] Serra | Image analysis and mathematical morphology[END_REF], as the difference between the origin image and its morphological erosion. Morphological erosion is one of the two fundamental operations in mathematical morphology [START_REF] Serra | Image analysis and mathematical morphology[END_REF], realizing the intuitive idea of eroding binary sets from their boundaries. rate by the ratio between the number of measures and the number of variables:

α = M N . (37)
Since belief propagation is not an exact algorithm, it is not guaranteed to reconstruct the exact original image. Nevertheless, we find experimentally that for a sufficient number of measures, the algorithm always reconstructs the exact image after a finite number of iterations. Once all spins have reached the correct orientation, we observe that their orientation does not change any more. For a given image, the evolution of the discrepancy between the segmentation of the magnetization (Eq. ( 33)) and the ground truth is plotted in Fig. 6 (a) against the number of iterations, for different measurement rates α. The number of errors decays roughly exponentially with time, but the decay rate increases with the number of measurements: the convergence towards the exact image is faster when more measurements are available. Also, a sharp transition is observed for a critical α, under which the number of errors does not reach zero, and eventually increases again at long times. The transition between the two regimes is hard to estimate accurately, since the convergence time seems to diverges when the transition is approached, as shown in Fig. 6 (b). We estimate the critical undersampling rate from the lowest number of angles for which exact reconstruction is reached before 400 iterations.

We have measured the critical undersampling rate α c for images with different sizes of the microstructure, hence different levels of sparsity. The critical undersampling rate is plotted against the boundary density ρ (Eq. ( 36)) in Fig. 6 (c). We observe a linear relationship

α c ≃ ρ. ( 38
)
This simple relationship corresponds to a very good recovery performance, since ρN can be seen as the number of unknowns needed to retrieve the image. Therefore, we only need a number of measurements M = αN comparable to the number of unknowns for an exact reconstruction of the image.

Numerical test in the noisy case

Now that we have established the phase diagram of our algorithm, we wish to assess its performance for noisy measurements. Different sources of noise may corrupt the measurements in X-ray tomography [START_REF] Tu | Noise simulation in cone beam ct imaging with parallel computing[END_REF]. Here we add a Gaussian white noise of fixed amplitude to the projections. We define the measure noise to signal ratio (NSR) by

NSR = σ L , ( 39
)
where σ is the standard deviation of the additive noise.

In the noisy case, it is not possible any more to stop the algorithm when all constraints are satisfied. Nevertheless, for a small value of the noise the algorithm still converges to an exact reconstruction (Fig. 7). When there is convergence to the exact solution, we observe that this solution is stable: all spins keep their orientation during further iterations of the algorithm. Convergence is reached after a number of iterations similar to the noise-free case. For a larger amplitude of the noise, a finite fraction of errors remains. We observe that the fraction of errors first decreases and reaches a minimum, then starts to slowly increase again. A good choice of the stop criterion is therefore important in order to optimize the quality of the reconstruction. Empirically, we found that the best number of iterations correlates well with the number of iterations needed to reach exact reconstruction in the small-noise case. It can also be detected by a change of slope in the decay of spins flipped between successive iterations (see the inset in Fig. 7).

Robustness to noise

In Fig. 8, we have plotted the reconstruction error for the same image against the measure noise to signal ratio, for two values of the measurement rate α = 1/4 and α = 1/10. An image with p = 14 (Fig. 5 left) was used, and J = 0.2 was used for all values of the NSR. We have also computed the reconstruction error for the convex minimization of Eq. ( 6). For a fair comparison, the best value of β minimizing the reconstruction error was computed using Brent's method (whereas no parameter was optimized for BP-tomo). We observe that the reconstruction quality is better for BPtomo than for the convex algorithm, up to a threshold above which the error is greater for BP-tomo. For the smaller undersampling rate α = 1/10, the values of the error reached at the threshold are too high for a satisfying reconstruction for several applications (most pixels lying on boundaries are wrongly labeled). Interestingly, there is a first regime for small noise, for which the fraction of errors is zero or very small for BP-tomo, while the error increases much faster for the convex algorithm. For large noise, however, the convex algorithm gives better results.

Fig. 9 shows that the reconstruction errors are mostly located on the boundaries between the two phases for the two algorithms (except for a few isolated errors for BPtomo), but that a larger fraction of the interfaces are correctly reconstructed for BPtomo. Fig. 9 also displays the continuous magnetization, as well as the non-segmented Stop criterion in the noisy case -Main axes: for α = 1/10, p = 14, fraction of error compared to ground truth vs. number of iterations n for a signal to noise ratio of 0.6 % (•), 1 % ( ) and 2 % (⋆). For a small amplitude of the noise, the number of errors converges to zero (SNR = 0.6%). For a large amplitude of the noise, the fraction of errors reaches a minimum after the number of iterations needed to reach an exact reconstruction for a smaller value of the noise (here, n = 15). Inset: number of flipped spins between two iterations vs. number of iterations. We observe that the minimum fraction or errors correlates well with a change of slope for the number of flipped spins. . Fraction of error N e /L 2 vs. normalized measure noise σ/L, for BP-tomo and the convex algorithm (TV). These numerical tests have been performed for two different undersampling rates α. For a small noise, BP-tomo always outperforms the convex algorithm. In particular, there is a finite interval of noise intensity for which the reconstruction is error-free. 10 -2 10 -1 10 0 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 Figure 10. Fraction or error N e /L 2 as a function of the measurement rate α, for BPtomo and the convex algorithm (TV) (inset in log-log coordinates). These numerical tests have been performed for a small noise σ/L = 0.002. For both algorithms, the number of errors decrease when the measurement rate is increased, but the evolution is very different. For BP-tomo, the evolution is very sharp: above a certain threshold, the reconstruction is error-free. However, below this threshold the number of errors displays a sharp transition and increases quickly to high values. The convex algorithm gives similar results as BP below the BP-tomo threshold, but does not have a "phase transition", so that it is outperformed by BP above the threshold. minimization of the convex algorithm. We see that the contours of the objects are delineated more accurately for BP-tomo than for the convex algorithm, especially for objects with concavities.

For a given value of the NSR, we observe in Fig. 8 that the reconstruction error increases when the measurement rate decreases. The effect of the measurement rate in the noisy case is studied in the next paragraph.

Influence of the number of measures

For the same image (p = 14, J = 0.2), we have fixed the SNR to 0.002 and computed the reconstruction error as a function as the measurement rate, for the two algorithms. Results are shown in Fig. 10. For BP-tomo, we observe a sharp transition between an exact reconstruction to a failed reconstruction when the measurement rate is decreased (a score of 0.5 would be obtained when labeling pixels at random). Therefore, the noisy case displays a transition between failure and success as in the noise-free case (Fig. 6). The transition is observed at the same measurement rate as the noise-free case. Approaching the transition from above, a finite error fraction appears (see the log-log inset in Fig. 10). In the noise-free case, exact reconstruction could be reached until the transition, but an increasing number of iterations was needed when approaching the transition (Fig. 6 (b)).

In contrast, the evolution of the error fraction when decreasing the measurement rate is much smoother for the convex algorithm, and does not display a phase transition in Fig. 10. As a result, the error fraction is larger for the convex algorithm above the BP-tomo transition. At low measurement rates, the error fraction is comparable for the two algorithms, but this regime corresponds to values of the error fraction that are not acceptable for most applications.

For a larger value of the noise, Fig. 11 shows that a transition between good and failed reconstruction is still observed. However, the fraction of errors is larger, therefore the transition is not as sharp as for smaller NSR. Comparing Fig. 10 and Fig. 11, we conclude that the reconstruction quality is better above the transition for a SNR of 0.01 than for the convex algorithm for a SNR five times smaller. undersampling rate to allow for a comfortable margin with respect to the critical noisefree undersampling rate. This will speed up the convergence rate (see Fig. 6 (b)) and reduce significantly the measurement error when the noise is important (compare the plots in Fig. 8 for α = 1/4 and α = 1/10 for 1% of noise and more). In other words, trying to approach the limit α c ≃ ρ is more interesting for theoretical than practical reasons. However, the good news is that one does not need to go very far from the transition α c for an efficient and accurate reconstruction, all the more for a small noise amplitude. A systematic evaluation of the reconstruction error as a function of α, ρ and the NSR is out of the scope of this article, but would be of interest for applications.

In future work, the recovery properties of our algorithm should be tested on real experimental data. In our prototype code, the BP algorithm is a few times (about 5 times) slower than the TV algorithm, principally because of the iterations needed to find H. However, the number of iterations needed to converge is much smaller for the BP algorithm, and we expect the speed of the algorithms to depend very much on the implementation. In order to reconstruct large 3-D volumes in a reasonable time, our prototype implementation of the algorithm can be improved on the numerical sidefor example with an implementation coded on the graphics processing unit (GPU)and probably also on the algorithmic side. For example, a different update scheme for the message passing might accelerate the convergence of the algorithm. For a better precision when the density of boundaries in the image is high (large ρ), the geometrical factors F µi (see Section 3) should be computed, as they are in usual algebraic reconstruction algorithms. For 3-D images, the additional knowledge that neighboring horizontal slices are likely to have the same pixel values can be easily implemented in the algorithm as well. In the same way, a time-regularization in a timeseries of several images could also be used [START_REF] Pustelnik | Parallel algorithm and hybrid regularization for dynamic pet reconstruction[END_REF][START_REF] Chaari | Multidimensional waveletbased regularized reconstruction for parallel acquisition in neuroimaging[END_REF], corresponding to a 3D + time regularization. This should lead to an even better recovery compared to the 2-D case.

We also plan to test the recovery properties in the multilabel case (q > 2). When q becomes large, it remains an open question to know how the discrete algorithm performs compared to the convex continuous algorithm.

Finally, for in-situ tomography [START_REF] Baruchel | Advances in synchrotron radiation microtomography[END_REF] the sample is rotating continuously to speed-up the acquisition rate, in the case of fast transformations. It would therefore be interesting to adapt our algorithm to the case of continuous acquisition, when one image taken on the detector integrates the projection of the sample over a finite angular sector.

Figure 2 .

 2 Figure 2. Construction of the factor graph for the tomographic reconstruction. Each tomographic line in the image (left) corresponds to a factor node in the corresponding graphical model (right). Each of the 1, . . . , M factor nodes imposes the statistical weight in eq. (8). Here, for instance, the variables denoted (i),(ii) and (iii) are all linked to the factor node number 3.

Figure 4 .

 4 Figure 4. Message-passing for the solving the sub-problem of estimating the marginal on each constraint. The left-moving message and the right-moving message should be computed on a one-dimensional line.

Figure 5 .

 5 Figure 5. Synthetic data generated for p = 14, 28, 38, and a size 256 × 256.

Figure 6 .

 6 Figure 6. (a) Decay of the number of errors w.r.t. ground truth vs. the number of iterations of BP-tomo, for p = 15, L = 256 and a number of angles varying from 66 (blue) to 20 (red). The same value of J = 0.2 was used for all cases. Note how the error reduction becomes slower when the number of projections is decreased. For smaller number of projections (not shown here), the fraction of errors saturates at an important level. (b) Number of iterations required to reach an exact reconstruction versus undersampling rate α, for p = 22 and J = 0.2. The number of iterations diverges when the transition between exact reconstruction and faulty reconstruction is approached. (c) Critical undersampling rate α c vs. boundary density ρ.

Figure 7 .

 7 Figure 7. Stop criterion in the noisy case -Main axes: for α = 1/10, p = 14, fraction of error compared to ground truth vs. number of iterations n for a signal to noise ratio of 0.6 % (•), 1 % ( ) and 2 % (⋆). For a small amplitude of the noise, the number of errors converges to zero (SNR = 0.6%). For a large amplitude of the noise, the fraction of errors reaches a minimum after the number of iterations needed to reach an exact reconstruction for a smaller value of the noise (here, n = 15). Inset: number of flipped spins between two iterations vs. number of iterations. We observe that the minimum fraction or errors correlates well with a change of slope for the number of flipped spins.

Figure 8

 8 Figure 8. Fraction of error N e /L 2 vs. normalized measure noise σ/L, for BP-tomo and the convex algorithm (TV). These numerical tests have been performed for two different undersampling rates α. For a small noise, BP-tomo always outperforms the convex algorithm. In particular, there is a finite interval of noise intensity for which the reconstruction is error-free.

Figure 9 .

 9 Figure 9.(a) BP-tomo, α = 1/10, p = 14, reconstructions for σ/L = 0.006, 0.01, 0.02, 0.03 (from left to right). The top row shows the magnetization of the pixels, and the bottom row the segmentation, with segmentation errors contoured in blue (resp. magenta) for wrong pixels in the x = -1 (resp. x = 1) phase. Note how the absolute value of the magnetization decreases when the measure noise increases, revealing a greater uncertainty. (b) TV (convex algorithm), α = 1/10, p = 14, reconstructions for x/L = 0.006, 0.01, 0.02, 0.03 (from left to right). The top row displays the result of the convex optimization, and the bottom row its segmentation. For a low noise amplitude, boundaries between domains are not as sharp as for BPtomo, resulting in more segmentation errors close to boundaries.

Figure 11 .

 11 Figure 11. Fraction or error N e /L 2 as a function of the measurement rate α, for BPtomo and two different values of the noise amplitude (inset in log-log coordinates). The transition is sharper for a small noise. In the inset, the dashed vertical line represents the empirical value of the noise-free critical undersampling rate.

Acknowledgments

The research leading to these results was supported by the ANR program "EDDAM" (ANR-11-BS09-027) and from the European Research Council under the European Union's 7 th Framework Programme (FP/2007-2013)/ERC Grant Agreement 307087-SPARCS. We gratefully acknowledge fruitful conversations with E. Chouzenoux, S. Roux, H. Talbot and G. Varoquaux.

Discussions and perspectives

The belief-propagation algorithm presented in this article is found to have excellent recovery properties for tomographic binary reconstruction. Indeed, in the noise-free case the exact original image can be reconstructed from a number of tomographic measures approximately equal to the number of pixels lying on the internal boundary of the objects, that can be viewed as the number of unknowns in our problem. We are not aware of other studies on discrete tomography where the empirical computation of the critical undersampling rate as a function of the data sparsity was performed. In Batenburg's DART algorithm [START_REF] Batenburg | DART: A practical reconstruction algorithm for discrete tomography[END_REF], a sharp transition between exact and faulty reconstruction is also observable for the different kinds of phantoms studied, but the relation between the complexity of the data and the critical undersampling rate was not described. Our search for recovery bounds is inspired by empirical and theoretical results obtained in the field of compressed sensing about phase transitions in recovery [START_REF] Donoho | Observed universality of phase transitions in high-dimensional geometry, with implications for modern data analysis and signal processing[END_REF][START_REF] Donoho | Precise undersampling theorems[END_REF][START_REF] Krzakala | Statistical-physics-based reconstruction in compressed sensing[END_REF][START_REF] Krzakala | Probabilistic reconstruction in compressed sensing: algorithms, phase diagrams, and threshold achieving matrices[END_REF]. Obtaining theoretical results for the discrete tomography problem is more difficult than for signals sparse in a known basis -our measure of image "sparsity" comes from set theory and mathematical morphology, and applies only to discrete images. Nevertheless, it is possible to obtain empirical bounds as we did in Section 4. For a more complete phase diagram of our algorithm, one should also investigate corrections depending on the size of the data N, that appear for example in the ℓ 0ℓ 1 Donoho-Tanner phase transition [START_REF] Donoho | Observed universality of phase transitions in high-dimensional geometry, with implications for modern data analysis and signal processing[END_REF]. A probability of exact recovery could also be obtained by computing the critical undersampling rate for a large number of images sharing the same sparsity [START_REF] Donoho | Observed universality of phase transitions in high-dimensional geometry, with implications for modern data analysis and signal processing[END_REF][START_REF] Donoho | Precise undersampling theorems[END_REF].

When Gaussian noise corrupts the measurements, we observe that an exact reconstruction can still be retrieved for small noise and a sufficient number of measurements. For a fixed non-zero level of noise, the fraction of segmentation errors grows when the number of measurements is decreased, but the fraction of errors is smaller than the value of the noise over signal ratio for a large window of undersampling rates. At the noise-free critical undersampling rate, the fraction of errors grows very fast when the number of measurements is decreased further. Above the noise-free critical undersampling rate and for moderate noise, we find that our algorithm systematically outperforms a convex algorithm implementing the convex relaxation of the binary tomography problem. The value of the noise corresponding to the crossover when the convex algorithm becomes better increases when the undersampling rates increases (i.e., when one moves away from the transition). At the crossover, the reconstruction error is high when the undersampling rate α is small -typically, most pixels lying on objects boundaries are wrong -, making the reconstruction unsuitable for further processing or interpretation in many cases. We therefore conclude that the belief-propagation algorithm is better in most cases where further processing of the images needs to be done.

For practical applications, if the measurement noise is known and the undersampling rate can be selected, we suggest that the operator should select the