This paper describes an algebraic reconstruction algorithm that uses total variation (TV) regularization for differential phase contrast computed tomography (DPC-CT) using a limited number of views. In order to overcome over-flattening inherent in TV regularization, a two-step reconstruction process is used: we first reconstruct tomographic images of gradient refractive index from differential projections with TV regularization; these images are then used to compute tomographic images of refractive index by solving the Poisson equation. We incorporate TV regularization in the reconstruction process because the distribution of gradient refractive index is much more flattened than the refractive index. Simulations of the proposed method demonstrate that it can achieve satisfactory image quality from a much smaller number of projections than is required by the Nyquist sampling theorem. We experimentally prove the feasibility of the proposed method using dark field imaging optics at PF-14C beamline at the Photon Factory, KEK. The differential phase contrast projection data was experimentally acquired from a biological sample and DPC-CT images were reconstructed. We show that far fewer projections are needed when the proposed algorithm is used.
© 2015 Optical Society of America
Phase contrast imaging (PCI) exhibits excellent soft tissue contrast and can delineate pathological conditions that may be missed by conventional X-ray imaging [1–7]. Four types of PCI method have been described in the literature: analyzer-based imaging, grating-based imaging, interferometry-based imaging, and inline holography [1, 8–11]. Among these, the first two methods have been widely used for tomographic reconstructions of biological samples because of their relative simplicity and ease with which projection images can be acquired. Since acquired projection data from the systems are sum of the refractive index along the incident beam direction, differentiated in a direction perpendicular to the incident direction, they are termed as differential projections. Tomographic image reconstruction from such projections is generically termed as the differential phase-contrast computed tomography (DPC-CT). DPC-CT can be divided into two main classes: in-plane CT  and out-of-plane CT  that correspond to reconstructed tomographic images representing the transverse and the longitudinal components of the three dimensional vector field of the gradient refractive index, respectively. Most DPC-CT systems are in-plane type because data processing for obtaining the three dimensional distribution of refractive index is relatively easier than that of the out-of-plane type. In this paper, we restrict our attention to in-plane type DPC-CT.
Reconstruction algorithms for in-plane CT can be further categorized into one-step or two-step reconstruction algorithms. A one-step algorithm directly reconstructs a map of refractive index decrement δ from the differential projections [14, 15]; here δ is related to the refractive index by the equation . A two-step algorithm first reconstructs a map of the gradient of the refractive index decrement ; this gradient map is then integrated to obtain a distribution of refractive index decrement . One-step algorithm has been extensively used because of its simplicity of data processing. In this paper, we will denote a distribution of refractive index decrement δ by a δ-map; similarly, a map of the gradient of the refractive index decrementwill be denoted by a gδ-map. Note that gδ-maps consist of two maps: and .
DPC-CT reconstruction algorithms, whether one-step or two-step, typically use the filtered backprojection (FBP) method. FBP algorithm has also been widely used in medical and industrial CT scanners because of the ease with which it can be implemented and parallelized. It is also based on sound mathematical theory and its properties have well studied in the literature. FBP algorithm, however, is derived from a continuous model while the measurements are performed in a discrete domain. As a result, in order to obtain an artifact-free reconstructed image, the CT measurements must satisfy the sampling theorem: , where and are the number of bins in a single projection and the number of projections, respectively . For example, if the number of detector elements in the detector is to 4000, as many as 6283 projections are required to satisfy the sampling theorem. In addition, if the system has low photon counts, long exposure times may be required in order to accumulated adequate number of counts in each projection. Large number of projections and long accumulation time for each projection result in long measurement times for FBP overall. It was, therefore, quite dramatic when Candes et al.  demonstrated that a 512x512 Shepp-Logan phantom can be satisfactorily reconstructed from merely 22 projections. Instead of FBP, these researchers used method called the total variation minimization.
Total variation (TV) regularization was introduced by Rudin et al. for denoising in  and has become one of the most popular and successful methodologies in image processing. It has been extensively used for tasks such as image restoration and image deconvolution because of its excellent edge preservative nature. However, the total variation regularization suffers from a well-known drawback: it has a tendency to over-flatten the target image and impart is a cartoon-like appearance. This method, therefore, is satisfactory for a piecewise constant image such as the Shepp-Logan phantom. However, for most biological images that have texture, TV regularization may impart an artificial look to the reconstructed image.
If we can minimize the over-flattening caused by TV regularization, while using its strong points—namely, its edge preservative nature and emancipation from the sampling theorem—we can reconstruct an image from limited number of views and greatly shorten the measurement time. For this purpose, leverage a key observation: the gδ-map for a typical image is considerably flatter than the δ-map for the same image. An example of this phenomenon is shown in Fig. 1. In this figure, Fig. 1(a) represents the phase image (i.e., the δ-map) for an in-vitro human coronary artery specimen while Figs. 1(b) and 1(c) represent the gδ-maps with the derivative performed with respect to x and y, respectively. As can be seen, the differential images have a relatively flat dynamic range. Therefore, application of TV regularization is more acceptable for the gδ-map than for the δ-map.
In this paper, we develop a two-step reconstruction algorithm demonstrate its effectiveness of using both simulated and experimental data. The experimental data was acquired from a biological specimen dark field imaging (DFI) optics, a form of analyzer-based phase contrast imaging technique.
DPC-CT reconstructs a three-dimensional image of refractive index decrement (δ) by separately reconstructing individual two-dimensional image planes (δ-maps) perpendicular to a rotational axis of a sample. A stack of these two-dimensional images provides the 3D refractive index map of the sample. Thus, the overall 3D reconstruction problem can be reduced to a set of 2D CT reconstruction problems.
Consider an ordinary CT geometry, that is, a pq-coordinate system which is obtained by rotating an xy-coordinate system around the rotational axis by θ, in the reconstruction plane fixed to the sample, such that the q-axis is oriented toward a direction of incident beam propagation. In  we derived the following equation that each differential projection should satisfy. From the ray equation, the fundamental equation in geometrical optics:Eq. (1). In , we obtained a δ-map, by first reconstructing gδ-map using the FBP method and then integrating the gδ-map.
We replace the FBP method with the algebraic reconstruction technique (ART) with TV regularization. Dividing Eq. (1) into the real and imaginary parts, we get:Eqs. (2) separately. Reconstruction from a limited number of views is expectable by incorporating TV regularization into the reconstruction process.
Denote digital images estimated for and as and , respectively; denote N-dimensional vectors representing and in the lexicographic ordering as and , respectively, where N = K × K, , and . Then, we discretize the Eqs. (2) according to  as follows:Eqs. (2), respectively; M is the total number of rays; the weighting factor , which is equal to the line segment of the jth pixel intercepted by the ith ray, represents the contribution of the jth pixel to the ith ray integral. Equations (3) can be expressed as and , where W is an M × N matrix such that is the (i, j)-th element; and are M-dimensional vectors such that and are the ith component of and , respectively.
We formulate the reconstruction for and with TV regularization as follows:20]. The algorithm obtains a solution from the minimization problems (4) by alternating the ART (Algebraic Reconstruction Technique) steps with TV minimization steps. In the ART step, the standard ART algorithm reconstructs an image by iterations . Then, the image reconstructed in the ART step is fed to the TV steps as an initial value. In the TV step, the objective functions of Eqs. (4) are minimized using the standard steepest descendent method with a derivatives of TV term , in which a sufficiently small constant step size was used.
Next, we will estimate a δ-map, , from the two tomographic images and estimated for gδ-maps, and , in the previous stage. If or is noise-free, we can obtain by integrating either or in a horizontal or vertical direction. However, in noisy conditions the integration deteriorates the quality of resultant image because of accumulation error. Here, we adopted a method to obtain by solving a Poisson equation, which was proposed Gasilov et al . In this research, we did not consider correlation between a reconstruction plane of interest and its adjacent reconstruction planes, while Gasilov et al. incorporated the correlation into the estimation process.
From the results in the reconstruction process, we obtain the following differential equations:Eqs. (5) with respect to x and y, respectively, and then adding them, we get:21]. Figure 2 shows the flow chart of the whole reconstruction procedure.
3. Experiments and results
We developed a simulator for DPC-CT reflecting physical parameters in the actual experiment based on geometrical optics. Each ray was refracted according to Snell’s law at each boundary encountered in the reconstruction plane. In order to evaluate the influence of number of projections on the image quality, initially we did not model measurement noise and absorption by the sample. A numerical phantom was a 20-mm-diameter circle including an 18 relatively small circles with a constant refractive index decrement of 1.0 × 10−7. The values of refractive index decrement in the regions inside the large circle and outside the small circles varied from 0.0 to 1.0 × 10−7 to simulate texture-like regions, as shown in Fig. 3; These values of refractive index decrements are typical for a biomedical sample at the incident energy of 30 keV; the value of refractive index decrement in air region was zero. The number of bins in a single projection was 440; the matrix size of the reconstructed image was 440 × 440; the pixel dimension was 0.05 × 0.05 mm2. The differential projections were generated at constant angular steps over a 180 degree span.
In order to investigate the effectiveness of our proposed two-step reconstruction method, we compared it to the one-step reconstruction method. We simulated two kinds of one-step methods: the first method reconstructed tomographic images by applying the FBP, the second method applied ART with TV regularization to the convolution of the differential projections with the signum function [14, 15]. The flow charts for these methods are schematically represented in Fig. 4. We refer to them as FBP and δ-ART + TV, respectively. We denote the proposed method by gδ-ART + TV. Comparison between gδ-ART + TV and FBP reveals the robustness for under-sampled projections, and comparison between gδ-ART + TV and δ-ART + TV reveals the advantage of using the differential images.
First, we investigated the convergence property of gδ-ART + TV. Figure 5 shows the results of 5 cases with different number of projections. The root mean square error (RMSE) between the numerical phantom and the reconstructed images is plotted against the number of iterations for different number of projections. As can be seen, the RMSE converges to values less than 10−8 over 40 iterations. The convergence is slow when only 18 and 36 projections are used. This exercise confirms that the proposed method exhibits a good convergence property.
Figure 6 shows reconstructed images for the 5 cases considered above with FBP, δ-ART + TV, and the gδ -ART + TV. Note that for 440 detector bins, as was the case here, the sampling theorem requires projections. From subjective evaluation of the image quality depicted in this figure, FBP offers a satisfactory reconstruction for 120 projections; However, for 60 projections, which is about 1/11 of the required value, considerable streak artifacts appear in the image because of view starvation. In contrast, δ-ART + TV and gδ-ART + TV produce acceptable image quality even for 60 projections. In fact, gδ -ART + TV can reconstruct an acceptable image even from 36 projections. Figure 7 shows zoomed regions that correspond to the region indicated by a square in Fig. 6. As can be seen, δ-ART + TV generates over-flattened regions, the so-called cartoon appearance. It can be concluded from these images that the quality of gδ-ART + TV is superior to that of δ-ART + TV.
Finally, we quantitatively evaluated the quality of reconstructed images. Figures 8 and 9 show the dependence of the structural similarity (SSIM)  and RMSE on the number of projections, respectively. SSIM and RMSE of FBP abruptly begin to decrease around 600 projections because the sampling theorem violated. On the other hand, δ-ART + TV and gδ-ART + TV maintain image quality even when only 120 projections are available. In addition, gδ-ART + TV keeps a value of 0.9 even with only 60 projections. The behavior of SSIM and RMSE is concordant with the subjective evaluation in Fig. 7.
We investigated the efficacy of the proposed algorithm using actual projection data obtained from a DPC-CT system based on the DFI optics (Fig. 10). This system uses a Laue case angular analyzer with the condition that forward diffraction intensity diminishes when an incident angle corresponds with its Bragg angle, i.e., dark field condition . A DPC-CT imaging system with DFI optics consists of a monochro-collimator, a rotational sample stage, an angular analyzer, and two CCD cameras, as shown in Fig. 10 . The monochromated and parallel x-rays are expanded into a volumetric parallel beam by the monochro-collimator and incident on a sample. Each ray is refracted by a sample. The refracted beam impinges on the analyzer, which produces two directional beams, i.e., forward diffraction and diffraction beams. These beams, whose intensities are individually modified according to their corresponding rocking curves, are simultaneously acquired by their corresponding CCD cameras. See [15, 25] for a detailed system description.
The sample was a human coronary artery exercised at autopsy from a patient with atherosclerosis. This specimen, 13 mm in diameter 40 mm in length, was put in a plastic container filled with formalin solution and 600 projections were acquired over 180 degrees. The measurement time of a single projection was 6s and the total measurement time was 1h. Since the number of detector elements in a horizontal direction was 4008, according to sampling theorem approximately projections should be obtained. While 600 projections are far fewer than the theoretical requirement, this is the maximum number that could be acquired with in the limited machine time that was available.
Figure 11 shows reconstructed images for different number of projections. Subjectively, gδ-ART + TV provides good image quality even for 200 projections, while the image is blurred for 100 projections. In contrast, for FBP and δ-ART + TV the images are blurred even for 600 projections. Figure 12 show zoomed images of the region indicated by a square in Fig. 11. These images indicate that gδ-ART + TV is superior to FBP and δ-ART + TV in robustness and image quality from a limited number of views.
Our results indicate that a two-step reconstruction algorithm with TV regularization is more effective than the one-step reconstruction with TV regularization. This advantage is primarily attributable to the fact that gδ-map are more flattened than the δ-map and TV regularization exhibits better performance in the differential image. A second reason for the superiority of the two-step reconstruction algorithm is a bit less intuitive and is explained below. The two-step reconstruction originates from the ray equation :Equation (1) was derived from the two equations in Eqs. (8) by using the paraxial approximation and then introducing the CT geometry . The one-step reconstruction algorithm, i.e., reconstruction after convolving the differential projections with the signum function, can be derived from only the lower equation without use of the upper equation . Therefore, the number of equations that should satisfy for the two-step reconstruction is double that for the one-step reconstruction. As the number of projections decreases, the two-step reconstruction is advantageous in the estimation accuracy because it imposes additional constraints that are not a part of the one-step reconstruction.
In this research, we employed classical ART and steepest descend method for optimization of reconstruction problem (4). Although in simulation the proposed method exhibited a good convergence property, we will be able to improve the convergence rate by using a more sophisticated optimization technique such as FISTA .
In the current research literature, l1 regularization is attracting considerable attention in signal and image processing communities. In this paper, we have not considered compressed sensing using l1 regularization because this research was originally motivated by the work of Candes . Not surprisingly, we could also reduce the views by replacing the TV regularization with l1 regularization. Compressed sensing requires sparseness. Also in this case the two-step approach will be more effective than the one-step approach, since gδ-maps are sparser than δ-map. A compressed sensing based approach is a topic of current on going research.
We obtained a three-dimensional CT image by forming a two-dimensional cross sectional images in individual reconstruction planes separately and then stacking them up. However, the adjacent reconstruction planes strongly correlate with each other . Using this property, one could potentially reduce the number of projections even more. Additional research is needed to fully explore this topic.
We proposed an algebraic reconstruction method with TV regularization for in-plane DPC-CT. In order to mitigate the effect of over-flattening inherent in TV regularization, the proposed algorithm first reconstructs the differential images with respect to x and y using TV regularization; the distribution of refractive index decrement is then obtained by solving the Poisson equation with the source term being the divergence of the differential images under the relevant boundary condition. From the simulation results, it was proved that the proposed method can reconstruct an adequate quality image from much smaller number of projections than mandated by the sampling theorem. Additionally, it was seen that the proposed two-step reconstruction avoids the over-flattening and exhibits a satisfactory performance compared with the one-step reconstruction. For the actual experimental data, the proposed algorithm provided adequate image quality with merely 200 projections or only about 3% of projections required by the sampling theorem.
This research was partially supported by a Grant-In-Aid for Scientific Research (#23602002, #26350494, #26750142) from the Japanese Ministry of Education Science and Culture, JGC-S scholarship foundation, Sumitomo Foundation, and Yazaki Memorial Foundation for Science and Technology. The experiment was performed under the auspices of KEK (2014G002).
References and links
2. M. N. Wernick, O. Wirjadi, D. Chapman, Z. Zhong, N. P. Galatsanos, Y. Yang, J. G. Brankov, O. Oltulu, M. A. Anastasio, and C. Muehleman, “Multiple-image radiography,” Phys. Med. Biol. 48(23), 3875–3895 (2003). [CrossRef] [PubMed]
3. M. J. Kitchen, K. M. Pavlov, S. B. Hooper, D. J. Vine, K. K. W. Siu, M. J. Wallace, M. L. L. Siew, N. Yagi, K. Uesugi, and R. A. Lewis, “Simultaneous acquisition of dual analyser-based phase contrast X-ray images for small animal imaging,” Eur. J. Radiol. 68(3Suppl), S49–S53 (2008). [CrossRef] [PubMed]
4. P. Coan, J. Mollenhauer, A. Wagner, C. Muehleman, and A. Bravin, “Analyzer-based imaging technique in tomography of cartilage and metal implants: a study at the ESRF,” Eur. J. Radiol. 68(3Suppl), S41–S48 (2008). [CrossRef] [PubMed]
5. S. Ichihara, M. Ando, E. Hashimoto, A. Maksimenko, H. Sugiyama, C. Ohbayashi, T. Yuasa, K. Yamasaki, Y. Arai, K. Mori, and T. Endo, “3-D reconstruction of high-grade ductal carcinoma in situ of the breast with casting type calcifications using refraction-based X-ray CT,” Virchows Arch. 452(1), 41–47 (2008). [CrossRef] [PubMed]
6. Y. Zhao, E. Brun, P. Coan, Z. Huang, A. Sztrókay, P. C. Diemoz, S. Liebhardt, A. Mittone, S. Gasilov, J. Miao, and A. Bravin, “High-resolution, low-dose phase contrast X-ray tomography for 3D diagnosis of human breast cancers,” Proc. Natl. Acad. Sci. U.S.A. 109(45), 18290–18294 (2012). [CrossRef] [PubMed]
7. M. Ando, N. Sunaguchi, Y. Wu, S. Do, Y. Sung, A. Louissaint, T. Yuasa, S. Ichihara, and R. Gupta, “Crystal analyser-based X-ray phase contrast imaging in the dark field: implementation and evaluation using excised tissue specimens,” Eur. Radiol. 24(2), 423–433 (2014). [CrossRef] [PubMed]
8. D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmür, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced x-ray imaging,” Phys. Med. Biol. 42(11), 2015–2025 (1997). [CrossRef] [PubMed]
9. P. Suortti, J. Keyrilainen, and W. Thomlinson, “Analyser-based x-ray imaging for biomedical research,” J. Phys. D Appl. Phys. 46(49), 494002 (2013). [CrossRef]
10. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of X-ray talbot interferometry,” Jpn. J. Appl. Phys. 42(2), L866–L868 (2003). [CrossRef]
11. A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of x-ray phase contrast microimaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. 66(12), 5486 (1995). [CrossRef]
12. A. Maksimenko, M. Ando, H. Sugiyama, and T. Yuasa, “Computed tomographic reconstruction based on x-ray refraction contrast,” Appl. Phys. Lett. 86(12), 124105 (2005). [CrossRef]
13. F. A. Dilmanian, Z. Zhong, B. Ren, X. Y. Wu, L. D. Chapman, I. Orion, and W. C. Thomlinson, “Computed tomography of x-ray index of refraction using the diffraction enhanced imaging method,” Phys. Med. Biol. 45(4), 933–946 (2000). [CrossRef] [PubMed]
15. N. Sunaguchi, T. Yuasa, Q. Huo, and M. Ando, “Convolution reconstruction algorithm for refraction-contrast computed tomography using a Laue-case analyzer for dark-field imaging,” Opt. Lett. 36(3), 391–393 (2011). [CrossRef] [PubMed]
16. A. C. Kak and M. Slaney, Principles of Computed Tomographic Imaging (IEEE, 1988).
17. E. J. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52(2), 489–509 (2006). [CrossRef]
18. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D 60(1–4), 259–268 (1992). [CrossRef]
19. T. Yuasa, A. Maksimenko, E. Hashimoto, H. Sugiyama, K. Hyodo, T. Akatsuka, and M. Ando, “Hard-x-ray region tomographic reconstruction of the refractive-index gradient vector field: Imaging principles and comparisons with diffraction-enhanced-imaging-based computed tomography,” Opt. Lett. 31(12), 1818–1820 (2006). [CrossRef] [PubMed]
20. G. H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” Med. Phys. 35(2), 660–663 (2008). [CrossRef] [PubMed]
21. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd Edition (Cambridge University, 1992).
22. S. Gasilov, A. Mittone, E. Brun, A. Bravin, S. Grandl, A. Mirone, and P. Coan, “Tomographic reconstruction of the refractive index with hard X-rays: an efficient method based on the gradient vector-field approach,” Opt. Express 22(5), 5216–5227 (2014). [CrossRef] [PubMed]
23. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. 13(4), 600–612 (2004). [CrossRef] [PubMed]
24. M. Ando, A. Maksimenko, H. Sugiyama, W. Pattanasiriwissawa, K. Hyodo, and C. Uyama, “Simple X-ray dark- and bright-field imaging using achromatic Laue optics,” Jpn. J. Appl. Phys. 41(2), 1016 (2002). [CrossRef]
25. N. Sunaguchi, T. Yuasa, Q. Huo, S. Ichihara, and M. Ando, “X-ray refraction-contrast computed tomography images using dark field imaging optics,” Appl. Phys. Lett. 97(15), 153701 (2010). [CrossRef]
26. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]