Abstract

It is an effective scheme to the phase retrieval for axial intensity derivative computing. In this paper, we demonstrate a method for estimating the axial intensity derivative and improving the calculation accuracy in the transport of intensity equation (TIE) from multiple intensity measurements. The method takes both the higher-order intensity derivatives and the noise into account, and minimizes the impact of detecting noise. The simulation results demonstrate that the proposed method can effectively reduce the error of intensity derivative computing.

©2012 Optical Society of America

1. Introduction

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

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

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

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

2. Multiple intensity measurements for phase retrieval

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

u(r)=[I(r)]1/2exp[iϕ(r)],
where r=(x,y) denotes the coordinate in the x-y plane, perpendicular to the z axis, I(r) and ϕ(r) are the intensity and phase distribution in the plane respectively. In this paper, the focal plane is assumed to be located atz=0.

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

kIz=(Iϕ),
where k=2π/λ is the wave number, and λ denotes the wavelength, ={x,y} is the two dimensional gradient operator in the x-y plane.

Assume the intensity distribution of m x-y planes located at z1,z2,zm is available, based on the Maclaurin formula, I(zi) (ignore the x-y dimension) can be expanded as

I(zi)=I(0)+I(0)zzi+k=2nzikk!kI(0)zk+zin+1(n+1)!n+1I(ξi)zn+1,
where kzkdenotes the kth order derivative, and ξi is some number located between 0 and zi. As a result, we get m equations. Multiply these equations with a number ai separately and sum their left and right hand side respectively, we can get the following expression

i=1maiI(zi)=I(0)i=1mai+I(0)zi=1maizi+k=2n1k!kI(0)zki=1maizik+1(n+1)!i=1maizin+1n+1I(ξi)zn+1.

Let

i=1mai=0,i=1maizi=1,i=1maizik=0k=2,3...n,
we can derive the derivative of intensity

I(0)z=i=1maiI(zi)1(n+1)!i=1maizin+1n+1I(ξi)zn+1.

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

[zI]c=i=1maiI(zi).

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

[zI]cI(0)z=i=1maiN(zi)+1(n+1)!i=1maizin+1n+1I(ξi)zn+1,
where N(zi) denotes the intensity noise at zi.Under the assumption of noise statistic independence, the mean square can be written as
ε2=D(N)+D(ξ),
where
D(N)=σ2i=1mai2,
D(ξ)=[1(n+1)!i=1maizin+1n+1I(ξi)zn+1]2,
and σ2denotes the variance of the noise. Because of the uncertainty of high-order item, we only evaluate the impact of detecting noise i.e. Equation (10). Our mission is to minimize it under the condition of Eq. (5).

Generally, Eq. (5) has another description form

[1111z1z2z3zmz12z22z32zm2z1nz2nz3nzmn]·[a1a2a3am]=[0100].

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

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

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

[AB]·[CD]=E,
where, B consists of the last n + 1 columns of the transformation matrix, which is nonsingular. AndC, D is divided to match A and B.

Then D can be written as

D=B-1EB-1AC.

As a result,

D(N)=σ2[i=1n+1ai2+j=1mn1(fjGjC)2].

To minimize it,

akD(N)=σ2[2akj=1n+12(fjGjC)gjk]=0k=1,2(mn1),
where fj denotes the jth element of the vector F=B-1E, and Gj is jth row vector of matrix G=B-1A. Obviously, Eq. (16) is equivalent to

(GTG+I)C=GTF.

Combining with Eq. (13), we get

[GTG+I0AB]·[CD]=[GTFE].

In Eq. (18), GTG+I and B are both nonsingular matrix, so the coefficient a1~amcan be derived easily.

3. Simulations

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

That is the 0th order Gauss mode

u00(x,y,z)=cw(z)exp(r2w2(z))exp{i[k(z+r22R)arctan(zf)]},
where

k=2πλr2=x2+y2w(z)=w01+(zf)2R=R(z)=z[1+(fz)2]f=πw02λw0=λfπ.

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

I(z)=c2w2(z),
where

w0=0.005m,λ=1550nm,f=πw02λ=50.6708m,c=0.01.

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

 figure: Fig. 1

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

Download Full Size | PPT Slide | PDF

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

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

 figure: Fig. 2

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

Download Full Size | PPT Slide | PDF

4. Conclusions

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

Acknowledgments

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

References and links

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

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

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

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

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

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

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

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

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

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

References

  • View by:
  • |
  • |
  • |

  1. M. Reed Teague, “Deterministic phase retrieval: a Green's function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983).
    [Crossref]
  2. D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. 80(12), 2586–2589 (1998).
    [Crossref]
  3. T. E. Gureyev and K. A. Nugent, “Rapid quantitative phase imaging using the transport of intensity equation,” Opt. Commun. 133(1-6), 339–346 (1997).
    [Crossref]
  4. L. J. Allen and M. P. Oxley, “Phase retrieval from series of images obtained by defocus variation,” Opt. Commun. 199(1-4), 65–75 (2001).
    [Crossref]
  5. T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation. II. orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A 13(8), 1670–1682 (1996).
    [Crossref]
  6. T. E. Gureyev, A. Roberts, and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation: matrix solution with use of Zernike polynomials,” J. Opt. Soc. Am. A 12(9), 1932–1942 (1995).
    [Crossref]
  7. M. Soto and E. Acosta, “Improved phase imaging from intensity measurements in multiple planes,” Appl. Opt. 46(33), 7978–7981 (2007).
    [Crossref] [PubMed]
  8. L. Waller, L. Tian, and G. Barbastathis, “Transport of intensity imaging with higher order derivatives,” Opt. Express 18(12), 12552–12561 (2010).
    [Crossref] [PubMed]
  9. B. D. Xue, S. L. Zheng, L. Y. Cui, X. Z. Bai, and F. G. Zhou, “Transport of intensity phase imaging from multiple intensities measured in unequally-spaced planes,” Opt. Express 19(21), 20244–20250 (2011).
    [Crossref] [PubMed]
  10. S. L. Zheng, B. D. Xue, W. F. Xue, X. Z. Bai, and F. G. Zhou, “Transport of intensity phase imaging from multiple noisy intensities measured in unequally-spaced planes,” Opt. Express 20(2), 972–985 (2012).
    [Crossref] [PubMed]

2012 (1)

2011 (1)

2010 (1)

2007 (1)

2001 (1)

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]

1998 (1)

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

1997 (1)

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]

1996 (1)

1995 (1)

1983 (1)

Acosta, E.

Allen, L. J.

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]

Bai, X. Z.

Barbastathis, G.

Cui, L. Y.

Gureyev, T. E.

Nugent, K. A.

D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. 80(12), 2586–2589 (1998).
[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]

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

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

Oxley, M. P.

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]

Paganin, D.

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

Reed Teague, M.

Roberts, A.

Soto, M.

Tian, L.

Waller, L.

Xue, B. D.

Xue, W. F.

Zheng, S. L.

Zhou, F. G.

Appl. Opt. (1)

J. Opt. Soc. Am. (1)

J. Opt. Soc. Am. A (2)

Opt. Commun. (2)

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]

Opt. Express (3)

Phys. Rev. Lett. (1)

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

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (2)

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

Equations (22)

Equations on this page are rendered with MathJax. Learn more.

u ( r ) = [ I ( r ) ] 1 / 2 exp [ i ϕ ( r ) ] ,
k I z = ( I ϕ ) ,
I ( z i ) = I ( 0 ) + I ( 0 ) z z i + k = 2 n z i k k ! k I ( 0 ) z k + z i n + 1 ( n + 1 ) ! n + 1 I ( ξ i ) z n + 1 ,
i = 1 m a i I ( z i ) = I ( 0 ) i = 1 m a i + I ( 0 ) z i = 1 m a i z i + k = 2 n 1 k ! k I ( 0 ) z k i = 1 m a i z i k + 1 ( n + 1 ) ! i = 1 m a i z i n + 1 n + 1 I ( ξ i ) z n + 1 .
i = 1 m a i = 0 , i = 1 m a i z i = 1 , i = 1 m a i z i k = 0 k = 2 , 3... n ,
I ( 0 ) z = i = 1 m a i I ( z i ) 1 ( n + 1 ) ! i = 1 m a i z i n + 1 n + 1 I ( ξ i ) z n + 1 .
[ z I ] c = i = 1 m a i I ( z i ) .
[ z I ] c I ( 0 ) z = i = 1 m a i N ( z i ) + 1 ( n + 1 ) ! i = 1 m a i z i n + 1 n + 1 I ( ξ i ) z n + 1 ,
ε 2 = D ( N ) + D ( ξ ) ,
D ( N ) = σ 2 i = 1 m a i 2 ,
D ( ξ ) = [ 1 ( n + 1 ) ! i = 1 m a i z i n + 1 n + 1 I ( ξ i ) z n + 1 ] 2 ,
[ 1 1 1 1 z 1 z 2 z 3 z m z 1 2 z 2 2 z 3 2 z m 2 z 1 n z 2 n z 3 n z m n ] · [ a 1 a 2 a 3 a m ] = [ 0 1 0 0 ] .
[ A B ] · [ C D ] = E ,
D = B - 1 E B - 1 A C .
D ( N ) = σ 2 [ i = 1 n + 1 a i 2 + j = 1 m n 1 ( f j G j C ) 2 ] .
a k D ( N ) = σ 2 [ 2 a k j = 1 n + 1 2 ( f j G j C ) g j k ] = 0 k = 1 , 2 ( m n 1 ) ,
( G T G + I ) C = G T F .
[ G T G + I 0 A B ] · [ C D ] = [ G T F E ] .
u 00 ( x , y , z ) = c w ( z ) exp ( r 2 w 2 ( z ) ) exp { i [ k ( z + r 2 2 R ) arctan ( z f ) ] } ,
k = 2 π λ r 2 = x 2 + y 2 w ( z ) = w 0 1 + ( z f ) 2 R = R ( z ) = z [ 1 + ( f z ) 2 ] f = π w 0 2 λ w 0 = λ f π .
I ( z ) = c 2 w 2 ( z ) ,
w 0 = 0.005 m , λ = 1550 n m , f = π w 0 2 λ = 50.6708 m , c = 0.01.

Metrics