1. Introduction

X-ray tomography has become essential in medical imaging, airport security, and industrial testing. The filtered back-projection algorithm (FBP) is the algorithm of choice for image reconstruction given its speed. However, irregular geometries of scanners and missing projection data produce artifacts and noise that make the reconstructions obtained with the FBP inadequate for diagnosis [1]. Furthermore, the numerous projections needed to obtain a high-quality reconstruction in medical X-ray computed tomography (CT) lead to high levels of radiation dose; thus, strategies to reduce the number of measurements are critical. Image estimation from incomplete data has been extensively researched since the introduction of CT [2]. Sparsity-promoting and total variation regularization algorithms have been recently used to alleviate the ill-posed inverse problem, obtaining better image reconstructions [3]. Most strategies, however, rely on optimizing the hardware settings by lowering the number of angles at which projections are taken [4]. These approaches, however, fall short when the number of measurements is further reduced leading to artifacts in the reconstructions. Brady et al. provide a study on how structured illumination can overcome these limitations while potentially reducing radiation exposure, decreasing the required number of detectors, and the acquisition time [5]. Particularly, the latter is achieved when using coded apertures in conjunction with multiplexed illumination, that is, the simultaneous illumination of the object using multiple sources. The acquisition time is reduced given that projections from more than one source may be acquired in a single step. Tomographic imaging using structured illumination was first introduced by Smith et al. in the 1980’s [6]. However, it was not further studied until the introduction of compressive sensing (CS) algorithms, which provided a strong mathematical framework to solve the resulting ill-posed problem. For instance, K. Choi et al. introduced coded aperture X-ray tomosynthesis, which outperforms sparse regularization due to the direct acquisition of compressive measurements [7]. Similarly, in [8], Kaganovsky et al. introduced subsampling strategies for transmission tomography in medical CT scanner geometries. These systems require coded apertures placed either, in front of the X-ray source to create structured X-ray bundles that interrogate the object, or at the detector to block the X-ray projections in some fashion. Interestingly, it was shown in [8] that randomly subsampling the detectors in each view angle presents the reconstructions with the lowest error. However, such strategy for coded aperture design relied on random coded apertures, which is suboptimal, as it does not exploit the structure of the sensing matrix of the system. Cuadros et al. proposed an optimization approach to design the structure of the coded apertures in a multiple source compressive X-ray tomosynthesis imaging system based on a uniform energy criteria on the voxels and detector elements, so that the object was uniformly sensed and the elements of the detector plane uniformly sensed the information [9, 10]. This approach for code optimization, however, cannot be applied directly to compressive X-ray CT in circular geometries since the characteristics of the underlying sensing matrix are radically different. Furthermore, the mathematical analysis of uniform sensing with respect to the restricted isometry property (RIP) and the coherence of the sensing matrix was not established. Besides the coded aperture optimization work for tomosynthesis in [9, 10], there has been only one attempt to design coded apertures for conventional X-ray CT systems. This approach introduced by Mojica et al. in [11] aims to optimize the coded apertures for a multiple shot CT system, where the coded apertures are modified multiple times per each view angle. Their formulation is based on reducing the mutual coherence of a high-resolution CT system matrix whose dimensions are constrained to be squared and invertible, which in most practical CT systems is not realistic [1]. Furthermore, their analysis aims at reducing the coherence of the CT system matrix without taking into account the representation basis of the object of interest, thus the coherence being minimized is not that of the sensing matrix used in the image reconstruction.

The contribution of this paper is the introduction of a rigorous coded aperture optimization framework for compressive X-ray CT. The proposed coded aperture design aims to minimize the coherence of the sensing matrix, as this parameter is often used as a metric to describe the effectiveness of compressed sensing projections [12, 13]. Furthermore, the coherence of the sensing matrix accounts for the sparse basis representation of the object. The design relies on the steepest descent algorithm in order to achieve incoherence for successful CT reconstruction from incomplete projection data. The cost function is based on the point spread function (PSF) of the system, introduced by Lustig et al. [14] as a measure of the sensing matrix quality based on the RIP and mutual coherence of the sensing matrix. Additionally, two regularization functions are introduced to obtain binary coded apertures and control the transmittance of the system while sensing all the voxels of the object uniformly. Given the set of coded measurements, CS reconstruction algorithms are used to recover the image of interest. The methods presented in this paper for compressive X-ray CT can be extended to other scanning geometries such as, projection radiography, cone-beam CT, and helical CT. Additionally, by adapting the cost function to the corresponding system matrix, the method can be extended to other imaging systems based on the X-ray transform. Such as, coded aperture snapshot spectral imaging (CASSI) [15], compressive X-ray tomosynthesis and spectral computed tomography (SCT).

This paper is organized as follows. In Section 2, the mathematical model of the coded aperture compressive X-ray CT system is proposed. Section 3 introduces the gradient descent optimization algorithm proposed for the design of the coded apertures as well as the regularization constraints proper of the system. In Section 4, a simulation study for fan beam compressive X-ray CT is presented to evaluate the performance of the optimized codes obtained using the proposed approach and compare them to the best random subsampling strategy introduced in [8]. The reconstructed images obtained from optimized coded aperture projections present higher peak-signal-to-noise ratio (PSNR) than the reconstructions obtained using random coded apertures. Additionally, reconstructions using real X-ray projection data are shown for both optimized and random coded apertures. The normalized absolute error of these reconstructions with respect to the reconstruction obtained when using the complete set of projections (no coded apertures) reveals that the proposed optimization framework leads to fewer artifacts in the reconstructed image. Finally, Section 5 presents the conclusions and future work.

2. Forward projection model

The X-ray transmission imaging model is given by the Beer-Lambert law for monochromatic X-rays $I=I_0\cdot e^{-{\displaystyle {\int }_0^{\infty }\mu }(\ell )d\ell }$, where I₀ is the intensity of a particular X-ray originated from the X-ray source passing through the object, I is the measured intensity in the detector, and μ(ℓ) is the linear attenuation coefficient at location ℓ [16]. If such X-ray source is located at position $\overrightarrow{s}$ and illuminates an object in direction $\widehat{\omega }$, the data function for the imaging model is given by: $y(\overrightarrow{s},\widehat{\omega })=-\ln (I/I_0)$. Thus, the Beer-Lambert law can be rewritten as $y(\overrightarrow{s},\widehat{\omega })={\displaystyle {\int }_0^{\infty }f(\overrightarrow{s}+\ell \widehat{\omega })d\ell }$, where f corresponds to the two-dimensional linear attenuation coefficient map. This continuous-to-continuous imaging model is known as the X-ray transform of the object. For a fan-beam X-ray CT system, the source location $\overrightarrow{s}$ varies continuously in a circular trajectory around the object, and the ray direction $\widehat{\omega }$ varies at each $\overrightarrow{s}$ in the plane of the source trajectory [17]. In real systems, only a discrete number of measurements are available. Furthermore, in order to use iterative image reconstruction algorithms to solve the inverse problem, the imaging model needs to be discretized. The discrete-to-discrete formulation of the problem leads to a finite linear system of equations of the form y = Hf, where $\mathbf{f}\in {ℝ}^{N^2}$ is a vectorized form of the attenuation coefficients of the N × N object to be reconstructed and each of the elements in $\mathbf{H}\in {ℝ}^{MP\times N^2}$ represents the intersection length of each ray with every pixel inside the image, with M the number of detector elements in the linear array for each of the P view angles. Specifically, the row index of the matrix H corresponds to a line integral and the column index of the matrix H corresponds to an image pixel [18]. The number of X-ray projections needed to reconstruct an image is determined by the Shannon-Nyquist sampling theorem [17, 18]. However, iterative reconstruction techniques may be used in conjunction with compressive sensing to reconstruct the image from fewer samples [3, 17, 19]. This approach uses ℓ₁ regularization to recover the object after subsampling the number of view angles or detectors. Although reducing the number of projection angles has been studied due to its straightforward implementation, reducing the number of detectors produces higher quality reconstructions for equivalent amounts of radiation as shown in [8, 20]. Coded aperture compressive X-ray CT was proposed by Kaganovsky et al. in [8] as an implementation of detector subsampling. It places coded apertures in front of the fan-beam source to modulate its energy producing a particular coded projection onto the detector strip. The coded apertures have the same number of elements as the detector array and the coded aperture pitch is fixed to obtain one to one correspondence with the detector elements. Different coded aperture patterns are used for each view angle position s_p with p = 1, ⋯, P as shown in Fig. 1(a), which depicts the coded aperture compressive X-ray CT system for a fan beam architecture.

 

Fig. 1 (a) Coded aperture compressive fan beam X-ray CT system. Different coded apertures are placed in front of the rotating source for each location s_p. Each row in the system matrix H describes the sensing for a particular detector element. The unblocked coded aperture elements select the rows of the matrix H associated with the detector elements that correspond to the unblocked coded aperture elements. (b) Sensing matrix $\mathbf{Q}\in {ℝ}^D{ }^{\times N^2}$, of a system with N² = 9, P = 4, D = 6 and M = 3. Q is formed by the rows selected by the unblocked elements of the coded apertures. Each column of the coded aperture matrix C is associated with a detector element. If the coded aperture element associated with the i^th detector element is 0 then c_i = 0. Otherwise, the vector c_i is a natural basis in ${ℝ}^D$.

Download Full Size | PDF

Each row of the CT system matrix, H, indexes a particular detector element at a given view angle, from the many detector arrays at all view angles in the system. Thus, as shown in Fig. 1(a) when a coded aperture is placed in front of the source, the coded aperture elements select the rows of the matrix H associated with the detector elements that correspond to the unblocking coded aperture elements. Mathematically, the coded aperture matrix C = [c₁, c₂, ⋯, c_MP] is defined as a binary matrix, where each D-long column vector c_i with i = 1, ⋯, MP, is associated with a particular detector element. If the coded aperture element associated with the i^th detector element is 0 (blocking element), the corresponding vector c_i is a zero vector. Otherwise c_i is a natural basis in R^D. Measurements are not multiplexed in the detector array; thus, the vectors composing the matrix C are disjoint. The resulting system matrix Q is formed by the rows of the matrix H that correspond to the detectors selected by the non-zero entries of matrix C, that is $\mathbf{Q}=\mathbf{CH}\in {ℝ}^D{ }^{\times N^2}$, where D is the number of unblocked X-rays. The construction of the matrix C from a set of coded apertures is shown in Fig. 1(b). Additionally, the sparse structure of the sensing matrix Q and the system matrix H for a real CT system using random coded apertures, with transmittance 12.5 %, is depicted in Fig. 2(a). The discretized set of line integrals $\mathbf{y}\in {ℝ}^D$, referred to as the sinogram is given by y = Qf. Each row of the sinogram represents the set of measurements taken with a particular detector at a particular view angle, Fig. 2(b) shows the spatially coded sinogram obtained by sampling 50% of the data required by the Shannon-Nyquist theorem for a 64 × 64 Shepp-Logan phantom. The coded aperture pattern placed in front of the source at the position s_p will be reflected in the p^th row of the sinogram. Figure 2(b) depicts a zoom of a particular row of the coded sinogram. The black elements denote the blocking elements of the coded aperture for that particular view angle, while the other detectors reflect the information of the corresponding spatial position of the un-coded sinogram depicted in Fig. 2(c). The coded aperture optimization problem reduces to finding the best pattern of blocked elements in the sinogram so that the sensing matrix has low coherence.

 

Fig. 2 (a) Sparse structures of the system matrix Q, the coded aperture matrix C, and the complete CT system matrix H for a fan beam coded aperture compressive X-ray CT system with P = 16 view angles, M = 16 detectors per view angle and a N × N = 8 × 8 image. The coded apertures have 12.5% transmittance, that is D = 32. Non-zero entries are depicted as black pixels for visualization purposes. (b) Under-sampled sinogram of 64 × 64 Shepp-Logan Phantom with P = 128, M = 128 detectors per view angle and 50% under-sampling. A zoom of a particular row of the coded sinogram is depicted. The black elements denote the blocking elements of the coded aperture for that particular view angle, while the other detectors reflect the information of the corresponding spatial position of the complete uncoded sinogram depicted in (c).

Download Full Size | PDF

In order to reconstruct the object f from the coded sinogram, the ill-conditioned problem can be solved using CS algorithms. In general, the function f can be recovered provided the function f is sufficiently sparse in some basis Ψ, and such basis is incoherent with the sensing matrix [13]. Let f be represented by: $\mathbf{f}=\mathbf{\Psi }z\in {ℝ}^N{ }^{{ }^2}$, where z is the sparse coefficient representation of the object, and Ψ is the basis representation. The sensing function can be rewritten as: y = QΨz = Az, where $\mathbf{A}=\mathbf{Q}\Psi \in {ℝ}^D{ }^{\times N^2}$ is the sensing matrix. The sparse representation of the object f can be obtained by solving the following nonlinear optimization problem [13]

3. Coded aperture optimization

A possible approach to arranging the blocking and unblocking elements of the coded apertures is to distribute them randomly as proposed in initial models of coded aperture compressive X-ray CT. Kaganovsky et al. performed an extensive study regarding random subsampling of both the view angles and the detectors in [8]. Four subsampling strategies were considered: (i) Uniform view (UV) subsampling, in which view angles were selected uniformly and the complete set of detectors was used; (ii) Random view (RV) subsampling, in which the view angles were selected randomly and the complete set of detectors was used; (iii) Uniform detector (UD) subsampling, in which all the view angles were used while the set of detectors in each view angle was uniformly sampled; and (iv) Random detector (RD) subsampling, in which all the view angles were used while the detectors were randomly selected for each view. RD subsampling was shown to yield the best performance. The best results in [8] pointed to the random detector (RD) strategy, where the subsampling is chosen at random. The sensing matrices of compressive X-ray CT, however, are sparse and highly structured as depicted in Fig. 2(a). As such, random coded apertures are not optimal [9, 10, 21]. Even though coded aperture design based on the RIP and the coherence of the sensing matrix has been proposed in [12, 22] for other imaging systems such as the CASSI and the colored CASSI; these results do not extend directly to CT given that the sensing matrices of the systems are different. While the system matrix in the CASSI system has well-structured diagonal elements [22], CT matrices do not have a regular structure as shown in Fig. 2(a). To illustrate the advantages of optimized coded apertures over the UD and UV sampling strategies proposed in [8], Fig. 3 depicts the reconstruction of a 128 × 128 Walnut phantom Fig. 3(a), made available by Hämäläinen et al. in [23], using: (b) Optimized coded apertures obtained using the proposed algorithm, (c) Uniform subsampling of line integrals (UD strategy [8]), and (d) Uniform subsampling of view angles (UV strategy [8]). In all scenarios, the number of line integrals is fixed to D = 16429 for a fair comparison. Note the PSNR of the reconstruction obtained when using optimized coded apertures Fig. 3(b) is more than 20 decibels (dBs) higher than the PSNR obtained using a uniformly distributed set of line integrals or a uniformly distributed set of view angles, Figs. 3(c) and 3(d) respectively. As it can be seen in the highlighted region of the original image and the reconstructions in Fig. 3, the reconstructions obtained by uniformly subsampling the measurements present numerous artifacts compared to using optimized coded apertures. This example illustrates that structured illumination leads to a better performance when reconstructing the object from highly under-sampled data.

 

Fig. 3 (a) 128 × 128 Walnut phantom [23] and reconstructions for different subsampling scenarios using: (b) Optimized coded apertures obtained using the proposed algorithm, and (c) UD, and (d) UV strategies [8]. Measurements were simulated for N = 128, M = P = 256 and D = 16429. A zoom of the highlighted region is depicted to emphasize the error in the reconstructions. Figure (b) was reconstructed using the gradient projection for sparse reconstruction (GPSR) algorithm, Figs. (c) and (d) were reconstructed using a TV regularization algorithm.

Download Full Size | PDF

The reconstruction of sparse or compressible signals from a very limited number of measurements relies on properties of the sensing matrix such as the RIP. The restricted isometry constant δ_s of the system matrix with normalized columns, $\tilde{\mathbf{A}}\in R^D{ }^{\times N^2}$, is defined as the smallest δ_s such that

3.1. Regularization

The proposed global thresholding approach to obtain binary values for the coded apertures is sub-optimal, and the transmittance of the coded apertures obtained by the optimization algorithm cannot be controlled. Thus, it is necessary to incorporate prior knowledge about the solution by constraining the solution space through the addition of appropriate regularization terms to the cost function. In general, the formulation can be described as follows:

3.1.1. Quadratic penalty

To obtain near-binary gray optimized coded apertures through the optimization process, the quadratic penalty is adopted as in [28, 29]. The quadratic penalty term is defined as V₁(b) = 4b^T (1 − b), such that for each coded aperture element V₁(b_i) = 1 − (2b_i − 1)². The penalty is maximum for b_i = 0.5 and minimum for 0 and 1 as shown in Fig. 4(a), that is the penalty increases as b_i moves away from the binary region. The gradient of V₁(b) is given by ∇V₁(b) = (−8b + 4), which can be used in conjunction with (13) while carrying out the steepest-descent iterations as before.

 

Fig. 4 (a) Quadratic penalty cost function encouraging solution in b_i ∈ {0, 1} since the minimum error is obtained for binary values. (b) Transmittance and uniformity penalty cost function for M = P = 2N = 16 and μ = 2, i.e. approximately 50 % transmittance. Note the minimum error is achieved for the desired average number of rays.

Download Full Size | PDF

3.1.2. Transmittance penalty

Uniform sensing of the object has been used as criteria for coded aperture design in previous works [9, 10]. Particularly, uniformity in the number of rays that sense a particular voxel across the region of interest (ROI) is desirable. To this end, a second penalty term is incorporated to obtain uniformity in the ROI, Ω_C, denoted by the circle of radius N inscribed in the N × N image while using a determined number of blocking coded aperture elements, therefore obtaining a coded aperture with a fixed transmittance. The quadratic penalty term is chosen as the variance of the sum of the length of the rays that measure each pixel, that is:

4. Simulation results

Three scenarios are considered to analyze the performance of the optimized coded apertures. Initially, a compressive X-ray CT configuration with a flat 1D detector strip composed by M = 64 elements, a fan-beam X-ray source, P = 64 view angles, and an object of interest f represented by a Shepp-Logan phantom of N × N = 32 × 32 pixels are used. Each of the pixels in the coded aperture corresponds to a particular detector element as detailed in Fig. 1(a). Therefore, the coded apertures placed in front of each of the sources are also composed by M = 64 elements. The ASTRA Tomography Toolbox (âĂIJAll Scale Tomographic Reconstruction AntwerpâĂİ) [30] was used to obtain the system matrix H as well as the projection measurements y of the X-ray source. The setup is adjusted to match the full sampling scenario described by Jorgensen et al. in [17]. That is, R₀ = 40 cm distance from the center of rotation to the X-ray source, T₀ = 80 cm source-to-detector-center distance, detector length 41.3 cm and M = P = 2N. After obtaining the system matrix using random and optimal coded apertures and the projections for the corresponding systems, the linear attenuation coefficient image is reconstructed by solving the minimization problem in (1) using a modified version of the gradient projection for sparse reconstruction (GPSR) algorithm. The modified GPSR includes a filtering step, as proposed by Mejia et al. in [31]. The linear attenuation coefficient CT image is represented on the Haar 2D Wavelet basis.

Using the algorithm described in Section 3 optimized coded apertures are obtained for different subsampling rates adjusting the regularization constants γ₁, γ₂ and the parameter D for each case. Note 50 % subsampling is equivalent to a 50% transmittance of the coded apertures in which case D = 0.5(MP), that is, the subsampling rate corresponds to the percentage of the full set of detectors used in each case. The elements of the random coded apertures used are random realizations of Bernoulli random variables, with levels of transmittance equivalent to the corresponding optimal coded aperture transmittance. The PSNR is used to compare the reconstructions obtained since it is suitable for comparing restoration results, as it does not depend strongly on the image intensity scaling. For a scenario with an image I and a reconstruction R of size N × N it is defined as $PSNR=10{\log }_{10}\left(\frac{Ma{x}_I^2}{MSE}\right)$, where Max_I is the maximum possible pixel value of the image I and the mean squared error (MSE), $MSE=\frac{1}{N^2}{\displaystyle {\sum }_{i=0}^{N-1}{\displaystyle {\sum }_{j=0}^{N-1}{[I(i,j)-R(i,j)]}^2}}$. Figure 5(a) shows the PSNR of the random and optimal coded apertures for multiple subsampling rates. For all the simulations in this paper, random coded apertures correspond to the RD subsampling strategy proposed by Kaganovsky et al. in [8]. The random selection of measurements was repeated 10 times, and the mean and standard error for this ensemble were computed and are represented in the graph by the solid line and shaded strips respectively. Note that the PSNR is higher when using optimized coded apertures for all the studied subsampling rates. Particularly, for 27.6%, 31.8%, and 34.8% subsampling rates the PSNR gain is higher than 15 dB. Furthermore, Fig. 5(b) shows that the error Γ of the system matrices obtained using optimized coded apertures is lower than the error obtained when using random coded apertures. This behavior is expected as Γ is the cost function minimized in the proposed optimization algorithm. For all the experiments the object used was a 32 × 32 Shepp-Logan phantom. A second simulation scenario is considered with the same geometrical configuration for the hardware, but changing the parameters M = P = 128 and N = 64. Figures 6(a) and 6(b) show the gradient descent iterations for the error function F associated with the PSF and the complete cost function J for subsampling rates of 47.3% and 20.3% respectively. As it can be seen, the error J and the Frobenius norm of the PSF (F) are minimized over time while maintaining the coded apertures in the solution space. Furthermore, the mutual coherence of the sensing matrix is calculated in each iteration for both subsampling rates and it can be seen in Figs. 6(c) and 6(d) that the algorithm converges to a lower mutual coherence value than the initial coded aperture pattern for both 47.3% and 20.3%. For the latter subsampling rate, reconstructions using a 64 × 64 Walnut phantom [23] Fig. 7(a) were obtained for: (b) Optimized and (c) Random coded apertures. Figures 7(d) and 7(e) depict the normalized absolute error for the respective cases. The error images were normalized for each case by dividing all the values by the maximum error of the 2 images (reconstruction error from random and optimized coded apertures); this process was performed in order to have the same error scale for both the absolute error of the random and optimized coded aperture cases. As it is expected from the 4 dB difference between the PSNR of the reconstructions obtained using optimized and random projections, the error images have higher values and more artifacts in the latter case. The singular value decomposition (SVD) of the system matrix has been reportedly used as a tool to compare different sampling strategies [8, 9]. The SVD of the matrix Q for the optimized and random cases is calculated and it is depicted in Fig. 7(f). For the latter, the mean for 10 different selections is obtained. Considering there is no prior information about the object under inspection, the measurement strategy that has more singular value components lying above certain noise level is considered to outperform the others since it would capture more orthogonal components of the object. It can be seen that the condition number (ratio of greatest singular value to the least non-zero singular value κ) of the optimal codes is better than the random coded apertures case by an order of magnitude. It can be concluded that the use of optimized codes leads to less ill-conditioned sensing matrices compared to using random codes, showing that the proposed optimization framework leads to a better conditioning of the forward operator. The maximum σ_max and minimum σ_min values are specified for the random and optimal cases and are depicted with a ‘*’ marker in Fig. 7(f).

 

Fig. 5 (a) PSNR of the reconstructions obtained using random and optimal coded apertures with different subsampling rates. For the random coded apertures, 10 realizations were computed and the graph depicts the mean and the standard error as a solid line and shaded area respectively. (b) Γ error for optimal and random coded apertures. The subsampling rate corresponds to the percentage of the full set of detectors used in each case.

Download Full Size | PDF

 

Fig. 6 Gradient descent iterations of the optimization function and regularization terms for 47.3% (a) and 20.3% (b) subsampling rates. Gradient descent iterations of the mutual coherence for 47.3% (c) and 20.3% (d) subsampling rates. The subsampling rate corresponds to the percentage of the full set of detectors used in each case. The parameters for the simulation were M = P = 2N = 128.

Download Full Size | PDF

 

Fig. 7 (a) 64 × 64 Walnut phantom [23]. Reconstructions using: (b) optimized and (c) random coded apertures with subsampling rate 20.3%. Absolute error images for the (d) optimized and (e) random cases show more artifacts for random projections.(f) SVD of the CT matrix with optimized codes and random codes. The condition number of both the random κ_r and optimized κ_o cases are specified. The maximum σ_max and minimum σ_min values are specified for the random and optimal cases. Optimized projections result in a less ill-conditioned system.

Download Full Size | PDF

To illustrate the versatility of the algorithm and the motivation for using a Haar wavelet in the majority of the simulations, reconstructions using a Symlet wavelet for optimized and random coded apertures are obtained. The PSNR for the optimized case is 122.28 dB for a subsampling rate of 32.43%. While the PSNR for a subsampling rate of 32.47% using the Haar wavelet is 133.51 dB, that is an 11 dB gain. In both cases, the reconstructions obtained using optimized coded apertures outperform the reconstructions obtained using random coded apertures by approximately 8 dBs. To illustrate the difference between the reconstructions obtained from random and optimized projections in each case, Figs. 8(a) and 8(b), depict the absolute error plots for the Haar wavelet and the Symlet wavelet respectively. In both figures, the left surface plot corresponds to the normalized absolute error of the reconstruction using optimized coded apertures and the right surface plot corresponds to the normalized absolute error for the reconstructions using random coded apertures. Note in both cases, the error is higher when using random coded apertures.

 

Fig. 8 Normalized absolute error plots for the reconstructions of the 64 × 64 Shepp-Logan phantom for subsampling rates of 32.47% using a Haar wavelet (a) and 32.4% using a Symlet wavelet (b) and M = P = 128.

Download Full Size | PDF

Figure 9 depicts reconstructions for the third scenario, in which real tomographic X-ray projections of a slice of a lotus root, made available by Bubba et al. in [32], are used. The optimization of the coded apertures for this system was performed using the system matrix and the sinogram provided by the authors corresponding to P = 120 view angles, M = 429 detectors and a discretized object of 128 × 128 pixels (N = 128). Figure 9(a) depicts the reconstruction of the lotus slice from real X-ray data using Landweber iterations without using coded apertures [32], that is, D = MP = 51480 and Figs. 9(b) and 9(c) depict the results obtained when using coded apertures with subsampling rate 20.6%. Note that the number of projection angles does not correspond to the optimal case (P = 2N) which results in more artifacts for the three cases. However, as shown before, the reconstruction results obtained when using optimized coded apertures have fewer artifacts than the results obtained using random coded apertures. Note that some of the components on the left-hand side of the pictures cannot be distinguished in Fig. 9(c) but they appear in Figs. 9(a) and 9(b). Figures 9(d) and 9(e) depict the normalized absolute error of the reconstructions using random and optimized coded apertures, respectively, compared to the reconstruction using all the available data Fig. 9(a). Note the reconstruction obtained using optimal coded apertures presents a lower error for the majority of the pixels in the image. Furthermore, the maximum error obtained for the optimized coded apertures case is 0.6 while in the random case is 1. As before, the SVD of the system matrix, as well as the condition number, were calculated for both the optimized and random coded aperture cases, the results are shown in Fig. 9(f). Note that the system matrix is less ill-conditioned when using optimized coded apertures. As mentioned in Section 3, to achieve uniform sensing, it is ideal that the number of rays that sense a particular pixel is approximately uniform across the ROI. Figure 10 shows the intensity graphs corresponding to the number of rays that measure certain pixel for different subsampling rates using both optimized and random coded apertures. The mean and standard deviation are calculated for each case in the circle of radius R = N which denotes the ROI for fan-beam X-ray CT systems. Note that for D = 7742 which corresponds to D > N² and D = 3327 which corresponds to D < N² the intensity graphs show that the optimized case results in a more uniform sensing of the object; particularly, note that the standard deviation is lower when using optimized coded apertures. Furthermore, the maximum radiation is obtained when using random coded apertures for both subsampling rates as shown by the hot-spots in the intensity graphs. This concentration is in general not ideal. It should be noted that the regularization term V₂ could be modified to control the radiation at different points. For example, to radiate specific parts of the object in applications with multiple ROIs.

 

Fig. 9 Reconstructions of a slice of a lotus root from real X-ray tomographic projections, using: (a) Landweber iterations without using coded apertures [32] and (b) optimized and (c) random coded apertures with subsampling rate 20.6%. Absolute error images for the (d) optimized and (e) random cases are depicted.(f) SVD of the CT matrix with optimized codes and random codes. The condition number of both the random κ_r and optimized κ_o cases are specified. The maximum σ_max and minimum σ_min values are specified for the random and optimal cases.

Download Full Size | PDF

 

Fig. 10 Intensity graphs of the number of rays that measure certain pixel for D = 7742, and D = 3327 i.e. 47.3% and 20.3% subsampling rates respectively, with the corresponding standard deviation, mean and maximum values. Hot Spots, with radiation higher than the 95% of the maximum intensity, are located in the random case for both subsampling rates and are shown as red dots in the figure.

Download Full Size | PDF

5. Conclusion

The analysis of the mutual coherence and the restricted isometry constant of the compressive X-ray CT system allows the formulation of an optimization algorithm based on the PSF of the system matrix to design the coded aperture patterns in a coded aperture compressive X-ray CT system. The performance of the proposed optimal codes has been compared to the random approach introduced in [8] using object-independent performance measures such as SVD, the mutual coherence, and the uniformity criteria. Additionally, reconstruction results for specific objects with synthetic and real data show that the coded apertures obtained through the proposed design lead to reconstructions with higher PSNR than the ones obtained using random coded apertures. This paper considers a fan-beam geometry and represents the object using Haar and Symlet wavelet bases, but the framework presented here is far more general and can be applied to other geometries and to other sparse representations of the object. Furthermore, the analysis can be extended to other imaging systems that make use of the X-ray transform, such as coded aperture snapshot spectral imaging (CASSI), compressive X-ray tomosynthesis and spectral CT, by adapting the cost function according to the corresponding system matrix. Different transmittance levels can be attained according to the constraints imposed by the system by changing a regularization constant in the proposed cost function and using the random stochastic gradient descent algorithm proposed in this work. Additionally, a regularization term can be added to the cost function to induce a more uniform sampling in multiple regions of interest. Fabrication of the coded apertures for circular tomography is under development. Calibration procedures will be used in order to mitigate the mismatching errors that might occur. Furthermore, the angular collimation produced by the implemented coded apertures can be accounted for in the system matrix H. The optimized codes would take into account this phenomenon.