Abstract
Discrete tomography refers to a class of reconstruction methods adapted to discrete-valued images. Many different approaches have been investigated to address the binary case, when a two-phase object is considered. This reconstruction problem is very important in medical or material applications where it is crucial to reduce the number of projections. In this paper, we address the problem of binary image reconstruction for X-ray CT imaging from a small number of projections. We propose a TV (Total Variation) regularization approach and compare the results obtained with or without an additional box convex constraint. The schemes are applied to a simple disk image and to more complex bone cross-sections for various noise levels. The minimization of the regularization functional is performed with the state-of-the-art ADMM (Alternate Direction Minimization Method) algorithm. The methods perform equally well on a simple disk image. The additional box convex constraints improves the reconstruction results for complex structures with fine details.
Introduction
Tomography reconstruction from a limited number of projections is an important problem in X-ray Computed Tomography. It is particularly crucial when imaging a moving organ, such as the beating heart or when the irradiation dose has to be reduced as for in vivo investigation of the bone micro-structure. Therefore, it is interesting to study the potential of new optimization schemes to reconstruct images from a limited number of projections. A number of reconstruction methods have been proposed to achieve low dose CT. These methods are generally iterative and rely on a sparsity prior which may be applied in the image domain or after a sparsifying transform such as a wavelet transform [1]. Recently, with the development of compressive sensing approaches, a number of algorithms based on TV regularization scheme have been investigated [2–6]. The TV regularization is well-known to preserve edge and it provides good reconstruction images with sparse view sampling.
Other approaches based on discrete tomography make the assumption that the image is a discrete valued function [7]. Such discrete tomography methods are very useful for X-ray tomography of industrial objects [8] or bone imaging [9, 11]. A number of methods have specifically been developed to address the binary case when a two-phase object is considered. In the framework of the linear imaging model for X-ray projection, the data are considered as the integrals of the attenuation coefficient. The binary tomography problem is associated with an under-determined linear system of equation Rf = p where R is the linear Radon projection operator, (f i ) 1≤i≤n the pixel values of the image with a binary constraint f = (f1, …, f n ) ∈ {0, 1} n and p the projection data. Only noisy projection values p δ are given for the noise level δ such that ∥p - p δ ∥ ≤ δ. Various approaches based on discrete algebraic reconstruction techniques [12, 15] have been proposed to solve this problem. Some methods minimize a functional that incorporates a data term and a binary constraint, withstochastic techniques [16] or convex analysis optimization [17, 18]. Markov random fields have also been much used [19]. Recently, a Belief Propagation reconstruction approach has been proposed [20].
The binary tomography problem is ill-posed and must be regularized. The TV regularization has been often used and this method gives good results [20]. This regulazation method improves the classical Tikhonov regularization which gives poor results for the reconstruction of non-smooth solutions [21]. It can be extended to the general discrete case in which the image to be reconstructed can take several discrete values. Yet, the TV regularization methods favor circular structures since it tend to minimize the perimeter of boundaries between the reconstructed regions. It is not clear if it is suited for more complex and elongated or tubular structures. These kinds of images are challenging due to the fine structures and the multi-scale details in the image. It is also necessary to handle the non-convex binary constraint. Some convexification methods have been proposed for segmentation tasks. These approaches can be used for the binary tomography problem but it is necessary to evaluate the improvement of the reconstruction results.
On the other hand, in order to evaluate the effectiveness of a regularization methods, it is necessary to choose an optimal regularization parameter. Several methods have been studied for the optimal selection of the regularization parameter for Tikhonovor TV regularization. Methods like the Morozov discrepancy principle and the Unbiased Predictive Risk Estimator are based on a knowledge of the noise level δ [22–25]. The L-curve method and the Generalized Cross-Validation can be used without the knowledge of this noise level [26–28].
The main objective of this work is to test TV regularization based reconstruction methods, without and with additional convex constraints, on a very simple disk image and on more complex images, obtained from CT bone cross-sections. The same methodology of choice of the regularization parameter will be used in all the cases, based on the Morozov discrepancy principle. The TV regularization functional will minimized by the state-of-the-art ADMM algorithm. Some comparison of the reconstruction results obtained with some state-of-the-art binary tomography will be also presented.
The outline of this paper is the following. After the introduction, the second section of this paper deals with the TV regularization based methods without or with convex constraints and with the ADMM minimization methodology. The numerical results obtained on a simple disk or on noisy bone CT cross-sections are reported and discussed in the last section. Concluding remarks are given in the last section.
TV regularization and ADMM approach
A common way to regularize the binary tomography problem is to construct a regularization functional E (f) with a data fidelity term that measures the consistency between the estimates and the measurements and a regularization term J (f) that imposes an a priori constraint on the solution. The data-fitting term is usually based on the L2 norm and the regularization functional can then be written as:
The parameter μ is the regularization parameter balancing the contribution of the two terms. The measured projection data p δ is the approximation of the correct data p, corresponding to the true solution f* with Rf* = p. The noisy data p δ are assumed to be corrupted by noise with a noise level δ, satisfying ∥p δ - p ∥ 2 ≤ δ
In the classical Tikhonov regularization, the regularization term is given by , where D is a differential operator. TV regularization was introduced by Rudin et al. [29] to solve the noise removal problem. It has been applied to various image processing problems [30–32] and it is very successful to preserve edges. Let Ω be a bounded open subset of , for an image belonging to the first-order Sobolev space, f ∈ H1 (Ω), this regularization is based on computing the L1 norm of the gradient:
The results obtained with the isotropic TV were very similar or slightly better than the ones achieved with the anisotropic TV for most cases investigated. We will present only the results obtained with the isotropic norm in this work.
The binary constraints leads to a non-convex inverse problem. Convexified models obtained by relaxation of the binary constraint have often been considered for segmentation tasks [33, 34]. In this study, we use the same type of methods and the function f to be reconstructed is thus allowed to take values continuously from [0, 1]. The convex constraints can be included in the regularization functional [31, 32]. The following optimization problems without and with convex constraints f ∈ C = [0, 1]
n
have been considered:
Various numerical methods have been used to solve the TV regularized deconvolution problem including partial differential equations, splitting or primal dual methods [35–37]. Results of extensive numerical experiments show that algorithms based on the ADMM and an augmented Lagrangian functional are among the state-of-the-art methods to minimize the TV regularization functional [30–32]. Algorithms based on the ADMM (SALSA and C-SALSA) have been proposed to solve a number of image processing tasks, such as image inpainting and deblurring [31, 32]. Our problem is thus formulated as a minimization problem of the ADMM form with a series of linear constraints. In order to include convex constraints, f ∈ C = [0, 1]
n
, the following augmented Lagrangian is considered:
The Lagrange multipliers (λ
i
) 1≤i≤n, λ
C
are vectors in and . For each pixel i, represents the first-order finite difference at pixel i in both horizontal and vertical directions, (g
i
) 1≤i≤n and h are the auxiliary unknowns corresponding to the gradient and the convex constraint. The ADMM algorithm searches for the saddle point of the augmented Lagrangian by iterating the following equations successively:
The h
k
update is:
The sequence which is generated by the ADMM algorithm converges to a Kuhn-Tucker point of problem (P2), , if (P2) has one. If (P2) does not have an optimal solution, then at least one of the sequences diverges. If the additional convex constraint is not included in the regularization functional, the additional unknowns h and λ C are not used and Equation 9 and Equation 12 are not considered.
In this section, we present the simulation details and the results obtained with the TV regularization without and with convex constraints.
Simulations details
In our experiments, the projection operator R is taken as the discrete approximation of the Radon transform, which is implemented on Matlab Image Toolbox. The TV regularization methods were applied to two small images of size 256 × 256, which are shown in Fig. 1, and to two big images of size 512 × 512 with low and high bone density, which are shown in Fig. 2 [10, 11].
The first small image is a simple disk image and the small bone images an experimental bone cross-sections reconstructed with 400 projections from FBP (Filtered Back Projections) algorithm and subsequently thresholded. The comparison of the reconstruction results obtained on the simple disk image and on the bone cross-section is useful to understand the effect of the complexity of the geometry. The large images are reconstructed with 729 projections. In these images, the pixel size is 15μm. These images are regarded as the ground-truth images.
In our simulations, the images were reconstructed from a limited number of projection angles M, with M = 20 and 50 Lower projection numbers such as M = 10 were not considered in this work since the reconstruction errors obtained on the boundaries are very large. For all images, the noisy projections p δ were obtained by adding a Gaussian noise with standard deviation σ to the raw projection datap.
This noise level δ can be estimated with δ2 = MN
r
σ2, where N
r
is the number of X-rays per projection. The σ values, the PPSNR values and the noise levels δ are summarized in Table 1 and Table 2 for the small and large images respectively. In order to evaluate the quality of reconstruction, the normalized error E (k) and misclassification rate MR (k) at the iteration k have been calculated as a function of the iteration number k:
The iterations were stopped when the regularization functional stagnates. At the end of the optimization process, E m will denote the minimum error for grey-level images f m obtained at the final iteration m. The final iteration index m is determined by the stopping condition: ∥fm+1 - f m ∥ 2/∥ f m ∥ 2 < 0.0001 .
In order to obtain the best reconstruction results, it is necessary to choose optimal regularization parameters. For the TV regularization method, there are two important parameters: the regularization parameter μ and Lagrangian parameter β. The β parameter controls the speed of convergence. This parameter was chosen to obtain the fastest decrease of the regularization functional.
The reconstructed image f
m
(μ) obtained at the end of the optimization process depends only on the regularization parameterμ. In order to find the best combination of these parameters, we have tested many values of μ and β. Our choice of the optimal μ is based on the Morozov discrepancy principle [23]. In the framework of this principle, the convergence of the regularization method is obtained if the final iterate, f
m
, satisfies the condition:
Finally, the difference map image f
diff
is used to evaluate the quality of binary images by visual inspection. It is defined as the difference between binary image at the end of the iterative algorithms and the ground-truth image f*:
For comparison, the iterative reconstruction algorithm DART [12] was applied to the data. The basic idea is to fix some pixels in the interior or exterior of the object to 0 or 1. A first reconstruction is computed as a starting point using the ART (Algebraic Reconstruction Technique) with 10000 iterations. Then a number of DART iterations are performed. The main loop of the algorithm is repeated until the data term is very close to the noise level δ. Without this stopping criterion, the algorithm diverges. In each inner iteration, 3 SART iterations were performed, updating only the pixels of the boundary region. The reconstructed image was smoothed at each step with a Gaussian filter. DART algorithm was implemented with the ASTRA tomography toolbox [9, 14]. The misclassification rate obtained at the end of the DART process is estimated from the binarized images. As compared to TV, DART is more an heuristic method with more adjustable parameters but is has been shown to give good results with enough SART iterations.
Small images
As an example of the reconstruction results obtained, the binary reconstructed images with 20 projection angles for TV regularization without and with box constraints are shown in Fig. 3 and Fig. 4 for the disk and the bone images respectively. We can see on the difference maps that the reconstruction errors are localized on the boundaries.
The evolution curves of data term (∥ Rf k - p δ ∥), normalized error (E (k)) and misclassification rate (MR (k)) with the iteration number k for small bone image with 20 projection angles are shown in Fig. 5. The evolution curves obtained with TV algorithm with convex constraints decrease faster and converge to smaller values than those without convex constraints.
Table 3, Table 4 and Table 5 summarize the minimum error E m and misclassification rate, MR m , for disk and bone image respectively, for each noise level with 20 and 50 projections obtained with the TV regularization without or with convex constraints.
The values obtained with the DART algorithm for the small bone image are also given. For the disk image, the reconstruction results obtained without convex constraints or with these constraints are very good and they are almost the same. The level lines of the image are circles and are well restored by the TV regularization term which tends to minimize their perimeter. When the noise level δ is high, TV without convex constraints performs better. On the contrary, for the bone image, with a more complex structure, the TV algorithm with convex constraints gives the minimum normalized error E m and misclassification rate MR m for all the noise levels investigated. When the noise level is low, the TV and TV box approaches perform better than DART. When the noise level is high, the DART method achieves better reconstruction results.
Large bone cross-section images
The TV algorithms without and with box constraints were also compared on big bone cross-sections images of size 512 × 512. The binary reconstructed images with TV regularization method without and with box constraints with σ = 3 and M = 20 projection angles are shown in Fig. 6 and Fig. 7 for low and high bone density images respectively. Similarly to the reconstruction results for small images, the reconstruction errors are localized on the boundaries.
The evolution curves of data term (||Rf k - p δ ||) , normalized error (E (k)) and misclassification rate (MR (k)) with the iteration number are shown in Fig. 8 for the big sparse bone image with 20 projection angles. The evolution curves for the big dense bone image are very similar. The evolution curves of the normalized error E (k) and misclassification rate MR (k) of TV algorithm with convex constraints decrease faster and give smaller values
than without convex constraints. It was not possible to find a regularization parameter μ and a Lagrange parameter β to make the data term , with ξ = 0.01, for TV algorithm without convex constraints because a good approximate solution can not be obtained with this regularization method. In our simulations, the smallest constant ξ which satisfies this relation is ξ = 0.5 . Table 6 and Table 7 summarize the minimum error Em and misclassification rate MRm obtained for sparse and dense bone cross-section images with 20 projections. The values achieved with the DART algorithm are given for comparison.
From these tables, we see that the TV with box constraints algorithm gives better reconstruction results than the algorithm without any constraints on both big sparse and dense bone cross-section images. The TV method performs poorly on both bone cross-section images with complex structures and elongated regions. For the two bone cross-sections, the DART method is the worst method.
These bone images are very challenging images due to the presence of fine structures. For a subset ⊂RN, it is well-known that the total variation of the characteristic function of the subset is equal to its perimeter [38], and the TV regularization is thus not the most efficient approach for complex morphologies. On the contrary, the TV regularization scheme with box constraints performs efficiently and it is able to extract many details and fine structures, when the parameter is chosen according to the Morozov principle. The reconstruction errors on the boundaries are much decreased.
Conclusion
In this paper, two TV based algorithms have been compared. The first algorithm is the classical TV regularization method. In the second approach, a convex constraint has been added to enforce the binary condition. The optimal solutions are obtained with the ADMM algorithm. The performance of these algorithms are compared under different levels of noise on two small images (a disk image and CT bone cross-section image) and two big images (sparse and dense CT bone cross-section images) with different number of projections. The regularization parameter are carrefully chosen with the Morozov discrepancy principle. Reconstruction errors are localized on boundaries for all kinds of images. Very similar normalized errors and misclassification rates are obtained on the simple disk image with the two schemes. For the more complex bone image, the reconstruction are much improved with the additional convex constraint. For the large images and the small bone cross-sections for low noise levels, the TV regularization approach with box constraints outperforms the DART method. In future work, the same comparison will be applied to real images. The method will be extended to 3D images and to real images. Some stochastic methods will be investigated to decrease the reconstruction errors on the boundaries.