1. Introduction

Optical holography can produce very realistic virtual images due to its capability to properly convey the depth cues that we use to interpret three-dimensional objects. Computational holography is the use of digital representations plus computational methods to carry out the holographic operations of construction and reconstruction; modeling the propagation of complex-valued wavefronts through three-dimensional space [1]. The large computational requirements of holographic simulations prohibit present-day existence of real-time holographic displays comparable in size to traditional two-dimensional displays.

Various approaches to computational hologram construction have been researched, of which the most prominant are interference-based [2, 3, 4], Fourier-based [5, 6, 7], and diffraction-specific [8, 9] hologram computation. This research addresses the Fourier-based approach. Fourier-based hologram computation is built around the fast Fourier transform (FFT), which provides efficient computational algorithms. The disadvantages of traditional Fourier techniques are that the sampling rates are fixed and that the sampled points of the hologram and object must both be centered along the axis of propagation [10, 11]. This research develops a fast computational method to construct and reconstruct holograms with arbitrary sampling rates and with a spatial shift away from the axis of propagation.

Section 2 gives an overview of the traditional discretization used in computing Fresnel diffraction. Section 3 introduces and derives a new form of the Fresnel diffraction integral that solves a problem with traditional Fourier-based hologram computation by allowing for an arbitrary spatial shift of the wavefront. We introduce and discuss the concept of tiling in Section 4, which allows for propagation of discrete wavefronts between sampled planes having different resolutions. In Section 5 we give a demonstration of our techniques, and in Section 6 we give an overview of our results and possible future enhancements.

2. Traditional discrete Fresnel diffraction

Traditional approaches that discretize the Fresnel diffraction integral to obtain a fast algorithm restrict the sampled points to be centered about the axis of propagation and constrain the sampling rates involved. The traditional discretization of the Fresnel diffraction integral begins with the following form of continuous Fresnel diffraction [10]:

where Ψ_d(x) represents the Fresnel diffraction pattern produced at a distance d by a one-dimensional object that is illuminated by a monochromatic plane wave of wavelength λ and amplitude ψ(x′). Extension of this relation to produce the Fresnel pattern of two-dimensional planar objects is straightforward, as is illustrated in Section 3. Let N be the sampled resolution in each domain, and let the sampling pitches be Δx and Δx′ in the x and x′ directions, respectively. The sampled Fresnel transform becomes:

where the integer indices m and p both range from -N/2-1 to N/2. If we impose the condition

then Eq. (2) can be written as

which has the structure of a one-dimensional, discrete Fourier transform,

The discrete function U_d(m) can be calculated by first taking the product of the input sampled function ψ(pΔx′) and the exponential function that is inside the braces under the summation, followed by an evaluation of the discrete Fourier transform of the one-dimensional array [effectively, g(p)] that results from this product, and finally multiplication by the leading terms outside the summation in Eq. (4).

The rather stringent condition imposed by Eq. (3) is traditionally made so that the sampled Fresnel transform can be computed using an FFT. However, a versatile method for generating discrete holographic images must allow for an independent selection of the various parameters that appear in Eq. (3), for example, an independent specification of the sampling pitches Δx and Δx′ and the sampled resolution N. Hence, it will be necessary to violate the condition imposed by Eq. (3) and simultaneously abandon the traditional Fourier-based approach to determining U_d(m).

3. Fresnel diffraction between shifted parallel planes

The techniques developed in this section overcome the limitations of traditional Fourier-based holography by allowing for arbitrary spatial separation from the axis of propagation and arbitrary pitch values, while also permitting an efficient computational algorithm.

We begin by defining a finite rectangular cartesian grid in the x-y and x′-y′ planes. Fig. 1 illustrates the geometric setup described. The x′-y′ plane is simply a shift of the x-y plane by a distance of d in the z direction. The x-y grid has one corner located at (x ₀ ,y ₀), with a resolution of M × N, and has x and y pitch values of Δx and Δy, respectively. The x′-y′ grid is defined similarly, having one corner located at (x′ ₀,y′ ₀), with a resolution of P × Q, and has x′ and y′ pitch values of Δx′ and Δy′, respectively. The coordinates of the set of indexed points for these two grids are determined by the following formulas:

where m, n, p, and q are integer indices with the following bounds:

 

Fig. 1. Shifted Fresnel diffraction geometry.

Download Full Size | PDF

Using the discrete functions U_d(m,n) and u(p,q) to represent the complex amplitude of the target and source wavefronts at these discrete coordinate locations, that is, defining

the two-dimensional analog of Eq. (2) yields the desired discrete shifted Fresnel transform, explicitly:

We note that if the number of image elements in the target plane matches the number in the source plane (i.e., if M = P and N = Q) this expression has the structure of a two-dimensional, scaled discrete Fourier transform, explicitly,

where, in our particular case, the scale parameters are s = (ΔxΔx′/λd) and t = (ΔyΔy′/λd). Note also that choosing the scale parameters s = 1/P= 1/M and t = 1/Q = 1/N results in the traditional Fourier transform, but we do not want to impose these stringent constraints here.

The technique that has been used to evaluate our scaled Fourier transform is based on the one developed by Bailey and Swarztrauber [12]. (We will not use their terminology of a “fractional” Fourier transform to avoid confusion with a different but more widely used definition of a fractional Fourier transform.) They noted that since (k- j)² = j ² + k ² -2jk, a one-dimensional scaled discrete Fourier transform can be rewritten as,

where,

Bailey and Swarztrauber [12] then pointed out that the final summation in Eq. (15) is in the form of a convolution. A discrete convolution can be computed by first restructuring the data into a circular convolution, which doubles the required storage space in each dimension, and then by applying the Convolution Theorem, which results in a computational cost proportional to three FFTs. This process is discussed in more detail in [12]. Given that the complexity of a two-dimensional FFT is O(N ² log₂ N), then the complexity of a two-dimensional scaled FFT by the method discussed here is

Note that the traditional Fresnel discretization stated in Eq. (2) can also make use of the fast scaled Fourier transform in order to avoid the restriction stated in Eq. (3). The derivation given above for shifted Fresnel diffraction was pursued to allow for arbitrary pitch values as well as a shift away from the axis of propagation.

The major restriction imposed by the fast shifted Fresnel transform, as presented, is that the number of discrete elements in the target image U_d(m,n) must be equal to the number of discrete elements in the source image u(p,q), that is, M = P and N = Q. This restriction will be overcome in the following section where we develop a tiling algorithm to handle source and target wavefronts with different size image arrays.

4. Rectangular tiling

The shifted Fresnel formula developed in Section 3 will be utilized in a technique we refer to as rectangular holographic tiling, or simply, rectangular tiling. This development is motivated by the fact that the imaging volume requires much less spatial resolution than the holographic fringe pattern. Human acuity is limited to approximately 0.175mm [8], meaning that two points of light closer than 0.175mm appear as one. Contrast this with holographic fringe patterns, which typically require spatial resolutions on the order of a micron. Thus, given a hologram and an imaging volume of roughly the same lateral size, the hologram will require a much higher spatial resolution. Thus fast numerical techniques allowing for the calculation of a diffraction pattern between vastly different spatial resolutions is needed.

Here we define rectangular tiling as the subdivision of a finite rectangular planar surface into smaller, possibly overlapping, rectangular planar surfaces referred to as tiles. These tiles form a set, and together their union covers the original planar surface. This restriction to rectangular tiling enables the use of the fast algorithm developed in Section 3. Our tiling construct allows for different resolutions for the imaging volume and the hologram, in either the construction or reconstruction of holograms. We begin with a source plane on which the values of a complex wavefront are known, and a target plane where we will calculate the values of the propagated wavefront. For hologram construction, the source plane is in the imaging volume and the target plane is the hologram; and vice versa for hologram reconstruction. As before, let u(p, q) be the discrete representation of the complex-valued wavefront ψ(x′,y′) at the source plane, and let U_d(m,n) be the discrete representation of the complex-valued wavefront ψ(x′, y′) at the target plane. If u and U_d have the same array dimensions in both x and y — i.e., if M = P and N= Q— then the computation makes direct use of the fast discrete shifted Fresnel transform (Eq. (12)). Otherwise, we will consider two cases: when M > P and N > Q; or when M < P and N < Q.

In the first case, when the source wavefront ψ(x′,y′) is sampled at a lower resolution than the target wavefront Ψ(x,y), we create a tiling of U_d such that each tile has the same array dimensions as that of u, as illustrated in Fig. 2. Then we sequentially compute the propagated wavefront of ψ for each tile. Since each tile has the same resolution as ψ, the fast shifted Fresnel transform may be used. In the second case, when the target wavefront Ψ(x,y) is sampled at a lower resolution than the source wavefront ψ(x′,y′), we use the superposition principle of light. We create a tiling of u such that each tile has the same array dimensions as that of U_d, as illustrated in Fig. 3. The wavefront Ψ is then determined from the linear superposition of the propagated wavefront of each tile of ψ.

 

Fig. 2. Lower resolution source wavefront for holographic tiling.

Download Full Size | PDF

Let us look at the complexity of the tiling algorithm in the first case just discussed. For simplicity, assume ψ has a resolution of Q × Q and Ψ has a resolution of N × N, such that 2^KQ = N for some non-negative integer K. The number of tiles generated to cover Ψ is (2^K)² and the complexity to transform each tile is O(Q ² log₂ Q). Therefore the total complexity is:

 

Fig. 3. Lower resolution target wavefront for holographic tiling.

Download Full Size | PDF

Without tiling, we would only be able to utilize the fast scaled Fourier transform if the array with the smaller dimensions, u, were enlarged to the dimensions of U_d and padded with zeros, in which case, as illustrated, the computational costs would have been O(N ² log₂ N). Hence we see that tiling reduces the computational costs.

The complexity of the tiling algorithm in the second case is derived similarly, except that the cost of the linear superpositioning is considered. Now assume that ψ has a resolution of N × N and Ψ has a resolution of Q × Q, such that 2^KQ = N for some non-negative integer K. The number of tiles generated to cover ψ is (2^K)², the complexity to transform each tile is O(Q ² log₂ Q), and the complexity of adding each tile’s contribution to the final wavefront at Ψ is O(Q ²). Therefore the total complexity is:

In considering holograms computed from multiple sources, as demonstrated in Section 5, the complexity increases linearly with the number of source wavefronts owing to the linear super-positioning property of holography.

5. Demonstration

In traditional Fourier-based hologram calculations, the digital arrays representing the discrete source and target wavefronts must have the same dimensions (i.e., the two wavefronts must be sampled at the same resolution) and they must both be centered along the axis of propagation. The demonstration that we present here essentially simulates on-axis holography, but with the computation of a hologram from three source images located at different depths, two of which are offset from the axis of propagation. The three two-dimensional source images that we have chosen to use (see Fig. 4) each have a sampling resolution of 256 × 256 and a size of 4 mm × 4 mm; hence, as itemized in Table 1, they each have a pitch of Δx′ = Δy′ = 15.6 μm. In terms of the Cartesian coordinates (x′ ₀, y′ ₀, d) defined in Fig. 1, Image A was located at (-2 mm, 2 mm, 50 cm), Image B was located at (0 mm, 0 mm, 52 cm), and Image C was located at (2 mm, -2 mm, 54 cm). A single hologram was calculated using a sampling resolution of 1024 × 1024 with a pitch Δx = Δy = 8 μm. This hologram, shown in Fig. 5, was computed as the superposition of the holographic fringe pattern that was produced by each of the three source images. The computation from each source image applied the tiling algorithm previously discussed, which in turn utilized the shifted Fresnel diffraction formula.

Table 1. Properties of images for tiling demonstration (λ = 633 nm; Δx = Δy = 8 μm).

View Table

 

Fig. 4. Source images used for demonstration having various depths.

Download Full Size | PDF

Three simulated reconstructions were carried out, each at a depth d of one of the original source images. The three reconstructed images are shown in Figures 6, 7, and 8. In each reconstruction, the proper image is in focus, located in the proper location and orientation, and the other two images are properly seen as out of focus. These reconstructions were also computed using the tiling algorithm and the fast shifted Fresnel transform.

A second demonstration was carried out to show the reduction of computational time of the tiling algorithm when compared to traditional techniques. Each source image was embedded in a 1024×1024 image in order to match the sampled resolution of the hologram, thus allowing computation of the hologram without use of the tiling algorithm. Both demonstrations were implemented using Java, executed on an Intel Dothan PC under Linux, representing a prototype implementation. The average computational time was 68 seconds for the first demonstration and 102 seconds for the second demonstration, representing an approximate decrease in computational time of 33 percent for the tiling algorithm. The hologram and reconstructions from the second demonstration, while not shown, were visually identical to those of the first demonstration.

Optical reconstruction of the demonstration holograms presented here was not available to us at the time of publication. However, reconstruction of earlier demonstration holograms which utilized both the shifted Fresnel diffraction and tiling algorithms has been carried out. These earlier holograms were optically reconstructed with help from a group at the University of Texas Southwestern Medical Center at Dallas[4]. The holographic display device used by this group consisted of a 17 μm pitch Texas Instruments digital micromirror device and a 633 nm laser. We are therefore confident that if the computed hologram shown here in Fig. 5 were projected through such a device, its optical reconstruction would match well with the computed reconstructions shown here in Figures 6–8.

 

Fig. 5. Hologram computed.

Download Full Size | PDF

 

Fig. 6. Reconstruction of image A.

Download Full Size | PDF

 

Fig. 7. Reconstruction of image B.

Download Full Size | PDF

 

Fig. 8. Reconstruction of image C.

Download Full Size | PDF

6. Conclusion

In this investigation, a discrete shifted Fresnel transform has been derived (Section 3) to overcome disadvantages of traditional Fourier-based computational holography techniques, namely that the sampling rates are fixed and the sampled points of the hologram and object must both be centered along the axis of propagation. This shifted Fresnel transform permits computation of Fresnel diffraction between two parallel planar surfaces. An efficient implementation, the fast shifted Fresnel transform, has been developed that utilizes a fast scaled Fourier transform for computation. To compute Fresnel diffraction between planar surfaces with different sampling resolutions, a tiling algorithm has been developed that uses the shifted Fresnel transform (Section 4).

The fast discrete shifted Fresnel transform developed in this research has complexity O(N log₂ N) for one-dimensional applications and complexity O(N ² log₂ N) for two-dimensional applications. This is more efficient than a “brute-force” approach that results in complexity, respectively, of O(N ²) or O(N ⁴). The fast calculation of our shifted Fresnel transform has relied on the computation of a scaled Fourier transform, which in turn has relied on applying the Convolution Theorem and computing three Fourier transforms. The computational time could be reduced even further if an alternative method of computing a scaled Fourier transform were found that required fewer than three executions of the FFT.

In the development of the fast shifted Fresnel transform, we have assumed that the source and target wavefronts were defined on parallel planes. But we have demonstrated that three-dimensional structure can be built through the superposition of holographic images that are computed from parallel planes at different depths. It would be beneficial to extend the formulations in this research to arbitrarily tilted finite rectangular planes as well as other geometric primitives, for example, cylinders or spheres. There has been much research in diffraction from tilted planes, and the interested reader may see [13, 14, 15, 16, 17].