## Abstract

Computer-Generated Holograms (CGHs) can be generated by superimposing zoneplates. A zoneplate is a grating that can concentrate an incident light into a point. Since a zoneplate has a circular symmetry, we reported an algorithm that rapidly generates a zoneplate by drawing concentric circles using computer graphic techniques. However, random memory access was required in the algorithm and resulted in degradation of the computational efficiency. In this study, we propose a fast CGH generation algorithm without random memory access using run-length encoding (RLE) based recurrence relation. As a result, we succeeded in improving the calculation time by 88 %, compared with that of the previous work.

© 2015 Optical Society of America

## 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]. However, there are many barriers to practical electro holography, one of which is the large amount of computational CGH calculations required.

Many methods of accelerating CGH calculations have been proposed [1–10]. The look-up table (LUT) is an effective approach to calculating CGHs [1]. As this method precalculates wavefronts of point light sources (PLSs), which are the minimum component of the PLS-based object model, CGHs can be calculated with fewer calculations. However, a large amount of memory space is required for storing the precalculated wavefronts.

Novel-LUT methods [4–6] have been proposed for reducing the required memory space; however, the amount of memory is still large when calculating a large 3D object or a high resolution CGH.

Our previous work also proposed a method for reducing the required memory space of LUT [7]. The method calculates a single line of wavefronts on a CGH plane and then copies the line circularly to create complete wavefronts. A CGH generated from the wavefront of one PLS can be regarded as a zoneplate, which can concentrate an incident light into a point. A CGH can be generated by superimposing zoneplates. As this method decompress a zoneplate from the single line utilizing its circular symmetry by computer graphics (CG) techniques, the required memory space for storing LUT is smaller than that of other works. However, the random memory access existing in the sequence of decompressing a zoneplate degrades the efficiency of CGH calculation.

In this study, we propose an effective LUT method without random memory access using a run-length encoding (RLE) approach. We developed an RLE-based recurrence relation that can obtain the run-length of the same values in the zoneplate. Using this recurrence relation, the same values can be written simultaneously, which improves the efficiency of CGH calculation.

The introduction of RLE approach to CGH calculation has already been reported in [5], which uses RLE for reducing the spatial redundancy of a 3D object itself, whereas our method is different because we use RLE for reducing the spatial redundancy of the zoneplate itself. As a result, we succeeded in improving memory access efficiency whilst keeping the small memory size of LUT.

In Section 2, we explain the proposed method. In Section 3, we present the results of the proposed method. In Section 4, we conclude 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 [8]. The CGH calculation is expressed by the following formula [1]:

*x*,

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

_{α}*x*) denotes the coordinates of a PLS in three dimensions,

_{j},y_{j},z_{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.

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

_{j}*j*-th PLS, and the phase is represented as ${\theta}_{j}({x}_{\alpha j},{y}_{\alpha j})=k{\{{x}_{\alpha j}^{2}+{y}_{\alpha j}^{2}+{z}_{j}^{2}\}}^{1/2}$. arg{·} is the operator for taking the argument of the complex amplitude. Strictly speaking, this equation is for kinoform, but the proposed method is also available for amplitude CGHs.

A CGH can be created by superimposing zoneplates [8], which have a concentric ring shape because Eq. (1) is a function of the distance between a certain pixel of a zoneplate and PLS. Our previous study utilizes the circular symmetry of a zoneplate to reduce the amount of memory required of LUT. In the study, a zoneplate can be calculated by two simple steps. Subsequently, the CGH can be synthesized by superimposing the zoneplates.

The first step calculates the complex amplitude of a single line of a zoneplate along the axial direction. The second step rolls the line to form a complete zoneplate using the CG technique to draw a discrete circle such as the Bresenham algorithm. The required memory size of LUT for storing the precalculated line is only proportional to the radius, which is smaller than that of other LUT methods. However, the efficiency of that calculation is degraded due to the random memory access in drawing the discrete circle.

There are two approaches similar to that of our previous work [9, 10]. S. Lee et al. proposed a similar approach [9]; however, their method requires precalculating the distribution of the distance from the center of a zoneplate in order to create a complete zoneplate, which increases the overhead of memory access. Z. Yang et al. also reported a similar approach [10]; however, their method requires the line to be calculated more than eight times oversampling in order to generate an accurate zoneplate.

#### 2.2. Zoneplate calculation using RLE approach

In order to improve the efficiency of CGH calculation, we developed a method for generating a zoneplate without random memory access using an RLE approach. Figure 1 shows the outline of our method.

The first step reads the precalculated line with a complex amplitude distribution from LUT. The line is shown at *y* = 0 in STEP 2 of Fig. 1. In STEP 1 of Fig. 1, the numbers in the boxes indicate indexes from the center of the zoneplate.

The second step copies the single line data to another area according to the sequence of the indexes in order to avoid the random memory access. For example, in calculating the line of *y* = 7, which is surrounded by a red frame in STEP 2, our method first copies the data indexed as 7 from the LUT three times, the data indexed as 8 from the LUT two times, and other copies are performed in the same way. The maximum index on each line is defined by (*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 an incident wave.

The properties of the sequence are as follows:

- The sequence can be grouped by the same index.
- The index of the adjacent group is incremented by 1.
- The index of the first group on each line corresponds to
*y*.

Therefore, we can obtain the sequence by calculating the length of each group. A complete zoneplate can be calculated by applying that process on each line repeatedly, as shown in STEP 3.

Figure 2 shows the method for calculating the length of groups. In Fig. 2, *n* is the group number, *Y* is the value of Axis *y* of a single line, *R _{n}* is the rounded value of

*r*, which is the real distance between the center of the zoneplate and a certain pixel.

_{n}*r*corresponds to the value in the group. The range of

_{n}*r*to be rounded to

_{n}*R*is [

_{n}*R*− 0.5 :

_{n}*R*+ 0.5).

_{n}*a*and

_{n}*b*are the intersection points between the line

_{n}*y*=

*Y*and circles whose radius are

*R*− 0.5 and

_{n}*R*+ 0.5. The same value of the complex amplitude distribution of a zoneplate runs between

_{n}*a*and

_{n}*b*.

_{n}The run-length of group *n* is defined as

*floor*(·) means round down and

*ceil*(·) means round out.

According to Fig. 2,

where,The relation of *a _{n}*,

*b*,

_{n}*d*is,

_{n}According to Eqs. (2), (6), (7) and (12), *l _{n}* is calculated by only three additions and one square-root calculation because

*b*and ${b}_{n}^{2}$ are equal to

_{n}*a*

_{n}_{+1}and ${a}_{n+1}^{2}$. In our method, both

*b*and ${b}_{n}^{2}$ are stored into registers temporarily after

_{n}*l*is obtained in order to use these values as

_{n}*a*

_{n}_{+1}and ${a}_{n+1}^{2}$. Therefore, the computational amount of generating the sequence is small.

#### 2.3. Discretization error compensation

The zoneplate calculated by our method is not exactly the same as that calculated by Eq. (1) because it includes discretization error [7].

Figure 3 describes why that error occurs. A certain phase *θ _{j}* on the CGH is calculated by

*θ*=

_{j}*k*(

*z*+

_{j}*pr*), where

*r*is the distance normalized by

*p*between the center of the zoneplate and a certain pixel on

*y*= 0. 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 the actual phase

*θ*+

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

_{ε}*r*+

*ε*, where

*θ*=

_{ε}*kpε*/

*z*,

_{j}*ε*=

*x*

^{2}+

*y*

^{2}

*− r*

^{2}. As our algorithm precaluclates the complex amplitude of radial direction and copies the values to another area, the discretization error

*θ*always exists.

_{ε}The complex amplitude *I*(*θ*_{0}) with the error is expressed as *I*(*θ*_{0} +*θ _{ε}*). Applying the additional theorem of the trigonometric function and approximation cosΔ ≈ 1, sinΔ ≈ Δ (when Δ ≪ 1), the compensated complex amplitudes are:

The discretization error *θ _{ε}* (

*x,r*) is obtained by the following recurrence relations:

Therefore, the additional cost of the error compensation is small.

## 3. Results and discussion

In this section, we evaluate our method in terms of image quality and computational speed. We examined the image quality by simulating reconstructed images of CGH using Fresnel diffraction calculation [11] at some distances around a desired focal position. For comparison, we also evaluate the CGH created by Eq. (1) in double precision and calculate the average signal-noise-ratio (SNR).

The object used for calculating CGH consisted of 11,646 PLSs located at *z _{j}* = 0.3 m and 0.05 m from CGH plane. The following parameters were set:

*λ*= 520 nm

*p*= 8

*μ*m, the resolution of CGH was 1,920 × 1,080 pixels. We used the following environment for the calculation: Microsoft Windows 7 Professional Service pack 1, Intel Core i7 4790K with 4 GHz (we used single core) and a memory of 8 GBytes, the compiler was Microsoft Visual C++ 2013 with/arch:AVX2, /fp:fast, /Ox, /Ob2, /Oi, /Ot compiler options.

Table 1 shows the results of our method in single and double precision and Fig. 4 shows reconstructed images obtained by a numerical simulation using Fresnel diffraction calculation [11] at 0.3 m and 0.05 m from GGHs: Fig. 4(a) and Fig. 4(d) are obtained by the direct calculation of Eq. (1), Fig. 4(b) and Fig. 4(e) are our previous method [7] and Fig. 4(c) and Fig. 4(f) are the proposed method with the compensation.

First, we describe the image quality. The SNR between the reconstructed image of direct calculation of Eq. (1) and that of proposed method are 30 dB over in each distances and each precisions. The human eye does not recognize the difference when the SNR is 30 dB over. Therefore, our method can obtain favorably reconstructed images.

Second, we describe the computational speed. Our method can obtain the CGH of the object about 30 times faster than the direct calculation of Eq. (1). Moreover, the computational time is reduced by 88 % compared to that of our previous method. Note that the calculation time of the proposed method with or without compensation is almost the same because the additional computational cost of the compensation is tiny. The speed-up ratio of 0.05 m is higher than that 0.3 m. In smaller distance, the radius of a zoneplate is short. Thus, the data size of storing the zoneplate in a buffer temporarily for overlapping the zoneplate is small enough to use memory band-width effectively. Therefore, the speed-up ratio of 0.05 m become higher than that of 0.3 m. On the other hand, our previous algorithm includes a random memory access in calculating a zoneplate, which causes a decreasing of the speed-up ratio in 0.05 m.

## 4. Conclusion

We proposed a fast CGH calculation method using an RLE based recurrence relation. Our method eliminates the need for random memory access in the sequence of creating zoneplates by the recurrence relation. As a result, we succeeded in reducing the calculation time by 88 % compared to that of our previous method without any degradation of image quality.

## Acknowledgments

This work is partially supported by JSPS KAKENHI Grant Numbers 25330125 and 25240015, and the Kayamori Foundation of Information Science Advancement and Yazaki Memorial Foundation for Science and Technology.

## References and links

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

**2. **N. Takada, T. Shimobaba, H. Nakayama, A. Shiraki, N. Okada, M. Oikawa, N. Masuda, and T. Ito, “Fast high-resolution computer-generated hologram computation using multiple graphics processing unit cluster system,” Appl. Opt. **51**, 7303–7307 (2012). [CrossRef] [PubMed]

**3. **H. Niwase, N. Takada, H. Araki, H. Nakayama, A. Sugiyama, T. Kakue, T. Shimobaba, and T. Ito, “Real-time spatiotemporal division multiplexing electroholography with a single graphics processing unit utilizing movie features,” Opt. Express **22**, pp.28052–28057 (2014). [CrossRef] [PubMed]

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

**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**, 1030–1041 (2009). [CrossRef] [PubMed]

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

**7. **T. Nishitsuji, T. Shimobaba, T. Kakue, N. Masuda, and T. Ito, “Fast calculation of computer-generated hologram using the circular symmetry of zone plates,” Opt. Express **20**, 27496–27502 (2012). [CrossRef] [PubMed]

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

**9. **S. Lee, H. C. Wey, D. K. Nam, D. S. Park, and C. Y. Kim, “Fast hologram pattern generation by radial symmetric interpolation,” Proc. SPIE 8498, Optics and Photonics for Information Processing VI, 84980O (2012).

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

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