1. Introduction
Hard X-ray imaging with a high spatial resolution is nowadays a powerful tool to investigate specimens in 2D or 3D in a non-destructive manner. For an object illuminated by partially coherent light sources, a simple and effective technique known as propagation-based phase contrast has a particular interest because of its simple imaging set-up. Additional optical elements are not required and the phase contrast images can be recorded by simply letting the X-ray beam propagate in free space after interaction with the sample [1] [2] (Figure 1).
Compared with attenuation-based imaging techniques, the main interest in X-ray phase imaging is the pos- sibility to study objects with either negligible absorption or dense objects with small density variations. More- over, in the hard X-ray region, the phase shift for low-Z elements improves the sensitivity with three orders of magnitude [3] , which makes this imaging modality attractive for biomedical imaging of soft tissues. The phase-contrast images do not yield directly the phase information and requires additional experimental set-ups, models and data analysis algorithm. The Fresnel diffraction intensity patterns set an ill-posed non-linear inverse problem. Phase retrieval and tomography can be coupled by a two-step process: first, the phase information is retrieved for all the projections, and secondly the three-dimensional tomographic reconstruction is performed on the retrieved phase images obtained for each angle of view (see Figure 2).
The most common algorithms for the phase retrieval problem for short propagation distances rely on the linearization of the Fresnel diffraction relationship [4] - [11] valid under restrictive assumptions. As far as phase tomography is concerned, several methods have been studied extensively. Langer et al. [9] have proposed to introduce the prior on the retrieved phase that the phase and the absorption are proportional. A single-distance phase retrieval method for a homogeneous object for a given ratio of the imaginary to the real part of the refractive index has been developed by Paganin [4] . This type of prior is valid for multi-material objects comprised of several homogeneous objects [7] [10] . A new inversion method where a prior phase estimate at each projection angle is obtained from a measured absorption index map evaluated with the intensity measured for a propagation distance 
m is described in [11] . This prior is introduced in the low-frequency range only. This method is an extension of the previous linear algorithm [8] including a Tikhonov regularization term to the tomographic case. Compressive sensing approaches have been also studied recently but they are restricted to small scale problems [12] [13] . The limitations of the approaches based on the linearization of the direct problem can be overcome by other methods which take into account the non-linearity of the phase problem. The phase retrieval problem is an inverse ill-posed problem, therefore regularization methods are necessary to obtain stable solutions less sensitive to noise. The non-linear contributions in the image contrast formation are non- negligible since large propagation distances and high spatial resolution are required. Consequently, the non- linearity of the phase-intensity relationship is a crucial aspect. New algorithms which take into account the non- linearity of the inverse problem for the radiographic case have been proposed recently [14] - [17] . These non- linear approaches are very promising and lead to a large decrease of the reconstruction errors.
![]()
Figure 1. Experimental set-up for propagation-based technique or in-line phase contrast imag- ing technique for a parallel X-ray beam. The incident field is assumed to have a degree of par- tial coherence and passes through a probed sample. Phase contrast images will be registered on the CCD-based detector for different distances 
in the Fresnel field.
In this paper, we consider the case of in-line phase tomography using different propagation distances. We extend to the tomographic case a modified non-linear algorithm proposed in [16] which is based on the Fréchet derivative of the intensity operator. As detailed in [16] the inverse ill-posed problem of the phase recovery is stabilized by a Tikhonov type regularization term with the square of the gradient phase term. In the following,
the Landweber iterative algorithm is modified by replacing this term 
with the phase term
. In this
paper, we first summarize this new multi-image non-linear ![]()
scheme, and then detail the results obtained on simulated images and for a tomographic reconstruction on a real multi-material 3D object.
2. Non-Linear Phase Retrieval Approach
2.1. Image Formation―The Direct Problem
For a parallel, partially coherent, monochromatic X-ray wave with wavelength![]()
, an object is characterized by the 3D complex refractive index [18] :
![]()
(1)
with ![]()
the spatial coordinates. The diffraction within the object is neglected due to the weak interaction of X-rays with matter and the exit wave can be described by the complex transmission function ![]()
at each projection angle ![]()
[18] :
![]()
(2)
where ![]()
is the spatial coordinates in the plane perpendicular to the propagation direction of the X- rays![]()
. The absorption ![]()
and phase shift ![]()
induced by the object for each projection angle ![]()
are considered as the projections of the absorption index ![]()
and of the refractive index decrement ![]()
respec- tively. The linear integral relationships used for the tomographic reconstruction, are given by
![]()
and![]()
.
The intensity distribution ![]()
for each distance ![]()
(where![]()
) and angle ![]()
can be related to ![]()
by the following expression:
![]()
(3)
where ![]()
denotes the coordinates in a plane perpendicular to the propagation direction ![]()
and ![]()
the 2D con-
volution of the complex transmission function ![]()
with the corresponding Fresnel propagator ![]()
at distance ![]()
and wavelength![]()
.
2.2. Non-Linear Inverse Problem-Phase Retrieval
As detailed in [16] [17] [19] , assuming that the phase ![]()
is defined on![]()
, an open subset of![]()
, the intensity operator ![]()
(Equation (3)) can be considered as a continuous and non-linear
function of the phase![]()
, which is a Fréchet differentiable in its domain [20] . Classical regularization algorithms in a Hilbert space can be used to study the non-linear phase retrieval problem [20] . The classical Landweber algorithm is based on the minimization of the regularization functional with its gradient [20] . The
previously proposed Landweber iterative approach [16] is based on the ![]()
norm of the phase gradient ![]()
as regularization term. In this study, better convergence results were obtained by replacing this term with the phase term![]()
.
The aim of this non-linear approach ![]()
is to minimize for each projection angle ![]()
the following cost functional:
![]()
(4)
where ![]()
is the measured noisy intensity at distance ![]()
(see Figure 1) for the projection angle ![]()
and ![]()
is the regularization parameter. This ill-posed problem is stabilized by a Tikhonov type regularization term with the square of the phase term ![]()
[17] . The regularization parameter ![]()
is chosen by trial-and-error for one projection angle and then fixed for all the projection angles.
The minimizer of the cost functional is calculated iteratively with a non-linear Landweber type scheme [21] :
![]()
(5)
The classical Landweber method is modified with a variable step ![]()
chosen using a dichotomy strategy. The algorithm in this form is a simplified version of the iterative Gauss-Newton method considered in [22] [23] . It was shown in [16] that at the point ![]()
the Fréchet derivative of the operator ![]()
is the linear operator defined as:
![]()
(6)
that can be obtained using the explicit calculation [16] :
![]()
(7)
where ![]()
denotes the convolution operator. Moreover, by using the standard definition of the scalar product in ![]()
spaces the adjoint operator ![]()
is given by [16] :
![]()
(8)
Thanks to the analytical expressions of the Fréchet derivative and of its adjoint, it was possible to decrease the computation time and to obtain a better convergence in the radiographic case. In order to take into account the intensity maps, better results were obtained when the propagation distances were used in a random way and not in a successive way during the iterations.
2.3. Initialization and Stopping Rules
It was shown that the non-linear algorithm improves the solution obtained with a linear algorithm in the radio- graphic case for simulated data [16] [17] . Yet, an initialization obtained with the mixed linear scheme is ne- cessary to obtain convergence of the non-linear method. In this work, for each projection angle, the ![]()
algo- rithm was initialized either with the phase retrieved without any a priori knowledge or with a fixed a priori value of the ratio ![]()
[8] [9] . In a second step, the tomographic reconstruction was performed from the whole set of phase maps![]()
, with a standard filter back projection (FBP), implemented at ESRF (PyHst) [24] .
For a projection angle![]()
, the iterations during the NL algorithm are terminated when:
![]()
(9)
where ![]()
is a parameter that was set at 0.01 by trial-and-error.
It is well known that the regularization parameter plays a crucial role in the convergence of the iterative regularization methods, therefore it has to be chosen carefully. In all the studies of this work, large and small values of the parameter leading to poor reconstruction results are first chosen. Then, the optimal value of the re- gularization parameter is gradually refined by trial-and-error with a decreasing interval. For a well-chosen para- meter, the errors on the intensity maps for all the three propagation distances are decreased. If the convergence is not achieved for all the propagation distances, the value of the regularization parameter is refined till the convergence is achieved.
2.4. Simulations and Data Acquisition
The new non-linear inversion method has been tested on simulated images and on experimental data for a multi- material object.
2.4.1. Simulation of the Image Formation
Two phantoms were defined in a deterministic fashion [25] , one for the absorption coefficient and one for the refractive index decrement. Figure 3(a) displays the 3D Shepp-Logan [24] , consisting of a series of ellipsoids on which the projections are based. Theoretical values for the absorption coefficient ![]()
and for the refractive index ![]()
of different materials at 24 keV were used (![]()
![]()
) in different regions (Table 1). Analytical projections were calculated in a parallel beam geometry with ![]()
pixels and the two resulting data sets were combined to form a complex representation of the wave exiting the object using Equation (2). Propa- gation in free-space was simulated using Equation (3) for three distances D = [0.035; 0.072; 0.222] m. The original phase contrast intensity maps were digitized to ![]()
pixels and corrupted with additive Gaussian
![]()
Table 1. Values of the absorption coefficient and refractive index at 24 keV for the materials used in the 3D phantom.
![]()
(a) (b)
Figure 3. (a) Ideal phase to be retrieved and (b) absorption image with PPSNR = 24 dB for strongly varying phase.
white noise with various peak-to-peak signal to noise ratios (PPSNR) of 24 dB. The peak-to-peak signal to noise ratio is defined by:
![]()
(10)
where ![]()
is the maximum signal amplitude and ![]()
is the maximum noise amplitude.
In order to asses the performance of the new regularization scheme, two types of objects were considered for short propagation distances and weak absorption, one with a slowly varying phase and another with a strongly varying phase.
2.4.2. Experimental Images
The experimental set-up used is equivalent to the one for the standard propagation based technique described in [9] , at the beam line ID19 at the European Synchrotron Radiation Facility (ESRF). The Fresnel diffraction in- tensity patterns for ![]()
projection angles were recorded using a FRELON CCD camera with ![]()
pixels for the energy ![]()
keV at four short distances ![]()
mm. The field of view was ![]()
mm for the given pixel size 0.68 μm. The multi-material object used is composed of 125 μm Aluminium![]()
, 200 μm Polyethylene Terephthalate ![]()
mono-filaments, 20 μm of Alumina ![]()
wires and 28 μm Polypropylene ![]()
fibres. Phase retrieval with the mixed approach was applied without any prior on ![]()
[8] (initialization (A)) and with ![]()
[9] corresponding to aluminium (initialization (B)).
3. Results
3.1. Simulated Data
The efficiency of the proposed new regularization scheme was analysed by comparing the numerical results obtained with the NL method with the phases retrieved with the CTF, TIE and mixed linear approach in the radiographic case. The four methods were tested for weakly and strongly varying phases, and for noise-free and noisy data.
Since the ideal reconstruction image is available, direct comparisons can be made. The method will be quantitatively evaluated by measuring the normalized mean square error (NMSE) using the ![]()
norm:
![]()
(11)
where ![]()
is the phase recovered at iteration ![]()
and ![]()
the ideal phase to be recovered.
The NMSE (Equation (11)) for all the methods are presented in Table 2. For the strongly varying phase without noise, the non-linear approach gives the most accurate results. For the weakly varying phase for noise- free data, the TIE method gives the best solution. On the other hand, for noisy simulated data with PPSNR = 24 dB, TIE yields the worst reconstructions. As shown in Table 2, the errors on the phase have been significantly reduced with our non-linear algorithm using as starting point the mixed phase map solution.
The evolution of the NMSE as a function of the iterations number is displayed in Figure 4 for the various cases investigated. In these plots, one iteration corresponds to a random cycle through the intensity images
![]()
Table 2. NMSE (%) values for different algorithms and objects.
obtained for the three distances. The initialization of the NL algorithm for these four situations was given by the linear mixed solution. These curves show that the proposed algorithm has good convergence properties. Very few iterates are necessary to obtain an improved stationary point.
3.2. Experimental Data for Non-Linear Phase Tomography
The reconstructed projections for the angle of view ![]()
retrieved with the mixed algorithms in these two cases are displayed in Figure 5(a) and in Figure 5(b), respectively. The non-linear phase map obtained for the initialization map given in Figure 5(a) is shown in Figure 5(c). Figure 5(d) displays the phase retrieved with NL with the starting point given by the linear solution displayed in Figure 5(b). Starting from these images the FBP is applied and the 3D refractive index decrement ![]()
is reconstructed.
The tomographic central slices of the refractive index decrement, in the case of the mixed algorithm with a standard Tikhonov regularization without any a priori knowledge on the ratio![]()
, is displayed in Figure 6(a). The corresponding central slice obtained with the non-linear approach initialized with this linear phase solution is shown in Figure 6(c). In order to have a quantitative estimate of the reconstruction errors, the theoretical values for ![]()
and the values estimated with the linear algorithm or with the non-linear approach are summarized in Table 3. In Table 4 the relative standard deviations (RSD[%]) and the normalized errors (NE[%]) for the four component materials have been measured for all reconstructions approaches. The RSD and the NE values were measured using:
![]()
Figure 5. Projection images corresponding to the angle of view 120˚ obtained after the phase retrieval step with the mixed algorithm (a) without a priori information [8] and (b) with ![]()
(Al) [9] . The projections obtained using NL initialized with these mixed solutions are displayed in (c) and (d) respectively. Gray-scale windows in (a), (c) is ![]()
and in (b), (d).
![]()
Table 3. Theoretical and measured values with different algorithms ![]()
(cm−1).
![]()
Figure 6. Tomographic central slice reconstructed with the mixed algorithm (a) without a priori information [8] and (b) with a priori information ![]()
(corresponding to aluminium) [9] . Corresponding central slice obtained with the non-linear algorithm initialized with the linear solution (c) without a priori (ini- tialization displayed in (a)) and (d) with a priori information ![]()
(initialization displayed in (b)).
![]()
Table 4. Values for relative standard deviation (RSD) and normalized error (NE) obtained with different algorithms.
![]()
(12)
and
![]()
(13)
where SD represents the standard deviation, ![]()
the theoretical value to be obtained and ![]()
the measurements (given in Table 3). The theoretical values were obtained using the tabulated values in the XOP software [26] . If the a priori ratio is not included in the initialization algorithm, the proposed approach reduces the total NE by![]()
. The reconstructed refractive index decrement obtained with the NL algorithm is better estimated for all the components of the sample (Table 3). In the case where the exact value of the a priori ratio is not known, which is the case for biomedical samples, this result shows that the non-linear algorithm is an interesting extension of the mixed approach.
The tomographic central slice obtained using the mixed approach with a prior value of the ratio ![]()
corresponding to aluminium is displayed in Figure 6(b), the corresponding non-linear reconstruction using this linear initialization is shown in Figure 6(d).
Comparing this reconstruction (Figure 6(b)) with the one where the a priori was not introduced (Figure 6(a)) in the mixed method, it can be observed that low-frequency noise artefacts are alleviated. The non-linear algo- rithm is very efficient to improve the uniformity of the reconstructed image as we can see in Figure 6. The non- linear solution retrieved using as starting point the linear solution with ![]()
provides more accurate re- constructions (Figure 6(d)). In this case, the NE [%] value corresponding to Al is overestimated, but the NE [%] for Al2O3 is reduced with 33.5% (Table 4). The overall reconstruction error is also decreased. The efficiency of the non-linear scheme depends thus strongly on the initial linear method used, and the best results are obtained when no assumption is made on the ratio![]()
. The total values of the normalized errors (Table 4) have been improved by the non-linear algorithm. The proposed approach reduces the global error in the reconstructed materials compared to the two linear initialization solutions. For most materials, the lowest error values are obtained by the non-linear algorithm. Nevertheless, most materials are underestimated (minus sign of NE in Table 4) which can also be related to some imperfection of the detector not taken into account in this study.
4. Discussion and Conclusion
In this paper, we have considered a non-linear phase retrieval method for phase tomography. The method has been evaluated quantitatively on simulated images and from experimental data acquired at three different propagation distances on a synchrotron X-ray micro-CT set-up. The proposed NL algorithm is achieving better results if no prior is introduced in the initialization solution. On the other hand, if the approach is initialized with the mixed solution including an a prior value, the improvement is not significant in terms of normalized errors. The proposed method decreases globally the reconstruction errors compared to the mixed algorithm applied with various priors [8] [9] . Then, the results suggest that the refractive index decrement for a non-homogeneous object can be retrieved more accurately in terms of global errors if the non-linearity of the phase problem is taken into account. Though the linear solution is necessary for the initialization of the algorithm, this approach is expected to open new possibilities for biomedical studies with phase tomographic imaging.