Sources and Coded Aperture Transmittance Analysis in Compressive Computed Tomography 1

Computed tomography (CT) allows the three-dimensional internal structure reconstruction of an object illuminated with X-ray light. In CT, a set of two-dimensional projections are taken to reconstruct the underlying object structure. The number of projections needed for sensing a CT scene is determined by the Nyquist limit. In some cases, the imposed projections number is excessive. Compressive sensing (CS) has emerged as a new sampling technique requiring fewer projections than those specified by the Nyquist criterion. Instead of measuring the samples directly, they are encoded before being integrated into the detector. This paper describes a CS system for CT based on coded apertures. An optimized value of transmittance and an aperture distribution are selected such that the quality of reconstruction is maximized. Simulations show that results in reconstruction with 50% of measurements are comparable with the traditional CT method based on Nyquist criterion. Similarly, results indicate that the PSNR of reconstructed images can be controlled according to the number of projections taken.


Introduction
Computed tomography (CT) is a technology established for the non-invasive acquisition of images from the internal structure of the objects in three dimensions (3D) [1].In CT, an object can be reconstructed from a set of two-dimensional projections, which are produced by illuminating the object with an X-ray source.CT is important for medical diagnosis and its applications address several areas of industry, material science, biophysics, among others [2]- [4].
The CT image reconstruction algorithms have been restricted by the Nyquist theorem [5]; thus, in some CT applications, the established number of projections is excessive since high X-ray dose could be destructive or carcinogenic for human beings [6].For this reason, it appears the interest to find acquisition methods for reducing the object exposure to the X-ray radiation, without sacrificing the quality of the images.
Furthermore, compressed sensing (CS) has recently emerged as a branch of signal processing.It is based on the fact that many signals in nature can be represented by a few coefficients in some representation basis [7].In CS, the samples are coded in order to reduce the data redundancy in a scene; these coded measurements are enough to reconstruct the signal with a comparable quality of the signal sampled following the Nyquist theorem [8].
Data reconstruction in traditional CT is performed from projections obtained in a detector that measures the X-ray attenuation after they pass through the object.To implement the CS theory in a CT configuration it is possible to include coded apertures that allow coding the measurements to take compressed samples.Coded apertures are two-dimensional arrays with patterns of opaque materials; sections that do not contain opaque material define the transmittance of the aperture, that is, the light fraction that can pierce them [7].In CT, this parameter can define the quantity of radiation that the object receives from each projection.Some works have focused on developing models for joining CS to CT [9]- [11].Nevertheless, the measurement architectures design to improve the quality of the reconstructions has not been detailed.In this work, a multisource CT system is designed from a measurement matrix model by using CT theory and the sampling principle of coded measurements, thus, the geometry of the architecture is described and an analysis of the parameters of coded apertures and number of X-ray sources are performed to ensure efficient tomography reconstruction.In order to specify the results, this study is organized as follows: Section 1 explains in detail the traditional acquisition model, compressed sensing theory, and the proposed model; Section 2 brings coded apertures design; Section 3 describes reconstruction algorithms; Section 4 shows simulations results, and finally, the discussion section is presented.

Acquisition Methods
Traditional CT method (known as transmission CT) is based on estimating the X-ray attenuation after the rays pass through the underlying object [1].In the last few years, some methods to measure changes in the X-ray intensity have been developed .These methods use multiple sources and optical instruments which modulate the beams to improve quality and to reduce the radiation dose.
Recently, systems based on the influence of the number of X-ray sources over CT reconstructions have been proposed.For instance the systems, Pioneering Dynamic Spatial Reconstruction of Mayo Clinic [12], Line Sources CT [13], Inverse-geometry CT, and Ultimately a rebirth of stationary CT [14] have multiple X-ray sources distributed in different ways.The usage of more sources allows reducing the data acquisition time, decreasing the radiation dose, and increasing the spatial resolution of images and reconstruction accuracy.Although developments of acquisition technologies have centered on detection systems, it is expected that the next advance focuses on the use of multiple X-ray sources [15].
On the other hand, compressive sensing technique has been used for transmission CT reconstruction, and thus Compressive CT (CCT) method has been established.CCT is investigated in multisource systems [10].The implementation of CS in CT seeks to reduce the number of measurements in the sampling process without risking the image quality.In the next subsections, the measure acquisition model in CT and CS is detailed and the CCT experimental system is also detailed.

Transmission CT
Figure 1 shows an object Q discretized when it is lightened with X-rays from a source S.This representation is composed of voxels that form Q 1 cross-sections or slices, each one with Q 2 × Q 3 superficial dimensions.For each voxel (with dimensions d x × d y × d z ) an attenuation µ is allocated.Considering a ray with intensity I 0 , which passes through the object with a non-homogeneous attenuation distribution μ(x), the ray intensity I(x) measured by the detector D depends on the distance x passed through and the attenuation μ(x) of each point in its trajectory.This can be modeled as: This phenomenon obeys Lambert Beer's law, and can be rewritten as: Eq. ( 2) is described as the fraction of transmitted light through the object and it can be assigned to a one-dimensional projection, received with regard to an angle of incidence.Reconstruction lies in estimating the distribution of the attenuation coefficients μ(x) from that information [1].

Compressive Sensing
Compressive sensing (CS) is a new theory to acquire and to reconstruct a signal efficiently by searching a sparse solution to an underdetermined system of linear equations [16].As opposed to the traditional signal acquisition process, CS allows sampling rates close to signal intrinsic information rate, which is much lower than Nyquist criteria.
The CS theory is based on two conditions principally: (1) the sparsity of digital signals and (2) the incoherence of the measurement matrix which depends on the sensing trajectory [17].
An image is sparse if the most of its elements are close or are equal to zero.Assuming a linear measurement process that calculates M« N internal products between f and a collection of vectors {φ j } j=1 M , as y i = <f, φ j >, then, where the set of y i projections form the vector y of M elements, Φ ∈ ℝ M × N is the measurement matrix, φ j T is the jth row of Φ, and f is the underlying signal.Taken into account the reconstruction of f from y, it is known that there exist infinite solutions for the equation (3), because the dimensionality of y is much lesser than the dimensionality of f.Thus, there are fewer equations than unknowns.
CS exploits the principle of most of the natural signals can be expressed in a suitable basis with a small number of coefficients.Sparsity is a key requirement for the application of CS theory.However, many signals of nature are not sparse, but they can become sparse adopting a transformation.For instance, a CT image has non-zero values in most of its pixels, as opposed to his wavelet coefficients; in these basis functions, the non-zero coefficients are sparse and they contain the most important information of the original image.
Mathematically, a discrete signal f ∈ ℝ n can be expressed as: where x is the coefficients sequence of f, and the basis Ψ is a matrix with columns ϕ 1 , ..., ϕ G .Clearly, f and x are equivalent representations of the same signal, f is a linear combination of barely F base vectors, F << G [17].
In accordance with equations ( 3) and ( 4): where Θ = ΦΨ ∈ ℝ M × N is the sensing matrix.The undetermined equation system in (5) causes the recuperation of x to be impossible without more information.Fortunately, in CS it is possible to solve (5) if we satisfy the condition that the measurement matrix Φ is incoherent with the sparse transformation Ψ.The incoherence means that object to be reconstructed having a sparse representation in Ψ cannot be sparse in the domain in which it was acquired.Coherence measures the highest correlation between the elements of Φ and Ψ.If Φ and Ψ contain correlated elements, coherence is high, otherwise it is low.
The condition of incoherence can preserve the information, since it requires that energy of the signal to be distributed for the entire detection domain.Each measurement has information of all the image components.Ideally, a measurement matrix ensures that relevant information on any compressible signal is not damaged by the reduced dimensionality from f to y [17], [18].

Compressive CT
A strategy to introduce CS theory in a CT configuration consists of including elements of the system that allow coding the measurements to get incoherence samples [19].These elements are coded apertures, due to their effects on the light.Coded apertures are two-dimensional arrays with patterns of opaque materials; sections that do not contain opaque material define their transmittance, i.e., the fraction of light that can pass through them.Figure 2 depicts the geometry of the sampling system for compressive CT used in this work.A system with multiple sources setting in an array including coded apertures T 1 and T 2 , which modulate the projected beams from the i-th source to a plane with multiple detectors that measure the attenuation generated by an object f.In this case, the coded apertures are located between the X-ray source and the object; therefore, they modulate the energy of conical beams producing a coded projection in the detector plane.The coded aperture elements of dimensions N × N, located between the i-th source S i and the object, are denoted as { } , where 0 blocks the X-ray beam and 1 allows it to cross .The matrix ) is defined.To generalize the measurement for several sources the matrix Φ is defined as the concatenation of the Φ i matrices, i.e., = 0 1 ... P 1 ; and the matrix Φ i is the sensing matrix related with the ith source.Defining the data cube arranging in the representation basis Ψ, the measurements that will be detected by the X-ray sensor can be modeled as: where P is the total number of sources, Φ is the sensing matrix, and Θ is the CS matrix for CT. Figure 3 illustrates a sampling matrix Φ for an object of two slices of Q 1 × Q 3 and a detector array of N 2 where each row represent the trajectory of a single X-ray and the grayscale values represent the different attenuation values.If measurement in a single projection is not sufficient for reconstruction, additional projections or shots per source are required, each with a distinct coded aperture.The number of shots can be expressed in terms of the compression ratio (Cr) and it is defined as: Where K is the number of shots.
The number of sources used in the CCT configuration may vary.Also, they are disposed in an array.Figure 4 shows some examples of source distributions.Figure 4a illustrates a distribution of 2 sources on the x and y-axes, Figure 4b shows a distribution of 3 sources on the x and y-axes.

Coded Aperture Design
Coded apertures are grids with patterns of opaque materials.Transmittance is the quantity of energy that passes through an object, for this case, it is the magnitude that expresses the quantity of X-ray light that is definitely projected over an object [20].The coded apertures have been developed and tested on X-ray systems based on fan beam geometry [21].The coded pattern is modeled in Matlab, a mold of the aperture is printed and is filled with Tungsten powder and it is sealed with epoxy resin.Source: authors' own presentation In a system based on coded apertures, the quality of reconstructions depends on the selection of the coded apertures used in the sampling.Coded apertures are traditionally used in multi-spectral signals sampling systems [22].Coded apertures used in such systems employ type Boolean, binary, grayscale, and Hadamard random codes.
Boolean patterns where the (x, y) element of the ith coded aperture, { }, have allowed getting the best results in spectral image reconstruction [20].
In this work, the Boolean pattern is employed, which not only encodes X-ray signals but also reduces the exposure of the object.The transmittance of a coded aperture is calculated as: Where N 2 represents the size of the coded aperture.Figure 5 illustrates different Boolean random matrixes representing coded apertures like the ones used in this work.Figure 5a shows one coded aperture with a transmittance of 30%, Figure 5b with a transmittance of 50%, and Figure 5c with a transmittance of 70%, this means that 30%, 50% and 70% of the elements, randomly distributed, of the apertures allow light crossing, respectively.Source: authors' own presentation

CT Reconstruction
The image reconstruction problem is to assign the suitable transmittance µ to each voxel in order to discretize the object.To realize such assignment, analytical or iterative methods are used.Analytical methods are the direct solution of a linear equation system; two examples of these methods are back projection and filtered back projection FBP, which is known to be very fast [1].Iterative methods include Algebraic Reconstruction Technique (ART), and iterative algorithms, such as the Simultaneous Iterative Reconstruction Technique (SIRT) [23], [24] used in this work for comparison.The SIRT consists of three phases, executed in an iterative fashion: (1) projection of the estimated object, (2) correction factor computation (updates), and (3) back projection updates of the estimated object.The SIRT algorithm is one of many methods to solve the system of linear equations ( 8) by minimizing: Where x represents the image, b represents the projections, and A represents the scanning process.
SIRT alternates forward and back projections.Its update equation is: Where C and R are diagonal matrices that contain the inverse of the sum of the columns and rows of the system matrix, and the transposed matrix A T back projects the projection images onto the reconstruction area.Given a ray sum, it describes which pixels are hit by that ray.
For these traditional reconstruction algorithms the number of projections must satisfy the Nyquist limit to avoid streaking artifacts.

CS Reconstruction
Given a set of measurements y, which depend on the source configuration, the reconstruction of the CT projections focuses on solving a linear underdeterminate equation system by estimating f as an optimization problem.For CT reconstruction it is necessary to use algorithms that can be adjusted to converge quickly on problems of this type; in this work the results obtained with the Gradient Projection for Sparse Reconstruction (GPSR) algorithm and the TwIST algorithm are presented.The GPSR method is an algorithm used for spectral image estimation with the assumption that the signal of interest is sparse or compressible in some basis.Then the reconstruction consists on recovering x such that the cost function is minimized as where Ψ is a sparse representation basis, Θ is the compressive sensing matrix, x is the sparse coefficients vector, the parameter τ is a regularization constant, and 1 and 2 correspond to the l 1 y l 2 norms, respectively.On the other side, the TwIST algorithm [25] is another framework used frequently to recover signals from compressed measurements.TwIST describes a data cube as the solution to the minimization problem: where the choices for the reguylarization function H TV (x) include, but are not limited to, the l 1 norm.Traditionally, TwIST uses the total variation (TV) regularizer H TV (x) given by Secondly, a real data cube is used from a lungs medical test; the data cube has 60 slices of 1024 × 1024 pixels.Figure 7 shows some cross sections of the real data cube.Thus, simulations to determine the effect of some parameters on the quality of reconstruction are performed.The first parameter to consider was the number of sources, for this experiment coded apertures with a transmittance of 50% were used, and the compression ratio was 0.5, established because that is the mean value.The average results obtained from the reconstruction with the GPRS algorithm are presented in Table 1.In the last column, the results of 200 iterations of the SIRT algorithm are included as a traditional reconstruction reference.As shown in Table 1, experiments were performed with different number and distribution of the sources.The best result in PSNR values is 31.98 [dB], with a distribution of 4 sources in the x-axis and 3 sources in the y-axis, a value that does not change significantly if the sources are increased in any axis due to the geometry and data cube size.These results show the effect of the number of sources on the image quality.It can be observed that 9 sources are sufficient to obtain high PSNR values with the system using a detector array of 64 × 64.Furthermore, it can be concluded that the results of the reconstructions with 50% of the measurements are comparable, and in some cases of higher quality, than those made with the number of measurements required for Nyquist.
For the second experiment random coded apertures were designed with different transmittance values; it is possible to find the relationship between the transmittance values of the coded apertures and the reconstruction quality.The geometry of 3 × 3 sources, 64 × 64 detector array, coded apertures, and the compression ratio of 0.5 were used for simulations.Figure 8 shows the results (in PSNR average values) obtained in ten repetitions for each specific transmittance value with random coded apertures.It shows that the simulations with the best PSNR were those corresponding to a transmittance value of 50% for the GPSR algorithm.Figure 9 shows the results of simulations performed for the TwIST algorithm were the best result was achieved with a transmittance value of 80%.However, it is lower than the GPSR results.Results of reconstructed cross sections of the synthetic data cube are presented in Figure 10. Figure 10a shows a cross section of the original data cube, Figure 10b shows a cross section reconstructed with the SIRT algorithm (traditional method), using a source array of 3 × 3 (PSNR = 30.01[dB]); Figure 10c shows a cross section reconstructed with the GPSR algorithm, using a source array of 3 × 3 and coded apertures with transmittance of 50% (PSNR = 31.22[dB]); and Figure 10d shows a cross section reconstructed with the TwIST algorithm, using a source array of 3 × 3 and coded apertures with transmittance of 80% (PSNR = 21.46 [dB]).It can be noted the quality of the reconstructed image with the GPSR algorithm is higher than the quality obtained by other algorithms.
In addition, in Figure 11 reconstructed cross sections from a real data cube of a human patient thorax are shown.The configuration used for this slice was a detector matrix of 128 × 128 and 9 sources, and the compression ratio was 0.33 as it is shown in Figure 4b. Figure 11a shows a cross section of the original data cube.Figure 11b shows a cross section reconstructed with SIRT algorithm (Nyquist), PSNR=33.37 [dB]. Figure 11c shows a cross section reconstructed with GPSR algorithm, transmittance of 50% and PSNR=35.07 [dB].Table 2 shows the results of simulations for different compression ratios with the GPSR algorithm.It can be seen that the quality of reconstruction with 4 and 9 sources increases while with 1 source remains unchanged.It can be observed in Table 2 that the image quality increases when more measurements are used.Source: authors' own presentation

Discussion
The compressive sensing technique allows compressing a signal at the acquisition step by using different sample patterns and the projections needed to recover data are less than the Nyquist rate.Some previous works have focused on developing models for joining CS to CT, nevertheless, it has not been detailed the measurement architectures design to improve the quality of the reconstructions.
We have proposed a multisource CT system for measuring coded projections physically by using coded apertures.It is demonstrated that physical coding can be used for data compression in CT.
Simulations indicate that the compressive CT architecture provides comparable results to those achieved with traditional CT architectures and the reconstruction algorithms used for CS require a fewer number of measurements than traditional algorithms used in CT to obtain comparable results.
In this work it was found that the number of sources, their distribution and the transmittance of coded apertures are important factors that largely define the quality of reconstructions.Additionally, with the result of simulations, it was found that the transmittance with the best results obtained is 50% and the corresponding average PSNR is 31.01[dB].

Figure 1 .
Figure 1.Diagram of a cone beam flat panel X-ray transmission system.S is an X-ray source, Q is an object, D is a detector array (flat panel detector), I 0 and I(x) are initial and measured intensities, respectively

Figure 2 .
Figure 2. CS system for CT based on coded apertures

Figure 3 .
Figure 3.The matrix Φ is presented for an object of two slices of Q 1 × Q 3 = 8 × 8. Grayscale pixels represent the different attenuation values in the trajectory of x-rays

Figure 4 .
Figure 4. Different source distributions.(a) 2 sources on the x and y-axes.(b) 3 Sources on the x and y-axes

Figure 7 .
Figure 7. Cross sections of human lungs

Figure 8 .
Figure 8. PSNR average for each transmittance value for GPSR

Figure 9 .
Figure 9. PSNR average for each transmittance value for TwIST

Figure 10 .
Figure 10.Comparison between a cross section of the original cube and a cross section reconstructed with 9 sources

Figure 11 .
Figure 11.Comparison between cross sections from a real data cube from a human thorax reconstructed from a configuration of 128 × 128 detectors and 3 × 3 sources

Table 1 .
PSNR average for systems with different number of sources

Table 2 .
PSNR average for systems with different number of sources and shots