## Abstract

A numerical reconstruction technique of digital holography based on angular spectrum diffraction by means of the ridge of Gabor wavelet transform (GWT) is presented. Appling the GWT, the object wave can be reconstructed by calculating the wavelet coefficients of the hologram at the ridge of the GWT automatically even if the spectrum of the virtual image is disturbed by the other spectrum. It provides a way to eliminate the effect of the zero-order and the twin-image terms without the spatial filtering. In particular, based on the angular spectrum theory, GWT is applied to the digital holographic phase-contrast microscopy on biological specimens. The theory, the results of a simulation and an experiment of an onion specimen are shown.

©2008 Optical Society of America

## 1. Introduction

Digital holography is applied to the analysis of the phase object because of its power of obtaining quantitative phase information from the hologram, such as the holographic phase-contrast microscopy[1–5]. A charge coupled device (CCD) camera is used to record a hologram onto a computer and numerical methods are subsequently applied to reconstruct the holographic image to enable direct access to both phase and amplitude information of the object wave in digital holography. The most popular numerical reconstruction methods include the well-known Fresnel diffraction integral method[6] and the angular spectrum method[7]. In the Fresnel diffraction integral method, because of the approximation of spherical Huygens wavelet by a parabolic surface, the calculation of the Fresnel diffraction integral can employ a Fourier transform. In order to filter out the zero-order term, the twin image term and the parasitic interferences, the process of the spatial filtering must be carried out. In the angular spectrum method, the angular spectrum of the object wave at the hologram plane is obtained by taking the spatial filtering and the object wave at any plane perpendicular propagation axis is reconstructed by taking the inverse Fourier transform. Therefore, it shows that the Fourier transform analysis and the spatial filtering are very important in these methods. In the in-line holography, applying phase-shifting technique[8] makes it possible to get rid of the twin image without the spatial filtering. But the phase-shifting technique demand more than one hologram, it is not fit for the dynamic analysis.

There are two main determining factors limit the quality of the reconstructed image by the process of the spatial filtering. First is the separation status of the spectrum of the virtual image from the other spectrum. Second are the shape and the size of the spatial filter selected. It is well known that the Fourier transform is a global operation with a poor spatial localization. Therefore, when there are some noises and parasitic interferences introduced into the hologram, the spectrum of the virtual image would be disturbed by the other spectrum. It brings difficulties to define the spatial filter because of the blurry boundary and non-regular distribution of the spectrum. Manual spatial filter should be made according to the human vision every time to obtain the correct spectrum of the virtual image[9]. However, in order to achieve the dynamic live cell analysis[4,9], one of the significant advantages of the digital holography, lots of holograms should be recorded, defining different manual spatial filters would consume plenty of time. Therefore, it is necessary to grope for a numerical reconstruction technique that can be performed automatically.

The wavelet transform, a tool excelling for its multiresolution and localization in the time- or space-frequency domains, is applied to the digital holographic phase-contrast microscopy in this paper. Liebling[10] et al. proposed a multiresolution wavelet method named Fresnelets for the digital holography analysis. Zhong and Weng[11,12] have presented an approach for the analysis of the spatial carrier-fringe patterns by means of the ridge of the wavelet transform for the three-dimension profilometry. Here the Gabor wavelet transform (GWT) is applied to the analysis of the digital hologram. By calculating the wavelet coefficients of the hologram at the ridge automatically, the object wave can be reconstructed. At the same time the effect of the zero-order and the twin-image terms are eliminated without the spatial filtering. Because the Fresnel approximation condition is not satisfied in holographic phase-contrast microscopy, the numerical reconstruction method here we present is based on the angular spectrum diffraction. A digital simulation and an experimental of an onion specimen are presented to demonstrate the validity of the principle.

## 2. Principle of the reconstruction technique employing the GWT

*A. GWT*

The continuous one-dimensional(1-D) wavelet transform(WT) is defined as follows:

where *a*>0 is the scale parameter related to the frequency; *b* is the shift parameter related to the position; *f(x)*is the signal to be analyzed;, *ψ _{a,b}(x)* is the analyzing wavelet obtained by shifting and scaling the mother wavelet

*ψ (x)*; * indicates the complex conjugate.

The Gabor wavelet[13] is employed as the mother wavelet. The Gabor wavelet function and its Fourier transform are given as:

where $\gamma =\pi \sqrt{2\u2044\mathrm{ln}2}$. The Gabor wavelet may be considered as a Gaussian window centered at *x*=0 and its Fourier transform centered at frequency *ω=2π*, as shown in Fig. 1. The function, *ψ _{a,b}(x)* is then centered around

*x=b*with a full width of half maximum Δ

*x*=2

*a*, and its Fourier transform is centered around

*ω*=2

*π/a*. The scale parameter

*a*controls the spatial and frequency resolution of the wavelet decomposition: A large

*a*corresponds to a long wave (high frequency resolution and low spatial resolution), where as a small

*a*corresponds to a short wave (low frequency resolution and high spatial resolution). The WT coefficients

*W(a,b)*using the Gabor wavelet describes the distribution of the signal

*f(x)*on the space-frequency plane. Namely, it can be considered that each of the

*W(a,b)*contains the information about the local frequency (

*ω=2π/a*) component of the signal at the position

*x=b*.

The reasons for the choice of the Gabor wavelet are as follows: (1) The Gabor wavelet consists of a sine wave and a cosine wave that have been modulated by a Gaussian function. (2) It can represent a complex sinusoidal function that correlates well with the sinusoidal characteristics of the interferogram. (3) The Gaussian function has the least spread in both domains of space and frequency.

Computing the 1-D WT, the modulus of the wavelet coefficients can be obtained by:

where *imag(W(a,b)*) and *real(W(a,b))* indicate the imaginary part and the real part of the *W(a,b)*. At the position *b*, the ridge of the wavelet transform[14] is defined as the location where the modulus |*W(a,b)*| reaches its local maximum along the scaling direction *a*.

*B. Reconstruction principle by means of the GWT*

The hologram is created by the interference, in off-axis geometry, between two coherent waves: on one side the wave of interest, called object wave *O(x, y)*, coming from the object, and on the other side a reference wave *R(x, y)* being plane:

where (*x, y*) are the coordinates of the hologram plane; *o(x, y)* and *R*
_{0} are the amplitude of the object and the reference waves; *ϕ(x, y)* is the phase of the object wave; *λ *is the wavelength; *θ* is the angle between the propagation direction of the object and the reference waves. The intensity of the hologram *I (x, y)* can be written as:

$$={\mid {R}_{0}\mid}^{2}+{\mid o(x,y)\mid}^{2}+o(x,y){R}_{0}\mathrm{exp}\left\{j\left[-2\pi \frac{\mathrm{sin}\theta}{\lambda}x+\varphi (x,y)\right]\right\}$$

$$+o(x,y){R}_{0}\mathrm{exp}\{-j\left[-2\pi \frac{\mathrm{sin}\theta}{\lambda}x+\varphi (x,y)\right]\}$$

The first and the second terms form the zero-order term, the third and the fourth terms are the virtual and the real image terms respectively. Take one row of the hologram to analyze. The 1-D intensity distribution of the hologram can be written as:

where *A(x)=|R _{0}|^{2}+|o(x)|^{2}* and $\phi \left(x\right)=-2\pi \frac{\mathrm{sin}\theta}{\lambda}x+\varphi \left(x\right)$. The phase

*φ(x)*can be expanded to a Taylor series around

*b*:

If the higher order terms of (*x-b*)can be neglected, the phase *φ(x)*can be expressed as:

The wavelet transform of *l(x)* becomes:

Setting *A(x)≈A, o(x)≈o*, and performing the above calculation by employing the Gabor wavelet as follows:

Because of the feature of *a*>0, the modulus |*W(a,b)*| reaches its maximum at:

Therefore, the wavelet coefficients at the ridge of the GWT, which is described as *W _{ridge}(b)*, become:

$$+\frac{\sqrt{\gamma}}{4\sqrt[4]{\pi}}{R}_{0}o\mathrm{exp}(-4{\gamma}^{2})\mathrm{exp}\{-j\left[-2\pi \frac{\mathrm{sin}\theta}{\lambda}b+\varphi \left(b\right)\right]\}$$

Here exp(-*γ ^{2}*)=exp[-

*π*

^{2}(2/ln 2)]≈0 and exp(-4

*γ*

^{2})=exp[-4

*π*

^{2}(2/ln 2)]≈0. Note that the position parameter

*b*relates to

*x*, so that Eqs. (14) can be rewritten as:

Multiplying the *W _{ridge}(x)* by an ideal wave corresponding to a replica of the reference wave, the reconstructed wave

*U*at the hologram plane is obtained as:

_{ridge}(x)A conclusion can be obtained from Eqs. (16) that the reconstructed wave at the ridge of the GWT is equal to the object wave at the hologram plane multiplied by a constant coefficient. Therefore the numerical reconstruction of the object wave can be carried out by means of the GWT. At the same time the effect of the zero-order and the twin-image terms is eliminated without the spatial filtering. It is noted that the conclusion is here obtained by assuming the higher order terms of (*x-b*)in Eqs. (9) being neglected. In the following two cases this approximation can be satisfied. First, if *ϕ′(x)* is of slow variation, it can be set that *ϕ″(x)≈ϕ‴(x)*≈…≈0, subsequently the higher order terms of (*x-b*) can be neglected. The phase of the object wave *ϕ(x)* corresponds to the optical path length (OPL) when a plane wave propagates through a phase object. Therefore, the slow variation of *ϕ′(x)* denotes that the first derivative of OPL is of slow variation or the first derivative of refraction-index distribution of the phase object is of slow variation in physics. Second, because the function, *ψ _{a,b}(x)* is concentrated around

*x=b*with a full width of half maximum m Δ

*x=2a*, the higher order terms of (

*x-b*) can be neglected when (

*x-b*)<1 in the full width.

*C. Angular spectrum method*

The angular spectrum reconstruction method has the significant advantage that it has no minimum reconstruction distance requirement as is the case for the Fresnel method. The wave front at the hologram plane is described as *E(x,y;0)*, and its angular spectrum *A(ξ,η;0)*=ℑ{*E(x, y;0)*} can be obtained by taking the Fourier transform. Here ℑ{} denotes the Fourier transform, and *ξ,η* are the corresponding the spatial frequencies of *x,y*. The angular spectrum of the reconstructed wave at the plane *z* perpendicular to the propagation axis, *A(ξ,η;z)*, can be calculated from *A(ξ,η;0)* as:

Subsequently the reconstructed wave can be obtained by taking the inverse Fourier transform as *E′(x′,y′;z)=ℑ-1{A(ξ,η;z)}*, where ℑ-1{} denotes the inverse Fourier transform.

## 3. Simulation

A computer simulation is performed to demonstrate the validity of the GWT algorithm. Fig. 2 (a) shows the amplitude of the phase object of 256×256 pixels employed in the simulation. The amplitude is assumed to be one in the region of the circle, and outside is zero. Fig. 2(b) shows the phase distribution, which can be calculated by:

Here (*m,n*) represents the position in pixel. The phase object is illuminated by a plane wave to produce an object wave, and the reference wave is assumed to be plane. Assuming that the wavelength is 632.8nm, the angle between the propagation direction of the object and the reference waves is 2.1°, and the propagation distance of the object wave is 0.09m. Figure 2(c) is the digital hologram generated by the computer. Figure 2(d) and (e) show the amplitude and the unwrapped phase of the reconstructed wave by employing the Single Fourier Transform Formulation(SFTF)[9]. Figure 2(f) shows the modulus of the wavelet coefficients at the 120^{th} row of the hologram. By detecting the maximum modulus value at each position, the ridge of the GWT can be determined and the corresponding wavelet coefficients can be obtained. Based on the angular spectrum method, the numerical reconstructed object wave is obtained at the same propagation distance. Figure 2(g) and (h) show the amplitude and the unwrapped phase of the reconstructed wave by the analysis of the GWT. Although the calculation time of the GWT algorithm (about 5 seconds) is longer than that of the SFTF algorithm (less than 1 second) obviously by a personal computer, the time consumed by defining the manual spatial filter in the SFTF method would be much longer.

Figure 2(i) shows the simulated phase and the reconstructed phase obtained by the SFTF and the GWT at the 120^{th} row. Figure 2(j) shows the deviation of the reconstructed phase data from the simulated phase data at the 120^{th} row within 0.0~0.08*rad* by the SFTF and that is within -0.08~0.04*rad* by the GWT respectively. The error distribution by the SFTF method is uniform corresponding to the poor spatial localization of the Fourier transform. In the GWT method, the error distribution shows the localization character. Figure 2(k) shows the first and the second derivative of the phase *ϕ(x)* at the 120^{th} row. It can be found that the error distribution curve by the GWT has an analogy with the *ϕ″(x)*curve. It shows that the phase errors analyzed by the GWT are coming from the higher order derivative of *ϕ′(x)*. It is interesting to mention that, although the approximation condition *ϕ″(x)≈ϕ‴(x)*≈…≈0 is not strictly satisfied in this example, the GWT method can still give a fairly accurate reconstruction. This is because the function, *ψ _{a,b}(x)* is concentrated around

*x=b*with a full width of half maximum Δ

*x*=2

*a*, which results in (

*x-b*)<1, the higher order terms of (

*x-b*) in Eqs. (9) can be neglected.

## 4. Experimental result

The digital holography experiment is performed employing the apparatus depicted in Fig. 3. The system is analogous to a Mach-Zehnder interferometer. A He-Ne laser beam (wavelength of 632.8nm) is separated into two beams by a beam splitter BS1. One beam serves as a reference wave (REF), and another beam illuminates the biological specimen (OBJ). A microscope objective MO1 (16x, 0.3NA) collects the object wave transmitted by the specimen and produces a magnified image of the specimen at a distance behind the CCD camera. Another MO2 (16x, 0.3NA) is placed in the reference arm, so that the two beams have matching wave front curvatures. The beam splitter BS2 placed in front of a CCD camera (MINTRON 22K9HC, 795×596 pixels, 8.33×8.33 µm2 per pixel) combines the two beams, and the holograms are recorded by the CCD camera with 8-bit gray scale output. The slight angle is introduced between the object and the reference beams by tilting the beam splitter BS2 for the off-axis holography.

An onion specimen placed between the slide and the cover slide is illuminated by the laser beam. Figure 4(a) shows the image of the onion specimen recorded by the CCD camera when the reference wave is turn off. Fig. 4(b) is the hologram with the onion specimen recorded by the CCD camera of 512×512 pixels. In order to compensate for phase aberrations[15,16] and image distortion[9], a reference hologram without the presence of any specimen is also recorded, as shown in Fig. 4(c). Figure 4(d) is the spectrum of the hologram with the specimen on the logarithmic coordinates. At first, the SFTF is applied to the analysis of the hologram. Considering the effect of the shape and the size of the spatial filter on the reconstructed wave, the hologram is numerically filtered with two square filters and a manual filter as shown in Fig. 4(d), where the small black box represents the square filter of size 15 pixels, the big black box represents the square filter of size 30 pixels and the white line window represents the manual spatial filter.

Figure 5(a) and (b) show the amplitude and the phase of the reconstructed wave by applying the SFTF with the square filter of size 15 pixels. Figure 5(c) and (d) show the results by employing the square filter of size 30 pixels. Figure 5(e) and (f) are that employing the manual filter. It shows that the amplitude and the phase are smooth when the square filter of size 15 pixels is employed, and lots of the high frequency information of the object is lost. When the square filter of size 30 pixels is employed, some high frequency information of the object can be reserved, whereas there is some rugged noises information introduced into the reconstructed wave. It can be concluded that the quality of the reconstructed wave varies with different spatial filters employed. The results by employing the manual spatial filter would be better, but it is not fit for the dynamic analysis.

Figure 6(a) shows the modulus of the wavelet coefficients at the 256^{th} row of the hologram. By detecting the ridge of the GWT and multiplying the wavelet coefficients at the ridge of the GWT by the numerical reference wave, the object wave is reconstructed. Figure 6(b) shows the spectrum of the reconstructed wave by applying the GWT. It shows that the high frequency information of the object wave is reserved, which might be filtered out by the spatial filter employed in the SFTF. Figures 6(c) and (d) show the amplitude and the phase modulo 2*π*, coded to 256 gray levels, of the reconstructed wave. The amplitude image shown in Fig. 6(c) seems to suffer from horizontal “ripples”. This “ripples” is likely to be an additional interference signal produced by the reflected wave on the slides, which appear in the image of the onion specimen shown in Fig. 4(a). Because 1-D WT is also carried out on the horizontal direction, it fails to eliminate this parasitic interference signal. It might be eliminated by employing a two-dimensional WT. After removal of the 2*π* ambiguity by a phase unwrapping process, the unwrapped phase is obtained, as shown in Fig. 6(e).

## 5. Conclusion

In this paper, the GWT algorithm is presented for the analysis of the digital holography to obtain the quantitative phase especially for the phase object. The theory of the GWT technique has been demonstrated in detail. The results of the digital simulation of a phase object and the experiment of an onion specimen show the conclusion that the numerical reconstruction of the phase object can be carried out by the ridge of the GWT. At the same time the effect of the zero-order term and the twin-image term are eliminated without the spatial filtering. One of the significant advantages of the technique here we presented is that it can be applied to the dynamic live biological specimen analysis in the digital holographic phase-contrast microscopy.

This research was supported by Chinese Natural Science Foundation under the contract of 60677019.

## References and links

**1. **D. Carl, B. Kemper, G. Wernicke, and G. von Bally, “Parameter-optimized digital holographic microscope for highresolution living-cell analysis,” Appl. Opt. **43**, 6536–6544 (2004). [CrossRef]

**2. **P. Marquet, B. Rappaz, P. J. Magistretti, E. Cuche, Y. Emery, T. Colomb, and C. Depeursinge, “Digital holographic microscopy: a noninvasive contrast imaging technique allowing quantitative visualization of living cells with subwavelength axial accuracy,” Opt. Lett. **30**, 468–470 (2005). [CrossRef] [PubMed]

**3. **C. J. Mann, L. Yu, C. Lo, and M. K. Kim, “High-resolution quantitative phase-contrast microscopy by digital holography,” Opt. Express **13**, 8693–8698 (2005). [CrossRef] [PubMed]

**4. **B. Kemper and G. von Bally, “Digital holographic microscopy for live cell applications and technical inspection,” Appl. Opt. **47**, A52–A61 (2008). [CrossRef] [PubMed]

**5. **P. Langehanenberg, B. Kemper, D. Dirksen, and G. von Bally, “Autofocusing in digital holographic phase contrast microscopy on pure phase objects for live cell imaging,” Appl. Opt. **47**, 176–182 (2008). [CrossRef]

**6. **J. W. Goodman, *Introduction to Fourier Optics* (McGraw-Hill,1996).

**7. **L. Yu and M. K. Kim, “Wavelength-scanning digital interference holography for tomographic three-dimensional imaging by use of the angular spectrum method,” Opt. Lett. **30**, 2092–2094 (2005). [CrossRef] [PubMed]

**8. **I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. **22**, 1268–1270 (1997). [CrossRef] [PubMed]

**9. **T. Colomb, J. Kühn, F. Charrière, and C. Depeursinge, “Total aberrations compensation in digital holographic microscopy with a reference conjugated hologram,” Opt. Express **14**, 4300–4306 (2006). [CrossRef] [PubMed]

**10. **M. Liebling, T. Blu, and M. Unser, “Fresnelets: new multiresolution wavelet bases for digital holography,” IEEE Trans. Image Proc. **11**, 1–14 (2002).

**11. **J. Zhong and J. Weng, “Spatial carrier-fringe pattern analysis by means of wavelet transform: wavelet transform profilomatry,” Appl. Opt. **43**, 4993–4998 (2004). [CrossRef] [PubMed]

**12. **J. Zhong and J. Weng, “Phase retrieval of optical fringe pattern form the ridge of a wavelet transform,” Opt. Lett. **30**, 2560–2562 (2005). [CrossRef] [PubMed]

**13. **H. Jeong, “Analysis of plate wave propagation in anisotropic laminates using a wavelet transform,” NDT E Inter. **34**, 185–190(2001). [CrossRef]

**14. **A. Cesar and K. Taeeeui, “Determination of strains from fringe patterns using space-frequency representations,” Opt. Eng. **42**, 3182–3193 (2003). [CrossRef]

**15. **F. Montfort, F. Charrière, and T. Colomb, “Purely numerical compensation for microscope objective phase curvature in digital holographic microscopy: influence of digital phase mask position,” J. Opt. Soc. Am. A **23**, 2944–2953 (2006). [CrossRef]

**16. **T. Colomb, E. Cuche, F. Charrière, J. Kühn, N. Aspert, F. Montfort, P. Marquet, and C. Depeursinge, “Automatic procedure for aberration compensation in digital holographic microscopy and applications to specimen shape compensation,” Appl. Opt. **45**, 851–863 (2006). [CrossRef] [PubMed]