Method for estimating the axial intensity derivative in the TIE with higher order intensity derivatives and noise suppression

Rui Bie; Xiu-Hua Yuan; Ming Zhao; Li Zhang

doi:10.1364/OE.20.008186

1. Introduction

Phase retrieval from intensity measurements is a noninterferometric technique suggested originally by Teague in 1983 for wavefront sensing [1]. Knowledge of the phase is helpful for correction of distortion produced by the turbulent atmosphere in some related fields such as free space optical communication (FSO) and astronomic observation.

Teague introduced the transport of intensity equation (TIE) which relate the object-plane to the axial derivative of the intensity distribution, $\partial I (x, y, z) / \partial z$ , in the Fresnel region to us, and if the intensity distribution has no zeros, the TIE has a unique solution (with an additive arbitrary constant). Many effective methods such as Fast Fourier Transform (FFT) [2, 3], the multigrid (MG) [4], the Zernike polynomial expansion methods [5, 6], and the Green’s function based method [1] have been proposed for solving it. And all these methods require two images that are defocused from each other in order for the axial derivation of intensity derivative. However, there are tradeoffs among the amount of defocus, the accuracy of the result and noise considerations.

In recent years, some algorithms for derivation of intensity derivative from intensity measurements in multiple planes are proposed: In 2007, Marcos Soto [7] evaluated the variance of the detecting noise and the axial second-order intensity derivative with multi-plane intensity distribution, based on which a method minimizing the error is proposed; and in 2010 Laura Waller [8] of MIT took the impact of higher order derivative into account and described it with a matrix transformation; a year later, Bindang Xue [9, 10] of Beihang University in China extended Laura Waller’s work to the situation with intensity distribution measured in unequally-spaced planes.

This paper proposed a novel method to evaluate intensity derivative, which takes both the detecting noise and higher order derivative into consideration. Based on Taylor expansions of the measured intensities, the axial intensity derivative in the TIE is expressed by a linear combination of the multiple intensity measurements, which is used to estimate the difference from the actual value. Numerical simulations are conducted to test it.

2. Multiple intensity measurements for phase retrieval

An optical plane wave $u (r)$ on the x-y plane can be described by optical intensity and phase as

u (r) = {[I (r)]}^{1 / 2} \exp [i ϕ (r)],

where

r = (x, y)

denotes the coordinate in the x-y plane, perpendicular to the z axis,

I (r)

and

ϕ (r)

are the intensity and phase distribution in the plane respectively. In this paper, the focal plane is assumed to be located at

z = 0

.

Under the paraxial approximation, the transport of intensity equation (TIE) can be derived as

- k \frac{\partial I}{\partial z} = \nabla \cdot (I \nabla ϕ),

where

k = 2 π / λ

is the wave number, and

λ

denotes the wavelength,

\nabla = {\frac{\partial}{\partial x}, \frac{\partial}{\partial y}}

is the two dimensional gradient operator in the x-y plane.

Assume the intensity distribution of m x-y planes located at $z_{1}, z_{2}, \dots z_{m}$ is available, based on the Maclaurin formula, $I (z_{i})$ (ignore the x-y dimension) can be expanded as

I (z_{i}) = I (0) + \frac{\partial I (0)}{\partial z} \cdot z_{i} + \sum_{k = 2}^{n} \frac{z_{i}^{k}}{k!} \frac{\partial^{k} I (0)}{\partial z^{k}} + \frac{z_{i}^{n + 1}}{(n + 1)!} \frac{\partial^{n + 1} I (ξ_{i})}{\partial z^{n + 1}},

where

\frac{\partial^{k}}{\partial z^{k}}

denotes the k^th order derivative, and

ξ_{i}

is some number located between 0 and

z_{i}

. As a result, we get m equations. Multiply these equations with a number

a_{i}

separately and sum their left and right hand side respectively, we can get the following expression

\sum_{i = 1}^{m} a_{i} I (z_{i}) = I (0) \cdot \sum_{i = 1}^{m} a_{i} + \frac{\partial I (0)}{\partial z} \cdot \sum_{i = 1}^{m} a_{i} z_{i} + \sum_{k = 2}^{n} \frac{1}{k!} \frac{\partial^{k} I (0)}{\partial z^{k}} \sum_{i = 1}^{m} a_{i} z_{i}^{k} + \frac{1}{(n + 1)!} \sum_{i = 1}^{m} a_{i} z_{i}^{n + 1} \frac{\partial^{n + 1} I (ξ_{i})}{\partial z^{n + 1}} .

Let

\begin{matrix} \sum_{i = 1}^{m} a_{i} = 0, & \sum_{i = 1}^{m} a_{i} z_{i} = 1, & \begin{matrix} \sum_{i = 1}^{m} a_{i} z_{i}^{k} = 0 & k = 2, 3... n \end{matrix} \end{matrix},

we can derive the derivative of intensity

\frac{\partial I (0)}{\partial z} = \sum_{i = 1}^{m} a_{i} I (z_{i}) - \frac{1}{(n + 1)!} \sum_{i = 1}^{m} a_{i} z_{i}^{n + 1} \frac{\partial^{n + 1} I (ξ_{i})}{\partial z^{n + 1}} .

Ignoring the second term on the right hand of Eq. (6),the intensity derivative can be calculated in practical system like

{[\partial_{z} I]}_{c} = \sum_{i = 1}^{m} a_{i} I (z_{i}) .

If the detecting noise is taken into consideration, the difference from Eq. (6) is

{[\partial_{z} I]}_{c} - \frac{\partial I (0)}{\partial z} = \sum_{i = 1}^{m} a_{i} N (z_{i}) + \frac{1}{(n + 1)!} \sum_{i = 1}^{m} a_{i} z_{i}^{n + 1} \frac{\partial^{n + 1} I (ξ_{i})}{\partial z^{n + 1}},

where

N (z_{i})

denotes the intensity noise at

z_{i}

.Under the assumption of noise statistic independence, the mean square can be written as

ε^{2} = D (N) + D (ξ),

where

D (N) = σ^{2} \sum_{i = 1}^{m} a_{i}^{2},

D (ξ) = {[\frac{1}{(n + 1)!} \sum_{i = 1}^{m} a_{i} z_{i}^{n + 1} \frac{\partial^{n + 1} I (ξ_{i})}{\partial z^{n + 1}}]}^{2},

and

σ^{2}

denotes the variance of the noise. Because of the uncertainty of high-order item, we only evaluate the impact of detecting noise i.e. Equation (10). Our mission is to minimize it under the condition of Eq. (5).

Generally, Eq. (5) has another description form

[\begin{matrix} 1 & 1 & 1 & \dots & 1 \\ z_{1} & z_{2} & z_{3} & \dots & z_{m} \\ z_{1}^{2} & z_{2}^{2} & z_{3}^{2} & \dots & z_{m}^{2} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ z_{1}^{n} & z_{2}^{n} & z_{3}^{n} & \dots & z_{m}^{n} \end{matrix}] \cdot [\begin{matrix} a_{1} \\ a_{2} \\ a_{3} \\ ⋮ \\ a_{m} \end{matrix}] = [\begin{matrix} 0 \\ 1 \\ 0 \\ ⋮ \\ 0 \end{matrix}] .

Obviously, the first matrix of left hand side called transformation matrix in this paper is a Vandermonde matrix.

(1) When $n + 1 > m$ , Eq. (12) are an overdetermined equations. Under this condition, it pointless to evaluate the n^th order derivative using intensity data of m planes.
(2) When $n + 1 = m$ , it is the case that Bindang Xue described, the coefficient $a_{1} ~ a_{m}$ can be calculated definitely from Eq. (12), when the detecting noise cannot be suppressed.
(3) When $n + 1 < m$ , Eq. (12) are an underdetermined equations, it is necessary to lay our attention on this situation which has not been considered.

We can divide the matrix of Eq. (12) as

[\begin{matrix} A & B \end{matrix}] \cdot [\begin{matrix} C \\ D \end{matrix}] = E,

where,

B

consists of the last n + 1 columns of the transformation matrix, which is nonsingular. And

C

,

D

is divided to match

A

and

B

.

Then $D$ can be written as

D = B^{- 1} E - B^{- 1} A C .

As a result,

D (N) = σ^{2} [\sum_{i = 1}^{n + 1} a_{i}^{2} + \sum_{j = 1}^{m - n - 1} {(f_{j} - G_{j} C)}^{2}] .

To minimize it,

\begin{matrix} \frac{\partial}{\partial a_{k}} D (N) = σ^{2} [2 a_{k} - \sum_{j = 1}^{n + 1} 2 (f_{j} - G_{j} C) g_{j k}] = 0 & k = 1, 2 \dots (m - n - 1) \end{matrix},

where

f_{j}

denotes the j^th element of the vector

F = B^{- 1} E

, and

G_{j}

is j^th row vector of matrix

G = B^{- 1} A

. Obviously, Eq. (16) is equivalent to

(G^{T} G + I) C = G^{T} F .

Combining with Eq. (13), we get

[\begin{matrix} G^{T} G + I & 0 \\ A & B \end{matrix}] \cdot [\begin{matrix} C \\ D \end{matrix}] = [\begin{matrix} G^{T} F \\ E \end{matrix}] .

In Eq. (18), $G^{T} G + I$ and $B$ are both nonsingular matrix, so the coefficient $a_{1} ~ a_{m}$ can be derived easily.

3. Simulations

We take the “Gauss Beam”, a widely used beam mode, to check the validity of the method described above.

That is the 0th order Gauss mode

u_{00} (x, y, z) = \frac{c}{w (z)} \exp (- \frac{r^{2}}{w^{2} (z)}) \exp {- i [k (z + \frac{r^{2}}{2 R}) - \arctan (\frac{z}{f})]},

where

\begin{matrix} k = \frac{2 π}{λ} & r^{2} = x^{2} + y^{2} & w (z) = w_{0} \sqrt{1 + {(\frac{z}{f})}^{2}} \\ R = R (z) = z [1 + {(\frac{f}{z})}^{2}] & f = \frac{π w_{0}^{2}}{λ} & w_{0} = \sqrt{\frac{λ f}{π}} \end{matrix} .

As a result, the intensity distribution of the axis can be written as

I (z) = \frac{c^{2}}{w^{2} (z)},

where

w_{0} = 0.005 m, λ = 1550 n m, f = \frac{π w_{0}^{2}}{λ} = 50.6708 m, c = 0.01.

The simulation error of Bindang Xue’s method (we call it high order derivative suppression method here) and the algorithm with noise suppression proposed in this paper are both shown in Fig. 1 .

Fig. 1 Axial intensity derivative error computed by noise suppression (described in this paper) and high order derivative suppression method (proposed by Bindang Xue) for different sampling number m: Each point is the average of 100 results. The blue lines show the results derived by noise suppression method for n = 4. (a)~(d) denotes the results for different sampling intervals $Δ z$ .

Download Full Size | PDF

When sampling interval $Δ z$ is large, the error raised by high order derivative impacts the final result significantly. Bindang Xue’s method plays a good role in this situation as shown in Fig. 1(a). But as it decreases, the influence of high order derivative vanishes gradually, while the detecting noise takes the dominated place step by step. As described in Fig. 1(a)–1(d) the performance of noise suppression method improves significantly.

Figure 2 shows the simulation error of noise suppression method for different n when m = 16. When n is small, the capability of noise suppression is strong, but not that of high order derivative suppression, which leads to big error from the actual value. As n increases, these two types of suppression reach balance point step by step. Under this condition, the error caused by the algorithm decreases, and finally reaches the lowest point. But when the growth of n is going on, the dominant position of error is taken by detecting noise gradually, which makes the error rise. From what is discussed above, we can safely get the conclusion that there exist an optimal n for a sampling number m, so it is necessary to search this point in practice.

Fig. 2 Axial intensity derivative error computed by noise suppression method for different n, the highest order derivative considered, when sampling number m = 16

Download Full Size | PDF

4. Conclusions

This paper proposed a novel method to derive the intensity derivative, which extends the work of Marcos Soto etc [7–10]. Including higher order intensity derivatives and detecting noise improves the accuracy of intensity derivative computing, which leads to enhancement of the algorithm performance based on TIE. By adopting correction for nonlinearities (higher order derivatives) and noise suppression, large stacks of data spanning longer distances can be use to estimate the derivative more precisely.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (No.61077058).

References and links

1. M. Reed Teague, “Deterministic phase retrieval: a Green's function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]

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

3. T. E. Gureyev and K. A. Nugent, “Rapid quantitative phase imaging using the transport of intensity equation,” Opt. Commun. 133(1-6), 339–346 (1997). [CrossRef]

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

5. T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation. II. orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A 13(8), 1670–1682 (1996). [CrossRef]

6. T. E. Gureyev, A. Roberts, and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation: matrix solution with use of Zernike polynomials,” J. Opt. Soc. Am. A 12(9), 1932–1942 (1995). [CrossRef]

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

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

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

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

Method for estimating the axial intensity derivative in the TIE with higher order intensity derivatives and noise suppression

Abstract

1. Introduction

2. Multiple intensity measurements for phase retrieval

3. Simulations

4. Conclusions

Acknowledgments

References and links

Cited By

Figures (2)

Equations (22)

Optics Express