The transport of intensity equation (TIE) is a two-dimensional second order elliptic partial differential equation that must be solved under appropriate boundary conditions. However, the boundary conditions are difficult to obtain in practice. The fast Fourier transform (FFT) based TIE solutions are widely adopted for its speed and simplicity. However, it implies periodic boundary conditions, which lead to significant boundary artifacts when the imposed assumption is violated. In this work, TIE phase retrieval is considered as an inhomogeneous Neumann boundary value problem with the boundary values experimentally measurable around a hard-edged aperture, without any assumption or prior knowledge about the test object and the setup. The analytic integral solution via Green’s function is given, as well as a fast numerical implementation for a rectangular region using the discrete cosine transform. This approach is applicable for the case of non-uniform intensity distribution with no extra effort to extract the boundary values from the intensity derivative signals. Its efficiency and robustness have been verified by several numerical simulations even when the objects are complex and the intensity measurements are noisy. This method promises to be an effective fast TIE solver for quantitative phase imaging applications.
© 2014 Optical Society of America
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  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 , X-ray diffraction , electron- microscopy  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 . 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 , 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  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  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 . 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.  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 axis, and let the complex amplitude of the object be . The derivative of intensity along the beam propagation direction, contains phase information that can be retrieved via TIE :Eq. (1), we obtain
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 with a piecewise smooth boundary . The intensity distribution is a continuous, non-negative function defined on the enclosure (consisting of the domain plus its boundary ) and is smooth and strictly positive in . The longitudinal intensity derivative is assumed to be a smooth function in . The phase is expected to be single-valued and smooth in . 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 . 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 boundary26]. 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 and , therefore if is a solution, then is also a solution for any constant . 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 .
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 conditionEq. (5) is known a priori, since the function 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 and , 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) becomesEqs. (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 unnecessaryEq. (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 ) 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 ) or periodic ones (the phase at the boundary repeats cyclically) as well. Gureyev and Nugent  suggested another way to eliminate the need of boundary conditions (LHS of Eq. (5)) by considering the case that the intensity inside the domain but strictly vanishes at the boundary (actually it corresponds to a special case of the homogeneous Neumann boundary conditions though not mentioned in ). Alternatively, without any additional requirement about the test object and experimental conditions, Volkov et al.  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, . 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, . 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 (): the intensity is smooth (but can be non-uniform) in but suddenly vanishes with a step-like discontinuity at the boundary, . Physically, it is equivalent to capturing intensity images through a hard-edged aperture at the object plane:2] as a plausible representation of the physical quantity concentrated on the boundary curve. The can be loosely thought of as the gradient magnitude of the aperture function or interpreted as the limit of a sequence of smooth functions. Mathematically, the function is an extension of the Dirac delta function at single point , which can be rigorously defined as a distribution such that for any test function that is smooth and compactly supportedEqs. (8) and (9) into the TIE (Eq. (2)), we obtainEq. (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  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 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 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 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 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 obtainEq. (6)), the LHS of the above equation should equal to zeroEq. (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 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 , an auxiliary function was introduced, which satisfies14, 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]Eq. (14). The Green's function representation of the solution to the Neumann boundary-value problem (Eqs. (14) and (16)) is given by30]. 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 and then integrate it over all space, we obtain:Eq. (19) with Eq. (17), we can rewrite the Green's function solution into a more compact area integral:Equation (20) is the central result of this section. Here we should reemphasize that the region 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 . However, this rule cannot be applied to our case since the integrand in Eq. (19) contains the generalized function . As shown in Eq. (10), the generalized function 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 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 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 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 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 equationEqs. (21) and (22)) isEq. (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 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 pixel image, will contain elements, which may require thousands of gigabyte memory to store. Besides, the direct numerical implementation of the integral requires expensive operations of 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 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], the eigen-functions form a complete and orthogonal basis over the domain . More specifically, for a rectangular domain, the eigen-functions can be represented as the Fourier cosine harmonicsEqs. (24)–(27) into Eq. (20), and interchanging the order of summation and integral, we obtainEq. (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 , and finally the function 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 is discretized into a finite rectangular grid cells of equal size. For convenience, in this subsection we still use the to represent the image traverse coordinate, but keep in mind they can only assume integral values, i.e., . The pixel coordinate corresponds the grid cell centered at the location in the continuous real space within the rectangular domain , where , is the cell interval (pixel size) in the and direction, respectively. Besides the intensity derivative , the cosine harmonics basis functions need to be discretized as well. As shown in Appendix B, the set of functions form an orthogonal basis over the discrete space while satisfies the homogeneous Neumann condition on the boundary. Correspondingly, the discretized form of Eq. (28) can be represented asEq. (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 , 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 time (one simple FFT-based fast DCT algorithm is given in Appendix B), the complexity of integral calculation is also of the order , which is a significant improvement over the 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  which uses the finite difference numerical method (using a five-point stencil approximating the Laplacian) to discretize the Poisson equation into a linear system.
So far, the only remaining problem is to calculate the source term 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 is obtained from Eqs. (53) and (55), the vector field can be easily calculated by point-wise multiplication. Finally, the divergence can be obtained throughEqs. (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 from the measured intensity images.
- Transform into cosine frequency domain by DCT (Eq. (43)).
- Reweight the DCT coefficients by (Eq. (27)).
- Transform the modified DCT coefficient back to spatial domain to get (Eq. (44)).
- Calculate .
- Transform into cosine frequency domain by DCT.
- Reweight the DCT coefficients by .
- 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 , 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.
Further, when the intensity is uniform in the defined region, Algorithm 1 can be further simplified as:
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 . 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.
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 . 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  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, is not defined at the zero point of the DCT space, therefore, we multiply by zero instead. Second, since the intensity 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 ). 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 . 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 ), 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 and with the distance of μm. The intensity derivative is estimated by central finite difference
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 where the lateral dimensions and 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, 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 equationFig. 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.
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 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.
To compare the proposed method with previously developed techniques, we applied the FFT-based method developed by Paganin and Nugent  and the FFT-based symmetrization method developed by Volkov et.al  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.
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.
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.
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 μ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%).
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 . The common countermeasures include regularization treatments , using more intensity measurements [21, 22], or simply increasing the defocus distance . By increasing the defocus distance 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 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 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 = 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 in the integral (Eqs. (20) and (23)) is only defined on the aperture region; we should choose as close to 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 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 is defined only on the closure , the gap (intensity zeros) region between 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 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 = 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 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 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.
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  (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 can be decomposed asEq. (14), it can be seen that the rotational term was simply ignored in Teague’s assumption , postulating that must be potential and irrotational. A detailed discussion about the validity and the inaccuracy resulting from “Teague’s assumption” can be found in . Using the identity , 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))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 .
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 defined on , it can be expended into an infinite but convergent series of cosine basis functions42]. The 1-D DCT of is defined as
Similar to the cosine transform of one dimension, for the two-dimensional (2-D) data sequence, their DCT and inverse transform are defined respectively as:
Appendix B: Calculate DCT based on FFT
The direct implementation of DCT involves four nested loops with complexity for an array, which is prohibitive in most applications. As a result, various fast and efficient DCT algorithms have been devised and studied extensively . However, here we only discuss one fast algorithm that utilizes the FFT to compute the DCT and its inverse  due to its simple implementation and relevance to one previous work . 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 asEq. (39), we get
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 mean the lth order derivative of the basis functions of the transform. Then lth order derivative of can be deduced as45]Eq. (51) or Eq. (52) along the direction, once for each :Eq. (27). We further denote the inverse version of Eq. (61) as :
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 and , 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 for an 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 extended from the extending the input N samples according to Eq. (45) with Fourier coefficients 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 as the following rules,Eq. (51)).Eq. (52)).Eqs. (63) and (65) make no difference if one wants to compute the derivative of use the DFT or FFT based method without considering that is an extended version of . Similar relations can be found in 2-D cases: to calculate , , and , one can just extend the input matrix both in and 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
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]
5. N. Streibl, “Phase imaging by the transport equation of intensity,” Opt. Commun. 49(1), 6–10 (1984). [CrossRef]
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]
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]
17. I. W. Han, “New method for estimating wavefront from curvature signal by curve fitting,” Opt. Eng. 34(4), 1232–1237 (1995). [CrossRef]
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]
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]
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).
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.
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.
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).
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]
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.