Evaluation of phase retrieval approaches in magnified X-ray phase nano computerized tomography applied to bone tissue

: X-ray phase contrast imaging offers higher sensitivity compared to conventional X-ray attenuation imaging and can be simply implemented by propagation when using a partially coherent synchrotron beam. We address the phase retrieval in in-line phase nano-CT using multiple propagation distances. We derive a method which extends Paganin’s single distance method and compare it to the contrast transfer function (CTF) approach in the case of a homogeneous object. The methods are applied to phase nano-CT data acquired at the voxel size of 30 nm (ID16A, ESRF, Grenoble, France). Our results show a gain in image quality in terms of the signal-to-noise ratio and spatial resolution when using four distances instead of one. The extended Paganin’s method followed by an iterative refinement step provides the best reconstructions while the homogeneous CTF method delivers quasi comparable results for our data, even without refinement step.


Introduction
Imaging the three-dimensional (3D) organization of bone structures has been a pioneering application in micro-computerized tomography (micro-CT) [1]. Conventional X-ray micro-CT based on laboratory X-ray sources is now routinely used to assess 3D bone micro-architecture. When coupled to synchrotron radiation (SR), with the possibility to use monochromatic or quasi-monochromatic sources, image quality is improved since beam hardening artifacts are avoided, which permits to obtain a quantitative map of the linear attenuation coefficient [2]. Over the past two decades, a number of third-generation SR sources have been actively developed and most of them have designed SR micro-CT setups.
An additional advantage of SR micro-CT is the facility to implement X-ray phase imaging with the high degree of spatial coherence. This technique is attracting more attention with the progresses in SR sources, detectors and computing facilities. Since its sensitivity can be more than three-orders of magnitude higher than attenuation-based imaging, it enables imaging of weakly absorbing samples such as soft tissue and of materials with small changes in attenuation [3,4]. Phase micro-CT has been used in many biomedical areas and has provided valuable images of small animals or biological samples at micrometric scale. For example, it was used to image small animals [5], brain tissue [6] and lung microstructure [7].
Phase contrast can be obtained experimentally by various techniques either on synchrotron radiation or laboratory sources (for a review see [8,9]). Here, we will focus on propagation-based phase contrast imaging [10] as it provides the highest spatial resolution. This technique can be used with one or several propagation distances. In previous works, the optimization of experimental parameters such as the source size and the detector resolution has been studied using a single propagation distance [11,12].
Phase retrieval can be expressed as solving a non-linear inverse problem. The corresponding direct problem of image formation in X-ray propagation based phase contrast imaging can be modeled by Fourier optics and Fresnel diffraction considering the weak interaction of X-rays with matter inside the object [13]. Under these circumstances, the intensity recorded by the detector can be expressed as the squared magnitude of the Fresnel transform of the transmittance function of the object.
While the phase retrieval inverse problem is intrinsically non-linear, various approaches have been proposed to linearize it. The first one is the transport of intensity equation (TIE) which linearizes the intensity by Taylor expansion of the propagator to the first order with respect to the propagation distance. The resulting partial derivative in the propagation direction is then approximated by a finite difference, which is only valid for small propagation distances and low spatial frequencies [14]. The other one is the contrast transfer function (CTF) which linearizes the intensity by expanding the transmittance function to the first order with respect to the amplitude and phase. This yields a linear relationship which is valid for weak absorption and slowly varying phase [15]. The TIE generally needs to take different images at two sufficiently close planes from the object, while the CTF is still limited to low absorbing objects. A mixed method which extends the validity to both long distances and strongly absorbing (but slowly varying) objects has been proposed with the asymptotical combination of CTF and TIE [16]. The different algorithms have been quantitatively compared in previous work on numerical and experimental phantoms [17]. Generally, phase retrieval in the Fresnel region suffers from sensitivity to noise in the low spatial frequency range yielding cloud-like artefacts due to a slowly varying background. This has prompted the introduction of a priori knowledge on the object in order to regularize the problem [18].
By adding an assumption on the homogeneity of the object directly in the TIE model, Paganin's method has been derived [19]. This method has proven popular since it only requires a single phase contrast image. It is essentially based on a Fourier low-pass filter which can be applied directly on the recorded projection image as a correction for the propagation effect. Recently, the effectiveness of the method for increasing signal-to-noise ratio while preserving spatial resolution was analyzed by modeling the properties of the noise with propagation [20]. It is essentially based on a Fourier low-pass filter which can be applied directly on the recorded projection image as a correction for the propagation effect. However, although it has been used in a number of applications, it requires setting the n δ β ratio of the imaged material, which is the ratio between the real and imaginary part of the refractive index of the object, a material and energy dependent constant [19]. For a given chemical composition, this ratio can be calculated from tabulated values, e.g. by the X-ray optics software XOP [21]. Besides, this method has also been extended to two-materials objects by modifying the filter [22] or by using 3D filtering based corrections in the object space [23].
In this paper, we focus on imaging bone tissue at the cellular scale by using a new magnified phase nano-CT set-up installed on the beamline ID16A at the ESRF. The beamline is capable of focusing the beam to the nanoscale while keeping a high photon flux. This setup permits to record one or several data sets at different propagation distances. However, acquiring data sets at different distances generates long scan time (about 1h per distance), and possibly more radiation-induced damage to the sample. Thus, we propose to investigate the gain of using several propagation distances instead of one in terms of quality of the reconstructed images, on one hand, and to evaluate phase retrieval algorithms for this type of relatively dense samples, on the other hand. To this aim, we first extend Paganin's single distance formula to multiple distance acquisitions. Then, for comparability, we describe the multi-distance CTF method for homogeneous objects with a known n δ β ratio. Further, we examine the relationships between these two methods by an analysis in the Fourier domain. The solution obtained with one of these two methods is then refined iteratively by a nonlinear conjugate gradient (NLCG) method. The two phase retrieval methods are applied to data sets of bone tissue samples acquired at the voxel size of 30 nm. In order to compare these different phase retrieval approaches, we use qualitative visual assessment and quantitative evaluation based on the SNR (Signal-to-Noise Ratio) and an estimate of spatial resolution. From a visual point of view, we show that the reconstructed images from the two multiple distances methods are sharper than those reconstructed from one propagation distance only. This is corroborated by the quantitative study, which showed an increase in estimated spatial resolution when using images recorded at several propagation distances.

Single distance phase retrieval methods
We briefly recall the notations of the forward problem and two solutions to the inverse phase retrieval problem using single distance data. An overview of single distance phase retrieval methods can be found in [24].

Model of image formation
The interaction between the 3D object and X-rays at a given wavelength λ can be modeled by the complex refractive index distribution [16], where n δ is the refractive index decrement, β the absorption index and ( ) , , x y z the 3D spatial coordinates.
If the propagation direction is z , and the illumination considered flat, the wave after passing through the object can be expressed by the transmittance function ( ) T x which can be written as: where ( ) a x is the amplitude, ( ) x the Fresnel propagator and k the distance index.
X-ray phase CT reconstructs the refractive index decrement n δ which is about threeorders of magnitude higher than the absorption index β , thus the technique offers higher sensitivity than conventional X-ray CT. The 3D phase CT images are obtained from phase maps retrieved at each projection angle as inputs of reconstruction instead of absorption images [15]. Therefore, it is necessary to estimate the phase shift ( ) ϕ x from recorded phase contrast images before tomographic reconstruction.

Paganin's method
The single distance Paganin's method can be derived from the TIE model by using an additional homogeneity constraint on the sample [19]. This is introduced directly in the contrast model. In the case of homogeneous composition, the phase map is proportional to the absorption map. The proportionality constant is given by the ratio between the real and imaginary part of the refractive index n δ β and is supposed to be known a priori. By using this constraint, the following retrieval formula is obtained: where Inc I is the incident intensity, and F (resp. 1 F − ) denote the 2D direct (resp. inverse) Fourier transform operator.

Homogeneous CTF method
The CTF method relies on a linearization of the direct problem, based on the first order Taylor expansion of the transmittance function (Eq. (2)), which can be written as: The CTF model at one distance can be expressed as [15]: Similarly to Paganin's method, the CTF method can be combined with the homogeneity assumption [25]. With the assumption of a known n δ β , the absorption is proportional to the phase: Substituting Eq. (7) to (6), the solution of the CTF method in the Fourier domain can be expressed as: Because of zero crossings in the denominator of Eq. (8), some frequencies of ϕ cannot be recovered. This generates artifacts in retrieved phase maps. This problem can be overcome by combining intensity information ( ) k D I  f at several propagation distances [26].

Multi-distance phase retrieval methods
In this section, we derive multiple distance versions of the previous methods: the extended Paganin's method and the homogeneous CTF method. Then we discuss the relationships between the two methods.

Multi-distance Paganin's method
To make the notations more concise, we introduce the three following variables: According to Eq. (9), Eq. (4) becomes: When considering different propagation distances, the phase retrieval problem consists in . This problem can be solved using an optimization method. We consider a regularized least squares minimization according to: The first term is the data fidelity term and the second is a Tikhonov regularization term, where α is a regularization parameter. The solution can be written as: The phase map estimated from acquisitions at different distances is finally obtained as:

Multi-distance homogeneous CTF method
Similarly, we introduce a new variable to express the CTF filter, The forward model of the homogeneous CTF can then be expressed as: Then the regularized least-squares solution in the case of homogeneous CTF and multidistance acquisition can be written as:

Relationships between the two approaches
Following Eq. (13) and (16), both methods can be interpreted as a linear combination of filtering of the recorded intensities in the Fourier domain. It may be noted that the frequency filters involved in both methods are radial. Let us first examine the single distance situation. The multi-distance expressions (13) and (16) The expression of the phase as a function of the filters differs in the two methods due to the natural logarithm term in the Paganin's expression.
If we consider the weak absorption case where , by the Taylor expansion at the first order, we get: Then calculating the Fourier transform of ( ) Substituting it into Eq. (4) yields, Finally, we obtain the solution of ( ) Under the assumptions of weak absorption and low spatial frequencies, Paganin's method and the homogeneous CTF method are identical. However, for higher spatial frequencies, the two filters essentially differ at locations corresponding to the zero crossings of the direct CTF model.
Let us now consider the case where information from several distances is used.
In this case, we can summarize the homogeneous CTF method as: (22) and the extended Paganin's method as: As previously, the homogeneous CTF method and the extended Paganin's method at multiple distances will be equivalent in the low frequency range when the absorption is weak, but will differ for higher frequencies.

Data acquisition
Data were acquired at the beamline ID16A, ESRF (Grenoble, France). This beamline combines coherent imaging techniques and X-ray fluorescence microscopy at a spatial resolution down to about 30 nm, enabling quantitative 3D characterization of morphology and elemental composition of specimens in their native state [27]. ID16A provides a highbrilliance beam at two specific energies (17 keV or 33.6 keV). The incoming parallel X-ray beam is monochromatized and focused into a focal spot by multilayer coated curved reflective optics (Kirkpatrick-Baez mirrors). After passing through the sample, the X-ray beam is converted into visible light by a scintillator and recorded by a lens-coupled FReLoN (Fast Readout Low Noise) camera developed by the instrument support group of ESRF with a binned CCD size of 2048 × 2048 pixels.
X-ray phase nano-CT data acquisition consisted in recording sets of angular projections at different sample positions between the focus and the detector. After being focused, the beam turns from parallel to divergent so that the images recorded at different distances contain not only phase contrast but also geometric magnifications. In practice, the position of the detector is fixed while the sample is moved downstream of the focus [28,29]. The magnification was chosen to obtain an object voxel size of 30 nm.
The samples were extracted from transverse cross-sections cut from mid-diaphysis in femurs from female cadavers (50-95 years old). Then, small cortical bone samples (0.4 × 0.4 × 3 mm 3 ) were prepared using a high precision low-speed circular saw. Each sample was positioned on the vacuum sample stage and rotated around a vertical rotation axis. In total, 2000 projections were acquired over a range of 180° at 4 different sample-to-detector distances. The acquisition time for a complete data set (i.e. with 4 propagation distances) was approximately 4 hours.

Data processing
First the images are so-called flat-field corrected by subtraction of the dark current image and division with the empty beam image recorded without sample. Next, the images at different propagation distances must be registered for each rotation angle. This is done after correction of the magnification k M that depends on the focus-to-sample distance k Z and the focus-todetector distance d Z according to the expression (Fig. 1): The images were rescaled at the highest magnification. Phase maps were estimated by using different retrieval methods. We used the equivalent propagation distance that can be expressed as: Table 1 lists some specific experimental parameters such as 4 different focus-to-sample distances, the focus-to-detector distance, the wavelength of the X-ray, the voxel size and Fresnel number k F which is calculated by: where A is the is the characteristic object size, λ the wavelength and k D the propagation distance. The Fresnel numbers are provided both for the osteocyte lacunae and the canaliculi by setting A as different typical sizes, 10 µm and 200 nm respectively.
The n δ β ratio was set to 645 according to the complex refractive index of cortical bone at 33.6 keV, which was calculated from tabulated values by the software XOP. To improve image quality we also considered refining the estimated phase maps by using a non-linear conjugate gradient method (NLCG) based on the direct nonlinear image formation model.
Finally, the extracted phase maps were used as the input to tomographic reconstruction to get the 3D structure of the sample. As there are a sufficient number of angular projections, the reconstruction can be performed by the Filtered Back-Projection (FBP) algorithm. We used the ESRF implementation PyHST2 for tomographic reconstruction [30].

Quantitative evaluation of the 3D reconstructions
To get a quantitative assessment of the different reconstruction methods, we computed the Signal-to-Noise Ratio (SNR) and an estimate of the spatial resolution. where w is the FWHM and 0.1 a = the position of the cut-off frequency. The estimation of the SNR and the spatial resolution was repeated at different locations to obtain average and standard deviation values. To calculate the spatial resolution, edges are selected in the middle slices of different lacunae from the reconstructed volumes.

Results
Figure 2(a) displays the filters of the CTF and Paganin's methods using a single distance given by Eq. (8) and (4); (c) (resp. (e)) shows the filters of the homogeneous CTF method given by Eq. (22) (resp. the extended Paganin's method given by Eq. (23)) using 1 to 4 different propagation distances omitting the image term; (b), (d) and (f) are zooms of the filters shown in (a), (c) and (e), respectively. We can see that the CTF method and Paganin's method behave similarly in the low frequency range when using a single or multiple distances. However, the filter corresponding to the single distance CTF method clearly contains numerous jump points with large amplitudes in the high frequency range, as opposed to Paganin's method which is a simple low-pass filter (Fig. 2(a)). When considering multidistance acquisitions, the amplitudes of the CTF filters are bound in the high frequency range, which shows that the use of multiple distances correctly handles the zero crossings in the direct CTF model. The homogeneous CTF filters again overlap with the extended Paganin's method filters in the low frequency range. However, while the filters in Paganin's method are essentially low-pass filters, the filters involved in the homogeneous CTF method (Fig. 2(d)) enable to retrieve the high frequency information during phase retrieval.
Minimum Intensity Projections (MIPs) of the 3D reconstructed volumes for a human cortical bone sample are shown in Figs. 3 and 4, respectively for the homogeneous CTF and the extended Paganin's method. In our case, the MIPs are calculated by selecting the minimum intensity value along the Z-axis for each XY coordinate. The MIPs permit to see lacunae (ellipsoidal-like black structures), and also to enhance canaliculi which are the very small channels connecting bone lacunae. Figures 3(a)-3(c) illustrate MIPs retrieved without refinement ("no it.") using 1, 2 and 4 distances respectively; Figs. 3(d)-3(f) display MIPs retrieved with 10 iterations refinement using 1, 2 and 4 distances respectively. Figure 4 shows the same MIPs, when the phase retrieval is performed using the extended Paganin's method.
First the impact of the iterative NLCG refinement is studied. Generally, the slices reconstructed and refined with 10 iterations look sharper and better contrasted than ones with no refinement. However, the gain in image quality is somehow different when considering the CTF or Paganin's method. For the CTF method, the non-linear refinement does not improve image quality very significantly, whatever the number of distances used (1, 2, or 4), as shown in Fig. 3. However, the non-linear refinement clearly improves the image quality for the extended Paganin's method. In Figs. 4(a)-4(c) the MIPs retrieved by the extended Paganin's method without iterations look blurred, which is in agreement with the low pass nature of the corresponding filters. It turns sharper with more high frequency information after the NLCG refinement as shown in Figs. 4(d)-4(f). Now, we consider the impact of using multiple distances rather than one on the reconstructed image quality. For the homogeneous CTF method, if we compare the single distance case (Figs. 3(a) and 3(d)) with the 4 distances case (Figs. 3(c) and 3(f)), it can be seen that the image using one propagation distance is more blurred and that the sharpness is improved when using 4 propagation distances, regardless the use of iterative refinement. For the extended Paganin's method without refinement, we observe that the difference between reconstructions using a single distance and multiple distances (see Figs. 4(a) and 4(c)) is barely visible. However after performing NLCG refinement, a better contrast of small channel structures as well as a decrease of noise can be observed when using multiple distances instead of one (see Figs. 4(d) and 4(f)). Additional slices corresponding to the results obtained in Figs Figures 4(a), 4(d), 4(c) and 4(f)). Table 2 presents the estimated SNR and spatial resolution in the reconstructed images for both the homogeneous CTF and extended Paganin's methods. These quantities are also displayed in Figs. 5(a) and 5(b) respectively. For the homogeneous CTF method, both the SNR and the spatial resolution get improved with the number of distances. In the single distance case, the iterative refinement improves the spatial resolution but decreases the SNR. For the extended Paganin's method without iterative refinement, the use of multiple distances does not yield significant differences. However, if we consider the non-linear refinement, the use of multiple distances obviously improves the SNR and to some extent the spatial resolution. Besides, the best results are obtained by the extended Paganin's and homogeneous CTF methods using 4 distances with 10 iterations.

Discussion
In this study, we compared and evaluated two phase retrieval methods for homogeneous objects based on the CTF and TIE approaches. We studied an extension of the single distance Paganin's method to multi-distance acquisitions and provided a similar formula for the multidistance CTF method using a homogeneity assumption on the sample. Moreover, assuming weak absorption, we derived relationships between the extended Paganin's and the homogeneous CTF methods. Considering the case of a single propagation distance, the homogeneous CTF and Paganin's filters are almost identical when restricted to low spatial frequencies and low absorption. Nevertheless, the filters differ at higher spatial frequency and particularly around the zero crossings of the direct CTF model where the homogeneous CTF filter reaches large values. When it comes to multi-distance acquisitions, the filters of both methods are still similar in the low frequency range. However, the amplitudes of the CTF filters are bound and different from zero which permits to retrieve the high frequency information while Paganin's method always acts as a low-pass filter.
These findings were confirmed by reconstructions performed on an experimental bone data set with a voxel size of 30 nm. From those reconstructions, we observed that the structures such as lacunae and canaliculi are similarly reconstructed in the images from the extended Paganin's and homogeneous CTF methods using 4 distances. The quantitative criteria, such as the SNR and the spatial resolution corroborate the visual observation. In this paper, we calculate the standard SNR which is used in the very large majority of studies and comparable to the literature. Other quality measures such as the recently proposed spatial SNR could be investigated to evaluate the quality of X-ray images [31].  The phase retrieval and tomographic reconstruction were performed sequentially. This has many practical advantages as only a large number of 2D problems have to be solved. Recently some methods have been proposed to combine the two steps together [32][33][34]. However, 3D approaches require extensive computer resources both in terms of computing time and memory. Therefore they are not practical for standard use on large data sets today, but they will be increasingly important in the future.
In previous works, Paganin's method was compared with the modified Bronnikov algorithm and was shown to provide better reconstructions on low-attenuating samples [35], but only one propagation distance was considered. Similarly, a good agreement was shown between the TIE and the CTF method in the low-spatial frequency range [25]. However, the work was limited to one propagation distance and not applied to CT reconstruction. Burvall et al. detailed the processing involved in seven phase retrieval methods found in the literature and made a comprehensive analysis of these methods [24]. The Paganin's and CTF methods considered here corresponding to their "single material" and "Fourier-Born" methods have the expected behavior. But again their analysis was performed for a single propagation distance and results were provided on 2D phase maps but not on 3D reconstructed images. Besides, Krenkel et al. proposed a Holo-TIE algorithm to address the phase retrieval in hard x-ray in-line holography and extended it from single-distance recordings to multiple-distance cases for achieving single-cell imaging [36].
In our experimental conditions, the best reconstructed images (i.e. with a sharper display of edge information and a better detection of small structures) were obtained after applying a non-linear conjugate gradient method to refine the results of the filtering based phase retrieval methods. Actually, iterative near-field methods have been recently developed a lot in phase contrast imaging. However, they are mostly based on a single distance and generally use iterative algorithms, such as Fienup's error reduction (ER) [37] or the Gerchberg-Saxton (GS) algorithms [38]. In this work, we used a NLCG algorithm, taking as initialization, the phase maps retrieved from the extended Paganin's or homogeneous CTF method. Our results indicate that this refinement critically improves the spatial resolution if the extended Paganin's method is used, while the homogeneous CTF alone is already close to the optimal result. Both initial linear estimates converge in the end to similar solutions. However, it is worth underlining that the improvement obtained with the non-linear refinement is marginal in the case of the homogeneous CTF method which is probably due to the fact that our samples approximately satisfy the weak absorption assumption. Thus, it would also be possible to use the multi-distance homogeneous CTF method without iterative refinement if computing time is an issue.

Conclusion
Phase retrieval which is a non-linear ill-posed inverse problem is a necessary step in the reconstruction of propagation-based SR X-ray phase nano-CT. The choice of the phase retrieval method and the best experimental conditions remained open questions in the multiple distance case. This has motivated our detailed study of two phase retrieval approaches in the case of multiple propagation distances from a theoretical and experimental point of view. Paganin's method using a single distance is one popular phase retrieval method due to its simple implementation. In this paper, we extended Paganin's method to multidistance acquisitions, and compared it to a homogeneous version of the CTF method for the reconstruction of bone samples at nanoscale spatial resolution. Our results showed that using 4 propagation distances generally improves image quality except in the case of the Paganin's method without non-linear refinement. When iterative refinement was used, the multidistance case provided the sharpest reconstruction with a higher spatial resolution and flatter background. In this case, the homogeneous CTF and the extended Paganin's method give comparable results. The advantage of using several propagation distances is that it provides a better coverage of the Fourier domain of the phase projections.
Our perspectives concern the quantitative analysis of the bone ultra-structure on which there are still very few data due to the limited 3D imaging modalities available for this purpose. In this respect, it is important to optimize image quality since the latter may have a considerable impact on the further quantification. While it is experimentally more demanding to use four propagation distances than one, our results showed that a clear gain is obtained in image quality. Future works will focus on the quantitative analysis of the data, from which we expect to obtain quantitative parameters, of the lacuno-canalicular network such as porosity and connectivity but also information about the extra cellular matrix.