1. Introduction

Computed tomography (CT) is a computational imaging technique used to reconstruct high-resolution three-dimensional (3D) volumetric images. It has been extensively used in medical/biological imaging, security inspection, non-destructive testing among others [1–3]. CT technology has continuously evolved since its introduction several decades ago. Most common CT systems, in either cone-beam CT or spiral CT architectures, utilize a circular scanning trajectory of the source-detector gantry around the subject under inspection. More recently, robotic CT imaging systems have been introduced as a means to improve targeting accuracy and to reduce radiation exposure to both patient and physician [4,5]. Unlike conventional CT systems, robotic CT systems can adjust their trajectories to maximize dose reduction and image reconstruction accuracy. A twin robotic X-ray system for projection and measurement is, for instance, shown in Fig. 1, where the source and detector are moved on a coordinated trajectory to ensure comfortable scan conditions and lower radiation dose. Given the set of projections and the geometry of the scanning trajectory, algebraic reconstruction techniques (ART) are used to reconstruct the volumetric images with limited-view projections [6].

 

Fig. 1. The schematic diagram of the twin robotic X-ray system with a static coded aperture. The source-detector pair scans the object with multi degree-of-freedom robotic arms along the pre-programmed trajectory. The plane of the source is designed parallel to that of the detector and thus the coded aperture has one-to-one correspondence to the pixels of the 2D detector.

Download Full Size | PDF

In order to further reduce the radiation dose, this work explores the use of coded illumination in robotic CT imaging systems. The use of binary “on-off” coded apertures placed in front of the source, blocking a significant portion of the X-rays transmitted to the subject, has been studied in conventional CT [7–10]. It has been shown that structured illumination in CT imaging reduces the radiation dose yet maintains the reconstruction image quality. Inspired by the recent work in [11,12], our work is based on coded illumination where only a single coded aperture is used in the robotic CT imaging system. The implementation of the proposed system with static coded projections can be obtained by following similar strategies as well as other feasible implementation of structured illumination [13,14]. The coded aperture is fixed and is static with respect to the source as depicted in Fig. 1. We thus explore the concept of coded illumination in robotic CT and explore novel scanning trajectories in robotic CT that can provide significant gains in reconstruction image quality. As in [11,12], our work focuses on the use of a single coded aperture in the robotic CT imaging gantry as depicted in Fig. 1. Our goal is then to design the scanning trajectory of the source-detector pair in a subject independent manner. The approach is to constrain the degrees of freedom in the gantry’s motion to only allow rotational scans with varying angles with respect to a fixed coordinate system. The plane of the source is designed parallel to that of the detector for maximum use of the X-ray projections, and the robotic source-detector pair scans the object with a non-zero pitch angle around a fixed axis. The irregular scanning trajectory driven by the robotic system significantly improves the sampling efficiency in a tilted orbit which is not possible with standard circular orbits in X-ray tomosynthesis.

The irregular projection structure and trajectories used in the proposed coded aperture robotic tomography system (CARTS) leads to complex structures for the sensing matrices at different view angles and positions. We thus introduce an efficient forward model to attain a diversity of view angles through rotational sensing by fixing the coded aperture to remain static. The two structures are equivalent and could be converted based on relative motion. As depicted in Fig. 2, the position and orientation of the X-ray source, the flat detector and the coded aperture are all fixed and the scans are performed by rotating the object. The object is placed on a plane rotating around the $r-$axis which has a pitch angle $\theta$ with the normal vector of the $x_1 - z_1$ plane. Again, this model is used to simplify the analysis of the robotic scanning system depicted in Fig. 1. Note that the proposed method is equivalent to the general cone-beam X-ray CT when $\theta = 0^{\circ }$ and a single coded aperture is used. The imaging performance in this $\theta =0^{\circ }$ case is poor as documented by Kaganovsky et. al. The proposed geometry has two advantages compared with the conventional one. First, the effect of source and coded aperture when $\theta >0^{\circ }$ at different view angles is replaced by the rotation of the object. Furthermore, the generation of multiple structure matrices at different view angles is replaced by that of the rotation matrix, leading to more efficient calculation. Second, the object uniformly rotates around the axis which is misaligned with the coordinate axes leading to more efficient sensing of the inner structure. In the forward model, the rotation matrix is numerically computed with two basic sub-matrices and the compressive reconstruction is performed at a sub-sampling rate. Due to the computational complexity, the existing coded aperture optimization methods cannot be used for large structured matrices in X-ray CT. This paper develops a more efficient algorithm based on the Pearson product-moment correlation coefficient (PPMCC) analysis, which can efficiently optimize the coded aperture associated with the high-dimensional sensing matrix that models the system in Fig. 1.

 

Fig. 2. Equivalent sensing mechanism of the CARTS system shown in Fig. 1. The 3D object illuminated by a stationary cone-beam X-ray source rotates around the $r-$axis with pitch angle $\theta$. The $p$th coded projection is measured by a stationary 2D detector, $p = 1, 2, \ldots, P$. A single block-unblock coded aperture is placed in front of the X-ray source and part of the X-ray radiation is blocked by the blocking elements.

Download Full Size | PDF

This paper is organized as follows. Section 2 introduces the forward model in matrix notation. The methods to generate the rotation matrix and optimize the coded aperture are presented in Section 3 and Section 4, respectively. The results using 3D simulated datasets and simulated projections of real datasets at different sampling rates are provided in Section 5. Conclusions and future work are presented in Section 6.

2. Coded aperture robotic tomography system

2.1 Forward model

As depicted in Fig. 2, a binary coded aperture is placed in front of a stationary X-ray source with a fixed orientation. The elements on the coded aperture have one-to-one correspondence to the pixels on a flat 2D detector. The 3D object rotates around the $r-$axis with pitch angle $\theta$ and part of the 3D object is illuminated by the coded X-ray patterns. The $p$th coded X-ray projection is measured by the 2D detector where $p = 1,2,\ldots,P$, and $P$ is the number of projections. Measurements corresponding to the blocking elements on the coded aperture are discarded. Then the linear attenuation map of the 3D object is reconstructed using the set of tomographic measurements and the projection geometries.

According to the Beer-Lambert law, the transmission model of a monochromatic X-ray path $l$ to a particular pixel on the detector is given by [15]

Due to the pixelated nature of the 2D detector, the 3D object is gridded into $N = N_x \times N_y \times N_z$ voxels where $N_x$, $N_y$ and $N_z$ are the dimensions along the $x_1$, $y_1$ and $z_1$ axes in Fig. 2(b), respectively. The 2D detector is gridded into $M = M_x \times M_y$ pixels where $M_x$ and $M_y$ are the dimensions along the $x_1$ and $y_1$ axes, respectively. Then, the discrete logarithmic measurements in matrix notation can be formulated as

The binary on-off coded aperture $T(x_0,y_0)$, $x_0 = 1, 2, \ldots, M_x$, $y_0 = 1, 2, \ldots, M_y$, is placed in front of the cone-beam X-ray source, and the coded aperture pitch is fixed to have one-to-one correspondence to the pixels of the 2D detector. The equivalent transmission function of the coded aperture can then be formulated as

Furthermore, sampling in X-ray CT can be considered as cutting the object using X-rays and thus more irregular voxels can lead to higher sampling efficiency in the absolute coordinate (3D free space). In our previous work, the number of intersections of X-rays in the object is theoretically demonstrated to determine the amount of information acquisition [16,17]. Limited by the projective geometry, the intersections of projections in X-ray CT are finite. However, misalignment between the voxels and the coordinates during the rotation can lead to more intersections and thus more details can be resolved. As depicted in Fig. 2, the 3D object rotates around the fixed $r-$axis with the pitch-roll angle representation $\theta$ and $\omega$, where $\theta$ is the angle of $r-$axis direction and $\omega$ is the magnitude of the rotation, respectively. Therefore, the coded measurements in the proposed CARTS system with the $r-$axis and the angle representation $\omega$ can be written as

An example of the sensing matrix with pitch angle $\theta = 0$ and $\theta = 45^{\circ }$ for $N_x = N_y = N_z = 8$ and $M_x = M_y = 32$, is depicted in Fig. 3(a) and (b), respectively. The number of unblocking elements on the coded aperture $D$ is 256 and thus the transmittance is 25%. It is seen that the sensing matrix with pitch angle $\theta = 0$ is highly structured and the regions marked in blue are not sensed. The distribution of elements in the sensing matrix with pitch angle $\theta = 45^{\circ }$ is more uniform and the corresponding regions marked in red are sensed. The singular value decomposition (SVD) analysis which has been reportedly used as a tool to compare the performance of coded apertures is depicted in Fig. 3(c). The measurement strategy with larger non-zero singular value components is considered to capture more orthogonal components of the object and thus, it leads to less ill-conditioned reconstruction. The condition numbers of sensing matrix with $\theta = 0$ and $\theta = 45^{\circ }$ are $1.03\times 10^{4}$ and $4.04\times 10^{3}$, respectively. It is seen that the misalignment implementation captures more orthogonal components.

 

Fig. 3. (a) Example of a sensing matrix for X-ray CT using rotational coding with pitch angle $\theta = 0$, $N_x = N_y = N_z = 8$ and $M_x = M_y = 32$. The coded aperture has 25% transmittance. (b) Example of a sensing matrix for X-ray CT using rotational coding with pitch angle $\theta = 45^{\circ }$, $N_x = N_y = N_z = 8$ and $M_x = M_y = 32$. The coded aperture has 25% transmittance. (c) The plot of the singular values as a function of component numbers. The highest and lowest singular values are highlighted in each case.

Download Full Size | PDF

It is difficult to reconstruct the linear attenuation coefficient map of the volumetric object from a single measurement. Multiple projections are measured at a series of $P$ rotations, and thus, the vectorized representation of the coded logarithmic measurements is given by

For simplicity, Eq. (7) can be re-written in matrix notation as

2.2 Rotation matrix

As shown in Eq. (8), the forward model of the proposed CARTS is determined by the coded aperture matrix $\widetilde {\textbf {C}}$, the structure matrix $\widetilde {\textbf {W}}$, and the rotation matrix $\widetilde {\textbf {R}}$. The closed-form solution of the coded aperture matrix $\widetilde {\textbf {C}}$ is provided in Eq. (9) and the structure matrix related to the geometric projections can be numerically calculated by using the open-access ASTRA Tomography Toolbox (“All Scale Tomographic Reconstruction Antwerp”) [18].

In this section, an approach to obtain the rotation matrix $\textbf {R}$ with the angle representation $(\theta,\omega )$ is proposed. Similar to the generation of the structure matrix, it is also difficult to obtain the closed-form solution of the volume portion of a rotated voxel associated with a particular element of the free space. Physically, the $(m,n)$th element in the rotation matrix $\textbf {R}$ represents the intersection volume of the $n$th voxel in the 3D object $\textbf {f}$ and the $m$th element in the space $\Omega$. Therefore, a numerical approach is proposed to estimate the intersection volume. As shown in Fig. 4, every voxel in the 3D object $f(x_1,y_1,z_1)$ is divided into $N_0\times N_0\times N_0$ smaller cubes and the coordinates of the centroids in the right-handed Cartesian coordinate system can be calculated as $x_{c}$, $y_{c}$ and $z_{c}$. Then, the number of smaller cubes for the 3D object is $NN_0^{3}$. Given a rotation axis $\boldsymbol {\mu } = (\mu _x, \mu _y, \mu _z)^{T}$ where $\mu _x^{2} + \mu _y^{2} + \mu _z^{3} = 1$ and a rotation angle $\delta$, the rotation transformation $\textbf {Q}_{\delta }(\boldsymbol {\mu },\delta )$ in the right-handed Cartesian coordinate system can be given by

 

Fig. 4. (a) Free space and the 3D object before rotation; (b) Free space and the 3D object after the rotation around the $r-$axis with angle representation $(\theta,\omega )$. Note that the coordinates are in the right-handed Cartesian coordinate system.

Download Full Size | PDF

Thus, the coordinates of the centroid $x_{c}$, $y_{c}$ and $z_{c}$ with the rotation axis $\boldsymbol {\mu }$ and the rotation angle $\delta$ can be calculated as

Since the projections are measured at a series of $P$ rotations, it is computationally inefficient to numerically estimate $P$ rotation matrices, $\textbf {R}(\textbf {r},\omega _i)$, $i = 1, 2, \ldots, P$. Here, the rotation represented by $\textbf {R}(\textbf {r},\omega )$ is decomposed as $\textbf {R}_{\theta }(\boldsymbol {\mu _1},\theta )\textbf {R}_{\omega }(\boldsymbol {\mu _2},\omega )$ where $\boldsymbol {\mu _1}$ and $\boldsymbol {\mu _2}$ are the normalized vectors along $z-$axis direction and $x-$ or $y-$axis direction, respectively. That is, $\boldsymbol {\mu _1} = (0,0,1)$, and $\boldsymbol {\mu _2} = (0,1,0)$ or $\boldsymbol {\mu _2} = (1,0,0)$. Besides, the $(i+1)$th rotation around the $r-$axis can be formulated by the $i$th rotation around the $r-$axis when the distribution of the angle, $\omega$, is uniform. That is, $\textbf {R}_{\omega }(\boldsymbol {\mu _2},\omega _{i+1}) = \textbf {R}_{\omega }(\boldsymbol {\mu _2},\omega _0)\textbf {R}_{\omega }(\boldsymbol {\mu _2},\omega _i)$, where $\textbf {R}_{\omega }(\boldsymbol {\mu _2},\omega _1)=\textbf {R}_{\omega }(\boldsymbol {\mu _2},\omega _0)$ is the initial rotation matrix with $\omega _0 = 2\pi / P$. Therefore, the $P$ rotation matrix can be calculated by $\textbf {R}_{\theta }(\boldsymbol {\mu _1},\theta )$ and $\textbf {R}_{\omega }(\boldsymbol {\mu _2},\omega _0)$. The numerical estimation aims at obtaining the rotation matrix is detailed in Algorithm 1. It is noted that the rotation matrix in this manuscript is pre-computed and stored with the CPU implementation. However, step 12 to step 18 in Algorithm 1 is a loop with large amount of independent computation requiring few memory. Then a GPU workstation with numerous cores can be used to significantly accelerate the algorithm and make an arbitrary acquisition scan possible.

2.3 Coded aperture optimization

Since the sensing matrix $\widetilde {\textbf {C}}\widetilde {\textbf {W}}\widetilde {\textbf {R}}$ is highly structured, the random coded aperture is not optimal in the reconstruction of CS. Recently, coded aperture optimization methods have been developed based on the RIP condition, the mutual coherence and information theory [10,16,17]. However, these methods cannot be directly used in the optimization of the proposed CARTS system due to two reasons. First, these methods are developed for the conventional CARTS systems where the coded apertures and source orientations at different view angles are different. However, the coded aperture and source orientation are stationary at a series of $P$ rotations in the proposed model. Second, the requirement of memory to calculate the Gram matrix in [10] or correlation matrix in [16] and [17] far exceeds the limit of regular computing platforms.

It has been demonstrated that the variables in the structure matrix represent the observations of the attenuation characteristics of X-ray projections and the coded apertures characterizes reducing the measurements of the structure matrix. The purpose of coded aperture design is to obtain the maximal non-overlappling information carried by unblocking X-ray projections. Thus the coded aperture optimization problem is analogous to the feature selection problem or the dimensionality reduction problem [17]. In [16], the PPMCC matrix is used as a metric to optimize the coded aperture and the cost function is given by

Since the system matrix in the proposed architecture is highly structured, the subset of columns cannot be chosen randomly. In practice, neighboring X-ray projections, with a much higher probability, illuminate a particular voxel at the same time and thus carry more overlapping information. Therefore, for a particular pixel on the coded aperture, the neighboring pixels in the diameter of 9 pixels are chosen as a subset. Then the PPMCC of the rows in the sensing matrix $\textbf {H}$ that correspond to the X-ray projections are calculated. The PPMCC of all pixels on the coded aperture with 9 neighboring pixels are calculated, and the average results are depicted in Fig. 5. Figure 5(a)-(d) illustrate the average coefficients at $\omega = 0^{\circ }$, $\omega = 90^{\circ }$, $\omega = 180^{\circ }$ and $\omega = 270^{\circ }$, respectively. It is seen that the PPMCC of the pixels with shorter distance to the neighboring pixel are much larger than that with longer distance. The PPMCC of the pixel pairs decreases much fast as the radius increases and the coefficients are almost zero at the diameter of 9 pixel. Therefore, the distance of the unblocking pixel on the coded aperture is used as the metric and the purpose of the optimization is to make the distance of any two neighboring unblocking elements as far away as possible. It is equivalent to the design where the unblocking elements are spread as homogeneously as possible and such distribution of binary pixels are blue noise coded apertures [19]. The direct binary search (DBS) method in our previous work can be directly used to efficiently obtain the blue noise coded aperture [20].

 

Fig. 5. Average PPMCC distribution of a particular element in the coded aperture and the neighboring elements in the radius of 9 elements at (a) $\omega = 0^{\circ }$, (b)$\omega = 90^{\circ }$, (c) $\omega = 180^{\circ }$ and (d) $\omega = 270^{\circ }$ for $M_x = 256$ and $M_y = 256$.

Download Full Size | PDF

2.4 Reconstruction

The discrete datacube for reconstruction often has a coarse resolution along the $z$ axial direction in that $N_z \ll min(N_x, N_y)$. With such resolution, the correlations between transverse slices are small and Discrete Cosine Transform (DCT) basis is used in $\ell _1$ norm to estimate the sparsity. Then, Kronecker CS reconstruction is often used in these cases [8,21]. In the cases where the resolution along the $z$ axial direction is equal to that of slices, Total Variation (TV) is often used as an implicit sparse basis (or dictionary) to minimize the sum of non-zero gradients. Then, the cost function for reconstruction is given by

3. Simulation results

3.1 Results of simulated datasets

In this section, we present simulations to illustrate the performance of the proposed CARTS system with an optimal coded aperture. The simulations are performed with a cone-beam X-ray source and a flat imager with $256 \times 256$ pixels. The coded aperture with $256 \times 256$ pixels, having one-to-one correspondence to the element on the detector, is placed in front of the X-ray source and the source orientation is disk-shaped in the center of the space $\Omega$. The side length of a pixel on the detector is 0.5 cm and thus the geometric length of the imaging system is 128 cm. The source-to-object and object-to-detector distances are 80 cm and 55 cm, respectively. The object rotates around the $r-$axis with pitch angle $\theta$ and the distribution of angle $\omega$ is uniform at the range of $0 \sim 2\pi$. The imaging system captures the coded X-ray projections at a series of 8 rotations and the 3D object is discretized into $64 \times 64 \times 64$ voxels of dimensions $0.5 \text {cm} \times 0.5 \text {cm} \times 0.5 \text {cm}$. That is, $M = 256^{2}$, $N = 64^{3}$ and $P = 8$. The discrete-to-discrete structure matrix is numerically obtained using the ASTRA Tomography Toolbox [18] and the discrete-to-discrete rotation matrix is calculated using the proposed algorithm in Section II-B. The blue noise coded aperture is generated utilizing the DBS algorithm proposed in [22]. The conventional CARTS where $\theta = 0^{\circ }$ is used to compared with the proposed geometry. The source-plane-to-object-plane and object-to-detector distances are 80 cm and 55 cm, respectively. The simulations are performed in a desktop architecture with an Intel Core i3 3.7GHz processor, 16G RAM, using Matlab 2015b. The object used in the simulations is generated by the 3D Shepp-Logan phantom [8]. To quantitatively evaluate the reconstructed images, the Peak Signal-to-Noise Ratio (PSNR) is used as

The pitch angle $\theta$ is determined numerically. An optimal first-order method for strongly convex TV regularization is used to reconstruct the 3D cubes with binary random coding and blue noise coding, respectively [23]. Since the X-ray projections through the $x_1-z_1$ plane and the $y_1-z_1$ plane are different, the simulations of the two cases are both performed. The results of the PSNR with 26% transmittance at the range of $0 \sim 90^{\circ }$ are shown in Table 1. Therefore, the average PSNR using optimized coding at $\theta = 0^{\circ }$, $15^{\circ }$, $30^{\circ }$, $45^{\circ }$, $60^{\circ }$, $75^{\circ }$ and $90^{\circ }$ are 29.2dB, 32.8dB, 38.1dB, 44.7dB, 38.1dB, 32.8dB and 29.2dB, respectively. The average PSNR using random coding at $\theta = 0^{\circ }$, $15^{\circ }$, $30^{\circ }$, $45^{\circ }$, $60^{\circ }$, $75^{\circ }$ and $90^{\circ }$ are 27.6dB, 31.3dB, 38.7dB, 42.7dB, 38.7dB, 31.3dB and 27.6dB, respectively. The strategy referred to as $\theta = 45^{\circ }$ yields the best performance. Thus, the pitch angle between the $r-$axis and the normal vector is $45^{\circ }$ in the following simulations.

Table 1. PSNR as a function of angle $\theta$ using optimized coding and random coding in the cases of $x_1-z_1$ plane and $y_1-z_1$ plane

View Table

The reconstructed tomographic images using random coding and optimized coding with 26% transmittance are depicted in Fig. 6. The execution time for the optimization of coded aperture with 26% transmittance is $\sim$198 seconds. That is, the sub-sampling rate in the CS reconstruction is 52%. As described in Section II-C, the proposed optimization method is not based on the structure matrix $\textbf {W}$ which is the computational burden of the other methods, the runtime of the optimization method is reduced by four orders of magnitude compared with the gradient approach utilized in [10] and one order of magnitude compared with the method based on information acquisition utilized in [16]. Note that the number of the elements in the structure matrix of the proposed method is more than $10^{11}$ while that of the methods in [10] and [16] are $\sim 10^{9}$. Furthermore, the memory requirements in the cases for the methods in [10] and [16] are more than 1TB. The convergence curves of $64\times 64\times 64$ reconstructed images using random coding and optimized coding are depicted in Fig. 7 (a) and (b), respectively.

 

Fig. 6. (a1) - (a5) are the reconstructed images with 26% transmittance of the 12th layer, 22th layer, 32th layer, 42th layer and 52th layer using the conventional system. (b1) - (b5) are the reconstructed images with 26% transmittance of the corresponding layers using random coding, respectively at pitch angle $\theta = 45^{\circ }$. (c1) - (c5) are the reconstructed images with 26% transmittance of the corresponding layers using optimized coding at pitch angle $\theta = 45^{\circ }$. (d1) - (d5) are the absolute error maps of (b1) - (b5) and (e1) - (e5) are the absolute error maps of (c1) - (c5), respectively.

Download Full Size | PDF

 

Fig. 7. The convergence curves of $64\times 64\times 64$ reconstructed images using (a) random coding and (b) optimized coding.

Download Full Size | PDF

The reconstructed images of the 12th layer, 22th layer, 32th layer, 42th layer and 52th layer, based on the conventional CARTS are shown in Fig. 6 (a1) - (a5), respectively. The corresponding PSNR are 17.3dB, 22.3dB, 19.4dB, 14.2dB and 8.6dB, respectively. It is seen that the structures of the upper layers ($12$th and $22$th layers) are basically reconstructed and the details are drowned out by artifacts. The structures of the lower layers ($42$th and $52$th layers) are not even reconstructed. The reconstructed images of the 12th layer, 22th layer, 32th layer, 42th layer and 52th layer, using stationary random coded aperture are shown in Fig. 6 (b1) - (b5), respectively. The corresponding PSNR are 35.2dB, 48.7dB, 52.0dB, 50.6dB and 49.8dB, respectively. Compared with the conventional system, the intensity distributions of the 3D Shepp-Logan phantom varying along the slices are reconstructed without many artifacts at both upper layers and lower layers. The PSNR gain is $\sim$30.9dB. Besides, the reconstructed images of the 12th layer, 22th layer, 32th layer, 42th layer and 52th layer, using stationary optimized coded aperture are shown in Fig. 6 (c1) - (c5), respectively. The corresponding PSNR are 37.3dB, 49.4dB, 53.4dB, 52.0dB and 50.6dB, respectively. Compared with random coding, the PSNR gains of using optimized coding are 2.1dB, 0.7dB, 1.4dB, 1.4dB and 0.8dB, respectively. The absolute error maps of using random coding and optimized coding are depicted in Fig. 6 (d1) - (d5) and (e1) - (e5), respectively. Note that the maximum value of the colorbar for a tomographic image is the maximum value of the reconstructed cubes which is 0.15 and the minimum value is 0. The red squares in the absolute error maps are used to highlight the difference in the reconstructed images using random coding and optimized coding. It is seen that the optimized coded aperture leads to fewer artifacts in the absolute error maps and the intensity of artifacts using optimized coding is less than that using random coding.

Furthermore, the plot of average PSNR of a $64\times 64\times 64$ 3D phantom as a function of transmittance for optimized coding and random coding is depicted in Fig. 8. It is seen that the proposed method performs much better than the conventional coded aperture X-ray tomosynthesis system. The $64 \times 64 \times 64$ 3D Shepp-Logan phantom cube can be reconstructed at high quality ($\gt$30dB) at only 16% transmittance. Besides, the reconstruction quality is improved by using the optimized coded aperture compared with that of random coded aperture. Note that the results of random coding are the average of 10 realizations in the simulations. The PSNR gain is up to 3dB at the range of 14% - 32% transmittance.

 

Fig. 8. The plot of average PSNR as a function of transmittance for optimized and random coding for the 3D Shepp-Logan phantom cube with $64 \times 64 \times 64$ voxels.

Download Full Size | PDF

3.2 Simulated results of simulated projections from real datasets

In this section, we present simulations using simulated projections of real datasets obtained in [21]. Experimental data was obtained at Chesapeake Testing Inc., with a Nikon metrology 225/450kV Vault CT scanning system with a 450kV micro-focus X-ray source capable of producing a spot-size down to 80$\mu$m. Multiple X-ray projections over 360 degrees around an object, in this particular case a vivofit watch, were acquired. These projection images were then reconstructed into a full 3D volumetric data set. The 3D data cube was re-sampled to the a cube of size $128 \times 128 \times 128$. The simulations are performed with a cone-beam X-ray source and a flat imager with $512 \times 512$ pixels. The coded aperture with $512 \times 512$ pixels is placed in front of the X-ray source. Note that the elements on the coded aperture have one-to-one correspondence to that on the detector and the source orientation is disk-shaped in the center of the space $\Omega$. The side length of a pixel on the detector is 0.5 cm and thus the geometric length of the imaging system is 256 cm. The source-to-object and object-to-detector distances are 160 cm and 110 cm, respectively. The object rotates around the $r-$axis with pitch angle $\theta = 45^{\circ }$ and the distribution of angle $\omega$ is uniform at the range of $0 \sim 2\pi$. The imaging system captures the coded X-ray projections at a series of 16 rotations and the 3D object is discretized into $128 \times 128\times 128$ voxels of dimensions $0.5 \text {cm} \times 0.5 \text {cm} \times 0.5 \text {cm}$. That is, $M = 512^{2}$, $N = 128^{3}$, $P = 16$ and the dimension of the sensing matrix is $4194304 \times 2097152$ (more than $8.8 \times 10^{12}$ elements in the matrix). It is noted that the content of the object out of the region of interest is truncated to avoid truncation artifacts. The simulations are performed in a desktop architecture with an Intel Core i3 3.7GHz processor, 64G RAM, using Matlab 2019b.

The reconstructed images using random coding and optimized coding with 12% transmittance are depicted in Fig. 9. The execution time for the coded aperture optimization with 12% transmittance is $\sim 5.5\times 10^{3}$ seconds. That is, the sub-sampling rate in the CS reconstruction is 24%. The convergence curves of $128\times 128\times 128$ reconstructed images using random coding and optimized coding are depicted in Fig. 10 (a) and (b), respectively. On the contrast, only two basic sub-matrices and a large structure matrix are required to be numerically calculated in the proposed method. It is seen that the intensity distributions of the 3D data cube varying along the slices are reconstructed without many artifacts using both random coding and optimized coding. The plot of average PSNR as a function of slices for optimized coding and random coding with $128\times 128\times 128$ voxels is depicted in Fig. 11. The PSNR gain is up to 4.8dB and the average PSNR gain for optimized coding is 1.5dB. It is seen that the proposed method outperforms the general CARTS system where $\theta = 0^{\circ }$ by using both optimized coding and random coding.

 

Fig. 9. The reconstructed tomographic images of the “watch” with $128 \times 128 \times 128$ voxels. The number “4”, “25”, “48”, “76”, “93” and “110” represent the $4$th, $25$th, $48$th, $76$th, $93$th and $110$th slices of the data cube, respectively.

Download Full Size | PDF

 

Fig. 10. The convergence curves of $128\times 128\times 128$ reconstructed images using (a) random coding and (b) optimized coding.

Download Full Size | PDF

 

Fig. 11. The plot of average PSNR as a function of slices for optimized coding and random coding with $128\times 128\times 128$ voxels. Conventional CARTS is shown at $\theta = 0^{\circ }$.

Download Full Size | PDF

4. Conclusion

Coded aperture robotic tomography system is introduced to obtain both flexible scanning and lower radiation dose. A static binary “on-off” coded aperture is placed in front of the X-ray source and the robotic source-detector pair scans the object with a pitch angle around a fixed axis. The degrees of freedom in the gantry’s motion are constrained to only rotational scans with respect to a fixed coordinate system. Therefore, the complex structured matrices of the robotic CT systems are simplified using an efficient rotational sensing model based on two basic sub-matrices. The stationary coded aperture is designed based on the PPMCC analysis and the DBS algorithm is used to obtain the optimized coding with one order of magnitude faster computation than state-of-the-art optimization algorithms. Computational experiments using simulated datasets show that the proposed method can reconstruct high-quality tomographic images (>30dB) at only 16% transmittance while that of conventional X-ray tomosynthesis is less than 15dB. Furthermore, the optimized coded aperture leads up to $\sim$3dB PSNR gains in the reconstruction quality compared with that attained by a random coded aperture.