## Abstract

We present a new prior for phase retrieval from X-ray Fresnel diffraction patterns. Fresnel diffraction patterns are achieved by letting a highly coherent X-ray beam propagate in free space after interaction with an object. Previously, either homogeneous or multi-material object assumptions have been used. The advantage of the homogeneous object assumption is that the prior can be introduced in the Radon domain. Heterogeneous object priors, on the other hand, have to be applied in the object domain. Here, we let the relationship between attenuation and refractive index vary as a function of the measured attenuation index. The method is evaluated using images acquired at beamline ID19 (ESRF, Grenoble, France) of a phantom where the prior is calculated by linear interpolation and of a healing bone obtained from a rat osteotomy model. It is shown that the ratio between attenuation and refractive index in bone for different levels of mineralization follows a power law. Reconstruction was performed using the mixed approach but is compatible with other, more advanced models. We achieve more precise reconstructions than previously reported in literature. We believe that the proposed method will find application in biomedical imaging problems where the object is strongly heterogeneous, such as bone healing and biomaterials engineering.

## 1. Introduction

In X-ray in-line phase tomography, the aim is to reconstruct the phase of a coherent X-ray beam, which can then be used to tomographically reconstruct the real part of the refractive index distribution in an object [1]. The main interest in this technique is the several orders of magnitude increase in sensitivity compared with standard, attenuation-based tomography. In the in-line technique, phase contrast is generated simply by letting the beam propagate in free space a relatively short distance after interaction with the object (figure 1) [2]. The measured intensity can be modelled in the Fresnel diffraction formalism; hence, the images recorded can be called Fresnel diffraction patterns. Using these diffraction patterns directly as input to a tomographic reconstruction algorithm yields an edge enhancement effect owing to the phase contrast and is called phase-contrast tomography. Reconstruction of the phase from diffraction patterns is called phase retrieval. Using the resulting phase maps to reconstruct the three-dimensional refractive index distribution is called phase tomography [3].

Several algorithms for phase retrieval from Fresnel diffraction patterns have been developed. Most of these rely on linearization of the Fresnel integral to yield efficient algorithms [1,4–8]. Nonlinear, iterative algorithms have also been proposed for the radiographic case [9–13], and recently also for tomography [14], with an eye to improve spatial resolution of the linear methods.

Owing to the weak transfer of low-frequency information by the Fresnel transform, phase retrieval from Fresnel diffraction patterns is inherently sensitive to noise in the low-frequency range. A typical situation is shown in figure 3*a*. This kind of noise makes image analysis, for example segmentation, very difficult without substantial post-treatment, which might compromise quantitative reconstruction of the refractive index. This issue has previously been addressed by introducing a homogeneous object assumption. This has been done by either setting the phase proportional to the absorption in the transport of intensity equation [5] or introducing a scaled version of the attenuation image as prior knowledge on the phase [15]. The first approach has the advantage of allowing for the use of a single image per projection angle, whereas the second requires a separate measurement of the attenuation (at zero propagation distance) but does not enforce strict proportionality between attenuation and phase. Both these approaches are based on two assumptions: that the object has homogeneous composition (that is a unique ratio between real and imaginary parts of the refractive index; density variations are allowed) and that the composition of the sample is known.

Recently, we proposed to instead introduce the prior knowledge in the object domain, thus allowing for heterogeneous samples where the ratio between attenuation and refractive indices is allowed to vary across the sample. As there are several unknowns entangled in the intensity images—the distance travelled in each material and the real and complex parts of the refractive index, respectively—such priors cannot be applied in the projections but have to be introduced in the object domain. Previously, a prior based on segmentation of the reconstruction of an attenuation scan to achieve an estimate of the object composition was presented [3]. This essentially allows for reconstruction of multi-material objects, i.e. heterogeneous objects consisting of several *a priori* known materials.

Here, we extend this approach to heterogeneous objects. Instead of thresholding the attenuation tomogram, knowledge of the object composition is used to find an empirical relationship between the attenuation and refractive index. The *a priori* phase maps are constructed by first reconstructing the attenuation index, applying the found empirical relationship to yield an initial guess of the refractive index, and finally forward project this estimate at each projection angle. The phase is then retrieved using a phase-retrieval algorithm, for example the mixed approach [6], using the *a priori* phase maps in the regularization term.

## 2. Material and methods

If we assume that imaging is done using a monochromatic, coherent, parallel X-ray beam, the optical properties of an object can be described by its complex refractive index
2.1where *β*(**r**) is the attenuation index distribution in the object, *δ*_{n}(**r**) the refractive index decrement distribution and **r**=(*r*,*s*,*t*) are the spatial coordinates in the object domain. The attenuation index *β* describes the damping of the X-ray beam in the object (owing to absorption and scattering), and the real part of *n*, 1−*δ*_{n}, is related to the phase velocity in the object, thus describing the phase shift induced by the object. Both *δ*_{n} and *β* depend on the material as well as the X-ray wavelength. The phase shift induced by a heterogeneous object can, analogously to the attenuation, be modelled as a straight-line projection. For the attenuation, we have
2.2and similarly for the phase shift
2.3where is the propagation direction of the X-rays. We can then describe the X-ray interaction with the object at each projection angle *θ* as a transmittance function
2.4where *a*_{θ}(**x**) is the amplitude modulation and *φ*_{θ}(**x**) the phase shift.

The beam propagation in free space over a relatively short distance *D* after interaction with the object can be described by the Fresnel transform [15]. If we assume uniform illumination, we have
2.5which is linear with respect to the wavefield. What can be measured, however, is the intensity *I*_{D}(**x**), so that
2.6where *γ*(**x**) accounts for noise and other imperfections in the image (these contributions are not necessarily additive or independent of the object, however). Moving the detector very close to the imaged object, effectively setting *D*=0, yields directly the attenuation image
2.7which is the quantity measured in standard tomography. In general, we can formulate the problem of reconstructing the phase shift from a series of images at different distances as the minimization problem
2.8where is the intensity calculated with some model, for example the mixed approach [6], and *φ*_{θ,0}(**x**) is an *a priori* known solution. Designing this *a priori* knowledge is the focus of this work. Previously, *φ*_{θ,0}(**x**) has been set to 0 [6], which stabilizes the problem but is sensitive to low-frequency noise (figure 3*a*). Setting *φ*_{θ,0}(**x**) proportional to *B*_{θ,0}(**x**) has also been reported [15], so that
2.9where *f*(**x**) is a low-pass filter so that the prior is introduced in the low-frequency range only. This constructs, at each projection angle, a prior phase estimate using the recorded attenuation image. Although this homogeneous object approach has been used successfully in practice [16], the homogeneity assumption can be restrictive (figure 3*b*,*c*). It would be desirable to let the *δ*_{n}/*β*-ratio vary in the sample, but this is not possible on the projection images, because the path length in each material is still unknown.

To allow for a heterogeneous object, the *a priori* knowledge has to be introduced in the object domain. Previously, this has been done by thresholding a reconstructed tomographic attenuation scan, effectively allowing for an object consisting of several discrete materials [3]. Here, we propose instead to set the *δ*_{n}/*β*-ratio directly as a function of *β*. This avoids the thresholding step, and lets the *δ*_{n}/*β*-ratio vary continuously across the object. The sample composition still has to be roughly known, however. Below, we describe the method and show the construction of the prior for two different types of objects.

### (a) Generation of *a priori* phase maps

The *a priori* phase maps are generated based on a reconstructed attenuation tomogram. Reconstructing *β*(**r**) from *I*_{θ,0}(**x**) yields
2.10where *Γ*(**r**) is noise and other imperfections owing to *γ*(**x**). The relationship between attenuation and refractive indices is introduced as a function *m*(**r**), so that the *a priori* estimate of *δ*_{n}(**r**) is obtained as
2.11The design of *m*(**r**) for two types of samples is given below. For phase retrieval, the prior estimate of the phase at each projection angle is obtained by forward projection as
2.12The prior is introduced in the low-frequency range only as before (equation (2.9)) by the application of a low-pass filter.

### (b) Design of prior object estimate

Here, we base the prior estimate on knowledge of the sample composition. We consider two types of objects: a test object constructed from fibres of Al, Al_{2}O_{3}, poly(ethylene terephthalate) (PET) and polypropylene (PP), and bone samples with adjacent soft and gradually mineralizing callus tissue (a typical situation in bone-healing studies). Theoretical values for the linear attenuation coefficient (*μ*), refraction coefficient (2*πδ*_{n}/*λ*) and *δ*_{n}/*β*-ratio calculated with the XOP software [17] are given in table 1. In the test object, no empirical relation can be found. In this case, we let *m*(**r**) vary as the linear interpolation between the theoretical values (figure 2*a*). In bone, on the other hand, we can find an empirical relationship for *m*(**r**). We assume a simplified composition of bone, where it consists of a fixed weight fraction of collagen and a variable balance of embedding material and hydroxyapatite (the mineral constituent of bone), from only embedding material in soft tissue to only mineral in fully mineralized bone. For the calculation, we used for hydroxyapatite the chemical formula Ca_{5}(PO_{4})_{3}(OH) and a density of 3.16 g cm^{−3}, for the embedding material (poly(methyl methacrylate) (PMMA)) C_{5}O_{2}H_{8} and 1.17 g cm^{−3}, and for collagen C_{2}H_{5}NOC_{5}H_{9}NOC_{5}H_{10}NO_{2} and 1.41 g cm^{−3}. Under this assumption, we find that the relationship between attenuation and refractive index follows a power law (figure 2*b*).

### (c) Imaging

The two samples were imaged at the ID19 beamline of the European Synchrotron Radiation Facility (ESRF), Grenoble, France. In both cases, a double Si crystal monochromator was used to achieve a high degree of monochromaticity. A FreLoN detector [18], consisting of a scintillating screen for X-ray to visible light conversion, conventional visible light microscope optics and a 2048×2048 pixels CCD, was used for imaging. The test object was imaged at 0.7 μm pixel size with an X-ray energy of 22.5 keV, using four propagation distances (*D*=[2, 10, 20, 45] mm) and 1500 projections over a 180^{°} rotation. A bone sample from a rat osteotomy model [19] was imaged at 5 μm pixel size with an X-ray energy of 27 keV, three propagation distances (*D*=[60, 300, 999] mm) and 2000 projections over a 360^{°} rotation.

## 3. Results

The test objects were reconstructed using the mixed approach [5], the mixed approach with the homogeneous object prior [14], Paganin's method for homogeneous objects [4] and the method proposed above. The *δ*/*β*-ratio was set to 367, corresponding to Al for the homogeneous object approaches. Reconstructed slices are shown in figure 3. The mixed-approach reconstruction (figure 3*a*) shows a very sharp reconstruction, but is contaminated with characteristic low-frequency noise. The homogeneous object approaches (figure 3*b*,*c*) clearly improve the low-frequency behaviour, but leave fringes around zones where the *δ*/*β*-ratio is selected too low (i.e. the plastics). With the proposed method, this deficiency is alleviated and strong area contrast is visible between the different materials.

Quantitative evaluation of the reconstructions is presented in table 2. Normalized error (NE) and relative standard deviation (RSD) were calculated as
3.1where *l*_{t} is the theoretical quantity from table 1, *l*_{m} the measured mean value and *s*_{m} the standard deviation in the corresponding material. Note that the phase reconstructions do not achieve as good an accuracy as the attenuation, but the proposed method reaches an order of magnitude improvement in precision. Resolution was measured by fitting an error function to a line profile on the Al–air interface, and then calculating the corresponding Gaussian full width at half-maximum (FWHM) based on the error-function-fitting parameters. All phase reconstructions add substantial blurring. Note however that no deconvolution of the detector response was included in the reconstructions, which can be easily added owing to the high signal-to-noise ratio.

Another aspect to note is that while Paganin's algorithm yields a better reconstruction on average, the mixed approach with the homogeneous object prior achieves better accuracy, precision and resolution close to the selected *δ*/*β*-ratio. Finally, with the proposed prior both accuracy and precision are improved compared with the homogeneous object approaches, apparently at the cost of some resolution. The loss of spatial resolution could be owing to the implementation of the forward projection. It should be noted that the resolution with the proposed prior is the same across the whole image, whereas in the homogeneous object approaches, it is degraded as we move away from the chosen *δ*/*β*-ratio. For example, setting the *δ*_{n}/*β*-ratio to correspond to PET in the homogeneous object approaches, then measuring resolution on the Al–air interface yields an FWHM of 46.7 μm with the mixed approach and 47.5 μm with Paganin's method. Reconstructed slices of the bone sample and corresponding histograms using Paganin's method with *δ*_{n}/*β*=516, corresponding to fully mineralized bone (*a*,*c*), and the proposed method (*b*,*d*) is shown in figure 4. Also here, the homogeneous object approach yields fringes on the edges of materials where the selected *δ*_{n}/*β*-ratio is too low, as in the embedding material. The proposed method, on the other hand, correctly reconstructs both soft and hard tissues. This can be seen on the density histograms of the shown slices (figure 4*c*,*d*). The local mass density can be estimated from the reconstructed refractive index as [20]
3.2with *λ* in angstroms. Qualitatively, it can be seen that the peak corresponding to the embedding material (PMMA) is shifted to higher density values with the proposed method compared with the homogeneous approach, and a peak is visible corresponding to the bone, which is not visible in the homogeneous object reconstruction. Quantitatively, we find a mean density in the PMMA of 1.14 g cm^{−3}, compared with the density given by the supplier of 1.19 g cm^{−3}; hence a relative error of −4.20%. Furthermore, in the bone (mainly mineralized callus visible) we find an average density of 1.43 g cm^{−3}, which is reasonable for the partially mineralized callus tissue visible in the reconstructed slice.

## 4. Conclusion

We presented a new prior term for in-line phase tomography. It extends the method for multi-material objects presented in Langer *et al.* [3] by allowing the ratio between refractive index decrement and attenuation index to vary in the object, thus avoiding the segmentation step. We showed that for certain classes of samples, for example bone, empirical relationships can be found between the attenuation and refractive indices. We evaluated the method on a test object and in the imaging of bone samples for bone-healing studies, and found that the proposed method improved reconstructions both qualitatively and quantitatively.

## Footnotes

One contribution of 16 to a Discussion Meeting Issue ‘Taking X-ray phase contrast imaging into mainstream applications’ and its satellite workshop ‘Real and reciprocal space X-ray imaging’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.