Open Access
Add to BibSonomy
Abstract
Compressive spectral X-ray imaging (CSXI) introduces a pixelated spectral modulator called K-edge coded aperture (KCA) in front of the X-ray source, which enables both, lower dosage to the subject, as well as the capability of spectral tomography while using low-cost integrating X-ray detectors. CSXI systems generally use hundreds of different spectral modulators, each with a distinct pattern to uniquely modulate the illumination at every view angle. In contrast, this paper introduces the use of a single and static coded aperture placed in a tomosynthesis gantry. The compressive system thus interrogates the subject with a fixed coded illumination pattern on all view angles. The advantages of the system are many including reduced cost and the feasibility of implementation. Given the reduced set of coded measurement and the limited spectral separation ability in the resulting architecture, the nonlinear inverse reconstruction problem results in a highly ill-posed problem. An efficient alternating minimization method with three-dimensional total variation regularization is developed for image reconstruction. Furthermore, rather than simply using a random pattern, the coded aperture is optimized under a uniform sensing criterion that shapes the spatial and spectral pattern of the coded aperture so as to minimize the overall radiation exposure placed on any volumetric area of the patient. This is of particular importance in medical imaging where patients at risk are recommended to have periodical X-ray tomosynthesis screenings. The coded aperture optimization is then posed as a binary programming problem solved by a gradient-based algorithm with equilibrium constraints. Numerical experiments show that spatial and spectral coding used in the proposed system to interrogate the subject not only reduces the radiation dose but it also improves the quality of image reconstruction. Gains close to 5dB in peak signal to noise ratio are observed in simulations. Furthermore, it is shown that the optimization of the KCA can effectively improve the uniformity of X-ray radiation compared to random KCA modulation, thus reducing the radiation dose throughout all volumetric sub-areas of the subject — an objective that is not possible with the use of random KCAs.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Tomosynthesis reconstructs multi-layer images of a three-dimensional (3D) object based on projection data acquired from a few incident angles, which allows 3D inspection with limited sensing geometry [1]. Since 2011, tomosynthesis has been applied in the 3D breast imaging [2]. Compared with traditional two-dimensional (2D) radiography, tomosynthesis avoids the superpositions of tissues at different depths, and can thus be applied to detect potential lung cancer nodes and breast cancer lesions [3–7]. A typical tomosynthesis system consists of a poly-energetic X-ray source and an energy integrating detector, which sums up the attenuation information of the object over the entire spectrum. The negative logarithm of the measured data is then used for the filtered backprojection (FBP) reconstruction [8]. Such a system however cannot resolve the energy-spectral attenuation information, which can be used to identify the material components in the object and to reduce the beam hardening artifacts [9–12]. Prototypes of spectral X-ray tomosynthesis systems typically use photon counting detectors, which can simultaneously obtain the measurement data with different energy levels [13,14]. However, the photon counting detector is rarely used in routine clinical examination due to its high cost and the effect of pulse pile-up at high X-ray-flux rates [15–17].
Recently, a cost-effective approach was proposed to acquire measurements containing the spectral attenuation information, which employs a set of balanced K-edge filters to reshape the X-ray spectra in sequential scans [20,21]. For instance, the spectral attenuation coefficient of silver (Ag) is shown in Fig. 1(a). It includes a sudden jump when the energy of incident X-ray photon is equal to the K-edge energy of the filtering material. As shown in Fig. 1(b), using two filters with adjacent atomic numbers and balanced thicknesses, a quasi-monochromatic spectrum can be obtained by subtracting the filtered intensity from each other. However, the radiation dose and the inspection time will be multiplied due to the sequential scans applied or a set of K-edge filters [22,23]. To overcome these limitations, a compressive spectral X-ray imaging model based on encoded illumination was proposed to collect the spectral attenuation data from a single scan [24,25]. The X-ray beams in compressive spectral X-ray imaging (CSXI) system are spatially and spectrally modulated by multiple coding masks consisting of randomly arranged pixelated K-edge filters, which are referred to as K-edge coded apertures (KCA). Subsequently, a multi-stage algorithm based on sinogram inpainting was proposed to reconstruct the energy-binned images [24], it first inpaints the sinograms of filtered intensity based on compressive sensing and then estimates the energy-binned sinograms that can be used as the inputs to reconstruct the energy-binned CT images [24]. Furthermore, these images can be decomposed as the images of given basis materials [25]. More recently, a nonlinear inversion approach based on the material decomposition model was proposed to directly reconstruct the energy-binned tomographic images from the CSXI measurements [26]. Despite the advanced reconstruction performance, the coded aperture used at each view angle in CSXI is different, which significantly increases the system complexity, cost, and calibration complexity. Recently, the method coined StaticCodeCT was proposed as a way to overcome the implementation complexity in conventional mono-energetic compressive tomography by using a single static coding mask for all views. The coding mask used by StaticCodeCT however is binary, which cannot modulate the spectrum of the X-ray source [27]. In addition, StaticCodeCT uses randomly patterned coded apertures for the blocking elements that reduce the radiation dose. Random coding generally leads to large variations on the radiation dose across the volume of the subject being inspected. This is undesirable as certain tissue volumes and cells of the patient can be over-exposed, hence suffer more from X-ray radiation damage, which increases the risk of cancer [28].
Fig. 1. Sketch of material’s K-edge attenuation. (a) Spectral attenuation coefficient of silver (Ag) that includes a sudden jump at 25.49 KeV [18]); (b) spectrum of a 50kV X-ray tube (black) [19], and the spectra filtered by a balanced K-edge filter pair of Palladium (Pd) (orange) and Ag (blue).
This paper introduces several contributions. First, it generalizes the static-code-CT architecture, from the conventional single-energy setting to the emerging modality of spectral tomography. By acquiring coded projections at different energy levels, the proposed technique can differentiate among various elements in the body based on their material density or atomic numbers. This capability is important in a number of medical and security applications. The coded aperture used is no longer a simple binary structure that either stops X-rays or allows their path unaltered. The proposed compressive spectral X-ray tomosynthesis (CSXT) approach uses a coded aperture made of K-edge materials that can spatially and spectrally modulate the poly-energetic cone-beam X-ray illumination of the system. The second contribution of the paper is to adopt the static-code-CT architecture, recently introduced in [27] for monochromatic CT, into spectral CT scanning where a static KCA is placed in front of the cone-beam X-ray source. The KCA synchronously rotates with the X-ray tube and detector during the scanning process. Unlike the case of static-code-CT which operates on a monochromatic energy model, the inverse image reconstruction problem in the proposed system becomes significantly more challenging as the problem is now nonlinear and the number of measurements obtained over a scan is proportionally much lower in the higher multi-dimensional setting. Consequently, the inverse imaging problem is severely ill-conditioned for which new reconstruction algorithms are required as is done in Section 3. The third contribution focuses on reducing the radiation dose in a scan. The goal is not only to reduce the overall radiation given to the volume of the subject interrogated [28], but also the minimization of high X-ray dose directed and concentrated at specific areas and tissues in the body during each scan. Dose reduction in concentrated volumetric areas of the subject reduces the obvious risks for women undergoing regular screening mammography [29]. Hence, rather than simply using a random pattern for the coded aperture that generally yields a non-uniform distribution of radiation on different portions of the subject, the coded aperture is optimized under a uniform sensing criterion that shapes the spatial and spectral pattern of the coded aperture so as to minimize the overall radiation exposure placed on any volumetric area of the patient. In order to reduce the overall radiation dose in the scan, opaque elements are also included in the KCA to block a portion of X-rays (see Fig. 2(b)).
Fig. 2. (a) Sketch of compressive spectral X-ray tomosynthesis based on the static KCA shown next to it, where different colors indicate different K-edge materials in the coding mask; (b) down-sampling in spatial domain by blocking a part of X-ray beams, where black squares indicate the blocked pixels.
Figure 2(a) illustrates the proposed system, where a static KCA is equipped in front of the cone-beam X-ray source and synchronously rotates with the X-ray tube during the scanning process. The KCA pixels are made of different K-edge materials, which filter the X-ray beams with different transmittance functions to sense the spectral attenuation information of the object. Similar to the StaticCodeCT proposed in [27], the KCA pattern used in CSXT does not change with the view angle, which effectively reduces the difficulty to fabricate and calibrate such systems. The forward projection imaging process of the CSXT system follows a nonlinear compressive sensing model. According to the theory of basis material decomposition [9,30], it is known that the spectral attenuation of the object is determined by its elements. The forward model is thus reformulated as the projection of the basis materials in the object. This reformulation benefits the spectral tomosynthesis reconstruction since the dimensionality of the variables can be significantly reduced. A 3D total variation regularization is introduced in the reconstruction process to preserve the edges and remove noise in the tomographic images [31–33]. To solve this regularized reconstruction problem, we generalize the alternating projection scheme to the nonlinear problem using the Gauss-Newton method [34–36].
To attain a uniform distribution of X-rays on the subject, we optimize the KCA using a uniform sensing criterion, which demands that each object voxel receives equal X-ray illumination in intensity coming from each type of filtered X-ray beam. We model the KCA optimization problem as an inequality constrained binary matrix optimization problem, which is then efficiently solved by a two-stage algorithm. Simulations under 50% dose setting show that the optimized KCA can improve the uniformity of X-ray radiation by 54.14% compared to the result using random KCA pattern. In addition, the optimization of KCA can also increase the peak-signal-to-noise-ratio (PSNR) of reconstructed images by 0.41dB.
The remainder of this paper is organized as follows. Section II describes the forward projection imaging model of the CSXT system with static KCA. Section III develops the nonlinear reconstruction algorithm for the proposed CSXT system. Section IV provides the optimization method of KCA. Section V presents the simulation results for the proposed methods. Conclusions and discussions are given in Section VI.
2. Forward projection model of CSXT
According to the spectral characterization of the Beer-Lambert law, as a mono-energetic X-ray beam with initial intensity $I_0(E)$ passes through the object, its intensity is attenuated as $I(E)=D(E)I_0(E)\exp \left (-\int _{l(r)}\mu \left (\vec {r},E\right )dr\right )$, where $D(E)$ is the spectral response function, $\mu \left (\vec {{r}},E\right )$ denotes the linear attenuation coefficient of the object, and $l(r)$ is the projection path that X-ray photons take to reach the detector element [37]. As shown in Fig. 2, a static KCA is placed in front of the X-ray source in CSXT, which plays a role in encoding the spectrum of the X-ray beam. The transmittance of the pixel in KCA is given by $f(E)=\exp \left (-t_f\mu _f(E)/\cos \psi \right )$, where $t_f$ and $\mu _f(E)$ are the thickness and linear attenuation coefficient of the K-edge material in this pixel, $\psi$ is the projection angle [24]. The filtered X-ray intensity acquired by the integrating detector is then given by $I=\int f(E)I(E)dE$. As previously mentioned, the linear attenuation coefficient of a K-edge filter has a sudden increase at its K-edge energy $E_f$. Suppose there are totally $F$ kinds of K-edge materials used for encoding, the X-ray spectrum can be divided into $K$ energy bins ($K=F+1$): $\left [E_{\min }, E_{1}\right ],\left [E_{1}, E_{2}\right ], \ldots$, and $\left [E_{F}, E_{\max }\right ]$, where $E_{\min }$ and $E_{\max }$ are the minimum and maximum photon energies of the X-ray source, respectively. Given that the width of each energy bin defined by the K-edge materials with nearly adjacent atomic numbers is small, the linear attenuation coefficients $\mu (\vec {r},E)$ and $\mu _f(E)$ at the $k\textrm {th}$ energy bin can thus be approximated by its effective attenuation coefficients denoted by $\bar \mu (\vec {r},E_k)$ and $\bar \mu _f(E_k)$, respectively. Subsequently, the measured intensity on detector can be rewritten as
Material decomposition is an important application of the spectral X-ray tomographic imaging [38,39]. Suppose the scanned object consists of $L$ known basis materials, the spectral attenuation of the object thus can be formulated as the weighted sum of the attenuation of the $L$ basis materials [40,41]:
3. Nonlinear reconstruction for CSXT
According to the forward mapping described in Eq. (5), the estimation of the tomographic images of the basis materials can be formulated as minimizing the unconstrained objective function $\frac {1}{2}\Vert \mathbf {y}-\mathcal {F}(\mathbf {X})\Vert ^2_2$. Considering the similarity of multi-layer tomographic images, this paper introduces the 3D total variation (TV) regularization per basis material to simultaneously preserve the edges and remove the noise in the images [31–33]. For a volumetric image denoted by the 3D tensor $\mathcal {X}\in \mathbb {R}^{N_x\times N_y\times N_z}$, the 3D TV norm is defined as [31]
Updating $\mathbf {X}$: Given $\mathbf {Z}$, the variable $\mathbf {X}$ can be updated using the Gauss-Newton (GN) method, which is a classical iterative method for nonlinear optimization [36]. It updates $\mathbf {X}$ by $\mathbf {X}^{(t+1)} = \mathbf {Z}^{(t)} + \alpha \mathbf {D}^{(t)}$, where $\mathbf {D}^{(t)}$ is the GN direction in the $t\textrm {th}$ iteration, and $\alpha$ is the step length. Let us denote $\mathcal {F}^{(t)}=\mathcal {F}\left (\mathbf {X}^{(t)}\right )$ and the residual $\mathbf {r}^{(t)}=\mathbf {y}-\mathcal {F}^{(t)}$ for simplicity. Then, the GN direction can be calculated by solving the following least square estimate problem:
Equation (11) is a high-dimensional linear problem, and thus cannot be efficiently solved by direct inversion. Here it is iteratively solved using the conjugate gradient (CG) method [36]. In the $\tau \textrm {th}$ CG iteration, the gradient of cost function Eq. (11) with respect to $\mathbf {D}^{(\tau )}_{(:,l)}$ is given by
Then, the descent direction $\mathbf {P}$ in the $\tau \textrm {th}$ iteration is determined by $\mathbf {P}^{(\tau )}=-\mathbf {G}^{(\tau )} + \beta ^{(\tau )}\mathbf {P}^{(\tau -1)}$, where $\beta ^{(\tau )}$ in the Flecther-Reeves formula is calculated by $\beta ^{(\tau )}_\textrm {FR} = {\Vert \mathbf {G}^{(\tau )}\Vert ^2_2}/{\Vert \mathbf {G}^{(\tau -1)}\Vert ^2_2}$ [42]. The step size $\alpha _\textrm {CG}^{(\tau )}$ of CG iteration can be obtained by a line-searching method along the descent direction, which is expressed as
Updating $\mathbf {Z}$: Given $\mathbf {X}$, the update of $\mathbf {Z}$ can be formulated as the following 3D TV denoising problem
This problem can be solved by applying the iterative clipping algorithm for each material in $\mathbf {X}$ [43].
4. KCA optimization for CSXT
As shown in Fig. 2(b), the proposed system includes opaque elements to block a partion of X-ray beam thereby reducing the overall X-ray dose in a scan. However, a random coded aperture generally yields nonuniform radiation dose on different portions of object. Some voxels are thus over-exposed and suffer more X-ray damage, which increases the risk for women undergoing regular screening mammography [28,29]. To avoid this, we propose to optimize the KCA based on the criterion of uniform sensing to minimize the overall radiation exposure placed on any volumetric area of the patient, i.e., every voxel is equally illuminated by each type of filtered X-ray beams. Other optimization criteria could be used but are left for future research [44,45]. The details of the optimization are presented as follows. As shown in Fig. 3, a KCA with $F$ types of K-edge materials can be formulated as a binary matrix $\mathbf {C}$ of size $F\times M_x M_y$, where the element $\mathbf {C}_{f,j}=1$ indicates that the $j\textrm {th}$ pixel in KCA is made of the $f\textrm {th}$ K-edge material. Evidently, the column-sum of $\mathbf {C}$ should be no more than 1 since only one type of K-edge material can be used in a specific KCA pixel. If all of the elements in one column of $\mathbf {C}$ are zero, the associated KCA pixel is blocked. Then, the radiation dose that different voxels in the object receive from the filtered X-ray beams can be measured by $\sum _{P=1}^P\mathbf {C}\mathbf {H}^p$. The average dose per voxel is $\bar \mu =D\sum ^P_{p=1}\mathbf {1}^\top _M\mathbf {H}^p\mathbf {1}_N/FNM$ when the total number of non-opaque pixels is limited to $D$, where $M=M_x M_y$. According to the criterion of uniform sensing, the optimization of KCA can be formulated as minimizing the variance of $\sum _{P=1}^P\mathbf {C}\mathbf {H}^p$. The problem is described as
Spatial coded aperture optimization: The subproblem of spatial coded aperture optimization can be derived from Eq. (15) by setting $F=1$, such that
Subsequently, Eq. (16) can be rewritten as
We solve Eq. (18) using the the Lagrangian multiplier algorithm [36]. The Lagrangian function for Eq. (18) is expressed as
Fig. 4. (a) Projection image of the patient's breast CT phantom on $0^\circ$ angular view; (b) original 30kV X-ray source spectrum, and the spectrum filtered with the material of Se, Sr, Nb, Ru, Ag, and In, respectively [18]
Spectral coded aperture optimization: Substituting the optimal spatial coded aperture $\mathbf {C}_\Lambda$ into Eq. (15), we obtain the problem for spectral coded aperture optimization, which is expressed as
Using Eq. (17) again to transform the binary constraint, the equivalent form of Eq. (20) is then given by
Similarly, this equality-constrained problem can be solved using the Lagrangian multiplier algorithm [36]. The Lagrangian function for Eq. (21) is given by
Fig. 5. Comparison of the material decomposition reconstruction for breast phantom. The first and fourth rows are the ground truth of adipose and the glandular in the tomographic images; the second and fifth rows are the proposed KCA-CSXT reconstruction results; the third and sixth rows are the conventional tomosynthesis reconstruction results using the proposed nonlinear algorithm. Columns from left to right correspond to the images at even layers.
5. Simulation results
This section will first studies the performance of the proposed CSXT imaging system. Since it is based on a single KCA, we also abbreviate it as KCA-CSXT. Then, we present the results of KCA optimization for uniform sensing. In the simulations, we use an X-ray tomosynthesis configuration with a flat 2D detector containing $M_x\times M_y=220\times 180$ pixels. The pitch size of each detector pixel is $1\textrm {mm}\times \textrm {mm}$. 11 projection images are taken from the angular range of $[-30^\circ, 30^\circ ]$ at $6^\circ$ intervals. The X-ray source rotates with the center of detector plane, and the distance from the X-ray source to the rotation center is 350mm. The distance from the object to the rotation center is 67.5mm. According to these parameters, we generate the system matrix $\mathbf {H}$ using the ASTRA software [47]. To evaluate the quality of the reconstructed images, we use the PSNR and the relative error (RE) as the metrics [37,48]. Given the reconstructed image $\hat {\mathbf {X}}\in \mathbb {R}^{N\times N}$ and the ground truth $\mathbf {X}$, they are respectively defined as $\operatorname {PSNR}=20\log _{10}\frac {\hat {\mathbf {X}}_\textrm {max}}{\Vert \mathbf {X}-\hat {\mathbf {X}}\Vert _\textrm {F}/N^2}$ and $\operatorname {RE} = \frac {\Vert \mathbf {X}-\hat {\mathbf {X}}\Vert _\textrm {F}}{\Vert {\mathbf {X}}\Vert _\textrm {F}}$.
5.1 Simulations on nonlinear reconstruction for CSXT
We illustrate the feasibility of the proposed KCA-CSXT method by applying it to mammographic imaging. Part of head Shepp-Logan phantom is extracted as the breast CT phantom used in our simulations [40], which is composed of the adipose and glandular (soft tissue), has a size of $128\textrm {mm}\times 128\textrm {mm}\times 50\textrm {mm}$, and is divided into $N_x\times N_y\times N_z=128\times 128\times 10$ voxels. Despite the nonsimilarity of the head phantom with the desired breast phantom in structure, the use of it is only to demonstrate the separation of those two materials in the reconstructed tomographic images. The projection image acquired at $0^\circ$ angular view is presented in Fig. 6(a). The tube voltage of the X-ray source is set to 30kV, whose spectrum is generated with the Spektr software [19]. To obtain the optimal reconstruction performance, a full sampling scheme shown in Fig. 2(a) is adopted in this simulation. Each pixel in KCA is randomly assigned to one of the six K-edge materials, including Selenium (Se), Strontium (Sr), Niobium (Nb), Ruthenium (Ru), Silver (Ag) and Indium (In). The corresponding K-edge energies are 12.66keV, 16.10keV, 18.99keV, 22.12keV, 25.51keV and 27.94keV. This choice of the K-edge materials allows the spectrum of the X-ray source to be divided into seven energy bins with approximately equal bandwidth of 2.5keV, which facilitates the uniform sensing of the spectral attenuation information throughout the entire spectrum. The thicknesses of the K-edge materials are respectively set as 20um, 20um, 10um, 10um, 10um and 10um to release approximately equal filtered X-ray intensity. Figure 6(b) shows the original and the filtered spectra, where the linear attenuation coefficients of the K-edge materials come from the National Institute of the Standards and Technology (NIST) X-ray attenuation database [18]. These spectra are discretized into seven energy bins determined by the materials’ K-edge energies. To obtain the spectral phantom, for each voxel in the breast CT phantom, we assign the spectral attenuation data according to its components. By substituting the discretized spectral attenuation of KCA and the spectral data cube of breast phantom into Eq. (5), we generate the KCA-CSXT measurement data used for simulation.
Fig. 6. Comparison of the reconstructed tomographic images for breast phantom. The first row is the ground truth of energy-binned tomographic images at the 6th energy bins; the second row is the KCA-CSXT reconstruction results using the proposed algorithm; the third row is the conventional tomosynthesis reconstruction results using the proposed nonlinear algorithm.
Figure 7 presents the reconstructed tomographic images of the two basis materials in the phantom. For comparison, we also present the reconstruction results by applying the nonlinear algorithm for the conventional tomosynthesis system without KCA, which is referred to as nonlinear tomosynthesis for simplicity. It should be pointed out that, different from KCA-CSXT, the blank intensity $\mathbf {W}$ in 4 for conventional tomosynthesis is only determined by the spectral response of the X-ray source and the detectors. The reconstruction results shown in Fig. 7 indicate that the KCA-CSXT system enables the separation of adipose and glandular tissue components in the reconstructed images. Meanwhile, in the reconstructed images for conventional tomosynthesis system, the component of adipose is reduced and misidentified as the glandular (soft tissue). Quantitative comparison on the average PSNR and RE values at all layers corresponding to these images are presented in Table 1. It is found that the material images reconstructed for KCA-CSXT have significantly lower reconstruction error when compared with the reconstruction results of conventional tomosynthesis.
Fig. 7. Quantitative comparison of nonlinear reconstruction for KCA-CSXT and conventional tomosynthesis. (a) The PSNR values of reconstructed tomographic images in all layers and all energy bins; (b) the RE values of reconstructed tomographic images in all layers and all energy bins.
Table 1. The average PSNR and RE values of the reconstructed basis material images.
The tomographic images of linear attenuation coefficients corresponding to the material images shown in Fig. 7 are presented in the second and third rows of Fig. 8. They are extracted from the $6\textrm {th}$ energy bin in the spectral phantom data cube. It is found that the high-quality tomographic images can be obtained from the proposed KCA-CSXT system. In addition, benefiting from the nonlinear reconstruction algorithm, the conventional tomosynthesis system also reconstructs the artifact-free CT images. However, it fails to separate the basis materials. Figure 9 presents the PSNR and RE values of the reconstructed tomographic images at all layers and energy bins in 3D surface form. It is found that the energy-binned tomographic images of KCA-CSXT have higher quality and lower error at most energy bins when compared with the results obtained from conventional tomosynthesis. The average PSNR value of reconstructed images at all layers and all energy bins is 38.48dB, which is 5.45dB higher when compared with the reconstruction results for conventional tomosynthesis system. The average RE value of these images is 1.60%, which is 1.76% lower when compared with the reconstruction results for conventional tomosynthesis system.
Fig. 8. Comparison of reconstruction with different settings for pixel binning in KCA. These images come from the sixth layer of the reconstructed 3D data cube.
Assume that the KCA is fixed in front of the X-ray source and the pixels of KCA are a one-by-one matching with the pixels on the detectors, the pixel sizes in KCA are then very small due to the projection amplification of the tomosynthesis system. To improve the manufacturability of KCA, we propose to bin the adjacent pixels together such that the binned pixels can be made of the same K-edge material. Simulations with different pixel binning areas of $2\times 2$, $4\times 4$, $8\times 8$ and $16 \times 16$ are implemented to demonstrate the feasibility. Table 2 presents the average PSNR and RE values of the reconstructed images at all layers. The corresponding images are shown in Fig. 10, which come from the sixth layer of the reconstructed data cube. It is found that the quality of reconstructed images gradually degrades as more pixels in KCA are binned together. Nevertheless, this degradation is acceptable since the pixel binning method facilitates the manufacturing of KCA.
Fig. 10. Frontal slice from the digital breast phantom. Each voxel in the phantom has been classified into one out four basis materials: adipose tissue, glandular tissue, skin, and air.
Table 2. The average PSNR and RE values of the reconstructed basis material images with different pixel binning settings.
The difficulty of manufacturing the KCA is also related to the required number of K-edge materials, We validate the performance of the proposed reconstruction method for the CSXT system with fewer types of K-edge materials. Additional two KCAs respectively made of [Sr, Nb, Ru, Ag] and [Nb, Ag] are used for simulations. Figure 11 presents the images that come from the sixth layer of the reconstructed 3D material data cube. Some images shown in Fig. 7 are presented again for comparison. Obviously, the CSXT system with zero type of K-edge material is equivalent to the conventional tomosynthesis system. It is found that the capability of material separation gradually declines as fewer types of K-edge materials are used in CSXT. We thus conclude that KCA should be made of as many types of materials as possible to obtain better reconstruction accuracy.
Fig. 11. Comparison of the material decomposition reconstruction for breast phantom derived from clinical breast images. The first and fourth rows are the ground truth of adipose and the glandular in the tomographic images; the second and fifth rows are the proposed KCA-CSXT reconstruction results; the third and sixth rows are the conventional tomosynthesis reconstruction results using the proposed nonlinear algorithm. Columns from left to right correspond to the images at even layers.
An additional digital breast phantom is derived from high-resolution clinical 3D breast images [49], which has been cropped and resized to adapt to the simulation setting. As shown in Fig. 12, each voxel in the phantom is classified into one of four main materials: glandular tissue, adipose tissue, skin, and air. Based on this segmentation, we generate the spectral phantom using the aforementioned method and then simulate the measurement data using the projection model. Figure 13 depicts the reconstructed tomographic images of the two basis material in the phantom. The images reconstructed from tomosynthesis system without KCA are also presented from comparison. It is found that the basis materials in KCA-CSXT images are with higher separation comparing with the conventional tomosynthesis images. Note that the region of skin in the phantom is reconstructed as glandular tissue since both of them are with the spectral attenuation characteristics similar to the soft tissue. Quantitative comparison presented in Table 3 also indicates that the material images reconstructed for KCA-CSXT have significantly lower reconstruction error when compared with the reconstruction results of conventional tomosynthesis.
Fig. 12. (a) The optimized spatial coded aperture obtained at the first stage of the proposed algorithm; (b) the optimized KCA using the proposed algorithm; (c) zoomed version of the part indicated by the red square in spatial coded aperture; (d) zoomed version of the part indicated by the red square in KCA.
Fig. 13. (a) X-ray radiation dose distribution in the fifth layer of the object when a random KCA is used; (b) X-ray radiation dose distribution in the fifth layer of the object when the optimized KCA is used.
Table 3. The average PSNR and RE values of the reconstructed basis material images for breast 3D phantom.
5.2 Simulations on optimization for KCA
In this subsection, we illustrate the performance of the proposed algorithm for optimizing the KCA in CSXT. The K-edge material set including Se, Sr, Nb, Ru, Ag and In is reused in this simulation. The transmittance of KCA is constrained within 50% to reduce the X-ray radiation dose. Fig. 12 presents both the optimized spatial coded aperture and the resultant KCA, which are respectively obtained at two stages of the proposed optimization algorithm. The zoomed images of the part indicated by the red squares in the coded apertures are also presented in Fig. 12. The black pixels in the spatial coded aperture represent the opaque pixels, while the white pixels represent the transparent pixels. The colorful pixels in KCA indicate that they are made of different K-edge materials shown in the legend. Given the optimized KCA $\mathbf {C}_\textrm {opt}$, we use the variance of $\sum _{p=1}^P\mathbf {C}_\textrm {opt}\mathbf {H}^p$ to measure the uniformity of radiation dose. For the six types of K-edge materials in the optimized KCA, the variances are 2.54, 2.54, 2.51, 2.50, 2.63 and 2.50 respectively. For comparison, we generate 10 different random KCAs and calculate the variances for the six types of K-edge materials again. The average variances for the random KCAs are 5.43, 5.68, 5.14, 5.52, 5.83 and 5.55 respectively, which means that the average uniformity is improved by 54.14% after optimization. Figure 13 shows the radiation dose distribution measured by the path length of X-ray beams in the fifth layer of the object, which comes from the X-ray beams filtered by the third type of K-edge material, i.e., the map of $\left (\mathcal {T}(\sum ^P_{p=1}{\mathbf {C}_{(3,:)}\mathbf {H}^p})\right )_{(:,:,5)}$. It is found that the optimized KCA provides a more uniform dose distribution in the object. In addition, using the CSXT with optimized KCA for 50% dose mammographic imaging, the average PSNR value of the reconstructed images is increased by 0.41dB compared to the use of a random KCA.
6. Conclusion and discussions
This paper develops a cost-effective spectral X-ray tomosynthesis system, which enables basis material decomposition in mammographic imaging by encoding the spectrum of X-ray beams with a static KCA. To reconstruct the tomographic images, a nonlinear alternating projection algorithm is proposed to solve the TV-regularized inverse imaging problem. Subsequently, an optimization algorithm of KCA is proposed to enable uniform distribution of X-ray radiation dose in the object. Simulations with two phantoms show that the proposed system enables the reconstruction of artifact-free tomographic images of both basis materials and linear attenuation coefficients of the object. In comparison, conventional tomosynthesis fails to identify the basis materials and merely reconstructs the attenuation images with inferior accuracy. Simulations also show that the uniformity of radiation dose is effectively improved when the optimized KCA is used. Besides, the image quality is slightly improved with the use of an optimized KCA. Future research will focus on developing novel criteria for KCA optimization, which enables high-efficiency sensing of the structural and spectral information in the object, thereby improves the nonlinear reconstruction accuracy of compressive spectral tomosynthesis.
Funding
National Key Research and Development Project under Grant (2019YFB2102300, 2019YFB2102301, 61936014); the Scientific Research Project of Shanghai Science and Technology Committee under Grant (19511103302); Fundamental Research Funds for the Central Universities; National Science Foundation (CIF 1717578); University Dissertation Award.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are available in [18,19].
References
1. D. G. Grant, “Tomosynthesis: a three-dimensional radiographic imaging technique,” IEEE Transactions on Biomed. Eng. pp. 20–28 (1972).
2. “Digital Breast Tomosynthesis (DBT) System,” https://www.accessdata.fda.gov/cdrh_docs/presentations/pghs/ Digital_Breast_Tomosynthesis_(DBT)_System.htm (2011). [Online; accessed 6-July-2021].
3. A. Ferrari, L. Bertolaccini, P. Solli, P. Salvia, and D. Scaradozzi, “Digital chest tomosynthesis: the 2017 updated review of an emerging application,” Ann. Transl. Med. 6(5), 91 (2018). [CrossRef]
4. L. T. Niklason, B. T. Christian, L. E. Niklason, D. B. Kopans, and R. F. Wirth, “Digital tomosynthesis in breast imaging,” Radiology 205(2), 399–406 (1997). [CrossRef]
5. J. T. Dobbins III and H. P. McAdams, “Chest tomosynthesis: technical principles and clinical update,” Eur. J. Radiol. 72(2), 244–251 (2009). [CrossRef]
6. A. Tingberg, “X-ray tomosynthesis: a review of its use for breast and chest imaging,” Radiat. Prot. Dosim. 139(1-3), 100–107 (2010). [CrossRef]
7. J. M. Park, E. A. Franken Jr, M. Garg, L. L. Fajardo, and L. T. Niklason, “Breast tomosynthesis: present considerations and future applications,” Radiographics 27(suppl_1), S231–S240 (2007). [CrossRef]
8. A. V. Bronnikov, “Cone-beam reconstruction by backprojection and filtering,” J. Opt. Soc. Am. A 17(11), 1993–2000 (2000). [CrossRef]
9. B. J. Heismann, B. T. Schmidt, and T. Flohr, Spectral computed tomography (SPIE, 2012).
10. D. Cabib, R. A. Buckwald, Z. Malik, Y. Garini, N. Katzir, and D. G. Soeknsen, “Spectral bio-imaging methods for biological research, medical diagnostics and therapy,” 1998. US Patent 5, 784, 162.
11. J. Chung, J. G. Nagy, and I. Sechopoulos, “Numerical algorithms for polyenergetic digital breast tomosynthesis reconstruction,” SIAM J. on Imaging Sci. 3(1), 133–152 (2010). [CrossRef]
12. W. R. Geiser, S. A. Einstein, and W.-T. Yang, “Artifacts in digital breast tomosynthesis,” AJR, Am. J. Roentgenol. 211(4), 926–932 (2018). [CrossRef]
13. F. F. Schmitzberger, E. M. Fallenberg, R. Lawaczeck, M. Hemmendorff, E. Moa, M. Danielsson, U. Bick, S. Diekmann, A. Pöllinger, F. J. Engelken, and F. Diekmann, “Development of low-dose photon-counting contrast-enhanced tomosynthesis with spectral imaging,” Radiology 259(2), 558–564 (2011). [CrossRef]
14. E. Fredenberg, M. Lundqvist, M. Åslund, M. Hemmendorff, B. Cederström, and M. Danielsson, “A photon-counting detector for dual-energy breast tomosynthesis,” in Medical Imaging 2009: Physics of Medical Imaging, vol. 7258 (International Society for Optics and Photonics, 2009).
15. K. Taguchi and J. S. Iwanczyk, “Vision 20/20: Single photon counting X-ray detectors in medical imaging,” Med. Phys. 40(10), 100901 (2013). [CrossRef]
16. J. Tanguay, H. K. Kim, and I. A. Cunningham, “The role of X-ray swank factor in energy-resolving photon-counting imaging,” Med. Phys. 37(12), 6205–6211 (2010). [CrossRef]
17. W. C. Barber, E. Nygard, J. S. Iwanczyk, M. Zhang, E. C. Frey, B. M. W. Tsui, J. C. Wessel, N. Malakhov, G. Wawrzyniak, N. E. Hartsough, T. Gandhi, and K. Taguchi, “Characterization of a novel photon counting detector for clinical CT: count rate, energy resolution, and noise performance,” in Medical Imaging 2009: Physics of Medical Imaging, vol. 7258 (2009), p. 725824.
18. M. Berger and S. Seltzer, “Radiation and biomolecular physics division, MPL NIST,” https://physics.nist.gov/PhysRefData/Xcom/html/xcom1-t.html (2008) [Online; accessed 19-July-2019].
19. J. Siewerdsen, A. Waese, D. Moseley, S. Richard, and D. Jaffray, “Spektr: A computational tool for X-ray spectral analysis and imaging system optimization,” Med. Phys. 31(11), 3057–3067 (2004). [CrossRef]
20. M. Saito, “Quasimonochromatic X-ray computed tomography by the balanced filter method using a conventional X-ray source: Quasimonochromatic X-ray CT by balanced filter method,” Med. Phys. 31(12), 3436–3443 (2004). [CrossRef]
21. Y. Rakvongthai, W. Worstell, G. El Fakhri, J. Bian, A. Lorsakul, and J. Ouyang, “Spectral CT using multiple balanced k-edge filters,” IEEE transactions on medical imaging 34(3), 740–747 (2015). [CrossRef]
22. P. Kirkpatrick, “On the theory and use of ross filters,” Rev. Sci. Instrum. 10(6), 186–191 (1939). [CrossRef]
23. P. Kirkpatrick, “Theory and use of ross filters. II,” Rev. Sci. Instrum. 15(9), 223–229 (1944). [CrossRef]
24. A. Cuadros, X. Ma, and G. R. Arce, “Compressive spectral X-ray tomography based on spatial and spectral coded illumination,” Opt. Express 27(8), 10745–10764 (2019). [CrossRef]
25. A. P. Cuadros, X. Ma, and G. R. Arce, “Compressive X-ray material decomposition using structured illumination,” in Developments in X-ray Tomography XII, vol. 11113 (International Society for Optics and Photonics, 2019), p. 111131H.
26. T. Zhang, S. Zhao, X. Ma, A. P. Cuadros, Q. Zhao, and G. R. Arce, “Nonlinear reconstruction of coded spectral X-ray CT based on material decomposition,” Opt. Express 29(13), 19319–19339 (2021). [CrossRef]
27. A. P. Cuadros, X. Ma, C. M. Restrepo, and G. R. Arce, “Staticcodect: single coded aperture tensorial X-ray CT,” Opt. Express 29(13), 20558–20576 (2021). [CrossRef]
28. R. Smith-Bindman, J. Lipson, R. Marcus, K. P. Kim, M. Mahesh, R. Gould, A. B. D. González, and D. L. Miglioretti, “Radiation dose associated with common computed tomography examinations and the associated lifetime attributable risk of cancer,” Arch. Intern. Med. 169(22), 2078–2086 (2009). [CrossRef]
29. R. E. Hendrick, “Radiation doses and cancer risks from breast imaging studies,” Radiology 257(1), 246–253 (2010). [CrossRef]
30. J. B. Weaver and A. L. Huddleston, “Attenuation coefficients of body tissues using principal-components analysis,” Med. Phys. 12(1), 40–45 (1985). [CrossRef]
31. M. Ertas, I. Yildirim, M. Kamasak, and A. Akan, “Digital breast tomosynthesis image reconstruction using 2D and 3D total variation minimization,” Biomed. engineering online 12(1), 112–212 (2013). [CrossRef]
32. L. Sun, Y. Zheng, and B. Jeon, “Hyperspectral restoration employing low rank and 3D total variation regularization,” in 2016 International Conference on Progress in Informatics and Computing (PIC) (IEEE, 2016), pp. 326–329.
33. A. Behrooz, H.-M. Zhou, A. A. Eftekhar, and A. Adibi, “Total variation regularization for 3D reconstruction in fluorescence tomography: experimental phantom studies,” Appl. Opt. 51(34), 8216–8227 (2012). [CrossRef]
34. X. Liao, H. Li, and L. Carin, “Generalized alternating projection for weighted-2, 1 minimization with applications to model-based compressive sensing,” SIAM J. on Imaging Sci. 7(2), 797–823 (2014). [CrossRef]
35. X. Yuan, “Generalized alternating projection based total variation minimization for compressive sensing,” in 2016 IEEE International Conference on Image Processing (ICIP) (IEEE, 2016), pp. 2539–2543.
36. J. Nocedal and S. Wright, Numerical optimization (Springer Science & Business Media, 2006).
37. L. Li, Z. Chen, W. Cong, and G. Wang, “Spectral CT modeling and reconstruction with hybrid detectors in dynamic-threshold-based counting and integrating modes,” IEEE transactions on medical imaging 34(3), 716–728 (2015). [CrossRef]
38. Y. Long and J. A. Fessler, “Multi-material decomposition using statistical image reconstruction for spectral CT,” IEEE transactions on medical imaging 33(8), 1614–1626 (2014). [CrossRef]
39. K. C. Zimmerman and T. G. Schmidt, “Experimental comparison of empirical material decomposition methods for spectral CT,” Phys. Med. Biol. 60(8), 3175–3191 (2015). [CrossRef]
40. G. Richard Hammerstein, D. W. Miller, D. R. White, M. Ellen Masterson, H. Q. Woodard, and J. S. Laughlin, “Absorbed radiation dose in mammography,” Radiology 130(2), 485–491 (1979). [CrossRef]
41. J. Hubbell, “Photon cross sections, attenuation coefficients, energy absorption coefficients,” Natl. Bureau Standards Rep. NSRDS-NBS29 (1969)
42. Y. Dai and Y.-x. Yuan, “Convergence properties of the fletcher-reeves method,” IMA J. Numer. Analysis 16(2), 155–164 (1996). [CrossRef]
43. A. Chambolle, “An algorithm for total variation minimization and applications,” Journal of Mathematical imaging and vision 20(1/2), 89–97 (2004). [CrossRef]
44. T. Mao, A. P. Cuadros, X. Ma, W. He, Q. Chen, and G. R. Arce, “Coded aperture optimization in X-ray tomography via sparse principal component analysis,” IEEE Transactions on Comput. Imaging 6, 73–86 (2020). [CrossRef]
45. T. Mao, A. P. Cuadros, X. Ma, W. He, Q. Chen, and G. R. Arce, “Fast optimization of coded apertures in X-ray computed tomography,” Opt. Express 26(19), 24461–24478 (2018). [CrossRef]
46. G. Yuan and B. Ghanem, “Binary optimization via mathematical programming with equilibrium constraints,” arXiv preprint arXiv:1608.04425 (2016).
47. W. van Aarle, W. J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. D. Beenhouwer, K. J. Batenburg, and J. Sijbers, “Fast and flexible X-ray tomography using the astra toolbox,” Opt. Express 24(22), 25129–25147 (2016). [CrossRef]
48. C. Lu, J. Feng, Y. Chen, W. Liu, Z. Lin, and S. Yan, “Tensor robust principal component analysis with a new tensor nuclear norm,” IEEE transactions on pattern analysis machine intelligence 42(4), 925–938 (2020). [CrossRef]
49. A. Sarno, G. Mettivier, F. di Franco, A. Varallo, K. Bliznakova, A. M. Hernandez, J. M. Boone, and P. Russo, “Dataset of patient-derived digital breast phantoms for in silico studies in breast computed tomography, digital breast tomosynthesis, and digital mammography,” Med. Phys. 48(5), 2682–2693 (2021). [CrossRef]
