Recording and recreating three-dimensional (3D) scenes can be done using holography as a useful tool. The data generation of digital holograms known as computer-generated holograms (CGHs) involves simulating the light waves in the process of optical holographic recording, which includes modeling and simulating phenomena such as reflection, diffraction, propagation, and interference of light in computers. CGHs can be transmitted to the spatial light modulator for display. The three commonly used methods for generating CGH are the point-based method^[1], the layer-based method^[2–4], and the polygon-based method^[5–7]. Due to their computational simplicity, point-based methods have been rapidly developed. Point-based CGHs become more efficient with the use of lookup table methods^[8], wave-front recording plane methods^[9], parallel computation in GPUs^[10], and field-programmable gate arrays^[11]. In layer-based methods, Yao et al.^[12] recently proposed a method of generating holograms using adaptive algorithms to reduce the depth error of reconstructed images. In polygon-based methods, Fu et al.^[7] proposed a complete analysis method using self-similarity segmentation along with a parallelogram approximation. They have also used texture mapping to render objects, thereby reducing a large amount of computational workload in high-precision texture rendering.

Polygon-based computer-generated holography has many advantages including compatibility with conventional computer graphics techniques and efficient representation of large 3D scenes. One of the polygon-based computer-generated holography methods is the traditional FFT-based method, which requires two FFTs and interpolation to obtain the spectrum of the polygon on the hologram plane^[13,14]. However, interpolation reduces the calculation speed of hologram generation. An alternative method within polygon-based computer-generated holography is the fully analytical approach, which computes the spectrum of the arbitrary triangular mesh using affine transform^[15–17]. There are multiple challenges associated with implementing a 3D display using computer-generated holography. Enhancing the fidelity of the reconstructed images presents a significant hurdle in this endeavor. When reconstructing a 3D scene, it is essential to account for the amplitude of the triangles, considering factors such as the illumination direction and the triangle’s normal vector^[16,17]. Employing uniform amplitude distribution within each triangle leads to what is known as flat shading. However, the realism of the reconstruction is hindered by the visibility of mesh boundaries caused by flat shading. Park et al. used a continuous shading method to provide continuously changing amplitudes within each triangle. However, the exaggeration of polygon edges still affects the fidelity of 3D reconstruction^[6]. On this basis, Yeom et al. proposed the use of discrete convolution to achieve the addition of object reflection highlights. However, due to the large computational complexity of discrete convolution, computational efficiency is reduced^[18]. Wang et al. proposed using the Blinn Phong reflection model, replacing diffuse and specular reflections with an approximately exponential function and implementing the addition of highlights in a fully analytical algorithm, but there are still cases of discontinuous polygonal edges^[19].

To solve the problem of edge discontinuity in polygon reconstruction for 3D objects, we propose a Phong shading approximation method. It involves dividing the (primitive) unit triangle into sub-triangles by segmentation in the fully analytical approach to determine the amplitude of each point within the sub-triangles, thereby achieving an approximation of Phong shading. The proposed method effectively eliminates bright line artifacts along the edges of polygonal meshes on object surfaces and ultimately produces superior surface approximations. The proposed method rigorously uses the Phong reflection model^[20] and has been validated through simulations and optical reconstructions of multiple 3D objects.

In Sec. 2, we summarize the theory of polygon-based method based on affine transformation. In Sec. 3, we discuss continuous shading. In Sec. 4, we derive a complete analytical form of Phong shading approximation using self-similarity segmentation and linear interpolation. In Sec. 5, we demonstrate numerical and optical reconstructions to verify the effectiveness of our proposed technique. Concluding remarks will be provided in Sec. 6.

In this section, we briefly summarize the important results of the fully analytic triangular mesh algorithm used for computer-generated holograms^{[16,17,21–23]}. In the global coordinate system $O\text{-}xyz$ as shown in Fig. 1, an $i$th arbitrarily tilted polygon of the 3D object is shown. Its spectrum on the hologram plane ($z=0$) is given by^[24–26]$$G_i(f_x,f_y)=\iint a_i[x,y;z(x,y)]\exp \{-j2\pi [f_{x}x+f_{y}y+f_{z}z-\frac{z(x,y)}{\lambda }]\}\mathrm{d}x\mathrm{d}y,$$(1)where $\lambda$ is the wavelength of light. $f_x$, $f_y$, and $f_z$ are the spatial frequencies along the $x$, $y$, and $z$ directions with $f_z=\sqrt{\frac{1}{{\lambda }^2}-f_x^2-f_y^2}$. $z(x,y)$ is the distance between a point on the tilted plane and the hologram plane, and $a_i$ is the amplitude of the $i$th arbitrarily tilted triangle given by $$a_i[x,y;z(x,y)]=\{\begin{array}{cc}1,& \text{inside triangle}\\ 0,& \text{outside triangle}\end{array}.$$(2)

The brightness factor for shading rendering can be determined by the amplitude of the triangle. In traditional methods, it is assumed that the amplitude of each triangle is uniform. Therefore, the 3D scene is represented by what is known as flat shading, which lacks realism upon hologram reconstruction.

In the fully analytic triangular mesh algorithm, the spectrum of the tilted triangle on the hologram plane is expressed in terms of the spectrum of a primitive triangle (unit right triangle) through affine transformation. If the affine coordinates of the three vertices of the primitive triangle on the $x_r{y}_r$ plane are (0,0), (1,0), and (0,1), and the coordinates of the vertex of the $i$th arbitrary triangle on the global coordinate system $O\text{-}xyz$ are $(x_{i1},y_{i1},z_{i1}),(x_{i2},y_{i2},z_{i2})$, and $(x_{i3},y_{i3},z_{i3})$, we have the following relation: $$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{c}(x_{i2}-x_{i1})(x_{i3}-x_{i1})x_{i1}\\ (x_{i2}-y_{i1})(y_{i3}-y_{i1})y_{i1}\\ (z_{i2}-z_{i1})(z_{i3}-z_{i1})z_{i1}\end{array}\right]\left[\begin{array}{c}{x}_r\\ y_r\\ 1\end{array}\right].$$

With this relation, Eq. (1) can be finally expressed as $$G_i(f_x,f_y)=\mathcal{J}·p·\iint a_r(x_r,y_r)\exp (-j2\pi f_{xi}^{\prime }{x}_r)\exp (-j2\pi f_{yi}^{\prime }{y}_r)\mathrm{d}{x}_r\mathrm{d}{y}_r.$$(3)

In Eq. (3), $a_r$ is the amplitude of the primitive triangle, $$\begin{gathered} a_r(x_r,y_r)=\{\begin{array}{cc}1,& \text{inside primitive triangle}\\ 0,& \text{outside primitive triangle}\end{array}, \\ \mathcal{J}=\left|\begin{array}{c}\frac{\partial x_r}{\partial x}\frac{\partial x_r}{\partial y}\\ \frac{\partial y_r}{\partial x}\frac{\partial y_r}{\partial y}\end{array}\right|=|x_{i2}{y}_{i3}+x_{i1}{y}_{i1}+x_{i1}{y}_{i2}-x_{i3}{y}_{i2}|, \\ {f}_{xi}^{\prime }=(x_{i2}-x_{i1})f_x+(y_{i2}-y_{i1})f_y+(z_{i2}-z_{i1})f_z-\frac{z_{i2}-z_{i1}}{\lambda }, \\ {f}_{yi}^{\prime }=(x_{i3}-x_{i1})f_x+(y_{i3}-y_{i1})f_y+(z_{i3}-z_{i1})f_z-\frac{z_{i3}-z_{i1}}{\lambda }, \\ p=\exp [-j2\pi (x_{i1}{f}_x+y_{i1}{f}_y+z_{i1}{f}_z-z_{i1}/\lambda )]. \end{gathered}$$

Now, the primitive triangle’s analytical spectrum $G_r(f_x^{\prime },f_y^{\prime })$ is given by $$\begin{gathered} G_r(f_x^{\prime },f_y^{\prime })={\iint }_{.}{a}_r(x_r,y_r)\exp [-j2\pi (x_r{f}_x^{\prime }+y_r{f}_y^{\prime })]\mathrm{d}{x}_r\mathrm{d}{y}_r \\ ={\int }_0^1{\int }_0^{-x_r+1}1·\exp [-j2\pi (x_r{f}_x^{\prime }+y_r{f}_y^{\prime })]\mathrm{d}{x}_r\mathrm{d}{y}_r=\frac{\exp (-j2\pi f_x^{\prime })-1}{4{\pi }^2{f}_x^{\prime }{f}_y^{\prime }}-\frac{\exp (-j2\pi f_x^{\prime })-\exp (-j2\pi f_y^{\prime })}{4{\pi }^2(f_x^{\prime }-f_y^{\prime })f_y^{\prime }}. \end{gathered}$$(4)

Comparing Eqs. (3) and (4), we have $$G_i(f_x,f_y)=\mathcal{J}·p·G_r(f_{xi}^{\prime },f_{yi}^{\prime }).$$(5)

In the polygon-based method, the total light field distribution $U_O(x,y)$ on the hologram plane is the superposition of the light fields emitted by all polygons (triangles). The total optical field distribution $U_O(x,y)$ can be obtained by inverse Fourier transform of the total spectrum obtained from all the arbitrary triangles $G_i(f_x,f_y)$: $$U_O(x,y)={\mathcal{F}}^{-1}\{\sum _{i=1}^n[\mathcal{J}·p·G_r(f_{xi}^{\prime },f_{yi}^{\prime })]\},$$(6)where ${\mathcal{F}}^{-1}$ represents the inverse Fourier transform and $n$ denotes the number of triangles of the 3D object.

The Chevreul illusion^[27] manifests as an optical effect where, between two adjacent triangles, if one triangle has a stronger contrast, it causes our visual system to perceive boundaries between these adjacent triangles. To address the Chevreul illusion^[27], Park et al. proposed a continuous shading model that interpolates the amplitudes within the primitive triangle based on the amplitudes of the three vertices of the primitive triangle. Therefore, Park et al.’s method is essentially Gouraud shading model^[6].

The coordinates of the three vertices of the primitive triangle on the $x_r{y}_r$ plane are ${\mathrm{V}}_1(0,0),{\mathrm{V}}_2(0,1)$, and ${\mathrm{V}}_3(1,0)$. ${\widehat{\mathit{n}}}_{v,1}$, ${\widehat{\mathit{n}}}_{v,2}$, and ${\widehat{\mathit{n}}}_{v,3}$ are the known three vertex unit normal vectors of the arbitrary triangle in $xyz$-coordinates, respectively, as shown in Fig. 2. Park et al. used the Phong reflection model to determine the amplitude $a_{\nu ,1}$, $a_{\nu ,2}$, and $a_{\nu ,3}$ of the three vertices of the primitive triangle according to $$a_{\nu ,i}=a_{\nu ,i,o}(k_a+k_d{\widehat{\mathit{n}}}_{v,i}·\widehat{\mathit{L}}),i=1,2,3,$$(7)where $a_{\nu ,i,o}$ and ${\widehat{\mathit{n}}}_{v,i}$ represent the reflectance and normal vector, respectively, and $\widehat{\mathit{L}}$ is the illumination vector. The ratio of ambient to diffusive reflection is dictated by the two coefficients, $k_a$ and $k_a$. Utilizing the amplitudes of $a_{\nu ,1}$, $a_{\nu ,2}$, and $a_{\nu ,3}$, the amplitude $g_r(x_r,y_r)$ of pixels $(x_r,y_r)$ within the primitive triangle is linearly interpolated as $$g_r(x_r,y_r)=\{\begin{array}{cc}(a_{\nu ,3}-a_{\nu ,1})x_r+(a_{\nu ,2}-a_{\nu ,1})y_r+a_{\nu ,1},& \text{inside triangle}\\ 0,& \text{outside triangle}\end{array}.$$(8)

Hence, $a_r(x_r,y_r)$ within the primitive triangle in Eq. (3) can now be given by $$a_r(x_r,y_r)=g_r(x_r,y_r).$$(9)

The angular spectrum of the primitive triangle $G_r(f_x^{\prime },f_y^{\prime })$ can be obtained as $$G_r(f_x^{\prime },f_y^{\prime })=(a_{\nu ,3}-a_{\nu ,1})D_1+(a_{\nu ,2}-a_{\nu ,1})D_2+a_{\nu ,1}{D}_3,$$(10)where $$\begin{gathered} D_1(f_x^{\prime },f_y^{\prime })={\int }_0^1{\int }_0^{1-x_r}{x}_r\exp [-j2\pi (x_r{f}_x^{\prime }+y_r{f}_y^{\prime })]\mathrm{d}{x}_r\mathrm{d}{y}_r \\ =\frac{1}{8{\pi }^3{f}_y^{\prime }}j\{\frac{e^{-2j\pi f_x^{\prime }}(-2j\pi f_x^{\prime }+e^{2j\pi f_x^{\prime }}-1)}{f_x{\prime }^2}+\frac{e^{-2j\pi (f_x^{\prime }+f_y^{\prime })}\{e^{2j\pi f_y^{\prime }}[1+2j\pi (f_x^{\prime }-f_y^{\prime })]-e^{2j\pi f_x^{\prime }}\}}{{(f_x^{\prime }-f_y^{\prime })}^2}\}, \end{gathered}$$(11)$$\begin{gathered} D_2(f_x^{\prime },f_y^{\prime })={\int }_0^1{\int }_0^{1-x_r}{y}_r\exp [-j2\pi (x_r{f}_x^{\prime }+y_r{f}_y^{\prime })]\mathrm{d}{x}_r\mathrm{d}{y}_r \\ =\frac{e^{-2j\pi (f_x^{\prime }+2f_y^{\prime })}[j{e}^{2j\pi (f_x^{\prime }+2f_y^{\prime })}{(f_x^{\prime }-f_y^{\prime })}^2]}{8{\pi }^3{f}_x^{\prime }{f}_y{\prime }^2(f_x^{\prime }-f_y^{\prime }{)}^2} \\ +\frac{f_x^{\prime }{e}^{2j\pi (f_x^{\prime }+f_y^{\prime })}[f_x^{\prime }(2\pi f_y^{\prime }-j)-2f_y^{\prime }(\pi f_y^{\prime }-j)]-j{e}^{4j\pi f_y^{\prime }}{f}_y^{\prime 2}}{8{\pi }^3{f}_x^{\prime }{f}_y^{\prime 2}(f_x^{\prime }-f_y^{\prime }{)}^2}, \end{gathered}$$(12)$$\begin{gathered} D_3(f_x^{\prime },f_y^{\prime })={\int }_0^1{\int }_0^{1-x_r}1·\exp [-j2\pi (x_r{f}_x^{\prime }+y_r{f}_y^{\prime })]\mathrm{d}{x}_r\mathrm{d}{y}_r \\ =\frac{\exp (-j2\pi f_x^{\prime })-1}{4{\pi }^2{f}_x^{\prime }{f}_y^{\prime }}-\frac{\exp (-j2\pi f_x^{\prime })-\exp (-j2\pi f_y^{\prime })}{4{\pi }^2(f_x^{\prime }-f_y^{\prime })f_y^{\prime }}. \end{gathered}$$(13)

By incorporating Eq. (10) into Eq. (6), Park et al. obtained a fully analytical spectrum with continuous shading. However, after numerical reconstruction, we have found significant discontinuities between the edges of the object surface polygons and bright line artifacts appear on the polygonal edges of the object surface reconstruction. To address this issue, we propose a Phong shading approximation method to make the surface of the object smoother, which will be detailed in Sec. 5 along with Park et al.’s results for comparison.

In the proposed method, we use a simple Phong reflection model that includes ambient reflection and diffuse reflection, which provides the amplitude $a$ of each surface point^[20]: $$a=k_a+k_d(\widehat{\mathit{n}}·\widehat{\mathit{L}}),$$(14)where $k_a$ and $k_d$ are the ambient and diffuse reflection coefficients, respectively, $\widehat{\mathit{n}}$ is the unit normal vector, and $\widehat{\mathit{L}}$ is the unit vector pointing toward the light source. This equation is reminiscent of Eq. (7).

In Fig. 3, we denote $(m,n)$ as an index value along the two unit lengths of the primitive triangle. The index value describes the position of the sub-triangle. Figure 3 also illustrates the coordinates of the vertices of $(m,n)$ sub-triangles (forward in green and reverse in orange). Each side of the unit triangle is divided into segments of length $1/M$, and $(m,n)$ is used as the index value along the two unit lengths. Therefore, we can easily obtain the vertex coordinates of $(m,n)$ sub-triangles (green forward and orange reverse, as shown in Fig. 3. ${\widehat{\mathit{n}}}_{\nu ,1}^{\prime }$, ${\widehat{\mathit{n}}}_{\nu ,2}^{\prime }$, and ${\widehat{\mathit{n}}}_{\nu ,3}^{\prime }$ are the vertex normals of the primitive triangle, as shown in Fig. 3. In our work, similar to Eq. (8), we use the vertex coordinates of the sub-triangle shown in Fig. 3 to obtain the vertex normals of the sub-triangle by linearly interpolating the vertex normals of the primitive triangle. The vertex normals ${\mathit{n}}_{\mathrm{F},1}$, ${\mathit{n}}_{\mathrm{F},2}$, and ${\mathit{n}}_{\mathrm{F},3}$ of the forward triangle $\mathrm{F}(m,n)$ and the vertex normals ${\mathit{n}}_{\mathrm{R},1}$, ${\mathit{n}}_{\mathrm{R},2}$, and ${\mathit{n}}_{\mathrm{R},3}$ of the reverse triangle $\mathrm{R}(m,n)$ are obtained by the following equations: $$\begin{gathered} {\mathit{n}}_{\mathrm{F},1}(m,n)=({\widehat{\mathit{n}}}_{\nu ,3}^{\mathrm{\prime }}-{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }})·\frac{m}{M}+({\widehat{\mathit{n}}}_{\nu ,2}^{\mathrm{\prime }}-{\widehat{\mathit{n}}}_{\nu ,1}^1)·\frac{\mathit{n}}{M}+{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }}, \\ {\mathit{n}}_{\mathrm{R},1}(m,n)=({\widehat{\mathit{n}}}_{\nu ,3}^{\mathrm{\prime }}-{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }})·\frac{m+1}{M}+({\widehat{\mathit{n}}}_{\nu ,2}^{\mathrm{\prime }}-{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }})·\frac{\mathit{n}+1}{M}+{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }}, \\ {\mathit{n}}_{\mathrm{F},2}(m,n)={\mathit{n}}_{\mathrm{R},2}(m,n)=({\widehat{\mathit{n}}}_{\nu ,3}^{\mathrm{\prime }}-{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }})·\frac{m}{M}+({\widehat{\mathit{n}}}_{\nu ,2}^{\mathrm{\prime }}-{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }})·\frac{\mathit{n}+1}{M}+{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }}, \\ {\mathit{n}}_{\mathrm{F},3}(m,n)={\mathit{n}}_{\mathrm{R},3}(m,n)=({\widehat{\mathit{n}}}_{\nu ,3}^{\mathrm{\prime }}-{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }})·\frac{m+1}{M}+({\widehat{\mathit{n}}}_{\nu ,2}^{\mathrm{\prime }}-{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }})·\frac{\mathit{n}}{M}+{\widehat{\mathit{n}}}_{\nu ,1}^{\mathrm{\prime }}. \end{gathered}$$(15)

We denote ${\widehat{\mathit{n}}}_{\mathrm{F},i}(m,n)$ and ${\widehat{\mathit{n}}}_{\mathrm{R},i}(m,n)$ as the normalized vertex normals of the forward and reverse sub-triangles, respectively, as follows: $$\begin{gathered} {\widehat{\mathit{n}}}_{\mathrm{F},i}(m,n)=\frac{{\mathit{n}}_{\mathrm{F},i}(m,n)}{\parallel {\mathit{n}}_{\mathrm{F},i}(m,n)\parallel },i=1,2,3, \\ {\widehat{\mathit{n}}}_{\mathrm{R},i}(m,n)=\frac{{\mathit{n}}_{\mathrm{R},i}(m,n)}{\parallel {\mathit{n}}_{\mathrm{R},i}(m,n)\parallel },i=1,2,3. \end{gathered}$$

Using the simple Phong reflection model given by Eq. (14), we can now obtain the amplitudes $a_{\mathrm{F},1}$, $a_{\mathrm{F},2}$, and $a_{\mathrm{F},3}$ of the vertices of the forward sub-triangles, $\mathrm{F}(m,n)$, and the amplitudes $a_{\mathrm{R},1}$, $a_{\mathrm{R},2}$, and $a_{\mathrm{R},3}$ of the vertices of the reverse sub-triangles, $\mathrm{R}(m,n)$ as follows: $$\begin{gathered} a_{\mathrm{F},1}(m,n)=K_a+K_d[{\widehat{\mathit{n}}}_{\mathrm{F},1}(m,n)·\widehat{\mathit{L}}], \\ {a}_{\mathrm{R},1}(m,n)=K_a+K_d[{\widehat{\mathit{n}}}_{\mathrm{R},1}(m,n)·\widehat{\mathit{L}}], \\ {a}_{\mathrm{F},2}(m,n)=a_{\mathrm{R},2}(m,n)=K_a+K_d[{\widehat{\mathit{n}}}_{\mathrm{F},2}(m,n)·\widehat{\mathit{L}}], \\ {a}_{\mathrm{F},3}(m,n)=a_{\mathrm{R},3}(m,n)=K_a+K_d[{\widehat{\mathit{n}}}_{\mathrm{F},3}(m,n)·\widehat{\mathit{L}}], \end{gathered}$$(16)where $a_{\mathrm{F},2}=a_{\mathrm{R},2}$ and $a_{\mathrm{F},3}=a_{\mathrm{R},3}$ are evident from Fig. 3. Finally, the amplitude $g_F(x_r,y_r,m,n)$ of pixels $(x_r,y_r)$ within the forward sub-triangles and the amplitude $g_{\mathrm{R}}(x_r,y_r,m,n)$ of pixels $(x_r,y_r)$ within the reverse sub-triangles are linearly interpolated as $$\begin{gathered} g_{\mathrm{F}}(x_r,y_r,m,n) \\ =\{\begin{array}{ll}[a_{\mathrm{F},3}(m,n)-a_{\mathrm{F},1}(m,n)]x_r(m,n)+[a_{\mathrm{F},2}(m,n)-a_{\mathrm{F},1}(m,n)]y_r(m,n)+a_{\mathrm{F},1}(m,n),& \text{inside of}\mathrm{F}(m,n)\\ 0,& \text{outside of}\mathrm{F}(m,n)\end{array}. \end{gathered}$$(17)$$\begin{gathered} g_{\mathrm{R}}(x_r,y_r,m,n) \\ =\{\begin{array}{ll}[a_{\mathrm{R},2}(m,n)-a_{\mathrm{R},1}(m,n)]x_r(m,n)+[a_{\mathrm{R},3}(m,n)-a_{\mathrm{R},1}(m,n)]y_r(m,n)+a_{\mathrm{R},1}(m,n),& \text{inside of}\mathrm{R}(m,n)\\ 0,& \text{outside of}\mathrm{R}(m,n)\end{array}. \end{gathered}$$(18)

In our method, the subdivision is crucial for approximating Phong shading, as only subdivision can obtain the vertex normals of the sub-triangles through interpolation and calculate the amplitudes within the sub-triangles. Our work combines continuous shading^[6] and subdivision to derive the amplitude within the sub-triangles, ensuring that the amplitude within the sub-triangles is nonuniform.

For general index values $(m,n)$, we can find the spectra of forward sub-triangles $\mathrm{F}(m,n)$ and Reverse sub-triangles $\mathrm{R}(m,n)$, using the integral limits shown in Fig. 6 of Fu et al.’s method^[7]: $$\begin{gathered} G_{\mathrm{F}}(f_x^{\prime },f_y^{\prime },m,n)={\int }_{\frac{m}{M}}^{\frac{m+1}{M}}{\int }_{\frac{n}{M}}^{-x_r+\frac{m+n+1}{M}}{g}_{\mathrm{F}}(x_r,y_r,m,n)·\exp [-j2\pi (x_r{f}_x^{\prime }+y_r{f}_y^{\prime })]\mathrm{d}{x}_r\mathrm{d}{y}_r \\ =\frac{1}{M}\{[a_{\mathrm{F},3}(m,n)-a_{\mathrm{F},1}(m,n)]·D_1(\frac{f_x^{\prime }}{M},\frac{f_y^{\prime }}{M})+[a_{\mathrm{F},2}(m,n)-a_{\mathrm{F},1}(m,n)] \\ ·D_2(\frac{f_x^{\prime }}{M},\frac{f_y^{\prime }}{M})+a_{\mathrm{F},1}(m,n)·D_3(\frac{f_x^{\prime }}{M},\frac{f_y^{\prime }}{M})\}\exp [-j2\pi \left(\frac{m}{M}\right)f_x^{\prime }-j2\pi \left(\frac{n}{M}\right)f_y^{\prime }], \end{gathered}$$(19)$$\begin{gathered} G_{\mathrm{R}}(f_x^{\prime },f_y^{\prime },m,n)={\int }_{\frac{m}{M}}^{\frac{m+1}{M}}{\int }_{-x_r+\frac{m+n+1}{M}}^{\frac{n+1}{M}}{g}_{\mathrm{R}}(x_r,y_r,m,n)·\exp [-j2\pi (x_r{f}_x^{\prime }+y_r{f}_y^{\prime })]\mathrm{d}{x}_r\mathrm{d}{y}_r \\ =\frac{1}{M}\{[a_{\mathrm{R},2}(m,n)-a_{\mathrm{R},1}(m,n)]·D_1(-\frac{f_x^{\prime }}{M},-\frac{f_y^{\prime }}{M})+[a_{\mathrm{R},3}(m,n)-a_{\mathrm{R},1}(m,n)] \\ ·D_2(-\frac{f_x^{\prime }}{M},-\frac{f_y^{\prime }}{M})+a_{\mathrm{R},1}(m,n)·D_3(-\frac{f_x^{\prime }}{M},-\frac{f_y^{\prime }}{M})\}\times \exp [-j2\pi \left(\frac{m-1}{M}\right)f_x^{\prime }-j2\pi \left(\frac{n+1}{M}\right)f_y^{\prime }]. \end{gathered}$$(20)

The calculation results for $D_1(f_x^{\prime },f_y^{\prime })$, $D_2(f_x^{\prime },f_y^{\prime })$, and $D_3(f_x^{\prime },f_y^{\prime })$ have been given by Eqs. (11)–(13). The sum of all the spectra of the forward and reverse sub-triangles gives the spectra of the primitive triangle: $$G(f_x^{\prime },f_y^{\prime })=\sum _{m=0}^{M-1}\sum _{n=0}^{M-1}{G}_{\mathrm{F}}(f_x^{\prime },f_y^{\prime },m,n)+\sum _{m=0}^{M-2}\sum _{n=0}^{M-2}{G}_{\mathrm{R}}(f_x^{\prime },f_y^{\prime },m,n),$$(21)where $G(f_x^{\prime },f_y^{\prime })$ is the spectrum of a primitive triangle with what we call an approximate Phong shading rendering method. Finally, we obtain the distribution of the object’s light field with approximate Phong shading rendering information using Eq. (6): $$U_O(x,y)={\mathcal{F}}^{-1}\{\sum _{i=1}^n[\mathcal{J}·p·G(f_{xi}^{\prime },f_{yi}^{\prime })]\}.$$(22)

We would also like to point out that Wang et al.^[28] recently proposed a method for simulating Phong shading using segmentation. A significant difference between our work and the method proposed by Wang et al. is that it solely employs a triangular subdivision method to simulate Phong shading, where the amplitude within the sub-triangle is determined by the amplitude at the barycenter of the sub-triangle, resulting in a uniform amplitude. In contrast, our work combines self-similar subdivision and continuous shading methods, where the amplitude within the sub-triangle is determined by interpolating the amplitudes of the three vertices of the sub-triangle, resulting in nonuniform amplitudes.

To validate the proposed method, holograms of 3D objects with different numbers of triangles were generated based on Eq. (22). We utilized various shadow models and conducted numerical and optical reconstructions to assess the results. For numerical reconstruction, we have generated the holograms with a resolution of $1024\text{pixel}\times 1024\text{pixel}$ with a pixel pitch of 8 µm and wavelength of 0.532 µm. The hardware includes CPU: Intel(R) Xeon(R) Gold 6226 R CPU at 2.90 GHz without using a GPU. For optical reconstruction, the model of the spatial light modulator (SLM) is HOLOEYE PLUTO2 (NIR-011) with a resolution of $1920\times 1080$ (Full HD 1080p) and with a pixel pitch of 8 µm. The active area is $15.36\mathrm{mm}\times 8.64\mathrm{mm}$. The wavelength of the light wave emitted by the laser is 0.532 µm. The image receiver is an MMRY UC900C charge-coupled device.

Two holograms are generated using the continuous shading method proposed by Park et al. and our Phong shading approximation method, using a sphere as a first example. Figure 4 shows the numerical reconstruction results. Regarding the hologram reconstruction produced by Park et al.’s continuous shading method, as shown in Fig. 4(a), it can be clearly seen that there is a noticeable discontinuity between the edges of the polygon. Figures 4(b) and 4(c) show the reconstructed images of our method with different subdivisions $M$. In the proposed method shown in Fig. 4(c), a smooth appearance of the surface can be observed with a larger $M$, successfully eliminating the inappropriate brightness of the polygon edges. The lighting parameters employed for the Phong reflection model are given as follows : $K_a=0.2$, $K_d=0.8$, and $\widehat{\mathit{L}}=(0,0,1)$.

In Fig. 5(a), reconstruction without shadows is depicted, employing uniform amplitude across all triangles. Figure 5(b) displays the hologram reconstruction utilizing flat shading, where the amplitude within each triangle of the 3D object is uniform. This uniform amplitude results in a uniform brightness within each reconstructed mesh surface, causing a sudden change in brightness between adjacent meshes. Figure 5(c) illustrates the hologram’s reconstruction employing Park et al.’s continuous shading methodology. The continuous shading effect makes the mesh boundaries visible, reducing the realism of the reconstruction. Figure 5(d) illustrates the reconstruction utilizing the proposed technique, which achieves better surface approximation and eliminates polygon edges.

Figure 6 shows the numerical reconstruction results of the teapot using Wang et al.’s method and the proposed method for different $M$ values. From Fig. 6(a), it is evident that the method proposed by Wang et al. produced noticeable grid edges on the object’s surface at $M=3$, achieving a smooth effect at $M=16$. In contrast, Fig. 6(b) shows that the proposed method already attained a smooth surface at $M=3$. In Fig. 7, we compare the calculation time of five different methods using a teapot as an example. From Figs. 6 and 7, we see that the proposed method achieves a smooth effect at $M=3$ with a computation time of 175.32 s, while Wang et al.’s method achieves a smooth effect at $M=16$ with a computation time of 1447.43 s. In comparison, although the overall computation time of our proposed method increases with $M$ and is longer than that of Wang et al.’s method, the advantage of our method lies in achieving a smooth effect in less time. According to Fig. 7, it is evident that the proposed method takes more time than Park et al.’s method, but we can select an appropriate $M$ within the perceptible smooth range of the human eye to achieve the most efficient computation.

Finally, Fig. 8 illustrates the reconstruction of multiple objects. Table 1 shows the size information of the 3D objects used in the experiment. There is a significant difference in realism between the results using the proposed method and those obtained using Park et al.’s method. Therefore, the experimental results shown in Figs. 4, 5, 6, and 8 confirm that the use of Phong shading approximation in the proposed method effectively achieves the better-expected approximation of surface curvature.

In this study, we propose a new Phong shading approximation method of computer-generated holography based on fully analytical triangle meshes, which provides a complete analytical formalism based on Phong shading. The method divides the primitive triangles into multiple sub-triangles (both forward and reverse) and linearly interpolates the vertex normals of the primitive triangles to obtain the vertex normals for the sub-triangles. The Phong reflection model is then used to determine the amplitudes at the vertices of the sub-triangles. Finally, the amplitudes at the vertices of the sub-triangles are linearly interpolated to achieve spatially varying amplitudes within the sub-triangles. During reconstruction, the Phong shading approximation method can be used to render 3D objects, which almost eliminates the bright line artifacts on the surface of 3D objects that occur during holographic reconstruction, thereby improving reconstruction quality compared to the continuous shading methods proposed by Park et al.^[6] We have conducted simulations and optical reconstructions to validate our proposed approach.