Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Fast calculation of computer-generated holograms based on 3-D Fourier spectrum for omnidirectional diffraction from a 3-D voxel-based object

Open Access Open Access

Abstract

We have derived the basic spectral relation between a 3-D object and its 2-D diffracted wavefront by interpreting the diffraction calculation in the 3-D Fourier domain. Information on the 3-D object, which is inherent in the diffracted wavefront, becomes clear by using this relation. After the derivation, a method for obtaining the Fourier spectrum that is required to synthesize a hologram with a realistic sampling number for visible light is described. Finally, to verify the validity and the practicality of the above-mentioned spectral relation, fast calculation of a series of wavefronts radially diffracted from a 3-D voxel-based object is demonstrated.

© 2012 Optical Society of America

1. Introduction

Since its invention by Gabor in 1948, holography [1] has been evolving and its applications have been proposed [2,3]. Among the applications, a 3-D display involving holography has been the object of many researchers’ attention as the most ideal display that meets all physiological and perceptual requirements. Computer-generated hologram (CGH), in which virtual objects in a computer can be reconstructed in real space, has recently been evolving rapidly along with the development of computers. This technique has the potential to be used for the realization of 3-D displays and virtual reality [4].

However, the enormous amount of diffraction calculation required in CGH is an obstacle in using CGH for realizing a 3-D display. Some effective methods such as the use of a GPU [5], the use of a look-up table [6], and the approximations of root operations [7] have been proposed to reduce the calculation time. The calculation method involving fast Fourier transforms (FFTs), in which the formula is expressed in the form of a Fourier transform, is the most commonly used method, and it drastically reduces the calculation time. However, while this method is applicable only under certain conditions, where an object plane and an observation plane should be parallel and their optical axes should be coincident, there are many cases where diffraction calculation in a direction different from the optical axis direction is required. Leseberg et al. [8] and Tommasi et al. [9] proposed methods involving the use of FFTs for diffraction calculation even in the case where an object plane and an observation plane are not parallel. A rigorous approach [10] involving these methods and the concept of angular spectrum [11] has also been reported.

However, the application of these methods involving FFTs is limited to 2-D planar objects, and the methods cannot be applied to a 3-D voxel-based object directly. Moreover, the spectral treatment in past studies has been confined to 2-D scenarios although diffraction between non-parallel planes involves light propagation in 3-D space. Thus, such spectral treatment requires complicated understanding of the calculation algorithm depending on geometrical alignment and its physical meaning, and limits the application potential.

Here, we have derived the basic spectral relation between a 3-D object and its diffracted wavefront by interpreting the diffraction calculation in the 3-D Fourier domain. This relation clearly indicates the information that is inherent in the diffracted wavefront on the 3-D object. In addition, this relation facilitates fast calculation of the diffraction of a 3-D object in an arbitrary direction; FFTs are used in the calculation. Unlike other methods used for the calculation of diffraction between non-parallel planes [810], in our method the spectral treatment is conducted not in the 2-D Fourier domain but in the 3-D Fourier domain. 3-D spectral approach makes it possible to graphically and intuitively understand coordinate transformations and resampling, which are inevitable for the fast calculation of diffraction between non-parallel planes; it is difficult to gain such an understanding through 2-D spectral approaches

In the following section, the basic spectral relation between a 3-D object and its diffracted wavefront is derived. It is found that the 2-D spectrum of the diffracted wavefront is completely identical to the spectrum on the surface of a hemisphere whose center is at the origin and whose radius is the reciprocal of the wavelength in the 3-D Fourier domain of the 3-D object. To apply the spectral relation to CGH designed for visible light, a method for obtaining the spectrum on a spherical surface with a feasible sampling number is also described. We numerically demonstrate fast calculation of a series of radially diffracted wavefronts.

2. Basic spectral relation between a 3-D object and its diffracted wavefront in 3-D Fourier domain

2.1. Derivation of the basic spectral relation

To elucidate the spectral relation, a virtual optical system and a coordinate system are introduced in Fig. 1. Without loss of generality, a 3-D object o(x0, y0, z0) can be assumed as the aggregation of point light sources isotropically emitting spherical wavefronts. Here, a spherical wave g(x0, y0, z0) is emitted from an individual point light source, and the superposition of spherical waves from all the point light sources is observed at point P(x,y,z) in 3-D space. Therefore, the wavefront f (x,y,z) diffracted from the object o(x0, y0, z0) can be expressed as follows:

f(x,y,z)=o(x0,y0,z0)g(xx0,yy0,zz0)dx0dy0dz0
=o(x,y,z)*g(x0,y0,z0),
where
g(x,y,z)=exp(i2πλx2+y2+z2)x2+y2+z2.

 figure: Fig. 1

Fig. 1 Schematic of the virtual optical system.

Download Full Size | PDF

The symbols * and λ represent 3-D convolution and the wavelength of the light sources, respectively. g(x,y,z), expressed by Eq. (3), represents a spherical wave. Because Eq. (3) is a solution of the scalar Helmholtz equation in 3-D free space, Eq. (1) is regarded as a rigorous diffraction integral. The 3-D Fourier spectrum G(u,v,w) of g(x,y,z) is given by [12]

G(u,v,w)=F[g(x,y,z)]=14π2(u2+v2+w21/λ2),
where (u,v,w) are the coordinates of (x,y,z) in the Fourier domain and F[·] denotes a 3-D Fourier transform operator. By converting Eq. (2) into the form of a Fourier transform, the following equation is obtained:
f(x,y,z)=14π2O(u,v,w)u2+v2+w21/λ2exp{i2π(ux+uy+wz)}dudvdw,
where O(u,v,w) is the 3-D Fourier spectrum of o(x,y,z). By considering the diffracted wavefront at the plane U(x,y,z = R), Eq. (5) can be written as
f(x,y)|z=R=14π2O(u,v,w)u2+v2+w21/λ2exp{i2π(ux+uy+wR)}dudvdw.
Next, to analytically calculate the integral in Eq. (6) with respect to the w axis, the integral is performed by considering the complex plane w = ξ + . The path of the line integral in the w plane is set to be a closed contour in which only the singularity w+=1/λ2u2v2 is included (blue solid line in Fig. 2). The other singularity w=1/λ2u2v2 is excluded since we take into account only the wavefront propagating in the +z0 direction. The radius of the hemicycle s is extended to infinity. By using the residue theorem, this integral is calculated by considering the contribution only from the singularity inside the closed contour:
f(x,y)|z=R=i4πO(u,v1/λ2u2v2)1/λ2u2v2exp(i2πR1/λ2u2v2)×exp{i2π(ux+uy)}dudv.

 figure: Fig. 2

Fig. 2 The path of the line integral (blue solid line).

Download Full Size | PDF

Finally, both sides of this equation are 2-D Fourier transformed.

F(u,v)|z=R=i4πO(u,v1/λ2u2v2)1/λ2u2v2exp(i2πR1/λ2u2v2)
This equation indicates the following important fact. When w(u,v)=1/λ2u2v2 is introduced, the equation u2 + v2 + w(u,v)2 = 1/λ2 is obtained. Therefore, Eq. (8) implies that the 2-D Fourier spectrum of the diffracted wavefront can be obtained by extracting the spectrum on the surface of the hemisphere whose center is located at the origin and whose radius is 1/λ, from the 3-D spectrum O(u,v,w) and by multiplying the weight factor and the phase components corresponding to the diffraction distance R. In other words, the spectrum on the hemispherical surface in the 3-D spectrum of the 3-D object inheres in the diffracted wavefront. This is a simple but important conclusion, and it clarifies the information contained in the diffracted wavefront, which is equivalent of the holographic information. Therefore, Eq. (8) can be regarded as a basic spectral relation for diffraction.

This conclusion and the procedure on which it is based are almost identical to those on diffraction tomography [12]. However, while in our study, the resulting Fourier spectrum is on a hemisphere whose center is located at the origin, in diffraction tomography, the resulting Fourier spectrum is on a hemisphere passing through the origin. This difference results from the fact that while in our study, it is not necessary to take into account the phase term of the incident light because a 3-D object is defined as self-luminous point light sources, in diffraction tomography, it is necessary to consider the phase term of the incident light because a series of diffracted wavefronts with varying incident angles of the incident light is considered. As will become apparent later, this difference has a special significance for our fast calculation method.

2.2. Illustration of object information inherent in diffracted wavefronts

A diagram is shown in Fig. 3 to understand Eq. (8) graphically and intuitively. As shown in Fig. 3, the spectrum of a diffracted wavefront in the +z0 direction is equivalent to the spectrum on a hemisphere with w as the axis of rotational symmetry. Similarly, the spectrum in the +x0 direction corresponds to a hemisphere with u as the axis of rotational symmetry, and its spectrum is partially shared with that of diffraction in the +z0 direction. This correspondence can be considered to hold for omnidirectional diffraction, including diffraction in the −z0 and −x0 directions. Therefore, it is possible to express diffracted wavefronts in all directions by a single spectrum on the spherical surface. In other words, it is concluded that the single spectrum on the spherical surface includes all information on diffraction. This simple result, which is easy to understand graphically and intuitively, is drastically different from that obtained from diffraction tomography and is the basis of the fast calculation of wavefronts using FFTs, as explained later.

 figure: Fig. 3

Fig. 3 Spectrum on the hemispherical surface corresponding to diffracted wavefronts in different directions.

Download Full Size | PDF

Needless to say, this result can be applied to a 2-D planar object by defining a 2-D object as a 3-D object distributed in a cross-sectional plane. This 3-D treatment of the 2-D object makes it easy to understand the calculation algorithm of the past studies on the diffraction between non-parallel planes [810] and its physical meaning. In such diffraction calculation, although the diffraction integral can be expressed in the form of Fourier transforms by the coordinate transformations corresponding to the inclination between non-parallel planes, the diffracted wavefront or its Fourier spectrum should be resampled at even intervals in the distorted coordinate system for the fast calculation using FFTs. However, these procedures by the 2-D spectral approach are somewhat difficult to understand and implement with computers. On the other hand, when such diffraction integral is considered with the 3-D spectral approach proposed here, these procedures are regarded as the simple extraction of the spectrum on the hemispherical surface, whose axis of rotational symmetry is corresponding to the inclination, and then the extracted spectrum is resampled in the orthogonal coordinate system perpendicular to the axis of rotational symmetry. Therefore, the 3-D spectral approach makes it possible to graphically and intuitively understand the diffraction calculation between non-parallel planes. Obviously, a 2-D FFT is enough to obtain the 3-D Fourier spectrum of such a 2-D planar object.

3. Acquisition of the spectrum on the spherical surface for visible light

The spectrum on a hemispherical surface with radius 1/λ is needed to calculate a wavefront diffracted from the 3-D object according to Eq. (8). It is however impractical to directly calculate a spectrum with a bandwidth of more than 2/λ in each direction by 3-D FFT for visible light. To solve this problem, we consider a 3-D object not as a continuously distributed object o(x,y,z) but as a 3-D object o′(x,y,z) sampled at intervals of δ = 1/W. Its 3-D Fourier spectrum O′(u,v,w) can be written as

O(u,v,w)=F[o(x,y,z)comb(xδ,yδ,zδ)]=1WO(u,v,w)*comb(uW,vW,wW),
where comb(x,y,z) is a comb function. When the spectral bandwidth of o(x,y,z) is limited to less than W in each direction, Eq. (9) implies that the original spectrum O(u,v,w) of the 3-D object o(x,y,z) is just repeated with period W in each direction. The full range of the spectrum is obtainable by the simple repetition of the spectrum O(u,v,w) with a bandwidth W. The sampled values of O(u,v,w) are obtained by the 3-D FFT of o(δnx, δny, δnz), where nx, ny and nz are integers. This process facilitates the application of the above-mentioned spectral relation with a feasible sampling number even to visible light.

4. Application of the method of fast calculation of a series of radially diffracted wavefronts

It has been found that obtaining the spectrum on the spherical surface in the 3-D Fourier domain is necessary to calculate wavefronts diffracted in any direction. After the acquisition of the spectrum on the spherical surface, the calculation procedure of the diffracted wavefronts consists of the extraction of the hemispherical spectrum corresponding to the diffraction direction, multiplication of the weight factor and phase term, projection of the spectrum onto the 2-D cross-sectional plane perpendicular to the axis of rotational symmetry, and inverse 2-D FFT. The time-consuming process among these is only the inverse 2-D FFT. Thus, fast calculation of a series of diffracted wavefronts in any direction is possible provided that the spectrum on the spherical surface has already been obtained. This method is extremely effective for the cases where the fast calculation of diffracted wavefronts in various directions is important, such as interactive CGH performed by varying the viewing angle, and polygonally arranged computer-generated holograms.

Here, as the first step in such applications, fast calculation of a series of radially diffracted wavefronts is demonstrated through a simulation by using the virtual optical system shown in Fig. 4 to show the validity of our method. As shown in Fig. 4, a globe on which the world map shown in Fig. 5 is spherically mapped is defined as a 3-D object o. The diffraction direction is identified by θ and ϕ, which correspond to the latitude and longitude of the globe, respectively. The diffraction distance R is the same as the radius of the spherical object. On the basis of this layout, the diffracted wavefronts are calculated by using FFTs by varying the diffraction direction θ and ϕ. The size and the sampling number of the 3-D computational area as the input are 256 × 256 × 256 and 2mm×2mm×2mm, respectively. The 3-D object o is defined inside this area and its radius R is 0.69mm. The wavelength is set to 632.8nm. Some simulation results obtained under these conditions are shown in Fig. 6. The central parts of Figs. 6(a) and 6(b) are in focus while the parts away from the center are out of focus. This is because the diffraction distance is set to the radius of the spherical object. The focused parts correspond to the diffraction direction; this can be seen by making a comparison with Fig. 5. On the other hand, since the diffraction direction of Fig. 6(c) is opposite to that in Fig. 6(b), Fig. 6(c) is the horizontal reversal of the defocused image of Fig. 6(b). From these results, it is found that our method can calculate diffracted wavefronts properly. Moreover, some results for a series of diffracted wavefronts obtained by continuously varying the diffraction direction are shown as images of movies in Fig. 7. Figures 7(a) ( Media 1) and (b) ( Media 2) show patterns for cases where 360 diffracted images have been calculated by varying ϕ and θ by 1°, respectively. Although the interference patterns between the wavefronts diffracted from the front surface and the rear surface are seen in the movies, only the central parts of the front surface are in focus. The calculation time required for the 360 diffracted images is only about 20s, while it takes about 32h to obtain one diffracted image when the diffraction integral of Eq. (1) is calculated directly without FFTs. This shows that our method involving the use of FFTs practically performs the calculation 2 × 106 times faster than the method involving the direct calculation of Eq. (1).

 figure: Fig. 4

Fig. 4 Schematic of the virtual optical system for the observation of radially diffracted wavefronts: (a) the perspective view and (b) top view.

Download Full Size | PDF

 figure: Fig. 5

Fig. 5 World map that is spherically mapped onto the object surface.

Download Full Size | PDF

 figure: Fig. 6

Fig. 6 Diffraction patterns in several directions: (a) θ = 46° and ϕ = 135°, (b) θ = 55° and ϕ = 0°, and (c) θ = −55° and ϕ = 180°.

Download Full Size | PDF

 figure: Fig. 7

Fig. 7 Images of movies showing 360 diffraction patterns; the patterns are obtained by varying the (a) diffraction angle θ ( Media 1) and (b) diffraction angle ϕ ( Media 2).

Download Full Size | PDF

5. Conclusion

We have derived the basic spectral relation between a 3-D object and its diffracted wavefront through rigorous calculation of the diffraction integral in the 3-D Fourier domain. It is found that the spectrum on the spherical surface includes all information on diffraction. This result is very simple and provides an intuitive understanding of diffraction calculation involving the use of FFTs. After the explanation of the procedure for obtaining the required spectrum for visible light, the application of our method to the fast calculation of a series of radially diffracted wavefronts is discussed, and the validity and high performance of the method are demonstrated.

Our method is extremely effective for interactive CGH performed by varying the viewing angle arbitrarily and polygonal arrangement of CGHs for the purpose of achieving a wide viewing angle, and it can also be used for obtaining CGHs from 3-D Fourier spectra obtained by X-ray CT, ultrasonic CT, MRI, and so on.

References and links

1. D. Gabor, “A new microscopic principle,” Nature 161, 777–778 (1948). [CrossRef]   [PubMed]  

2. L H. Lin and H. L. Beauchamp, “Write-read-erase in situ optical memory using thermoplastic holograms,” Appl. Opt. 9, 2088–2092 (1970). [CrossRef]   [PubMed]  

3. X. Zhang, E. Dalsgaard, S. Liu, H. Lai, and J. Chen, “Concealed holographic coding for security applications by using a moiré technique,” Appl. Opt. 36, 8096–8097 (1997). [CrossRef]  

4. A. W. Lohmann and D. P. Paris, “Binary Fraunhofer holograms, generated by computer,” Appl. Opt. 6, 1739–1748 (1967). [CrossRef]   [PubMed]  

5. T. Shimobaba, T. Ito, N. Masuda, Y. Ichihashi, and N. Takada, “Fast calculation of computer-generated-hologram on AMD HD5000 series GPU and OpenCL,” Opt. Express 18, 9955–9960 (2010). [CrossRef]   [PubMed]  

6. M. Lucente, “Interactive computation of hologram using a look-up table,” J. Electron. Imaging 2, 28–34 (1993). [CrossRef]  

7. T. Yamaguchi, T. Fujii, and H. Yoshikawa, “Fast calculation method for computer-generated cylindrical holograms,” Appl. Opt. 47, D63–D70 (2008). [CrossRef]   [PubMed]  

8. D. Leseberg and C. Frère, “Computer-generated holograms of 3-D objects composed of tilted planar segments,” Appl. Opt. 27, 3020–3024 (1988). [CrossRef]   [PubMed]  

9. T. Tommasi and B. Bianco, “Computer-generated holograms of tilted planes by a spatial frequency approach,” J. Opt. Soc. Am. A 10, 299–305 (1993). [CrossRef]  

10. K. Matsushima, H. Schimmel, and F. Wyrowski, “Fast calculation method for optical diffraction on tilted planes by use of the angular spectrum of plane waves,” J. Opt. Soc. Am. A 20, 1755–1762 (2003). [CrossRef]  

11. J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996).

12. A. C. Kak and M. Slaney, “Tomographic Imaging with Diffracting Sources,” in Principles of Computerized Tomographic Imaging, R. F. Cotellessa, J. K. Aggarwal, and G. Wade, eds. (Institute of Electrical and Electronics Engineers, 1988), pp. 203–273.

Supplementary Material (2)

Media 1: MPG (3730 KB)     
Media 2: MPG (3730 KB)     

Cited By

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

Alert me when this article is cited.


Figures (7)

Fig. 1
Fig. 1 Schematic of the virtual optical system.
Fig. 2
Fig. 2 The path of the line integral (blue solid line).
Fig. 3
Fig. 3 Spectrum on the hemispherical surface corresponding to diffracted wavefronts in different directions.
Fig. 4
Fig. 4 Schematic of the virtual optical system for the observation of radially diffracted wavefronts: (a) the perspective view and (b) top view.
Fig. 5
Fig. 5 World map that is spherically mapped onto the object surface.
Fig. 6
Fig. 6 Diffraction patterns in several directions: (a) θ = 46° and ϕ = 135°, (b) θ = 55° and ϕ = 0°, and (c) θ = −55° and ϕ = 180°.
Fig. 7
Fig. 7 Images of movies showing 360 diffraction patterns; the patterns are obtained by varying the (a) diffraction angle θ ( Media 1) and (b) diffraction angle ϕ ( Media 2).

Equations (9)

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

f ( x , y , z ) = o ( x 0 , y 0 , z 0 ) g ( x x 0 , y y 0 , z z 0 ) d x 0 d y 0 d z 0
= o ( x , y , z ) * g ( x 0 , y 0 , z 0 ) ,
g ( x , y , z ) = exp ( i 2 π λ x 2 + y 2 + z 2 ) x 2 + y 2 + z 2 .
G ( u , v , w ) = F [ g ( x , y , z ) ] = 1 4 π 2 ( u 2 + v 2 + w 2 1 / λ 2 ) ,
f ( x , y , z ) = 1 4 π 2 O ( u , v , w ) u 2 + v 2 + w 2 1 / λ 2 exp { i 2 π ( u x + u y + w z ) } d u d v d w ,
f ( x , y ) | z = R = 1 4 π 2 O ( u , v , w ) u 2 + v 2 + w 2 1 / λ 2 exp { i 2 π ( u x + u y + w R ) } d u d v d w .
f ( x , y ) | z = R = i 4 π O ( u , v 1 / λ 2 u 2 v 2 ) 1 / λ 2 u 2 v 2 exp ( i 2 π R 1 / λ 2 u 2 v 2 ) × exp { i 2 π ( u x + u y ) } d u d v .
F ( u , v ) | z = R = i 4 π O ( u , v 1 / λ 2 u 2 v 2 ) 1 / λ 2 u 2 v 2 exp ( i 2 π R 1 / λ 2 u 2 v 2 )
O ( u , v , w ) = F [ o ( x , y , z ) comb ( x δ , y δ , z δ ) ] = 1 W O ( u , v , w ) * comb ( u W , v W , w W ) ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.