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

doi:10.1364/OE.20.000972

« Show journal navigation

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

Shiling Zheng, Bindang Xue, Wenfang Xue, Xiangzhi Bai, and Fugen Zhou »View Author Affiliations

Optics Express, Vol. 20, Issue 2, pp. 972-985 (2012)
http://dx.doi.org/10.1364/OE.20.000972

View Full Text Article

Acrobat PDF (1355 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. TIE with unequally-spaced noisy intensity measurements
3 Accuracy analysis of the estimated intensity derivative
4. Simulations
5. Conclusions
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (21)
Cited By
Figures (10)
Metrics

Abstract

The noise problem is generally inevitable for phase retrieval by solving the transport of intensity equation (TIE). The noise effect can be alleviated by using multiple intensities to estimate the axial intensity derivative in the TIE. In this study, a method is proposed for estimating the intensity derivative by using multiple unevenly-spaced noisy measurements. The noise-minimized intensity derivative is approximated by a linear combination of the intensity data, in which the coefficients are obtained by solving a constrained optimization problem. The performance of the method is investigated by both the error analysis and the numerical simulations, and the results show that the method can reduce the noise effect on the retrieved phase. In addition, guidelines for the choice of the number of the intensity planes are given.

1. Introduction

Phase imaging is an important technique for studying phase objects in biological or material research. TIE [1

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

, 2

F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt. 27(7), 1223–1225 (1988). [CrossRef] [PubMed]

] provides a convenient and deterministic way for quantitative phase reconstruction from intensity measurements. It relates the axial intensity derivative to the object-plane phase by a second order elliptic partial differential equation, and the phase can be reconstructed uniquely (within an additive constant) by solving the equation with appropriate boundary conditions [3

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]

–5

J. Frank, S. Altmeyer, and G. Wernicke, “Non-interferometric, non-iterative phase retrieval by Green’s functions,” J. Opt. Soc. Am. A 27(10), 2244–2251 (2010). [CrossRef] [PubMed]

]. Many effective numerical methods have been proposed to solve the TIE, such as the Green’s function based method [1

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

, 5

J. Frank, S. Altmeyer, and G. Wernicke, “Non-interferometric, non-iterative phase retrieval by Green’s functions,” J. Opt. Soc. Am. A 27(10), 2244–2251 (2010). [CrossRef] [PubMed]

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

], fast Fourier transform (FFT) method [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

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]

], and multigrid method [8

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]

, 9

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]

]. Therefore, TIE has been applied in a variety of applications including adaptive optics [9

, 10

C. Roddier and F. Roddier, “Wave-front reconstruction from defocused images and the testing of ground-based optical telescopes,” J. Opt. Soc. Am. A 10(11), 2277–2287 (1993). [CrossRef]

], microscopy [11

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

, 12

J. Frank, J. Matrisch, J. Horstmann, S. Altmeyer, and G. Wernicke, “Refractive index determination of transparent samples by noniterative phase retrieval,” Appl. Opt. 50(4), 427–433 (2011). [CrossRef] [PubMed]

] and X-ray optics [13

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]

, 14

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

] etc. Although different methods may retrieve the phase distribution with different precisions, a fundamental factor which determines the accuracy of the retrieved phase is the precision of the axial intensity derivative. The intensity derivative is generally calculated from two defocused intensity measurements by the finite difference method. However, the main problem inherent in this method is that it ignores the nonlinear effect of the higher-order intensity derivatives and is vulnerable to the noise effect [15

D. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent, “Quantitative phase-amplitude microscopy. III. The effects of noise,” J. Microsc. 214(1), 51–61 (2004). [CrossRef] [PubMed]

]. Multiple intensity measurements based methods [16

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

–18

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

] can mitigate the nonlinear effect by estimating the higher-order intensity derivatives and thus can improve the accuracy of the retrieved phase. However, these methods still suffer from the noise problem. Although the noise can be reduced by multiple intensity measurements weighted average in one plane [16

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

] or by de-noising techniques [17

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

, 19

L. Waller, M. Tsang, S. Ponda, S. Y. Yang, and G. Barbastathis, “Phase and amplitude imaging from noisy images by Kalman filtering,” Opt. Express 19(3), 2805–2814 (2011). [CrossRef] [PubMed]

], this will add inconvenience and complication to the calculation. In this case, the noise-suppressed methods [20

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

] are more preferred to estimate the intensity derivative from multiple intensity measurements. In addition, the intensity derivative is mostly computed from multiple intensities measured in equally-spaced planes. These methods, however, cannot apply in arbitrarily-spaced intensity measurements. Recently, the authors have developed a method which uses multiple unevenly sampled noise-free intensities to obtain a higher-order approximation to the TIE [18

]. However, this method did not consider the noise problem.

In this study, the noisy intensity measurements based TIE is proposed. In this equation, the intensity derivative is approximated by a linear combination of the multiple unequally-spaced noisy intensities. The weight coefficients are derived by minimizing the mean square noise error of the estimated and the ideal intensity derivatives, and then expressed by explicit formulas. This approach can be regarded as a generalization of the method presented in [20

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

], which adopted multiple equally-spaced intensity images to estimate the intensity derivative. In this work, the noise-reduction property of the method is analyzed, and the error caused by the higher-order intensity derivatives is also discussed to provide guidelines for the choice of the newly used intensities. Numerical simulations are conducted to test the method using multiple noisy intensities contaminated by different levels of Gaussian noise and Poisson noise.

2. TIE with unequally-spaced noisy intensity measurements

Suppose a coherent monochromatic plane wave with a wavelength

λ

propagates along the

Z

axis in the coordinate system

(r, z)

, where

r = (x, y)

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

Z

axis. With the paraxial approximation, TIE can be derived from the parabolic wave equation [1

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

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

(1)

where

k = 2 π / λ

is the wave number,

\nabla = (\partial x, \partial y)

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

I_{z} (r)

and

ϕ_{z} (r)

are the intensity and the phase distributions respectively.

Estimating the axial intensity derivative is essential to retrieving a precise phase by solving Eq. (1) [16

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

]. Let us assume that the unevenly sampled intensities

I_{z_{j}}

and

I_{z_{- j}}

, with the defocused distances

z_{j}

and

z_{- j}

(j = 1, 2, ..., M)

, are given by

I_{z_{j}} (r) = {\bar{I}}_{z_{j}} (r) + n_{z_{j}} (r) .

(2)

I_{z_{- j}} (r) = {\bar{I}}_{z_{- j}} (r) + n_{z_{- j}} (r) .

(3)

where,

{\bar{I}}_{z_{j}} (r)

and

{\bar{I}}_{z_{- j}} (r)

represent the ideal intensities.

n_{z_{j}}

n_{z_{- j}}

are supposed to be zero-mean, additive Gaussian noise with variance

σ^{2}

, independent of the intensities. It should be mentioned that

z_{j}

and

z_{- j}

(j = 1, 2, ..., M)

can be arbitrary values at which TIE is valid. The smaller of them, the closer of the intensity planes to the focal plane.

Then, the intensity derivative in the focal plane

z = 0

, i.e.

\partial I_{0} (r) = \partial I_{z} (r) / \partial z |_{z = 0}

, can be estimated by:

\partial {\hat{I}}_{0} (r) = \sum_{j = 1}^{M} \frac{a_{j}}{z_{j} - z_{- j}} (I_{z_{j}} (r) - I_{z_{- j}} (r)) .

(4)

It should be mentioned that

z_{- j}

can be negative or just be

z_{0}

, which means that we can estimate the axial intensity derivative independently from either the double-sided or single-sided diffraction intensities. The coefficients

a_{j}, j = 1, 2, ..., M

must minimize the mean square error between the estimated and the real intensity derivatives

ε^{2} (r) = E {{[\partial {\hat{I}}_{0} (r) - \partial I_{0} (r)]}^{2}},

(5)

where

ε^{2} (r)

represents the mean square error,

E {•}

denotes the mean value. In this study,

ε^{2} (r)

is divided into two part, one is caused by the noise and the other is caused by the higher order intensity derivatives. To solve Eq. (5) for the coefficients, Taylor expanding the ideal intensity

{\bar{I}}_{z} (r)

z = 0

and with the Taylor theorem, we can get

{\bar{I}}_{z_{j}} (r) = {\bar{I}}_{0} (r) + z_{j} \partial I_{0} (r) + z_{j}^{2} \partial^{2} I_{c_{j}} (r) / 2, c_{j} \in [0, z_{j}] .

(6)

{\bar{I}}_{z_{- j}} (r) = {\bar{I}}_{0} (r) + z_{- j} \partial I_{0} (r) + z_{- j}^{2} \partial^{2} I_{c_{- j}} (r) / 2, c_{- j} \in [z_{- j}, 0] .

(7)

where

\partial^{2} I_{c_{\pm j}} (r) = {[\partial^{2} I_{z} (r) / \partial z^{2}]}_{z = c_{\pm j}}

. Substituting Eq. (6) and Eq. (7) into Eq. (4), we can obtain

\partial {\hat{I}}_{0} (r) = \sum_{j = 1}^{M} \frac{a_{j}}{z_{j} - z_{- j}} (n_{z_{j}} (r) - n_{z_{- j}} (r)) + \partial I_{0} (r) \sum_{j = 1}^{M} a_{j} + \sum_{j = 1}^{M} \frac{a_{j}}{z_{j} - z_{- j}} (\frac{z_{j}^{2}}{2} \partial^{2} I_{c_{j}} (r) - \frac{z_{- j}^{2}}{2} \partial^{2} I_{c_{- j}} (r)) .

(8)

Then, substituting Eq. (8) into Eq. (5) we can get

\sum_{j = 1}^{M} a_{j} - 1 = 0.

(9)

ε^{2} (r) = ε_{n}^{2} (r) + ε_{h}^{2} (r) is minimum,

(10)

with

ε_{n}^{2} (r) = 2 σ^{2} \sum_{j = 1}^{M} \frac{a_{j}^{2}}{{(z_{j} - z_{- j})}^{2}},

(11)

and

ε_{h}^{2} (r) = \frac{1}{4} {\sum_{j = 1}^{M} \frac{a_{j}}{z_{j} - z_{- j}} [z_{j}^{2} \partial^{2} I_{c_{j}} (r) - z_{- j}^{2} \partial^{2} I_{c_{- j}} (r)]}^{2} .

(12)

where

ε_{n}^{2} (r)

and

ε_{h}^{2} (r)

represent the error caused by noise and the higher order intensity derivatives respectively. When the noise is not considered, i.e.

ε_{n}^{2} (r)

is zero, the coefficients

a_{j}

j = 1, 2, ..., M

can be calculated by the method introduced in [18

]. It calculates the higher-order intensity derivatives by solving the linear equations with Vandermonde matrix as the coefficient matrix.

In this study, our objective is to minimize the noise effects on the estimated intensity derivative. Thus, from Eq. (11) we obtain:

\sum_{j = 1}^{M} \frac{a_{j}^{2}}{{(z_{j} - z_{- j})}^{2}} = minimum .

(13)

Equation (9) and Eq. (13) can construct a nonlinear optimization problem with an equality constraint. With the Lagrange multiplier method, the coefficients can be calculated by

a_{j} = \frac{{(z_{j} - z_{- j})}^{2}}{\sum_{j = 1}^{M} {(z_{j} - z_{- j})}^{2}} .

(14)

Substituting Eq. (14) into Eq. (4) and then substituting the obtained equation into Eq. (1), we can get the generalized TIE with unequally-spaced noisy intensity measurements as:

\nabla \cdot [I_{z} (r) \nabla ϕ_{z} (r)] = - k \sum_{j = 1}^{M} s_{j} (I_{z_{j}} (r) - I_{z_{- j}} (r)),

(15)

where the weight coefficients are

s_{j} = \frac{(z_{j} - z_{- j})}{\sum_{j = 1}^{M} {(z_{j} - z_{- j})}^{2}} .

(16)

Equation (15) can be regarded as a generalization of Eq. (3) in [20

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

] from equally-spaced intensities to unequally-spaced intensities. It can be verified easily that when the measured intensities are equally-spaced, i.e.:

z_{j} = j Δ z, z_{- j} = - j Δ z, j = 1, 2, ...

(17)

The coefficients for the intensity measurements become:

s_{j} = \frac{3 j}{M (M + 1) (2 M + 1)}, j = \pm 1, \pm 2, ...,

(18)

which correspond to that in Eq. (9) in [20

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

In the next section, we will analyze the effects of the defocused values on

ε_{h}^{2} (r)

and

ε_{h}^{2} (r)

. In addition, we will discuss the relationship between the number of the intensity planes and the accuracy of the estimated intensity derivative.

3 Accuracy analysis of the estimated intensity derivative

Substituting Eq. (14) into Eq. (11) and Eq. (12), we have:

ε_{n}^{2} (r) = 2 σ^{2} / \sum_{j = 1}^{M} {(z_{j} - z_{- j})}^{2},

(19)

ε_{h}^{2} (r) = \frac{1}{4} {\sum_{j = 1}^{M} \frac{z_{j} - z_{- j}}{\sum_{j = 1}^{M} {(z_{j} - z_{- j})}^{2}} [z_{j}^{2} \partial^{2} I_{c_{j}} (r) - z_{- j}^{2} \partial^{2} I_{c_{- j}} (r)]}^{2} .

(20)

As shown in Eq. (19), the noise error becomes smaller with the increase of

M

. However,

ε_{h}^{2} (r)

is indeterminable as

M

increases. This is because, firstly, the second-order intensity derivatives cannot be precisely known in advance, and secondly, the sign of

(z_{j} - z_{- j}), j = 1, 2, ..., M

is uncertain when

z_{j}

and

z_{- j}

are arbitrary values. Consequently, we cannot make sure whether

ε_{h}^{2} (r)

is decreasing by increasing the number of the intensity planes. In this study, we will provide a general idea of it by analyzing

ε_{h}^{2} (r)

with some prescribed assumptions.

Let us assume that

z_{j}

is positive and

z_{- j}

is negative. The lower bound and the upper bound of

ε_{h}^{2} (r)

are adopted to analyze its effects, and they can be derived from Eq. (20) as:

ε_{h (l o)}^{2} \leq ε_{h}^{2} (r) \leq ε_{h (u p)}^{2},

(21)

with

\begin{array}{l} ε_{h (u p)}^{2} = U^{2} Z_{u p}^{2} (M) / 4, Z_{u p}^{2} (M) = \frac{{[\sum_{j = 1}^{M} (z_{j} - z_{- j}) (z_{j}^{2} + z_{- j}^{2})]}^{2}}{{[\sum_{j = 1}^{M} {(z_{j} - z_{- j})}^{2}]}^{2}}, \\ ε_{h (l o)}^{2} = 0, \end{array}

and

U = \max_{r, c} (| \partial^{2} I_{c} (r) |), c \in [z_{\min}, z_{\max}], z_{\min} \leq z_{j}, z_{- j} \leq z_{\max}, j = 1, 2, ..., M .

where

[z_{\min}, z_{\max}]

belongs to the interval where TIE is valid and the defocused values of the used intensity measurements always lie in this interval.

The accuracy of the estimated intensity derivative can be improved by using more intensity planes [16

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

, 20

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

], so we express

ε^{2} (r)

ε_{M + 1}^{2} (r) < ε_{M}^{2} (r) \Rightarrow ε_{h_(M + 1)}^{2} (r) - ε_{h_M}^{2} (r) < ε_{n_M}^{2} (r) - ε_{n_(M + 1)}^{2} (r),

(22)

where

ε_{M}^{2} (r)

ε_{n_M}^{2} (r)

and

ε_{h_M}^{2} (r)

denote respectively the total error, the noise error and the error caused by the higher intensity derivatives with

M

intensity planes. While,

ε_{M + 1}^{2} (r)

ε_{n_(M + 1)}^{2} (r)

and

ε_{h_(M + 1)}^{2} (r)

are those with

M + 1

intensity planes.

To further analyze

ε^{2} (r)

, we here assume

ε_{h}^{2} (r)

approximates to

ε_{h (u p)}^{2}

. Then substituting

ε_{h}^{2} (r)

ε_{h (u p)}^{2}

so that Eq. (22) can be expressed as

ε_{h (u p)_(M + 1)}^{2} (r) - ε_{h (u p)_M}^{2} (r) < ε_{n_M}^{2} (r) - ε_{n_(M + 1)}^{2} (r) .

(23)

where

ε_{h (u p)_M}^{2} (r)

ε_{h (u p)}^{2}

and

ε_{h (u p)_(M + 1)}^{2} (r)

are respectively the upper bounds of

ε_{h_M}^{2} (r)

and

ε_{h_(M + 1)}^{2} (r)

In Eq. (23), if

ε_{h (u p)_(M + 1)}^{2} (r) - ε_{h (u p)_M}^{2} (r) < 0

, i.e.

ε_{h (u p)}^{2}

is a decreasing function of

M

, The inequality is always true. From the expression of in Eq. (21), we know that this can be achieved by properly choosing the defocused values to make

Z_{u p}^{2} (M)

decrease. As an example, we can choose the defocused values of the already and newly used intensities as:

{\begin{matrix} z_{M + 1} < z_{j + 1} < z_{j}, \\ z_{M + 1} \leq | z_{- (j + 1)} | < | z_{- j} |, \\ | z_{- (M + 1)} | < | z_{- (j + 1)} |, \\ | z_{- (M + 1)} | \leq z_{j + 1} . \end{matrix} (j = 1, 2, ..., M - 1) .

(24)

This means that the newly used intensities with the defocused distances

z_{M + 1}

and

z_{- (M + 1)}

are closer to the focal plane than the already used intensities with defocused values

z_{j}

and

z_{- j}

(j = 1, 2, ..., M)

. Then we have:

\begin{array}{l} (z_{M + 1}^{2} + z_{- (M + 1)}^{2}) (z_{j + 1} - z_{- (j + 1)}) - (z_{j + 1}^{2} + z_{- (j + 1)}^{2}) (z_{M + 1} - z_{- (M + 1)}) \\ = z_{M + 1} z_{j + 1} (z_{M + 1} - z_{j + 1}) - z_{M + 1} z_{- (j + 1)} (z_{M + 1} + z_{- (j + 1)}) \\ + z_{- (M + 1)} z_{j + 1} (z_{- (M + 1)} + z_{j + 1}) - z_{- (M + 1)} z_{- (j + 1)} (z_{- (M + 1)} - z_{- (j + 1)}) < 0 \end{array}

(25)

Consequently, we can get:

\begin{array}{l} Z_{u p} (M + 1) - Z_{u p} (M) \\ = \frac{\sum_{j = 1}^{M + 1} (z_{j} - z_{- j}) (z_{j}^{2} + z_{- j}^{2})}{\sum_{j = 1}^{M + 1} {(z_{j} - z_{- j})}^{2}} - \frac{\sum_{j = 1}^{M} (z_{j} - z_{- j}) (z_{j}^{2} + z_{- j}^{2})}{\sum_{j = 1}^{M} {(z_{j} - z_{- j})}^{2}} \\ = \frac{\sum_{j = 1}^{M} (z_{M + 1} - z_{- (M + 1)}) (z_{j} - z_{- j}) [(z_{M + 1}^{2} + z_{- (M + 1)}^{2}) (z_{j} - z_{- j}) - (z_{j}^{2} + z_{- j}^{2}) (z_{M + 1} - z_{- (M + 1)})]}{\sum_{j = 1}^{M} {(z_{j} - z_{- j})}^{2} \sum_{j = 1}^{M + 1} {(z_{j} - z_{- j})}^{2}} < 0. \end{array}

(26)

Therefore,

Z_{u p}^{2} (M)

decreases with

M

increasing.

However, If

ε_{h (u p)_(M + 1)}^{2} (r) - ε_{h (u p)_M}^{2} (r) > 0

, i.e.,

ε_{h (u p)}^{2}

is an increasing function of

M

, The inequality in Eq. (23) may not be always true. For example, when the defocused values satisfy

{\begin{matrix} z_{M + 1} > z_{j + 1} > z_{j}, \\ z_{M + 1} \geq | z_{- (j + 1)} | > | z_{- j} |, \\ | z_{- (M + 1)} | > | z_{- (j + 1)} |, \\ | z_{- (M + 1)} | \geq z_{j + 1} . \end{matrix} \begin{matrix}  \end{matrix} (j = 1, 2, ..., M - 1),

(27)

Z_{u p}^{2} (M)

increases with

M

. This is just a contrary situation to the example in the last paragraph, in which the newly used intensities are farther to the focal plane than the already used intensities. In this case, when

ε_{n}^{2} (r)

is the dominant factor in

ε^{2} (r)

, it can be significantly reduced by using more intensity planes. Whereas as the reducing of

ε_{n}^{2} (r)

and the increasing of

ε_{h (u p)}^{2}

ε_{h (u p)}^{2}

may becomes dominant when the number of the used intensity planes increases to a critical value of

M_{c}

. The value of

M_{c}

is determined by the defocused values and the noise levels. It cannot be precisely determined in advance.

In summary,

ε_{n}^{2} (r)

can be reduced by using more intensity measurements, and

ε_{h}^{2} (r)

can be reduced when

Z_{u p}^{2} (M)

decreases.

4. Simulations

Numerical simulations are conducted to test the derived formulas for the intensity derivative estimation. Figure 1 shows the simulated phase profile with phase value between 0 and 0.5

π

radians (0 denotes the darkest pixels and 0.5

π

refers to the brightest pixels). The dimensions of the image are 512 × 512 pixels with the assigned pixel size 4

μ m

. The radiation wavelength is chosen to be 550nm. The intensity distribution in the input plane is uniform.

Fig. 1 Simulated phase distribution in the

z = 0

plane.

The intensities are sampled as two-dimensional images, so that the normalized mean square error (NMSE) is used to quantitatively evaluate the accuracy of the estimated intensity derivative, which is:

E^{2} = \sum_{x = 1}^{X} \sum_{y = 1}^{Y} {[\partial {\hat{I}}_{0} (x, y) - \partial I_{0} (x, y)]}^{2} / \sum_{x = 1}^{X} \sum_{y = 1}^{Y} {[\partial I_{0} (x, y)]}^{2} .

(28)

where

{\hat{I}}_{0} (x, y)

and

I_{0} (x, y

) denote respectively the measured intensity and the ideal intensity. X × Y denotes the number of the sampled pixels in the images.

The simulated diffraction intensities in the different planes are calculated by the angular spectrum method [21

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

]. Therefore, the ideal intensity derivative can be derived as (see Appendix A):

\partial I_{0} = - 2 Im [u_{0}^{*} (r) F^{- 1} {F {u_{0} (r)} k {(1 - λ^{2} g^{2})}^{1 / 2}}] .

(29)

Where

F {•}

and

F^{- 1} {•}

denote respectively the Fourier and inverse Fourier transforms.

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

denotes the spatial frequency, conjugate to

r

. The superscript

*

stands for the complex conjugate. The TIE is solved by the fast Fourier transform (FFT) method subject to the periodic boundary conditions [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

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]

, 14

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

], and the root-mean-square-error (RMSE) is employed to quantify the accuracy of the retrieved phases:

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

(30)

where

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

and

θ_{0} (x, y)

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

(x, y)

In this study, thirty diffraction intensity images are used o test the method with random defocused distances between

- 28 μ m

and

29 μ m

.The positive defocused values of

z_{j}

are assigned to be

29,27,26,23,22,20,19,17,16,12,11,10,8,5 μ m

and

3 μ m

. The negative defocus values of

z_{- j}

are assigned to be

-28,-25,-23,-21,-18,-16,-15,-13,-10,-9,-7,-6,-3 μ m

and

- 2 μ m

. The maximum and minimum sampled spaces are respectively assigned to be

4 μ m

and

1 μ m

. These defocused values can make

Z_{u p} (M)

diminish with

M

increases. Intensities with Gaussian noise, Poisson noise and the mix of them will be simulated and tested respectively. For each numerical experiment, the method for each noise levels will be carried out for 60 times, and the averages of the corresponding NMSE and RMSE will be calculated.

4.1 Estimate the intensity derivative from multiple unequally-spaced intensity measurements with Gaussian noise

We assume that the intensities are contaminated with the independent Gaussian white noise and four noise levels in the images are simulated. The corresponding signal-to-noise-ratio (SNR) of the images is 70db, 60db, 55db and 50db, respectively. Figure 2 shows an example of the diffraction images in the

z = 19 μ m

plane with different SNR.

Fig. 2 Simulated noisy diffraction intensities in the

z = 19 μ m

plane with different SNRs of (a)

70 d b

; (b)

60 d b

; (c)

55 d b

; (d)

50 d b

Figure 3 shows the averaged NMSEs as a function of

M

with different SNR. Figure 4 shows the retrieved phases with different number of noisy intensity measurements and with different SNR. In Table 1 , NMSEs and RMSEs are compared with different number of noisy intensity measurements and with different levels of SNR.

Fig. 3 The average NMSEs of 60 times for each noise level are plotted as a function of M and SNR. The noise in the intensities is assumed to be Gaussian white noise.

Fig. 4 Retrieved phases from

M = 1

(the first row) and

M = 15

(the second row) Gaussian-noise damaged intensity measurements with different SNRs of 70db, 60db, 55db and 50db (from left to right). The RMSEs for the phase images in the first row are respectively 0.042, 0.180, 0.330 and 0.458, and for those in the second row are respectively 0.017, 0.043, 0.139 and 0.169.

Table 1 Comparison of the errors between the estimated axial intensity derivatives and the reconstructed phases by solving TIE from different numbers of intensities with different levels of Gaussian noise.

SNR(db)	M	NMSE	RMSE
70	1	0.006	0.043
70	15	0.002	0.019
60	1	0.019	0.148
60	15	0.004	0.063
55	1	0.051	0.243
55	15	0.010	0.109
50	1	0.150	0.460
50	15	0.028	0.184

As shown in Fig. 3, the estimation accuracy is improved with the increase of

M

and the NMSEs decrease monotonically. The NMSEs curves also show that the method is effective for phase retrieval by multiple intensities with much noise, since the NMSE for SNR of 50db is reduced from 0.150 to 0.028 (see Table 1), and the precision of the estimated derivative is improved about 81%. In Table 1, the NMSEs are diminished from 0.006, 0.019, and 0.051 to 0.002, 0.004, and 0.010, with an improvement of about 67%, 79%, and 80% respectively for SNR of 70db, 60db, 55db. However, the convergence speed of the NMSEs becomes slower as

M

increases in Fig. 3. This can be explained by Eq. (19), which shows that the intensities with smaller defocused values have less impact on the noise error, so that the noise error reduces slower or even become stagnation as the defocused values decrease.

In Table 1, RMSEs show that the accuracy for the phases retrieved from more intensity measurements are much higher. The visual impression in Fig. 4 is that the cloudy effect in the retrieved phases induced by the measurement noise is largely removed when multiple intensities are adopted.

4.2 Estimate the intensity derivative from multiple unequally-spaced intensity measurements with Poisson noise

In this simulation, we assume that the intensities are contaminated with Poisson noise, and the noise is simulated respectively with an average of 100,000, 80,000, 32,000, and 10,000 photons per pixel in the diffraction image. We denote the average photons per pixel in the images as

γ

. The less number of the average photons means that the images are measured under the lower-light conditions. In this case, the Poisson-noise problem of the images is more severe.

Figure 5 shows the noisy diffraction image in the

z = 19 μ m

plane. As shown, the intensity images with smaller

γ

are with stronger noise. Figure 6 shows the averaged NMSEs as a function of

M

and

γ

. As shown, the trends of the curves are similar to those in the Gaussian-noise case and the NMSEs for all the four cases are monotonically decreasing as

M

increases.

Fig. 5 Simulated noisy diffraction intensities in the

z = 19 μ m

plane with Poisson noise. The average photons for each pixel of the images are respectively: (a) 100,000; (b) 80,000; (c) 32,000; and (d) 10,000.

Fig. 6 The NMSEs between the estimated intensity derivative and the ideal intensity derivative are plotted as a function of M and the average photons per pixel in the diffraction images. The noise in the intensities are assumed to be Poisson noise.

Quantitative comparisons are listed in Table 2 . As shown, RMSEs are smaller with more intensity measurements. This indicates that the accuracy of the reconstructed phases from more intensity measurements is much higher. Figure 7 shows the corresponding retrieved phases of the four cases above mentioned. The first row listed the retrieved phases retrieved from intensities with

M = 1

and the second row displayed those from intensities with

M = 15

. Visually, the phases are much better when multiple intensities are used and more close to the simulated phase.

Table 2 Comparison of the errors for the estimated axial intensity derivatives and the reconstructed phases by solving TIE from different numbers of intensities with different levels of Poisson noise.

$γ$	M	NMSE	RMSE
$1 \times 10^{5}$	1	0.150	0.454
$1 \times 10^{5}$	15	0.028	0.120
$8 \times 10^{4}$	1	0.186	0.491
$8 \times 10^{4}$	15	0.034	0.209
$3.2 \times 10^{4}$	1	0.459	0.771
$3.2 \times 10^{4}$	15	0.083	0.342
$1 \times 10^{4}$	1	1.459	1.358
$1 \times 10^{4}$	15	0.262	0.561

Fig. 7 Retrieved phases from

M = 1

(the first row) and

M = 15

(the second row) Poisson-noise damaged intensity measurements with average photons per pixel of 100,000, 80,000, 32,000, and 10,000 (from left to right). The phase images in the first row have the RMSEs of 0.542, 0.422, 0.730, and 1.682 respectively, and those in the second row have the RMSEs of 0.223, 0.223, 0.300, and 0.750 respectively.

4.3 Estimate the intensity derivative from multiple unequally-spaced intensity measurements with mixed noise

In this section, we assume that the intensities are contaminated with the mixed noise, i.e. a combination of the Poisson noise and white Gaussian noise. The noisy images are simulated by firstly adding the Poisson noise to the intensities and then the independent white Gaussian noise to the obtained Poisson-noise corrupted intensities. We denote the standard deviation of the Gaussian noise as

σ

. Four cases of the diffraction intensities are tested: (1)

γ = 500, 000

σ = 0.001

; (2)

γ = 500, 000

σ = 0.01

; (3)

γ = 32, 000

σ = 0.001

; (4)

γ = 32, 000

σ = 0.01

. Figure 8 shows the noisy intensity images in the

z = 19 μ m

plane. The SNRs for the intensity images in the four cases are respectively 55db, 39db, 44db and 38db. The curves for the averaged NMSEs and the retrieved phases for the four cases are respectively shown in Fig. 9 and Fig. 10 . The conclusions are all in accordance with those obtained in Sections 4.1 and 4.2. When

M = 1

, the RMSEs for the four retrieved phases are respectively 0.231, 1.366, 0.6271 and 1.207. While the RMSEs with

M = 15

are much smaller, which are 0.075, 0.677, 0.263 and 0.802.

Fig. 8 Simulated diffraction intensities corrupted with mixed noise in the

z = 19 μ m

plane. (a) case 1; (b) case 2; (c) case 3; and (d) case 4.

Fig. 9 The NMSEs between the estimated intensity derivative and the ideal intensity derivative are plotted as a function of M for the four cases (as indicated in the main text). The noise in the intensities are assumed to be mixed noise.

Fig. 10 Retrieved phases from

M = 1

(the first row) and

M = 15

(the second row) mixed-noise damaged intensity measurements for the four cases (please refer to the main text for details) (from left to right).

5. Conclusions

We presented a method for estimating the axial intensity derivative in the TIE from multiple noisy unequally-spaced intensity measurements. An important characteristic of the method is that it could minimize the noise effect on the retrieved phase by using multiple unequally-spaced intensities and this has been proven by the mathematical error analysis. When multiple intensities corrupted with Gaussian noise, Poisson noise and mixed noise, the simulation results illustrate that the accuracy of the retrieved phase is advanced by our method.

Appendices

Appendix A

The propagation of the wave function

u_{z} (r)

in the

z \neq 0

plane from the

z = 0

plane can be calculated by the angular spectrum method:

\begin{array}{l} u_{z} (r) = F^{- 1} {F {u_{0} (r)} H_{z} (g)}, \\ H_{z} (g) = \exp (j k z {(1 - λ^{2} g^{2})}^{1 / 2}) . \end{array}

(A1)

where

F {•}

and

F^{- 1} {•}

denote respectively the Fourier and inverse Fourier transforms.

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

is the spatial frequency, conjugate to

r

. So, the analytical expression of the first-order intensity derivatives with respect to

z

can be obtained as:

\partial I_{z} = \partial_{z} (u_{z} (r) u_{z}^{*} (r)) = 2 Re [u_{z}^{*} (r) \partial_{z} u_{z} (r)] = 2 Re [u_{z}^{*} (r) F^{- 1} {F {u_{0} (r)} j k {(1 - λ^{2} g^{2})}^{1 / 2} H_{z} (g)}] .

(A2)

Let

z = 0

, Eq. (A2) turns to:

\partial I_{0} = - 2 Im [u_{0}^{*} (r) F^{- 1} {F {u_{0} (r)} k {(1 - λ^{2} g^{2})}^{1 / 2}}] .

(A3)

Acknowledgments

This work is partly supported by a grant from the National Natural Science Foundation of China (No.60502018, 61171005), the Innovation Foundations of BUAA for PhD Students, the Fundamental Research Funds for the Central Universities and the Aeronautical Special Founds 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.	F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt. 27(7), 1223–1225 (1988). [CrossRef] [PubMed]
3.	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]
4.	V. V. Volkov, Y. Zhu, and M. De Graef, “A new symmetrized solution for phase retrieval using the transport of intensity equation,” Micron 33(5), 411–416 (2002). [CrossRef] [PubMed]
5.	J. Frank, S. Altmeyer, and G. Wernicke, “Non-interferometric, non-iterative phase retrieval by Green’s functions,” J. Opt. Soc. Am. A 27(10), 2244–2251 (2010). [CrossRef] [PubMed]
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. 12(9), 1932–1941 (1995). [CrossRef]
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.	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]
9.	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]
10.	C. Roddier and F. Roddier, “Wave-front reconstruction from defocused images and the testing of ground-based optical telescopes,” J. Opt. Soc. Am. A 10(11), 2277–2287 (1993). [CrossRef]
11.	A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett. 23(11), 817–819 (1998). [CrossRef] [PubMed]
12.	J. Frank, J. Matrisch, J. Horstmann, S. Altmeyer, and G. Wernicke, “Refractive index determination of transparent samples by noniterative phase retrieval,” Appl. Opt. 50(4), 427–433 (2011). [CrossRef] [PubMed]
13.	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]
14.	D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. 80(12), 2586–2589 (1998). [CrossRef]
15.	D. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent, “Quantitative phase-amplitude microscopy. III. The effects of noise,” J. Microsc. 214(1), 51–61 (2004). [CrossRef] [PubMed]
16.	L. Waller, L. Tian, and G. Barbastathis, “Transport of intensity imaging with higher order derivatives,” Opt. Express 18(12), 12552–12561 (2010). [CrossRef] [PubMed]
17.	W. X. Cong and G. Wang, “Higher-order phase shift reconstruction approach,” Med. Phys. 37(10), 5238–5242 (2010). [CrossRef] [PubMed]
18.	B. Xue, S. Zheng, L. Cui, X. Bai, and F. Zhou, “Transport of intensity phase imaging from multiple intensities measured in unequally-spaced planes,” Opt. Express 19(21), 20244–20250 (2011). [CrossRef] [PubMed]
19.	L. Waller, M. Tsang, S. Ponda, S. Y. Yang, and G. Barbastathis, “Phase and amplitude imaging from noisy images by Kalman filtering,” Opt. Express 19(3), 2805–2814 (2011). [CrossRef] [PubMed]
20.	M. Soto and E. Acosta, “Improved phase imaging from intensity measurements in multiple planes,” Appl. Opt. 46(33), 7978–7981 (2007). [CrossRef] [PubMed]
21.	J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996), pp. 55–61.

OCIS Codes
(000.3860) General : Mathematical methods in physics
(100.3010) Image processing : Image reconstruction techniques
(100.5070) Image processing : Phase retrieval

ToC Category:
Image Processing

History
Original Manuscript: November 28, 2011
Revised Manuscript: December 21, 2011
Manuscript Accepted: December 21, 2011
Published: January 4, 2012

Citation
Shiling Zheng, Bindang Xue, Wenfang Xue, Xiangzhi Bai, and Fugen Zhou, "Transport of intensity phase imaging from multiple noisy intensities measured in unequally-spaced planes," Opt. Express 20, 972-985 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-2-972

Sort: Author | Year | Journal | Reset

References

M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am.73(11), 1434–1441 (1983). [CrossRef]
F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt.27(7), 1223–1225 (1988). [CrossRef] [PubMed]
T. E. Gureyev, A. Roberts, and K. A. Nugent, “Partially coherent fields, the transport-of-intensity equation, and phase uniqueness,” J. Opt. Soc. Am. A12(9), 1942–1946 (1995). [CrossRef]
V. V. Volkov, Y. Zhu, and M. De Graef, “A new symmetrized solution for phase retrieval using the transport of intensity equation,” Micron33(5), 411–416 (2002). [CrossRef] [PubMed]
J. Frank, S. Altmeyer, and G. Wernicke, “Non-interferometric, non-iterative phase retrieval by Green’s functions,” J. Opt. Soc. Am. A27(10), 2244–2251 (2010). [CrossRef] [PubMed]
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]
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]
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]
S. V. Pinhasi, R. Alimi, L. Perelmutter, and S. Eliezer, “Topography retrieval using different solutions of the transport intensity equation,” J. Opt. Soc. Am. A27(10), 2285–2292 (2010). [CrossRef] [PubMed]
C. Roddier and F. Roddier, “Wave-front reconstruction from defocused images and the testing of ground-based optical telescopes,” J. Opt. Soc. Am. A10(11), 2277–2287 (1993). [CrossRef]
A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett.23(11), 817–819 (1998). [CrossRef] [PubMed]
J. Frank, J. Matrisch, J. Horstmann, S. Altmeyer, and G. Wernicke, “Refractive index determination of transparent samples by noniterative phase retrieval,” Appl. Opt.50(4), 427–433 (2011). [CrossRef] [PubMed]
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]
D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett.80(12), 2586–2589 (1998). [CrossRef]
D. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent, “Quantitative phase-amplitude microscopy. III. The effects of noise,” J. Microsc.214(1), 51–61 (2004). [CrossRef] [PubMed]
L. Waller, L. Tian, and G. Barbastathis, “Transport of intensity imaging with higher order derivatives,” Opt. Express18(12), 12552–12561 (2010). [CrossRef] [PubMed]
W. X. Cong and G. Wang, “Higher-order phase shift reconstruction approach,” Med. Phys.37(10), 5238–5242 (2010). [CrossRef] [PubMed]
B. Xue, S. Zheng, L. Cui, X. Bai, and F. Zhou, “Transport of intensity phase imaging from multiple intensities measured in unequally-spaced planes,” Opt. Express19(21), 20244–20250 (2011). [CrossRef] [PubMed]
L. Waller, M. Tsang, S. Ponda, S. Y. Yang, and G. Barbastathis, “Phase and amplitude imaging from noisy images by Kalman filtering,” Opt. Express19(3), 2805–2814 (2011). [CrossRef] [PubMed]
M. Soto and E. Acosta, “Improved phase imaging from intensity measurements in multiple planes,” Appl. Opt.46(33), 7978–7981 (2007). [CrossRef] [PubMed]
J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996), pp. 55–61.

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper