Abstract
X-ray grating interferometry offers a novel framework for the study of weakly absorbing samples. Three kinds of information, that is, the attenuation, differential phase contrast (DPC), and dark-field images, can be obtained after a single scanning, providing additional and complementary information to the conventional attenuation image. Phase shifts of X-rays are measured by the DPC method; hence, DPC-CT reconstructs refraction indexes rather than attenuation coefficients. In this work, we propose an explicit filtering based low-dose differential phase reconstruction algorithm, which enables reconstruction from reduced scanning without artifacts. The algorithm adopts a differential algebraic reconstruction technique (DART) with the explicit filtering based sparse regularization rather than the commonly used total variation (TV) method. Both the numerical simulation and the biological sample experiment demonstrate the feasibility of the proposed algorithm.
1. Introduction
X-ray grating interferometry [1–4] with conventional X-ray tubes develops rapidly in recent years and is becoming the most promising technology among various phase contrast imaging methods for clinical applications. Three kinds of information, that is, the attenuation, differential phase contrast (DPC), and dark-field images, can be obtained through one single scanning, and the latter two images provide additional and complementary information to the conventional attenuation image. The DPC method measures phase shifts of X-rays by obtaining the line integral of the directional derivatives of refractive index decrements (δ) [5], that is, the refraction angle of the beam. The refraction index reconstructed afterwards is 1000 times larger than the absorption index.
The phase-stepping approach of the grating interferometry, which requires a number of images to retrieve information, significantly increases the examine time and the dose delivered to the patient [6]. The problem becomes even more severe for DPC-CT because of the requirement of multiangle scanning. Therefore, reducing the number of projections, the exposure time, and the delivered dose is of great value. And that is why the low dose DPC reconstruction algorithm is proposed.
As mentioned above, the reconstruction problem for DPC-CT is to obtain the refraction index from the refraction angle data. The analytical method, the filtered backprojection (FBP) algorithm with the Hilbert transform, was first applied [7, 8]. Afterwards, several iterative algorithms, such as the maximum likelihood (ML) algorithm [9] and the differential algebraic reconstruction technique (DART) [10], were proposed. However, these algorithms rely on the completeness of data and the large number of projections.
The recently proposed compressed sensing (CS) theory [11] makes image reconstruction from incomplete data possible. Essentially, it illustrates that if the image is sparse in a domain which has small coherence with the sampling domain, according to the Shannon/Nyquist sampling theorem, fewer projections can almost accurately recover the images. A typical image reconstruction method exploits TV as the sparse regularization [12] (from CS measurements). Applications in both the absorption imaging [12] and the DPC imaging [10] have been implemented. Instead of the implicit regularization coming from the penalty, another sparse regularization method based on explicit filtering is proposed, which exploits spatially adaptive filters sensitive to image features and details [13]. However, no similar algorithms for DPC imaging have been suggested so far, which is just the problem to be solved in this paper.
In this work, we propose an explicit filtering based low-dose differential phase reconstruction algorithm. The algorithm combines the DART iterative algorithm and the explicit filtering based CS method. It has the potential to exactly reconstruct the refractive index distribution using few-view projections, thus reducing the exposure time and the delivered dose, making DPC-CT closer to clinical applications. The feasibility of the low dose reconstruction algorithm is verified by both the numerical simulation and the biological sample experiments.
2. Methods
2.1. Grating-Based Imaging
Figure 1 illustrates the schematic diagram of a typical grating interferometry. Two kinds of apparatuses are shown, the Talbot effect based interferometry with coherent source, that is, Figure 1(a), and the Talbot-Lau effect based interferometry with incoherent source, that is, Figure 1(b). The first grating G1 creates its self-image through Talbot-Lau effect or classical optics in the position of G2 where Moire fringes occur. The source grating G0 splits the source into an array of line sources, enabling the use of the large-focal-spot X-ray tube, that is, the incoherent source. The phase-stepping approach is adopted for image acquisition, capturing a series of raw images at every step of one of the gratings along the transverse direction, obtaining the intensity oscillation curve, Figure 1(c). The changes of the background oscillation curves determine three kinds of information, namely, the attenuation image, the DPC image, and the dark-field image.
(a)
(b)
(c)
To analyze the changes quantitatively, the oscillation curve for each pixel is expressed by the Fourier expansion series: where is the amplitude coefficient, is the corresponding phase coefficient, and is the period of G2. Then, the attenuation, dark-field, and DPC images are given by , , and , respectively, where the superscripts () and () denote the values with the sample in place and as a reference without, respectively, and is the visibility of the oscillation curve.
2.2. Differential Phase-Contrast Reconstruction Algorithm
The DPC image measured by the grating interferometer is the refraction angle, which is related to the phase shift of the sample:
In other words, the grating interferometer obtains the line integral of the directional derivatives of refraction index decrements. Therefore, the reconstruction problem of DPC-CT can be expressed by where is the refraction angle projection, is the refraction index decrement of the samples, is the path of X-ray beam in the medium and is the perpendicular direction to .
By contrast, the projection of the conventional X-ray transmission comes from the linear integration of attenuation coefficient; hence, the reconstruction problem can be expressed by where is the intensity projection and is the linear attenuation coefficient.
Such difference in mathematical expressions of the projections requires different reconstruction algorithms.
Inspired by the widely used algebraic reconstruction technique (ART) [14], Wang et al. proposed a differential algebraic reconstruction technique (DART) by discretizing the projection process of the differential phase contrast imaging into a linear partial derivative matrix [10]: where denotes the net interpolation coefficient corresponding to each pixel.
Therefore, the forward projection process of differential phase contrast imaging can be expressed by where is named as the linearly partial-derivative matrix. Equation 3 can be used to reconstruct the refractive index by algorithms similar to ART directly shown as 4: where is the image vector in the th iteration; presents the forward projection process.
However, the DART algorithm relies on the completeness of data [10]. It is incompetent at the ill-posed reconstruction problems, such as the cases of few-view or limited-angle projections.
2.3. Explicit Filtering Based Low-Dose Differential Phase-Contrast Reconstruction Algorithm
In this part, we propose an explicit filtering based low-dose differential phase-contrast reconstruction algorithm based on the above DART algorithm for DPC reconstruction.
Generally, the compressed sensing DPC-CT reconstruction method can be summarized as
Here, the fitness function , which tries to match the estimation to the data, is accomplished by the DART algorithm above, where the operation represents the -2 norm. The regularization function expresses some priori known property of the unknown object, which usually defines a sparse representation of after a specific transformation. The regularization parameter balances the fitness function with the regularization function. Different regularization methods are used to find the solution of the mathematic model under different constraints of minimizing the complexity of reconstructed signals in different representations. The typical TV-type regularization method uses the total variations as the sparse transformation, which is a parametric regression technique. Another solution is to replace the parametric regression by spatially adaptive filters, which is sensitive to image features and details. In this paper, the BM3D filter is adopted, that is, DART_BM3D algorithm.
The BM3D algorithm is based on an enhanced sparse representation in transform domain [15]. The enhancement of the sparsity is achieved by grouping similar 2D fragments of the image into 3D data arrays. And then the collaborative filtering with the hard thresholding is carried out, which enhances the similarity between the blocks while at the same time preserves even the finest details with their essential unique features shared by the jointly filtered 2D fragments. The implementation of BM3D filter is shown as Figure 2.
BM3D filter matches the similar block regions in the image (block matching, BM) using the similarity defined as where is the numerical metric of the image, is the side length of block, and is the difference between two blocks.
Blocks with small differences can be grouped as a set . For blocks in the same group, a hard threshold is used to assign zeros for small value pixels to generate a weight for blocks: where is the number of nonzero pixels. Then, the filtered image can be obtained as a similarity-weighted average:
The process of our reconstruction algorithm is shown in Figure 3 and the steps are described as follows.(A)Initialization. Creating the linearly partial-derivative matrix for few-view scanning based on the geometry parameters. An initial guess of the reconstruction image is given and .(B)DART reconstruction: where DART represents the reconstruction method in 4.(C)BM3D filter: where BM3D represents the method shown as in Figure 2.(D)Add excitation noise into and go back to step (B).
Gaussian noise is added in the unobserved portions, that is, the missing angles caused by the few-view or limited-angle scanning, in frequency domain, which works as a random generator of the missing components in the spectrum: where represents the Fourier transform operation and is the inverse Fourier transform operation, is the standard deviation of the Gaussian noise which is determined empirically by the noise of the projections caused by the quantum noise and electronics noise.
3. Experiments and Results
The proposed explicit filtering based low-dose differential phase-contrast reconstruction algorithm was validated by both the numerical simulation and biological experiments.
3.1. Numerical Simulation
The Shepp-Logan phantom with the resolution of 256 × 256 was used in the study. The phantom is shown as in Figure 4(a) and the reconstructed results are shown as in Figures 4(b)–4(g). All images are shown in the same display window .
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figures 4(b) and 4(c) are the reconstructed results with the FBP and DART (500 iterations) methods, respectively, with 180 views within 180° with 1° angular interval, which verified the effectiveness of the FBP and DART methods in case of complete data. Figures 4(d) and 4(e) are the reconstructed results with the FBP and DART, respectively, methods with 10 views within 180° with 18° angular interval, with an angular downsampling factor of 18. The results demonstrate the dependency on the completeness of the data of the FBP and DART (2000 iterations) methods.
Figure 4(f) is the reconstructed result with the typical TV-type compressed sensing method (2000 iterations). The method is an effective reconstruction method dealing with incomplete data. Figure 4(g) is the reconstructed results with proposed explicit filtering based compressed sensing methods (2000 iterations), which is also skilled in dealing with incomplete data. The reconstructed results with the compressed sensing methods are almost artifact-free and they are in high accordance with the phantom in both appearance and values, while the results of FBP and DART algorithms show severe streaking artifacts due to the downsampling. Furthermore, profiles as shown in Figure 4(h) show that the proposed method is better than the TV-type method in preserving details.
3.2. Biological Sample Experiment
The biological sample results have been used to test the proposed reconstruction method. The experiments were performed at the TOMCAT beamline using a two-grating interferometer operated at 25 KeV and in the 3rd Talbot distance at the Swiss Light Source of the Paul Scherrer Institute in Switzerland. The pitch of the phase grating was 3.981 m with a height of m. The corresponding values for the second grating (gold absorber grating) were m and m. The sample was a guinea pig eyeball packed in a plastic pipe and was scanned within 180° with equivalent angular interval of 1°. For each view, an eight-step phase stepping process was adopted and the refraction angular projections were retrieved by information retrieving algorithm.
The reconstruction results are shown in Figure 5. Figures 5(a) and 5(b) are the results of FBP algorithm and DART algorithm with 180 views, respectively. Since the reconstruction image is more complex than the phantom in the previous section, an angular down-sampling factor of 9 was adopted. Figures 5(c) and 5(d) is the results of FBP algorithm and DART algorithm with 20 views, respectively. The results show severe streaking artifacts, especially the FBP result. As shown in Figures 5(e) and 5(f), the results of the two compressed sensing methods are nearly artifact-free. Both the results of the two compressed sensing methods show great image quality improvement compared with DART and FBP in few-view reconstruction. And the red ellipses in Figures 5(e) and 5(f) illustrate the better capabilities of preserving details of the proposed explicit filtering based method than the TV-type method, which can also be seen from Figure 6, which shows the profiles of the 188th line in the red circles on the right.
(a)
(b)
(c)
(d)
(e)
(f)
4. Conclusions
In this paper, an explicit filtering based low-dose differential phase-contrast reconstruction algorithm is proposed. The algorithm is an application of the compressed sensing theory and has the potential to accurately reconstruct the distribution of the refractive index with few-view projections.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors thank Marco Stampanoni for the help with biological sample experiments. The work was supported by a grant from the National Natural Science Foundation of China (NNSFC), Grants 11235007, and a Tsinghua University Independent Research Project Grant, “Research on Key Technologies and CT Reconstruction Methods of Multienergy X-Ray imaging”.