Transport of Intensity phase imaging by intensity spectrum fitting of exponentially spaced defocus planes.

: We propose an alternative method for solving the Transport of Intensity equation (TIE) from a stack of through–focus intensity images taken by a microscope or lensless imager. Our method enables quantitative phase and amplitude imaging with improved accuracy and reduced data capture, while also being computationally efﬁcient and robust to noise. We use prior knowledge of how intensity varies with propagation in the spatial frequency domain in order to constrain a ﬁtting algorithm [Gaussian process (GP) regression] for estimating the axial intensity derivative. Solving the problem in the frequency domain inspires an efﬁcient measurement scheme which captures images at exponentially spaced focal steps, signiﬁcantly reducing the number of images required. Low–frequency artifacts that plague traditional TIE methods can be suppressed without an excessive number of captured images. We validate our technique experimentally by recovering the phase of human cheek cells in a brightﬁeld microscope.


Introduction
Quantitative phase imaging has found useful applications in biology, surface metrology and X-ray imaging [1,2]. Methods that use a series of through-focus intensity images (e.g., [3][4][5][6][7][8][9]) are especially popular due to their experimental simplicity. In-focus intensity images contain no phase information; however, defocus introduces phase contrast. In fact, any imaging system with a complex transfer function will provide some phase contrast. These images can then be inverted to recover phase and amplitude quantitatively. In defocus-based methods, the experimental procedure involves simply moving the sample (or camera) axially while capturing multiple images through-focus (see Fig. 1), or using any of the recently proposed schemes for simultaneous multi-plane capture [10][11][12][13].
Recovering phase (and amplitude) from a series of defocused images is a nonlinear problem because intensity is bilinear with phase and amplitude. One of the most successful methods for this inversion is the iterative method, which takes a nonlinear convex optimization approach [14, z stage Experimetal setup: Nikon TE300 Through-focus intensity stack z camera Fig. 1. Generalized experimental setup for an imaging system (e.g. a microscope) which captures intensity images at a range of axial defocus distances in order to recover phase. 15]. Though the problem is non-convex, results usually converge, especially when multiple images are used at various focal planes [5]. In another approach, where a direct solution is desired, the problem is linearized for small defocus distances using the Transport of Intensity Equation (TIE) [3,4]. The TIE relates phase and amplitude to variations of intensity along the optical axis z [3]: where I(x, y, 0) is the intensity at focus, φ (x, y) is the phase, λ is the spectrally-weighted mean wavelength of illumination and ∇ ⊥ denotes the gradient operator in lateral directions (x, y) only. Besides being simple, the TIE method achieves phase accuracy equivalent to that of interferometric methods, but it is much less sensitive to coherence in the illumination [16]. Numerical solvers require no phase unwrapping and are computationally efficient, employing FFT-based inversion solvers [17]. Numerical solutions for the TIE can be understood by considering the case of a pure-phase object, where I(x, y, 0) is constant [18]. In this situation, the right hand side of Eq. (1) simplifies to a second derivative (Laplacian). Phase can then be recovered from a Laplacian inversion where Φ(u, v) is the Fourier transform of φ (x, y), F(u, v) is the Fourier transform of the measured first derivative ∂ I(x,y,z) ∂ z z=0 scaled by 2π/λ , and u, v are the spatial frequency variables.
When I(x, y, 0) is not constant, two Laplacian inversions are required [3]. Note that the denominator of Eq. (2) goes towards zero as the spatial frequency goes towards zero, and a small regularization parameter must be added in order to avoid division-by-zero instability. This means that the DC phase term is lost and low frequency noise is amplified, resulting in phase errors that give cloudy phase results. We discuss here strategies for mitigating these problems. Limitations of the TIE method mainly stem from noise and nonlinearity in the intensity derivative estimate, both of which we address here. Nonlinearity error comes from the estimation of the intensity derivative, ∂ I(x,y,z) (fitting to noise) and the optimal order to use will be phase dependent, making it a difficult parameter to optimize. The Savitzky-Golay differentiation filter (SGDF) TIE [20] was proposed to solve this trade-off between the order of fitting and noise in higher order TIE. It recovers phase images with different orders of polynomial fitting, then combines these phase images into a final reconstructed phase with band-pass filters. Limitations are that the fitting becomes unstable when the order is large [21] and the derivation assumes equally spaced focus steps, which we show here to be non-ideal.
We propose here a new TIE phase recovery method in which we also perform throughfocus fitting of the data, except we do the fitting in Fourier space, rather than in real space. The advantage of this method can be understood from Fig. 2. On the left, we show a stack of through-focus intensity images, with a few example plots of the intensity vs. z behavior. These are the curves that would be fit to polynomials in higher order TIE [7]. It is easy to see that they do not conform to any particular functional form. On the right side of Fig. 2, we take 2D Fourier transforms of the intensity images at each of the defocus steps (showing here only the real part), then plot the intensity spectrum vs. z plots for a few example spatial frequencies. In contrast to the real space plots, these curves display a distinct sinusoidal behavior which can be predicted by the Contrast Transfer Function [22][23][24]. By using this a priori information of how light propagates in spatial frequency space, we can simultaneously deal with nonlinearity and noise in order to obtain a better derivative estimate. Here, we propose to fit the spectral data with Gaussian process (GP) regression [25]. GP regression does not require that the inputs are equally spaced and can suppress the unwanted high-frequency components in the fitted function by using the squared exponential covariance function (see details in Appendix A). The curves in Fig. 2 follow sinusoidal trajectories, but actual behavior may deviate from this model due to breaking of the small phase approximation or unmodeled coherence effects. GP regression is able to flexibly fit the data to variations in the model using the Contrast Transfer Function as a priori information to set a threshold for suppressing the unwanted high frequency components in the fitted function. In this way, the over-fitting to noise is avoided without knowing or fitting to the precise function of a sinusoid.

Intensity
Intensity spectrum  In studying the frequency behavior for through-focus images, we find that equally spaced defocus steps in z are not ideal. Previous methods generally use images equally spaced in z, though higher order TIE can accommodate arbitrary z steps [26]. In a recent study [27], an exponentially-growing measurement scheme is proposed to deal with the regularization problem. Here, we show that exponentially spaced defocus steps can provide a large reduction in the number of images required for a high quality phase result. We suggest an efficient scheme for choosing the z plane distances (within the limits of the motion stage), from the perspective of efficiently transferring phase information to intensity images, and show how our method can achieve high quality phase results with far fewer images than equally spaced schemes.
The rest of the paper is structured as follows. In Section 2, we apply GP regression to TIE phase recovery and propose an exponential spacing measurement scheme. In Section 3, we compare the error performance of the proposed algorithm with traditional TIE methods. In Section 4, we discuss the experimental results. We offer concluding remarks in Section 5.

Intensity spectrum fitting
We consider the intensity derivative estimation problem from the spatial frequency domain. When the phase φ (x, y) is small, the variations of intensity over z can be approximated in frequency space by the Contrast Transfer Function [22-24]: where δ (u, v) denotes Dirac delta function, and I (u, v, z), U(u, v), and Φ(u, v) are Fourier transforms of I(x, y, z), − 1 2 ln I(x, y, 0), and φ (x, y), respectively. For each spatial frequency (u m , v n ), I (u m , v n , z) follows a sinusoidal pattern that oscillates with frequency πλ (u 2 m + v 2 n ). Thus, instead of fitting the intensity curves in real space to polynomials (as in higher order TIE [7]), we use the sinusoids in the intensity spectra through-focus as prior in the fitting.

TIE using Gaussian process regression
To fit the data to our model, we employ Gaussian process regression. We estimate the derivative of the intensity spectrum with respect to defocus ∂ I (u,v,z) which are the 2D Fourier transforms of the measured intensity images. For each spatial fre- can be obtained from the continuous function fitted from the discrete data points I (u m , v n , z 1 ), ..., I (u m , v n , z N ). GP regression is a suitable fitting method that allows us to fit the discrete points I (u m , v n , z 1 ), ..., I (u m , v n , z N ) to a continuous sinusoidal function (see details in Appendix A). In the fitted function, frequency components of the axial oscillation that are higher than a threshold s c can be suppressed by initializing appropriate hyper-parameters in the regression (see details in Appendix B). From the Contrast Transfer Function (Eq. (3)), I (u m , v n , z), as a function of z, does not contain frequency components higher than πλ (u 2 m + v 2 n ). Therefore, the threshold can be set just above πλ (u 2 m + v 2 n ) in order to eliminate high-frequency noise while fitting over the inputs I (u m , v n , z 1 ), ..., I (u m , v n , z N ). Although the Contrast Transfer Function is derived under the weak phase assumption, it only serves as a guideline to pick a frequency threshold, above which variations are treated as noise. Thus, large phase variations can still be recovered, beyond the small phase approximation. Furthermore, by using GP regression, other prior knowledge about the sinusoidal pattern could be incorporated into the regressionfor example, coherence effects of the illumination could be modeled in the Contrast Transfer For each frequency v m , v n : Initialize the hyper-parameters in GP regression, in order to suppress the frequency components higher than s c . (4) Output is obtained by calculating through Eq. (17).
Function as a product term, since reduced spatial coherence blurs defocused images and thus reduces the amplitude of the oscillations away from the focus plane. The resulting algorithm of GP TIE is summarized in Table 1. The scaling factor γ in Step (2) of the table allows a trade-off between accuracy and noise filtering. Since the computational complexity of GP regression is proportional to the cubic of the number of measurements [25], the algorithm has a computational complexity of O(N(N z ) 3 ), where N is the number of pixels in one image, and N z is the number of measurements. This algorithm is very amenable to parallel processing and can be implemented with Graphics Processing Units (GPUs) for near real-time performance, in conjunction with GPU-based TIE solvers [28]. Further, GP regression does not require the measuring positions z 1 , z 2 , ..., z N to be equally spaced, so any set of measurement positions along the axial dimension can be used as input. We show in the next session why this is not only convenient, but also leads to better phase recovery.

Exponential spacing measurement scheme
In looking at this problem from Fourier space, we see that equally spaced z steps are not the ideal measurement scheme. The choice of z distances is crucial, since it defines how much phase information from the object is transferred into intensity contrast in the defocused images [29]. Each z distance measurement can be thought of as having its own phase to intensity transfer function, with its own set of preferred spatial frequencies. The low-frequency information of phase is particularly poorly represented, as seen in the Laplacian inversion of Eq. (2). To better capture this low frequency phase information, images with large z are required, ideally out to the maximum range of the axial motion stage. However, the high-frequency components favor a small z and are important for recovering fine detail in the phase result. We would like to capture images both at the minimum z possible as well as at the maximum z possible. For linearly spaced schemes, this will require excessive data capture in between. To avoid this, we introduce an unequally spaced measurement scheme that balances these concerns. Since diffraction goes as (∆x) 2 /λ z [18], with ∆x relating to the object size, we expect the ideal spacing to follow a nonlinear trajectory.
To derive an ideal z sampling scheme, we start with Eq. (3) and extract the phase transfer function, g(u, v, z): This equation tells us how much phase information for a given spatial frequency is converted into intensity contrast, and thus relates to the signal-to-noise ratio (SNR) of the measurement. For a given position z, the phase transfer function is a sine function of πλ (u 2 + v 2 ) that will be maximized for particular spatial frequencies, as shown in Fig. 3. Since each object contains a different mix of spatial frequencies in its phase information, the theoretically optimal measurement planes are object-dependent [29], so it is not possible to determine them without knowing the object itself. Thus, we aim not for the absolute optimal set of measurements, but rather for the optimal coverage of sensitivity across all spatial frequencies that pass through the imaging system, treating each equally. Consider the goal of selecting a set of z planes such that the maximum phase measurement sensitivity (set by the value of the transfer function) is at least α for the largest range of spatial frequencies possible. In the following, we demonstrate an exponential spacing measurement scheme that achieves this goal with the least number of measurement required. The highest frequency that can be recovered is set by the diffraction limit as NA/λ , where NA is the numerical aperture of the imaging system. This will define the minimum defocus distance from which we can capture relevant information. We define f = πλ (u 2 +v 2 ) and f ≤ πλ ( NA λ ) 2 . First, we select the minimum defocus distance z 1 such that the sensitivity g(u, v, z 1 ) is α at frequencies corresponding to the maximum frequency allowed f 1 = πλ ( NA λ ) 2 , the solution of which is z 1 = π−arcsin(α) π(NA) 2 λ . Defocus steps smaller than z 1 will not provide useful information and are thus unnecessary, though in practice the minimum z step size may be set by the axial motion stage. As shown by the blue curve in Fig. 3, the sensitivity remains above α until f 2 = arcsin(α)π(NA) 2 λ [π−arcsin(α)] . Next, we select the second defocus distance z 2 = β z 1 , where β = π−arcsin(α) arcsin(α) > 1, such that the sensitivity g(u, v, z 2 ) is α at f 2 , and will remain at least α for a range of frequencies, as illustrated by the green curve in Fig. 3. By induction, the optimal measurement scheme that satisfies a minimum phase measurement sensitivity α should satisfy the following exponential relation for the defocus distances The exponential spacing implies that a large z can be reached with far fewer measurements as compared to the equal-spacing measurement schemes. Large z images are crucial for transferring low-frequency phase information, and the exponential spacing scheme enables us to reach this without taking an excessive number of measurements or trading off high-frequency information. The larger the maximum z, the better the low-frequency result, and we desire to keep the minimum z as small as possible in order to maintain diffraction-limited phase resolution, both extremes being limited by the motion stage axial range. Accuracy can be traded off against number of images through the choice of α. Larger values results in more accurate phase retrieval, at the cost of requiring more images. π π π π λ λ λ λ (u 2 +v 2 ) Phase transfer function Fig. 3. Rationale for exponential spacing measurement scheme. Plot shows phase transfer functions for exponentially spaced z steps, which ensure a minimum sensitivity of α across a range of frequencies. g(u, v, z 1 ), g(u, v, z 2 ), and g(u, v, z 3 ) are the phase transfer functions at z 1 , z 2 , and z 3 , respectively. The minimum sensitivity plot shows the frequencies which are transferred at the sensitivity higher than α by choosing z 1 , z 2 , and z 3 . Larger z brings more low-frequency sensitivity.

Simulations
We compare the performance of various phase recovery methods that use equal spacing and exponential spacing. In the simulation, the illumination wavelength is set as 632.8nm, and each image has 100 × 100 pixels (pixel size 2µm × 2µm). Equally Spaced Stack has 9 intensity images simulated with a constant defocus step size of 5µm. Exponentially Spaced Stack also has 9 images; however, they are exponentially spaced, at z positions around focus of ±5µm, ±20µm, ±80µm and ±320µm (see Fig. 4). Although both data sets use the same number of images, the exponentially spaced data set contains more information about the low-frequency phase information because it has a higher z range. In order to assess the error performance, the intensity images of the equally and exponentially spaced stacks are corrupted by white noise with SNR ranging from 18.5 to 8 dB (noise variance from 0 to 0.02). Fig. 5 shows the average mean square error (MSE) of the recovered phase over 50 trials as SNR decreases. For the Equally Spaced Stack, higher order TIE performs significantly worse than SGDF TIE, and GP TIE is slightly better than SGDF TIE. For the Exponentially Spaced Stack, we show the results for two possible choices of Higher order TIE: m = 9, which performs better in low noise, and m = 5, which performs better in high noise. GP TIE with exponential spacing clearly exhibits the lowest MSE. This can be explained by the fact that the exponentially spaced data contains more low-frequency phase content in the measurements than the equally spaced data, and there is no trade-off between noise and nonlinearity. Figure 6 gives an example of the recovered phase at SNR of 11.1dB.

Experimental results
We tested our method experimentally in a microscope (magnification 20x, NA=0.5) with filtered white light illumination (center wavelength 650nm, 10nm bandwidth). Data Set 1 comprises 129 images of human cheek cells, equally spaced by a constant small step size dz = 4µm over a large defocus range [−252µm to 252µm] (Fig. 8(a)). In Fig. 7, the GP fitted intensity spatial frequency variations along the propagation direction z for 3 different (u, v) values are shown. Both the measured and fitted curves follow nearly sinusoidal patterns, as predicted by Eq. (3). Note that the delay of each sinusoid is dependent on the absorption of the object at the corresponding spatial frequency. When the frequency (u, v) is high, the intensity spectrum variations diminish for large defocus distance (see plot of (u 60 , v 60 ) in Fig. 7). This is due to the effect of partially coherent illumination, which will be the subject of future work. With our exponentially spaced measurement scheme, GP TIE requires fewer images to be captured in order to obtain a high quality phase result. To demonstrate this, we compare the recovered phase by GP TIE from different subsets of Data Set 1 (Fig. 8) by sampling the image stack along z using various strategies. Figure 8 z-planes. First, we show the best possible phase result, when all 129 images are used. With equally spaced planes, we have two possibilities for reducing the number of images: we can either reduce the defocus range (keeping the step size small) or increase the step size (reducing the defocus range). If we reduce the defocus range, the recovered phase becomes susceptible to low-frequency noise (see the middle of Fig. 8(b)). This is due to the fact that the lowfrequency information of phase is not well captured at small defocus distances. If we instead increase the step size (keeping the defocus range large), high-frequency components are lost due to nonlinearity (see the right of Fig. 8(b)). In order to accurately capture both high and low-frequency information with the same reduced number of images, we need nonlinearly spaced measurements. Near the focus, the step size should be small, yet a large focus range can still be covered with only a few measurements in our exponentially spaced scheme. In Fig. 8(c), we extract images from Data Set 1 according to the exponential spacing scheme described in Section 2.3. As can be seen from Fig. 8(c), the phase results are free of low-frequency noise and also have high resolution. Even after further reducing the dataset to only 9 exponentially spaced images, we obtain a similar result to that with all 129 images. Thus, we have reduced the data capture requirement by more than an order of magnitude, without sacrificing quality.
Having shown that exponential spacing is advantageous, we now compare GP TIE with Higher order TIE, both of which are capable of using unequally spaced data. Data Set 2 comprises 9 images of human cheek cells, exponentially spaced about focus at ±5.7µm, ±11.4µm, ±22.8µm, and ±45.6µm ( Fig. 9(a)). The imaging system has 10x magnification and NA of 0.5. Each image has 350 × 360 pixels of size 0.62µm × 0.62µm. Figure 9(b) shows the recovered phase of Data Set 2 obtained by Higher order TIE and GP TIE. We show the results for higher order TIE with the order of polynomial fitting m equal to 2, 3, and 4. The phase images of m = 2 and 3 have small low-frequency noise but appear blurred. The phase of m = 4 has strong contrast in some regions, but contains low-frequency noise. In contrast, GP TIE does not suffer from this tradeoff, so the phase recovered has less low-frequency noise and also exhibits high contrast. We can see details inside cells in the phase image recovered by GP TIE.

Conclusions
In this paper, we proposed a TIE phase recovery method that exploits prior knowledge of spatial frequency variations of the intensity in order to yield more accurate phase images from exponentially spaced images through-focus. Our proposed method estimates the intensity derivative by means of GP regression. It is robust and stable with noise, while incorporating the a priori knowledge of the sinusoidal patterns predicted when fitting spectrum variations over z.
We derived an exponentially spaced measurement scheme which guarantees a minimum phase measurement sensitivity across the full range of available spatial frequencies with the least possible number of images required. In future work, we will extend the current technique to add more a priori information, such as partial coherence of the illumination and large phase effects.
With the freedom to choose nonlinearly spaced measurement positions, we can develop new strategies to optimize measurements if any prior knowledge of the phase spectrum is known. It should be noted that improving the TIE result by incorporating prior knowledge has been previously considered in [30][31][32]. The difference is that the previous approach considers priors over the phase [30,31] or refractive index distribution [32], whereas here we use a prior on intensity spectral evolution to improve the intensity derivative estimate. We expect the approach considered here should have wide application as it not only can be directly apply to the special cases in [30][31][32], but also to more general situations without constraints on the phase distribution. If any reader is interested in this algorithm, open source code can be obtained by emailing the authors or visiting the website either www.laurawaller.com or www.dauwels.com.