Boundary-artifact-free phase retrieval with the transport of intensity equation: fast solution with use of discrete cosine transform

Chao Zuo; Qian Chen; Anand Asundi

doi:10.1364/OE.22.009220

1. Introduction

Phase is an important component of an optical wavefield providing information of the refractive index, optical thickness, and the topology of the specimen. Phase retrieval is a central problem in many areas of physics and optics since the phase of a wavefield is not accessible directly. In contrast, it is relatively easy to measure the intensity of the wavefield; therefore it is tempting to seek methods for phase reconstruction from intensity measurements. The commonly known transport of intensity equation (TIE) proposed by Teague [1] is a non-interferometric method for the optical phase retrieval. This equation provides a relationship between intensity and phase of a light wave in the near Fresnel regime. With two or more intensity images taken at distinct planes orthogonal to the optical axis, solving the TIE directly gives the phase distribution, in a non-iterative manner. In the past few decades, TIE has found a variety of applications in adaptive optics [2], X-ray diffraction [3], electron- microscopy [4] and optical quantitative phase imaging [5, 6].

TIE is a second order elliptic partial differential equation for the phase function, and solving this equation does not appear to be difficult. Mathematically, the TIE phase retrieval is essentially a boundary value problem – seeking the solution of the partial differential equation subject to certain boundary conditions. The well-posed nature of this problem is well acknowledged in 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 [7]. Several method have been proposed to solve the TIE for the unknown phase: the Green’s function method [1, 8], the multi-grid method [9, 10], the Zernike polynomial expansion method [11, 12], and the fast Fourier transform (FFT) method [9, 13–15]. However, one important practical problem pervading these TIE solvers is that the associated boundary conditions are difficult to measure or to know a priori. For example, in Teague’s Green’s function solution [1], one needs to know the phase function at the region boundary as the Dirichlet boundary conditions, which is rather difficult without resorting to other means. The use of Neumann boundary conditions, studied by Roddier [2, 16] and later by Woods and Greenaway [8] has the advantage that the boundary values can be found from the intensity measurements at the pupil boundary. However, the results make one major assumption - that the intensity distribution should be uniform within the region [2, 16]. Moreover, distinguishing the boundary signal from the interior intensity derivative has been known to cause serious difficulties [11, 17]. Gureyev and Nugent [12] proposed to solve the TIE under non-classical boundary conditions – assuming the intensity distribution is non-zero inside the domain but drops to zero at the boundary. They proved that under such condition, TIE has a unique solution up to an additive constant. They then, employed orthogonal polynomial expansion (Zernike polynomial expansion for the circular region [11, 12] and fast Fourier transform (FFT) for the rectangular region [12, 13]) to solve the TIE without imposing any separate boundary conditions. However, the soft-edged intensity profile is hard to realize experimentally, and one may face the difficulty on how to define the correct domain as the intensity tends to zero as it approaches the boundary.

Another important issue among these TIE solvers is that only very few of them are sufficiently fast for reasonable sized image. The FFT-based methods [9, 13–15] are probably the most popular TIE solvers especially in the field of quantitative phase imaging [18–23] for its simplicity, efficiency, and applicability to rectangular domain. Though Gureyev and Nugent [12, 13] claimed that no boundary conditions are required for getting a unique (up to an arbitrary additive constant) solution of the TIE, the FFT-based methods have been shown to be quite sensitive to the boundary artifacts, in particular when the actual phase distribution covers the boundary of the image [15, 24, 25]. The root of the problem lies in the fact that the FFT solver implying periodic boundary conditions due to the cyclic nature of the discrete Fourier transform, which is often in contradiction with the actual situations. Throughout literature [18–23], a common way to avoid such boundary error is to consider the object to be isolated and placed in the center of the image [24]. However, this configuration does not reflect the general experimental conditions, and is not practical when the measured object is larger than the camera field of view. An alternative FFT-based method proposed by Volkov et al. [15] extends the measured intensity by a mirror padding scheme, which provides homogenous Neumann boundary conditions rather than the periodic ones. However, it still does not fundamentally solve the boundary artifacts problem since the boundary conditions are formed purely by a mathematical trick which is generally not physically grounded.

In this paper, we will solve the TIE with experimentally measurable inhomogeneous Neumann boundary conditions, through an efficient numerical algorithm using the fast discrete cosine transform (DCT). This approach not only addresses the above mentioned difficulties such as obtaining boundary conditions under non-uniform intensity distribution, distinguishing the boundary signal from the interior intensity derivative, and defining the correct domain for the solution, but also retains the major advantages of the FFT-based methods – simple and fast for the rectangular domain. Rather than conventional differential equation methods where the TIE is discretized directly, we derive the analytic integral solution via the Green’s function method, and then the DCT is applied to calculate the integral solution numerically. All the efforts are made to simplify the numerical calculation process as well as the enforcement of the experimentally obtained Neumann boundary conditions. Its good performance is demonstrated by several numerical examples.

2. Transport of intensity equation, boundary conditions, and uniqueness of the solution

Assume that a paraxial beam is propagating along the $z$ axis, and let the complex amplitude of the object be $\sqrt{I (r)} \exp [i k ϕ (r)]$ . The derivative of intensity along the beam propagation direction, $z$ contains phase information that can be retrieved via TIE [1]:

- k \frac{\partial I (r)}{\partial z} = \nabla \cdot [I (r) \nabla ϕ (r)],

where

k

is the wave number

2 π / λ

,

r

is the position vector representing the spatial coordinates

(x, y)

.

I (r)

is the intensity, located without loss of generality at the plane

z = 0

.

\nabla

is the gradient operator over

r

which is normal to the beam propagation direction. Expanding the right hand side (RHS) of Eq. (1), we obtain

- k \frac{\partial I}{\partial z} = \nabla I \cdot \nabla ϕ + I \nabla^{2} ϕ .

Note here the

r

dependence has been omitted for notation simplicity. It is seen that the TIE links the longitudinal intensity derivative with the slope and curvature of the wavefront that produces the changes in intensity as the wavefront propagates.

2.1 Boundary conditions and uniqueness of the solution

The TIE is an elliptic partial differential equation for the phase function, $ϕ$ . We assume the region governed by the TIE to be a general open and bounded domain $Ω \subset ℝ^{2}$ with a piecewise smooth boundary $\partial Ω$ . The intensity distribution $I$ is a continuous, non-negative function defined on the enclosure $\bar{Ω}$ (consisting of the domain $Ω$ plus its boundary $\partial Ω$ ) and is smooth and strictly positive in $Ω$ . The longitudinal intensity derivative $\partial I / \partial z$ is assumed to be a smooth function in $Ω$ . The phase $ϕ$ is expected to be single-valued and smooth in $\bar{Ω}$ . Before deriving a solution to the TIE, the issues concerning its solvability (well-posedness) and the uniqueness of the solution must be addressed. The classic theory of elliptic partial differential equation states that the uniqueness of the TIE solution subject to certain boundary conditions [26]. We shall consider two classes of possible boundary conditions that could be applied to the solution of TIE: (1) Dirichlet boundary conditions. The values of phase function $ϕ$ are specified on the domain boundary

{ϕ |}_{\partial Ω} = g .

(2) Neumann boundary conditions. The product of

I

and the normal derivative of

ϕ

, are specified on the domain boundary

{I \frac{\partial ϕ}{\partial n} |}_{\partial Ω} = g .

Here

g

is a smooth function on the boundary

\partial Ω

,

\partial ϕ / \partial n

is the outward normal derivative. For the case of Dirichlet boundary conditions, the solution to the TIE always exists and is unique [26]. The case of the Neumann boundary conditions demands special attention because a solution may or may not exist (depending on whether the compatibility condition (Eq. (5)) is satisfied [7, 27]). Furthermore, it is clear that

ϕ

= constant satisfies the homogeneous problem with

- k \partial I / \partial z = 0

and

g = 0

, therefore if

ϕ

is a solution, then

ϕ + C

is also a solution for any constant

C

. In other words, the solution of the TIE subject to Neumann boundary conditions (Eq. (4)) should be unique up to an arbitrary additive constant, assuming that the solution exists in the first place. This constant is not essential for the phase retrieval problem. The proof of the uniqueness theorem is straightforward: we can first assume that two qualified solutions have been exhibited and then prove they can only differ from one additive constant using the maximum principle for elliptic equations [28].

2.2 Compatibility condition and energy conservation law

For the Neumann boundary problem, there is an important compatibility condition that must be satisfied in order that a solution exists. Starting from the TIE (Eq. (1)), we integrate both sides of the equation over $Ω$ and apply the divergence theorem to get the compatibility condition

\oint_{\partial Ω} I (r) \frac{\partial ϕ (r)}{\partial n} d s = \iint_{Ω} - k \frac{\partial I (r)}{\partial z} d r,

where

s

is a parameterization of

\partial Ω

. The RHS of Eq. (5) is known a priori, since the function

\partial I / \partial z

can be experimentally obtained. The left hand side (LHS) of Eq. (5) is determined solely by the Neumann boundary conditions (Eq. (4)). Thus, mathematically all these possible Neumann boundary conditions must respect Eq. (5) in order for there to be a nontrivial solution to the TIE. Generally, this compatibility condition is only considered to be a necessary condition for the solvability of the TIE. But with the smoothness assumptions imposed on

I

and

g

, the sufficiency of the condition can be justified. In other words, once Eq. (5) is satisfied, a solution to the Neumann boundary value problem (Eqs. (1) and (4)) must exist (and use the uniqueness theorem, this solution is unique up to an overall additive constant).

It is more interesting to examine the physical picture described by Eq. (5): actually, it can be thought as an expression of local energy conservation law - the loss of energy (intensity) inside the region arising from energy flow across the region boundary. If we extend $Ω$ to the unbounded space, the contour integral vanishes and Eq. (5) becomes

\iint_{ℝ^{2}} \frac{\partial I (r)}{\partial z} d r = 0.

It represents the law of energy conservation for unlimited free space. In essence, the energy conservation expressed by Eqs. (5) and (6) is a universal law of physics, and their validity relies only on the paraxial approximation.

2.3 Solving the TIE without any boundary values

Since solving the TIE requires boundary conditions either the Dirichlet boundary value or the normal derivative of the phase (Neumann boundary conditions) to ensure the uniqueness of the solution, and both kinds of boundary data are not easy to acquire experimentally, researchers are looking for ways to solve the TIE without taking any boundary measurements. Coincidently, all the efforts aim to find some way to nullify the LHS of Eq. (5), making boundary conditions unnecessary

\iint_{Ω} \frac{\partial I (r)}{\partial z} d r = 0.

Generally, Eq. (7) cannot be really valid on a bounded region unless the contour integral (LHS of Eq. (5)) vanishes. Physically, it is equivalent to the special case that there is no overall energy transfer across the region boundary.

One simple and common way to satisfy this is to let the test object be centrally isolated in the field of view. Then the zero phase changes at the boundary (corresponding to homogeneous Neumann boundary conditions ${I \partial ϕ / \partial n |}_{\partial Ω} = 0$ ) can be assumed thus the LHS of Eq. (5) vanishes. In fact, in this case one can define not only the Neumann boundary conditions, but any uniform Dirichlet boundary conditions (constant phase at the boundary ${ϕ |}_{\partial Ω} = C$ ) or periodic ones (the phase at the boundary repeats cyclically) as well. Gureyev and Nugent [12] suggested another way to eliminate the need of boundary conditions (LHS of Eq. (5)) by considering the case that the intensity $I > 0$ inside the domain $Ω$ but strictly vanishes at the boundary (actually it corresponds to a special case of the homogeneous Neumann boundary conditions ${I \partial ϕ / \partial n |}_{\partial Ω} = 0$ though not mentioned in [12]). Alternatively, without any additional requirement about the test object and experimental conditions, Volkov et al. [15] proposed a pure mathematical trick to nullify the energy flow across the boundary through appropriate symmetrization of input images. However, it assumes there is no energy dissipation through the image boundary for any class of images and objects, which is generally not physically grounded. To summarize former discussion, without boundary value measurements does not mean that we can solve the TIE without imposing any boundary conditions, or more exactly, we have to confine our test object or experimental configuration to some implicit boundary conditions. Obviously, the assumption of the test object to be isolated is apparently limiting and does not reflect general cases. Furthermore, the assumption is impractical when the object is larger than the image field of view. On the other hand, intensity measurement satisfying Gureyev and Nugent’s conditions tends to be difficult as well. One needs to use a special-designed soft-edged aperture or use specific illumination with intensity vanishing at the boundary, $\partial Ω$ . Furthermore, defining the boundary of the domain tends to be very difficult in practice since very-near zero and at zero may seem almost the same but can be entirely different vis-a-vis the well-posedness of the TIE.

3. Solution to transport of intensity under non-uniform intensity

As suggested in last section, because of its physical origin, solving the TIE essentially boils down to boundary value problems – to find a solution in $Ω$ with conditions imposed on its boundary, $\partial Ω$ . However, the major obstacle for solving the TIE is that the associated boundary conditions are difficult to measure or know a priori. Therefore, in this section, a simple approach to get the boundary conditions and construct an analytic solution to the TIE using Green’s function method is suggested.

3.1 Boundary values generation with a hard-edged aperture

Let us consider a particular intensity distribution in the object plane ( $z = 0$ ): the intensity $I$ is smooth (but can be non-uniform) in $Ω$ but suddenly vanishes with a step-like discontinuity at the boundary, $\partial Ω$ . Physically, it is equivalent to capturing intensity images through a hard-edged aperture at the object plane:

I = A_{Ω} I_{0} = {\begin{cases} I_{0} r \in \bar{Ω} \\ 0 o t h e r s \end{cases},

where

I_{0}

is the intensity when there is no aperture;

A_{Ω}

is the aperture function (

A_{Ω}

= 1 when

r \in \bar{Ω}

,

A_{Ω}

= 0 when

r \notin \bar{Ω}

). Since the intensity

I

is smooth inside

Ω

, it is reasonable to assume that its gradient on the boundary,

\partial Ω

is dominated by the component in the boundary normal direction, which has the impulse value that equals the negative of the intensity boundary value and points outward from the domain:

\nabla I = {\begin{cases} - I δ_{\partial Ω} n r \in \partial Ω \\ \nabla I_{0}, r \in Ω \\ 0 others \end{cases},

where

δ_{\partial Ω}

is the Dirac delta function around the aperture edge

\partial Ω

and

n

is an outward-pointing unit vector, normal to

\partial Ω

. Strictly speaking, the gradient

\nabla I

does not exist in calculus since

I

is discontinuous at

\partial Ω

. So we call

\nabla I

here as a generalized gradient. Here we adopted the notation

δ_{\partial Ω}

that first used by Roddier [2] as a plausible representation of the physical quantity concentrated on the boundary curve. The

δ_{\partial Ω}

can be loosely thought of as the gradient magnitude of the aperture function

δ_{\partial Ω} = | \nabla A_{Ω} |

or interpreted as the limit of a sequence of smooth functions. Mathematically, the function

δ_{\partial Ω}

is an extension of the Dirac delta function

δ_{r}

at single point

r

, which can be rigorously defined as a distribution such that for any test function

f (r)

that is smooth and compactly supported

\iint_{ℝ^{2}} f (r) δ_{\partial Ω} d r = \oint_{\partial Ω} f (r) d s .

Substituting Eqs. (8) and (9) into the TIE (Eq. (2)), we obtain

- k \frac{\partial I}{\partial z} = A_{Ω} (I_{0} \nabla^{2} ϕ + \nabla I_{0} \cdot \nabla ϕ) - I \frac{\partial ϕ}{\partial n} δ_{\partial Ω} .

The first term on RHS is the intensity variation inside the domain due to the phase slope and curvature as if the aperture is not present. The second term is a delta-function-like a sharply peaked signal at the aperture boundary, which provides the exact Neumann boundary conditions for the TIE (Eq. (4)). Since the axial intensity derivative in LHS is experimentally measurable, and the two RHS terms do not overlap in space, there is enough information to solve the TIE with no other measurement specially for the boundary values (Teague [1] suggested to measure the phase values on the domain boundary for the Dirichlet boundary conditions (Eq. (3)) using Hartmann sensors). However, from a purely mathematical viewpoint, the new form of the TIE defined in Eq. (11) is not rigorous because the function

δ_{\partial Ω}

belongs to the class of generalized functions, which is meaningful only as a part of an integral expression. Thus, this new form serves merely as an apparatus for describing the physical quantities (longitudinal intensity derivative) in a convenient and adequate manner. When writing the unusual form of the TIE as in Eq. (11), we always keep in mind that

δ_{\partial Ω}

will eventually become a part of an integral, and produce valid results that can be used.

Let us put aside these complications for a while and assumes we already have the separated longitudinal intensity derivative inside $Ω$ and the amplitude of the $δ_{\partial Ω}$ along the boundary line through some peculiar way. When applying these measured quantities to solve the TIE (us the interior signal as the source term of the TIE (LHS of Eq. (1)) and take the amplitude of delta-function term in Eq. (11) as the function $g$ in Eq. (4)), we must make sure the Neumann boundary value problem is well-posed and yields a unique solution to the problem. Integrating both sides of Eq. (11) over all space and using Eq. (10), we obtain

\iint_{ℝ^{2}} - k \frac{\partial I (r)}{\partial z} d r = \iint_{\bar{Ω}} - k \frac{\partial I (r)}{\partial z} d r = \iint_{Ω} - k \frac{\partial I (r)}{\partial z} d r - \oint_{\partial Ω} I (r) \frac{\partial ϕ (r)}{\partial n} d s .

According to the energy conservation law for unbounded space (Eq. (6)), the LHS of the above equation should equal to zero

\iint_{\bar{Ω}} - k \frac{\partial I (r)}{\partial z} d r = \iint_{Ω} - k \frac{\partial I (r)}{\partial z} d r - \oint_{\partial Ω} I (r) \frac{\partial ϕ (r)}{\partial n} d s = 0,

which is the same expression as the compatibility condition, Eq. (5) in Section 2.2. This leads to the inescapable conclusion that the solution to this Neumann boundary problem always exists and is unique up to an arbitrary additive constant. Note a similar approach was suggested by Roddier [2, 16] to obtain the normal derivative of the phase function at the pupil edge. Such an approach relies on the major assumption that the intensity has to be uniform inside the region

Ω

, which is reasonable in adaptive optics. However, our derivation represents a more general case and is valid for non-uniform intensity distribution - which is quite common in many practically important applications such as quantitative phase imaging.

3.2 Analytic solution to the TIE with the Green’s function approach

In this subsection, we derive a special analytic solution to the TIE with the Neumann boundary conditions found from the intensity derivative boundary based on Green’s function approach. The reason why we do not use numerical methods directly is that these methods generally require the intensity derivative inside the aperture and the Neumann boundary value (delta-function term at the aperture edge) being treated separately. More specifically, in order to get the Neumann boundary value, one has to find some appropriate ways to extract the amplitude of the delta-function around the aperture edge. However, as we mentioned previously, the Eq. (11) is an only an idealized representation which is not suitable for direct application. The delta function $δ_{\partial Ω}$ is only defined on infinitely small support and takes infinite values, reflecting the fact that a physical quantity cannot be measured along one line. Although in reality the boundary signal assumes finite values and can be measured by a camera, it occupies a finite area, bringing about a series of non-trivial issues like how to define the boundary region, and how to separate the boundary values from the interior signals [11, 17]. In the following, we will introduce a way to bypass these difficulties.

Firstly, as suggested by Teague [1], an auxiliary function $ψ$ was introduced, which satisfies

\nabla ψ = I \nabla ϕ

This physical meaning of the vector field

I \nabla ϕ

can be interpreted as the Poynting vector [14, 29], which is characterized by the scalar potential

ψ

. With this auxiliary function

ψ

, the elliptic TIE can be converted down to a standard Poisson equation [9, 13, 23]

- k \frac{\partial I}{\partial z} = \nabla^{2} ψ,

with the associated Neumann boundary conditions unchanged

{\frac{\partial ψ}{\partial n} |}_{\partial Ω} = {I \frac{\partial ϕ}{\partial n} |}_{\partial Ω} .

Once

ψ

is solved, the phase function can be uniquely determined up to an arbitrary additive constant by integrating Eq. (14). The Green's function representation of the solution to the Neumann boundary-value problem (Eqs. (14) and (16)) is given by

ψ (r) = \iint_{Ω} (I \nabla^{2} ϕ + \nabla I \cdot \nabla ϕ) G (r, r^{'}) d r^{'} - \oint_{\partial Ω} G (r, r^{'}) I (r^{'}) \frac{\partial ϕ (r^{'})}{\partial n} d s^{'},

where

G (r, r^{'})

is the Neumann Green’s function, which is the solution of the Neumann boundary value problem:

\nabla^{2} G (r, r^{'}) = δ (r - r^{'}) - A^{- 1}, \frac{\partial G (r, r^{'})}{\partial n} |_{Ω} = 0,

where

A

is the area of the region

Ω

[30]. As we can see, the Green's function method represents the solution of the partial differential equation by one area and one curve integral of the Green's function. The two terms stands for the two contributory sources to the solution: the intensity derivative inside the aperture and the boundary signal at the aperture edge. It seems that the issue on how to distinguish the boundary value from the interior signal has not yet been solved. However, if we multiply both sides of Eq. (11) by the Green’s function

G (r, r^{'})

and then integrate it over all space, we obtain:

\begin{array}{l} \iint_{ℝ^{2}} - k \frac{\partial I (r^{'})}{\partial z} G (r, r^{'}) d r^{'} = \iint_{\bar{Ω}} - k \frac{\partial I (r^{'})}{\partial z} G (r, r^{'}) d r^{'} \\ = \iint_{Ω} [A_{Ω} (I_{0} \nabla^{2} ϕ + \nabla I_{0} \cdot \nabla ϕ) - I_{\partial Ω} \frac{\partial ϕ}{\partial n} δ_{\partial Ω}] G (r, r^{'}) d r^{'} \\ = \iint_{Ω} (I \nabla^{2} ϕ + \nabla I \cdot \nabla ϕ) G (r, r^{'}) d r^{'} - \oint_{\partial Ω} G (r, r^{'}) I (r^{'}) \frac{\partial ϕ (r^{'})}{\partial n} d s^{'} . \end{array}

Comparing Eq. (19) with Eq. (17), we can rewrite the Green's function solution into a more compact area integral:

ψ (r) = \iint_{\bar{Ω}} - k \frac{\partial I (r^{'})}{\partial z} G (r, r^{'}) d r^{'} .

Equation (20) is the central result of this section. Here we should reemphasize that the region

\bar{Ω}

in Eq. (20) should be closed so that all the boundary values in the intensity derivative signal is covered by this area integral. By doing so, the solution of the Poisson equation can be simply represented as only one area (convolution) integral over the aperture region so that the two terms (interior signal and boundary value) can be considered as one entity, instead of separately. Before going further, there are two important points need to be clarified. First, as we know for the normal integrable functions, there makes no difference between the integral over the interior

Ω

and its closure

\bar{Ω}

. However, this rule cannot be applied to our case since the integrand in Eq. (19) contains the generalized function

δ_{\partial Ω}

. As shown in Eq. (10), the generalized function

δ_{\partial Ω}

defines a distribution mapping from function to the value of its curve integral, which is used here to greatly simplify the original Green’s function solution (Eq. (17)). Second, the new form of solution (Eq. (20)) does not go against the fact that the solution of the partial differential equation must depend on the boundary conditions, since the experimentally measured boundary value information has already been encapsulated in the area integral. It can be seen that once the Green’s function

G (r, r^{'})

is determined, the solution of the Poisson equation (Eq. (15)) can be obtained from Eq. (20) with single area integral via Green's function; and we emphasize the Green's function here is solely determined by the shape of the region

Ω

, and independent on experimental data.

Once $ψ$ is obtained by the area integral (Eq. (20)), the phase function $ϕ$ can be reconstructed from its determined gradient $\nabla ϕ = I^{- 1} \nabla ψ$ uniquely up to an additive constant. Mathematically, phase integration is a well-defined problem and can be simply achieved by a line-by-line accumulation procedure. However, the phase integration is ill-posed (and not well defined) as the vector field $I^{- 1} \nabla ψ$ is rarely integrable (see Sec. 6.2 for details). A common approach is to find the scalar function whose gradients are closest to the determined vector field $I^{- 1} \nabla ψ$ in a least squares sense. It can be shown by calculus of variations [31, 32], the solution to this least squares error problem is equivalent to that of finding the solution of the following Poisson equation

\nabla \cdot (I^{- 1} \nabla ψ) = \nabla^{2} ϕ .

with homogeneous Neumann boundary conditions

{\frac{\partial ϕ}{\partial n} |}_{\partial \bar{Ω}} = 0.

The Green's function solution of the Neumann boundary problem (Eqs. (21) and (22)) is

ϕ (r) = \iint_{\bar{Ω}} \nabla \cdot (I^{- 1} \nabla ψ) G (r, r^{'}) d r^{'},

where

G (r, r^{'})

is the Neumann Green’s function, which is exactly the same one as in Eq. (20). As can be seen, Eq. (20) and (23) give analytic convolution integral solution to the TIE with the Green’s function, and the object phase distribution can be obtained by calculating the two area integrals numerically. There are two major benefits of our approach. First, it explicitly considers the case of non-uniform intensity distribution. Second, the delta-function terms at the aperture edge (second terms of RHS of Eqs. (11) and (17)) are enclosed in Eq. (20), without requiring special-purpose detection scheme to separate these boundary values from the interior signals. The well-posedness and uniqueness of the solution is guaranteed by the energy conservation law (Eq. (13)), which can be used to define the integral domain

\bar{Ω}

and check the consistency of experimentally measured intensity data.

4. Fast numerical implementation with use of discrete cosine transform

Evidently, the analytic solution to TIE (Eqs. (20) and (23)) is not very convenient for numerical implementation since the Green's function is a four dimensional function - the direct calculation is memory and computationally demanding. For an $M \times N$ pixel image, $G (r, r^{'})$ will contain $M^{2} N^{2}$ elements, which may require thousands of gigabyte memory to store. Besides, the direct numerical implementation of the integral requires expensive operations of $O (M^{2} N^{2})$ complexity. Therefore, in this section, we focus on developing an efficient numerical implementation of the solution. More specifically, we consider the case of a rectangular domain since most detectors (like the CCD camera) are of rectangular geometry, thus in order to maximally fill the camera sensor area, using a rectangular aperture is a natural choice. It is found that the Green's function under Neumann boundary conditions with a rectangular region can be represented as an infinite-series expansion in terms of Fourier cosine harmonics, leading to a fast and efficient numerical implementation of the TIE solution with use of fast DCT. It should be noted that in practice this approach can be easily extended to the domains with arbitrary shapes (e.g. the circular or annular regions) as long as the corresponding eigen-functions of Laplacian can be obtained.

4.1 Eigen-functions expansion of the Green's function

Although the analytic representation of the Green's function of Poisson equation under Neumann boundary conditions for a rectangular domain $[0, a] \times [0, b]$ can be found in many mathematical textbooks [30, 33], here we prefer to expand it into an infinite series to gain some insight to the efficient numerical implementation. Consider the Green's function in Eq. (20) associated with the Neumann condition problem, which obeys Eq. (18), we can expand it in terms of eigen-functions of Laplacian [33]

G (r, r^{'}) = \sum_{m = 0}^{\infty} \sum_{n = 0}^{\infty} \frac{Ψ_{m, n} (r) Ψ_{m, n} (r^{'})}{λ_{m, n}} .

Here

Ψ_{m, n}

are eigen-functions of Laplacian satisfying the homogeneous Neumann conditions on the boundary:

\nabla^{2} Ψ_{m, n} = λ_{m, n} Ψ_{m, n}

,

\partial Ψ_{m, n} / \partial n |_{Ω} = 0

and

λ_{m, n}

are the associated eigenvalues. According to Sturm-Liouville theory [33], the eigen-functions form a complete and orthogonal basis over the domain

\bar{Ω}

. More specifically, for a rectangular domain, the eigen-functions can be represented as the Fourier cosine harmonics

Ψ_{m, n} (x, y) = a_{m, n} \cos (\frac{m π}{a} x) \cos (\frac{n π}{b} y),

with the normalization coefficients

a_{m, n} = {\begin{cases} \frac{2}{\sqrt{a b}}, m, n \neq 0 \\ \frac{\sqrt{2}}{\sqrt{a b}}, m o r n = 0 \end{cases},

and the associated eigenvalues

λ_{m, n} = - π^{2} (\frac{m^{2}}{a^{2}} + \frac{n^{2}}{b^{2}}) .

Substituting Eqs. (24)–(27) into Eq. (20), and interchanging the order of summation and integral, we obtain

ψ (r) = \sum_{m = 0}^{\infty} \sum_{n = 0}^{\infty} λ_{m, n}^{- 1} [a_{m, n}^{2} \iint_{\bar{Ω}} - k \frac{\partial I (r^{'})}{\partial z} \cos (\frac{m π}{a} x^{'}) \cos (\frac{n π}{b} y^{'}) d r^{'}] \cos (\frac{m π}{a} x) \cos (\frac{n π}{b} y) .

Although it seems to be more complex, Eq. (28) explicitly interprets the convolution integral as a three-step process. As it is shown in Appendix A, the expression inside the square brackets on the RHS is by definition the coefficient of the Fourier cosine series expansion of the axial intensity derivative signal. These coefficients are then multiplied with the factor

λ_{m, n}^{- 1}

, and finally the function

ψ (r)

is recomposed by the Fourier cosine basis with these modified coefficients. The same approach can be applied to compute Eq. (23) to obtain the phase function. The major advantage of using the Fourier cosine expansion is the availability of the fast algorithm, which considerably increases the speed and efficiency of all the numerical computations. We will discuss it in more detail in next subsection.

4.2 Numerical implementation with fast discrete cosine transform

In practice, the intensity distributions are always recorded by a pixelated detector - the continuous function $\partial I / \partial z$ is discretized into a finite $M \times N$ rectangular grid cells of equal size. For convenience, in this subsection we still use the $r = (x, y)$ to represent the image traverse coordinate, but keep in mind they can only assume integral values, i.e., $x = 0, 1, ..., M - 1;$ $y = 0, 1, ..., N - 1$ . The pixel coordinate $(x, y)$ corresponds the grid cell centered at the location $(Δ x (0.5 + x), Δ y (0.5 + y))$ in the continuous real space within the rectangular domain $[0, a] \times [0, b]$ , where $Δ x = a / M$ , $Δ y = b / N$ is the cell interval (pixel size) in the $x$ and $y$ direction, respectively. Besides the intensity derivative $\partial I / \partial z$ , the cosine harmonics basis functions need to be discretized as well. As shown in Appendix B, the set of functions $\cos [m π (x + 0.5) / M] \cos [n π (y + 0.5) / N]$ form an orthogonal basis over the $M \times N$ discrete space while satisfies the homogeneous Neumann condition on the boundary. Correspondingly, the discretized form of Eq. (28) can be represented as

S_{m, n} = A_{m} A_{n} \sum_{x^{'} = 0}^{M - 1} \sum_{y^{'} = 0}^{N - 1} - k \frac{\partial I (r^{'})}{\partial z} \cos [\frac{m π}{M} (x^{'} + 0.5)] \cos [\frac{n π}{N} (y^{'} + 0.5)];

ψ (r) = \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} A_{m} A_{n} λ_{m, n}^{- 1} S_{m, n} \cos [\frac{m π}{M} (x + 0.5)] \cos [\frac{n π}{N} (y + 0.5)] .

The normalization coefficients

A_{m} A_{n}

for the discrete cosine basis functions are defined in Eq. (40). It is shown that the two equations can be considered as the forward and backward (inverse) DCT (see the definition of DCT in Appendix A). Thus, with use of DCT, the numerical implementation of Eq. (20) is quite straightforward: perform DCT of the measured axial intensity derivative, reweight its DCT coefficients with the factor

λ_{m, n}^{- 1}

, and then perform inverse DCT,

ψ

can be obtained. The phase distribution

ϕ

can be obtained by Eq. (23) in a similar way. Since the DCT and inverse DCT can be calculated efficiently using some fast algorithm or using the FFT in

O (M N \log (M N))

time (one simple FFT-based fast DCT algorithm is given in Appendix B), the complexity of integral calculation is also of the order

O (M N \log (M N))

, which is a significant improvement over the

O (M^{2} N^{2})

complexity of the direct computation. Meanwhile, the memory requirement is substantially reduced since there is no need to store the four dimensional Green's function. Note that although the DCT has been already used as a fast solver for the discrete Poisson equation with Neumann boundary conditions, our method is quite different from the common one in the textbook [34] which uses the finite difference numerical method (using a five-point stencil approximating the Laplacian) to discretize the Poisson equation into a

M N \times M N

linear system.

So far, the only remaining problem is to calculate the source term $\nabla \cdot (I^{- 1} \nabla ψ)$ in Eq. (23). As $ψ$ can be expanded into DCT basis vectors, its derivative can be efficiently calculated through DCT-based differential operators (all the relevant operators can be found in Appendix C). Once the gradient field $\nabla ψ = (\partial_{x - D C T} ψ /, \partial_{y - D C T} ψ)$ is obtained from Eqs. (53) and (55), the vector field $I^{- 1} \nabla ψ = (I^{- 1} \partial_{x} ψ, I^{- 1} \partial_{y} ψ)$ can be easily calculated by point-wise multiplication. Finally, the divergence $\nabla \cdot (I^{- 1} \nabla ψ)$ can be obtained through

\nabla \cdot (I^{- 1} \nabla ψ) = \nabla_{D C T} \cdot (I^{- 1} \nabla_{D C T} ψ) = \partial_{x - D C T} (I^{- 1} \partial_{x - D C T} ψ) + \partial_{y - D C T} (I^{- 1} \partial_{y - D C T} ψ),

with two more DCT-based derivative operations. As we mentioned in Appendix C, the basis functions of the obtained derivative are no longer cosine functions, which means we cannot reconstruct the modified coefficients via inverse DCT. One can either directly calculate it with the Eqs. (53) and (55), or use some fast algorithms (as shown in Appendix D, the DCT-based derivatives can be efficiently calculated through FFT as well).

Combining the above analysis, the complete numerical implementation of the TIE solution can be summarized as the following steps:

Calculate the intensity derivative term $- k \partial I / \partial z$ from the measured intensity images.
Transform $- k \partial I / \partial z$ into cosine frequency domain by DCT (Eq. (43)).
Reweight the DCT coefficients by $λ_{m, n}^{- 1}$ (Eq. (27)).
Transform the modified DCT coefficient back to spatial domain to get $ψ$ (Eq. (44)).
Calculate $\nabla ψ$ based on DCT (Eqs. (53), (55), and (58)).
Calculate $I^{- 1} \nabla ψ$ .
Calculate $\nabla \cdot (I^{- 1} \nabla ψ)$ (Eqs. (53), (55), and (59)).
Transform $\nabla \cdot (I^{- 1} \nabla ψ)$ into cosine frequency domain by DCT.
Reweight the DCT coefficients by $λ_{m, n}^{- 1}$ .
The phase can be determined up to an additive constant by transforming the modified DCT coefficient back to spatial domain. Steps (8)-(10) are analogous to steps (2)-(4).

It should be noted that all the source intensity images should be strictly defined on the rectangular-shaped area $\bar{Ω}$ , which includes both the aperture boundary and the region inside it. However, it is evident that the above steps have redundant operations: the auxiliary function $ψ$ is the intermediate step in the solution and therefore should not be express explicitly. Taking out all the redundancies from the above steps and adopting the simple notation we used in Appendix C, the final algorithm is given in Algorithm 1.

Algorithm 1: Solution to TIE using the DCT-based Green's function method under non-uniform intensity distribution ( $I \neq constant$ ). (This gives identical results to Algorithm 2 when $I = constant$ ).

View Table | View all tables in this article

Further, when the intensity is uniform in the defined region, Algorithm 1 can be further simplified as:

Algorithm 2: Solution to TIE using the DCT-based Green's function method under non-uniform intensity distribution ( $I = constant$ ).

View Table | View all tables in this article

In this case, the TIE can be directly simplified as a Poisson equation without the need of the auxiliary function $ψ$ , which corresponds to the case discussed in [8,16]. However, our method involves only one forward and one inverse DCT, which substantially reduces the computation and memory demand compared with the classical Green's function method [8]. As shown in Appendices B and D, the DCT as well as all DCT-based derivative operators involved can be easily calculated through FFT with even extension scheme. These methods are quite helpful to simplify the numerical implementation of the proposed algorithms since FFT is available in almost all computer packages, but fast algorithms for DCT and some other transforms like Eq. (54) are not. The two alternatives to Algorithm 1 and 2 based on FFT and even extension are given in Algorithm 3 and 4, respectively.

Algorithm 3: Alternative to Algorithm 1 based on FFT and even extension (this gives the exactly the same solution as Algorithm 1).

View Table | View all tables in this article

Algorithm 4: Alternative to Algorithm 2 based on FFT and even extension] (this gives the exactly the same solution as Algorithm 2).

View Table | View all tables in this article

The meaning of the bar (‘-’) and angle brackets notations (‘< >’) used here can be found in Appendix D. At the first glance, the Algorithms 3 and 4 seem to be very similar with the previous work using FFT and even symmetrization [15]. However, one should note the key differences: our algorithm only applies to the experimental data captured with the presence of an aperture, which means the axial derivative of intensity generally contains a delta-function term at the boundary. It was derived from the Green's function method with the inhomogeneous Neumann boundary value measured at the aperture edge without imposing any implicit or prior assumptions on test objects. We found the solution to the Poisson equation can be expressed by cosine series expansion and then we use DCT as a fast numerical method for the Green's function solution. The FFT and even symmetrization schemes merely serve as alternative ways to perform DCT and its related computations. While previous work [15] used even symmetrization to nullify the contour integral (corresponding to the second terms of RHS of Eq. (17)) by making a silent assumption for the test object to be even symmetry, which is usually unrealistic especially when the object extends out of the image boundary, leading to non-optimal or even erroneous results.

Before proceeding further, several points on practical implementation of our algorithm need to be clarified: First, $λ_{m, n}^{- 1}$ is not defined at the zero point of the DCT space, therefore, we multiply by zero instead. Second, since the intensity $I$ appears in the denominator, one should avoid dividing too small values to guarantee the stability of solution by setting a certain threshold value (e.g. 1% of the maximum value of $I$ ). Finally, one should correctly define the domain for each captured intensity image before carrying out the algorithms. It is trivial when the small defocus distance requirement is fulfilled (Eq. (11) is valid) - one may simply choose the pixels within aperture region $Ω$ as well as those around the aperture boundary $\partial Ω$ . However, when the defocus distance is not sufficiently small (increasing the defocus distance is a common way to improve the SNR of the intensity derivative estimates [35]), some boundary signal will either go into the aperture zone or escape outside from the aperture edge due to the diffraction effect (Eq. (11) will no longer well hold). In this case, defining the proper domain seems problematic, and we will examine this point in Section 5.4.

5. Simulations and comparisons

Simulations were performed to test the proposed method. We simulated three challenging situations for TIE phase retrieval. In these simulations, the complex field was defined on a 256 × 256 grid with 2μm × 2μm pixel size. The wavelength is 632.8nm and the object is assumed to be larger than the camera field of view (200 × 200 pixels) which represents a more general situation than the common assumption in the literatures that the object is isolated and placed in the center of the camera field of view. Angular spectrum numerical propagation was used to generate two defocused images $I_{+}$ and $I_{-}$ with the distance of $Δ z = \pm 10$ μm. The intensity derivative $\partial I / \partial z$ is estimated by central finite difference

\frac{\partial I (r)}{\partial z} \approx \frac{I_{+} (r) - I_{-} (r)}{2 Δ z} .

The metric relative root-mean-square error (RMSE) is 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 RMSE evaluation.

5.1 Pure phase object with uniform intensity distribution

First we focus on the case that the test object is a pure phase function with unitary intensity. As shown in Fig. 1(a), the phase function is predefined as $ϕ (x, y) = 10 x^{2} - 10 y^{2} - 0.7 x + 2 y + 0. 82$ where the lateral dimensions $x$ and $y$ are limited within a range from −0.2 mm to 0.2 mm in the camera field of view. Since the object is a pure phase one, $I$ is constant when no aperture exists, and can be pulled out of the divergence operator [21, 36, 37] from the TIE (Eq. (1)), resulting in a standard Poisson equation

- k \frac{\partial I}{\partial z} = I \nabla^{2} ϕ .

This Poisson equation states that the intensity derivative is proportional to the Laplacian of the phase. Since the simulated phase function has zero Laplacian, the difference between two defocus images

I_{+}

and

I_{-}

provided a null signal as expected (Fig. 1(b)), therefore reconstructing the phase using these intensity measurements can never be achieved in this situation [38, 39]. This once again reveals the dependence of the TIE solution on the boundary conditions. By introducing an aperture (aperture size 180 × 180 pixels) located at the center of the object plane, intensity derivative signal can (only) be observed around the aperture boundary (Fig. 1(c)), which is consistent with our theoretical analysis (Eq. (11)). With use of the proposed method (Algorithm 2), the phase function was retrieved accurately with the RMSE 0.91% (Figs. 1(d) and 1(e)). Since only one DCT and inverse DCT is needed to retrieve the phase, the algorithm can be performed extremely fast. The computation time in this example was only 16 ms using a 2.5 GHz laptop in the MATLAB environment. In contrast, the direct computation of the convolution integral (Eq. (20)) requires a large amount of memory to store the Green’s function, therefore limits the maximum size of the phase matrix allowed to only 50 × 50 pixels. Even for such low resolution data, the direct computation took almost 12 seconds.

Fig. 1 Simulation results for the phase function with zero Laplacian. (a) True phase distribution. (b) Longitudinal intensity derivative distribution in the absence of the aperture. (c) Longitudinal intensity derivative distribution with the aperture. (d) Reconstructed phase with the proposed algorithm. (e) Residual error map (RMSE 0.91%). The white dashed box in (a) outlines the aperture edge.

Download Full Size | PDF

5.2 Phase retrieval with non-uniform intensity distribution

To test the proposed algorithm in the case of non-uniform intensity distribution, we further consider the same phase map but with the non-uniform intensity distribution $I (x, y) = \exp [- (x^{2} + y^{2}) / 2 σ^{2}]$ with $σ$ = 0.18 mm (Fig. 2(a)). In this case, the intensity derivative signal obtained is non-zero in the absence of the aperture, but still quite weak, as shown in Fig. 2(b). When the aperture was introduced, boundary signal can be observed at the aperture edge, which is much stronger compared with the derivative signal inside the aperture (Fig. 2(c)). Again, the phase distribution was faithfully retrieved by the proposed method (Algorithm 1) (Fig. 2(d)), with the residual error plotted in Fig. 2(e). The RMSE of the reconstructed phase is 2.23%, which is slightly higher than the uniform intensity case. The increased error may result from the “Teague’s assumption” (see Section 6.2 for details). The runtime of the algorithm in this example was 66 ms, which is nearly four times as the uniform intensity case (Algorithm 2). However, if one simply assumes constant intensity and applies Algorithm 2 to the same dataset when the intensity is in fact non-uniform, the resultant phase error will be much higher, as shown in Fig. 2(f). This further verifies that the intensity variations cannot be neglected even its distribution is relatively smooth.

Fig. 2 Reconstruction results for the zero-Laplacian phase function with non-uniform intensity distribution. (a) Simulated non-uniform intensity distribution. (b) Longitudinal intensity derivative distribution without the aperture. (c) Longitudinal intensity derivative distribution with the aperture. (d) Reconstructed phase with the proposed algorithm (using Algorithm 1 or 3). (e) Residual error map. (f) Residual error when the intensity distribution is assumed to be uniform (using Algorithm 2 or 4). The white dashed box in (a) outlines the aperture edge.

Download Full Size | PDF

To compare the proposed method with previously developed techniques, we applied the FFT-based method developed by Paganin and Nugent [14] and the FFT-based symmetrization method developed by Volkov et.al [15] to the same dataset but without the presence of the aperture. Besides, for the FFT-based method, we also tested it with the same dataset as our DCT-based algorithm in the presence of the aperture. The RMSEs of different methods are summarized in the first row (Dataset 1) of Table 1.It can be seen that only the proposed algorithm reconstructed the phase reliably, while all other methods failed to retrieve the phase function with zero Laplacian correctly, leading to very large RMSE values.

Table 1. RMSEs of the Phases Reconstructed by Different Solutions to the TIE

View Table | View all tables in this article

5.3 Complex phase and intensity distribution

The next simulation focuses on more complex phase map and intensity distribution. We consider the phase profile of letters ‘TIE’ (Fig. 3(a)) with the non-uniform intensity distribution of letters ‘COLE’ (Fig. 3(b)). Figures 3(c) and 3(d) plot the longitudinal intensity derivative distribution without/with the presence of the aperture, respectively. The two derivative images look quite similar at first glance. However, a closer inspection on the enlarged area (Figs. 3(e) and 3(f)) reveals that the aperture generated additional boundary signals at the aperture edge which cannot be observed for the no aperture case.

Fig. 3 Simulation result of complex intensity and phase distribution. (a) True phase map. (b) Simulated intensity map. (c) Longitudinal intensity derivative without the aperture. (d) Longitudinal intensity derivative with the aperture. (e) Enlarged region corresponding to the boxed area of (c). (f) Enlarged region corresponding to the boxed area of (d). The white dashed box in (a) outlines the aperture edge, which defines the integral region $\bar{Ω}$ shown in (f).

Download Full Size | PDF

Figure 4 gives the pictorial comparison of different solutions and the RMSE values are presented in the second row (Dataset 2) of Table 1. It can be seen that the conventional FFT method could not retrieve the phase correctly no matter with (Figs. 4(b) and 4(f)) or without (Figs. 4(c) and 4(d)) the aperture since it simply assumes periodic boundary conditions, which is of course not true in this situation. Furthermore, the boundary error of the even symmetrization method was also significant in the absence of the aperture (Figs. 4(d) and 4(h)). In contrast, the proposed method provided a faithful result with much higher accuracy (RMSE 2.75%), proving that boundary errors can only be eliminated by the proposed DCT-based method with the practical boundary value information obtained on the aperture edges.

Fig. 4 Comparison of different TIE solutions. The first row shows the recovered phase by (a) the proposed algorithm with the aperture, (b) FFT-based algorithm with the aperture, (c) FFT-based algorithm without the aperture, (d) even symmetrization method without the aperture. The second row (e)-(f) shows the corresponding residual errors for different methods.

Download Full Size | PDF

5.4 Phase retrieval under noise and large defocus distance

In the last simulation, we tested the accuracy and stability of the proposed method in the presence of noise. Different levels of pseudo-random Gaussian noise with standard deviations (STD) of 0.0001, 0.001 and 0.01 were generated and superimposed on each intensity image of the data set from Section 5.3. To clearly demonstrate the noise levels, we show the calculated noisy intensity derivatives in Figs. 5(a)–5(c) (the defocus distance is $Δ z = \pm 10$ μm, same as the previous simulations). The reconstructed phases and corresponding residual errors are shown in Figs. 5(d)–5(f) and Figs. 5(g)–5(i), respectively. The RMSE values of the reconstructed phases are further given in the first row (defocus distances 10μm) of Table 2It can be seen that when the noise level was low, the current method could still recover the phase fairly accurately, with a slight increase in RMSE, from 2.23% to 3.14%. However, with increase of noise level, low-frequency artifacts can be clearly observed from the recovered phases. Especially when the noise STD reached 0.01, the original phase signal was almost submerged by the cloud-like artifacts, leading to significant error (RMSE 213.44%).

Fig. 5 Phase reconstruction results by the proposed method at different noise levels when defocus distance is 10μm. The first row shows the longitudinal intensity derivative distribution with noise stand deviation (a) 0.0001, (b) 0.001, and (c) 0.01. The second row shows the reconstructed phases at different noise levels (d)-(f). The third row shows the corresponding residual error at different noise levels (g)-(i).

Download Full Size | PDF

Table 2. RMSEs of the Phases Reconstructed by the Proposed Algorithm Under Different Noise Levels and at Difference Defocus Distances

View Table | View all tables in this article

The high sensitivity to low-frequency phase error is a common and notorious problem in the TIE because the inverse Laplacian operator is ill-conditioned at near zero-frequency [35]. The common countermeasures include regularization treatments [4], using more intensity measurements [21, 22], or simply increasing the defocus distance [35]. By increasing the defocus distance $Δ z$ from 10μm to 100μm, the derivative SNR is almost increased by 10 times since the image noise is amplified by the inverse of defocus distance $Δ z$ in the derivative estimates, as verified by comparing the two intensity derivative distributions shown in Figs. 5(b) and 6(a).However, on the other hand, increasing $Δ z$ also compromised the linearity assumption in the finite-difference derivative approximation (Eq. (32)), leading to non-linearity errors [21, 22]. From the enlarged region (Fig. 6(b)), we can clearly see the estimated derivative became blurred compared with the more accurate one obtained with $Δ z$ = 10μm (Fig. 3(f)). Moreover, some boundary signals escaped from the aperture edge and went into the exterior region, suggesting Eq. (11) is no longer well hold. In this case, a new integration region $Ω^{+}$ should be defined to make sure all the escaped signals can be enclosed in to fulfill the energy conservation law (Eq. (13)). Obviously, the choice of region $Ω^{+}$ satisfying Eq. (13) is not unique and one can even extend $Ω^{+}$ to the unbounded space. However, since the Green's function $G (r, r^{'})$ in the integral (Eqs. (20) and (23)) is only defined on the aperture region; we should choose $Ω^{+}$ as close to $\bar{Ω}$ as possible to minimize this region mismatch. Therefore, a simple searching process was used to quickly define the appropriate integration region $Ω^{+}$ by checking the integral value (Eq. (13)) starting from the initial region $Ω^{+} = \bar{Ω}$ and then gradually expanded outwards until it is sufficiently close to zero. The new integral region $Ω^{+}$ (190 × 190 pixels) found by this procedure is indicated by the solid black line in Fig. 6(b), which is slightly larger than the original aperture region (180 × 180 pixels). Since the intensity $I (r)$ is defined only on the closure $\bar{Ω}$ , the gap (intensity zeros) region between $\bar{Ω}$ and $Ω^{+}$ can duplicate the intensity values of the outermost pixels within the aperture region (or simply restrict the domain of phase integration in Eq. (23) to $\bar{Ω}$ only). The reconstructed phase and the residual error are shown in Fig. 6(c) and 6(d), respectively, which demonstrates a great improvement in phase accuracy compared with the result obtained when $Δ z$ = 10μm (Fig. 5(e)), with the RMSE value decreasing from 29.38% to 9.30%. But for the high-noise-level condition (noise STD 0.01), we need to further increase the defocus distance $Δ z$ to 500μm. In this case, although the SNR of the derivative signal was improved greatly compared with Fig. 5(c), the nonlinear effect was significant – the estimated derivative became quite blurry, as shown from the enlarged region (Fig. 6(f)). It is obvious that we cannot expect the result to be very accurate for such a situation. Followed by a similar procedure, we found the integration region $Ω^{'}$ in this case should be extended to 198 × 198 pixels indicated by the solid black line in Fig. 6(f). Once again, the reconstructed phase demonstrates a significant reduction in low-frequency artifacts, with the RMSE value dropping from 213.44% to 17.92%. However, the blurring effects can be easily perceived in the reconstructed phase and residual error map, which means the non-linearity error was quite dominant and the $\partial I / \partial z$ obtained by Eq. (32) was no longer a valid estimate of axial intensity derivative. However, even under such conditions, the RMSE of our algorithm is still much lower than all the other algorithms under ideal, noise-free conditions. These results confirm the stability of our algorithm with respect to the noise, and its validity under large defocus distances and nonlinearity errors.

Fig. 6 Phase reconstruction results by the proposed method when the defocus distance is not small. The first row shows the (a) longitudinal intensity derivative distribution, (b) enlarged region corresponding to the boxed area of (a), (c) reconstructed phase map and (d) residual error when the noise standard deviation is 10⁻³, and the defocus distance is 100μm. The second row (e)-(h) shows same figures when the noise standard deviation is 10⁻², and the defocus distance is 500μm.

Download Full Size | PDF

6. Discussions

6.1 Relation between FFT-based method

The main idea of the present method is to consider the solution of the TIE as an inhomogeneous Neumann boundary value problem - with the Neumann boundary value experimentally measurable around a hard-edged aperture. The validity of boundary value is guaranteed by the energy conservation law so that our method represent a more general case and can be applied to any complex object without imposing any implicit or prior assumptions. Conversely, as we discussed in Sec. 2.3, if we are able to make some safe assumptions on the test object or experimental configuration such as constant phase at the boundary (which are a quite common assumptions to ensure the validity the FFT-based methods), the boundary measurements and the aperture become unnecessary because there would be no boundary signal at the aperture edge (see Eq. (9)). In this case, one can define not only the Neumann boundary, but any uniform Dirichlet boundary conditions or periodic ones as well. The eigen-functions of the Laplacian can then be more general Fourier harmonics [33] (can be cosine, sine, or exponential harmonics), and thus there will be no essential difference between the proposed DCT-based solution and the traditional FFT-based solution.

6.2 Inaccuracy results from “Teague’s assumption”

It has to be mentioned that the TIE solution presented here is based on the “Teague’s assumption” that there exists an auxiliary function $ψ$ satisfying Eq. (14). However, as some previous work suggested, Teague’s assumption may not be always valid [9, 40]. According to the Helmholtz’s theorem, the vector field $I \nabla ϕ$ can be decomposed as

I \nabla ϕ = \nabla ψ + r o t η,

where

r o t η

is the vector field

r o t η = (\partial η / \partial y, \partial η / \partial x)

and

η

is a scalar function. Compared with Eq. (14), it can be seen that the rotational term

r o t η

was simply ignored in Teague’s assumption [9], postulating that

I \nabla ϕ

must be potential and irrotational. A detailed discussion about the validity and the inaccuracy resulting from “Teague’s assumption” can be found in [40]. Using the identity

\nabla \cdot r o t η = 0

, the first Poisson equation (Eq. (15)) is still valid and this approximation does not affect accurate reconstruction of function

ψ

. However, the missing rotational term does affect the solution of the second Poisson equation (Eq. (21))

\nabla^{2} ϕ = \nabla \cdot (I^{- 1} \nabla ψ) + \nabla \cdot (I^{- 1} r o t η)

It should be noted that the second term

\nabla \cdot (I^{- 1} r o t η)

may be generally non-zero thus error will be introduced by simply ignoring it. However, in general the contribution of the rotational term

r o t η

is quite small [40]. As is shown in our simulations, inaccuracies due to “Teague’s assumption” are much smaller than the error introduced by incorrect boundary conditions. Thus this assumption can be regarded as a reasonable one in our work. Besides, the small phase error can be further compensated by an iterative process considering the rotational term as a perturbation term [41].

7. Conclusions

In conclusion, we have presented a method for TIE phase retrieval with use of DCT. The new method has several substantial advantages: First, it clearly defines the TIE as an inhomogeneous Neumann boundary value problem with the boundary value experimentally measurable around a hard-edged aperture. Second, the method is valid for the case of non-uniform intensity distribution. Third, the measured intensity data are treated as one entity, without requiring special-purpose detection scheme to separate boundary values from the interior signals. Finally, for rectangular region, the DCT based numerical implementation is quite simple, fast, and of low memory usage. Its good performance has been confirmed by several numerical simulations even when the objects are located at the image borders. Those precise results are clearly essential when dealing with quantitative phase imaging, which was the main driver for this research. The method has already been tested experimentally with the application to micro-optics characterization and found to be robust to many source of experimental error. The experimental results will be published in a subsequent paper.

Appendix A: Cosine series expansion and discrete cosine transform

The preceding Appendixes served to motivate the development of DCT as an efficient tool for solving the two Poisson equations, i.e., calculating the two area integrals of Eqs. (20) and (23) related to the TIE solution. Considering a band limited continuous signal $f (x)$ defined on $[0, a]$ , it can be expended into an infinite but convergent series of cosine basis functions

f (x) = \sum_{k = 0}^{\infty} F_{k} \cos (\frac{π k}{a} x) .

The coefficient for this expansion can be represented as

F_{k} = a_{k} \int_{0}^{a} f (x) \cos (\frac{π k}{a} x) d x,

with the normalization factors

a_{k} = {\begin{cases} \frac{2}{a}, k \neq 0 \\ \frac{1}{a}, k = 0 \end{cases} .

It should be noted the basis functions of cosine series expansion

Ψ_{k} (x) = \cos (π k x / a)

is the eigen-function of Laplacian with the Neumann boundary conditions in one-dimension (1-D) case. Without loss of generality, we sample

f (x)

at a uniform interval to produce data sequence

{f (n), n = 0, 1, ..., N - 1}

. The cosine series expansion then becomes DCT accordingly [42]. The 1-D DCT of

f (n)

is defined as

F (k) = A_{k} \sum_{n = 0}^{N - 1} f (n) \cos [\frac{π k}{N} (n + 0.5)]; k = 0, 1, ..., N - 1.

A_{k} = {\begin{cases} \frac{\sqrt{2}}{\sqrt{N}}, k \neq 0 \\ \frac{1}{\sqrt{N}}, k = 0 \end{cases} .

Its inverse transform is

x (n) = \sum_{k = 0}^{N - 1} A_{k} X (k) \cos [\frac{π k}{N} (n + 0.5)]; n = 0, 1, ..., N - 1.

Note that correct to a scaling factor

A_{k}

, the forward and backward (inverse) transforms have identical transformation kernels. The DCT basis vectors are

Ψ_{k} (n) = {A_{k} \cos [\frac{π k}{N} (n + 0.5)]}_{n = 1, 2, ..., N - 1}; k = 0, 1, ..., N - 1.

DCT is an orthonormal transform since the half point shift of coordinate making the cosine base function to be an orthogonal, complete set over the N-dimensional vector space and meanwhile fulfills the Neumann at the region boundary.

Similar to the cosine transform of one dimension, for the two-dimensional (2-D) data sequence ${f (x, y), x = 0, 1, ..., N - 1, y = 0, 1, ..., M - 1}$ , their DCT and inverse transform are defined respectively as:

F (m, n) = A_{m} A_{n} \sum_{x = 0}^{M - 1} \sum_{y = 0}^{N - 1} f (x, y) \cos [\frac{π m}{M} (x + 0.5)] \cos [\frac{π n}{N} (y + 0.5)];

f (x, y) = \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} A_{m} A_{n} F (m, n) \cos [\frac{π m}{M} (x + 0.5)] \cos [\frac{π n}{N} (y + 0.5)] .

Note that the transformation kernels of 2-D DCT are separable, so that the 2-D transform is equivalent to a 1-D DCT performed along a single dimension followed by another 1-D DCT in the other dimension.

Appendix B: Calculate DCT based on FFT

The direct implementation of DCT involves four nested loops with complexity $O (M^{2} N^{2})$ for an $M \times N$ array, which is prohibitive in most applications. As a result, various fast and efficient DCT algorithms have been devised and studied extensively [43]. However, here we only discuss one fast algorithm that utilizes the FFT to compute the DCT and its inverse [44] due to its simple implementation and relevance to one previous work [15]. The fast computation procedure consists of extending the input sequence of N samples to an 2N-point sequence with even symmetry, taking an 2N-point DFT, and saving N terms from it. The even extension is defined as

\bar{f} (n) = {\begin{cases} f (n) n = 0, 1, ..., N - 1 \\ f (2 N - 1 - n) n = N, N + 1, ..., 2 N - 1 \end{cases} .

The 2N-point DFT vector of

\bar{f} (n)

is

\bar{F} (k) = \frac{1}{\sqrt{2 N}} \sum_{n = 0}^{2 N - 1} \bar{f} (n) e^{- j (\frac{2 π k n}{2 N})} = \frac{1}{\sqrt{2 N}} e^{j (\frac{π k}{N})} \sum_{n = 0}^{N - 1} f (n) \cos [\frac{π k}{N} (n + 0.5)], k = 0, 1, ..., 2 N - 1.

Comparing it with the DCT coefficients

F (k)

defined in Eq. (39), we get

F (k) = \frac{A_{k}}{\sqrt{2 N}} e^{- j (\frac{π k}{2 N})} \bar{F} (k); k = 0, 1, ..., N - 1.

Similarly, it can be deduced that for the inverse DCT, one can recover the 2N-point DFT coefficients

\bar{F} (k)

from the DCT coefficients

F (k)

as follows:

\bar{F} (k) = {\begin{cases} \frac{A_{k}}{\sqrt{2 N}} e^{j (\frac{π k}{2 N})} F (k); k = 0, 1, ..., N - 1 \\ 0; k = 0 \\ - \frac{A_{k}}{\sqrt{2 N}} e^{j (\frac{π k}{2 N})} F (2 N - k); k = N + 1, ..., 2 N - 1 \end{cases} .

Therefore, computing the 2N-point inverse DFT, the reconstructed data by inverse DCT can be obtained in the first N elements of the inverse DFT result.

f (n) = \bar{f} (n); n = 0, 1, ..., N - 1,

For simplicity, we denote the even extension and the corresponding interception operation as the bar (‘-’) and angle brackets (‘< >’) respectively. This algorithm is presented here for the 1-D DCT, nevertheless it can be easily extended to the 2-D case (by extending the input matrix both in

x

and

y

coordinates) due to the separability of the 2D-FFT and 2D-DCT.

Appendix C: Derivative calculation based on DCT

In this section, expressions are derived for computing derivative operations on a signal, from the DCT coefficients of the original signal’s samples. Applying the derivative operation to both sides of Eq. (36), and adopting the notation $Ψ_{k}^{(l)} (x)$ mean the lth order derivative of the basis functions of the transform. Then lth order derivative of $f$ can be deduced as

f^{(l)} (x) = \sum_{k = 0}^{\infty} a_{k} F_{k} Ψ_{k}^{(l)} (x) = \sum_{k = 0}^{\infty} a_{k} F_{k} \cos^{(l)} (\frac{π k}{a} x) .

Similarly, for discrete data let us adopt the notation

f^{(l)} (n)

to mean the discrete derivative of

f^{(l)} (x)

at

x = n

, and we can use it as an approximation of the original derivative. Then the 1st and 2nd order derivative calculated with DCT can be represented as [45]

f^{'} (n) = - \sum_{k = 0}^{N - 1} A_{n} (\frac{π k}{a}) F (k) \sin [\frac{π k}{N} (n + 0.5)];

f^{″} (n) = - \sum_{k = 0}^{N - 1} A_{n} {(\frac{π k}{a})}^{2} F (k) \cos [\frac{π k}{N} (n + 0.5)];

where

F (k)

is the DCT coefficients of

f (n)

. Higher order derivatives can be analogized. It is seen that for the 2nd derivative, one only needs to reweight the DCT coefficients with a multiplicative term

- {(π k / a)}^{2}

and then reconstruct it by inverse DCT. However, for 1st derivative, not only the multiplicative term

- (π k / a)

is induced, but also the basis of the obtained derivative are changed. So one could not reconstruct the modified coefficients via inverse DCT. Analogous rules apply to rectangular domains in two dimensions. To compute a single derivative like

\partial f (x, y) / \partial x

or

\partial^{2} f (x, y) / \partial x^{2}

, we can apply Eq. (51) or Eq. (52) along the

x

direction, once for each

y

:

\frac{\partial f (x, y)}{\partial x} = - \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} A_{m} A_{n} (\frac{π m}{a}) F (m, n) \sin [\frac{π m}{M} (x + 0.5)] \cos [\frac{π n}{N} (y + 0.5)];

\frac{\partial^{2} f (x, y)}{\partial x^{2}} = - \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} A_{m} A_{n} {(\frac{π m}{a})}^{2} F (m, n) \cos [\frac{π m}{M} (x + 0.5)] \cos [\frac{π n}{N} (y + 0.5)] .

Similar equations applies for derivatives in the

y

- direction

\frac{\partial f (x, y)}{\partial y} = - \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} A_{m} A_{n} (\frac{π n}{b}) F (m, n) \cos [\frac{π m}{M} (x + 0.5)] \sin [\frac{π n}{N} (y + 0.5)];

\frac{\partial^{2} f (x, y)}{\partial y^{2}} = - \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} A_{m} A_{n} {(\frac{π n}{b})}^{2} F (m, n) \cos [\frac{π m}{M} (x + 0.5)] \cos [\frac{π n}{N} (y + 0.5)] .

We denote the above formulations in more convenient forms:

Eq . (C4) - \partial_{x - D C T} f

,

Eq . (C5) - \partial_{x - D C T}^{2} f

,

Eq . (C6) - \partial_{y - D C T} f

,

Eq . (C7) - \partial_{y - D C T}^{2} f

. With these basic elements, we can easily describe the Hamilton

\nabla

with DCT-based operators

\nabla_{D C T} = (\partial_{x - D C T}, \partial_{y - D C T}) .

By applying product to scalars or dot product to vectors,

\nabla_{D C T}

may denote the gradient of a scalar field and the divergence of a vector field:

\nabla_{D C T} f = (\partial_{x - D C T} f, \partial_{y - D C T} f),

\nabla_{D C T} \cdot f = \partial_{x - D C T} f_{x} + \partial_{y - D C T} f_{y} .

It is more interesting to consider the Laplacian operator

\begin{array}{l} \nabla_{D C T}^{2} f (x, y) = \partial_{x - D C T}^{2} f + \partial_{y - D C T}^{2} f = \\ - \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} A_{m} A_{n} [{(\frac{π m}{a})}^{2} + {(\frac{π n}{b})}^{2}] F (m, n) \cos [\frac{π m}{M} (x + 0.5)] \cos [\frac{π n}{N} (y + 0.5)], \end{array}

from which it is easy to find that the multiplicative term above is exactly the Laplacian eigenvalues shown in Eq. (27). We further denote the inverse version of Eq. (61) as

\nabla_{^{D C T}}^{^{- 2}} f

:

\nabla_{^{D C T}}^{^{- 2}} f = - \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} A_{m} A_{n} {[{(\frac{π m}{a})}^{2} + {(\frac{π n}{b})}^{2}]}^{- 1} F (m, n) \cos [\frac{π m}{M} (x + 0.5)] \cos [\frac{π n}{N} (y + 0.5)] .

Or more compactly,

\nabla_{^{D C T}}^{^{- 2}} f = D C T^{- 1} λ_{m, n}^{- 1} D C T (f) .

Appendix D: Calculate DCT-based derivatives using FFT

In appendix C, we present several formulations for differentiation differential operators based on DCT. It has been found that for operators related to only even order derivatives, such as $\partial_{x - D C T}^{2} f$ and $\nabla_{^{D C T}}^{^{- 2}} f$ , the reconstruction is a simple inverse DCT. However, this rule does not work for odd order derivatives. The direct calculation of Eqs. (55) and (56) requires expensive computation $O (M^{2} N^{2})$ for an $M \times N$ matrix. In Appendix B, we demonstrate a handy fast algorithm to compute DCT based on 2N point FFT. Accordingly, here we present the corresponding fast algorithms to calculate DCT-based differential operators using FFT.

We again consider the 2N-point sequence $\bar{f} (n)$ extended from the extending the input N samples according to Eq. (45) with Fourier coefficients $\bar{F} (k), k = 0, 1, ..., 2 N - 1$ defined in Eq. (46). Using the relationship between DCT coefficients and the 2N-point DFT coefficients (Eqs. (47) and (48)), we can deduce that if we modify the $\bar{F} (k)$ as the following rules,

\bar{F} (k) = {\begin{cases} \frac{2 π j k}{a} F (k); k = 0, 1, ..., N \\ \frac{2 π j}{a} (k - 2 N) F (k); k = N, ..., 2 N - 1 \end{cases} .

Computing the 2N-point inverse FFT, then the first N elements of the reconstruction result corresponds to the sampled derivative calculated by DCT-based method (Eq. (51)).

f^{'} (n) = {\bar{f}}^{'} (n); n = 0, 1, ..., N - 1.

If the

\bar{F} (k)

is modified as

\bar{F} (k) = {\begin{cases} - {(\frac{2 π k}{a})}^{2} F (k); k = 0, 1, ..., N - 1 \\ - {[\frac{2 π (k - 2 N)}{a}]}^{2} F (k); k = N, ..., 2 N - 1 \end{cases},

then the first N elements in the reconstruction result by FFT corresponds to the sampled second order derivative calculated by DCT-based method (Eq. (52)).

f^{″} (n) = \bar{f^{″}} (n); n = 0, 1, ..., N - 1.

It is not hard to find the weighting coefficients in Eqs. (63) and (65) make no difference if one wants to compute the derivative of

\bar{f} (n)

use the DFT or FFT based method without considering that

\bar{f} (n)

is an extended version of

f (n)

[46]. Similar relations can be found in 2-D cases: to calculate

\partial_{x - D C T} f

,

\partial_{y - D C T} f

, and

\nabla_{^{D C T}}^{^{- 2}} f

, one can just extend the input matrix

f

both in

x

and

y

coordinates symmetrically, then perform the same procedure as in DFT or FFT based derivative calculation methods, and finally intercept the part corresponding to its original size from the reconstructed result. By denoting the symmetrical extension and the corresponding interception operation as the bar (‘-’) and angle brackets (‘< >’) respectively, the fast algorithms for some 2-D differential operators can be presented as

\partial_{x - D C T} f = 〈 \partial_{x - F F T} \bar{f} 〉,

\partial_{y - D C T} f = 〈 \partial_{y - F F T} \bar{f} 〉,

\nabla_{^{D C T}}^{^{- 2}} f = 〈 \nabla_{^{F F T}}^{^{- 2}} \bar{f} 〉,

where

\partial_{x - F F T} f

,

\partial_{y - F F T} f

, and

\nabla_{^{F F T}}^{^{- 2}} f

represent their counterpart implemented based on FFT.

Acknowledgments

We are grateful to the anonymous reviewers for their constructive suggestions. This project was supported by the Research Fund for the Doctoral Program of Ministry of Education of China (No. 20123219110016), 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. F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt. 27(7), 1223–1225 (1988). [CrossRef] [PubMed]

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

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

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

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

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

8. S. C. Woods and A. H. Greenaway, “Wave-front sensing by use of a Green’s function solution to the intensity transport equation,” J. Opt. Soc. Am. A 20(3), 508–512 (2003). [CrossRef] [PubMed]

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

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

11. T. Gureyev, A. Roberts, and K. 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]

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

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

14. D. Paganin and K. A. Nugent, “Noninterferometric Phase Imaging with Partially Coherent Light,” Phys. Rev. Lett. 80(12), 2586–2589 (1998). [CrossRef]

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

16. F. Roddier, “Wavefront sensing and the irradiance transport equation,” Appl. Opt. 29(10), 1402–1403 (1990). [CrossRef] [PubMed]

17. I. W. Han, “New method for estimating wavefront from curvature signal by curve fitting,” Opt. Eng. 34(4), 1232–1237 (1995). [CrossRef]

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

19. E. D. Barone-Nugent, A. Barty, and K. A. Nugent, “Quantitative phase-amplitude microscopy I: optical microscopy,” J. Microsc. 206(3), 194–203 (2002). [CrossRef] [PubMed]

20. S. S. Kou, L. Waller, G. Barbastathis, and C. J. R. Sheppard, “Transport-of-intensity approach to differential interference contrast (TI-DIC) microscopy for quantitative phase imaging,” Opt. Lett. 35(3), 447–449 (2010). [CrossRef] [PubMed]

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

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

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

24. 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).

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

26. D. A. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Second Order (Springer, 2001), Vol. 224.

27. K. Bhamra, Partial Differential Equations: An Introductory Treatment With Applications (PHI, 2010).

28. R. Courant and D. Hilbert, “Potential theory and elliptic differential equations,” in Methods of Mathematical Physics (Wiley-VCH Verlag GmbH, 2008), pp. 240–406.

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

30. D. W. Trim, Applied Partial Differential Equations (PWS-Kent, 1990).

31. A. Agrawal, R. Raskar, and R. Chellappa, “What is the range of surface reconstructions from a gradient field?” in Computer Vision–ECCV 2006 (Springer, 2006), pp. 578–591.

32. A. Talmi and E. N. Ribak, “Wavefront reconstruction from its gradients,” J. Opt. Soc. Am. A 23(2), 288–297 (2006). [CrossRef] [PubMed]

33. F. Morse, Methods of Theoretical Physics (1981), Vols. 1 and 2.

34. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran 77: The Art of Scientific Computing (Cambridge University, 1992).

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

36. C. Dorrer and J. D. Zuegel, “Optical testing using the transport-of-intensity equation,” Opt. Express 15(12), 7165–7175 (2007). [CrossRef] [PubMed]

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

38. C. Campbell, “Wave-front sensing by use of a Green’s function solution to the intensity transport equation: comment,” J. Opt. Soc. Am. A 24(8), 2480–2482 (2007). [CrossRef] [PubMed]

39. S. C. Woods, H. I. Campbell, and A. H. Greenaway, “Wave-front sensing by use of a Green's function solution to the intensity transport equation: reply to comment,” J. Opt. Soc. Am. A 24(8), 2482–2484 (2007). [CrossRef]

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

41. V. V. Volkov and Y. Zhu, “Lorentz phase microscopy of magnetic materials,” Ultramicroscopy 98(2-4), 271–281 (2004). [CrossRef] [PubMed]

42. N. Ahmed, T. Natarajan, and K. R. Rao, “Discrete cosine transform,” IEEE Trans. Comput. C-23(1), 90–93 (1974). [CrossRef]

43. E. Feig and S. Winograd, “Fast algorithms for the discrete cosine transform,” IEEE Trans. Signal Process. 40(9), 2174–2193 (1992). [CrossRef]

44. J. Tribolet and R. E. Crochiere, “Frequency domain coding of speech,” IEEE Trans. Acoust. Speech Signal Process. 27(5), 512–530 (1979). [CrossRef]

45. R. Reeves and K. Kubik, “Shift, scaling and derivative properties for the discrete cosine transform,” Signal Process. 86(7), 1597–1603 (2006). [CrossRef]

46. R. N. Bracewell, The Fourier Transform and Its Applications (McGraw-Hill, 1986), Vol. 3.

Dataset	Proposed (Algorithm 1)	FFT-based method (with aperture)	FFT-based method	FFT with odd symmetrization	FFT with even symmetrization
1	2.23%	77.80%	89.54%	71.34%	96.18%
2	2.75%	19.25%	131.49%	33.83%	493.51%

Dataset	Proposed (Algorithm 1)	FFT-based method (with aperture)	FFT-based method	FFT with odd symmetrization	FFT with even symmetrization
1	2.23%	77.80%	89.54%	71.34%	96.18%
2	2.75%	19.25%	131.49%	33.83%	493.51%

Boundary-artifact-free phase retrieval with the transport of intensity equation: fast solution with use of discrete cosine transform

Abstract

1. Introduction

2. Transport of intensity equation, boundary conditions, and uniqueness of the solution

2.1 Boundary conditions and uniqueness of the solution

2.2 Compatibility condition and energy conservation law

2.3 Solving the TIE without any boundary values

3. Solution to transport of intensity under non-uniform intensity

3.1 Boundary values generation with a hard-edged aperture

3.2 Analytic solution to the TIE with the Green’s function approach

4. Fast numerical implementation with use of discrete cosine transform

4.1 Eigen-functions expansion of the Green's function

4.2 Numerical implementation with fast discrete cosine transform

5. Simulations and comparisons

5.1 Pure phase object with uniform intensity distribution

5.2 Phase retrieval with non-uniform intensity distribution

5.3 Complex phase and intensity distribution

5.4 Phase retrieval under noise and large defocus distance

6. Discussions

6.1 Relation between FFT-based method

6.2 Inaccuracy results from “Teague’s assumption”

7. Conclusions

Appendix A: Cosine series expansion and discrete cosine transform

Appendix B: Calculate DCT based on FFT

Appendix C: Derivative calculation based on DCT

Appendix D: Calculate DCT-based derivatives using FFT

Acknowledgments

References and links

Cited By

Figures (6)

Tables (6)

Equations (69)

Optics Express

Defocus distances(μm)	Low noise level (noise STD 10⁻⁴)	medium noise level (noise STD 10⁻³)	High noise level (noise STD 10⁻²)
10	3.14%	29.38%	213.44%
100	8.88%	9.30%	38.75%
500	17.21%	17.32%	17.92%