|Computers, Materials & Continua |
Binary Tomography Reconstruction with Limited-Data by a Convex Level-Set Method
1Information and Systems, University of Tsukuba, Tsukuba, Ibaraki, 305-8573, Japan
2Department of Mathematics, Faculty of Science, Sohag University, Sohag, 82524, Egypt
*Corresponding Author: Haytham A. Ali. Email: email@example.com
Received: 03 March 2022; Accepted: 20 April 2022
Abstract: This paper proposes a new level-set-based shape recovery approach that can be applied to a wide range of binary tomography reconstructions. In this technique, we derive generic evolution equations for shape reconstruction in terms of the underlying level-set parameters. We show that using the appropriate basis function to parameterize the level-set function results in an optimization problem with a small number of parameters, which overcomes many of the problems associated with the traditional level-set approach. More concretely, in this paper, we use Gaussian functions as a basis function placed at sparse grid points to represent the parametric level-set function and provide more flexibility in the binary representation of the reconstructed image. In addition, we suggest a convex optimization method that can overcome the problem of the local minimum of the cost function by successfully recovering the coefficients of the basis function. Finally, we illustrate the performance of the proposed method using synthetic images and real X-ray CT projection data. We show that the proposed reconstruction method compares favorably to various state-of-the-art reconstruction techniques for limited-data tomography, and it is also relatively stable in the presence of modest amounts of noise. Furthermore, the shape representation using a compact Gaussian radial basis function works well.
Keywords: Binary tomography; parametric level-set method; inverse problem; shape recovery; Gaussian function; convex optimization
Computed tomography (CT) is a non-invasive imaging technique that reconstructs images of non-accessible or non-visible objects from a set of projection data [1,2]. Computed tomography has a wide range of applications in diverse fields, such as medicine, science, industrial non-destructive testing, and electron tomography [3–8]. For many of these applications, it is extremely desirable to minimize the number of projections data measured for the object, because minimizing the amount of collected data leads to low-dose imaging and faster imaging. For example, in electron tomography, only a small number of projections, i.e., less than 100, are measured to limit the acquisition time and avoid damaging the sample .
Unfortunately, reducing the number of projection data results in artifacts in the reconstructed image. There exist two classes of reconstruction algorithms in CT. The first class of reconstruction algorithms, known as analytical reconstruction methods, such as Filtered Back Projection (FBP) , requires a large number of projection data acquired throughout the whole angular range to achieve acceptable reconstructions. The second class of algorithms is known as iterative reconstruction methods, such as the Simultaneous Iterative Reconstruction Technique (SIRT) , which allows for the incorporation of prior knowledge about the object into the reconstruction. Thanks to the use of prior knowledge, high-quality reconstructions may be achieved from a relatively small number of projection data or other limited projection data.
Recently, reconstruction algorithms using sparsity of magnitude of intensity gradient as prior knowledge have been used with the name of Total-Variation (TV) to handle sparse-view and limited-angle CT reconstruction . Finally, to further improve image quality or reduce the number of necessary projection data, using the constraint that image intensity is binary has been investigated with the name of discrete (binary) tomography [13–15]. For example, Discrete Algebraic Reconstruction Technique (DART) algorithm for discrete tomography has been developed, which provides high-quality reconstructions with a very simple modification of Algebraic Reconstruction Technique (ART) method [16–18]. Despite the fact that DART or its variations have been effectively used for electron tomography [19,20], micro-CT [21,22], and magnetic resonance imaging (MRI) , the use of this technique in practical applications is still limited because it suffers from a longer computation time.
To reduce computing time or to increase reconstruction quality for a given computation time, we consider a parametric level-set function that may be defined using a collection of basic functions. This parametric level-set function is likely to behave well since there is no need for re-initialization, unlike the classical level-set method, and we can easily derive the corresponding evolution equation based on the underlying parameters. According to Kilmer et al. , the polynomial basis can be applied to diffuse optical tomography by expressing the level set function as a fixed-order polynomial, and it evolves by updating the polynomial coefficients at every iteration. Additionally, the parametric level-set functions have been used in applications related to image processing. In , Gelas et al. used compactly supported radial basis functions (RBFs) as a basis function to reduce the dimension of the problem through the parametric representation and computational cost due to the sparsity provided by this type of basis function. Moreover, Bernard et al.  parameterized the level-set function by using B-splines as the basis function.
In this paper, we propose a novel representation of the parametric level-set (PALS) for shape-based inverse problems. There are some similarities to the previous methods in our approach, such as the use of a basis function to represent the parametric level-set function . However, we chose a different representation model, in which a Gaussian function is used as the basis function. Furthermore, unlike previous works that used basis functions to represent the parametric level-set, our representation has several advantages, including high accuracy, easy applicability to problems of a high dimension, and the intrinsic smoothness of Gaussian basis functions, which guarantees the smoothness of the solution.
Besides the above advantages, the proposed technique is based totally on convex optimization, so it has a consistent performance with a small number of projection data by solving the issue of the cost function’s local minimum. Moreover, the number of parameters in the level-set function is often far smaller than the number of image pixels. Finally, the results of numerical experiments with synthetic images and real X-ray CT projection data show that the proposed method compares favorably to other state-of-the-art reconstruction techniques, and that it is also relatively stable in the presence of modest amounts of noise.
This paper is organized as follows. First, we describe how to formulate the shape-based reconstruction problems in tomography dealt with in this paper. Second, we present the mathematical basis of the parametric level-set technique, and various practical issues are discussed. Then, we give the experimental results of our method for some examples and real data, comparing it to other possible reconstruction methods. Finally, we conclude the paper.
2 Mathematical Formulation
In this study, we consider the discrete tomography reconstruction problem, which may be expressed mathematically as follows. To begin, any image is defined on grid of pixels, which is expressed as . From the image, we can calculate the projection data by taking a ray-sum along a set of straight lines measuring the projection data. The resulting projection data is represented by using the following linear system of equations.
denotes the value of each pixel in the image, denotes the weighted sum of pixel values along the ray as shown in Fig. 1a, and A denotes the projection matrix.
The image reconstruction problem is an inverse problem in which the image u is recovered from the projection data y using the matrix A, i.e.,
In many cases of CT imaging, Eq. (1) is underdetermined because the number of measure projection data is substantially smaller than the number of unknowns N. Also, due to the noise in the projection data or errors in the projection model, the system of equations is inconsistent, such that the solution may be far from the truth. To address these issues, the reconstruction problem is typically expressed as a least-squares problem as
For example, to solve this optimization problem, the following Newton-like algorithm can be used to find the image u .
where denotes the step-size parameter, denotes an approximation of the Hessian matrix of f calculated at , and the gradient of f is given by
3 Basic Principle of Level-Set Method
In the areas of geometric inverse problems, such as shape optimization, image segmentation, and interface tracking, level-set methods have received a lot of attention due to their ability to adapt to topological changes in the region structures. In , Sethian and Osher introduced the classical level-set methods for solving the Hamiltonian-Jacobi equation, which is also known as a level-set equation.
where denotes the level-set function and F denotes the normal velocity and this velocity is often calculated from the gradient of the cost function with respect to the parameters of the model [30,31].
3.1 Shape-Based Approach
In a large class of shape-based inverse problems [31–33], the model of interest is composed of two classes, namely foreground and background. Then, our objective is to determine the boundary between these classes and the characteristics that define the values of u in each class. Therefore, with any closed domain with boundary , is modeled as
where denotes the value inside D and denotes the value outside D. Therefore, the corresponding model is represented as follows.
where is the indicator function of D. From Eq. (3), cannot be differentiated on the boundary unlike what was stated earlier and finding the boundary in this case is the main objective. To solve this issue, the level-set function is known to be particularly useful. The level-set approach represents the region boundary by by using the zero contour of a higher-dimensional function . The function is called the level-set function [34,35], and the relation between and D is expressed as
By using Eq. (4), the indicator function can be represented as , where H is the Heaviside function defined by As a result, can be represented by
3.2 Parametric Level-Set
In many existing shape-based approaches, the level-set function is selected as a signed distance function . As an alternative to this classical approach of level-set functions, consider that the level-set function is not only a function of x but also a function of parameters . Based on this parameterization, the level-set function can be defined as a parametric level-set (PALS), and we can modify Eq. (4) as follows
where c is a positive small value. Now, we can say that for some specified α, can explicitly determine the level-set function for the entire domain . By updating α, the required evolution of is performed to solve the underlying inverse problems. Based on Eq. (6), the model in Eq. (5) is expressed as
A PALS function is an example of a basis expansion using known basis functions and unknown weights. As a result of this statement, the level-set function may be represented parametrically by using the set of RBFs as follows .
where N denotes the number of RBFs, and are unknown parameters of that determine the weight attached to each The meaning of Eq. (8) is that is represented by a linear combination of basis functions taken with the weights There are several basis set choices available, including Gaussian, multi-quadric, poly-harmonic splines, and thin spline polynomial. Here, the Gaussian radial-basis function (GRBF) is used, and it is expressed as
where is the width of Gaussian, is the center of GRBF, and is the Euclidean norm.
4 Shape Representation and Algorithms
In this section, we describe the method for reconstructing the image using the parametric level-set and how it evolves when parameters and GRBFs are updated. Consider the level-set represented in Fig. 2. In Fig. 2, we illustrate how a smooth shape discretized on a grid of can be reconstructed.
Next, we discuss how to select the Gaussian width parameter of GRBF in the PALS reconstruction. For this purpose, we rewrite the GRBF basis function of Eq. (9) in the alternative form as
where Obviously, by changing the parameter value , i.e., the width of GRBF, representation power of the method will be changed significantly. The spreading parameter is very important. In particular, has a large effect on the recovery accuracy (representation power) when the number of GRBFs placed in the 2-D plane is small. In our experiments, we discretized our model on a cartesian grid with pixels, and the nodes to place the GRBFs were selected as a sparse cartesian grid with pixels, i.e., N = 784. The PALS reconstruction process was initialized with a very rough representation with N = 10 GRBFs and . As shown in Fig. 2, this model is finally discretized on a Cartesian grid. Therefore, by using Eq. (7) in Eq. (2), the discretized reconstruction problem for the parametric level-set method can be expressed as
To calculate the gradient of Eq. (10), taking the derivative of Eq. (7) with respect to gives
where is Dirac delta function [38,39]. Using the chain rule, we have
In Eq. (12), B is the matrix constructed by the basis vectors of GRBF, i.e., matrix converting the parameters into the level-set function . Using Eq. (12), the gradient of cost function f is obtained as follows.
where the residual operator is defined by
In our implementation, we solve the minimization problem of Eq. (10) with respect to the level-set parameters by using the quasi-Newton method with the L-BFGS method to approximate the Hessian . The iteration formula is expressed as
The final reconstruction algorithm is summarized in Algorithm below.
From Eq. (5), during the optimization process, we need to compute Heaviside function from the level-set function . As a smooth approximation of the Heaviside function, we can use the following regularization of .
However, Eq. (15) is not a good choice due to the following reason. Considering the nature of parametric level-set function, Eq. (7) shows that for every and for any . According to Fig. 3, using the regularization of the Heaviside function means that the resulting parametric level-set function by using GRBFs takes a relatively large positive value in the region D and a rather large negative value outside the region D. To avoid this problem, we use the following regularization of Heaviside function .
A final remark on solving the under-determined inverse problems is as follows. Our original cost function to solve the under-determined CT reconstruction is the standard least-squares error of Eq. (2). In the setup of sparse-view CT and limited-angle CT, however, due to the lack of sufficient projection data, we normally need to include a regularization term into the cost function. In the proposed method, instead of using the regularization term, we have used the deterministic constraint that the image can be expressed by the simple model of Eq. (7). Since the dimension to represent the image in this model is essentially very small, it works similarly to the regularization. In this model-based approach, the detailed form of the cost function for image reconstruction cannot be explicitly written, unlike the case of using the regularization term.
5 Numerical Experiments
In this section, we describe experimental results for image reconstruction in sparse-view CT and limited-angle CT. The experiments were performed by using synthetic images and real X-ray CT projection data of a carved cheese slice . We compared the proposed method to some of the state-of-the-art reconstruction methods for binary tomography. In Section 5.1, we will discuss the details of the compared methods.
For the experiments with synthetic images, we used four images shown in Fig. 4. All these binary images consist of pixels. We considered standard parallel-beam geometry, in which every detector consists of pixels. For the case of sparse-view CT, the number of projection data uniformly sampled over the angular range was set to and Our aim is to decrease scanning time in sparse-view CT by reducing the number of angles. For the case of limited-angles CT, we considered four different angular coverages and . Finally, we evaluated the performance of the proposed method by adding different levels of Gaussian noise to the projection data, where the noise level was 0.1%, 5%, and 10%. Up to the 5% noise case, the proposed method could stably reconstruct a reasonable image.
5.1 Other Reconstruction Methods Used for Comparison
There are several reconstruction methods for tomography. We are focusing on the following four methods, which have been developed for binary tomography aiming at reconstructing nice images with limited projection data.
First: Total-variation (TV) method based on solving the following optimization problem.
denotes the regularization parameter, and D is the matrix that captures the horizontal and vertical discrete gradients. The above problem was solved using the Chambolle-Pock method in .
Second: Dual Problem (DP) method proposed by Kadu and Van Leeuwen , which is based on solving the following so-called Largrangian dual problem.
where denotes an optimal dual variable. This method is a convex formulation based on Lagrangian dual of the binary constrained least-squares problem, and they solved this problem by using L-BFGS Quasi-Newton method [43,44].
Third: Batenburg and Sijbers proposed the Discrete Algebraic Reconstruction Technique (DART) in . Although DART works well for a large class of binary images by exploiting the prior knowledge on gray levels in each region of the image, our proposed method is still better and can compare favorably with this method.
Fourth: The parametric level-set (PALS) representation proposed by Kadu et al. . In their work, they used -smooth Wendland’s function as the radial-basis function (RBF), and the corresponding PALS function is expressed as
where denotes each RBF, L denotes the number of RBFs, is the center to place RBF, and is a scaling parameter.
5.2 Sparse-View CT Case
In Fig. 5, we show reconstructed images for Example 1 when the number of projection data is 4. From the results, the proposed method provided an image that is very close to the ground truth compared to the other methods.
Fig. 6 shows the performance of our proposed method when the number of projection data is changed between and . We find that our method can reconstruct a reasonable image with as few as projection data.
In addition to the synthetic images, we applied the proposed method to real X-ray CT projection data of a carved cheese slice . In Fig. 7, we show the measured sinogram and the corresponding Filtered Back Projection (FBP) reconstruction from projection data measured over the angular range . The object can be observed to be approximately binary, consisting of only two gray levels, representing cheese region and air region. Now, we will show the performance of the proposed method with this data. Fig. 8 shows the comparison between the proposed method and the other methods. In Fig. 9, we show object boundaries (red contours) reconstructed from and projection data. We observe that the proposed method provided nice images with a small number of projection data.
5.3 Limited-Angle CT Case
In the case of limited-angle CT, the problem becomes much more difficult, but the proposed method still works well, and the reconstructed images were very close to the true images. For example, Fig. 10 shows the images by the proposed method and the other reconstruction methods for Example 3. In this experiment, the number of projection data was only and the angular range was set to .
In Fig. 11, we show how reconstructions by the proposed method vary by decreasing the angular range for Example 4. In CT fields, it is well-known for a long time that it is difficult to reconstruct images with sufficient accuracy from limited-angle projection data. However, as shown in Fig. 11, we could reconstruct almost correct images even when the angular range was . This is clearly thanks to the use of binary constraint.
For the real data, in Figs.12 and 13, we show reconstructed images by the proposed method when the angular range was limited to , where the number of projection data was and . We observe that it is still possible to correctly identify the existence of letters C and T.
5.4 Evaluation of Effect of Noise
We evaluated the performance of the proposed method in the presence of three different levels of noise in the projection data. We added Gaussian noise with ratio 0.1%, 5%, and 10% to the projection data, from which image reconstruction was performed. In Fig. 14, we show reconstructed images for Example 1 and Example 2. We observe that the proposed method is stable in the presence of modest amounts of noise but degrades gradually as the noise level increases.
In this paper, we proposed a novel approach using a parametric level-set method for image reconstruction in binary tomography. The basic formulation of the problem is summarized as follows. Shapes of objects are represented via a level-set function, which is represented by using a Gaussian radial basis function (GRBF). We are using an appropriate parametrization of the level-set function, which can reduce the problem dimension and essentially regularizes the problem in a successful way. Based on this modeling, the level-set function is more likely to behave well, so it doesn’t need re-initialization like the traditional methods. Furthermore, unlike the use of traditional level-set methods, the problem to be solved for image reconstruction becomes a tractable convex optimization. We solved the optimization problem using a Newton-type method because the number of underlying parameters in a parametric level-set approach is usually much smaller than the number of pixels that we finally would like to get the image by discretizing the level-set function. We tested the proposed method on synthetic images and real data. According to our experiments, the proposed approach compares favorably with some state-of-the-art reconstruction methods, and it is also relatively stable in the presence of modest amount of noise.
Acknowledgement: Real projection data used in this study is available from the members of the X-ray lab at the Faculty of Science, University of Helsinki, Finland upon request http://www.fips.fi/dataset.php. We greatly appreciate them to provide us with the data. We also thank four anonymous reviewers for their comments, which significantly improved the manuscript.
Funding Statement: This work was supported by JST-CREST Grant Number JPMJCR1765, Japan.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
|This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.|