REVIEW OF SCIENTIFIC INSTRUMENTS 85, 013704 (2014)
Improved algorithm for processing grating-based phase contrast interferometry image sets Shashidhara Marathe,1,a) Lahsen Assoufid,1,b) Xianghui Xiao,1,c) Kyungmin Ham,2,d) Warren W. Johnson,3,e) and Leslie G. Butler4,f) 1
Advanced Photon Source, Argonne National Laboratory, 9700 S. Cass Ave., Argonne, Illinois 60439, USA CAMD, LSU, Baton Rouge, Louisiana 70806, USA 3 Department of Physics and Astronomy, LSU, Baton Rouge, Louisiana 70803, USA 4 Department of Chemistry, LSU, Baton Rouge, Louisiana 70803, USA 2
(Received 10 December 2013; accepted 19 December 2013; published online 17 January 2014) Grating-based X-ray and neutron interferometry tomography using phase-stepping methods generates large data sets. An improved algorithm is presented for solving for the parameters to calculate transmissions, differential phase contrast, and dark-field images. The method takes advantage of the vectorization inherent in high-level languages such as Mathematica and MATLAB and can solve a 16 × 1k × 1k data set in less than a second. In addition, the algorithm can function with partial data sets. This is demonstrated with processing of a 16-step grating data set with partial use of the original data chosen without any restriction. Also, we have calculated the reduced chi-square for the fit and notice the effect of grating support structural elements upon the differential phase contrast image and have explored expanded basis set representations to mitigate the impact. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4861199] I. INTRODUCTION
Grating-based interferometry has been implemented at many X-ray and neutron imaging beamlines worldwide to perform phase contrast imaging and tomography with applications in materials science, biomedical science, and beam wavefront and optics characterization.1–10 These measurement tools generate large data sets that are processed using specialized computer codes. A variety of algorithms have been developed for this purpose. A discrete Fourier transform algorithm1 is quite popular, is fast, yet can suffer from spectral leakage when the grating steps are not uniformly spaced within a period. A Levenberg-Marquardt algorithm can solve non-uniformly spaced data, yet is slow. In this paper, we describe the application of an algorithm from the gravity wave toolbox that provides fast and robust processing. The toolbox was originally developed to rapidly process the large data volume that is typically generated by the interferometers used to search for gravity waves caused by events such as the collision of two black holes.11 The X-ray interferometry data were collected with a prototype experimental setup at the Advanced Photon Source. A grating-based phase-stepping interferometry data set consists of a set of pixelated images acquired as a function of a grating stepping. The triangular interferometric signal13, 14 is approximated as a sinusoid for which we wish to determine the amplitude and phase; the sinusoid period is usually known from the grating motion. The speed of this algorithm is based on pre-computation of the sinusoidal functions and a) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected] d) Electronic mail: [email protected] e) Electronic mail: [email protected] f) Electronic mail: [email protected]
0034-6748/2014/85(1)/013704/6/$30.00
vectorization of the pixel data. The amplitude and phase are then calculated with standard matrix manipulations. Two samples are studied: The first sample is made of polystyrene spheres supported by a Kapton film; reconstructions are performed with complete and partial data sets. The second sample is a calcium-carbonate-rich shell of a foraminifera, a one-cell protist, which shows a highly structured, millimeter-sized prolate object in a tomographic study. The reconstructed absorption and differential phase images demonstrate the applicability of this new algorithm. II. THEORY
The interferometry data are a set of X-ray counts cgp , where g = 1, . . . , M is the index that identifies the exposure, with one exposure at each grating displacement xg ; and p = 1, . . . , N is the index that identifies the pth pixel. Typically M ∼ 7 to 16 and N ∼ millions. We have adapted a fitting algorithm from mathematical physics (section 14-6 in Mathews and Walker15 ) that transforms the fitting operation into a simple matrix problem. We fit each set of exposures pixel by pixel, cgp , to each pixel’s expected dependence on grating position xg : 2π cˆgp = a1p + ap sin xg + φp (1) pgrat 2π xg ap cos(φp ) ≡ [1] a1p + sin pgrat 2π xg ap sin(φp ) (2) + cos pgrat ≡
3
Bgμ aμp
(3)
μ=1
85, 013704-1
© 2014 AIP Publishing LLC
013704-2
Marathe et al.
Rev. Sci. Instrum. 85, 013704 (2014)
a2p ≡ ap cos(φp ), a3p ≡ ap sin(φp ),
(4)
where ap and φ p are the amplitude and phase of the interference term for the pth pixel, and pgrat is the period of the translated grating. The interference term is expressed in both polar form (Eq. (1)) and Cartesian form (Eq. (4)). The latter is used to define the M × 3 matrix B (Eq. (3)) that represents the three fitting functions shown in brackets in Eq. (2): constant, sine, and cosine. The aμp is one element in the coefficient matrix a, which is a 3 × N matrix of amplitudes, or weights, of the three fitting functions for the N pixels. Later, we will reshape a into a 3D matrix with dimensionality [rows, columns, 3]. The best fit can be chosen to be the set of coefficients aμp that minimize the deviation-squared, Dp , for each pixel, defined by M (cgp − cˆgp )2 . Dp ≡
(5)
g=1
To minimize Dp , we need only calculate the derivatives of each deviation-squared with respect to each component of a, set all deviations to zero, and solve the resulting matrix equations. The closed form expression for the coefficient matrix a is found to be a = G · c,
(6)
G = (BT · B)−1 · BT ,
(7)
where
and where the superscript T indicates the matrix transpose. So the optimization problem is reduced to multiplying a 3 × M fixed matrix G into the M × N data matrix c, to find the 3 × N coefficient matrix a. With an efficient matrix manipulation program like MATLAB, the multiplication is extremely fast, of the order of 1 s for 1k × 1k images times 16 grating steps. By inspection, we recover the polar coefficients from the Cartesian in the usual way: 2 2 + a3p , φp = tan−1 (a3p /a2p ). (8) ap = a2p We have extended this analysis to use knowledge of the ex2 pected uncertainty or variance of each datum σgp at the price of increased computing time. The model equation, Eq. (4), can be extended to include anharmonic terms to describe the grating effects more accurately and to include more sophisticated grating motion trajectories. The Appendix has more details about the derivation of Eq. (7). In the stepped-grating interferometry experiment, two data sets are measured: one without the sample in the beam (reference data set) and the other with the sample in the beam (sample data set): transmission =
a1p (sample) , a1p (reference)
III. EXPERIMENTAL
A two-grating interferometer was assembled at the Advanced Photon Source 2-BM-B beamline (Fig. 1).16 The stepped-grating system was installed 23 m downstream of the double multilayer monochromator, which provides a monochromatic beam at 25 keV with 1% bandwidth. The phase grating, G1 , was optimized for π -phase shift at 25 keV with a period of 4.8 μm; the analyzer grating, G2 , had a period of 2.4 μm. The gratings were produced by Microworks (Karlsruhe, Germany) and to support the high-aspect ratio structures, both gratings included support structures: the G1 phase grating was fabricated with a broken-line structure and the G2 analyzer included a bridge-structure connecting adjacent lines. Experiments were performed with a sample-to-phase grating distance of 270 mm; a much shorter distance was desired but not possible due to mechanical interference between stacks of positioning stages in the test setup. The analyzer grating was then positioned at the first fractional Talbot distance (m = 1) 58 mm from the phase grating. Both phase and analyzer gratings were mounted on dual-tilt stages. Due to a smaller source size along the vertical plane, the X-ray radiation has its greatest coherence along the lab vertical; both gratings were positioned with the gratings aligned with the lab horizontal. The tomography sample rotation was about the lab vertical axis. The stepping scan motion along the lab vertical axis was performed with a piezoelectric-based positioner stage with a 200-μm range and sub-nanometer resolution. A 100-μm thick LuAG:Ce scintillator was imaged with a 10 × optical magnification lens and a Coolsnap HQ2 CCD with a 1040 × 1392 array of 6.45-μm-square pixels; the small effective pixel size of 0.645 μm was used to accommodate another experiment performed in the test setup. The exposure time for each interferogram was 600 ms, and 16 interferograms were measured across 4/3 period of the analyzer grating structure. Reference data were collected before and after the tomography data; 601 projections were acquired in the angle range of 0◦ to 180◦ . Also, background images (X-rays off) were collected to account for CCD dark-current and the offset signal.12 The tomography data were acquired over several hours. A comparison of the two reference data sets shows some drift as seen in the values of φ p (reference). Test regions in the air around the tomography sample were used to correct the
(9)
differential phase contrast = φp (sample) − φp (reference) , (10) dark-field =
ap (sample) /ap (reference) , a1p (sample) /a1p (reference)
where the polar coefficients are defined in Eq. (1).
(11) FIG. 1. A schematic side-view of the experimental setup for the two-grating system.
013704-3
Marathe et al.
φ p (reference) values as used in Eq. (10) and this correction was successful over about 40% of the field of view. IV. RESULTS AND DISCUSSION
To assess the new algorithm, we first examine 2D images of polystyrene spheres supported on a Kapton film and imaged at the m = 1 (58 mm) fractional Talbot distance. The 16grating steps cover 4/3-period of the analyzer structure. A discrete Fourier transform method is best used with 10-grating steps, a full period, while Levenberg-Marquardt and the new algorithm can use any or all of the grating steps. In addition, the gratings were fabricated with support structures that leak into the differential phase contrast (DPC) image. The new algorithm is tested with several different selections of grating steps to find an optimum pattern for image processing. The computational time for Eq. (6) for the sample data set was about 1 s in either Mathematica v9 or MATLAB 2012b running on a common workstation computer. This is about three orders of magnitude faster than fitting Eq. (1) with a Levenberg-Marquardt routine and about the same as for the discrete Fourier transform method. The differential phase contrast image of one 50-μm-diameter sphere, Fig. 2(a), is shown after baseline correction (the mean DPC of a 50 × 50 pixel region near the sphere) to remove the phase shift from the Kapton film support. The sensitivity of the phase contrast experiment for light atom samples such as polystyrene when imaged at 25 keV is remarkable. The line probe through the sphere, Fig. 2(b), shows the sphere edges and center. Two pix-
Rev. Sci. Instrum. 85, 013704 (2014)
els along the line probe are selected for further examination; the pixel at {row=104,column=70}, Fig. 2(c), represents a typical fit while the pixel at {96,70}, Fig. 2(f), has a large χν2 that is attributed to the support structure in the G1 , G2 gratings. The support structure is visually apparent in the differential phase image, Fig. 2(a), as well as in the transmission and dark-field images (not shown). To illustrate the robustness and flexibility of the new algorithm, the data from the 7th grating step have been eliminated from the processing. For comparison, the discrete Fourier transform method is applied to the 10-grating step data set, also with omission of the 7th grating step. The results in Fig. 3 show a poor fit to the data points, more noise in the DPC image, and much higher χν2 values. The results of several fitting methods and grating step selections are summarized in Table I. The Levenberg-Marquardt and new algorithm are nearly equivalent for use with a full period of interferometry data or with under or oversampling. The discrete Fourier transform method performs worse than the other methods when used with a partial data set. The new algorithm gives roughly equivalent results for grating step patterns of [1-10], [1-6,8-10], and [2-11] for a single pixel as listed in Table I and for DPC images (not shown). The [1-16] grating step pattern does yield a visually better image as will be discussed next. We are at an early stage for developing optimal methods for interferometry using gratings with support structures. In Fig. 4, we compare the new algorithm with 16-grating steps (4/3-period) to the discrete Fourier transform method with
FIG. 2. The new algorithm is used to analyze a cropped image from the study of polystyrene spheres supported on a Kapton film. The data from the 7th stepped grating position have been omitted from the processing. The image origin is in the upper left corner at {1,1} and the field of view is 200 × 200 pixels. (a) Differential phase contrast (DPC) image (Eq. (10)) for a 50-μm-diameter sphere and (d) corresponding plot of χν2 (sample) (from Eq. (1)). (b) The DPC line probe for column 70. (e) The χν2 (sample) line probe plot shows a wide range of values; two representative pixels (asterisks) are shown in (c) and (f); the poor fit in (f) is attributed to the grating support structures. For (c) and (f), the solid trace, •, is the reference interferometry data set and the dashed trace, ◦, is the interferometry data set with the sample in position.
013704-4
Marathe et al.
Rev. Sci. Instrum. 85, 013704 (2014)
FIG. 3. ((a)–(f)) A discrete Fourier transform method is used to analyze the same data set as shown in Fig. 2. Again, the data from the 7th grating step have been eliminated from the processing.
TABLE I. Comparison of three methods for pixel {104,70}. Grating steps 1-10 1-6, 8-10 1-6, 8-10 1-6, 8-10 1-16 2-11 a b
Algorithma
Transmission
DPC
Errorb (%)
Dark-field
χν2 (sample)
LM, FT, Basis LM Basis FT Basis Basis
0.9345 0.9331 0.9331 0.9233 0.9376 0.9353
−0.1743 −0.1641 −0.1641 −0.1457 −0.1862 −0.1765
... 5.88 5.88 16.4 6.83 1.24
1.0819 1.0904 1.0904 1.1530 1.0985 1.0612
0.63 0.75 0.75 5.42 1.34 0.17
LM: Levenberg-Marquardt, FT: discrete Fourier transform, and Basis: new algorithm in Sec II. Error with respect to Levenberg-Marquardt fit for grating steps 1-10.
FIG. 4. Comparison of algorithms with best possible sampling for this data set. ((a) and (b)) New algorithm with 16-grating steps (4/3-period). ((c) and (d)) Discrete Fourier transform with 10-grating steps (1-period). Image artifacts from the G1 , G2 grating support structures are diminished with the new algorithm and a larger grating step pattern.
013704-5
Marathe et al.
Rev. Sci. Instrum. 85, 013704 (2014)
FIG. 5. Slices from a tomographic volume. (a) A filtered back-projection reconstruction of a sinogram constructed from transmission images (Eq. (9)). (b) An ART reconstruction of the differential phase contrast images (Eq. (10)) followed by two-level segmentation to emphasize positive and negative regions of φ values. (c) A filtered back-projection reconstruction of the differential phase contrast images using the imaginary reconstruction kernel described by Pfeiffer and co-workers.18
10-grating steps (1 period). The additional data does yield an improved image. This suggests that a procedure exists for matching support structure design with a basis set, not necessarily sinusoidal, and corresponding grating step structure, not necessarily equally spaced. The new algorithm is tested with phase contrast tomography. The foraminifera sample is nominally calcium carbonate; the expected uniformity is seen in the X-ray transmission image in Fig. 5(a). Interestingly, the differential phase contrast image shows a significantly different structure, especially when the differences are highlighted in a two-level segmentation as shown in Fig. 5(b) with white regions for positive φ values and gray for negative φ values. The differential phase contrast sinogram was reconstructed with algebraic reconstruction tomography.17 Figure 5(c) shows an alternative reconstruction based on the imaginary kernel proposed by Pfeiffer and co-workers.18
knowledged. This material is based upon work supported by the National Science Foundation under the NSF EPSCoR Cooperative Agreement No. EPS-1003897 with additional support from the Louisiana Board of Regents.
V. CONCLUSIONS
We can expand cˆg as
A fast, robust, easily implemented, and flexible algorithm is presented. The algorithm offers a route to tuning data acquisition so as to reduce the artifacts from phase and analyzer grating support structures. ACKNOWLEDGMENTS
We thank Professor Gabor Herman, Joanna Klukowska, and Professor Stephen Shipman for enlightening conversations about reconstruction and phase analysis. Use of the Advanced Photon Source and Center for Nanoscale Materials, Office of Science User Facilities operated for the U.S. Department of Energy (DOE), Office of Science, by Argonne National Laboratory, was supported by the U.S. DOE under Contract No. DE-AC02-06CH11357. The support of National Science Foundation (NSF) CHE-0910937 is gratefully ac-
APPENDIX: A DETAILED DERIVATION
Let us return to the development of Eqs. (5) to (7). First, let us simplify the discussion by considering a single pixel. Then, Eq. (5) reduces to the scalar D for the deviation, a data ¯ˆ The data and fit vectors have length ¯ and a fit vector c. vector c, M, the number of grating steps. The deviation is defined as before: D=
M (cg − cˆg )2 .
(A1)
g=1
cˆg =
3
Bgμ aμ ,
(A2)
μ=1
where B is an M × 3 matrix of basis vectors spanning an M-dimensional vector space. Recall c¯ˆ is a fit vector in the M-dimensional vector space with three degrees of freedom, given by the coefficients a¯ and the basis vectors in B. We wish to minimize the distance, Eq. (A1), between the data vector c¯ ˆ¯ Using the calculus of variations, we difand the fit vector c. ferentiate Eq. (A1) with respect to the three coefficients aν to get ∂ cˆg ∂D = 2(cg − cˆg )1 , ∂aν ∂aν g=1 M
(A3)
013704-6
Marathe et al.
Rev. Sci. Instrum. 85, 013704 (2014)
where
The reduced chi-square statistic is the chi-square divided by the degrees of freedom:
∂ cˆg = Bgμ δμν , ∂aν μ=1 3
(A4) χν2 =
and where δ μν is the Dirac delta function. Substituting Eq. (A4) into Eq. (A3), we get an important intermediate result: ∂D = ∂aν
M
2(cg − cˆg )1 · Bgν .
(A5)
g=1
∂D Then, for ν = 1. . . 3, we set ∂a = 0 since we are seeking ν ¯ the closest alignment of the fit vector c¯ˆ with the data vector c. Then, Eq. (A5) becomes
0=
M
cg Bgν −
g=1
0=
M
cg Bgν
M
cˆg Bgν ,
(A6)
g=1
⎛ ⎞ M 3 ⎝ − Bgμ aμ ⎠ Bgν .
g=1
g=1
(A7)
μ=1
¯ via a sequence of matrix multiNow, we solve for aμ (as a) plication, transpose, and inverse operations: ⎛ ⎞ M 3 M ⎝ Bgμ aμ ⎠ Bgν = cg Bgν , g=1
M
μ=1
⎛ T ⎝ Bνg
(A8)
g=1
3
⎞ Bgμ aμ ⎠ =
M
T Bνg cg ,
(A9)
⎛ ⎞ 3 M M T T ⎝ Bνg Bgμ ⎠ aμ = Bνg cg ,
(A10)
g=1
μ=1
μ=1
g=1
g=1
g=1
3
(BT · B)νμ aμ = BT · c¯ ν ,
(A11)
μ=1
¯ ν = (BT · c) ¯ ν, [(BT · B) · a]
(A12)
¯ (BT · B) · a¯ = BT · c,
(A13)
¯ a¯ = (BT · B)−1 · BT · c.
(A14)
Recall Eqs. (6) and (7) for all pixels: ¯ a¯ = G · c,
(A15)
G = (BT · B)−1 · BT .
(A16)
M 1 (cg − cˆg )2 ν g=1 σg2
(A17)
(cg − cˆg )2 1 , (M − 3 − 1) g=1 cg
(A18)
M
=
where we have assumed that the standard deviations are given by the square root of the number of counts. 1 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). 2 A. Momose, A. Fujii, H. Kadowaki, and H. Jinnai, “Three-dimensional observation of polymer blend by X-ray phase tomography,” Macromolecules 38, 7197–7200 (2005). 3 A. Momose, “Recent advances in X-ray phase imaging,” Jpn. J. Appl. Phys., Part 1 44, 6355–6367 (2005). 4 F. Pfeiffer, C. Grunzweig, O. Bunk, G. Frei, E. Lehmann, and C. David, “Neutron phase imaging and tomography,” Phys. Rev. Lett. 96, 215505 (2006). 5 F. Pfeiffer, C. Kottler, O. Bunk, and C. David, “Hard X-ray phase tomography with low-brilliance sources,” Phys. Rev. Lett. 98, 108105 (2007). 6 P. P. Zhu, K. Zhang, Z. L. Wang, Y. J. Liu, X. S. Liu, Z. Y. Wu, S. A. McDonald, F. Marone, and M. Stampanoni, “Low-dose, simple, and fast grating-based X-ray phase-contrast imaging,” Proc. Natl. Acad. Sci. U.S.A. 107, 13576–13581 (2010). 7 S. W. Lee, Y. K. Jun, and O. Y. Kwon, “A neutron dark-field imaging experiment with a neutron grating interferometer at a thermal neutron beam line at HANARO,” J. Korean Phys. Soc. 58, 730–734 (2011). 8 J. Herzen, T. Donath, F. Beckmann, M. Ogurreck, C. David, J. Mohr, F. Pfeiffer, and A. Schreyer, “X-ray grating interferometer for materialsscience imaging at a low-coherent wiggler source,” Rev. Sci. Instrum. 82, 113711 (2011). 9 S. Rutishauser, L. Samoylova, J. Krzywinski, O. Bunk, J. Gruenert, H. Sinn, M. Cammarata, D. M. Fritz, and C. David, “Exploring the wavefront of hard X-ray free-electron laser radiation,” Nat. Commun. 3, 947 (2012). 10 I. Zanette, M. Bech, A. Rack, G. Le Duc, P. Tafforeau, C. David, J. Mohr, F. Pfeiffer, and T. Weitkamp, “Trimodal low-dose X-ray tomography,” Proc. Natl. Acad. Sci. U.S.A. 109, 10199–10204 (2012). 11 B. P. Abbott, R. Abbott et al., “The Laser Interferometer GravitationalWave Observatory,” Rep. Prog. Phys. 72, 076901 (2009). 12 S. M. Gruner, M. W. Tate, and E. F. Eikenberry, “Charge-coupled device area x-ray detectors,” Rev. Sci. Instrum. 73, 2815–2842 (2002). 13 T. Koehler, G. Martens, U. van Stevendaal, and E. Roessl, “Non-scatter contributions to the dark-field signal in differential phase contrast imaging,” AIP Conf. Proc. 1466, 205–210 (2012). 14 R. M. Baur, “Development and application of a grating interferometer at the Cornell High Energy Synchrotron Source,” Ph.D. dissertation (Cornell University, 2013). 15 J. Mathews and R. L. Walker, Mathematical Methods of Physics, (W. A. Benjamin, 1965). 16 S. Marathe, X. Xiao, M. J. Wojcik, R. Divan, L. G. Butler, K. Ham, K. Fezzaa, M. Erdmann, H. H. Wen, W.-K. Lee, A. T. Macrander, F. De Carlo, D. C. Mancini, and L. Assoufid, “Development of grating-based x-ray Talbot interferometry at the Advanced Photon Source,” AIP Conf. Proc. 1466, 249–254 (2012). 17 J. Klukowska, R. Davidi, and G. T. Herman, “SNARK09—A software package for reconstruction of 2D images from 1D projections,” Comput. Methods Programs Biomed. 110, 424–440 (2013). 18 F. Pfeiffer, O. Bunk, C. Kottler, and C. David, “Tomographic reconstruction of three-dimensional objects from hard X-ray differential phase contrast projection images,” Nucl. Instrum. Methods Phys. Res., Sect. A 580, 925– 928 (2007).
Review of Scientific Instruments is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/rsio/rsicr.jsp
Improved algorithm for processing grating-based phase contrast interferometry image sets Shashidhara Marathe,1,a) Lahsen Assoufid,1,b) Xianghui Xiao,1,c) Kyungmin Ham,2,d) Warren W. Johnson,3,e) and Leslie G. Butler4,f) 1
Advanced Photon Source, Argonne National Laboratory, 9700 S. Cass Ave., Argonne, Illinois 60439, USA CAMD, LSU, Baton Rouge, Louisiana 70806, USA 3 Department of Physics and Astronomy, LSU, Baton Rouge, Louisiana 70803, USA 4 Department of Chemistry, LSU, Baton Rouge, Louisiana 70803, USA 2
(Received 10 December 2013; accepted 19 December 2013; published online 17 January 2014) Grating-based X-ray and neutron interferometry tomography using phase-stepping methods generates large data sets. An improved algorithm is presented for solving for the parameters to calculate transmissions, differential phase contrast, and dark-field images. The method takes advantage of the vectorization inherent in high-level languages such as Mathematica and MATLAB and can solve a 16 × 1k × 1k data set in less than a second. In addition, the algorithm can function with partial data sets. This is demonstrated with processing of a 16-step grating data set with partial use of the original data chosen without any restriction. Also, we have calculated the reduced chi-square for the fit and notice the effect of grating support structural elements upon the differential phase contrast image and have explored expanded basis set representations to mitigate the impact. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4861199] I. INTRODUCTION
Grating-based interferometry has been implemented at many X-ray and neutron imaging beamlines worldwide to perform phase contrast imaging and tomography with applications in materials science, biomedical science, and beam wavefront and optics characterization.1–10 These measurement tools generate large data sets that are processed using specialized computer codes. A variety of algorithms have been developed for this purpose. A discrete Fourier transform algorithm1 is quite popular, is fast, yet can suffer from spectral leakage when the grating steps are not uniformly spaced within a period. A Levenberg-Marquardt algorithm can solve non-uniformly spaced data, yet is slow. In this paper, we describe the application of an algorithm from the gravity wave toolbox that provides fast and robust processing. The toolbox was originally developed to rapidly process the large data volume that is typically generated by the interferometers used to search for gravity waves caused by events such as the collision of two black holes.11 The X-ray interferometry data were collected with a prototype experimental setup at the Advanced Photon Source. A grating-based phase-stepping interferometry data set consists of a set of pixelated images acquired as a function of a grating stepping. The triangular interferometric signal13, 14 is approximated as a sinusoid for which we wish to determine the amplitude and phase; the sinusoid period is usually known from the grating motion. The speed of this algorithm is based on pre-computation of the sinusoidal functions and a) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected] d) Electronic mail: [email protected] e) Electronic mail: [email protected] f) Electronic mail: [email protected]
0034-6748/2014/85(1)/013704/6/$30.00
vectorization of the pixel data. The amplitude and phase are then calculated with standard matrix manipulations. Two samples are studied: The first sample is made of polystyrene spheres supported by a Kapton film; reconstructions are performed with complete and partial data sets. The second sample is a calcium-carbonate-rich shell of a foraminifera, a one-cell protist, which shows a highly structured, millimeter-sized prolate object in a tomographic study. The reconstructed absorption and differential phase images demonstrate the applicability of this new algorithm. II. THEORY
The interferometry data are a set of X-ray counts cgp , where g = 1, . . . , M is the index that identifies the exposure, with one exposure at each grating displacement xg ; and p = 1, . . . , N is the index that identifies the pth pixel. Typically M ∼ 7 to 16 and N ∼ millions. We have adapted a fitting algorithm from mathematical physics (section 14-6 in Mathews and Walker15 ) that transforms the fitting operation into a simple matrix problem. We fit each set of exposures pixel by pixel, cgp , to each pixel’s expected dependence on grating position xg : 2π cˆgp = a1p + ap sin xg + φp (1) pgrat 2π xg ap cos(φp ) ≡ [1] a1p + sin pgrat 2π xg ap sin(φp ) (2) + cos pgrat ≡
3
Bgμ aμp
(3)
μ=1
85, 013704-1
© 2014 AIP Publishing LLC
013704-2
Marathe et al.
Rev. Sci. Instrum. 85, 013704 (2014)
a2p ≡ ap cos(φp ), a3p ≡ ap sin(φp ),
(4)
where ap and φ p are the amplitude and phase of the interference term for the pth pixel, and pgrat is the period of the translated grating. The interference term is expressed in both polar form (Eq. (1)) and Cartesian form (Eq. (4)). The latter is used to define the M × 3 matrix B (Eq. (3)) that represents the three fitting functions shown in brackets in Eq. (2): constant, sine, and cosine. The aμp is one element in the coefficient matrix a, which is a 3 × N matrix of amplitudes, or weights, of the three fitting functions for the N pixels. Later, we will reshape a into a 3D matrix with dimensionality [rows, columns, 3]. The best fit can be chosen to be the set of coefficients aμp that minimize the deviation-squared, Dp , for each pixel, defined by M (cgp − cˆgp )2 . Dp ≡
(5)
g=1
To minimize Dp , we need only calculate the derivatives of each deviation-squared with respect to each component of a, set all deviations to zero, and solve the resulting matrix equations. The closed form expression for the coefficient matrix a is found to be a = G · c,
(6)
G = (BT · B)−1 · BT ,
(7)
where
and where the superscript T indicates the matrix transpose. So the optimization problem is reduced to multiplying a 3 × M fixed matrix G into the M × N data matrix c, to find the 3 × N coefficient matrix a. With an efficient matrix manipulation program like MATLAB, the multiplication is extremely fast, of the order of 1 s for 1k × 1k images times 16 grating steps. By inspection, we recover the polar coefficients from the Cartesian in the usual way: 2 2 + a3p , φp = tan−1 (a3p /a2p ). (8) ap = a2p We have extended this analysis to use knowledge of the ex2 pected uncertainty or variance of each datum σgp at the price of increased computing time. The model equation, Eq. (4), can be extended to include anharmonic terms to describe the grating effects more accurately and to include more sophisticated grating motion trajectories. The Appendix has more details about the derivation of Eq. (7). In the stepped-grating interferometry experiment, two data sets are measured: one without the sample in the beam (reference data set) and the other with the sample in the beam (sample data set): transmission =
a1p (sample) , a1p (reference)
III. EXPERIMENTAL
A two-grating interferometer was assembled at the Advanced Photon Source 2-BM-B beamline (Fig. 1).16 The stepped-grating system was installed 23 m downstream of the double multilayer monochromator, which provides a monochromatic beam at 25 keV with 1% bandwidth. The phase grating, G1 , was optimized for π -phase shift at 25 keV with a period of 4.8 μm; the analyzer grating, G2 , had a period of 2.4 μm. The gratings were produced by Microworks (Karlsruhe, Germany) and to support the high-aspect ratio structures, both gratings included support structures: the G1 phase grating was fabricated with a broken-line structure and the G2 analyzer included a bridge-structure connecting adjacent lines. Experiments were performed with a sample-to-phase grating distance of 270 mm; a much shorter distance was desired but not possible due to mechanical interference between stacks of positioning stages in the test setup. The analyzer grating was then positioned at the first fractional Talbot distance (m = 1) 58 mm from the phase grating. Both phase and analyzer gratings were mounted on dual-tilt stages. Due to a smaller source size along the vertical plane, the X-ray radiation has its greatest coherence along the lab vertical; both gratings were positioned with the gratings aligned with the lab horizontal. The tomography sample rotation was about the lab vertical axis. The stepping scan motion along the lab vertical axis was performed with a piezoelectric-based positioner stage with a 200-μm range and sub-nanometer resolution. A 100-μm thick LuAG:Ce scintillator was imaged with a 10 × optical magnification lens and a Coolsnap HQ2 CCD with a 1040 × 1392 array of 6.45-μm-square pixels; the small effective pixel size of 0.645 μm was used to accommodate another experiment performed in the test setup. The exposure time for each interferogram was 600 ms, and 16 interferograms were measured across 4/3 period of the analyzer grating structure. Reference data were collected before and after the tomography data; 601 projections were acquired in the angle range of 0◦ to 180◦ . Also, background images (X-rays off) were collected to account for CCD dark-current and the offset signal.12 The tomography data were acquired over several hours. A comparison of the two reference data sets shows some drift as seen in the values of φ p (reference). Test regions in the air around the tomography sample were used to correct the
(9)
differential phase contrast = φp (sample) − φp (reference) , (10) dark-field =
ap (sample) /ap (reference) , a1p (sample) /a1p (reference)
where the polar coefficients are defined in Eq. (1).
(11) FIG. 1. A schematic side-view of the experimental setup for the two-grating system.
013704-3
Marathe et al.
φ p (reference) values as used in Eq. (10) and this correction was successful over about 40% of the field of view. IV. RESULTS AND DISCUSSION
To assess the new algorithm, we first examine 2D images of polystyrene spheres supported on a Kapton film and imaged at the m = 1 (58 mm) fractional Talbot distance. The 16grating steps cover 4/3-period of the analyzer structure. A discrete Fourier transform method is best used with 10-grating steps, a full period, while Levenberg-Marquardt and the new algorithm can use any or all of the grating steps. In addition, the gratings were fabricated with support structures that leak into the differential phase contrast (DPC) image. The new algorithm is tested with several different selections of grating steps to find an optimum pattern for image processing. The computational time for Eq. (6) for the sample data set was about 1 s in either Mathematica v9 or MATLAB 2012b running on a common workstation computer. This is about three orders of magnitude faster than fitting Eq. (1) with a Levenberg-Marquardt routine and about the same as for the discrete Fourier transform method. The differential phase contrast image of one 50-μm-diameter sphere, Fig. 2(a), is shown after baseline correction (the mean DPC of a 50 × 50 pixel region near the sphere) to remove the phase shift from the Kapton film support. The sensitivity of the phase contrast experiment for light atom samples such as polystyrene when imaged at 25 keV is remarkable. The line probe through the sphere, Fig. 2(b), shows the sphere edges and center. Two pix-
Rev. Sci. Instrum. 85, 013704 (2014)
els along the line probe are selected for further examination; the pixel at {row=104,column=70}, Fig. 2(c), represents a typical fit while the pixel at {96,70}, Fig. 2(f), has a large χν2 that is attributed to the support structure in the G1 , G2 gratings. The support structure is visually apparent in the differential phase image, Fig. 2(a), as well as in the transmission and dark-field images (not shown). To illustrate the robustness and flexibility of the new algorithm, the data from the 7th grating step have been eliminated from the processing. For comparison, the discrete Fourier transform method is applied to the 10-grating step data set, also with omission of the 7th grating step. The results in Fig. 3 show a poor fit to the data points, more noise in the DPC image, and much higher χν2 values. The results of several fitting methods and grating step selections are summarized in Table I. The Levenberg-Marquardt and new algorithm are nearly equivalent for use with a full period of interferometry data or with under or oversampling. The discrete Fourier transform method performs worse than the other methods when used with a partial data set. The new algorithm gives roughly equivalent results for grating step patterns of [1-10], [1-6,8-10], and [2-11] for a single pixel as listed in Table I and for DPC images (not shown). The [1-16] grating step pattern does yield a visually better image as will be discussed next. We are at an early stage for developing optimal methods for interferometry using gratings with support structures. In Fig. 4, we compare the new algorithm with 16-grating steps (4/3-period) to the discrete Fourier transform method with
FIG. 2. The new algorithm is used to analyze a cropped image from the study of polystyrene spheres supported on a Kapton film. The data from the 7th stepped grating position have been omitted from the processing. The image origin is in the upper left corner at {1,1} and the field of view is 200 × 200 pixels. (a) Differential phase contrast (DPC) image (Eq. (10)) for a 50-μm-diameter sphere and (d) corresponding plot of χν2 (sample) (from Eq. (1)). (b) The DPC line probe for column 70. (e) The χν2 (sample) line probe plot shows a wide range of values; two representative pixels (asterisks) are shown in (c) and (f); the poor fit in (f) is attributed to the grating support structures. For (c) and (f), the solid trace, •, is the reference interferometry data set and the dashed trace, ◦, is the interferometry data set with the sample in position.
013704-4
Marathe et al.
Rev. Sci. Instrum. 85, 013704 (2014)
FIG. 3. ((a)–(f)) A discrete Fourier transform method is used to analyze the same data set as shown in Fig. 2. Again, the data from the 7th grating step have been eliminated from the processing.
TABLE I. Comparison of three methods for pixel {104,70}. Grating steps 1-10 1-6, 8-10 1-6, 8-10 1-6, 8-10 1-16 2-11 a b
Algorithma
Transmission
DPC
Errorb (%)
Dark-field
χν2 (sample)
LM, FT, Basis LM Basis FT Basis Basis
0.9345 0.9331 0.9331 0.9233 0.9376 0.9353
−0.1743 −0.1641 −0.1641 −0.1457 −0.1862 −0.1765
... 5.88 5.88 16.4 6.83 1.24
1.0819 1.0904 1.0904 1.1530 1.0985 1.0612
0.63 0.75 0.75 5.42 1.34 0.17
LM: Levenberg-Marquardt, FT: discrete Fourier transform, and Basis: new algorithm in Sec II. Error with respect to Levenberg-Marquardt fit for grating steps 1-10.
FIG. 4. Comparison of algorithms with best possible sampling for this data set. ((a) and (b)) New algorithm with 16-grating steps (4/3-period). ((c) and (d)) Discrete Fourier transform with 10-grating steps (1-period). Image artifacts from the G1 , G2 grating support structures are diminished with the new algorithm and a larger grating step pattern.
013704-5
Marathe et al.
Rev. Sci. Instrum. 85, 013704 (2014)
FIG. 5. Slices from a tomographic volume. (a) A filtered back-projection reconstruction of a sinogram constructed from transmission images (Eq. (9)). (b) An ART reconstruction of the differential phase contrast images (Eq. (10)) followed by two-level segmentation to emphasize positive and negative regions of φ values. (c) A filtered back-projection reconstruction of the differential phase contrast images using the imaginary reconstruction kernel described by Pfeiffer and co-workers.18
10-grating steps (1 period). The additional data does yield an improved image. This suggests that a procedure exists for matching support structure design with a basis set, not necessarily sinusoidal, and corresponding grating step structure, not necessarily equally spaced. The new algorithm is tested with phase contrast tomography. The foraminifera sample is nominally calcium carbonate; the expected uniformity is seen in the X-ray transmission image in Fig. 5(a). Interestingly, the differential phase contrast image shows a significantly different structure, especially when the differences are highlighted in a two-level segmentation as shown in Fig. 5(b) with white regions for positive φ values and gray for negative φ values. The differential phase contrast sinogram was reconstructed with algebraic reconstruction tomography.17 Figure 5(c) shows an alternative reconstruction based on the imaginary kernel proposed by Pfeiffer and co-workers.18
knowledged. This material is based upon work supported by the National Science Foundation under the NSF EPSCoR Cooperative Agreement No. EPS-1003897 with additional support from the Louisiana Board of Regents.
V. CONCLUSIONS
We can expand cˆg as
A fast, robust, easily implemented, and flexible algorithm is presented. The algorithm offers a route to tuning data acquisition so as to reduce the artifacts from phase and analyzer grating support structures. ACKNOWLEDGMENTS
We thank Professor Gabor Herman, Joanna Klukowska, and Professor Stephen Shipman for enlightening conversations about reconstruction and phase analysis. Use of the Advanced Photon Source and Center for Nanoscale Materials, Office of Science User Facilities operated for the U.S. Department of Energy (DOE), Office of Science, by Argonne National Laboratory, was supported by the U.S. DOE under Contract No. DE-AC02-06CH11357. The support of National Science Foundation (NSF) CHE-0910937 is gratefully ac-
APPENDIX: A DETAILED DERIVATION
Let us return to the development of Eqs. (5) to (7). First, let us simplify the discussion by considering a single pixel. Then, Eq. (5) reduces to the scalar D for the deviation, a data ¯ˆ The data and fit vectors have length ¯ and a fit vector c. vector c, M, the number of grating steps. The deviation is defined as before: D=
M (cg − cˆg )2 .
(A1)
g=1
cˆg =
3
Bgμ aμ ,
(A2)
μ=1
where B is an M × 3 matrix of basis vectors spanning an M-dimensional vector space. Recall c¯ˆ is a fit vector in the M-dimensional vector space with three degrees of freedom, given by the coefficients a¯ and the basis vectors in B. We wish to minimize the distance, Eq. (A1), between the data vector c¯ ˆ¯ Using the calculus of variations, we difand the fit vector c. ferentiate Eq. (A1) with respect to the three coefficients aν to get ∂ cˆg ∂D = 2(cg − cˆg )1 , ∂aν ∂aν g=1 M
(A3)
013704-6
Marathe et al.
Rev. Sci. Instrum. 85, 013704 (2014)
where
The reduced chi-square statistic is the chi-square divided by the degrees of freedom:
∂ cˆg = Bgμ δμν , ∂aν μ=1 3
(A4) χν2 =
and where δ μν is the Dirac delta function. Substituting Eq. (A4) into Eq. (A3), we get an important intermediate result: ∂D = ∂aν
M
2(cg − cˆg )1 · Bgν .
(A5)
g=1
∂D Then, for ν = 1. . . 3, we set ∂a = 0 since we are seeking ν ¯ the closest alignment of the fit vector c¯ˆ with the data vector c. Then, Eq. (A5) becomes
0=
M
cg Bgν −
g=1
0=
M
cg Bgν
M
cˆg Bgν ,
(A6)
g=1
⎛ ⎞ M 3 ⎝ − Bgμ aμ ⎠ Bgν .
g=1
g=1
(A7)
μ=1
¯ via a sequence of matrix multiNow, we solve for aμ (as a) plication, transpose, and inverse operations: ⎛ ⎞ M 3 M ⎝ Bgμ aμ ⎠ Bgν = cg Bgν , g=1
M
μ=1
⎛ T ⎝ Bνg
(A8)
g=1
3
⎞ Bgμ aμ ⎠ =
M
T Bνg cg ,
(A9)
⎛ ⎞ 3 M M T T ⎝ Bνg Bgμ ⎠ aμ = Bνg cg ,
(A10)
g=1
μ=1
μ=1
g=1
g=1
g=1
3
(BT · B)νμ aμ = BT · c¯ ν ,
(A11)
μ=1
¯ ν = (BT · c) ¯ ν, [(BT · B) · a]
(A12)
¯ (BT · B) · a¯ = BT · c,
(A13)
¯ a¯ = (BT · B)−1 · BT · c.
(A14)
Recall Eqs. (6) and (7) for all pixels: ¯ a¯ = G · c,
(A15)
G = (BT · B)−1 · BT .
(A16)
M 1 (cg − cˆg )2 ν g=1 σg2
(A17)
(cg − cˆg )2 1 , (M − 3 − 1) g=1 cg
(A18)
M
=
where we have assumed that the standard deviations are given by the square root of the number of counts. 1 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). 2 A. Momose, A. Fujii, H. Kadowaki, and H. Jinnai, “Three-dimensional observation of polymer blend by X-ray phase tomography,” Macromolecules 38, 7197–7200 (2005). 3 A. Momose, “Recent advances in X-ray phase imaging,” Jpn. J. Appl. Phys., Part 1 44, 6355–6367 (2005). 4 F. Pfeiffer, C. Grunzweig, O. Bunk, G. Frei, E. Lehmann, and C. David, “Neutron phase imaging and tomography,” Phys. Rev. Lett. 96, 215505 (2006). 5 F. Pfeiffer, C. Kottler, O. Bunk, and C. David, “Hard X-ray phase tomography with low-brilliance sources,” Phys. Rev. Lett. 98, 108105 (2007). 6 P. P. Zhu, K. Zhang, Z. L. Wang, Y. J. Liu, X. S. Liu, Z. Y. Wu, S. A. McDonald, F. Marone, and M. Stampanoni, “Low-dose, simple, and fast grating-based X-ray phase-contrast imaging,” Proc. Natl. Acad. Sci. U.S.A. 107, 13576–13581 (2010). 7 S. W. Lee, Y. K. Jun, and O. Y. Kwon, “A neutron dark-field imaging experiment with a neutron grating interferometer at a thermal neutron beam line at HANARO,” J. Korean Phys. Soc. 58, 730–734 (2011). 8 J. Herzen, T. Donath, F. Beckmann, M. Ogurreck, C. David, J. Mohr, F. Pfeiffer, and A. Schreyer, “X-ray grating interferometer for materialsscience imaging at a low-coherent wiggler source,” Rev. Sci. Instrum. 82, 113711 (2011). 9 S. Rutishauser, L. Samoylova, J. Krzywinski, O. Bunk, J. Gruenert, H. Sinn, M. Cammarata, D. M. Fritz, and C. David, “Exploring the wavefront of hard X-ray free-electron laser radiation,” Nat. Commun. 3, 947 (2012). 10 I. Zanette, M. Bech, A. Rack, G. Le Duc, P. Tafforeau, C. David, J. Mohr, F. Pfeiffer, and T. Weitkamp, “Trimodal low-dose X-ray tomography,” Proc. Natl. Acad. Sci. U.S.A. 109, 10199–10204 (2012). 11 B. P. Abbott, R. Abbott et al., “The Laser Interferometer GravitationalWave Observatory,” Rep. Prog. Phys. 72, 076901 (2009). 12 S. M. Gruner, M. W. Tate, and E. F. Eikenberry, “Charge-coupled device area x-ray detectors,” Rev. Sci. Instrum. 73, 2815–2842 (2002). 13 T. Koehler, G. Martens, U. van Stevendaal, and E. Roessl, “Non-scatter contributions to the dark-field signal in differential phase contrast imaging,” AIP Conf. Proc. 1466, 205–210 (2012). 14 R. M. Baur, “Development and application of a grating interferometer at the Cornell High Energy Synchrotron Source,” Ph.D. dissertation (Cornell University, 2013). 15 J. Mathews and R. L. Walker, Mathematical Methods of Physics, (W. A. Benjamin, 1965). 16 S. Marathe, X. Xiao, M. J. Wojcik, R. Divan, L. G. Butler, K. Ham, K. Fezzaa, M. Erdmann, H. H. Wen, W.-K. Lee, A. T. Macrander, F. De Carlo, D. C. Mancini, and L. Assoufid, “Development of grating-based x-ray Talbot interferometry at the Advanced Photon Source,” AIP Conf. Proc. 1466, 249–254 (2012). 17 J. Klukowska, R. Davidi, and G. T. Herman, “SNARK09—A software package for reconstruction of 2D images from 1D projections,” Comput. Methods Programs Biomed. 110, 424–440 (2013). 18 F. Pfeiffer, O. Bunk, C. Kottler, and C. David, “Tomographic reconstruction of three-dimensional objects from hard X-ray differential phase contrast projection images,” Nucl. Instrum. Methods Phys. Res., Sect. A 580, 925– 928 (2007).
Review of Scientific Instruments is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/rsio/rsicr.jsp
