Edge illumination X-ray phase tomography of multi-material samples using a single-image phase retrieval algorithm

In this paper we present a single-image phase retrieval algorithm for multi-material samples, developed for the edge illumination (EI) X-ray phase contrast imaging method. The theoretical derivation is provided, along with any assumptions made. The algorithm is evaluated quantitatively using both simulated and experimental results from a computed tomography (CT) scan using the EI laboratory implementation. Qualitative CT results are provided for a biological sample containing both bone and soft-tissue. Using a single EI image per projection and knowledge of the complex refractive index, the algorithm can accurately retrieve the interface between a given pair of materials. A composite CT slice can be created by splicing together multiple CT reconstructions, each retrieved for a different pair of materials. Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. OCIS codes: (340.7430) X-ray coded apertures; (340.7440) X-ray imaging; (100.5070) Phase retrieval; (100.6950) Tomographic image processing. References and links 1. S. W. Wilkins, Y. I. Nesterets, T. E. Gureyev, S. C. Mayo, A. Pogany, and A. W. Stevenson, “On the evolution and relative merits of hard X-ray phase-contrast imaging methods,” Phil. Trans. R. Soc. 372, 20130021 (2014). 2. A. Bravin, P. Coan, and P. Suortti, “X-ray phase-contrast imaging: from pre-clinical applications towards clinics,” Phys. Med. Biol. 58, R1–35 (2013). 3. 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, 5486–5492 (1995). 4. S. W. Wilkins, T. E Gureyev, D. Gao, A. Pogany and A. W. Stevenson, “Phase-contrast imaging using polychromatic hard x-rays,” Nature 384, 335–338 (1996). 5. S. C. Mayo, A. W. Stevenson, and S. W. Wilkins, “In-line phase-contrast x-ray imaging and tomography for materials science,” Materials 5, 937–965 (2012). 6. M. Topperwien, M. Krenkel, D. Vincenz, F. Stober, A. M. Oelschlegel, J. Goldschmidt and T. Salditt, “ Threedimensional mouse brain cytoarchitecture revealed by laboratory-based x-ray phase-contrast tomography,” Sci. Rep. 7, 42847 (2017). 7. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai and Y. Suzuki, “ Demonstration of x-ray talbot interferometry,” Jpn. J. Appl. Phys. 42, L866–L868 (2003). 8. F. Pfeiffer, T. Weitkamp, O. Bunk and C. David, “Phase retrieval and differential phase-contrast imaging with lowbrilliance x-ray sources,” Nature Phys. 2, 258–261 (2006). 9. F. Pfeiffer, C. Kottler, O. Bunk and C. David, “Hard x-ray phase tomography with low-brilliance sources,” Phys. Rev. Lett. 98, 108105 (2007). 10. K. S. Morgan, D. M. Paganin and K. K.W. Siu, “Quantitative single-exposure x-ray phase contrast imaging using a single attenuation grid,” Opt. Exp. 19, 19781–19789 (2011). 11. D. Macindoe, M. J. Kitchen, S. C. Irvine, A. Fouras and K. S. Morgan, “Requirements for dynamical differential phase contrast x-ray imaging with a laboratory source,” Phys. Med. Biol. 61, 8720–8735 (2016). 12. M. Endrizzi, F. A. Vittoria, G. Kallon, D. Basta, P. C. Diemoz, A. Vincenzi, P. Delogu, R. Bellazzini and A. Olivo, “Achromatic approach to phase-based multi-modal imaging with conventional x-ray sources,” Opt. Exp. 23, 16473– 16480 (2015). 13. A. Olivo, F. Arfelli, G. Cantatore, R. Longo, R. H. Menk, S. Pani, M. Prest, P. Poropat, L. Rigon, G. Tromba, E. Vallazza and E. Castelli, “An innovative digital imaging set-up allowing a low-dose approach to phase contrast Vol. 25, No. 10 | 15 May 2017 | OPTICS EXPRESS 11984 #286760 https://doi.org/10.1364/OE.25.011984 Journal © 2017 Received 15 Feb 2017; revised 5 Apr 2017; accepted 6 Apr 2017; published 12 May 2017 applications in the medical field,” Med. Phys. 28, 1610–1619 (2001). 14. A. Olivo and R. Speller, “Modelling of a novel x-ray phase contrast imaging technique based on coded apertures,” Phys. Med. Biol. 52, 6555–6573 (2007). 15. P. C. Diemoz, M. Endrizzi, C. E. Zapata, Z. D. Pesic, C. Rau, A. Bravin, I. K. Robinson, and A. Olivo, “X-ray phase-contrast imaging with nanoradian angular resolution,” Phys. Rev. Lett. 110, 138105 (2013). 16. C. K. Hagen, P. R. T. Munro, M. Endrizzi, P. C. Diemoz, and A. Olivo, “Low-dose phase contrast tomography with conventional x-ray sources,” Med. Phys. 41, 070701 (2014). 17. A. Olivo and R. Speller, “A coded-aperture technique allowing x-ray phase contrast imaging with conventional sources,” Appl. Phys. Lett. 91, 074106 (2007). 18. P. C. Diemoz and A. Olivo, “On the origin of contrast in edge illumination x-ray phase-contrast imaging,” Opt. Exp. 22, 28199–28214 (2014). 19. A. Olivo and R. Speller, “Image formation principles in coded-aperture based x-ray phase contrast imaging,” Phys. Med. Biol. 53, 6461–6474 (2008). 20. P. R. T. Munro, K. Ignatyev, R. D. Speller, and A. Olivo, “Phase and absorption retrieval using incoherent x-ray sources,” Proc. Natl. Acad. Sci. 109, 13922–13927 (2012). 21. P. C. Diemoz, F. A. Vittoria, C. K. Hagen, M. Endrizzi, P. Coan, E. Brun, U. H. Wagner, C. Rau, I. K. Robinson, A. Bravin, and A. Olivo, “Single-image phase retrieval using an edge illumination x-ray phase-contrast imaging setup,” J. Synchrotron Radiat. 22, 1072–1077 (2015). 22. D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Miller, and S. W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206, 33–40 (2002). 23. P. C. Diemoz, C. K. Hagen, M. Endrizzi, M. Minuti, R. Bellazzini, L. Urbani, P. De Coppi, and A. Olivo, “Singleshot X-ray phase-contrast computed tomography with non-microfocal laboratory sources,” Phys. Rev. Appl. (2017, in press). 24. M. A. Beltran, D. M. Paganin, K. Uesugi, and M. J. Kitchen, “2D and 3D X-ray phase retrieval of multi-material objects using a single defocus distance,” Opt. Express 18, 6423–36 (2010). 25. M. A. Beltran, D. M. Paganin, K. K. W. Siu, A. Fouras, S. B. Hooper, D. H. Reser, and M. J. Kitchen, “Interfacespecific x-ray phase retrieval tomography of complex biological organs,” Phys. Med. Biol. 56, 7353–7369 (2011). 26. T. E. Gureyev, Y. I. Nesterets, A. W. Stevenson, P. R. Miller, A. Pogany, and S. W. Wilkins, “Some simple rules for contrast, signal-to-noise and resolution in in-line x-ray phase-contrast imaging,” Opt. Express 16, 3223–3241 (2008). 27. P. Diemoz, F. Vittoria, and A. Olivo, “Spatial resolution of edge illumination x-ray phase-contrast imaging,” Opt. Express 22, 15514–29 (2014). 28. M. Reed Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. 73, 1434 (1983). 29. D. M. Paganin, Coherent X-Ray Optics, (Oxford University, 2009). 30. P. Munro and A. Olivo, “X-ray phase-contrast imaging with polychromatic sources and the concept of effective energy,” Phys. Rev. A 87, 053838 (2013). 31. 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). 32. International Commission on Radiation Units and Measurements, “Tissue substitutes in radiation dosimetry and measurement,” ICRU Report 44 (1989). 33. C. K. Hagen, P. C. Diemoz, M. Endrizzi, and A. Olivo, “The effect of the spatial sampling rate on quantitative phase information extracted from planar and tomographic edge illumination x-ray phase contrast images,” J. Phys. D. Appl. Phys. 47, 455401 (2014). 34. R. Bellazzini, G. Spandre, A. Brez, M. Minuti, M. Pinchera, and P. Mozzo, “Chromatic X-ray imaging with a fine pitch CdTe sensor coupled to a large area photon counting pixel ASIC,” J. Instrum. 8, C02028 (2013). 35. M. Endrizzi, P. C. Diemoz, C. K. Hagen, F. A. Vittoria, P. R. T. Munro, L. Rigon, D. Dreossi, F. Arfelli, F. C. M. Lopez, R. Longo, M. Marenzana, P. Delogu, A. Vincenzi, L. D. Ruvo, G. Spandre, A. Brez, R. Bellazzini, and A. Olivo, “Edge-illumination X-ray phase contrast imaging: matching the imaging method to the detector technology,” J. Instrum. 9, C11004 (2014). 36. P. R. T. Munro, K. Ignatyev, R. D. Speller and A. Olivo, “Design of a novel phase contrast x-ray imaging system for mammography,” Phys. Med. Biol. 55, 4169–4185 (2010). 37. N. Otsu, “A threshold selection method from gray-level histograms,” IEEE Trans. Syst. Man. Cybern. 9, 62–66 (1979). 38. A. Horng, E. Brun, A. Mittone, S. Gasilov, L. Weber, T. Geith, S. Adam-neumair, S. D. Auweter, A. Bravin, M. F. Reiser, and P. Coan, “Cartilage and soft tissue imaging using x-rays,” Invest. Radiol. 49, 627–634 (2014). 39. K. Ignatyev, P. R. T. Munro, R. D. Speller, and A. Olivo, “Effects of signal diffusion on x-ray phase contrast images,” Rev. Sci. Instrum. 82, 073702 (2011). Vol. 25, No. 10 | 15 May 2017 | OPTICS EXPRESS 11985


Introduction
X-ray phase contrast tomography is a technique capable of reconstructing the phase map of an object, often leading to an increased image contrast when compared with conventional, absorption-based X-ray computed tomography (CT), with the added benefit of a potential reduction in radiation dose [1].This makes it a desirable tool for the investigation of low-absorbing samples, which naturally suffer from weak contrast in conventional CT, and particularly those constrained by dose limits, e.g.biological samples.Many phase contrast methods require a highly coherent X-ray source and are therefore limited to synchrotrons, however, several techniques have been implemented with laboratories sources [2].In propagation-based imaging (PBI), although the technique relies on spatial coherence to detect interference patterns [3], the development of microfocus sources enabled its implementation outside synchrotron facilities [4,5].While microfocus sources are mostly associated with reduced flux, therefore leading to long exposures, promising new technologies such as the liquid-metal-jet source can overcome this limitation [6].In the grating interferometry method, phase sensitivity is achieved by employing two gratings to create and detect the Talbot self-imaging effect [7].Spatial coherence is therefore required, however, the technique can be made compatible for use with extended laboratory sources, provided that a third, "source" grating is used [8,9].Another technique, the single grid method, employs an absorption (or phase) grid to create a reference beam, which is later compared to the detected beam in the presence of a sample, in order to extract information about the sample's complex refractive index [10].While the method was first designed for use with synchrotron radiation, a recent study explored the possibility of adapting it for use with microfocus, polychromatic sources to perform dynamical imaging [11].The edge illumination (EI) technique, considered in this study, is a non-interferometric, achromatic imaging method [12], which can be used with either synchrotron radiation or with polychromatic, extended sources (with spot size up to ~100 μm) [13][14][15][16].
The working principles of EI have been discussed in detail previously [17,18]; in general, in the laboratory, phase sensitivity is achieved by employing sets of apertures (masks), placed behind and in front of the sample.While the first mask divides the incoming, divergent beam into physically separated beamlets, the second mask is placed in contact with the detector, creating insensitive regions between neighboring pixels.By slightly misaligning the masks such that only part of each beamlet reaches the active area of the pixel, refraction due to the sample can be detected, since the latter will deflect the beamlet towards or away from the pixel, therefore leading to an increase or decrease in the number of detected photons.[18,19] Until recently, in EI, at least two frames per projection were required to extract an object's phase information, usually acquired at different positions of the sample-mask [20].Diemoz et al. have recently developed a phase retrieval algorithm where a single image per projection is sufficient [21].In a similar way to the algorithm developed by Paganin et al. for PBI [22], Diemoz's algorithm for EI relies on several assumptions, the main one being that the sample is quasi-homogeneous (i.e. it has a constant ratio between the real and imaginary parts of its complex refractive index).Using projection data acquired with synchrotron radiation, the algorithm has been shown to provide quantitative results when applied to a homogeneous sample, and enhanced image quality when imaging complex biological samples [21].
Translation of this single-image retrieval algorithm to CT scans would have significant impact, as it would considerably reduce scan times and simplify the acquisition sequence, as well as possibly lead to a dose reduction.Indeed, recent developments have made the algorithm suitable for use on CT data acquired in a laboratory environment, with an extended and polychromatic X-ray source [23].While high image quality was observed for a range of samples which can be considered approximately homogeneous (e.g.rat heart), many biological samples of interest are not quasi-homogeneous (e.g.contain both soft-tissue and bone), and applying the algorithm to such samples will result in under or over-retrieved characteristic features (dark/bright fringes, blurring).For example, in a paper by Diemoz et al., this can be seen as blurring around a chalk detail embedded in a phantom containing various plastic materials [23].
In this paper, we present an extension to the retrieval algorithm developed by Diemoz et al., which enables EI phase retrieval of multi-material samples.Our new algorithm follows the work of Beltran et al. who extended Paganin's single-image retrieval algorithm for PBI, for the case of a multi-material object [24,25].While still requiring only a single image per projection, our retrieval algorithm can correctly reconstruct the interface between any two materials (when one is encased by the other), by tuning the relevant input values.Using a "splicing" method presented here, when a sample is made of various materials, the algorithm can be applied multiple times, once for each pair of materials, and a composite CT slice can then be obtained which will simultaneously display all materials adequately.

Theory
In this section we will derive the theory for EI single-image phase retrieval for samples made of multiple materials.A schematic of a laboratory EI system is shown in figure 1.A divergent X-ray beam is divided into beamlets by a mask made of alternating absorbing and transmitting lines, placed immediately before the sample.A second mask is placed before the detector, with a corresponding and magnified period and aperture size.When both masks are misaligned with respect to each other, refraction of the beamlets caused by the sample in the direction perpendicular to the masks' apertures is converted into detectable intensity changes.These intensity changes are proportional to the first derivative of the phase, φ [20].Image formation in the direction parallel to the masks' apertures can be treated as PBI, where the signal recorded is proportional to the second derivative of the phase [26,27].It should be noted however that in our case, the PBI signal is very weak due to the use of an extended laboratory source.Assuming a near-field regime, and that the sample's attenuation and phase are varying slowly, the signal can be expressed by the transport-of-intensity equation [26,28].
In their work, Diemoz et al. have shown that the normalized signal recorded by an EI setup is given by [21]: where T = exp(−2k β(x, y, z)dz) and φ = −k δ(x, y, z)dz are the transmission and phase shift caused by the sample, respectively.The X-ray wave-number is represented by k and the complex refractive index is n = 1 − δ + i β.The illumination curve, C(y e ), describes the measured intensity as a function of the relative masks' displacement, y e , without a sample in the field-of-view.The sample-to-detector distance is denoted by z 2 , and LSF y is the detector's line spread function along the direction perpendicular to the masks' apertures.Differentiation with respect to x or y is represented by ∇ x,y and * indicates convolution.When considering a homogeneous sample which is sufficiently thin as to satisfy the projection approximation [29], the transmission and phase at the sample's exit surface can be expressed as a function of the sample's thickness, t(x, y): φ(x, y) = −kδt(x, y) .
Diemoz et al. have shown that by substituting Eqs. ( 2) and (3) into Eq.( 1) and defining J EI = z 2 C (y e )C −1 (y e ), an expression for the sample's projected thickness can be obtained: where μ = 2k β is the linear attenuation coefficient, F and F −1 represent the two dimensional Fourier transform and its inverse, k x = 2π f x and k y = 2π f y where f x and f y are the Fourier space coordinates, and the system's modulation transfer function along y is given by MT F y (k y ) = F{LSF y }.When a polychromatic beam is used, both δ and μ should be evaluated at their effective energies [23,30].Therefore, when a homogeneous sample of known material is imaged with monochromatic radiation, the algorithm will retrieve quantitative values for the object's thickness.While for the case of inhomogeneous samples the retrieved results would not be quantitative, often high image quality can still be achieved by tuning the input δ and μ values following an initial guess.For simplicity, here we will derive the multi-material retrieval algorithm assuming monochromatic radiation.Here, we adapt the derivation first presented for PBI by Beltran et al. for a ternary object (composed of two materials and voids) [24], to the case of EI.Consider an object made of two materials denoted by 1 & j, where material j is fully encased by material 1.The corresponding refractive indices are n 1 = 1 − δ 1 + i β 1 and n j = 1 − δ j + i β j , with β 1 = μ 1 /2k and β j = μ j /2k.For monochromatic radiation, Eq. ( 1) describes the detected normalized signal.Assuming the sample satisfies the projection approximation, the transmission and phase shift at the sample's exit plane are given by: (5) where t 1 (x, y) and t j (x, y) are the projected thicknesses of materials 1 and j in the direction of wave propagation, z.Let us now introduce the total projected thickness, defined as A(x, y) = t 1 (x, y) + t j (x, y).This allows us to express the transmission and derivatives of the phase as: where dependencies on x and y have been discarded for notation simplicity.By substituting Eqs. ( 7) and ( 8) into Eq.( 1), a solution for t j can be obtained.Eq. ( 7) is in fact the first term on the right-hand-side (RHS) of Eq. ( 1).The second term on the RHS becomes: where we have used the identity: and made the assumption that the projected thickness of the encasing material is slowly varying (i.e. that the first term on the RHS of Eq. ( 8) is negligible compared to the second term), allowing us to ignore the terms containing its spatial derivatives.This assumption was previously made by Beltran et al. [24], and while it is violated in certain cases (e.g.near the edges of the encasing material, where phase gradients due to the latter are non-negligible), errors resulting from these violations are localized and should not hinder the retrieval of t j as long as the latter is not in contact, or in the immediate vicinity of another interface.
By taking the two-dimensional Fourier transform of Eq. ( 13) and making use of the Fourier derivative theorem, we get the following expression: An expression for t j can then be obtained by inverse transforming both sides of Eq. ( 14) and rearranging: For a given experimental setup, Eq. ( 15) can be readily implemented to obtain quantitative results of the projected thickness of material j, provided that the total projected thickness, A, and both materials' refractive indices are known.The total projected thickness can be easily estimated in common situations where the sample is placed in a cylindrical, full container (as is the case in both examples provided below).However, in other cases where there are voids to consider, A(x, y) can be estimated by forward-projecting a CT slice reconstructed by using Eq. ( 4) for material 1, and setting a threshold to separate object from void, as was previously suggested by Beltran et al. [24].
For samples containing multiple materials j = 2, 3, ...N, Eq. ( 15) can be applied N − 1 times, each time adjusting Δμ and γ such that they correspond to the pair of materials of interest.As long as each material j is fully encased by material 1, the interface between the two materials will be correctly reconstructed.Independent N CT slices can be obtained by applying a conventional filtered-backprojection (FBP) algorithm to sinograms of t 1 and t j=2,3,... N .However, each of these slices will provide an accurate reconstruction of the corresponding pairs of materials, while other interfaces will be either under or over-retrieved (appearing either as residual fringes or as blurred interfaces).A composite CT slice, with all interfaces sharply reconstructed, can then be obtained by splicing together these CT slices corresponding to all different materials.The splicing procedure consists of extracting from each slice the area containing the pair of materials focused on in the retrieval, and digitally inserting these into one composite slice, after the adjustment of background offsets.
Note that while Eq. ( 15) was derived for the case of monochromatic radiation, an analogous expression could be developed for polychromatic beams.In a similar way to Diemoz's work [23], since both δ and β vary with energy, a polychromatic extension would involve expressing the measured detector signal as a weighted sum of all its monochromatic components, taking into account the source spectrum and the detector's energy response.Therefore, as long as there is no significant beam hardening by the sample, Eq. ( 15) can also be applied to data acquired with a polychromatic beam, if the input values for δ and β are calculated at their effective energies [30].

Quantitative validation
The novelty of the proposed algorithm lies in its ability to enable the imaging of samples containing materials with significantly different refractive index properties.In order to test its validity, a simulation study was performed, followed by the collection of experimental data for direct comparison.As a test sample, a water filled cylinder containing aluminium and low-density polyethylene (LDPE) rods was chosen.
A wave-optics simulation, previously developed to replicate the experimental EI setup [31], was used with a photon energy of E = 17.5 keV (corresponding to the k-alpha line of the molybdenum target in the used laboratory X-ray tube).A numerical phantom was created, which consisted of a vertically-aligned cylinder of water (diameter = 1.85 cm) containing an aluminium rod (diameter = 2 mm) and a LDPE rod (diameter = 3.9 mm).Corresponding refractive index values were taken from the ICRU 44 database [32] and are reported in table 1. Mask parameters were chosen to match those used in the laboratory, where each mask is 150 μm thick and is made of gold.The masks' periods and aperture sizes were 48 μm and 12 μm, respectively, for the sample-mask, and 62 μm and 15 μm for the detector-mask.The masks' apertures were aligned with the vertical direction.The detector pixel size was 62 × 62 μm 2 .The source-to-sample distance was z 1 = 1.6 m and the sample-to-detector distance was z 2 = 0.4 m.For simplicity, in the simulation we assumed that there is no gap between the sample-mask and sample, and similarly no gap between the detector-mask and detector; in truth, a small (~1-2 cm) gap is present which is not expected to have any noticeable effect on the images.An illumination curve was simulated to obtain C(y e ) and C −1 (y e ).A CT scan was then simulated, consisting of 2400 views over 360 degrees (i.e.angular step = 0.15 degrees).The relative displacement of the sample-mask with respect to the detector-mask was set as y e = 10 μm.To increase the spatial sampling rate, at each angular view 4 images were taken at different sub-pixel positions of the sample and were later recombined to form a high-resolution image.This process is often referred to as "dithering" in EI, and is not required unless an increase in spatial resolution over that determined by the detector pixel size is sought [33].The phase retrieval part comprised of three steps.First, Eq. ( 4) was applied to the raw projections using the refractive index values of water.The second and third steps consisted of applying Eq. ( 15) to raw projections, once for the water-aluminium interface, and then for the water-LDPE interface, using the values listed in table 1.In this case, A(x, y) was calculated as the projection of a circle with a diameter corresponding to that of the water cylinder.Three separate sinograms were created, corresponding to the different retrieval steps.CT slices of each of these phase-retrieved sinograms were reconstructed by means of a FBP algorithm.
For experimental data collection, a phantom was created which matched the simulated one (apart from the fact that a hollow plastic cylinder of unknown material, with inner and external diameters of 1.75 cm and 1.85 cm, respectively, was filled with water for practical reasons).As a source, the Rigaku MicroMax 007 HF rotating anode (molybdenum) X-ray tube (Rigaku Corporation, Japan) with a focal spot of approximately 70 μm was used.The source was operated at 25 mA and 40 kVp, with a 30 μm Molybdenum filter.A Pixirad single photon counting, energyresolving detector [34,35], with a pixel size of 62 × 62 μm 2 was used and placed 2.07 m away from the source.The sample and detector masks were placed at 1.6 m and 1.97 m downstream of the source, respectively, following studies from the method's early days in which distances between the source, masks and detector were optimized [14,36].Both masks were fabricated by electroplating gold strips onto a graphite substrate (Creatv Microtech Inc., Potomac, MD, USA).The masks' period and aperture sizes were the same as used in the simulation, apart from the period of the detector-mask which was 59 μm, to compensate for the limited magnification arising from the small gap between the detector-mask and detector.An illumination curve scan was performed, followed by a CT scan of the phantom with the same parameters reported for the simulation.Flat-field images (i.e.images in the absence of the sample) were acquired at each angular view, and were later used to normalize the sample images.The exposure time per projection was 4 s.
For phase retrieval, Eqs. ( 4) and ( 15) were applied to the normalized data.Here, since the source spectrum is polychromatic and the complex refractive index varies with energy, the refractive index values used in the retrieval were estimated.In general, these can be estimated using the concept of effective energies if the source spectrum, detector energy response and sample composition and geometry are known [30].In the case discussed in this paper, one can consider that an optimal CT reconstruction has been achieved when two conditions are satisfied.The first one is the maximization of the sharpness of the interface between the pair of materials under question.The second is that the reconstructed material j has a mean value of 1 in the CT slice.In order to simultaneously fulfill these conditions, CT reconstruction is carried out multiple times while varying the refractive index values.The chosen values used in the retrieval algorithm of the experimental data are presented in table 2, along with the effective energies these correspond to for each material.These effective energies are reasonable considering the source spectrum produced by the molybdenum target.The difference between the values obtained for water and LDPE compared to aluminium can be explained by taking into account the beam hardening caused by the higher absorption of the aluminium rod.Phase retrieved sinograms were created, followed by a CT reconstruction using FBP.For a quantitative evaluation of the proposed algorithm, simulated and experimental projected thicknesses of each material were plotted together and are shown in Fig. 2. The plots demonstrate the quantitativeness of the algorithm: all simulated profiles retrieve the nominal thickness values accurately (retrieving diameters of 1.85 cm, 2 mm and 3.9 mm for water, aluminium and LDPE, respectively).As expected, in each case the profiles are quantitative only for the interface between the two materials targeted by the algorithm through selection of the relative input parameters.A good agreement can also be noted between the experimental and simulated profiles for each material of interest.Discrepancies in e.g. the retrieved thickness of water, in the location where the aluminium rod is embedded, arise from the difference in (effective) energies between the monochromatic simulation and the experimental polychromatic spectrum.Another inconsistency in retrieved values can be observed between simulated and experimental profiles of both aluminium and LDPE, just near the edges of the cylinder.This inconsistency arises from the fact that while a water cylinder was implemented in the simulation, the experimental phantom consisted of a plastic cylinder filled with water.Finally, by looking at the retrieved projected thickness of both aluminium and LDPE, it can be seen that the algorithm breaks down at the edges of the outer cylinder; this can be expected since in that case our assumption of a slowly varying encasing material thickness is violated.However, this affects the retrieved values only locally, at the edges of the encasing material, and values at the interface of interest are not polluted.
In order to present all interfaces correctly in a single CT slice, the three different slices were spliced together.This step was performed using the open source software ImageJ, for both the experimental and simulated slices.While the segmentation part of the splicing procedure could be automated using methods such as Otsu thresholding [37,38] or by developing expressions for the "bleed-width" in a similar fashion to Beltran's work [24], for this proof-of-concept study the segmentation was done manually, by visually estimating the blurring bleed-width.
Individual simulated CT slices, each displaying reconstructions of the water-air (a), aluminium-water (b) and LDPE-water (c) interfaces are shown in Fig. 3, along with the simulated and experimental spliced slices.Each slice in panels (a-c) demonstrates a sharp reconstruction of the material pair of interest, while residual phase-contrast fringes and blurred interfaces of other materials are highlighted by arrows.In order to achieve a quantitative spliced slice preserving across the phantom the δ values used in the retrieval, each slice reconstructed from t j had to be manipulated according to: slice * t j = [slice t j × (δ j − δ 1 )] + δ 1 , before digitally inserting its segmented area of interest into the corresponding region of slice * t 1 = slice t 1 × δ 1 .The same splicing procedure was applied to both simulation and experimental slices, however, the δ values used for quantitative splicing were different.For each case, the δ values applied during splicing were the same as the ones previously used in the retrieval algorithm, hence some differences exist (most notably for aluminium, see tables 1 and 2).This can be appreciated in Fig. 3(f) where profiles across both spliced slices are plotted together for comparison.Ignoring the discrepancies which have already been discussed above (due to the plastic container and the polychromaticity of the spectrum), both profiles demonstrate a sharp, quantitative reconstruction of all parts of the phantom.These results confirm that the algorithm is indeed quantitative for the case of monochromatic radiation, and can be applied when a polychromatic spectrum is used, if the concept of effective energies is adopted [30].

Qualitative evaluation
Following validation by a phantom study, the proposed retrieval algorithm and splicing method were applied to experimental data of a complex biological sample.The imaged sample was a chicken bone surrounded by soft tissue, placed in a plastic cylinder of approximately 8 mm diameter.The sample was "fresh", i.e. no additional sample preparation was applied (e.g.formalin fixation).The same source used for the phantom scan was used, operated at 40 kVp and 25 mA, with no additional filters.The detector used for this scan was a Hamamatsu C9732DK flat panel detector with a passive-pixel CMOS sensor (Hamamatsu, Japan), and was positioned 2 m away from the source.The detector pixel size was 50 × 50 μm 2 , however the effective pixel size in the x-direction was 100 μm due to the line-skipping design of the detector-mask, where every second detector pixel column is completely covered by the detector-mask.This mask design is used to reduce the negative effect of pixel cross-talk on the EI signal [39].The sample and detector masks were placed at 1.6 m and 1.96 m downstream the source, respectively.The sample-mask period and aperture size were 79 μm and 10 μm, respectively, while the detectormask period was 98 μm with a 17 μm aperture size.Following an illumination curve scan, a CT scan was performed with 720 views over 360 degrees (i.e.0.5 degrees step).To increase the spatial resolution, 8 projections were taken at different sub-pixel positions of the sample, at each angular view.The exposure time per projection was 3 s and the relative position of the sample- mask with respect to the detector-mask was set as y e = 8 μm.Raw projections were corrected for dark current and flat-field inhomogeneities prior to the phase retrieval step, which used Eq. ( 4) once (for the cylinder-air interface) and Eq. ( 15) twice (where the interfaces between soft tissue types, e.g.fat-muscle, and between soft tissue and bone were considered).For the retrieval of the bone-soft tissue interface, refractive index values of bone and water (commonly used as a tissue-equivalent material) at an estimated effective energy of 22 keV were used [32].Here, the effective energy was estimated according to the theoretical source spectrum and a predicted linear detector energy response, since the detector used in this setup was an integrating detector with an energy response function which has yet to be fully characterized (unlike in the phantom scan reported in section 3, where a photon-counting detector was employed).Moreover, for the retrieval of the interface between soft tissue types, since the specific tissue types were not known, and the complex sample geometry in this case did not strictly satisfy the algorithm's condition for a material fully encased by another, a purely qualitative evaluation was performed.Therefore, to maximize the contrast between the different soft tissues, an iterative approach was employed, where Δμ and γ were varied until the contrast was maximized.The reconstructed slices are shown in Fig. 4. On the left hand side, the three separate reconstructions are shown, each rendering a specific interface while exhibiting the expected blurring\fringes at other interfaces.A spliced, composite CT slice is presented in the larger image on the right, in which all different materials are adequately reconstructed.As previously noted, although the algorithm was derived for the case of a material j fully encased within material 1, this is not strictly the case for the chicken bone sample.However, as already observed by Beltran, a locally incorrect choice of refractive index values will only affect the part of the interface for which the choice is not ideal [24].Since quantitative results were not expected, the composite slice was created by windowing each retrieved slice differently, in order to achieve the best visualization.As before, the splicing procedure was performed using ImageJ and manual segmentation.The spliced image shown in Fig. 4 demonstrates that this method could be used to obtain high-quality, albeit not quantitative, images of complex samples with unknown composition, using a single image per projection.

Discussion and conclusion
We have derived an algorithm capable of quantitatively retrieving the projected thickness of multi-material samples illuminated by monochromatic radiation, using a single input image acquired with an EI setup.A CT phantom simulation study validated the algorithm's accuracy, while results from an experimental scan provided good agreement, within the limits of polychromaticity.A "splicing" method was then presented as means for constructing a composite CT slice where all parts of a sample are adequately reconstructed, thereby eliminating artefacts previously arising from locally incorrect choices of refractive index values for multi-material samples, during the phase retrieval step.For strictly quantitative results, one material must be fully encased by the other, and the refractive indices of both materials and the total projected thickness must be known.Therefore, for quantitative imaging, certain applications (e.g. in the field of materials science) could benefit from using the proposed method on data acquired using monochromatic synchrotron radiation.Although we have yet to apply the algorithm to synchrotron radiation data, we expect it to provide more accurate quantitative CT reconstructions, as was demonstrated via simulation.However, we have shown that enhanced image quality and visualization can be achieved for complex samples of unknown materials, imaged in a laboratory environment.
It should be noted that, if the requirement for quantitativeness is relinquished and when splic-

Fig. 1 .
Fig. 1.A top-view schematic of the laboratory implementation of an EI setup.The sample's rotation axis is aligned with the y-direction.

Fig. 2 .
Fig. 2. Retrieved projected thicknesses for water, aluminium and LDPE, using values reported in tables 1 and 2. The red arrows point at the material of focus in each case.The experimental profiles were averaged over 5 pixel rows.

Fig. 3 .
Fig. 3. Simulated CT slices focusing on the water-air (a), aluminium-water (b) and LDPEwater interfaces, using E=17.5 keV.The white arrows indicate artefacts arising from locally incorrect choices of refractive index values.Simulated (d) and experimental (e) spliced slices along with a plot of profiles through them (f), in the position indicated by the red line in (d).The difference in the retrieved values of the aluminium rod in (f, see arrow) is expected and is due to the difference in effective energies (see tables 1 and 2).

Fig. 4 .
Fig. 4. Axial slices of a chicken bone sample reconstructed using values optimized for the cylinder-air interface (a), intra-soft tissue contrast (b) and bone-soft tissue interface (c).Panel (d) shows a composite slice obtained by splicing slices (a-c).