Non-linear Phase Tomography Based on Fréchet Derivative

Phase imaging coupled to micro-tomography acquisition has emerged as a powerful tool to investigate specimens in a non-destructive manner. While the intensity data can be acquired and recorded , the phase information of the signal has to be " retrieved " from the data modulus only. Phase retrieval is an ill-posed non-linear problem and regularization techniques including a priori knowledge are necessary to obtain stable solutions. Several linear phase recovery methods have been proposed and it is expected that some limitations resulting from the linearization of the direct problem will be overcome by taking into account the non-linearity of the phase problem. To achieve this goal, we propose and evaluate a non-linear algorithm for in-line phase micro-tomo-graphy based on an iterative Landweber method with an analytic calculation of the Fréchet derivative of the phase-intensity relationship and of its adjoint. The algorithm was applied in the projection space using as initialization the linear mixed solution. The efficacy of the regularization scheme was evaluated on simulated objects with a slowly and a strongly varying phase. Experimental data were also acquired at ESRF using a propagation-based X-ray imaging technique for the given pixel size 0.68 μm. Two regularization scheme were considered: first the initialization was obtained without any prior on the ratio of the real and imaginary parts of the complex refrac-tive index and secondly a constant a priori value was assumed on r δ β. The tomographic central slices of the refractive index decrement were compared and numerical evaluation was performed. The non-linear method globally decreases the reconstruction errors compared to the linear algorithm and is achieving better reconstruction results if no prior is introduced in the initialization solution. For in-line phase micro-tomography, this non-linear approach is a new and interesting method in biomedical studies where the exact value of the a priori ratio is not known.


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 possibility to study objects with either negligible absorption or dense objects with small density variations.Moreover, 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 1 0 D = 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 nonnegligible since large propagation distances and high spatial resolution are required.Consequently, the nonlinearity of the phase-intensity relationship is a crucial aspect.New algorithms which take into account the nonlinearity of the inverse problem for the radiographic case have been proposed recently [14]- [17].These nonlinear approaches are very promising and lead to a large decrease of the reconstruction errors. ).For each angle θ , the phase map is retrieved using the 4 phase contrast images.Starting from these phase maps, the filter back projection is applied and the 3D refractive index decrement reconstruction is obtained.
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 ϕ .In this paper, we first summarize this new multi-image non-linear ( ) NL scheme, and then detail the results obtained on simulated images and for a tomographic reconstruction on a real multi-material 3D object.

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]: with ( ) are considered as the projections of the absorption index β and of the refractive index decrement r δ respec- tively.The linear integral relationships used for the tomographic reconstruction, are given by The intensity distribution where  denotes the coordinates in a plane perpendicular to the propagation direction z and * the 2D con- volution of the complex transmission function

Non-Linear Inverse Problem-Phase Retrieval
As detailed in [16] [17] [19], assuming that the phase ( ) : 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 2 L norm of the phase gradient The aim of this non-linear approach ( ) NL is to minimize for each projection angle θ the following cost functional:  is the measured noisy intensity at distance n D (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]: .
The classical Landweber method is modified with a variable step ,k θ τ 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 that can be obtained using the explicit calculation [16]: where * denotes the convolution operator.Moreover, by using the standard definition of the scalar product in 2 L spaces the adjoint operator ( ) is given by [16]: 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.

Initialization and Stopping Rules
It was shown that the non-linear algorithm improves the solution obtained with a linear algorithm in the radiographic case for simulated data [16] [17].Yet, an initialization obtained with the mixed linear scheme is necessary to obtain convergence of the non-linear method.In this work, for each projection angle, the NL algorithm was initialized either with the phase retrieved without any a priori knowledge or with a fixed a priori value of the ratio r δ β [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: , 1 , , 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 regularization parameter is gradually refined by trial-and-error with a decreasing interval.For a well-chosen parameter, 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.

Simulations and Data Acquisition
The new non-linear inversion method has been tested on simulated images and on experimental data for a multimaterial object.

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 r δ and for the refractive index β of different materials at 24 keV were used ( 0.5166 λ = Å ) in different regions (Table 1).Analytical projections were calculated in a parallel beam geometry with 2048 2048× pixels and the two resulting data sets were combined to form a complex representation of the wave exiting the object using Equation (2).Propagation in free-space was simulated using Equation ( 3  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: where max f is the maximum signal amplitude and max n 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.

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 intensity patterns for 1500 projection angles were recorded using a FRELON CCD camera with 2048 2048× pixels for the energy 22.5 keV at four short distances Al O wires and 28 μm Polypropylene ( ) PP fibres.Phase retrieval with the mixed approach was applied without any prior on r δ β [8] (initialization (A)) and with 367 r δ β = [9] corresponding to aluminium (initialization (B)).

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 ( ) where k ϕ is the phase recovered at iteration k 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 noisefree 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  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.

Experimental Data for Non-Linear Phase Tomography
The reconstructed projections for the angle of view 120 θ =  retrieved with the mixed algorithms in these two cases are displayed in Figure 5(a 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 r δ β , 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 2π r δ λ 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: 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 64% .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 r δ β 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 algorithm is very efficient to improve the uniformity of the reconstructed image as we can see in Figure 6.The nonlinear solution retrieved using as starting point the linear solution with 367 r δ β = provides more accurate reconstructions (Figure 6(d)).In this case, the NE [%] value corresponding to Al is overestimated, but the NE [%] for Al 2 O 3 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 r δ β .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.

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.

Figure 1 .
Figure 1.Experimental set-up for propagation-based technique or in-line phase contrast imaging technique for a parallel X-ray beam.The incident field is assumed to have a degree of partial coherence and passes through a probed sample.Phase contrast images will be registered on the CCD-based detector for different distances n D in the Fresnel field.

Figure 2 .
Figure 2. Principles of phase tomography.For each sample-to-detector distance ( 1 D (blue dashed line), 2 D (red dashed line), 3 D (green dashed line) and 4 D (purple dashed line)) 2D phase contrast images for Shepp Logan phantom are acquired.For each projection distance, the sample is rotated over minimum 180  and different 2D projection angles are considered (three angles are displayed 0 0 θ =  , by the following expression: . In this study, better convergence results were obtained by replacing this term with the phase

Figure 3 .
Figure 3. (a) Ideal phase to be retrieved and (b) absorption image with PPSNR = 24 dB for strongly varying phase.
field of view was 1.4 mm for the given pixel size 0.68 μm.The multi-material object used is composed of 125 μm Aluminium ( )

Figure 4 .
Figure 4. Normalized mean square error for the phase versus iteration number.(a) Strongly varying phase without noise; (b) Weakly varying phase without noise; (c) Strongly varying phase with PPSNR = 24 dB; (d) Weakly varying phase with PPSNR = 24 dB.
) 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 r δ is reconstructed.

Figure 5 .
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 367 r δ β = (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 [ 30 30] −

Figure 6 .
Figure 6.Tomographic central slice reconstructed with the mixed algorithm (a) without a priori information [8] and (b) with a priori information 367 r δ β = (corresponding to aluminium) [9].Corresponding central slice obtained with the non-linear algorithm initialized with the linear solution (c) without a priori (initialization displayed in (a)) and (d) with a priori information 367 r δ β =

Table 1 .
Values of the absorption coefficient and refractive index at 24 keV for the materials used in the 3D phantom.

Table 3 .
Theoretical and measured values with different algorithms 2π r

Table 4 .
Values for relative standard deviation (RSD) and normalized error (NE) obtained with different algorithms.