Partially coherent phase imaging with simultaneous source recovery

: We propose a new method for phase retrieval that uses partially coherent illumination created by any arbitrary source shape in K¨ohler geometry. Using a stack of defocused intensity images, we recover not only the phase and amplitude of the sample, but also an estimate of the unknown source shape, which describes the spatial coherence of the illumination. Our algorithm uses a Kalman ﬁltering approach which is fast, accurate and robust to noise. The method is experimentally simple and ﬂexible, so should ﬁnd use in optical, electron, X-ray and other phase imaging systems which employ partially coherent light. We provide an experimental demonstration in an optical microscope with various condenser apertures.


Introduction
The 'phase problem' can be summarized as the need to recover a sample's complex-field, ψ(x, y), from measurements of intensity, I(x, y) = |F{ψ(x, y)}| 2 , in an optical system described by the operator F{·}. At the sample, intensity images contain no information about the phase of ψ(x, y). However, with a suitable choice of optical system, F, the sample's phase variations cause intensity variations at the camera. Thus, the first goal of phase imaging is to design an optical system for transferring phase into intensity information (the forward problem). In fact, any complex transfer function will induce phase contrast; however, the ideal choice of F involves tradeoffs between experimental complexity, cost, accuracy, and the sample itself [1]. Various systems have been proposed, broadly categorized as interferometric [2], grating-based [3] and defocus-based methods [4,5]. Here, we choose defocus as our F, because it is convenient and experimentally simple, even for lensless imagers or non-traditional wavelengths (e.g. X-ray [6], TEM [7]). The transfer function for defocus is well-studied [8,9] and can be tuned by varying defocus distance, z. Though we use here a stack of images at varying defocus ( Fig. 1), we note that our algorithm is general and can be adapted to any appropriate optical system.
The second goal of phase imaging is to invert the measurements in order to recover complex field (the inverse problem). Being a nonlinear problem, one can employ nonlinear optimization (e.g. iterative methods [4]). These methods are simple to implement using convex procedures, but the problem is non-convex, so they are not proven to converge. To solve this, one can lift the problem into higher dimensions where it is linear [10], but the computational requirements become prohibitive. Alternatively, the problem can be linearized by assuming a weak object [1,11] or a small defocus distance (via the Transport of Intensity Equation (TIE) [5]). The TIE gives a direct solution for phase as long as defocus is small. Unfortunately, images taken with small defocus have limited phase contrast, causing noise instabilities which are only partially mitigated by more z steps [12][13][14] or non-linear spacing [15].
In our work, we use a Kalman filtering approach [16], which provides the near-optimal phase solution, even in severe noise. The method can be thought of as a hybrid between iterative methods and TIE, in that it uses a nonlinear forward model (as in iterative methods), but updates the phase estimate according to a linearization at each step (as in TIE). Thus, we get the benefits of using a large range of defocus distances, with good convergence. Because the Kalman filter stores an estimate of the noise covariance for all pixels, it provides pixel-wise adaptive filtering, which significantly improves noise performance. Unfortunately, this means that one must store and manipulate a 4D covariance matrix of size N 2 (where N is the total number of pixels), so memory requirements are prohibitive. We recently solved this problem by proving that the covariance matrix should be sparse [17], reducing both memory and computational complexity by many orders of magnitude. Note that the sparsity of the covariance matrix does not imply sparsity of the object [18][19][20]. The Kalman filter also provides the ability to easily modify the forward model for any complex transfer function. The result is a general method that can achieve nanometer level phase sensitivity, converges reliably and is robust to noise.
This work focuses on an important extension of our Kalman filter method which enables partially spatially coherent illumination [21]. Coherence has been discussed extensively in the context of TIE, since it plays a role in defocus. The TIE method is often cited as being robust to partial coherence [22], but this is only true if the defocus distance is small and the source is circular [23,24]. When non-circular sources are used, the phase result is incorrect [25], so coherence effects must be explicitly incorporated [11,[26][27][28]. Further, as propagation distances increase beyond the 'small defocus' regime, partial coherence becomes more problematic [24,[29][30][31]. We seek here a general algorithm which can incorporate illumination coherence efficiently across a large range of defocus, in order to achieve better phase recovery.
We restrict our analysis to imaging systems with incoherent sources placed in the Fourier plane of the sample (Köhler geometry). In this configuration, illumination is spatially invariant and coherence is controlled by source shape (via the Van Cittert-Zernike theorem). Smaller sources provide more spatial coherence, usually at the cost of reduced photon flux, lower resolution and unwanted speckle. In some cases, partially coherent light is unavoidable. For example, in lithography applications [32], current trends towards source-mask optimization (SMO) [33,34] mean that non-traditional source shapes are the norm. Thus, incorporating partial coherence may be not only desirable, but necessary, for phase retrieval.
Here, we propose a modification to our Kalman filter method in which coherence effects are modeled as a simple convolution, allowing fast and accurate phase recovery with arbitrary source shape. Further, by adding an extra optimization procedure and some assumptions on the source, we show that it is possible to also recover the source size (coherence) simultaneously, meaning that one need not know the illumination coherence a priori.

Theory
Our method is based on the Kalman filter presented in [16,17], which we refer to here as the 'coherent Kalman filter'. The inherent assumption of spatial coherence made in this previous work breaks down in the case of large or non-circular sources, causing phase errors (Fig. 2). Here, we modify our algorithm to account for coherence explicitly. We note that the phase we solve for is well-defined as that of the sample, not the illumination [35].
Generally, modeling spatial coherence for a 2D beam requires a 4D function, which is difficult to know or measure a priori [36]. However, if the source is in Fourier space (Köhler configuration), the Van Cittert-Zernike theorem shows that we can model coherence with a much simpler 2D function related to the source shape. Upon defocus, each source point creates the same coherent diffraction pattern, but shifted according to the source position. Neglecting pupil filtering effects, we can therefore treat the forward model in the same way as the coherent case, but with an extra convolution step for each defocus distance. This convolution model was first introduced for phase recovery in [21] and later used in other methods [37,38].
Specifically, the convolution model describes the partially coherent intensity I(x, y, z) at defocus distance z as a 2D convolution between the intensity that would have been measured with coherent illumination |A(x, y, z)| 2 and a scaled source distribution S(x, y) [24,39]: where f is the focal length of the condenser lens. Based on this model, we formulate a statespace description and develop a sparse Kalman filtering algorithm for phase recovery.

Partially coherent phase recovery with known source shape
First, we describe our method for estimating the complex field at focus A(x, y, z = 0) given a known source shape S(x, y). Our dataset consists of multiple intensity images, I(x, y, z n ), taken at different defocus, where z n is the intensity measurement position along z axis. In order to fit into the framework of our state-space model, we first discretize the continuous model. We rasterize the 2D continuous x −y plane A(x, y, z n ) into a column vector a n and denote the Fourier transform of A(x, y, z n ) with b n . Propagation in Fourier domain is written as: where K is the discrete Fourier transform (DFT) matrix, and H n is a diagonal matrix describing wave propagation by a distance z n − z n−1 .
Assuming that intensity measurements incur Gaussian noise [17], Eq. (1) is discretized as: where the vector I n is the discretized form of I(x, y, z n ), |a n | 2 is for discretizing |A(x, y, z n )| 2 , C n is a circulant matrix describing the convolution by S(− f z n x, − f z n y), and v is Gaussian vector noise in the measurement with zero mean and covariance R n = diag(C n |a n | 2 ). Note that each entry in the vector |a n | 2 is the square of the absolute value of the corresponding entry in a n . From the eigenvalue decomposition property of circulant matrices, C n = K H S n K. Here S n is a diagonal matrix with its diagonal entries equal to the Fourier transform of the scaled source S(− f z n x, − f z n y). We now derive a complex extended Kalman filter, by linearizing the nonlinear observation model [40]: whereâ n is the state predicted from the previous n − 1 observations, and Eq. (5) is the first order Taylor series expansion of Eq. (4) with respect toâ n . Here a * n denotes the complex conjugate of a n , and diag(a * n ) is a diagonal matrix with its corresponding diagonal entries equal to the elements in the vector a * n . Summarizing, the augmented model is: observation: where v is white Gaussian noise with covariance R n , R n = diag(C n |â n | 2 ), and J n = C n diag(â * n )K H .
We modify the state space Eq. (6) - (7), to accommodate the sparse version in [17]. The resulting partially coherent phase retrieval algorithm is summarized in Table 1. The matrices Q 0 and P 0 are initialized as diagonal matrices, so Q n and P n remain as diagonal matrices during the update in Eq. (10)- (12). Because the matrices H n , S n , Q n , and P n Eqs. (10)-(12) are diagonal, the computational complexity of matrix multiplication in the algorithm is reduced to be linear in the size N of the image. Comparing to O(N 3 ) in the standard extended Kalman filter, the computational complexity of the sparse algorithm is O (N log N). Thus, the algorithm is computationally efficient and can run in nearly real-time, even with complicated source shapes. (1) Initialization of b 0 , Q 0 and P 0 .
2.2. Partially coherent phase recovery with unknown source shape The algorithm described in Table 1 assumes a known source shape, S(x, y). Any errors in measurement or estimation of the source will result in errors in the recovered phase. To remedy this, we extend our method to the case where the source shape is unknown. Because the sample phase information is imprinted on each of the coherent modes (source points), the resulting through-focus intensity stack contains sufficient information for recovering both phase and source shape. One could, for example, use blind deconvolution to recover the coherent intensity stack from the partially coherent stack, then recover phase from the coherent stack.
Here, we aim to solve for both simultaneously within a single algorithm. Consider the simplified case of source size estimation for a known complex field a 0 and source S(x, y) with parameters θ (for example, θ is the radius of a circular source). The problem is to estimate θ from the measured partially coherent intensity I n . The estimation of θ can be formulated as a nonlinear optimization problem: where |a n | 2 can be obtained by defocusing a 0 . The solution finds an optimal θ to minimize the squared error between measurements I n and the estimation C n (θ )|a n | 2 .
To estimate both the complex field a 0 and θ simultaneously, we propose an iterative update to a (i) 0 and θ (i) , where i denotes the iteration number. The algorithm for unknown source shape can be realized by iterating through the following steps: (i) For the first iteration i = 0, initialize a (i=0) 0 and θ (i=0) .
(ii) With the current estimated parameters θ (i) , recover the complex field a This iterative method provides an approximation to the optimal solutions for a 0 and θ in terms of the least squares error between the measurement I n and the estimated intensity C n (θ )|a n |. By iterating between Step (ii) and Step (iii), the estimation of the complex-field and source shape parameters are updated gradually, giving an approximation of optimal value when the iterations stop. Although the Kalman filter is an optimal estimator when the state-space is linear [41,42] and noise is Gaussian, our method in Table 1 cannot be proven to be the optimal estimator for the extended Kalman filter case [43,44]. However, it performs well in practice and has good convergence properties, even with unknown source shapes.

Experimental results
We demonstrate our method experimentally in a Nikon TE300 microscope (Fig. 1) with broadband illumination centered around wavelength 550 nm. We capture intensity focal stacks of cheek cell samples, using different condenser apertures to vary the source shape (see Fig. 2, top row). The first three shapes were obtained by progressively opening up the condenser aperture to reduce spatial coherence, and the annular ring was obtained by choosing the Ph1 condenser. In order to validate our algorithm, images of the source plane were captured at the side port of the microscope. The measured numerical apertures (NA) of the first three circles are given in Table 2, with the two numbers for Ph1 representing the NA of the outer and inner circles. For each shape, the sample was defocused symmetrically about focus at 9 z-planes from 1µm to 16µm, exponentially spaced as in [15]. The second row of Fig. 2 shows the phase images recovered by our previous method, the coherent Kalman filter [17]. In the case of Size 1, the light source can be approximately viewed as coherent since the source size is small. The coherent Kalman filter thus recovers an accurate phase result. As the size of light source increases, however, the level of coherence reduces and the recovered phase becomes more blurred, degrading the phase result. Next, we test our partially coherent Kalman filter with known (measured) source shapes, showing that it can recover an accurate phase result for each case. Its recovered phase images do not suffer from blurring and contain fine details for all source sizes and shapes, as shown in the third row of Fig. 2. Note that the phase recovered from the larger sources has slightly better resolution than the coherent case, due to the partial coherence of the illumination. Finally, we show results for our algorithm with both phase and source size recovery. The algorithm successfully recovers a source NA similar to the measured values (see Table 2) and we can see in Fig. 3 that the estimated NA converges after only a few iterations. The phase recovered with the unknown source algorithm also shows good fidelity (Fig. 2, bottom row). In fact, the phase with the unknown source algorithm contains somewhat better details than those with known source shapes, suggesting that the recovered source shapes are more accurate than the measured ones. This suggests that the unknown source algorithm may be used (with the measured source as an initial guess), even when the source shape is known.
In our experiments, we find that the annular aperture provides the best phase result for several reasons. It achieves the improved spatial resolution that comes with using partially coherent illumination, but does not lose contrast through-focus due to blurring, as in the case of the larger circular aperture. Another advantage of using partially coherent illumination is that images do not suffer from speckle noise or out-of-focus artifacts, as the coherent situation. An interesting topic for future work might be to derive the optimal source shape for a given experiment.

Conclusion
We have demonstrated new Kalman filtering methods to recover phase and illumination coherence from a stack of defocused intensity images. Our algorithm can either incorporate a known source of arbitrary shape, or estimate the source size computationally. This embedded calibration step is particularly helpful for making our algorithm more general and easy to implement.
We expect it will find use not only in biological microscopes, but also in X-ray and TEM imaging, where source coherence is difficult to control and may change between experiments [45]. The ability to work with partially coherent illumination, along with the simplicity of phasefrom-defocus techniques, eases constraints on experiments, opening up new application areas.