PARTIALLY COHERENT PHASE RECOVERY

Abstract

A system and method for incorporating partially coherent illumination models into the problem of phase and amplitude retrieval from a stack of intensity images. The recovery of phase could be realized by many methods, including Kalman filters or other nonlinear optimization algorithms that provide least squares error between the measurement and estimation.

Claims

1. A system for phase recovery of an intensity image with partially coherent illumination, the system comprising: (a) a computer processor; and (b) a memory storing instructions executable on the processor; (c) the instructions, when executed, performing the steps comprising: (i) mathematically modeling partial coherence with respect to the intensity image; (ii) formulating a state-space Kalman filter from the partial coherence model; and (iii) applying a standard Kalman filter to the state-space Kalman filter to estimate a phase with respect to the intensity image.

2. A system as recited in claim 1, wherein the image comprises a stack of defocused intensity images; and wherein a propagation transfer function is applied between captured intensity images.

3. A system as recited in claim 1, further comprising: applying a sparse model to the Kalman filter to reduce computational complexity.

4. A system as recited in claim 3, wherein the sparse model optimizes according to:
.sub.n=K.sup.H{circumflex over (b)}.sub.n
Q.sub.n={circumflex over (Q)}.sub.n({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*+({circumflex over (P)}.sub.n)*+qS.sup.2).sup.1({circumflex over (Q)}.sub.n+({circumflex over (P)}.sub.n)*)
P.sub.n={circumflex over (P)}.sub.n({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*+({circumflex over (P)}.sub.n)*+qS.sup.2).sup.1({circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*)
b.sub.n={circumflex over (b)}.sub.n+(Q.sub.n+P.sub.n)Kdiag(.sub.n)C.sub.nR.sup.1(I.sub.n|.sub.n|.sup.2); where .sub.n is the state predicted from the previous n1 observations, I.sub.n and b.sub.n are discretized column vectors, C.sub.n describes convolution of scaled source intensity distribution, and K is the discrete Fourier transform (DFT) matrix.

5. A system as recited in claim 1, wherein the instructions are configured for estimating a phase with respect to a partially coherent intensity image generated from one or more of: optical microscopes, electron imaging systems, tabletop X-ray imaging systems, synchrotron imaging systems, and lithography aerial imaging systems.

6. A system as recited in claim 1, wherein the instructions are configured for estimating phase for an image having an arbitrary shape of Khler illumination.

7. A system as recited in claim 1, wherein the image comprises a stack of intensity images; and wherein phase or amplitude filters or other linear optical systems are applied as a transfer function between captured intensity images.

8. A method for phase recovery of an intensity image with partially coherent illumination, the method comprising: (a) mathematically modeling partial coherence with respect to the intensity image; (b) formulating a state-space Kalman filter from the partial coherence model; and (c) applying a standard Kalman filter to the state-space Kalman filter to estimate a phase with respect to the intensity image; (d) wherein said method is performed by executing programming on at least one computer processor, said programming residing on a non-transitory medium readable by the computer processor.

9. A method as recited in claim 8, wherein the image comprises a stack of defocused intensity images; and wherein a propagation transfer function is applied between captured intensity images.

10. A method as recited in claim 9, further comprising: applying a sparse model to the Kalman filter to reduce computational complexity.

11. A method as recited in claim 10, wherein the sparse model optimizes according to:
.sub.n=K.sup.H{circumflex over (b)}.sub.n
Q.sub.n={circumflex over (Q)}.sub.n({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*+({circumflex over (P)}.sub.n)*+qS.sup.2).sup.1({circumflex over (Q)}.sub.n+({circumflex over (P)}.sub.n)*)
P.sub.n={circumflex over (P)}.sub.n({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*+({circumflex over (P)}.sub.n)*+qS.sup.2).sup.1({circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*)
b.sub.n={circumflex over (b)}.sub.n+(Q.sub.n+P.sub.n)Kdiag(.sub.n)C.sub.nR.sup.1(I.sub.n|.sub.n|.sup.2); where .sub.n is the state predicted from the previous n1 observations, I.sub.n and b.sub.n are discretized column vectors, C.sub.n describes convolution of scaled source intensity distribution, and K is the discrete Fourier transform (DFT) matrix.

12. A method as recited in claim 8, wherein the method is configured for estimating a phase with respect to a partially coherent intensity image generated from one or more of: optical microscopes, electron imaging systems, tabletop X-ray imaging systems, synchrotron imaging systems, and lithography aerial imaging systems.

13. A method as recited in claim 8, wherein the method is configured for estimating phase for an image having an arbitrary shape of Khler illumination.

14. A system as recited in claim 8, wherein the image comprises a stack of intensity images; and wherein phase or amplitude filters or other linear optical systems are applied as a transfer function between captured intensity images.

15. A system for phase recovery of an intensity image with partially coherent illumination, the system comprising: (a) a computer processor; and (b) a memory storing instructions executable on the processor; (c) the instructions, when executed, performing the steps comprising: (i) mathematically modeling partial coherence with respect to one or more intensity image having partially coherent illumination; (ii) applying a nonlinear least square error function to fit over the one or more measured intensity images; and (iii) estimating a phase with respect to the one or more intensity images.

16. A system as recited in claim 15, wherein applying a nonlinear least square error function comprises finding an optimal estimation to minimize an average square error between measurement and predicted measurement (cost function) with respect to the one or more intensity images.

17. A system as recited in claim 16, wherein applying a nonlinear least square error function comprises: minimizing the cost function using nonlinear least square error method; and deriving a gradient and Hessian.

18. A system as recited in claim 16, further comprising: computing a multiplication of the Hessian with vector with Fast Fourier transform (FFT).

19. A system as recited in claim 15, wherein one or more of the following algorithms are used to solve for said phase: gradient descent, L-BFGS, Newton's, Gauss-Newton, Levenberg-Marquardt, or other nonlinear conjugate gradient algorithms.

20. A system as recited in claim 17, wherein the gradient and Hessian are computed according to: $\begin{array}{c}{F}^{}\left(b,\stackrel{\_}{b}\right)=\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.\left[\begin{array}{c}-A_n^H.Math.\mathrm{diag}\left(A_n.Math.b\right).Math.C_n\\ -A_n^T.Math.\mathrm{diag}\left({\stackrel{\_}{A}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\end{array}\right]\left[I_n-C_n^H.Math.{.Math.A_n.Math.b.Math.}^2\right];\\ F^{}\left(b,\stackrel{\_}{b}\right)\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.\left[\begin{array}{cc}\begin{array}{c}{\stackrel{\_}{H}}_n.Math.K.Math..Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K.Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.K^H.Math.H_n\end{array}& \begin{array}{c}{\stackrel{\_}{H}}_n.Math.K.Math..Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.K.Math..Math.{\stackrel{\_}{H}}_n\end{array}\\ \begin{array}{c}{H}_n.Math.K^H.Math..Math.\mathrm{diag}\left(K.Math..Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K.Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.K^H.Math.H_n\end{array}& \begin{array}{c}{H}_n.Math.K^H.Math..Math.\mathrm{diag}\left(K.Math..Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.K.Math..Math.{\stackrel{\_}{H}}_n\end{array}\end{array}\right].\end{array}$

21. A system as recited in claim 20, further comprising solving the linear equations system: $F^{}\left(b,\stackrel{\_}{b}\right)\left[\begin{array}{c}x\\ \stackrel{\_}{x}\end{array}\right]=-F^{}\left(b,\stackrel{\_}{b}\right).$

22. A system as recited in claim 21, wherein the linear equation is solved as a function of a conjugate gradient method iteratively by computing the gradient and the multiplication of the Hessian with a vector with FFT.

Description

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWING(S)

[0026] The technology described herein will be more fully understood by reference to the following drawings which are for illustrative purposes only:

[0027] FIG. 1 shows a schematic diagram of the image processing system of the present invention implemented with a 4f microscope.

[0028] FIG. 2 shows a high-level flow diagram for a method of obtaining the state-space Kalman filter in accordance with the present description.

[0029] FIG. 3A shows a stack of simulated noisy intensity images with partially coherent illumination.

[0030] FIG. 3B shows an image phase recovered by sparse ACEKF which ignores partial coherence.

[0031] FIG. 3C shows an image phase recovered by use of the Kalman filtering method of the present description which includes partial coherence.

[0032] FIG. 4A shows biological intensity images through focus, with partially coherent illumination.

[0033] FIG. 4B shows an image phase recovered by sparse ACEKF assuming fully coherent illumination.

[0034] FIG. 4C shows an image phase recovered by Kalman filtering with knowledge of the source shape.

[0035] FIG. 5 shows a series of images for comparison of experimental results by coherent Kalman filter and partially coherent Kalman filter. In the top row are source shapes taken by a 2f system at the microscope's side port. The middle row shows recovered phase by the coherent Kalman filter. The bottom row shows recovered phase by the partially coherent Kalman filter of the present description.

DETAILED DESCRIPTION

1. System Configuration

[0036] Referring to FIG. 1 a schematic diagram of the image processing system 10 of the present invention is shown implemented in a Khler system configuration with a 4f microscope 30.

[0037] When an object is illuminated by a coherent source, the optical field at z is described with A(x, y, z)=|A(x, y, z)|e.sup.i(x, y, z), where |A(x, y, z)| is the amplitude, and (x, y, z) is the phase. Free space wave propagation is modeled with wave propagation Eq. 1:

[00001]$\begin{array}{cc}\frac{A\left(x,y,z\right)}{z}=\frac{i.Math..Math.}{4.Math..Math.}.Math.{}_{}^2.Math.A\left(x,y,z\right),& \mathrm{Eq}..Math.1\end{array}$

where is the wavelength of the illumination, and .sub. is the gradient operator in the lateral (x, y) dimensions.

[0038] To include partially coherent illumination into the Kalman filter model of the present description, an incoherent extended source is assumed in the Khler system configuration. The intensity I(x, y, z) defocused by z can be written as a convolution between the intensity from coherent illumination |A(x, y, z)|.sup.2 and a scaled source intensity distribution S(x, y):

[00002]$\begin{array}{cc}I\left(x,y,z\right)={.Math.A\left(x,y,z\right).Math.}^2.Math.S\left(-\frac{f}{z}.Math.x,-\frac{f}{z}.Math.y\right),& \mathrm{Eq}..Math.2\end{array}$

where f is the focal length of the condenser lens 40 of microscope 30. Light from incoherent white light source 32 is filtered by a narrow-band color filter 34, and is placed after collimator 38 at the front focal plane of the condenser 40. The spatial coherence of the illumination is adjusted by the size of the iris 36 placed immediately behind the color filter 34. A 4f system images the object onto the camera 44 and a through-focus intensity stack 46 is taken by moving the sample (at object plane 48) along the optical axis (or moving the focal plane 42 of the camera or objective, e.g. with an x-y stage (not shown)).

[0039] Image processing of the stack 46 of intensity images for phase 68 and amplitude 66 retrieval is performed via computer 12 and application software 18. Application software 18 comprises instructions, including algorithms of the present invention (such as method 50 of FIG. 2), stored in memory 16 and executable on processor 14 for receiving intensity image stack images 46 and transforming the data to an output image 20.

[0040] While FIG. 1 is detailed with reference to an optical microscope 30, it is appreciated that the application software may be implemented on one or more of: electron or tabletop X-ray imaging systems (medical), synchotron imaging systems, and lithography aerial imaging systems, and phase imaging in both wafer and mask inspection and metrology, among other applications.

2. Recursively Updating the Estimate by Sparse Kalman Filter

[0041] The optical complex field A(x, y, z) is discretized as a raster-scanned complex column vector a.sub.n, and the intensity measurement I(x, y, z) is similarly discretized as column vector I.sub.n. The 2-D Fourier transform of A(x, y, z) is denoted by b(u, v, z). The column vector b.sub.n is again raster-scanned from b(u, v, z), and hence can be expressed as b.sub.n=Ka.sub.n, where K is the discrete Fourier transform (DFT) matrix.

[0042] The relation between two images with distance z in Fourier domain can be written as Eq. 3:

b.sub.n=Hb.sub.n-1.Eq. 3

[0043] Supposing the measurement of the intensity of the camera is modeled with Poisson noise, the Poisson noise with Gaussian distribution can be modeled as:

I.sub.n=C.sub.n|a.sub.n|.sup.2+v,Eq. 4

where v is a Gaussian vector with zero mean and covariance R=diag(C.sub.n|a.sub.n|.sup.2). The matrix C.sub.n describes the convolution of the scaled source intensity distribution S(x, y) at z.sub.n. Here a.sub.n* denotes the complex conjugate of a.sub.n, and diag(a.sub.n*) is a diagonal matrix with its corresponding diagonal entries equal to the elements in the vector a.sub.n*.

[0044] The nonlinear observation model in Eq. 4 is linearized as:

I.sub.n=C.sub.n|.sub.n|.sup.2+C.sub.ndiag(.sub.n*)(a.sub.n.sub.n)+C.sub.ndiag(.sub.n)(.sub.n.sub.n*)+v,Eq. 5

where .sub.n is the state predicted from the previous n1 observations, and Eq. 5 is the first order Taylor series expansion of Eq. 4 with respect to .sub.n.

[0045] Summarizing, the augmented state space model is given as:

[00003]$\begin{array}{cc}.Math.\mathrm{State}.Math.\text{:}.Math.\left[\begin{array}{c}{b}_n\\ b_n^{*}\end{array}\right]=\left[\begin{array}{cc}H& 0\\ 0& H^{*}\end{array}\right]\left[\begin{array}{c}{b}_{n-1}\\ b_{n-1}^{*}\end{array}\right],& \mathrm{Eq}..Math.6\\ \mathrm{Observation}.Math.\text{:}.Math..Math.I_n=\left[\begin{array}{cc}{J}_n& J_n^{*}\end{array}\right]\left[\begin{array}{c}{b}_n\\ b_n^{*}\end{array}\right]-.Math..Math.C_n.Math.{.Math.{\hat{a}}_n.Math.}^2+v,\mathrm{with}.Math..Math.v~\left(0,R\right),& \mathrm{Eq}..Math.7\end{array}$

where

R=diag(C.sub.n|a.sub.n|.sup.2), and J.sub.n=C.sub.ndiag(.sub.n*)K.sup.H.Eq. 8

[0046] FIG. 2 shows a high-level flow diagram for a method 50 of obtaining the state-space Kalman filter in accordance with the present description. First, the process starts with the experimental setup as step 52 (see setup shown in FIG. 1). Next, the partial coherence is modeled mathematically at step 54 and a state-space Kalman filter is formulated from the physical model at step 56. After the state-space Kalman filter, the standard Kalman filter is used to estimate the phase at step 58.

[0047] However, the computational complexity of this estimation is (N.sup.3) where N is the total number of the pixels in the image. A sparse model may be adopted at step 56 to reduce the computational complexity to (N log N). The resulting sparse algorithm is summarized in Table 2. The diagonal entries of diagonal matrix S in Table 2 are equal to Fourier transform of the scaled source

[00004]$S\left(-\frac{f}{z}.Math.x,-\frac{f}{z}.Math.y\right).$

In order to compensate the error of the first order approximation, a forward and backward sweep of the Kalman filter improves the recovered result iteratively. The computational complexity increases linearly with number of iterations.

3. Nonlinear Least Square Error Method

[0048] It is appreciated that the Kalman filter algorithm detailed above is only one of several methods that may be used to solve phase from intensity images of partially coherent illumination. The Kalman filter method processes the intensity images in a recursive manner. Basically, the optimization problem could be formulated as nonlinear least square problem, minimizing a cost function over the estimation:

[00005]$\begin{array}{cc}\underset{b}{\min }.Math.F=\frac{1}{2}.Math.\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.{\left(I_n-C_n.Math.{.Math.K^H.Math.H_n.Math.b.Math.}^2\right)}^H.Math.\left(I_n-C_n.Math.{.Math.K^H.Math.H_n.Math.b.Math.}^2\right),& \mathrm{Eq}..Math.9\end{array}$

where b is a vector denoting Fourier transform of the unknown complex field at z=0, and H.sub.n is a diagonal matrix denote wave propagation for the plane of the unknown complex field to the plane of the measurement I.sub.n, and N.sub.im is the total number of the intensity images. Minimizing the cost function (the error between measurement and prediction) gives an optimal solution of b.

[0049] The nonlinear least square problem may be solved with a Gradient descent method, Newton's method, Gauss-Newton method, or similar methods. These methods involve calculating the gradient and Hessian (approximation of Hessian for Gauss-Newton method) of the cost function. Since the b is complex variable, we use an augmented complex variable [b b].sup.T so that the theory of CR calculus can be used to calculate the gradient and Hessian derivation. Gradient, F(b, b), and approximation of Hessian F(b, b) are derived as follows:

[00006]$\begin{array}{cc}{F}^{}\left(b,\stackrel{\_}{b}\right)=\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.\left[\begin{array}{c}-{\stackrel{\_}{H}}_n.Math.K_{\mathrm{diag}}\left(K^H.Math.H_n.Math.b\right).Math.C_n^H\\ -H_n.Math.K_{\mathrm{diag}}^H\left(K.Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\end{array}\right]\left[I_n-C_n.Math.{.Math.K^H.Math.H_n.Math.b.Math.}^2\right],& \mathrm{Eq}..Math.10\end{array}$

[00007]$\begin{array}{cc}{F}^{}\left(b,\stackrel{\_}{b}\right)\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.\left[\begin{array}{cc}\begin{array}{c}{\stackrel{\_}{H}}_n.Math.K.Math..Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K.Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.K^H.Math.H_n\end{array}& \begin{array}{c}{\stackrel{\_}{H}}_n.Math.K.Math..Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.K.Math..Math.{\stackrel{\_}{H}}_n\end{array}\\ \begin{array}{c}{H}_n.Math.K^H.Math..Math.\mathrm{diag}\left(K.Math..Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K.Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.K^H.Math.H_n\end{array}& \begin{array}{c}{H}_n.Math.K^H.Math..Math.\mathrm{diag}\left(K.Math..Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.K.Math..Math.{\stackrel{\_}{H}}_n\end{array}\end{array}\right].& \mathrm{Eq}..Math.11\end{array}$

[0050] With the Gradient, one may use a Gradient descent method to solve for the nonlinear optimization problem. However, gradient descent method is slow in convergence. Improvement of convergence usually needs the second order derivative of the cost functionHessian or approximation of Hessian. Newton's method, Gauss-Newton method and Levenberg-Marquardt algorithm use the second order derivative information to obtain a search direction for the estimation update, but require solving a big linear equations system. The linear equations system is related to Gradient and Hessian, expressed as:

[00008]$\begin{array}{cc}{F}^{}\left(b,\stackrel{\_}{b}\right)\left[\begin{array}{c}x\\ \stackrel{\_}{x}\end{array}\right]=-F^{}\left(b,\stackrel{\_}{b}\right),& \mathrm{Eq}..Math.12\end{array}$

where the solution of the linear equations [x x].sup.T is the search direction. Since the Hessian matrix of our problem (size N.sup.2) is large, it is difficult to compute and invert the big matrix or solve linear equations of the Hessian matrix. An alternative choice is use L-BFGS method which approximates the Hessian and inversion of Hessian from the gradient of the previous estimations. The drawback of L-BFGS method is that the convergence is highly reliant on the number of estimations used to approximate the Hessian.

[0051] The present disclosure incorporates a new method to solve the linear equations system iteratively. The matrix F(b, b) includes DFT matrix K and circulant matrix C.sub.n, which denotes convolution. So the multiplication of the matrix F(b, b) with a vector can be computed with FFT, which has the computational complexity to (N log N). The efficiency of calculating the multiplication could be exploited in solving the linear equations with conjugate gradient method, because the conjugate gradient needs to compute the multiplication during the iteration which acts as the bottleneck of the computational complexity. Not only is the computation sped up by adapting FFT, but also the storage requirement is alleviated because we don't have to compute the matrix F(b, b) explicitly. It is also noted that the gradient F(b, b) could also be computed with FFT. And hence, the search direction is solved with a conjugate method which is fast and efficient. After the search direction is obtained, one can use a line search method to find a proper step size to update the estimation of current step with the search direction. The computational time can be further reduced by parallelizing the computing in the gradient calculation and the multiplication of Hessian matrix with a vector.

4. Results

[0052] The experimental setup used for demonstration is shown FIG. 1. An incoherent white light source 32, filtered by a narrow-band color filter 34 with center wavelength 650 nm, is placed at the front focal plane of the condenser 40. The spatial coherence of the illumination is adjusted by the size of an iris placed immediately behind the color filter. A 4f system 30 images the object onto the camera 44 and a through focus intensity stack 46 is taken by moving the sample along the optical axis (or moving the camera or objective).

[0053] The efficacy of the system and method of the present description are first demonstrated by simulating a phase and amplitude object illuminated by a circular incoherent source 32 (2 mm in diameter, and the focal length of the condenser 40 is 10 mm). An intensity stack 46, partially shown in FIG. 3A, comprised of 101 images that were generated by defocusing the object symmetrically about the focus with 10 mm step size. The data was further corrupted by white Gaussian noise with a variance of 0.0015.

[0054] To demonstrate the effect of partial coherence, the data was first processed using a fully coherent model, with results shown the output image 60 in FIG. 3B. Significant blur present in the coherent phase reconstruction result is caused by the intensity smearing due to off-center source points, which cause the illumination to be less coherent and the model to break down.

[0055] The phase reconstruction by incorporating partially coherent illumination, on the other hand, successfully eliminates the blurring artifact, as shown in the output image 62 FIG. 3C.

[0056] Next, the system and methods of the present description were evaluated experimentally with a biological sample and a microscope. A cheek cell sample was defocused symmetrically about the focus at 81 z-planes ranging from 2.5 mm to 2.5 mm, shown as images 70 in FIG. 4A. Each image contains 701625 pixels. The circular iris had a 2 mm diameter and the light is collimated by a 100 mm focal length condenser. FIG. 4B shows image phase 72 recovered by sparse ACEKF assuming fully coherent illumination. FIG. 4C shows an image phase 74 recovered by Kalman filtering with partial coherence.

[0057] As seen in FIG. 4B and FIG. 4C, the phase recovered by the partially coherent Kalman filtering method of the present description has higher contrast and shows much more detail inside the cell, giving a much improved result.

[0058] Referring now to FIG. 5, the method of the present description was further evaluated with different levels of partial coherence (different source shapes) on the microscope. The temporal coherence is not considered here since the white light is filtered by a color filter of 550 nm. Intensity focal stacks were taken with different source shapes to vary the partial coherence of the illumination. For the first three shapes, the sizes were varied by opening the condenser aperture. The ring shape in the fourth case was obtained by choosing the annular diagram of the Phi condenser. Note that the objective lens is not a Phi type, so the images are still brightfield mode. The images of the source shapes were taken by adding a 2f system at the side port of the microscope. For each shape, the same cheek cell sample was defocused symmetrically about the focus at 131 z-planes ranging from 32.5 m to 32.5 m. Each image contains 800 by 800 pixels.

[0059] FIG. 5 compares the recovered phase images using the coherent Kalman filter and the partially coherent Kalman filter of the present disclosure. As the coherence decreases (aperture opens up), the phase images from the coherent Kalman filter become more blurred. However, the recovered phase images from the partially coherent Kalman filter still have high contrast and contain fine details.

5. Alternative Configurations

[0060] The above experimental demonstrations involve sources of increasing radius circles and intensity images with propagation as the transfer function. In fact, the method of the present disclosure is more general and can be used with any size or shape of source (dipole, quadropole, etc.), simply by modifying the shape in the convolution of Eq. 2. Further, the transfer function between consecutive images, H, can be of any complex form, without changing the algorithm or method. Propagation is just one example, as it is the most common scheme, but one could instead place filters in the imaging pathway successively to collect various images with different types of phase contrast.

[0061] The nonlinear least square error method is derived to recover phase from partially coherent intensity images. However, the algorithm could be adapted into similar situations. For example, the cost function takes form:

[00009]$\begin{array}{cc}\underset{b}{\min }.Math.F=\frac{1}{2}.Math.\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.{\left(I_n-C_n.Math.{.Math.A_n.Math.b.Math.}^2\right)}^H.Math.\left(I_n-C_n.Math.{.Math.A_n.Math.b.Math.}^2\right)& \mathrm{Eq}..Math.13\end{array}$

where b is the unknown, and the matrix A denotes the transfer function. The Hessian (or its approximation) and gradient could take the same form:

[00010]$\begin{array}{cc}{F}^{}\left(b,\stackrel{\_}{b}\right)=\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.\left[\begin{array}{c}-A_n^H.Math.\mathrm{diag}\left(A_n.Math.b\right).Math.C_n^H\\ -A_n^T.Math.\mathrm{diag}\left({\stackrel{\_}{A}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\end{array}\right]\left[I_n-C_n.Math.{.Math.A_n.Math.b.Math.}^2\right],& \mathrm{Eq}..Math.14\\ F^{}\left(b,\stackrel{\_}{b}\right)\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.\left[\begin{array}{cc}\begin{array}{c}{A}_n^H.Math.\mathrm{diag}\left(A_n.Math.b\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left({\stackrel{\_}{A}}_n.Math.\stackrel{\_}{b}\right).Math.A_n\end{array}& \begin{array}{c}{A}_n^H.Math.\mathrm{diag}\left(A_n.Math.b\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(A_n.Math.b\right).Math.{\stackrel{\_}{A}}_n\end{array}\\ \begin{array}{c}{A}_n^T.Math.\mathrm{diag}\left({\stackrel{\_}{A}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left({\stackrel{\_}{A}}_n.Math.\stackrel{\_}{b}\right).Math.A_n\end{array}& \begin{array}{c}{A}_n^T.Math.\mathrm{diag}\left({\stackrel{\_}{A}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(A_n.Math.b\right).Math.{\stackrel{\_}{A}}_n\end{array}\end{array}\right].& \mathrm{Eq}..Math.15\end{array}$

[0062] The same idea could be easily adapted to solve the big linear equations related to F(b, b) when the multiplication of F(b, b) with a vector can be computed efficiently, for example with FFT. For example, in case phase recovery from coherent intensity, the matrices C.sub.n are identity matrix and the matrix A.sub.n=K.sup.HH.sub.n. And hence the multiplication of F(b, b) with a vector could be computed with FFT. In the case of ptychography or Fourier ptychography, the multiplication of F(b, b) with a vector could also be computed efficiently. Specifically, it should be noted that in some cases the Hessian, instead of its approximation, could be derived explicitly, the error in the cost function could be weighted, and an extra regularization over the estimation could be added to the cost function. In all of these situations, the same idea could be used as long as the multiplication of a big matrix with a vector could be computed efficiently.

[0063] Another possible variation for performing partially coherent phase recovery is that extra steps may be taken to estimate the sources. The source shape's size or form is controlled by several limited number of parameters. The source estimation step could then aim to estimate these parameters. The source estimation step could also recover arbitrary source shape. Therefore, the phase can possibly be recovered from the measured intensity images without providing information of the source shape, and source shape could be recovered simultaneously.

6. Conclusion

[0064] Table 1 compares the phase recovery method of the present disclosure with phase contrast microscopy and digital holography. As seen in Table 1, significant unexpected improvements are gained from the systems and methods of the present disclosure.

[0065] Previous coherent Kalman filtering methods assume a coherent illumination in their models. The method of the present disclosure extends the model to the important case of partially coherent illumination, and recovers phase of an object when it is illuminated with partially coherent illumination. The result is a valid construct for commercial optical microscopes, electron, tabletop X-ray imaging systems (medical), synchrotron imaging systems, and lithography aerial imaging systems, among other applications.

[0066] The system and method of the present disclosure is capable of phase recovery for arbitrary shapes of Khler illumination. In lithography, dipole or multiple illumination is widely used and phase is a well-known diagnostic test for mask and wafer defects. This new method will be crucial to phase imaging in both wafer and mask inspection and metrology.

[0067] The system and method of the present disclosure can be applied to the case of strong noise, such as low photon counts caused by high-throughput imaging. The Kalman filtering technique is robust to noise because the Kalman filter gives a near-optimal estimation from a series of noisy measurements.

[0068] The system and method of the present disclosure is valid for any variety of complex transfer functions between the captured intensity images (propagation (defocus) was used for exemplary purposes only as a simple example). Alternatively, phase or amplitude filters or other linear optical systems can replace the propagation transfer function with their appropriate transfer function to achieve a phase result with the same method.

[0069] The system and method of the present disclosure may further be refined to be computationally efficient by using a convolution model. It can be programmed to run on a Graphics Processing Unit (GPU) to achieve greater speed-up in processing times.

[0070] Embodiments of the present technology may be described with reference to flowchart illustrations of methods and systems according to embodiments of the technology, and/or algorithms, formulae, or other computational depictions, which may also be implemented as computer program products. In this regard, each block or step of a flowchart, and combinations of blocks (and/or steps) in a flowchart, algorithm, formula, or computational depiction can be implemented by various means, such as hardware, firmware, and/or software including one or more computer program instructions embodied in computer-readable program code logic. As will be appreciated, any such computer program instructions may be loaded onto a computer, including without limitation a general purpose computer or special purpose computer, or other programmable processing apparatus to produce a machine, such that the computer program instructions which execute on the computer or other programmable processing apparatus create means for implementing the functions specified in the block(s) of the flowchart(s).

[0071] Accordingly, blocks of the flowcharts, algorithms, formulae, or computational depictions support combinations of means for performing the specified functions, combinations of steps for performing the specified functions, and computer program instructions, such as embodied in computer-readable program code logic means, for performing the specified functions. It will also be understood that each block of the flowchart illustrations, algorithms, formulae, or computational depictions and combinations thereof described herein, can be implemented by special purpose hardware-based computer systems which perform the specified functions or steps, or combinations of special purpose hardware and computer-readable program code logic means.

[0072] Furthermore, these computer program instructions, such as embodied in computer-readable program code logic, may also be stored in a computer-readable memory that can direct a computer or other programmable processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the block(s) of the flowchart(s). The computer program instructions may also be loaded onto a computer or other programmable processing apparatus to cause a series of operational steps to be performed on the computer or other programmable processing apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable processing apparatus provide steps for implementing the functions specified in the block(s) of the flowchart(s), algorithm(s), formula (e), or computational depiction(s).

[0073] It will further be appreciated that the terms programming or program executable as used herein refer to one or more instructions that can be executed by a processor to perform a function as described herein. The instructions can be embodied in software, in firmware, or in a combination of software and firmware. The instructions can be stored local to the device in non-transitory media, or can be stored remotely such as on a server, or all or a portion of the instructions can be stored locally and remotely. Instructions stored remotely can be downloaded (pushed) to the device by user initiation, or automatically based on one or more factors. It will further be appreciated that as used herein, that the terms processor, computer processor, central processing unit (CPU), and computer are used synonymously to denote a device capable of executing the instructions and communicating with input/output interfaces and/or peripheral devices.

[0074] From the description herein, it will be appreciated that that the present disclosure encompasses multiple embodiments which include, but are not limited to, the following:

[0075] 1. A system for phase recovery of an intensity image with partially coherent illumination, the system comprising: (a) a computer processor; and (b) a memory storing instructions executable on the processor; (c) the instructions, when executed, performing the steps comprising: (i) mathematically modeling partial coherence with respect to the intensity image; (ii) formulating a state-space Kalman filter from the partial coherence model; and (iii) applying a standard Kalman filter to the state-space Kalman filter to estimate a phase with respect to the intensity image.

[0076] 2. The system of any preceding embodiment, wherein the image comprises a stack of defocused intensity images; and wherein a propagation transfer function is applied between captured intensity images

[0077] 3. The system of any preceding embodiment, further comprising: applying a sparse model to the Kalman filter to reduce computational complexity.

[0078] 4. The system of any preceding embodiment, wherein the sparse model optimizes according to:

.sub.n=K.sup.H{circumflex over (b)}.sub.n

Q.sub.n={circumflex over (Q)}.sub.n({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*+({circumflex over (P)}.sub.n)*+qS.sup.2).sup.1({circumflex over (Q)}.sub.n+({circumflex over (P)}.sub.n)*)

P.sub.n={circumflex over (P)}.sub.n({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*+({circumflex over (P)}.sub.n)*+qS.sup.2).sup.1({circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*)

b.sub.n={circumflex over (b)}.sub.n+(Q.sub.n+P.sub.n)Kdiag(.sub.n)C.sub.nR.sup.1(I.sub.n|.sub.n|.sup.2);

where .sub.n is the state predicted from the previous n1 observations, I.sub.n and b.sub.n are discretized column vectors, C.sub.n describes convolution of scaled source intensity distribution, and K is the discrete Fourier transform (DFT) matrix.

[0079] 5. The system of any preceding embodiment, wherein the instructions are configured for estimating a phase with respect to a partially coherent intensity image generated from one or more of: optical microscopes, electron imaging systems, tabletop X-ray imaging systems, synchrotron imaging systems, and lithography aerial imaging systems.

[0080] 6. The system of any preceding embodiment, wherein the instructions are configured for estimating phase for an image having an arbitrary shape of Khler illumination.

[0081] 7. The system of any preceding embodiment, wherein the image comprises a stack of intensity images; and wherein phase or amplitude filters or other linear optical systems are applied as a transfer function between captured intensity images.

[0082] 8. A method for phase recovery of an intensity image with partially coherent illumination, the method comprising: (a) mathematically modeling partial coherence with respect to the intensity image; (b) formulating a state-space Kalman filter from the partial coherence model; and (c) applying a standard Kalman filter to the state-space Kalman filter to estimate a phase with respect to the intensity image; (d) wherein said method is performed by executing programming on at least one computer processor, said programming residing on a non-transitory medium readable by the computer processor.

[0083] 9. The method of any preceding embodiment, wherein the image comprises a stack of defocused intensity images; and wherein a propagation transfer function is applied between captured intensity images.

[0084] 10. The method of any preceding embodiment, further comprising: applying a sparse model to the Kalman filter to reduce computational complexity.

[0085] 11. The method of any preceding embodiment, wherein the sparse model optimizes according to:

.sub.n=K.sup.H{circumflex over (b)}.sub.n

Q.sub.n={circumflex over (Q)}.sub.n({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*+({circumflex over (P)}.sub.n)*+qS.sup.2).sup.1({circumflex over (Q)}.sub.n+({circumflex over (P)}.sub.n)*)

P.sub.n={circumflex over (P)}.sub.n({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n+{circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*+({circumflex over (P)}.sub.n)*+qS.sup.2).sup.1({circumflex over (P)}.sub.n+({circumflex over (Q)}.sub.n)*)

b.sub.n={circumflex over (b)}.sub.n+(Q.sub.n+P.sub.n)Kdiag(.sub.n)C.sub.nR.sup.1(I.sub.n|.sub.n|.sup.2);

where .sub.n is the state predicted from the previous n1 observations, I.sub.n and b.sub.n are discretized column vectors, C.sub.n describes convolution of scaled source intensity distribution, and K is the discrete Fourier transform (DFT) matrix.

[0086] 12. The method of any preceding embodiment, wherein the method is configured for estimating a phase with respect to a partially coherent intensity image generated from one or more of: optical microscopes, electron imaging systems, tabletop X-ray imaging systems, synchrotron imaging systems, and lithography aerial imaging systems.

[0087] 13. The method of any preceding embodiment, wherein the method is configured for estimating phase for an image having an arbitrary shape of Khler illumination.

[0088] 14. The system of any preceding embodiment, wherein the image comprises a stack of intensity images; and wherein phase or amplitude filters or other linear optical systems are applied as a transfer function between captured intensity images.

[0089] 15. A system for phase recovery of an intensity image with partially coherent illumination, the system comprising: (a) a computer processor; and (b) a memory storing instructions executable on the processor; (c) the instructions, when executed, performing the steps comprising: (i) mathematically modeling partial coherence with respect to one or more intensity image having partially coherent illumination; (ii) applying a nonlinear least square error function to fit over the one or more measured intensity images; and (iii) estimating a phase with respect to the one or more intensity images.

[0090] 16. The system of any preceding embodiment, wherein applying a nonlinear least square error function comprises finding an optimal estimation to minimize an average square error between measurement and predicted measurement (cost function) with respect to the one or more intensity images.

[0091] 17. The system of any preceding embodiment, wherein applying a nonlinear least square error function comprises: minimizing the cost function using nonlinear least square error method; and deriving a gradient and Hessian.

[0092] 18. The system of any preceding embodiment, further comprising: computing a multiplication of the Hessian with vector with Fast Fourier transform (FFT).

[0093] 19. The system of any preceding embodiment, wherein one or more of the following algorithms are used to solve for said phase: gradient descent, L-BFGS, Newton's, Gauss-Newton, Levenberg-Marquardt, or other nonlinear conjugate gradient algorithms.

[0094] 20. The system of any preceding embodiment, wherein the gradient and Hessian are computed according to:

[00011]$\begin{array}{c}{F}^{}\left(b,\stackrel{\_}{b}\right)=\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.\left[\begin{array}{c}-A_n^H.Math.\mathrm{diag}\left(A_n.Math.b\right).Math.C_n^H\\ -A_n^T.Math.\mathrm{diag}\left({\stackrel{\_}{A}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\end{array}\right]\left[I_n-C_n.Math.{.Math.A_n.Math.b.Math.}^2\right];\\ F^{}\left(b,\stackrel{\_}{b}\right)\underset{n=1}{\overset{N_{\mathrm{im}}}{.Math.}}.Math.\left[\begin{array}{cc}\begin{array}{c}{\stackrel{\_}{H}}_n.Math.K.Math..Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K.Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.K^H.Math.H_n\end{array}& \begin{array}{c}{\stackrel{\_}{H}}_n.Math.K.Math..Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.K.Math..Math.{\stackrel{\_}{H}}_n\end{array}\\ \begin{array}{c}{H}_n.Math.K^H.Math..Math.\mathrm{diag}\left(K.Math..Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K.Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.K^H.Math.H_n\end{array}& \begin{array}{c}{H}_n.Math.K^H.Math..Math.\mathrm{diag}\left(K.Math..Math.{\stackrel{\_}{H}}_n.Math.\stackrel{\_}{b}\right).Math.C_n^H\\ C_n.Math.\mathrm{diag}\left(K^H.Math.H_n.Math.b\right).Math.K.Math..Math.{\stackrel{\_}{H}}_n\end{array}\end{array}\right].\end{array}$

[0095] 21. The system of any preceding embodiment, further comprising solving the linear equations system:

[00012]$F^{}\left(b,\stackrel{\_}{b}\right)\left[\begin{array}{c}x\\ \stackrel{\_}{x}\end{array}\right]=-F^{}\left(b,\stackrel{\_}{b}\right).$

[0096] 22. The system of any preceding embodiment, wherein the linear equation is solved as a function of a conjugate gradient method iteratively by computing the gradient and the multiplication of the Hessian with a vector with FFT.

[0097] Although the description herein contains many details, these should not be construed as limiting the scope of the disclosure but as merely providing illustrations of some of the presently preferred embodiments. Therefore, it will be appreciated that the scope of the disclosure fully encompasses other embodiments which may become obvious to those skilled in the art.

[0098] In the claims, reference to an element in the singular is not intended to mean one and only one unless explicitly so stated, but rather one or more. All structural, chemical, and functional equivalents to the elements of the disclosed embodiments that are known to those of ordinary skill in the art are expressly incorporated herein by reference and are intended to be encompassed by the present claims. Furthermore, no element, component, or method step in the present disclosure is intended to be dedicated to the public regardless of whether the element, component, or method step is explicitly recited in the claims. No claim element herein is to be construed as a means plus function element unless the element is expressly recited using the phrase means for. No claim element herein is to be construed as a step plus function element unless the element is expressly recited using the phrase step for.

TABLE-US-00001 TABLE 1 Advantages Over Phase Contrast Microscopy And Digital Holography Partially Phase Contrast Digital Coherent Phase Property Microscopy Holography Recovery Digital processing Difficult Yes Yes Phase wrapping No Yes No problem Quantitative No Yes Yes Noise resilient No No Yes Experimental setup Complicated Complicated Simple Real-time imaging Yes Potentially Yes

TABLE-US-00002 TABLE 2 Sparse Augmented Complex Extended Kalman Filter For Estimating A Wave Field 1. Initialization b.sub.0 , Q.sub.0 and P.sub.0 . 2. Prediction: {circumflex over (b)}.sub.n = Hb.sub.n1;{circumflex over (Q)}.sub.n = Q.sub.n1; {circumflex over (P)}.sub.n = HP.sub.n1H. 3. Update: .sub.n = K.sup.H{circumflex over (b)}.sub.n Q.sub.n = {circumflex over (Q)}.sub.n ({circumflex over (Q)}.sub.n + {circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n + {circumflex over (P)}.sub.n + ({circumflex over (Q)}.sub.n)* + ({circumflex over (P)}.sub.n)* + qS.sup.2).sup.1 ({circumflex over (Q)}.sub.n + ({circumflex over (P)}.sub.n)*) P.sub.n = {circumflex over (P)}.sub.n ({circumflex over (Q)}.sub.n + {circumflex over (P)}.sub.n)({circumflex over (Q)}.sub.n + {circumflex over (P)}.sub.n + ({circumflex over (Q)}.sub.n)* + ({circumflex over (P)}.sub.n)* + qS.sup.2).sup.1 ({circumflex over (P)}.sub.n + ({circumflex over (Q)}.sub.n)*) b.sub.n = {circumflex over (b)}.sub.n + (Q.sub.n + P.sub.n)Kdiag(.sub.n)C.sub.nR.sup.1 (I.sub.n |.sub.n |.sup.2).