1. Introduction

The vast potential of holography in true three-dimensional (3D) displays has attracted remarkable research attention [1,2]. Nevertheless, for many years, the calculation burden of computer-generated hologram (CGH) has hindered its application in real-time display. Many methods have been proposed to accelerate the computation [3,4]. Such methods include coherent ray tracing (CRT) [5], look-up tables (LUT) [6–10], and wave-front plane(WRP) [11,12] based on point-cloud distribution, fast Fourier transform (FFT) [13–15] and analytical frequency spectrum (AFS) [16–22] based on polygon-unit distribution, and angular-spectrum transfer [23–26] based on layer distribution. Some of these methods combine hardware, such as graphics processing units (GPUs) and field programmable gate arrays (FPGAs) for computing in parallel [27,28]. Moreover, some methods improve the quality of reconstruction, such as Gerchber–Saxton (GS) algorithm [29,30], further combine improved coding methods, such as double phase-amplitude coding [31,32], to display complex amplitude on a spatial light modulator (SLM). In recent reports, neural network based on deep learning has been employed to build consumer-grade holographic generation frameworks, which have significantly improved reconstruction quality and calculation speed [33,34].

Admittedly, polygon-based methods have their distinct potential in some special applications, such as relatively smooth 3D models (e.g., the human face) and simple 3D projection without the high requirement for details like head-up display (HUD) [35]. The amount of calculation in polygon-based methods could be reduced compared with that of the point cloud–based methods. Nevertheless, expressions calculated using polygon-based methods are not as simple and flexible as the point cloud–based methods. The conventional FFT-based polygon methods have shown the possibilities and correctness of using polygon division to generate hologram, but it is time-consuming because of the FFTs and interpolations required [13–15]. Ahrenberg et. al. [17], Pan et. al. [21] and Zhang et.al. [22] employed the polygon-based method in full analytical spectra expressions using affine transformation, which greatly saves much time for calculating spectra. Liu et. al. [19] developed the full analytical expressions in the spatial domain. To further improve polygon-based hologram, herein, we extract an idea of LUT used in the point-based method and implement it in the polygon-based method, which has been proved to improve calculation efficiency.

LUT methods are used to eliminate redundant computations for the point spread function. The diffractive field of a point with different z depth or x–, and y– coordinates comprises many identical operations, such as exponent, square-root, and several multiplications [8]. All repetitive operations are precalculated, and the data are stored in a memory or hard disk as LUTs, this is referred as offline calculation. Then, loading and look up data from LUTs calculate online diffractive field distribution at each hologram point; however, the LUT methods require a large memory size. Several methods based on LUT can reduce memory size from GigaByte-level [7] to KiloByte-level [8,10]. In addition, there are some repetitive operations in the full-analytical frequency spectrum expression of polygon-based methods [21,22]. The online calculation load must exhibit lightweight to eliminate redundant computations. Herein, we propose a novel analytical algorithm using LUT for polygon-based hologram called the analytical frequency spectrum-LUT method (AFS-LUT).

There are two primary steps to generate a polygon-based hologram: solving affine transformation matrix and computing spectra by analytical expression. Zhang’s method [22] employs a different idea from that of Pan’s method [21] to solve the affine matrix, bypassing the solution of the pseudo-inverse matrix. However, herein, we adopt Pan’s method to combine LUT due to its relative simplicity. Meanwhile, we improve the way of solving the affine transformation matrix to make it lighter. The most time-consuming part is computing analytical expressions, which are much more complicated than PSF used in point-cloud methods as they contain many multiplications, exponents, and even divisions.

The method proposed herein mainly pre-calculates AFSs and stores them in LUTs, thus eliminating the need for exponential calculations and reducing the number of multiplications in the online calculation. Experimental results show that the calculation speed is increased by almost twice compared to that of the full analytical calculation method [17,19,21,22]. However, AFSs have much redundant data, which would cause the need for large memory. Therefore, we performed a principal component analysis (PCA) on the LUTs, compressing much memory [36], which is employed in the point-based method [37]. Comparing the image quality between the proposed and existing methods using the peak signal-to-oise ratio (SNR) and structural similarity index measure (SSIM), we confirmed that the image obtained using the proposed method has the better image quality without visible artifact.

2. Proposed method

2.1 Solving an affine matrix simply

In Fig. 1(a), a triangle $\Gamma$ with vertexes $\left (x_{i},y_{i},z_{i}\right )$, $i=1,2,3$, and a primitive triangle $\Delta$ with vertexes $(x_{0i}, y_{0i}, z_{0i})$, $i=1,2,3$, show an affine relationship $\mathbf {T}: (x_{i},y_{i},z_{i})\mapsto (x_{0i}, y_{0i}, z_{0i})$ [21,22], which can be represented by

The equation above can be simplified as $\left [\mathbf {x},\mathbf {y},\mathbf {z}\right ]^\intercal = \mathbf {T} [\mathbf {x}_{0}, \mathbf {y}_{0}, \mathbf {z}_{0},1]^\intercal$, where $\mathbf {x}, \mathbf {y}, \mathbf {z}$ and $\mathbf {x}_{0}, \mathbf {y}_{0}, \mathbf {z}_{0}$ are the coordinate vectors of the triangle $\Gamma$ and the triangle $\Delta$, respectively; the superscript $^\intercal$ denotes the transpose operation and $\mathbf {T}=[T_{ij}]_{3\times 4}$. In Fig. 1(a), three vertexes of the primitive triangle $(x_{0i}, y_{0i}, z_{0i})$ are $(1, 0, 0), (1, 1, 0)$, and $(0, 0, 0)$, respectively, so defining the non-square matrix $[\mathbf {x}_{0}, \mathbf {y}_{0}, \mathbf {z}_{0},1]^\intercal =\mathbf {A}=\left [\begin {smallmatrix}1 & 1 & 0\\0 & 1 & 0\\0 & 0 & 0\\1 & 1 & 1\end {smallmatrix}\right ]$. Because $\textrm {rank}(\mathbf {A})=3$, $\mathbf {A}$ is a full column rank matrix, there must exist a pseudo-inverse matrix $\mathbf {A}^{\dagger }$ satisfying the criteria of $\mathbf {A^{\dagger }A=\mathbf {I_3}}$, where the superscript $^\dagger$ denotes the pseudo-inverse operation of matrix. Through the pseudo-inverse matrix solution method cited in Ref. [38] that $\mathbf {A}^{\dagger }=(\mathbf {A}^{\intercal }\mathbf {A})^{-1}\mathbf {A}^{\intercal }$, we can pre-compute that $\left ([\mathbf {x}_{0}, \mathbf {y}_{0}, \mathbf {z}_{0},1]^\intercal \right ) ^{\dagger } =\mathbf {A}^{\dagger }= \left [\begin {smallmatrix}1 & -1 & 0 & 0\\0 & 1 & 0 & 0\\-1 & 0 & 0 & 1\end {smallmatrix}\right ]$. Therefore, the affine transform matrix can be solved by $\mathbf {T} =[\mathbf {x},\mathbf {y},\mathbf {z}]^\intercal \left ([\mathbf {x}_{0}, \mathbf {y}_{0}, \mathbf {z}_{0},1]^\intercal \right ) ^{\dagger }$, which is expanded as

 

Fig. 1. (a) Affine transformation of an arbitrary triangle $\Gamma$ to a primitive triangle $\Delta$. (b) Physical scene of a plane wave diffraction passing through a titled triangle aperture $\Gamma$.

Download Full Size | PDF

The triangle mesh is assumed to be a transmission aperture with the amplitude transmittance function $a(x,y,z)$, where $(x,y,z)$ are the coordinate points within the triangle $\Gamma$, as shown in Fig. 1(a). A planar wave illuminating to the triangle $\Gamma$ can be expressed as $\exp [j2\pi (x\cos \alpha +y\cos \beta )/\lambda )]$, where $\lambda$ is wavelength, and $j=\sqrt {-1}$. $\alpha$ and $\beta$ are the angles with regard to the $x-$ and $y-$ axis, respectively, as shown in Fig. 1(b). After the plane wave passes through the aperture, the complex amplitude immediately behind the triangle $\Gamma$ is represented by $E_0= \exp [j2\pi (x\cos \alpha +y\cos \beta )/\lambda )]\cdot a(x,y,z)\exp (j2\pi z\cos \gamma /\lambda )$, where $\cos \gamma =\sqrt {1-(\cos ^2\alpha +\cos ^2\beta )}$ and $\gamma$ is the angle with regard to $z-$axis. In this study, the planar wave illuminating to the triangle $\Gamma$ is assumed to propagate along the positive direction of $z-$axis toward the hologram plane(located at the $z=0$), which indicates $\cos \alpha =\cos \beta =0$, as shown by the black dashed line in Fig. 1(b). Thereby, $E_0=a(x,y,z)\cdot \exp (j2\pi z/\lambda )$, where $z=T_{31}x_{0}+T_{32}y_{0}+T_{33}z_{0}+T_{34}$ according to Eq. (1), $(x_0,y_0,z_0)$ represent the coordinate points within the triangle $\Delta$. Using angular spectrum $H=\exp (-j2\pi w_0z)$, where $w_0=(1/\lambda ^2-u_0^2-v_0^2)^{1/2}$, as transfer function to perform diffraction calculation, the frequency spectrum distribution of the triangle $\Gamma$ in the hologram plane can be expressed as

It is worth noting that, different from the other methods, e.g., Ref. [14,19,21] and [22], the Eq. (5) indicates that the zero-frequency of the new frequency coordinates $(u', v')$ lies in the center of the sampling frequency coordinates $(u_0,v_0)$ without artificially subtracting a frequency offset whose value is $\frac {T_{31}}{\lambda }$ and $\frac {T_{32}}{\lambda }$ in Ref. [21]. This is because if the light waves passing through the triangle $E_0$ is not defined, it will be defaulted to propagate along the normal direction of triangle, as shown by the red dashed line in Fig. 1(b). The triangular spectrum therefore produces the frequency offset due to the triangle $\Gamma$ being tilted in general, and thus, it may exceed the frequency sampling limit of the hologram plane, i.e., $(u_0)_{max}$, as shown by the red spectrum in Fig. 1(b). Here we can learn that subtracting the factors of $\frac {T_{31}}{\lambda }$ and $\frac {T_{32}}{\lambda }$ is physically represented as light propagating along the $z-$axis instead of its normal axis.

2.2 Solving analytical frequency spectrum by LUT

Because $w_0$ can be expanded by the Taylor series as $w_0\approx 1/\lambda -(\lambda /2)(u^2_0+v^2_0)$, Eq. (5) can be further expressed as

In addition, solving the integral for Eq. (4), we can obtain AFS as follows:

In the AFS expression, the main workload is the last term (when $u'$ and $v'$ and $u'+v'$ are not equal to $0$). Therefore, we will focus on the last term to implement the LUT method. For the last term,

2.2.1 For the $\exp (-j2\pi u')$ and $\exp [-j2\pi (u'+v')]$ term

Firstly, for the $\exp (-j2\pi u')$ terms, $u'$ can be deduced from Eq. (6) in a perfect square expression as

Thus, we define the following four tables for storing pre-calculated data,

For the $\exp [-j2\pi (u'+v')]$ term in Eq. (8), $(u'+v')$ can similarly be expanded as

From Eq. (2), $T_{11}+T_{12}=x_3-x_1$, $T_{21}+T_{22}=y_3-y_1$ and $T_{31}+T_{32}=z_3-z_1$, which are in the same range as $T_{11}$, $T_{21}$ and $T_{31}$. Therefore, the three of them are used as input arguments ($\delta x$, $\delta y$ and $\delta z$) for TU1, TV1, TU2, and TV2 tables, respectively. Then, the $\exp [-j2\pi (u'+v')]$ term can be solved by

Eqs. (11) and (13) indicate that both $\exp (-j2\pi u')$ and $\exp [-j2\pi (u'+v')]$ can be obtained from the same set of tables, TU1, TV1, TU2, and TV2. Figure 2 illustrates the process for computing F1 online using TU1 and TV1 (computing F1 if $T_{31}=0$ using TU2 and TV2 is straightforward and not shown in Fig. 2). Then, we can obtain the final $[$F1$]_{n\times m}$ by multiplying with the constant factor $C_1$. $[$F1$']_{n\times m}$ can be computed following the same procedure.

 

Fig. 2. Schematic diagram of the holographic plane spectrum obtained using TV1 and TU1.

Download Full Size | PDF

2.2.2 $\exp (-j2\pi w')$ term

Next, we analyze the $\exp (-j2\pi w')$ term in Eq. (8). Based on Eq. (5), similar to Eq. (9), $w'$ in $\exp (-j2\pi w')$ can be expanded as

Similarly, we define two LUTs as follow,

However, considering the practical situation of holographic calculations, the distance from the 3D object to the hologram plane is generally much greater than the hologram size in the xy-direction. Thus, $\frac {T_{14}}{T_{34}}$=$\frac {x_1}{z_1}$, is a finite range of values because the denominator $z_1$ is much greater than the numerator $x_1$, and $\frac {T_{14}}{T_{34}} \frac {1}{\lambda }$ will be scaled up because $\frac {1}{\lambda }$ is a value in thousands of millimeters, thus reducing the requirement of $\frac {T_{14}}{T_{34}} \frac {1}{\lambda }$ for precision. As a result, $\frac {T_{14}}{T_{34}} \frac {1}{\lambda }$ is treated as a single variable $s_1$, where $s_1\in \left [ \frac {(T_{14})_{min}}{(T_{34})_{max}}\frac {1}{\lambda }, \frac {(T_{14})_{max}}{(T_{34})_{min}}\frac {1}{\lambda }\right ]$. Similarly, $\frac {T_{24}}{T_{34}} \frac {1}{\lambda }$ is named $s_2$. Further analysis shows that the vector $[u_0-s_1]_{1\times i}$ in table TU0 and the vector $[v_0-s_2]_{k\times 1}$ in table TV0 can be discretely sampled with the period defined as $S_p$. $i$ and $k$ are the sampling amounts, which are expressed as

The sampling period $S_p$ is also regarded as the precision of $s_1$ and $s_2$. In general, $S_p$ is set to the order of magnitude of $\frac {T_{14}}{T_{34}} \frac {1}{\lambda }$ divided by $10^3$, which has been experimentally confirmed to produce no significant error. As an example, herein, the 3D object size is set to be comparable to SLM, i.e., $x\in [−8,8]mm$, and the position on the $z$-axis is set between 450 and 500$mm$. Thus, $(s_1)_{min}\approx$ −30, $(s_1)_{max}\approx$ 33, its order of magnitude is $10^{1}$, then ${S_p}=10^1/10^3=0.01$.

In Eq. (16), we define $l$ as the sampling number of $T_{34}$. Therefore, the 2D table $[$TU0$]_{i\times l}$ requires $(i \times l\times 2\times 8)$ bytes of memory, and $[$TV0$]_{k\times l}$ requires $(k\times l\times 2\times 8)$ bytes of memory, which are much smaller than the memory of TU1 and TV1 because both are 2D tables.

When computing Eq. (16) online, as shown in Fig. 3, $s_1$ and $s_2$ are first obtained according to $T_{14}$ , $T_{24}$ and $T_{34}$ , respectively. Next, the starting position and step length($S_l=s_1/S_p$ or $s_2/S_p$) in LU0 and LV0 are computed. Then, considering the data from TV0 and TU0 based on $S$, the matrix $[$F0$]_{m\times n}$ with the same size as the hologram is obtained by multiplying the column vector $[$TV0$]_{m\times 1}$ with the row vector $[$TU0$]_{1\times n}$ and the constant $C_{0}$.

 

Fig. 3. Schematic diagram of the holographic plane spectrum obtained using TV0 and TU0.

Download Full Size | PDF

2.3 PCA to compress the memory of tables

As mentioned in Section 2.2, LUTs require large memory, especially TU1 and TV1. Herein, PCA is used to further compress the LUT size. PCA extracts the principal components from the data set and removes data redundancy, maintaining almost the features of the original data [36]. The outline of PCA and how to apply it to the proposed method are as follows:

Take a data matrix $[\mathbf {X}]_{n\times p}$ as an example that has already shifted the sample mean of each column to zero, and $n$ and $p$ are the numbers of samples and variables, respectively. First, we solve the covariance matrix between different variables in matrix $\mathbf {X}$, named $[\mathbf {K}]_{p\times p}$, given by

We define $\mathbf {V}_{i}$ as the $i^{\,th}$ column eigenvector with corresponding eigenvalue $D_{ii}$, where $D_{ii}$ is the diagonal element of the diagonal matrix $\mathbf {D}$. Next, the corresponding eigenvectors are reordered in the order of eigenvalues $D_{ii}$ from largest to smallest. Set $W_i=\frac {D_{ii}}{\sum _{i=1}^{p}D_{ii}}$, it means that the $i^{\,th}$ eigenvector accounts for the weight $W_i$ of the original matrix $\mathbf {X}$. There must exist a $k$ that satisfies that $\sum _{i=1}^{k}W_i \approx 1$, which $k\ll p$. This implies that only the first $k$ eigenvectors in $\mathbf {V}$ approximate the features of $\mathbf {X}$. Finally, let $[\mathbf {C}]_{p\times k}=[\mathbf {V}_{1}, \mathbf {V}_{2},\dots, \mathbf {V}_{k}]$ as the coefficient of the principal component of $\mathbf {X}$, which the principal component is defined as $\mathbf {P}$ obtained by

As shown in Eq. (20), since $k\ll p$, the principal component matrix $\mathbf {P}$ and the coefficient matrix $\mathbf {C}$ require only small memory. Thus $\mathbf {P}$ and $\mathbf {C}$ are stored as LUTs instead of the original matrix $\mathbf {X}$. When calculating online, the original data $\mathbf {X}$ can be approximately retrieved by

Notable, $\mathbf {C}$ is an orthogonal matrix, so $\mathbf {C}^{-1}=\mathbf {C}^{^\intercal }$. Finally, we can compress all LUTs used in Eqs. (11) and (16) by the PCA.

3. Implementation of the method

3.1 Algorithm

The proposed method was implemented in the MATLAB environment. This sub-section presents the pseudo–code of the proposed algorithm, AFS-LUT, based on the parameters defined in Table 1.

Table 1. Definition of parameter symbols in the algorithm.

View Table | View all tables in this article

In Table 1, the precision of $dx$, $dy$, and $dz$ are determined by the physical size of the object. In general, the object size is treated the same as the hologram size. In the $z$–direction, the polygon-based method does not require as many layers as in the point-based method. Here, the sampling interval $dz$ is set as $0.04$mm at a depth range of $20$mm, which is sufficient to describe the object in details without changing its shape. The pseudo–codes for offline and online calculations are listed in Table 2.

Table 2. Pseudo MATLAB-like codes of the proposed algorithm, AFS-LUT.

View Table | View all tables in this article

3.2 Results and performance evaluation

To verify the proposed method, here, we present the calculated results for 3D objects with several triangles. The program developed using MATLAB 2021a was implemented on a personal computer configured with CPU: AMD Ryzen 5-3600 @3.59GHz RAM: 32GB. Parallel computing in GPU(Nvidia GeForce RTX 3070) was implemented by calling the CUDA functions in MATLAB.

Figure 4 shows two 3D models composed of different numbers triangles. Considering occlusion, the actual polygons involved in the calculation are 3084 triangles for a teapot and 27,619 triangles for a human head. The parameters used for the calculation are listed in Table 1. To evaluate the performance of the proposed method, Table 3 lists the evaluation criteria in comparison with the AFS method [19,21,22] for the calculated result of the teapot. Compared to the AFS method, the proposed AFS-LUT method achieves higher computational efficiency mainly because the exponential operation is eliminated and extra multiplication operations do not increase the load. On the other hand, three types of data tables require a large amount of memory because all possible values of the exponential term are stored. Here, the PCA technique is used to compress the redundant data about four times without much loss in the image quality as revealed by PSNR and SSIM. Based on the parameters in Table 1, Table 4 lists the actual memory required for each LUTs and compares them with that required by PCA. LU2 and LV2 are not compressed by PCA due to very small memory.

 

Fig. 4. 3D models composed of triangles and the corresponding numerical reconstructed image. (a) 3084 triangles of the teapot, (b) the reconstructed image of teapot at 520mm, (c) 27, 619 triangles of the human’s head, (d) the reconstructed image of human’s head at 500mm.

Download Full Size | PDF

Table 3. Performance comparison among methods for calculating "teapot" model.

View Table | View all tables in this article

Table 4. Memory of each table compared with PCA compression(Unit in MB)."P" and "C" denote the principal component and its coefficient matrices.

View Table | View all tables in this article

Figure 5(a) depicts the performance difference of each method in calculating different numbers of triangles. Parallel calculation implemented using "gpu" functions in MATLAB greatly improves the calculating speed. Based on this, the proposed method, AFS-LUT, accelerates further about twice compared with the AFS method in GPU [17,21,22]. The specific time consumed is listed at the bottom of Fig. 5(b). Notably, as the number of triangles increases, the AFS-LUT with PCA method performs faster than the AFS method, because the time spent to restore the LUTs is constant. The more the number of triangles, the lower the proportion of time to restore the LUTs. In addition, the required memory does not change as the number of polygons increases, as shown in the right axes of Fig. 5(a). At the top of Fig. 5(b), the highly consistent PSNR and SSIM values indicate that no artifacts are produced even with PCA.

 

Fig. 5. Performance comparison among methods for calculating five 3D objects. (a) Elapsed time by different methods (left axis) and memory required with and without PCA (right axis). (b) Top: PSNR (left axis) and SSIM (right axis) with and without PCA compared to AFS; Bottom: Specific time elapsed corresponding to (a).

Download Full Size | PDF

3.3 Experimental reconstruction

The holograms calculated by the proposed method are optically reconstructed by a reflective liquid crystal on silicon(LCoS) of Holoeye GAEA-2 phase-only SLM with 4160$\times$2464 pixels. The maximum angle of diffraction by the SLM is $(\theta )_{max}=\arcsin {(\frac {\lambda }{2p}})\approx 4^{\circ }$, where $p=3.74\mu m$ is the pixel pitch of the SLM and $\lambda =532nm$. By loading a virtual off-axis grating as the reference light, the complex amplitude hologram, $H_{\mu }$, can be represented by

The experimental optics framework is shown in Fig. 6(a). Due to the property of the SLM, the phase of only horizontally polarized light (p light) is modulated by the liquid crystal molecules of the SLM, and the vertically polarized light (s light) will be the background noise for the reflected light. Two polarizers are used to weaken the background light of SLM, thus increasing the contrast. The "Lens" in the figure is designed to produce a frequency spectrum in the focal plane, and the aperture in the figure is used to filter out the zero-order light spectra and to pass through the reconstructed image. The reconstructed images are captured by a camera (Sony $\alpha$7R-II) with a zoom lens (vari-focal: 24–240 mm), as shown in Figs. 6(b) and (c).

 

Fig. 6. (a) Schematic diagram of optical experiment reconstructions. (b) and (c) reconstruction images in the optical experiment.

Download Full Size | PDF

4. Conclusion

In this study, we investigated algorithms for generating holograms using the polygon-based method. We first propose a new method called AFS-LUT, which employs LUTs. Then, to solve the problem of the large memory required in the LUTs, PCA was used to compress the memory from GB- to MB-level. By computing different numbers of polygons and comparing the calculation efficiencies, the AFS-LUT method was confirmed to improve the calculation speed by twice. Considering PSNR and SSIM, we found that the image reconstructed by the proposed method has almost no artifacts. Further, we implemented parallel calculations in GPU by evoking the build-in functions in MATLAB, which could increase the calculation speed several tens of times compared with the calculations in the CPU. The proposed method was verified by reconstructing in simulation and optical experiments.