Occlusion handling in computer-generated holography is of vast importance as it enhances depth information by presenting correct motion parallax of the 3D scene within the viewing angle. In this paper, we propose a computationally efficient occlusion handling technique based on a fully analytic mesh based computer generated holography. The proposed technique uses angular spectrum convolution that renders exact occlusion while preserving all other aspects of the fully analytic mesh based computer generated holography. The proposed method is computationally efficient as only a single convolution operation is required for each mesh without numerical propagation between the meshes. The proposed method is also exact as it performs the occlusion processing in the tilted mesh plane, being free from artifacts coming from orthographic spatial masking. The proposed method can be applied to the self and the mutual occlusions between the objects in the 3D scene. The computer simulated results show the feasibility of the proposed method.
© 2017 Optical Society of America
Computer generated hologram (CGH) synthesizes the complex optical field of a 3D scene in a numerical way. The 3D scene is represented using different primitives such as a point [1–10], a light ray [11–15], a layer [16–18], and a triangular mesh [19–27], and the complex field for individual primitives are calculated and aggregated to form the complex field for the entire scene. The unique feature of the CGH over still images is that it provides not a single image at a specific viewpoint but a continuous motion parallax within a viewing angle. In order to fully acquire its benefit, the occlusion should be incorporated properly within the viewing angle as it is the most important factor in enhancing the depth information by producing correct motion parallax.
Several methods have been proposed to deal with the occlusion in the CGH. In the point cloud based methods where the 3D scene is sampled into a number of self-luminous points [1–10], the occlusion processing can be done by applying a visibility test that works in an iterative way of checking each point in the rear object to be blocked by a front object [5–8]. For better reconstruction quality and visual representation, the object must be sampled at a very high-density, which increases the computational cost significantly. Although the visibility check can be reduced by grouping pixels in the hologram plane [9-10], the computational cost is still high.
In the light ray based CGHs [11–15], the occlusion is handled in the light ray or stereoscopic view space before they are converted into the complex wave field. As the occlusion handling in the light ray or the stereoscopic view space has been well developed in computer graphics field, the occlusion can be implemented relatively easily in this type of CGHs. However, the light ray based CGHs have the limitation in the point that they reproduce not the exact wavefront of the 3D scene but the collection of the light rays or stereoscopic views which only approximate the wavefront.
In the layer based methods [16–18] and the triangular mesh based methods [25–27], the spatial masking techniques have been proposed for the occlusion processing. In the basic spatial masking techniques, the object space is divided into many planes parallel to the hologram. For each plane, the wavefield from the rear space is numerically propagated to the current plane and masked by the object slice or the orthographic silhouette of the triangular meshes associated with the current plane [16, 26]. The computational load of this method can be reduced by using the silhouette aperture instead of the silhouette mask and limiting the processing area . However, as the masking is performed not in the tilted mesh plane but in the plane parallel to the hologram, these techniques show artifacts at oblique observation angles. Although the spatial masking in the tilted mesh plane has been proposed , it requires much more computations.
In this paper, we propose a novel occlusion processing technique based on the fully-analytic triangular mesh based CGH. The fully-analytic triangular mesh based CGH technique [19–24] of this paper is differentiated from the previous triangular mesh based CGH techniques [25–27] which have been used for the spatial masking, in the way that it calculates the global angular spectrum of each triangular mesh. In the previous techniques, the local angular spectrum is first obtained in the local spatial frequency grid by taking the discrete Fourier transform in the tilted local plane of the triangle, then the global angular spectrum is obtained by resampling and interpolating the local angular spectrum, considering geometrical rotation and translation between the global hologram plane and the local plane. To the contrary, in the fully-analytic technique which is used in this paper, the global angular spectrum is obtained directly from the analytic formula of the angular spectrum of the reference triangle, without performing the discrete Fourier transform in the local plane. Therefore no resampling or interpolation of the local angular spectrum is required in the fully analytic technique, which ensures exact complex field calculation for the given sampling grid in the hologram plane.
The proposed method performs the occlusion processing not in the spatial domain like previous methods but in the angular spectrum domain. In the proposed method, the global angular spectrum of the rear meshes is convoluted with the global angular spectrum of the current mesh to determine the angular spectrum part that should be occluded by the current mesh. The convolution operation is performed in the hologram plane, not in the individual local mesh plane. Therefore, the numerical propagations between the silhouette mask planes are not required, which reduces the computational cost of the proposed method. The global angular spectrum convolution achieves the blocking of the rear wave field in the tilted local mesh plane, not in the plane parallel to the hologram. Therefore the proposed method is exact, being free from oblique angle artifacts of the previous silhouette methods.
In the following section, we briefly review the principle of the fully-analytic triangular mesh based CGH. In next sections, we present the principle of the proposed method and the necessary algorithms for its implementation. Finally, we show the numerical simulation results to verify the proposed method.
2. Brief introduction to fully analytic mesh-based CGH
The global angular spectrums of individual triangular meshes are calculated and added to the hologram plane using the analytic formula of the reference triangle. Finally, the aggregated global angular spectrum is Fourier transformed in the hologram plane to give the complex wave field.
2.1 Single plane carrier wave
Suppose that the (x,y,z) is the global coordinate where z = 0 is the hologram plane. The (xl,yl,zl) is the local coordinate where the triangle lies in zl = 0 plane and one of the vertices of the triangle coincides with the local coordinates origin (xl,yl,zl) = (0,0,0). The global and local coordinates are related by the 3 × 3 rotation matrix R and the 3 × 1 shift vector c such that rxl,yl,zl = Rrx,y,z + c where rx,y,z = [x,y,z]T and rxl,yl,zl = [xl,yl,zl]T are 3 × 1 global and local position vectors, respectively. A plane carrier wave u(rx,y,z) = aexp[j2πνx,y,zTrx,y,z] is illuminating the triangle where νx,y,z = [νx, νy, νz]T with νz = (1/λ2-νx2-νy2)1/2 is the 3 × 1 global spatial frequency vector and a is the initial complex amplitude of the plane carrier wave. The global angular spectrum AS(fx,y; νx,y; a) of the mesh in the hologram plane is given byEq. (1), it is assumed that we use the reference triangle having one of its three vertices at the coordinate origin. The 3 × 1 vectors fx,y,z = [fx, fy, fz]T, fxl,yl,zl = [fxl, fyl, fzl]T are spatial frequencies in the global and the local plane where fz = (1/λ2-fx2-fy2)1/2, fzl = (1/λ2-fxl2-fyl2)1/2, and fxl,yl,zl = Rfx,y,z. The fx,y is a global 2 × 1 spatial frequency vector defined by fx,y = [fx, fy]T. The νx,y = [νx, νy]T is the 2 × 1 spatial frequency vector of the plane carrier wave. The 2 × 2 matrix A is the transformation matrix relating the local triangle gl(rxl,yl) = gl([xl,yl]T) and the reference triangle gr by gl(rxl,yl) = gr(Arxl,yl). The 3 × 1 vector rox,y,z is the position vector in the global coordinates of the triangle vertex which is located at the local coordinates origin. In Eq. (1), the first term B(fx,y; νx,y) is given by
2.2 Multiple plane carrier waves
The angular spectrum AS(fx,y; νx,y; a) calculated using Eqs. (1) and (2) is concentrated around the carrier wave spatial frequency νx,y, giving narrow viewing angle around the carrier wave direction. In order to enlarge the viewing angle as required in the case of the diffuse surface, the angular spectrum is calculated for multiple carrier waves through accumulating the angular spectrums of the individual carrier waves byEq. (3) requires heavy computation, its reduction technique using angular spectrum convolution has been reported . In , it has been shown that B(fx,y;νx,y) can be approximated to21] indicates that this approximation causes slight mismatch of the spatial frequency grid on the reference angular spectrum ASr at large inclination of the triangular mesh with respect to the hologram plane and at large viewing angle. However, its impact on the final angular spectrum or complex field is negligible as the mismatch only happens where the reference angular spectrum ASr has small amplitude, which justifies this approximation . With this approximation, the angular spectrum for multiple carrier waves given in Eq. (3) can now be obtained from Eq. (1) by21].
Finally, the angular spectrums of the individual triangular meshes are aggregated and the final complex field in the hologram plane is then obtained by
3. Proposed method
Figure 2 illustrates the proposed occlusion handling method. First, all triangular meshes in the 3D scene are sorted in the depth order from rear to front. The angular spectrum calculation and the occlusion handling is performed sequentially from the first (rear) triangle to the last (front) triangle. Suppose that the angular spectrum of the wave field from the first n triangles is given by AS1:n(fx,y). The angular spectrum including the next, i.e. (n + 1)-th triangle can be written by
The AS1:n,n+1(fx,y) is the angular spectrum of the rear wave field part which passes through the (n + 1)-th triangle area. We can think that the AS1:n,n+1(fx,y) is the angular spectrum of the clear aperture which has the same shape as the (n + 1)-th triangle and is illuminated by the rear wave field. The rear wave field can be represented using the collection of the plane waves and their complex amplitudes in the hologram plane are the angular spectrum of the rear wave field AS1:n(fx,y), which is already given. Therefore, the AS1:n,n+1(fx,y) can be obtained using Eq. (3) by
The angular spectrum convolution technique explained in the previous section 2.2 can be applied here as well. From Eqs. (5) and (6), the AS1:n,n+1(fx,y) can be written byEq. (11) and all variables in Eq. (12) correspond to the (n + 1)-th triangle. Therefore, the angular spectrum of the rear wave field part which should be occluded by the (n + 1)-th triangle can simply be obtained by taking a single convolution between the (n + 1)-th triangle term and the angular spectrum of the rear wave field.
The computational load can be more reduced by using the fact that many terms are shared in AS1:n,n+1(fx,y) and ASn+1(fx,y). In Eq. (8), ASn+1(fx,y) - AS1:n,n+1(fx,y) can be combined intoEq. (8), i.e. the blocking of the occluded rear wave field and the addition of the new mesh, can be obtained by a single convolution operation.
Note that the additional processing for the occlusion handling over the conventional fully-analytic mesh based CGH with target reflectance distribution is only the replacement of Dn+1(fx,y) with Dn+1(fx,y) - D1:n(fx,y) in Eq. (13). Unlike previous spatial masking techniques, numerical propagations between the mesh planes are not required and all processing is performed in the angular spectrum domain in the hologram plane. These two features make the proposed method computationally efficient.
The whole procedure maintains the fully-analytic framework as it requires no resampling in the local coordinates. The proposed angular spectrum convolution achieves the occlusion in the local plane which is arbitrarily slanted with respect to the hologram plane. Therefore the proposed method is free from artifacts of the conventional orthographic silhouette method which performs the spatial masking in the transverse plane.
Note that Eq. (13) is valid when the reference triangle is a clear aperture having unity inside and zero outside. Therefore, in the techniques using the reference triangle of non-uniform amplitude, for example, the continuous shading technique , the simplification using Eq. (13) cannot be applied. In these cases, however, AS1:n,n+1(fx,y) and ASn+1(fx,y) can still be calculated using single convolution, respectively. Therefore, the computational cost increase is manageable.
The sequential processing for individual meshes in the depth order is presented in this section. The depth ordering of the meshes, however, is sometimes ambiguous due to complicated arrangement of the meshes. In these cases, we can group the meshes in similar depth and apply the same rear wave field to them. Although the occlusion between the meshes in the same group is not realized in this case, it is rarely a problem as the meshes in the group have the similar depth, not occluding each other within large viewing angle.
4. Simulation Results
The proposed method was verified using numerical simulations. In all numerical simulations, the complex valued holograms H(rx,y) are used without additional encoding process to amplitude only or phase only holograms. Figure 3 shows the numerical results of each step of the proposed method. In Fig. 3, two triangular meshes at different depths, i.e. 20 mm and 24 mm from the hologram plane are used as a scene. The hologram of the two triangles was calculated over 3000 × 3000 regular sampling grid with 2 × 2 um sampling pitch for 532 nm wavelength. For a clear illustration of the angular spectrum, only a single plane carrier wave normal to the hologram plane is used as the carrier wave. In Fig. 3, the upper row shows the amplitude part of the angular spectrums in the hologram plane and the lower row shows their numerical propagation results to the front triangle plane at z = 24 mm. The wave field from the rear triangle whose numerical reconstruction at the front triangle plane is blurred as shown in Fig. 3(a) is effectively blocked by the front triangle as shown in Figs. 3(b) and 3(c). By adding the wave field for the front triangle shown in Fig. 3(d), the final hologram with correct occlusion is synthesized as shown in Fig. 3(e). Note that the conventional method which simply adds AS1(fx,y) and AS2(fx,y) shows overlapping of two triangles as shown in Fig. 3(f). Figure 4 shows another example of a mug object having 3450 triangular meshes. Similarly, the overlapping between the triangular meshes in the rear and front surface of the mug is removed by the proposed method.
Figure 5 shows numerical observations of the two-triangle hologram from different directions. The parameters of the two triangles and the hologram are the same as those of Fig. 3. But in this case, we use multiple plane carrier waves with random phases to represent the diffuse surface. For numerical observation of the hologram, we place a lens at different transverse positions in z = 36 mm plane and the intensity images in the conjugate plane of the front triangle are calculated. As shown in Fig. 5 and the associated movie, it can be confirmed that the proposed method realizes continuous and correct occlusion within the viewing angle of the hologram.
Figure 6 is another numerical observation from different directions of the hologram of two triangles. In this case, the front triangle is tilted such that it is not parallel to the hologram plane. Figure 6 shows the result when we apply the spatial masking method where the wave field from the rear triangle is blocked by the orthographic projection of the tilted front triangle following the conventional silhouette method and the result of the proposed method. From Fig. 6, it can be seen that the overlapping and black line artifacts of the orthographic spatial masking approach are corrected in the proposed method as highlighted by the white box. This result ensures the proposed method handles the occlusion in the tilted plane properly without any artifact.
Figures 7 and 8 show the numerical observations of the complex torus objects having 1600 triangular meshes in the scene. The hologram is calculated over 3000 × 3000 sampling grid with 2 μm sampling pitch. The wavelength is 532 nm. As mentioned in the last paragraph of section 3, the meshes of the scene are divided into 16 groups according to their depths and the same rear wave field is used for every mesh in the same group to avoid depth order ambiguity between the meshes of the similar depths. Figure 7 shows the case when only a single plane carrier wave is used to highlight the overlapping artifact due to lack of the occlusion processing. Figure 8 shows numerical observations of the proposed method with multiple carrier waves. It can be seen that the occlusion is realized properly within the full viewing angle of the hologram.
Finally, the computation time for the numerical simulations is shown in Table 1, comparing the conventinoal method without occlusion processing, the sihllouette switch back method , and the proposed method. All calculations were done using Matlab software for diffusive polygon surfaces. In order to reduce the computation time, the computation window reduction technique [26,27] was applied to our code implementation. The analytic calculation of the angular spectrum and the discrete fast Fourier transform (FFT) are performed not in the full 3K × 3K window (or even larger window with zero padding in case of FFT to avoid possible edge artifact originating from the periodicity of the FFT) but in the small computation window around each triangular mesh considering the diffraction angle . The size of the reduced computation window was kept the same for all three methods under comparison.
Note that when the computation window reduction technique is not applied, all intermediate data AS1:n(fx,y) is kept in the angular spectrum domain and the proposed occlusion processing method requires 1 convolution operation with Dn+1(fx,y) - D1:n(fx,y) in Eq. (13), or 3 FFTs in full 3K × 3K window for each mesh to process the occlusion and the reflectance distribution. Since the conventional method without the proposed occlusion processing also requires 1 convolution with Dn+1(fx,y) for the reflectance distribution processing, there is no significant increase in the computaional load in the proposed method. By applying the computation window reduction technique, however, the intermediate data should be kept in the spatial domain H1:n(rx,y) for extracting the appropriate computation window around the mesh. This requires 2 additional FFTs for each mesh in the proposed method, i.e. one for the preparation of D1:n(fx,y) from H1:n(rx,y) and one for the transformation of the resultant AS1:n + 1(fx,y) to H1:n + 1(rx,y). In the case of the conventional method without the proposed occlusion processing, only one additional FFT is required for each mesh to transform the angular spectrum of each mesh AS(fx,y) to its spatial complex field H(rx,y). Therefore, the compution overhead by the proposed method can be slightly larger in the computation window reduction technique.
Also note that the original silhouette switch back method  is based on the non-analytic mesh based CGH where the local angular spectrum of each mesh is obtained by the FFT. But in our implementation, we applied the sihllouette switch back method to the fully-analytic mesh based CGH like other two methods under comparison, i.e. the conventional method without occlusion processing and the proposed method. Therefore in our implementaion, the proposed method and the silhouette switch back method  generate the angular specturm of each mesh in the same way, but processes the occlusion differently.
From Table 1, it can be seen that the increase in computation time due to the occlusion processing of the proposed method is small but not negligible. This increase in the computation time is believed to originate from slightly larger computatioanl overhead when the computation window reduction technique is applied to the proposed method as explained above. But the computation time of the proposed method is still smaller than the silhouette switch back method  in our implementation. It is because the silhouette switch back method uses additional numerical propagations to obtain the rear complex field that passes through the aperture of the current mesh, which increases required number of the FFTs.
In this paper, we proposed an angular spectrum convolution technique for accurate self or mutual occlusion processing in a 3D scene based on the fully analytic triangular mesh based CGH. The method involves convolution between the angular spectrum of the rear wave field and the current triangular aperture. This convolution happens in hologram plane with global spatial frequency grid and thus it requires no additional resampling, maintaining the benefit of fully analytic mesh based CGH. The proposed method outperforms the orthographic spatial masking technique in terms of computation cost and occlusion accuracy at oblique observation angles. The principle of the proposed method was verified by numerical reconstructions, ascertaining full parallax with correct occlusion handling in a 3D scene.
Basic Science Research Program, National Research Foundation of Korea (NRF) [NRF-2017R1A2B2011084]; ITRC support program, MSIT [IITP-2017-2015-0-00448];
References and links
1. M. Lucente, “Interactive computation of holograms using a look-up table,” J. Electron. Imaging 2(1), 28–34 (1993).
2. B. R. Brown and A. W. Lohmann, “Complex spatial filtering with binary masks,” Appl. Opt. 5(6), 967–969 (1966). [PubMed]
3. H. Zhang, Q. Tan, and G. Jin, “Holographic display system of a three- dimensional image with distortion-free magnification and zero-order elimination,” Opt. Eng. 51(7), 075801 (2012).
4. Y. Ogihara and Y. Sakamoto, “Fast calculation method of a CGH for a patch model using a point-based method,” Appl. Opt. 54(1), A76–A83 (2015). [PubMed]
5. J. S. Underkoffler, “Occlusion processing and smooth surface shading for fully computed synthetic holography,” Proc. SPIE 3011, 19–30 (1997).
6. R. H.-Y. Chen and T. D. Wilkinson, “Computer generated hologram with geometric occlusion using GPU-accelerated depth buffer rasterization for three-dimensional display,” Appl. Opt. 48(21), 4246–4255 (2009). [PubMed]
7. H. Zhang, N. Collings, J. Chen, B. Crossland, D. Chu, and J. Xie, “Full parallax three-dimensional display with occlusion effect using computer generated hologram,” Opt. Eng. 50(7), 074003 (2011).
8. H. Zhang, Q. Tan, and G. Jin, “Full parallax three-dimensional computer generated hologram with occlusion effect using ray casting technique,” J. Phys. Conf. Ser. 415, 012048 (2013).
9. R. H.-Y. Chen and T. D. Wilkinson, “Computer generated hologram from point cloud using graphics processor,” Appl. Opt. 48(36), 6841–6850 (2009). [PubMed]
10. I. Hanák, M. Janda, and V. Skala, “Detail-driven digital hologram generation,” Vis. Comput. J. 26(2), 83–96 (2010).
11. T. Yatagai, “Stereoscopic approach to 3-D display using computer-generated holograms,” Appl. Opt. 15(11), 2722–2729 (1976). [PubMed]
12. M. Yamaguchi, H. Hoshino, T. Honda, and N. Ohyama, “Phase-added stereogram: calculation of hologram using computer graphics technique,” Proc. SPIE 1914, 25–33 (1993).
13. Q. Y. J. Smithwick, J. Baradas, D. E. Smalley, and V. M. Bove, “Real-time shader rendering of holographic stereograms,” Proc. SPIE 7233, 723302 (2009).
14. K. Wakunami and M. Yamaguchi, “Calculation for computer generated hologram using ray-sampling plane,” Opt. Express 19(10), 9086–9101 (2011). [PubMed]
15. K. Wakunami, H. Yamashita, and M. Yamaguchi, “Occlusion culling for computer generated hologram based on ray-wavefront conversion,” Opt. Express 21(19), 21811–21822 (2013). [PubMed]
16. H. Zhang, L. Cao, and G. Jin, “Computer-generated hologram with occlusion effect using layer-based processing,” Appl. Opt. 56(13), F138–F143 (2017). [PubMed]
17. N. Okada, T. Shimobaba, Y. Ichihashi, R. Oi, K. Yamamoto, M. Oikawa, T. Kakue, N. Masuda, and T. Ito, “Band-limited double-step Fresnel diffraction and its application to computer-generated holograms,” Opt. Express 21(7), 9192–9197 (2013). [PubMed]
18. Y. Zhao, L. Cao, H. Zhang, D. Kong, and G. Jin, “Accurate calculation of computer-generated holograms using angular-spectrum layer-oriented method,” Opt. Express 23(20), 25440–25449 (2015). [PubMed]
19. L. Ahrenberg, P. Benzie, M. Magnor, and J. Watson, “Computer generated holograms from three dimensional meshes using an analytic light transport model,” Appl. Opt. 47(10), 1567–1574 (2008). [PubMed]
20. Y.-M. Ji, H. Yeom, and J.-H. Park, “Efficient texture mapping by adaptive mesh division in mesh-based computer generated hologram,” Opt. Express 24(24), 28154–28169 (2016). [PubMed]
21. H.-J. Yeom and J.-H. Park, “Calculation of reflectance distribution using angular spectrum convolution in mesh-based computer generated hologram,” Opt. Express 24(17), 19801–19813 (2016). [PubMed]
22. H. Kim, J. Hahn, and B. Lee, “Mathematical modeling of triangle-mesh-modeled three-dimensional surface objects for digital holography,” Appl. Opt. 47(19), D117–D127 (2008). [PubMed]
23. J.-H. Park, H.-J. Yeom, H.-J. Kim, H. Zhang, B. Li, Y.-M. Ji, and S.-H. Kim, “Removal of line artifacts on mesh boundary in computer generated hologram by mesh phase matching,” Opt. Express 23(6), 8006–8013 (2015). [PubMed]
24. J.-H. Park, S.-B. Kim, H.-J. Yeom, H.-J. Kim, H. Zhang, B. Li, Y.-M. Ji, S.-H. Kim, and S.-B. Ko, “Continuous shading and its fast update in fully analytic triangular-mesh-based computer generated hologram,” Opt. Express 23(26), 33893–33901 (2015). [PubMed]
25. K. Matsushima, “Exact hidden-surface removal in digitally synthetic full-parallax holograms,” Proc. SPIE 5742, 25–32 (2005).
26. K. Matsushima and S. Nakahara, “Extremely high-definition full-parallax computer-generated hologram created by the polygon-based method,” Appl. Opt. 48(34), H54–H63 (2009). [PubMed]
27. K. Matsushima, M. Nakamura, and S. Nakahara, “Silhouette method for hidden surface removal in computer holography and its acceleration using the switch-back technique,” Opt. Express 22(20), 24450–24465 (2014). [PubMed]