2.1

The limited-angle CT reconstruction problem

Assume that the size of reconstruction image u is M×N. The vector is a concatenate form along the columns of u, and u_i describes the ith entry of , i = 1, 2,…, J. The CT reconstruction problem could be formulated as solving a linear system

where A ∈ ℝ^I×J is the system matrix and is a vector of length I = V×D which represents the acquired projection data. V and D denote the number of projection views and the number of detector cells, respectively accounts for any measurement bias and additive noise.

For limited-angle data, the linear system (1) with I ≪ J is severely ill-posed, therefore, images reconstructed by conventional reconstruction algorithms will introduce streak and blurring artifacts. This is demonstrated in Fig.(1). Without loss of generality, here, we consider the fan-beam scanning with limited-angle range The rectangle phantom and the scanning configuration are shown in Fig.(1a) and Fig.(1b), respectively. Under this scanning configuration, according to the theory of visible and invisible boundaries, the edges close to vertical are visible and can be easily reconstructed, while for the edges close to be horizontal will be invisible and cannot be recovered well by conventional reconstruction algorithms like FBP or SART, as shown in Fig.(1c).

Figure 1:

Illustration of the theory of visible and invisible boundaries. (a) The rectangle phantom; (b) the scanning configuration; (c) the limited-angle reconstruction for the scanning angular range .

2.2

The DTV model

The DTV algorithm¹² is to solve the following minimization problem

where ∇x and ∇y represent matrices size of J × J, corresponding to the discrete x-direction and y-direction gradient operators, respectively, and t_x and t_y are two scalars, specifying the allowed total variations along the x-direction and y-direction, respectively. Since the model (2) is convex, the CP algorithm could be employed to compute a global minimizer. It should be noted that for limited-angle problems, the matrix A has a very large kernel space, so that the model (2) could possess multiple global minimizers. In this case, different parameterizations or different initializations could lead to different solutions.

2.3

The visible edge aware DTV model (VEA-DTV)

As mentioned earlier, to better utilize the visible edges prior, we reformulate the DTV model (2) as the following one

The parameter ϵ controls the noise-level of reconstructed image, which has a clear physical meaning.¹³ It’s easy to check that

so the proposed VEA-DTV model (3) is theoretically equivalent to the DTV model (2). However, since the two formulations are not the same, when applying the CP algorithm, the resulting solving algorithms would be different. As mentioned earlier, for limited-angle problems, the models (3) and (2) are not strictly convex and since the system matrix A has a large kernel space, each of the two models admits multiple solutions, in which case different algorithms might reach different global minimizers. So, starting from the formulation (3), the CP algorithm might compute a solution different from that of the DTV algorithm. This is also the case when comparing AEDS and DTV. The model AEDS(l₁, l₁) coincides exactly with that of DTV, since condition (4) is met. However, since AEDS and DTV employ different minimization algorithms, their performance could be different. In fact, the alternating minimization algorithm adopted by AEDS takes constant step-sizes, according to the framework of incremental methods,¹⁴ it only converges to a neighbourhood of the optimum. On the other hand, the CP algorithm can be proved to converges to a saddle point corresponding to a optimum.

2.4

Numerical algorithm

The CP algorithm is adapted to develop an iterative algorithm for solving (3) by

where indicator functions are defined as:

Then, the min-max formulation of (5) is given by

Where

Applying the proximal point algorithm to solve (6) and taking an additional extrapolation step, we thus obtain the VEA-DTV algorithm described in Algorithm 1, where ||·||₂ is computed by the power method suggested in,¹³1_J ϵ ℝ ^J denotes the constant vector with all elements set to 1, operator sgn(·) returns the sign of a real number, and l₁Ball_a(·) projects a vector onto the l₁ball with radius of a. The symbol neg(·) represents the negative thresholding function, i.e. projects any positive elements of its argument to zero.

3.1

Inverse crime test

The inverse crime occurs when employing the same forward reconstruction model to generate, as well as to invert, synthetic data. To avoid the inverse crime, analytic projection data are acquired in CTSim. Both noise-free and noisy projection data are tested. The scanning angular range is set to . Poisson noise with incidence intensity I₀ = 1.5 ×10⁵ is added to the analytic projection data.

The results are shown in Fig.2. From left to right, the columns 1 and 2 show the SART (10 iterations) reconstructions, with full data and limited data, respectively, and the columns 3 and 4 show the reconstructions of DTV and VEA-DTV, respectively. As Fig.1 has demonstrated, in the limited-angle reconstructions, the invisible edges are too blurred to be recognized. The first row and second row show the noise-free and noisy reconstructions, respectively. We can easily observe that for the noise-free case, DTV and VEA-DTV achieves similar high quality reconstructions, while for the noisy case, VEA-DTV demonstrates superior results. Distortions and blurring could be easily recognized in the DTV reconstructions, especially at the right bottom part. For the proposed VEA-DTV, blurring has been completely removed, and just small local distortions could be identified along the diagonal of the big parallel gram. Same conclusion could be drawn from the quantitative measures listed in Table 1.

Figure 2:

The analytic rectangle phantom. From left to right, the images are reconstructed by full-angle SART, SART, DTV, VEA-DTV, respectively. The first row shows the reconstructions without noise, while the second row shows the results with added Poisson noise, with incidence intensity I₀ = 1.5 ×10⁵ The display window is set to [0, 0.5].

Table 1:

PSNR, SSIM and NRMSE for the analytic rectangle phantom.

3.2

Invisible edges recovery capability test

One rhombus phantom with tilt angle of 5 degrees is constructed in CTSim. Its projection data without noise are also acquired in CTSim. Both of them are analytic, which are used to test VEA-DTV’s capacity to recover invisible edges. Since the boundary between the two triangles are completely invisible and not distributed along the axes, it is quite challenging to recover it.

The results are shown in Fig.3. From top to bottom, the images are reconstructed with angular ranges 130, 120 and 110 degrees, respectively. From left to right, the columns 2, 3, and 4 show the reconstructions by SART, DTV and VEA-DTV, respectively. As shown in the second column, the SART method fails to recover the invisible edge. For both DTV and VEA-DTV, the quality of the reconstructions decreases with reducing angular ranges, as demonstrated in the last two columns of Fig.3. When the scanning angular range is 130 degrees, both DTV and VEA-DTV recover the invisible edge nearly perfectly. However, when reducing the angular ranges, the performance of DTV deteriorates quickly, while the proposed VEA-DTV could demonstrate certain resistances to such changes. Same conclusion could be drawn from the quantitative measures listed in Table 2.

Figure 3:

The rhombus phantom. From left to right, the images are reconstruction results from full-angle SART, SART, DTV, STV, respectively. From up to bottom, each row shows the reconstructions with different angular ranges. The display window is set to [0, 0.18].

Table 2:

PSNR, SSIM and NRMSE for the rhombus phantom.