## Abstract

The noise problem is generally inevitable for phase retrieval by solving the transport of intensity equation (TIE). The noise effect can be alleviated by using multiple intensities to estimate the axial intensity derivative in the TIE. In this study, a method is proposed for estimating the intensity derivative by using multiple unevenly-spaced noisy measurements. The noise-minimized intensity derivative is approximated by a linear combination of the intensity data, in which the coefficients are obtained by solving a constrained optimization problem. The performance of the method is investigated by both the error analysis and the numerical simulations, and the results show that the method can reduce the noise effect on the retrieved phase. In addition, guidelines for the choice of the number of the intensity planes are given.

© 2012 OSA

## 1. Introduction

Phase imaging is an important technique for studying phase objects in biological or material research. TIE [1, 2] provides a convenient and deterministic way for quantitative phase reconstruction from intensity measurements. It relates the axial intensity derivative to the object-plane phase by a second order elliptic partial differential equation, and the phase can be reconstructed uniquely (within an additive constant) by solving the equation with appropriate boundary conditions [3–5]. Many effective numerical methods have been proposed to solve the TIE, such as the Green’s function based method [1, 5], Zernike polynomial expansion method [6], fast Fourier transform (FFT) method [7, 8], and multigrid method [8, 9]. Therefore, TIE has been applied in a variety of applications including adaptive optics [9, 10], microscopy [11, 12] and X-ray optics [13, 14] etc. Although different methods may retrieve the phase distribution with different precisions, a fundamental factor which determines the accuracy of the retrieved phase is the precision of the axial intensity derivative. The intensity derivative is generally calculated from two defocused intensity measurements by the finite difference method. However, the main problem inherent in this method is that it ignores the nonlinear effect of the higher-order intensity derivatives and is vulnerable to the noise effect [15]. Multiple intensity measurements based methods [16–18] can mitigate the nonlinear effect by estimating the higher-order intensity derivatives and thus can improve the accuracy of the retrieved phase. However, these methods still suffer from the noise problem. Although the noise can be reduced by multiple intensity measurements weighted average in one plane [16] or by de-noising techniques [17, 19], this will add inconvenience and complication to the calculation. In this case, the noise-suppressed methods [20] are more preferred to estimate the intensity derivative from multiple intensity measurements. In addition, the intensity derivative is mostly computed from multiple intensities measured in equally-spaced planes. These methods, however, cannot apply in arbitrarily-spaced intensity measurements. Recently, the authors have developed a method which uses multiple unevenly sampled noise-free intensities to obtain a higher-order approximation to the TIE [18]. However, this method did not consider the noise problem.

In this study, the noisy intensity measurements based TIE is proposed. In this equation, the intensity derivative is approximated by a linear combination of the multiple unequally-spaced noisy intensities. The weight coefficients are derived by minimizing the mean square noise error of the estimated and the ideal intensity derivatives, and then expressed by explicit formulas. This approach can be regarded as a generalization of the method presented in [20], which adopted multiple equally-spaced intensity images to estimate the intensity derivative. In this work, the noise-reduction property of the method is analyzed, and the error caused by the higher-order intensity derivatives is also discussed to provide guidelines for the choice of the newly used intensities. Numerical simulations are conducted to test the method using multiple noisy intensities contaminated by different levels of Gaussian noise and Poisson noise.

## 2. TIE with unequally-spaced noisy intensity measurements

Suppose a coherent monochromatic plane wave with a wavelength $\lambda $propagates along the $Z$ axis in the coordinate system $(r,z)$, where $r=(x,y)$ is the spatial coordinates in the x-y plane, orthogonal to the $Z$ axis. With the paraxial approximation, TIE can be derived from the parabolic wave equation [1]:

where $k=2\pi /\lambda $ is the wave number, $\nabla =\left(\partial x,\partial y\right)$ is the two-dimensional differential operator in the x-y plane. ${I}_{z}(r)$ and ${\varphi}_{z}(r)$ are the intensity and the phase distributions respectively.Estimating the axial intensity derivative is essential to retrieving a precise phase by solving Eq. (1) [16]. Let us assume that the unevenly sampled intensities ${I}_{{z}_{j}}$and ${I}_{{z}_{-j}}$, with the defocused distances ${z}_{j}$and ${z}_{-j}$, $(j=1,2,\mathrm{...},M)$, are given by

where, ${\overline{I}}_{{z}_{j}}(r)$and ${\overline{I}}_{{z}_{-j}}(r)$represent the ideal intensities. ${n}_{{z}_{j}}$, ${n}_{{z}_{-j}}$are supposed to be zero-mean, additive Gaussian noise with variance${\sigma}^{2}$, independent of the intensities. It should be mentioned that ${z}_{j}$and ${z}_{-j}$, $(j=1,2,\mathrm{...},M)$ can be arbitrary values at which TIE is valid. The smaller of them, the closer of the intensity planes to the focal plane.Then, the intensity derivative in the focal plane$z=0$, i.e.$\partial {I}_{0}(r)=\partial {I}_{z}(r)/\partial z|{}_{z=0}$, can be estimated by:

It should be mentioned that ${z}_{-j}$ can be negative or just be ${z}_{0}$, which means that we can estimate the axial intensity derivative independently from either the double-sided or single-sided diffraction intensities. The coefficients ${a}_{j},j=1,2,\mathrm{...},M$ must minimize the mean square error between the estimated and the real intensity derivatives

Then, substituting Eq. (8) into Eq. (5) we can get

withIn this study, our objective is to minimize the noise effects on the estimated intensity derivative. Thus, from Eq. (11) we obtain:

Equation (9) and Eq. (13) can construct a nonlinear optimization problem with an equality constraint. With the Lagrange multiplier method, the coefficients can be calculated by

Substituting Eq. (14) into Eq. (4) and then substituting the obtained equation into Eq. (1), we can get the generalized TIE with unequally-spaced noisy intensity measurements as:

Equation (15) can be regarded as a generalization of Eq. (3) in [20] from equally-spaced intensities to unequally-spaced intensities. It can be verified easily that when the measured intensities are equally-spaced, i.e.:

The coefficients for the intensity measurements become:

which correspond to that in Eq. (9) in [20].In the next section, we will analyze the effects of the defocused values on ${\epsilon}_{h}{}^{2}(r)$and${\epsilon}_{h}{}^{2}(r)$. In addition, we will discuss the relationship between the number of the intensity planes and the accuracy of the estimated intensity derivative.

## 3 Accuracy analysis of the estimated intensity derivative

Substituting Eq. (14) into Eq. (11) and Eq. (12), we have:

As shown in Eq. (19), the noise error becomes smaller with the increase of $M$. However, ${\epsilon}_{h}{}^{2}(r)$ is indeterminable as $M$increases. This is because, firstly, the second-order intensity derivatives cannot be precisely known in advance, and secondly, the sign of $\left({z}_{j}-{z}_{-j}\right),j=1,2,\mathrm{...},M$is uncertain when ${z}_{j}$and${z}_{-j}$are arbitrary values. Consequently, we cannot make sure whether ${\epsilon}_{h}{}^{2}(r)$is decreasing by increasing the number of the intensity planes. In this study, we will provide a general idea of it by analyzing ${\epsilon}_{h}{}^{2}(r)$with some prescribed assumptions.

Let us assume that ${z}_{j}$ is positive and ${z}_{-j}$is negative. The lower bound and the upper bound of ${\epsilon}_{h}{}^{2}(r)$ are adopted to analyze its effects, and they can be derived from Eq. (20) as:

withThe accuracy of the estimated intensity derivative can be improved by using more intensity planes [16, 20], so we express ${\epsilon}^{2}(r)$ as

To further analyze${\epsilon}^{2}(r)$, we here assume ${\epsilon}_{h}{}^{2}(r)$ approximates to ${\epsilon}_{h(up)}{}^{2}$. Then substituting${\epsilon}_{h}{}^{2}(r)$ by ${\epsilon}_{h(up)}{}^{2}$so that Eq. (22) can be expressed as

In Eq. (23), if ${\epsilon}_{h(up)\_(M+1)}{}^{2}(r)-{\epsilon}_{h(up)\_M}{}^{2}(r)<0$, i.e. ${\epsilon}_{h(up)}{}^{2}$ is a decreasing function of $M$, The inequality is always true. From the expression of in Eq. (21), we know that this can be achieved by properly choosing the defocused values to make ${Z}_{up}{}^{2}(M)$decrease. As an example, we can choose the defocused values of the already and newly used intensities as:

This means that the newly used intensities with the defocused distances ${z}_{M+1}$ and ${z}_{-\left(M+1\right)}$ are closer to the focal plane than the already used intensities with defocused values ${z}_{j}$and ${z}_{-j}$, $(j=1,2,\mathrm{...},M)$. Then we have:

Consequently, we can get:

Therefore, ${Z}_{up}{}^{2}(M)$ decreases with $M$increasing.

However, If ${\epsilon}_{h(up)\_(M+1)}{}^{2}(r)-{\epsilon}_{h(up)\_M}{}^{2}(r)>0$, i.e., ${\epsilon}_{h(up)}{}^{2}$is an increasing function of $M$, The inequality in Eq. (23) may not be always true. For example, when the defocused values satisfy

In summary, ${\epsilon}_{n}{}^{2}(r)$ can be reduced by using more intensity measurements, and ${\epsilon}_{h}{}^{2}(r)$can be reduced when ${Z}_{up}{}^{2}(M)$decreases.

## 4. Simulations

Numerical simulations are conducted to test the derived formulas for the intensity derivative estimation. Figure 1 shows the simulated phase profile with phase value between 0 and 0.5$\pi $ radians (0 denotes the darkest pixels and 0.5$\pi $ refers to the brightest pixels). The dimensions of the image are 512 × 512 pixels with the assigned pixel size 4$\mu m$. The radiation wavelength is chosen to be 550nm. The intensity distribution in the input plane is uniform.

The intensities are sampled as two-dimensional images, so that the normalized mean square error (NMSE) is used to quantitatively evaluate the accuracy of the estimated intensity derivative, which is:

The simulated diffraction intensities in the different planes are calculated by the angular spectrum method [21]. Therefore, the ideal intensity derivative can be derived as (see Appendix A):

In this study, thirty diffraction intensity images are used o test the method with random defocused distances between $-28\mu m$and$29\mu m$.The positive defocused values of${z}_{j}$ are assigned to be$\text{29,27,26,23,22,20,19,17,16,12,11,10,8,5}\mu \text{m}$ and$3\mu m$. The negative defocus values of ${z}_{-j}$ are assigned to be $\text{-28,-25,-23,-21,-18,-16,-15,-13,-10,-9,-7,-6,-3}\mu \text{m}$ and $\text{-}2\mu m$. The maximum and minimum sampled spaces are respectively assigned to be$4\mu m$and$1\mu m$. These defocused values can make ${Z}_{up}(M)$diminish with $M$ increases. Intensities with Gaussian noise, Poisson noise and the mix of them will be simulated and tested respectively. For each numerical experiment, the method for each noise levels will be carried out for 60 times, and the averages of the corresponding NMSE and RMSE will be calculated.

#### 4.1 Estimate the intensity derivative from multiple unequally-spaced intensity measurements with Gaussian noise

We assume that the intensities are contaminated with the independent Gaussian white noise and four noise levels in the images are simulated. The corresponding signal-to-noise-ratio (SNR) of the images is 70db, 60db, 55db and 50db, respectively. Figure 2 shows an example of the diffraction images in the $z=19\mu m$ plane with different SNR.

Figure 3 shows the averaged NMSEs as a function of $M$ with different SNR. Figure 4 shows the retrieved phases with different number of noisy intensity measurements and with different SNR. In Table 1 , NMSEs and RMSEs are compared with different number of noisy intensity measurements and with different levels of SNR.

As shown in Fig. 3, the estimation accuracy is improved with the increase of $M$ and the NMSEs decrease monotonically. The NMSEs curves also show that the method is effective for phase retrieval by multiple intensities with much noise, since the NMSE for SNR of 50db is reduced from 0.150 to 0.028 (see Table 1), and the precision of the estimated derivative is improved about 81%. In Table 1, the NMSEs are diminished from 0.006, 0.019, and 0.051 to 0.002, 0.004, and 0.010, with an improvement of about 67%, 79%, and 80% respectively for SNR of 70db, 60db, 55db. However, the convergence speed of the NMSEs becomes slower as $M$ increases in Fig. 3. This can be explained by Eq. (19), which shows that the intensities with smaller defocused values have less impact on the noise error, so that the noise error reduces slower or even become stagnation as the defocused values decrease.

In Table 1, RMSEs show that the accuracy for the phases retrieved from more intensity measurements are much higher. The visual impression in Fig. 4 is that the cloudy effect in the retrieved phases induced by the measurement noise is largely removed when multiple intensities are adopted.

#### 4.2 Estimate the intensity derivative from multiple unequally-spaced intensity measurements with Poisson noise

In this simulation, we assume that the intensities are contaminated with Poisson noise, and the noise is simulated respectively with an average of 100,000, 80,000, 32,000, and 10,000 photons per pixel in the diffraction image. We denote the average photons per pixel in the images as$\gamma $. The less number of the average photons means that the images are measured under the lower-light conditions. In this case, the Poisson-noise problem of the images is more severe.

Figure 5 shows the noisy diffraction image in the $z=19\mu m$ plane. As shown, the intensity images with smaller $\gamma $ are with stronger noise. Figure 6 shows the averaged NMSEs as a function of $M$ and $\gamma $. As shown, the trends of the curves are similar to those in the Gaussian-noise case and the NMSEs for all the four cases are monotonically decreasing as $M$increases.

Quantitative comparisons are listed in Table 2 . As shown, RMSEs are smaller with more intensity measurements. This indicates that the accuracy of the reconstructed phases from more intensity measurements is much higher. Figure 7 shows the corresponding retrieved phases of the four cases above mentioned. The first row listed the retrieved phases retrieved from intensities with $M=1$ and the second row displayed those from intensities with $M=15$. Visually, the phases are much better when multiple intensities are used and more close to the simulated phase.

#### 4.3 Estimate the intensity derivative from multiple unequally-spaced intensity measurements with mixed noise

In this section, we assume that the intensities are contaminated with the mixed noise, i.e. a combination of the Poisson noise and white Gaussian noise. The noisy images are simulated by firstly adding the Poisson noise to the intensities and then the independent white Gaussian noise to the obtained Poisson-noise corrupted intensities. We denote the standard deviation of the Gaussian noise as$\sigma $. Four cases of the diffraction intensities are tested: (1) $\gamma =500,000$,$\sigma =0.001$; (2) $\gamma =500,000$,$\sigma =0.01$; (3) $\gamma =32,000$,$\sigma =0.001$; (4) $\gamma =32,000$,$\sigma =0.01$. Figure 8 shows the noisy intensity images in the $z=19\mu m$plane. The SNRs for the intensity images in the four cases are respectively 55db, 39db, 44db and 38db. The curves for the averaged NMSEs and the retrieved phases for the four cases are respectively shown in Fig. 9 and Fig. 10 . The conclusions are all in accordance with those obtained in Sections 4.1 and 4.2. When $M=1$, the RMSEs for the four retrieved phases are respectively 0.231, 1.366, 0.6271 and 1.207. While the RMSEs with $M=15$ are much smaller, which are 0.075, 0.677, 0.263 and 0.802.

## 5. Conclusions

We presented a method for estimating the axial intensity derivative in the TIE from multiple noisy unequally-spaced intensity measurements. An important characteristic of the method is that it could minimize the noise effect on the retrieved phase by using multiple unequally-spaced intensities and this has been proven by the mathematical error analysis. When multiple intensities corrupted with Gaussian noise, Poisson noise and mixed noise, the simulation results illustrate that the accuracy of the retrieved phase is advanced by our method.

## Appendix A

The propagation of the wave function ${u}_{z}(r)$ in the $z\ne 0$plane from the $z=0$plane can be calculated by the angular spectrum method:

where $F\{\u2022\}$ and ${F}^{-1}\{\u2022\}$ denote respectively the Fourier and inverse Fourier transforms. $g=({g}_{x},{g}_{y})$is the spatial frequency, conjugate to $r$. So, the analytical expression of the first-order intensity derivatives with respect to $z$ can be obtained as:

Let $z=0$, Eq. (A2) turns to:

## Acknowledgments

This work is partly supported by a grant from the National Natural Science Foundation of China (No.60502018, 61171005), the Innovation Foundations of BUAA for PhD Students, the Fundamental Research Funds for the Central Universities and the Aeronautical Special Founds of AVIC(No.CXY2010BH02).

## References and links

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

**2. **F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt. **27**(7), 1223–1225 (1988). [CrossRef] [PubMed]

**3. **T. E. Gureyev, A. Roberts, and K. A. Nugent, “Partially coherent fields, the transport-of-intensity equation, and phase uniqueness,” J. Opt. Soc. Am. A **12**(9), 1942–1946 (1995). [CrossRef]

**4. **V. V. Volkov, Y. Zhu, and M. De Graef, “A new symmetrized solution for phase retrieval using the transport of intensity equation,” Micron **33**(5), 411–416 (2002). [CrossRef] [PubMed]

**5. **J. Frank, S. Altmeyer, and G. Wernicke, “Non-interferometric, non-iterative phase retrieval by Green’s functions,” J. Opt. Soc. Am. A **27**(10), 2244–2251 (2010). [CrossRef] [PubMed]

**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. **12**(9), 1932–1941 (1995). [CrossRef]

**7. **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]

**8. **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]

**9. **S. V. Pinhasi, R. Alimi, L. Perelmutter, and S. Eliezer, “Topography retrieval using different solutions of the transport intensity equation,” J. Opt. Soc. Am. A **27**(10), 2285–2292 (2010). [CrossRef] [PubMed]

**10. **C. Roddier and F. Roddier, “Wave-front reconstruction from defocused images and the testing of ground-based optical telescopes,” J. Opt. Soc. Am. A **10**(11), 2277–2287 (1993). [CrossRef]

**11. **A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett. **23**(11), 817–819 (1998). [CrossRef] [PubMed]

**12. **J. Frank, J. Matrisch, J. Horstmann, S. Altmeyer, and G. Wernicke, “Refractive index determination of transparent samples by noniterative phase retrieval,” Appl. Opt. **50**(4), 427–433 (2011). [CrossRef] [PubMed]

**13. **K. A. Nugent, T. E. Gureyev, D. J. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard X rays,” Phys. Rev. Lett. **77**(14), 2961–2964 (1996). [CrossRef] [PubMed]

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

**15. **D. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent, “Quantitative phase-amplitude microscopy. III. The effects of noise,” J. Microsc. **214**(1), 51–61 (2004). [CrossRef] [PubMed]

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

**17. **W. X. Cong and G. Wang, “Higher-order phase shift reconstruction approach,” Med. Phys. **37**(10), 5238–5242 (2010). [CrossRef] [PubMed]

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

**19. **L. Waller, M. Tsang, S. Ponda, S. Y. Yang, and G. Barbastathis, “Phase and amplitude imaging from noisy images by Kalman filtering,” Opt. Express **19**(3), 2805–2814 (2011). [CrossRef] [PubMed]

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

**21. **J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996), pp. 55–61.