Transport of intensity phase imaging from multiple intensities measured in unequally-spaced planes

Bindang Xue; Shiling Zheng; Linyan Cui; Xiangzhi Bai; Fugen Zhou

doi:10.1364/OE.19.020244

1. Introduction

The TIE can relate the object-plane phase to the axial derivative of the intensity distribution in the Fresnel region [1], and if the intensity distribution has no zeros, the TIE has a unique solution (within an additive constant) for coherent and partially coherent fields [2]. Many effective numerical methods have been proposed for solving the TIE. Therefore, the TIE has been widely used for the recovery of optical phase from multiple intensity measurements, and phase recovery has been an essential imaging technique in many fields, e.g., optics [3], electron [4], neutrons [5] and x-ray microscopy [6].

In the numerical methods for solving the TIE, such as Fast Fourier Transform (FFT) [7,8], the multigrid (MG) [9,10], the Zernike polynomial expansion methods [11,12], and the Green’s function based method [13], the axial intensity derivative is commonly approximated by the finite difference method from two intensity measurements. However, two intensity measurements only achieve a maximum of second-order precision with respect to the separation of the adjacent intensity measurements if the intensity distributions do not contain noise [14]. To improve the accuracy of the retrieved phase, multiple intensity measurements can be used to minimize the effect of the noise on the retrieved phase [14] or estimate the higher order axial intensity derivatives in the TIE [15,16]. The previous work mostly assumed that the intensity measurements are taken in the equally-spaced planes. However, the measurement errors exist objectively in practical application. Therefore, the intensity measurements are not really in the equally-spaced planes.

In this study, a generalization of the TIE is proposed, which uses multiple unevenly sampled intensity measurements for phase retrieval. Based on the higher-order Taylor expansions of the measured intensities, the axial intensity derivative in the TIE is expressed by a linear combination of the multiple unevenly sampled intensity measurements. Numerical simulations are conducted to test the method.

2. Unequally-spaced intensity measurements for phase retrieval

An optical plane wave $u_{z} (r)$ in the x-y plane propagating along the $Z$ axis can be expressed as:

u_{z} (r) = I_{z}^{1 / 2} (r) \exp (i ϕ_{z} (r)),

where

r = (x, y)

is the spatial coordinates in the x-y plane, perpendicular to the

Z

axis.

I_{z} (r)

and

ϕ_{z} (r)

are respectively the intensity and the phase distributions. The focal plane is located in the

z = 0

plane. With the paraxial approximation, the transport of intensity equation (TIE) states that [1]:

\nabla \cdot [I_{z} (r) \nabla ϕ_{z} (r)] = - k \partial I_{z} (r) / \partial z,

where _{$k = 2 π / λ$} is the wave number,

λ

is the wavelength, and_{$\nabla = (\partial x, \partial y)$} denotes the two-dimensional gradient operator in the x-y plane.

Assume that intensities measured in $z = {0, z_{1}, z_{2} \dots, z_{M}}$ planes are available for solving the TIE. Expanding $I_{z} (r)$ with Taylor series at $z = 0$ as:

I_{z} (r) - I_{0} (r) = z \partial I_{0} (r) + \sum_{n = 2}^{\infty} z^{n} \partial^{n} I_{0} (r) / n!

where

n!

denotes the factorial of

n

(

n = 2, 3, ..., M

).

\partial^{n} I_{0} (r)

is the

n^{t h}

order axial intensity derivative evaluated at the point

z = 0

. In the traditional TIE, the second term on the right side of Eq. (3) is ignored, which produce a linear Taylor approximation respect to

z

. To improve the accuracy, we should estimate and compensate the higher-order intensity derivative in Eq. (3). This can be done by constructing a linear system from Eq. (3) as:

[\begin{matrix} z_{1} & z_{1}^{2} & \dots & z_{1}^{M} \\ z_{2} & z_{2}^{2} & \dots & z_{2}^{M} \\ ⋮ & ⋮ & \dots & ⋮ \\ z_{M} & z_{M}^{2} & \dots & z_{M}^{M} \end{matrix}] [\begin{matrix} \partial I_{0} (r) \\ \partial^{2} I_{0} (r) / 2! \\ ⋮ \\ \partial^{M} I_{0} (r) / M! \end{matrix}] = [\begin{matrix} I_{z_{1}} (r) - I_{0} (r) \\ I_{z_{2}} (r) - I_{0} (r) \\ ⋮ \\ I_{z_{M}} (r) - I_{0} (r) \end{matrix}] .

When the intensities are equally-spaced and symmetrically distributed relative to the focal plane, Eq. (4) is equivalent to the equation in [15] (Eq. (7) in [15]). When $M = 4$ , Eq. (4) is reduced to the Eq. (7) in [16].

To solve Eq. (4), we denote the coefficient matrix in the left hand of it as $A$ , and the inverse matrix of $A$ as $B$ . Since $A$ is an $M^{t h}$ order Vandermonde matrix, the specific expression of $B$ is [17]:

B = {[b]}_{M}, b_{i j} = \sum_{\begin{array}{l} 1 \leq k_{1} < ... < k_{M - i} \leq M \\ k_{1}, ..., k_{M - i} \neq j \end{array}} {(- 1)}^{i - 1} z_{k_{1}} ... z_{k_{M - i}} / [z_{j} \prod_{\begin{array}{l} 1 \leq k \leq M \\ k \neq j \end{array}} (z_{k} - z_{j})] .

The solution to Eq. (4) can be obtained by inverting the matrix $A$ :

[\begin{matrix} \partial I_{0} (r) \\ \partial^{2} I_{0} (r) / 2! \\ ⋮ \\ \partial^{M} I_{0} (r) / M! \end{matrix}] = B [\begin{matrix} I_{z_{1}} (r) - I_{0} (r) \\ I_{z_{2}} (r) - I_{0} (r) \\ ⋮ \\ I_{z_{M}} (r) - I_{0} (r) \end{matrix}] .

Substituting Eq. (6) into Eq. (3), the

M^{t h}

order approximation to the intensity derivatives can be expressed with the form:

\partial I_{0} (r) = \sum_{m = 1}^{M} c_{m} [I_{z_{m}} (r) - I_{0} (r)] .

Where the weight coefficients have the form

c_{m} = b_{1 j} = \sum_{\begin{array}{l} 1 \leq k_{1} < ... < k_{M - 1} \leq M \\ k_{1}, ..., k_{M - 1} \neq m \end{array}} z_{k_{1}} ... z_{k_{M - 1}} / [z_{m} \prod_{\begin{array}{l} 1 \leq k \leq M \\ k \neq m \end{array}} (z_{k} - z_{m})], m = j .

Substituting Eq. (7) into Eq. (2), we can obtain a new generalization of the TIE with high order intensity derivatives:

\nabla \cdot [I_{0} (r) \nabla ϕ_{0} (r)] = - k \sum_{m = 1}^{M} c_{m} [I_{z_{m}} (r) - I_{0} (r)] .

In this study, the forward formula is used to estimate the intensity derivative. Table 1 shows the difference between the weight coefficients $c_{m}, (m = 1, ..., 4)$ with unequally-spaced planes and that with equally-spaced planes [15]. If the intensity measurement planes are equally-spaced, i.e. $z_{m} = m Δ z$ , the proposed expressions of the coefficients are reduced to that in [15].

Table 1. Comparison of the weight coefficients for the given order of accuracy.

View Table

3. Simulations

Numerical simulation experiments are conducted to demonstrate the feasibility of the proposed phase retrieval approach using intensity measurements in multiple unequally-spaced planes. A simulated optical fiber is chosen as the true phase profile (see Fig. 1 ).The dimensions of the image are 512 × 512 pixels with an assigned pixel size of 4 $μ m$ . The phase value is between 0 and 0.5 $π$ radians. The radiation wavelength is chosen to be 550nm. The intensity distribution in the input plane is uniform. The root-mean-square-error (RMSE) is employed to quantify the accuracy of the retrieved phase:

\begin{array}{l} R M S E = \sqrt{\sum_{x, y} {(θ_{0}^{'} (x, y) - θ_{0} (x, y) - a)}^{2} / (X \times Y)} \\ a = \sum_{x, y} (θ_{0}^{'} (x, y) - θ_{0} (x, y)) / (X \times Y), \end{array}

where

θ_{0}^{'} (x, y)

and

θ_{0} (x, y)

denote the calculated phase distribution and the true phase distribution at pixel

(x, y)

respectively.

X \times Y

is the sampled pixels of the phase image.

Fig. 1 Simulated phase in the focal plane.

Download Full Size | PDF

The diffraction images are calculated by the angular spectrum method [9,18] as follows:

\begin{array}{l} I_{z_{m}} = u_{z_{m}} (r) u_{z_{m}}^{*} (r), \begin{matrix}  \end{matrix} u_{z_{m}} (r) = F^{- 1} {F {u_{0} (r)} H_{z_{m}} (g)}, \\ H_{z_{m}} (g) = \exp (i k z_{m} {(1 - λ^{2} g^{2})}^{1 / 2}) . \end{array}

Where

F {•}

and

F^{- 1} {•}

respectively denote the Fourier and inverse Fourier transforms.

g = (g_{x}, g_{y})

is the spatial frequency, conjugate to

r

. The superscript

*

denotes conjugation.

z_{m} (m = 1, ..., M)

denotes the defocused distance. Figure 2 shows an example of the simulated diffraction images, which can be used to solve the TIE for phase retrieval. It should be mentioned that the diffraction images are sampled in unequally-spaced planes.

Fig. 2 An example of the simulated intensity focal stack. The defocused distances for the intensities are respectively (from left to right): $0, 4.02 μ m, 8.03 μ m, 12.00 μ m, 16.00 μ m, 20.03 μ m, 24.07 μ m$ .

Download Full Size | PDF

We carry out the phase retrieval simulations by solving the TIE with the FFT method [8]. The first simulation assumes that the intensities are measured in the planes $z_{m} = m Δ z + ε_{m}$ , $(m = 1, 2, ..., M) .$ $ε_{m}$ is a random disturbance between 0 and $0.1 μ m$ , which denotes the measurement errors of the defocused distance. This numerical experiment includes the following two cases:

Case1: Ignoring $ε_{m}$ , and the equally-spaced formula [15] is used to solve the TIE.

Case2: Considering $ε_{m}$ , and the unequally-spaced formula is used to solve the TIE.

Figure 3 shows the recovered phases and the corresponding phase errors for case 2. In this numerical experiment, $Δ z = 4 μ m$ and $M = 6$ , the defocused distances for the six intensities are respectively assigned to be $4.02, 8.03, 12.00, 16.00, 20.03 μ m$ and $24.07 μ m$ . The visual impression of the reconstructed phases is in good agreement with the true phase, and from the phase error images we can see that the visual differences between the reconstructed phases and the true phase tend to be smaller with the increase of the number of intensity measurements.

Fig. 3 Recovered phases (radians) and phase errors (radians) for case 2: The first row are listed the recovered phases for $M = 1, 2, 6$ . The second row are listed the corresponding phase errors relative to the true phase.

Download Full Size | PDF

In the other simulation experiment, the intensities are measured in $M$ random planes between $[z_{1}, z_{M}]$ . Here we choose $z_{1} = 4 μ m$ and $z_{M} = 16 μ m$ . The intervals between the adjacent defocused planes are randomly chosen. In this simulation, the defocused distances of the six intensities are respectively set to be $4, 7, 9, 13, 15 μ m$ and $16 μ m$ . We name the second simulation as case 3.

Figure 4 shows the RMSE as a function of $M$ on a logarithmic scale. Figure 4(a) shows that the defocused distance errors can effect the oscillation of the RMSE curve, which corresponds to case 1. In this case, $ε_{m}$ is ignored and the equally-spaced formula [15] is used to solve the TIE. Because the defocused distance errors are ignored, the accuracy of phase reconstruction cannot be improved by including higher order intensity derivatives. In case 2, we have considered the measurement errors, and the unequally-spaced formulas are used for phase retrieval. Consequently, the accuracy of the retrieved phases can be improved with the increase of the number of intensity measurements. The difference of the RMSE between case 2 and case 3 shown in Fig. 4(b) is due to the defocused distances. Since the noise is not considered, the smaller the defocused distances, the lower the RMSE are. Therefore, the RMSE in case 2 is larger than that in case 3.

Fig. 4 RMSEs as a function of the number of intensity measurements with a logarithmic scale. (a) case 1 and 2; (b) case 2 and case 3.

Download Full Size | PDF

4. Conclusions

We demonstrate a phase retrieval method based on the TIE from multiple unequally-spaced intensities, which helps relax the limitation of equally-spaced intensity measurements. Therefore, the proposed phase imaging method is more convenient for practical application. The numerical simulation results illustrate that the measurement error of the defocused distances effect the accuracy of phase retrieval and the proposed method can eliminate this influence. It can also improve the accuracy of the phase reconstruction with higher-order intensity derivative components just as the method with equally-spaced intensity measurements. The method has a potential use in the optical measurement system [13] for retrieving the phase with high accuracy in real time. However, the method presented in this paper does not consider the noise effect. More work will be done to simultaneously minimize the effects of the noise and the higher-order intensity derivatives on the retrieved phase from multiple unequally-spaced intensity measurements in the future.

Acknowledgments

This work is partly supported by a grant from the National Natural Science Foundation of China (No.60502018, No.61171005), the Innovation Fund and the Scholarship Award of BUAA for PhD Students, the Fundamental Research Fund for the Central Universities and the Aeronautical Special Fund of AVIC (No.CXY2010BH02).

References and links

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

2. T. E. Gureyev, A. Roberts, and K. A. Nugent, “Partially coherent fields, the transport-of-intensity equation, and phase uniqueness,” J. Opt. Soc. Am. A 12(9), 1942–1946 (1995). [CrossRef]

3. A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett. 23(11), 817–819 (1998). [CrossRef] [PubMed]

4. K. Ishizuka and B. Allman, “Phase measurement of atomic resolution image using transport of intensity equation,” J. Electron Microsc. (Tokyo) 54(3), 191–197 (2005). [CrossRef] [PubMed]

5. B. E. Allman, P. J. McMahon, K. A. Nugent, D. Paganin, D. L. Jacobson, M. Arif, and S. A. Werner, “Phase radiography with neutrons,” Nature 408(6809), 158–159 (2000). [CrossRef] [PubMed]

6. K. A. Nugent, T. E. Gureyev, D. J. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard X rays,” Phys. Rev. Lett. 77(14), 2961–2964 (1996). [CrossRef] [PubMed]

7. 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]

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

9. 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]

10. W. Xiao, M. Heng, and Z. Dazun, “Phase retrieval based on intensity transport equation,” Acta Opt. Sin. 27, 2117–2122 (2007).

11. 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. 12(9), 1932–1941 (1995). [CrossRef]

12. S. V. Pinhasi, R. Alimi, L. Perelmutter, and S. Eliezer, “Topography retrieval using different solutions of the transport intensity equation,” J. Opt. Soc. Am. A 27(10), 2285–2292 (2010). [CrossRef] [PubMed]

13. J. Frank, G. Wernicke, J. Matrisch, S. Wette, J. Beneke, and S. Altmeyer, “Quantitative determination of the optical properties of phase objects by using a real-time phase retrieval technique,” Proc. SPIE 8082, 80820N, 80820N-9 (2011). [CrossRef]

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

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

16. W. X. Cong and G. Wang, “Higher-order phase shift reconstruction approach,” Med. Phys. 37(10), 5238–5242 (2010). [CrossRef] [PubMed]

17. A. E. Knuth, The Art of Computer Programming: Volume 1: Fundamental Algorithms (Addison-Wesley Professional, 1997), pp. 38.

18. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, New York, 1996), pp. 55–61.

Intensity planes	Intensity difference	2nd order coefficients	3rd order coefficients	4th order coefficients
Unequally spaced	$I_{z_{1}} - I_{0}$	$c_{1} = \frac{z_{2}}{z_{1} (z_{2} - z_{1})}$	$c_{1} = \frac{z_{2} z_{3}}{z_{1} (z_{2} - z_{1}) (z_{3} - z_{1})}$	$c_{1} = \frac{z_{2} z_{3} z_{4}}{z_{1} (z_{2} - z_{1}) (z_{3} - z_{1}) (z_{4} - z_{1})}$
	$I_{z_{2}} - I_{0}$	$c_{2} = \frac{z_{1}}{z_{2} (z_{1} - z_{2})}$	$c_{2} = \frac{z_{1} z_{3}}{z_{2} (z_{1} - z_{2}) (z_{3} - z_{2})}$	$c_{2} = \frac{z_{1} z_{3} z_{4}}{z_{2} (z_{1} - z_{2}) (z_{3} - z_{2}) (z_{4} - z_{2})}$
	$I_{z_{3}} - I_{0}$		$c_{3} = \frac{z_{1} z_{2}}{z_{3} (z_{1} - z_{3}) (z_{2} - z_{3})}$	$c_{3} = \frac{z_{1} z_{2} z_{4}}{z_{3} (z_{1} - z_{3}) (z_{2} - z_{3}) (z_{4} - z_{3})}$
	$I_{z_{4}} - I_{0}$			$c_{4} = \frac{z_{1} z_{2} z_{3}}{z_{4} (z_{1} - z_{4}) (z_{2} - z_{4}) (z_{3} - z_{4})}$
Equally spaced	$I_{Δ z} - I_{0}$	$c_{1} = \frac{2}{Δ z}$	$c_{1} = \frac{3}{Δ z}$	$c_{1} = \frac{4}{Δ z}$
	$I_{2 Δ z} - I_{0}$	$c_{2} = - \frac{1}{2 Δ z}$	$c_{2} = - \frac{3}{2 Δ z}$	$c_{2} = - \frac{3}{Δ z}$
	$I_{3 Δ z} - I_{0}$		$c_{3} = \frac{1}{3 Δ z}$	$c_{3} = \frac{4}{3 Δ z}$
	$I_{4 Δ z} - I_{0}$			$c_{4} = - \frac{1}{4 Δ z}$

Transport of intensity phase imaging from multiple intensities measured in unequally-spaced planes

Abstract

1. Introduction

2. Unequally-spaced intensity measurements for phase retrieval

3. Simulations

4. Conclusions

Acknowledgments

References and links

Cited By

Figures (4)

Tables (1)

Equations (11)

Optics Express