Received 16 September 2014; revised 21 October 2014; accepted 3 November 2014

ABSTRACT

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 refractive index and secondly a constant a priori value was assumed on. 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.

Keywords:

Phase Retrieval, In-Line Phase Tomography, Inverse Problems, Non-Linear Problem, Non-Linear Optimization, Fréchet Derivative, Coherent Imaging, Fresnel Diffraction, Phase Contrast, X-Ray Imaging

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.

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

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

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:

(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 Al₂O₃ 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.

References