# Non linear phase retrieval from fresnel diffraction patterns using the frechet derivative

B. Sixou, V. Davidoiu, M. Langer, F. Peyrin

► To cite this version:

B. Sixou, V. Davidoiu, M. Langer, F. Peyrin. Non linear phase retrieval from fresnel diffraction patterns using the frechet derivative. IEEE International Symposium on Biomedical Imaging: From Nano to Macro, Mar 2011, Chicago, United States. pp.1370-1373, 10.1109/ISBI.2011.5872655. hal-01825259

HAL Id: hal-01825259

https://hal.science/hal-01825259

Submitted on 14 Apr 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# NON LINEAR PHASE RETRIEVAL FROM FRESNEL DIFFRACTION PATTERNS USING THE FRECHET DERIVATIVE

B. Sixou*, V. Davidoiu*, M. Langer*†, F. Peyrin*†

* Université de Lyon, CREATIS-LRMN ; CNRS UMR5220 ; Inserm U630 ; INSA-Lyon ; Université Lyon 1, F-69621 Villeurbanne Cedex, France

† European Synchrotron Radiation Facility, 6 rue Jules Horowitz, F-38043 Grenoble Cedex, France

## ABSTRACT

Several methods of phase retrieval for in line phase tomography have already been investigated based on the linearization of the relation between the phase shift induced by the object and the diffracted intensity. They use the Transport Intensity Equation (TIE) or the Contrast Transfer Function (CTF), or mixed approaches. In this work, we present a non linear iterative approach using the Frechet derivative of the intensity recorded at a few number of propagation distances. The inverse problem is regularized with the smoothing $L_2$ norm of the phase gradient. The evaluation of the method was performed using a simple phase map, both with and without noise. Our approach outperforms the linear methods on simulated noisy data up to high noise levels.

**Index terms** – X-ray Imaging, Coherent Imaging, Phase Contrast, Phase Retrieval, Nonlinear Optimization

## 1. INTRODUCTION

Phase contrast imaging is a new X-ray imaging modality enhancing the visibility of low absorption objects and weakly contrasted structures. Phase contrast can be achieved by propagation when using coherent X-ray beams, for example extracted from synchrotron radiation. The relationship between the phase shift induced by a sample and the intensity recorded at some sample-to-detector distance can be modeled by the Fresnel diffraction theory. This quantitative relationship can be used to calculate the phase shift induced by the object from intensity radiographs recorded at different propagation distances, a process called phase retrieval. When coupled to X-ray tomography, it is possible to recover a three dimensional map of the refractive index decrement.

Several phase retrieval methods have been proposed, using the Transport Intensity Equation [1, 2, 3], the Contrast Transfer Function (CTF) [4], the mixed approach between the two former methods [5], and a TIE-based approach that combine the phase retrieval step with the inverse Randon transform [6]. A common problem in phase retrieval is low frequency artifacts. To minimize this problem, regularization techniques based on a homogeneous prior on the object have

previously been proposed [7]. However, this method is only appropriate for some classes of objects. Further, all previous approaches are based on linearized relationships between the phase and the intensity, thus involving some approximations in the direct problem $\mathcal{I}(\varphi)$ where $\varphi$ is the phase shift map.

In this work, we investigate the resolution of the non linear inverse problem, which has so far not been considered. We propose a new iterative method based on the Frechet derivative of the $\mathcal{I}(\varphi)$ relationship. Further, the phase retrieval problem considered in this work is an ill-posed problem since the inverse operator $\mathcal{I}^{-1}(\varphi)$ is not a bounded. Thus some regularization is required. As detailed in the following, we search for a smooth solution and a continuous phase map. Failure of convergence and were observed when this a priori information on the smoothness of the solution was not included. First results are presented from simulations of noise-free and noisy data.

## 2. THE DIRECT PROBLEM OF IMAGE FORMATION

We consider an object illuminated with partially coherent X-rays of wavelength $\lambda$. Neglecting the diffraction within the object, the interaction of X-rays with matter can be described by a transmittance function $T$ of the coordinates $x = (x, y)$ in a plane perpendicular to the propagation direction $z$.

$$T(\mathbf{x}) = \exp[-B(\mathbf{x}) + i\varphi(\mathbf{x})] = a(\mathbf{x}) \exp[i\varphi(\mathbf{x})] \quad (1)$$

The absorption, $a(\mathbf{x})$, and phase shift $\varphi(\mathbf{x})$ induced by the object can be considered as projections of the absorption and refraction index respectively. The Fresnel diffracted intensity at a distance $D$ is given by the squared modulus of the exit wave:

$$I_D(\mathbf{x}) = |T(\mathbf{x}) * P_D(\mathbf{x})|^2 \quad (2)$$

where * denotes the convolution of the transmitance with the Fresnel propagator

$$P_D(\mathbf{x}) = \frac{1}{i\lambda D} \exp\left(i\frac{\pi}{\lambda D}|\mathbf{x}|^2\right) \quad (3)$$

$D$ being the propagation distance along $z$. The convolution can be efficiently calculated in the Fourier domain. The real
and imaginary part of the complex refractive index can be re-
constructed by recording the intensity at different angles and
with tomographic reconstruction algorithms. In this work, we
consider a single projection. An attenuation image can be ob-
tained by moving the detector to D = 0.

3. THE NON LINEAR INVERSE PROBLEM

**3.1. Initialization**

In order to save computing time, the algorithm is initialized
with the Mixed approach [5]. The derivation of this relation-
ship relies on the hypothesis that the absorption and phase are
slowly varying. In the fourier domain, the linearized problem
may be expressed as:

$$
\begin{equation}
\begin{aligned}
\tilde{I}_D(\mathbf{f}) = \tilde{I}_D^{\varphi=0}(\mathbf{f}) &+ 2 \sin(\pi \lambda D |\mathbf{f}|^2) \mathcal{F}(I_0 \varphi)(\mathbf{f}) \\
& + \frac{\lambda D}{2\pi} \cos(\pi \lambda D |\mathbf{f}|^2) \mathcal{F}[\nabla(\varphi \nabla I_0)](\mathbf{f})
\end{aligned}
\tag{4}
\end{equation}
$$

where $I_D^{\varphi=0}(\mathbf{f})$ is the intensity at distance $D$ if the phase was zero. To include several distances, a least squares minimization procedure is used. Assuming that the absorption is given, the phase retrieval problem can be stated with the phase-attenuation product $\psi(\mathbf{x}) = I_0(\mathbf{x})\varphi(\mathbf{x})$ as

$$
\psi(\tilde{\mathbf{f}}) = \arg\min_{\tilde{\psi}} \sum_{D} |A_D(\mathbf{f})\tilde{\psi}(\mathbf{f}) - \tilde{I}_D(\mathbf{f}) - \tilde{I}_D^{\varphi=0}(\mathbf{f}) - \Delta_D(\mathbf{f})|^2 \quad (5)
$$

with

$$
A_D(\mathbf{f}) = 2 \sin(\pi \lambda D |\mathbf{f}|^2) \quad (6)
$$

and

$$
\Delta_D(\mathbf{f}) = \frac{\lambda D}{2\pi} \cos(\pi \lambda D |\mathbf{f}|^2) \mathcal{F}\{\nabla[\psi \nabla \ln(I_0)]\}(\mathbf{f}) \quad (7)
$$

Several regularization scheme have been tested to solve
this linear inverse problem, like classical quadratic Tikhonov
regularization and wavelet shrinkage [5, 8, 7]. In this work,
we used a regularization method based on the L₂ norm of the
phase gradient.

**3.2. Frechet derivative of the intensity**

We assume that the operator $\mathcal{I}_D(\varphi)$ is a non linear Frechet
differentiable operator in its domain. In the following, we will
consider the phase has a bounded support $\Omega$ and that its do-
main is the functional space [9]:

$$
W_{\diamond}^{2,2}(\Omega) = \{\varphi \in W^{2,2}(\Omega), \frac{\partial \varphi}{\partial \vec{n}} = 0\} \quad (8)
$$

where $\frac{\partial \varphi}{\partial \vec{n}}$ is the normal derivative of the phase. In the case of the intensity-phase relationship, the Frechet derivative of the operator $\mathcal{I}_D(\varphi)$ at the point $\varphi_k$ is the linear operator defined by the relationship:

$$
\mathcal{I}_D(\varphi_k + \varepsilon) = \mathcal{I}_D(\varphi_k) + G_k(\varepsilon) + O(\varepsilon^2) \quad (9)
$$

The linear operator $\mathcal{I}'_D(\varphi_k)(\varepsilon) = G_k(\varepsilon)$ can be given explicitly as:

$$
G_k(\varepsilon) = \{-[ia\varepsilon \exp(-i\varphi)] * P_D\}\{[(a \exp(i\varphi)] * \overline{P_D}) \\
+ \([a \exp(-i\varphi)] * P_D\)\{[ia\varepsilon \exp(i\varphi)] * \overline{P_D}\} \quad (10)
$$

where * denotes the convolution operator.

**3.3. A Landeweber type iterative method**

As a first approach to retrieve the phase from the intensity
measurements and to obtain a rather smooth solution, we con-
sider the problem of minimizing a simple Tikhonov's func-
tional of the following form:

$$
J_{\alpha}(\varphi) = \frac{1}{2} \| \mathcal{I}_D(\varphi) - I_D \|_{L_2(\Omega)}^2 + \frac{\alpha}{2} \| \nabla \varphi \|_{L_2(\Omega)}^2 \quad (11)
$$

where $\alpha$ is a regularizing parameter. The optimality condition
is then:

$$
\mathcal{I}'_D(\varphi)^*[\mathcal{I}_D(\varphi) - I_D] - \alpha\Delta\varphi = 0 \quad (12)
$$

where $\mathcal{I}'(\varphi)^*$ is the adjoint of the Frechet derivative of the intensity. This optimality condition defines the descent direction of our steepest descent iterative method and the next iterate is obtained with:

$$
\varphi_{k+1} = \varphi_k - \tau_k \nabla J_{\alpha_k}(\varphi_k). \quad (13)
$$

Starting from the current phase estimate $\varphi_k$ at the iteration $k$, a linear search procedure is introduced with a variable step $\tau_k$ yielding the following modification of the standard Landweber method:

$$
\varphi_{k+1} = \varphi_k - \tau_k \{ \mathcal{I}'_D(\varphi_k)^* [\mathcal{I}_D(\varphi_k) - I_D] - \alpha_k \Delta \varphi_k \} \quad (14)
$$

The step length parameter $\tau_k$ is chosen in order to minimize the Tikhonov's functionnal along the descent direction:

$$
\tau_k = \operatorname{argmin} J_{\alpha_k}(\varphi_k - \tau \delta_k), \quad (15)
$$

where $\delta_k = \mathcal{I}'_D(\varphi_k)^* (\mathcal{I}_D(\varphi_k) - I_D) - \alpha_k \Delta \varphi_k$ is the descent direction. An approximate value is obtained with a dichotomy strategy. We have chosen a geometrically decreasing sequence of regularization parameters $\alpha_k$. The initial value of the regularizing parameter is chosen such that the two terms of the Tikhonov's functionnal are of the same order of magnitude. The regularization parameter $\alpha_k$ is decreased when the Tikhonov's functionnal $J$ stagnates. The next iterate is accepted if both the $L_2$ norm of the gradient and the discrepancy term decrease.

The described algorithm involves a single distance. In order to include the intensity maps obtained for the different distances, we consider a cyclic iteration over several distances, one iteration being performed for each distance. The use of several distances should improve the reconstruction since it allows a better coverage of the frequency domain and it improves the statistics.
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Original phase map to be retrieved.</figcaption></figure>

## 4. SIMULATIONS

The diffracted intensity was simulated as described in [5]. Two phantoms were defined, one for the absorption coefficient and one for the phase shift based on experimental values of the absorption coefficient and refractive index of a test material at 24 keV X-ray energy. Propagation in free space was simulated in the Fourier domain using Eq. 2 for the eight propagation distances 0.2 m, 0.4 m, 0.6 m, 0.8 m, 1.2 m, 1.4 m and 1.6 m. The images obtained were downsampled to $N = 74 \times 74$ pixels. The original phase map to be retrieved is displayed in Figure 1. The phase contrast images were corrupted with additive Gaussian white noise with peak to peak signal to noise ratios (PPSNR) of 20 dB, 0 dB and no noise.

As a first approach, in order to calculate the gradient of the potentially noisy intensity function, we have applied the implicit filtering method described by Kelley et al. [10]. In its simplest form, implicit filtering is the steepest descent algorithm with finite difference gradients, where the difference increment varies as the iteration progresses. Because the gradient is only an approximation, the computed descent direction may fail to be a descent direction, and the line search may fail. In this case the difference increment is reduced. The derivatives in the gradient $\mathcal{I}'_D$ are thus approximated by centered differences formulae. This finite difference method requires many evaluations of the Tikhonov's functional. The Frechet

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Fresnel diffraction pattern at propagation distance $D=1.4$ m.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Reconstructed phase for our Landweber type method and noise free simulations.</figcaption></figure>

derivative $G_k$ calculated with the finite difference method at the point $\varphi_k$ is a thus matrix of $\mathbb{R}^{N \times N}$.

The mixed approximation of the forward problem (Eq. 4) is used as the starting point of our reconstructions. This initial guess is close enough to the ideal result to ensure convergence of our algorithm. The decrease of the phase error reflects the improvements over traditional linear methods. In the iteration, we have chosen to not introduce any priori information on the phase values. The iteration was stopped when the data term $||\mathcal{I}_D(\varphi_k) - \mathcal{I}_D||$ was equal to the noise level.

## 5. RESULTS

A noise-free Fresnel diffraction pattern for $D = 1.4$ m is shown in Figure 2. The phase maps obtained with our method are displayed in Figure 3 and Figure 4, without noise and with noise at 20dB. The cyclic finite difference Landweber type method has been tested on noise free and noisy data. If $\varphi^*$ is the phase to be recovered, the normalized least square error of $L_2(\Omega)$ norm $||\varphi_k - \varphi^*||/||\varphi^*||$ is used to measure the quality of the phase recovery. The results obtained with noise free and noisy data are displayed in Figure 5 as a function of the number of iterations. Our iterative phase retrieval algorithm retrieves the phase better than the linear Mixed approach. The errors on the phase have been significantly reduced for the noise free simulations.

## 6. DISCUSSION AND CONCLUSION

In this work we have presented a new non linear iterative method of phase retrieval for in line phase tomography based on the Frechet derivative of the phase intensity relationship. As a first approach to obtain a smooth solution, a regularization term with the norm $L_2$ of the gradient is introduced. This method uses the intensity radiographs obtained for different propagation distances in a cyclic iteration. The Frechet derivative is evaluated with a finite difference method with a phase increment adapted to the noise level.

The evaluation was performed using a numerical phan-
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Reconstructed phase for our Landweber type method and 20 dB PPSNR.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Evolution of the mean square error for the phase for the noise free (i) and the noisy data with PPSNR=20 dB (ii).</figcaption></figure>

tom, used to simulate phase contrast tomography data, both with and without noise. For the simulated data, the normalized mean square error was measured. Our approach outperforms the mixed approach and it performs well on simulated noisy data. Thus, if the accuracy of the reconstruction is the primary goal, especially in a quantitative analysis, this method could be very interesting.

There are many perspectives to this work. A particular difficulty of this problem is to take into account several propagation distances. We applied a cyclic iteration over the distances, whereas other iteration schemes could be utilized. Since the long computing time could limit the applicability of the present method, further improvements will be proposed to speed up the algorithm. In addition, a quadratic regularization term has been used. The drawbacks of the regularization function $||∇φ||_2$ are well-known. An isotropic smoothing effect is obtained and the boundaries are not well preserved. The noise is suppressed but the high values of the gradient are too penalized on the edge. The gradient $∇φ$ may be replaced by a non linear function of $φ$ corresponding to non linear or anisotropic diffusion [9] and other norms will be investigated. This type of regularization will also be considered in future works. Finally, it is expected to test the method on experimental data acquired at the ESRF.

## 7. REFERENCES

[1] D.M. Paganin, *Coherent X-Ray Optics*, Oxford University Press, New York, 2006.

[2] S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. W. Stevenson, “Quantitative phase imaging using hard x-rays,” *Phys. Rev. Lett.*, vol. 77, pp. 2961–2964, 1996.

[3] L. J. Allen and M. P. Oxley, “Phase retrieval from series of images obtained by defocus variation,” *Opt. Comm.*, vol. 199, pp. 65–75, 2001.

[4] P. Cloetens, R. Barrett, J. Baruchel, J. P. Guigay, and M. Schlenker, “Phase objects in synchrotron radiation hard x-ray imaging,” *J. Phys. D*, vol. 29, pp. 133–146, 1996.

[5] M. Langer, P. Cloetens, J. P. Guigay, and F. Peyrin, “Quantitative comparison of direct phase retrieval algorithms in in-line phase tomography,” *Medical Physics*, vol. 35, pp. 4556–4565, 2008.

[6] A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” *J. Opt. Soc. Am. A*, vol. 19, pp. 472–480, 2002.

[7] M. Langer, P. Cloetens, and F. Peyrin, “Regularization of phase retrieval with phase attenuation duality prior for 3d holotomography,” *IEEE Trans. Image Process.*, vol. 19, pp. 2425–2436, 2010.

[8] M. Langer, P. Cloetens, and F. Peyrin, “Fourier-wavelet regularization of phase retrieval in x-ray in line phase tomography,” *J. Opt. Soc. Am. A*, vol. 28, pp. 1877–1882, 2009.

[9] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen, *Variational Methods in Imaging*, Springer Verlag, New York, 2008.

[10] C. T. Kelley, *Iterative methods for optimization*, SIAM, Philadelphia, 1999.