Open Access
Add to BibSonomyAbstract
Numerical wave-optical simulations of X-ray differential phase-contrast imaging using grating interferometry require the oversampling of gratings and object structures in the range of few micrometers. Consequently, fields of view of few millimeters already use large amounts of a computer’s main memory to store the propagating wave front, limiting the scope of the investigations to only small-scale problems. In this study, we apply an approximation to the Fresnel-Kirchhoff diffraction theory to overcome these restrictions by dividing the two-dimensional wave front up into 1D lines, which are processed separately. The approach enables simulations with samples of clinically relevant dimensions by significantly reducing the memory footprint and the execution time and, thus, allows the qualitative comparison of different setup configurations. We analyze advantages as well as limitations and present the simulation of a virtual mammography phantom of several centimeters of size.
© 2014 Optical Society of America
1. Introduction
X-ray differential phase-contrast (DPC) imaging is a new modality that provides information about a sample not obtainable by conventional absorption-based X-ray imaging. DPC analyzes the phase-shift that a sample imposes on a wave front and enhances contrast in weakly attenuating materials such as soft tissue [1–5]. Therefore, it has great potential e.g. in medical diagnostics for the detection of breast cancer [6–8].
The DPC signal, together with the absorption and the dark-field signal - a third contrast, which yields information about the scattering by a sample and already showed promising results e.g. in lung imaging [9–11] and mammography [12] -, can be obtained within a single scan using a grating interferometer (GI, Fig. 1) [2]. The setup comprises a highly coherent X-ray source, a detector and two gratings (a phase grating and an absorption grating). The phase grating (G1) causes a phase shift in the incoming wave front, which reproduces in form of an interference pattern at certain distances downstream the propagation direction (z-axis). As the periods of gratings and pattern are in the order of a few micrometers, it is not possible to directly resolve the interference pattern with clinically used detectors. Therefore, the intensity variation is transformed into a stepping curve by moving an analyzer grating (G2) perpendicular to the beam direction and the grating lines (i.e. along the x-direction) before the signals are extracted e.g. by Fourier analysis. By placing a second absorption grating near the source it is possible to relax the requirements for high coherence, which is necessary for the interference effect, and to use this imaging modality at conventional laboratory X-ray sources (Talbot-Lau interferometer [9]).
Fig. 1 The setup for X-ray phase-contrast imaging as assumed in the presented simulations consists of an X-ray source (represented by a plane wave front), a detector and a Talbot interferometer: The phase grating (G1) imposes a phase shift on the incoming wave front, which reproduces as intensity variation at certain distances (Talbot distances dT) downstream of the beam path. As the interference pattern has a period smaller than the pixel size of commonly used detectors an analyzer grating (G2) is stepped over it to convert the pattern into an intensity curve, from which the signals can be extracted.
In order to investigate DPC and dark-field contrast and create a better understanding of the processes leading to the signal formation, simulations are an inexpensive and practicable tool as compared to changes in real table-top setups. In particular, wave-optical simulations based on the Fresnel approximation to the scalar Fresnel-Kirchhoff diffraction theory have proven to provide valuable information, as they are able to reproduce the involved physical effects like interference of waves and subsequent Fresnel diffraction phenomena [2, 3, 9, 13]. For the simulation of a DPC image the discrete wave front can be modeled as a two-dimensional array of complex values, which holds the information about amplitude and phase at the according position. In order to avoid aliasing caused by an insufficiently high sampling rate a wave-front spacing of tens to hundreds of nanometers has proven to be suitable. These small dimensions result from the micron-sized periods of the gratings and the need to resolve small sample features. Consequently, fields of view (FOV) of few millimeters already lead to an enormous usage of a computer’s main memory (random access memory, RAM), which makes qualitative comparisons between different setup configurations based on large samples with dimensions of several centimeters not feasible.
In this study we try to overcome restrictions caused by the array size and present an approach that reduces the memory footprint by dividing up the two-dimensional wave front into lines which are processed separately and in parallel. It is thereby possible to execute wave-optical simulations of large FOVs while still mimicking the propagation effects and the interaction with matter of sample and gratings. Furthermore, we evaluate the results obtained with this approach and compare it to standard 2D simulations in order to assess its performance.
2. Materials and methods
We developed a wave-optical simulation framework for the numerical quantitative investigation of X-ray phase-contrast and dark-field imaging. A detailed description can be found in [14]. The object-oriented code is written in Python and partly in C++ for increased speed. It allows the adjustment of all relevant parameters of the components of a DPC-GI setup including source spectra, gratings, detectors and samples.
The wave front is propagated according to the Fresnel approximation to the scalar diffraction theory [15]
The theory describes the wave propagation as a three-dimensional process by a complex scalar x–y wave phasor moving along the z-axis, with U(x, y, 0) being the wave field at the source position, λ the wave length of the monochromatic radiation and k the absolute value of the wave vector .The phase shift Φ imposed by the sample on the wave front leads to a refraction angle
which has components both in x- and y-direction. The Talbot interferometer - as shown in Fig. 1 - is characterized by two aligned gratings with the grating bars oriented in y-direction. As the analysis of the signal is based on measuring the relative displacement of the self-images of the one-dimensionally structured phase grating G1, the GI is only sensitive to phase shifts in x-direction: Consequently, shifts in y-direction cannot be resolved by the setup and, therefore, do not contribute to the measurement signal.The presented simulation approach makes use of the directional sensitivity of the GI by splitting up the 2D array of the wave front into one-dimensional lines in x-direction, which are propagated independently through the virtual setup and recombined at the detector. The propagation formula (1) (for a plane wave) is then simplified to a two-dimensional problem at the constant position y0 of the corresponding array line:
This reduces the memory needed for the wave front by a factor equal to the number of lines in y-direction. The simulation can be executed in parallel by creating different processes that take care of separate lines and distributing them to the available CPUs. The parallelization will decrease the run-time at the cost of RAM usage, which increases linearly with the number of processes in the current implementation.The simulation results presented here were obtained using plane-wave monochromatic illumination at an energy of 45.7 keV and gratings with periods of 5 μm and duty cycle of 0.5 (i.e. the rectangular bars cover 50% of the grating period). The samples were simulated in projection approximation [16] and assumed 1 cm in front of G1. The phase grating was chosen to be π/2 phase-shifting and comprised 8 μm Ni lines on a 500 μm Si wafer. The analyzer grating, which was stepped over ten steps, was positioned at a distance of 46.07 cm to G1, which equals the first fractional Talbot distance, and was simulated with 150 μm gold lines on a 500 μm Si wafer. At the detector the interaction of the wave front with a 450 μm thick silicon scintillation layer was assumed and the resulting intensities (calculated from the squared absolute wave field) are presented at the resolution of the wave front (0.25 μm, Fig. 2) or were binned into square detector pixels of 25 μm (Fig. 3) or 100 μm (Fig. 4, wave-front resolution 0.17 μm), respectively. The DPC signal was calculated using Fourier analysis. The wave-front spacing was chosen based on an oversampling factor of the grating pitch, i.e. how many wave-front pixels represent the width of one grating period. The factors were 20 - as validated by Malecki et al. [14] - or 30, respectively.
Fig. 2 Simulation of an aluminum oxide sphere. Propagation of the X-rays leads to Fresnel fringes at the projected edges of the sphere (diameter 800 μm), which become more pronounced with increasing distance from the sample (a). At certain distances behind the object (here, first fractional Talbot distance of a π/2 phase grating at 45.7 keV) an interference pattern is formed. Magnifications b to e show the effects at the height of the sphere center. Simulations including the grating interferometer with a one-dimensional (d) or a two-dimensional (c) wave front, respectively, deliver almost identical results including the Fresnel fringes visible without the gratings (b). Subtraction of c and d proves that the lines of the interference pattern are shifted in the same way (e, absolute difference). The 1D approach cannot reproduce the propagation effect between wave front lines and therefore does not present the Fresnel fringes alternating in y-direction (f,i), in contrast to the conventional simulation (g). The subtraction image, however, demonstrates that the pattern lines of both approaches coincide well (h).
Fig. 3 Differential phase image of an aluminum oxide sphere (DPC signal given in rad). The sphere (diameter 800 μm) leads to a phase shift in two opposite directions. The results of both the conventional simulation with a 2D wave front (a) and the 1D approach (b) look very similar. The subtraction image (c, absolute difference) and line plot through the sphere center (e) show that the absolute error is close to zero at the projected edges that are parallel to the grating lines (y-direction). It becomes larger for increasing influence of the y-gradient on the phase shift (d). Towards the top and the bottom of the sphere the error vanishes because the phase shift is close to zero and the fringes cancel out.
Fig. 4 Simulation of a mammography phantom. An anthropomorphic breast software phantom by Pokrajac et al. [17] (a) was used to demonstrate the advantage of the 1D approach. The three-dimensional model mimics 450 ml breast compressed to 5cm thickness in medio-lateral position consisting of adipose tissue (light gray), fibrogladular tissue (white), Cooper’s ligaments (dark gray), skin and a medium. The representation had a voxel size of 200 μm and covered an area of 20 × 7 cm2 in the sagittal plane. It was simulated in projection approximation with a monochromatic plane wave at 45.7 keV in water. The differential phase image (b, DPC signal given in rad) is sampled at 100 μm detector pixel size. The memory footprint of the running simulation of around 500 MByte was very moderate compared to the theoretical value of the full 2D wave front (approx. factor 16300, overhead for 2D case not included).
To test the approach a digital mammography phantom [17], mimicking compressed female breast of a volume of approximately 450 ml, was simulated in the same setup configuration. Four different tissue types are contained: skin, adipose tissue, dense fibroglandular tissue and Cooper’s ligaments. Material compositions and densities were taken from [18, 19]. In order to avoid artifacts from the surface structure of the skin the first and last voxel layer in z-direction was removed. 41 voxels were added in x-direction to simulate the background. With a resolution of 350 × 1009 × 248 pixels (x-, y-, z-axis) and a voxel side length of 200 μm the wave front covering the whole sample had dimensions of 1212000 lines with 420000 pixels. The phantom was upsampled to the resolution of the wave front by two-dimensional linear interpolation.
3. Results
We compared DPC-GI simulations with the 1D wave-front approach to the propagation of a two-dimensional wave front. Fig. 2 shows the resulting intensities as they would be detected by an ideal detector with a quantum efficiency of 100 percent. The effect of an aluminum oxide sphere (diameter 800 μm) 47 cm behind the sample is visible through the characteristic Fresnel fringes at the projected edges of the sphere (Fig. 2(a,b,f)). If a grating interferometer is introduced, the phase shift caused by G1 induces an interference pattern at (in this case) the first fractional Talbot distance and the lines get shifted in x-direction by the sample (Fig. 2(c,d,g,i)). At the height of the sphere center both images show the same shift and also the Fresnel fringes; the subtraction image demonstrates that the difference of both cases is very small (Fig. 2(b–e)). The Fresnel fringes along y disappear for the 1D case, however, with increasing influence of the sample’s y-gradient on the refraction angle (Fig. 2(f–i)). Propagation effects in y-direction are not reproduced by the 1D approach, because the single wave-front lines are processed separately and independently from each other.
A differential phase-contrast image (Fig. 3) can be obtained by the stepping of G2 and the analysis of the stepping curve [2]. If the simulation of the same sphere is conducted assuming a detector with a pixel size of 25 μm both DPC images look very similar. At the y-position of the sphere center the error is very small as the phase shift is solely along the x-axis. Discrepancies can only be observed on the projected edges and become more pronounced at the upper and lower region of the sphere, where the influence of the gradient in y-direction increases. At the x-position of the sphere center the phase shift is zero and both simulations deliver the same result.
The memory usage is highly reduced by the 1D processing, which enables the simulations of large FOVs and samples (Fig. 4). In the example of the breast phantom, the calculated RAM usage of a 2D wave front with 128-bit complex values would be 8.14 TByte. With reduced accuracy (64-bit complex values) and 10-times oversampling the resulting wave front would still use up 452 GByte of memory. However, these values do not take into account the array representing the phantom - which is usually upsampled to the resolution of the wave front and, thus, further increases the memory footprint by a factor of 2 - and the overhead caused by the memory management of Python. Therefore, the standard approach with a two-dimensional wave front would be impracticable even on many currently used high-end compute servers. In contrast, the linewise parallelized simulation performed at approx. 500 MByte per CPU core, including also the representation of the mammography phantom. The resulting DPC projection image is presented in Fig. 4(b).
As the simulation of the full phantom with a 2D wave front was not possible because of its size, ten smaller parts of 8×8 mm2 at random positions were compared. The maximum absolute difference of the DPC signals that was observed is 0.000387 rad. This equals a relative error of 0.236% at that point (absolute difference between 1D and 2D wave-front simulation divided by the DPC signal of the 2D simulation). The discrepancies were found at edges not perpendicular or parallel to the grating bars. The same procedure was repeated with a finer sampled version of the phantom (50×50×50 μm3 voxel size). Here, the values were larger (0.00217 rad and 0.66%) but still reasonably small. The difference between the two phantom resolutions is caused by the influence of the linear interpolation on the edge steepness, when the volume is upsampled to the resolution of the wave front.
Exploiting the approach even more by propagating only one wave-front line per sample layer leads to a drastic reduction of the execution time: running on 25 cores on the utilized compute server (4 AMD Opteron 6274 2.2 GHz 16-core CPUs, 512 GByte RAM, Ubuntu Linux 13.04 64-bit) the monochromatic simulation of the phantom takes less than three minutes including the flat-field and processing, a polychromatic simulation using a tungsten spectrum with 73 energy bins based on an adapted version of Spektr [20] is finished in less than 2 hours 20 minutes (48 cores).
4. Discussion
With this study we demonstrated an approach to overcome the limitations that computer hardware imposes on wave-optical GI simulations. Processing of a 1D wave front relaxes the requirements for RAM size to roughly the square root of the value for the entire two-dimensional wave front (assuming square shape). Thus, it is possible to simulate much larger samples than with the regular two-dimensional simulation at a given resolution. Alternatively, for a given sample the grid spacing of the wave front can be refined to improve results by a higher sampling rate [21].
The influence of the y-component on the wave propagation is in this case not accounted for as the single lines of the wave front are processed independently from another. Therefore, results differ slightly from propagating a 2D wave front, as presented in Figs. 2 and 3. These differences are caused by intensity variations at the edges, as shown in Fig. 2 (Fresnel fringes), in combination with the binning of the wave front at the detector.
An exact quantification of the error is difficult as it depends on a large number of factors, such as e.g. size and material of the sample, its position relative to the detector pixel, or the propagation distance. As demonstrated in Fig. 3 the orientation of the structure plays an important role. For samples that are constant in the y-direction, for example a cylinder oriented along the y-axis, the 1D approach is able to perfectly reproduce the result of the 2D wave front as the phase is only shifted in x-direction. Moreover, if the same cylinder is tilted by 90 degrees the error will again vanish as the GI is not sensitive to shifts in y-direction and the DPC signal itself is zero. For all orientations between zero and 90 degrees an error is introduced. In order to keep the influence of this effect reasonably small one detector pixel should average over several periods of the interference pattern.
Simulations of a flat-field - i.e. simulations without a sample in the GI, which are needed to obtain the DPC signal - or of samples with zero gradient in y-direction, consequently, highly profit from the 1D approach in terms of speed without any reduction of result accuracy. Samples with a more variable structure, like the breast phantom, show deviations from the simulation of a 2D wave front at edges that are neither perpendicular nor parallel to the grating bars. Although the differences remained small, in practice it may be advisable to use the approach to generate an overview image and validate quantitative values by simulating small regions of interest that contain sharp tilted edges with a 2D wave front.
For any sample the execution time can be significantly reduced by parallel calculation of the single lines. Also, phantoms that are described as a discretized volume (Fig. 4) can be processed much faster if they are upsampled to the resolution of the wave front only in x-direction instead of both directions. In the presented case this led to a further speed-up of a factor of about 1000.
The dark-field signal of the sphere in Fig. 3 yielded similar results as the comparison of the DPC signals. However, a relevant sample containing small scatterers like in [14] was not explicitly tested in this study. In the future, the simulation will also be expanded to fan-beam geometry. Here, the approach also applies if the opening angle lies in the x–z-plane.
5. Conclusion
We presented an approach to significantly reduce the hardware requirements of wave-optical simulations of X-ray phase-contrast imaging using grating interferometry. Splitting the two-dimensional wave front up into single lines and processing them separately allows the investigation of samples and fields of view of clinically relevant sizes and the qualitative comparison of different setup configurations. We also discussed the limitations of the approach regarding orientation of the sample structure relative to the direction of the grating bars. Especially for flat-field simulations and samples with a non-zero gradient only in x-direction, a significant speed-up can be achieved without any alteration of the result.
Acknowledgments
J.W. would like to thank the Graduate School of Information Science in Health (GSISH) for the financial support of his studies. M.S. would like to thank the TUM Graduate School. We gratefully acknowledge Dr. Predrag Bakic for providing the software breast phantom. This work has been partly funded by the Bayerisches Staatsministerium für Wirtschaft, Infrastruktur, Verkehr und Technologie under contract number 13 30 / 89265 / 5 / 12. We acknowledge financial support through the DFG Cluster of Excellence Munich-Centre for Advanced Photonics (MAP), the DFG Gottfried Wilhelm Leibniz program and the European Research Council (ERC, FP7, StG 240142).
References and links
1. A. Momose, “Phase-sensitive imaging and phase tomography using x-ray interferometers,” Opt. Express 11, 2303–2314 (2003). [CrossRef] [PubMed]
2. 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). [CrossRef] [PubMed]
3. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance x-ray sources,” Nature Phys. 2, 258–261 (2006). [CrossRef]
4. F. Pfeiffer, O. Bunk, C. David, M. Bech, G. L. Duc, A. Bravin, and P. Cloetens, “High-resolution brain tumor visualization using three-dimensional x-ray phase contrast tomography,” Phys. Med. Biol. 52, 6923 (2007). [CrossRef] [PubMed]
5. T. Donath, F. Pfeiffer, O. Bunk, C. Grünzweig, E. Hempel, S. Popescu, P. Vock, and C. David, “Toward clinical x-ray phase-contrast ct: Demonstration of enhanced soft-tissue contrast in human specimen,” Invest. Radiol. 45, 445–452 (2010). [PubMed]
6. M. Stampanoni, Z. Wang, T. Thüring, C. David, E. Roessl, M. Trippel, R. A. Kubik-Huch, G. Singer, M. K. Hohl, and N. Hauser, “The first analysis and clinical evaluation of native breast tissue using differential phase-contrast mammography,” Invest. Radiol. 46, 801–806 (2011). [CrossRef] [PubMed]
7. S. Grandl, M. Willner, J. Herzen, D. Mayr, S. D. Auweter, A. Hipp, F. Pfeiffer, M. Reiser, and K. Hellerhoff, “Evaluation of phase-contrast CT of breast tissue at conventional X-ray sources - presentation of selected findings,” Z. Med. Phys. 23, 212–221 (2013). [CrossRef] [PubMed]
8. A. Sztrókay, J. Herzen, S. D. Auweter, S. Liebhardt, D. Mayr, M. Willner, D. Hahn, I. Zanette, T. Weitkamp, K. Hellerhoff, F. Pfeiffer, M. F. Reiser, and F. Bamberg, “Assessment of grating-based X-ray phase-contrast CT for differentiation of invasive ductal carcinoma and ductal carcinoma in situ in an experimental ex vivo set-up,” Eur. Radiol. 2, 381–387 (2012).
9. F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. F. Eikenberry, C. Bronnimann, C. Grünzweig, and C. David, “Hard-x-ray dark-field imaging using a grating interferometer,” Nat. Mater. 7, 134–137 (2008). [CrossRef] [PubMed]
10. S. Schleede, F. G. Meinel, M. Bech, J. Herzen, K. Achterhold, G. Potdevin, A. Malecki, S. Adam-Neumair, S. F. Thieme, F. Bamberg, K. Nikolaou, A. Bohla, A. Yildirim, R. Loewen, M. Gifford, R. Ruth, O. Eickelberg, M. Reiser, and F. Pfeiffer, “Emphysema diagnosis using x-ray dark-field imaging at a laser-driven compact synchrotron light source,” Proc. Natl. Acad. Sci. U.S.A. 109, 17880–17885 (2012). [CrossRef] [PubMed]
11. A. Yaroshenko, F. G. Meinel, M. Bech, A. Tapfer, A. Velroyen, S. Schleede, S. Auweter, A. Bohla, A. O. Yildirim, K. Nikolaou, F. Bamberg, O. Eickelberg, M. F. Reiser, and F. Pfeiffer, “Pulmonary emphysema diagnosis with a preclinical small-animal x-ray dark-field scatter-contrast scanner,” Radiology 269, 427–433 (2013). [CrossRef] [PubMed]
12. T. Michel, J. Rieger, G. Anton, F. Bayer, M. W. Beckmann, J. Durst, P. a. Fasching, W. Haas, A. Hartmann, G. Pelzer, M. Radicke, C. Rauh, A. Ritter, P. Sievers, R. Schulz-Wendtland, M. Uder, D. L. Wachter, T. Weber, E. Wenkel, and A. Zang, “On a dark-field signal generated by micrometer-sized calcifications in phase-contrast mammography,” Phys. Med. Biol. 58, 2713–2732 (2013). [CrossRef] [PubMed]
13. P. Bartl, J. Durst, W. Haas, T. Michel, A. Ritter, T. Weber, and G. Anton, “Simulation of x-ray phase-contrast imaging using grating-interferometry,” in “Nuclear Science Symposium Conference Record (NSS/MIC), 2009 IEEE,” (2009), pp. 3577–3580.
14. A. Malecki, G. Potdevin, and F. Pfeiffer, “Quantitative wave-optical numerical analysis of the dark-field signal in grating-based X-ray interferometry,” EPL 99, 48001 (2012). [CrossRef]
15. J. Goodman, Introduction to Fourier Optics (Roberts & Company Publishers, Greenwood Village CO, USA, 2005), 3rd ed.
16. D. M. Paganin, Coherent x-ray optics, Oxford series on synchrotron radiation (Oxford University Press, Oxford, 2006). [CrossRef]
17. D. D. Pokrajac, A. D. A. Maidment, and P. R. Bakic, “Optimized generation of high resolution breast anthropomorphic software phantoms,” Med. Phys. 39, 2290–2302 (2012). [CrossRef] [PubMed]
18. International Commission on Radiation Units and Measurements, Tissue substitutes in radiation dosimetry and measurement. ICRU Report 44 (International Commission on Radiation Units and Measurements, Bethesda MD, USA, 1989).
19. H. Q. Woodard and D. R. White, “The composition of body tissues,” Br. J. Radiol. 59, 1209–1219 (1986). [CrossRef] [PubMed]
20. J. H. Siewerdsen, a. M. Waese, D. J. Moseley, S. Richard, and D. a. Jaffray, “Spektr: A computational tool for x-ray spectral analysis and imaging system optimization,” Med. Phys. 31, 3057 (2004). [CrossRef] [PubMed]
21. F. A. Vittoria, P. C. Diemoz, M. Endrizzi, L. Rigon, F. C. Lopez, D. Dreossi, P. R. T. Munro, and A. Olivo, “Strategies for efficient and fast wave optics simulation of coded-aperture and other x-ray phase-contrast imaging methods,” Appl. Opt. 52, 6940–6947 (2013). [PubMed]
