Fast CGH computation using S-LUT on GPU

1. Introduction

Holographic 3D display uses diffraction effect from holograms to reconstruct 3D objects. It is a very promising true 3D display technology, which has capability to provide wide color gamut and full depth cues, without using any glasses [1]. However, generation of holograms has heavy computational load [2]. To solve this problem needs fast algorithm and powerful computer hardware.

There are different kinds of algorithms for computer generated hologram (CGH) [2,3]. Fourier Ping-Pong algorithm was used to create Fourier holograms for holographic 3D display [4,5]. In this algorithm, one hologram-size fast Fourier transform operation is required for each slice of object, thus it is slow for objects with many depth layers. Lucente [6] developed diffraction specific algorithm and significantly reduced CGH computational load. But it is more suitable for horizontal parallax only (HPO) holograms and display system with high horizontal resolution. Coherent ray trace (CRT) algorithm simulates light propagation from all sampling object points to all hologram pixels [7]. It can achieve very high quality image reconstruction with full parallax, while keeping low memory consumption. However, because of direct many-to-many convolution, CRT algorithm is very time consuming. Look-up table (LUT) was used to reduce the computational load of CRT algorithm by pre-calculating possible results and storing them in a table for in-line look-up [8]. As the amount of possible results is huge, the size of look-up table is big. Kim et al. [9] implemented a novel look-up table (N-LUT) method to reduce memory usage but it is still in the order of gigabyte (GB).

All CGH algorithms need to trade off among reconstruction quality, computational speed and memory usage. For the holographic 3D display system developed in Data Storage Institute (DSI) [10], we need to compute full parallax holograms with fast speed on PC, even for big holograms containing complicated 3D object specifications. Thus, we have developed a new algorithm using split look-up tables (S-LUT), which can meet the requirements of high speed, low memory usage and full parallax reconstruction.

Using powerful computer hardware is another way to increase CGH computational speed. In the past, one used super computers / clusters for this task due to its heavy computational load [2,3]. Special-purpose computational chip was also designed to increase the speed [11]. Recently there is a trend to use GPU for CGH computation [7,12,13]. In this paper, we implement CGH algorithms on GPU in order to further increase the computational speed. But instead of using GPU’s graphic pipeline / OpenGL [7,12,13], we use C + + -like parallel GPU application programming interface (API), Compute Unified Device Architecture (CUDA) [14], as CUDA is designed especially for direct GPU computing without limitations by graphic pipeline / shader.

2. Conventional CRT and LUT algorithms

Conventional CRT algorithm directly simulates physical light propagation process from each object point in 3D space to each hologram pixel. Its mathematical formula can be written as in Eq. (1) [7,13]:

Algorithm steps of CRT follow Eq. (1), and can be listed as:

From the algorithm described above, the most inner loop ① runs $XYN$ times. The total computational complexity is in $O\left(XYN\right)$, where $O\left(\right)$ is the big O notation. For each run of loop ①, operations include one cosine, one square root, five additions and five multiplications. These operations introduce heavy computational load, especially cosine and square root operations, and make the algorithm run very slow. On the other hand, CRT algorithm has some advantages such as no additional memory requirement other than input 3D object specification and output hologram, as well as high reconstruction quality.

In order to solve slow speed problem of CRT, LUT algorithm [8,9] had been developed. In LUT, all possible results of the cosine function in Eq. (1) are off-line pre-computed and stored in a table; later whenever a result is needed for in-line computation, it will be directly read out from the table. Original LUT uses HPO, and reduces CGH computational load by one order [8] while keeping affordable memory requirement. But using HPO is considered as downgrading of reconstruction quality, as HPO lacks vertical parallax. If using LUT for full parallax hologram computation, memory requirement will dramatically increase.

The mathematical detail of LUT can be described in Eq. (2):

Converted into algorithm, it becomes:

Although the most inner loop ② of in-line computation part also runs $XYN$ times. The total computational complexity is also in $O\left(XYN\right)$. Each time only one multiplication and one addition operations are performed. The rest of computational operations are done by the off-line pre-computation in LUT. Therefore, the computational speed of LUT algorithm is much faster than CRT. Considering the memory requirement of LUT, $x_h$ values are from 0 to $d{x}_{h}X$ with stepping $d{x}_h$, and $x_j$ values are from 0 to $d{x}_{o}W$ with stepping $d{x}_o$, where $d{x}_h$ is the pixel size of SLM, $d{x}_o$ the sampling interval of object space in horizontal direction and W the number of object space sampling points along horizontal direction. For typical holographic display system, $d{x}_o\gg d{x}_h$, and in our display system W and X are in the same order of magnitude. Thus, the minimum of $\Delta x$ is $-d{x}_{o}W$ and the maximum is $d{x}_{h}X$, and the stepping of $\Delta x$ is $d{x}_h$. As a result, the table size in $\Delta x$ dimension is $\left(d{x}_{o}W+d{x}_{h}X\right)/d{x}_h=d{x}_{o}W/d{x}_h+X$. Similarly, the table size in $\Delta y$ dimension is $d{y}_{o}H/d{y}_h+Y$, where $d{y}_o$ is the sampling interval of object space, H the number of object space sampling points, $d{y}_h$ the pixel size of SLM, all in vertical direction. In total, the table size is $\left(d{x}_{o}W/d{x}_h+X\right)\left(d{y}_{o}H/d{y}_h+Y\right)D$, where D is the number of object space sampling points along depth direction (depth layers). Normally $d{x}_h$ and $d{y}_h$ are in the order of micrometer while $d{x}_o$ and $d{y}_o$ are in the order of sub-millimeter, so the table size can easily boost and uses up all available memory space of a PC. This limits the application of LUT algorithm although it is faster than CRT. However the reconstruction quality of LUT is as good as that of CRT because their computational results are mathematically identical.

3. A novel S-LUT algorithm

In order to solve the speed problem of CRT and the memory problem of LUT without sacrificing reconstruction quality, we propose a new algorithm using look-up tables for horizontal and vertical directions, separately. We name it split look-up table (S-LUT) algorithm.

Recalled from Eq. (1), the resulted hologram can be described as:

As in our display system, reconstruction size is small as compared to object-hologram distance, i.e. $\left|\Delta x\right|\ll z_j$ and $\left|\Delta y\right|\ll z_j$, using Fresnel approximation, Eq. (3) can be approximated as Eq. (4):

Because wavelength is much smaller than the object-hologram distance, i.e. $\lambda \ll z_j$, small change of $z_j$ in the order of λ will not affect the reconstruction quality. Therefore $z_j$ can be adjusted to the nearest multiple of λ. Under this condition, $2\pi z_j/\lambda$ is multiple of $2\pi$, and one more $2\pi z_j/\lambda$ can be added into Eq. (4) without changing the value of $I\left(x_h,y_h\right)$. As a result, Eq. (4) becomes Eq. (5):

Using inverse Fresnel approximation on Eq. (5) results in Eq. (6):

Rewrite Eq. (6) in complex form as Eq. (7):

Let us define horizontal light modulation factor as$H(\Delta x,z_j)=e^{i\left(\genfrac{}{}{0.1ex}{}{2\pi }{\lambda }\sqrt{\Delta x^2+z_j{}^2}\right)}$ and vertical light modulation factor as $V(\Delta y,z_j)=e^{i\left(\genfrac{}{}{0.1ex}{}{2\pi }{\lambda }\sqrt{\Delta y^2+z_j{}^2}\right)}$, then Eq. (7) can be rewritten as Eq. (8):

For n object points falling on the same vertical line $\left(x_j,z_j\right)$, $H\left(\Delta x,z_j\right)$ is all the same, thus the hologram resulted from these n object points can be calculated by Eq. (9):

As ${\displaystyle \colorbox[rgb]{}{${\sum }_{j=0}^{n-1}{a}_{j}V\left(\Delta y,z_j\right)$}}$ is independent of $\Delta x$, Eq. (9) can be broken down into two steps as Eq. (10) and Eq. (11):

The final hologram $I_c\left(x_h,y_h\right)$ is the sum of all contributions from object points on different vertical lines and can be computed by Eq. (12):

Therefore our new algorithm S-LUT can be listed as:

The list of object points falling on the same vertical line is generated by the program in input processing during in-line computation. The n object points with non-zero intensity are marked down and stored in an array for each vertical line together with respective y coordinates.

In in-line computation of S-LUT algorithm, for a vertical line containing n object points with non-zero intensity, loop ③ runs $nY$ times and loop ④ runs $XY$ times. Therefore, total in-line computational complexity is in $O\left(nY+XY\right)$, and during each step of both loop ③ and loop ④, only one multiplication and one addition operations need to be performed. Similar to memory analysis of LUT, the table size of $H\left(\Delta x,z_j\right)$ in $\Delta x$ dimension is $d{x}_{o}W/d{x}_h+X$, and that of $V\left(\Delta y,z_j\right)$ in $\Delta y$ dimension is $d{y}_{o}H/d{y}_h+Y$, while the size in $z_j$ dimension is D. As a result, the total size of two tables is $\left(d{x}_{o}W/d{x}_h+X+d{y}_{o}H/d{y}_h+Y\right)D$. Reconstruction quality of S-LUT is as good as CRT, if all conditions, $\left|\Delta x\right|\ll z_j$, $\left|\Delta y\right|\ll z_j$, and $z_j$ being multiple of λ, are met.

The comparison in computational complexity, operation and memory usage for CRT, LUT and S-LUT algorithms is listed in Table 1 , for n object points falling on the same vertical line.

Table 1. Complexity, operation and memory usage comparison among algorithms

View Table | View all tables in this article

From Table 1, it can be seen that both computational complexity and memory usage of S-LUT are one order less than that of LUT, i.e. S-LUT runs much faster than LUT, while using much less memory.

4. Implementation of algorithms on GPU

By using CUDA, mathematical operations can be performed on GPU cores in parallel, without mapping them to graphic pipeline. As generally GPU has more cores than CPU (30 computational cores for nVidia GTX285 GPU [15] and 4 cores for Intel core i7 965 CPU [16]), parallel performance of GPU is better than that of CPU, although individual core speed of GPU is slower than that of CPU (680 MHz for MSI GTX 285 OC [17], 3.2 GHz for Intel core i7 965 [16]). Multiple GPUs can be used in a PC, which provides more GPU cores to speed up CGH computation. Another advantage of GPU over CPU is that each GPU has its own memory interface to its dedicated memory chips. As a result, memory for different GPUs can read/write simultaneously. Currently only global memory on GPU is used to store the pre-computed data of $H\left(\Delta x,z_j\right)$ and $V\left(\Delta y,z_j\right)$. In the future, shared memory on GPU may be used for further computational speed optimization. However CPU only has one interface to system main memory and data for CPU cores need to be queued up, which may create computational bottleneck.

GPU also has some disadvantages. High speed computation on GPU is mainly for single precision calculation. Error rate of GPU and its memory is higher than that of CPU and main memory. Therefore, direct computation on GPU without error correction may not be suitable for high precision scientific computation. As hologram has high error tolerance, GPU still can meet the requirement of CGH computation, especially when resulted holograms need to be binarised before they are lunched to SLM for reconstruction.

Using a multi-GPU PC for CGH computation needs to consider how to distribute and balance workload among GPUs. In our implementation, the program will distribute loads to each CUDA-enabled GPU in unit of vertical line in object space, in such a way that computation time on each GPU is as close to each other as possible.

Detailed implementation steps of S-LUT algorithm on GPU are listed below, where m is the number of CUDA-enabled GPUs in the PC:

From steps listed above, it is noticed that additional memory needs to be allocated for workload $L_k$ and holograms $I_c\left(x_h,y_h\right)|{}_{GPUk}$ of individual GPUs in main memory. This increases memory requirement by $WHDm+XYm$ as compared to that running on CPU.

5. Experimental results

Our computational experiments are done on a PC with specifications listed in Table 2 . Our program uses only single CPU core and total 3 × 30 GPU cores. Some parameters used in CGH computation are listed in Table 3 .

Table 2. Specifications of PC

View Table | View all tables in this article

Table 3. CGH computation parameters

View Table | View all tables in this article

The time for calculation of $H\left(\Delta x,z_j\right)$ and $V\left(\Delta y,z_j\right)$ is shown in Fig. 1 . It can be seen that the time needed for this off-line calculation is proportional to the number of depth layers.

 

Fig. 1 Time for calculation of $H\left(\Delta x,z_j\right)$ and $V\left(\Delta y,z_j\right)$for S-LUT.

Download Full Size | PPT Slide | PDF

The comparison of computational speed among different algorithms implemented on both CPU and GPU is shown in Fig. 2 . The detail of computation time of S-LUT on GPU is shown in Fig. 3 . The object points are randomly distributed in the object space (Table 3) with 21 depth layers.

 

Fig. 2 Comparison of computational speed among algorithms on CPU and GPU.

Download Full Size | PPT Slide | PDF

 

Fig. 3 Computation time of S-LUT on GPU.

Download Full Size | PPT Slide | PDF

From CPU computational results shown in Fig. 2, it is clear that S-LUT is much faster than CRT and LUT, especially when the number of object points increases. LUT also shows faster computation as compared with CRT. Running time of S-LUT depends mainly on the number of vertical lines rather than the number of object points in object space as described in Table 1. As a result, running time for S-LUT keeps almost constant when the number of object points goes beyond a certain value, which in our test cases is 10k. Beyond this value, almost all vertical lines in object space have at least one object point and need to be computed for different number of object points above 10k.

From GPU computational results shown in Fig. 2, the speed of CRT on GPU increases drastically as compared with that on CPU. CUDA automatically optimizes the execution of kernel code running on GPUs to cover read / write time of GPU memory. This optimization effect is more obvious for larger number of object points for the cases of CRT and LUT on GPU. Therefore, their computation time does not have linear relationship with the number of object points. From Fig. 2 and Fig. 3, the speed of S-LUT on GPU is about 100 times faster than S-LUT on CPU and CRT on GPU. It is shown that S-LUT is the fastest among all algorithms. However, LUT is even slower than CRT on GPU. This is because LUT reads huge amount of data from large area in memory on GPU, and different GPU cores read data simultaneously, which causes a bottleneck in memory bandwidth. As for S-LUT, there is no such a bottleneck because it only reads small amount of data from small area in memory, as shown in Table 1.

As shown in Fig. 3, S-LUT on GPU achieves 0.3 second running time for less than 1k object points and 0.5 second for more than 10k object points. For the number of object points larger than 40k, S-LUT on GPU is faster than LUT on CPU by at least 700 times.

The comparison of memory usage is shown in Fig. 4 , in which base memory refers to the memory used to store input 3D object specification and final output hologram.

 

Fig. 4 Comparison of memory usage among algorithms on CPU and GPU.

Download Full Size | PPT Slide | PDF

LUT on both CPU and GPU uses much more memory space than CRT and S-LUT. It reaches 12GB physical memory limitation of our PC at 90 depth layers. With virtual memory support, LUT can handle up to 200 layers by using 27GB memory, but the computation as well as the whole PC slows down, and further increase in the number of depth layers stops the OS from responding. S-LUT on CPU uses very small amount of memory and only needs 1.2GB for 20k depth layers. S-LUT on GPU needs roughly twice of base memory because tables of object points are created to optimize overall performance of algorithm. However it still keeps memory usage at 7.5GB for 10k depth layers, which are far more than enough for normal CGH computation.

In order to test memory limitations, we purposely increase resolutions of SLM and object space to four times in both horizontal and vertical directions, i.e. 4096 × 3072 for hologram and 1200 × 1200 for object space, and the comparison is shown in Fig. 5 . Under this condition, LUT on both CPU and GPU stops at 8 depth layers with 17.5GB memory usage; S-LUT on CPU stops at 1490 depth layers with 474MB memory for tables $H\left(\Delta x,z_j\right)$ and $V\left(\Delta y,z_j\right)$ because at 1490 depth layers base memory occupies available memory space of 8GB. S-LUT on GPU uses 12GB at 1000 depth layers, and stops at 1300 depth layers with 15GB memory usages. This shows that S-LUT on GPU can meet most of application requirements within main memory limitation.

 

Fig. 5 Comparison of memory usage for high resolution SLM and object space.

Download Full Size | PPT Slide | PDF

Based on the results presented above, S-LUT on GPU has shown the fastest speed while maintaining low memory usage. It is most efficient in addressing the balance between computational speed and memory usage among all the algorithms investigated in this paper.

Figure 6 shows objects reconstructed from CGHs computed by CRT on CPU and S-LUT on GPU with our holographic 3D display system, in which the SLM device has 1920 × 1080 pixels, and its pixel size is 10.8 µm in both horizontal and vertical directions. The pictures were captured by a 2D SLR camera. Random-phase noise reduction technique [18] is adopted to improve reconstruction quality for all the objects. By naked eyes or even camera, no difference in reconstruction is observed between above two algorithms for both 2D (resolution chart, 372k object points) and 3D (teapot, 100 depth layers, 37k object points) objects. It is shown that S-LUT algorithm implemented on GPU can preserve full parallax and high reconstruction quality of original CRT algorithm.

 

Fig. 6 Reconstructed objects: (a) 2D resolution chart by CRT on CPU, (b) 2D resolution chart by S-LUT on GPU, (c) 3D teapot by CRT on CPU, and (d) 3D teapot by S-LUT on GPU.

Download Full Size | PPT Slide | PDF

Figure 7 shows a snapshot from a short video clip of a rotating 3D globe (150 depth layers, 62k object points). It is displayed at video rate by our holographic 3D display system [10]. The holograms used to reconstruct this dynamic 3D globe are computed by S-LUT on GPU.

 

Fig. 7 (Media 1) Video clip of a rotating 3D globe reconstructed from holograms computed with S-LUT on GPU.

Download Full Size | PPT Slide | PDF

6. Conclusion

We have developed a new algorithm by using split look-up tables (S-LUT) and successfully implemented it on GPU using CUDA. Comparisons among different algorithms implemented on both CPU and GPU have also been done in terms of computational speed, memory usage as well as reconstruction quality. S-LUT on GPU can achieve a speed faster than LUT on CPU by at least 700 times for the number of object points larger than 40k under the specified system specifications. It is also about 100 times faster than CRT on GPU. The memory usage of S-LUT on GPU is much less than that of LUT, which enables fast CGH computation for complicated 3D objects. The reconstruction of objects from holograms computed with S-LUT on GPU has shown full parallax and the same image quality as that with CRT on CPU. S-LUT algorithm implemented on GPU is very promising to meet the requirements of high speed, low memory usage and high quality full-parallax reconstruction for real-time and interactive holographic 3D display in the future.