Bounds on the quality of reconstructed images in binary tomography

Bounds on the quality of reconstructed images in binary tomography K. Joost Batenburga,c, Wagner Fortesa,b, Lajos Hajdud,e, Robert Tijdemanb aCentrum Wiskunde &Informatica, Science Park 123, 1098 XG Amsterdam, The Netherlands bMathematical Institute, Leiden University, The Netherlands cVision Lab, University of Antwerp, Belgium dInstitute of Mathematics, University of Debrecen, Hungary eNumber Theory Research Group of the Hungarian Academy of Sciences, Debrecen, Hungary Abstract Binary tomography deals with the problem of reconstructing a binary image from its projections. In particular, there is a focus on highly underdetermined reconstruction problems for which many solutions may exist. In such cases, it is important to have a quality measure for the reconstruction with respect to the unknown original image. In this article, we derive a series of upper bounds that can be used to guarantee the quality of a reconstructed binary image. The bounds limit the number of pixels that can be incorrect in the reconstructed image with respect to the original image. We provide several versions of these bounds, ranging from bounds on the difference between any two binary solutions of a tomography problem to bounds on the difference between approximate solutions and the original object. The bounds are evaluated experimentally for a range of test images, based on simulated projection data. Keywords: binary tomography, discrete tomography, error bounds, stability 1. Introduction Tomography is a technique for reconstructing an image of an object from a series of projections of this object, acquired from a range of viewing angles. The projection images are typically recorded using a scanning device, which can employ various types of beams (e.g., Xrays, neutrons, electrons) that traverse the object, after which a detector measures the result of the beam-object interaction. Provided that a large number of high-quality projections are available, sampled from a full range of angles, an accurate reconstruction of the object can be computed using a tomographic reconstruction algorithm [15, 19]. In practice, the set of angles for which projections are acquired is often limited. Due to dose constraints, it can be desirable to record as few projections images as possible, while still attaining sufficient image quality. Also, the angular range can be restricted by the particular scanning setup, such as in electron tomography, where the shape of the sample holder limits the Email addresses: [email protected] (K. Joost Batenburg), [email protected] (Wagner Fortes), [email protected] (Lajos Hajdu), [email protected] (Robert Tijdeman) Preprint submitted to Discrete Applied Mathematics September 10, 2012 angular range of the projections [20]. The resulting image reconstruction problems, based on just a small number of projections, are known as limited data problems. For tomographic reconstruction from severely limited data, classical algorithms based on analytical inversion of the Radon transform, such as the Filtered Backprojection algorithm, often yield low quality reconstructions that are polluted by strong artefacts. In such cases, it makes sense to exploit available prior knowledge of the unknown object. Incorporation of this knowledge in the reconstruction algorithm can potentially result in a reduction of the required number of projections, increased accuracy of the reconstruction, or an improved ability to deal with noisy projection data. A prior that has received much attention recently concerns the sparsity of the image, or of its gradient, which is exploited in the field of Compressed Sensing [12, 13, 23, 24]. A related, but more strict type of prior knowledge is exploited in Discrete Tomography, which focuses on the reconstruction of images that consist of a small, discrete set of grey values [16, 17]. The actual set of grey levels is typically assumed to be known in advance. Here, we focus specifically on the reconstruction of binary images, which consist of just two grey levels, 0 and 1. Several reconstruction algorithms have been proposed in literature for reconstructing binary images from a small number of projections, often even less than 6 [5, 8, 17, 22]. These algorithms do not work well for all types of binary images, yet they have demonstrated the ability to accurately reconstruct various phantom images. Despite the strong constraint imposed on the grey values in discrete tomography, many valid solutions of the reconstruction problem can exist, all corresponding to the same set of projections. If the projections are obtained by performing measurements on some unknown ground truth object, the reconstruction can then deviate substantially from the true object. As a consequence, there is a need for an upper bound on the difference between binary solutions of the reconstruction problem. As the ground truth is a solution by itself, this would also yield a bound on the reconstruction error with respect to the ground truth. A related problem in discrete tomography is the so-called stability problem, which deals with the question of how the reconstruction changes if the projections are slightly perturbed. For images represented on a discrete grid, both lower and upper bounds have been obtained for the magnitude of such changes [1, 2, 3, 11, 28]. For the case of binary image reconstruction from just two projections, horizontal and vertical, sharp bounds on the difference between binary images having the same projections have been obtained by Alpers and Brunetti [2] and extended by Van Dalen [26, 27]. In this article, we present a series of bounds which are highly general. Our bounds can be computed for any set of projections and in different geometrical settings, for lattice images as well as discretized continuous images. A key idea in deriving these bounds is an observation first made by Hajdu and Tijdeman in [14], concerning the fact that all binary solutions of the reconstruction problem must lie on a certain hypersphere, of which both the center and radius can be computed. The center of this hypersphere, which we call the central reconstruction, is the shortest real-valued solution of the tomography problem. This hypersphere construction leads directly to a simple bound on the distance between any two binary solutions, based on the triangle inequality. Stronger bounds can be derived by focusing on the distance between binary solutions of the tomography problem and the binary image that is obtained by rounding the central reconstruction. We derive several bounds that combine properties of the real-valued solution with combinatorial properties that are satisfied by the binary solutions. In particular, the fact that the sum of the pixel values in the unknown image is fully determined by the projection data can be used to improve the error bounds for binary images. The structure of this article is as follows. Our notation is introduced in Section 2, after which 2 a general class of reconstruction problems is introduced. We then discuss two specific examples of such problems, based on the grid model and strip model for computing the projections. In Section 3, the central reconstruction is introduced, which plays an important role in all results that follow. We discuss two different strategies for its computation. Section 4 contains the main results of this paper: a series of quality bounds for binary solutions of the tomography problem. It is divided in three parts: in Section 4.1, a general bound is derived on the difference between two binary images having a given set of projections. Section 4.2 deals with bounds that are based on properties of the binary images that are obtained by rounding the central reconstruction. These bounds are subsequently refined in Section 4.3 by including the knowledge of the total number of 1’s in any binary solution, which can be determined from the projection data. Section 5 presents a series of simulation experiments and their results. From these results, the practical value of the proposed bounds can be evaluated for different types of images. The results are further discussed in Section 6. Section 7 concludes the paper. This article is an extended version of the conference paper [7]. Compared to the conference paper, several new bounds have been included in the present article, in particular (i) bounds on the difference between any binary image, not necessarily being a solution of the reconstruction problem, and the unknown original image, and (ii) bounds that incorporate knowledge of the total number of 1’s in the unknown original image. Various other bounds have also been refined. The experimental section has been extended with an example, depicting the different images and concepts involved in the various bounds. A series of experiments on bounding the error of binary reconstructions computed by the DART algorithm, a heuristic method that is not guaranteed to yield an exact solution [8], but yields high quality reconstructions in practice, is also included. 2. Notation and model Throughout the discrete tomography literature, several imaging models have been considered. In the grid model, an image is formed by assigning a value to each point in a regular grid. In the case of binary images, each point is assigned a value of either 0 or 1. Here, we consider square grids of the form A={(i,j)∈Z2: 1 ≤i,j≤s}for s∈N,s≥1; see Fig. 1(a). For the grid model, we refer to the points in Aas pixels. A binary image defined on Acan be represented by a map A→ {0,1}. A projection of an image fis formed by considering the set of parallel lines through one or more grid points in a certain direction (a,b)∈Z2, with a≥0 and (a,b) coprime, and summing the values of the points on each line. For a line given by the equation ax −by =t(t∈Z), the line projection p is defined as p=X (x,y)∈A:ax−by=t f(x,y). The grid model can be used to model nanocrystals, that consist of discrete atoms positioned in a regular grid [18, 25]. In many tomography applications, a continuous representation of the object is more realistic, as there is no intrinsic grid structure. In such cases, the unknown image is typically approximated by an image defined on a discrete pixel grid, using square pixels. A common model for computing the projections of such a pixelized image is the strip model [19, section 7.4.1], [30]. In the strip model, a projection is computed by considering a set of parallel strips in a given direction and and for each strip computing the weighted sum of all the pixels which intersect that strip with a weight equal to the intersection area of the strip and the pixel. 3 (a) Grid model (b) Strip model Figure 1: Two different projection models. We now define some general notation. An image is represented by a vector x=(xi)∈Rn. We refer to the entries of xas pixels, which correspond to unit squares in the strip model and to points in the grid model. The derivation of our main results does not depend on the particular projection model. Throughout this paper we assume that all images are square, consisting of c rows and ccolumns, where n=c2. A binary image corresponds with a vector ¯x∈ {0,1}n. For a given set of kprojection directions, the projection map maps an image xto a vector p∈ Rmof projection data, where mdenotes the total number of line measurements. As the projection map is a linear transformation, it can be represented by a matrix W=(wi j)∈Rm×n, called the projection matrix. Entry wi j represents the weight of the contribution of xjto projected line i. Note that for the grid model the projection matrix is a binary matrix, while for the strip model its entries are real values in [0,1]. The projection matrix Wand vector pcan be decomposed into k blocks as W=            W1 . . . Wk             ,p=            p1 . . . pk             ,(1) where each block Wd(d=1,...,k) represents the projection map for a single direction and each block pdrepresents the corresponding projection data. From this point on, we assume that the projection matrix has the property that Pm i=1wi j =k for all j=1,...,n. This property is certainly satisfied for the grid model, as every xjis counted with weight 1 for exactly one line in each projection direction. The property is also satisfied for the strip projection model, as the total pixel weight for each projection angle is equal to the area of a pixel, which is 1. For most other projection models commonly used in tomography, such as the line model, where the weight of a pixel is determined by the length of its intersection with a line, this property is approximately satisfied, but not always exactly. The general reconstruction problem consists of finding a solution of the system Wx =p for given projection data p, i.e., to find an image that has the given projections. In binary tomography, one seeks a binary solution of the system. For a given projection matrix Wand given projection data p, let SW(p)={x∈Rn:Wx =p}, the set of all real-valued solutions corresponding with the projection data, and let ¯ SW(p)=SW(p)∩ {0,1}n, the set of binary solutions of the system. As the main goal of incorporating prior knowledge of the binary grey levels in the reconstruction is to reduce the number of required projections, we focus on the case where mis small with respect to n, such that the real-valued reconstruction problem is severely underdetermined. Despite the strong constraint that each pixel value ximust belong to the set {0,1}, the binary reconstruction problem often does not have a unique solution. Instead of posing the uniqueness 4 problem as a yes/no question, we focus in this paper on the number of different pixels between different binary solutions. For example, if the reconstruction problem has no unique solution, but all pairs of solutions have at most 4 different pixels, then one can say that any such solution must be close to the original object from which the projections have been obtained, even if the exact set of differences cannot be determined. For any two vectors ¯ u,¯ v∈ {0,1}n, define the difference set D(¯ u,¯ v)={i: ¯ui,¯vi}and the number of differences d(¯ u,¯ v)=#D(¯ u,¯ v), where the symbol # denotes the cardinality operator for a finite set. Note that d(¯ u,¯ v)=k¯ u−¯ vk1. 3. The central reconstruction As the projection matrix is typically not a square matrix, and also does not have full rank, it does not have an inverse. Recall that the Moore-Penrose pseudo inverse of an m×nmatrix Ais an n×mmatrix A†, which can be uniquely characterized by the two geometric conditions A†b⊥ N(A) and (I−AA†)b⊥ R(A) for all b∈Rm, where N(A) is the nullspace of Aand R(A) is the range of A, [10, page 15]. Let x∗=W†p. Then x∗has the property (see Chapter 3 of [9]) that it is the real-valued solution of minimal Euclidean norm of the system Wx =p, provided that the latter system is solvable. We call x∗the central reconstruction of p. The central reconstruction plays an important role in the bounds we derive for the binary reconstruction problem. We will show in the next section that all binary solutions of the system have equal distance to x∗, so that one can consider the central reconstruction as lying “in the middle” of all binary solutions. As all bounds presented in this paper depend on x∗, accurate computation of x∗is necessary to compute the corresponding difference bounds. One approach to computing the central reconstruction of a consistent system Wx =pis to use the QR decomposition of WT. We will only sketch the computation here and refer to [4] for details. For clarity of presentation, we assume that Whas full row rank. In fact, this assumption is not satisfied for tomography, and the extended QR decomposition should be used. The QR decomposition factorizes the matrix WTinto an orthogonal matrix Qand an uppertriangular matrix Rof full column rank, as WT=Q R 0!. The central reconstruction is then given by x∗=Q(RT)−1p, which can be computed efficiently by first solving the system RTy=pfor yby back substitution, and then computing x∗=Qy. However, due to the size of the matrix W, calculation of the QR decomposition is usually unpractical for large images. For the case m<n, the QR decomposition requires O(n3) operations. Moreover, the n×nmatrix Qis typically dense, requiring a vast amount of computer memory. As an alternative, an iterative method for solving the system Wx =p, called CGLS (Conjugate Gradient Least Squares), can be used [21]. The CGLS algorithm can effectively exploit the sparse structure of the projection matrix to reduce the required computation time, and does not require storage of large, dense matrices. Apart from numerical errors, applying CGLS to the system Wx =presults, after convergence, in the computation of W†p, while not computing the matrix W†explicitly (see also [29]). For all experiments in Section 5, the time required to compute the central reconstruction varied from a few seconds up to two minutes on a standard PC, depending on the number of projections and image size. 5 4. Quality bounds for binary solutions In all the results in the following subsections, we consider a fixed system Wx =pcorresponding to a binary tomography problem, and refer to the central reconstruction of this system as x∗. As a substantial number of bounds will be given throughout this paper, we introduce the following notation that will be further defined in the remainder of the paper: •The expressions a(i) (i=1,2,3,4) will represent bounds on the number of pixel differences between any two binary solutions of the reconstruction problem. •The expressions b(i) (i=1,2,3,4) will represent bounds on the number of pixel difference between a certain given binary image (not necessarily a solution) and any binary solution of the reconstruction problem. •The expressions c(i) (i=1,2) will represent bounds on the number of pixel differences between the rounded central reconstruction ¯ rand any binary solution. The bounds within each class a,b, and crepresent upper bounds for the same distance measure and can therefore be compared. 4.1. Elementary bounds based on the central reconstruction In this subsection, a first set of bounds are derived. They follow from the fact that the Euclidean distance between the central reconstruction and any binary solution of the reconstruction problem can be determined from the projections. We start by noticing that the Euclidean norm of any binary solution of the tomography problem is determined by the projection data: Lemma 1. Let ¯ x∈¯ SW(p). Then, k¯ xk2 2=kpk1 k. Proof. By the definition of the `1-norm, kpk1=Pm i=1|pi|=Pm i=1pi, since pi≥0 (i=1,...,n). Also, m X i=1 pi= m X i=1        n X j=1 wi j ¯xj        = n X j=1       m X i=1 wi j      ¯xj= n X j=1 k¯xj,(2) and therefore kpk1=kPn j=1¯xj. As ¯ x∈ {0,1}n, we have k¯ xk2 2=k¯ xk1=Pn j=1¯xj=kpk1 k. The following lemma illustrates the importance of the central reconstruction, the shortest real-valued solution in SW(p), by showing that the binary solutions are the shortest among all integer solutions of the system. Lemma 2. Let ¯ x∈¯ SW(p)and y∈SW(p)∩Zn. Then k¯ xk2≤ kyk2, with equality if and only if y∈¯ SW(p). Proof. Note that the statement is proved in [14], see Problem 2 and the subsequent paragraph. However, for the convenience of the reader we give the proof here. We have || ¯ x||2 2= n X i=1 ¯x2 i= n X i=1 ¯xi= n X i=1 yi=kpk1 k.(3) 6 Observing that n X i=1 yi≤ n X i=1 y2 i=||y||2 2,(4) with equality if and only if yis binary, yields the result. Lemma 3. Let ¯ x∈¯ SW(p). Then k¯ x−x∗k2=qkpk1 k− kx∗k2 2. Proof. From the definition of x∗we have (¯ x−x∗)∈ N(W), and x∗⊥(¯ x−x∗). Applying Pythagoras’ Theorem and Lemma 1 yields k¯ x−x∗k2 2=kpk1 k− kx∗k2 2.(5) Define R=qkpk1 k− kx∗k2 2. We will use this constant throughout the remainder of this article, and refer to Ras the central radius. According to Lemma 3, any binary solution of the reconstruction problem is on the hypersphere centered in x∗with radius R. Supposing the existence of at least two different binary solutions, Lemma 3 allows us to derive an upper bound for the number of pixel differences between those solutions. Theorem 4. Let ¯ x,¯ y∈¯ SW(p)and put a(1) =4R2. Then d(¯ x,¯ y)≤a(1). Proof. According to Lemma 3, we have k¯ x−x∗k2=k¯ y−x∗k2=R. Therefore, k¯ x−¯ yk2≤ k ¯ x−x∗k2+k¯ y−x∗k2=2R. As ¯ xand ¯ yare binary, we have d(¯ x,¯ y)=k¯ x−¯ yk1=k¯ x−¯ yk2 2. Using the triangle inequality, a simple bound can also be given for the distance between any binary image and a solution of the reconstruction problem, as follows: Corollary 5. Let ¯ v∈ {0,1}nbe a given binary image and put b(1) =(R+k¯ v−x∗k2)2. Then for any binary image ¯ x∈¯ SW(p), we have d(¯ x,¯ v)≤b(1). 4.2. Bounds based on rounding the central reconstruction The fact that all elements of ¯ SW(p) have equal distance to the central reconstruction x∗, combined with the facts that binary solutions are the shortest solutions among all integer solutions (Lemma 2) and that x∗is the shortest real-valued solution, suggests that binary solutions can often be found near x∗. It is therefore natural to consider the image that is obtained by rounding each entry of x∗to the nearest binary value. In this section, we will derive several bounds based on the number of differences between a binary solution of the reconstruction problem and a binary image obtained by rounding x∗. For α∈R, let bin(α)=min(|α|,|1−α|). Put T=qPn i=1bin2(x∗ i), i.e., the Euclidean distance from x∗to the nearest binary vector. We will use this constant throughout this paper and refer to Tas the central rounding distance. 7 Corollary 6. If R <T, then ¯ SW(p)=∅. If R =T, then all solutions in ¯ SW(p)can be obtained by rounding the values in x∗to the nearest binary values, and variations are only possible for the entries i where x∗ i=1 2. Let T={¯ v∈ {0,1}n:k¯ v−x∗k2=T}and let ¯ r∈ T , i.e., ¯ ris among the binary vectors that are nearest to x∗in the Euclidean sense. If R>Tand R−Tis small, it is possible to say that a fraction of the rounded values are correct, i.e., to provide an upper bound on the number of pixel differences between any solution in ¯ SW(p) and ¯ r. In most cases we can not say which rounded values are correct. Lemma 7. Let ¯ r∈ T and let ¯ v∈ {0,1}nbe any binary vector. Then k¯ v−x∗k2 2=T2+Pi∈D(¯ v,¯ r)|2x∗ i−1|. Proof. We have the following identities: k¯ v−x∗k2 2=k¯ v−¯ r+¯ r−x∗k2 2 =k¯ r−x∗k2 2+2h¯ r−x∗,¯ v−¯ ri+h¯ v−¯ r,¯ v−¯ ri =T2+2h¯ r−¯ v,x∗i+k¯ vk2 2− k¯ rk2 2 =T2+2 n X i=1 (¯ri−¯vi)x∗ i+ n X i=1 ¯vi− n X i=1 ¯ri =T2+ n X i=1 (¯ri−¯vi)(2x∗ i−1) =T2+X i∈D(¯ v,¯ r) |2x∗ i−1|. Lemma 7 can be interpreted as follows: consider the set of entries where ¯ rand ¯ vare different. If we transform ¯ rinto ¯ vby performing a sequence of single-entry changes (either from 0 to 1, or from 1 to 0), each time an entry iof ¯ ris changed the squared Euclidean distance from the current vector to x∗increases by si=|2x∗ i−1|. Let πbe a permutation of {1,...,n}such that sπ(1) ≤sπ(2) ≤. . . ≤sπ(n), which can be obtained by sorting the entries siin increasing order. Corollary 8. Let ¯ r∈ T and let ¯ v∈ {0,1}nbe any binary vector. Then k¯ v−x∗k2 2≥T2+P` i=1sπ(i), where `=d(¯ r,¯ v). Proof. According to Lemma 7 we have k¯ v−x∗k2 2=T2+X i∈D(¯ r,¯ v) si≥T2+ ` X i=1 sπ(i). As the Euclidean distance from x∗to any ¯ x∈¯ SW(p) is R, a bound can now be derived on the maximal number of pixels in ¯ rthat must be changed to move from ¯ rto ¯ x. 8 Theorem 9. Let ¯ r∈ T ,¯ x∈¯ SW(p). Put c(1) =max        0≤`≤n: ` X i=1 sπ(i)≤R2−T2       . Then d(¯ x,¯ r)≤c(1). Proof. As ¯ x∈¯ SW(p), we have k¯ x−x∗k2 2=R2.Applying Lemma 7, we find that R2−T2=X i∈D(¯ x,¯ r) si≥ d(¯ x,¯ r) X i=1 sπ(i), which implies that d(¯ x,¯ r)≤c(1). The proof of Theorem 9 can be interpreted as follows: consider the set of entries where ¯ r and ¯ xare different. If we transform ¯ rinto ¯ xby performing a sequence of single-entry changes (either from 0 to 1, or from 1 to 0), each time an entry iof ¯ ris changed the squared Euclidean distance from the current vector to x∗increases by si=|2x∗ i−1|. As all binary solutions of the reconstruction problem are on a hypersphere centered in x∗with radius R, we know that once we have crossed the boundary of this hypersphere, a binary solution can no longer be obtained by changing the values of additional entries that have not yet been changed. An upper bound on the number of differences between ¯ rand ¯ xcan be obtained by counting the number of steps required to cross the hypersphere, each time choosing a pixel which results in the minimal increase of the distance to x∗. The following two Corollaries follow directly from Theorem 9: Corollary 10. Let ¯ r∈ T ,¯ x,¯ y∈¯ SW(p)and let a(2) =2c(1) with c(1) defined as in Theorem 9. Then d(¯ x,¯ y)≤a(2). Corollary 11. Let ¯ r∈ T and let ¯ v∈ {0,1}nbe a given binary image and let b(2) =c(1) +d(¯ r,¯ v) with c(1) defined as in Theorem 9. Then d(¯ x,¯ v)≤b(2). In fact, the bound from Corollary 10 can be sharpened by noting that we can assume that the sets D(¯ r,¯ x) and D(¯ r,¯ y) are disjoint, as entries that occur in both sets do not contribute to the number of differences between ¯ xand ¯ y. This observation leads to the following Theorem: Theorem 12. Let ¯ r∈ T ,¯ x,¯ y∈¯ SW(p). Put a(3) =max        0≤`≤n: ` X i=1 sπ(i)≤2(R2−T2)       . Then d(¯ x,¯ y)≤a(3). Proof. Define ˆ xby ˆxi=¯riif ¯xi=¯yi, and ˆxi=¯xiotherwise. Define ˆ yanalogously. Then d(ˆ x,ˆ y)=d(¯ x,¯ y), kˆ x−x∗k2 2≤R2, and kˆ y−y∗k2 2≤R2. Hence, 2R2≥ k ˆ x−x∗k2 2+kˆ y−y∗k2 2=2T2+ Pi∈D(¯ r,ˆ x)si+Pi∈D(¯ r,ˆ y)si.As D(¯ r,ˆ x) and D(¯ r,ˆ y) are disjoint, we have 2R2−2T2≥Pd(¯ r,ˆ x)+d(¯ r,ˆ y) i=1sπ(i). This implies that d(¯ x,¯ y)=d(ˆ x,ˆ y)≤d(¯ r,ˆ x)+d(¯ r,ˆ y)≤a(3). A similar bound can be derived for the case where a particular binary image ¯ v, not necessarily a solution of the tomography problem, is given. For this, we transform ¯ vinto ¯ rand then we perform a sequence of single-entry changes in ¯ rwith the exclusion of the pixels that differ between ¯ vand ¯ rbecause they have already been counted as wrong pixels of ¯ v. 9 wrong with respect to the unknown phantom image. 6. Discussion of the results Despite the facts that the four phantoms have strong differences in shape and morphology, and that the grid and strip models are quite different, the results shown in Figs. 6–9 are surprisingly consistent throughout all experiments. Most of the bounds become smaller as the number of projection directions is increased but monotonicity is not a property of all the bounds presented in Section 4. From the difference between the bounds presented in Section 4.1 and the bounds based on the rounded central reconstruction, we see that in most cases the phantom ¯ xis substantially closer to ¯ rthan to x∗. In Figs. 6 and 8, it can be observed that the true fraction of pixel differences between the phantom image ¯ xand the rounded central reconstruction ¯ r, denoted by Ec, is often approximated quite well by the bound c, in particular for the grid model. This indicates that with respect to ¯ r, the bounds presented in this article can be quite sharp. In Fig. 8(a), parts of the graphs for the bounds aand c, for more than 6 projections, are missing. In fact, in this case all of them are zero, such that they cannot be displayed in the logarithmic scale. This illustrates that our theorems for bounding the distance between any two binary solutions can be used to prove uniqueness of a binary solution, even when the corresponding real-valued system of equations is underdetermined. In most of the experiments, the graphs for a(2) and a(4) are almost indistinguishable and the bound a(3) is smaller than a(4) but it can happen, as in Fig. 9(a), that this is not the case. In Fig. 9(a) the bound a(4) proves uniqueness from 6 projection directions or more while the bound a(3) does not for 12, 14 and 16 projections. Although a(3) is the sharpest bound in most of the experiments, this is not true in general and all bounds should be used to obtain the best result. The bounds a(2), a(3) and a(4) use similar principles and their computation can be performed simultaneously. 7. Outlook and conclusions In this article, we have presented a range of general bounds on the accuracy of reconstructions in binary tomography, with respect to the unknown original object. The bounds can be computed within reasonable time and give guarantees on the number of pixels that can be different between any two binary solutions of the tomography problem, on the difference between an image obtained by rounding the central reconstruction and any binary solution, and on the difference between any binary image and any binary solution. The experimental results show that by using these bounds, one can prove that the number of differences between binary reconstructions must be very small, even when the corresponding real-valued system of equations is severely underdetermined. In order to make these bound practically useful, our results will have to be extended to deal with noisy projection data, which we will incorporate in future research. Acknowledgements L.H. was supported by the OTKA grants K67580 and K75566, and by the T ´ AMOP 4.2.1./B09/1/KONV-2010-0007 project. The project is implemented through the New Hungary Devel16 opment Plan, cofinanced by the European Social Fund and the European Regional Development Fund. W.F. acknowledges support from the Erasmus Mundus program of the European Union. K.J.B. was supported by the Netherlands Organisation for Scientific Research (NWO), programme 639.072.005. [1] Alpers, A., 2003. Instability and Stability in Discrete Tomography. Ph.D. thesis. Technische Universit¨ at M¨ unchen. Shaker Verlag, ISBN 3-8322-2355-X. [2] Alpers, A., Brunetti, S., 2007. Stability results for the reconstruction of binary pictures from two projections. Image and Vision Computing 25, 1599–1608. [3] Alpers, A., Gritzmann, P., 2006. On stability, error correction, and noise compensation in discrete tomography. SIAM Journal on Discrete Mathematics 20, 227–239. [4] Batenburg, K.J., 2003. Analysis and optimization of an algorithm for discrete tomography. Electr. Notes Discrete Math. 12, 35–46. [5] Batenburg, K.J., 2008. A network flow algorithm for reconstructing binary images from continuous X-rays. J. Math. Im. Vision 30, 231–248. [6] Batenburg, K.J., Bals, S., Sijbers, J., K¨ ubel, C., Midgley, P.A., Hernandez, J.C., Kaiser, U., Encina, E.R., Coronado, E.A., Van Tendeloo, G., 2009. 3D imaging of nanomaterials by discrete tomography. Ultramicroscopy 109, 730– 740. [7] Batenburg, K.J., Fortes, W., Hajdu, L., Tijdeman, R., 2011. Bounds on the difference between reconstructions in binary tomography, in: Proceedings of the 16th IAPR international conference on Discrete geometry for computer imagery, Springer-Verlag, Berlin, Heidelberg. pp. 369–380. [8] Batenburg, K.J., Sijbers, J., 2011. DART: a practical reconstruction algorithm for discrete tomography. IEEE Trans. Image Processing 20, 2542–2553. [9] Ben-israel, A., Greville, T.N.E., 2002. Generalized inverses: Theory and applications. Canadian Math. Soc. [10] Bj¨ orck, Å., 1996. Numerical methods for least square problems. SIAM, Link¨ oping University, Sweden. [11] Brunetti, S., Daurat, A., 2005. Stability in Discrete Tomography: some positive results. Discrete Appl. Math. 147, 207–226. [12] Chen, G.H., Tang, J., Leng, S., 2008. Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets. Medical Physics 35, 660–663. [13] Donoho, D., 2006. Compressed sensing. IEEE Trans. Information Theory 52, 1289–1306. [14] Hajdu, L., Tijdeman, R., 2001. Algebraic aspects of discrete tomography. J. Reine Angew. Math. 534, 119–128. [15] Herman, G.T., 2009. Fundamentals of Computerized Tomography: Image reconstruction from projections. Springer. [16] Herman, G.T., Kuba, A. (Eds.), 1999. Discrete Tomography: Foundations, Algorithms and Applications. Birkh¨ auser, Boston. [17] Herman, G.T., Kuba, A. (Eds.), 2007. Advances in Discrete Tomography and its Applications. Birkh¨ auser, Boston. [18] Jinschek, J.R., Batenburg, K.J., Calderon, H.A., Kilaas, R., Radmilovic, V., Kisielowski, C., 2007. 3-D reconstruction of the atomic positions in a simulated gold nanocrystal based on discrete tomography. Ultramicroscopy 108(6), 589–604. [19] Kak, A.C., Slaney, M., 2001. Principles of Computerized Tomographic Imaging. SIAM. [20] Midgley, P.A., Dunin-Borkowski, R.E., 2009. Electron tomography and holography in materials science. Nature Materials 8, 271–280. [21] Saad, Y., 2003. Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia, PA, USA. [22] Sch¨ ule, T., Schn¨ orr, C., Weber, S., Hornegger, J., 2005. Discrete tomography by convex-concave regularization and D.C. programming. Discr. Appl. Math 151, 229–243. [23] Sidky, E.Y., Anastasio, M.A., Pan, X., 2010. Image reconstruction exploiting object sparsity in boundary-enhanced X-ray phase-contrast tomography. Optics Express 18, 10404–10422. [24] Sidky, E.Y., Kao, C.M., Pan, X., 2006. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT. J. X-ray Sci. Tech. 14, 119–139. [25] Van Aert, S., Batenburg, K.J., Rossell, M.D., Erni, R., Van Tendeloo, G., 2011. Three-dimensional atomic imaging of crystalline nanoparticles. Nature 470, 374–377. [26] Van Dalen, B., 2009a. On the difference between solutions of discrete tomography problems. Journal of Combinatorics and Number Theory 1, 15–29. [27] Van Dalen, B., 2009b. On the difference between solutions of discrete tomography problems II. Pure Mathematics and Applications 20, 103–112. [28] Van Dalen, B., 2009c. Stability results for uniquely determined sets from two directions in discrete tomography. Discrete Mathematics 309, 3905–3916. [29] Van der Sluis, A., Van der Vorst, H.A., 1990. SIRT and CG-type methods for the iterative solution of sparse linear least-squares problems. Linear Algebra Appl. 130, 257–302. 17 [30] Zhua, J., Li, X., Ye, Y., Wang, G., 2008. Analysis on the strip-based projection model for discrete tomography. Discrete Appl. Math. 156, 2359–2367. 1234567 10−3 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (a) Phantom 1, 32×32 0 5 10 15 20 25 30 10−4 10−3 10−2 10−1 Number of angles Fraction of pixels a c Ec (b) Phantom 1, 128×128 0 50 100 150 200 10−6 10−5 10−4 10−3 10−2 10−1 Number of angles Fraction of pixels a c Ec (c) Phantom 1, 512×512 0 2 4 6 8 10 12 10−3 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (d) Phantom 2, 32×32 0 10 20 30 40 50 60 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (e) Phantom 2, 128×128 0 50 100 150 200 10−5 10−4 10−3 10−2 10−1 Number of angles Fraction of pixels a c Ec (f) Phantom 2, 512×512 0 5 10 15 10−3 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (g) Phantom 3, 32×32 0 10 20 30 40 50 60 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (h) Phantom 3, 128×128 0 50 100 150 200 10−6 10−4 10−2 100 Number of angles Fraction of pixels a c Ec (i) Phantom 3, 512×512 0 5 10 15 10−3 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (j) Phantom 4, 32×32 0 20 40 60 80 10−5 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (k) Phantom 4, 128×128 0 50 100 150 200 10−6 10−4 10−2 100 Number of angles Fraction of pixels a c Ec (l) Phantom 4, 512×512 Figure 6: Grid model: computed bounds as a function of the number of projection directions. 18 0 5 10 15 20 10−3 10−2 10−1 100 101 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (a) Phantom 1, 32×32 0 20 40 60 80 10−5 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (b) Phantom 1, 128×128 0 50 100 150 200 10−5 10−4 10−3 10−2 10−1 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (c) Phantom 1, 512×512 0 5 10 15 20 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (d) Phantom 2, 32×32 0 20 40 60 80 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (e) Phantom 2, 128×128 0 50 100 150 200 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (f) Phantom 2, 512×512 0 5 10 15 20 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (g) Phantom 3, 32×32 0 20 40 60 80 10−5 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (h) Phantom 3, 128×128 0 50 100 150 200 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (i) Phantom 3, 512×512 0 5 10 15 20 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (j) Phantom 4, 32×32 0 20 40 60 80 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (k) Phantom 4, 128×128 0 50 100 150 200 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (l) Phantom 4, 512×512 Figure 7: Grid model: computed bounds as a function of the number of projection directions. Note that the graphs for a(2) and a(4) have strong overlap and can hardly be distinguished. 19 0 2 4 6 8 10−3 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (a) Phantom 1, 32×32 0 20 40 60 80 10−5 10−4 10−3 10−2 10−1 Number of angles Fraction of pixels a c Ec (b) Phantom 1, 128×128 0 50 100 150 200 10−6 10−5 10−4 10−3 10−2 10−1 Number of angles Fraction of pixels a c Ec (c) Phantom 1, 512×512 0 5 10 15 20 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (d) Phantom 4, 32×32 0 20 40 60 80 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (e) Phantom 4, 128×128 0 50 100 150 200 10−3 10−2 10−1 100 Number of angles Fraction of pixels a c Ec (f) Phantom 4, 512×512 Figure 8: Strip model: computed bounds as a function of the number of projection directions for Phantoms 1 and 4. 0 5 10 15 20 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (a) Phantom 1, 32×32 0 20 40 60 80 10−4 10−3 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (b) Phantom 1, 128×128 0 50 100 150 200 10−4 10−3 10−2 10−1 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (c) Phantom 1, 512×512 0 5 10 15 20 10−0.7 10−0.5 10−0.3 10−0.1 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (d) Phantom 4, 32×32 0 20 40 60 80 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (e) Phantom 4, 128×128 0 50 100 150 200 10−2 10−1 100 Number of angles Fraction of pixels a(1) a(2) a(3) a(4) (f) Phantom 4, 512×512 Figure 9: Strip model: computed bounds as a function of the number of projection directions for Phantoms 1 and 4. Note that the graphs for a(2) and a(4) have strong overlap and can hardly be distinguished. 20