Cylindrical angular spectrum using Fourier coefficients of point light source and its application to fast hologram calculation

Seungtaik Oh; Il Kwon Jeong

doi:10.1364/OE.23.029555

1. Introduction

In recent years, computer-generated hologram(CGH) which is a digitally synthesized hologram by a numerical simulation of a diffracted wave field, has been one of the most attractive fields in holography. One advantage of CGH is that one can generate a digital hologram for a 3D virtual object. Generally, the object wave recorded on a CGH is optically reconstructed by a spatial light modulator(SLM).

One important issue in CGH is the narrow field of view. In order to widen the field of view or to provide a wider range of object view in reconstruction, the size of a CGH needs to be enlarged and at the same time, the pixel size of the CGH should be decreased to record the object wave correctly by the Nyquist sampling condition. As an extreme case, if one wants to record the object wave in 180°, then the pixel size of a CGH should go to the half of the wavelength and a significantly large number of pixels are needed.

A cylindrical hologram can be an effective solution to record and reconstruct the object wave in 360° [1, 2]. Cylindrically-arranged planar CGHs were used to increase the field of view [3, 4]. Cylindrical CGHs to capture the object wave in 360° have been studied by many authors. Yamaguchi et al. [5] proposed a fast calculation method for a cylindrical CGH based on segmentation and lookup table. Jackin et al. [6] introduced another Fourier transform-based high speed calculation method which can be understood as an extension of the angular spectrum method to the cylindrical coordinates system. In particular, Sando et al. [7,8] calculated the 3D Fourier spectrum to record the object wave in all directions and the method had been developed to a fast cylindrical CGH generation method based on 1D convolution and Fourier transform [9].

In this paper we are going to introduce a new simple analytic formula for the Fourier coefficients of the 3D field distribution of a point light source to generate a cylindrical angular spectrum which can capture a 3D Fourier spectrum of the object wave in all horizontal directions. Conceptually, the cylindrical angular spectrum can be understood as a cylindrical version of the omnidirectional spectral method in [7, 8]. Our Fourier coefficient formula for a point light source is derived from an intuitive observation about a fundamental property of point light source and a surface integral expression over a sphere in the 3D Fourier space for the 3D field distribution of light. This derivation can provide a more natural and clearer understanding of capturing the angular spectrum of the object wave in all directions since the Fourier coefficients can be directly calculated for frequency vectors lying on the entire surface of a sphere in the 3D Fourier space. In contrast, the previous method in [7, 8] computed a basic frequency spectrum defined only on the hemisphere by using a diffraction formula to a plane, z = z₀ and merged some overlapping basic frequency spectrums to obtain a single spectrum on the entire sphere. Therefore, our Fourier coefficient formula is more frequency-oriented and allows us to generate a cylindrical angular spectrum in a more efficient manner without any spectrum merging and expensive frequency vector rotations on a hemisphere which is required in the previous method to align a basic spectrum to the reference propagation axis.

We will also propose a generalized frequency-based occlusion culling process tightly coupled with the generation of a cylindrical angular spectrum. Our occlusion culling is based on an intersection computation between a ray and the object geometry and so it is more effective for arbitrary complex objects which cannot be handled by the previous method in [8].

One important advantage of the cylindrical angular spectrum is that one can synthesize a planar CGH along any given horizontal direction at a high speed. The procedural steps to synthesize a planar CGH from a cylindrical angular spectrum are explained in details. We will finally present an experimental interactive hologram result based on realtime synthesis of a CGH from a cylindrical angular spectrum. Since the diffraction to a target plane can be calculated by a single operation of the fast Fourier transform(FFT), the CGH synthesis can be done almost in realtime.

2. Cylindrical angular spectrum generation

2.1. Fourier coefficients of the 3D field generated by a point source

Let us first consider a field Ũ(x) of wavelength λ generated by a point source at the origin in the 3D space. Since Ũ(x) is a 3D monochromatic field, Ũ(x) can be written as a surface integral over the sphere of radius 1/λ denoted by S_1/λ (see [10]):

\tilde{U} (x) = \tilde{U} (x, y, z) = \int_{k \in ℝ^{3}} δ_{S_{1 / λ}} (k) A (k) exp (j 2 π k \cdot x) d k

= \int_{k \in S_{1 / λ}} A (k) exp (j 2 π k \cdot x) d S,

where the δ_{S_1/λ}(·) is the delta function over S_1/λ and k ∈ S_1/λ is a frequency vector. Note that Eq. (1) is a triple integral over the full 3D Fourier space and Eq. (2) is a surface integral over S_1/λ in the 3D Fourier space. Our intuitive observation is that the light generated by a point source radiates uniformly in all directions with a uniform spectrum of propagating plane waves, which implies that all the Fourier coefficients need to be constant for Ũ(x) and we may write A(k) ≡ 1 in Eq. (2).

If we translate Ũ(x) by shifting the source to a position p, the translated field U_p(x) := Ũ(x − p) is written as

U_{p} (x) = \int_{k \in S_{1 / λ}} A_{p} (k) exp (j 2 π k \cdot x) d S,

where

A_{p} (k) = exp (- j 2 π k \cdot p) .

Therefore, our analytic formula of the Fourier coefficient of an arbitrary point light source with respect to a frequency vector is given by Eq. (4).

Assume a given object is sampled by points {p_i} and each point p_i has a complex amplitude Ū_{p_i}. Then, the Fourier coefficient A(k) for the object wave U (x) is obtained by adding up Ū_{p_i} A_{p_i} (k) for all sampling points:

A (k) = \sum_{i} {\bar{U}}_{p_{i}} A_{p_{i}} (k) = \sum_{i} {\bar{U}}_{p_{i}} exp (- j 2 π k \cdot p_{i}) .

Now, we define cylindrical angular spectrum by the set of pairs of a frequency vector and the corresponding Fourier coefficient {(k_i, A(k_i))}, where k_i ∈ S_1/λ and a uniform scaling of k_i lies on a cylinder in the 3D Fourier space.

In order to generate a cylindrical angular spectrum for U(x), we consider a cylindrical grid structure in the 3D Fourier space(see Fig. 1). The grid structure has uniform grid spacings, that is, a horizontal angular spacing Δθ and a vertical spacing Δβ. Each grid point and grid value refer to a frequency vector and the corresponding Fourier coefficient, respectively. A grid point v with a 2D index (i, j) is given by $v = (\frac{1}{λ} cos (i Δ θ), - \frac{Δ β N_{y}}{2} + j Δ β, γ = \frac{1}{λ} sin (i Δ θ))$ , where N_y is the number of the vertical number of pixels in the cylindrical grid structure, and the associated frequency vector k with the grid point v is expressed as $k = \frac{1}{λ} \frac{v}{‖ v ‖}$ . A cylindrical angular spectrum is generated by computing the Fourier coefficient A(k) from Eq. (5) for each frequency vector associated with a grid point of the cylindrical grid structure.

Fig. 1 Cylindrical grid structure in the 3D Fourier space for cylindrical angular spectrum.

Download Full Size | PDF

Since the object wave is propagating to the outward direction, we impose a constraint on the frequency vectors involved in the generation of the cylindrical angular spectrum:

k \cdot n_{p} > 0,

where n_p is the outward normal vector of the object at p.

2.2. Ray-based occlusion culling in the 3D Fourier space

Occlusion culling is one of the most important factors to enhance the quality of a hologram. We introduce a generalized occlusion culling method in the 3D Fourier space based on identification of a frequency vector with a ray as seen in [8]. It can be assumed that each object sampling point emits a set of rays whose origin and direction are the position of the sampling point and a frequency vector, respectively. Then, we are going to apply a ray-object intersection test for a occlusion culling in the 3D Fourier space. More precisely, if a ray given by an object point p and a frequency vector k intersects with the object, the corresponding Fourier coefficient A_p(k) is ignored in the generation of a cylindrical angular spectrum by setting

A_{p} (k) = 0 .

An additional surface representation of the object such as a triangular mesh or an algebraic equation is needed for a proper ray-object intersection test. In this work, we take a triangular mesh for its convenience and popularity. Then, our ray-object intersection test can be reduced to a ray-triangle intersection test. If the given ray intersects with a triangle of the mesh, the ray intersects with the object. If the ray does not intersect with any triangle of the mesh, the ray does not intersect with the object. The ray-triangle intersection test is a fundamental problem of computer graphics [11]. The detailed explanation about how to determine if a ray intersects with a triangle is given in the appendix.

For comparison with the previous work, it is noted that the frontward point selection algorithm described in [8] might be ambiguous and only effective for a simple object such as a sphere. For a spherical object, the frontward points can be selected simply by the hemisphere determined by the ray direction. However, this simple selection method does not work for general cases. In order to select the frontward points in a general complex object, one needs to determine if two interesting points are on the same ray associated with the propagating direction, but such condition is hard or ambiguous to determine in practical situations, and moreover comparisons for all points can be very inefficient in computation. Thus, the ray-object intersection computation for an occlusion culling is inevitable and most effective for general complex objects.

The ray-object intersection computation is a very expensive operation and it is common to use an acceleration structure to speed up the computation. Since the intersection computation has a nice structure for parallel processing, one can take a maximum efficiency by making use of a high performance parallel computing structure such as modern programmable graphics processing unit(GPU).

3. Fourier-based CGH synthesis from cylindrical angular spectrum

Now it is assumed that we have a cylindrical angular spectrum for the object wave. We will explain here how to compute the diffraction of the object wave to a target plane based on our cylindrical angular spectrum, or equivalently, how to synthesize a planar CGH at an arbitrary horizontal view angle from our cylindrical angular spectrum. Let Δp be the pixel size of a CGH to be synthesized and let Δϕ = sin⁻¹(λ/(2Δp)) be the diffraction angle determined by Δp. Our CGH synthesis consists of several steps as follows:

Identify a patch of the cylindrical angular spectrum required to synthesize a CGH.
Rotate the angular spectrum patch around the β-axis so that the center of the rotated patch exists on the +γ-axis.
Project the rotated angular spectrum patch to the αβ-plane.
Propagate the projected angular spectrum by a distance given by the hologram plane.
Apply the inverse Fourier transform to obtain the desired CGH.

To explain the steps in more detail, let c and θ be the center of the hologram plane and the view angle, respectively, as shown in Fig. 2(a).

Fig. 2 CGH synthesis at view angle θ and distance d(a) and rotation of cylindrical angular spectrum patch(b). P colored in gray is an angular spectrum patch determined by the given CGH and P′ is the rotated angular spectrum patch whose center is aligned to the +γ-axis. k is a frequency vector in P and k′ is the rotated frequency vector in P′.

Download Full Size | PDF

In step 1, we identify a curved rectangular patch P of the cylindrical angular spectrum to generate a planar CGH. The patch is represented by a 2D index array of grid points [i₀, i₁] × [j₀, j₁], where i and j represent the horizontal and vertical indices of the cylindrical grid structure, respectively. Let (i_c, j_c) be the 2D index of the closest grid point to c and let Δi and Δj be the numbers of grid points such that Δθ × Δi ≈ Δϕ and $Δ β \times Δ j \approx \frac{1}{λ} {tan}^{- 1} (Δ ϕ)$ . Then, we set i₀ = i_c − Δi, i₁ = i_c + Δi and j₀ = j_c − Δj, j₁ = j_c + Δj. In our case, Δϕ is relatively small and so P occupies a quite small portion of the cylindrical grid structure.

In step 2, each frequency vector in the angular spectrum patch P is rotated by −θ with respect to the β-axis so that the center of the rotated angular spectrum P′ exists on the +γ-axis as shown in Fig. 2(b). Note that the Fourier coefficient data is unchanged during the rotation and the rotation computation can be done at a high speed since the number of the frequency vectors in the patch is small.

In step 3, the rotated angular spectrum patch P′ is projected to the αβ-plane. Each frequency vector k = (α, β, γ) in P′ is projected to k′ = (α, β, 0). At the same time, each Fourier coefficient is scaled by 1/γ due to the Jacobian of the projection map from S_1/λ to the αβ-plane [10]:

A^{'} (k^{'}) = \frac{1}{γ} A (k),

where A and A′ are the Fourier coefficient and the projected Fourier coefficient, respectively. In fact, the end points of the projected frequency vectors, generally, do not match exactly with the grid points of a uniform grid on the αβ-plane which is required to synthesize the desired CGH. So, by finding the projected Fourier coefficient of the nearest projected frequency vector or interpolating the projected Fourier coefficients of some neighboring projected frequency vectors, an angular spectrum A′_z=0(α, β) defined on a uniform grid of the αβ-plane is obtained. The sampling condition implies that if the CGH has N × N pixels and the entire domain of the uniform grid is [−α_max, α_max] × [−β_max, β_min], then α_max = β_max = 1/(2Δp) and the grid spacings are Δα = 2α_max/N and Δβ = Δα. It is also noted that for a correct CGH synthesis, the grid structure determined by the projected frequency vectors needs to be as fine as the uniform grid structure determined by the CGH to be synthesized:

\tilde{Δ α} \leq Δ α

and

\tilde{Δ β} \leq Δ β

, where

\tilde{Δ α}

and

\tilde{Δ β}

are the average grid spacings of the projected frequency vectors.

In step 4, the angular spectrum defined on the uniform grid of the αβ-plane is propagated by the distance d = ‖c‖ following the well known angular spectrum propagation formula [10]:

A_{z = d}^{'} (α, β) = A_{z = 0}^{'} (α, β) exp (j 2 π d \sqrt{\frac{1}{λ^{2}} - α^{2} - β^{2}}) .

In step 5, the desired CGH is obtained by applying the inverse Fourier transform to the angular spectrum A′_z=d associated with the hologram plane.

4. Results

In this section we verify our formula Eq. (4) for the Fourier coefficient of a point light source and provide an experimental result of an interactive CGH synthesis with a user key-input based on our cylindrical angular spectrum and planar CGH synthesis algorithm.

4.1. Verification of Fourier coefficient formula for point light source

We compute a field generated by a point source at (0, 0, 1) on the xy-plane based on Eq. (4). Note that Eq. (4) is the Fourier coefficients defined on the sphere S_1/λ and in order to compute an angular spectrum A(α, β) on the αβ-plane, a factor 1/γ needs to be applied to Eq. (4)(see Eq. (8)). The field on the xy-plane is obtained by applying an inverse FFT to A(α, β). The bipolar intensity pattern of the field on the xy-plane is the same as a zone plate pattern by a single spherical wave as shown in Fig. 3(a) and the point light is correctly restored at z = 1 by numerical reconstruction. In addition, some focus analyses are performed for further verification. A focus estimation algorithm based on frequency analysis in [12] computes the focus distance at z = 1.00124. A focus measure-based analysis in [13] also shows that the field has the focus at z = 1 as illustrated in Fig. 3(b).

Fig. 3 Bipolar intensity pattern at z = 0(a) and focus analysis(b) for the field generated by a point source at (0, 0, 1). In the graph of (b), the peak point represents the focus distance.

Download Full Size | PDF

4.2. Example of interactive CGH

First, the wavelength is set to λ = 633nm. Construct a cylindrical grid structure in the 3D Fourier space as seen in Fig. 1. The horizontal angular and the vertical spacings of the grid structure are set to Δθ = 3.86353 × 10⁻⁵(rad) and Δβ = 61.0352(1/m) and then, we have 162,628 × 2,048 grid points in the cylindrical grid structure.

With the cylindrical grid structure, a cylindrical angular spectrum is generated for a bunny mesh model with 34,834 points. Note that we have 162,628 × 2,048 frequency vectors where to calculate the Fourier coefficients. For each frequency vector k, the Fourier coefficient A(k) is calculated by adding up the elementary Fourier coefficients based on Eq. (5) over some object points with respect to which the ray direction given by k is outward(see Eq. (6)). In adding up, for occlusion culling, the elementary Fourier coefficient is set to zero (see Eq. (7)) if the ray determined by k and an object point intersects with the object mesh.

In our implementation, we make use of a free GPU-based acceleration structure OptiX provided by NVIDIA for the ray-object intersection computation which is the core process of our frequency-based occlusion culling. With the GPU-based approach, approximately 30% performance improvement has been achieved in the generation of a cylindrical angular spectrum compared to a CPU-based approach. Even though the occlusion culling performance has been improved thanks to the GPU-based acceleration approach, the total computation for the cylindrical angular spectrum generation takes a while, about 22 hours, which is a limitation of the current work. An occlusion culling result for a simple test example is shown in Fig. 4.

Fig. 4 A result of frequency-based occlusion culling. Before occlusion culling(a) and after occlusion culling(b). The front object(left lower) is focused for both cases.

Download Full Size | PDF

The CGH to be synthesized is set to have 2,048 × 2,048 pixels and 8μm pixel size. The horizontal view angle is changed by a user key input and so holographic image can be rotated by pressing keys to increase or decrease the horizontal view angle. Given a horizontal view angle, the Fourier-based CGH synthesis process described in Section 3 is performed to generate a CGH. In this example, the distance to the CGH is fixed to 20cm. For any view angle, a CGH is generated almost in realtime. A more detailed timing information for our CGH synthesis is presented in Table 1. The most expensive step in the CGH synthesis is data copying between CPU and GPU, and we have optimized the data copying with pinned memories. The numerical reconstruction results for horizontal view angles −20°, 70°, 160° and 250° are shown in Fig. 5. For a horizontal view angle −20°, numerical reconstruction results at different focus distances are presented as well in Fig. 6.

Table 1. Computation time to synthesize a CGH with 2,048 × 2048, pixels from a cylindrical angular spectrum with 162,628 × 2,048 grid points for a bunny model with 34,834 points.

View Table

Fig. 5 Numerical reconstructions for view angle −20°, 70°, 160°, 250° from cylindrical angular spectrum(see Visualization 1 and Visualization 2 for 360° and interactive results).

Download Full Size | PDF

Fig. 6 Numerical reconstructions with different focus distance for view angle −20°. The chest(a) and the nose(b) are focused, respectively.

Download Full Size | PDF

5. Conclusion

We have proposed a new simple analytic formula for the Fourier coefficients of the 3D distribution of a point light source to generate a cylindrical angular spectrum. Our Fourier coefficient formula is based on an intuitive observation that a point light radiates uniformly in all directions with a uniform spectrum. Our formula is defined over all frequency vectors lying on the entire sphere in the 3D Fourier space and so it is more natural and computationally more efficient for all around recording of the object wave than the previous formula which is defined on the hemisphere and is indirectly derived from a diffraction formula to a plane, z = z₀. Our formula has been also verified with a various supportive data. In addition, a generalized ray-based occlusion culling method in the 3D Fourier space has been introduced for arbitrary complex objects. As a result, an interactive CGH synthesis with a key-input to rotate the holographic image horizontally is presented and a 2048×2048 CGH for a 34,834 points model can be synthesized almost in real time(≈ 32ms) at an arbitrary horizontal view angle.

Appendix

We will explain here how to determine if a ray intersects a triangle in ℝ³. Let r₀ and k be the origin and the direction of the ray, respectively. Let p₀, p₁ and p₂ be the positions of the three vertices of the triangle. First, we compute the intersection point x between the ray and the plane T determined by the triangle. One may write x = r₀ + tk for t ∈ ℝ. Let n be the normal of the triangle. Then, we have $t = - \frac{r_{0} \cdot p_{0}}{k \cdot n}$ since (x − p₀) · n = 0. Second, we need to determine if x is inside the triangle. To do this, let v = p₁ − p₀ and w = p₂ − p₀. Since {v, w} is a basis for the plane T, we can write x − p₀ = av + bw. The coefficients a and b are computed by

\begin{array}{l} [\begin{array}{l} a \\ b \end{array} & = & {[\begin{matrix} v \cdot v & v \cdot w \\ v \cdot w & w \cdot w \end{matrix}}^{- 1} & [\begin{matrix} (x - p_{0}) \cdot v \\ (x - p_{0}) \cdot w \end{matrix} . \end{array}

Finally, x is inside the triangle, or equivalently, the ray intersects with the triangle if and only if a + b ≤ 1, a ≥ 0, b ≥ 0.

Acknowledgments

The authors would like to thank all anonymous reviewers for their valuable comments and suggestions. This work was supported by the ‘Cross-Ministry Giga KOREA Project’ of the Ministry of Science, ICT and Future Planning, South Korea[ GK15C0100, Development of Interactiveand RealisticMassive Giga-ContentTechnology].

References and links

1. T. Jeong, “Cylindrical holography and some proposed applications,” J. Opt. Soc. Am. 57, 1396–1398 (1967). [CrossRef]

2. O. Soares and J. Fernandes, “Cylindrical hologram of 360° field of view,” Appl. Opt. 21, 3194–3196 (1982). [CrossRef] [PubMed]

3. J. Hahn, H. Kim, Y. Lim, G. Park, and B. Lee, “Wide viewing angle dynamic holographic stereogram with a curved array of spatial light modulators,” Opt. Express 16, 12372–12386 (2008). [CrossRef] [PubMed]

4. F. Yaraş, H. Kang, and L. Onural, “Circular holographic video display system,” Opt. Express 19, 9147–9156 (2011). [CrossRef]

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

6. B. Jackin and T. Yatagai, “Fast calculation method for computer-generated cylindrical hologram based on wave propagationin spectral domain,” Opt. Express 18, 25546–25555 (2010). [CrossRef] [PubMed]

7. Y. Sando, D. Barada, and T. Yatagai, “Fast calculation of computer-generated holograms based on 3-d fourier spectrum for omnidirectional diffraction from a 3-d voxel-based object,” Opt. Express 20, 20962–20969 (2012). [CrossRef] [PubMed]

8. Y. Sando, D. Barada, and T. Yatagai, “Hidden surface removal of computer-generated holograms for arbitrary diffraction directions,” Appl. Opt. 52, 4871–4876 (2013). [CrossRef] [PubMed]

9. Y. Sando, D. Barada, B. Jackin, and T. Yatagai, “Fast calculation method for computer-generated cylindrical holograms based on the three-dimensional fourier spectrum,” Opt. Lett. 38, 5172–5175 (2013). [CrossRef] [PubMed]

10. L. Onural, “Exact solution for scalar diffraction between tilted and translated planes using impulse functions over a surface,” J. Opt. Soc. Am. A 28, 290–295 (2011). [CrossRef]

11. P. Shirley and R. Morley, Realistic Ray Tracing, 2nd ed. (A K Peters, 2003).

12. S. Oh, C.-Y. Hwang, I. K. Jeong, S.-K. Lee, and J.-H. Park, “Fast focus estimation using frequency analysis in digital holography,” Opt. Express 22, 28926–28933 (2014). [CrossRef] [PubMed]

13. P. Memmolo, C. Distante, M. Paturzo, A. Finizio, P. Ferraro, and B. Javidi, “Automatic focusing in digital holography and its application to stretched holograms,” Opt. Lett. 36, 1945–1947 (2011). [CrossRef] [PubMed]

Name	Description
Visualization 1: MOV (7908 KB)	360 degree reconstruction result
Visualization 2: MOV (27950 KB)	realtime capture video

Procedure	CPU	GPU(GTX690)

Identify angular spectrum patch and initialize	11 ms	11.00 ms
Rotate angular spectrum patch	4 ms	1.90 ms
Project rotated angular spectrum patch	82 ms	0.93 ms
Propagate angular spectrum to a target plane	30 ms	0.71 ms
Inverse FFT to obtain the final CGH	242 ms	2.51 ms

Copy memory between CPU and GPU	-	16 ms

Total time for CGH synthesis	369 ms	32.05 ms

Cylindrical angular spectrum using Fourier coefficients of point light source and its application to fast hologram calculation

Abstract

1. Introduction

2. Cylindrical angular spectrum generation

2.1. Fourier coefficients of the 3D field generated by a point source

2.2. Ray-based occlusion culling in the 3D Fourier space

3. Fourier-based CGH synthesis from cylindrical angular spectrum

4. Results

4.1. Verification of Fourier coefficient formula for point light source

4.2. Example of interactive CGH

5. Conclusion

Appendix

Acknowledgments

References and links

Supplementary Material (2)

Cited By

Figures (6)

Tables (1)

Equations (10)

Optics Express