Interior tomography with curvelet-based regularization

2.1 Interior CT problem

When x-rays pass through an object, the radiation intensity is attenuated. If we assume that the incoming x-ray photons are monoenergetic, the attenuation follows the well-known Beer–Lambert law. A relationship can be obtained between the linear attenuation coefficients and the projection data p = - ln ( I / - I 0 ) = ∫ L f ( x ) dx(1)

where p is the projection data, I₀ is the incoming x-ray intensity, I is the outgoing x-ray intensity, f is the linear attenuation coefficient of the object, and L represents the x-ray path. The CT image reconstruction problem consists of reconstructing the linear attenuation coefficient map f (x) from the measured projection data.

The first generation CT involves parallel-beam geometry. Without loss of generality, here we assume parallel-beam geometry because the projection data collected from narrow or general fan-beam geometry can be converted to parallel-beam projections via rebinning. In a circular parallel-beam CT scan, the object is scanned from different views. In each view, all the x-ray paths are parallel. In general, a parallel-beam CT system can be modeled by the well-known 2D Radon transform R, which integrates a function f on 2D real space over lines [2] Rf ( θ , s ) = R θ f ( s ) = ∫ - ∞ ∞ f ( s θ → + t θ → ⊥ ) dt(2)

where θ denotes the scanning view angle, θ → = ( cos θ , sin θ ) and θ → ⊥ = ( - sin θ , cos θ ). Assuming that the linear attenuation coefficient f is supported on the unit disk centered at the origin, the interior problem is characterized by the measurements in the range |s| ≤ a < 1. The interior problem does not have a unique solution [2]. In this paper, we denote the truncated Radon transform data as R _I with |s| ≤ a < 1. For a more detailed description of the interior problem, please refer to [1].

The curvelet transform is a multiscale directional transform that allows an almost optimal nonadaptive sparse representation of objects with edges [19]. For most natural images, curvelet frame approximation is not as effective as wavelet orthonormal bases. The reason is that the former has the redundancy factor of at least five, and most natural images include certain specific fine structures that are more irregular than piecewise ℂ 2functions (here ℂ 2denotes the space of two times continuously differentiable functions) [20]. However, Candès and Donoho showed that simple thresholding of curvelet coefficients yields nearly optimal approximations for piecewise ℂ 2 images, as opposed to complex triangulation optimization [21].

In the two dimensional real space ℝ 2, let us consider the spatial variable x = (x₁, x₂) and the polar coordinates (r, ω) in the frequency-domain. A curvelet frame can be generated by one basic curvelet ψ_j,0,0 at the scale 2^-j (the number of wedges at this scale is N _j ) using dilations, rotations and translations. ψ_j,0,0 can be defined via its r- and ω-dependence in polar coordinates in the frequency domain. Let ψ ˆ be the Fourier transform of ψ [17], then ψ ˆ j , 0 , 0 ( r , ω ) = 2 - 3 j / 4 W ( 2 - j r ) V ( 2 ⌈ j / 2 ⌉ + 1 ω / - π )(3)

where W (r) is a “radial window” and V (ω) is an “angular window”. Here denotes the smallest integer greater than or equal to x. The windows are smooth, nonnegative and real-valued, supported on r ∈ (1/2, 2)and ω ∈ (-1, 1), and they obey the following admissibility conditions: ∑ j = - ∞ ∞ W 2 ( 2 j r ) = 1 , r ∈ ( 3 / 4 , 3 / 2 ) ; ∑ l = - ∞ ∞ V 2 ( ω - l ) = 1 , ω ∈ ( - 1 / 2 , 1 / 2 ) .

There are many windows satisfying the above condition, and examples (e.g. scaled Meyer windows) can be found in [19].

The curvelet ψ_j,l,k at the scale 2^-j with scale-dependent rotation angle α_j,l = lπ2^-⌈j/2⌉-1, -2^⌈j/2⌉+1 ≤ l ≤ 2^⌈j/2⌉+1 and translational parameters k = ( k 1 , k 2 ) T ∈ ℤ 2 (ℤ denotes the set of integers) can be defined by ψ j , l , k ( x ) = ψ j , 0 , 0 ( R ˜ α j , l ( x - x k j , l ) )(4)

where x k j , l = R ˜ α j , l - 1 ( k 1 2 - j , k 2 2 - j / 2 ) Tis the center of the curvelet, and R ˜ α j , l is the rotation matrix defined by R ˜ α = [ cos α sin α - sin α cos α ], R ˜ α j , l - 1 = R ˜ α j , l T = R ˜ - α j , l.

Because the union of the supports ∪ j , l , k supp ( ψ ˆ j , l , k ) (here supp (·) denotes the support of a function) does not cover the whole frequency domain, we need to define a low-pass coarse scale element [17, 18] ψ ˆ - 1 , 0 , 0 ( r , ω ) = W 0 ( r ) , W 0 2 ( r ) : = 1 - ∑ j = 0 ∞ W 2 ( 2 - j r ) .(5)

The system of curvelets [17, 19] { ψ - 1 , 0 , k : k ∈ ℤ 2 } ∪ { ψ j , l , k : j ∈ N j , - 2 ⌈ j / 2 ⌉ + 1 ≤ l ≤ 2 ⌈ j / 2 ⌉ + 1 , k ∈ ℤ 2 } satisfies the tight frame property. An arbitrary square-integrable function f ∈ L 2 ( ℝ 2 ) can be represented by a curvelet series f = ∑ n 〈 f , ψ n 〉 ψ n , with n = n ( j , l , k ) ,(6)

and the Parseval’s relationship ∑ n | 〈 f , ψ n 〉 | 2 = ∥ f ∥ L 2 ( ℝ 2 ) 2 holds. The components in the set {c _n  : =〈 f, ψ _n  〉 } are called the curvelet coefficients of f.

Based on the CS theory, it was proved that an internal ROI can be reconstructed uniquely and stably if the ROI is piecewise constant/polynomial. As a result, total variation (TV) minimization [10], high-order TV (HOT) minimization [22], and weighted l _p -norm of wavelet coefficients minimization [13] are coupled with data discrepancy minimization to solve the interior problem. Most frequently, CS inspired algorithms first construct an objective function and then minimize it by an iterative method. The objective function usually contains a regularization term, and the choice of this term is one of the key steps. Frequently, one uses a frame that is capable of representing an image in a sparse fashion with a small number of nonzero coefficients, and the regularization term is a norm of these coefficients.

In 2006, Candès et al. [17] proposed two fast discrete curvelet transforms (FDCT). The FDCT and its inverse were used for image reconstruction from noisy [23] and incomplete tomographic data [18]. Considering the discrete version of the reconstruction problem, the unknown object f can be represented as a finite linear combination of curvelets, i.e. f = ∑ n = 1 N 〈 f , ψ n 〉 ψ n .(7)

In the Fourier transform domain, the curvelets are localized with respect to both radial and angular variables (see Fig. 1). This feature has been used for limited angle CT image reconstruction [18]. Compared with the TV minimization method in reference [10], where a piecewise constant image model is assumed, the curvelet transform uses a more practically appealing framework to express the details of an image. Let l _θ denote the radial line in polar coordinates corresponding to the value of the angular parameter θ. Then, we have

R θ ψ n ≡ 0 if supp ( ψ ˆ n ) ∩ l θ = ∅ ,(8)

which implies that the curvelet ψ _n is invisible for the angle θ in the projection domain. In reference [18], to reduce the dimension of the problem, the A-CSR algorithm was constructed that uses only visible curvelet coefficients. The experimental results in [18] demonstrate that the A-CSR algorithm can achieve similar performance as the curvelet sparse regularization, even though the latter uses all curvelet coefficients.

For the interior problem, the truncated projection data p can be expressed as

p = R I f = R I ∑ n = 1 N 〈 f , ψ n 〉 ψ n = ∑ n = 1 N c n R I ψ n .(9)

Denote ρ _n  =∥ R _I ψ _n  ∥/  - ∥ Rψ _n  ∥. Then, this quantity can be used to measure the visibility of the curvelet ψ _n from the projection data. If ρ _n  = 0, the curvelet ψ _n is invisible in the projection domain. If ρ _n  = 1, then the curvelet ψ _n is completely visible in the projection domain. In the intermediate case 0 < ρ _n  < 1, the curvelet is partially visible. The visibility ρ _n also measures data completeness with respect to the curvelet ψ _n . Let S : ={ n|ρ _n  > τ  } (τ ∈ (0, 1)) denote the set of indices of those curvelets whose visibility value is greater than τ. We split the curvelet coefficients in (7) into two sets according to S. Therefore, the unknown object f can be represented as f = ∑ n ∈ S 〈 f , ψ n 〉 ψ n + ∑ n ∉ S 〈 f , ψ n 〉 ψ n .(10)

It is obvious that, when reconstructing f, the curvelet coefficients in the set S are more reliable than other coefficients. Based on the properties of the curvelet frame and inspired by the CS theory and the aforementioned considerations, a curvelet-based regularization model is constructed as follows arg min f { μ 1 ∑ n ∈ S | 〈 f , ψ n 〉 | + μ 2 ∑ n ∉ S | 〈 f , ψ n 〉 | } , s . t . ∥ Af - p ∥ 2 2 = 0 ,(11)

where A = (a_i,j) is a matrix representing the discretized truncated Radon transform R _I , a_i,j is the contribution of the jth pixel to the ith ray-sum, and μ₁ and μ₂ are the regularization parameters. A closely related idea was proposed in [13], where the authors assign weights to basis functions (wavelets) depending on whether the support of a basis function is a subset of the ROI, partially overlaps the ROI, or is outside the ROI. In contrast, here the weight of a basis function is determined by its visibility from the data.

To solve the problem (11), we develop a heuristic iterative algorithm inspired by two alternative minimization procedures whose convergence was proved by Csisz and Tusnády [24] and Daubechies et al. [25]. In [25], each iteration step corresponds to a soft-thresholding filtering of the standard Landweber iteration obtained by adding the backprojected residual error to the previous iteration. This proposed iterative algorithm can also be viewed as an instance of the framework of projection onto convex sets (POCS). The iterative algorithm we propose can be summarized as follows:

S1: Initialization: k = 0, f^(k) = 0, g = 0, t _k  = 1.

Repeating the main loop:

S2: Incorporating the data fidelity term by using the ordered subsets simultaneous algebraic reconstruction technique (OS-SART) [26] { f ˜ ( k ) = f ( k ) f ˜ j ( k + 1 ) = f ˜ j ( k ) + ∑ i ∈ V m a i , j ∑ n ∈ V m a n , j p i - 〈 A i , f ˜ j ( k ) 〉 ∑ l a i , l , ( j = 1 , 2 , … , N ) ,(12)

where i ∈ V _m denotes the ith ray in the mth subset. The reconstructed intermediate image in this substep is denoted f ^OSSART . In [25], the system matrix A is assumed to be normalized for the convenience of theoretical analysis. In CT image reconstruction applications, the adjoint operator A^* can be simplified as the transpose matrix A ^T (backprojection), and the system matrix A usually is not normalized. Therefore, an additional scale factor is required for the control parameter to guarantee convergence. In the OS-SART algorithm, the residual error is normalized before backprojection, and this is equivalent to the normalization of the system matrix. In addition, the ordered subsets method is used to accelerate the convergence, and the non-negativity constraint is used to improve image quality: f j OSSART = 0 , if f j OSSART < 0 , ( j = 1 , 2 , . . . , N ) .

S3: Incorporating the regularization term by updating the result with a thresholding function[25, 27] f TH = ∑ n ∈ S T μ 1 ( 〈 f OSSART , ψ n 〉 ) ψ n + ∑ n ∉ S T μ 2 ( 〈 f OSSART , ψ n 〉 ) ψ n ,(13)

where T _μ is the thresholding function defined by T μ ( x ) = { x + μ if x ≤ - μ 0 if | x | < μ . x - μ if x ≥ μ .(14)

The non-negativity constraint is also used: f j TH = 0 , if f j TH < 0 , ( j = 1 , 2 , . . . , N ) .

S4: Employing the fast weighting to accelerate the convergence [28]. t k + 1 = 1 + 1 + 4 t k 2 2 ;            h = f ^TH ; (h = f ^OSSART , if images are reconstructed by the OS-SART only) f ( k + 1 ) = h + t k - 1 t k + 1 ( h - g ) ;           g = h;

Non-negativity constraint: f j ( k + 1 ) = 0 , if f j ( k + 1 ) < 0 , ( j = 1 , 2 , . . . , N ) . k = k + 1;

Until the stopping criteria are satisfied.

In this work, we refer to the above steps as a curvelet frame based regularization method (CFRM) for short. In the CFRM, τ, μ₁ and μ₂ are three key parameters. Because the curvelets in S are more reliable than others, different values should be assigned to the regularization parameters μ₁ and μ₂. They should be chosen to balance the data fidelity term and the regularization term. There are some general strategies for choosing values of regularization parameters [29]. However, it is difficult to choose “optimal” regularization parameters, and these strategies cannot be used in our model. As far as we know, most regularization coefficients in regularization based CT image reconstruction algorithms are chosen empirically. Therefore, in this work, we follow the most common approach and select these parameters in an empirical fashion.

Figure 2 shows the image of the Shepp-Logan phantom. In the experiment, the parameters τ, μ₁ and μ₂ were considered separately. At first, some values were chosen randomly, and we found that μ₁ = 0.00005 and μ₂ = 0.05 is a relatively good combination. As shown in Table 2, to study the effect of each parameter, three groups of experiments (I, II and III) were designed. The last column in Table 2 shows the minimum RMSE and the corresponding iteration number. The RMSEs vs. iteration numbers were shown in Fig. 3. When only the OS-SART is used, it converges after 43 iterations, and the minimum RMSE is 0.0119. We can see that, in most cases (except for the last one), the CFRM outperforms the OS-SART. Although the performance of CFRM is affected by three parameters, we can see that the choices of τ and μ₂ are more flexible than the choice of μ₁. This is because μ₁ directly affects the set of visible curvelets, which are more important to the ROI than the other curvelets. Overall, μ₁ = 0.00005, μ₂ = 0.005 and τ = 0.3 was found to be a fairly good combination.

To study the convergences of the CFRM (μ₁ = 0.00005, μ₂ = 0.005 and τ = 0.3) and OS-SART, 30000 iterations were used. The objective function and RMSEs vs. iteration number are shown in Fig. 4. Compared with OS-SART, the CFRM converged to a solution with a smaller RMSE. Figure 5 shows the reconstructed results using CFRM (μ₁ = 0.00005, μ₂ = 0.005 and τ = 0.3) after 202 iterations. The CFRM result matches well with the original Shepp-Logan phantom inside the ROI.

To test the stability of the CFRM against data noise, Poisson noise was added to the simulated projections of the Shepp-Logan phantom, assuming 10⁴ photons per detector element. We repeated the aforementioned experiments. Figure 6 shows the reconstructed image using CFRM (μ₁ = 0.00005, μ₂ = 0.005 and τ = 0.3) after 30 iterations. Although the result contains some noise, the gray level is very close to the ground truth of the original phantom inside the ROI.

3.2 Thorax image phantom experiment

The thorax image phantom was reconstructed from a clinical dataset, which was in the fan-beam geometry with 888 detector elements. To save the computational cost, the 256×256 image matrix was first reconstructed by the OS-SART algorithm and normalized as shown in Fig. 7. Then, the phantom image was used to generate parallel-beam projections. In order to investigate the robustness of the three parameters for different phantoms, the same set of parameters μ₁ = 0.00005, μ₂ = 0.005 and τ = 0.3 in section 3.1 were used in this experiment. The minimum RMSE of the CFRM is 0.0078 (after 895 iterations). Similar to Fig. 5, the results are shown in Fig. 8. We can see that the ROI is well reconstructed by the CFRM, and the parameters are still a good choice for the new phantom.

3.3 Modified Shepp-Logan phantom experiment

Recently, Klann et al. developed a wavelet method with a weighted sparsity penalty for interior tomography [13]. In this work, a weighted l _p -penalty is used, and the weight of a wavelet depends on its location relative to the ROI. For the choice of weights, they consider three types of wavelet basis functions: 1) whose support either contains the ROI or is contained in the ROI. The weights are equal to 1 for these basis functions. 2) whose support does not have overlap with the ROI. A fixed value for the corresponding weights is chosen. 3) whose support has an overlap with the ROI, but not of the first type. The weight interpolates between 1 and the value chosen for the outside weight. Klann et al. use a modified Shepp-Logan phantom [10] to test their algorithm. In order to compare the performance of the proposed CFRM to the method in [13], in this experiment the same modified Shepp-Logan phantom and the associate parameters are used to generate projections. The same parameters μ₁ = 0.00005, μ₂ = 0.005 and τ = 0.3 were used for the CFRM and the reconstructed result after 550 iterations is shown in Fig. 9. While the reconstructed ROI by Klann et al. (see Fig. 6 in [13]) still has visible shift, the CFRM result matches the original phantom very well inside the ROI. Since Klann et al. have showed that their method outperforms the TV based method [13], it is indirectly verified that the proposed method also outperforms the TV based method.