1. Introduction

Holography is a technique to record and reconstruct a three-dimensional scene on a 2D plane (i.e. a hologram) based on optical diffraction and interference. In the early stage, holograms can only be acquired optically and recorded on photographic films as hardcopies. With the development of electronic and computer technology, a hologram can be generated numerically from a virtual 3D object by mathematically simulating the optical wave propagation on computers, known as Computer Generated Hologram (CGH). A CGH can be displayed on a spatial light modulator (SLM) to reconstruct the 3D object scene represented by the CGH. CGH can eliminate the difficulty of complicated instrumental setups in optical recording and have advantages of convenient generation, processing and storage. On the downside, the heavy amount of computation in CGH generation is a major challenge. In many dynamic holographic applications, such as 3D movie and 3D television, real-time CGH generation at high speed is favorable. To fulfill such requirement, it is necessary to develop fast CGH generation algorithms to simply the computation and accelerate the computing speed.

In CGH computation, point light source model is commonly adopted to model the virtual 3D object. The 3D object is assumed to be composed of many self-illuminated points and each point contributes an elementary hologram, known as Fresnel zone plate. The entire CGH can be obtained by overlapping the Fresnel zone plates from each individual point source.

For CGH computation based on point light source model, look-up table approaches are proposed to effectively reduce the computation time. The look-up table (LUT) approach is first designed by pre-computing all possible elementary CGH patterns off-line and storing them in a huge look-up Table [1]. For each object point, the corresponding elementary fringe pattern can be directly retrieved from the look-up table when a hologram is calculated. This method can accelerate the computing speed significantly but the drawback is tremendous memory usage.

To address this problem, a novel look-up table (N-LUT) [2] approach is proposed to reduce the data amount of LUT. In the original LUT method [1], an elementary hologram pattern is included in the look-up table for every point located at different horizontal and vertical positions of each depth plane. In N-LUT method, only the CGH pattern of center point on each depth plane is pre-computed. CGH patterns of the other points on the same plane is obtained by shifting the pattern of central point. N-LUT table method can significantly lower memory usage at the expense of extra shifting operations.

Despite the success, the total data size of N-LUT table is still quite large when there are many depth planes in the object scene. The memory size of N-LUT method can be further reduced by methods such as rearrangement of horizontal, vertical and longitudinal light modulation factors [3,4], trigonometric decomposition [5], single line representation in circular symmetry [6–8]. Among these methods, compressing N-LUT to a single line based on circular symmetry (referred to as “radial symmetric interpolation method” in this paper) [6–8] can generally achieve the highest compression ratio of memory storage size. In this paper, we notice that there is still data redundancy after the look-up table is already compressed by radial symmetric interpolation method. The size of look-up table can be further compressed by as much as 10 to 20 times with our proposed mini look-up table approach based on Principal Component Analysis (PCA).

2. Look-up table methods for fast CGH calculation

In point light source model, it is assumed that a three-dimensional (3D) object is composed of N points and the jth point has an intensity of $a_j$ located at the position of $(x_j,y_j,z_j)$, where $x_j$ and $y_j$ are horizontal and vertical position and $z_j$ is the depth position. The complex hologram $H(x,y)$ located at $z=0$ can be computed by Eqs. (1) and (2)

The basic ideas of N-LUT method [2] and radial symmetric interpolation method [6–8] are briefly demonstrated below. In N-LUT method, the elementary hologram (Fresnel zone plate) of the central point (when $x_j=0$ and $y_j=0$) at each depth plane $z_k$ (it is supposed that there are totally M depth planes, $1\le k\le M$) is pre-calculated as a look-up table $F{C}_k(x,y)$, given by Eqs. (3) and (4).

The Fresnel diffraction pattern of an arbitrary object point $D(x_j,y_j,z_j)$ can be computed by a multiplication of intensity $a_j(x_j,y_j,z_j)$ with a shifted $F{C}_k(x,y)$ when $z_j=z_k$. For example, the elementary hologram contributed by the point $a_j(x_j,y_j,z_j)$ can be synthesized by Eq. (5)

Suppose the size of one elementary hologram is $X\times Y$, there are totally M depth planes, the total size of look-up table is $X\times Y\times M$.

In radial symmetric interpolation method, the representation of 2D elementary hologram (Fresnel zone plate) at each depth plane is simplified to be a one single line due to its circular symmetry property, given by Eqs. (6)-(8)

3. Look-up table compression based on Principal Component Analysis

It is assumed the total number of depth planes is M and the number of r values is Q in each radius line at different depth planes. The entire look-up table FC in radial symmetric interpolation method, given by Eq. (9), can be considered as a two-dimensional (2D) complex number array of size Q × M. Each row in the array stores a set of radius line data for the Fresnel zone plate at the corresponding depth. Alternatively, FC can be regarded as two 2D real number arrays of size Q × M, $F{C}_{real}$ for real part and $F{C}_{imag}$ for imaginary part of the complex table correspondingly, given by Eq. (10).

One such 2D array (Re(FC) or Im(FC)) can be visualized as an intensity image, shown in Fig. 1. It can be observed that similar image patterns exist across different parts of the image. The radius lines at varying depth are highly correlated with each other instead of independent random lines. This implies there is still much data redundancy in the entire look-up table of radial symmetric interpolation method and data compression can be performed. Common data compression methods [10] such as Discrete Cosine Transform (DCT, adopted in JPEG image compression standard) and Discrete Wavelet Transform (DWT, adopted in JPEG 2000 image compression standard) can be employed to compress the lookup table data. Both DWT and DCT methods are capable of compressing natural images such as photographs with high compression ratio and satisfactory image quality. However, the texture pattern of the lookup table in this paper has unique data characteristics and general compression methods such as DWT and DCT do not give optimal result. Principal Component Analysis (PCA) method [9] is proposed to reduce the data size of look-up table in radial symmetric interpolation method in this paper, which demonstrates better performance than DWT and DCT methods, shown in Section 4.1. The entire flowchart of our proposed highly compressed mini look-up table incorporated with radial symmetric interpolation method for CGH generation is shown in Fig. 2.

 

Fig. 1 One portion of the visualized intensity image representing the real part of a look-up table in radial symmetric interpolation method

Download Full Size | PDF

 

Fig. 2 Flowchart of proposed highly compressed look-up table incorporated with radial symmetric interpolation method for CGH generation

Download Full Size | PDF

In our proposed scheme, a 2D rectangular array Re(FC) or Im(FC) with size Q × M is first partitioned into non-overlapping rectangular blocks of size a × b. Each block $S_w$ has $L=a\times b$ pixels and each pixel is denoted by $S_{wl}(1\le w\le W,1\le l\le L)$. There are totally a number of $W=(Q\times M)/L$ blocks in the 2D array. An example is illustrated in Fig. 3, where Q = 6, M = 4, a = 3, b = 2, L = 6 and W = 4.

 

Fig. 3 Partition of a 2D look-up table (size 6 x 4) into 4 non-overlapping rectangular blocks: (a) before partition; (b) after partition.

Download Full Size | PDF

Principal component analysis (PCA) [9] is a statistical method that employs an orthogonal transformation to convert a set of observations of correlated variables into a set of values of linearly uncorrelated variables named principal components. In our case, each block is regarded as one variable and the original 2D array contains W variables totally. Since the data blocks have similarity and duplication, they are highly statistically correlated. There are L pixels in each block, which indicates that each variable has L observations. PCA is capable of converting W correlated variables into $W_{pca}$ principal components and each principal component is referred to as $C_p(1\le p\le W_{pca})$, based on L observations of each variable. Each principal component has the same number of observations as original vectors. $C_p$ is composed of L intensity values $C_{pl}(1\le l\le L)$ corresponding to the L observations in arbitrary original data block $S_w$.

$W_{pca}$ can be significantly smaller than W but the $W_{pca}$ principal components can account for most of the variability in the original data (all W variables). All the data in original 2D array can be approximately represented by a linear combination of the $W_{pca}$ principal components, which is a data compression process. The $W_{pca}$ principal components can be derived by different methods such as matrix eigenvalue decomposition [9]. For brevity, the details shall not be discussed in this paper. The weighting parameters of linear combination can also be obtained when the $W_{pca}$ principal components are derived.

In the decompression process to recover original data, each observation $S_{wl}$ in each variable $S_w$ can be represented by a multiplication of corresponding set of weighting parameters ${\mu }_{wp}(1\le w\le W,1\le p\le W_{pca})$ with the elements in corresponding principal components, plus the average value of each corresponding observations in all variables $f_l(1\le l\le L)$, given by Eq. (11) and (12).

In PCA compression, the original data size is $W\times L(=Q\times M)$ and the compressed data size is $L+W{W}_{pca}+W_{pca}L$, since $f_l$, ${\mu }_{wp}$ and $C_{pl}$ values all need to be stored in the compressed look-up table to recover the original data. It shall be noted that the values of $W_{pca}$ and L are both adjustable, which will lead to different compression ratios.

With PCA compression methods described above, a 2D rectangular array Re(FC) or Im(FC) can be compressed off-line and stored in a highly compressed min look-up table. An arbitrary element in the original uncompressed 2D array can be decompressed and retrieved by Eq. (12) from the mini look-up table. In CGH calculation, the elements belonging to the corresponding radius line will be retrieved when the calculation is performed at a certain depth plane. The data belonging to the radius lines at other depths will remain compressed in the mini look-up table.

4. Numerical simulation experiment results

It is assumed that the 3D object points are located in the depth range of [0.2m 0.5m] from the hologram plane. The wavelength of illumination light is 632.8nm and the hologram pixel size is 6.5μm. The size of one hologram fringe pattern is 1024 × 1024 and the length of the radius line is set to be 768, which is close to the half-diagonal length of the pattern. The number of depth planes (or layers) within [0.2m 0.5m] is set to be 768. It shall be noted that the number of depth planes can be any arbitrary number, not necessarily equal to the length of radius line. Each pixel of hologram fringe pattern is represented by 8 bytes, 4 bytes for real part and 4 bytes for imaginary part.

For radial symmetric interpolation method, the size of 2D lookup table will be 768 × 768 (Q = 768, M = 768). In our proposed PCA based compression scheme, the partitioned block size is set to be $a=b=16$ ($L=256$). The number of principal components is set to different values ($W_{pca}=22$, $W_{pca}=17$, and $W_{pca}=12$). The memory storage size of LUT and compression ratios in these three cases compared to conventional radial symmetric interpolation method are calculated and shown in Table 1.

Table 1. Comparison of SSIM and memory storage size performance of different CGH calculation methods.

View Table | View all tables in this article

4.1 Reconstructed image quality at different depth planes

First, a numerical simulation experiment is conducted to evaluate the image quality of holographically reconstructed images when our prosed method is applied. The reconstructed image quality is measured by the Structural Similarity Index (SSIM) [11] with reference to the images reconstructed from holograms generated by direct analytic method. SSIM is a generally accepted image quality metric by assessing three visual characteristics including luminance, contrast and structure in the form of one single value between 0 and 1 (1 is best quality). A “pepper” test image of size 256 x 256, shown in Fig. 4(a) is placed at 0.2m, 0.3m, 0.4m and 0.5m respectively from the hologram plane. The CGH is calculated with direct analytic calculation base on Eqs. (1) and (2), radial symmetric interpolation method and proposed mini lookup table method with different $W_{pca}$ values respectively for four different distances. The average SSIM of the reconstructed images at four different distances are illustrated Table 1. As an example, the reconstructed images at the distance of 0.3m are illustrated in Fig. 4.

 

Fig. 4 (a) Original “pepper” image; (b-f) Reconstructed images from CGH calculated by (b) direct analytic method; (c) Radial symmetric interpolation method; (d) Proposed mini look-up table method with $W_{pca}=22$; (e) Proposed mini look-up table method with $W_{pca}=17$; (f) Proposed mini look-up table method with $W_{pca}=12$.

Download Full Size | PDF

It can be observed from Table 1 that the memory storage size of look-up table can be reduced from 4.5MB to 200-500 KB with proposed method, compared with conventional radial symmetric interpolation method. When the number of principal components $W_{pca}$ is decreased, a higher compression ratio can be achieved and more reconstructed image quality degradation is induced at the same time. However, it is evident that the holographically reconstructed image, shown in Fig. 4(f), with proposed mini look-up table method when $W_{pca}=12$ (compression ratio = 19.0413), is still very close to the one with conventional radial symmetric interpolation method, shown in Fig. 4(c). The quality degradation is trivial even when the memory storage size is highly compressed. With more error tolerance in the reconstructed image, the memory usage is likely to be further reduced to less than 200KB.

Here, we also evaluate the SSIM performance of reconstructed images when the lookup table is compressed by Discrete Cosine Transform (DCT) and Discrete Wavelet Transform (DWT) method. In these two methods, the compression is achieved by discarding the small value coefficients under a certain threshold after the transform. By adjusting the threshold, the compression ratio in each method can be set close ( ± 0.1) to the three ratios in Table 1 and their reconstructed image quality is compared with PCA method, shown in Table 2.

Table 2. Comparison of computation speed performance of different CGH calculation methods.

View Table | View all tables in this article

The results indicate that our proposed PCA method outperforms both DCT and DWT method since a higher SSIM can be achieved at the same compression ratio.

Another comparative study is conducted on the reconstructed image quality when the number of depth planes varies in the lookup table. The SSIM value versus compression ratio curve when the number of depth planes is 128, 256, 512 and 768 is plotted respectively in Fig. 5. The results indicate that the compression capability of our proposed method will be reduced when depth planes are sparser, which is due to the fact that different parts in the lookup table are less correlated when the intervals between depth planes are increased. But our proposed compression method can still give an acceptable result. For example, a SSIM of more than 0.9 at a compression ratio of 10 can still be maintained when the number of depth planes is recued to 128.

 

Fig. 5 SSIM value of reconstructed image versus compression ratio curve when the number of depth planes varies.

Download Full Size | PDF

4.2 Computation speed performance

Next, to evaluate the computation speed performance of proposed method in general situations, a 3D object is assumed to consist of N random points located within the depth range [0.2m 0.5m]. Their horizontal and vertical coordinates are distributed in a random manner as well. Then a CGH is generated from such a point cloud by direct analytic method, radial symmetric interpolation method and proposed mini lookup table method respectively. The total calculation time is summarized in Table 3.

Table 3. Comparison of computation speed performance of different CGH calculation methods.

View Table | View all tables in this article

It can be observed that the computation time of proposed mini look-up table method is almost the same as radial symmetric interpolation method for different number of objects points. This indicates that the extra computation cost in the look-up table decompression step is trivial and negligible in the entire CGH computation process. Both proposed method and original radial symmetric interpolation method are much faster than direct analytic method.

4.3 Reconstruction of 3D object model

Finally, we calculate holograms for real 3D object models. In the first experiment, a hologram is calculated from a simple 3D scene consisting of four characters “A”, “B”, “C” and “D”, placed at 0.2m, 0.22m, 0.24m and 0.26m respectively from the hologram plane. Both direct analytic method and our proposed mini lookup table method are employed for hologram generation. The generated hologram is numerically reconstructed at each focusing distance where only one character is in focus and the remaining three are defocused. The results are illustrated in Fig. 6. It can be observed that the reconstructed images, from hologram generated by our proposed method when the lookup table is highly compressed (compression ratio = 19.0413), can focus correctly. The performance of our proposed method is very close to direct analytic method.

 

Fig. 6 Numerical reconstruction results for a 3D scene containing four characters placed at different distances from a hologram generated by direct analytic method ((a)-(d)) and our proposed method (compression ratio = 19.0413) ((e)-(h)).

Download Full Size | PDF

In second experiment, a “Elephant” 3D model consisting of 24955 points is employed to calculate the hologram. The calculation time by direct analytic method, radial symmetric interpolation method and proposed mini lookup table method are 3338.2 seconds, 115 seconds and 115.1 seconds respectively. The hologram generated with proposed method (compression ratio = 19.0413) is numerically reconstructed at different distances and it can be observed that different body parts of the elephant are focused and defocused at different distances correspondingly (the focused part is labeled with red line), shown in Fig. 7.

 

Fig. 7 Numerically reconstructed images from a hologram generated by our proposed method (compression ratio = 19.0413) at different distances, showing the focus change.

Download Full Size | PDF

The SSIM difference between the reconstructed images at the same distance by direct analytic method and by our proposed method is 0.9165 (compression ratio = 19.0413), 0.9626 (compression ratio = 13.4737) and 0.9997 (compression ratio = 10.4253) respectively. The results reveal that our proposed method can maintain good image quality for 3D holographic reconstruction. In addition to numerical reconstruction, the generated hologram is optically reconstructed by a self-assembled phase-only spatial light modulator (SLM). Numerical and optical reconstruction by direct analytic method and our proposed method is illustrated in Fig. 8. Due to the SLM device limitation in our laboratory, the optically reconstructed image quality is degraded compared to numerical reconstructed result. However, it can be observed that the reconstructed images by direct analytic method and our proposed method are very close to each other using the same device.

 

Fig. 8 Numerically reconstructed 3D image: (a) by direct analytic method; (b) by our proposed method (compression ratio = 19.0413); Optically reconstructed 3D image: (c) by direct analytic method; (d) by our proposed method (compression ratio = 19.0413).

Download Full Size | PDF

The results in Section 3.1, 3.2 and 3.3 reveal that our proposed mini look-up table method can reduce the memory storage size of 2D look-up table in conventional radial symmetric interpolation method substantially by 10 to 20 times or even more. In the meanwhile, there is minor extra cost in computation time and reconstructed image quality. In practice, the calculation of CGH can be implemented with CPU and GPU on a personal computer or a supercomputer. Alternatively, it can be implemented by microprocessor chips such as FPGA [12,13]. The memory reduction of look-up table by our method is meaningful for both cases. In the former case, a single reduction of memory size at MBytes and KBytes level can be very trivial but it will be significant when a large number of holograms are calculated simultaneously and tiny memory reductions are accumulated to a huge effect. In the latter case, the memory size of microprocessor chips is usually very limited and our work on mini lookup table can directly lower the difficulty of system design.

5. Conclusion

In fast calculation of Computer Generated Hologram (CGH), look-up table methods can significantly enhance the computation efficiency. On the downside, the memory usage of look-up tables is tremendous. To address this problem, look-up tables can be compressed in different ways and radial symmetric interpolation is an effective method with high compression ratio. In this paper, the size of look-up table in radial symmetric interpolation method can still be further compressed to 1/10-1/20 or even less, by our proposed mini look-up table approach based on Principal Component Analysis (PCA). Since PCA compression can cause information loss and the decompression takes some computation time, we investigate the possible speed slowdown and reconstructed image quality degradation, compared to conventional radial symmetric interpolation method. The results reveal that our proposed method will only introduce negligible loss in computation efficiency and trivial quality degradation. As a consequence, the compression of look-up table is performed at almost no extra cost. Our proposed method can achieve the smallest look-up table size so far in the category of look-up table based fast CGH calculation methods.