An image reconstruction model regularized by edge-preserving diffusion and smoothing for limited-angle computed tomography

Limited-angle computed tomography is a very challenging problem in applications. Due to a high degree of ill-posedness, conventional reconstruction algorithms will introduce blurring along the directions perpendicular to the missing projection lines, as well as streak artifacts when applied on limited-angle data. Various models and algorithms have been proposed to improve the reconstruction quality by incorporating priors, among which the total variation, i.e. l1 norm of gradient, and l0 norm of the gradient are the most popular ones. These models and algorithms partially solve the blurring problem under certain situations. However, the fundamental difficulty remains. In this paper, we propose a reconstruction model for limited-angle computed tomography, which incorporates two regularization terms that play the role of edge-preserving diffusion and smoothing along the x-direction and y -direction respectively. Then, an alternating minimization algorithm is proposed to solve the model approximately. The proposed model is inspired by the theory of visible and invisible singularities of limited-angle data, developed by Quinto et al. By incorporating visible singularities as priors into an iterative procedure, the proposed algorithm could produce promising results and outperforms state-of-the-art algorithms for certain limited-angle computed tomography applications. Extensive experiments on both simulated data and real data are provided to validate our model and algorithm.


Introduction
As a non-destructive imaging technology, computed tomography (CT) has been widely applied in many fields like industrial inspection, medical diagnosis, etc. Conventional CT needs full angular scanning data for reconstruction. This requirement, however, cannot be satisfied in certain important applications. For example, in medical examinations, like breast tomosynthesis [34], C-arm neuro imaging [15], and image guided surgery [2], they all acquire limited-angle projection data. In micro-CT applications, when the object under scanning possesses a special shape, like long and thin plates, then the scanning might be confined to a certain angular range, or the acquired data for certain angles might be useless due to the heavy absorption of the x-rays, which also leads to limited-angle data.
Limited-angle CT is a highly ill-conditioned problem. It is well known that incomplete data leads to heavy artifacts in the images reconstructed by conventional algorithms, like the ART [20], SART [1] and FBP [16,33]. In fact, limited-angle CT has been studied as early as in the 1970s [33]. Since then, its ill-posedness and artifact characterization have been researched persistently. In [33,47], Smith and Natterer have discussed the uniqueness and stability of the limited-angle CT problem. By utilizing the tool of singular value decomposition, it was proved that the reconstruction should be highly unstable with small scanning angular ranges (e.g. less than 2π/3). In [31] and [36], it has been shown that in numerical inversion of the exterior Radon transform, boundaries tangent to the projection lines in the data set should be easy to reconstruct, while boundaries not tangent to projection lines in the data set should be blurred and difficult to reconstruct. This property was further investigated in [18,37], where the framework of micro-local analysis was utilized to characterize the behavior of conventional reconstruction algorithms when applied to the limited-angle CT. It was further shown that from the limited-angle data (principle 1 in [37]): (a) if a boundary is tangent to a line in the limited data set, then that boundary should be easy to reconstruct. Such boundaries are called visible boundaries (from this limited data); (b) if a boundary is not tangent to any line in the limited data set, then that boundary should be difficult to reconstruct. Such boundaries are called invisible boundaries (from this limited data). Besides, it is also verified that for the limited-angle problem (principle 2 in [37]), streak artifacts can occur along the ending lines of the limited angular range and tangent to the boundary of some features in the object.
With the development of theoretical analysis for the limited-angle CT, various modified algorithms that aim to reconstruct better images from incomplete data have been proposed. Early methods try to restore the complete projection data by some kind of 'extrapolation' on the limited-angle data, with smoothness prior or global properties of the projection data like the H-L consistency [32,53]. Recent methods try to incorporate various priors, including geometrical shape, distribution of edges or densities (gray values) about the image, into iterative reconstruction algorithms [35].
For certain medical and industrial applications, the CT images can be approximated well by piecewise constant functions. So image gradients should possess the sparsity property. Based on this prior, many total variation (TV)-based models and algorithms have been proposed.
The first work in this line can be found in [7], which is motivated by the compressed sensing theory and deals with parallel beams. The first practical model for fan-beam and cone-beam reconstruction within this category seems to be introduced in [46]. After that, modifications and improvements include the adaptive steepest descent-projection onto a convex sets method [45], prior image constrained compressed sensing method [10], soft-threshold filtering approach [22], and improved TV-based image reconstruction method [39], etc. More recent work utilizing the sparsity in the curvelet and shearlet tranform domain can be found in [17,38]. These methods, however, do not really take the properties of limited-angle CT into consideration, i.e. the TV does not encode any information about the configuration of the limited-angle scanning, or the visible and invisible features from the limited-angle data.
The anisotropic total variation (ATV) reconstruction approach, proposed in [26] and then refined in [11], explicitly considers the property of reconstructed images. The idea is based on the fact that when applying conventional iterative methods, e.g. ART or SART, on limited-angle data, its edge-recovery ability is angle-dependent. So the authors proposed to use the direction-aware ATV rather than the isotropic TV to serve as the regularization term. It was shown that the isotropic TV was not suitable for limited-angle reconstruction. This idea marks a breakthrough for limited-angle CT, which was further explored in [21] for exterior tomography.
Besides the image TV, i.e. the l 1 norm of the image gradient, the l 0 norm of the image gradient, which was proposed in [54], has also been applied to limited-angle reconstruction. In [56] and [55], it was demonstrated that the l 0 norm could perform better than the l 1 norm, because the l 0 norm is better at preserving edges. However, the properties of limited-angle data were not considered in [56] and [55]. Another interesting work is the reweighted anisotropic TV approach [51]. It is known that reweighted TV might approximate the l 0 norm of gradient during iterations. For the reweighted anisotropic TV, however, it is not clear where it converges to.
When the image could be approximated well by a piecewise constant function, its gray value levels should be sparsely distributed. So rather than using the sparsity of the image gradient, the DART method proposed in [4] makes use of the sparsity of the image gray levels, which achieves promising experimental results. A segmentation procedure needs to be incorporated in each iteration of the DART method, and prior information like the number of different gray levels and good approximations of them needs to be known in advance. Besides, the segmentation error during iterations might lead to unacceptable results. A recent improvement of the DART method can be found in [57], which does not need to know the gray levels in advance. A similar idea that utilizes the Mumford-Shah image segmentation model for limited data reconstruction with an application to electron tomography can be found in [28].
If the structure or shape information is available, then they could also be utilized to regularize the limited-angle reconstruction. With the assumption that the CAD file or other description file is at hand, Schorr et al proposed a reconstruction algorithm which took the edge distribution information provided by the CAD file as priors [41][42][43]. However, this kind of information is rarely available. So Liu et al [30] proposed to use optical scanning method to acquire the outer surface of the scanning object, which was then utilized as a constraint for the reconstructed image, i.e. the image values were set to zero outside the region defined by the surface. However, this method cannot remove the interior blurring and artifacts on the image.
Overall, theoretical analysis has provided deeper understanding of the limited-angle reconstruction problem. However, the theoretical results have not been taken to its full advantage to develop corresponding reconstruction algorithms. In this paper, based on the key property of limited-angle CT that when applying conventional reconstruction algorithms like the SART on limited-angle data, image edges should be rather accurately recovered along certain directions, while blurred along some other directions, we try to build an optimization model which performs anisotropic regularizations, i.e. edge-preserving diffusion along a certain direction, and edge-preserving smoothing along the corresponding perpendicular direction. This separated regularization strategy results in more flexibility in choosing regularizers and leads to better reconstructions over the state-of-the-art algorithms for certain limited-angle computed tomography applications.
The remainder of this paper is organized as follows. The proposed model and its solution algorithm are described in section 2, and section 3 is dedicated to experiments on simulated data as well as real data. Discussions clarifying several aspects of the proposed algorithm will be carried out in section 4. We conclude our paper in section 5.

Image projection and its discretization
For simplicity, we consider the two-dimensional limited-angle parallel-beam CT. It should be no problem to extend our model and algorithm to a fan-beam CT and 3D cone-beam CT. The scanning configuration for the limited-angle parallel-beam CT is illustrated in figure 1. In the Oxy coordinates system, let x = (x, y), u the location of the detector cells, and β the rotation angle, v (β) = (cos β, sin β), v ⊥ (β) = (− sin β, cos β), and f ( x) denote the linear attenuation distribution of the object to be scanned. Then parallel-beam projection data for f ( x) can be expressed as (1) (1) refers to the typical limited-angle CT problem.
The problem (1) can be discretized as follows. With respect to the field of view, f ( x) is firstly discretized as a two-dimensional digital image { f i,j |i = 1, 2, ..., K 1 , j = 1, 2, ..., K 2 }, where i and j denote the row and column indices respectively, while K 1 and K 2 denote the number of rows and columns, respectively. Let and then equation (1) can be approximated as where τ denotes the transposition operation, then the limited-angle imaging problem can be expressed as where the system matrix A = (a m,n ) is usually called the projection matrix or forward projection operator.

The proposed optimization model
In the following, we will construct our optimization model by utilizing the properties of the images reconstructed from limited-angle data. We found that the locations and the gray values of the visible edges in the reconstructed images are more reliable than those of other points. Figure 2 gives an example to show this phenomenon, where figure 2(a) is the Shepp-Logan phantom and figure 2(b) is the image reconstructed with the SART method from 90 degree data (β ∈ [π/4, 3π/4]). In figure 2(b), obvious blurring and distortion can be seen approximately along the y -direction while most edges approximately aligned to the y -direction are recovered faithfully, and streak artifacts appear along the directions of the ending lines of the angular range.
Since the SART method recovers the visible edges rather accurately, then one can fix them and diffuse their values to other areas to obtain a better approximation to the ideal image. Indeed, if the ideal image is piecewise constant or piecewise smooth, which is a reasonable assumption in real applications, then with a linear interpolation, good approximation can be obtained for non-edge points. Unfortunately, due to artifacts and noise, direct interpolation might be problematic. So we need a better strategy for edge-preserving diffusion. This motivates us to propose the following optimization model: where W is a diagonal matrix that defines a weighted norm [25], λ is a weighting parameter, and E h performs certain sparsifying transform along the x-direction. Note that the y -direction regularization is missing. This is because sparsifying transform along the y -direction is much less meaningful due to blurring and distortion. It needs to be pointed out that alternative norms other than the W-norm could have been specified in the above formulation. The W-norm is chosen because the SART method is utilized for reconstructions, which minimizes a weighted objective function defined by the W-norm [25].
If the ideal image is piecewise constant, then when approaching it, the y -direction should also demonstrate the property of the piecewise constant. Another issue is that the onedimensional regularization along the x-direction might introduce horizontal striping artifacts, since the x-direction regularization is done row by row, which might break the correlation between the image rows. So we need another regularizer to suppress such possible artifacts. This motivates us to add the y -direction regularization back and arrive at where E h and S v denote the horizontal x-direction transform and vertical y -direction transform operators, respectively, while φ and ψ define two metrics, such that the minimization of φ E h f and ψ S v f has the effect of edge-preserving diffusion and edge-preserving smoothing, respectively. The parameters λ 1 and λ 2 are weighting parameters to adjust the effect of φ E h f and ψ S v f in (8).
A typical choice for E h is the one-dimensional gradient operator, and φ(·) can be chosen as the l 1 or l 0 norm. To simplify the expressions, with a bit of abuse of notations, we would like to denote the l 1 norm of gradient as TV l1 , and l 0 norm of gradient as TV l0 , respectively. It turns out that TV l0 might do a better job, because it behaves better than TV l1 at preserving edges. If one also chooses TV l0 regularization along the y -direction, then we have where sgn denotes the sign function with the property sgn(0) = 0. Besides, we define f 0,j f K1,j , f i,0 f i,K2 to enforce the periodic boundary condition which is needed due to fast Fourier transform (FFT) shall be utilized to solve the Subproblem-2 and Subproblem-3 described below when the l 0 norm is involved.
It is necessary to point out that ψ S v f aims at edge-preserving smoothing, rather than diffusion. One can choose other regularizers for the y -direction smoothing like TV l1 , and even the median filter could do the job well, which will be demonstrated in the experiments section. On the other hand, the choice for φ E h f is vital for successful reconstructions. While alternative choices are possible, in this paper, we found that choosing φ E h f = TV l0 worked the best and produced much better images than TV l1 .

The alternating edge-preserving diffusion and smoothing algorithm
The optimization problem (8) might be non-convex, and one can only expect local minimizers. So the solution strategy might put a heavy influence on the computed solution. In our practice, we adopt a simple alternating minimization technique to solve (8). Let f (k) denote the solution after the kth iterations, and λ 0 be a positive constant, then the updating for f (k+1) can be split into three subproblems: Subproblem-1: Subproblem-2: Subproblem-3: The above splitting strategy is motivated physically, rather than mathematically. Denote the ideal image as u , which is assumed to be piecewise constant. Starting from some initial guess u 0 , our algorithm will first perform image reconstruction to get u 1/3 . According to the theory of visible and invisible boundaries, some of the edges will be faithfully reconstructed in u 1/3 . Then a piecewise constant approximation to u 1/3 is constructed, based on the visible edges. This is done by solving the Subproblem-2 to get u 2/3 which is a more piecewise constant along the rows, followed by solving the Subproblem-3 to get u 1 , which is a more piecewise constant along the columns. To respect the fact the gradient of u along the x-direction is much more faithful than its gradient along the y -direction, the weighting parameters λ 1 and λ 2 should be set properly, and the best ratio for λ 1 /λ 2 is problem dependent. Now u 1 is a piecewise constant approximation to u 0 , and one can expect that u 1 is also more close to u than u 0 , since u is piecewise constant. Then next iteration starts, which should compute a solution even closer to u , and so on.
Even with such a physically meaningful explanation, the proposed algorithm, however, lacks a strict mathematical justification. We cannot find an existing theoretical framework to accommodate it. This will be investigated further in the next subsection.
The Subproblem-1 can be solved approximately by applying one iteration (all the rays are traversed once) of the SART method with f (k) as the initial guess. In our tests, by employing the random ordering technique, this strategy works well. Note that traditional SART is to minimize the second term of the objective function (11). The rationale behind the application of SART for approximating the solution of Subproblem-1 can be justified as follows. Let the total number of views be V, then the second term can be written as where A i denotes the projection operator (matrix) for the ith view of scanning, p i denotes the projection data (vector) for the ith view, and W i denotes the diagonal matrix whose diagonal elements are assembled from the rows of A i . By introducing A 0 = I f as the identity matrix that maps f to itself and p 0 = f (k) , which could be thought of being the projection operator and the projection data corresponding to an additional view. Let W 0 = 1 √ λ0 I , then the minimization problem (11) can be written as a weighted least square one which perfectly fits to the SART method for pursuing the solution. By employing the random ordering technique, just one iteration of SART could compute good approximations. If one fixes A 0 as the first view, one iteration of SART corresponds to using f (k) as the initial guess, and then applying one iteration of SART on the remaining projection data.
In practice, the SART method is usually combined with a hard-thresholding procedure to impose the non-negativity constraint (prior). In this case, the objective function should be changed accordingly. We deal with this situation by changing the objective function of Subproblem-1 (so also the original objective function) to where and R + denotes the set of all vectors (the same dimension as f ) with non-negative components. An approximate solution of the above equation can be obtained in two steps: the first step is the one iteration of SART method described above, which gives an intermediate result denoted as f , then f (k+1/3) is obtained by thresholding f . One can also think that the original objective function consists of four terms, and the updating of f (k+1) is decomposed into four subproblems. Please refer to algorithm 1 for details. The Subproblem-2 can actually be decoupled into a series of one-dimensional problems, which are then solved by utilizing the algorithm proposed in [54], with a little extra work to tailor its two-dimensional code for solving one-dimensional problems. It can be briefly described as follows. Firstly, introducing an auxiliary variable g = {g i,j |i = 1, 2, ..., K 1 , j = 1, 2, ..., K 2 }, and let g n = g i,j , where n = (i − 1) K 2 + j. Then the Subproblem-2 can be transformed to , λ is used to control the similarity between D h f and g. Then the alternating minimization method is employed to minimize f and g separately in an iterative manner. Just as described in [54], λ is increased by multiplying a constant κ > 1 during the iterations, until λ λ max , where the iteration stops. In our experiments, we always set λ max = 10 4 and κ = 5.
The solution of Subproblem-3 depends on the choice of the regularization term. If TV l0 is chosen, then its solution algorithm is similar to Subproblem-2, and if TV l1 is chosen, i.e.
ψ S v f = D v f 1 , then the split-Bregman algorithm [19] or Chambolle's algorithm [8] can be employed. In our experiments, the split-Bregman algorithm is adopted. Using a median filter to play the role of regularization needs more explanation. We know that the solution of Subproblem-3 f (k+1) should be a smoother version of f (k+2/3) with the property of edge preserving. Intuitively, this functionality could be fulfilled by other edge-preserving filters. In fact, this idea has been presented in the literature [50] as early as in 2013, or might be even earlier. Recent studies have extended this idea to many other areas, bearing the name 'plug-andplay'. For example, subproblems evolved in applying ADMM [5] could be replaced by more popular off-the-shelf denoisers. According to the reports, by playing with the 'plug-and-play' technique, better results could be obtained in spite of its ad-hoc nature. The problem with the 'plug-and-play' is that theoretical analysis is very difficult to carry out. Indeed, replacing a energy minimization problem with a smoothing filter invalidates the optimization interpretation of the whole algorithm immediately. Still, some efforts dedicated to the convergence analysis by assuming certain properties of the denoising operators can be found in [9,40,49]. The 'plug-and-play' technique provides a way to extend the functionality or improve the performance of existing algorithms. Replacing the Subproblem-3 with a median filter is exactly in the spirit of 'plug-and-play'. The radius of the template used in the median filter could play the role of the weighting parameter λ 2 . Of course, other edge-preserving smoothing filters also deserve a try. A median filter is tested because of its simplicity and efficiency.
We would like to name the above algorithm as alternating edge-preserving diffusion and smoothing (AEDS). It should be mentioned that after submitting this paper, several other members of our research group also dived into limited-angle reconstruction. They tried the regularization of zero norm of second-order derivatives to see if it could loose the piecewise constant assumption. The algorithm proposed in this paper was borrowed to develop a similar algorithm to solve their optimization model. The result was then presented in the fifth international conference on image formation in x-ray computed tomography.
To clarify different diffusion and smoothing operators, we parameterize AEDS with appropriate arguments, e.g. AEDS(l 0 , l 1 ) means that TV l0 and TV l1 are employed as the diffusion operator and smoothing operator respectively. If an argument is not given, it means that the corresponding regularization is not performed, e.g. AEDS(l 0 ) means that only TV l0 along x-direction is performed. For the median filtering, argument mf will be adopted.
The whole algorithm can be summarized as algorithm 1 described below. Note that the function symbol SART f , p, ω denotes the operation of applying one iteration of SART method using f as the initial guess, p as the projection data, and ω as the relaxation parameter. Another function symbol shrink is defined as When x is a vector, then shrink performs the thresholding operation in an element-wise way.
Step 1: Step 2: Solving optimization problem: Solving optimization problem: Using the split-Bregman algorithm to solve the optimization problem: Filtering f (k+2/3) to arrive at f (k+1) by applying the median filter along the vertical direction, with a problem-dependent template size. end if end for

Convergence behavior of the AEDS algorithm
The proposed algorithm AEDS aims to solve the minimization problem (8). Since the l 0 norm of gradient is always utilized, one can only expect local minimizers from the AEDS algorithm. Unfortunately, we cannot provide a theoretical characterization of the behavior of the AEDS algorithm at the current stage. The difficulties come from two aspects. One is that the functional consists of three terms, while most of the theoretical analysis for alternating optimization algorithms is done for two terms splitting. The other one is that the l 0 norm of gradient lacks proper regularity properties, e.g. smoothness and convexity, which are usually required for carrying out theoretical analysis.
Existing approaches that could be utilized to solve problem (8) include ADMM [5], primaldual [46] and proximal splitting methods [13]. The basic form of these methods handles only two terms splitting. There are extensions to deal with multiple term splittings, but very few convergence results have been established. In [48], ADMM is utilized to minimize a similar functional like AEDS(l 0 , l 0 ) does, and sequence convergence is proved. However, this convergence result is established at the price of not updating some of the Lagrangian multipliers. Besides, it says nothing about whether or not the sequence converges to local minimizers. By extending 'projection' to 'proximal operators', the proximal splitting methods [13] have been shown to be more general to take ADMM, projected gradient method [29], Douglas-Rachford splitting [14] etc as special cases. However, it still focuses on two terms splitting. Its extension to multiple terms splitting with guaranteed convergence is actually pulled back to the two terms case by reformulation, which leads to the parallel proximal algorithm [12].
The proposed AEDS algorithm can find its roots in proximal operators. In a sense, AEDS is an extension to the backward-backward algorithm [3], which can be applied to functionals consisting of only non-smooth terms. Let Γ 0 (R N ) denote lower semicontinuous convex functions mapping R N to (−∞, ∞), then the proximal operator Prox g (·) is defined as follows: admits a unique solution, which is denoted by Prox g (x). Suppose that we are trying to find a solution for the problem where ζ > 0 is a weighting constant. Starting from some initial point x 0 , the backward-backward algorithm performs the following iterations: for n = 0, 1, . . . , It has been shown in [3] that, under certain conditions, the sequence (x n , y n ) converges to the minimizer of (16). One could think that the backward-backward algorithm is to find a minimizer for the problem To make (17) easier to solve, we introduce another variable y = x, and then use the penalty method to deal with the constraint y = x. So problem (16) can be thought of being a relaxation to (17). Then let us consider the following problem: To easier its solving, one could also introduce variables y and z, and requires that y = x, z = x. Then, given x 0 , a direct extension of the backward-backward algorithm gives the following iteration: Here, the weighting parameters ζ i , i = 1, 2, 3 can be thought of playing the similar role as ζ in (16). If one substitutes u(x), v(x) and w(x) with the corresponding fidelity term and regularization terms in (8), the AEDS algorithm is then immediately identified. It should be pointed out that, at each iteration step, the AEDS algorithm computes only approximate solutions (e.g. local minimizers rather than the required global ones) for the proximal operators, which could be thought of being another extension to the backward-backward algorithm.
Unfortunately, at the current stage, we cannot provide a convergence analysis for the iteration (19). In fact, a rather thorough investigation of the literature reveals that no existing theory or framework could handle the proposed algorithm. So in the next section, numerical experiments will be instead carried out to demonstrate the convergence behavior of AEDS.

Experiments
In this section, we will first try to validate the generality of the proposed model (8) and the AEDS algorithm. Experiments on simulated discrete and analytic phantoms will be carried out to demonstrate the behaviors of AEDS, such as sensitivity to noise and stability against model perturbations. Then experiments on simulated as well as real data will be carried out to study the effectiveness of AEDS against state-of-the-art algorithms like the TV l1 + SART, TV l0 + SART [56], DART [4] and the Potts model [48]. For brevity, the abovementioned four competing algorithms shall be named TV l1 , DART, TV l0 , and Potts model, respectively. It is easy to recognize their meanings referring to either regularizer or algorithm from the context.
For simulated discrete phantoms, the projections are obtained by applying a fast version of the Siddon's method [44], i.e. the incremental ray casting method [24]. For analytic phantoms, the projections are obtained by utilizing the open source computed tomography simulator-CTSim (www.ctsim.org). Since the forward projection calculation in CTSim is analytic, the inverse crime [27] should have been removed. For real data, even though widely used, the linear imaging model (6) is not accurate enough to model the real forward projection process. So experiments on real data would help to further demonstrate AEDS for their tolerance to imaging model perturbations in real applications.
In all the experiments, if noisy data are needed, Poisson noise specified by the incident intensity, denoted as I 0 , would always be added to the raw data, i.e.
where the symbols log and poissrnd denote the Matlab functions for logarithm transform and Poisson random numbers generation, respectively, while p and p noisy denote the noise-free and noisy projection data, respectively.
Note that the scanning configuration illustrated in figure 1 is adopted in all the experiments with simulated data. So, roughly speaking, vertical edges are visible, and horizontal edges are invisible. The geometrical parameters listed in table 1 is used to acquire the raw data, if not otherwise specified. To start the proposed algorithm, SART with 10 iterations is applied on the limited data to compute an initial guess. This is done in all the experiments.
A C++ implementation of the proposed AEDS algorithm and the competing algorithms (except the Potts model) can be downloaded from https://github.com/jinqiuXu/ AEDS-Limited-angle-reconstruction-

Validation of the proposed AEDS algorithm
In this subsection, experiments will be carried out to demonstrate various aspects of the proposed model and the AEDS algorithm. Mechanism verification will be first performed, which aims to provide numerical evidence on why and how AEDS works, e.g. utilizing visible edges as prior to recover invisible edges. Then sequence convergence will be examined by studying the energy curves and absolute changes between successive iterations. Inverse crime tests will be then performed to clarify whether or not inverse crime [27] is a big deal for AEDS, which helps to demonstrate AEDS's stability against model perturbations. Finally, AEDS will be further tested on specially designed phantoms, which tries to clarify when and where AEDS works.
3.1.1. Mechanism verification. Two phantoms are used to test how the proposed model (8) and the AEDS algorithm perform. The first phantom is the classical Shepp-Logan phantom, as shown in figure 2(a), while the other one, as shown in figure 3, is a rectangular phantom with several elliptic regions, disks, squares and polygon regions. Generally speaking, the rectangular phantom is more challenging for limited-angle imaging, due to its high length/width aspect ratio. Note that since the discrete phantoms and their reconstructed images are of the same size, and the same forward projection operator is utilized for both data simulation and reconstruction, so the so-called inverse crime might exist.     Figure 5 shows the reconstructions of AEDS(l 0 ), AEDS(l 0 , mf ), AEDS(l 0 , l 1 ) and AEDS(l 0 , l 0 ), respectively. The parameters used by the algorithms are listed in table 2, where ω refers to the relaxation parameter of SART, r denotes the diameter of the template used in the median filtering and µ balances the fidelity and regularization involved in the split-Bregman algorithm [19] for solving the TV l1 problem. For the AEDS(l 0 ) algorithm, most distortion and blurring introduced by the SART method (as shown in figure 4(b)) have been removed, which indicates that the TV l0 regularization along the x-direction would help to recover most of the features of the true image. However, some horizontal striping artifacts are observable. We think that these artifacts are introduced because the regularization along different image rows is completely uncoupled, and no respect is given to the correlation or smoothness property between the image rows. So, regularization along the y -direction is needed to suppress such artifacts. Note that edge-preserving property is also needed for the y -direction regularization, otherwise it will smooth out reconstructed true edges. Three regularizers for the y -direction smoothing have been tested, i.e. the median filtering, TV l1 , and TV l0 , and the reconstructed images are shown in the first row of figure 5. One can see that both the TV l0 and TV l1 regularizers can effectively remove the striping artifacts and almost perfect reconstructions have been obtained. To look into more details, the residual images, which are defined as the absolute differences between the reconstructed images and the true images, are also computed and shown in the second row of figure 5. Clearly, without y -direction regularization, large errors will be introduced, as shown by the residual image of AEDS(l 0 ). With edge-preserving smoothing along the y -direction incorporated, high-quality reconstructions are obtained, as shown by the residual images of AEDS(l 0 , l 1 ) and AEDS(l 0 , l 0 ). Another observation is that, the residual image of AEDS(l 0 , l 0 ) is quite close to zero, which indicates that the result of AEDS(l 0 , l 0 ) is of the highest quality among the four algorithms.
To compare the results quantitatively, we have computed the peak signal to noise ratio (PSNR) [23], mean square error (MSE), and structural similarity (SSIM) [52] indices of the reconstructed images, which are shown in table 3. To compute the above indices, reference images are always needed. The three algorithms with y -direction regularization achieve higher PSNR than AEDS(l 0 ). In addition, the MSE and SSIM indices are also consistently better, which agree with the residual images. According to the indices, AEDS (l 0 , l 0 ) performs better than AEDS (l 0 , l 1 ) in this test.

Experiment 2.
Our second experiment is performed on the rectangular phantom size of 512 × 150, with scanning angular range [π/6, 5π/6]. Poisson noise with incident intensity I 0 = 1.5 × 10 6 is added to the scanning data. Figures 6(a) and (b) show the full-angle and limited-angle reconstructions, respectively, by applying 10 iterations of SART. The reconstructed images are size of 512 × 512. However, for a better display, they are clipped to their original size 512 × 150. Figure 7 shows the reconstruction results of the proposed algorithms, with the reconstruction parameters listed in table 4. It can be seen that all the image edges are reconstructed correctly with the proposed AEDS(l 0 , * ) algorithms. A similar conclusion can be drawn that regularization along the x-direction helps to recover the features of the true image, while regularization along the y -direction helps to suppress the underlying horizontal striping artifacts.
The computed quantitative measures, PSNR, MSE, and SSIM indices, are shown in table 5. Among the four AEDS(l 0 , * ) algorithms, AEDS (l 0 , l 1 ) performs the best. This is different from the previous experiment.
From the experimental results, one can see that AEDS (l 0 , l 1 ) performs a bit better for the more challenging rectangular phantom. Considering space limitations, in the following experiments, we show only the results of either AEDS (l 0 , l 0 ) or AEDS (l 0 , l 1 ).

Convergence test of AEDS.
We choose AEDS (l 0 , l 0 ) as the representative to explore the convergence behavior of the proposed AEDS algorithms. In fact, AEDS (l 0 , l 1 ) is also tested, and its convergence behavior is quite similar. The energy functional is computed to check how it evolves with the iterations. Besides the energies, the absolute increments between successive iterations are also computed to check how the computed image sequences behave with the iterations. To this end, we define d (k) = f k − f (k−1) 2 2 , which will be computed and illustrated. The energy functional corresponding to AEDS (l 0 , l 0 ) reads To ease the computations, the W-norm has been replaced by l 2 norm in computing the energies. This replacement should not lead to contrary conclusions since the two norms are equivalent. Besides, direct computation of D h f 0 and D v f 0 are not practically possible due to numerical errors. Instead, they are approximated by where the nth component of g h and g v is computed as where λ max denotes a big positive number, which has been set to 10 4 in all our experiments. The data used for plotting are obtained from experiment 1 performed in section 3.1.1. The energy and the increment curves are plotted in figure 8. Despite of some local small oscillations, all the curves keep decreasing with the iteration numbers. After 800 iterations, both the energy and increment curves reach stable states.
The global decreasing property of the energy curves suggest that the image sequence f k produced from AEDS (l 0 , l 0 ) is to minimize the energy functional. The decreasing behavior of d(k), even though does not imply convergence, indicates that the sequence is approaching a steady state, especially considering the fact that the oscillations diminishes with the iterations.

Inverse crime test of AEDS.
Projection simulation from discrete phantoms, especially when the phantoms and the reconstructed images are of the same size, may lead to an inverse crime which helps to achieve over-positive results. To avoid an inverse crime, we employ CTSim to acquire the projection data analytically, i.e. the forward projection is calculated analytically on an analytic phantom. To provide quantitative measures, reference image is needed. This is done by applying the SART method (10 iterations) on the full data acquired from the analytic phantom.
The algorithm variant AEDS (l 0 , l 1 ) is chosen as a representative of AEDS algorithms to perform the tests against the competing algorithms: DART, TV l1 , TV l0 and the Potts model. To further alleviate the concern of the inverse crime, both noise-free and noisy projection data are simulated to test whether or not inverse crime could be identified.
In this test, the projection data is acquired with an analytic rectangular phantom constructed in CTSim. The scanning angular range is 100 degrees, i.e. [2π/9, 7π/9]. Poisson noise with I 0 = 1.5 × 10 5 is added to the scanning data, which is also much higher than that with previous experiment.
The full-angle reconstructions (SART, 10 iterations) for noise-free and noisy data are shown in figures 9(a) and (b), respectively, and the limited-angle reconstructions (SART, 10 iterations) for noise-free and noisy data are shown in figures 9(c) and (d), respectively. Similar to the previous experiment, the full-angle reconstructions are degraded by streak artifacts and noise. In the limited-angle reconstructions, all the edges which are supposed to be invisible are completely blurred and can not be identified, i.e. the invisible edges are completely missing.
AEDS(l 0 , l 1 ) is chosen to compare with the competing algorithms, and as stated before, the reconstructions shown in figures 9(c) and (d) serve as the initial guesses for reconstructions with noise-free data and noisy data, respectively. The parameters used by AEDS(l 0 , l 1 ) and the competing algorithms are listed in table 6. Figure 10 shows the reconstructed images of the analytic phantom. Clearly, the reconstructions of DART, TV l1 and TV l0 suffer from severely distortions, even for noise-free data. For the Potts model, local distortions are still clearly seen at the right bottom part, and noise leads to more distortions. When checking the results of   AEDS(l 0 , l 1 ) shown in the last column of figure 10, two observations could be made. One is that the reconstruction quality is much higher than the competing algorithms, e.g. no obvious distortions could be observed. Another one is that in terms of distortions the reconstruction is less affected by noise.
When checking the quality measures of PSNR, MSE and SSIM, as listed in tables 7 and 8 for the noise-free and noisy data, respectively, one can see that the indices agree with the reconstructions well. Potts model and the proposed algorithm perform the best, while AEDS(l 0 , l 1 ) wins the first position with a clear margin.

Limitation test of AEDS.
Experiments will be performed to test AEDS's ability to recover invisible edges when the visible edges are not distributed in a favorable way. One rhombus phantom with tilt angle of 5 degrees is constructed in CTSim to acquire the projection data. Both the phantom and the forward projection are analytic, and no noise is added to the data. This is a very challenging problem since the inverse crime has been removed, and  most importantly, the horizontal diffusion does not match the invisible edge due to the nonzero tilt angle. The reconstructed image is size of 512 × 512. AEDS(l 0 , l 1 ) is chosen as the testing algorithm. Figure 11(a) shows the rasterized image (with resolution 2048 × 2048) of the designed phantom by using the CTSim software, while the image reconstructed by applying the SART method (10 iterations) on the full data is shown in figure 11(b). The parameters adopted by AEDS(l 0 , l 1 ) are listed in table 9.
The reconstruction results are shown in figure 12. As shown in the first row, the SART reconstructions with angular ranges from 120 degrees to 150 degrees fail to restore the invisible edge. For the AEDS(l 0 , l 1 ) algorithm, the reconstruction quality demonstrates dependence on the scanning angular range. For the 120 degree case, the invisible edge is not approximated well. It seems that the algorithm tries to match it with a horizontal edge. When the scanning angular range is increased to 130 degrees, however, the invisible edge has been successfully recovered. If one continues to increase the scanning angular range, the reconstruction quality can be improved, though slowly. Careful examination on the results of AEDS(l 0 , l 1 ) reveals that there exists one horizontal strip artifact, i.e. the reconstructed gray values are a bit higher than expected, just above the recovered invisible edge, though quite weak. This artifact can be attributed to the horizontal diffusion process of the AEDS algorithm.

Performance test of AEDS against state-of-the-art algorithms
In this subsection, we will test how the proposed model and algorithm perform under various situations against state-of-the-art algorithms like the TV l1 + SART, DART, TV l0 + SART and the Potts model, so as to verify their effectiveness and robustness. Experiments will be done on simulated data as well as real data. As explained before, the abovementioned four competing algorithms will be named TV l1 , DART, TV l0 , and Potts model, respectively.
The TV l0 method regularizes the conventional two-dimensional l 0 norm of gradients, while both the Potts model and the proposed AEDS(l 0 , l 0 ) regularize the images with the l 0 norm of directional derivatives. The difference between them lies in their weighting strategies for the directional derivatives of the images. By assigning weights to different directions, the l 0 terms in the Potts model aim to approximate the Euclidean length. The proposed AEDS(l 0 , l 0 ) considers only the directional derivatives along the x-direction and the y -direction, and the weights are assigned according to the properties of the image reconstructed by SART from limited-angle data. The best weights would be problem-dependent. For all the algorithms, SART is always applied for reconstruction. The only parameter involved with SART is its relaxation parameter ω. For the DART method, parameter ρ defines the image gray value levels needed by the segmentation procedure during each DART iteration. For the TV l1 algorithm, just one parameter µ, which balances the fidelity and the regularization and keeps the same meaning as in [19] needs to be tuned. For the TV l0 method, the parameters λ l0 corresponds to the parameters λ defined in Yu's algorithm [56]. For the Potts model, the code from https://github.com/mstorath/Pottslab is utilized for its implementation, and the parameters involved are the selected directions and the weight γ for balancing the regularization term, which are defined in [48]. The directions are actually used by the Potts model will be denoted by vectors like (1, 0), (0, 1) and (1, 1), etc, while the parameter γ is problem dependent.
For the TV l1 , DART, TV l0 , AEDS(l 0 , l 0 ) and AEDS(l 0 , l 1 ) algorithms, 1000 iterations are performed for simulated data, since they always demonstrate convergence before 1000 iterations. The Potts model algorithm usually needs more iterations for convergence. The iteration number shall be specified with each experiment. For real data, the four competing algorithms, i.e. TV l1 , DART, TV l0 and AEDS(l 0 , l 0 ), would run 2000 iterations to guarantee convergence.
The initial guess f (0) for all tested algorithms is set to be the reconstruction result by applying the SART algorithm with zero initial guess and 10 iteration cycles, on the limited-angle data.

Shepp-Logan phantom.
Experiments will be done first on the Shepp-Logan phantom, which makes it easier to compare the performance of different algorithms. Two sets of projection data, i.e. 120 degree noise free and noisy data, are tested. The noise level is set to I 0 = 1.0 × 10 7 . Figures 4(a) and (b) provide reference images for full-angle and limited-angle reconstructions, respectively. Streak artifacts and edge distortions along vertical directions, especially at the top and bottom parts, are clearly seen. On the other hand, visible edges are abundant and well distributed, which indicates that the reconstruction should be relatively easy. Various parameters for the abovementioned five algorithms with different configurations are listed in table 10. The parameters have subjected to fine-tuning in terms of PSNR. Figure 13 shows the reconstruction results of the proposed algorithm and the four competing algorithms. It can be seen that both DART and TV l1 could effectively suppress the streak artifacts, but fail to remove the image blurring at the top and the bottom parts of the phantom. In addition, the blurring becomes heavier with the noisy data. Besides, DART appears to be more sensitive to noise. This is understandable since the DART method considers only the gray value levels of the image, and the local spacial smoothness prior is ignored. For this simple phantom, all the TV l0 , the Potts model, and the proposed method produce high quality images. The image edges and structures have been well recovered, even from noisy projections.
When we check the quantitative measures shown in tables 11 and 12, however, it turns out that the Potts model performs quite well for both noise-free and noisy data tests, and the performance of AEDS(l 0 , l 0 ) is comparable with the Potts model.

The rectangular phantom.
To further demonstrate the behavior of the proposed algorithm, tests on the more challenging rectangular phantom will be performed against the competing algorithms. Two sets of projection data, i.e. 120 degree noise free and noisy data, are tested. The noise level is increased to I 0 = 5 × 10 5 . The full-angle reconstructions (by SART, 10 iterations) for noisy data is shown in figure 14(a), while the corresponding noisy limitedangle reconstruction is shown in figure 14(b). AEDS(l 0 , l 1 ) is chosen to compare with other algorithms, since experiments show that it performs better than AEDS(l 0 , l 0 ) for dealing with high noise levels. The parameters for the reconstruction algorithms are listed in table 13. Figure 15 shows the reconstruction results of the rectangular phantom. We can see that DART, TV l1 and TV l0 introduce various artifacts and the reconstructed images are severely distorted. The Potts model restores some of the structures of the phantom, but failed to reconstruct the horizontal edges of the phantom even for noise-free data. On the other hand, the proposed AEDS(l 0 , l 1 ) algorithm produces nearly perfect reconstruction. TV l1 and TV l0 are not aware of the properties of limited-angle CT, and the x-direction regularization is tangled with the y -direction regularization. We think that this entanglement is the key reason that accounts for the bad results. The Potts model applies different weights to different directional derivatives. But the purpose is to approximate the Euclidean length instead of encoding the prior of the reconstructed image or the properties of the scanning configuration. While our model, by assigning proper weights to different directional derivatives, has taken the properties of the scanning configuration into consideration, it thus gets the best reconstruction results.
The above analysis and conclusion can be further validated by comparing the quantitative indices like PSNR, MSE, and SSIM, which are computed and listed in tables 14 and 15.

Experiments on real data.
To test the potential capability of our method in practical applications, experiments are carried out on a real flat object, as shown in figure 16. This is a typical application of limited-angle tomography, and the reconstruction is very challenging, since there are very few visible edges. The difficulty also comes from the large ratio between length and width. Due to the long and thin structure, the edge information needs to travel from one end to the other end to successfully recover the invisible edges. To make this process easy,       For this test, the complete projection data was firstly acquired by full-angular scanning. Then the SART algorithm was applied on this complete data to construct a reference image. The limited-angle data were extracted from the complete data. The angular range is [π/6, 5π/6]. The reconstructed images are size of 512 × 512. However, for better display, they are clipped to the size of 512 × 150. The full angle reconstruction (SART, 10 iterations) is shown in figure 17(a). The SART method (10 iterations) is also applied on the limited-angle data, and the reconstruction is shown in figure 17(b).  The parameters for the considered algorithms are listed in table 17. Reconstruction results corresponding to the solutions after 2000 iterations for the TV l1 , DART, TV l0 , and the proposed AEDS(l 0 , l 0 ) algorithms are shown in figure 18. All the competing algorithms reconstruct more or less correctly the visible edges along the x-direction. However, only the AEDS(l 0 , l 0 ) algorithm could recover most of the missing edges well. Even some artifacts are still noticeable, the global structures have been recovered well. The PSNR, MSE and SSIM indices computed and listed in table 18 further demonstrate the superiority of AEDS(l 0 , l 0 ) over the competing ones. As pointed out earlier, we have employed an open source implementation of the Potts model regularization, which works only on parallel-beam data. So the Potts model is not tested on this real data.
The PSNR, MSE and SSIM indices are computed and listed in table 18.

Discussions
From the experiments section, we have shown that the proposed AEDS algorithm is robust against noise and model perturbations. One can tell that the object shape has also substantial   effect on the reconstructions. The long plate shape can reduce significantly the effectiveness of the TV l1 , DART and TV l0 algorithms. On the other hand, the proposed AEDS(l 0 , l 0 ) is still robust against the object shape changing, and high quality reconstruction can still be obtained. Experiments show that AEDS(l 0 , l 1 ) could perform better than AEDS(l 0 , l 0 ) for very noisy data. This might be because that the l 1 norm is better in suppressing noise, though weaker than the l 0 norm in preserving edges. In this section, we want to further clarify several aspects of the proposed model and algorithm. DART TV l1 TV l0 AEDS(l 0 , l 0 )

Comparison with the Potts model
In [48], the authors proposed a model which seems to incorporate similar regularization terms as our model (8). They started from the variational formulation of the Potts model arg min which was then discretized as where f and p are the corresponding discretized vector of f and p , ∇ ps denotes finite difference along some direction p s , and S is the number of directions considered. Then the minimization problem (25) was relaxed by introducing additional variables and augmented Lagrangian was employed to enforce the constraint. The relaxed problem was then solved by ADMM. Although the algorithm is similar to the proposed AEDS(l 0 , l 0 ) if one considers only the x and y direction regularizations, they are essentially quite different.
Firstly, the authors aim to solve (24) rather than (25), i.e. problem (25) is just an approximation to problem (24) on discrete grids. Note that (24) corresponds exactly to the TV l0 algorithm. So the weighting parameters ω s were not meant to respect the properties of the limited-angle data, they were chosen to approximate the continuous counterpart of ∇f 0 . In a sense, the Potts model was utilized to deal with sparse angle imaging problem, rather than the limited-angle CT. For our AEDS(l 0 , l 0 ) algorithm, the parameters λ 1 and λ 2 play different roles. λ 1 is chosen to show respect to the visible edges of limited-angle projection data, while the λ 2 is chosen to smooth out possible stripping artifacts as well as to promote piecewise constantness along the y -direction.
Secondly, the proposed AEDS algorithms are innovative, i.e. different from the popular ADMM which has been used in [48]. Since the objective functions are non-convex, different algorithms might result in quite different reconstructions. The proposed AEDS algorithm possesses a favorable physical background which might be helpful in explaining its behaviors, and its superior performance in practical applications.

The piecewise constant assumption
All the algorithms addressed so far place an underlying assumption on the true images to be reconstructed: the ideal image should be piecewise constant. This is an implicit assumption for both the TV l1 and TV l0 norms. The DART method, however, requires this assumption indirectly, i.e. it requires that the gray value distribution concentrates on several levels. This might lead to vulnerability to noise, which has been demonstrated in the experiments section.
For real projection data, due to various reasons like beam hardening [6], the solution of the discrete linear imaging model (6) cannot be regarded as a piecewise constant function, and piecewise smooth might be a more appropriate assumption. In this case, the competing algorithms fail to reconstruct quality images. On the other hand, the proposed AEDS(l 0 , l 0 ) algorithm still faithfully recovers high-quality approximation to the true image. From figure 18 we can see that the recovered image is very close to a piecewise constant function, which approximates the reference image shown in figure 17(a).

Parameter selection
For the AEDS algorithm to work properly, the strength of diffusion and smoothing has to be balanced well. This raises the problem of parameter tuning. Take the AEDS(l 0 , l 0 ) for example, the two parameters λ 1 and λ 2 have to be properly set, which are actually problemdependent. A rule of thumb is that diffusion should be stronger than smoothing, so usually λ 1 > λ 2 . When the horizontal and vertical regularizers are of different type, however, we do not have a guideline at hand. Other parameters related to the SART method might also contribute a lot to the reconstruction quality. The information provided by the regularization terms (diffusion and smoothing) and the information provided by the SART method (from the projection data) must be 'blended' with a reasonable ratio. This is also a concern for all the competing algorithms. In our experiments, for relative fair comparison, all the algorithms have subjected to parameter tuning by the technique of error and trial to select the best ones according to PSNR.

Performance versus distribution of visible edges
By the virtue of the proposed algorithm, the x-direction diffusion is completely decoupled from the y -direction smoothing. By solving Subproblem-2, the information of vertical edges is propagated along the x-direction. This is just the behavior that being devised to utilize the visible edges. However, there might be a side-effect of this strategy: the algorithm favors perfect horizontal edge restoration (vertical edges are visible). This phenomenon can be observed in the results of the limitation test shown in figure 12 of section 3.1.4, where, for the 120 degrees case, the proposed algorithm computes a horizontal edge to approximate the edge with 5 degree tilt angle. So, the performance of the proposed algorithm depends on the distribution of visible edges (so as the invisible ones), in terms of both shape and direction.
This distribution, however, has dependence on the coordinate system that we build to describe the scanning system. One could rotate the Oxy coordinates system to change the distribution of visible edges (so as the invisible ones) to reveal potential performance of the algorithm. Taking the real data test for example, the plate phantom has been placed almost horizontal such that the invisible edges coincide well with the propagating direction, so as to easy edge information propagation. If the phantom is tilted a small degree, e.g. 2 degrees, the reconstruction quality would have been degraded.
This should not raise the concern that the proposed algorithm could only recover perfect horizontal invisible edges. It just indicates that if the true edges are close to horizontal, then the proposed algorithms has the tendency to approximate them by horizontal edges. Note that solving Subproblem-2 is just one step of the proposed algorithm, and edge information propagation along the vertical direction is also performed similarly by solving Subproblem-3. The algorithm's ability for recovering invisible tilt edges have been validated by other experiments, such as the ones for inverse crime tests. In fact, for the limitation test of section 3.1.4, when increasing the angular range from 120 degree to 130 degree, the 5 degree edge is already reconstructed very well.

Conclusions
To remove the blurring and artifacts introduced by the limited-angle CT, we proposed an optim ization model that aims to take the properties of limited-angle data to its full advantage. The x-ray dependent features are roughly combined into a 'global feature', i.e. the vertical edges are visible, while the horizontal edges are invisible, according to the specified scanning configuration. This global feature is utilized to devise two regularization terms in an optimization model. Then the alternating minimization technique is employed to tackle the optimization problem, which results in the alternating edge-preserving diffusion and smoothing algorithm (AEDS). Various experiments on simulated and real data suggest that the proposed AEDS algorithm outperforms state-of-the-art algorithms, e.g. TV l1 , DART, Potts model and TV l0 for practical applications.
Generally speaking, the performance of AEDS depends on the distribution of visible edges. The more visible the edges and the more evenly distributed they are, the higher the chance of reconstructing high-quality images.
For now, AEDS seems slow since usually about 1000-2000 iterations are needed. We will investigate its acceleration techniques in our future work.