Abstract

For applications in the domain of digital holographic microscopy, we present a fast algorithm to propagate scalar wave fields from a small source area to an extended, parallel target area of coarser sampling pitch, using the first Rayleigh-Sommerfeld diffraction formula. Our algorithm can take full advantage of the fast Fourier transform by decomposing the convolution kernel of the propagation into several convolution kernel patches. Using partial overlapping of the patches together with a soft blending function, the Fourier spectrum of these patches can be reduced to a low number of significant components, which can be stored in a compact sparse array structure. This allows for rapid evaluation of the partial convolution results by skipping over negligible components through the Fourier domain pointwise multiplication and direct mapping of the remaining multiplication results into a Fourier domain representation of the coarsly sampled target patch. The algorithm has been verified experimentally at a numerical aperture of 0.62, not showing any significant resolution limitations.

© 2010 Optical Society of America

Full Article  |  PDF Article

References

  • View by:
  • |
  • |
  • |

  1. M. Gu, Advanced Optical Imaging Theory, vol. 75 of Optical Sciences (Springer, Berlin, Germany, 1999).
  2. P. P. Vaidyanathan, "Generalizations of the Sampling Theorem: Seven Decades After Nyquist," IEEE Trans. Circ. Syst. I Fundam. Theory Appl. 48, 1094-1109 (2001).
    [CrossRef]
  3. F. Shen, and A. Wang, "Fast-Fourier-transform based numerical integration method for the Rayleigh-Sommerfeld diffraction formula," Appl. Opt. 45(6), 1102-1110 (2006).
    [CrossRef] [PubMed]
  4. D. Gabor, "A new microscopic principle," Nature 161(4098), 777 (1948).
    [CrossRef] [PubMed]
  5. U. Schnars, "Direct phase determination in hologram interferometry with use of digitally recorded holograms," J. Opt. Soc. Am. A 11(7), 2011-2015 (1994).
    [CrossRef]
  6. H. J. Kreuzer, "Low energy electron point source microscopy," Micron 26(6), 503-509 (1995).
    [CrossRef]
  7. H. J. Kreuzer, "Holographic microscope and method of hologram reconstruction," US 6,411,406 B1, June 25, 2002 (Canadian Patent CA2376395).
  8. R. W. Gerchberg, and W. O. Saxton, "A practical algorithm for the determination of phase from image and diffraction plane pictures," Optik (Stuttg.) 35, 227-246 (1972).
  9. G.-Z. Yang, B.-Z. Dong, B.-Y. Gu, J.-Y. Zhuang, and O. K. Ersoy, "Gerchberg-Saxton and Yang-Gu algorithms for phase retrieval in a nonunitary transform system: a comparison," Appl. Opt. 33(2), 209-218 (1994).
    [CrossRef] [PubMed]
  10. A. Grjasnow, A. Wuttig, and R. Riesenberg, "Phase resolving microscopy by multi-plane diffraction detection," J. Microsc. 231(1), 115-123 (2008).
    [CrossRef] [PubMed]
  11. G. C. Sherman, "Application of the Convolution Theorem to Rayleigh’s Integral Formulas," J. Opt. Soc. Am. 57(4), 546-547 (1967) (Proof: Rayleigh-Sommerfeld-integral ) (CV) ( equals angular spectrum ) (AS).
    [CrossRef] [PubMed]
  12. J. Li, Z. Peng, and Y. Fu, "Diffraction transfer function and its calculation of classic diffraction formula," Opt. Commun. 280(2), 243-248 (2007).
    [CrossRef]
  13. J. W. Goodman, Introduction to Fourier Optics, Electrical and Computer Engineering, 2nd ed. (McGraw-Hill Companies, Inc., USA, 1996).
  14. M. Kanka, R. Riesenberg, and H. J. Kreuzer, "Reconstruction of high-resolution holographic microscopic images," Opt. Lett. 34(8), 1162-1164 (2009).
    [CrossRef] [PubMed]
  15. M. Kanka, A. Wuttig, C. Graulig, and R. Riesenberg, "Fast exact scalar propagation for an in-line holographic microscopy on the diffraction limit," Opt. Lett. 35(2), 217-219 (2010).
    [CrossRef] [PubMed]
  16. V. Nascov, and P. C. Logofătu, "Fast computation algorithm for the Rayleigh-Sommerfeld diffraction formula using a type of scaled convolution," Appl. Opt. 48(22), 4310-4319 (2009).
    [CrossRef] [PubMed]
  17. M. Paturzo, P. Memmolo, L. Miccio, A. Finizio, P. Ferraro, A. Tulino, and B. Javidi, "Numerical multiplexing and demultiplexing of digital holographic information for remote reconstruction in amplitude and phase," Opt. Lett. 33(22), 2629-2631 (2008).
    [CrossRef] [PubMed]
  18. M. Frigo, and S. G. Johnson, "The Design and Implementation of FFTW3," in Proc. IEEE, vol. 93, pp. 216-231 (2005).
    [CrossRef]

2010

2009

2008

2007

J. Li, Z. Peng, and Y. Fu, "Diffraction transfer function and its calculation of classic diffraction formula," Opt. Commun. 280(2), 243-248 (2007).
[CrossRef]

2006

2001

P. P. Vaidyanathan, "Generalizations of the Sampling Theorem: Seven Decades After Nyquist," IEEE Trans. Circ. Syst. I Fundam. Theory Appl. 48, 1094-1109 (2001).
[CrossRef]

1995

H. J. Kreuzer, "Low energy electron point source microscopy," Micron 26(6), 503-509 (1995).
[CrossRef]

1994

1972

R. W. Gerchberg, and W. O. Saxton, "A practical algorithm for the determination of phase from image and diffraction plane pictures," Optik (Stuttg.) 35, 227-246 (1972).

1967

1948

D. Gabor, "A new microscopic principle," Nature 161(4098), 777 (1948).
[CrossRef] [PubMed]

Dong, B.-Z.

Ersoy, O. K.

Ferraro, P.

Finizio, A.

Fu, Y.

J. Li, Z. Peng, and Y. Fu, "Diffraction transfer function and its calculation of classic diffraction formula," Opt. Commun. 280(2), 243-248 (2007).
[CrossRef]

Gabor, D.

D. Gabor, "A new microscopic principle," Nature 161(4098), 777 (1948).
[CrossRef] [PubMed]

Gerchberg, R. W.

R. W. Gerchberg, and W. O. Saxton, "A practical algorithm for the determination of phase from image and diffraction plane pictures," Optik (Stuttg.) 35, 227-246 (1972).

Graulig, C.

Grjasnow, A.

A. Grjasnow, A. Wuttig, and R. Riesenberg, "Phase resolving microscopy by multi-plane diffraction detection," J. Microsc. 231(1), 115-123 (2008).
[CrossRef] [PubMed]

Gu, B.-Y.

Javidi, B.

Kanka, M.

Kreuzer, H. J.

Li, J.

J. Li, Z. Peng, and Y. Fu, "Diffraction transfer function and its calculation of classic diffraction formula," Opt. Commun. 280(2), 243-248 (2007).
[CrossRef]

Logofatu, P. C.

Memmolo, P.

Miccio, L.

Nascov, V.

Paturzo, M.

Peng, Z.

J. Li, Z. Peng, and Y. Fu, "Diffraction transfer function and its calculation of classic diffraction formula," Opt. Commun. 280(2), 243-248 (2007).
[CrossRef]

Riesenberg, R.

Saxton, W. O.

R. W. Gerchberg, and W. O. Saxton, "A practical algorithm for the determination of phase from image and diffraction plane pictures," Optik (Stuttg.) 35, 227-246 (1972).

Schnars, U.

Shen, F.

Sherman, G. C.

Tulino, A.

Vaidyanathan, P. P.

P. P. Vaidyanathan, "Generalizations of the Sampling Theorem: Seven Decades After Nyquist," IEEE Trans. Circ. Syst. I Fundam. Theory Appl. 48, 1094-1109 (2001).
[CrossRef]

Wang, A.

Wuttig, A.

Yang, G.-Z.

Zhuang, J.-Y.

Appl. Opt.

IEEE Trans. Circ. Syst. I Fundam. Theory Appl.

P. P. Vaidyanathan, "Generalizations of the Sampling Theorem: Seven Decades After Nyquist," IEEE Trans. Circ. Syst. I Fundam. Theory Appl. 48, 1094-1109 (2001).
[CrossRef]

J. Microsc.

A. Grjasnow, A. Wuttig, and R. Riesenberg, "Phase resolving microscopy by multi-plane diffraction detection," J. Microsc. 231(1), 115-123 (2008).
[CrossRef] [PubMed]

J. Opt. Soc. Am.

J. Opt. Soc. Am. A

Micron

H. J. Kreuzer, "Low energy electron point source microscopy," Micron 26(6), 503-509 (1995).
[CrossRef]

Nature

D. Gabor, "A new microscopic principle," Nature 161(4098), 777 (1948).
[CrossRef] [PubMed]

Opt. Commun.

J. Li, Z. Peng, and Y. Fu, "Diffraction transfer function and its calculation of classic diffraction formula," Opt. Commun. 280(2), 243-248 (2007).
[CrossRef]

Opt. Lett.

Optik (Stuttg.)

R. W. Gerchberg, and W. O. Saxton, "A practical algorithm for the determination of phase from image and diffraction plane pictures," Optik (Stuttg.) 35, 227-246 (1972).

Other

M. Frigo, and S. G. Johnson, "The Design and Implementation of FFTW3," in Proc. IEEE, vol. 93, pp. 216-231 (2005).
[CrossRef]

H. J. Kreuzer, "Holographic microscope and method of hologram reconstruction," US 6,411,406 B1, June 25, 2002 (Canadian Patent CA2376395).

J. W. Goodman, Introduction to Fourier Optics, Electrical and Computer Engineering, 2nd ed. (McGraw-Hill Companies, Inc., USA, 1996).

M. Gu, Advanced Optical Imaging Theory, vol. 75 of Optical Sciences (Springer, Berlin, Germany, 1999).

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 (5)

Fig. 1
Fig. 1

The main aim of the algorithm is the computational scalar wave field propagation of a small, densely sampled source field to a large, coarsely sampled detector image using an unaltered Rayleigh-Sommerfeld kernel. The sampling of the detector image will here differ from the sampling of the source field by an integer factor T. The propagation is divided into several patches with overlapping target regions on the detector.

Fig. 2
Fig. 2

Geometrical arrangement of propagator kernel patches around an arbitrarily chosen kernel patch hμν; Different patches are arranged with a pitch of SN pixels, leading to an overlap between adjacent patches of w = NS pixels. The sum of all patches gives the total propagator kernel h. The blending function shown here is linked to the particular decomposition of h and will be described in the next section. The bitmap overlay for hμν is figurative.

Fig. 3
Fig. 3

Effect of soft blending of a propagator kernel patch; The left column shows the spatial amplitude distribution of an arbitrarily chosen (off-axis) propagator kernel patch without (a) and with (b) application of a soft blending function. Uncolored (grey and white) parts mean a phase angle of 0 or π, while red and blue color overlays stand for phase angles of π/2 and (3/2)π, respectively. The same color encoding is used for the Fourier domain representation of the patches shown on the right side. The figure visualizes the strong reduction of relevant Fourier components through the application of a soft-blending.

Fig. 4
Fig. 4

Schematic flow chart of the patch propagation algorithm; The point-wise multiplication and Fourier domain resampling are carried out in the packed domain, causing the main acceleration of the algorithm.

Fig. 5
Fig. 5

The validity of the patched propagation was investigated by propagating an object wave of a real inline holography sample containing clusters of 1.06 μm polymer beads (a) to a detector using the described algorithm and then back to the object plane (b), using the Fourier domain tile superposition propagation algorithm [15]. The images show a 40 × 76 μm2 detail of the object plane. The distance between object and detector plane was about 3.8 mm. A truncation threshold of 10−8 was used for the forward propagation. No perceivable differences can be seen, indicating that the full information of the sample is retained through the forward-backward iteration cycle. This is verified by the 5000 times intensified difference intensity image (c), showing only minor, non-systematic differences. The wavelength of the measurement was 661 nm, indicating a realized numerical aperture of at least 0.62. Image (d) shows the 2 billion times intensified difference between the round-tripped image (b) and a round-tripped image generated without truncation in the forward propagation.

Tables (1)

Tables Icon

Table 1 Some corner marks for propagation times, average propagator fill factors, and total memory requirements for different detector sizes, propagation distances and patch sizes (N). All other parameters were held constant: wavelength λ =500 nm; sampling raster Δx=125nm; detector pixel pitch T Δx=4 μm; overlap between the patches: N/8; truncation threshold for propagator Fourier components: 10−8. Computer: Intel Motherboard DP35DP with Processor Core Duo E6750 at 2.66 GHz, 2GiB RAM. FFT using the FFTW library [18] and float (32 bit) arithmetics; Propagator patches were pre-computed using double precision (64 bit) arithmetics; Compiler: MSVC++ 2005; The given memory requirements include also the memory for the input and output fields. Propagator patches were reused up to eight times, when linked to other patches through a mirroring operation. No parallelization of the patch propagation was done

Equations (17)

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

u ( x , y ) = d x d y h ( x x , y y , z ) o ( x , y ) ,
o m , n = o ( x 0 + Δ x m , y 0 + Δ y n )
u p , q = u ( x 0 + Δ x p , y 0 + Δ y q ) ,
u p , q = m n h ( p Δ x m Δ x , q Δ y n Δ y , z ) o m , n .
u p , q = m n h p m , q n o m , n { h * o } p , q ,
h p , q = Δ x Δ y 2 π z r p , q 2 ( 1 r p , q i k ) exp ( i k r p , q z | z | ) ,
r p , q = [ ( Δ x p ) 2 + ( Δ y q ) 2 + z 2 ] 1 / 2 ,
Δ x Δ y ,
u = 1 [ ( h ) ( o ) 1 [ H O ] .
[ ( u ) ] s , t U s , t = p = 0 N 1 q = 0 N 1 u p , q exp ( 2 π i s p + t q N ) .
h p , q = μ ν h p S μ , q S ν μ ν , ( 0 { p S μ , q S ν } < N ) ,
u p , q = μ ν ( m = 0 N 1 n = 0 N 1 h p S μ m , q S ν n μ ν o m , n ) = μ ν { h μ ν * o } p S μ , q S ν .
b ( d ) = sin 2 ( π 2 d w ) ,
u p = 1 2 N t = 0 2 N 1 U t exp ( 2 π i t p 2 N ) .
u ^ j u T j + s = 1 2 N t = 0 2 N 1 U t exp ( 2 π i t ( T j + s ) 2 N ) = 1 2 N l = 0 2 N 1 U ^ l e x p ( 2 π i j l 2 N ) ,
U ^ l 1 T m = 0 T 1 U l + 2 N m exp ( 2 π i s ( l + 2 N m ) 2 N ) .
u ^ j = [ 1 U ^ ] j

Metrics