Discrete tomography is concerned with the recovery of discrete images (i.e., images whose pixels take on a small number of prescribed grey values) from a few of their projections (i.e., sums of the pixel values along various directions). J. Steiner first introduced the problem of discrete tomography in 1838 . Early work on the subject mostly deals with the mathematical analysis, combinatorics, and geometry. Since the 1970s, the development of algorithms for discrete tomography has become an active area of research as well. It has found applications in image processing, atomic-resolution electron microscopy, medicine and material science (for example,[2, 3, 4, 5, 6, 7, 8]).
I-a Mathematical formulation
The discrete tomography problem may be mathematically formulated as follows. We represent an image by a grid of pixels taking values . The projections are linear combinations of the pixels along
different (lattice) directions. We denote the linear transformation from image to projection data by
where denotes the value of the image in the cell, is the (weighted) sum of the image along the ray and gives the length of the ray in the cell.
The goal is to find a solution to this system of equations with the constraint that , i.e.,
When the system of equations does not have a unique solution, finding one that only takes values in has been shown to be to be an NP-hard problem when .
Due to the presence of noise, the system of equations may not have a solution, and the problem is sometimes formulated as a constrained least-squares problem
Obviously, the constraints are non-convex and solving (1) exactly is not trivial. Next, we briefly discuss some existing approaches for solving it.
I-B Literature Review
Methods for solving (1) can be roughly divided into four classes: algebraic methods, stochastic sampling methods, (convex) relaxation and (heuristic/greedy) combinatorial approaches.
The algebraic methods exploit the algebraic structure of the problem and may give some insight into the (non-) uniqueness and the required number of projections [10, 11]. While theoretically very elegant, these methods are not readily generalised to realistic projection models and noisy data.
The stochastic sampling methods12, 13, 14]. These methods are very flexible but may require a prohibitive number of samples when applied to large-scale datasets.
Relaxation methods are based on some form convex or non-convex relaxation of the constraint. This allows for a natural extension of existing variational formulations and iterative algorithms [15, 16, 17, 18]. While these approaches can often be implemented efficiently and apply to large-scale problems, getting them to converge to the correct binary solution can be challenging.
algorithms, finally, combine ideas from combinatorial optimization and iterative methods. Such methods are often efficient and known to perform well in practice[19, 20].
A more extensive overview of various methods for binary tomography and variants thereof (e.g., with more than two grey levels) are discussed in .
I-C Contributions and outline
We propose a novel, convex, reformulation for discrete tomography with two grey values (often referred to as binary tomography). Starting from the constrained least-squares problem (1) we derive a corresponding Lagrange dual problem, which is convex by construction. Solving this problem yields an image with pixel values in . Setting the remaining zero-valued pixels to or generates a feasible solution of (1) but not necessarily an optimal one. In this sense, our approach is a relaxation. Exhaustive enumeration of small-scale () images with few directions () show that if the problem has a unique solution, then solving the dual problem yields the correct solution. When there are multiple solutions, the dual approach finds the common elements of the solutions, leaving the remaining pixels undefined (zero-valued). We conjecture that this holds for larger and as well. This implies that we can only expect to usefully solve problem instances that allow a unique solution and characterize the non-uniqueness when there are a few solutions. For practical applications, the most relevant setting is where the equations alone do not permit a unique solution, but the constrained problem does. Otherwise, more measurements or prior information about the object would be needed in order to usefully image it. With well-chosen numerical experiments on synthetic and real data, we show that our new approach is competitive for practical applications in X-ray tomography as well.
The outline of the paper is as follows. We first give an intuitive derivation of the dual problem for invertible before presenting the main results for general . We then discuss two methods for solving the resulting convex problem. Then, we offer the numerical results on small-scale binary problems to support our conjecture. Numerical results on numerical phantoms and real data are presented in Section IV. Finally, we conclude the paper in Section V.
Ii Dual Problem
For the purpose of the derivation, we assume that the problem has been rescaled such that the desired pixel values are . The least-squares binary tomography problem can then be formulated as:
where is an auxiliary variable and where denotes the elementwise signum function. In our analysis, we consider the strict signum function, i.e., is not defined. The Lagrangian for this problem is given by :
The variable is the Lagrange multiplier associated with the equality constraint . We refer to this variable as the dual variable in the remainder of the paper. We define the dual function corresponding to the Lagrangian (3) as:
The primal problem (2) can be expressed in terms of this dual objective as
As the dual function is always concave , this provides a way to define a convex reformulation of the original problem. We should note two important aspects of duality theory here: i) we are not guaranteed in general that maximizing yields a solution to the primal problem; ii) the reformulation is only computationally useful if we can efficiently evaluate . The conditions under which the dual problem yields a solution to the primal problem are known as Slater’s conditions  and are difficult to check in general unless the primal problem is convex. We will later show, by example, that the dual problem does not always
solve the primal problem. Classifying under which conditions wecan solve the primary problem via its dual is beyond the scope of this paper.
It turns out we can obtain a closed-form expression for . Before presenting the general form of , we first present a detailed derivation for invertible to provide some insight.
The Lagrangian is separable in terms of and . Hence, we can represent the dual function as the sum of two functions, and .
First, we consider . Setting the gradient to zero we find the unique minimizer:
Substituting back in the expression and re-arranging some terms we arrive at the following expression for :
where is a weighted -norm.
Next, we consider . Note that this function is separable in terms of
The function achieves its smallest value for when . This solution is not unique of course, but that does not matter as we are only interested in the sign of . When the function takes on value regardless of the value of . We thus find
Hence, the dual function for the Lagrangian in (3) takes the following explicit form:
and leads to the following dual problem to (2).
This optimization problem is famously known as the least absolute shrinkage and selection operator (LASSO) 
in the statistics literature. It tries to find a sparse vector in image space by minimizing the distance to the back projection of the data in a scalednorm. The primal solution can be synthesized from the solution of the dual problem via .
It is important to note at this point that the solution of the dual problem only determines those elements of the primal problem, , for which
. The remaining degrees of freedom inneed to be determined by alternative means. The resulting solution is at least a feasible solution of the primal problem, but not necessarily the optimal one.
To gain some insight into the behaviour of the dual objective, consider a one-dimensional example with :
The solution to this problem is given by . The corresponding dual problem is
the solution of which is given by . Hence, for , the solution of the dual problem yields the desired solution. For , however, the dual problem yields in which case the primal solution is not well-defined. We will see in a later section that the when using certain iterative methods to solve the dual problem, the iterations will naturally approach the solution from the correct side, so that the sign of the approximate solution may still be useful.
Ii-B Main results
We state the main results below. The proofs for these statements are provided in the Appendix section.
The dual objective of (2) for general is given by
where denotes the pseudo-inverse and is the row-space of (i.e., the range of ). This leads to the following dual problem
The dual problem (12) can be restated as
and the primal solution is recovered through .
For and full rank, we have and the formulation (13) simplifies to
This form implicitly handles the constraints on the search space of in Theorem 1. It allows us to use the functional form for matrix thereby reducing the storage and increasing the computational speed to find an optimal dual variable .
The dual problem for a binary tomography problem with grey levels is given by:
where is an asymmetric one-norm. The primal solution is obtained using
where denotes the Heaviside function.
We summarize the procedure in algorithm 1 for finding the optimal solution via solving the dual problem.
Ii-C Solving the dual problem
where and the soft thresholding operator is applied component-wise to its input. We can interpret this algorithm as minimizing subsequent approximations of the problem, as illustrated in figure 1.
An interesting note is that, when starting from , the first iteration yields a thresholded version of . As such, the proposed formulation is a natural extension of a naive segmentation approach and allows for segmentation in a data-consistent manner.
If is invertible we have and it seems more natural to solve (14) instead. Due to the appearance of in the one-norm, it is no longer straightforward to apply a proximal gradient method. A possible strategy is to replace the one-norm with a smooth approximation of it, such as . As illustrated in figure 2, this will slightly shift the minimum of the problem. Since we are ultimately only using the sign of the solution, this may not be a problem. The resulting objective is smooth and can be solved using any gradient-descent algorithm.
We also note that splitting methods can be used to solve (14). For example, the alternating direction method of multipliers (ADMM)  and/or split-Bregman method . Another class of method that can solve (14) are the primal-dual methods (e.g., Arrow-Hurwicz primal-dual algorithm , Chambolle-Pock algorithm ). These methods rely on the proximal operators of functions and iterate towards finding the saddle point of the problem. If the proximal operators are simple, these are computationally faster than the splitting methods.
The dual problem (15) for binary tomography problem with grey levels is also solved using proximal gradient method. We provide the proximal operator for an asymmetric one-norm in the following theorem.
The proximal operator for an asymmetric one-norm function
with , is given by
where , and is an asymmetric soft-thresholding function
See the Appendix. ∎
Iii Numerical experiments - Binary tomography
To illustrate the behaviour of the dual approach, we consider the simple setting of reconstructing an image from its sums along lattice directions (here restricted to the horizontal, vertical and two diagonal directions, so ). For the problem is known to be NP-hard. For small we can simply enumerate all possible images, find all solutions in each case by a brute-force search and compare these to the solution obtained by the dual approach. To this end, we solved the resulting dual problem (13) using the CVX package in Matlab . This yields an approximate solution, and we set elements in the numerically computed dual solution smaller than to zero. We then compare the obtained primal solution, which has values to the solution(s) of the binary tomography problem. From performing these computations for and we conclude the following:
If the problem has a unique solution then the dual approach retrieves it.
If the problem has multiple solutions then the dual approach retrieves the intersection of all solutions. The remaining pixels in the dual solution are undetermined (have value zero).
An example is shown in 3.
A summary of these results is presented in table I. The table shows the number of cases with a unique solution where the dual approach gave the correct solution and, in case of multiple solutions, the number of cases where the dual approach correctly determined the intersection of all solutions). In a few instances with multiple solutions, CVX failed to provide an accurate solution (denoted with in the table).
Based on these experiments, we conjecture that there is a subclass of the described binary tomography problem that is not NP-hard. We should note that, as grows the number of cases that have a unique solution grows smaller unless grows accordingly. It has been established that binary images that are h,v,d-convex111For h,v,d-convexity one uses the usual definition of convexity of a set but one considers only line segments in the horizontal, vertical and diagonal directions. can be reconstructed from their horizontal, vertical and diagonal projections in polynomial time [32, 19]. However, we can construct images that are not h,v,d-convex but still permit a unique solution, see figure 4. Such images are also retrieved using our dual approach.
Iv Numerical experiments - X-ray tomography
In this section, we present numerical results for limited-angle X-ray tomography on a few numerical phantoms and an experimental X-ray dataset. First, we describe the phantoms and the performance measures used to compare our proposed dual approach to several state-of-the-art iterative reconstruction techniques. We conclude this section with results on an experimental dataset. All experiments are performed using Matlab in conjunction with the ASTRA toolbox .
For the synthetic tests, we consider four phantoms shown in figure 5. All the phantoms are binary images of size 128 128 pixels. The detector has 128 pixels, and the distance between the adjacent detectors is the same as the pixel size of the phantoms. We consider a parallel beam geometry for the acquisition of the tomographic data in all the simulation experiments.
We perform three different tests to check the robustness of the proposed method. First, we consider the problem of sparse projection data. For all the phantoms, we first start with 45 projections at angles ranging from to and subsequently reduce the number of angles. This setup is also known as sparse sampling, where the aim is to reduce the scan time by decreasing the number of angles.
Next, we consider a limited angle scenario. Such situation usually arises in practice due to the limitations of the setup. For the test, we acquire projections in the range for . Reconstruction of limited angle data is known to lead to so-called streak artifacts in the reconstructed image. Strategies to mitigate these streak artifacts include the use of regularization method with some prior information. We will experiment how the discrete tomography can lead to the removal of these artifacts.
Finally, we test the performance of the proposed method in the presence of noise. We consider an additive Gaussian noise in these experiments. We measure the performance of our approach for tomographic data with an SNR of dB.
Iv-C Comparison with other reconstruction methods
There exist a vast amount of reconstruction methods for tomography. Here, we consider the following three:
: The total-variation method leads to an optimization problem described below:
where is a corresponding regularization parameter, and matrix captures the discrete gradient in both directions. We use the Chambolle-Pock method  to solve the above optimization problem with non-negativity constraints on the pixel values. For each case, we perform 200 iterations. To avoid slow convergence, we scale the matrices and to have unit matrix norm. The regularization parameter is selected using Morozov’s discrepancy principle using the correct noise level. We segment the TV-reconstructed image through Otsu’s thresholding algorithm.
: We use the method described in  for DART on the above binary images. The grayvalues are taken to be the same as true grayvalues. We perform 20 ARM iterations initially before performing 40 DART iterations. In each DART iteration, we do 3 ARM iterations. We use the segmented image as a result of the DART iterations to perform the further analysis of the method.
Iv-D Performance Measures
In order to evaluate the performance of reconstruction methods, we use the following two criteria
The root-mean-square error
measures how well the forward-projected reconstructed image matches the projection data. This measure is useful in practice, as it does not require knowledge of the ground truth. If the RMS value is close to the noise level of the data, the reconstruction is considered as a good reconstruction.
The Jaccard index
measures the similarity between the reconstructed image and the ground truth in a discrete sense. If JI has high value (close to 1), the reconstruction is considered good. Although this measure cannot generally be applied on real datasets, it is a handy measure to compare the various reconstruction methods on synthetic examples.
Iv-E Experimental data setup
We use the experimental X-ray projection data of a carved cheese slice provided by the FIPS group . Figure 6 shows a high-resolution filtered back-projection reconstruction of the data. The cheese contains the letters C and T and the object is (approximately) binary with two grey levels corresponding to calcium-containing organic compounds of cheese and air. The dataset consists of projection data with three different resolutions (, , ) and the corresponding projection matrix modelling the linear operation of X-ray transform. We perform two sets of experiments: (1) Sparse sampling with 15 angles ranging from to and limited-angle using 15 projections from to .
Iv-F Sparse projections test
Figure 7 presents the reconstruction results from various methods for phantoms 1 and 2. The tomographic data are generated for ten equidistant projection angles from to . The reconstruction results show the difference between the reconstructed image and the ground truth. It is evident that, compared to the other methods, the proposed method reconstructions are very close to the ground truth. The results from LSQR are the worst as it does not incorporate any prior information about the model. The TV method also leads to artifacts as it includes the partial information about the model. The DART and DP are very close to each other. For all the phantoms, we tabulate the data misfit (RMS) and Jaccard index (JI) in Table II.
We also show reconstructions with the proposed approach for a varying number of projection angles in Figure 8. We note that the problem becomes harder to solve as the number of projection angles gets smaller. Hence, we may also expect the reconstruction to become poorer. We see that the proposed approach can reconstruct almost correctly with as few as ten projection angles.
Iv-G Limited angle test
Figure 9 shows the results of phantom 2 with various reconstruction methods for limited angle tomography ( 10 equispaced angles in the range to ). It is visible that the reconstructions from DART and the proposed method are very close to true images of the phantoms. The values for data misfit and Jaccard index for all the tests with each of the synthetic phantoms are tabulated in Table III.
We also look at how the reconstructions with the proposed method varies with limiting the angle (see Figure 10). As the angle gets limited, the reconstruction problem gets difficult. The proposed method can reconstruct almost perfectly with angle limited to .
Iv-H Noisy projection test
This test aims to check the robustness of the proposed method to noise in the data. We perform four experiments with varying levels of additive Gaussian noise in the data. Figure 11 shows the results on phantoms 1 and 2 for increasing noise level. We see that the reconstruction is stable against a moderate amount of noise and degrades gradually as the noise level increases.
Iv-I Real data test
We look at the results of reconstructions from the proposed method for two sets of experiments at various resolutions and compare them with the reconstructions from LSQR and TV. Since the ground truth image is not available, we compare these reconstructions visually.
In order to apply DART and DP, we first need to estimate the grey values of the object. The object, a thin slice of cheese, consists of two materials; the organic compound of the cheese, which is we assume to be homogeneous, and air. For air, the grey value is zero. We estimate the grey value of the organic compound of cheese from the histogram of an FBP reconstruction provided with the data. Figure12 represents the histogram. We obtain a value of 0.00696 for this compound.
We first consider the reconstructions from sparse angular sampling. We have a tomographic data from 15 projections spanning from to . The tests are performed on three different resolutions: , and . Figure 13 presents the results of the reconstructions with LSQR, TV and DP for these resolutions. The DP reconstruction is discrete and correctly identifies the letters C and T with also a little hole at the left side of C. Although LSQR reconstruction is poor for , it improves with the resolution. We still see the mild streak artifacts in these reconstructions. The TV reconstruction removes these streak artifacts but fails to identify the homogeneous cheese slice correctly.
In the second test, we limit the projection angles to . Figure 14 shows the results of the reconstructions from LSQR, TV, and DP for three different resolutions. We see that the reconstructions improve with increment in the resolution. LSQR reconstructions have severe streak artifacts, which are the characteristics of the limited data tomography. TV and DP reconstructions do not possess these artifacts. TV reconstruction can capture the shape of the cheese, but it blurs out the carved parts C and T. DP reconstructs the shape of cheese quite accurately and has C and T are also identified.
We presented a novel convex formulation for binary tomography. The problem is primarily a generalized LASSO problem that can be solved efficiently when the system matrix has full row rank or full column rank. Solving the dual problem is not guaranteed to give the optimal solution, but can at least be used to construct a feasible solution. In a complete enumeration of small binary test cases (images of pixels for ) we observed that if the problem has a unique solution, then the proposed dual approach finds it. In case the problem has multiple solutions, the dual approach finds the part that is common in all solutions. Based on these experiments we conjecture that this holds in the general case (beyond the small test images). Of course, verifying beforehand if the problem has a unique solution may not be possible.
We test the proposed method on numerical phantoms and real data, showing that the method compares favourably to state-of-the-art reconstruction techniques (Total Variation, DART). The proposed method is also reasonably stable against a moderate amount of noise.
We currently assume the grey levels are known apriori. Extension to multiple (i.e., more than 2) unknown grey levels is possible in the same framework but will be left for future work. To make the method more robust against noise additional regularization may be added.
Appendix A Proofs
A-a Theorem 1
In (5), the has a closed-form expression for general . To see this, let us first denote
We are interested in the infimum value of this function with repsect to . To obtain this, we set the gradient of with respect to to zero
Since is a general matrix, may be rank-deficient. Hence, the optimal value only exists if is in the range of (same as the row space of ) and it is given by
where denotes the Moore-Penrose pseudo-inverse of the matrix. Substituting this value in (17), we get the following:
where denotes the row-space of .
The above dual objective leads to the following maximization problem with respect to the dual variable
As we are only interested in the maximum value of the dual objective, the space of can be constrained to the range of . This is valid as the dual objective is for the outside the range of . Hence, the maximization problem reduced to the following minimization problem:
A-B Corollary 1
Since the search space for the dual variable is constrained to the range of , we can express this variable as , where . Substituting in (12), we get
We can simplify the weighted norm. Let’s consider the economical form of the singular value decomposition (SVD) of
The matrices and are orthogonal matrices, i.e. , where
is an identity matrix. The diagonal matrixhas singular values for on the diagonal, where is the rank of the matrix. It can be easily shown that the pseudo-inverse of the matrix is represented by . Hence, weighted norm reduced to
Also, we note that the pseudo-inverse of matrix is given by . Hence we have . Since , the above weighted norm can be represented in terms of as follows:
The dual problem (18) is now given by
The optimal solution to the above problem is denoted by . Correspondingly, the primal solution related to the dual optimal is
A-C Theorem 2
The primal problem for binary tomography problem with grey levels can be stated as:
where denotes the Heaviside function and is an auxiliary variable. Such problem admits a Lagrangian
where is a Lagrangian multiplier (also known as dual variable) corresponding to the equality constraint. This gives rise to a dual function
Since we already know (refer to (7)), we require the explicit form for . For its computation, we use the componentwise property of the Heaviside function to separate the infimum.
Since the range of Heaviside function is only two values, namely , we get the simple form for :
This infimal value is attained at . Now the dual problem reads
We note that the last two terms in the dual objective can be compactly represented by
where is known as an asymmetric one-norm. The optimal point of the problem (19) is denoted by and the corresponding primal optimal is retrieved using
A-D Theorem 3
The minimization problem for the proximal operator of an asymmetric one-norm function reads
where is a parameter. Since the function is convex, we get the following from the first-order optimality condition :
where is an optimal point of (20), and is a sub-differential of function at . This sub-differential is
Now coming back to the first-order optimality condition in (21), we get the explicit form for optimal solution :
We recognize this function as an asymmetric soft-thresholding function. ∎
This work is part of the Industrial Partnership Programme (IPP) ‘Computational sciences for energy research’ of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO). This research programme is co-financed by Shell Global Solutions International B.V. The second author is financially supported by the Netherlands Organisation for Scientific Research (NWO) as part of research programme 613.009.032.
-  J. Steiner, “Einfache beweise der isoperimetrischen hauptsätze.” Journal für die reine und angewandte Mathematik, vol. 18, pp. 281–296, 1838.
-  J.-M. Carazo, C. Sorzano, E. Rietzel, R. Schröder, and R. Marabini, “Discrete tomography in electron microscopy,” in Discrete Tomography. Springer, 1999, pp. 405–416.
-  B. M. Carvalho, G. T. Herman, S. Matej, C. Salzberg, and E. Vardi, “Binary tomography for triplane cardiography,” in Biennial International Conference on Information Processing in Medical Imaging. Springer, 1999, pp. 29–41.
-  M. Servieres, J. Guédon, and N. Normand, “A discrete tomography approach to pet reconstruction,” in Fully 3D Reconstruction In Radiology and Nuclear Medicine, no. 1, 2003, pp. 1–4.
-  B. Sharif and B. Sharif, “Discrete tomography in discrete deconvolution: Deconvolution of binary images using ryser’s algorithm,” Electronic Notes in Discrete Mathematics, vol. 20, pp. 555–571, 2005.
-  M. Balaskó, A. Kuba, A. Nagy, Z. Kiss, L. Rodek, and L. Ruskó, “Neutron-, gamma-and x-ray three-dimensional computed tomography at the budapest research reactor site,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 542, no. 1-3, pp. 22–27, 2005.
-  S. Van Aert, K. J. Batenburg, M. D. Rossell, R. Erni, and G. Van Tendeloo, “Three-dimensional atomic imaging of crystalline nanoparticles,” Nature, vol. 470, no. 7334, p. 374, 2011.
-  K. J. Batenburg, S. Bals, J. Sijbers, C. Kübel, P. Midgley, J. Hernandez, U. Kaiser, E. Encina, E. Coronado, and G. Van Tendeloo, “3d imaging of nanomaterials by discrete tomography,” Ultramicroscopy, vol. 109, no. 6, pp. 730–740, 2009.
-  R. J. Gardner, P. Gritzmann, and D. Prangenberg, “On the computational complexity of reconstructing lattice sets from their X-rays,” Discrete Mathematics, vol. 202, no. 1–3, pp. 45–71, 1999. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0012365X98003471
-  A. E. Yagle, “An algebraic solution for discrete tomography,” in Discrete Tomography. Springer, 1999, pp. 265–284.
-  L. Hajdu and R. Tijdeman, “Algebraic aspects of discrete tomography,” Journal fur die Reine und Angewandte Mathematik, vol. 534, pp. 119–128, 2001.
-  S. Matej, A. Vardi, G. T. Herman, and E. Vardi, “Binary tomography using gibbs priors,” in Discrete Tomography. Springer, 1999, pp. 191–212.
-  M. T. Chan, G. T. Herman, and E. Levitan, “Probabilistic modeling of discrete images,” in Discrete Tomography. Springer, 1999, pp. 213–235.
-  T. Frese, C. A. Bouman, and K. Sauer, “Multiscale bayesian methods for discrete tomography,” in Discrete Tomography. Springer, 1999, pp. 237–264.
-  Y. Censor and S. Matej, “Binary steering of nonbinary iterative algorithms,” in Discrete tomography. Springer, 1999, pp. 285–296.
-  Y. Vardi and C.-H. Zhang, “Reconstruction of binary images via the em algorithm,” in Discrete Tomography. Springer, 1999, pp. 297–316.
-  T. Capricelli and P. Combettes, “A convex programming algorithm for noisy discrete tomography,” in Advances in Discrete Tomography and Its Applications. Springer, 2007, pp. 207–226.
-  T. Schüle, C. Schnörr, S. Weber, and J. Hornegger, “Discrete tomography by convex–concave regularization and dc programming,” Discrete Applied Mathematics, vol. 151, no. 1-3, pp. 229–243, 2005.
-  E. Barcucci, S. Brunetti, A. Del Lungo, and M. Nivat, “Reconstruction of lattice sets from their horizontal, vertical and diagonal X-rays,” Discrete Mathematics, vol. 241, no. 1-3, pp. 65–78, 2001.
-  K. J. Batenburg and J. Sijbers, “Dart: a practical reconstruction algorithm for discrete tomography,” IEEE Transactions on Image Processing, vol. 20, no. 9, pp. 2542–2553, 2011.
-  G. T. Herman and A. Kuba, Discrete tomography: Foundations, algorithms, and applications. Springer Science & Business Media, 2012.
-  R. Rockafellar, “Convex analysis,” 1970.
-  M. Slater, “Lagrange multipliers revisited,” Cowles Foundation for Research in Economics, Yale University, Tech. Rep., 1959.
-  R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society. Series B (Methodological), pp. 267–288, 1996.
-  I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, vol. 57, no. 11, pp. 1413–1457, 2004.
-  A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM journal on imaging sciences, vol. 2, no. 1, pp. 183–202, 2009.
S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed
optimization and statistical learning via the alternating direction method of
Foundations and Trends® in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011.
-  T. Goldstein and S. Osher, “The split bregman method for l1-regularized problems,” SIAM journal on imaging sciences, vol. 2, no. 2, pp. 323–343, 2009.
K. J. Arrow, L. Hurwicz, and H. Uzawa, “Studies in linear and non-linear programming,” 1958.
-  A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” Journal of mathematical imaging and vision, vol. 40, no. 1, pp. 120–145, 2011.
-  M. Grant, S. Boyd, and Y. Ye, “Cvx: Matlab software for disciplined convex programming,” 2008.
-  E. Barcucci, A. Del Lungo, M. Nivat, and R. Pinzani, “Reconstructing convex polyominoes from horizontal and vertical projections,” Theoretical Computer Science, vol. 155, no. 2, pp. 321–347, 1996.
-  W. van Aarle, W. J. Palenstijn, J. De Beenhouwer, T. Altantzis, S. Bals, K. J. Batenburg, and J. Sijbers, “The astra toolbox: A platform for advanced algorithm development in electron tomography,” Ultramicroscopy, vol. 157, pp. 35–47, 2015.
-  C. C. Paige and M. A. Saunders, “Lsqr: An algorithm for sparse linear equations and sparse least squares,” ACM Transactions on Mathematical Software (TOMS), vol. 8, no. 1, pp. 43–71, 1982.
-  N. Otsu, “A threshold selection method from gray-level histograms,” IEEE transactions on systems, man, and cybernetics, vol. 9, no. 1, pp. 62–66, 1979.
-  D. C. Liu and J. Nocedal, “On the limited memory bfgs method for large scale optimization,” Mathematical programming, vol. 45, no. 1-3, pp. 503–528, 1989.
-  T. A. Bubba, M. Juvonen, J. Lehtonen, M. März, A. Meaney, Z. Purisha, and S. Siltanen, “Tomographic x-ray data of carved cheese,” arXiv preprint arXiv:1705.05732, 2017.