## Abstract

The transport of intensity equation (TIE) has long been recognized as a *quantitative* method for phase retrieval and phase contrast imaging. However, it is shown that the most widely accepted fast Fourier transform (FFT) based solutions do *not* provide an exact solution to the TIE in general. The root of the problem lies in the so-called “Teague’s assumption” that the transverse flux is considered to be a conservative field, which cannot be satisfied for a general object. In this work, we present the theoretical analysis of the phase discrepancy owing to the Teague’s assumption, and derive the necessary and sufficient conditions for the FFT-based solution to coincide with the exact phase. An iterative algorithm is then proposed aiming to compensate such phase discrepancy in a simple yet effective manner.

© 2014 Optical Society of America

## 1. Introduction

Transport of intensity equation (TIE) is a powerful tool for phase retrieval due to its non-interferometric, deterministic nature [1]. This method has recently received revived interest, owing to its potential to provide wavefront metrology with simple and inexpensive setups, without requiring external interferometric reference beam. By solving the TIE, the quantitative phase can be directly recovered with only a minimum of two intensity measurements at distinct planes orthogonal to the optical axis. The TIE relaxes the stringent beam-coherence requirements for interferometry, making it possible to realize quantitative phase imaging over a wide range of light- and electron-beams [2, 3]. In the past few decades, the TIE has found a variety of applications in adaptive optics [4], X-ray diffraction [5], electron microscopy [6], optical quantitative phase imaging [7] and computational light field imaging [8, 9]. Several optical configurations for realizing dynamic TIE phase imaging have also been developed with use of the chromatic aberrations [10], tunable lens [11], and spatial light modulator [12].

Despite these great improvements and extensive applications, the literature survey conducted by the authors reveals that one important issue regarding the solution of the TIE has not yet been well addressed in the practice of the TIE-based methods. The well-posedness and the uniqueness of the TIE solution are well studied in the literature - with proper boundary conditions, one may find a unique (apart from an arbitrary additive constant) solution to the TIE if the intensity distribution is strictly positive [13]. Nevertheless, the standard procedure for solving the TIE includes the introduction of the Teague’s auxiliary function [1] that transforms the TIE into a classical two-dimensional (2D) Poisson equation [14, 15]. This treatment makes the TIE conveniently solvable by Fourier transform-based methods for periodic or some simplified homogeneous boundary conditions [2, 16], or the discrete cosine transform based method for inhomogeneous Neumann boundary conditions [15]. However, it is important to remark that the Teague’s auxiliary function does not always exist in practical situations since the transverse energy flux may not be conservative, and consequently it would produce results that would not adequately match the exact solution [17, 18]. This discrepancy has not yet been recognized as a significant problem since previously the TIE was mainly applied for weak absorbing specimens, like optical components and biological samples. Although in such cases the methods using Teague’s auxiliary function are not exactly accurate, they provide a very good approximation for the exact solution to the TIE because the transverse energy flux is nearly irrotational. However, as will be shown in this paper, the phase error is often nonnegligible and may even be relatively large when the measured sample exhibits a strong absorption. This obviously poses a significant problem in the practice of the TIE phase retrieval, which has been generally acknowledged to be quantitative. To this end, this work aims to present the theoretical analysis of the phase discrepancy owing to the Teague’s auxiliary function as well as a simple iterative algorithm to compensate such phase errors. It addresses a long-lasting issue regarding the foundations of the widely applied TIE-based phase retrieval methods and extends the scope of their applicability.

## 2. Transport of intensity equation, partial coherence, and geometrical vector flux

The TIE was originally derived by Teague [1] from the stationary paraxial wave equation with the assumption of a purely coherent field. Its validity for partially coherent fields was then explored and confirmed by properly extending the physical meaning of phase, connecting it to the current density or the energy flow (Poynting vector) [2, 9]. Without loss of generality, we consider a quasi-monochromatic, stationary and ergodic scalar field $U\left(x\right)$, where $x$ is the 2D spatial vector. To characterize such statistical field, we introduce the Wigner distribution function (WDF) in 4D phase space $\left(x,u\right)$ [19], defined as

Though Eqs. (6) and (7) suggest that the transverse flux for a fully coherent field takes the form

one should always treat it with care because the $\nabla \varphi \left(x\right)$ remains undefined where the amplitude goes to zero. Since the zeros of intensity are usually associated with the presence of phase vortices, there is no way to define a $\varphi \left(x\right)$ so that it will be differentiable everywhere, which in turn leads to the fact that $\nabla \varphi \left(x\right)$ in Eq. (9) cannot be considered to be simply the gradient of some scalar potential. According to the Helmholtz decomposition theorem for vector fields [23], it is therefore more properly to modify the $\nabla \varphi \left(x\right)$ in Eq. (9) as the sum of the gradient of a scalar potential function ${\varphi}_{s}\left(x\right)$ and the curl of a vector potential function ${\varphi}_{v}\left(x\right)$*z*. Therefore, as the wave propagates, both of the two kinds of phase will manifest themselves, contributing to the longitudinal intensity variations. However, the vortex cores are associated with intensity zeros that break the uniqueness theorems of the TIE [13], making conventional TIE solvers fail in such scenario [2]. Although some recent work extends the applicability of the TIE to arbitrary phase functions by explicitly taking into account the effect of the vorticity [24, 25], using intensity zeros to identify singularities tends to be difficult since

*very-near*zero and

*at*zero may seem almost the same but can be entirely different

*vis-a-vis*the presence of a phase singularity point. In the remainder of this paper, we exclude the consideration of the vector phase ($\varphi \left(x\right)={\varphi}_{s}\left(x\right)$) by assuming

*a priori*that the phase is single-valued and continuous. This assumption is generally reasonable since the TIE phase measurements are habitually performed directly in (or very close to) the object plane (or its conjugated plane), where the phase function is single-valued and well-defined by the projected optical path length and the refractive index of the sample.

## 3. Theoretical analysis of the phase discrepancy

Putting the phase vorticity aside, we now turn to the problem of phase discrepancy in well-acknowledged TIE solvers. In the literature, the TIE is usually solved under the so-called “Teague’s assumption” that the transverse flux is conservative so that an a scalar potential (auxiliary function) $\psi $ exists that satisfying

Despite its widespread use, the “solution” given by Paganin and Nugent’s method [Eq. (14)] is only an approximation to the exact solution, and there has been some research regarding its validity and accuracy [17, 18]. In this work, we refer this error as the phase discrepancy and take it to be equal to the difference between the true phase $\varphi $ and the reconstructed phase ${\varphi}_{0}$ by Eq. (14) [corresponding to Eq. (13) of [17]]. It should be plain to see that the phase discrepancy originates from the Teague’s assumption, since as we know from the Helmholtz decomposition theorem that the transverse flux should be decomposed in terms of the gradient of a scalar potential $\psi $ and the curl of a vector potential $\eta $

Next we examine how this missing rotational term will affect the phase reconstruction. When inserting the correct decomposition Eq. (15) into the TIE [Eq. (6)] we obtain the exact same equation as Eq. (12) because the rotational term $\nabla \times \eta \left(x\right)$ is divergence free. Consequently, the curl-free component of the transverse flux ${j}_{x}\left(x\right)$ can be reconstructed by solving this Poisson equation. The second Poisson equation, however, is different from Eq. (13), obtained by moving the $I\left(x\right)$ to the RHS of Eq. (15) and then taking the divergence on both sides

There are a few special cases with clear physical meanings when the above conditions are established. The first two trivial cases are uniform intensity (Teague’s assumption is unnecessary) and constant phase (physically trivial). Another interesting special case is that $\nabla I\left(x\right)\times \nabla \varphi \left(x\right)=0$, which is also the necessary and sufficient conditions for the validity of Teague’s assumption [since $\eta \left(x\right)$ is defined by Eq. (17)]. This means $\nabla I\left(x\right)$ and $\nabla \varphi \left(x\right)$ must be parallel, or equivalently, $\nabla \varphi \left(x\right)=c\left(x\right)\nabla I\left(x\right)$ for some scalar function $c\left(x\right)$. When $c\left(x\right)$ is taken as $\gamma {I}^{-1}\left(x\right)$ such that $\nabla \varphi \left(x\right)=\gamma \nabla \mathrm{ln}I\left(x\right)$, we have$\varphi \left(x\right)=\gamma \mathrm{ln}I\left(x\right)$. This is an equivalent definition of the Lambert–Beer's law of absorption [26] which states there is a logarithmic dependence between the transmission of light and the density or path length through the material. Due to the inherent correlation between the intensity and phase, we could usually expect $\nabla I\left(x\right)$ and $\nabla \varphi \left(x\right)$ to be nearly parallel for a typical transmissive object. More physical interpretations and examples regarding the validity of Teague’s assumption can be found in [17]. However, generally speaking, the condition given by Eq. (20) cannot be satisfied for an arbitrary intensity and phase distribution, and therefore the Paganin and Nugent’s method would not produce results that adequately match the exact solution of TIE.

The effects of the phase discrepancy have not yet been recognized as a significant problem since previously the TIE was mainly applied for weak absorbing specimens, like optical components, biological soft tissues, and cells. For these kind of objects, the phase discrepancy is negligible [second term on the RHS of Eq. (18) is close to zero] because the total intensity variation is small ${\Vert {\scriptscriptstyle \frac{\nabla I\left(x\right)}{I\left(x\right)}}\Vert}_{2}\approx 0$, where ${\Vert \cdot \Vert}_{2}$is the normalized *L _{2}* norm, defined as

However, for a strongly absorbing specimen (like a reflective object, whose intensity and phase are generally uncorrelated), the effects of the phase discrepancy can be easily observed. Figure 1 shows a realistic example of the error arising from a simulated reflective object (a piezo-attenuated circular micro-membrane diaphragm), reconstructed using the Paganin and Nugent’s method. The phase and intensity of the object is shown in Figs. 1(a) and 1(b) respectively. The minimum value of the intensity is only 0.05 due to the strong absorption of the object. The simulation is designed to replicate the real measurement from the data obtained by a digital holographic microscope. Figure 1(c) shows the retrieved phase by the Paganin and Nugent’s solution. Note the phase recovered by the TIE is continuous without 2π jumps, and the phase shown in Fig. 1(c) is digitally rewrapped into the interval [-π, π) for better comparison. The phase discrepancy (the difference between the ground truth phase [Fig. 1(a)] with the reconstructed phase [Fig. 1(c)] is shown in Fig. 1(d). The relative root-mean-square error (RMSE) of the reconstruction is 6.76%. This example obviously indicates that the phase discrepancy can be a practical problem when dealing with absorbing objects with use of the TIE, which has been generally acknowledged to be quantitative.

## 4. Iterative algorithm for phase discrepancy compensation

We now proceed to develop a simple Picard-type iterative algorithm [27] to compensate the phase discrepancy owing to Teague’s assumption. Assuming we have already got an inaccurate solution ${\varphi}_{0}\left(x\right)$ by Paganin and Nugent’s method, the Helmholtz decomposition theorem states that

To verify the convergence of the proposed method, it is sufficient to prove that the sequence of functions $\Delta {\varphi}_{n}$ converges uniformly to zero. In practice, the rigorous proof of the convergence is not possible without some *a priori* information. Here we will take a conservative stance and assume that the Paganin and Nugent’s method is a *valid* (though not accurate) solver for the TIE, which means the phase discrepancy should *at least smaller than* the exact solution in norm: ${\Vert \varphi -{\varphi}_{0}\Vert}_{2}\le L{\Vert \varphi \Vert}_{2},$with constant$\text{}L1$. By analogy, we can derive ${\Vert \varphi -{\varphi}_{n+1}\Vert}_{2}\le L{\Vert \varphi -{\varphi}_{n}\Vert}_{2}$, which means the phase discrepancy can only decrease at each iteration. Besides, since ${\Vert \varphi -{\varphi}_{n+1}\Vert}_{2}\le {L}^{n}{\Vert \varphi -{\varphi}_{0}\Vert}_{2}$, and ${L}^{n}\to 0$as $n\to \infty $, the compensated phase function is shown to be a Cauchy sequence and thus it converges to the exact solution. In the ideal noise-free condition, the above iterative algorithm usually converges very rapidly and the error decreases exponentially with the iterations. Within two to four iterations the phase discrepancy can be reduced to a negligible level and the exact solution to the TIE can be thus obtained. However, in a real measurement the correction term can never reach a very low level because the noise, finite sampling effect, and the nonlinearity error in the finite difference estimation of the axial intensity derivative will introduce additional errors. Thus, in practice, the iteration should stop when there is no significant decrease in the RMS of the correction term $\Delta {\varphi}_{n}$ as the iteration number increases. Besides, the maximum iteration number can also be set as an additional termination condition to prevent unnecessary ineffective iterations.

## 5. Simulation results

#### 5.1 Smooth wavefront

Several simulations were performed to test the proposed method. First, we adopted the “counter example” suggested in [17] to verify whether the proposed approach can effectively compensate the phase discrepancy even with sufficiently large scale. In the simulation, the complex field was defined on an 256 × 256 grid with 0.025μm × 0.025μm pixel size. The wavelength is 632.8nm (HeNe Laser). The intensity derivative $\partial I/\partial z$ is estimated by central finite difference between two defocused images ${I}_{+}$ and ${I}_{-}$ with the distance of $\Delta z=\pm 10$μm generated by the angular spectrum numerical propagation. As shown in Fig. 3, the phase and intensity functions are defined as

andwhere${a}_{0}={b}_{0}={10}^{-2}{\text{\mu m}}^{-2}$ and ${a}_{1}={b}_{2}={0.25}^{-2}\times {\text{10}}^{-8}{\text{\mu m}}^{-8}$. The metric RMSE was used to quantify phase retrieval accuracy. Since the TIE can only determine the phase uniquely up to an additive constant, the piston term (the mean value) of each result was removed before the RMSE evaluation.To illustrated the effects of the Teague’s assumption in the phase reconstruction, the
Helmholtz decomposition was performed based on the transverse flux field
$I\nabla \varphi $ (neglecting the constant –*k*) and the
results are shown in Fig. 4.The butterfly-like transverse flux field can be decomposed into a gradient field and a
rotational field. The Teague’s assumption simply neglects the rotational field, which
obviously is comparable in norm to the gradient field. As a result, a large phase discrepancy
(RMSE 9.53%) appeared when the Paganin and Nugent’s solution was used [Eq. (14)], as demonstrated in Fig. 5(a).With the proposed compensation process after the first iteration, the phase RMSE reduced
significantly to around 0.99%, as shown in Fig. 5(b).
The error was even down to a negligible level - around 0.016% after the third iteration, as
shown in Figs. 5(c). The evolution of the RMSE values
versus the iteration times is shown as one-line profile in Fig.
5(d), which clearly illustrates the rapid convergence of the proposed approach.

In a practical TIE measurement system, noise exists on each captured intensity images. To test
the performance of the proposed method in the presence of noise, different levels of
pseudo-random Gaussian noise with standard deviations (STD) of 10^{−5},
10^{−4}, and 10^{−3} was generated and superimposed on each
intensity image. In such scenario, the phase error was not only resulted from the
Teague’s assumption but also from the influence of the noise. For the low-noise-level
case, the underlying saddle-like structure in the phase error of the Paganin and
Nugent’s solution before compensation can be clearly perceived [Fig. 6(a)].When our method was applied, the majority amount of phase discrepancy could be
compensated after three iterations, with the residual phase error shown in Fig. 6(b). The iterative process successfully removed the large waviness due
to the phase discrepancy, leaving only the cloud-like random errors [28, 29] from the intensity noise.
However, with the increase of the noise level, the noise effect became dominant. The initial
phase errors were almost submerged by the cloud-like artifacts, showing no evident sign of the
phase discrepancy (Fig. 7).Their corresponding RMSE curves only decreased slightly and could never go down to very
small values due to the nonnegligible noise.

To verify whether our method can effectively remove the phase discrepancy in the presence of significant noise, we also show the correction term by our method (which is the difference between the initial solution and the final refined solution after 4 iterations) in Fig. 7. For the medium noise level, the correction term closely resembles the phase discrepancy in the noise-free case. While for the high noise level, the correction term was distorted and no longer matches the noise-free phase discrepancy, which can be expected since the noise has already changed the basic structure of intensity distribution fundamentally. These results clearly confirm the validity and the effectiveness of the compensation algorithm in the presence of noise.

#### 5.2 Complex phase and intensity distribution

In the previous simulation, the complex field was constructed with very simple and smooth
mathematical functions. Next, we focus on more complex phase and intensity distributions. We
consider the complex intensity distribution of letters ‘COLE’ (Fig. 8(a)) and the phase profile of letters
‘TIE’ (Fig. 3(b)). The value range of the intensity is [0.05, 1]. For such complex object with a lot of high
spatial frequency details, the effect of the finite sampling in both spatial
(*x-y*) and axial (*z*) dimensions must be given additional
attentions.

For the *x-y* dimension, the Shannon theorem indicates that for the correct
acquisition of a given 2D signal, the sample frequency must be at least twice as high as the
highest signal frequency of interest [30]. For the axial
direction, the situation is a bit more complicated since the axial intensity derivative
$\partial I/\partial z$ in TIE is calculated via finite differences from two (or more)
defocused intensities. When the separation between the two measurement planes was sufficiently
small ($\Delta z=\pm 10$μm), the effect of the finite sampling was not evident.
Figure 8(c) shows the phase discrepancy of the Paganin
and Nugent’s solution before compensation (RMSE 15.91%). After five iterations of
compensation, the majority amount of phase discrepancy was removed, with the residual phase
error shown in Fig. 8(d). However, for relatively larger
$\Delta z$ (usually required for high signal-to-noise measurements [29, 31, 32]), the breakdown of the linear approximation that
underlies the finite difference approximation induces nonlinear errors, and the estimate
$\partial I/\partial z$ does not accurately represent the axial intensity derivative.
Typically, this results in a loss of resolution of the retrieved phase map. It can be clearly
seen from Fig. 9 that object edge information (high
frequency details) appears in the phase error maps when a larger defocus distance is used
($\Delta z=\pm 100$ or $\pm 500$μm). When our compensation method was applied, phase discrepancies were effectively
compensated, leaving only the edge error resulting from the nonlinearity errors. In all the
three cases, the compensation procedure converges after 4 to 5 iterations and does not show
significant improvement with more subsequent iterations. These results verify the stability of
our algorithm under large defocus distances and nonlinearity errors.

#### 5.3 A realistic example and some practical considerations

Finally, let us turn to the realistic MEMS sample shown in Fig. 1. The value range of the intensity is [0.05, 1] and the RMSE of the Paganin and Nugent’s solution before compensation is 6.76%. Once again, our approach converges to the exact solution (blue curve) with very small residual error (RMSE 0.119%), and the final compensated phase closely matches the ground-truth image (Fig. 10).

Before concluding, several important issues in the practice of our approach need to be clarified. Although no rigorous proof of convergence of our approach has been devised, the process did converge in all above cases. This seems reasonable because the iterative process is essentially self-correcting. In the presence of noise and nonlinearity error, our algorithm effectively extracts and removes the phase discrepancy without noticeable disturbance by other types of errors. However, as we mentioned at the end of Section 4, the stability and convergence of our method depend heavily on the accuracy of the Paganin and Nugent’s solution. Since the intensity $I$ appears in the denominator [Eq. (14)], the Paganin and Nugent’s solution becomes quite unstable when the minimum intensity is close to zero. Therefore, caution is required when very small values appear in the intensity distribution. On the other hand, these small intensity values also increase the phase discrepancy accordingly since the total intensity variation ${\Vert {\scriptscriptstyle \frac{\nabla I\left(x\right)}{I\left(x\right)}}\Vert}_{2}$ can become much larger. For example, if we rescale the intensity range of the MEMS sample to [0.001, 1], the iterative process will be out of control after two iterations (red curve). The singular point creates significant artifacts that degrade the whole phase reconstruction prevailingly. In this limiting case, we should avoid dividing too small values to guarantee the stability of solution by setting a certain intensity threshold (e.g. 1% of the maximum value of $I$) when evaluating the term in the square brackets in Eq. (14). When the intensity threshold (0.01) is used, our method becomes convergent again and finally reaches a reasonable accuracy (RMSE 2.02%).

## 5. Conclusions

To conclude, this work has presented the theoretical analysis of the phase discrepancy in the most widely applied FFT-based solution of the TIE owing to the Teague’s assumption. The necessary and sufficient conditions for the solution to coincide with the exact phase have been derived. We also have proposed a simple and novel iterative algorithm to compensate the phase discrepancy due to the Teague’s assumption when such conditions cannot be satisfied for a general object. Several numerical results have been presented showing the effectiveness of the proposed algorithm.

A clear advantage of the proposed method resides in the very rapid convergence of the iterative phase corrections. The phase error drops exponentially with iterations and shows no sign of stagnation for a wide variety of test cases. Moreover, the convergence is insensitive to the intensity noise and the separation between the input image planes. Probably the most limiting characteristic of the method is its instability to the intensity singularity, a feature inherited from the FFT-based TIE solvers. This problem can be countered by setting a certain intensity threshold. Since generally only 2 to 4 iterations are sufficient to compensate most phase discrepancy and thanks to the efficient FFT-based operations, the iterative nature of our method does not preclude its application to high-speed, real-time phase imaging. It should be noted that although the phase discrepancy problem can be avoided by using some other TIE solvers that do not involve the use of Teague’s assumption, like the multigrid method [14, 33], our approach can be still considered as an attractive solution for its better spatial resolution (FFT-based solvers do not use finite difference to approximate the spatial derivatives) [34], simplified treatment for boundary conditions [15], and easy implementation.

## Acknowledgments

This project was supported by the Research Fund for the Doctoral Program of Ministry of Education of China (No. 20123219110016), and the National Natural Science Foundation of China (No. 61271332). C. Zuo acknowledges the financial support from China Scholarship Council (No. 201206840009).

## 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. **A. M. Zysk, R. W. Schoonover, P. S. Carney, and M. A. Anastasio, “Transport of intensity and spectrum for partially coherent fields,” Opt. Lett. **35**(13), 2239–2241 (2010). [CrossRef] [PubMed]

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

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

**6. **K. Ishizuka and B. Allman, “Phase measurement of atomic resolution image using transport of intensity equation,” J. Electron Microsc. (Tokyo) **54**(3), 191–197 (2005). [CrossRef] [PubMed]

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

**8. **A. Orth and K. B. Crozier, “Light field moment imaging,” Opt. Lett. **38**(15), 2666–2668 (2013). [CrossRef] [PubMed]

**9. **C. Zuo, Q. Chen, and A. Asundi, “Light field moment imaging: comment,” Opt. Lett. **39**(3), 654 (2014). [CrossRef] [PubMed]

**10. **L. Waller, S. S. Kou, C. J. R. Sheppard, and G. Barbastathis, “Phase from chromatic aberrations,” Opt. Express **18**(22), 22817–22825 (2010). [CrossRef] [PubMed]

**11. **C. Zuo, Q. Chen, W. Qu, and A. Asundi, “High-speed transport-of-intensity phase microscopy with an electrically tunable lens,” Opt. Express **21**(20), 24060–24075 (2013). [CrossRef] [PubMed]

**12. **C. Zuo, Q. Chen, W. Qu, and A. Asundi, “Noninterferometric single-shot quantitative phase microscopy,” Opt. Lett. **38**(18), 3538–3541 (2013). [CrossRef] [PubMed]

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

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

**15. **C. Zuo, Q. Chen, and A. Asundi, “Boundary-artifact-free phase retrieval with the transport of intensity equation: fast solution with use of discrete cosine transform,” Opt. Express **22**(8), 9220–9244 (2014). [CrossRef] [PubMed]

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

**17. **J. A. Schmalz, T. E. Gureyev, D. M. Paganin, and K. M. Pavlov, “Phase retrieval using radiation and matter-wave fields: Validity of Teague's method for solution of the transport-of-intensity equation,” Phys. Rev. A **84**(2), 023808 (2011). [CrossRef]

**18. **J. A. Ferrari, G. A. Ayubi, J. L. Flores, and C. D. Perciante, “Transport of intensity equation: Validity limits of the usually accepted solution,” Opt. Commun. **318**, 133–136 (2014). [CrossRef]

**19. **M. Testorf, B. Hennelly, and J. Ojeda-Castañeda, *Phase-Space Optics: Fundamentals and Applications: Fundamentals and Applications: Fundamentals and Applications* (McGraw Hill Professional, 2009).

**20. **R. Winston and W. T. Welford, “Geometrical vector flux and some new nonimaging concentrators,” J. Opt. Soc. Am. **69**(4), 532–536 (1979). [CrossRef]

**21. **M. J. Bastiaans, “Application of the Wigner distribution function to partially coherent light,” J. Opt. Soc. Am. A **3**(8), 1227–1238 (1986). [CrossRef]

**22. **B. Boashash, “Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals,” P IEEE **80**(4), 520–538 (1992). [CrossRef]

**23. **H. Jeffreys and B. Jeffreys, *Methods of Mathematical Physics* (Cambridge University, 1999).

**24. **L. J. Allen, H. M. Faulkner, K. A. Nugent, M. P. Oxley, and D. Paganin, “Phase retrieval from images in the presence of first-order vortices,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **63**(3), 037602 (2001). [CrossRef] [PubMed]

**25. **A. Lubk, G. Guzzinati, F. Börrnert, and J. Verbeeck, “Transport of Intensity Phase Retrieval of Arbitrary Wave Fields Including Vortices,” Phys. Rev. Lett. **111**(17), 173902 (2013). [CrossRef] [PubMed]

**26. **S. R. Crouch and J. D. Ingle, *Spectrochemical Analysis* (Prentice Hall, 1988).

**27. **V. Berinde, “Approximating fixed points of weak contractions using the Picard iteration,” in *Nonlinear Analysis Forum* (2004), 43–54.

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

**29. **C. Zuo, Q. Chen, Y. Yu, and A. Asundi, “Transport-of-intensity phase imaging using Savitzky-Golay differentiation filter--theory and applications,” Opt. Express **21**(5), 5346–5362 (2013). [CrossRef] [PubMed]

**30. **A. J. Jerri, “The Shannon sampling theorem—Its various extensions and applications: A tutorial review,” P IEEE **65**(11), 1565–1596 (1977). [CrossRef]

**31. **J. Martinez-Carranza, K. Falaggis, and T. Kozacki, “Optimum measurement criteria for the axial derivative intensity used in transport of intensity-equation-based solvers,” Opt. Lett. **39**(2), 182–185 (2014). [CrossRef] [PubMed]

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

**33. **T. Gureyev, C. Raven, A. Snigirev, I. Snigireva, and S. Wilkins, “Hard x-ray quantitative non-interferometric phase-contrast microscopy,” J. Phys. D Appl. Phys. **32**(5), 563–567 (1999). [CrossRef]

**34. **J. Martinez-Carranza, K. Falaggis, T. Kozacki, and M. Kujawinska, “Effect of imposed boundary conditions on the accuracy of transport of intensity equation based solvers,” Proc. SPIE **8789**, 87890N (2013).