Occlusion Culling for Wide-Angle Computer-Generated Holograms Using Phase Added Stereogram Technique

: A computer-generated hologram (CGH) allows synthetizing view of 3D scene of real or virtual objects. Additionally, CGH with wide-angle view offers the possibility of having a 3D experience for large objects. An important feature to consider in the calculation of CGHs is occlusion between surfaces because it provides correct perception of encoded 3D scenes. Although there is a vast family of occlusion culling algorithms, none of these, at the best of our knowledge, consider occlusion when calculating CGHs with wide-angle view. For that reason, in this work we propose an occlusion culling algorithm for wide-angle CGHs that uses the Fourier-type phase added stereogram (PAS). It is shown that segmentation properties of the PAS can be used for setting efﬁcient conditions for occlusion culling of hidden areas. The method is efﬁcient because it enables processing of dense cloud of points. The investigated case has 24 million of point sources. Moreover, quality of the occluded wide-angle CGHs is tested by two propagation methods. The ﬁrst propagation technique quantiﬁes quality of point reproduction of calculated CGH, while the second method enables the quality assessment of the occlusion culling operation over an object of complex shape. Finally, the applicability of proposed occlusion PAS algorithm is tested by synthetizing wide-angle CGHs that are numerically and optically reconstructed. key point of occlusion culling when using PAS is the allocation of frequency information of the sources that are closer to the hologram plane ﬁrst, which are labeled as front sources. These sources contribute the proper object information to the CGH. Other point sources of the same frequency but farther from the hologram plane are labeled as back sources, which must be disregarded. Thus, sorting of the cloud in respect to the z -axis allows identifying points from the nearest to the farthest in respect to the hologram plane.


Introduction
The aim of 3D display technology is to reproduce the depth perception of surrounding world. An important feature that must be considered in this technology is the ability to reproduce image providing a complete 3D experience to the viewers [1], i.e., full-parallax. Among several solutions that have been proposed in the last decades for developing 3D display technology, digital holography constitutes the best-known framework due to the ability of truly reconstructing the 3D wavefront [2]. The possibility of complex wavefront reproduction makes digital holography popular choice in several areas as microscopy [3], metrology [4], non-destructive sensing [5], security [6], design of diffractive optical elements [7]. In the realm of display technology, digital holography enables accurate reproduction of physiological cues [8,9], which is important for photorealistic 3-D images. However, digital holographic systems cannot capture scenes for wide-angle viewing display. With the development of computing technology, holograms can be registered without the necessity of a physical holographic setup. Computer-generated hologram (CGH) does not have limitations of physical systems [10] such as optical aberrations, extension, noise, or sensitivity to vibrations. Nevertheless, CGHs have limitations for providing the truly 3D experience. The first limitation is related with the pixel pitch of the spatial light modulators (SLMs) that are used for displaying CGHs. To have a full-parallax, pixel pitches below one micrometer must be used. High-end SLM technologies can offer SLMs with pixel pitch up to 3.74 µm and 4K resolution. Devices in research stage have shown that the pixel pitch The hologram is a set of dark and bright fringes that encode amplitude and phase information of the object. When using a digital camera for registering H, a digital hologram (DH) is obtained. DH can be processed numerically for reconstructing the encoded object information. However, capturing of DHs requires interferometric systems that are sensitive to optical aberrations, environmental and mechanical vibrations, or noise.
With the development of computing technology, holograms can be registered without the necessity of a physical holographic setup. When synthetizing CGH, only the object wave field O is calculated. Solutions for calculating this wavefield have been studied for more than three decades [30,[34][35][36][37]. The most popular approach consists in decomposing the 3D object into a set of illuminating point sources that are distributed in space. Each point source emits a spherical wavefront that propagates to the hologram plane. Superposition of all propagated wavefronts generate the CGH as follows where D p = (x 2 − x p ) 2 + (y 2 − y p ) 2 + z 2 p , (x p , y p , z p ) are the spatial position and a p amplitudes of the point sources, respectively. Nevertheless, this solution is time consuming [26,28]. Efficient techniques based on propagation have been proposed for calculating O [21,25,36,38,39]. With the aim to give the viewer a complete 3-D experience with fullparallax, wide-angle CGHs are required. The scheme of geometry for recording wide-angle CGHs is presented in Figure 1. The distance between hologram plane (x 1 , y 1 ) and object plane (x 2 , y 2 ) is z. Field of View (FoV) in a CGH is calculated as FoV = 2sin −1 (λ/2∆), where λ and ∆ are the wavelength and pixel pitch of the hologram, respectively. Thus, wide-angle CGHs are obtained by employing small pixel pitch at the hologram plane. For example, for given λ = 0.64 µm and ∆ = 0.62 µm the FoV ≈ 62 • . Large FoV allows synthetizing CGHs of objects with spatial extension B x2 = B y2 = 2z tan(FoV). These spatial extensions are much larger than the CGH size; for example, when using z = 1000 mm, B x2 = B y2 = 1209 mm while an a 4K CGH (4120 × 2060) will have the dimensions B x1 = 2.6 mm and B y2 = 1.3 mm.

Synthetizing CGH with the PAS
A hologram H is obtained as the result of measuring the intensity that results from the linear superposition of a diffracted object wavefield O and a reference wavefield R. The detected intensity can be mathematically expressed as follows 2 2 * * .
The hologram is a set of dark and bright fringes that encode amplitude and phase information of the object. When using a digital camera for registering H, a digital hologram (DH) is obtained. DH can be processed numerically for reconstructing the encoded object information. However, capturing of DHs requires interferometric systems that are sensitive to optical aberrations, environmental and mechanical vibrations, or noise.
With the development of computing technology, holograms can be registered without the necessity of a physical holographic setup. When synthetizing CGH, only the object wave field O is calculated. Solutions for calculating this wavefield have been studied for more than three decades [30,[34][35][36][37]. The most popular approach consists in decomposing the 3D object into a set of illuminating point sources that are distributed in space. Each point source emits a spherical wavefront that propagates to the hologram plane. Superposition of all propagated wavefronts generate the CGH as follows are the spatial position and ap amplitudes of the point sources, respectively. Nevertheless, this solution is time consuming [26,28]. Efficient techniques based on propagation have been proposed for calculating O [21,25,36,38,39]. With the aim to give the viewer a complete 3-D experience with full-parallax, wide-angle CGHs are required. The scheme of geometry for recording wide-angle CGHs is presented in Figure 1. The distance between hologram plane (x1, y1) and object plane (x2, y2) is z. Field of View (FoV) in a CGH is calculated as FoV = 2sin −1 (λ/2Δ), where λ and Δ are the wavelength and pixel pitch of the hologram, respectively. Thus, wideangle CGHs are obtained by employing small pixel pitch at the hologram plane. For example, for given λ = 0.64 µm and Δ = 0.62 µm the FoV ≈ 62°. Large FoV allows synthetizing CGHs of objects with spatial extension Bx2 = By2 = 2z tan(FoV). These spatial extensions are much larger than the CGH size; for example, when using z = 1000 mm, Bx2 = By2 = 1209 mm while an a 4K CGH (4120 × 2060) will have the dimensions Bx1 = 2.6 mm and By2 = 1.3 mm.  Unfortunately, propagation-based techniques are not an optimal approach for synthetizing CGHs of wide-angle, since they do not consider geometry, where the object is much larger than the hologram [25]. This geometrical feature is challenging for rigorous propagation algorithms. Alternatives for efficient synthesis of CGHs are PAS-based solutions. These solutions do not have the restrictions of propagation-based methods since they consider the basic idea of summing point source responses at the hologram plane [28,31,32]. The PAS algorithm divides the CGH into smaller segments of size N s × N s . The number of segments within the CGH is N × M, where N is the number of segments in the horizontal direction and M in the vertical direction. At each segment, PAS approximates the spherical wavefront that comes from the p-th point source of the object with a local tilted plane wave [31], as shown in Figure 2a. The direction of the plane wave is calculated according to the value of the local instantaneous frequency at the center of the segment (x c n , y c m ), where n = 1, . . . , N, and m = 1, . . . , M indicates segment position over the horizontal and vertical direction, respectively. Since the frequency information of each segment is coded according to the local approximation of the tilted plane, the PAS can calculate the diffracted field of objects larger than the hologram size. The term O for a given segment (m,n) is calculated with the PAS as follows where D pmn is the distance from p-th point source to the center of the segment (m,n), and are the local instantaneous frequencies. The angles θ xpnm , θ ypnm , θ rxnm , and θ rynm are the incident angles of object and reference beams on x-axis and y-axis in relation to the center of the segment, respectively. Notably, the hologram in the corresponding segment is a superposition of plane waves. When the reference beam has a constant incident angle, spatial frequencies depend on the object beam [31].
they consider the basic idea of summing point source responses at the hologram plane [28,31,32]. The PAS algorithm divides the CGH into smaller segments of size Ns × Ns. The number of segments within the CGH is N × M, where N is the number of segments in the horizontal direction and M in the vertical direction. At each segment, PAS approximates the spherical wavefront that comes from the p-th point source of the object with a local tilted plane wave [31], as shown in Figure 2a. The direction of the plane wave is calculated  according to the value of the local instantaneous frequency at the center of the segment   ( )   ,   c  c  n  m x y , where n = 1, …, N, and m = 1, …, M indicates segment position over the horizontal and vertical direction, respectively. Since the frequency information of each segment is coded according to the local approximation of the tilted plane, the PAS can calculate the diffracted field of objects larger than the hologram size. The term O for a given segment (m,n) is calculated with the PAS as follows ( ) where Dpmn is the distance from p-th point source to the center of the segment (m,n), and sin sin , are the local instantaneous frequencies. The angles θxpnm, θypnm, θrxnm, and θrynm are the incident angles of object and reference beams on x-axis and y-axis in relation to the center of the segment, respectively. Notably, the hologram in the corresponding segment is a superposition of plane waves. When the reference beam has a constant incident angle, spatial frequencies depend on the object beam [31].  Direct computation of Equation (3) can still be slow for a dense cloud of points. Nevertheless, this process can be significantly accelerated when employing the FT. When representing Equation (3) in frequency domain, the plane waves in the segment are transformed into a 2D distribution of Dirac deltas. The positions of these deltas need to be set to a discrete value that is specified by the segment size. Hence, Equations (4) and (5) are rewritten as:f where ∆f = 1/N s ∆ is the frequency sampling, and · is the rounding operation. When employing Equations (6) and (7), the term O R is expressed in the frequency domain aŝ After placing all the Dirac deltas for the 3D point distribution into the segment, the inverse FT is applied, and then the complex field for a given segment is obtained. This procedure is repeated for all segments for obtaining the CGH, as shown in Figure 2b.

Occlusion Culling of Back Points with the PAS
Representation of 3D object for point-based approach offers a simple and flexible way of depicting its shape. Moreover, this representation allows associating the sources with additional information such as colors, textures, and normals [40]. Each point source generates a complex field from the scene towards the hologram plane, and the final CGH is the sum of all corresponding spherical waves. In the point-based representation, point sources cannot occlude one another, unless they accidentally fall along the same direction from the viewpoint. The wavefields from back and front points are encoded in the CGH. This prevents correct reproduction of depth cues of the 3D object in holographic image. In order to eliminate the back points, occlusion techniques are needed. Occlusion techniques can rely on geometric considerations [40], light shielding [17], or hybrid [15,23]. Since the CGH must provide motion parallax, the chosen occlusion technique must consider that the visibility of object changes according to the movement of the viewer.
Occlusion can be introduced by employing the segmentation property of the PAS [29]. As shown in Equation (3), the term O for a given segment (m,n) is a superposition of plane waves having directional wave vectors calculated according to the center of the segment. When applying the FT to the segment, the coordinates of 3D point sources are mapped into the sets of 2D Dirac deltas in frequency domain. The frequency support of this segment is given by ∆ −1 . Since frequency allocation depends on the size of the segment, the gaps between adjacent Dirac deltas can be controlled in the frequency domain by the value of N s . Figure 3 illustrates this effect for three point sources when using different segments sizes. The point sources x p , and x p+1 are adjacent front points, while x m is a back point. The angles of corresponding plane waves are θ p , θ p+1 , and θ m , respectively. In this example, the angles fulfill the relation θ p < θ m < θ p+1 . Figure 3a presents the case of a small segment. In this case, the wavefront information for all point sources is located in the same frequency pixel sincef p =f p+1 =f m . For carrying out the occlusion, information related with the point source x m needs to be disregarded because it has the farthest separation from the hologram plane. When enlarging the segment, as shown in Figure 3b, wavefront information of front points x p and x p+1 is allocated in frequency pixelsf p andf p+1 , respectively. In this case, the information related to the source x m will be placed to the nearest neighborf p+1 . Similarly to Figure 3a, the information related to the source x m needs to be disregarded for the occlusion culling. Finally, for large segment the wavefront information of each point source is stored in independent pixels, as presented in Figure 3c. As seen in Figure 3, the key point of occlusion culling when using PAS is the allocation of frequency information of the sources that are closer to the hologram plane first, which are labeled as front sources. These sources contribute the proper object information to the CGH. Other point sources of the same frequency but farther from the hologram plane are labeled as back sources, which must be disregarded. Thus, sorting of the cloud in respect to the z-axis allows identifying points from the nearest to the farthest in respect to the hologram plane.
the nearest neighbor 1p f + . Similarly to Figure 3a, the information related to the source xm needs to be disregarded for the occlusion culling. Finally, for large segment the wavefront information of each point source is stored in independent pixels, as presented in Figure  3c. As seen in Figure 3, the key point of occlusion culling when using PAS is the allocation of frequency information of the sources that are closer to the hologram plane first, which are labeled as front sources. These sources contribute the proper object information to the CGH. Other point sources of the same frequency but farther from the hologram plane are labeled as back sources, which must be disregarded. Thus, sorting of the cloud in respect to the z-axis allows identifying points from the nearest to the farthest in respect to the hologram plane. With the aim to calculate the maximum segment size that enables occlusion, let's consider a set of points that are distributed in the plane (x2,z). Then, adjacent points xp and xp+1 and an arbitrary hologram segment with center at c n x are chosen. The difference between corresponding local frequencies is given by Note that θ(p+1)nx = θpnx + Δθ, where Δθ is the angular difference between adjacent sources. To eliminate any gap between adjacent frequencies, the right-hand side of Equation (9) must be equal or smaller than the frequency sampling of the segment. This can be expressed as The denominator on the right-hand side of Equation (10) has been obtained by using the sum-to-product identity. When meeting this condition, the two wavefronts corresponding to the discrete frequencies pxn  With the aim to calculate the maximum segment size that enables occlusion, let's consider a set of points that are distributed in the plane (x 2 ,z). Then, adjacent points x p and x p+1 and an arbitrary hologram segment with center at x c n are chosen. The difference between corresponding local frequencies is given by Note that θ (p+1)nx = θ pnx + ∆θ, where ∆θ is the angular difference between adjacent sources. To eliminate any gap between adjacent frequencies, the right-hand side of Equation (9) must be equal or smaller than the frequency sampling of the segment. This can be expressed as The denominator on the right-hand side of Equation (10) has been obtained by using the sum-to-product identity. When meeting this condition, the two wavefronts corresponding to the discrete frequencies f pxn and f (p+1)xn are occupying two adjacent discrete locations in the frequency matrix. When considering a back point x m , which view angle θ m lies between θ p and θ p+1 , the holographic information associated with the frequency f mxn can be accommodated into the nearest neighbor f pxn or f (p+1)xn . Nevertheless, the contribution from f mxn can be disregarded since the frequency positions are already occupied by the front point sources, as shown in Figure 3b. Finally, maximum segment size can be found by solving Equation (10) for N s , which yields Notably, Equation (11) allows calculating the segment size that enables occlusion culling of back points. Once the occlusion culling of back points is complete, inverse FT is applied, and thus, the hologram segment with the proper depth perspective is obtained.
It is important to differentiate two cases in the occlusion culling. The first case can be exemplified by Figure 3a. In this case, a small segment is employed for the occlusion. As depicted in this figure, the holographic information of two sources (or more) is allocated in the same frequency pixel. This means that resolution between the sources is lost. Thus, the smaller the segment is, the smaller the resolution in the CGH. On the other hand, when the chosen segment size is closer to N s , wavefront information of the sources is placed in different frequency pixels, as presented in Figure 3b, which allows preserving resolution in the reconstructed image. Finally, it is important to note that position of the segment shifts allocation of the discrete frequencies but not the condition imposed by Equation (11). For segments that are too far from the center of the CGH, the shifting causes discrete frequencies to fall outside from the frequency support ∆ −1 .
Visualization of the occlusion culling for a segment in the frequency domain is presented in Figure 4. The employed object is a point cloud of a gargoyle figurine with dimensions 1055 mm (height) × 971 mm (width) × 485 mm (depth) composed of approximately 24M points with a density of 9 points/mm 2 . The distance from the center of the object to the hologram plane is one meter. Employed features allow calculating an average angular spacing of 0.038 • . Employed features allow calculating an average angular spacing of 0.035 • . When employing, ∆ = 0.62 µm and λ = 640 nm, which enables maximum angle view of 62.5 • , it is found that N s must be smaller than 1550 pixels. Thus, segment size is selected as N s = 1024. The size of the calculated CGH is chosen as 4096 × 4096 pixels. In terms of the PAS technique, the CGH consists of 4 segments of size N s × N s , and thus N, M = 4. Figure 4 shows the amplitudes of the segment m = 1 and n = 1 in the frequency domain. Figure 4a depicts the cloud of points without occlusion. Image shows that back and front points contribute to the holographic information of the segment. Thus, the depth cues are not correctly displayed. Figure 4b shows the cloud of points when applying the proposed occlusion culling based PAS. After the occlusion culling of back points, only the front of the gargoyle is seen, which is the correct observable surface. Later, the cloud of points is rotated 180 • and PAS with culling occlusion is applied again. Hence, the back of the gargoyle is obtained. The result is shown in Figure 4c. These examples show that our method can eliminate unwanted back points by using the segmentation properties of the PAS.
in the same frequency pixel. This means that resolution between the sources is lost. Thus, the smaller the segment is, the smaller the resolution in the CGH. On the other hand, when the chosen segment size is closer to Ns, wavefront information of the sources is placed in different frequency pixels, as presented in Figure 3b, which allows preserving resolution in the reconstructed image. Finally, it is important to note that position of the segment shifts allocation of the discrete frequencies but not the condition imposed by Equation (11). For segments that are too far from the center of the CGH, the shifting causes discrete frequencies to fall outside from the frequency support Δ −1 .
Visualization of the occlusion culling for a segment in the frequency domain is presented in Figure 4. The employed object is a point cloud of a gargoyle figurine with dimensions 1055 mm (height) × 971 mm (width) × 485 mm (depth) composed of approximately 24M points with a density of 9 points/mm 2 . The distance from the center of the object to the hologram plane is one meter. Employed features allow calculating an average angular spacing of 0.038°. Employed features allow calculating an average angular spacing of 0.035°. When employing, Δ = 0.62 µm and λ = 640 nm, which enables maximum angle view of 62.5°, it is found that Ns must be smaller than 1550 pixels. Thus, segment size is selected as Ns = 1024. The size of the calculated CGH is chosen as 4096 × 4096 pixels. In terms of the PAS technique, the CGH consists of 4 segments of size Ns × Ns, and thus N, M = 4. Figure 4 shows the amplitudes of the segment m = 1 and n = 1 in the frequency domain. Figure 4a depicts the cloud of points without occlusion. Image shows that back and front points contribute to the holographic information of the segment. Thus, the depth cues are not correctly displayed. Figure 4b shows the cloud of points when applying the proposed occlusion culling based PAS. After the occlusion culling of back points, only the front of the gargoyle is seen, which is the correct observable surface. Later, the cloud of points is rotated 180° and PAS with culling occlusion is applied again. Hence, the back of the gargoyle is obtained. The result is shown in Figure 4c. These examples show that our method can eliminate unwanted back points by using the segmentation properties of the PAS. In its basic configuration, the PAS algorithm introduces a truncation error when discretizing the instantaneous local frequencies (Equations (6) and (7)) [32], which decreases the reconstruction quality of the CGH. This truncation error depends on the segment size and can be easily reduced by using a large segment. In the general case of CGH calculation, arbitrarily large segment cannot be selected since the wavefront at the hologram plane may not be well represented. For the proposed occlusion culling PAS, the restriction of the segment becomes stricter due to the limitation imposed with Equation (11). Accurate PAS (A-PAS) can be employed for solving the discretization problem while fulfilling the N s constraint and sampling theory [31] at the same time. The A-PAS considers that the segment is enlarged in the discrete frequency domain by an integer factor r. In this way, the frequency sampling is increased r-times, which allows a better approximation of the discrete values of the frequencies to their respective continuous values. After applying the inverse FT to the enlarged segment, the original segment is taken. The A-PAS is implemented in the occluded PAS algorithm by considering that the enlarged segment in frequency domain is not larger than N s . This means that we must find a hologram segment of size N so and extension factor r such that N s > N so r. When selecting proper N so and r, the accuracy of the PAS algorithm can be reached while preserving the occlusion culling of unwanted back points.

Numerical Reconstruction Methods for Wide View-Angle CGHs
With the aim to measure accuracy and quality of the holographic image encoded in the CGH, propagation methods [33,41,42] must be employed. In this section, we investigate two propagation techniques. The first, which is based on the direct integration method [33] is used to assess the accuracy of the CGH calculated by the PAS. The second, that is a generalization of the MFFT-AS technique [10], allows for reconstructing the whole object and visualizing quality of the PAS based occlusion culling.

Off Axis Direct Integration Method
Recovering the full complex field from wide-angle CGH with classical propagation techniques requires a huge amount of zero padding in x and y directions. Thus, an enormous amount of computational power is required for storing and operating such information. Moreover, the time that is necessary for obtaining the reconstruction will be on the scale of hours or even days. Therefore, it is proposed to employ the integration method [33] for reconstructing a very small off-axis window of the entire complex field. This is done with the aim to carry out a fast quantification of point reproduction quality of calculated CGH.
Off-axis wavefield recovery can be carried out using the Fourier-based solver of the Rayleigh-Sommerfeld integral [43]: where x 2 = x 1 − x oa , y 2 = y 1 − y oa , and (x oa , y oa ) are the coordinates of transversal shift, u in is the complex field from the hologram, u out is the reconstructed field, the subindexes 1 and 2 point out the input and output planes, respectively, and where D = x 2 + y 2 + z 2 . Application of the FFT in Equation (12) generates a circular convolution affecting accuracy of reconstruction. Influence of circular convolution is avoided by increasing twice size of the input matrix by using zero-padding. When solving Equation (12) with these considerations, the desired field is given by the N x × N y lower right submatrix, where N x is the size of the input field. This method preserves the pixel size between input at output planes. Thus, the size of the integrated window is x 2 ∈ (−N x ∆/2, N x ∆/2), y 2 ∈ (−N x ∆/2, N x ∆/2). This off-axis direct integration (OADI) method enables fast wavefield reconstruction from any point of view with small amount of memory [33,41,44] but for very small window.

Multi Fast Fourier Transform Angular Spectrum Method
For the reconstruction of full holographic data from the wide-angle CGH, a generalization of the MFFT-AS method is proposed [10]. The MFFT-AS was initially employed for the reconstruction of Fourier horizontal parallax-DHs (FHPO-DHs) [45], which are small in vertical and large in horizontal directions, respectively. This method applies (q+1)N x × N y FFTs instead of two large FFTs from the classical angular spectrum (AS) method, where N x , Photonics 2021, 8, 298 9 of 16 N y are the dimensions of the FHPO-DH and q = B y2 /B y1 is the compression factor for y direction (B y1 is vertical dimension of the hologram). A critical step in the MFFT-AS is the summation in form of "tiles" in y direction. This makes possible to significantly reduce the computational load. Hence, reconstruction of FHPO-DH with AS is possible. In this work, we extend the MFFT AS application to two dimensions, which allows obtaining accurate reconstruction of full parallax CGH.
Implementation of the MFFT-AS can be described as follows: let O be the 2D object wave of N x × N y samples and spatial dimensions B x1 × B y1 . We wish to find its diffracted field at plane (x 2 , y 2 ). The MFFT-AS approach accomplishes this task by calculation of p × q FFTs as where p = B x2 /B x1 is the compression factor in horizontal direction, lp+α, jq+β represent the frequency coordinates separated by 1/pB x1 and 1/qB y1 , respectively. Equation (14) is multiplied by a shifted propagation kernel, which yields It is worth noting that the Equation (15) has size of N x p × N y q; however, this large matrix is never allocated. Sequentially, calculated vectors of size N x × N y with Equation (15) are assembled as matrices of the same length N x × N y by next part of the MFFT-AS method. This element of the MFFT-AS is based on the property of FFT that was proposed in references [46] and allows efficient calculations of high-resolution focused wave of small spatial extend. The property enables reduction of bandwidth of FFT without the need of computing FFT of the full-sized vector. To explain the property, let assume a 1D signal A of length qN y . The method starts with summing up "tiles" as A s (n) = ∑ q−1 t=0 A(tq + n), where n ∈ 0, . . . , N y − 1 and then computes FFT. Notably, the A s is a vector, which contains q-th frequency of A s . In the 2D MFFT-AS method, the tiling is applied in both directions and the introduced property enables reduction of the sampling rate of the calculated diffracted field. Thus, the result is calculated using at first 'tiles' summation method, and then FT as Notably, full field reconstruction with the 2D MFFT-AS AS requires q × p + 1 FFTs of size N x × N y , and pixel size at the output plane is determined by ∆x 2 = p ∆x 1 , and ∆y 2 = q ∆y 1 .

Comparison of Full and Small off Axis Window Reconstruction
The strategies based on off axis window and full image reconstruction are here compared from the computational point of view. For this, the set of holograms for reconstruction distance z = 1000 mm, pixel count 4096 × 4096, and several pixel pitches were generated and reconstructed. Table 1 presents the image size and the obtained computation times for both reconstruction methods when using selected pixel pitches. Notably, the MFFT-AS enables to obtain image sizes ranging from one quarter of meter up to one meter. However, the price to pay is a large processing time. For example, reconstruction of the smallest pixel pitch (∆ = 0.6 µm) requires almost 12 h of processing. On the other hand, the OADI method provides results in a few seconds for reconstructing the complex field but of window size of a few millimeters.

Diffraction Efficiency Based on Assessment of CGH Quality
In this section, assessment of CGH quality based on diffraction efficiency measure is carried out. This is done with the aim to quantify quality of point reproduction of calculated CGH. Reconstruction of wide-angle CGHs require, see first row of Table 1, a large amount of memory and processing time. Nevertheless, for measuring the diffraction efficiency of CGHs of this work, we retrieve a small region where a point source is reconstructed. Hence, the OADI method is suitable for measuring the diffraction efficiency of CGHs calculated with the PAS algorithm.
The first diffraction efficiency-based test is evaluated with two CGHs of ten point sources distributed linearly along the diagonal from the center of the image P 0 (0, 0, z) up to corner point P 10 = (−0.475B x2 , −0.475B y2, z), where z = 1000 mm. Pixel pitch of the CGHs is ∆ = 0.62 µm. Theoretical CGH, which is reference one, is calculated using Equation (2), while the second CGH with the PAS algorithm employing N so = 256 and r = 4. Calculated CGHs are processed with Equation (12) and the corresponding off-axis coordinate. As a result, the single point source with its nearest neighborhood is reconstructed. Figure 5a,b present the reconstruction of point P 10 of theoretical and PAS CGHs, respectively. Visually, reconstructed point sources have identical profiles and almost the same peak value (P-V). Figure 5c presents the accuracy test in relation to point source position along the diagonal of image space. Maximum amplitude values are normalized with respect to the on-axis point P 0 from theoretical CGH. Regardless, the CGH is made with PAS method or with Equation (2), the linear loss of amplitude due to point source position is present. Moreover, it is clearly seen in the plot that accuracy of reconstructed point sources is strongly related to the selection of the parameter r, e.g., for r ≥ 3, the maximum amplitude of each point starts to behave stable and is close to the result of theoretical CGH. Accordingly, CGHs of high accuracy made with PAS can be achieved by using a value N so and large r. Nevertheless, the constraint provided by Equation (11) must be fulfilled all the time if occlusion culling is desired. For this reason, we perform a test that is aimed to measure the accuracy of calculated CGHs with the PAS with different parameters N so and r. As a measure of accuracy, the ratio of maximum amplitude of reconstructed point P 10 between PAS and theoretical CGH is used. Corresponding results are shown in Figure 6a. The obtained ratios demonstrate that when using original segment size (r = 1), the diffraction efficiency is low and unstable, which means that calculated CGH with PAS is not optimal. When using large values of r (r > 2), the maximum amplitude increases its value independently of the N so , which is pointed out by large dash box into the plot.
r. As a measure of accuracy, the ratio of maximum amplitude of reconstructed point P10 between PAS and theoretical CGH is used. Corresponding results are shown in Figure 6a. The obtained ratios demonstrate that when using original segment size (r = 1), the diffraction efficiency is low and unstable, which means that calculated CGH with PAS is not optimal. When using large values of r (r > 2), the maximum amplitude increases its value independently of the Nso, which is pointed out by large dash box into the plot.  The second test is carried out for measuring the computation time of calculated CGH with a point cloud of one million points. The cloud has a similar point density to the gargoyle figurine presented in Section 3. Different Nso and r values without and with occlusion culling routine are employed. The results are presented in Figure 6b. Note that PAS calculation of a segment consists of two parts: allocation of holographic information according to the local instantaneous discrete frequencies and performing FFT. Time consumption for allocation depends on the number of points, and thus the larger the cloud is, the larger the computational effort. Time consumption is particularly demanding when the number of points within the cloud is counted in millions. In comparison, the computational effort of the FFT only depends on the size of the frequency matrix (Nso × r), where the holographic information is stored. Since the computational complexity of FFT is given between PAS and theoretical CGH is used. Corresponding results are shown in Figure 6a. The obtained ratios demonstrate that when using original segment size (r = 1), the diffraction efficiency is low and unstable, which means that calculated CGH with PAS is not optimal. When using large values of r (r > 2), the maximum amplitude increases its value independently of the Nso, which is pointed out by large dash box into the plot.  The second test is carried out for measuring the computation time of calculated CGH with a point cloud of one million points. The cloud has a similar point density to the gargoyle figurine presented in Section 3. Different Nso and r values without and with occlusion culling routine are employed. The results are presented in Figure 6b. Note that PAS calculation of a segment consists of two parts: allocation of holographic information according to the local instantaneous discrete frequencies and performing FFT. Time consumption for allocation depends on the number of points, and thus the larger the cloud is, the larger the computational effort. Time consumption is particularly demanding when the number of points within the cloud is counted in millions. In comparison, the computational effort of the FFT only depends on the size of the frequency matrix (Nso × r), where the holographic information is stored. Since the computational complexity of FFT is given The second test is carried out for measuring the computation time of calculated CGH with a point cloud of one million points. The cloud has a similar point density to the gargoyle figurine presented in Section 3. Different N so and r values without and with occlusion culling routine are employed. The results are presented in Figure 6b. Note that PAS calculation of a segment consists of two parts: allocation of holographic information according to the local instantaneous discrete frequencies and performing FFT. Time consumption for allocation depends on the number of points, and thus the larger the cloud is, the larger the computational effort. Time consumption is particularly demanding when the number of points within the cloud is counted in millions. In comparison, the computational effort of the FFT only depends on the size of the frequency matrix (N so × r), where the holographic information is stored. Since the computational complexity of FFT is given by (N so × r) × log(N so × r), then a small product of N so and r will result in a short computation time, which is almost negligible, while a large product between N so and r, will increase it. However, this FFT-related time is still shorter than for point processing. A benefit of using large segments N so is that less segments are computed for hologram generation, which speeds up the CGH computation. When measuring time for the occlusion culling, it is seen that the occlusion PAS is faster than normal PAS. For example, generation of the CGH without occlusion for N so = 32 and r = 3, and N so = 512 and r = 3 took around 674 and 4.6 s, respectively. On the other hand, calculation of the CGH with occlusion for previous parameters took around 297 s and 3.9 s, respectively. Notably, for segments of size N so × r that do not fulfill Equation (11) the growth of computation time is higher for the occlusion culling than for the computation without occlusion. This is because of the finer frequency sampling of the segment, which means that back points are allowed to store their holographic information. Therefore, the number of processed points grows, and thus the computation time as well.
As discussed in former paragraph, selection of N so and r for ensuring high accuracy and fast execution time must consider the constraint imposed by Equation (11). However, Equation (11) relies in the angular spacing ∆θ. The parameter ∆θ depends on the geometry of the object, the dimension of the hologram, and the separation between object and hologram planes. Due to the mentioned features, the estimation of the angular spacing is a non-trivial issue. Here, it is decided to estimate the local angular spacing with the fitting of a Gaussian distribution to the computed local density of points. Also, considering the standard deviations, we estimate the quasi-minimal angular spacing value for the points clouds that statistically work properly for over 99% of the object surface. Table 2 presents the connection between the average angular spacing of the point cloud and the proper segment size N so as it is stated in Equation (11). Four point clouds with the same geometry from gargoyle figure from Section 3 are employed. These clouds have 200 thousand, 1 million, 5 million, and 24 million of points. Investigation of these exemplary point clouds presents the proper choice of the segment size due to the point cloud angular spacing and their respective CGH time calculation.

Reconstruction of Wide-Angle CGHs Calculated with the Occlusion Culling Based-PAS Algorithm
In this section, numerical and optical reconstructions of wide-angle CGHs obtained with the occlusion culling based PAS are carried out. PAS algorithm is implemented in Matlab software, which is installed on a PC equipped with Intel i9-9900K processor, 32 GB of RAM, and run without support of GPU. Selected size of the CGHs is 4160 × 2160 pixels, which corresponds to the size of the 4K SLM Holoeye GAEA 2,0, pixel pitch ∆ SLM = 3.74 µm. CGHs are calculated from the cloud of 24 million points defined in Section 3 and PAS parameters: N so = 512, r = 3, and ∆ = 0.62 µm, which according to Table 2 fulfills N s limit. Calculation time for each CGH is around 190 s, which can be reduced up to 30 s when using the Matlab parallel toolbox. When comparing this time with Table 1 from reference [29], it can be seen that our approach improves the calculation time by a factor of 3000. Figure 7a-c show numerical reconstructions for three perspectives angles: 0 • , 45 • , and 180 • of CGH calculated with the 2D MFFT-AS, respectively. Notably, the application of the occlusion culling based PAS algorithm allows eliminating unwanted back points, and hence, reconstruction of the object with proper depth cues. It is worth noting that SFI method is not able to provide any of the presented results since it is a method designed under the constraints of the paraxial regime. Video S1 presents the numerical reconstruction of the gargoyle view angles from 0 • to 360 • . The cloud of points was rotated 8 • between each frame. the occlusion culling based PAS algorithm allows eliminating unwanted back points, and hence, reconstruction of the object with proper depth cues. It is worth noting that SFI method is not able to provide any of the presented results since it is a method designed under the constraints of the paraxial regime. Visualization 1 presents the numerical reconstruction of the gargoyle view angles from 0° to 360°. The cloud of points was rotated 8° between each frame. Optical reconstruction of CGH is carried out in the display setup presented in Figure  8. Illumination laser source (λR = 640 nm) hits half-wave plate (HWP), which controls polarization state. Next, the beam passes through microscopic objective (MO) and a pinhole (PH), generating a spherical wavefront. A collimator lens C (Fc = 300 mm) is placed in front of the pinhole to obtain an extended plane wave. Plane wave illumination is directed to the phase only 4K SLM by a non-polarized beam splitter. SLM is imaged by an    3 ]. Since the dimension of the cloud is reduced, the average angular spacing of the point cloud decreases to 0.030°. For this value, the NS must be smaller than 1900 pixels. Hence, former PAS parameters fulfill accuracy and speed for CGH calculation as well. Obtained results prove that the occlusion culling based-PAS can eliminate unwanted back points for each angle. Visualization 2 presents Optical reconstruction of CGH is carried out in the display setup presented in Figure 8. Illumination laser source (λ R = 640 nm) hits half-wave plate (HWP), which controls polarization state. Next, the beam passes through microscopic objective (MO) and a pinhole (PH), generating a spherical wavefront. A collimator lens C (F c = 300 mm) is placed in front of the pinhole to obtain an extended plane wave. Plane wave illumination is directed to the phase only 4K SLM by a non-polarized beam splitter. SLM is imaged by an afocal demagnifying system composed of an objective (F 1 = 200 mm, Ø 1 = 50 mm) and an eyepiece (F 2 = 33 mm, Ø 2 = 25 mm). At the viewing window (VW) plane pixel pitch of the image of SLM is ∆' SLM = F 2 /F 1 × ∆ SLM = 0.62 µm. The small pixel pitch allows for a wide viewing angle (62.5 • ). However, the physical aperture of the eyepiece limits it to 36 • . method is not able to provide any of the presented results since it is a metho under the constraints of the paraxial regime. Visualization 1 presents the nume struction of the gargoyle view angles from 0° to 360°. The cloud of points was between each frame.    3 ]. Since the dimension of the duced, the average angular spacing of the point cloud decreases to 0.030°. For the NS must be smaller than 1900 pixels. Hence, former PAS parameters fulfi and speed for CGH calculation as well. Obtained results prove that the occlus based-PAS can eliminate unwanted back points for each angle. Visualization   3 ]. Since the dimension of the cloud is reduced, the average angular spacing of the point cloud decreases to 0.030 • . For this value, the N S must be smaller than 1900 pixels. Hence, former PAS parameters fulfill accuracy and speed for CGH calculation as well. Obtained results prove that the occlusion culling based-PAS can eliminate unwanted back points for each angle. Video S2 presents the optical reconstruction of the gargoyle from 0 • to 360 • . Here, the cloud of points was rotated 1 • per frame.
OR PEER REVIEW 14 of 16 the optical reconstruction of the gargoyle from 0° to 360°. Here, the cloud of points was rotated 1° per frame.

Conclusions
In this work, we proposed a simple, effective and fast occlusion culling strategy for wide-angle view holograms. Our approach takes advantage of PAS algorithm, which is based on mapping of 3D object information from spatial to frequency domain. It is shown that optimal segment size that enables occlusion of back points can be calculated according to the average angular spacing of the point cloud. Moreover, sorting of the cloud of points in the z-axis allows ordering the point sources from the nearest to the farthest respect to the hologram plane. In this way, sources that are closer to the hologram plane are labeled as front ones and allocate their holographic information into the frequency matrix first. Other point sources of the same frequency but farther from the hologram plane are labeled as back ones, and thus, their information into the frequency matrix is omitted. Moreover, two methods for numerical reconstruction of wide-angle view CGH were studied. The first was the OADI method, which allows reconstructing a small window from the whole image space. In this way, efficient numerical technique is obtained, which enables measuring the diffraction efficiency of the PAS algorithm. The second was the 2D MFFT-AS, which allows reconstructing the full image space. It was employed for reconstructing the whole object and visualize the quality of the occlusion culling that was carried out with the PAS. Finally, optical reconstruction of occluded CGHs in wide-angle holographic display verified the validity of our occlusion culling strategy. Notably, our methodology enables fast computation time for obtaining occluded CGHs from very large point clouds.  Data Availability Statement: Data underlying the results presented in this paper are publicly available at this time but may be obtained from the authors upon request.

Acknowledgments:
The authors would like to thank to the Warsaw University of Technology for the statutory funds provided to this research.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
In this work, we proposed a simple, effective and fast occlusion culling strategy for wide-angle view holograms. Our approach takes advantage of PAS algorithm, which is based on mapping of 3D object information from spatial to frequency domain. It is shown that optimal segment size that enables occlusion of back points can be calculated according to the average angular spacing of the point cloud. Moreover, sorting of the cloud of points in the z-axis allows ordering the point sources from the nearest to the farthest respect to the hologram plane. In this way, sources that are closer to the hologram plane are labeled as front ones and allocate their holographic information into the frequency matrix first. Other point sources of the same frequency but farther from the hologram plane are labeled as back ones, and thus, their information into the frequency matrix is omitted. Moreover, two methods for numerical reconstruction of wide-angle view CGH were studied. The first was the OADI method, which allows reconstructing a small window from the whole image space. In this way, efficient numerical technique is obtained, which enables measuring the diffraction efficiency of the PAS algorithm. The second was the 2D MFFT-AS, which allows reconstructing the full image space. It was employed for reconstructing the whole object and visualize the quality of the occlusion culling that was carried out with the PAS. Finally, optical reconstruction of occluded CGHs in wide-angle holographic display verified the validity of our occlusion culling strategy. Notably, our methodology enables fast computation time for obtaining occluded CGHs from very large point clouds.  Data Availability Statement: Data underlying the results presented in this paper are publicly available at this time but may be obtained from the authors upon request.