Scaled diffraction calculation between tilted planes using nonuniform fast Fourier transform

Chenliang Chang; Jun Xia; Jun Wu; Wei Lei; Yi Xie; Mingwu Kang; Qiuzhi Zhang

doi:10.1364/OE.22.017331

1. Introduction

Holography is the most prominent technology for three-dimensional display due to its capability for reconstructing all of the information of an object [1]. One of the important works for holographic display is the calculation of computer generated hologram (CGH). The conventional method for CGH calculation is the point-based method, in which the 3D object is sampled into a large amount of points and the CGH is obtained by diffraction calculation from all of the points [2]. Although some methods for improving the calculation speed have been proposed, the point-based method is still time consuming. In order to speed-up the calculation, researchers have developed the polygon-based method. The object is represented by a set of planes rather than points, and the CGH is obtained by calculating the diffraction from each plane, in which the fast Fourier transform (FFT) is used to accelerate the calculation [3–5]. Based on this theory, several advanced methods have been proposed in the past decade. For example, L. Ahrenberg proposed a method to compute the CGH by using an analytical light transport model to avoid calculation FFT for every plane [6]. Y. Pan proposed a fast polygon-based method based on two-dimensional Fourier analysis of 3D affine transformation to reconstruct the 3D scene with the solid effect [7].

For the polygon-based method, since a three-dimensional object is divided into a few planes which are not parallel to the hologram plane, the key task is the diffraction calculation between the tilted plane and the hologram plane. D. Leseberg first proposed a method to calculate diffraction of objects on tilted planar using a coordinate transformation and a subsequent multiplication of a quadratic phase [8]. K. Matsushima and S. D. Nicola proposed a method to calculate the diffraction between tilted planes based on the angular spectrum and the coordinate rotation in the Fourier domain [9, 10]. The angular spectrum based method uses the coordinate rotation to transform the Fourier spectrum from the tilted plane to the reference plane, which get nonlinear sampling rates on the reference plane. Therefore an interpolation is implemented on the reference plane. Moreover, by using the FFT to calculate the Fourier spectrum, the sampling rates on both of the spatial domain and the frequency domain is restricted by the Nyquist theory, which means that one can only calculate the diffraction from tilted plane to hologram plane with the fixed sampling rates, so the size of the three-dimensional object is restricted.

In this paper, we proposed a simple method to calculate the diffraction from a tilted plane with variable sampling rates to the hologram plane. The nonuniform fast Fourier transform (NUFFT) algorithm is used to calculate the diffraction instead of FFT. The use of NUFFT has two advantages. Firstly, it can directly calculate the spectrum values on the tilted plane with a nonuniform sampled spectrum coordinate, so the interpolation operation after coordinate rotation is not necessary. Secondly, the NUFFT can calculate a scaled Fourier transformation which can easily realize an image scaling.

2. Method

2.1 Calculate diffraction from tilted plane using NUFFT

In order to calculate the diffraction from a tilted plane, we first set up a reference plane which is paralleled to the hologram plane. As shown in Fig. 1(a), the reference plane and the tilted plane share the same origins. Let us assume that f_h (x_h, y_h) represents the light distribution on the hologram plane and f_t (x_t, y_t) represents the light distribution on the tilted plane respectively.

Fig. 1 Schematic diagram of diffraction calculation of tilted plane. (a) Optical locations of tilted plane, reference plane and hologram plane. (b) Conventional calculation between tilted plane and reference plane. (c) NUFFT based method of calculation between tilted plane and reference plane.

Download Full Size | PDF

The flow chart of the calculation can be seen in Fig. 1(b). We assume the uniform sampling on the reference plane. The relationship in Fourier spectrum domain between these two planes, which can be expressed as

[\begin{matrix} u_{t} \\ v_{t} \\ w_{t} (u_{t}, v_{t}) \end{matrix}] = T \cdot [\begin{matrix} u_{r} \\ v_{r} \\ w_{r} (u_{r}, v_{r}) \end{matrix}] = [\begin{matrix} a_{1} & a_{2} & a_{3} \\ a_{4} & a_{5} & a_{6} \\ a_{7} & a_{8} & a_{9} \end{matrix}] \cdot [\begin{matrix} u_{r} \\ v_{r} \\ w_{r} (u_{r}, v_{r}) \end{matrix}]

Where T is the transformation matrix and (u_t, v_t) is the coordinate in spectrum domain on the tilted plane. The parameters in matrix T is determined by the angle of plane rotation [9]. According to Eq. (1), the spectrum coordinate of tilted plane can be further rewritten as

\begin{array}{l} u_{t} = α (u_{r}, v_{r}) = a_{1} u_{r} + a_{2} v_{r} + a_{3} w_{r} (u_{r}, v_{r}) \\ v_{t} = β (u_{r}, v_{r}) = a_{4} u_{r} + a_{5} v_{r} + a_{6} w_{r} (u_{r}, v_{r}) \end{array}

Where, w_r = (λ⁻²-u_r²-v_r²)^1/2. Equation (2) is a nonlinear function which will get nonuniform sampling rates on the tilted plane as shown in Fig. 1(b). Conventionally, an interpolation is implemented in order to get uniform sampled G_t (u_t, v_t) [see the dotted area in Fig. 1(b)]. Here we propose to use the nonuniform fast Fourier transform (NUFFT) to directly calculate G_t (u_t, v_t).

NUFFT is an algorithm to calculate the Fourier transform when the sampling rate in either the spatial domain or the frequency domain is nonuniform [11–13]. Generally the NUFFT can be categorized into two types: the first type is the calculation from nonuniform sampled points to uniform sampled points (the type 1), and the second type is the calculation from uniform sampled points to nonuniform sampled points (the type 2). In our situation, the calculation from uniform sampled plane f_t (x_t, y_t) to nonuniform sampled spectrum G_t (u_t, v_t) is type 2. So by using the NUFFT algorithm, the calculation of spectrum on the tilted plane can be expressed as:

\begin{matrix} G_{t} (u_{t}, v_{t}) = G_{t} [m_{u} Δ (u_{t}), m_{v} Δ (v_{t})] \\ = N U F F T 2 [f_{t} (m_{x} Δ x, m_{y} Δ y)] \\ = N U F F T 2 [f_{t} (x_{t}, y_{t})] \end{matrix}

Where NUFFT2 represents the second type of NUFFT. m_u, m_v, m_x and m_y are integers with range of [-N/2, N/2] where N is the sampling resolution in each domain. Δ(u_t) and Δ(v_t) are the sampling rates in the spectrum domain of tilted plane which are nonuniformly sampled. The brackets means the sampling rate is associated with the spectrum coordinate. The sampling rates on the tilted plane Δx and Δy are constant and uniform. Once the spectrum on the tilted plane is obtained, the spectrum on the reference plane can be easily obtained by

G_{r} (u_{r}, v_{r}) = G_{r} (α^{- 1} (u_{t}, v_{t}), β^{- 1} (u_{t}, v_{t})) = G_{t} (u_{t}, v_{t})

Where α⁻¹ and β⁻¹ represent the inverse function. Then the spectrum on the hologram plane can be calculated by the propagation of plane wave as [14]

G_{h} (u_{h}, v_{h}) = G_{r} (u_{r}, v_{r}) \cdot \exp (2 π i z \sqrt{\frac{1}{λ^{2}} - u_{r}^{2} - v_{r}^{2}})

Where, z is the propagation distance between the reference plane and the hologram plane, λ is the wavelength of the light. The NUFFT is used to directly calculate the nonuniform sampled spectrum on the tilted plane so the interpolation in the conventional method is not necessary. In the following part, we will introduce the three-dimensional scaling that is a byproduct of the NUFFT.

2.2 Scaled diffraction calculation between tilted planes based on NUFFT

For simplicity, the following analysis is presented in the one dimensional form. Figure 2(a) shows the relationship of sampling rates in the conventional angular spectrum based calculation. We assume that the sampling rate on the hologram is ΔH, the sampling rate on tilted plane is ΔT. The bandwidth of the spectrum on tilted plane is L_t = 1/ΔT according to the Nyquist theory. After the coordinate rotation in spectrum domain, the bandwidth of the spectrum on the reference plane is L_r = (1/cosθ)L_t, here θ is the tilting angle between the tilted plane and reference plane. Then the spectrum is propagated to the hologram plane with the same bandwidth L_h = L_r. Therefore the sampling rate of the hologram is ΔH = 1/L_h = cosθ ΔT. This indicates that in the conventional method the sampling rates are constraint, which means once the sampling rate on the hologram is known (this is usually the sampling rate of the spatial light modulator), the size of the image on the tilted plane is determined by ΔT = ΔH/cosθ.

Fig. 2 Schematic diagrams of scaled diffraction calculation between tilted planes. (a) Conventional calculation process based on angular spectrum and interpolation between tilted planes with the constraint sampling rates. (b) Scaled diffraction calculation based on NUFFT from tilted plane with variable sampling rate when R<1. (c) Scaled diffraction calculation based on NUFFT from tilted plane with variable sampling rate when R>1.

Download Full Size | PDF

In order to calculate the image on tilted plane with different image size, we used the NUFFT instead of FFT to overcome the limitation imposed by sampling theory. According to Eq. (5), the distribution on the hologram using NUFFT can be written as

\begin{matrix} f_{h} (x_{h}) = N U F F T 1 {G_{r} (u_{r}) \cdot \exp (2 π i z \sqrt{\frac{1}{λ^{2}} - u_{r}^{2}})} \\ = N U F F T 1 {G_{t} (u_{t}) \cdot \exp (2 π i z \sqrt{\frac{1}{λ^{2}} - u_{r}^{2}})} \\ = N U F F T 1 {N U F F T 2 [f_{t} (x_{t})] \cdot \exp (2 π i z \sqrt{\frac{1}{λ^{2}} - u_{r}^{2}})} \end{matrix}

The second and third equal sign in Eq. (6) are deduced according to Eq. (4) and Eq. (3) respectively. Let us assume that the sampling rate we want to calculate on the tilted plane is ΔT’ = RΔT (R≠1), in which R is a scaling factor. And we first consider the situation R<1. In this paper, we used Greengard and Lee’s NUFFT, which is based on the oversampled gridding (for more details, see [13]). The calculation of the type 1 NUFFT is

F (n) = N U F F T 1 [f (x_{m})] = \sum_{m} f (x_{m}) \exp (- i n x_{m})

Where n is integral value with uniform sampling, x_m = mΔ(x_m) is nonuniform sampled and its range must be confined in the interval of [-π, π]. Similarly, the calculation of type 2 of the NUFFT is

F (x_{n}) = N U F F T 2 [f (m)] = \sum_{m} f (m) \exp (- i x_{n} m)

Where m is integral value with uniform sampling and x_n = nΔ(x_n) is nonuniform sampled. The range of x_n is also in the interval of [-π, π]. The expression of the spectrum G_t of the tilted plane is

\begin{matrix} G_{t} (u_{n}) = \sum f_{t} (x_{m}) \cdot \exp (- i 2 π m R Δ T u_{n}) \\ = \sum f_{t} (x_{t}) \cdot \exp (\frac{- 2 i π R m u_{n}}{L_{t}}) \end{matrix}

Where m and n is the index number of sampling points in spatial and frequency domains respectively, u_n = nΔ(u_n) is the coordinate of spectrum which is nonuniform sampled. If R<1, we have u_n’ = (2πRu_n/L_t)-π which is in the interval of [-π, π], so Eq. (9) can be calculated by NUFFT2 according to Eq. (8). Therefore, a variable sampling rate ΔT’ = RΔT is achieved. The detailed calculation process is shown in Fig. 2(b), the red points on tilted plane means the scale effect of sampling rate change by R<1. It should be noted that here the NUFFT2 performs two functions, one is the calculation of scale factor R and another is the direct calculation of nonuniform sampled spectrum.

When R>1, the NUFFT2 [Eq. (8)] cannot be used to calculate Eq. (9), because u_n’ = (2πRu_n/L_t)-π exceed the interval of [-π, π]. So here NUFFT2 only perform the calculation of nonuniform sampled spectrum. And due to Nyquist sampling theory, the bandwidth on the spectrum domain is changed to L_t’ = (1/R)L_t after the rotation. The bandwidth of spectrum on the hologram plane is also changed to L_h’ = (1/R)L_h, so now the distribution on hologram plane can be written as

\begin{matrix} f_{h} (m) = \sum G_{h} (u_{h}) \cdot \exp (i 2 π m Δ H n Δ u_{h}') \\ = \sum G_{h} (u_{h}) \cdot \exp (\frac{i 2 π m Δ H n L_{h}'}{N}) \\ = \sum G_{h} (u_{h}) \cdot \exp (\frac{i 2 π m Δ H n L_{h}}{R N}) \\ = \sum G_{h} (u_{h}) \cdot \exp (\frac{i 2 π m n}{R N}) \end{matrix}

Now if we define n’ = -(2πn/RN)-π, then n’ is in the interval of [-π, π] because of R>1, so Eq. (10) can be calculated by NUFFT1 according to Eq. (7). The whole calculation process is shown in Fig. 2(c). The red points means the sampling rate is changed compared to the conventional angular spectrum based method. It can be seen that due to the scale of the sampling on tilted plane (R>1), the bandwidth of spectrum in all of the three planes is also scaled, and then by using the NUFFT1, the hologram with the sampling rate ΔH is calculated.

In conclusion, the whole calculation process can be simply generalized as three steps: 1. the calculation from tilted plane to its spectrum domain by NUFFT2, here the tilted plane is uniformly sampled and its spectrum domain is nonuniformly sampled, 2. the propagation from the spectrum domain of the tilted plane to the spectrum domain of the hologram plane, 3. the calculation from spectrum domain to the hologram plane by NUFFT1 (R>1) or FFT (R<1), here both of the hologram plane and its spectrum domain is uniformly sampled. Figure 3 shows the 2D diagram of coordinate grids on each domain in the calculation. The samplings on both of the tilted and hologram plane are uniform and the sampling on the spectrum domain of tilted plane is distorted due to the nonlinear coordinate rotation. The calculation from tilted plane to its spectrum is by NUFFT. The scaled diffraction calculation from tilted plane with sampling rates ΔT’ = RΔT can also be achieved by using NUFFT at the same time.

Fig. 3 Schematic diagram of coordinate grids on each domain.

Download Full Size | PDF

2.3 Speed comparison

In the conventional angular spectrum based method, when rotating from tilted plane to parallel plane, interpolation is always required. In our method, the interpolation operation is avoided by using NUFFT algorithm, which simplifies the calculation procedure. We compared the calculation time between our method and the angular spectrum based method. Figure 4 shows that the ratios of the calculation time between our method and the angular spectrum based method are in the range of 70%-90% when the test image size is varied from 500 × 500 to 1500 × 1500, which demonstrates that the calculation time of our method is less than the angular spectrum based method.

Fig. 4 Comparison of calculation time.

Download Full Size | PDF

3. Experiment and results

3.1 Numerical simulation

We first performed a numerical simulation experiment. A 512 × 512 image on the tilted plane is composed of four letters “A”, “B”, “C” and “D” and it is expanded to 1024 × 1024 by zero padding in order to avoid the aliasing. The hologram plane is at the distance of z = 200mm from the reference plane. The sampling rate of the hologram is ΔH = 8μm which is the same with the spatial light modulator (SLM) used in our optical experiment. The hologram is calculated from the image with different sampling rates ΔT’ = RΔT by using our proposed method. For each scale factor R, the hologram is calculated with three tilting angles 30°, 45° and 60°.

In order to show the scale effect, the numerical reconstructions on the tilted plane are all calculated by using the angular spectrum based method, that means all of the holograms are reconstructed on the plane of the same size (ΔT × 1024). The simulation results are presented in Fig. 5.The three columns in Fig. 5 represent the reconstructions of different angles. The images in first row are the reconstructions with the scale factor R = 1. The second and third rows are the reconstruction of the holograms which are calculated from the tilted plane with the scale factors R = 0.5 and R = 1.5 respectively. Due to the different sampling rates of the original images, the zoom-out and zoom-in effect can be seen in the reconstructed images, which proves the feasibility of our proposed scaled diffraction.

Fig. 5 Simulation results of the reconstruction on tilted plane from CGH.

Download Full Size | PDF

Next, the holograms of polygon-based object with different sampling rates are simulated. Here, a simple object composed by three polygons is simulated [Fig. 6(a)]. The front polygon is paralleled with the hologram plane and the other two polygons are tilted. Each polygon is calculated to the hologram plane based on the scaled diffraction described above with the same scale factor R, and the hologram is synthesized by adding all of the three distributions from the polygons. The distance from the front polygon to the hologram is 200mm. The complex hologram is encoded into the phase-only hologram (kinoform) by enforcing the amplitude of the complex distribution to unity. The holograms are calculated in this way with scale factor R = 0.6, R = 0.8, R = 1 and R = 1.6. All of the holograms are reconstructed on parallel planes by using the angular spectrum method. The results of reconstruction on the plane of z = 200mm are shown in Figs. 6(b)-6(e). The scale effect of zoom-in and zoom-out can be seen from Figs. 6(b), 6(c) and 6(e) compared with Fig. 6(d) which is the normal size (R = 1). The depth of field for a three-dimensional object is obvious when the scale factor is R = 1.6 in Fig. 6(e). Figure 6(f) is the reconstruction on the plane of 210mm in which both of the edge and letter “D” is clear while the letter “A” and the front polygon is blurred.

Fig. 6 Reconstructions of the polygon-based object. (a) The object (b) R = 0.6. (c) R = 0.8. (d) R = 1. (e) R = 1.6. (e) and (f) are reconstructions at front and back plane respectively.

Download Full Size | PDF

3.2 Optical experiment

Figure 7(a) shows the optical setup of the reconstruction. The 532nm laser is collimated by the lens and then illuminates the SLM. The reconstructed images are captured directly by the CMOS of the camera. The SLM we used is holoeye Pluto with 1920 × 1080 resolutions and 8μm pixel pitch. We calculated the holograms of the polygon-based object with different scale factors. We used the multi-random phase method to suppress the speckle noise caused by using laser. Figures 7(b)-7(e) are the reconstructions of the object with different scale factors R = 0.6, 0.8, 1, and 1.6 respectively. It is clearly that the size of the reconstruction in space is different with different scale factors, which means the display of polygon-based object with scaled size can be realized by using our proposed method. The depth of field of the three-dimensional object can be seen by comparing Figs. 7(e) and 7(f). The supplemented Media (Media 1) shows the dynamic process by moving the image plane from front to back.

Fig. 7 Optical reconstructions of polygon-based object. (a) Optical setup for holographic display. (b)-(e) Reconstructions of polygon-based object with different scale factors. (f) Reconstruction of object on the back focus plane (Media 1).

Download Full Size | PDF

4. Conclusion

We proposed a method of calculating the scaled diffraction between tilted planes. The diffraction from the tilted plane with variable sampling rates to the hologram is calculated by using the nonuniform fast Fourier transform. This scaled diffraction calculation includes two NUFFTs (or one NUFFT and one FFT) and one propagation operation, in which the conventional interpolation is not necessary. It overcomes the limitation of the sampling rates imposed by the FFT in conventional angular spectrum based method. The numerical and optical experiment exhibit the application of our method in calculating the hologram of three-dimensional object. This proposed method paves the way to the application of zoom-able holographic display of polygon-based objects.

Acknowledgments

This work was supported by the Major State Basic Research Development Program of China (2013CB328803 and 2013CB328804), a grant from the National High Technology Research and Development Program of China (2012AA03A302 and 2013AA013904), and the Aeronautical Science Foundation of China (No. 20125169021).

References and links

1. C. Slinger, C. Cameron, and M. Stanley, “Computer-generated holography as a generic display technology,” Computer 38(8), 46–53 (2005). [CrossRef]

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

3. K. Matsushima, “Computer-generated holograms for three-dimensional surface objects with shade and texture,” Appl. Opt. 44(22), 4607–4614 (2005). [CrossRef] [PubMed]

4. K. Matsushima and S. Nakahara, “Extremely high-definition full-parallax computer-generated hologram created by the polygon-based method,” Appl. Opt. 48(34), H54–H63 (2009). [CrossRef] [PubMed]

5. Y. Z. Liu, J. W. Dong, Y. Y. Pu, B. C. Chen, H. X. He, and H. Z. Wang, “High-speed full analytical holographic computations for true-life scenes,” Opt. Express 18(4), 3345–3351 (2010). [CrossRef] [PubMed]

6. L. Ahrenberg, P. Benzie, M. Magnor, and J. Watson, “Computer generated holograms from three dimensional meshes using an analytic light transport model,” Appl. Opt. 47(10), 1567–1574 (2008). [CrossRef] [PubMed]

7. Y. Pan, Y. Wang, J. Liu, X. Li, and J. Jia, “Fast polygon-based method for calculating computer-generated holograms in three-dimensional display,” Appl. Opt. 52(1), A290–A299 (2013). [CrossRef] [PubMed]

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

9. 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(9), 1755–1762 (2003). [CrossRef] [PubMed]

10. S. De Nicola, A. Finizio, G. Pierattini, P. Ferraro, and D. Alfieri, “Angular spectrum method with correction of anamorphism for numerical reconstruction of digital holograms on tilted planes,” Opt. Express 13(24), 9935–9940 (2005). [CrossRef] [PubMed]

11. A. Dutt and V. Rokhlin, “Fast Fourier transforms for nonequispaced data,” SIAM J. Sci. Comput. (USA) 14(6), 1368–1393 (1993). [CrossRef]

12. Q. H. Liu and N. Nguyen, “An accurate algorithm for nonuniform fast Fourier transforms (NUFFTs),” IEEE Microw. Guid. Wave Lett. 8(1), 18–20 (1998). [CrossRef]

13. L. Greengard and J. Y. Lee, “Accelerating the nonuniform fast Fourier transform,” SIAM Rev. 46(3), 443–454 (2004). [CrossRef]

14. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company Publishers, 2004).

Scaled diffraction calculation between tilted planes using nonuniform fast Fourier transform

Abstract

1. Introduction

2. Method

2.1 Calculate diffraction from tilted plane using NUFFT

2.2 Scaled diffraction calculation between tilted planes based on NUFFT

2.3 Speed comparison

3. Experiment and results

3.1 Numerical simulation

3.2 Optical experiment

4. Conclusion

Acknowledgments

References and links

Supplementary Material (1)

Cited By

Figures (7)

Equations (10)

Optics Express