Implementation of a fast method for high resolution phase contrast tomography

We report the implementation of a method which can yield the 3D distribution of the phase (refractive index) of a weakly absorbing object from a single tomographic data set. In order to reduce the residual absorption artifact (due to the fact that only one data set is used) the original algorithm presented by A. V. Bronnikov is amended by adding in the filter a factor found by using a semi empirical approach. The quality of the reconstruction is largely sufficient for optimal segmentation and further postprocessing even though the filter correction is based on assumption of constant absorption. This one step approach allows keeping radiation dose to the minimum. Spatial resolution is comparable to the conventional absorption based technique. The performance of the method is validated by using an established phase contrast technique. ©2006 Optical Society of America OCIS codes: (340.7440) X-ray Imaging; (170.6960) Tomography (100.5070) Phase retrieval References and links 1. U. Bonse and F. Bush, “X-ray computed microtomography using synchrotron radiation,” Prog. Biophys. Mol. Biol., 65, 133-169 (1996). 2. B. A. Dowd, G. H. Campbell, R. B. Marr, V. V. Nagarkar, S. V. Tipnis, L. Axe and D. P. Siddons, “Developments in synchrotron x-ray computed microtomography at the National Synchrotron Light Source,” in Developments in X-Ray Tomography II, U. Bonse, ed., Proc. SPIE 3772, 224-236 (1999). 3. M. Stampanoni, G. Borchert, P. Wyss, R. Abela, B. Patterson, S. Hunt, D. Vermeulen and P. Rüegsegger, “High resolution X-ray detector for synchrotron-based microtomography,” Nucl Instrum Methods A 491, 291-301 (2002). 4. T. Weitkamp, C. Raven, and A. Snigirev, “An imaging and microtomography facility at the ESRF beamline ID 22,” in Developments in X-Ray Tomography II, U. Bonse, ed., Proc. SPIE 3772, 311-317 (1999). 5. U. Bonse and M. Hart, “An x-ray interferometer,” Appl. Phys. Lett. 6, 155-156 (1965). 6. U. Bonse, and F. Beckmann, “Multiple-beam X-ray interferometry for phase-contrast microtomography,” Journal of Synchrotron Radiat. 8, 1-5 (2001). 7. A. Momose, “Demonstration of phase-contrast x-ray computed-tomography using an x-ray interferometer,” Nuclear Instruments and Methods A 352, 622-628 (1995). 8. T. Weitkamp, A. Diaz, C. David, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “X-ray phase imaging with a grating interferometer,” Opt. Express 13, 6296-6304 (2005). 9. M. R. Teague, “Deterministic phase retrieval a Green-function solution,” J. of the Opt. Soc. Am. 73, 1434– 1441 (1983). 10. K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard x rays,” Phys. Rev. Lett. 77, 2961-2964 (1996). 11. L. J. Allen, W. McBride, N. L. O’Leary, and M. P. Oxley, “Exit wave reconstruction at atomic resolution,” Ultramicroscopy 100, 91-104 (2004). 12. W. K. Hsieh, F. R. Chen, J. J .Kai, A. I. Kirkland, “Resolution extension and exit wave reconstruction in complex HREM,” Ultramicroscopy 98, 99-114 (2004). 13. S. C. Mayo, P. R. Miller, S. W. Wilkins, T. J. Davis, D. Gao, T. E. Gureyev, D. Paganin, D. J. Parry, A. Pogany and A. W. Stevenson, “Quantitative X-ray projection microscopy: phase-contrast and multi-spectral imaging,” J. Microsc. 207, 79-96 (2002). 14. T. E. Gureyev, A. Pogany, D. M. Paganin and S. W. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region,” Optics Commun. 231, 53-70 (2004). 15. A. G. Peele, F. DeCarlo, P. J. Mahon, B. B. Dhal and K. A. Nugent, “X-ray phase contrast tomography with a bending magnet source,” Rev. Sci. Instrum. 76, 083707 (2005). #72197 $15.00 USD Received 20 June 2006; revised 21 August 2006; accepted 23 August 2006 (C) 2006 OSA 4 September 2006 / Vol. 14, No. 18 / OPTICS EXPRESS 8103 16. T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation II. Orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A 13, 1670-1682 (1996). 17. P. Cloetens, W. Ludwig, J. Baruchel, D. Van Dyck, J. Van Landuyt, J. P. Guigay and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. Lett. 75, 2912-2914 (1999). 18. U. Bonse, F. Beckmann, M. Bartscher, T. Biermann, F. Busch, and O. Günnewig, “Phase contrast x-ray tomography using synchrotron radiation” in Developments in X-Ray Tomography, U. Bonse, ed., Proc. SPIE 3149, 108-119 (1997). 19. A. V. Bronnikov, “Reconstruction formulas in phase-contrast tomography,” Optics Communications 171, 239-244 (1999). 20. A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. of the Opt. Soc. of Am. A 19, 472-480 (2002). 21. A. Groso, M. Stampanoni, R. Abela, P. Schneider, S. Linga, and R. Müller, “Phase contrast tomography: an alternative approach,” Appl. Phys. Lett. 88, 214104 (2006). 22. The computation of the two-dimensional Radon transform and backprojection over the angle ω are combined in a single step of the algorithm [20]. 23. J. M. Cowley, Diffraction Physics (North-Holland, 1995), Chap. 3.4. 24. A. Pogany, D. Gao, and S. W. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instr. 68, 2774-2782 (1997).


Introduction
Absorption based x-ray computed tomography (CT) is nowadays a standard technique with spatial resolution around one micron [1][2][3][4]. On the other hand, a wide range of samples studied in materials science, biology and medicine show very weak absorption contrast, yet produce significant phase shifts of the x-ray beam. The use of phase information can therefore provide new and otherwise not accessible information. Phase contrast imaging techniques have two key advantages: first, light elements -showing poor contrast in absorption radiography -can be easily detected and second, phase-contrast radiography helps to reduce the radiation dose deposited on the object under investigation.
As far as phase tomography (quantitative 3D reconstruction of the phase or the refractive index from 2D phase images) is concerned, several attempts were made by using interferometric [5][6][7][8] and non-interferometric phase retrieval methods [9][10][11][12][13][14][15][16][17]. All these methods are based on a two step approach: first, the projections of the phase are determined in the form of Radon projections and then the object function, i.e. the refractive index δ is reconstructed applying a conventional filtered backprojection algorithm. Since far from absorption edges δ is linearly related to the electron density ρ, which is -except for hydrogen rich materials -proportional to the mass density, the resulting reconstruction represents approximately the distribution of the mass density in the sample [18].
The reconstruction algorithm suggested by Bronnikov in his theoretical works [19,20] presents an alternative approach which requires no intermediate step of 2D phase retrieval and provides a direct 3D reconstruction of the refractive index from the intensity distributions (radiographic projections at different angles) measured in a single plane (pure phase object) and from the intensity distributions measured in two planes in the case of mixed phaseamplitude object.
In a previous work [21] we have simulated and validated the quantitative phase tomography method described above. We have adapted the original algorithm in order to get 3D reconstruction of the refractive index from the intensity distributions measured in a single plane also in the case of mixed weakly absorbing-phase object. When using this one step approach in phase reconstruction in 3D the experimental setup can be kept extremely simple and the radiation damage is kept to the minimum, which is of special importance for biological specimens. Besides, phase retrieval as well as data alignment (necessary in case of measurements in multiple planes) are superfluous.
In this work we describe in detail the filter adaptation and the means of choosing the filter parameters. The method is applied to real experimental data and also validated by the Differential Phase Contrast Imaging Technique [8].

Theoretical background
When a monochromatic plane wave propagates along the positive z-axis and impinges upon a thin mixed phase-amplitude object (characterized by the linear absorption coefficient ) , , ( z y x μ and the real part of the index of refraction ) , , ( z y x δ ), the intensity distribution at distance z = d and angle of rotation θ is approximated by the following expression, also known as Transport of Intensity equation [9][10][11][12][13][14][15][16]: is the phase function of the object. This is valid in the near field Fresnel region, where d << a 2 /λ (a is the transversal size of the smallest structure in the object and λ is the x-ray wavelength). For a mixed phase amplitude object with weak and almost (1) simplifies further to [19,20]: . The goal of quantitative phase tomography is to reformulate Eq. (2) to , applying then the 3D Radon transform (denoted by the symbol ∧) and calculating finally the second derivative with respect to the variable ω ω cos sin y x s + = one gets: Expression (3) is a theorem which states that from the 2D Radon transform of the measured value g, one can directly find the 3D Radon transform of δ. An explicit inversion formula for the 3D Radon transform was given by Radon and Lorentz already at the beginning of the previous century, and one therefore has a direct solution to the reconstruction problem for quantitative phase contrast CT [22]: where the stars indicate a 2D convolution and This convolution integral can be computed in the Fourier domain by taking the two dimensional Fourier transform. In the Fourier domain, Eq. (5) has the low-pass filter form given by [20]: Expressed in this (filtered backprojection) form [Eq. (4)], this reconstruction approach becomes very interesting since in case of a pure phase object ( 1 ) it allows recovering the refractive index in 3D from only one single tomographic data set: This is a significant improvement since (a) the experimental setup can be kept extremely simple (it is actually the same as for standard, absorption based tomography) and (b) the radiation damage is kept to the minimum, which is of special importance for biological specimens and for potential future in vivo studies. In addition, time consuming data alignment routines as well as phase retrieval processes become superfluous.

Results and analysis
In a previous work [21], we have validated the theory by implementing the aforesaid reconstruction algorithm and successfully recovered the refractive index of the simulated object. In the second step, the original algorithm [19,20] has been used to reconstruct phase tomograms from experimental data sets, obtained at the Tomography Station of the Materials Science beamline at the Swiss Light Source [3]. The reconstructed images using one distance data set have shown artifacts due to residual absorption (see section 3.1). We have found a solution to this problem by modifying the original algorithm (section 3.2).

Experimental results
As already pointed out, the experiment setup is equivalent to the one for standard absorption tomography. Nevertheless, two additional conditions have to be satisfied: i) The sample-detector distance (z) is now increased from the minimum, (z=0) used in absorption tomography, to the near field Fresnel region (z << a 2 /λ).
ii) Photon energy is increased to maximum attainable (with sufficient flux) at the experimental station (25 keV) in order to satisfy, as much as possible, the weak and almost homogeneous After traversing the sample, the x-rays are converted into visible light by a 20 microns thick YAG scintillator and then magnified (10× or 20×) by an optical microscope and detected with a 14 bit, 2048 × 2048 pixel CCD camera, giving an effective pixel size of 1.4 and 0.7 microns respectively.
As a test, we investigated the internal microstructure of 350 microns thin wood samples. Even when low photon energy (10 keV) is used there is a weak absorption contrast and only the edges of internal structures are visible [see Fig. 1(a)].
In order to apply the phase tomography method in a pure phase object approximation, a single data set (I θ,z=d (x,y)) of 501 angular projections at sample detector distance of 32 mm and 25 keV photon energy was acquired. The nonuniform illumination intensity was normalized by dividing through reference (sample-out of the beam) images. For the practical implementation of the filter in Eq. (6) a sufficiently small sampling interval has to be used in the frequency domain. This can be achieved by zero padding the discrete data of the experimental function g θ (x,y). This is a standard technique to avoid the wraparound error in digital signal processing. Fig. 1(b) presents a slice through the sample reconstructed using Eq. (6) and Eq. (4).

Analysis
The reconstructed slice shown in Fig. 1(b) is severely corrupted by the residual absorption artifact (see cap in the center). This means that, even though the absorption level for the (wood) sample is calculated to be in the range of only 2 %, pure phase object condition is not satisfied and an intensity measurement at 'zero' distance (I θ,z=0 (x,y)) is also required. We have adopted another approach -the reconstruction artifacts are corrected by amending the original method [20] leading to what we named Modified version of Bronnikov's Algorithm (MBA).
The correction consists of adding, in the denominator of the low-pass filter given in Eq.
(6), an absorption dependent correction factor α exp . Consequently, the new filter has the following form: The values of α exp to be used are found by using a semi empirical (simulations-experiment) approach. First, we performed simulations using a mathematical phantom made up of three spheres of different size inside a cube. Two cases, a pure phase object (β=0, where β is the imaginary part of the refractive index) and mixed, weakly absorbing phase object (β of the order of 10 -9 ) were studied. The phase function φ was calculated analytically from the knowledge of δ and sets of phase contrast projection data (at different rotation angles) were generated using the Fresnel propagator [23,24] with λ = 1 Å, sample detector distance z = 25 mm, 185 2 pixels detector array and an isotropic pixel size of 1 ) the phase information (refractive index) has been successfully reconstructed with errors around 1%. For mixed phase-amplitude object, reconstruction using two distances data sets ( leads to similar result. If only a data set of angular projections at one distance is used for the reconstruction of a simulated mixed object, one gets a cap similar to the one shown in Fig.  1(b). The right correction factor in case of simulations (α simul ) has been obtained by the iterative minimization of the sum square error (SSE) defined as: The achieved values for SSE were of the order of 10 -4 . The same iterative procedure was applied for different values (in the range of 1-10 %) of absorption and the corresponding correction factors are shown in Fig. 2. For a specific experimental data set, the correction factor α exp to be used in the filter of Eq. (8) is then calculated in the following way: Average absorption/transmission level of the object is measured/calculated based on a 'zero' distance intensity projection. The corresponding α simul is found from the graph depicted in Fig.   2.
Based on α simul , the correction factor α exp is found according to the following expression: Where pixel simul and array simul are the pixel size and the detector array size (in pixels) used for simulations and pixel exp and array exp are the corresponding values for a specific experimental data set. Applying the MBA method, we reconstructed the same data set, for the wood sample. The result is presented in Fig. 3. It appears clearly that the absorption artifact has completely disappeared.

Validation
In order to validate our approach, we have compared our results with the ones obtained using an established phase contrast method, the Differential Phase Contrast (DPC) Technique [8].
To make a direct evaluation, experiments using the two approaches have been performed on the same samples and with the same detector effective pixel size (1.4 microns). The used photon energies were different, optimized for the individual techniques (17.5 keV for the DPC and 25 keV for the MBA). performed without taking detector and beam characteristics into account. However, as well shown in Fig. 3(a) and Fig. 4(b), the quality of the reconstruction is largely sufficient to perform optimal segmentation and further post processing.
The MBA was validated using the Differential Phase Contrast Imaging technique [8]. For the particular example shown in the Fig. 4 the contrast obtained using MBA is even better. In general, the two methods are complementary. The full 3D phase imaging approach (MBA) is experimentally simple, very fast and it is ideally suited for small objects when resolution around 1 micron is needed. Deposited dose is equivalent (or lower, since higher photon energy is used) to the one deposited in conventional tomography. The MBA spatial resolution is limited, equivalently to standard absorption tomography, by the detector resolution to around 1 micron [3]. The DPC is more demanding in terms of instrumentation and acquisition time (therefore also deposited dose) but it is more sensitive and can be scaled up to large fields of view. Spatial resolution is limited to two period of analyzer grating [8] e.g. to 2 microns with present fabrication technology.
Future developments will focus on the fine tuning of the filtering step and on the streamlining of the whole process in order to allow high-resolution and high throughput imaging of soft biological and materials science specimens.