## Abstract

A rapid calculation method of Fresnel computer-generated-hologram (CGH) using look-up table and wavefront-recording plane (WRP) methods toward three-dimensional (3D) display is presented. The method consists of two steps: the first step is the calculation of a WRP that is placed between a 3D object and a CGH. In the second step, we obtain an amplitude-type or phase-type CGH to execute diffraction calculation from the WRP to the CGH. The first step of the previous WRP method was difficult to calculate in real-time due to the calculation cost. In this paper, in order to obtain greater acceleration, we apply a look-up table method to the first step. In addition, we use a graphics processing unit in the second step. The total computational complexity is dramatically reduced in comparison with conventional CGH calculations. We show optical reconstructions from a 2,048 × 2,048 phase-type CGH generated by about 3 × 10^{4} object points over 10 frames per second.

© 2010 OSA

## 1. Introduction

Computer-generated-holograms (CGH) are widely used in the field of optics, especially as a three-dimensional (3D) display technique [1, 2]. However, the computational time involved in CGH obstructs the realization of a practical 3D display system using the CGH. Many methods have been proposed for accelerating the computational time toward real-time calculation. Current CGH calculations are mainly categorized into two approaches: polygon-based and point cloud approaches. In this paper, we focus on the acceleration of the point cloud approach, which treats a 3D object as a composite of multiple point light sources, and generates a CGH from a 3D object more flexibly than polygon-based approaches.

Look-up table (LUT) methods [3, 4, 5] pre-calculate light intensity distributions on a CGH of point light sources, then, we reference the LUT during the CGH calculation. An image hologram approach carried out real-time 3D color reconstruction because the approach placing a 3D object close to a CGH calculates a small region on the CGH in place of calculating the whole region on the CGH [6]. In addition, the method accomplishes greater acceleration using a LUT. LUT methods [3, 4, 5] calculate a CGH faster than without using LUT methods; unfortunately, they required large memory for the light intensity distributions, and cannot decrease the computational complexity, which is proportional to the total number of point light sources and the total number of CGH pixels. The image hologram approach reduces the computational complexity and memory amount for the light intensity distributions; unfortunately, the method cannot reconstruct a 3D object in the Fresnel region.

In the paper, a rapid calculation method of Fresnel CGH using LUT and the wave-recording plane (WRP) methods is presented. The WRP method reduces the computational complexity of CGH [7]. The method consists of two steps: the first step is the calculation of a virtual plane, which is referred to as WRP, which is placed between a 3D object and a CGH. In the second step, we obtain an amplitude-type or phase-type CGH to execute diffraction calculation from the WRP to the CGH. The first step of the previous WRP method was difficult to calculate in real-time due to the calculation cost. In this paper, in order to obtain greater acceleration, we apply a look-up table (LUT) method to the first step. In addition, we use a graphics processing unit (GPU) in the second step. The total computational complexity is dramatically reduced in comparison with conventional CGH calculations. We show optical reconstructions from a 2,048 × 2,048 phase-type CGH generated by about 3 × 10^{4} object points in real-time.

## 2. Improved wave-recording plane method by LUT method

The WRP method inspired by Ref.[6] reduces the computational complexity of CGH by two-step calculation. Figure 1(a) shows an outline of the algorithm.

In the first step, we record the complex amplitude on the WRP of a 3D object using a raytracing between the 3D object and the WRP. We record the complex amplitude *u _{w}* of all points of a 3D object using the following equation [3]:

where,
${R}_{\mathrm{wj}}=\sqrt{{({x}_{w}-{x}_{j})}^{2}+{({y}_{w}-{y}_{j})}^{2}+{d}_{j}^{2}}\approx {d}_{j}+\frac{({({x}_{w}-{x}_{j})}^{2}+{({y}_{w}-{y}_{j})}^{2})}{\left(2{d}_{j}\right)}$
is the distance between a 3D object point and a coordinate (*x _{w},y_{w}*) on the WRP, (

*x*) and

_{j},y_{j},z_{j}*A*are the coordinates whose origin is on the CGH and intensity of

_{j}*j*-th object point,

*λ*is the wavelength of the reference light, and

*N*is the total number of 3D object points,

*d*=

_{j}*z*−

_{j}*z*is the perpendicular distance between the WRP and a 3D object point, and $i=\sqrt{-1}$ . When placing the WRP close to the 3D object, the object light traverses small regions (the radius is

*W*) on the WRP. Therefore, the computational amount of Eq.(1) is very small.

_{j}The second step calculates the complex amplitude *u*(*x*,*y*) on the CGH using the diffraction calculation between the WRP and CGH. Since the WRP has amplitude and phase information for a 3D object light, the diffraction calculation is equivalent to the case of directly calculating the complex amplitude on the CGH from the 3D object. We use the Fresnel diffraction:

$$\phantom{\rule{.9em}{0ex}}\phantom{\rule{.9em}{0ex}}\phantom{\rule{.7em}{0ex}}\phantom{\rule{.3em}{0ex}}\phantom{\rule{.2em}{0ex}}=\frac{\mathrm{exp}\left(i\frac{2\pi}{\lambda}z\right)}{i\lambda z}{\mathcal{F}}^{-1}\left[\mathcal{F}\left[{u}_{w}(x,y)\right]\xb7\mathcal{F}\left[h(x,y)\right]\right]$$

where, *𝓕*[·] and *𝓕*
^{−1}[·] are the Fourier and inverse Fourier operators, *z* is the perpendicular distance between the WRP and the CGH, and the impulse response
$h(x,y)=\mathrm{exp}\left(i\frac{\pi}{\lambda z}({x}^{2}+{y}^{2})\right).$
.

We analytically obtain the Fourier transform of *h*(*x,y*), so that the numerical calculation of Eq.(2) is accelerated by using two fast Fourier transform (FFT) operations. In Eq.(2), *𝓕*[*u _{w}*(

*x,y*)] is accelerated by the FFT and

*𝓕*

^{−1}[·] is accelerated by the inverse FFT.

Lastly, we calculate the CGH. If we calculate an amplitude-type CGH, we calculate the light intensity *I*(*x,y*) = ∣*u*(*x,y*)+*R*(*x,y*)^{∣2}, where, *R*(*x,y*) is a reference light. On the other hand, if we calculate a phase-type CGH, we calculate θ(*x,y*)=*arg*{*u*(*x,y*)}, where, the operator *arg*{·} calculates the argument of *u*(*x,y*).

#### 2.1. Improvement of the WRP method

As shown in Fig.1(a), the maximum diffraction angle *θ _{o}* for reconstructing a 3D object light from the CGH is the following equation: sin

*θ*=

_{o}*λ*/(2

*p*), where,

*p*is the sampling pitch on the CGH. We need to set the spread angle for a 3D object point under the angle

*θ*. The radius

_{o}*W*on the WRP is the following equation:

_{j}where, the approximation is good when *θ _{o}* is small.

The first step of the previous WRP method calculated Eq.(1) on circular regions decided by the radius *W _{j}*. During the calculation, we need to judge either the circular regions or other regions. In this paper, in order to accelerate the judgment, we adopt the rectangular region inscribed in the circular region as shown in Fig.1(b). The length of a side on the rectangular region,

*L*is,

_{j}It is possible to omit a condition branch for the judgment of the circular regions in the programming. We could not confirm the deterioration on optical or numerical reconstructions using the rectangular region instead of the circular one.

In addition, we apply a LUT method to the first step. The previous WRP method directly calculated Eq.(1); however, the calculation time was not accomplished as a real-time calculation. In order to apply the first step to the LUT, we store the pre-calculated distributions of the complex amplitude Eq.(1) corresponding to *d _{j}* into the LUT. During the calculation of the first step, we accumulate the distributions of the complex amplitude by referencing the LUT using the parameter,

*x*and,

_{j},y_{j}*d*

_{j}The amount of memory of the LUT, which is proportional to *d _{j}*, is decided by the following equation:

where, the operators *max*{·} and *min*{·} take a maximum and minimum value in the argument. In Fig.1(c), the areas of rectangular regions increase proportional to *d _{j}*, thus, the LUT has a pyramid-like structure. For example, in case of a 3D object with the depth range 0.5

*mm*≤

*d*≤10.5

_{j}*mm*, the sampling pitch along the depth 0.5

*mm*,

*p*= 8

*µm*and

*λ*= 532

*nm*, the amount of memory is only about 220 Kbytes.

We consider the computational amount of the first step. Because *L _{j}* depends on

*d*, let us treat the average of

_{j}*L*for all points of a 3D object: $\overline{L}=\frac{1}{N}{\Sigma}_{j}^{N}{L}_{j}$ . The average of the rectangular area on the WRP is expressed as

_{j}*L̄*

^{2}. Therefore, the computational amount in the first step is 2

*αNL̄*

^{2}, where, “2” means the calculations of the real and imaginary parts, and

*α*is the number of arithmetic operations (e.g. additions and multiplications) in the summation of Eq.(1).

Next, the computational amount in the second step, which includes two FFT operations, is about 2*βN*
^{2}
_{x}log*N _{x}*, where the CGH size is

*N*×

_{x}*N*,

_{x}*β*is the number of arithmetic operations (e.g. additions and multiplications) in the FFT. Note that this computational amount does not depend on the number of 3D object points.

Therefore, the total computational amount is about 2*αNL̄*
^{2} + 2*βN*
^{2}
_{x}log*N _{x}*. And, if

*N*is very large, the total computational amount is expressed as 2

*αNL̄*

^{2}+ 2

*βN*

^{2}

_{x}log

*N*≈ 2

_{x}*αNL̄*

^{2}. When placing the WRP close to the 3D object,

*L̄*is much smaller than

*N*, therefore, we can see that the computational amount of the proposed method is much smaller than that of conventional ray-tracing algorithms [3, 4], which is

_{x}*αNN*

^{2}

_{x}.

In addition, in order to obtain more acceleration, we use a multi-thread programming by OpenMP [8] in the first step, and GPU computing in the second step. In recent years, several research have been reported on GPU-based CGH calculation [9, 10, 11, 12, 13, 14]. In this paper, for the Fresnel diffraction, we used our GPU-basedWave Optics (GWO) library [15, 16]. The library is a numerical calculation library for the diffraction calculations using a GPU. The library and a sample code can be downloaded from Ref.[17].

As shown in Eq.(2), the Fresnel diffraction can accelerate the computational time using FFT and inverse FFT; however, recent central processing units (CPUs) do not have sufficient computational power for real-time calculation. Therefore, we use GPU instead of CPU.

## 3. Results

We compare the calculation times when using a conventional algorithm without LUT [3] implemented on a CPU and GPU, the previous WRP method [7] implemented on a CPU, and, the improved WRP method.

We used Microsoft Windows XP Professional Service Pack 3 as the operating system, a computer with an Intel Core i7 920 processor (2.8 GHz), memory of 3G bytes, and Microsoft Visual C++ 2005 with Intel C++ compiler version 11.1. We used FFTW3.2.1 library [18] for the FFT operations on the CPU. In addition, we used one Nvidia GeForce GTX260 GPU board, and CUDA 2.3 for GPU programming.

The calculation conditions are as follows: the range of *d _{j}* is 2

*mm*≤

*d*≤ 12

_{j}*mm*,

*z*= 300

*mm*,

*λ*= 532

*nm*,

*p*= 8

*µm*and

*N*×

_{x}*N*= 2,048 × 2,048. We calculated CGHs from two 3D objects: “Nanami”, which is NHK (Japan Broadcasting Corp.) character, and “Earth”, consisting of 4,596 and 30,492 points, respectively.

_{y}Table 1 shows the CGH calculation times using each method. All of the calculations on the CPU used eight CPU threads. As shown in the table, for little points such as “Nanami”, the previous and improved WRP methods are ineffective since the FFT operations become a dominant factor. Conversely, for many object points such as “Earth”, the previous and improved WRP methods are very effective. In addition, the calculation time required for the improved WRP method on the CPU is faster than that for the conventional algorithm on the GPU.

Table 2 shows the CGH calculation times of the step 1 and 2 in the previous and improved WRP methods. As shown in the table, by adopting the LUT and the rectangular region, the calculation time of Step 1 of the improved WRP method is over four times faster than that of the previous method. And, by adopting the GPU, the calculation time of Step 2 of the improved WRP method is over 80 times faster than that of the previous method.

Figure 2 shows optical reconstruction movies from phase-type CGHs (Kinoforms) generated by the improved WRP method. Note that the current improved WRP method cannot treat an image occlusion technique [9]. This is our next work. We used a green semiconductor laser with the wavelength of 532 nm and a phase-modulated liquid crystal display (LCD) panel made by Holoeye, HEO-1080, with the pixel pitch of 8*µm* and 1,920 × 1,080 pixels. We displayed a portion of 2,048 × 2,048 pixels CGH on the LCD panel. The size of all the reconstructed images is approximately 16*mm*(*Width*) × 16*mm*(*Height*) × 10*mm*(*Depth*).

Comparing numerical reconstruction from a CGH generated by the conventional algorithm with that from a CGH generated by the improved WRP algorithm, we find a blur in micrometerorder in the improved WRP algorithm; however, we do not almost observe the blur.

## 4. Conclusion

In conclusion, we proposed an improved WRP method using LUT and GPU computing. The improved WRP method calculates a CGH faster than the previous method. Next we will apply an occlusion techinique to the improved WRP method. We thank NHK for use of the character. This work is supported by the Ministry of Internal Affairs and Communications, Strategic Information and Communications R&D Promotion Programme (SCOPE), 2009.

## References and links

**1. **C. Slinger, C. Cameron, and M. Stanley
M, “Computer-Generated Holography as a Generic Display Technology,” Computer **38**, 46–53 (2005). [CrossRef]

**2. **S. A. Benton et al., Holographic Imaging, (Wiley-Interscience, 2008).

**3. **M. Lucente, “Interactive Computation of holograms using a Look-up Table,” J. Electron. Imaging **2**, 28–34 (1993). [CrossRef]

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

**6. **H. Yoshikawa, T. Yamaguchi, and R. Kitayama, “Real-Time Generation of Full color Image Hologram with Compact Distance Look-up Table,” OSA Topical Meeting on Digital Holography and Three-Dimensional Imaging 2009, DWC4 (2009).

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

**9. **M. Lucente and T. A. Galyean, “Rendering Interactive Holographic Images,” Proc. of SIGGRAPH 95 387–394 (1995).

**10. **N. Masuda, T. Ito, T. Tanaka, A. Shiraki, and T. Sugie, “Computer generated holography using a graphics processing unit,” Opt. Express **14**, 587–592 (2006). [CrossRef] [PubMed]

**11. **L. Ahrenberg, P. Benzie, M. Magnor, and J. Watson, “Computer generated holography using parallel commodity graphics hardware,” Opt. Express **14**, 7636–7641 (2006). [CrossRef] [PubMed]

**12. **H. Kang, F. Yaras, and L. Onural, “Graphics processing unit accelerated computation of digital holograms,” Appl. Opt. **48**, H137–H143 (2009). [CrossRef] [PubMed]

**13. **Y. Pan, X. Xu, S. Solanki, X. Liang, R. Bin, A. Tanjung, C. Tan, and T. C. Chong, “Fast CGH computation using S-LUT on GPU,” Opt. Express **17**, 18543–18555 (2009). [CrossRef]

**14. **T. Shimobaba, T. Ito, N. Y. Ichihashi, and N. Takada “Fast calculation of computer-generated-hologram on AMD HD5000 series GPU and OpenCL,” Opt. Express bf 18, 9955–9960 (2010). [CrossRef]

**15. **T. Shimobaba, et.al., “Numerical calculation library for diffraction integrals using the graphic processing unit : the GWO library,” J. Opt. A: Pure Appl. Opt. **10**, 075308 (5pp) (2008). [CrossRef]

**16. **T. Shimobaba, Y. Sato, J. Miura, M. Takenouchi, and T. Ito, “Real-time digital holographic microscopy using the graphic processing unit,” Opt. Express **16**, 11776–11781 (2008). [CrossRef] [PubMed]