A convex formulation for Discrete Tomography

07/24/2018 ∙ by Ajinkya Kadu, et al. ∙ 0

Discrete tomography is concerned with the recovery of binary images from a few of their projections (i.e., sums of the pixel values along various directions). To reconstruct an image from noisy projection data, one can pose it as a constrained least-squares problem. As the constraints are non-convex, many approaches for solving it rely on either relaxing the constraints or heuristics. In this paper we propose a novel convex formulation, based on the Lagrange dual of the constrained least-squares problem. The resulting problem is a generalized LASSO problem which can be solved efficiently. It is a relaxation in the sense that it can only be guaranteed to give a feasible solution; not necessarily the optimal one. In exhaustive experiments on small images (2x2, 3x3, 4x4) we find, however, that if the problem has a unique solution, our dual approach finds it. In the case of multiple solutions, our approach finds the commonalities between the solutions. Further experiments on realistic numerical phantoms and an experiments X-ray dataset show that our method compares favourably to alternative approaches, including Total Variation and DART.



There are no comments yet.


page 5

page 6

page 7

page 8

page 9

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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, 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 [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 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 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 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.

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].

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:

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 :


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 [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 closed-form expression for . Before presenting the general form of , we first present a detailed derivation for invertible to provide some insight.

Ii-a Invertible

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) [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 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.

Theorem 1.

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

Remark 1.

In case and has full rank, is invertible and the general form (12) simplifies to (9).

Corollary 1.

The dual problem (12) can be restated as


and the primal solution is recovered through .

Remark 2.

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 .

Theorem 2.

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.

0:  ,
1:  if  then
2:     if  then
4:        return  
5:     else
7:        return  
8:     end if
9:  else
11:     return  
12:  end if
Algorithm 1 Dual problem for various cases

Ii-C Solving the dual problem

When is invertible, the dual formulation (12) can be readily solved using a proximal gradient algorithm ([25, 26]):


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.

Figure 1: Plot of the dual function (gray line) corresponding the the primal objective for (left) and (right) and its approximations (red line) at .

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) [27] and/or split-Bregman method [28]. Another class of method that can solve (14) are the primal-dual methods (e.g., Arrow-Hurwicz primal-dual algorithm [29], Chambolle-Pock 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 one-norm in the following theorem.

Theorem 3.

The proximal operator for an asymmetric one-norm function

with , is given by

where , and is an asymmetric soft-thresholding function


See the Appendix. ∎

Figure 2: Plot of the dual function (gray line) corresponding the the primal objective for (left) and (right) and its smooth approximation (red line) using with .

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 [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
Table I: Summary of complete enumeration experiments.

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).

Figure 3: Example for and (using directons , and ). Two images with the same projections are shown in the top row while the intersection and the results obtained by the pseudo-inverse and the dual problem are shown in the bottom row.

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.

Figure 4: Example of a binary image that is not h,v,d-convex but does permit a unique solution.

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 [33].

Iv-a 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.

Figure 5: Phantom images used in the simulation experiments. (a) Phantom 1, (b) Phantom 2, (c) Phantom 3, (d) Phantom 4.

Iv-B 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 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:


: Least squares QR method described in [34]. We perform a total of 1000 iterations with a tolerance of . We segment the resulting reconstruction using Otsu’s thresholding algorithm [35].


: 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 [30] 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 [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.


: We solve the dual formulation (13) with a smooth approximation of the -norm (as discussed in section II ) using L-BFGS method [36] with a maximum of 500 iterations.

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 [37]. 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 .

Figure 6: The high resolution filtered back-projection reconstruction of the carved cheese from 360 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.

(a) LSQR
(b) TV
(c) DART
(d) DP
Figure 7: Limited projection test I for Phantom 1. Performance of various reconstruction methods with 10 projection angles from 0 to . The results correspond to difference between reconstructed image with true image. Blue pixels denote the missing pixels while red pixels denote the overestimated pixels.
(a) 45
(b) 20
(c) 10
(d) 5
Figure 8: Limited projection test II for Phantom 3. Performance of proposed method vs number of projections. Blue pixels denote the missing pixels while red pixels denote the overestimated pixels.
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
Table II: Limited projection test performance measures

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 .

(a) LSQR
(b) TV
(c) DART
(d) DP
Figure 9: Limited angle test I for Phantom 2. Performance of various reconstruction methods with 10 projection angles from 0 to . Blue pixels denote the missing pixels while red pixels denote the overestimated pixels.
Figure 10: Limited angle test II for Phantom 4. Performance of proposed method vs maximum angle. Blue pixels denote the missing pixels while red pixels denote the overestimated pixels.
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
Table III: Limited angle test performance measures

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.

(a) 50 dB
(b) 40 dB
(c) 30 dB
(d) 20 dB
(e) 50 dB
(f) 40 dB
(g) 30 dB
(h) 20 dB
Figure 11: Noisy Projection test on phantoms 1 and 2. Performance of proposed method vs signal-to-noise ratio. Blue pixels denote the missing pixels while red pixels denote the overestimated pixels.

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. Figure

12 represents the histogram. We obtain a value of 0.00696 for this compound.

Figure 12: Histogram of filtered backprojection image of the carved cheese.

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.

(a) LSQR
(b) TV
(c) DP
(d) LSQR
(e) TV
(f) DP
(g) LSQR
(h) TV
(i) DP
Figure 13: Real Data Test I - Sparse projection tomography. Performance of various methods with different resolutions. Top row corresponds to pixels. Middle row corresponds to pixels. Bottom row corresponds to pixels. Figure below each image denote the histogram.

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.

(a) LSQR
(b) TV
(c) DP
(d) LSQR
(e) TV
(f) DP
(g) LSQR
(h) TV
(i) DP
Figure 14: Real Data Test II - Limited angle tomography. Performance of various methods with different resolutions. Top row corresponds to pixels. Middle row corresponds to pixels. Bottom row corresponds to pixels. Figure below each image denote the histogram.

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 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 .

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:

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 matrix

has 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:

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


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 [22]:


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.


  • [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-, 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.
  • [7] 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.
  • [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 X-rays,” 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. 1-3, 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 X-rays,” Discrete Mathematics, vol. 241, no. 1-3, 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 shrinkage-thresholding 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 l1-regularized 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 non-linear programming,” 1958.

  • [30] 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.
  • [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 gray-level 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. 1-3, pp. 503–528, 1989.
  • [37] 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.