Source and coded aperture joint optimization for compressive X-ray tomosynthesis

Xu Ma; Qile Zhao; Angela Cuadros; Tianyi Mao; Gonzalo R. Arce

doi:10.1364/OE.27.006640

1. Introduction

X-ray computed tomography (CT) has been extensively used in clinical diagnosis, security inspection, as well as industrial examination, detection, analysis and so on [1–6]. Traditional CT makes a complete set of measurements around the object, and reconstructs the content of the object using computer-processed algorithms, such as the filtered back-projection algorithm (FBP) [1, 7]. However, numerous measurements are generally required in X-ray CT leading to high radiation exposure, long acquisition time, and numerous angle views that increase the geometrical complexity. In some applications, the geometry of the object under inspection prevents a full set of measurements around the object. As an alternative technique to X-ray CT in such scenarios, X-ray tomosynthesis makes an incomplete set of projections from limited angle views, from which the three-dimensional (3D) object can be reconstructed based on iterative algorithms [8–10]. Tomosynthesis has several merits including lower levels of radiation dose, fewer measurements, and simpler sensing geometry [8, 9, 11]. Radiation dose is, in particular, a critical concern in practical clinical applications [12, 13]. The incomplete measurement strategy of tomosynthesis is beneficial to reduce radiation, but inevitably leads to an ill-posed inverse reconstruction problem [7, 14]. Thus, sparsity priors and total variation regularizations have been used to improve the image reconstruction performance in tomosynthesis imaging [9, 15–17]. In the past, Kaganovsky et al. studied and compared four different subsampling measurement strategies in X-ray tomography imaging process: the uniform view subsampling, random view subsampling, uniform detector subsampling, and random detector subsampling [18]. It is proven that the random detector subsampling performs best, which can be physically implemented by the coded aperture illumination used in this paper [18, 19]. A noise analysis based on a singular decomposition analysis was provided in [18].

Choi and Brady recently introduced the concept of compressive X-ray tomosynthesis, to alleviate the tradeoff between the small number of measurements and the ill-posedness of the inverse reconstruction problem [20]. As shown in Fig. 1, compressive tomosynthesis exposes the 3D object by a set of X-ray sources distributed on a source plane. A coded aperture is placed under each of the X-ray sources to modulate the light rays transmitted through the object. Typically, binary coded apertures with random patterns are used to block a portion of X-rays in order to reduce the radiation dose. The coded apertures in X-ray systems can be implemented by 3D printed Tungsten masks [21–23]. Then, the coded illumination is attenuated by the 3D object and then measured by an integrating two-dimensional (2D) detector array. Based on the sparsity assumption of the objects of interest, the 2D measurements received by the detector can be used to recover the 3D object using compressive sensing (CS) reconstruction algorithms [24–27]. While random coded projections are adequate in general CS scenarios, it has been shown that they are suboptimal in cases where the sensing matrices are highly structured, as in this case. Recently, Cuadros, et al. developed a coded aperture optimization method to improve the reconstruction performance of tomosynthesis [27]. Using the direct binary search algorithm, optimized coded apertures were obtained to achieve uniform sensing and complementary codes for multiple shots. More recently, the concept of code optimization was extended to the general tomography imaging system, where a gradient-based numerical optimization and several regularization approaches were utilized to optimize the coded apertures [28].

Fig. 1 Sketch of compressive tomosynthesis system, which exposes the 3D object by a set of X-ray sources distributed on a source plane. A coded aperture is placed under each source to modulate the light rays. Then, the coded illumination is attenuated by the 3D object and then measured by an integrating 2D detector array.

Download Full Size | PDF

To date, the optimization of compressive X-ray tomosynthesis has been limited to the design of coded aperture structures. However, the locations where the X-ray sources are placed, their radiation angles, in addition to the coded apertures can be shown to have a significant impact on the attained imaging quality. This paper thus focuses on the jointly optimization of all three imaging components, which have an influence of the sensing matrix. Thus, the source locations and the incident angles will be optimized determining the configuration of the X-ray sources. The coded aperture optimization associated with each source is synergistically optimized to obtain their transmissivities and structures. Based on this concept, this paper develops a gradient-based source and coded aperture joint optimization (SCO) method for compressive tomosynthesis systems to break through the limitations of current code optimization algorithms. A set of X-ray sources are deployed on the source plane, and these sources are turned on successively to obtain multiple projection measurements on the detector. In order to solve for the tomosynthesis problem using a linear sensing model, at any moment only one radiation source is allowed. In addition, the number of sources is set as constant to fix the emitted radiation dose. Given the positions of the object and detector, the locations of X-ray sources determine the sensing geometry and the system matrix of the tomosynthesis architecture. On the other hand, the incident angles of the X-ray beams determine the orientation of the corresponding source. The incident angles determine the directions from which the object is exposed, and which portions of the object are projected. Finally, the coded apertures are used to modulate the structure illuminations and reduce the radiation dose.

The source locations, incident angles, and the coded apertures are jointly optimized based on the restricted isometry property (RIP) condition and the incoherence of the sensing matrix, both well-established methods in CS theory [29–32]. In particular, the cost function is set as the Frobenius norm of the difference between the point spread function (PSF) of the system and the identity matrix, and the gradient-based algorithm is used to solve for the optimization problem. During the optimization, variables of the source locations and coded apertures are first relaxed to continuous ranges, and then the optimized variables are quantized into discrete domains. A set of regularization methods are used to reduce the quantization error, and to obtain the uniform sensing on the object. Based on the sparsity assumption, the gradient projection for sparse reconstruction (GPSR) algorithm is used to reconstruct the 3D object from the compressive measurements under the optimized settings [33]. Simulations prove that the proposed SCO method can effectively improve the reconstruction accuracy, and outperform the individual optimization of coded apertures or incident angles. In the following, the forward imaging model of X-ray tomosynthesis is described in Section 2. The proposed SCO method is developed in Section 3. The simulations and analysis are provided in Section 4. The conclusions and future work are presented in Section 5.

2. Forward imaging model of X-ray tomosynthesis

According to the Beer-Lambert law, the X-ray transmission imaging model through an object is given by

I = I_{0} \cdot \exp [- \int_{0}^{\infty} μ (r) \cdot d r],

where I₀ is the original intensity of the X-ray, µ(r) is the linear attenuation coefficient of the object at location r, and I is the measured intensity on the detector after it is transmitted through an object [1]. The transmission imaging model can be transferred to the form:

- \ln (I / I_{0}) = \int_{0}^{\infty} μ (r) \cdot d r

. Let f represent the linear attenuation coefficient map of the 3D object. Suppose

\vec{s}

and

\hat{θ}

represent the location of X-ray source and the direction to illuminate the object, respectively. Then, the imaging model can be rewritten as

- \ln (I / I_{0}) ≜ y (\vec{s}, \hat{θ}) = \int_{0}^{\infty} f (\vec{s} + r \hat{θ}) \cdot d r .

The continuous form of the imaging model in Eq. (2) is referred to as the X-ray transform of the object. In order to apply numerical algorithms to solve for the object reconstruction problem, the imaging model is discretized as follows. First, the 3D object is gridded into N = N_x × N_y × N_z voxels, where N_x, N_y and N_z are the dimensions of the object along the x, y and z axes, respectively. Let f ∈ ℝ^{N ×}¹ be the raster-scanned vector representation of the 3D attenuation coefficients corresponding to an object. The detector is gridded into M = M_x × M_y pixels, where M_x and M_y are the dimensions along the x and y axes. Let y ∈ ℝ^M×¹ be the raster-scanned vector representing the sequential measurements on the detector. Then, the discrete form of the imaging model is formulated as

y = Hf,

where H ∈ ℝ^M×N is termed as the system matrix. The (m, n)th element of H represents the volume portion of the nth voxel that is irradiated by the X-ray associated with the mth detector pixel. The structure of H is determined by the sensing geometric relationship between the source, object and detector in the 3D free space. Specifically, each row of H corresponds to the information collected by a certain detector pixel, while each column of H corresponds to the weight coefficients of a particular voxel [27].

As shown in Fig. 1, the proposed X-ray tomosynthesis system deploys multiple X-ray sources on the source plane above the object. Each source emits an X-ray cone-beam to expose the underlying object. In this paper, the source plane is gridded into an N_s × N_s pixel array referred to as the source pattern. The source pattern is denoted by a binary matrix $J \in ℝ^{N_{s} \times N_{s}}$ , where J_i means the ith pixel on J. If an X-ray source is placed at pixel J_i, then J_i = 1, otherwise J_i = 0. Suppose there are P X-ray sources distributed on the source plane, which means J includes P one-valued pixels, or ║J║₀ = P where ║║₀ denotes the l₀-norm. Considering the imaging process of the pth source point, the corresponding measurements on the detector are y_p = H_pf, where H_p is the system matrix corresponding to the pth source.

In the optimization of the tomosynthesis imager, three kinds of optimization freedoms are taken into account, i.e., the source locations, source orientations (or incident directions), and the coded apertures. As shown in Fig. 2(a), a coded aperture is placed in front of each X-ray source to reduce the radiation dose and modulate the projection pattern on the detector. The coded apertures for different sources are not necessarily the same. It is noted that the thickness of the X-ray masks is a very important factor for the implementation of the coded apertures. The thickness of the masks depends on the energy of the X-ray source and the materials of the masks. Given the X-ray source energy used and the mask materials, we can calculate the necessary thickness of the masks to block the X-ray radiation. It is natural to ask if we could attenuate the radiation by directly reducing the dose of X-ray sources but use full-sampling to interrogate the object without coded apertures. However, in the presence of noise this method will degrade the signal-to-noise ratio (SNR) of the measurement on each pixel of the detector. Without using the coded aperture illumination, we cannot optimize the sensing matrix of the X-ray tomosynthesis systems to further improve the reconstruction performance. Mathematically, the coded aperture associated to the pth source is represented by an M_x × M_y binary matrix T_p, where the elements with values 0 and 1 represent the blocked and unblocked pixels, respectively. The dimension of the coded aperture is the same as that of detector, and thus the pixels on the coded aperture have one-to-one correspondence to the pixels on the detector. All of the pixel values of the coded apertures can be flipped during the optimization process to modulate the cone-beam projection. On the other hand, the source locations are represented by the source pattern J, all pixels of which can be turned on or turned off. As illustrated in Fig. 2(b), turning off a point source on J is equivalent to blocking all pixels on the corresponding coded aperture. In addition, the orientation of each source can be modified to obtain the optimal incident angles of X-rays. Suppose the exposed area on the coded aperture by the X-ray source is a circular region, which is shown as the red shaded area in Fig. 2(c). Consequently, there are no X-rays passing through the coded aperture pixels outside the exposed area. Thus, the pixels outside the exposed area can be treated as blocked pixels. Rotating the source is equivalent to shifting the exposed area on the coded aperture, which is shown by the blue circle in Fig. 2(c).

Fig. 2 The optimization freedom of the proposed compressive tomosynthesis system: the optimization of (a) the coded aperture, (b) source pattern, and (c) source orientation.

Download Full Size | PDF

The principal idea of co-optimization is to integrate all of the optimization degree of freedom into the variables of coded apertures. Then, the equivalent transmission function of the coded aperture for the pth source can be formulated as

{\tilde{T}}_{p} = \prod_{p} ⊙ Λ_{p} ⊙ T_{p},

where ⊙ is the element-by-element multiplication, Π_p, Λ_p and T_p are the M_x × M_y transmission functions corresponding to the source pattern, source orientation and coded aperture, respectively. All elements in the matrix Π_p have the same value. If the pth source is turned on, Π_p is set to be a one-valued matrix, otherwise Π_p is a zero-valued matrix. Λ_p is a circular window function to represent the exposed area on the coded aperture. The central coordinate of the pth circular function is (x_p, y_p). The pixels of Λ_p within the circular function have the value of 1, while other pixels have the value of 0. T_p is a binary matrix representing the coded aperture. Intuitively,

{\tilde{T}}_{p}

represents the overlap between the exposed area and the coded aperture when the pth source is turned on. Otherwise,

{\tilde{T}}_{p}

is a zero matrix if the source is off.

For the equivalent coded aperture ${\tilde{T}}_{p}$ , we define a code matrix C_p ∈ ℝ^M×M (M = M_x · M_y), which is a diagonal matrix with the diagonal elements equal to the vectorized representation of ${\tilde{T}}_{p}$ . Figure 3 illustrates the generation of the code matrix. In Fig. 3(a), the X-ray source exposes the right-bottom corner of the coded aperture within the red dashed circle. Then, we stack the pixels of the coded aperture into a vector and place it along the diagonal of C_p. Notice that only the unblocked pixels in the exposed areas are counted. In Fig. 3(b), the X-ray source is turned off, thus all pixels in ${\tilde{T}}_{p}$ are equivalently blocked and C_p becomes a zero-valued matrix. Involving the code matrix, the transmission imaging model under the pth source is modified as

y_{p} = C_{p} H_{p} f .

According to Eq. (5), the code matrix C_p selects the rows of system matrix H_p associated with the unblocked detector elements within the exposed area.

Fig. 3 Formation of the code matrix C_p based on the source pattern, source orientation and coded aperture.

Download Full Size | PDF

Taking into account all of the P sources, the measurements on the detector are given by

[\begin{matrix} y_{1} \\ y_{2} \\ ⋮ \\ y_{P} \end{matrix}] = [\begin{matrix} C_{1} & 0 & \dots & 0 \\ 0 & C_{2} & \dots & 0 \\ ⋮ & ⋱ & ⋱ & ⋮ \\ 0 & 0 & \dots & C_{P} \end{matrix}] [\begin{matrix} H_{1} \\ H_{2} \\ ⋮ \\ H_{P} \end{matrix}] f .

For simplicity, Eq. (6) can be written as

y = CHf = C H Ψ θ,

where

y = {[y_{1}^{T}, y_{2}^{T}, \dots, y_{P}^{T}]}^{T}

, C = diag (C₁, C₂,…, C_P), and

H = {[H_{1}^{T}, H_{2}^{T}, \dots, H_{P}^{T}]}^{T}

. Assume f can be sparsely represented on a basis Ψ, which is referred to as the sparse basis. Then, the sparse coefficient vector is denoted by θ, where only a few elements in θ are non-zero or significant [34].

In the following, we define $W = H Ψ = {[w_{1}^{T}, w_{2}^{T}, \dots, w_{M P}^{T}]}^{T}$ , where w_i ∈ ℝ¹^×N (N = N_x · N_y · N_z) is the ith row of W. In addition, we define C in Eq. (7) as [c₁, c₂,…, c_MP ], where c_i ∈ ℝ^MP×¹ represents the i the column of C. It is noted that each column of C, i.e. c_i, contains at most one unblocked pixel. That is because C₁, C₂,…, and C_P in Eq. (6) are all diagonal matrices, thus C is a diagonal matrix, each column of which includes only one unblocked pixel at most. This idea has been illustrated in Fig. 3. If the corresponding element in ${\tilde{T}}_{p}$ is 0, then c_i is a zero vector. Otherwise, c_i is a natural basis in R^M . Then, the sensing matrix of the tomosynthesis system is defined as

A = C H Ψ = CW .

The (m, n)th element in A is

{[A]}_{m, n} = \sum_{i = 1}^{M P} c_{i} (m) w_{i} (n),

where c_i(m) is the mth element in the vector c_i, and w_i (n) is the nth element in the vector w_i.

3. Joint optimization of source and coded aperture

3.1. Cost function of the SCO problem

In order to reconstruct the object f from the compressive measurements, we need to solve for the ill-conditioned problem in Eq. (7). Suppose f can be sparsely represented on the basis Ψ, such that the basis is incoherent with the projection matrix CH. Then, the sparse set of coefficients of f can be recovered by solving for the following optimization problem [34]:

\hat{θ} = \underset{θ}{\arg \min} {‖ y - A θ ‖}_{2}^{2} + γ {‖ θ ‖}_{1},

where ║·║₁ and ║·║₂ respectively represent the l₁ and l₂ norms, and γ is the weight of regularization term.

According to CS theory, the successful reconstruction of coefficients representing the signal f from Eq. (10) is guaranteed by the RIP condition of the sensing matrix. Define $\tilde{A} \in ℝ^{M \times N}$ as the sensing matrix with normalized columns. Define x ∈ ℂ^{N ×}¹ as an s-sparse signal, which means that x contains s non-zero elements. For any s-sparse signal x, the restricted isometry constant δ_s of $\tilde{A}$ is defined as the smallest δ_s satisfying the following inequality

(1 - δ_{s}) {‖ x ‖}_{2}^{2} \leq {‖ \tilde{A} x ‖}_{2}^{2} \leq (1 + δ_{s}) {‖ x ‖}_{2}^{2} .

The matrix $\tilde{A}$ is said to satisfy the RIP condition if δ_s is small for reasonably large s. According to [35], δ_s can be calculate as

δ_{s} = \max_{S \subset [N], | S | \leq s} {‖ {\tilde{A}}_{S}^{H} {\tilde{A}}_{S} - I ‖}_{2},

where [N] means the set of integers from 1 to N, that is [N] ≜ 1, 2, ⋯, N. |S| is the cardinality of S.

{\tilde{A}}_{S}

represents the sub-matrix of

\tilde{A}

that only consists of the columns indexed by S,

{\tilde{A}}_{S}^{H}

represents the conjugate transpose of

{\tilde{A}}_{S}

, and I is the identity matrix. As shown in [13], δ_s is bounded by the Frobenius norm of the difference between the PSF of the system and the identity matrix, such that

δ_{s} \leq {‖ {\tilde{A}}^{H} \tilde{A} - I ‖}_{F},

where

{\tilde{A}}^{H} \tilde{A}

is the PSF of the system. In order to satisfy the RIP condition, we need to minimize the upper bound of δ_s. Thus, the cost function of the optimization problem is set as

F = {‖ {\tilde{A}}^{H} \tilde{A} - I ‖}_{F} .

According to Eqs. (4), (8), and (14), the optimization problem is reduced to searching for the optimal matrices Π_p, Λ_p and T_p for all sources to minimize the cost function, i.e.,

({\hat{Π}}_{p}, {\hat{Λ}}_{p}, {\hat{T}}_{p}) = \underset{Π_{p}, Λ_{p}, T_{p}}{\arg \min} {‖ {\tilde{A}}^{H} \tilde{A} - I ‖}_{F},

where Π_p, Λ_p and T_p have been defined in Eq. (4).

Based on Eq. (9), the inner product between the mth and nth columns of the normalized sensing matrix $\tilde{A}$ is given by:

{[{\tilde{A}}^{T} \tilde{A}]}_{m, n} = \frac{\sum_{i, j = 1}^{M P} 〈 c_{i}, c_{j} 〉 w_{i} (m) w_{j} (n)}{{[\sum_{i, j = 1}^{M P} 〈 c_{i}, c_{j} 〉 w_{i} (m) w_{j} (m)]}^{1 / 2} {[\sum_{i, j = 1}^{M P} 〈 c_{i}, c_{j} 〉 w_{i} (n) w_{j} (n)]}^{1 / 2}},

where (·, ·) is the inner product between two vectors, and w(m) and w(n) are the mth and nth elements in w, respectively. Since C is a diagonal matrix, the inner product 〈c_i, c_j 〉 = 0 if i ≠ j. Therefore, Eq. (16) can be simplified as:

{[{\tilde{A}}^{T} \tilde{A}]}_{m, n} = \frac{\sum_{i = 1}^{M P} b_{i} w_{i} (m) w_{i} (n)}{{[\sum_{i = 1}^{M P} b_{i} w_{i} (m) w_{i} (m)]}^{1 / 2} {[\sum_{i = 1}^{M P} b_{i} w_{i} (n) w_{i} (n)]}^{1 / 2}},

where b_i = 〈c_i, c_i〉. Additionally, the cost function F in [14] is simplified as:

F = \sum_{m = 1}^{N} \sum_{n = 1}^{N} {| \frac{\sum_{i = 1}^{M P} b_{i} w_{i} (m) w_{i} (n)}{{[\sum_{i = 1}^{M P} b_{i} w_{i} (m) w_{i} (m)]}^{1 / 2} {[\sum_{i = 1}^{M P} b_{i} w_{i} (n) w_{i} (n)]}^{1 / 2}} - I_{m, n} |}^{2},

where I_m,n is the (m, n)th element in the identity matrix. As shown in Appendix A, the cost function can be further simplified as follows

F = \sum_{m < n} {\frac{d (m, n)}{{[d (m, m)]}^{1 / 2} {[d (n, n)]}^{1 / 2}}}^{2},

where

d (m, n) = \sum_{i = 1}^{M P} b_{i} R_{i}^{m n}, and R_{i}^{m n} = w_{i} (m) w_{i} (n) .

Assume i = u · M + v, where u = 1, 2,…, P indexes the source, and v = 1, 2,…, M indexes the measurements for a specific source. Let b_i = π_u · λ_i · t_i, where π_u, λ_i, t_i are all binary parameters corresponding to the source pattern, source orientation and coded aperture. In particular, π_u = 1 means the uth source is turned on, otherwise it is turned off. λ_i = 1 means the vth pixel on the coded aperture for the uth source is within the exposed area of the cone-beam projection, otherwise the pixel is out of the exposed area due to the oblique orientation of the source. t_i, where i = u · M + v, is the value of the vth pixel on the coded aperture associated to the uth source. In order to relax Eq. (14) from a combinatorial optimization problem into a continuous one, we use the following parameter transformations:

π_{u} = \frac{1 + \tanh α_{u}}{2}, and t_{i} = \frac{1 + \cos θ_{i}}{2},

where tanh is the hyperbolic tangent function. Suppose (x_u, y_u) is the central coordinate of the exposed area on the coded aperture associated to the uth source, r_u is the radius of the exposed area, and (x_i, y_i) is the coordinate of the ith pixel on the coded aperture. If

{(x_{u} - x_{i})}^{2} + {(y_{u} - y_{i})}^{2} < r_{u}^{2}

, then the pixel (x_i, y_i) is within the exposed area and λ_i = 1, otherwise λ_i = 0. Thus, we can calculate λ_i as

λ_{i} = Γ [r_{u}^{2} - {(x_{u} - x_{i})}^{2} - {(y_{u} - y_{i})}^{2}],

where Γ(x) is the hard threshold function, if x > 0 then Γ(x) = 1, otherwise Γ(x) = 0. In order to make the cost function differentiable, we use the sigmoid function to replace the hard threshold function as follows:

\begin{matrix} λ_{i} = sigmoid (x_{i}, y_{i}) \\ = \frac{1}{1 + \exp {- a_{λ} [- {(x_{u} - x_{i})}^{2} - {(y_{u} - y_{i})}^{2} + r_{u}^{2}]}} \\ = \frac{1}{1 + \exp [a_{λ} {(x_{u} - x_{i})}^{2} + a_{λ} {(y_{u} - y_{i})}^{2} - a_{λ} r_{u}^{2}]}, \end{matrix}

where a_λ is the parameter indicating the steepness of the sigmoid function. In the following simulations, a_λ = 1 and r_u is a predefined parameter to constrain the size of the exposed area. Thus, the optimization problem is reformulated as

({\hat{α}}_{u}, {\hat{θ}}_{i}, {\hat{x}}_{u}, {\hat{y}}_{u}) = \underset{α_{u}, θ_{i}, x_{u}, y_{u}}{\arg \min} {‖ {\tilde{A}}^{H} \tilde{A} - I ‖}_{F} .

It is noted that the forward imaging model of the X-ray CT systems with fan-beam illuminations has a similar mathematical form as Eq. (7). In the past, the gradient-based algorithm has been used to optimize the coded apertures in the fan-beam X-ray CT systems [28]. Based on the idea in this paper, the proposed SCO method can be extended to fan-beam X-ray CT systems. In addition, the proposed SCO method can be also applied to other type of realistic X-ray tomosynthesis systems using a single pair of X-ray source and coded aperture that moves around the object to collect different measurements. The coded aperture could be translated or rotated to modify the coding pattern when the X-ray source is moved to different locations. However, these topics are out of the scope of this paper and will be studied and discussed in future work.

3.2. Regularization methods

In order to attain more stable and robust solutions for the optimization problem in Eq. (24), regularization methods are used to constrain the solution space and improve the optimization performance. In particular, we adopt three kinds of penalty terms called quadratic penalty, source penalty and transmittance penalty. Using the parameter transformation in Eq. (21), the definition domain of π_u and t_i is relaxed to [0, 1]. However, π_u and t_i have to be binarized after optimization to obtain binary source patterns and binary coded apertures. To reduce the quantization error, the quadratic penalty is adopted to obtain near-binary optimized source pattern and coded apertures with gray-scaled pixel values [36, 37]. The quadratic penalty term is defined as

R_{Q} = \sum_{u = 1}^{P} [1 - {(2 π_{u} - 1)}^{2}] + \sum_{i = 1}^{M P} [1 - {(2 t_{i} - 1)}^{2}] .

As the values of π_u and t_i approach 0 and 1 respectively, the penalty term R_Q will approach 0, otherwise, the value of R_Q will increase and reach the maximum penalty when π_u = t_i = 0.5. Thus, the penalty term in Eq. (25) constrains the solution spaces of π_u and t_i close to 0 and 1, thus benefiting a reduction in the quantization error.

Similarly, the source penalty is incorporated to constrain the number of X-ray sources to be constant. Assume there are P X-ray sources deployed on the source plane. Then, the source penalty is chosen as the square of the difference between the source number and the constant P, that is:

R_{S} = {[\sum_{u = 1}^{P} sigmoid (π_{u}) - P]}^{2}, where sigmoid (π_{u}) = \frac{1}{1 + \exp [- a_{s} (π_{u} - t_{s})]} .

In the above equation, a_s is the steepness index and t_s is the threshold value. In the following simulations, a_s = 20 and t_s = 0.5. In Eq. (26), the sigmoid function indicates all of the sources with intensities exceeding the threshold t_s. Thus, the source penalty enforces the number of sources equal to P.

Finally, a uniform sensing of the object has been shown to improve the reconstruction performance [27, 38]. In the design of coded apertures, the numbers of X-rays that sense all voxels are desirable to be constant to achieve uniform sensing. To this end, the transmittance penalty term is chosen as the variance of the sum of intersections of X-ray with every voxel, that is:

R_{T} = \frac{1}{N} \sum_{r = 1}^{N} {[\sum_{i = 1}^{M P} \sqrt{b_{i}} \cdot H (i, r)] - μ}^{2} = \frac{1}{N} \sum_{r = 1}^{N} {[\sum_{i = 1}^{M P} π_{u} \cdot λ_{i} \cdot t_{i} \cdot H (i, r)] - μ}^{2},

where N = N_x · N_y · N_z, and µ is the desired transmittance for each voxel. H(i, r) is the (i, r)th element in H, which means the intersection length of the ith X-ray with the rth voxel.

\sum_{i = 1}^{M P} \sqrt{b_{i}} \cdot H (i, r)

is the sum of the intersection lengths of all X-rays in the rth voxel. The second equality in Eq. (27) is held since b_i is binary, thus

\sqrt{b_{i}} = b_{i} = π_{u} \cdot λ_{i} \cdot t_{i}

.

3.3. Gradient-based optimization algorithm

Taking into account the penalty terms in Section 3.2, the SCO problem can be rewritten as following:

({\hat{α}}_{u}, {\hat{θ}}_{i}, {\hat{x}}_{u}, {\hat{y}}_{u}) = \underset{α_{u}, θ_{i} x_{u}, y_{u}}{\arg \min} (F + γ_{Q} R_{Q} + γ_{S} R_{S} + γ_{T} R_{T}),

where F, R_Q, R_S and R_T are described in Eqs. (19), (25), (26), and (27), respectively. γ_Q, γ_S and γ_T are the weight coefficients of the penalty terms. Prior to the optimization in Eq. (28), the optimization variables have to be initialized. If the uth source is turned on, α_u is initialized as ln3, otherwise α_u is initialized as ln(−3). θ_i is initialized as 4π/5 and π/5 for the blocked and unblocked pixels on the coded apertures, respectively. In the following simulations, the source pattern is gridded into 3 × 3 pixels, each pixel represents one source location. In the initialization stage, the central source pixel and the four source pixels at corners are turned on. During the optimization, there are always 5 sources turned on. For the exposed area of all X-ray sources, the initial origin coordinate (x_u, y_u) is set right in the center of the coded aperture. The radius of the exposed area spanned by the cone-beam projection is set to be one fourth of the side length of coded aperture. The size of the source location grid as well as the radius of exposed area per source can be both modified without loss of generality.

After the initialization, the steepest descent algorithm is used to solve for the optimization problem. Let ω represent any optimization variable α_u, θ_i, x_u, or y_u. In the kth iteration, the variable is updated as

ω^{k} = ω^{k - 1} - step \cdot (\nabla F + γ_{Q} \nabla R_{Q} + γ_{S} \nabla R_{S} + γ_{T} \nabla R_{T}),

where ∇F, ∇R_Q, ∇R_S and ∇R_T are the gradients of the cost function and penalty terms with respect to the variable ω^k−¹. The derivation of ∇F is provided in Appendix B, and the derivations of ∇R_Q, ∇R_S and ∇R_T are provided in Appendix C. After the optimization algorithm is terminated, the variables are quantized to obtain the binary source pattern and binary coded apertures.

4. Simulation and analysis

This section illustrates the simulation results of the joint optimization algorithm, and compares it to the approaches where the source orientations and coded aperture structures are optimized individually. The simulations use a compressive X-ray tomosynthesis configuration with a flat 2D detector plane with N₁ × N₂ = 64 × 64 elements. Above the detector, 9 cone-beam X-ray sources (P = 9) are placed uniformly in a 3 × 3 array, where only 5 sources are turned on as shown in Fig. 1. A 64 × 64 coded aperture is placed under each X-ray source, where the pitch of the coded aperture is adjusted to obtain one-to-one correspondence with the detector elements. An object under detection is placed between the sources and detector with dimension of N_x × N_y × N_z = 32 × 32 × 4. The ASTRA Tomography Toolbox (“All Scale Tomographic Reconstruction Antwerp”) [39] is used to obtain the system matrices H_p and the projection measurements y_p for all of the X-ray sources. The weight coefficients of all penalty terms are set to be 1 × 10⁻⁴.

Figure 4 illustrates the simulation results based on the conditions mentioned above. From left to right, it shows the images of the four layers in the object along the z-axis. The colors from black to white represent the pixel values from 0 to 1.25. The first row shows the original images of the object, which are generated by the 3D Shepp-Logan phantom [40]. The second row illustrates the reconstructed images with the initial simulation settings without any optimization. In particular, the conditions of initial source pattern, source orientation and coded apertures have been described in Section 3.3.

Fig. 4 The simulation results of different optimization methods based on the proposed compressive tomosynthesis system. From left to right, are the images of the four layers in the object along the z-axis. From top to bottom, are the original images, reconstructed images using the initial settings, reconstructed images under the individual optimization of the source orientation, individual optimization of the coded apertures, and the joint optimization of the source and coded apertures.

Download Full Size | PDF

The GPSR is used to reconstruct the object by solving for the ill-posed problem in Eq. (7) aided by sparsity regularization [33]. The sparse basis Ψ in Eq. (7) is defined as the Kronecker product of a two-dimensional wavelet transform basis in the x-y plane and a one-dimensional discrete Fourier transform (DCT) basis along the z-axis [41]. Hereafter, the peak signal-to-noise ratio (PSNR) is used as the metric to evaluate the reconstruction results as it does not depend strongly on the image intensity scaling. Given the original image I ∈ ℝ^N×N and its reconstruction R ∈ ℝ^N×N , the PSNR is defined as $PSNR = 10 \log_{10} (I_{\max}^{2} / MSE)$ , where I_max is the maximum possible pixel value within the image I, and MSE is the mean square error of the reconstructed image defined as $MSE = \frac{1}{N^{2}} \sum_{i = 1}^{N} \sum_{j = 1}^{N} {[I (i, j) - R (i, j)]}^{2}$ . Using the original images as the benchmark, the PSNR of all reconstructed images are presented below, and the average PSNR of the reconstructed images without optimization is 24.92dB. In order to compare the details of the reconstruction images from different algorithms, Fig. 5 illustrates the magnified regions within the red boxes in Fig. 4.

Fig. 5 Magnified regions within the red boxes in Fig. 4.

Download Full Size | PDF

The third rows in Figs. 4 and 5 illustrate the reconstructed images under the individual optimization of source orientations, where the source pattern and coded apertures are not changed. In the optimization process, the step size to update the orientation variables is 0.01. As mentioned above, optimizing source orientations is equivalent to optimizing the locations of the exposed areas on the coded apertures. Figure 6 shows the exposed areas before and after optimizing the source orientations. Figure 6(a) gives the initial exposed areas corresponding to the 9 X-ray sources. After optimization, the exposed areas are changed from Fig. 6(a) to the positions in Fig. 6(b). The difference between the initial and optimized exposed areas is shown in Fig. 6(c). Notice that in this simulation, only the X-ray sources in the center and at the four corners are turned on. Thus, only the exposed areas associated with the opened sources are modified. From Fig. 4, it is observed that the optimization on the source orientation will improve the reconstruction PSNR to 26.20 on average. It is also worth noting that when the exposed areas are magnified, the impact of the source orientations will be diminished.

Fig. 6 Individual optimization of source orientations: (a) initial exposed areas on the coded apertures, (b) optimized exposed areas on the coded apertures, and (c) the difference between the initial and optimized exposed areas.

Download Full Size | PDF

The fourth rows in Figs. 4 and 5 illustrate the reconstructed images under the individual optimization of coded apertures by fixing the source pattern and orientations. The step length to update coded apertures is 0.1. The code optimization seems more powerful than the orientation optimization in improving the reconstruction accuracy. That is because the degree of optimization freedom in coded apertures is much higher than that of the source orientations. Figure 7 shows the coded apertures before and after optimization, as well as their differences. Again, only the coded apertures associated with the opened sources are modified, while other four coded apertures are fixed.

Fig. 7 Individual optimization of coded apertures: (a) initial coded apertures, (b) optimized coded apertures, and the (c) the difference between the initial and optimized coded apertures.

Download Full Size | PDF

The fifth rows in Figs. 4 and 5 illustrate the reconstructed images under the proposed joint optimization method. In this simulation, the source pattern, source orientations and coded apertures are simultaneously optimized, which further increases the degree of optimization freedom. The step lengths to update the source pattern, source orientations and coded apertures are 0.01, 0.01 and 0.1, respectively. Figures 8(a) and 8(b) show the change of the source pattern, where white and black circles represent the opened and closed X-ray sources, respectively. Figures 8(c) and 8(d) illustrate the initial exposed areas and the optimized exposed areas due to the change of source orientations. Figure 8(e) is the difference between the initial and optimized exposed areas. Figures 8(f)-8(h) show the initial and optimized coded apertures, as well as their differences. Compared to the individual optimization on the source orientations or coded apertures, SCO approach can improve the average PSNR of reconstruction images by 6.3dB and 4.5dB, respectively. Figure 9 illustrates the error patterns of the magnified regions within the red boxes in Fig. 4. The error pattern is defined as the normalized absolute difference between the reconstructed image and the original image. From blue to red colors represent the range of [0,1]. It is shown that the SCO approach leads to less pattern errors than the reconstruction methods based on initial setting, orientation optimization or coded aperture optimization.

Fig. 8 Joint optimization of source and coded apertures. Top to bottom rows show the source patterns, source orientations and coded apertures before and after optimization.

Download Full Size | PDF

Fig. 9 Error patterns of the magnified regions within the red boxes in Fig. 4.

Download Full Size | PDF

Figure 10 compares the convergence properties of different optimization methods mentioned above. Figure 10(a) shows the convergence curves of the gray-scaled cost functions calculated using the optimization variables without quantization. On the other hand, Fig. 10(b) shows the convergence curves of the binarized cost functions, which are calculated using the optimization variables after quantization. That means the coded apertures and source patterns have to be binary. It is found that the orientation optimization oscillates more than others probably because of the local minima. In contrast, the coded aperture optimization and joint optimization converge much more smoothly. Due to the length limitation of the paper, we will provide experimental results of this work based on other object data sets in our future work.

Fig. 10 The convergence curves of (a) the gray-scaled cost function and (b) the binarized cost function.

Download Full Size | PDF

In previous works, the singular value decomposition (SVD) analysis has been used to compare the conditions of the tomosynthesis or tomography sensing matrix obtained by different optimization algorithms [27, 42]. Figure 11 illustrates the first 850 normalized singular values of the sensing matrices obtained by the individual optimization on source orientations (black dotted), individual optimization on coded apertures (red dashed), and the joint optimization method (green solid). Here all of the singular values are normalized by the largest one. It proves that better measurement strategy should have more SVD values lying beyond a given noise level, since it will capture more orthogonal components of the underlying object. It is observed that the proposed joint optimization approach outperforms other methods based on this criterion. In addition, the condition number is used to measure how ill-conditioned the system matrix is [43]. Here, the condition number κ is defined as the ratio between the greatest singular value and the least nonzero singular value in Fig. 11. It can be seen that the condition number for the joint optimization method (κ = 3.6140) is less ill-conditioned compared to the individual coded aperture optimization (κ = 4.5296) or individual source orientation optimization (κ = 5.2999). That means the joint optimization method further improves the condition of the tomosynthesis sensing matrix, thus benefits to reduce the reconstruction artifacts.

Fig. 11 Normalized singular values of the tomosynthesis system matrix after the individual optimization on source orientations (black dotted), individual optimization on coded apertures (red dashed), and the joint optimization on source and coded apertures (green solid).

Download Full Size | PDF

5. Conclusion

A source and coded aperture joint optimization method is developed for X-ray tomosynthesis system to effectively improve the reconstruction performance of 3D objects. The forward imaging model of the proposed compressive tomosynthesis system is described. Then, the proposed SCO problem is formulated according to the RIP condition of the sensing matrix. Based on parameter transformations, the SCO problem is relaxed to a continuous optimization problem, and then a gradient-based algorithm is used to iteratively access the optimal combination of source locations, source orientations and coded apertures. A set of regularization constraints are used to reduce the quantization error and obtain uniform sensing. The proposed SCO method proves to be more powerful in improving the reconstruction performance over either coded aperture structure, or angle of incidence individual optimization methods. In the future, the proposed SCO method will be generalized to different source geometries with larger source grids, as well as fan-beam X-ray CT systems or X-ray tomosynthesis systems with different sensing geometries. In addition, the proposed SCO method will be verified based on an experimental X-ray tomosynthesis testbed.

A. Appendix A: derivation of Eq. (19)

According to Eq. (14), we have

F = \sum_{m = 1}^{N} \sum_{n = 1}^{N} {| {[{\tilde{A}}^{T} \tilde{A}]}_{m, n} - I_{m, n} |}^{2} = \sum_{m \neq n} {({[{\tilde{A}}^{T} \tilde{A}]}_{m, n} - 0)}^{2} + \sum_{m = n} {({[{\tilde{A}}^{T} \tilde{A}]}_{m, n} - 1)}^{2} .

Since the columns of $\tilde{A}$ are normalized, ${[{\tilde{A}}^{T} \tilde{A}]}_{m, n} = 1$ if m = n. Thus, the second term in Eq. (30) is equal to 0. In addition, considering ${\tilde{A}}^{T} \tilde{A}$ is a symmetric matrix, minimizing the cost function in Eq. (30) is equivalent to minimizing the cost function in the following

\begin{array}{l} F = \sum_{m < n} {({[{\tilde{A}}^{T} \tilde{A}]}_{m, n})}^{2} \\ = {\sum_{m < n} {\frac{\sum_{i, j = 1}^{M P} 〈 c_{i}, c_{j} 〉 w_{i} (m) w_{j} (n)}{{[\sum_{i, j = 1}^{M P} 〈 c_{i}, c_{j} 〉 w_{i} (m) w_{j} (m)]}^{1 / 2} {[\sum_{i, j = 1}^{M P} 〈 c_{i}, c_{j} 〉 w_{i} (n) w_{j} (n)]}^{1 / 2}}}}^{2} \\ = \sum_{m < n} {\frac{d (m, n)}{{[d (m, m)]}^{1 / 2} {[d (n, n)]}^{1 / 2}}}^{2}, \end{array}

where

d (m, n) = \sum_{i = 1}^{M P} b_{i} R_{i}^{m n} and R_{i}^{m n} = w_{i} (m) w_{i} (n) .

B. Appendix B: derivation of gradient of cost function

From Eqs. (19) and (20), we have

F = \sum_{m < n} {[\frac{\sum_{i = 1}^{M P} b_{i} R_{i}^{m n}}{{(\sum_{i = 1}^{M P} b_{i} R_{i}^{m m})}^{1 / 2} {(\sum_{i = 1}^{M P} b_{i} R_{i}^{n n})}^{1 / 2}}]}^{2} = \sum_{m < n} \frac{{(\sum_{i = 1}^{M P} b_{i} R_{i}^{m n})}^{2}}{(\sum_{i = 1}^{M P} b_{i} R_{i}^{m m}) (\sum_{i = 1}^{M P} b_{i} R_{i}^{n n})} .

Then, we can calculate the partial derivative of F to b_i as following:

\begin{matrix} \frac{\partial F}{\partial b_{i}} = \sum_{m < n} \frac{2 d (m, n) R_{i}^{m n} d (m, m) d (n, n) - d {(m, n)}^{2} [R_{i}^{m m} d (n, n) + R_{i}^{n n} (m, m)]}{d {(m, m)}^{2} d {(n, n)}^{2}} \\ = \sum_{m < n} [\frac{2 d (m, n) R_{i}^{m n}}{d (m, m) d (n, n)} - \frac{d (m, n) R_{i}^{m m}}{d {(m, m)}^{2} d (n, n)} - \frac{d {(m, n)}^{2} R_{i}^{n n}}{d (m, m) d {(n, n)}^{2}}] . \end{matrix}

Since i = u · M + v, b_i can be also written as b_uv. Then, the derivative of F to the optimization variables can be calculated as follows:

\frac{\partial F}{\partial α_{u}} = \sum_{v = 1}^{M} \frac{\partial F}{\partial b_{u v}} \cdot \frac{\partial b_{u v}}{\partial α_{u}} = \frac{1}{2} (1 - \tanh^{2} α_{u}) \cdot \sum_{v = 1}^{M} [\frac{\partial F}{\partial b_{u v}} \cdot λ_{i} \cdot t_{i}],

\frac{\partial F}{\partial θ_{i}} = \frac{\partial F}{\partial b_{i}} \cdot \frac{\partial b_{i}}{\partial θ_{i}} = - \frac{1}{2} \sin θ_{i} \cdot [\frac{\partial F}{\partial b_{i}} \cdot π_{u} \cdot λ_{i}],

\frac{\partial F}{\partial x_{u}} = \sum_{v = 1}^{M} \frac{\partial F}{\partial b_{u v}} \cdot \frac{\partial b_{u v}}{\partial x_{u}} = - 2 a \sum_{v = 1}^{M} [\frac{\partial F}{\partial b_{u v}} \cdot π_{u} \cdot t_{i} \cdot λ_{i} \cdot (1 - λ_{i}) \cdot (x_{u} - x_{i})],

and

\frac{\partial F}{\partial y_{u}} = \sum_{v = 1}^{M} \frac{\partial F}{\partial b_{u v}} \cdot \frac{\partial b_{u v}}{\partial y_{u}} = - 2 a \sum_{v = 1}^{M} [\frac{\partial F}{\partial b_{u v}} \cdot π_{u} \cdot t_{i} \cdot λ_{i} \cdot (1 - λ_{i}) \cdot (y_{u} - y_{i})] .

C. Appendix C: derivation of gradients of penalty terms

The derivatives of the quadratic penalty to the optimization variables are:

\frac{\partial R_{Q}}{\partial α_{u}} = (1 - \tanh^{2} α_{u}) \cdot (2 - 4 π_{u}),

\frac{\partial R_{Q}}{\partial θ_{i}} = - \frac{1}{2} \sin θ_{i} \cdot (- 8 t_{i} + 4) = \sin θ_{i} \cdot (4 t_{i} - 2), and \frac{\partial R_{Q}}{\partial x_{u}} = \frac{\partial R_{Q}}{\partial y_{u}} = 0 .

The derivatives of the source penalty to the optimization variables are:

\frac{\partial R_{S}}{\partial α_{u}} = a_{s} [\sum_{u = 1}^{P} sigmoid (π_{u}) - P] \cdot sigmoid (π_{u}) \cdot [1 - sigmoid (π_{u})] \cdot (1 - \tanh^{2} α_{u}),

and ∂R_S/∂θ_i = ∂R_S/∂x_u = ∂R_S/∂y_u = 0. The derivatives of the transmittance penalty to the optimization variables are derived as following:

\frac{\partial R_{T}}{\partial π_{u}} = \frac{2}{N} \sum_{r = 1}^{N} {[(\sum_{i = 1}^{M P} π_{u} \cdot λ_{i} \cdot t_{i} \cdot H (i, r)) - μ] \cdot (\sum_{v = 1}^{M} λ_{i} \cdot t_{i} \cdot H (i, r))} .

Thus, we have

\frac{\partial R_{T}}{\partial π_{u}} = \frac{1 - \tanh^{2} α_{u}}{N} \cdot \sum_{r = 1}^{N} {[(\sum_{i = 1}^{M P} π_{u} \cdot λ_{i} \cdot t_{i} \cdot H (i, r)) - μ] \cdot (\sum_{v = 1}^{M} λ_{i} \cdot t_{i} \cdot H (i, r))} .

In addition,

\frac{\partial R_{T}}{\partial t_{i}} = \frac{2}{N} \sum_{r = 1}^{N} {[(\sum_{i = 1}^{M P} π_{u} \cdot λ_{i} \cdot t_{i} \cdot H (i, r)) - μ] \cdot [π_{u} \cdot λ_{i} \cdot H (i, r)]} .

Thus, we have

\frac{\partial R_{T}}{\partial θ_{i}} = \frac{\sin θ_{i}}{N} \sum_{r = 1}^{N} {[(\sum_{i = 1}^{M P} π_{u} \cdot λ_{i} \cdot t_{i} \cdot H (i, r)) - μ] \cdot [π_{u} \cdot λ_{i} \cdot H (i, r)]} .

Finally,

\frac{\partial R_{T}}{\partial λ_{i}} = \frac{2}{N} \sum_{r = 1}^{N} {[(\sum_{i = 1}^{M P} π_{u} \cdot λ_{i} \cdot t_{i} \cdot H (i, r)) - μ] \cdot [π_{u} \cdot t_{i} \cdot H (i, r)]} .

Thus, we have

\frac{\partial R_{T}}{\partial x_{u}} = - \frac{4 a}{N} \sum_{r = 1}^{N} {[(\sum_{i = 1}^{M P} π_{u} \cdot λ_{i} \cdot t_{i} \cdot H (i, r)) - μ] \cdot \sum_{v = 1}^{M} [π_{u} \cdot t_{i} \cdot H (i, r)] \cdot λ_{i} \cdot (1 - λ_{i}) \cdot (x_{u} - x_{i})},

and

\frac{\partial R_{T}}{\partial y_{u}} = - \frac{4 a}{N} \sum_{r = 1}^{N} {[(\sum_{i = 1}^{M P} π_{u} \cdot λ_{i} \cdot t_{i} \cdot H (i, r)) - μ] \cdot \sum_{v = 1}^{M} [π_{u} \cdot t_{i} \cdot H (i, r)] \cdot λ_{i} \cdot (1 - λ_{i}) \cdot (y_{u} - x_{i})} .

Funding

National Science Foundation (NSF) (CIF 1717578); Nokia Foundation and Fulbright Finland under the 2017 Fulbright-Nokia Distinguished Chair Award; Fundamental Research Funds for the Central Universities (2018CX01025); Chinese Scholarship Council (201706035012, 201706840090).

References

1. A. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (Society for Industrial and Applied Mathematics, 2011).

2. J. T. Dobbins III and D. J. Godfrey, “Digital x-ray tomosynthesis: current state of the art and clinical potential,” Phys. Med. Biol. 48, R65–R106 (2003). [CrossRef] [PubMed]

3. F. Elarnaut, J. P. O. Evans, D. Downes, A. J. Dicken, S. X. Godber, and K. D. Rogers, “Sporadic absorption tomography using a conical shell x-ray beam,” Opt. Express 25, 33029–33042 (2017). [CrossRef]

4. H. Gao, L. Zhang, Z. Chen, Y. Xing, J. Cheng, and Y. Yang, “Application of x-ray ct to liquid security inspection: system analysis and beam hardening correction,” Nucl. Instrum. Meth. A 579, 395–399 (2007). [CrossRef]

5. H. Gao, L. Zhang, Z. Chen, Y. Xing, H. Xue, and J. Cheng, “Straight-line-trajectory-based x-ray tomographic imaging for security inspections: system design, image reconstruction and preliminary results,” IEEE Trans. Nucl. Sci. 60, 3955–3968 (2013). [CrossRef]

6. T. M. Buzug, Computed Tomography: from Photon Statistics to Modern Cone-Beam CT (Springer, 2008).

7. F. Natterer, The Mathematics of Computerized Tomography (Society for Industrial and Applied Mathematics, 1986).

8. I. Reiser and S. Glick, Tomosynthesis Imaging (Taylor and Francis, 2014).

9. J. S. Jørgensen, E. Y. Sidky, and X. Pan, “Quantifying admissible undersampling for sparsity-exploiting iterative image reconstruction in x-ray ct,” IEEE Trans. Med. Imag. 32, 460–473 (2013). [CrossRef]

10. V. Krishnaswamy, K. E. Michaelsen, B. W. Pogue, S. P. Poplack, I. Shaw, K. Defrietas, K. Brooks, and K. D. Paulsen, “A digital x-ray tomosynthesis coupled near infrared spectral tomography system for dual-modality breast imaging,” Opt. Express 20, 19125–19136 (2012). [CrossRef] [PubMed]

11. H. Kudo, T. Suzuki, and E. A. Rashed, “Image reconstruction for sparse-view ct and interior ct-introduction to compressed sensing and differentiated backprojection,” Quant. Imag. Med. Surg. 3, 147–161 (2013).

12. R. Smith-Bindman, J. Lipson, R. Marcus, K. Kim, M. Mahesh, R. Gould, A. González, and D. L. Miglioretti, “Radiation dose associated with common computed tomography examinations and the associated lifetime attributable risk of cancer,” Arch. Internal Med. 169, 2078–2086 (2009). [CrossRef]

13. W. Hou and C. Zhang, “Analysis of compressed sensing based ct reconstruction with low radiation,” in Proc. ISPACS, 291–296 (2014).

14. A. A. Samarskii and P. N. Vabishchevich, Numerical Methods for Solving Inverse Problems of Mathematical Physics (Walter de Gruyter, 2007).

15. K. Hämäläinen, A. Kallonen, V. Kolehmainen, M. Lassas, K. Niinimäki, and S. Siltanen, “Sparse tomography,” SIAM J. Sci. Comput. 35, 644–665 (2013). [CrossRef]

16. X. Cui, L. Mili, and H. Yu,“Sparse-prior-basedprojectiondistanceoptimizationmethodforjointct-mrireconstruction,” IEEE Access 5, 20099–20110 (2017). [CrossRef]

17. E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Phys. Med. Biol. 53, 4777–4807 (2008). [CrossRef] [PubMed]

18. Y. Kaganovsky, D. Li, A. Holmgren, H. Jeon, K. P. MacCabe, D. G. Politte, J. A. O’Sullivan, L. Carin, and D. J. Brady, “Compressed sampling strategies for tomography,” J. Opt. Soc. Am. A 31, 1369–1394 (2014). [CrossRef]

19. D. J. Brady, A. Mrozack, K. Maccabe, and P. Llull, “Compressive tomography,” Adv. Opt. Photonics 7, 756–813 (2015). [CrossRef]

20. K. Choi and D. J. Brady, “Coded aperture computed tomography,” in Proc. SPIE7468, 74680B (2009).

21. A.A.M. Muñoz, A. Vella, M.J.F. Healy, D.W. Lane, I. Jupp, and D. Lockley,“ 3D-printedcodedaperturesforx-ray backscatter radiography,” in Proc. SPIE10393, 103930F (2017).

22. A. Holmgren, “Coding strategies for x-ray tomography,” Dissertation, Duke University (2016).

23. Wolfmet, “Precision tungsten collimators from wolmet,” https://static.mimaterials.com/SLMflyerweb.pdf.

24. X. Yang, R. Hofmann, R. Dapp, T. Kamp, T. Rolo, X. Xiao, J. Moosmann, J. Kashef, and R. Stotzka, “TV-based conjugate gradient method and discrete l-curve for few-view ct reconstruction of x-ray in vivo data,” Opt. Express 23, 5368–5387 (2015). [CrossRef] [PubMed]

25. M. A. Davenport, M. F. Duarte, Y. C. Eldar, and G. Kutyniok, “Introduction to compressed sensing,” in Compressed Sensing: Theory and Applications (Cambridge University, 2011).

26. M. Elad, “Optimized projections for compressed sensing,” IEEE Trans. Signal Process. 55, 5695–5702 (2007). [CrossRef]

27. A. P. Cuadros, C. Peitsch, H. Arguello, and G. R. Arce, “Coded aperture optimization for compressive x-ray tomosynthesis,” Opt. Express 23, 32788–32802 (2015). [CrossRef] [PubMed]

28. A. P. Cuadros and G. R. Arce, “Coded aperture optimization in compressive x-ray tomography: a gradient descent approach,” Opt. Express 25, 23833–23849 (2017). [CrossRef] [PubMed]

29. E. Candés, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52, 489–509 (2006). [CrossRef]

30. D. L. Donoho, “Compressive sensing,” IEEE Trans. Inf. Theory 52, 1289–1306 (2006). [CrossRef]

31. H. Arguello and G. R. Arce, “Colored coded aperture design by concentration of measure in compressive spectral imaging,” IEEE Trans. Image Process. 23, 1896–1908 (2014). [CrossRef] [PubMed]

32. J. M. Duarte-Carvajalino and G. Sapiro, “Learning to sense sparse signals: simultaneous sensing matrix and sparsifying dictionary optimization,” IEEE Trans. Image Process. 18, 1395–1408 (2009). [CrossRef] [PubMed]

33. M. Figueiredo, R. D. Nowak, and S. J. Wright, “Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems,” IEEE J. Sel. Topics Signal Process. 1, 586–597 (2007). [CrossRef]

34. E. Candès and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25, 21–30 (2008). [CrossRef]

35. H. Rauhut, “Compressive sensing and structured random matrices,” Radon Series Comp. Appl. Math XX, 1–94 (2011).

36. X. Ma and G. R. Arce, “Binary mask optimization for inverse lithography with partially coherent illumination,” J. Opt. Soc. Am. A 25, 2960–2970 (2008). [CrossRef]

37. X. Ma and G. R. Arce, Computational Lithography (John Wiley & Sons, Inc, 2010). [CrossRef]

38. A. P. Cuadros, K. Wang, C. Peitch, H. Arguello, and G. R. Arce, “Coded aperture design for compressive x-ray tomosynthesis,” in Imaging and Applied Optics, OSA, CW2F.2 (2015).

39. W. Xu, F. Xu, M. Jones, B. Keszthelyi, J. Sedat, D. Agard, and K. Mueller, “High-performance iterative electron tomography reconstruction with long-object compensation using graphics processing units,” J. Struct. Biol. 171, 142–153 (2010). [CrossRef] [PubMed]

40. V. Cunha, C. Catherine, and M. Schabel, “3D shepp-logan phantom,” http://www.mathworks.com/matlabcentral/fileexchange/9416-3d-shepp-logan-phantom (2006).

41. G. R. Arce, D. J. Brady, L. Carin, H. Arguello, and D. S. Kittle, “Compressive coded aperture spectral imaging: an introduction,” IEEE Signal Process. Mag. 31, 105–115 (2014). [CrossRef]

42. T. Mao, A. P. Cuadros, X. Ma, W. He, Q. Chen, and G. R. Arce, “Fast optimization of coded apertures in x-ray computed tomography,” Opt. Express 26, 24461–24478 (2018). [CrossRef] [PubMed]

43. D. J. Brady, Optical Imaging and Spectroscopy (Wiley, 2009). [CrossRef]

Source and coded aperture joint optimization for compressive X-ray tomosynthesis

Abstract

1. Introduction

2. Forward imaging model of X-ray tomosynthesis

3. Joint optimization of source and coded aperture

3.1. Cost function of the SCO problem

3.2. Regularization methods

3.3. Gradient-based optimization algorithm

4. Simulation and analysis

5. Conclusion

A. Appendix A: derivation of Eq. (19)

B. Appendix B: derivation of gradient of cost function

C. Appendix C: derivation of gradients of penalty terms

Funding

References

Cited By

Figures (11)

Equations (48)

Optics Express