## Abstract

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

## 1. Introduction

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 [12] and out-of-plane CT [13] 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 $n=1-\delta $. A two-step algorithm first reconstructs a map of the gradient of the refractive index decrement $\nabla \delta $; this gradient map is then integrated to obtain a distribution of refractive index decrement [12]. 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 decrement$\nabla \delta $will be denoted by a g*δ*-map. Note that g*δ*-maps consist of two maps: $\frac{\partial \delta}{\partial x}$ and $\frac{\partial \delta}{\partial y}$.

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: ${N}_{p}\approx \pi /2\cdot {N}_{b}$, where ${N}_{b}$ and ${N}_{p}$ are the number of bins in a single projection and the number of projections, respectively [16]. 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. [17] 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 [18] 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.

## 2. Method

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 [19] we derived the following equation that each differential projection should satisfy. From the ray equation, the fundamental equation in geometrical optics:

*l*denotes a line, $x\mathrm{sin}\theta +y\mathrm{cos}\theta =p$, corresponding to an incident beam trajectory; ${r}_{l}=\left(p\mathrm{cos}\theta -q\mathrm{sin}\theta ,p\mathrm{sin}\theta +q\mathrm{cos}\theta \right)$; $\nabla \delta \left({r}_{l}\right)=\left(\frac{\partial \delta \left({r}_{l}\right)}{\partial x},\frac{\partial \delta \left({r}_{l}\right)}{\partial y}\right)$; $\varphi ({r}_{l})$ is a phase of the vector $\nabla \delta $ at ${r}_{l}$; $\Delta \alpha \left(p\right)$ is an in-plane angular deviation at the beam position

*p*. Noting that the integrand in the left hand side is a polar coordinate form of $\left(\raisebox{1ex}{$\partial $}\!\left/ \!\raisebox{-1ex}{$\partial x$}\right.+i\raisebox{1ex}{$\partial $}\!\left/ \!\raisebox{-1ex}{$\partial y$}\right.\right)\delta $, the left hand side represents a line integral of a complex-valued representation of $\nabla \delta $ along the

*q*-axis, i.e., the beam direction. In other words, it corresponds to a Radon transform of $\nabla \delta $. On the other hand, the values in the right hand side can be known by the measurement. Therefore, we can obtain g

*δ*-maps, $\nabla \delta (r)$, by solving the inverse problem of Eq. (1). In [12], 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:

*l*,

*i.e.*, the beam propagation, respectively. Thus, we can obtain the g

*δ*-maps of $\left(\frac{\partial \delta \left(r\right)}{\partial x},\frac{\partial \delta \left(r\right)}{\partial y}\right)$ by reconstructing 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 $\frac{\partial \delta \left(r\right)}{\partial x}$ and $\frac{\partial \delta \left(r\right)}{\partial y}$ as $\xi \left(m,n\right)$ and $\eta \left(m,n\right)$ $\left(1\le m,n\le K\right)$, respectively; denote *N*-dimensional vectors representing $\xi \left(m,n\right)$ and $\eta \left(m,n\right)$ in the lexicographic ordering as $\xi ={({\xi}_{1},{\xi}_{2},\cdots ,{\xi}_{N})}^{t}$ and $\eta ={({\eta}_{1},{\eta}_{2},\cdots ,{\eta}_{N})}^{t}$, respectively, where *N* = *K* × *K*, ${\xi}_{mK+n}=\xi \left(m,n\right)$, and ${\eta}_{mK+n}=\eta \left(m,n\right)$. Then, we discretize the Eqs. (2) according to [16] as follows:

*M*is the total number of rays; the weighting factor ${w}_{ij}$, which is equal to the line segment of the

*j*th pixel intercepted by the

*i*th ray, represents the contribution of the

*j*th pixel to the

*i*th ray integral. Equations (3) can be expressed as $W\text{\hspace{0.17em}}\xi =h$and $W\text{\hspace{0.17em}}\eta =v$, where

*W*is an

*M*×

*N*matrix such that ${w}_{ij}$ is the (

*i*,

*j*)-th element; $h$ and $v$ are M-dimensional vectors such that ${h}_{i}$ and ${v}_{i}$ are the

*i*th component of $h$ and $v$, respectively.

We formulate the reconstruction for $\text{\xi}$ and $\text{\eta}$ with TV regularization as follows:

*λ*controls the relative weights of the data fidelity and regularization terms. In order to obtain reconstructed differential images, we apply the method proposed by [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 [16]. 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 [21], in which a sufficiently small constant step size was used.

Next, we will estimate a *δ*-map, $\delta (r)$, from the two tomographic images $\text{\xi}$ and $\text{\eta}$ estimated for g*δ*-maps, $\frac{\partial \delta \left(r\right)}{\partial x}$ and $\frac{\partial \delta \left(r\right)}{\partial y}$, in the previous stage. If $\text{\xi}$ or $\text{\eta}$ is noise-free, we can obtain $\delta (r)$ by integrating either $\text{\xi}$ or $\text{\eta}$ 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 $\delta (r)$ by solving a Poisson equation, which was proposed Gasilov et al [22]. 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:

*x*and

*y*, respectively, and then adding them, we get:

## 3. Experiments and results

#### 3.1 Simulation

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 mm^{2}. 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 $691\text{\hspace{0.17em}}\left(\approx \raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$2$}\right.\times 440\right)$ 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) [23] 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.

#### 3.2 Experiment

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 [24]. 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 [25]. 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 $6296\text{\hspace{0.17em}}\left(\approx \raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$2$}\right.\times 4008\right)$ 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.

## 4. Discussion

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 [19]:

*s*is an arc length along the ray; $r=r\left(s\right)$ is the ray trajectory; $t\left(r\right)$ is a unit vector tangential to the ray trajectory. From the equation, the following two equations are derived:

*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 [15]. Therefore, the number of equations that $\delta (r)$ 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 [26].

In the current research literature, *l*_{1} regularization is attracting considerable attention in signal and image processing communities. In this paper, we have not considered compressed sensing using *l*_{1} regularization because this research was originally motivated by the work of Candes [17]. Not surprisingly, we could also reduce the views by replacing the TV regularization with *l*_{1} 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 [22]. Using this property, one could potentially reduce the number of projections even more. Additional research is needed to fully explore this topic.

## 5. Conclusion

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.

## Acknowledgment

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

**1. **A. Momose, T. Takeda, Y. Itai, and K. Hirano, “Phase-contrast X-ray computed tomography for observing biological soft tissues,” Nat. Med. **2**(4), 473–475 (1996). [CrossRef] [PubMed]

**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]

**14. **G. W. Faris and R. L. Byer, “Three-dimensional beam-deflection optical tomography of a supersonic jet,” Appl. Opt. **27**(24), 5202–5212 (1988). [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]