Open Access
Add to BibSonomyAbstract
In Digital Holography there are applications where computing a few samples of a wavefield is sufficient to retrieve an image of the region of interest. In such cases, the sampling rate achieved by the direct and the spectral methods of the discrete Fresnel transform could be excessive. A few algorithmic methods have been proposed to numerically compute samples of propagated wavefields while allowing down-sampling control. Nevertheless, all of them require the computation of at least two 2D discrete Fourier transforms which increases the computational load. Here, we propose the use of an aliasing operator and a single discrete Fourier transform to achieve an efficient method to down-sample the wavefields obtained by the Fresnel transform.
©2012 Optical Society of America
1. Introduction
Digital Holography (DH) requires to numerically compute samples of a wavefield propagated from an interference pattern (formed by the light scattered from the objects in a scene and the light of a reference wave) recorded by a digital camera (CCD or CMOS).
The discrete Fresnel transform (DFreT) has been widely used in DH to numerically compute arrays of samples of propagated wavefields. The DFreT can be directly implemented using the so called direct method, in which case a single two-dimensional discrete Fourier transform (2D DFT) is required. In the direct method, the sampling rate (inverse of the sampling period) of the recovered wavefield depends on the aperture of the recording camera, the recovery distance d, and the recovery wavelength λ [1–3]. The DFreT can also be implemented with the spectral method where a convolution approach requiring the computation of at least two 2D DFTs is used [1–3]. In the spectral method, the sampling rate remains equal to that of the digital camera independently of d and λ [1–3].
The modification of the sampling rate in the propagated wavefield is of interest in several applications, as highlighted in [3–9]. This modification can be interpreted as a digital zoom operation. In [3], the existing zoom algorithms were classified according to the number of the required chirp multiplications and 2D DFTs. The modified spectral method [3], the double direct method [4], the chirp z transform method [5], the Rhodes light tube method [6] and the convolution approach with magnification method [7], all require three chirp multiplications and at least two 2D DFTs. In [8], a method that requires the computation of more than two 2D DFTs was proposed, which consists in computing first the angular spectrum method and then a form of direct method to obtain samples of the propagated wavefield over a tilted plane. In [9], a method based on zero padding the input array (interference pattern) and then using the direct method to compute the required sampled wavefield was proposed. Such a method can only be used to increase the sampling rate with regards to that obtained with the direct method. An advantage of the latter method is that it requires computing only a single 2D DFT, although when a very large sampling rate should be obtained this 2D DFT could be more time-consuming than the 2D DFTs required by the methods reported in [3–8].
The appropriate sampling rate depends on the frequency content of the propagated wavefield to be sampled. In DH, this frequency content is due to the recorded scene and the impulse response of the recording system [10]. The numerically reconstructed wavefield obtained by the Fresnel transform directly applied to the interference pattern (without any sort of filtering or processing) is composed of a DC term (the reference wave multiplied by the addition of the reference wave intensity and the object wave intensity), a real image and a virtual image of the recorded scene [1,2]. If an off-axis set-up is used during the recording process, these components are recovered spatially separated [2,3,11]. In such a case, and if a conjugate reference wave is used for reconstruction, the undistorted real image would contain frequency components mainly given by the original recorded scene.
When the scene to be recorded is static, the recording set-up can be adjusted or modified to maximize the frequency content in the captured interference pattern. Nevertheless, there are applications that need to record changing scenes and thus the recording set-up cannot be modified, for instance in digital holographic particle tracking velocimetry [12]. In that sort of applications, scenes with smooth objects could be recorded and the region of interest in the propagated wavefield would contain mainly low frequency components. In such case, the sampling rate obtained by the direct or the spectral method could result excessive. Any of the algorithmic methods in [3–8] can be used to reduce the recovery sampling rate, but with the cost of requiring the computation of at least two 2D DFTs. The method in [9] that requires a single 2D DFT cannot be used to reduce the sampling rate, but only to increase it.
A common practice in DH to reduce the sampling rate is to use the direct method on a truncated input array (a segment of the interference pattern recorded by the digital camera). Nevertheless, this method is equivalent to use a smaller aperture than that of the original digital camera and thus the resolution of the recovery system and the reconstructed wavefield is modified [10]. Further, if smooth objects are being recorded, truncation could eliminate relevant information required to recover a complete image of the objects. The latter happens since the light scattered from smooth objects spread out slowly (near field propagation), and thus impinges only on small regions of the digital camera during the recording process.
In this paper, we propose a method that uses an aliasing operator and a single 2D DFT to down-sample the propagated wavefield under the sampling rate obtained by the direct method of the DFreT. First, a formulation of the DFreT is described. Then, the use of an aliasing operator is proposed to introduce down-sampling control without modifying the resolution of the recording system. Finally, the usefulness of our proposal is demonstrated for digital holography when off-axis set-ups are used in the recording process.
2. The discrete Fresnel transform
In the next description of the DFreT we consider that it is directly applied over the product of the interference pattern recorded by the digital camera and the conjugate complex reference wave. Here, this product is referred to as hr(m,n). The size and sampling rate of the interference pattern are given by the characteristics of the digital camera. We consider that the digital camera (CCD or CMOS) is an array of Nx × Ny pixels with sampling periods of Δx in the x-axis and Δy in the y-axis, thus the physical size of the interference pattern is (NxΔx) × (NyΔy) = Lx × Ly. Given hr(m,n) for 0≤m≤Nx-1 and 0≤n≤Ny-1, the DFreT can be used to compute samples of the wavefield in the ξ-η plane propagated at a distance d from the interference pattern in the x-y plane. The computed array of samples is a window of Nξ × Nη elements with sampling periods of Δξ in the ξ-axis, and Δη in the η-axis, thus the physical size of the sampled wavefield is (NξΔξ) × (NηΔη) = Lξ × Lη. Here, the computed array of samples is referred to as Γ(k,l), and defined for 0≤k≤Nξ-1 and 0≤l≤Nη-1.
The DFreT can be implemented with the direct method as
where DFT2D is a 2D discrete Fourier transform, andThe terms Nx/2, Ny/2, Nξ/2 and Nη/2 in Eq. (2) and Eq. (3), and the complex exponential exp[iπ(m-Nx/2 + n-Ny/2)] in Eq. (2) are introduced to align the symmetry center of the complex exponentials to the center of processed arrays.
For the direct method of the DFreT, it is assumed that Nξ = Nx and Nη = Ny. In this case, the sampling periods Δξ and Δη, and in consequence the sampling rates (1/Δξ and 1/Δη), are given by
The sampling rates produced by Eq. (4) achieve an optimal sampling considering the maximum physical resolution that the propagated wavefield can reach [1,2]. Nevertheless, when high frequency components are not required to be recovered but only low frequency components, it could be preferred to obtain sampled wavefields with smaller sampling rates than those obtained by the direct method.
3. Down-sampling the Fresnel transform
By analyzing Eq. (4), it can be inferred that the sampling rates can be modified by changing the size of the array to be transformed. Expressions for the size of the new array are obtained by working out Nξ and Nη from Eq. (4). Thus, the number of elements in the new array should be given by
In Eq. (5), the desired sampling rates 1/Δξ and 1/Δη should be defined. It is important to highlight that since Nξ and Nη should be integers, the sampling rates can only vary in discrete steps defined by 1/Δξ = aΔx/λd and 1/Δη = bΔy/λd for positive integer values of a and b. A direct alternative of down-sampling (Nξ<Nx and/or Nη<Ny) is truncating the input array (interference pattern); in this case, the wavefield is consequently modified. Truncation could be a good alternative when rough objects are recorded since their information is spread out all over the recording medium (far field propagation), but a bad one for smooth objects since their information is spread out only in a portion of the recording medium (near field propagation). In fact, information is lost in both cases when truncating. In the former case, truncation eliminates high frequency information, while in the latter case, truncation additionally eliminates spatial information. The latter effect will be shown later.
Here, we propose a novel method based on the relationship between the discrete-time Fourier transform (FT) and the DFT to obtain an aliasing operator that can be used to down-sample the Fresnel transform without lost of information.
3.1. Aliasing operator
For the sake of brevity, the following analysis is carried out for the one-dimensional case. According to the analysis presented in [13], samples of the FT of an aperiodic sequence f(m) can be thought of as DFT coefficients of a periodic sequence g(m) obtained through summing periodic replicas of f(m). These samples are equally spaced in frequency in a range defined by 0≤ω<2π. This relationship is expressed as:
whereand Ns defines the number of required samples of the FT and, in consequence, the length of the DFT (also of g(m)). The different replicas of f(m) in Eq. (7) are given by the shift index r. Thus, aliasing (superposition of shifted replicas of a signal) can be introduced in the input domain to modify the sampling in the Fourier domain. Ns does not need to be equal to the length of the aperiodic signal f(m), here defined as Ni. If Ns<Ni, there are some shifted replicas of f(m) that contribute to the summation in Eq. (7). An example of how g(m) is generated from different replicas of f(m) is shown in Fig. 1 .
Fig. 1 Generation of g(m) from f(m): (a) example function f(m); (b) g(m) generated from the addition of shifted replicas of f(m).
Figure 1 shows that since f(m) is finite with length Ni, just a few of its replicas are required to obtain g(m). In fact, due to the different contributions of the shifted replicas, all the elements in f(m) are used to generate g(m). Thus, considering the finite length of f(m), g(m) can be computed as
where ⌊ ⌋ is the floor operator. The maximum value of the upper bound of the summation in Eq. (8) is obtained when the fraction m/Ns does not surpass the fractional part of (Ni-1)/Ns. Hence, the maximum value of the upper bound occurs whenConversely, when m>[(Ni-1) mod Ns], the upper bound takes its minimum value. Thus, the computation of g(m) can be carried out by
Equation (10) can be viewed as a general aliasing operator since Ni and Ns can take any integer value that fulfill Ns<Ni. Thus, our proposal consists in using the aliasing operator in Eq. (10) prior to computing the DFT, which will result in equally-spaced samples of the FT.
In Fig. 2 , an example of sampling the FT by using the proposed down-sampling method and truncation are shown; in this the 1D DFT of the zero padded input sequence and the 1D DFT directly applied on the input sequence are used as references. It can be observed in Fig. 2 that the values obtained with the proposed down-sampling method correspond to precise samples of the continuous frequency spectrum of the original input sequence. On the other hand, when the 1D DFT is computed using a truncated input sequence, the obtained values do not correspond to samples of the continuous frequency spectrum of the original input sequence. The aforementioned, since the elimination or modification of elements from the input sequence changes the output frequency spectrum.
Fig. 2 Samples of the FT of f(m) = [1,2,3,4,4,4,3,2,1] (0≤m≤8) obtained by using the proposed down-sampling method with Ns = 7 (Ns<Ni), the 1D DFT of a truncated input sequence with Ns = 7 (Ns<Ni Trunc), the 1D DFT of the original input sequence (Ns = Ni), and the 1D DFT of a zero padded input sequence (Ns>>Ni): (a) magnitude of the frequency spectrum; (b) phase of the frequency spectrum.
The relationship in Eq. (6) is not fulfilled when the symmetry shift Ni/2 is included in the complex exponential of the FT. For this case, the FT can be computed as
Thus, the relationship between the FT and the DFT is given by
Comparing Eq. (6) and Eq. (12), it can be noted that an extra complex exponential is required when the FT is used with a symmetry shift in the input domain. The new input sequence g(m) should be computed as specified by Eq. (10) from the original input sequence f(m).
The aliasing operator in Eq. (10) can be extended to the 2D case. First, Eq. (10) should be applied to each column of the input array generating an intermediate array, and then to each row of the previously obtained array. The elements of the 2D DFT computed with the new input array would correspond to the required samples of the Fourier transform.
3.2. Using the aliasing operator to down-sample the Fresnel transform
The values of Nξ and Nη obtained with Eq. (5) correspond to the number of equally-spaced samples of the Fresnel transform in each column and each row, respectively, that have to be computed to achieve the required sampling rate. Thus, down-sampling is required when Nξ<Nx and/or Nη<Ny, where Nx is the length of the columns and Ny of the rows in the digital camera. In order to use the aliasing operator, the array on which this should be applied has to be identified. This array is f(m,n) as expressed in Eq. (2). Thus, the aliasing operator should be applied over f(m,n) to generate the new array g(m,n) to be transformed by the 2D DFT. Since f(m,n) has columns of length Nx and rows of length Ny, the aliasing operator should be applied over the columns if Nξ<Nx and over the rows if Nη<Ny to generate g(m,n). Then, the DFreT should be implemented as
where p(k,l) is defined in Eq. (3), and the term exp[iπ(kNx/Nξ + lNy/Nη)] is introduced as indicated by Eq. (12) to correctly recover the phase information since the symmetry shifts Nx/2 and Ny/2 are used in the input domain of the DFreT (Eq. (2)).It should be noted that due to the use of the aliasing operator, our proposal is a non reversible method (f(m,n) cannot be recovered from Γ(k,l)). Thus, in order to obtain recoveries with different sampling rates, different aliasing operations over f(m,n) should be carried out.
4. Experimental results
We show sampled wavefields obtained from two off-axis digital holograms to validate the usefulness of our proposal. The reference waves are collimated beams with wavelength λrec and inclination angle θ that fulfill the Nyquist criterion. The registered objects are smooth transmission-slides. A synthetic interference pattern (H1) with Nx × Ny = 512 × 512 and Δx = Δy = 6.4 μm is generated using a reference wave with λrec = 532 nm and θ = −0.02078 rad, and a spoke target as object. A second interference pattern (H2) is recorded on a CCD sensor of Nx × Ny = 768 × 1024 and Δx = Δy = 4.7 μm using a reference wave with λrec = 632 nm and θ = 0.0286 rad, and a USAF test chart as object. For the interference pattern H1, the object was placed at a distance d = 17 cm, and for H2, it was placed at a distance d = 22.7cm. The required DFTs are computed with the FFTW [14] that is one of the most efficient algorithms to compute DFTs of any size. In Fig. 3 , the interference patterns H1 and H2 are shown. It can be observed that since the recorded objects are smooth transmission-slides, the scattered light from the objects does not spread out all over the recording sensor as stated previously.
Fig. 3 Interference patterns: (a) H1 from a synthetic spoke target as object; (b) H2 from a real USAF test target as object.
In order to recover undistorted real images, conjugate reference waves are used. In Fig. 4 , sampled in-focus real-images computed with the direct method of the DFreT are shown, where Δξ = Δη = 27.6 μm in Fig. 4(a), and Δξ = 39.75 μm and Δη = 29.81 μm in Fig. 4(b). In all the following recoveries the complete wavefields are shown together with sub-images of the region of interest.
Figure 5 and Fig. 6 show sampled in-focus real images computed with the direct method of the DFreT applied to truncated interference patterns and with the proposed down-sampling method, respectively. In Fig. 5(a) and 5(b), it can be observed that since the registered objects are smooth transmission-slides, the sampled wavefields computed from truncated interference patterns are incomplete versions of the required wavefields. As shown in Fig. 6(a) and 6(b), this problem is not present when utilizing the proposed down-sampling method since an aliasing operation is used instead of truncation. Even the small phase tilt shown in Fig. 4(b) is correctly recovered with the proposed method as observed in Fig. 6(b).
Fig. 5 Sampled in-focus real-images using truncated interference patterns for Δξ = Δη = 51.13 μm: (a) from truncated H1 with Nξ = Nη = 276; (b) from truncated H2 with Nξ = Nη = 597.
Fig. 6 Sampled in-focus real-images using the proposed down-sampling method for Δξ = Δη = 51.13 μm: (a) from H1 with Nξ = Nη = 276; (b) from H2 with Nξ = Nη = 597.
The recoveries in Fig. 6(a) and Fig. 6(b) show some down-sampled wavefields, but with heuristically chosen sampling rates. The differences between the recoveries in Fig. 4 and Fig. 6 obey the used sampling rate; in fact, it can be observed that there are well defined regions with the performed down-sampling, while it is also obvious there are others regions that require a larger sampling rate (for instance to recover the center of the spoke target image and the small numbers and bars in the USAF test target image). Note that the sampling rate cannot be larger than that obtained with the direct method (Fig. 4). A step forward from this proposal, could be the deduction of a criterion to select a suitable sampling rate, maybe based on the frequency bandwidth of the registered object. Due to the exact relationship expressed in Eq. (12), it can be affirmed that the proposed method yields precise samples of the propagated wavefield and no quantification is required to measure its performance.
5. Conclusions
We have proposed a method that requires a single 2D DFT to down-sample the Fresnel transform. This method consists in applying an aliasing operation in the space domain previous to compute the 2D DFT. Due to the use of the aliasing operator, the resolution of the recorded system is conserved resulting in precise samples of the propagated wavefield, although a non reversible method results. The size of the used 2D DFT is smaller than those required by the previously reported methods. It was demonstrated that the proposed method can be used in digital holography even to retrieve images of smooth objects. Hence, this proposal is an effective, precise and efficient method to down-sample wavefields computed with the discrete Fresnel transform.
References and links
1. T. M. Kreis, M. Adams, and W. P. P. Jüptner, “Methods of digital holography: a comparison,” Proc. SPIE 3098, 224–233 (1997). [CrossRef]
2. U. Schnars and W. Jueptner, Digital Holography (Springer, 2005).
3. B. M. Hennelly, D. P. Kelly, D. S. Monaghan, and N. Pandey, “Zoom Algorithms for Digital Holography,” in Information Optics and Photonics: Algorithms, Systems, and Applications, T. Fournel and B. Javidi, eds. (Springer, 2010), pp. 187–204.
4. X. Deng, B. Bihari, J. Gan, F. Zhao, and R. T. Chen, “Fast algorithm for chirp transforms with zooming-in ability and its applications,” J. Opt. Soc. Am. A 17(4), 762–771 (2000). [CrossRef] [PubMed]
5. F. Zhang, I. Yamaguchi, and L. P. Yaroslavsky, “Algorithm for reconstruction of digital holograms with adjustable magnification,” Opt. Lett. 29(14), 1668–1670 (2004). [CrossRef] [PubMed]
6. W. T. Rhodes, “Light tubes, Wigner diagrams and optical signal propagation simulation,” in Optical Information Processing: A Tribute to Adolf Lohmann, H. J. Caulfield, ed. (SPIE Press, 2002), pp. 343–356.
7. J. C. Li, P. Tankam, Z. J. Peng, and P. Picart, “Digital holographic reconstruction of large objects using a convolution approach and adjustable magnification,” Opt. Lett. 34(5), 572–574 (2009). [CrossRef] [PubMed]
8. L. Yu and M. K. Kim, “Pixel resolution control in numerical reconstruction of digital holography,” Opt. Lett. 31(7), 897–899 (2006). [CrossRef] [PubMed]
9. P. Ferraro, S. De Nicola, G. Coppola, A. Finizio, D. Alfieri, and G. Pierattini, “Controlling image size as a function of distance and wavelength in Fresnel-transform reconstruction of digital holograms,” Opt. Lett. 29(8), 854–856 (2004). [CrossRef] [PubMed]
10. D. P. Kelly, B. M. Hennelly, N. Pandey, T. J. Naughton, and W. T. Rhodes, “Resolution limits in practical digital holographic systems,” Opt. Eng. 48(9), 095801 (2009). [CrossRef]
11. E. N. Leith and J. Upatnieks, “Reconstructed wavefronts and communication theory,” J. Opt. Soc. Am. 52(10), 1123–1130 (1962). [CrossRef]
12. S. Satake, T. Kunugi, K. Sato, and T. Ito, “Digital Holographic Particle Tracking Velocimetry for 3-D Transient Flow around an Obstacle in a Narrow Channel,” Opt. Rev. 11, 162–164 (2004).
13. A. V. Oppenheim, R. W. Schafer, and J. R. Buck, Discrete-time signal processing (Prentice Hall, 1999).
14. M. Frigo and S. G. Johnson, “The FFTW web page,” (2007), http://www.fftw.org/.
