Acceleration of the calculation speed of computer-generated holograms using the sparsity of the holographic fringe pattern for a 3 D object

In computer-generated hologram (CGH) calculations, a diffraction pattern needs to be calculated from all points of a 3-D object, which requires a heavy computational cost. In this paper, we propose a novel fast computer-generated hologram calculation method using sparse fast Fourier transform. The proposed method consists of two steps. First, the sparse dominant signals of CGHs are measured by calculating a wavefront on a virtual plane between the object and the CGH plane. Second, the wavefront on CGH plane is calculated by using the measured sparsity with sparse Fresnel diffraction. Experimental results proved that the proposed method is much faster than existing works while it preserving the visual quality. © 2016 Optical Society of America OCIS codes: (090.0090) Holography; (090.1760) Computer holography. References and links 1. J. Weng, T. Shimobaba, N. Okada, H. Nakayama, M. Oikawa, N. Masuda, and T. Ito, “Generation of real-time large computer generated hologram using wavefront recording method,” Opt. Express 20(4), 4018–4023 (2012). 2. C. Slinger, C. Cameron, and M. Stanley, “Computer-generated holography as a generic display technology,” Computer 38(8), 46–53 (2005). 3. M. Lucente, “Interactive computation of holograms using a look-up table,” J. Electron. Imaging 2(1), 28–34 (1993). 4. S.-C. Kim and E.-S. Kim, “Effective generation of digital holograms of three-dimensional objects using a novel look-up table method,” Appl. Opt. 47(19), D55–D62 (2008). 5. S.-C. Kim and E.-S. Kim, “Fast computation of hologram patterns of a 3D object using run-length encoding and novel look-up table methods,” Appl. Opt. 48(6), 1030–1041 (2009). 6. S.-C. Kim, J.-M. Kim, and E.-S. Kim, “Effective memory reduction of the novel look-up table with onedimensional sub-principle fringe patterns in computer-generated holograms,” Opt. Express 20(11), 12021–12034 (2012). 7. T. Nishitsuji, T. Shimobaba, T. Kakue, and T. Ito, “Fast calculation of computer-generated hologram using runlength encoding based recurrence relation,” Opt. Express 23(8), 9852–9857 (2015). 8. K. Matsushima and M. Takai, “Recurrence formulas for fast creation of synthetic three-dimensional holograms,” Appl. Opt. 39(35), 6587–6594 (2000). 9. H. Yoshikawa, S. Iwase, and T. Oneda, “Fast computation of Fresnel holograms employing difference,” Proc. SPIE 3956, 48–55 (2000). 10. T. Shimobaba, N. Masuda, and T. Ito, “Simple and fast calculation algorithm for computer-generated hologram with wavefront recording plane,” Opt. Lett. 34(20), 3133–3135 (2009). 11. A. Symeonidou, D. Blinder, A. Munteanu, and P. Schelkens, “Computer-generated holograms by multiple wavefront recording plane method with occlusion culling,” Opt. Express 23(17), 22149–22161 (2015). 12. D. Arai, T. Shimobaba, K. Murano, Y. Endo, R. Hirayama, D. Hiyama, T. Kakue, and T. Ito, “Acceleration of computer-generated holograms using tilted wavefront recording plane method,” Opt. Express 23(2), 1740–1747 (2015). 13. T. Tommasi and B. Bianco, “Computer-generated holograms of tilted planes by a spatial frequency approach,” J. Opt. Soc. Am. A 10(2), 299–305 (1993). 14. 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). 15. N. Okada, T. Shimobaba, Y. Ichihashi, R. Oi, K. Yamamoto, M. Oikawa, T. Kakue, N. Masuda, and T. Ito, “Band-limited double-step Fresnel diffraction and its application to computer-generated holograms,” Opt. Express 21(7), 9192–9197 (2013). Vol. 24, No. 22 | 31 Oct 2016 | OPTICS EXPRESS 25317


Introduction
Holographic 3-D is a promising technology in realistic 3-D display, which could provide fully satisfying 3-D perception without shutter/polarized 3D glasses in the 3-D stereoscopy [1].Computer-generated holograms (CGHs) are attracting increasing interest for the holographic display.CGHs are generated by numerically calculating a holographic fringe pattern on a computer simulation.CGHs have attractive features, such as not requiring specialized holographic recording materials and the physical manifestation of synthetic objects [2].
There are various CGH calculation methods for 3-D object, which are point light source (PLS) based methods [3][4][5][6][7][8][9][10][11][12], polygon based methods [13,14], RGB-D images based method [15], and so on.The PLS-based methods deal with 3-D object points as an aggregate of PLS.By using a ray tracing algorithm, they can generate CGHs more flexible than polygon-based method [10].However, in ray tracing process, a heavy computation is required to calculate the complex amplitude of all 3-D object points at every pixel on the CGH plane.As a result, the computational complexity is known as one of the challenging problems [10].
To accelerate the computational speed of CGH calculation, several methods were proposed.In [3], a look-up table (LUT) based approach was proposed.It could reduce the computational cost by pre-calculating and storing holographic fringe patterns.Kim et al. proposed novel LUT based approaches to reduce the required memory space for storing the pre-calculated values [4][5][6].However, a large amount of memory and additional operations (e.g., file read or array access) are still required for high resolution CGHs [7].In [8] and [9], recurrence based approaches were proposed using recurrence relations.However, approximation errors could be accumulated and propagated.Recently, wavefront recording plane (WRP) based approaches were proposed [10][11][12].In [10], a computational complexity was reduced by calculating the wavefront (i.e., complex amplitude) on the WRP followed by Fresnel diffraction calculation with fast Fourier transform (FFT).In [11] and [12], multi-WRP and tilted WRP methods were proposed.
As aforementioned, a holographic fringe pattern can be represented by FFT based light propagation such as Fresnel diffraction or angular spectrum.Most of coefficients of the Fourier representation are very small or zero in image, video, and audio signals [16,17].Without the loss of generality, the holographic fringe patterns are sparse in the Fourier domain as well [18].Based on the sparsity property of holographic fringe pattern, compressive holography methods were proposed that reconstruct original data from a small set of sampled signals [18,19].To reduce redundant calculations, the sparsity property can play key role in the CGH calculation as well as the compressive holography.The existing CGH calculation methods do not consider the sparsity property of CGHs yet.
In this paper, we propose a novel fast CGH calculation method using the sparsity of holographic fringe patterns.The main contributions of this paper are the following.2) The sparsity property of holographic fringe pattern is calculated in the wavefront on a virtual plane.In this paper, the wavefront on the virtual plane is effectively obtained by clustering 3-D object points.Multiple wavefronts of corresponding clustered points are calculated in parallel.Using the wavefront on the virtual plane, the sparsity of the holographic fringe pattern is measured.With the measured sparsity, sparse fast Fourier transform (sFFT) is applied to calculate the wavefront on the CGH plane using Fresnel diffraction.The computational complexity of sFFT is about O(k log k), which is lower than that of FFT, O(N log N).Note that N indicates the entire signals (e.g., N-dimensional entire CGH pixels) and k indicates the number of sparse dominant signals.
The remainder of this paper is organized as follows.In Section 2, the sparsity of the CGHs is investigated.In Section 3, the proposed fast CGH calculation method is described.Section 4 presents the experiments and results to evaluate the performance of the proposed method.Finally, conclusions are drawn in Section 5.

Sparsity of CGH and the visual quality
In this section, we observe the visual quality of CGH signals.CGH signals are sampled with two different sampling scheme.The visual qualities of corresponding numerical reconstructions are observed.The details of the experimental parameters (e.g., wavelength, sampling pitch, distance between object and CGH plane, CGH resolution, etc.) are the same as the experimental conditions described in Table 2 in Section 4. Figure 1 shows examples of numerically reconstructed results by sampled CGH signals, where 10% of the signals have been used.Figure 1(a) and 1(b) show a 3-D object point cloud and its holographic fringe pattern (i.e., CGH) generated by ray tracing, respectively.Figure 1(c) shows the result reconstructed by randomly sampled signals.Figure 1(d) shows the result reconstructed from the top 10% magnitude of CGH signals.As shown in the Fig. 1(c) and 1(d), sampled sparse signals could provide visually acceptable results.In particular, magnitude based sampling could provide better visual quality than random sampling with the same sampling ratio.This observation is consistent with the previous works [18].As a result, we observe that top magnitude samples of the CGHs signals have an important effect on the visual quality.In addition, we measure the visual quality depending on the amount of samples in CGHs. Figure 2 shows average peak signal to ratio (PSNR) [dB] according to sampling ratio for four data sets (please refer to Section 4 for details).The PSNR is calculated between numerical reconstructed results (i.e., numerically reconstructed 2-D image as shown in Fig. 1(c) and 1(d)) from the entire CGH signals and sampled CGH signals.The numerical reconstructed results from the entire CGH signals are considered as a reference to measure the PSNR.To obtain the numerical reconstructed results from CGH of 3-D data, the direct calculation using ray tracing is used in our experiments [20].As shown in the Fig. 2, in the magnitude based sampling case, 30 dB PSNR can be achieved with around 5% dominant signals (i.e., top 5% magnitude of CGH signals).PSNR over 30 dB provides viewers with enough visual image quality [21].
We observe that the CGH signals are sparse enough so that about top 5 -10% dominant signals in the CGHs could provide feasible visual quality of the numerical reconstruction.Based on this observation, we propose a fast CGH calculation considering the sparsity property of holographic fringe pattern in the following section.

Proposed fast CGH calculation method
Figure 3 shows the overview of the proposed CGH calculation for 3-D object based on the sparsity of holographic fringe pattern.The proposed method consists of two steps.First, the proposed method calculates the wavefronts of object points clustered on a virtual plane which is located between the object and the CGH plane.The sparsity of the holographic fringe pattern is measured by the number of large coefficients (red dots in Fig. 3) on the virtual plane.Second, the CGH is calculated by propagating the wavefront on the virtual plane to the CGH plane using sparse Fresnel diffraction with the measured sparsity.The details of the proposed fast CGH calculation are illustrated in next subsections.

Multiple ray tracing for calculating wavefront on the virtual plane
Sparse Fresnel diffraction using sFFT for calculating wavefront on the CGH plane Fig. 3. Overview of the proposed fast CGH calculation for 3-D objects.In the first step, the wavefront on the virtual plane is calculated using multiple ray tracing for each object point cluster.In the second step, the wavefront on the CGH plane is calculated using sparse Fresnel diffraction with sFFT.Note that the red dots indicate a small number of dominant signals on the virtual plane.For sFFT, the number of dominant signals (i.e, sparsity) is measured based on the magnitude.z 1 represents the distance between the object and the virtual plane.z 2 represents the distance between the virtual plane and the CGH plane.

Calculation of the wavefront on the virtual plane using multiple ray tracing
In this section, we present the fast calculation of the complex amplitude (i.e., wavefront) on the virtual plane.Since the virtual plane is placed close to the object, point lights traverse limited areas [10].As a result, the computational complexity of each point light is reduced by employing the virtual plane.However, a heavy computational cost is still required to calculate the wavefronts of a large number of object points.In this paper, the computational time of the wavefront is reduced by simultaneously calculating multiple complex amplitudes on the virtual plane.To that end, 3-D object points are divided into S clusters.Let C t denote the t-th cluster (t = 1, …, S) and N t denote the number of points belonging to the t-th cluster.Let u VP denote the wavefront on the virtual plane.The wavefront on the virtual plane, u VP , is calculated by integrating multiple complex amplitudes of clustered points.It can be written as ), in the t-th cluster (C t ).A t i indicates the intensity of i-th object point in the t-th cluster (C t ).k represents the wave number, k = λ/2π.λ indicates the wave length of the reference light.
Figure 4 shows the proposed fast calculation of the wavefront on the virtual plane using multiple ray tracing.As shown in Fig. 4, the number of object points within a cluster is much smaller than that of all 3-D object points.The areas traversed by point lights are small on the virtual plane.By calculating the wavefronts of point clusters in parallel, the computational speed of the ray tracing can be significantly improved for 3-D object.In this paper, the multiple wavefronts on the virtual plane are calculated by CPU multi threads in openMP.The ray tracing can be parallelized with unclustered 3-D points.However, a 3-D object data set has thousands of points.If a large number of CPU threads are available, additional costs such as memory allocation may significantly increase for the point-grained parallel computing.On the other hand, the cluster-grained parallel computing can make effective use of the limited number of CPU threads in practical (the number of clusters can be much less than the number of 3-D points).By assigning each cluster to each CPU thread for cluster-based parallel computing, the calculation times of wavefronts of all 3-D points can be effectively improved.The 3-D object points are simply divided into clusters according to the order of point number.Fig. 4. Proposed fast calculation of wavefront on the virtual plane using multiple ray tracing.The red dots and the green dots represent the object points belonging to the first and second clusters, C 1 and C 2 , respectively.The blue dots represent the object points belonging to the S-th cluster, C S .Notably, multiple wavefronts of each cluster on the virtual plane are calculated using ray tracing in parallel.Note that the radius of small area traversed by i-th point light in the C t is defined as W t i = |z t i |tanθ = |z t i |tan(sin −1 (λ/2p)), as reported in [10].z t i represents the distance between i-th point in the C t and the virtual plane.p is the sampling pitch, which is 8.5 μm in this paper.
In this paper, the virtual plane plays a key role in obtaining the signal characteristics of the holographic fringe pattern as well as to reduce the computational cost of ray tracing.The wavefront on the virtual plane is followed by Fresnel diffraction with sFFT for CGH calculation.The details of the proposed Fresnel diffraction with sFFT are described in the next session.

Calculation of the wavefront on CGH plane using sparse Fresnel diffraction with sFFT
To calculate the wavefront on the CGH plane from the virtual plane, Fourier transform based propagation such as Fresnel diffraction has been widely used [3][4][5][6][7][8][9][10][11][12].Fourier transform based propagation reduces computational time of CGH calculation by using FFT.However, it is not enough to reduce the computational complexity for calculating CGH with a high resolution.In this paper, we propose a fast CGH calculation using sFFT in Fourier based propagation.
Recently, sFFT methods have been proposed for reducing the computational complexity of FFT [16,17].The main idea of sFFT is to sample a small number of signals (with the sparsity of signals).By calculating the FFT with sparse signals, sFFT could achieve a low computational complexity.However, it requires the locations and values of dominant spare signals in the frequency domain to avoid data loss or incorrect FFT calculation.By performing filtering or a permutation scheme, sFFT increases the probability of capturing dominant signals in Fourier domain and outperforms the FFT [22].
As described in Section 2, sparse signals of CGHs could provide feasible visual quality of numerical reconstruction.A sparse Fresnel diffraction with sFFT can accelerate the calculation of the wavefront on the CGH plane.The main idea of the proposed method is to measure the sparsity of the holographic fringe pattern from the wavefront recorded on the virtual plane.With the measured sparsity, the CGH is rapidly generated using sFFT.The signal characteristics of CGH are similar to those of the wavefront on virtual plane since the diffraction calculation on the virtual plane is equivalent to CGH calculation from the 3-D object [10].Figure 5 shows the proposed fast CGH calculation using sparse Fresnel diffraction with sFFT.The sparsity of the holographic fringe pattern is measured by finding the number of dominant signals on the virtual plane [23].In this paper, dominant signals are found as top magnitude signals, i.e., whose magnitude is greater than a threshold.The wavefront on the virtual plane and its sparsity (the number of dominant signals) are inputs of sFFT.In the sFFT, input signals are randomly permuted and sampled referring to the sparsity of wavefront on virtual plane.Then, the FFT is performed with the permuted and sampled input signals.
Finally, the wavefront on CGH plane is obtained by sparse Fresnel diffraction using sFFT.Let u(ξ,η) denote the wavefront on the CGH plane.It can be written as where  [•] and  −1 [•] represent the sFFT and inverse sFFT operators.(ξ,η) is the pixel position of the CGH in real domain.z 2 represents the perpendicular distance between the virtual plane and the CGH plane.h(ξ,η) is the impulse response function, which is equal to

Data sets
To evaluate the performance of the proposed fast CGH calculation for 3-D object, four publicly available 3-D point cloud data sets were utilized.Among them, three data sets were collected from the Berkeley instance recognition data set (BigBIRD) [24]: Baby toy, Soft soap, and Syrup.In addition, Bunny was collected from the Stanford 3D scanning repository [25].Table 1 shows the data sets used in our experiment.For example, the number of object points for bunny is 35,947 and it is a synthetic object.

Experimental setup
For our computing environment to calculate the CGHs from the 3-D point cloud data sets, we used Microsoft Windows 7 Professional Service pack 1, Intel Core i7-4770 CPU @ 3.40 GHz and a 32 GBytes memory, Microsoft Visual Studio 2013 and openMP.Table 2 illustrates the CGH calculation conditions.We used the FFTW 3.2.1 library [26] for the FFT operation and the sFFT algorithm of [17].

Performance evaluation results for visual quality and computational speed
Experiments were performed in terms of visual quality and computational time in order to evaluate the performance of the proposed method.To demonstrate the visual quality, we compared the numerically reconstructed results from the CGHs generated by five CGH calculation methods, which were a ray tracing method, an LUT based method [3], a recurrence based method [8] (with stride 2 pixels), a WRP based method [10], and the proposed method.Notably, the results by the ray tracing method were considered as ground truth [7].In addition, we compared the PSNR between the numerical reconstructed results of ray tracing and other methods in order to quantitatively evaluate the performance of the visual quality of the CGHs.The PSNR value is increased as the difference between the CGH calculation method and the ground truth is decreased.Table 3 shows the computational time of the CGH calculation for each data set.We measured the average calculation times when a 3-D object is rotated 360 degrees.As seen in Table 3, the proposed method shows a low computational complexity for the CGH of 3-D objects.The computational times of the ray tracing method (i.e., ground truth) and the LUT based method were rapidly increased with the large number of object points.The recurrence based method was faster than the ray tracing.But it still required a lot of computational time.The WRP based method was faster than other three methods.The proposed method was at least 15 times faster than the WRP based method while providing visually plausible results (please see Fig. 7 and Table 4).For example, the computational times of calculating the wavefront on the virtual plane and the CGH plane are 0.08 s and 0.07 s, respectively, for Baby toy in the proposed method.
Figure 7 shows the reconstructed images obtained by a numerical reconstruction from each CGH calculation method.In the Fig. 7, the first row shows the numerical reconstructed results from the CGHs generated by ray tracing algorithm (i.e., ground truth).The second and third rows show the numerical reconstructed results from the CGHs generated by the LUT based and recurrence based methods, respectively.The fourth and fifth rows show the results from the CGHs calculated by the WRP based method and the proposed method, respectively.As shown in Fig. 7(b), the visual results by the LUT based method were very similar to the ground truth because pre-calculated values were just loaded and used.In Fig. 7(c), the visual result of the recurrence based method for Bunny had some distortions while visual results for other data sets provided visually plausible results.While a small number of object points are sparsely distributed in the Baby toy, Soft soap, and Syrup, a large number of points (35,947) stand close together in the Bunny.Therefore, in the Bunny case, it is easy to be affected by the approximation errors by the stride (2 pixels = 17 μm).As shown in Fig. 7(e), the proposed method provided visually plausible results for four data sets regardless of the number of points, similar to Fig. 7(d).Furthermore, to handle the specular objects (specular surfaces), polygon-based methods can be used [29,30] because they can deal with the 3-D object (3-D surface) with both diffuse and specular reflectance on each individual polygon.On the other hand, the proposed method is a point-based approach in order to calculate CGH of 3-D point cloud data.To handle the specular object using the proposed method, the conversion of specular object to point cloud can be one of alternatives.For example, the specular object data is converted into the point cloud data.Then, the proposed method can be applied to calculate CGH of the point cloud data.The specular object can be reconstructed by a rendering algorithm for specular surfaces [31].

Conclusions
This paper presented a fast CGH calculation method for 3-D objects.In the proposed CGH calculation framework, the sparsity property of holographic fringe pattern was considered.To measure the signal characteristics of CGHs, wavefront on a virtual plane was calculated.By dividing all object points into sub clusters and calculating the multiple wavefronts of the clusters in parallel, we improved the computational speed.In particular, based on the wavefront recorded on the virtual plane, we introduced a novel fast CGH calculation using sparse Fresnel diffraction with sFFT.In our experiment, we achieved at least 15 times faster than other existing works while preserving the visual quality.The numerical reconstructed results showed that the proposed CGH calculation method could achieve good visual quality.In addition, we verified the performance of the proposed method using an objective quality metric, PSNR.

1 )
We observe the visual quality of the numerical reconstruction and a small (sparse) set of CGH signals.The experiments show that the CGHs of 3-D object are signals which are sparse enough so that a 3-D object can be reconstructed with a small number of dominant CGH signals.Based on this observation, we propose a fast CGH calculation using the sparsity of the CGHs.

Fig. 1 .
Fig. 1.Examples of numerical reconstructed images of sampled CGH signals, where 10% of the signals were used.(a) 3-D object point cloud for Bunny.(b) Holographic fringe pattern (i.e., CGH) (c) The numerical reconstructed result from randomly selected 10% of CGH signals.(d) The numerical reconstructed result from the top 10% magnitude of CGH signals.

Fig. 2 .
Fig. 2. The quality of the numerical reconstruction according to the sampling ratio.The visual quality is measured by PSNR.

Fig. 5 .
Fig. 5. Schematic diagram of the proposed fast CGH calculation using sparse Fresnel diffraction with sFFT.The wavefront on the CGH plane is calculated from the wavefront on a virtual plane using sFFT.To that end, the sparsity of the wavefront is measured by counting the dominant signals based on the magnitude.
In the proposed method, the selection process is performed in Fourier domain of CGHs to investigate the magnitude threshold for about top 10% of the CGH fringe patterns of 3-D point cloud data in advance.In our experiment, the threshold value (Th) is 0.35 with which about top 10% of signals of CGH of 3-D point cloud data sets can be extracted.By applying the threshold to the magnitudes of wavefront on the virtual plane, the sparsity of the CGH can be estimated.The sparsity (k) estimated by the threshold (Th = 0.35) depends on the signal distributions of each 3-D point cloud data.With the wavefront on the virtual plane and estimated sparsity value (k) of CGH signals, k dominant CGH signals can be generated by sparse Fresnel diffraction.The threshold parameter and the sparsity of each 3-D point cloud data set are described in Section 4.2. Figure 6 shows an example of the distributions of the CGH fringe pattern for Bunny.As shown in Fig. 6(a), the CGH fringe pattern in Fourier domain has a few dominant signals and a lot of small or zero signals.Figure 6(b) shows the sparse distribution of the CGH fringe pattern after the selection process.

Fig. 6 .
Fig. 6.Examples of the distributions of the CGH fringe pattern for Bunny.(a) Original signal distribution of the CGH fringe pattern in Fourier domain.(b) Sparse distribution of the CGH fringe pattern in Fourier domain after the selection process.

Table 2 .
CGH calculation conditions in our experiments.Parameters Values Wave length (λ) 633 [nm] Sampling pitch (p) 8.5 [μm] Distance between the object and the virtual plane (z 1 ) 0.0500 [m] Distance between the virtual plane and CGH plane (z 2 ) 0.1844 [m] Distance between the object and the CGH plane (z = z 1 + z 2 ) 0.2344 [m] Resolution of CGH 1024x1024 [pixels] The number of object point clusters (S

Fig. 7 .
Fig. 7. Visual results of the numerical reconstruction from the CGHs generated by four existing methods and the proposed method for each data set.(a) Results of the ray tracing, (b) Results of the LUT based method [3], (c) Results of the recurrence based method [8], (d) Results of the WRP based method [10], (e) Results of the proposed method.

Fig. 8 .
Fig. 8.The quality and the computation time (step 2) as a function of k for Bunny.