Interior tomography with curvelet-based regularization

Abstract

The interior problem, i.e. reconstruction from local truncated projections in computed tomography (CT), is common in practical applications. However, its solution is non-unique in a general unconstrained setting. To solve the interior problem uniquely and stably, in recent years both the prior knowledge- and compressive sensing (CS)-based methods have been developed. Those theoretically exact solutions for the interior problem are called interior tomography. Along this direction, we propose here a new CS-based method for the interior problem based on the curvelet transform. A curvelet is localized in both radial and angular directions in the frequency domain. A two-dimensional (2D) image can be represented in a curvelet frame. We employ the curvelet transform coefficients to regularize the interior problem and obtain a curvelet frame based regularization method (CFRM) for interior tomography. The curvelet coefficients of the reconstructed image are split into two sets according to their visibility from the interior data, and different regularization parameters are used for these two sets. We also presents the results of numerical experiments, which demonstrate the feasibility of the proposed approach.

Keywords

Computed tomography (CT)local reconstruction interior tomography radon transform curvelet transform regularization method

1 Introduction

While the classical computed tomography (CT) theory targets exact reconstruction of the whole cross-section or of an entire object from complete projections, practical applications are often focused on a much smaller internal region-of-interest (ROI) and result in an interior problem, i.e. reconstructing from local truncated projections. Because in this ROI-focused interior problem only the ROI is illuminated, it can help to handle large objects, minimize radiation dose, suppress scattering artifacts, enhance temporal resolution, reduce system cost, and increase scanner throughput [1].

In our opinion, development of local reconstruction algorithms for solving the interior problem can be divided into the following three stages. Stage (I) is approximate reconstruction. For a long time, approximate local reconstruction methods were widely employed because the interior problem does not have a unique solution [2]. Stage (II) is knowledge based interior tomography. The interior problem is exactly and stably solved by assuming a known sub-region inside the ROI [3 –9]. This knowledge-aided theoretically exact solution to the interior problem is called interior tomography. Although the CT numbers of certain sub-regions such as air in a trachea and blood in an aorta indeed can be assumed known, to obtain precise knowledge of a sub-region can be difficult in many cases such as in studies of rare fossils or certain biomedical structures [10]. Stage (III) is compressive sensing (CS) based interior tomography. Inspired by the CS theory, it was proved that an internal ROI can be reconstructed uniquely and stably if the ROI is piecewise constant/polynomial [10 –13]. Interior tomography reconstruction algorithms can be classified into two categories: iterative reconstruction methods [4 , 14] and singular value decomposition (SVD) based differentiated back-projection (DBP) approaches [6 , 16].

In recent years, it was reported that two-dimensional (2D) square integrable functions can be represented by a linear combination of curvelets, which are localized in both radial and angular directions in the frequency domain [17]. With a limited scanning angle, if a curvelet’s support in the frequency domain is outside the visible region, its Radon transform always equals zero. In [18], this property was used to construct an adapted curvelet sparse regularization (A-CSR) algorithm to solve the limited angle problem of tomography. For the interior problems, data truncation affects the visibility of curvelets in the projection domain. Curvelets with higher visibility will be more reliable and more accurate when reconstructing the ROI than curvelets with lower visibility. This motivates us to develop a new iterative reconstruction algorithm for the interior problem by incorporating a curvelet-based regularization. In this pilot study, the coefficients of the curvelet expansion of the reconstructed image are split into two sets according to their visibility from the interior data, and their L₁-norms are used to regularize the interior reconstruction. This is in the spirit of reference [13], where a different strategy was used to split the wavelets into three sets with different regularization parameters.

The rest of this paper is organized as follows. In the next section, the interior problem is defined and the curvelet frame based reconstruction method is presented. In section 3, extensive numerical simulations are performed with different datasets and parameters. Relevant issues are discussed and the paper is concluded in the last section.

2 Method

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}) = \int_{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) = \int_{- \infty}^{\infty} f (s \vec{θ} + t {\vec{θ}}^{⊥}) dt$ (2)

where θ denotes the scanning view angle, $\vec{θ} = (cos θ, sin θ)$ and ${\vec{θ}}^{⊥} = (- 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].

2.2 Curvelet frame

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 $ℂ^{2}$ functions (here $ℂ^{2}$ denotes 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 $\hat{ψ}$ be the Fourier transform of ψ [17], then ${\hat{ψ}}_{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: $\sum_{j = - \infty}^{\infty} W^{2} (2^{j} r) = 1, r \in (3 / 4, 3 / 2); \sum_{l = - \infty}^{\infty} V^{2} (ω - l) = 1, ω \in (- 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} \in ℤ^{2}$ ( $ℤ$ denotes the set of integers) can be defined by $ψ_{j, l, k} (x) = ψ_{j, 0, 0} ({\tilde{R}}_{α_{j, l}} (x - x_{k}^{j, l}))$ (4)

where $x_{k}^{j, l} = {\tilde{R}}_{α_{j, l}}^{- 1} {(k_{1} 2^{- j}, k_{2} 2^{- j / 2})}^{T}$ is the center of the curvelet, and ${\tilde{R}}_{α_{j, l}}$ is the rotation matrix defined by ${\tilde{R}}_{α} = [\begin{matrix} cos α sin α \\ - sin α cos α \end{matrix}]$ , ${\tilde{R}}_{α_{j, l}}^{- 1} = {\tilde{R}}_{α_{j, l}}^{T} = {\tilde{R}}_{- α_{j, l}}$ .

Because the union of the supports $\underset{j, l, k}{\cup} supp ({\hat{ψ}}_{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] ${\hat{ψ}}_{- 1, 0, 0} (r, ω) = W_{0} (r), W_{0}^{2} (r) : = 1 - \sum_{j = 0}^{\infty} W^{2} (2^{- j} r) .$ (5)

The system of curvelets [17, 19] ${ψ_{- 1, 0, k} : k \in ℤ^{2}} \cup {ψ_{j, l, k} : j \in N_{j}, - 2^{⌈ j / 2 ⌉ + 1} \leq l \leq 2^{⌈ j / 2 ⌉ + 1},$ $k \in ℤ^{2}}$ satisfies the tight frame property. An arbitrary square-integrable function $f \in L^{2} (ℝ^{2})$ can be represented by a curvelet series $f = \sum_{n} 〈 f, ψ_{n} 〉 ψ_{n}, with n = n (j, l, k),$ (6)

and the Parseval’s relationship $\sum_{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.

2.3 Proposed algorithm

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 = \sum_{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} \equiv 0 if supp ({\hat{ψ}}_{n}) \cap l_{θ} = \emptyset,$ (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} \sum_{n = 1}^{N} 〈 f, ψ_{n} 〉 ψ_{n} = \sum_{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 = \sum_{n \in S} 〈 f, ψ_{n} 〉 ψ_{n} + \sum_{n \notin 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 $\underset{f}{arg min} {μ_{1} \sum_{n \in S} | 〈 f, ψ_{n} 〉 | + μ_{2} \sum_{n \notin 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] ${\begin{matrix} {\tilde{f}}^{(k)} = f^{(k)} \\ {\tilde{f}}_{j}^{(k + 1)} = {\tilde{f}}_{j}^{(k)} + \sum_{i \in V_{m}} \frac{a_{i, j}}{\sum_{n \in V_{m}} a_{n, j}} \frac{p_{i} - 〈 A_{i}, {\tilde{f}}_{j}^{(k)} 〉}{\sum_{l} a_{i, l}}, (j = 1, 2, \dots, N), \end{matrix}$ (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} = \sum_{n \in S} T_{μ_{1}} (〈 f^{OSSART}, ψ_{n} 〉) ψ_{n} + \sum_{n \notin S} T_{μ_{2}} (〈 f^{OSSART}, ψ_{n} 〉) ψ_{n},$ (13)

where T_μ is the thresholding function defined by $T_{μ} (x) = {\begin{matrix} x + μ & if x \leq - μ \\ 0 & if | x | < μ . \\ x - μ & if x \geq μ . \end{matrix}$ (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} = \frac{1 + \sqrt{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 + \frac{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.

3 Results

To verify the feasibility of the proposed local reconstruction method for solving the interior problem, we implemented it using Matlab and performed extensive numerical simulations. A discrete curvelet transform package (downloaded from www.curvelet.org) was used for the curvelet transform and its inverse transform. In our experiments, three test objects were used: 1) Shepp-Logan phantom which was generated by using the Matlab function “phantom(256)”. Poisson noise was also added into the projections to simulate realistic situations. 2) Thorax image phantom which was reconstructed from a clinical dataset [30]. 3) Modified Shepp-Logan phantom [10]. The experiment setups are summarized in Table 1. All noise-free projections were generated using a ray-driven model. For each view, the projection is truncated on two sides to simulate the interior problem.

The root mean square error (RMSE) was adopted to quantitatively evaluate the reconstructed image quality $RMSE (f^{ROI}, {\hat{f}}^{ROI}) = \sqrt{\frac{1}{J} \sum_{j = 1}^{J} {(f_{j}^{ROI} - {\hat{f}}_{j}^{ROI})}^{2}},$ (15)

where ${\hat{f}}^{ROI}$ represents the ROI of a reconstructed image, and f^ROI is the ROI of the reference image. Because the goal is to reconstruct the interior ROI, the RMSE is calculated only for the interior ROI. To investigate the effect of the regularization term, images were also reconstructed by the OS-SART only, that is, the CFRM without the thresholding filtering sub-step. 16 ordered subsets were used in the OS-SART step. To eliminate the effect of model mismatch, the same ray driven model was used in both projection and backprojection routines.

3.1 Shepp-Logan phantom experiment

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.

4 Discussion and conclusions

When the CFRM is used to solve the interior problem, choosing the values of parameters and the number of iterations is important. For clinical applications, the targets are human beings which have structural similarity. We can train those parameters (including the number of iterations) on a large number of sample cases. As we can see from the experimental results, the RMSE changes gradually when the parameters change. Although it is difficult to find the “optimal” result, it is not hard to find a satisfactory one.

Based on the work of Daubechies et al. [25], we analyzed the SART-type reconstruction framework from a limited number of projections with a sparsity constraint [31]. By introducing two diagonal matrices to normalize the CT system matrix, the convergence of SART-type reconstruction was established by following the same steps as in [25]. Because there is only one regularization parameter and term in [31] and there are two regularization parameters and terms in this paper, the theoretical analysis of convergence in [31] cannot be directly applied. However, it is our hypothesis that the convergence of the proposed algorithm can be established with additional efforts to appropriately extend the work of [31]. Meanwhile, the proposed algorithm is not identical to the one in [31]. While a steepest gradient method was implemented to find the optimal regularization parameter in each iteration (see equations 11–15 in [31]) for the fastest convergence, in this paper we did not implemented this steepest gradient method for the regularization parameter. In fact, it can be viewed as the original version of [25] with a SART-type reconstruction for the data fidelity term. Regarding the step size in the SART, a default value “1” is employed.

Our MATLAB codes were developed on a PC with Xeon CPU E5-2620 and a 16GB RAM. Currently, the computational cost is a major issue for the proposed CFRM. In our experiment with the Shepp-Logan phantom, the average computational time for each iteration is 2.79 seconds. To make the CFRM practical, the projection/backprojection, soft thresholding, curvelet transform and inverse curvelet transform steps used in our implementation need to be optimized and accelerated. The graphics processing unit (GPU) and other hardware-based high-performance computing technologies can be adopted. Given the fact that this is one of the first-of-its-kind papers along this direction, the software and hardware acceleration for the CFRM are beyond the scope of this paper. Meanwhile, we need to convert the CFRM for parallel-beam geometry to fan-beam or cone-beam geometry.

In recent years, many wavelet constructions (e.g. contourlet, surfacelet, shearlet) have been developed and used for different image processing tasks. Collectively, these wavelets are called X-lets. The relationship of curvelets to other X-lets was briefly discussed in [19]. Although it is still an open problem to prove the uniqueness and stability, we believe that other X-lets can also be applied for interior tomography. However, because different X-lets correspond to different image sparsity, these X-lets may be valid for different imaging objects satisfying different sparse models.

In conclusion, we proposed a curvelet-based regularization method for interior tomography. Numerical experiments with different datasets demonstrate the feasibility of our approach. Because the CFRM is a heuristic method without a rigorous proof, in the future we will study more sophisticated optimization algorithms to make it more suitable for practical applications.

Footnotes

Acknowledgments

B.D. Liu was sponsored in part by the Instrument Developing Project of the Chinese Academy of Sciences YZ201410, in part by IHEP-CAS Scientific Research Foundation 2013IHEPYJRC801 and in part by SRF for ROCS, SEM. A. Katsevich was supported in part by NSF Grant DMS No. 1211164. H.Y. Yu was supported in part by NSF Grant DMS No. 1619550.

References

Wang

, Yu

, The meaning of interior tomography, Physics in Medicine and Biology 58 (2013), R161.

Natterer

, The mathematics of computerized tomography, Society for Industrial and Applied Mathematics, Philadelphia, 2001.

, Yu

, Wang

, Exact interior reconstruction with cone-beam CT, International Journal of Biomedical Imaging 2007 (2007), 1–6.

, Yu

, Wei

, Wang

, A general local reconstruction approach based on a truncated Hilbert transform, International Journal of Biomedical Imaging 2007 (2007), Article ID: 63634, 63638 pages.

, Yu

H.Y.

, Wang

, Exact interior reconstruction from truncated limited-angle projection data, International Journal of Biomedical Imaging 2008 (2008), Article ID: 427989, 427986 Pages.

H.Y.

, Ye

Y.B.

, Wang

, Interior reconstruction using the truncated Hilbert transform via singular value decomposition, Journal of X-Ray Science and Technology 16 (2008), 243–251.

Courdurier

, Noo

, Defrise

, Kudo

, Solving the interior problem of computed tomography using a priori knowledge, Inverse Problems 24 (2008), Article ID 065001, 065027 pages.

Kudo

, Courdurier

, Noo

, Defrise

, Tiny a priori knowledge solves the interior problem in computed tomography, Physics in Medicine and Biology 53 (2008), 2207–2231.

Gong

, Liu

, Yu

, Liu

, Zhou

, Kan

, He

, Cao

, Interior tomographic imaging of mouse heart in a carbon nanotube micro-CT, Journal of X-ray Science and Technology 24 (2016), 549–563.

10.

, Wang

, Compressed sensing based interior tomography, Physics in Medicine and Biology 54 (2009), 2791–2805.

11.

, Yang

, Jiang

, Wang

, Supplemental analysis on compressed sensing based interior tomography, Phys Med Biol 54 (2009), N432–N425.

12.

Han

, Yu

, Wang

, A total variation minimization theorem for compressed sensing based tomography, International Journal of Biomedical Imaging 2009 (2009), Articel ID:125871, 125873 pages.

13.

Klann

, Quinto

E.T.

, Ramlau

, Wavelet methods for a weighted sparsity penalty for region of interest tomography, Inverse Problems 31 (2015).

14.

, Mou

, Wang

, Sieren

, Hoffman

E.A.

, Yu

, Statistical Interior Tomography, IEEE Transactions on Medical Imaging 30 (2011), 1116–1128.

15.

Katsevich

, Singular value decomposition for the truncated Hilbert transform, Inverse Problems 26 (2010), 115011.

16.

Jin

, Katsevich

, Yu

, Wang

, Li

, Chen

, Interior Tomography With Continuous Singular Value Decomposition, IEEE Transactions on Medical Imaging 31 (2012), 2108–2119.

17.

Candes

, Demanet

, Donoho

, Ying

L.X.

, Fast discrete curvelet transforms, Multiscale Modeling & Simulation 5 (2006), 861–899.

18.

Frikel

, Sparse regularization in limited angle tomography, Applied and Computational Harmonic Analysis 34 (2013), 117–141.

19.

, Plonka

, The Curvelet Transform, IEEE Signal Processing Magazine 27 (2010), 118–133.

20.

Mallat

, A wavelet tour of signal processing, Academic press, 1999.

21.

Candes

E.J.

, Donoho

D.L.

, New tight frames of curvelets and optimal representations of objects with piecewise C-2 singularities, Communications On Pure And Applied Mathematics 57 (2004), 219–266.

22.

Yang

, Yu

, Jiang

and Wang

, High-order total variation minimization for interior tomography, Inverse Problems 26 (2010).

23.

Candes

E.J.

, Donoho

D.L.

, Curvelets and reconstruction of images from noisy radon data, Wavelet Applications in Signal and Image Processing Viii Pts 1 and 2 4119 (2000), 108–117.

24.

Csisz

, Tusnády

, Information geometry and alternating minimization procedures, Statistics and decisions (1984).

25.

Daubechies

, Defrise

, De Mol

, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Communications On Pure And Applied Mathematics 57 (2004), 1413–1457.

26.

Kak

A.C.

, Slaney

, Principles of computerized tomographic imaging, Society for Industrial and Applied Mathematics, Philadelphia, 2001.

27.

, Wang

, A soft-threshold filtering approach for reconstruction from a limited number of projections, Physics in Medicine and Biology 55 (2010), 3905–3916.

28.

Beck

, Teboulle

, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, Siam Journal on Imaging Sciences 2 (2009), 183–202.

29.

Scherzer

, Handbook of mathematical methods in imaging Springer, New York; London, 2011.

30.

, Wang

, Hsieh

, Entrikin

D.W.

, Ellis

, Liu

, Carr

J.J.

, Compressive Sensing–Based Interior Tomography: Preliminary Clinical Application, Journal Of Computer Assisted Tomography 35 (2011), 762.

31.

, Wang

, SART-type image reconstruction from a limited number of projections with the sparsity constraint, International Journal of Biomedical Imaging 2010 (2010), Article ID: 934847, 934849 pages.