Real-Time Diffraction Field Calculation Methods for Computer-Generated Holograms

Holographic three-dimensional television systems provide a natural 3D visualization. Fast calculation of the diffraction field from a three-dimensional object is essential to achieve video rate. In the literature, there are myriads of fast algorithms for diffraction field calculation from three-dimensional objects, but most of them omit the pixelated structure of the dynamic display devices which are used in the reconstruction process. In this chapter, the look-up table-based fast algorithm for diffraction field calculation from a three-dimensional object for a pixelated dynamic display device is presented. Real-time diffraction field calculations are obtained by running the algorithm in parallel on a graphical processing unit. Performance of the algorithm is evaluated in terms of computation time of the diffraction field and the normalized mean square error on the reconstructed object. To have optimization on the required memory space for the look-up table, two different sampling policies along the longitudinal axis are implemented. Uniform sampling policy along the longitudinal axis provides better error performance than nonuniform sampling policy. Furthermore, optical experiments are performed, and it is observed that both numerical and optical reconstructions are similar to each other. Hence, the proposed method provides successful results.


Introduction
Holography is the only visualization technique that satisfies all the depth cues [1][2][3]. Therefore, it gives a natural three-dimensional (3D) visualization. Holography is based on capturing the diffracted optical waves from an object and regenerating those waves again by illuminating the recording media [1][2][3][4][5][6]. Captured optical waves provide a significant amount of information related to the object such as surface profile, depth, and refractive index of the object. Hence, holography has a myriad of applications. For instance, holograms can be used as optical elements like prisms, lenses, and mirrors [7,8]. Also, parallel optical computing is possible when holograms are employed [9,10]. Furthermore, holograms are useful in metrology [11][12][13] and microscopic imaging to visualize very small objects like cells and bacterias [14,15]. Another application of holography is related to nondestructive testing [16][17][18]. Nevertheless, major application of holography is related to 3D visualization, and it is used in education [19,20], dentistry [21,22], gaming [23], demonstration of cultural heritage [24], and more.
Holography setups can be assembled by using different configurations depending on the application. In optical holography setups, holographic patterns are stored on high-resolution holographic films [25,26] and some type of crystals [27]. However, in some of the applications, we need to process the captured holographic patterns by numerical methods. Then, digital sensing devices are employed as a capturing device. Those types of setups are called as digital holography, and it has a vast amount of applications especially in nondestructive testing and microscopy. In [28], digital holography-based measurement method of 3D displacement is presented. Observed material is illuminated from four different directions sequentially; then they are combined to improve the resolution in the order of 10 nm. As a nondestructive testing method, digital holography is used in the analysis of cortical bone quality and strength impact in [29]. Furthermore, a method based on digital holography is implemented for detecting and measuring effect of moisture on the hygroscopic shrinkage strain on wood [30]. Another application of digital holography is in precise and accurate measurement of the initial displacement of the canine and molar in human maxilla [31]. By using subpixel registration and fusion algorithms, an improvement of profile measurements and expanding the field of view (FOV) in continuous-wave terahertz reflective digital holography is achieved [32]. A comprehensive review of denoising methods on phase retrieval from digital holograms in terms of signal-to-noise ratio (SNR) and computation time is presented in [33]. Removal of phase distortions by using principal component analysis (PCA) method is given in [34].
Holography is a versatile tool for visualization, measurement, and testing. In optical and digital holography methods, we need some optical sensing elements like polymers and digital devices to capture the diffracted field from the object. However, in computer-generated holography (CGH), diffraction field calculations are performed by using numerical methods and signal processing algorithms [4][5][6]35]. Then, we can obtain the hologram from the calculated diffraction field and use it to drive dynamic display devices such as spatial light modulators (SLMs). After that, illumination of the SLM with a coherent light source will provide an optical reconstruction of the original object. When CGHs are calculated sequentially and used in driving SLMs, then we can have a holographic 3D television (H3DTV) as a product. An overview on holographic displays is presented in [36]. Generally, coherent light sources are used in H3DTV systems, and those light sources can generate speckle noise in the reconstructions. Low computational method for improving image quality and decreasing the speckle noise in CGH is proposed in [37]. Diffraction field calculations as in CGH are also used in other 3D display systems to improve the resolution of reconstructed objects. For instance, in integral imaging-based 3D display system, distortions on the elemental images are corrected by using holographic functional screen [38].
In diffraction field calculation from a 3D object, we have to generate a synthetic 3D object. There are plenty of ways for generating a synthetic 3D object in a computer. For instance, we can form a 3D object by using a set of point light sources which are distributed over the space. Those types of objects are called as point cloud objects. To calculate the diffraction field from a point cloud object, we superpose the diffraction fields emitted by each point light source [35,[39][40][41][42][43][44]. Another 3D object generation method is based on stitching small planar patches. As in the process of diffraction field calculation from point cloud objects, once again the diffracted fields from each patch are superposed to obtain the diffraction field of the object [45][46][47][48][49][50][51][52][53][54][55]. The third method which can be used in the generation of synthetic 3D object is based on having multiple two-dimensional (2D) cross sections of the object along the longitudinal axis. Then, superposition of diffracted fields from those 2D cross sections will give the diffraction field of the 3D object [56][57][58][59][60]. A detailed summary on CGHs in terms of resolution, field of view, eye relief, and optical setups for different 3D object generation methods can be seen in [61,62].
CGHs of the objects should be calculated rapidly to obtain H3DTV systems. Hence, fast methods such as fast Fourier transform (FFT) and look-up table (LUT)based methods are utilized in CGH calculations. In [39,52], algorithms which are based on FFT are used for decreasing the calculation time of CGH. Precomputed LUTs are also used for achieving fast calculations in CGH calculations [2,39,41,42,63,64]. Another way to achieve fast calculation in CGH is based on segmentation of diffraction field from point light sources [43,44]. Parallel processing of diffraction field calculation provides further improvements on the computation time. Graphical processing units (GPUs) are special hardware to run parallel calculations. Thus, they are one of the most convenient hardware for H3DTV systems [40,44,65,66]. Time-division method can also be used in the calculation of CGHs for layered 3D objects to achieve fast computations [67].
Imposing some approximations in the diffraction field calculations provides to decrease the computational complexity, and it paves the way to obtain fast diffraction field calculations. In the meantime, we have to improve the quality of the reconstructed object. An accurate calculation method of diffraction field which is based on angular spectrum decomposition is explained in [68]. Furthermore, diffraction field calculation methods for SLMs with pixelated structure are presented in [69][70][71]. However, the computational complexities of those methods are too high to have real-time diffraction field calculations. As a result of this, the algorithms presented in [72][73][74][75] are proposed as a solution to both computation time and quality in the reconstructed object in H3DTV. Further computational time improvements can be obtained by utilizing a LUT which is optimized for parallel processing on a GPU to achieve real-time calculations. Moreover, the pixel structure of the employed SLM in the reconstruction process is taken into account in forming LUT. Calculated LUT has one-dimensional (1D) kernels to decrease the allocated memory space.

Calculation of diffraction pattern used in driving SLM with pixelated structure
In CGH, it is possible to obtain 3D reconstructions of both synthetic and real objects. By employing dynamic display devices like SLMs in the reconstruction process, we can have H3DTV systems. To drive SLMs, we have to calculate diffraction fields from 3D objects by using numerical analysis methods and signal processing techniques. Calculation of diffraction field depends on the 3D object generation method. In this work, we assumed that 3D objects are represented as point clouds, because it is one of the simplest methods in 3D object generation. The diffraction field of the 3D object is calculated by superposing the diffraction fields emitted from the points that form the 3D object.
Superposition on diffraction field calculation from a point cloud object over a planar surface can be expressed as where ψ r 0 ð Þ and ψ r l ð Þ are the diffraction fields over SLM and diffraction field at l th sample point of the 3D object, respectively. Surface of SLM is represented by the position vector r 0 ¼ x; y; 0 ½ , and the sampling points of 3D object are shown by r l ¼ x l ; y l ; z l Â Ã . We assume that Fresnel approximation is valid and the term h F r ð Þ denotes diffracted field on the SLM from a point light source expressed as where r ¼ x; y; z ½ , k is the wave number, and λ is the wavelength of the light source used in illumination of the object.
Scaled and superposed diffraction fields from point light sources provide the diffraction field of the 3D object, and its phase component is used for driving the SLM. Then, entire surface of the SLM is illuminated by a plane wave. After that, the reflected optical wave from the surface of the SLM generates an optical replica of the 3D object. Most of the off-the-shelf SLMs have square pixels with very high filling factors like 93% [76]. Hence, the filling factor in the simulated SLM is approximated as 100%. The pixel structure of the simulated SLM is illustrated in Figure 1.
Simulation of optical setup can be improved when the pixelated structure of the SLM is taken into consideration. For that purpose, we have to perform surface integration over each pixel area on the SLM. It is assumed that gray value over each pixel area has a constant value. The diffraction field over SLM can be found as where n and m stand for indices of SLM along x-and y-axes, respectively. It is also possible to represent Eq. (3) by scaling and superposing 2D kernels K α l , 2D , where P r l ð Þ ¼ Àjψ r l ð Þe jkz l and α l ¼ 1 ffiffiffiffi λz l p 2D kernel K α l , 2D can be decomposed into 1D kernels as where x l and y l refer to locations of l th point light source, used in generation of 3D object, along x-and y-axes, respectively. Each 1D kernel K α l , 1D can be represented as and its elements can be calculated as where ζ l, n ¼ x n Àx l √λ z l . The operators C Á ð Þ and S Á ð Þ stand for cosine and sine Fresnel integrals, respectively [5,6], and they are calculated as Numerical evaluation of cosine and sine Fresnel integrals given in Eq. (8) is calculated by adaptive Lobatto quadrature [77].
In the standard algorithm, diffraction field of each point is obtained by evaluating Eq. (8). Then, superposition of those fields is performed to obtain CGH. As a result of this, computational complexity of diffraction field calculation is too high to have real-time applications. As a solution to the computation time problem, we present a fast algorithm to calculate 2D kernel, K α l , 2D , based on LUT and parallel processing.

Proposed algorithm for fast calculation of CGH
Fast computation of diffraction field and improved quality of reconstructed 3D object are essential issues in H3DTV. As a solution to those problems, we propose a method based on calculation of 2D kernels K α l , 2D without evaluating sine and cosine Fresnel integrals. To achieve fast calculation, precomputed LUT is utilized, and the diffraction field of the 3D object can be obtained by scaling and superposing the 2D kernels K α l , 2D asψ whereψ 2D, z¼0 denotes the estimated diffraction field of the 3D object on the SLM andK α l , 2D is the 2D kernel which denotes the diffraction field of l th point of the 3D object on SLM. 2D kernelK α l , 2D is calculated by multiplying 1D kernels K α l , 1D from LUT as shown in Eq. (5). Each 1D kernel K α l , 1D represents the diffraction field on SLM from specific depth along longitudinal axis. A simple arithmetic operation is used for speeding up data fetching from the LUT. As result of this, total computation time of the diffraction field can be improved in terms of data fetching. By increasing the number of precomputed 1D kernels in LUT, we can achieve better diffraction field estimations for proposed method, but it causes to allocate more memory space. Hence, we apply different sampling policies along longitudinal axis to optimize memory space allocation. In the first sampling policy, uniform sampling along longitudinal axis is performed. In the second sampling policy, we sample the parameter α l ¼ 1 ffiffiffiffi ffi λz l p uniformly. Thus, we have nonuniform sampling along the longitudinal axis.

Simulation results
Performance assessment of the proposed diffraction field calculation method is obtained by implementing different scenarios, but a few of them are presented to give an insight to the reader. Two major performance evaluation criteria are taken into account: total computation time of the CGH and the normalized mean square error (NMSE) on the reconstructed object. NMSE on the reconstructed object can be calculated as where ψ 2D, z¼z 0 n; m ð Þandψ 2D, z¼z 0 n; m ð Þdenote reconstructed objects at z ¼ z 0 plane from the diffraction field calculated by the standard and the proposed algorithms, respectively. Simulated scenario for a CGH is illustrated in Figure 2.
First, a 3D point cloud object is generated in computer environment. The generated 3D object has 3144 points which are distributed over the space. The volume occupied by the object has x e ¼ 2:8mm, y e ¼ 4:1mm, and z e ¼ 4:1mm extensions along x-, y-, and z-axes. There is a distance between the object and the screen, and it is taken as z 0 ¼ 61:6mm. We assume that simulated SLM has 100% fill factor and pitch distance X s is taken as 8 μm. Also, the simulated SLM has 512 pixels along both xand y-axes, respectively. We assume that green laser is employed for illumination purpose; hence the wavelength is taken as 532 nm.
The proposed algorithm is implemented by using two platforms: MATLAB and Visual C++. To have shorter computation time for diffraction fields, we utilize CUDA libraries and parallel computation power of GPU. The assembled computer system has i5-2500 CPU at 3.3 GHz, 4GB RAM, and a GTX-680 GPU to run the algorithm. Operating system of the computer is chosen as 64-bit Windows 7. Generally, off-the-shelf SLMs have pixelated structure, and phase parts of the calculated diffraction fields are used for driving the SLM. When the pixelated structure of SLM is not taken into account in CGH calculations, it is not easy to differentiate focused and unfocused parts of the reconstructed 3D objects. An illustration of such a result can be seen in Figure 3a. As a result of the similarity in focused and unfocused parts, the quality of the reconstructed object is decreased significantly. On the contrary, the difference between focused and unfocused parts in the reconstructed 3D object is clear when the proposed method is used in diffraction field calculation. Those results can be seen easily in Figure 3b.
Furthermore, numerical and optical reconstructions are very similar to each other, and that similarity in the reconstructions can be seen in Figure 4.
To calculate the diffraction field in the standard method, we need to perform cosine and sine Fresnel integrals for each pixel on SLM and for each point light source in 3D object. As a result of this, computational complexity of the standard method is extremely high, and CGH is calculated at 2701:10 s. Significant improvement in computation time can be achieved when the proposed algorithm is  employed in CGH calculation. When we use LUT-based method for the same scenario which is mentioned above, we need 8:15 s to calculate the CGH. Further improvement in computation time can be obtained if the presented algorithm is implemented in parallel on a GPU. Although, significant gain on the computation time of CGH is obtained by using LUT, there will be negligible amount of error on the reconstructed objects, because of having finite number of kernels and the quantization effect along the longitudinal axis. The performance of the proposed method is summarized in Table 1.
By increasing the number of kernels in LUT, we can improve error performance of the algorithm without having any extra computational load, but there is an increase in the size of the required memory. As a result of this, installed memory space may not be enough to perform the diffraction field calculations. To overcome memory allocation problem, we use another sampling policy in generation of LUT. Two different sampling policies along the longitudinal axis are proposed. The first sampling policy is based on uniform sampling of longitudinal axis. The second sampling policy is related to uniform sampling of α l . Hence, there will be nonuniform sampling along the longitudinal axis. Tables 2 and 3 summarize performances of the sampling policies in terms of NMSE and required memory allocation by the precomputed LUT. As it can be seen from Tables 2 and 3, when the size  of the LUT is fixed, uniform sampling policy along longitudinal axis provides better NMSE performance than nonuniform one. In terms of calculated numerical errors, there should be a significant amount of deviation between reconstructed objects from CGHs obtained by standard and proposed method, but it is not easy to differentiate the reconstructions visually. Illustrations of numerically reconstructed objects by using both methods are shown in Figure 5a and b, respectively. To see the difference between to reconstructions, we subtract two reconstructions from each other and then take the magnitude of that difference. Then, we scale difference image linearly between 0 and 255 to improve the visibility of insignificant deviations. Those deviations can be seen in Figure 5c. Most of the deviations are in the unfocused region, and those deviations will not decrease the quality of the reconstruction. As a result of this, the proposed algorithm provides successful results. LUT is formed by uniform sampling of α l parameter. Each element in 1D kernels is represented by four bytes. Table 3.
Performance of the proposed algorithm according to the number of kernels used in LUT, NMSE, and allocated memory space. Performance assessment of the presented algorithm is tested by optical reconstructions as well. For that purpose, we assembled an optical setup which is shown in Figure 6. Green laser with λ ¼ 532 nm is used as a coherent light source, and HoloEye Pluto phase-only SLM is employed as a dynamic display device. A couple of optically reconstructed objects are shown in Figure 7.

Conclusions
Two major problems in H3DTV systems can be called as decreasing the computation time of CGH and improving the quality of the reconstructed object. Using  fast algorithms in diffraction field calculations will be helpful to decrease the computation time, but most of those fast algorithms impose some approximations that decrease the quality of the reconstructed object. In this work, we propose a diffraction field calculation algorithm that paves the way to achieve real-time calculations of diffraction fields from point cloud objects. In the meantime, the quality of the reconstructed objects is improved by taking into account the pixelated structure of SLM. Also, the proposed method can be run in parallel on a GPU. Performed numerical and optical experiments provide similar results. The proposed method utilizes precomputed LUT to decrease the computational load. To store the precomputed LUT, we need significant amount of memory allocation, and optimization of the occupied memory space is obtained by having two different sampling policies along the longitudinal axis. In the first sampling policy, LUT is formed by having uniform sampling along longitudinal axis. In the second one, nonuniform sampling is applied. When we fix size of the LUT, better NMSE performance is obtained by uniform sampling policy. As a result of this, when we use uniform sampling policy in computation of LUT, we need to allocate less amount of memory to store it.