Sparse ACEKF for phase reconstruction.

: We propose a novel low-complexity recursive ﬁlter to efﬁciently recover quantitative phase from a series of noisy intensity images taken through focus. We ﬁrst transform the wave propagation equation and nonlinear observation model (intensity measurement) into a complex augmented state space model. From the state space model, we derive a sparse augmented complex extended Kalman ﬁlter (ACEKF) to infer the complex optical ﬁeld (amplitude and phase), and ﬁnd that it converges under mild conditions. Our proposed method has a computational complexity of N z N log N and storage requirement of O ( N ) , compared with the original ACEKF method, which has a computational complexity of O ( N z N 3 ) and storage requirement of O ( N 2 ) , where N z is the number of images and N is the number of pixels in each image. Thus, it is efﬁcient, robust and recursive, and may be feasible for real-time phase recovery applications with high resolution images.


Introduction
When light propagates, its amplitude and phase evolve according to the wave equation. These amplitude and phase perturbations contain important information about the complex-field in the original (focal) plane. Intensity can be measured easily; however, phase cannot be measured directly and it needs to be reconstructed from intensity images with computational methods. Most solutions for phase recovery do not explicitly account for noise in the measurement, and thus typical methods are not applicable in high noise situations, such as in X-ray imaging or photon-starved situations. Recently, Kalman filters have been used to estimate phase from very noisy intensity images (signal-to-noise ratio less than one) [1]; however the method is very computationally intensive and impractical for images with large pixel counts (e.g. 1 megapixel or greater). In this paper, we solve these problems by developing efficient inference methods for recovering the phase from a set of defocus intensity images, which can be captured either sequentially by moving camera along the optical axis, or all at once by using various focal stack collection methods [2,3]. Such methods can be useful for applications in medical imaging, neuroscience, and materials science, where low noise images are difficult to obtain. Furthermore, since our proposed method is recursive (not iterative), it has the ability to estimate phase in real-time during the measurement sequence, progressively improving the estimate as each new image is captured.
Traditional methods for recovering phase involve complicated interferometric setups, so there is a significant experimental advantage to methods such as ours, which only require a set of intensity images taken with different amounts of defocus. One of the originally proposed algorithms, the Gerchberg-Saxton(GS) method [4,5], uses only two images and treats the problem as convex (which it is generally not), iterating back and forth between two domains to reduce error. In [6], a generalization of the GS method is given for cases of multiple defocus images. These methods tend to be strongly sensitive to the noise in the last image, although techniques exist for averaging at each iteration. Direct methods generally linearize the problem, either by assuming that the object is weak or that the propagation distance is small. It allows one to exploit the Transport of Intensity Equation (TIE) [7], which reveals relationship between phase and the first derivative of intensity with respect to the optical axis. The first derivative of intensity is usually approximated by finite difference method [8,9], and it is not robust to severe noise. A few statistical approaches have been proposed as well. An approximation to the maximum likelihood estimator is derived in [10,11]. However, it easily gets stuck in local maxima, and sometimes leads to poor results.
In [1], an augmented complex extended Kalman filter (ACEKF) was used to solve for phase in the presence of significant noise corruption. The Kalman filter is a well-known recursive algorithm for estimating dynamically changing quantities from a set of noisy measurements. In the linear case, with Gaussian noise statistics, the Kalman filter is the optimal estimation method. However, intensity is a nonlinear measure of complex field, and noise in images generally follows Poisson statistics, so the solution will not be provably optimal. Extended Kalman filters add a linearization step in order to account for nonlinear measurements. Further, we use a complex EKF in order to deal with complex numbers. The resulting ACEKF outperforms many previous algorithms, at the cost of very high computational complexity of O(N z N 3 ) and storage requirement of O(N 2 ). The memory storage requirements are particularly limiting -for example, an image with pixel counts on the order of megapixels will require terabytes of data to be stored and processed -thus, the algorithm is infeasible for large images and impractical for real-time applications. In [12], a diagonalized complex extended Kalman filter (diagonalized CEKF) was proposed to alleviate those issues, without jeopardizing the reconstruction accuracy. The diagonalized CEKF is iterative: it needs to cycle through the set of intensity images repeatedly, yielding a more accurate phase reconstruction after each cycle. The computational complexity increases with each cycle.
The sparse ACEKF proposed here uses the same augmented state space model as ACEKF. With the assumption that the phase is small, we derive theorems to simplify the update equations of ACEKF. It eventually results in a low-complexity noise-robust phase reconstruction algorithm. The proposed method has a computational complexity of N z N log N, while ACEKF requires N z N 3 . However, our method still achieves better phase reconstruction than both ACEKF and the diagonalized CEKF.
This paper is organized as follows. In the next section, we turn the propagation and observation models from scalar wave optics into an augmented state space model in signal processing. In Section 3, we derive the proposed sparse ACEKF algorithm. In Section 4, we analyze the convergence and stability of the new ACEKF. In Section 5, we compare results of different methods using synthetic and experimental data. We offer concluding remarks in Section 6.

Problem description and state space model of the optical field
We aim to estimate the 2D complex-field A(x, y, z 0 ) at the focal plane z 0 , from a sequence of noisy intensity images I(x, y, z) captured at various distance z 0 , z 1 , z 2 , ..., z N z . We assume a linear medium with homogeneous refractive index and coherent (laser) illumination, such that the complex-field at z 0 fully determines the complex-field at all other planes. The complex optical field at z is A(x, y, z) = |A(x, y, z)|e iφ (x,y,z) , where |A(x, y, z)| is the amplitude, and φ (x, y, z) is the phase. Propagation is modeled by the homogeneous paraxial wave equation: where λ is the wavelength of the illumination, and ∇ ⊥ is the gradient operator in the lateral (x, y) dimensions. The noisy measurements I(x, y, z) usually adhere to a (continuous) Poisson distribution: where γ is the photon count detected by the camera. The measurement at each pixel I(x, y, z) is assumed statistically independent of any other pixel (conditioned on the optical field A(x, y, z)).
We discretize the optical field A(x, y, z) as a raster-scanned complex column vector a n , and similarly discretize the intensity measurement I(x, y, z) as column vector I n . We denote by b(u, v, z) the 2-D Fourier transform of A(x, y, z). The column vector b n is again raster-scanned from b(u, v, z), and hence can be expressed as b n = Ka n , where K is the discrete Fourier trans-form matrix. Since K is unitary, we can write KK H = K H K = U (with normalization), where U is the identity matrix or unit matrix and K H denotes the hermitian of K.
Let us assume the distance between two consecutive planes is a constant ∆z. Then, we can define the propagation matrix from z n−1 to z n as [13]: where L x and L y are the width and height of the image, respectively. The relation between two images with distance ∆z in Fourier domain can be written as: In order to fit within the Kalman filter framework, we approximate the Poisson observation model as in Eq. (2) with a Gaussian distribution of the same mean and covariance. In particular, we consider the approximate observation model: where v is a Gaussian vector with zero mean and covariance R = γ diag(a * n )diag(a 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 . The nonlinear observation model in Eq. (5) is linearized as [14]: whereâ n is the state predicted from the previous n − 1 observations, and Eq. (6) is the first order Taylor series expansion of Eq. (5) with respect toâ n . Summarizing, the augmented state space model is given as: observation: where We adopt the augmented state space model because the ACEKF with the state augmented outperforms CEKF in both estimation error and convergence [14]. If the state is not augmented, such as in CEKF algorithm, the covariance E[b n b H n ] is not sufficient for a complete description of the complex state. Since the state in ACEKF is augmented from together provide a more comprehensive statistical description of the complex state.

State estimation by sparse augmented complex extended Kalman filter
The state covariance matrix of the augmented state has the form: From the update equations of ACEKF [1,15], we have the following steps: 3. Update: The size of S Q n or S P n is N 2 , where N is the total number of the pixels in the image. Thus, the covariance matrix dominates the memory requirements. The inversion of the covariance matrix has a computational complexity of O(N 3 ) in each step. Both the storage requirement and computational burden make the above update algorithm impractical for real applications.
Here we make some mild assumptions and derivations, resulting in a low-complexity algorithm with reduced storage requirement. After some derivations (see more in Appendices A and B), we prove Lemma 1 as well as Theorems 1 and 2. Theorem 1 is derived from the ACEKF update Eqs. (11)- (14). It shows that if the covariance matrix at step n − 1 has the form S Q n−1 = Q n−1 and S P n−1 = P n−1 E, where Q n−1 and P n−1 are diagonal, then the covariance matrix at step n has the same form as that of step n − 1 (see Eqs. (23)-(24)). Therefore, once the first covariance matrix is initialized as S Q 0 = Q 0 and S P 0 = P 0 E, then the subsequent matrices will keep the same form. Since the matrices in Theorem 1 are diagonal, the computational complexity in Theorem 1 is O(N). Theorem 2, derived from the Eqs. (15)- (16), provides a new Kalman gain and update formula using the diagonal matrices Q n and P n . Lemma 1. If a matrix M is diagonal and its diagonal entries are rotationally symmetric in 2-D, then EME = M, where E = KK T , and K is the Discrete Fourier Transform Matrix.
where Q n−1 and P n−1 are diagonal, then we derive from the ACEKF update Eqs.
P n = HP n−1 H (20) Update: S Q n = Q n (23) where Q n and P n are diagonal, and q = 1 γ . Theorem 2. The Kalman gain and update formula for the state are (1) Initialization b 0 , Q 0 and P 0 . Q n =Q n − (Q n +P n )(Q n +P n + (Q n ) * + (P n ) * + qI) −1 (Q n + (P n ) * ) (31) The resulting algorithm, referred to as the sparse ACEKF, is summarized in Table 1. Matrices Q n and P n are diagonal, and hence they can be stored as two vectors. The storage burden of Eqs. (13)- (14) in the update step is reduced from N 2 to N. The inverse of J n in Eq. (33) can be computed by Fast Fourier Transform (FFT), which has a computational complexity of O(N log(N). Since Q n and P n are diagonal, the matrix multiplications and inversions in Eqs.

Analysis of the convergence of the augmented complex extended Kalman filter
The convergence of ACEKF, for the complex state space model, has not been well analyzed in the literature. However, the convergence and stability of the extended Kalman filter (EKF) for the real state is well studied [16][17][18]. Thus we first convert the ACEKF state space model into its equivalent dual form as a real state space model (EKF), then we use existing methods to analyze the convergence of the dual form.
We follow [14] to derive the dual real form of the ACEKF state space model (7)- (8). Define where U is identity matrix. When multiplied by M, a complex variable is converted to its corresponding real form: where Re(b n ) and Im(b n ) are the real and imaginary parts of b n , respectively. The complex model (7)- (8) becomes a real value model: observation: where By applying the theory of EKF convergence in [16] to our real value model, we obtain Theorem 3 which provides the convergence and stability of the model. 1. There exist positive numbersā,c,s, s, v > 0, such that for every n ≥ 0: (43)

2.
A n is non-singular for every n. In Theorem 3, . denotes the spectral norm of matrices. Assuming the initial estimation error and the noise are small enough, let us check each condition of Theorem 3. From the definition in Eq. (3), the propagation matrix H is non-singular and determined by the propagation distance ∆z. Since M is full rank, the matrix A n given in Eq. (38) is bounded and non-singular. Therefore, the inequality (40) and the second condition of Theorem 3 are satisfied. The linearization in Eq. (6) is a first-order truncation of Taylor series expansion of a quadratic function. Therefore, the linearization has a bounded error and the real observation matrix C n of Eq. (37) is bounded as well, as required by the inequality (41) and the third condition of Theorem 3. The matrix R is the covariance matrix of the Poisson noise in Eq. (5), and it is bounded as the amplitude of the state is bounded.
The inequality (42) requires the covariance matrix to be bounded for every n. In [16], it is shown that the boundedness of the covariance matrix is closely related to observability and detectability properties of the EKF model. However, it is difficult to prove the observability of the model directly, because A n and C n are based on the observed data at each step. We calculate the rank of the observability matrix and the bound on the covariance matrix by testing our algorithm on the synthetic data (see Data Set 2 in section 5). We initialize the covariance matrix Q 0 and P 0 in Table 1 with 30 different values. The matrix Q 0 is initialized as qU (U is identity matrix), with q ranging from 70 to 7000, and P 0 is initialized as pU, with p ranging from 50 to 5000. Figure 1 shows that the minimum singular value of the observability matrix is larger than zero, so the observability matrix is full rank at each step. Figures 2 and 3 show the maximum and minimum singular value of the covariance matrix at each step with different initializations. The covariance matrix is bounded, as required by inequality (42).  Fig. 1. Minimum singular value of observability matrix with 30 different initializations of the covariance matrix Q 0 and P 0 . The matrix Q 0 is initialized as qU, with q ranging from 70 to 7000, and P 0 is initialized as pU, with p ranging from 50 to 5000.

Experimental results
We consider three sets of data to assess the performance the augmented Kalman filter. Data Set 1 consists of 100 images of size 100 × 100 pixels artificially generated to simulate a complex field propagating from focus in 0.5 µm steps over a distance of 50 µm with illumination wavelength of 532 nm. Pixels are corrupted by Poisson noise so that, on average, each pixel detects γ = 0.998 photons. Data Set 2 comprises 50 images of size 150 × 150 pixels acquired by a microscope. The wavelength was again 532 nm, and the defocused intensity images were captured by moving the camera axially with a step size of 2 µm over a distance of 100 µm. Data Set 3 has 101 images of size 492 × 656 pixels acquired by a microscope. The light source is partially coherent, and filtered by a narrow-band color filter of center wavelength 633nm. The images were captured by moving the camera axially with a step size of 2 µm. Figure 4 shows the images of synthetic data (Data Set 1) and experimental data (Data Set 2 and Data Set 3).   N log N), and it takes 0.40 seconds to process the 100 (full) images, thus giving a speedup factor of 35000x. As illustrated in Table 2, the computational complexity of the diagonalized CEKF is lower than that of ACEKF. However, the latter yields better results in terms of phase error. In order to reduce the error of the diagonalized CEKF, forward and backward sweeps (iterations) are applied in [12]. However, the iteration increases the computational complexity linearly, and makes the method no longer recursive. The sparse ACEKF method has an intensity error of 0.0071, and a phase error of 0.0143 (radian). Compared with the diagonalized CEKF, the sparse ACEKF has the same computational complexity and storage requirement, but returns more accurate images.

Synthetic data
We assess the quality of reconstruction by root mean square error (RMSE). The proposed sparse ACEKF has an error near to that of ACEKF, while the recovered phase and intensity images of the sparse ACEKF in Fig. 5 might look better. The images recovered by ACEKF exhibit a block effect as straight lines crossing the images, whereas the result of sparse ACEKF is free of such block effect. The sparse ACEKF has a much lower complexity, which avoids the need of dividing the images into independent blocks. The images recovered by ACEKF and the diagonalized CEKF contain traces of phase in the intensity images due to errors. However, the trace of phase is mostly removed in the estimated intensity image from the sparse ACEKF.

Experimental data
In Fig. 6, we compare the estimated intensity and phase images of Data Set 2 using the ACEKF, the diagonalized CEKF, and the sparse ACEKF. Stripes in the phase image recovered by the diagonalized CEKF look darker, while the stripes in the recovered phase image of the sparse ACEKF method have stronger contrast. In Fig. 7(a), we show the recovered phase of the Data Set 3 by ACEKF, the diagonalized CEKF, and the sparse ACEKF. The real depth of the sample in Data Set 3 is around 75 ± 5nm. The proposed method takes 20.24 seconds to process 101 images of size 492 × 656. However, the ACEKF method takes 54.15 hours and each image is separated into 117 pieces of 50 × 50 blocks. In Fig. 7(b), we compare the depth across the black line of the recovered phase in Fig. 7(a). The sparse ACEKF method shows a result much closer to the true value, compared to the ACEKF and the diagonalized CEKF. However, the reconstructed phase image obtained by sparse ACEKF contains low-frequency noise, especially at the edges. This shadow effect may be attributed to the partial coherence of the light, which is not incorporated into our model (we assume full coherence). The low frequency issue may be alleviated through a nonlinear diffusion filter, which we will explore in future work.

Conclusions
The proposed method efficiently recovers phase and amplitude from a series of noisy defocused images. Although a constant defocus step size is used for the through-focal series in this paper, it can be generalized to cases of non-constant step-size. The method can directly deal with large data sets and high-resolution images because of its low computational complexity and storage requirement. It can process the images recursively in real-time after the data has been captured, or even during the capture sequence. For example, with further work, this method could form the basis for an adaptive phase imaging method, in which the current estimate of the phase informs the choice of the next measurement plane (defocus distance). Since the optimal measurement planes will depend on the object itself [19], real-time estimation during capture could lead to optimization not only of the processing, but also optimization of the measurement scheme itself. Furthermore, due to the scalability of the wave equations and the simplicity of the measurement technique, this method could find use in phase imaging beyond optical wavelengths (for example, X-ray or neutron imaging), where high-quality images are difficult to obtain and noise is significant and unavoidable.