I Introduction
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 [1]. 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, atomicresolution electron microscopy, medicine and material science (for example,[2, 3, 4, 5, 6, 7, 8]).
Ia 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 NPhard problem when [9].
Due to the presence of noise, the system of equations may not have a solution, and the problem is sometimes formulated as a constrained leastsquares problem
(1) 
Obviously, the constraints are nonconvex and solving (1) exactly is not trivial. Next, we briefly discuss some existing approaches for solving it.
IB 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 methods
typically construct a probability density function on the space of discrete images, allowing one to sample images and use Markov Chain Monte Carlo type methods to find a solution
[12, 13, 14]. These methods are very flexible but may require a prohibitive number of samples when applied to largescale datasets.Relaxation methods are based on some form convex or nonconvex 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 largescale problems, getting them to converge to the correct binary solution can be challenging.
The heuristic
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 [21].
IC 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 leastsquares 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 zerovalued 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 smallscale () 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 (zerovalued). 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 nonuniqueness 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 wellchosen numerical experiments on synthetic and real data, we show that our new approach is competitive for practical applications in Xray 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 smallscale 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 leastsquares binary tomography problem can then be formulated as:
(2)  
subject to 
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 :
(3) 
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
(4) 
As the dual function is always concave [22], 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 [23] 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 we
can solve the primary problem via its dual is beyond the scope of this paper.It turns out we can obtain a closedform expression for . Before presenting the general form of , we first present a detailed derivation for invertible to provide some insight.
Iia Invertible
The Lagrangian is separable in terms of and . Hence, we can represent the dual function as the sum of two functions, and .
(5) 
First, we consider . Setting the gradient to zero we find the unique minimizer:
(6) 
Substituting back in the expression and rearranging some terms we arrive at the following expression for :
(7) 
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
(8) 
Hence, the dual function for the Lagrangian in (3) takes the following explicit form:
and leads to the following dual problem to (2).
(9) 
This optimization problem is famously known as the least absolute shrinkage and selection operator (LASSO) [24]
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 scaled
norm. 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 in
need 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 onedimensional example with :
(10) 
The solution to this problem is given by . The corresponding dual problem is
(11) 
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 welldefined. 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.
IiB Main results
We state the main results below. The proofs for these statements are provided in the Appendix section.
Theorem 1.
The dual objective of (2) for general is given by
where denotes the pseudoinverse and is the rowspace of (i.e., the range of ). This leads to the following dual problem
(12) 
Corollary 1.
Remark 2.
For and full rank, we have and the formulation (13) simplifies to
(14) 
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 .
Theorem 2.
The dual problem for a binary tomography problem with grey levels is given by:
(15) 
where is an asymmetric onenorm. 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.
IiC Solving the dual problem
When is invertible, the dual formulation (12) can be readily solved using a proximal gradient algorithm ([25, 26]):
(16) 
where and the soft thresholding operator is applied componentwise 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 dataconsistent manner.
If is invertible we have and it seems more natural to solve (14) instead. Due to the appearance of in the onenorm, it is no longer straightforward to apply a proximal gradient method. A possible strategy is to replace the onenorm 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 gradientdescent algorithm.
We also note that splitting methods can be used to solve (14). For example, the alternating direction method of multipliers (ADMM) [27] and/or splitBregman method [28]. Another class of method that can solve (14) are the primaldual methods (e.g., ArrowHurwicz primaldual algorithm [29], ChambollePock algorithm [30]). 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 onenorm in the following theorem.
Theorem 3.
The proximal operator for an asymmetric onenorm function
with , is given by
where , and is an asymmetric softthresholding function
Proof.
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 NPhard. For small we can simply enumerate all possible images, find all solutions in each case by a bruteforce 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 [31]. 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.
total  unique  multiple  

2  16  14/14  2/2 
3  512  230/230  282/282 
4  65536  6902/6902  58541/58634 
2  16  16/16  0/0 
3  512  496/496  16/16 
4  65536  54272/54272  10813/11264 
2  16  16/16  0/0 
3  512  512/512  0/0 
4  65536  65024/65024  512/512 
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 NPhard. 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,dconvex^{1}^{1}1For h,v,dconvexity 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,dconvex but still permit a unique solution, see figure 4. Such images are also retrieved using our dual approach.
Iv Numerical experiments  Xray tomography
In this section, we present numerical results for limitedangle Xray tomography on a few numerical phantoms and an experimental Xray dataset. First, we describe the phantoms and the performance measures used to compare our proposed dual approach to several stateoftheart 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 [33].
Iva Phantoms
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.
IvB Tests
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 socalled 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.
IvC Comparison with other reconstruction methods
There exist a vast amount of reconstruction methods for tomography. Here, we consider the following three:
 LSQR
 TV

: The totalvariation 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 ChambollePock method [30] to solve the above optimization problem with nonnegativity 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 TVreconstructed image through Otsu’s thresholding algorithm.
 DART

: We use the method described in [20] 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.
 DP
IvD Performance Measures
In order to evaluate the performance of reconstruction methods, we use the following two criteria
 RMS

The rootmeansquare error
measures how well the forwardprojected 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.
 JI

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.
IvE Experimental data setup
We use the experimental Xray projection data of a carved cheese slice provided by the FIPS group [37]. Figure 6 shows a highresolution filtered backprojection reconstruction of the data. The cheese contains the letters C and T and the object is (approximately) binary with two grey levels corresponding to calciumcontaining 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 Xray transform. We perform two sets of experiments: (1) Sparse sampling with 15 angles ranging from to and limitedangle using 15 projections from to .
IvF 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.
Test  Model  LSQR  TV  DART  DP  

RMS  JI  RMS  JI  RMS  JI  RMS  JI  
45  P1  31.44  0.9967  19.20  0.9997  31.54  0.9967  17.56  1.0000 
P2  26.20  0.9959  11.98  1.0000  36.40  0.9901  11.98  1.0000  
P3  48.07  0.9899  24.53  0.9981  52.42  0.9877  16.37  1.0000  
P4  11.35  0.9978  5.01  1.0000  13.61  0.9956  5.01  1.0000  
20  P1  54.30  0.9781  34.50  0.9888  24.06  0.9954  12.07  1.0000 
P2  60.51  0.9577  12.75  0.9986  16.99  0.9961  8.91  1.0000  
P3  54.22  0.9760  26.98  0.9938  40.18  0.9885  11.07  1.0000  
P4  20.17  0.9927  3.35  1.0000  6.84  0.9978  3.35  1.0000  
10  P1  122.25  0.9312  62.41  0.9615  22.39  0.9914  8.29  0.9999 
P2  254.42  0.7601  80.98  0.9110  17.02  0.9904  11.20  0.9967  
P3  73.76  0.9485  38.59  0.9811  28.74  0.9897  7.43  1.0000  
P4  42.66  0.9700  5.42  0.9985  5.88  0.9971  2.36  1.0000  
5  P1  269.68  0.8241  77.51  0.8694  23.65  0.9734  52.76  0.9065 
P2  370.80  0.6322  158.46  0.7129  38.79  0.7613  72.82  0.7391  
P3  79.98  0.9089  64.88  0.9335  22.49  0.9846  21.49  0.9757  
P4  79.69  0.8637  22.85  0.9581  4.65  0.9956  1.68  1.0000 
IvG 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 .
Test  Model  LSQR  TV  DART  DP  

RMS  JI  RMS  JI  RMS  JI  RMS  JI  
150  P1  158.27  0.9157  101.94  0.9389  23.61  0.9910  5.41  1.0000 
P2  215.94  0.7949  79.64  0.9127  16.39  0.9917  4.99  1.0000  
P3  82.41  0.9447  53.29  0.9673  28.54  .9899  5.32  0.9998  
P4  54.90  0.9555  4.65  0.9985  5.33  0.9971  9.45  0.9920  
120  P1  199.67  0.8836  141.67  0.9145  18.68  0.9927  6.66  0.9996 
P2  189.32  0.7780  146.46  0.8304  22.00  0.9798  13.30  0.9926  
P3  96.15  0.9291  70.45  0.9499  27.4086  0.9856  3.40  1.0000  
P4  68.23  0.9170  9.11  0.9934  7.06  0.9949  0.75  1.0000  
105  P1  213.17  0.8695  164.62  0.8965  20.10  0.9924  7.67  0.9993 
P2  231.76  0.7641  182.08  0.8118  21.70  0.9798  16.97  0.9846  
P3  105.54  0.9227  79.85  0.9431  28.27  0.9846  3.65  1.0000  
P4  84.50  0.8965  8.05  0.9942  7.64  0.9949  0.93  1.0000  
90  P1  258.58  0.8493  293.04  0.8517  19.09  0.9925  17.96  0.9920 
P2  205.27  0.7533  192.45  0.8094  22.72  0.9819  18.19  0.9849  
P3  127.09  0.9054  99.23  0.9350  30.99  0.9834  4.04  1.0000  
P4  128.16  0.8288  27.09  0.9672  7.79  0.9949  7.19  0.9949 
IvH 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.
IvI 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. Figure
12 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.









V Conclusion
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 stateoftheart 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
Aa Theorem 1
Proof.
In (5), the has a closedform expression for general . To see this, let us first denote
(17) 
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 rankdeficient. 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 MoorePenrose pseudoinverse of the matrix. Substituting this value in (17), we get the following:
where denotes the rowspace of .
Now we return to the dual objective in equation (5). Substituting the explicit forms for from above and from equation (8), we get the expression for the dual objective:
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:
∎
AB Corollary 1
Proof.
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
(18) 
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 matrix
has singular values for on the diagonal, where is the rank of the matrix. It can be easily shown that the pseudoinverse of the matrix is represented by . Hence, weighted norm reduced toAlso, we note that the pseudoinverse 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
∎
AC Theorem 2
Proof.
The primal problem for binary tomography problem with grey levels can be stated as:
subject to 
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
(19) 
We note that the last two terms in the dual objective can be compactly represented by
where is known as an asymmetric onenorm. The optimal point of the problem (19) is denoted by and the corresponding primal optimal is retrieved using
∎
AD Theorem 3
Proof.
The minimization problem for the proximal operator of an asymmetric onenorm function reads
(20) 
where is a parameter. Since the function is convex, we get the following from the firstorder optimality condition [22]:
(21) 
where is an optimal point of (20), and is a subdifferential of function at . This subdifferential is
Now coming back to the firstorder optimality condition in (21), we get the explicit form for optimal solution :
We recognize this function as an asymmetric softthresholding function. ∎
Acknowledgment
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 cofinanced 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.
References
 [1] J. Steiner, “Einfache beweise der isoperimetrischen hauptsätze.” Journal für die reine und angewandte Mathematik, vol. 18, pp. 281–296, 1838.
 [2] 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.
 [3] 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.
 [4] 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.
 [5] 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.
 [6] M. Balaskó, A. Kuba, A. Nagy, Z. Kiss, L. Rodek, and L. Ruskó, “Neutron, gammaand xray threedimensional 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. 13, pp. 22–27, 2005.
 [7] S. Van Aert, K. J. Batenburg, M. D. Rossell, R. Erni, and G. Van Tendeloo, “Threedimensional atomic imaging of crystalline nanoparticles,” Nature, vol. 470, no. 7334, p. 374, 2011.
 [8] 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.
 [9] R. J. Gardner, P. Gritzmann, and D. Prangenberg, “On the computational complexity of reconstructing lattice sets from their Xrays,” Discrete Mathematics, vol. 202, no. 1–3, pp. 45–71, 1999. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0012365X98003471
 [10] A. E. Yagle, “An algebraic solution for discrete tomography,” in Discrete Tomography. Springer, 1999, pp. 265–284.
 [11] L. Hajdu and R. Tijdeman, “Algebraic aspects of discrete tomography,” Journal fur die Reine und Angewandte Mathematik, vol. 534, pp. 119–128, 2001.
 [12] S. Matej, A. Vardi, G. T. Herman, and E. Vardi, “Binary tomography using gibbs priors,” in Discrete Tomography. Springer, 1999, pp. 191–212.
 [13] M. T. Chan, G. T. Herman, and E. Levitan, “Probabilistic modeling of discrete images,” in Discrete Tomography. Springer, 1999, pp. 213–235.
 [14] T. Frese, C. A. Bouman, and K. Sauer, “Multiscale bayesian methods for discrete tomography,” in Discrete Tomography. Springer, 1999, pp. 237–264.
 [15] Y. Censor and S. Matej, “Binary steering of nonbinary iterative algorithms,” in Discrete tomography. Springer, 1999, pp. 285–296.
 [16] Y. Vardi and C.H. Zhang, “Reconstruction of binary images via the em algorithm,” in Discrete Tomography. Springer, 1999, pp. 297–316.
 [17] 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.
 [18] 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. 13, pp. 229–243, 2005.
 [19] E. Barcucci, S. Brunetti, A. Del Lungo, and M. Nivat, “Reconstruction of lattice sets from their horizontal, vertical and diagonal Xrays,” Discrete Mathematics, vol. 241, no. 13, pp. 65–78, 2001.
 [20] 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.
 [21] G. T. Herman and A. Kuba, Discrete tomography: Foundations, algorithms, and applications. Springer Science & Business Media, 2012.
 [22] R. Rockafellar, “Convex analysis,” 1970.
 [23] M. Slater, “Lagrange multipliers revisited,” Cowles Foundation for Research in Economics, Yale University, Tech. Rep., 1959.
 [24] R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society. Series B (Methodological), pp. 267–288, 1996.
 [25] 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.
 [26] A. Beck and M. Teboulle, “A fast iterative shrinkagethresholding algorithm for linear inverse problems,” SIAM journal on imaging sciences, vol. 2, no. 1, pp. 183–202, 2009.

[27]
S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed
optimization and statistical learning via the alternating direction method of
multipliers,”
Foundations and Trends® in Machine Learning
, vol. 3, no. 1, pp. 1–122, 2011.  [28] T. Goldstein and S. Osher, “The split bregman method for l1regularized problems,” SIAM journal on imaging sciences, vol. 2, no. 2, pp. 323–343, 2009.

[29]
K. J. Arrow, L. Hurwicz, and H. Uzawa, “Studies in linear and nonlinear programming,” 1958.
 [30] A. Chambolle and T. Pock, “A firstorder primaldual algorithm for convex problems with applications to imaging,” Journal of mathematical imaging and vision, vol. 40, no. 1, pp. 120–145, 2011.
 [31] M. Grant, S. Boyd, and Y. Ye, “Cvx: Matlab software for disciplined convex programming,” 2008.
 [32] 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.
 [33] 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.
 [34] 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.
 [35] N. Otsu, “A threshold selection method from graylevel histograms,” IEEE transactions on systems, man, and cybernetics, vol. 9, no. 1, pp. 62–66, 1979.
 [36] D. C. Liu and J. Nocedal, “On the limited memory bfgs method for large scale optimization,” Mathematical programming, vol. 45, no. 13, pp. 503–528, 1989.
 [37] T. A. Bubba, M. Juvonen, J. Lehtonen, M. März, A. Meaney, Z. Purisha, and S. Siltanen, “Tomographic xray data of carved cheese,” arXiv preprint arXiv:1705.05732, 2017.
Comments
There are no comments yet.