## Abstract

It is an effective scheme to the phase retrieval for axial intensity derivative computing. In this paper, we demonstrate a method for estimating the axial intensity derivative and improving the calculation accuracy in the transport of intensity equation (TIE) from multiple intensity measurements. The method takes both the higher-order intensity derivatives and the noise into account, and minimizes the impact of detecting noise. The simulation results demonstrate that the proposed method can effectively reduce the error of intensity derivative computing.

©2012 Optical Society of America

## 1. Introduction

Phase retrieval from intensity measurements is a noninterferometric technique suggested originally by Teague in 1983 for wavefront sensing [1]. Knowledge of the phase is helpful for correction of distortion produced by the turbulent atmosphere in some related fields such as free space optical communication (FSO) and astronomic observation.

Teague introduced the transport of intensity equation (TIE) which relate the object-plane to the axial derivative of the intensity distribution,$\partial I(x,y,z)/\partial z$, in the Fresnel region to us, and if the intensity distribution has no zeros, the TIE has a unique solution (with an additive arbitrary constant). Many effective methods such as Fast Fourier Transform (FFT) [2, 3], the multigrid (MG) [4], the Zernike polynomial expansion methods [5, 6], and the Green’s function based method [1] have been proposed for solving it. And all these methods require two images that are defocused from each other in order for the axial derivation of intensity derivative. However, there are tradeoffs among the amount of defocus, the accuracy of the result and noise considerations.

In recent years, some algorithms for derivation of intensity derivative from intensity measurements in multiple planes are proposed: In 2007, Marcos Soto [7] evaluated the variance of the detecting noise and the axial second-order intensity derivative with multi-plane intensity distribution, based on which a method minimizing the error is proposed; and in 2010 Laura Waller [8] of MIT took the impact of higher order derivative into account and described it with a matrix transformation; a year later, Bindang Xue [9, 10] of Beihang University in China extended Laura Waller’s work to the situation with intensity distribution measured in unequally-spaced planes.

This paper proposed a novel method to evaluate intensity derivative, which takes both the detecting noise and higher order derivative into consideration. Based on Taylor expansions of the measured intensities, the axial intensity derivative in the TIE is expressed by a linear combination of the multiple intensity measurements, which is used to estimate the difference from the actual value. Numerical simulations are conducted to test it.

## 2. Multiple intensity measurements for phase retrieval

An optical plane wave $u(r)$ on the x-y plane can be described by optical intensity and phase as

where $r=(x,y)$ denotes the coordinate in the x-y plane, perpendicular to the z axis, $I(r)$ and $\varphi (r)$ are the intensity and phase distribution in the plane respectively. In this paper, the focal plane is assumed to be located at$z=0$.Under the paraxial approximation, the transport of intensity equation (TIE) can be derived as

where $k=2\pi /\lambda $ is the wave number, and $\lambda $ denotes the wavelength, $\nabla =\left\{\frac{\partial}{\partial x},\frac{\partial}{\partial y}\right\}$ is the two dimensional gradient operator in the x-y plane.Assume the intensity distribution of *m* x-y planes located at ${z}_{1},{z}_{2},\cdots {z}_{m}$ is available, based on the Maclaurin formula, $I({z}_{i})$ (ignore the x-y dimension) can be expanded as

^{th}order derivative, and ${\xi}_{i}$ is some number located between 0 and ${z}_{i}$. As a result, we get m equations. Multiply these equations with a number ${a}_{i}$ separately and sum their left and right hand side respectively, we can get the following expression

Let

Ignoring the second term on the right hand of Eq. (6),the intensity derivative can be calculated in practical system like

If the detecting noise is taken into consideration, the difference from Eq. (6) is

Generally, Eq. (5) has another description form

Obviously, the first matrix of left hand side called transformation matrix in this paper is a Vandermonde matrix.

- (1) When$n+1>m$, Eq. (12) are an overdetermined equations. Under this condition, it pointless to evaluate the n
^{th}order derivative using intensity data of m planes. - (2) When$n+1=m$, it is the case that Bindang Xue described, the coefficient ${a}_{1}~{a}_{m}$ can be calculated definitely from Eq. (12), when the detecting noise cannot be suppressed.
- (3) When$n+1<m$, Eq. (12) are an underdetermined equations, it is necessary to lay our attention on this situation which has not been considered.

We can divide the matrix of Eq. (12) as

Then $D$ can be written as

As a result,

To minimize it,

^{th}element of the vector $F={B}^{-1}E$, and ${G}_{j}$ is j

^{th}row vector of matrix $G={B}^{-1}A$. Obviously, Eq. (16) is equivalent to

Combining with Eq. (13), we get

In Eq. (18), ${G}^{T}G+I$ and $B$ are both nonsingular matrix, so the coefficient ${a}_{1}~{a}_{m}$can be derived easily.

## 3. Simulations

We take the “Gauss Beam”, a widely used beam mode, to check the validity of the method described above.

That is the 0th order Gauss mode

As a result, the intensity distribution of the axis can be written as

whereThe simulation error of Bindang Xue’s method (we call it high order derivative suppression method here) and the algorithm with noise suppression proposed in this paper are both shown in Fig. 1 .

When sampling interval $\Delta z$is large, the error raised by high order derivative impacts the final result significantly. Bindang Xue’s method plays a good role in this situation as shown in Fig. 1(a). But as it decreases, the influence of high order derivative vanishes gradually, while the detecting noise takes the dominated place step by step. As described in Fig. 1(a)–1(d) the performance of noise suppression method improves significantly.

Figure 2 shows the simulation error of noise suppression method for different n when m = 16. When n is small, the capability of noise suppression is strong, but not that of high order derivative suppression, which leads to big error from the actual value. As n increases, these two types of suppression reach balance point step by step. Under this condition, the error caused by the algorithm decreases, and finally reaches the lowest point. But when the growth of n is going on, the dominant position of error is taken by detecting noise gradually, which makes the error rise. From what is discussed above, we can safely get the conclusion that there exist an optimal n for a sampling number m, so it is necessary to search this point in practice.

## 4. Conclusions

This paper proposed a novel method to derive the intensity derivative, which extends the work of Marcos Soto etc [7–10]. Including higher order intensity derivatives and detecting noise improves the accuracy of intensity derivative computing, which leads to enhancement of the algorithm performance based on TIE. By adopting correction for nonlinearities (higher order derivatives) and noise suppression, large stacks of data spanning longer distances can be use to estimate the derivative more precisely.

## Acknowledgments

This work is supported by the National Natural Science Foundation of China (No.61077058).

## References and links

**1. **M. Reed Teague, “Deterministic phase retrieval: a Green's function solution,” J. Opt. Soc. Am. **73**(11), 1434–1441 (1983). [CrossRef]

**2. **D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. **80**(12), 2586–2589 (1998). [CrossRef]

**3. **T. E. Gureyev and K. A. Nugent, “Rapid quantitative phase imaging using the transport of intensity equation,” Opt. Commun. **133**(1-6), 339–346 (1997). [CrossRef]

**4. **L. J. Allen and M. P. Oxley, “Phase retrieval from series of images obtained by defocus variation,” Opt. Commun. **199**(1-4), 65–75 (2001). [CrossRef]

**5. **T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation. II. orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A **13**(8), 1670–1682 (1996). [CrossRef]

**6. **T. E. Gureyev, A. Roberts, and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation: matrix solution with use of Zernike polynomials,” J. Opt. Soc. Am. A **12**(9), 1932–1942 (1995). [CrossRef]

**7. **M. Soto and E. Acosta, “Improved phase imaging from intensity measurements in multiple planes,” Appl. Opt. **46**(33), 7978–7981 (2007). [CrossRef] [PubMed]

**8. **L. Waller, L. Tian, and G. Barbastathis, “Transport of intensity imaging with higher order derivatives,” Opt. Express **18**(12), 12552–12561 (2010). [CrossRef] [PubMed]

**9. **B. D. Xue, S. L. Zheng, L. Y. Cui, X. Z. Bai, and F. G. Zhou, “Transport of intensity phase imaging from multiple intensities measured in unequally-spaced planes,” Opt. Express **19**(21), 20244–20250 (2011). [CrossRef] [PubMed]

**10. **S. L. Zheng, B. D. Xue, W. F. Xue, X. Z. Bai, and F. G. Zhou, “Transport of intensity phase imaging from multiple noisy intensities measured in unequally-spaced planes,” Opt. Express **20**(2), 972–985 (2012). [CrossRef] [PubMed]