## Abstract

We investigate the measurement of a thin sample’s optical thickness using the transport of intensity equation (TIE) and demonstrate a version of the TIE, valid for partially coherent illumination, that allows the measurement of a sample’s optical path length by the removal of illumination effects.

© 2013 OSA

## 1. Introduction

The space–varying component of a monochromatic, coherent, scalar field *U* may be represented in terms of its spatial distribution of intensity *I*, which is directly measurable by a detector, and its phase *ϕ*, which is not. Phase retrieval is of interest, particularly in imaging systems, because the intensity image carries information about the attenuation of light through the sample, while phase carries information about the optical path length through the sample. In microscopy of unstained biological samples or wavefront characterization, for example, much of the useful information may be contained in phase. Techniques to quantitatively recover the phase of a coherent field are often based upon interferometry of the field of interest with a second field, rendering the phase differences between the two fields visible as an intensity modulation. An alternative set of techniques relies on the interference of the field with itself upon propagation. One such propagation–based method relies upon guessing the initial unknown phase, and uses repeated iterations of computational propagation between two (or more) planes and the enforcement of prior knowledge about the field in order to reduce the error in the estimated phase [1]. Iterative solutions are widely used, but can require many iterations to converge, and do not always converge to the correct phase. Here we investigate a different class of propagation–based techniques, which rely on intensity measurements in closely spaced planes to provide a linear equation relating those intensities and the underlying phase. Phase recovery can be formulated as the inversion of the transport of intensity equation (TIE) [2].

Consider a paraxial, monochromatic, coherent beam, propagating along the *z* axis. This field may be represented over space in the complex representation, *U*:

**x**represents position transverse to the optical axis and

*I*and

*ϕ*are real–valued intensity and spatially–varying phase. The physics of paraxial propagation under this representation are governed by the paraxial wave equation:

*D*Laplacian operator over

**x**.

An alternative formulation of Eq. (3) in terms of propagation of the real–valued intensity *I* and phase *ϕ* may be derived by applying Eq (3) to Eq. (2), and separating the real and imaginary parts. The cost of eliminating the complex–valued field *U* is that propagation is described by a pair of coupled, nonlinear differential equations:

**is the 2D gradient operator over**

_{x}**x**, ⊗ denotes the dyadic product between two vectors, and the notation

**v**(

_{x}**x**,

*z*) = ∇

_{x}*ϕ*(

**x**,

*z*)/

*k*has been adopted for notational simplicity. This pair of equations represents a coupled set of continuity equations expressing the conservation of intensity

*I*and transverse flux density

*I*

**v**upon propagation, and were initially introduced in quantum mechanics [3, 4] [where the probability density and phase of a particle’s wavefunction satisfy analogous forms to Eqs. (4)]. Physically, these equations express the propagation of intensity along flux lines [Eq. (4a)] (flux being the analogue of Poynting vector for scalar fields) [5], as well as the propagation of the flux lines themselves [Eq. (4b)]. The term in parentheses on the right–hand side of Eq. (4b) represents the so–called “quantum potential,” which becomes important near caustics of the field and prevents these flux lines from crossing and represents a diffraction term in optics. If this term is neglected, the second equation represents the eikonal equation of geometrical optics. By substituting

*k*

**v**= ∇

_{x}*ϕ*, Eq. (4b) can be simplified to a phase rather than flux density transport equation,

*∂ϕ/∂z*, ([see, for example, Eq. (3) of Ref. [6]).

In TIE imaging, one measures changes of intensity using a detector that can be displaced along the optical axis and considers regions away from caustics and propagation distances small enough that Eq. (4b) can be ignored. Moreover, for a coherent field **v**(**x**, *z*) = ∇_{x}*ϕ*(**x**, *z*)/*k* so that Eq. (4a) becomes the TIE:

*ϕ*to changes in intensity

*∂I/∂z*through a single, linear differential equation. The derivative along

*z*can be approximated by finite differences between intensity measurements at two or more closely spaced planes [7–10], allowing this equation to be computationally inverted for phase. TIE–based approaches are of practical interest because of their ease of implementation in optical imaging systems as well as in cases where optical elements are difficult to manufacture, such as X–ray [11, 12], electron [13] and neutron beam imaging [14].

In many cases, the concern is recovering the optical properties of a sample being illuminated rather than the phase of the incident field. Consider a thin sample, located without loss of generality at the plane *z* = 0, characterized through a transmission function

*T*and

*ψ*are real–valued. Under the thin–sample approximation, the total phase accrued upon propagation through the object is equal to the product of the optial path length (OPL) through the object with the wave number of the illumination,

*ψ*=

*k*OPL. If the sample consists of a uniform material of known index of refraction

*n*, its thickness is given by OPL/

*n*. When the sample is illuminated coherently, the complex field immediately after passing through the sample is given by $U\left(\mathbf{x},0\right)=\sqrt{T\left(\mathbf{x}\right)I\left(\mathbf{x},0\right)}\text{exp}\left\{\text{i}\left[\varphi \left(\mathbf{x},0\right)+\psi \left(\mathbf{x}\right)\right]\right\}$. If the object plane is relayed onto a detector,

*∂*|

*U*|

^{2}/

*∂z*can be estimated by defocusing the detector, and the total phase

*ϕ*(

**x**, 0) +

*ψ*(

**x**) can be retrieved from the TIE. The sample’s OPL can be recovered if

*ϕ*(

**x**, 0) is known. In this manuscript, we consider the case in which partially coherent light illuminates a thin sample and show that although a partially coherent field does not have a well–defined 2D phase in the sense of Eq. (2), it is possible to characterize the incident illumination such that a TIE–type measurement yields

*ψ*for a thin sample.

## 2. The TIE under illumination of any state of coherence

While a monochromatic, coherent, paraxial field may be represented in terms of a complex–valued function over a plane, a partially coherent field contains statistical fluctuations and is characterized by second–order correlations between pairs of points in that plane [15]. For simplicity we consider a paraxial field of central frequency whose second–order correlations are described in the frequency domain by the cross-spectral density, *W*(**x**_{1}, **x**_{2}, *z*, *ν*), over a plane of constant *z*, where *ν* denotes frequency. We further restrict ourselves to quasi-monochromatic fields or a single frequency component of a polychromatic field and suppress explicit reference to *ω* in what follows. The cross–spectral density may be represented as *W*(**x**_{1}, **x**_{2}, *z*) = 〈*U*^{*}(**x**_{1}, *z*)*U*(**x**_{2}, *z*)〉, where 〈·〉 denotes the ensemble average over a statistical ensemble of realiziations of the field, {*U*(**x**, *z*)}, each of which satisfies a paraxial wave equation, Eq. (3). The spectral density, *S*(**x**, *z*) = *W*(**x**, **x**, *z*), represents the power incident upon a detector placed perpendicular to the *z* axis, such that *S* is an analogue to the coherent intensity for a quasimonochromatic field. However, for fixed *z*, *W* is complex–valued and defined over the 4D space spanned by the pair of planar 2D position vectors, **x**_{1,2}. The phase of *W* is therefore not a 2D function of a local variable, but rather describes the correlations of the field at all pairs of points in that plane. One can still measure a focal stack of spectral densities and apply Eq. (5) in order to reconstruct a 2D function in this plane, although this “phase” clearly does not provide the phase of *W*. Paganin and Nugent have interpreted this as a scalar phase whose gradient is related to the normalized transverse energy flux vector of the field [5]. This interpretation of “phase” as a scalar–valued quantity describing local energy flow was further studied by Gureyev *et. al.* as a generalized eikonal [16] (although one should note that Gureyev’s generalized eikonal approximately satisfies the eikonal equation in the short wavelength limit only for highly coherent fields [17]). Recently, Zysk *et. al.* pointed out that since a partially coherent field may be decomposed into a sum of fully coherent modes, each of which is fully incoherent with respect to the other modes, TIE recovers a phase that is a weighted superposition of the phases of these modes [18]. However, the phases of the individual modes cannot be determined from a single TIE measurement. Finally Gureyev *et. al* have also considered case of a thin sample placed in a partially coherent beam, derived a result for the intensity after propagation, and suggested that the sample’s thickness could be characterized if the illumination were first measured in the absence of the sample [19]. In what follows, we consider the relationship between the quantity reconstructed by the TIE and the optical properties of a thin sample placed in the beam’s path and demonstrate that under small displacements of the detector, the effects of the illumination may be removed provided its flux vector is first characterized. It should be noted that if the goal is to reconstruct the 4D phase of *W*, a similar technique known as phase space tomography allows the reconstruction of this phase in an initial plane from a stack of many defocused intensity measurements (though since space is 3D, the extra measurement dimension is generally provided by introducing an astigmatic lens), but requires many more intensity measurements [20–22].

The cross–spectral density of a paraxial, partially coherent field satisfies a pair of paraxial wave equations in both transverse variables:

*j*= 1, the lower when

*j*= 2 and ${\nabla}_{{\mathbf{x}}_{j}}^{2}$ denotes a 2D Laplacian with respect to

**x**

*. By averaging Eqs. (7) for*

_{j}*j*= 1, 2, changing variables to mean and difference coordinates,

**x**

_{1,2}=

**x**∓

**x′**/2, and invoking the relationship ${\nabla}_{{\mathbf{x}}_{2}}^{2}-{\nabla}_{{\mathbf{x}}_{1}}^{2}=2{\nabla}_{\mathbf{x}}\cdot {\nabla}_{{\mathbf{x}}^{\prime}}$, a transport equation for

*W*is obtained:

*S*(

**x**,

*z*) =

*W*(

**x**,

**x**,

*z*), the partially coherent TIE (PC–TIE) can be defined by taking the limit of Eq. (8) as

**x′**→

**0**, yielding

**F**is related to

*W*by the expression

*S*over the transverse plane, where

**F**represents a transverse spectral density flux vector. In other words,

**F**expresses the fact that for spectral density to change at a given transverse position

**x**as the field propagates along the axis, energy from neighboring points must flow into or out of that point. This transverse flow of energy is accounted for by

**F**, the transverse energy flux density [23]. Assuming there are no points of null intensity present, it is possible to describe the transverse flux as the product of the spectral density and a transverse vector flow:

**F**(

**x**,

*z*) =

*S*(

**x**,

*z*)

**v**(

**x**,

*z*). It is further possible to express

**v**, through a Helmholtz decomposition, as the sum of a divergence–free component and a curl–free component. Since the TIE measures ∇

**·**

_{x}**F**, the divergence–free component is undetectable in this technique, so it is sufficient to express

**v**as a gradient,

**v**(

**x**,

*z*) = ∇

_{x}*ϕ*(

**x**,

*z*).

If *∂S/∂z* and *S*(**x**, *z*) are used in place of *∂I/∂z* and *I*(**x**, *z*) in the coherent TIE, Eq. (5), a meaningful scalar–valued phase may be reconstructed. The “phase,” *ϕ*, retrieved by using the TIE on a partially coherent field is related to the flux vector by

**v**represents the flow vector of spectral density. In other words, the phase recovered by the coherent TIE is a scalar–valued function whose gradient yields the transverse spectral density flow as the field propagates, and only in the case of a coherent field does this represent the phase of the complex–valued function describing the field itself. This is the same as the notion of the scalar phase introduced by Paganin and Nugent [5].

Except in certain extremely simplified states of coherence, knowledge of **F** and *S* is insufficient to characterize *W*. However, it is possible to obtain the optical thickness of samples using Eq. (9) under certain illumination conditions. Consider a field with cross–spectral density *W*_{inc}(**x**_{1}, **x**_{2}, *z*) and spectral density *S*_{inc}(**x**, *z*) passing through a thin sample described by Eq. (6). Immediately after the sample, the cross–spectral density is given by *τ*^{*}(**x**_{1})*τ*(**x**_{2})*W*_{inc}(**x**_{1}, **x**_{2}, *z*), with spectral density *S*_{tot} = *T*(**x**)*S*_{inc}(**x**, *z*) and the subscript “inc” denotes properties of the incident illumination. The PC–TIE for the total spectral density in the detector plane in this case is

*z*is infinitesimal. In reality, one is limited to estimating

*∂S/∂z*from finite displacements in

*z*, where the physics embodied by Eq. (7) can more accurately be represented by a pair of Fresnel propagation integrals. A more detailed derivation of the PC–TIE, examining the validitity of the TIE approximation used to reduce the propagation integrals to Eq. (12) is given in Appendix A.

Notice that if **F**_{inc} = **0**, *i.e.* the incident illumination is symmetric with respect to the optical axis, Eq. (12) simplifies to

**F**

_{inc}=

**0**include the case of coherent, plane wave illumination of a sample (for which

*S*is also uniform), and also the case of imaging about the waist of a partially coherent Gaussian Schell–model beam or using symmetric Köhler illumination.

The difference in form between Eqs. (12) and (13) is the presence of the term ∇** _{x}** · [

*T*(

**x**)

**F**

_{inc}(

**x**)]. This characterizes the intensity change along the axis due to interaction of the transmissivity

*T*, of the sample and flux of the incident illumination. The term on the right–hand side of Eq. (12) represents the flux induced by the sample,

*T*∇

_{x}*ψ*/

*k*, carrying the spectral density of the illumination. If

**F**

_{inc}≠

**0**, then

*ψ*may still be recovered provided that

**F**

_{inc}and

*S*

_{inc}, properties of the illumination, are known. This may be achieved by first measuring

*∂S*

_{inc}/

*∂z*with no sample in place, to characterize the scalar phase of

**F**

_{inc}, and then solving Eq. (12). Since the PC–TIE is valid for illumination in any state of coherence, it also includes the case of the coherent limit through the substitution

**F**

_{inc}(

**x**,

*z*) =

*k*

^{−1}

*I*

_{inc}(

**x**,

*z*)∇

_{x}*ϕ*

_{inc}(

**x**,

*z*), where

*I*

_{inc}and

*ϕ*

_{inc}are the intensity and phase of the coherent incident illumination in the sample plane, respectively.

Another special case of the PC–TIE is when the sample is nearly pure–phase, *i.e. T* ≈ *T*_{0} = constant, where Eq. (12) reduces to the simpler form

*∂S/∂z*measurements, with and without the sample in place.

## 3. Simulation results

Consider the arrangement shown in Fig. 1, assuming a monochromatic, planar source located one focal length from the first, collimating lens. Here we assume a perfect lens of focal length *f* and neglect the effect of finite apertures for simplicity. The field immediately after the collimating lens consists of plane waves from different points on the source. If the source is delta–correlated, which is a good approximation for LED or incandescent illumination passed through a diffuser, this system produces uniform Köhler illumination after the first lens. A single point on the source of intensity *I*_{0}, located at transverse position **x**_{0}, will produce a cross–spectral density of the form:

**d**is a result of the assumption of an idealized lens of infinite aperture, though Eq. (15) provides a good approximation within a small transverse region located close to the lens.

Since the source points are completely uncorrelated with respect to each other, each contribution of the form in Eq. (15) is a coherent mode of the field after the lens, whose weight is given by *I*_{0}(**x**_{0}), such that the total cross–spectral density is simply the integral over contributions from all points on the source, which corresponds to a Fourier transform over the source’s intensity distribution,

Here, we consider a planar, delta–correlated source that is a disk of radius *R* and uniform intensity *I*_{0} such that the cross–spectral density after the lens is given by

**u**in the sample which satisfy For a fully incoherent source, this requirement on spatial frequency has an intuitive explanation: each point in the source plane produces a separate coherent intensity pattern of the sample. This image will arrive at the detector with a transverse displacement of

*λz*

**u**with respect to the position of the sample. Since the source produces a set of mutually incoherent plane waves, the intensity pattern on the detector will be the superposition of these displaced TIE patterns,

*i.e.*a convolution of the coherent TIE result with the source intensity distribution. This convolution will degrade low spatial frequencies less than high spatial frequencies. In particular, high spatial frequencies in the measurement for which the TIE approximation is not valid will lead to inaccurate phase reconstructions, so that suppression of those spatial frequencies may be desirable. Notice that for this disk source,

**F**

_{inc}=

**0**due to symmetry. In order to produce a partially coherent field with non–zero flux, a field modulation mask may be placed in the beam prior to the sample being imaged, as illustrated in Fig. 1.

In order to test the validity of Eq. (12), we simulated the system illustrated in Fig. 1. The source is a disk of radius *R* = 250 *μ*m emitting monochromatic light with wavelength of 620 nm. The collimating lens is taken to have a focal length of *f* = 5 cm. We use a 512×512 pixel detector with 2.2 *μ*m pixels. The field modulation mask consists of a transparency with five vertical bars of OPL of 0.3 *μ*m. The field is numerically propagated (using Fresnel propagation integrals) a distance of 200 *μ*m to the in–focus plane and to an additional pair of defocused image planes at ±100 *μ*m from this plane. This is performed with and without a sample in place. Although, strictly speaking, such a system does not obey the TIE approximation for all spatial frequencies present in the image, we used a test object for which the majority of spatial frequency content does satisfy the TIE approximation. We added minor Poisson noise (SNR≈ 45 dB) to each image, such that we did not have to worry about significant denoising of the resulting images (considerably more noise is present in the experimental results presented in the following section). In order to solve both the TIE and PC–TIE, an FFT–based Poisson solver was used for direct inversion of the linear equations [26]. Since the TIE is insensitive to an additive contant offset in phase, we subtract the mean value from each OPL and phase presented here in order to compare results.

Results are illustrated in Fig. 2. In the top row, Figs. 2(a)–2(b) show the sample’s OPL and transmittance, while Fig. 2(c) shows the OPL of the mask used to modulate the illumination. The difference between defocused images Δ*S*_{tot} = *S*_{tot}(**x**, *z*) − *S*_{tot}(**x**, −*z*) is illustrated in Fig. 2 (d). The bottom row illustrates reconstructed optical thickness using a variety of techniques. In Fig. 2(e), Eq. (5) is used to construct OPL from the image, which clearly combines information about the illumination and sample. Figure 2(f) is the reconstructed scalar phase of the illumination. Figures 2(g) and 2(h) illustrate two ways of compensating for this background. In Fig. 2(g), Eq. (14) was used, which clearly does not compensate well for the illumination, leaving artifacts in regions of high transmittance variation of the sample such as the feet. In Fig. 2(h), the the PC–TIE is employed, using the previously characterized background flux in Eq. (12). Because the PC–TIE includes a term compensating for the interaction of the illumination’s flux with the sample’s transmittance, this method more effectively removes artifacts from regions in which both of these quantities are non-negligible, such as the feet. The remaining discrepancies between Figs. 2(h) and 2(a) are likely due to noise as well as spatial frequency components for which the TIE approximation does not hold. Notice the presence of low–frequency noise in all reconstructed images, since we only used extremely weak Tikhonov regularization in the inversion process to reduce the noise.

## 4. Experimental results for a pure–phase object

The PC–TIE was also tested experimentally, using the same arrangement in Fig. 1. An LED of central wavelength 620 nm is placed before a circular pinhole of 500 *μ*m. We use the previously described Köhler illumination system to create a quasi–uniform spectral density distribution after the first lens, and then illuminate either one or two objects. The first object, an MIT logo, is used as the field modulation mask, while the second, the MIT beaver, is the sample being imaged. Both samples were etched onto a glass microscope slide (*n* ≈ 1.5) to a depth of approximately 100 nm. Since the samples are of known refractive index, we reconstruct thickness rather than OPL. After the sample, the light passes through a 4f system to produce an in-focus image at unity magnification on a detector. To provide a baseline measurement for comparison, measurements were first taken with only the phase sample in place and quasi–uniform illumination [*S*(**x**) ≈ constant, **F**(**x**) ≈ **0**]. In order to test the PC–TIE, Eq. (12), measurements were taken with only the field modulating mask in place, and with both sample and mask in place. For each case, an underfocused, overfocused and in–focus images were taken with axial displacement in steps of 250 *μ*m.

In–focus intensity images for the sample alone, the illumination modulation mask alone, and both sample and mask in place are shown in Figs. 3(a)–3(c), respectively. For display, all intensities are normalized so that the maximum value of the in–focus intensity in Fig. 3(a) is unity. All reconstructions are invariant to this normalization. The difference of defocused intensities with both the mask and sample in place is shown in Fig. 3(d). This exhibits significantly more noise than the simulation, which results in low spatial frequency artifacts in the reconstruction. In order to alleviate the worst of these artifacts, we adopt Tikhonov regularization to suppress the low–frequency components in the reconstruction, although alternative techniques may prove more effective in the reduction of low–frequency noise [27, 28]. Results of the baseline reconstruction of thickness using Eq. (13) with quasi–uniform illumination are shown in Fig. 3(e). Figure 3(f) illustrates scalar phase reconstruction with only the modulation mask in place. The result of applying Eq. (5) to the measurement to retrieve thickness, without accounting for the illumination, is shown in Fig. 3(g), which clearly mixes illumination flux with sample OPL. The thickness reconstructed from Eq. (14), assuming a pure–phase sample, is shown in Fig. 3(h), while the result of applying Eq. (12) is shown in Fig. 3(i). The difference between Figs. 3(h) and 3(i) is minimal, due to the fact that the in–focus intensity of the sample exhibits minimal variation, *i.e.* is it nearly pure–phase, so that the term ∇** _{x}** · [

*T*(

**x**)

**F**

_{inc}(

**x**,

*z*)] that differentiates Eqs. (12) and (14) is negligible.

## 5. Concluding remarks

The results demonstrated here clarify the meaning of the TIE as a “phase” measurement technique when partially coherent illumination is used. Comparing images in closely spaced planes along the optical axis characterizes the divergence of a flux vector describing the transverse flow of spectral density between those planes. This flux vector’s curl–free component (under a Helmholtz decomposition) may be described by a scalar–valued phase. Understanding this distinction allows the TIE to be formulated as an OPL retrieval problem for samples under partially coherent illumination. The sample’s OPL causes a deviation of the flux vectors which can be measured by an axial stack of spectral density measurements. By compensating for the flux of the illumination, as expressed in Eq. (12), this data can be used to retrieve the sample’s OPL.

In cases where the illumination does not satisfy the TIE approximation described in Appendix A, a more general formula of the form Eq. (21) is required. For certain types of illumination, *e.g.* delta–correlated sources, this equation may be significantly simplified, such that a simple linear relationship still holds between *∂S*_{tot}/*∂z* and the sample’s OPL [19, 29]. This may create opportunities for properly chosen source distributions to enhance the performance of TIE systems, which will be considered in future work.

## A. Derivation of the TIE from finite propagation distances

In order to explore the limits in which the PC–TIE is valid, we consider in this appendix the propagation integral expressing the spectral density of a quasimonochromatic field of any state of coherence after passing through a thin sample. The results are similar, although simplified, to those given by Gureyev for more general sources [19]. The spectral density *S* produced at a propagation distance *z* due to a field with cross–spectral density *W* described over an initial plane at *z* = 0 is described by a pair of Fresnel propagation interals acting on the initial *W*,

*z*= 0, such that the cross–spectral density immediately after the sample is given by

*W*

_{inc}represents the incident field’s cross–spectral density at the sample plane, and $\sqrt{T}\text{exp}\left(\text{i}\psi \right)$ is the sample’s transmission function. By inserting this form of the spectral density into Eq. (19) and making the change of variables

**x′**→

*λz*

**u′**to simplify the results, the spectral density after propagation can be written as

The derivation of the PC–TIE proceeds from here under the assumption that all functions of *λz***u′** can be well–approximated by first–order Taylor expansions in *λz***u′**, in which case

**F**

_{inc}, is defined in terms of

*W*

_{inc}by Eq. (10). Factors of

**u′**in the integrand may be removed by expressing them as gradient operators on the exponential, since

By taking a Fourier transform of Eq. (21) from **x** to spatial frequency **u**, and performing the integration over **u′** in closed form, we can examine the validity of the TIE approximation in terms of the spatial frequency content of the image formed at the detector plane,

*Ŝ*is the 2D Fourier transform of

*S*measured by the defocused detector from transverse position

**x**to spatial frequency

**u**. Equation (24) results again from a first–order Taylor expansion with respect to

*λz*

**u**. In theory, one could always pick displacement

*z*small enough so that this expansion holds for all spatial frequencies present in the system. However, notice that the OPL of the sample is present through the Taylor expansion of

*ψ*in only one term on the right–hand side of Eq. (22), while the expansion of

*W*

_{inc}is present through

**F**

_{inc}in a separate term. If the limit on

*z*for the expansion of

*W*

_{inc}is significantly more restrictive than that of

*ψ*, the term containing

**F**

_{inc}will likely be significantly larger than that containing

*ψ*, making the OPL of the sample difficult to detect, especially in the presence of noise. If a larger propagation distance is used, the TIE reconstruction will no longer be accurate for higher spatial frequencies since

*W*is not accurately represented by this expansion. To accurately reconstruct those spatial frequencies in the sample’s OPL, a more accurate expression for propagated intensity based on Eq. (25) must be employed.

## Acknowledgments

This work was supported by the Department of Homeland Security, Science and Technology Directorate through contract HSHQDC-11-C-00083, the National Research Foundation (NRF) of Singapore through the Singapore-MIT Alliance for Research and Technology (SMART) Center for Environmental Sensing and Modeling (CENSAM) and the Chevron Energy Technology Company. The authors would like to thank Chih-Hung Hsieh and the Microsystems Technology Laboratories at the Massachusetts Institute of Technology for fabricating the test objects used in experiments and Laura Waller for useful discussions.

## References and links

**1. **J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. **21**, 2758–2769 (1982) [CrossRef] [PubMed] .

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

**3. **E. Madelung, “Quantentheorie in hydrodynamische Form,”Z. für Phys. **40**, 322–326 (1926) [CrossRef] .

**4. **D. Bohm, “A suggested interpretation of the quantum theory in terms of “hidden” variables. I,” Phys. Rev. **85**, 166–179 (1952) [CrossRef] .

**5. **D. Paganin and K. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. **80**, 2586–2589 (1998) [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**, 1932–1941 (1995) [CrossRef] .

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

**8. **T. C. Petersen, V. Keast, K. Johnson, and S. Duvall, “TEM based phase retrieval of p-n junction wafers using the transport of intensity equation,” Phil. Mag. **87**, 3565–3578 (2007) [CrossRef] .

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

**10. **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**, 5346–5362 (2013) [CrossRef] [PubMed] .

**11. **S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. W. Stevenson, “Phase-contrast imaging using polychromatic hard X-rays,” Nature **384**, 335–338 (1996) [CrossRef] .

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

**13. **T. C. Petersen, V. J. Keast, and D. M. Paganin, “Quantitative TEM-based phase retrieval of MgO nano-cubes using the transport of intensity equation,” Ultramicroscopy **108**, 805–815 (2008) [CrossRef] [PubMed] .

**14. **B. E. Allman, P. J. McMahon, K. A. Nugent, D. Paganin, D. L. Jacobson, M. Arif, and S. A. Werner, “Phase radiography with neutrons,” Nature **408**, 158–159 (2000) [CrossRef] [PubMed] .

**15. **L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge University, 1995). Ch. 4 [CrossRef] .

**16. **T. E. Gureyev, D. M. Paganin, A. W. Stevenson, S. C. Mayo, and S. W. Wilkins, “Generalized eikonal of partially coherent beams and its use in quantitative imaging,” Phys. Rev. Lett. **93**, 068103 (2004) [CrossRef] [PubMed] .

**17. **A. M. Zysk, P. S. Carney, and J. C. Schotland, “Eikonal method for calculation of coherence functions,” Phys. Rev. Lett. **95**, 043904 (2005) [CrossRef] [PubMed] .

**18. **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**, 2239–2241 (2010) [CrossRef] [PubMed] .

**19. **T. E. Gureyev, Y. I. Nesterets, D. M. Paganin, A. Pogany, and S. W. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region. 2. Partially coherent illumination,” Opt. Commun. **259**, 569–580 (2006) [CrossRef] .

**20. **K. A. Nugent, “Wave field determination using three–dimensional intensity information,” Phys. Rev. Lett. **68**, 2261–2264 (1992) [CrossRef] [PubMed] .

**21. **M. G. Raymer, M. Beck, and D. F. McAlister, “Complex wave-field reconstruction using phase-space tomography,” Phys. Rev. Lett. **72**, 1137–1140 (1994) [CrossRef] [PubMed] .

**22. **L. Tian, J. Lee, S. B. Oh, and G. Barbastathis, “Experimental compressive phase space tomography,” Opt. Express **20**, 8296–8308 (2012) [CrossRef] [PubMed] .

**23. **L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge University, 1995), Sec. 5.7.1 [CrossRef] .

**24. **N. Streibl, “Phase imaging by the transport equation of intensity,” Opt. Commun. **49**, 6–10 (1984) [CrossRef] .

**25. **L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge University, 1995), pp. 188–193.

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

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

**28. **L. Tian, J. C. Petruccelli, and G. Barbastathis, “Nonlinear diffusion regularization for transport of intensity phase imaging,” Opt. Lett. **37**, 4131–4133 (2012) [CrossRef] [PubMed] .

**29. **J. Petruccelli, L. Tian, and G. Barbastathis, “Source diversity for transport of intensity phase imaging,” to appear in to *Computational Optical Sensing and Imaging* (Optical Society of America, 2013), June 2013, paper CTu2C.3.