## Abstract

Computer-Generated Holograms (CGHs) can be generated from three-dimensional objects composed of point light sources by overlapping zone plates. A zone plate is a grating that can focus an incident wave and it has circular symmetry shape. In this study, we propose a fast CGH generating algorithm using the circular symmetry of zone plates and computer graphics techniques. We evaluated the proposed method by numerical simulation.

© 2012 OSA

## 1. Introduction

Electro-holography is a promising three-dimensional (3D) display technology, using computer-generated holograms (CGHs), which is a method of digitally creating holograms [1]. For CGH generation, an object and reference lights are interfered by simulating light propagation on a computer. Then, the CGH is displayed on a spatial light modulator (SLM) such as a liquid crystal display (LCD) in order to reconstruct the 3D object recorded on the CGH.

There are three methods of creating CGH: point light source (PLS) based approach, polygon based approach and arbitrary shape surface based approach. In this study, we focus on the PLS approach, which is the simplest but most flexible approach because the PLS approach is suitable for parallel processing and hardware implementations such as field programmable gate arrays (FPGAs) and graphic processing units (GPUs).

Many methods of accelerating the CGH calculation have been proposed. The recurrence formulas [2–4] are an effective approach to improving the conventional method [1], which can substitute the square root in the conventional calculation for addition. As this approach is suitable for implementing on parallel computers or devices, some studies have already been reported using GPUs or FPGAs [5, 6].

Look-up-table (LUT) methods are also an effective approach from a different aspect [1, 7, 8]. As a CGH can be constructed by overlapping each wavefront of PLSs [9], this approach pre-calculates these wavefronts and stores them in a LUT. In order to generate a CGH using LUT methods, the wavefronts are read from the LUT and accumulated on the CGH. After calculating the LUT, this approach can calculate a CGH without any complex computation; however, a large amount of memory space for the LUT is required for storing the pre-calculated wavefronts. The N-LUT [8] method for reducing memory space has also been proposed; however, the amount of memory is still large when calculating a large 3D object or a high resolution CGH.

In this study, we propose a simple and efficient method of generating a CGH. A CGH generated from one PLS can be regarded as a zone plate. A zone plate is a grating that can focus an incident wave. As mentioned above, a CGH can be generated by overlapping zone plates. This method is capable of calculating CGH at high speed by generating and overlapping the zone plates quickly using the circular symmetry shape of the zone plate and computer graphics (CG) technique.

In Section 2, we explain the proposed method. In Section 3, we present the results of the proposed method and discuss the implications. Section 4 concludes this study.

## 2. Theory

#### 2.1. Conventional CGH calculation from PLS-based 3D object

A CGH is calculated by accumulating the wavefronts from all PLSs. The CGH calculation is expressed by the following formula [1]:

*x*,

_{α}*y*) denotes the coordinates of a CGH plane, (

_{α}*x*,

_{j}*y*,

_{j}*z*) denotes the coordinates of a PLS in three dimensions,

_{j}*x*=

_{α}_{j}*x*−

_{α}*x*,

_{j}*y*=

_{α}_{j}*y*−

_{α}*y*,

_{j}*i*is the imaginary unit,

*ϕ*is the argument of the object wave,

*k*is the wave number of reference light,

*N*is the total number of the PLSs,

*A*is the amplitude of

_{j}*j*-th PLS.

*I*represents the complex amplitude distribution on the CGH for

_{j}*j*-th PLS, and the phase is represented as ${\theta}_{j}\left({x}_{\alpha j},{y}_{\alpha j}\right)=k{\left\{{x}_{\alpha j}^{2}+{y}_{\alpha j}^{2}+{z}_{j}^{2}\right\}}^{1/2}$. arg{·} is the operator for taking the argument of the complex amplitude.

A CGH can be created by overlapping zone plates [9]. A zone plate has concentric ring shape because Eq. (1) is a function of the distance between the zone plate and PLS. In this study, we propose the fast algorithm for CGH using the property of the zone plate.

#### 2.2. Fast calculation of CGH using the circular symmetry of zone plate and LUT

An outline of our method is shown in Fig. 1. Utilizing the circular symmetry shape of the zone plate, it can be calculated by simple steps (STEP 1 and STEP 2). Subsequently, the CGH can be synthesized by STEP 3 in this figure.

The first step calculates a single line with complex amplitude of a zone plate along the radial direction. The line is calculated by the following recurrence formulas of Eqs. (2) and (3) [3]. We derived the recurrence formulas from the relation of the adjacent phases. The length of the line is defined by *R* = (*z _{j}*/

*p*) × tan(sin

^{−1}(

*λ*/2

*p*)) where

*p*is the pixel pitch of the device for displaying a CGH and

*λ*is the wave length of a incident wave. These recurrence formulas can calculate an adjacent phase on the CGH only in two additional operations faster than the direct calculation of

*θ*(

_{j}*x*,

_{α}_{j}*y*)

_{α}_{j}*γ*= 2

*π*

*/*

*λ*

*z*. As the start position of the formula is set as the center of the zone plate, the first terms of the formulas become

_{j}*θ*(0, 0) = 2

_{j}*π*

*z*

_{j}*/*

*λ*,

*δ*(0) =

_{j}*π*

*/*

*λ*

*z*, which is simpler than Ref. [3].

_{j}The second step rolls the line in order to form a complete zone plate, then we store it in a zone plate table, which is a temporary LUT for storing the zone plate. In this study, a CG technique for drawing a discrete circle is used in order to roll the line.

Fast drawing methods for a discrete circle have been developed. Bresenham’s algorithm [10] is a fundamental algorithm for drawing a discrete circle; however, it cannot fill all points inside the circle if we draw Bresenham’s circle by varying the integer radius [11]. Instead of Bresenham’s algorithm, we used the “arithmetical circle algorithm” [11], which solves this problem by slightly improving on Bresenham’s algorithm (Fig. 2). Similarly to Bresenham’s algorithm, this algorithm can draw a discrete circle by selecting the nearest pixels to be a circle that has integer radius. In this algorithm, the difference between the distance between a drawing pixel (*x _{α}*,

*y*) and the integer radius

_{α}*r*is defined as $\varepsilon ={x}_{\alpha}^{2}+{y}_{\alpha}^{2}-{r}^{2}$. Evaluating the

*ε*, the algorithm selects the nearest pixels to draw the circle. Due to

*ε*can be easily updated by simple addition or subtraction as shown in Fig. 2, it can draw a discrete circle rapidly. Therefore, we can obtain a zone plate from the line of the radius quickly, and it is stored in the zone plate table.

In the third step, one zone plate generated by STEP 1 and STEP 2 corresponding to PLSs with the same depth are translated according to (*x _{j}*,

*y*) and then overlapped on the CGH. It is unnecessary for the proposed method to retain the wavefronts for all PLSs in the table. The proposed method retains only the wavefronts for the same depth, so that the required memory for the zone plate table is small. For example, the required memory for the zone plate table is 90 MBytes when the parameters are

_{j}*λ*= 520 nm,

*z*= 0.3 m,

_{j}*p*= 8

*μ*m (in the parameters, the radius of the zone plate is 1,220 pixels). After overlapping zone plates corresponding to the same depth, we recalculate the zone plate table for the other depth. Note that the zone plates are treated as complex amplitude distribution when overlapping them.

Most recently, a similar method was proposed by Z. Yang et al. [12]; however, this approach needed to calculate the line more than eight times oversampling in order to create an accurate zone plate.

#### 2.3. Discretization error compensation

In this section, we describe the discretization error of the zone plate and the method for the error compensation. Figure 3 explains why the discretization error occurs. A certain phase *θ _{j}* on the CGH is calculated by

*θ*=

_{j}*k*(

*z*+

_{j}*pr*). For example, the phase of a coordinate (0,6) on the line as shown in Fig. 3 is represented as

*θ*=

_{o}*k*(

*z*+

_{j}*pr*)|

_{r}_{=6}. Considering the phase of coordinate (5,4) in Fig. 3, we express an actual phase

*θ*+

_{o}*θ*due to the distance from the center becomes

_{ε}*r*+

*ε*, where

*θ*=

_{ε}*kpε/z*. However, copying the phase

_{j}*θ*to the coordinate (5,4) by the “arithmetical circle algorithm” [11] is not an exact phase because the phase error between the actual and copied phases has

_{o}*θ*.

_{ε}The complex amplitudes *I*(*θ _{o}*) along the circle by the “arithmetical circle algorithm” is expressed by

*Re*[

*I*(

*θ*)] = cos(

_{o}*θ*),

_{o}*Im*[

*I*(

*θ*)] = sin(

_{o}*θ*). For the error compensation, applying the additional theorem of the trigonometric function and approximation cosΔ ≈ 1, sinΔ ≈ Δ (when Δ ≪ 1), the compensated complex amplitudes:

_{o}*ε*, which is calculated in Arithmetical circle algorithm and cos

*θ*and sin

_{o}*θ*have been already calculated.

_{o}## 3. Results and Discussion

In this section, we evaluate our method in terms of the image quality and computational speed. First, we examined the image quality. It is important to evaluate the focus property of a single zone plate because the CGH consists of many zone plates. We calculated three zone plates created by the direct calculation of Eq. (1), the proposed methods with or without the compensation.

In order to evaluate the focus property of the zone plates quantitatively, we analyze the zone plates by calculating Fresnel diffraction calculation [13]. Figure 4 shows the intensities of the diffracted light from the zone plates. Parameters are set as *λ* = 520 nm, *z _{j}* = 0.3 m,

*p*= 8

*μ*m, the CGH resolution is 1,920 × 1,080 pixels. We used the following environment: the operating system is Microsoft Windows 7 Professional Service pack 1, the CPU is AMD A8-3850 with 2.89 GHz (we used single core) and a memory of 4 GByte, the compilers are Microsoft Visual C++ 2010, the computation precision is set as double. In this figure, the light intensities are normalized at the maximum of the direct calculation of Eq. (1). Compared with the result by the direct calculation of Eq. (1), which is just focused at 0.3 m, the maximum light intensity of our proposed method without compensation decays by about 20% and the focal point is shifted by 400

*μ*m. Applying the error compensation mentioned in subsection 2.3, the maximum light intensity is improved by about 10%; however, the focal point is shifted by 488

*μ*m. Nevertheless, the problem of shifting focal point is not a serious problem because it can be solved easily by changing the depth of PLSs preliminarily.

Next, we estimate the reconstructed images of CGHs generated by the direct calculation of Eq. (1), the proposed method with the compensation, wavefront-recording plane (WRP) method [14–18] and WRP method combining the proposed method.

Our method can also be effective for the WRP method, which is an the effective way of calculating PLS-based CGH. The WRP method can calculate a CGH by two steps: the first step is the calculation of the wavefronts of PLSs on the virtual plane, which is referred to as WRP, which is placed between a 3D object and CGH. In the first step, the complex amplitude on a WRP is calculated by superimposing all of the wavefronts from all PLSs. We apply the proposed method to the first step. Here, the perpendicular distance between a 3D object and WRP is set as 5 cm. In the second step, we calculate the wave propagation from the WRP to CGH by fast Fourier transform (FFT) based diffraction calculation.

Figure 5 shows reconstructed images from each CGH: (a) is obtained by the direct calculation of Eq. (1), (b) is obtained by our method and (c) is obtained by the WRP method combining the proposed method with the compensation. The insets in Fig. 5 show magnified images of the legs. We measured the signal to noise ratio (SNR) between three images. In general, the human eye does not recognize the difference between reconstructed images when the SNR is 30 dB over. The SNR between (a) and (b) is 34.4 dB, (a) and (c) is 32.6 dB. Therefore, our method can obtain favorably reconstructed images.

Lastly, we describe the computational speed. Table 1 shows a comparison of calculation times between conventional methods and proposed methods of CGHs generated from a 3D object that consists of 11,646 PLSs. Our method can calculate the CGH 42.4 times faster than the direct calculation of Eq. (1), and 349.3 times faster than by the WRP method with the proposed method. In the WRP methods, the calculation times of the second step are 2.5 sec respectively. Therefore, we can reduce the calculation cost of the first step of the WRP method. Note that the calculation time of the proposed method with or without error compensation are almost the same because the additional computational cost of the compensation is tiny.

## 4. Conclusion

We proposed a fast algorithm for creating CGH using the properties of zone plate and CG technique. The proposed method can calculate CGH about 350 times faster than the direct calculation of Eq. (1) on a CPU. As our method is based on a simple CG technique, it will be useful implemented on GPUs or FPGAs. Moreover, it will contribute to developing 3D displays.

## Acknowledgments

The present research was supported in part by Support Center for Advanced Telecommunications Technology Research Foundation (SCAT).

## References and links

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

**2. **K. Matsushima and M. Takai, “Recurrence formulas for fast creation of synthetic three-dimensional holograms,” Appl. Opt. **39**, 6587–6594 (2000). [CrossRef]

**3. **T. Shimobaba and T. Ito, “An efficient computational method suitable for hardware of computer-generated hologram with phase computation by addition,” Comput. Phys. Commun. **138**, 44–52 (2001). [CrossRef]

**4. **H. Yoshikawa, “Fast Computation of Fresnel Holograms Employing Difference,” Opt. Rev. **8**, 331–335 (2001). [CrossRef]

**5. **Y. Ichihashi, H. Nakayama, T. Ito, N. Masuda, T. Shimobaba, A. Shiraki, and T. Sugie, “HORN-6 special-purpose clustered computing system for electroholography,” Opt. Express **17**, 13895–13902 (2009). [CrossRef] [PubMed]

**6. **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]

**7. **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**, D55–D62 (2008). [CrossRef] [PubMed]

**8. **S. -C. Kim, J. -M. Kim, and E. -S. Kim, “Effective memory reduction of the novel look-up table with one-dimensional sub-principle fringe pattern in computer-generated holograms,” Opt. Express **20**, 12021–12034 (2012). [CrossRef] [PubMed]

**9. **L. Mertz and N. O. Young, “Fresnel Transformation of Images,” in *Preceedings, International Conference on Optical Instruments and Techniques*, K. J. Habell, Ed. (Chapman & Hall, London), 305–310 (1961).

**10. **J. Bresenham, “A linear algorithm for incremental digital display of circular arcs,” Commun. ACM **20**, 100–106 (1977). [CrossRef]

**11. **E. Andres, “Discrete circles, rings and spheres,” Comput. Graphics **18**, 695–706 (1994). [CrossRef]

**12. **Z. Yang, Q. Fan, Y. Zhang, J. Liu, and J. Zhou, “A new method for producing computer generated holograms,” J. Opt. **14**, 095702 (2012). [CrossRef]

**13. **J. W. Goodman, *Introduction to Fourier Optics*Roberts & Company (2004).

**14. **T. Shimobaba, N. Masuda, and T. Ito, “Simple and fast calculation algorithm for computer-generated hologram with wavefront recording plane,” Opt. Lett. **20**, 3133–3135 (2009). [CrossRef]

**15. **T. Shimobaba, T. Ito, N. Masuda, Y. Ichihashi, and N. Takada, “Rapid calculation algorithm of Fresnel computer-generated-hologram using look-up table and wavefront-recording plane methods for three-dimensional display,” Opt. Express **18**, 19504–19509 (2010). [CrossRef] [PubMed]

**16. **P. Tsang, W. -K. Cheung, T. -C. Poon, and C. Zhou, “Holographic video at 40 frames per second for 4-million object points,” Opt. Express **19**, 15205–15211 (2011). [CrossRef] [PubMed]

**17. **J. Weng, T. Shimobaba, N Okada, H. Nakayma, M. Oikawa, N. Masuda, and T. Ito, “Generation of real-time large computer generated hologram using wavefront recording method,” Opt. Express **20**, 4018–4023 (2012). [CrossRef] [PubMed]

**18. **P. W. M. Tsang, K. W. K Cheung, and T.-C. Poon, “Real-time relighting of digital holograms based on wavefront recording plane method,” Opt. Express **20**, 5962–5967 (2012). [CrossRef] [PubMed]