OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 20, Iss. 2 — Jan. 16, 2012
  • pp: 972–985
« 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)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

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.

© 2012 OSA

1. Introduction

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

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

, 2

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

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

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

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

, 5

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

, 19

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

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

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]

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

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

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

]:
·[Iz(r)ϕz(r)]=kIz(r)/z,
(1)
where k=2π/λ is the wave number, =(x,y) is the two-dimensional differential operator in the x-y plane. Iz(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

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 Izjand Izj, with the defocused distances zjand zj, (j=1,2,...,M), are given by
Izj(r)=I¯zj(r)+nzj(r).
(2)
Izj(r)=I¯zj(r)+nzj(r).
(3)
where, I¯zj(r)and I¯zj(r)represent the ideal intensities. nzj, nzjare supposed to be zero-mean, additive Gaussian noise with varianceσ2, independent of the intensities. It should be mentioned that zjand zj, (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 planez=0, i.e.I0(r)=Iz(r)/z|z=0, can be estimated by:

I^0(r)=j=1Majzjzj(Izj(r)Izj(r)).
(4)

It should be mentioned that zj can be negative or just be z0, which means that we can estimate the axial intensity derivative independently from either the double-sided or single-sided diffraction intensities. The coefficients aj,j=1,2,...,M must minimize the mean square error between the estimated and the real intensity derivatives
ε2(r)=E{[I^0(r)I0(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 I¯z(r)at z=0 and with the Taylor theorem, we can get
I¯zj(r)=I¯0(r)+zjI0(r)+zj22Icj(r)/2,cj[0,zj].
(6)
I¯zj(r)=I¯0(r)+zjI0(r)+zj22Icj(r)/2,cj[zj,0].
(7)
where 2Ic±j(r)=[2Iz(r)/z2]z=c±j. Substituting Eq. (6) and Eq. (7) into Eq. (4), we can obtain

I^0(r)=j=1Majzjzj(nzj(r)nzj(r))+I0(r)j=1Maj+j=1Majzjzj(zj222Icj(r)zj222Icj(r)).
(8)

Then, substituting Eq. (8) into Eq. (5) we can get
j=1Maj1=0.
(9)
ε2(r)=εn2(r)+εh2(r)isminimum,
(10)
with
εn2(r)=2σ2j=1Maj2(zjzj)2,
(11)
and
εh2(r)=14{j=1Majzjzj[zj22Icj(r)zj22Icj(r)]}2.
(12)
whereεn2(r)and εh2(r) represent the error caused by noise and the higher order intensity derivatives respectively. When the noise is not considered, i.e. εn2(r)is zero, the coefficientsaj, j=1,2,...,M can be calculated by the method introduced in [18

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]

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

j=1Maj2(zjzj)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

aj=(zjzj)2j=1M(zjzj)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:
·[Iz(r)ϕz(r)]=kj=1Msj(Izj(r)Izj(r)),
(15)
where the weight coefficients are

sj=(zjzj)j=1M(zjzj)2.
(16)

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

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

zj=jΔz,zj=jΔz,j=1,2,...
(17)

The coefficients for the intensity measurements become:
sj=3jM(M+1)(2M+1),j=±1,±2,...,
(18)
which correspond to that in Eq. (9) in [20

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

].

3 Accuracy analysis of the estimated intensity derivative

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

εn2(r)=2σ2/j=1M(zjzj)2,
(19)
εh2(r)=14{j=1Mzjzjj=1M(zjzj)2[zj22Icj(r)zj22Icj(r)]}2.
(20)

As shown in Eq. (19), the noise error becomes smaller with the increase of M. However, εh2(r) is indeterminable as Mincreases. This is because, firstly, the second-order intensity derivatives cannot be precisely known in advance, and secondly, the sign of (zjzj),j=1,2,...,Mis uncertain when zjandzjare arbitrary values. Consequently, we cannot make sure whether εh2(r)is decreasing by increasing the number of the intensity planes. In this study, we will provide a general idea of it by analyzing εh2(r)with some prescribed assumptions.

Let us assume that zj is positive and zjis negative. The lower bound and the upper bound of εh2(r) are adopted to analyze its effects, and they can be derived from Eq. (20) as:
εh(lo)2εh2(r)εh(up)2,
(21)
with
εh(up)2=U2Zup2(M)/4,Zup2(M)=[j=1M(zjzj)(zj2+zj2)]2[j=1M(zjzj)2]2,εh(lo)2=0,
and
U=maxr,c(|2Ic(r)|),c[zmin,zmax],zminzj,zjzmax,j=1,2,...,M.
where [zmin,zmax]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

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

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) as
εM+12(r)<εM2(r)εh_(M+1)2(r)εh_M2(r)<εn_M2(r)εn_(M+1)2(r),
(22)
where εM2(r),εn_M2(r)and εh_M2(r)denote respectively the total error, the noise error and the error caused by the higher intensity derivatives with Mintensity planes. While, εM+12(r),εn_(M+1)2(r)and εh_(M+1)2(r) are those with M+1intensity planes.

To further analyzeε2(r), we here assume εh2(r) approximates to εh(up)2. Then substitutingεh2(r) by εh(up)2so that Eq. (22) can be expressed as
εh(up)_(M+1)2(r)εh(up)_M2(r)<εn_M2(r)εn_(M+1)2(r).
(23)
where εh(up)_M2(r)εh(up)2 and εh(up)_(M+1)2(r)are respectively the upper bounds of εh_M2(r)and εh_(M+1)2(r).

In Eq. (23), if εh(up)_(M+1)2(r)εh(up)_M2(r)<0, i.e. εh(up)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 Zup2(M)decrease. As an example, we can choose the defocused values of the already and newly used intensities as:

{zM+1<zj+1<zj,zM+1|z(j+1)|<|zj|,|z(M+1)|<|z(j+1)|,|z(M+1)|zj+1.(j=1,2,...,M1).
(24)

This means that the newly used intensities with the defocused distances zM+1 and z(M+1) are closer to the focal plane than the already used intensities with defocused values zjand zj, (j=1,2,...,M). Then we have:

(zM+12+z(M+1)2)(zj+1z(j+1))(zj+12+z(j+1)2)(zM+1z(M+1))=zM+1zj+1(zM+1zj+1)zM+1z(j+1)(zM+1+z(j+1))+z(M+1)zj+1(z(M+1)+zj+1)z(M+1)z(j+1)(z(M+1)z(j+1))<0
(25)

Consequently, we can get:

Zup(M+1)Zup(M)=j=1M+1(zjzj)(zj2+zj2)j=1M+1(zjzj)2j=1M(zjzj)(zj2+zj2)j=1M(zjzj)2=j=1M(zM+1z(M+1))(zjzj)[(zM+12+z(M+1)2)(zjzj)(zj2+zj2)(zM+1z(M+1))]j=1M(zjzj)2j=1M+1(zjzj)2<0.
(26)

Therefore, Zup2(M) decreases with Mincreasing.

However, If εh(up)_(M+1)2(r)εh(up)_M2(r)>0, i.e., εh(up)2is an increasing function of M, The inequality in Eq. (23) may not be always true. For example, when the defocused values satisfy
{zM+1>zj+1>zj,zM+1|z(j+1)|>|zj|,|z(M+1)|>|z(j+1)|,|z(M+1)|zj+1.(j=1,2,...,M1),
(27)
Zup2(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 εn2(r) is the dominant factor inε2(r), it can be significantly reduced by using more intensity planes. Whereas as the reducing of εn2(r)and the increasing of εh(up)2, εh(up)2may becomes dominant when the number of the used intensity planes increases to a critical value of Mc. The value of Mcis determined by the defocused values and the noise levels. It cannot be precisely determined in advance.

In summary, εn2(r) can be reduced by using more intensity measurements, and εh2(r)can be reduced when Zup2(M)decreases.

4. Simulations

Numerical simulations are conducted to test the derived formulas for the intensity derivative estimation. Figure 1
Fig. 1 Simulated phase distribution in the z=0plane.
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.

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:
E2=x=1Xy=1Y[I^0(x,y)I0(x,y)]2/x=1Xy=1Y[I0(x,y)]2.
(28)
whereI^0(x,y)and I0(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

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):
I0=2Im[u0*(r)F1{F{u0(r)}k(1λ2g2)1/2}].
(29)
Where F{} and F1{} denote respectively the Fourier and inverse Fourier transforms. g=(gx,gy)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

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

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

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:
RMSE=x=1Xy=1Y[θ0(x,y)θ0(x,y)μ]2/(X×Y),μ=x=1Xy=1Y[θ0(x,y)θ0(x,y)]/(X×Y)
(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μmand29μm.The positive defocused values ofzj are assigned to be29,27,26,23,22,20,19,17,16,12,11,10,8,5μm and3μm. The negative defocus values of zj 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 be4μmand1μm. These defocused values can make Zup(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
Fig. 2 Simulated noisy diffraction intensities in the z=19μmplane with different SNRs of (a)70db; (b)60db; (c)55db; (d)50db.
shows an example of the diffraction images in the z=19μm plane with different SNR.

Figure 3
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.
shows the averaged NMSEs as a function of M with different SNR. Figure 4
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.
shows the retrieved phases with different number of noisy intensity measurements and with different SNR. In Table 1

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.

table-icon
View This Table
| View All Tables
, NMSEs and RMSEs are compared with different number of noisy intensity measurements and with different levels of SNR.

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.

Quantitative comparisons are listed in Table 2

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.

table-icon
View This Table
| View All Tables
. 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
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.
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.

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
Fig. 8 Simulated diffraction intensities corrupted with mixed noise in the z=19μmplane. (a) case 1; (b) case 2; (c) case 3; and (d) case 4.
shows the noisy intensity images in the z=19μmplane. 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
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.
and Fig. 10
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).
. 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.

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.

Appendix A

The propagation of the wave function uz(r) in the z0plane from the z=0plane can be calculated by the angular spectrum method:

uz(r)=F1{F{u0(r)}Hz(g)},Hz(g)=exp(jkz(1λ2g2)1/2).
(A1)

where F{} and F1{} denote respectively the Fourier and inverse Fourier transforms. g=(gx,gy)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:

Iz=z(uz(r)uz*(r))=2Re[uz*(r)zuz(r)]=2Re[uz*(r)F1{F{u0(r)}jk(1λ2g2)1/2Hz(g)}].
(A2)

Let z=0, Eq. (A2) turns to:

I0=2Im[u0*(r)F1{F{u0(r)}k(1λ2g2)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

  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. A12(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,” Micron33(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. A27(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. A27(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. A10(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. Express18(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. Express19(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. Express19(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.

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.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited