## Abstract

In this work, an optimum frequency combination (OFC) method is proposed to reconstruct high quality phase information of the complex light field, which is really valuable for many objects such as optical elements and cells. It is shown that the difference image between two symmetrical separated, larger defocused planes contains a lot of lower frequency components of the phase distribution and the higher frequency components can be easily observed in the difference image between two nearly focused planes. Based on the phase transfer function (PTF), our method combines different frequency components with high Signal-to-Noise Ratio (SNR) together to estimate a more accurate frequency spectrum of the object’s phase distribution without any complicated linear or nonlinear regression. Then, we can directly reconstruct a high-quality phase map through inverse Fourier transform. What’s more, in order to compensate the phase discrepancy resulted from strong absorption in the intensity, an iterative compensation algorithm is proposed. Both the simulation and experimental results demonstrate that our iterative OFC (IOFC) method can give a computationally efficient and noise-robust phase reconstruction for absorptive phase objects with higher accuracy and fewer defocus planes.

© 2015 Optical Society of America

## 1. Introduction

Phase information cannot be obtained directly, yet this information is often extremely valuable. Traditional imaging systems only record the intensity of light, however the phase of complex light field contains a lot of important information for many classes of objects such as optical elements, biological soft tissues, and cells [1]. Therefore, quantitative phase imaging (QPI) increasingly plays an important role in many fields such as surface measurement, X-ray imaging [2], electron-beam microscopy [3], neutron radiography [4], and quantitative phase microscopy (QPM) [5–9 ].

Numerous techniques have been developed for the phase retrieval problem, including digital holography [10], phase-shifting interferometry [11–13 ], intensity-based iterative methods [14, 15], transport of intensity equation (TIE) based methods [16, 17], and methods based on Ptychography [18–20 ]. Among these different quantitative phase retrieval methods, the TIE-based method [16, 17] has the advantages of being non-interferometric [21], not needing phase unwrapping [22], applicable with partially coherent beams [23–25 ], and without significant hardware modification and complicated computation [26, 27]. Moreover, by using recently developed multi-plane simultaneously capture schemes instead of simply moving the sample axially [28, 29], TIE method can capture a series of through-focus intensity images in real-time. In the TIE method, the phase of the object is extracted from the first derivative of intensity along the optical axis that occur due to propagation. Therefore, an accurate estimate of intensity derivative should be obtained in order to solve a Poisson equation [30–32 ]. However, because the Poisson inversion process significantly amplifies low spatial frequency information, the TIE has an intrinsic problem of amplifying low-frequency noise, which results in cloudy phase results [33].

Specifically, the cloudy noise problem of TIE phase reconstruction has been found to be dependent on a couple of factors, for example the distances between planes, the number of planes, the absorption of objects, and the SNR level of captured intensity images. Simply increasing the number of planes or utilizing multi-frame de-noise algorithm could improve the reconstruction accuracy certainly, but neither of them is a time-efficient strategy. Therefore, many techniques have been proposed to improve the reconstruction accuracy in TIE. Soto and Acosta [34] proposed a noise-reduction finite difference formula to minimize the noise effect. Afterwards, two well-known high-order derivative estimation methods [35, 36] were respectively introduced by L. Waller and R. Bie to further suppress the cloudy noise in retrieved results. Successively, it was proved that these different derivative estimation methods can be generalized by Savitzky-Golay differentiation filter (SGDF), and then an optimal frequency selection (OFS) method was reported in [37] by analyzing the frequency response of the SGDF and the ideal derivative filter in TIE. In addition, these methods have been extended to the situation in which the intensity is measured in unequally-spaced planes [38, 39]. However, a drawback of these techniques is the lack of analyses for the optimum separation between measurement planes. Therefore, an optimum selection scheme of value and number of the sample-to-detector distances and the photon energy has been carried out to optimize the phase contrast using hard x rays [40]. Lately, two optimum plane selection criteria have been proposed to minimize the noise effect [41, 42]. Each of them introduced an optimum plane selection scheme based on mathematical optimization models and both of them improve the reconstruction quality notably. However, their optimization performance depends on some priori knowledge, such as the noise level and the characteristics of the experimental data. Zhong et al. estimated the derivative of the intensity spectrum by using Gaussian process (GP) regression with an exponential spacing measurement scheme, decreasing the number of required images significantly [43]. Although this algorithm is amenable to parallel processing and can be implemented with Graphics Processing Units (GPUs) for near real-time performance, it is not easy to be accomplished due to the computational complexity of GP regression.

In this paper, we introduce a quantitative phase retrieval algorithm called OFC method for reconstructing high quality phase information. It is shown that the difference image between two symmetrical separated, lager defocused planes contains a lot of lower frequency components of the phase distribution and the higher frequency components can be easily observed in the difference image between two nearly focused planes. Based on this fact, utilizing a combined PTF and an exponential spacing measurement scheme, our method can estimate a more accurate frequency spectrum of the object’s phase distribution, by combining different frequency components with high SNR together without any complicated linear or nonlinear regression. In addition, we introduce an IOFC approach with compensative iteration process to compensate the phase discrepancy resulted from the strong absorption of objects. Both the simulation and experimental results demonstrate that this IOFC method can give a computationally efficient and noise-robust phase reconstruction for absorptive phase objects with higher accuracy and fewer defocus planes.

The remainder of this paper is structured as follows: in Section 2.1 we derive the phase transfer function for TIE and analyze its characteristic. Based on it, in Section 2.2 we then design an OFC phase transfer function, combining all the frequency components of the phase information together without any losing or excessive information. Furthermore, in Section 2.3 we propose the IOFC compensation algorithm to get higher reconstruction accuracy for different objects. Numerical simulation and experimental results are presented in Section 3 and Section 4 respectively. At last, conclusions are summarized in Section 5.

## 2. Principle

#### 2.1. Phase transfer function

Firstly, let us consider the intensity derivative estimation problem in the spatial frequency space. Note this intensity difference is the basis for numerous phase retrieval methods based on TIE or contrast transfer function (CTF) [33, 37, 40, 44]. When the phase *φ*(*x,y*) is small and the intensity image in focus *I*(*x,y*,0) is almost a constant *I*
_{0}, the intensity for a given plane at the z position can be approximate in frequency space by the CTF as [40, 44]

*U*(

*u,v*), and Φ(

*u,v*) are Fourier transforms of $I\left(x,y,z\right),-\frac{1}{2}\mathit{ln}\left(\frac{I\left(x,y,0\right)}{{I}_{0}}\right)$, and

*φ*(

*x,y*), respectively.

*δ*(

*u,v*) denotes Dirac delta function. Therefore, the intensity differences between images at

*z*and −

*z*can be deduced in frequency domain as Eq. (2):

For each defocus distance *z*,
$\widehat{I}\left(u,v,z\right)-\widehat{I}\left(u,v,-z\right)$ varies as a sinusoidal function with frequency (*u*
^{2}+*v*
^{2}). Therefore, to derive an ideal sampling scheme, we extract the phase transfer function *G*(*u,v,z*) from Eq. (2):

For a given position *z*, the phase transfer function is a sinusoidal function of *πλz*(*u*
^{2} + *v*
^{2}) which would be close to zero for particular frequencies, as shown in Figs. 1(a) and 1(b). Figure 1(a) presents the one-dimension PTF curve taken as a function of (*u*
^{2}+*v*
^{2}) and Figure 1(b) shows the 2D representation of |*sin*[*πλz*(*u*
^{2}+*v*
^{2})]|. In Fig. 1(b), black rings denote the range of frequencies where the phase transfer function is close to zero while white rings denote the frequency area where the phase information could be reliably transformed into intensity difference by PTF.

Since in both TIE and our approach, the phase is calculated based on the intensity difference image, their quality is curial to the final phase retrieval result. So, at first, we analyze the influence on the intensity difference image between two defocused planes *I*(*u,v,z*) and *I*(*u,v*,−*z*) in frequency domain which is exerted by the PTF and the noise level by simulating two different situations. For a typical pure-phase object, when the captured images *I*(*u,v,z*) and *I*(*u,v*,−*z*) are free of noise, the frequency spectrum of their difference image
$\widehat{I}\left(u,v,z\right)-\widehat{I}\left(u,v,-z\right)$ is shown in Fig. 1(c). After manually adding white Gaussian noise in *I*(*u,v,z*) and *I*(*u,v,z*), another spectrum can be obtained [Fig. 1(d)]. In Fig. 1(c), blue rings represent the frequency range where objects phase information is lost and the rest purple area denotes the phase components which are transformed into intensity difference. However, in Fig. 1(d), those blue rings are covered by the frequency spectrum resulted from the noise in the captured images *I*(*u,v,z*) and *I*(*u,v*,−*z*). This illustrates that, in an actual experiment, the frequency components within these blue rings have much more information of noise and almost contain no phase information. Furthermore, the distribution of those blue rings is decided by the defocus distance *z* in Eq. (3). So the choice of *z* defines how much phase information from the object is transferred into intensity contrast in the defocused images and how many frequency components around some particular frequencies are lost.

#### 2.2. Exponential spacing measurement scheme

To address this problem in Fourier domain, we note that in order to capture low-frequency phase information, images with large defocus distance are required. Meanwhile, the high-frequency components are transferred in those images with small defocus distance, which are important for recovering accuracy of the result. Therefore, to recover a high quality result, we should capture images both at small defocus distance as well as at large defocus distance.

Here, we utilize an exponential spacing measurement scheme, which has been reported in [43], to build OFC phase transfer function. Aiming to reconstruct phase distribution by combining different frequency components with high SNR, we set a sensitivity threshold *α* to minimize the range of frequencies where the PTF is close to 0. It is known that, the highest frequency
${f}_{1}=\pi \lambda {\left(\frac{NA}{\lambda}\right)}^{2}$ that can be recovered is set by the diffraction limit as *NA*, which is the numerical aperture of the microscope. So, the minimum defocus distance *z*
_{1} is selected as
${z}_{1}=\lambda \frac{\pi -\mathrm{arcsin}(\alpha )}{\pi {\left(NA\right)}^{2}}$ to make sure that *G*(*u,v,z*
_{1}) is *α* at the maximum frequency *f*
_{1}, as presented by the yellow curve intersecting with the purple dot-line in Fig. 2(a). Then, noticing that *G*(*u,v,z*
_{1}) falls below *α* from *f*
_{2} to 0, we select the second defocus distance *z*
_{2} = *βz*
_{1}, where
$\beta =\frac{\pi -\mathrm{arcsin}(\alpha )}{\mathrm{arcsin}(\alpha )}>1$. So that *G*(*u,v,z*
_{2}) is *α* at *f*
_{2}, and will remain at least *α* for a range of frequencies to *f*
_{3}, as illustrated by the green curve in Fig. 2(a).

Similarly, subsequent defocus distances should satisfy the following exponential relation,

Using this exponential spacing measurement scheme, we can build OFC phase transfer function to give a computationally efficient and noise-robust phase reconstruction with higher accuracy and fewer defocus planes. As shown by the red dotted curve in Fig. 2(a), the OFC phase transfer function

which combines those band-pass filtered parts of each phase transfer function*G*(

*u,v,z*) together, satisfies a minimum phase measurement sensitivity

_{k}*α*. Here, $\mathcal{B}\left\{\mathrm{...}\right\}$ denotes the bandpass filter for each phase transfer function

*G*(

*u,v,z*). In other words, the most parts of the proposed OFC phase transfer function are kept above the purple dot-line in Fig. 2(a) and the range of frequencies where OFC falls below

_{k}*α*is minimized. Figure 2(b) shows the 2D representation of the OFC phase transfer function. Comparing with Fig. 1(b), Fig. 2(b) has fewer frequency components where the OFC phase transfer function is close to zero. Similarly, a combined Fourier spectrum of intensity distribution differences could be presented as

Comparing with Fig. 1(c), this combined Fourier spectrum, shown in Fig. 2(c), has no obvious rings where the phase information is missing. So, by using this exponential spacing measurement scheme and the OFC phase transfer function, we can combine all the frequency components with high SNR together and then we can reconstruct the phase distribution with higher accuracy and fewer defocus planes. Mathematically, the Fourier spectrum of phase distribution Φ(*u,v*) could be obtained according to Eq. (2):

At last, we can efficiently reconstruct a high-quality phase distribution *φ*(*x,y*) through inverse Fourier transform.

#### 2.3. Iterative compensation algorithm

We have demonstrated that we can reconstruct the phase distribution of an object by using Eq. (7) when the phase is small and the intensity image in focus is almost a constant. But, as mentioned before, the reconstruction accuracy is affected by the absorption of objects because in a general situation the phase could be large and the intensity image in focus is spatially varying. Figure 3 presents the simulation results of two different situations.

As shown in Figs. 3(a)–3(d), when the intensity image in focus is almost a constant and the phase is small, the reconstructed phase map is very close to the ideal phase map and the difference map between them can be hardly observed. But when the intensity image in focus is spatial varying or the phase is large, the discrepancy can be easily observed by using OFC method, as shown in Figs. 3(e)–3(h). So we propose an iterative compensation algorithm to address this problem.

Firstly, considering a strong-absorptive phase object, its phase distribution can be retrieved by OFC-TIE method. But the retrieved phase map contains a lot of phase discrepancy since the intensity image in focus is spatial varying between [0.1,0.9]. Therefore, we derived a mathematical expression for the discrepancy *R* which need to be compensated in Eq. (14) (see details in Appendix A). Considering the amplitude of
$\mathcal{F}\left\{R\right\}$ is much smaller than
$OF{C}_{\widehat{I}}$, an iterative phase compensation algorithm for OFC-TIE is proposed with a compensation process as descripted below. The block diagram of the entire IOFC-TIE algorithm is summarized in Fig. 4.

Step 1: assuming that the target is a weak pure-phase object. So the Fourier spectrum of phase distribution Φ_{0}(*u,v*) could be obtained easily from a stack of captured images according to Eq. (7). Then, using inverse Fourier transform, the phase distribution *φ*
_{0}, which actually includes phase error, can be obtained.

Step 2: set this incorrect phase map *φ*
_{0} as the initial phase map
${{\phi}^{\prime}}_{j}$ and start iterative compensation algorithm. Here *j* = 0.

Step 3: iterative OFC compensation process invovles four sub-steps from Step 3a to Step 3d.

- Step 3b: figure out an updated phase distribution ${{\phi}^{\prime}}_{j+1}$ according to Eq. (17).
- Step 3c: calculate the difference phase map $\mathrm{\Delta}{{\phi}^{\prime}}_{j}$ between ${{\phi}^{\prime}}_{j}$ (before compensation) and ${{\phi}^{\prime}}_{j+1}$ (after compensation). Here $\mathrm{\Delta}{{\phi}^{\prime}}_{j}={{\phi}^{\prime}}_{j+1}-{{\phi}^{\prime}}_{j}$.
- Step 3d: set a threshold
*ε*to terminate the iteration process. If the difference map $\mathrm{\Delta}{{\phi}^{\prime}}_{j}$ is quite large, send ${{\phi}^{\prime}}_{j+1}$ go back into Step 3a and*j*=*j*+ 1. Otherwise, if $\left|\mathrm{\Delta}{{\phi}^{\prime}}_{j}\right|<\epsilon $, which means the iteration converged, then terminate iteration algorithm and set ${{\phi}^{\prime}}_{j+1}$ as the output.

Step 4: set the output of the iterative OFC compensation process
${{\phi}^{\prime}}_{j+1}$ as the correct phase map of the object *φ*.

Unlike the mixed transfer function approach for TIE (MTF-TIE) based on the mixed contrast transfer model [44] and the iterative compensation algorithm for TIE (IC-TIE) compensating the phase discrepancy by utilizing the Teagues assumption [45], our IOFC-TIE algorithm is proposed based on the OFC phase transfer function. This is the fundamental difference between IOFC-TIE and other iterative phase compensation algorithms. In order to evaluate the efficiencies of IC-TIE and IOFC-TIE, both of them will be tested afterwards.

## 3. Simulations

#### 3.1. Pure phase object

To confirm and characterize the performance of various phase recovery methods that use equal spacing and exponential spacing schemes, we simulated the propagation of a complex field of a pure phase profile of the letters ‘TIE’, shown in Fig. 6(a), which was defined on a grid with 256 × 256 pixels (pixel size of 2*µm* × 2*µm*). In the simulation, the illumination wavelength is set as 632.8nm. Equally Spaced Stack 1 has 129 intensity images, simulated with a constant defocus step size of 2*µm* over a large defocus range [−128*µm* to 128*µm*]. Then, we extract some images from Equally Spaced Stack 1 to build up other three stacks. Equally Spaced Stack 2 has 17 intensity images, simulated with a constant defocus step size of 2*µm* over a small defocus range [−16*µm* to 16*µm*]. Equally Spaced Stack 3 also has 17 intensity images, but simulated with a larger constant defocus step size of 16*µm* over a large defocus range [−128*µm* to 128*µm*]. Exponentially Spaced Stack has 15 images and they are exponentially spaced, at z positions around focus of *±*2*µm*, *±*4*µm*, ±8*µm*, ±16*µm*, *±*32*µm*, *±*64*µm*, and ±128*µm*. In order to evaluate the accuracy all the simulated intensity images are corrupted by additive white Gaussian noise with standard deviation ranging from 0 to 0.2. The Poisson noise case has been also tested on the same data. The results did not differ significantly from the case of Gaussian and hence not discussed further in this work.

Figure 5 shows the root mean square error (RMSE) of the recovered phase by using three different phase retrieval algorithms as the standard deviation of noise increases. OFS-TIE method is employed with the equally spaced stacks while GP-TIE and OFC-TIE methods are employed with the exponentially spaced stack. For the equally spaced stacks, by analyzing the ascending tendency of each curve in Fig. 5, it can be concluded that small defocusing distance performs obviously worse than large defocusing distance with same number of sampling planes, and the stack which has 129 measurement planes performs best. This can be explained by the fact that the Equally Spaced Stack 2 contains more high-frequency phase content in the measurements while the Equally Spaced Stack 3 contains more low-frequency phase content. So the phase recovered from Equally Spaced Stack 2 is contaminated by low-frequency noise while the phase recovered from Equally Spaced Stack 3 is perceptibly blurred. Since the Equally Spaced Stack 1 has 129 images, which containing all the excessive frequency components of the object, the phase recovered from it clearly exhibits the lowest RMSE expectedly as the standard deviation of noise increases. For the Exponentially Spaced Stack, it only requires 15 images and it almost performs as good as the equally spaced stack which has 129 measurement planes. This is because that the exponentially spaced stack contains both the low-frequency components and the high-frequency components of the phase information, and there is no trade-off between noise and nonlinearity. In addition, comparing with GP-TIE, OFC-TIE achieves better performance as shown in Fig. 5. Figure 6 gives an example of the recovered phase at standard deviation of 0.01. It seems that the phase retrieved from the Equally Spaced Stack 1 and Equally Spaced Stack 3 appears almost same. This is due to the fact that low-frequency information is generally sufficient to present the outline of a target. In order to highlight the details under cover, the enlargement part of Fig. 6(e) is presented and the blurring edge is quite evident comparing with the accurately retrieved phase map in enlargement part of Fig. 6(c). What’s more, it only takes OFC-TIE about 0.028s to finish the phase reconstruction with 256 × 256 pixels. But it takes GP-TIE 0.32s, about 11 times longer than OFC-TIE, to reconstruct the phase distribution in the same situation.

#### 3.2. Complex phase and intensity distribution

In the previous simulation, the complex field was constructed with a weak pure-phase object. Next, to test the performance of IOFC-TIE algorithm for absorptive phase objects, we consider the phase profile of letters ‘TIE’ in Fig. 7(a) and the complex intensity distribution of a ‘digital eye’ logo in Fig. 7(b). The value range of the intensity is [0.1, 0.9] and the noise standard deviation is 0.04. Figures 7(c), 7(d), and 7(e) show the phase distribution of the OFS-TIE method with different equally spaced stacks. Since the intensity distribution in-focus varying spatially between [0.1, 0.9], the phase error mainly emerges from where the intensity is close to zero, as shown in Fig. 7(c). Figures 7(f) and 7(g) show the phase distribution of the GP-TIE, OFC-TIE methods before compensation. The reconstruction quality of the result in Fig. 7(g) is little better than the result in Fig. 7(f). But both of them contain a lot of low-frequency noise and amount of phase discrepancy which are resulted from strong absorption. Therefore, two iterative phase compensation algorithms are tested.

Figure 8(a) presents the comparison of IC-TIE algorithm and IOFC-TIE algorithm for their compensation efficiencies. Utilizing IC-TIE method, after five iterations of compensation, the majority amount of phase discrepancy, which was caused by the absorption, was removed as shown in Fig. 8(c). But the iterative IC-TIE algorithm converge to a result which still contains amount of obvious phase discrepancy where the intensity is very small. However, utilizing IOFC-TIE compensation algorithm, after five iterations of compensation, the phase discrepancy was removed more radically, as shown in Fig. 8(e), and the compensated phase distribution almost equals to the true phase distribution in Fig. 7(a).

These results demonstrate that, comparing with GP-TIE algorithm, OFC-TIE algorithm could achieve better result before compensation. Moreover, comparing with IC-TIE algorithm, IOFC-TIE compensation algorithm could significantly reduce the phase error and converging to an accurate phase result with no more than five iteration times, which demonstrating its high efficiency.

## 4. Experiments

The simulation results have presented that IOFC requires fewer sampling planes to obtain a high quality phase result, with exponential spacing measurement scheme. To demonstrate it experimentally, we compare the recovered phase map by using OFS-TIE, GP-TIE and IOFC-TIE with different sampling strategies. First we test these methods under coherent illumination. The experimental test setup is presented in Fig. 9(a). The test source is a beam from a He-Ne laser (wavelength 632.8nm) that has been expanded and collimated, then passed through the object under test. The object was imaged onto a monochrome CCD camera with a standard 4 *f* system - two lenses of focal length *f* = 25*mm* separated by the distance 2*f*, and the distance from the object to the first lens is *f*. The camera is set on a translation stage in order to modify the virtual distance Δ*z* between the image of the surface and the camera. The phase object under test is a geometry pattern etched on polymethyl methacrylate (PMMA) substrate, which is also shown in Fig. 9(a). Figure 9(b) illustrates the tested phase sample captured in a DIC Microscope. The DIC microscope rendered a pseudo three-dimensional relief shading intensity which related to the phase gradient of the object. It is evidently shown that the test object includes phase ramps of different heights, and more importantly, it contains some regular periodic vertical line structures. It can be estimated that the period of the line structures is 25.9*µm*. Therefore, the corresponding spatial frequency is 38.61*cycles/mm*, which was sufficient to be resolved by the 4*f* system.

We shall be interested in the behavior and limitations of the phase retrieval algorithm across a range of defocus distance. For this purpose, we acquired five sets of images at defocus distances 40*µm*, 120*µm*, 600*µm*, 1040*µm*, and 1280*µm*, respectively. Their corresponding difference images are shown in Fig. 10. The left column shows three images for each defocus distance. Left two are the captured defocusing images and each of them are either less defocused (left up) or further defocused (left down) on top of the defocus distances mentioned, while the image on the right shows the difference between left two images. The contrast ratio of the right image represents the amount of phase information which has been transformed into the intensity difference with light propagating. Since in both TIE and our approach, the phase is calculated based on the intensity difference image, their quality is curial to the final phase retrieval result. The corresponding theoretical PTF and TIE PTFs are shown in the right column. It can be seen that when the defocus value was small (40*µm* or 120*µm*), the phase contrast effect was not evident. This is because the PTF in most of the object frequency range is near zero, as shown in Figs. 10(a) and 10(b). In this case, the intensity difference image was barely usable for phase retrieval. With the defocus value increased, the quality of the intensity difference improved gradually. When the defocus distance reached 600*µm*, the contrast of the period structure was the best, which is in accordance with the corresponding PTF curves plotted in Fig. 10(f). When the defocus value continued increasing to about 1040*µm*, the regular vertical line structures almost disappeared. This is because the spatial frequency of the structure fall into the vicinity of zero crosses of PTF in Fig. 10(h). When the defocus value increase again to about 1280*µm*, the structure reappeared and showing inverse contrast accordingly. Besides, the rapid changing PTF curve in high frequency part results in the evident ripple effect around the edge of the large phase changes. All these results confirmed that the estimated model presented in Section 2 is sufficiently accurate.

Then, in order to reconstruct the phase images, a Data Set containing 129 images was captured, equally spaced by a constant small step size *dz* = 20*µm* over a large defocus range [−1280*µm* to 1280*µm*]. Some data samples from the intensity stack are shown in Fig. 11(a). The intensity distribution in-focus, varying between [0.41, 0.50], illustrating the phase object under test is a weak-absorptive object. First, the best phase result in Fig. 11(b) is recovered by using OFS-TIE method with Equally Spaced Stack 1, which involves all 129 images. Then, we extract some images from Equally Spaced Stack 1 to make up other three stacks. Equally Spaced Stack 2 has 17 intensity images, simulated with a constant defocus step size of 20*µm* over a small defocus range [−160*µm* to 160*µm*]. Equally Spaced Stack 3 also has 17 intensity images, but simulated with a larger constant defocus step size of 160*µm* over a large defocus range [−1280*µm* to 1280*µm*]. Exponentially Spaced Stack also has 15 images and they are exponentially spaced, at z positions around focus of *±*20*µm*, *±*40*µm*, *±*80*µm*, *±*160*µm*, *±*320*µm*, *±*640*µm*, and *±*1280*µm*. Using Equally Spaced Stack 2, if we recover phase with OFS-TIE method, the result becomes susceptible to low-frequency noise, as shown in Fig. 11(c). This is due to the fact that the low-frequency information of phase is not well captured at small defocus distances. If we use Equally Spaced Stack 3, increasing the step size, high-frequency components are lost and low-frequency noise still exists because of the intensity varying in-focus, shown in Fig. 11(d). In order to accurately capture both high and low-frequency information with the same reduced number of images, we utilized exponentially spaced measurements. Unfortunately, as can be seen from Fig. 11(e), the phase result still involves low-frequency noise by using GP-TIE. Because only if the captured images are free of noise, the excessive 112 images in Equally Spaced Stack 1 is completely redundant comparing with the 17 images involved in Exponentially Spaced Stack. But when the captured images are suffer from heavy noise, those excessive images in Equally Spaced Stack 1 will restrain the reconstructed cloudy noise in Fig. 11(b). Same conclusion can be drawn from Fig. 5. Except for the situation when the noise standard deviation is 0, OFS-TIE method with Equally Spaced Stack 1 performs much better than the GP-TIE method with Exponentially Spaced Stack as noise level increasing. However, comparing with the phase map recovered by using GP-TIE algorithm shown in Fig. 11(e), the result recovered from OFC-TIE method in Fig. 11(f) involves less cloudy phase error since OFC-TIE method only use the frequency components with high SNR to retrieve the phase information. Moreover, it should be emphasized that OFS-TIE method performs better than GP-TIE and OFC-TIE methods only for the pure-phase objects or weak-absorptive objects with an redundant image dataset such as Equally Spaced Stack 1. With fewer defocus planes, OFC-TIE gives a more reliable result.

To better assess the accuracy of the phase measurements, we use a confocal microscope with a 50× objective (NA=0.8) to independently measure the surface profile of a small area corresponding to the magnified area shown in Fig. 11(b). Fig. 12(a) shows the inverse height distribution and one line profile was extracted for better quantitatively evaluation in Fig. 12(b). The same line from Figs. 11(b)–11(f) were also extracted, plotted in Fig. 12(c), converted to inverse height (refractive index *n _{PMMA}* = 1.49). The offsets of the five curves were adjusted to the same level. The result obtained by confocal microscope gave 136nm for the height of the large step and 42.1nm for the average peak-to-valley height of the small periodic structure. The OFS-TIE method produced a step height of 138nm and underestimated the height of the small variations 41.2nm when using Equally Spaced Stack 1. For OFS-TIE method using Equally Spaced Stack 2, the step height extremely reduced from 160nm to about 80nm due to the low-frequency noise. When using OFS-TIE method with Equally Spaced Stack 3, the step height varies stably between 136nm and 152nm but the average peak-to-valley height of the small periodic structure reduced to 15.8nm, which means a lot of high-frequency information is missing. Comparing with GP-TIE, OFC-TIE method provided a more faithful result with 139nm step height and the height of the small variations 39.2nm with Exponentially Spaced Stack. Thus, we have reduced the data capture requirement and avoided the negative impacts resulting from intensity absorption, without sacrificing reconstruction quality.

In addition, we also tested our approach for measuring a sample of unstained *SMMC*−7721 human hepatocellular carcinoma cells, using bright-field microscope (magnification 20*×*, NA = 0.45) with filtered white light illumination (center wavelength 525nm). Then, a Data Set containing 129 images was captured, equally spaced by a constant small step size *dz* = 1*µm* over a defocus range [−64*µm* to 64*µm*]. Similarly, three equally spaced stacks and one exponentially spaced stack are extracted from the data set. The reconstructed results are presented in Fig. 13 and again they demonstrate that IOFC-TIE method can achieve high-quality phase reconstruction somewhat better than other methods. Note that three cloudy phase maps are recovered by OFS-TIE as shown in Figs. 13(b)–13(d). This is due to the fact that the intensity distribution in-focus varies enormously between [0.12, 0.97]. In addition, although the recovered phase map is free of low-frequency noise using GP-TIE, the details of the cells are blurred as shown in Fig. 13(e). However, comparing with GP-TIE, IOFC-TIE method can recover a high-quality phase distribution with more detailed information and the result is also free of low-frequency noise, as shown in Fig. 13(f).

## 5. Conclusion

This paper has demonstrated both theoretically and experimentally that a high-quality, and noise-robust phase reconstruction can be obtained efficiently by an IOFC-TIE method based on the PTF with fewer defocus planes. Different from those equally spacing measurement strategy, our proposed IOFC-TIE method utilizing an exponential spacing measurement scheme, decreasing the number of required images significantly. Comparing with the developed GP-TIE algorithm, which involves a complex GP regression procedure, OFC-TIE method is simpler and much better not only for its accuracy, but also for its efficiency. Furthermore, even for an absorptive phase object, its phase distribution still can be obtained using IOFC-TIE method in a high quality.

Although the IOFC-TIE improves the noise resistibility and the phase resolution of TIE significantly, its performance also limited by the absorption of the object. When the minimum value of the intensity in-focus is lower than 0.05, our iterative compensation algorithm may not converge to the correct phase distribution. Moreover, the conventional TIE method has an advantage of being non-interferometric as mentioned before. However, when the illumination source loses coherence, reconstruction error will occur in the retrieved phase image by employing IOFC-TIE algorithm since the PTF no longer varies as an ideal sinusoidal function with frequency (*u*
^{2} +*v*
^{2}) for partially coherent illumination [23–25
]. Therefore, the proposed IOFC-TIE algorithm could achieve its best performance only if a coherent illumination is adopted. This seems to be another limitation of our approach and it will be the subject of future work.

## Appendix A: derivation of OFC-TIE

At the beginning, assume that the complex amplitude of an object in focus is

*w*

_{0}is the complex amplitude of an object,

*A*and

*φ*present the amplitude and the phase of the object. Here we defined

*A*=

^{0}*Acos*(

*φ*) and

*φ′*=

*tan*(

*φ*). So the Fourier spectrum of

*w*

_{0}is

*z*is

*z*,

*λ*,

*u*, and

*v*present the position of the defocus plane, the wavelength and the frequency coordinates respectively. Here we define $\omega =kz(1-\sqrt{1-{\lambda}^{2}({u}^{2}+{v}^{2})})$. Next, the intensity at the defocus plane can be presented as

*z*and −

*z*

Equation (12) seems quite different from TIE. But when employing weak phase approximation, paraxial approximation and nearly focused approximation, TIE can be derived from Eq. (12) easily.

In order to point out difference between Eq. (12) and Eq. (2), we assume that amplitude *A′* is made up by two components, the average part
${{A}^{\prime}}_{m}$ and the residual part
${{A}^{\prime}}_{r}$. The average part is a constant and the residual part varies spatially. Then we can get

*R*denotes the phase discrepancy which need to be compensated. So Eq. (13) can be derived into

Then, from the OFC phase transfer function Eq. (5), Eq. (15) can be reduced to

At last, we can reconstruct the correct phase distribution *φ*(*x,y*) through inverse Fourier transform and inverse tangent function

## Acknowledgments

This work was supported by the National Natural Science Fund of China ( 11574152, 61505081), ‘ Six Talent Peaks’ project ( 2015-DZXX-009, Jiangsu Province, China) and ‘ 333 Engineering’ research project ( BRA2015294, Jiangsu Province, China), Fundamental Research Funds for the Central Universities ( 30915011318), and Open Research Fund of Jiangsu Key Laboratory of Spectral Imaging & Intelligent Sense ( 3092014012200417). C. Zuo thanks the support of the ‘Zijin Star’ program of Nanjing University of Science and Technology. The authors would also like to thank Jingshan Zhong for the open source code of GP-TIE.

## References and links

**1. **G. Popescu, *Quantitative Phase Imaging of Cells and Tissues* (McGraw Hill Professional, 2011).

**2. **K. A. Nugent, “Coherent methods in the X-ray sciences,” Adv. Phys. **59**, 1–99 (2010). [CrossRef]

**3. **S. Bajt, A. Barty, K. A. Nugent, M. McCartney, M. Wall, and D. Paganin, “Quantitative phase-sensitive imaging in a transmission electron microscope,” Ultramicroscopy **83**, 67–73 (2000). [CrossRef] [PubMed]

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

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

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

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

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

**9. **C. Zuo, J. Sun, J. Zhang, Y. Hu, and Q. Chen, “Lensless phase microscopy and diffraction tomography with multi-angle and multi-wavelength illuminations using a LED matrix,” Opt. Express **23**, 14314–14328 (2015). [CrossRef] [PubMed]

**10. **P. Marquet, B. Rappaz, P. J. Magistretti, E. Cuche, Y. Emery, T. Colomb, and C. Depeursinge, “Digital holographic microscopy: a noninvasive contrastimaging technique allowing quantitative visualization of living cells with subwavelength axial accuracy,” Opt. Lett. **30**, 468–470 (2005). [CrossRef] [PubMed]

**11. **R. K. Crane, “Interference phase measurement,” Appl. Opt. **8**, 538 (1969).

**12. **O. Y. Kwon, “Multichannel phase-shifted interferometer,” Opt. Lett. **9**, 59–61 (1984). [CrossRef] [PubMed]

**13. **T. Weitkamp, A. Diaz, C. David, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “X-ray phase imaging with a grating interferometer,” Opt. Express **13**, 6296–6304 (2005). [CrossRef] [PubMed]

**14. **R. Gerchberg and W. Saxton, “Phase determination for image and diffraction plane pictures in the electron microscope,” Optik **34**, 275–284 (1971).

**15. **J. R. Fienup, “Reconstruction of an object from the modulus of its Fourier transform,” Opt. Lett. **3**, 27–29 (1978). [CrossRef] [PubMed]

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

**17. **L. J. Allen and M. P. Oxley, “Phase retrieval from series of images obtained by defocus variation,” Opt. Commun. **199**, 65–75 (2001). [CrossRef]

**18. **J. M. Rodenburg and H. M. Faulkner, “A phase retrieval algorithm for shifting illumination,” Appl. Phys. Lett. **85**, 4795–4797 (2004). [CrossRef]

**19. **P. Thibault, M. Dierolf, O. Bunk, A. Menzel, and F. Pfeiffer, “Probe retrieval in ptychographic coherent diffractive imaging,” Ultramicroscopy **109**, 338–343 (2009). [CrossRef] [PubMed]

**20. **G. Zheng, R. Horstmeyer, and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nat. Photonics **7**, 739–745 (2013). [CrossRef]

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

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

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

**24. **A. M. Zysk, R. W. Schoonover, P. S. Carney, and M. A. Anastasio, “Transport of intensity and spectrum for partially coherent fields,” Opt. Lett. **35**, 2239–2241 (2010). [CrossRef] [PubMed]

**25. **M. H. Jenkins, J. M. Long, and T. K. Gaylord, “Multifilter phase imaging with partially coherent light,” Appl. Opt. **53**, D29–D39 (2014). [CrossRef] [PubMed]

**26. **S. S. Gorthi and E. Schonbrun, “Phase imaging flow cytometry using a focus-stack collecting microscope,” Opt. Lett. **37**, 707–709 (2012). [CrossRef] [PubMed]

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

**28. **P. M. Blanchard and A. H. Greenaway, “Simultaneous multiplane imaging with a distorted diffraction grating,” Appl. Opt. **38**, 6692–6699 (1999). [CrossRef]

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

**30. **M. Beleggia, M. A. Schofield, V. V. Volkov, and Y. Zhu, “On the transport of intensity technique for phase retrieval,” Ultramicroscopy **102**, 37–49 (2004). [CrossRef] [PubMed]

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

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

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

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

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

**36. **R. Bie, X. H. Yuan, M. Zhao, and L. Zhang, “Method for estimating the axial intensity derivative in the TIE with higher order intensity derivatives and noise suppression,” Opt. Express **20**, 8186–8191 (2012). [CrossRef] [PubMed]

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

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

**39. **S. Zheng, B. Xue, W. Xue, X. Bai, and F. Zhou, “Transport of intensity phase imaging from multiple noisy intensities measured in unequally-spaced planes,” Opt. Express **20**, 972–985 (2012). [CrossRef] [PubMed]

**40. **S. Zabler, P. Cloetens, J. P. Guigay, J. Baruchel, and M. Schlenker, “Optimization of phase contrast imaging using hard x rays,” Rev. Sci. Instrum. **76**, 073705 (2005). [CrossRef]

**41. **K. Falaggis, T. Kozacki, and M. Kujawinska, “Optimum plane selection criteria for single-beam phase retrieval techniques based on the contrast transfer function,” Opt. Lett. **39**, 30–33 (2014). [CrossRef]

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

**43. **J. Zhong, R. A. Claus, J. Dauwels, L. Tian, and L. Waller, “Transport of Intensity phase imaging by intensity spectrum fitting of exponentially spaced defocus planes,” Opt. Express **22**, 10661–10674 (2014). [CrossRef]

**44. **J. P. Guigay, M. Langer, R. Boistel, and P. Cloetens, “Mixed transfer function and transport of intensity approach for phase retrieval in the Fresnel region,” Opt. Lett. **32**, 1617–1619 (2007). [CrossRef] [PubMed]

**45. **C. Zuo, Q. Chen, L. Huang, and A. Asundi, “Phase discrepancy analysis and compensation for fast Fourier transform based solution of the transport of intensity equation,” Opt. Express **22**, 17172–17186 (2014). [CrossRef] [PubMed]