# Real-Time Calculation of 3D Computer-Generated Holograms by Phase Space Warping

Antoine Lagrange, Antonin Gilles, Patrick Gioia, Vincent Savaux, Kevin Heggarty, Bruno Fracasso

► To cite this version:

Antoine Lagrange, Antonin Gilles, Patrick Gioia, Vincent Savaux, Kevin Heggarty, et al. Real-Time Calculation of 3D Computer-Generated Holograms by Phase Space Warping. 2026. hal-05524264

HAL Id: hal-05524264

https://hal.science/hal-05524264v1

Preprint submitted on 23 Feb 2026

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Real-Time Calculation of 3D Computer-Generated Holograms by
Phase Space Warping

Antoine Lagrange, Antonin Gilles, Patrick Gioia, Vincent Savaux, Kevin Heggarty, Bruno Fracasso

***Abstract***—Modeling light propagation from a 3D structure to a plane in a medium with a constant refractive index, according to scalar diffraction theory, is essential in many optical fields, particularly in computer-generated holography for 3D visualization. However, this process is computationally expensive, often scaling with the geometry complexity of the 3D scene or requiring compromises in computational rigor. To overcome this issue, we propose to model this phenomenon in phase space, which represents both spatial and spatial frequency information in a joint domain. By decomposing wavefields into sets of Gaussian beams, our formulation provides a more precise control over light propagation, significantly reducing the dependency on scene geometry. We implemented our method in CUDA, leveraging acceleration with GPU. Our experiments reported up to a 26.1-fold speedup compared to the conventional layer-based methods, while maintaining a similarity of 35.1 dB in PSNR.

***Index Terms***—Computer-Generated Holograms, Phase Space, Gabor Transform, Gaussian beams, Diffraction Modelling

## I. INTRODUCTION

MANY fields of research in optics, such as computer-generated holography, rely on scalar diffraction theory to digitally model the propagation of light in a medium with a constant refractive index. More precisely, a central result of this theory is the Rayleigh-Sommerfeld diffraction formula, which expresses any wavefield point in space in terms of the values of the same wavefield on a given surface [1]. However, direct evaluation of this integral is computationally expensive. To improve efficiency, the angular spectrum method propagates the wavefield from an initial planar surface to a parallel plane at a normal distance *z* using only two Fast Fourier Transforms (FFTs). Under paraxial approximation, the computation time can be further reduced by using the Fresnel diffraction integral. While these methods provide an efficient means to model light propagation from a plane, extensions to 3D structures with complex geometries typically entail a substantial increase in computational cost or a loss of modeling rigor. Overcoming this limitation would be highly beneficial for computer-generated holography, particularly in 3D visualization applications, where the trade-off between computation time and visual quality remains prohibitive for practical use.

In recent years, deep learning methods have been extensively studied for hologram calculation in 3D visualization

applications [2], including both the modeling of light propagation and the quantization required to meet display constraints. Although these optimization processes are particularly well suited for minimizing the loss of information induced by quantization [3], modeling light propagation itself is not an inverse problem. Hence, an efficient analytical method should be preferred to avoid the intrinsic issues of deep learning such as training cost and spectral bias.

Such analytical formulations include primitive-based methods that calculate Computer-Generated Holograms (CGHs) by first dividing the 3D model into primitives, such as points [4] or polygons [5]. Light is then propagated from each primitive to a hologram plane using scalar diffraction theory. The resulting hologram corresponds to the coherent sum of all these contributions. Despite rigorously modeling light propagation, these methods have an algorithmic complexity that scales with scene geometry, making the computational cost prohibitive for complex 3D models. Approaches such as split-Lohmann [6] or Holographic Stereogram (HS) methods avoid this scaling but compromise rigorous light propagation, limiting depth perception. Recent methods address this issue, yet most either do not provide all depth cues due to the division of the hologram into hogels [7], or fail to meet real-time requirements [8].

To model diffraction both efficiently and rigorously, we have developed an analytical method formulated in phase space. For light waves propagating along the *z*-axis, this joint spatial and spatial frequency space consists of four dimensions: the spatial components representing position in the transverse plane, and the spatial frequency coordinates representing the wave vector components in the transverse directions. Hence, in this domain, a more precise control of light wavefront is possible, since wavefields are decomposed into beams that can be manipulated individually. Previous work on this topic was limited to diffraction between two planes, including rigorous calculation [9], [10] and paraxial approximation [11], [12]. In this paper, we extend these methods to model light propagation from a 3D model to a plane in a medium with a constant refractive index.

Specifically, we make the following contributions:

* We generalize previous formulas that describe diffraction between two planes in phase space, including a covariance formula, to handle 3D scenes.

* We develop methods for working around the uncertainty principle to increase the accuracy with which light can be controlled by relying on a depth-dependent linear phase and a band-limited phase as well as a layered depth map in phase space.

Antoine Lagrange, Antonin Gilles and Vincent Savaux are with Research and Technology Institute b<>com, 35510 Cesson-Sévigné, France. Patrick Gioia is with Orange Innovation, 35510 Cesson-Sévigné, France. Antoine Lagrange, Kevin Heggarty and Bruno Fracasso are with IMT Atlantique, Technopole Brest-Iroise, CS 83818, 29285 Brest, France.

Manuscript received April 19, 2021; revised August 16, 2021.
* We implement our modeling method in CUDA to demonstrate practical usage, including the Gabor transform.

* We assess the effect of interpolation over fast varying coefficients in phase space, which is crucial for our modeling method to be efficient.

* We evaluate the speedup and similarity of our modeling approach compared to the layer-based approach.

## II. RELATED WORK

### A. *Primitive-based Methods*

In primitive-based methods, the 3D scene is decomposed into a set of primitives from which light propagation is computed according to scalar diffraction theory. The hologram is then obtained by summing all diffraction contributions. Depending on the chosen type of primitive, different trade-offs can be achieved between computational cost and the accuracy of scene representation. However, regardless of the approach, the computation time scales linearly with the number of primitives and thus with the complexity of the scene.

The first types of primitives introduced were points [13], polygons [14], and layers [15]. Since these initial formulations, improvements have been proposed to further mitigate the trade-off between computational cost and visual quality. Recent methods often exploit the sparse representation of the wavefield in the spatial frequency domain or phase space [16], [17], where most information is concentrated in a small number of significant coefficients. As a result, fewer operations are required to compute a hologram while preserving most of the relevant information. Other approaches focus on enhancing visual realism by more accurately modeling optical effects [18], lighting [19], or occlusions [20].

More recently, two formulations introducing new primitive types, lines [21] and Gaussians [22], have been proposed. The first method efficiently models light propagation from 3D wireframes by integrating point spread functions along line paths. The second relies on 2D Gaussian splatting, where a 3D scene is represented as a set of scaled and rotated 2D Gaussians. Light propagation is modeled from tilted planes similarly to polygon-based approaches, although the intensity is weighted by Gaussian functions. View-dependent effects are handled through an alpha blending strategy inspired by conventional computer graphics.

### B. *Light-Field Methods*

The first approach to exploit light-field information for modeling light propagation from a 3D scene is the HS technique [23]. In this method, the hologram is divided into spatial blocks called hogels, within which the light wave emitted from a single object point is approximated as a plane wave. Compared to primitive-based methods, the HS technique has an algorithmic complexity independent of scene geometry and inherently handles view-dependent effects, which reduces computation time. However, the formulas from scalar diffraction theory are only approximated, which limits depth perception.

To address this limitation, subsequent methods have aimed to reproduce scalar diffraction theory more faithfully by incorporating depth information into the calculations while keeping computational costs low. Phase-Added Stereogram (PAS) techniques [24] introduce a quadratic depth-dependent phase factor to account for wavefront curvature and to reduce incoherencies between hogels. Other strategies use point-based [7] or layer-based [25] approaches to compute each hogel, but computation time often increases with both the extent and the discretization of the depth range. Despite these improvements, all methods that decompose the hologram into hogels face an inherent trade-off between spatial and spatial-frequency resolution.

To suppress this trade-off, non-hogel methods [8], [26], that invert the light-field with a phase space transform, were proposed. However, light-fields only provide amplitude information, while phase is lost, which makes this approach ill-posed. To retrieve part of the phase information, Padmanaban et al. [27] proposed to multiply the light-field by the PAS phase factor. However, this method struggle to meet real-time requirements. Another method [28] addresses the trade-off by reproducing an additive compressive light-field display via a computationally expensive optimization strategy.

### C. *Single-Step Methods*

To develop approaches whose algorithmic complexity is independent of scene geometry, single-step methods concentrate depth-specific information at separate locations in the spatial or spatial-frequency domain. All depth layers are thus processed independently yet simultaneously. However, as the number of depth layers increases, these methods approximate scalar diffraction formulas more coarsely, degrading depth perception.

The first single-step method [29] relies on non-uniform sampling of the depth layers and requires depth-specific information to be localized in the hologram plane. However, this condition implies a reduced depth range for the 3D scene. More recently, the split-Lohmann method [6], which digitally reproduces a multifocal display constituted by a Lohmann lens, was proposed. This method was then extended to Fresnel holography [30] and enhanced to eliminate lateral shift artifacts [31]. In split-Lohmann methods, depth-specific information is concentrated in different spatial-frequency areas. However, to guarantee sufficient localization in this domain, constraints are imposed on the depth range of the 3D scene.

## III. PRELIMINARIES ON PHASE SPACE

### A. *Physical Interpretation*

In the context of wave propagation, phase space consists of two spatial dimensions and two spatial frequency dimensions, enabling spatially local description of spatial frequencies. Each Phase Space Coefficient (PSC) then indicates the contribution of a wavefield in a given direction and at a given spatial position.

To give more insight about the physical interpretation of phase space, the PSCs of a point source are analyzed in different configurations, see [32] for an in-depth study. Without any loss of generality, a 1D cross-section of a wavefield is
<figure><img src="image_2.png" /><figcaption>Fig. 1. Point source PSCs. Each image represents the PSCs of a 1D cross-section of a hologram generated from a single point source located on the optical axis. In the left image, the point source lies in the same plane as the wavefield, thus all spatial frequencies are associated with the same spatial position. For the two other images, the point source is located at a smaller z coordinate than the wavefield plane. Hence, the PSCs spread over the spatial domain, with different behavior in the rigorous and paraxial cases.</figcaption></figure>

analyzed. Therefore, phase space is reduced to two dimensions, where the spatial and spatial frequency coordinates are denoted $\tau_{\xi}$ and $\nu_{\xi}$, respectively. If the point source is located in the hologram plane at the origin, the PSCs of the hologram correspond to the left image in Figure 1. The information is concentrated at the origin along the spatial axis and spread over the entire spatial frequency extent. This represents a point source which is localized in space and emits light in every direction.

Now consider that the point source is located on the optical axis and that light propagates in the positive z direction over a distance d to the hologram. The corresponding PSCs, without and with paraxial approximation, are shown in the right-hand part of Figure 1. In both cases, the point spread function becomes more widely distributed across the hologram plane. In addition, since light is emitted from a single point, each spatial position is solely associated to a given spatial frequency interval. To better understand the distribution of these PSCs, the local spatial frequencies of the corresponding holograms are computed from the kernel functions of the Rayleigh-Sommerfeld diffraction formula and of the Fresnel diffraction integral. In the rigorous case, the local spatial frequencies are expressed as

$$f_{\xi}(r) = \frac{1}{2\pi} \frac{\partial kr}{\partial \xi} = \frac{\xi - x}{\lambda r}, \qquad (1)$$

where $\xi$ is the hologram coordinate, $x$ is the point source coordinate, and

$$r = \sqrt{(\xi - x)^2 + d^2}.$$

Similarly, with the paraxial approximation, the local spatial frequencies are

$$f_{\xi}(r) = \frac{1}{2\pi} \frac{\partial k \left( d + \frac{(\xi - x)^2}{2d} \right)}{\partial \xi} = \frac{\xi - x}{\lambda d}. \qquad (2)$$

These results explain why in the rigorous case the local spatial frequencies converge towards a given value, while with the paraxial approximation the local spatial frequencies rise linearly. Nevertheless, the local spatial frequency only represents the average spatial frequency at a given position, whereas the wave nature of light inherently produces a band of spatial frequencies at each spatial point.

## B. Discrete Gabor Transform

The PSCs of a wavefield can be computed using various phase space transforms, each employing distinct basis functions for decomposition. For instance, the Gabor transform is a specific case of the Short Time Fourier Transform (STFT) in which the window function is a Gaussian. Thus, the Heisenberg uncertainty principle is minimized, achieving the best possible joint localization in both the spatial and the spatial frequency domains. In the discrete case, this transform is expressed as follows

$$\mathcal{V}_u(\mathbf{w}_x) = \sum_{x,y} u(x,y) \overline{\psi(x - \tau_x, y - \tau_y)} e^{-j(x\nu_x + y\nu_y)}, \quad (3)$$

where $\psi$ is the Gaussian function, and $\mathcal{V}_u$ is the PSCs of $u$, where the coordinates in the spatial domain and phase space are $(x, y)$ and $(\tau_x, \tau_y, \nu_x, \nu_y)$, respectively. For brevity, the notation $\mathbf{w}_x = (\tau_x, \tau_y, \nu_x, \nu_y)$ is used in the following.

The Gabor transform relies on the concept of frame, which is a generalization of a basis in a Hilbert space $\mathcal{H}$. While a basis provides a unique decomposition of any element of this space, a frame allows for possibly redundant decompositions, still ensuring stable reconstruction. A formal introduction to frame theory is provided in [33]. For each $\mathbf{w}_x$, the space-frequency shift operator is defined as follows

$$\pi(\mathbf{w}_x)\psi(x, y) = e^{j2\pi(x\nu_x + y\nu_y)}\psi(x - \tau_x, y - \tau_y). \quad (4)$$

The Gabor system is then

$$\mathcal{G}_{\psi} = \{\pi(\mathbf{w}_x)\psi | \mathbf{w}_x \in \Lambda\}, \qquad (5)$$

where $\Lambda$ is a 4D lattice. This system is called a Gabor frame for $\mathcal{H}$ if there exist constants $A, B > 0$ that satisfy

$$\forall u \in \mathcal{H}, \quad A \|u\|^2 \le \sum_{\mathbf{w}_x \in \Lambda} |\langle u, \pi(\mathbf{w}_x)\psi \rangle|^2 \le B \|u\|^2. \quad (6)$$

If this inequality holds, which can be ensured by choosing a sufficiently dense lattice, the dual frame $\mathcal{G}_{\tilde{\psi}}$ is defined such that

$$\begin{aligned} \forall u \in \mathcal{H}, \quad & u = \sum_{\mathbf{w}_x \in \Lambda} \langle u, \pi(\mathbf{w}_x)\psi \rangle \pi(\mathbf{w}_x)\tilde{\psi} \\ & = \sum_{\mathbf{w}_x \in \Lambda} \langle u, \pi(\mathbf{w}_x)\tilde{\psi} \rangle \pi(\mathbf{w}_x)\psi. \end{aligned} \qquad (7)$$
If the direct Gabor transform is performed with $G_ψ$, then $G_ψ$ is employed to reconstruct $u$, and vice versa. The redundancy $R$ of a Gabor frame quantifies overcompleteness: the more overcomplete the frame, the higher the redundancy. This quantity is defined as the ratio of the number of lattice points to the spatial domain dimension.

A fast Gabor transform is essential for modeling efficiently diffraction in phase space without compromising the accuracy of the calculations. Hence, we implemented an optimized Gabor transform based on two key elements. First, we compute the transform on a regular lattice following the algorithm introduced in [34], which allows us to exploit separability along the spatial axes. Secondly, we developed a parallel implementation with CUDA to leverage hardware acceleration.

### C. 2D Diffraction

Since $u$ can be expressed in terms of the Gabor frame as in Equation (7), the effect of light propagation is studied through the warping of this family of functions. More formally, understanding how a propagation operator $\mathcal{T}$, such as the angular spectrum or the Fresnel diffraction integral, acts on the Gabor frame is essential.

The Egorov theorem, supported by a solid theoretical background, can be used to investigate this warping as in [10]. Instead, we adopt a simpler method based on a covariance formula derived for any translation and rotation between two planes [9], but here restricted to a translation along the optical axis

$$ \mathcal{T}\pi(\mathbf{w}_x)\psi \approx e^{j2\pi d\sqrt{\frac{1}{\lambda_x^2} - \nu_x^2 - \nu_y^2}} \pi \circ \chi(\mathbf{w}_x)\psi, \quad (8) $$

where $d$ is the normal distance between the two parallel planes. The symplectomorphism $\chi$ associated to the propagation operator $\mathcal{T}$ is a transformation of spatial and spatial frequency coordinates that preserves the symplectic structure and Hamiltonian flow, see [35] for more information on this topic. This function is expressed as follows

$$ \begin{pmatrix} \tau_{\xi} \\ \tau_{\nu} \\ \nu_{\xi} \\ \nu_{\nu} \end{pmatrix} = \chi \begin{pmatrix} \tau_{x} \\ \tau_{y} \\ \nu_{x} \\ \nu_{y} \end{pmatrix} = \begin{pmatrix} \tau_{x} + \frac{\nu_{x} d}{\sqrt{\frac{1}{\lambda_{x}^{2}} - \nu_{x}^{2} - \nu_{y}^{2}}} \\ \tau_{y} + \frac{\nu_{y} d}{\sqrt{\frac{1}{\lambda_{x}^{2}} - \nu_{x}^{2} - \nu_{y}^{2}}} \\ \nu_{x} \\ \nu_{y} \end{pmatrix}. \quad (9) $$

The warping of the Gabor frame is defined by both a phase space mapping corresponding to the symplectomorphism $\chi$ and a multiplication by the phase factor, $e^{j2\pi d\sqrt{\frac{1}{\lambda_x^2} - \nu_x^2 - \nu_y^2}}$. These operations can be interpreted physically as the symplectomorphism accounting for the displacements of light beams upon propagation and the phase factor accumulating the phase on these paths.

After warping phase space according to Equation (8), the separability of the inverse Gabor transform can no longer be exploited because the warped coefficients lie on an irregular lattice. To address this issue, the coefficients are first arranged in the $w_\xi$-space on a regular grid, which is then subjected to an inverse coordinate mapping. The output coefficients are subsequently obtained by interpolating over the regular lattice in $w_x$-space, as illustrated in the first row of Figure 2.

## IV. 3D DIFFRACTION IN PHASE SPACE

### A. Phase Space Warping

Previous methods only model diffraction in phase space between two planes, we propose to extend these formulations to account for light propagation from multiple depth layers, i.e. to reproduce a layer-based approach. First, PSCs are computed from an all-in-focus wavefield, including all depth layers, with a Gabor transform. Then, phase space is warped by considering that the normal distance over which light beams propagate varies according to the spatial area from which each beam originates. This phenomenon is accounted for by replacing the constant depth in Equations (8) and (9) by a Phase Space Depth Map (PSDM) $d_{w_\xi}$. Finally, the hologram is obtained by applying an inverse Gabor transform on the resulting Gabor spectrum.

The PSDM is captured in $w_\xi$-space to compute inverse mappings, which are followed by interpolation in $w_x$-space, see the second row of Figure 2. This process is carried out by using a ray-casting technique. Virtual cameras are positioned in the hologram plane at each position $(\tau_\xi, \tau_\eta)$ and cast rays in directions related to the spatial frequencies $(\nu_\xi, \nu_\eta)$. The depths at which these rays intersect a layer then constitute $d_{w_\xi}$.

Nevertheless, due to the uncertainty principle, the effective support of the Gaussian functions in the spatial and spatial frequency domains cannot both be arbitrarily small. This principle causes two main limitations that impair the accuracy of diffraction modeling. First, the window in Equation (3) may spread over several depth layers when computing the PSCs of the all-in-focus wavefield. Hence, objects should be well-separated along the $x$ and $y$ axes to avoid information from distinct depths to be merged in phase space. Second, the Gaussian functions used in the inverse Gabor transform also spread in both the spatial and spatial frequency domains. Hence, light beams emitted from different depth layers may contribute to the same coefficient in the $w_\xi$-space. Consequently, when considering a single depth value at each position in $w_\xi$-space, portions of the wavefield are neglected, which results in an incorrect propagation of light. Two workarounds are presented in the next section to address these issues.

### B. Nonzero Effective Support

The two issues identified in the previous section concern the extent of the effective support, in the spatial and spatial frequency domains, of the functions constituting the Gabor and dual frames. To address these problems, this extent is first analyzed in the spatial domain. Since $\psi$ is $L^2$-normalized, the ratio of the energy contained by the effective support on the total energy of the Gaussian is defined as follows

$$ \begin{align} \rho_x &= \frac{2}{\alpha} \int_0^{2\pi} \int_0^{l_x} re^{-\frac{2\pi r^2}{\alpha}} dr d\theta \tag{10} \\ &= 1 - e^{-\frac{2\pi l_x^2}{\alpha}}, \nonumber \end{align} $$

where the $l_x$ defines the integration area, $\rho_x \in [0, 1]$ corresponds to a percentage, and $\alpha$ controls the spread of the
<figure><img src="image_3.png" /><figcaption>Fig. 2. **Phase space warping.** In the first row, phase space is warped assuming a single depth layer, thus the displacement of the PSCs along the spatial axes is only a function of the spatial frequency components. In the second row, coefficients at given spatial frequency coordinates undergo different mappings according to the depth of the different layers. Both rows illustrate the inverse mapping and the interpolation used to compute direct and inverse Gabor transforms on regular grids.</figcaption></figure>

Gaussian function. A similar analysis can be conducted on the Fourier transform of the Gaussian window

$$ \hat{g}(f_x, f_y) = \sqrt{2\alpha e^{-\pi\alpha(f_x^2 + f_y^2)}} \qquad (11) $$

to obtain the extent of the effective support in the spatial frequency domain. Equation (10) is rewritten as follows

$$ \rho_{f_x} = 1 - e^{-2\pi\alpha l_{f_x}^2} \qquad (12) $$

where $l_{f_x}$ and $\rho_{f_x}$ are the counterparts of $l_x$ and $\rho_x$ in the spatial frequency domain, respectively. When the same lattice is used in $w_x$-space and $w_\xi$-space, the effective supports of the Gaussian windows in the spatial and spatial frequency domains associated to the hologram have identical extents as in Equations (10) and (12). For the direct Gabor transform, a dual window is employed. With increasing redundancy of the Gabor frame, the effective support of the dual window approaches that of the Gaussian. Therefore, Equation (12) is also used in $w_\xi$-space.

As seen in Section IV-A, at a given spatial position in $w_x$-space, if the samples from different depth layers are not well-separated along the $x$ and $y$ axes of the all-in-focus wavefield, information specific to distinct layer is merged in phase space. To enable depth discrimination when this condition is not satisfied, information from different depth layers is embedded in distinct spatial frequency intervals as illustrated in Figure 3. For this purpose, the phase of the all-in-focus wavefield is defined as the product of a band-limited random phase and a depth-dependent linear phase.

The former restricts the spread of information, while the latter localizes the spatial frequency ranges. Although most of the energy associated to a depth layer is contained in the corresponding interval, the band-limited random phase does not strictly bound the spread of the information. As a result, superpositions between different information clusters may create notable artifacts in the reconstructions especially as the spatial density of samples belonging to different depth layers increases. To limit this issue, information from depth layers that are close along the optical axis are mapped to distant regions along the spatial frequency axes. The cut-off of high spatial frequencies also limits the diffusion of light and therefore the perception of depth. Hence, the spread of the information clusters is subject to a trade-off between visual quality and 3D perception. Besides, the shifts induced by the depth-dependent linear phase have to be compensated during the warping of phase space to prevent any errors in the reconstructions.

As described in Section IV-A, collecting a unique depth value for each coefficient in $w_\xi$-space omits part of the wavefield in the hologram plane. To address this issue, several virtual cameras are placed in the effective spatial support associated to a given $(\tau_\xi, \tau_\eta)$, whose size is defined in Equation (10). Similarly, each of these cameras cast multiple rays over the spatial frequency extent determined by Equation (12). The depths at which rays intersect layers are gathered into a set from which duplicates are removed. The PSDM $d_{w_\xi}$ then consists of the remaining elements and has multiple layers.
<figure><img src="image_4.png" /><figcaption>**Fig. 3. Spatial frequency embeddings.** When Gaussian functions from the Gabor frame spread over samples belonging to different depth layers in the all-in-focus wavefield, information associated to distinct depth is merged in phase space. The depth-dependent linear phase and the band-limited random phase create information clusters in distinct regions along the spatial frequency axes to separate the information specific to each depth layers.</figcaption></figure>

Afterwards, for each of these depth values at a given $w_\xi$, an inverse mapping is applied followed by the use of the transfer function. The corresponding coefficient is then the sum of the resulting values. Consequently, the computation time grows slowly as the number of depth layers in the 3D scene increases. Nevertheless, since the multiple inverse mappings and phase modulations do not necessitate the calculation of several Gabor transforms, this relationship is computationally far less significant than in primitive-based methods.

# V. EXPERIMENTS

## A. Set-up

We compared our phase space warping (PSW) method with the conventional layer-based approach [15] through digital experiments. In the layer-based approach, the depth layers are summed in the spatial frequency domain after being multiplied by the transfer function. All methods were implemented as C++ Unity plugins on a GPU with CUDA. The PC system included an Intel Core i9-9900X operating at 3.5GHz, 32GB RAM and a NVIDIA GeForce RTX 2080 Ti graphics card. The calculation parameters of the PSW method are shown in Table I. The first digital set of parameters was used to assess the effect of interpolation with the Scene 1, see Figure 4. The other set was employed to compare the PSW method with the conventional layer-based approach. In this experiment, there were two different types of scenes. In Scene 2, the samples of different depth layers are well-separated along the *x* and *y* axes, whereas in the Scene 3 this assumption is not satisfied, as illustrated in Figure 5. All experiments were conducted with different redundancy *R* of the Gabor frame to determine the effect of this parameter.

## B. Interpolation Effects

As stated in [9], the coefficients in $w_x$-space oscillate rapidly and interpolating over this lattice may not capture significant information. The following experiment provides a first answer to this concern with Scene 1 consisting of four objects located at different depths. The samples in different

<table><caption>**TABLE I**
**EXPERIMENT PARAMETERS.** THE FIRST SET OF PARAMETERS WAS USED TO ASSESS THE EFFECTS OF INTERPOLATION AND THE LAST TO COMPARE THE PSW METHOD WITH THE LAYER-BASED APPROACH.</caption><thead><tr><td>Parameters</td><td>Variables</td><td>Interpolation</td><td>Analysis</td></tr></thead><tbody><tr><td>Pixel pitch</td><td>δ<sub>ξ</sub> = δ<sub>η</sub></td><td>4.5 μm</td><td>4.5 μm</td></tr><tr><td>Wavelength</td><td>λ</td><td>635 nm</td><td>635 nm</td></tr><tr><td>Pixel count</td><td>N<sub>ξ</sub> × N<sub>η</sub></td><td>2048 × 2048 px</td><td>2048 × 2048 px</td></tr><tr><td>Minimal depth</td><td>d<sub>min</sub></td><td>0.1 cm</td><td>0.1 cm</td></tr><tr><td>Maximal depth</td><td>d<sub>max</sub></td><td>1.375 cm</td><td>2.65 cm</td></tr></tbody></table>

depth layers are well-separated along the *x* and *y* axes to focus on the study of the interpolation effects in phase space. Hence, to diffuse light during the CGH calculation, each depth layer was multiplied by a random diffuser instead of a band-limited random phase. Similarly, no depth-dependent linear phase was employed. A linear interpolation method was compared with an approach that eliminates the need for interpolation. Without interpolation, the method consists in first computing the PSCs of the all-in-focus wavefield on an irregular lattice. The resulting coefficients are warped onto a regular grid and then subject to an inverse Gabor transform. The reconstructions were computed by propagating complex amplitude holograms to the 3D scene with the angular spectrum method and by normalizing the resulting signal. Hence, the PSW method was assessed without considering quantization. The reconstruction quality was measured as a PSNR where the reference is a reconstruction obtained from a CGH computed with the layer-based method.

The reconstruction obtained from a CGH computed according to the PSW method with the linear interpolation and a redundancy of 1 suffers from huge artifacts, as illustrated in Figure 4. Nevertheless, visual quality is significantly enhanced with higher redundancies, especially with a redundancy of 16. On the other hand, computing CGH without any interpolation produces reconstructions of high visual quality even with low redundancy. However, the calculation time for a redundancy of 1 is 1 h 9 min and 32 s. The calculations even became intractable for a redundancy of 16, and were stopped after 60 h. Consequently, the interpolation approach with a high redundancy achieve high visual quality at a lower computational cost than the method that does not use interpolation, as indicated by the results shown in Section V-C and Figure 4.

## C. Calculation Analysis

The computation time of the layer-based and the PSW methods were compared by calculating color CGHs from Scenes 2 and 3. Both scenes were sliced uniformly along the optical axis to define the position of the depth layers. Scene 2 is constituted of multiple objects placed at different depths: the closer is in the top left corner and the farthest at the bottom right corner. With this layout, the objects are well-separated along the *x* and *y* axes of the all-in-focus wavefield. Hence, to diffuse light during the CGH calculation, each depth layer was multiplied by a random diffuser instead of a band-limited random phase. Similarly, no depth-dependent linear phase was employed. Scene 3 represents a lake with multiple trees where the samples of the depth layers are not separated
<figure><img src="image_5.png" /><figcaption>Fig. 4. **Digital reconstructions.** Digital reconstructions obtained from holograms computed with the PSW method are compared as a function of the interpolation approach. Each row is associated with an interpolation strategy: no interpolation and linear, respectively. Each column correspond to a given redundancy from 1 to 16.</figcaption></figure>

Therefore, information from different depth layers is embedded in distinct spatial frequency interval with a band-limited random phase and a depth-dependent linear phase. The calculation time includes neither the loading and saving operations nor the computation of the PSDM that can leverage conventional ray casting techniques.

<table><caption>TABLE II
CALCULATION TIME. THE CALCULATION TIMES FOR COLOR HOLOGRAMS USING THE LAYER-BASED APPROACH AND THE PSW METHOD ARE COMPARED AS FUNCTIONS OF THE NUMBER OF DEPTH LAYERS USED TO REPRESENT THE 3D SCENE.</caption><thead><tr><th rowspan="3">Layers</th><th rowspan="3">Layer-based (ms)</th><th colspan="4">PSW (ms)</th></tr><tr><th colspan="2">Scene 2</th><th colspan="2">Scene 3</th></tr><tr><th>R = 4</th><th>R = 16</th><th>R = 4</th><th>R = 16</th></tr></thead><tbody><tr><td>8</td><td>71.56</td><td>19.49</td><td>61.09</td><td>20.62</td><td>66.42</td></tr><tr><td>16</td><td>135.64</td><td>19.54</td><td>61.91</td><td>21.08</td><td>69.24</td></tr><tr><td>32</td><td>262.52</td><td>19.58</td><td>62.22</td><td>22.19</td><td>73.17</td></tr><tr><td>64</td><td>513.13</td><td>19.63</td><td>62.33</td><td>24.18</td><td>84.44</td></tr></tbody></table>

The results shown in Table II prove that even with a few depth layers the PSW method is faster than the layer-based approach and up to 26.1 times faster. As predicted by the theoretical analysis, the calculation time as a function of the number of layers rises much faster with the layer-based approach than with the PSW method. Indeed, the number of Fourier transforms and depth layers increases jointly in the layer-based approach, whereas in the PSW method only two

Gabor transforms are computed, independently of the number of layers.

Considering Scene 2, at given coordinates $(\tau_x, \tau_y)$ in the all-in-focus wavefield, all samples belong to the same depth layer. Consequently, at a given position of the PSDM, there are only few different values even when the number of depth layers increases. The number of inverse warpings needed during the CGH calculation is then low and rises slowly as a function of the number of layers. This trend is emphasized in Table II, where the calculation time with Scene 2 is lower and increases at a slower rate than that with Scene 3. In addition, in both experiments the computation time is higher with a redundancy of 16 than with $R = 4$ since in this configuration the Gabor transform and the warpings of phase space are computationally more expensive.

## D. Quality Assessment

The visual quality provided by the PSW method was assessed with the Scenes 2 and 3, constituted of 64 depth layers, using the same computation process as in Section V-C. The PSW method was assessed without considering quantization. Hence, the reconstructions were computed by propagating the complex amplitude holograms to the 3D scene with the angular spectrum method and by normalizing the resulting signal. The reconstruction quality was measured as a PSNR where
<figure><img src="image_6.png" /><figcaption>Fig. 5. **Digital reconstructions.** Digital reconstructions obtained from holograms computed with the layer-based and the PSW methods are compared. The first and second row are associated with reconstructions from Scene 2 and 3. Columns correspond to the following methods: layer-based, PSW with *R* = 4 and PSW with *R* = 16.</figcaption></figure>

the reference was a reconstruction obtained from a complex amplitude CGH computed with the layer-based approach.

In Scene 2, information from distinct depth layers is merged in both the spatial and spatial frequency domains after propagation, due the use of a random diffuser. Nevertheless, the PSW method still achieves high visual quality as illustrated in the first row of Figure 5. This result is in agreement with the theoretical expectations because phase space is a domain in which information specific to each depth layers is often concentrated in distinct regions. This statement is especially confirmed in Scene 2 because all the samples at a given spatial position ($\tau_x, \tau_y$) belong to the same depth layer. This layout is leveraged by the PSW method to propagate light waves from a 3D scene without considering multiple primitives individually. It is worth noting that the visual quality increases with the redundancy.

Scene 3 is less suited for the PSW method because samples belonging to distinct depth layers are located at the same spatial position ($\tau_x, \tau_y$). Thus, unlike in Scene 2, information specific to different depth layers is not inherently concentrated in distinct regions of phase space. Nevertheless, by combining the band-limited random phase with the depth-dependent linear phase, this information is separated and acceptable quality is achieved as illustrated in the last row of Figure 5. The PSW method with a redundancy of 4 provides a PSNR of 29 dB and visual quality improves with a redundancy of 16. Note that the

<figure><img src="image_7.png" /><figcaption>Fig. 6. **Digital reconstructions.** Digital reconstructions obtained from holograms computed on Scene 3 with the PSW methods without using the workarounds presented in Section IV-B are shown. The first and second images are associated with redundancy of 4 and 16, respectively.</figcaption></figure>

workarounds presented in Section IV-B are essential to obtain these results. Without these techniques, reconstructions suffer from huge artifacts, as illustrated in Figure 6.

## VI. DISCUSSION

In this article, we propose to model diffraction from a 3D structure in phase space to provide both efficient and rigorous calculation, which is essential in many fields of optics such as computer-generated holography for 3D visualization. For this purpose, we have extended previous formulations [9], [12] that
were restricted to light propagation between two planes. Our method is up to 26.1 times faster than the conventional layer-based approach [15] while maintaining a 35.1 dB similarity in PSNR. Additionally, by multiplying the number of depth layers by 8, the computational load of the layer-based method grows linearly, whereas with our method, the computation time only increases by 0.12 ms, reducing dependency on scene geometry. This feature stems from the wavefield decomposition as a set of Gaussian beams in phase space that enable a finer control of light propagation than in the spatial or the spatial frequency domain alone.

We have also investigated the quality loss induced by interpolation over PSCs which is essential to determine whether fast and rigorous diffraction modeling is possible. Our results indicate that a linear interpolation over a Gabor spectrum computed with a redundancy of at least 4 is sufficient for reproducing with high fidelity (37.9 dB) rigorous light propagation from a 3D scene. Additional study could either investigate more advanced interpolation strategies or fast implementations of Gabor transform on irregular lattices to reduce the loss of information at low redundancy while maintaining computation time low.

Despite the improvements of our method, a major limitation is the use of a band-limited random phase when the spatial density of samples belonging to different depth layers increases, restricting the diffusion of light and thus 3D perception. To tackle this issue, the geometry of the local spatial area from which Gaussian beams are emitted should be accounted for during the warping of phase space. For instance, we can consider a beam propagating from a plane that is tilted relative to the hologram plane. In this situation, with our method, the spatial support of the Gaussian function gathers samples belonging to different depth layers requiring the use of a band-limited random phase and a depth-dependent linear phase. Instead, we could warp phase space to consider that a single Gaussian beam is emitted from a tilted plane. Future work could also address view-dependent effects, such as occlusions and lighting, which can easily be handled in phase space.

## VII. CONCLUSION

We proposed to model in phase space light diffraction from a 3D structure in a medium of constant refractive index by extending previous formulations restricted to 2D planes. Our digital experiment demonstrated up to 26.1-fold speedup while maintaining a similarity of 35.1 dB in PSNR compared to the conventional layer-based method, due to a near-constant scaling relative to the scene geometry. These results stem from the wavefield decomposition into a set of Gaussian beams providing a fine control of light propagation. This study has a particular impact for application in computer-generated holography for 3D visualization to enable both efficient calculation and high visual realism.

## ACKNOWLEDGMENTS

The 3D model associated to Scene 3 was downloaded from the Unity Asset Store.

## REFERENCES

[1] J. W. Goodman, *Foundations of Scalar Diffraction Theory*, 3rd ed. Roberts and Company Publishers, 2005, pp. 32–62.

[2] X. Yu, H. Zhang, Z. Zhao, X. Fan, S. Hu, Z. Li, W. Chen, D. Li, S. Shi, W. Xiong, and H. Gao, “On the use of deep learning for computer-generated holography,” *iScience*, vol. 28, no. 5, p. 112507, May 2025.

[3] C. Zhong, X. Sang, B. Yan, H. Li, D. Chen, X. Qin, S. Chen, and X. Ye, “Real-Time High-Quality Computer-Generated Hologram Using Complex-Valued Convolutional Neural Network,” *IEEE Transactions on Visualization and Computer Graphics*, vol. 30, no. 7, pp. 3709–3718, Jul. 2024.

[4] P. W. M. Tsang, T.-C. Poon, and Y. M. Wu, “Review of fast methods for point-based computer-generated holography [Invited],” *Photonics Research*, vol. 6, no. 9, pp. 837–846, Sep. 2018.

[5] Y. Zhang, H. Fan, F. Wang, F. Wang, X. Gu, X. Qian, and T.-C. Poon, “Polygon-based computer-generated holography: a review of fundamentals and recent progress [Invited],” *Applied Optics*, vol. 61, no. 5, pp. B363–B374, Feb. 2022.

[6] C. Chang, X. Ding, D. Wang, Z. Ren, B. Dai, Q. Wang, S. Zhuang, and D. Zhang, “Split Lohmann computer holography: fast generation of 3D hologram in single-step diffraction calculation,” *Advanced Photonics Nexus*, vol. 3, no. 3, p. 036001, Mar. 2024.

[7] J.-P. Liu and S.-L. Lu, “Fast calculation of high-definition depth-added computer-generated holographic stereogram by spectrum-domain lookup table [Invited],” *Applied Optics*, vol. 60, no. 4, pp. A104–A110, Feb. 2021.

[8] K. Min, D. Min, and J.-H. Park, “Wigner inverse transform based computer generated hologram for large object at far field from its perspective light field,” *Optics Communications*, vol. 532, p. 129229, Apr. 2023.

[9] P. Gioia, A. Gilles, A. El Rhammad, and S. Vū Ngoc, “Phase Space Formulation of Light Propagation on Tilted Planes,” *Photonics*, vol. 11, no. 11, p. 1034, Nov. 2024.

[10] P. Gioia and S. V. Ngoc, “Accurate semiclassical analysis of light propagation on tilted hyperplanes,” 2025.

[11] P. Gioia, A. Gilles, A. Lagrange, A. El Rhammad, and S. Vu-Ngoc, “Linear canonical transformations in phase space: the Gabor frames approach,” in *Optics, Photonics, and Digital Technologies for Imaging Applications VIII*, vol. 12998, Jun. 2024, p. 129980I.

[12] D. Blinder, T. Birnbaum, and P. Schelkens, “Efficient numerical Fresnel diffraction with Gabor frames,” *Photonics Research*, vol. 13, no. 2, pp. 330–339, Feb. 2025.

[13] M. E. Lucente, “Interactive computation of holograms using a look-up table,” *Journal of Electronic Imaging*, vol. 2, no. 1, pp. 28–34, Jan. 1993.

[14] K. Matsushima, H. Schimmel, and F. Wyrowski, “Fast calculation method for optical diffraction on tilted planes by use of the angular spectrum of plane waves,” *Journal of the Optical Society of America A*, vol. 20, no. 9, pp. 1755–1762, Sep. 2003.

[15] Y. Zhao, L. Cao, H. Zhang, D. Kong, and G. Jin, “Accurate calculation of computer-generated holograms using angular-spectrum layer-oriented method,” *Optics Express*, vol. 23, no. 20, pp. 25440–25449, Oct. 2015.

[16] F. Wang, T. Ito, and T. Shimobaba, “High-speed rendering pipeline for polygon-based holograms,” *Photonics Research*, vol. 11, no. 2, pp. 313–328, Feb. 2023.

[17] B. Moussa Fares, H.-s. Kim, P. Schelkens, and D. Blinder, “Sparse point cloud computer-generated holography with the Gabor transform,” *Optics Express*, vol. 33, no. 16, pp. 34303–34319, Aug. 2025.

[18] H. Nishi and K. Matsushima, “Rendering of transparent objects in large-scale full-parallax polygon-based computer holography,” *Optics Express*, vol. 33, no. 1, pp. 704–716, Jan. 2025.

[19] K. Watanabe, K. Yamauchi, and Y. Sakamoto, “Realistic rendering method for specular reflections with continuous motion parallax in a computer-generated hologram,” *Applied Optics*, vol. 63, no. 7, pp. B126–B133, Mar. 2024.

[20] X. Duan and Z. Le, “Inter-layer edge supplementation in computer-generated multi-depth 3D holograms,” *Optics Express*, vol. 33, no. 10, pp. 21453–21463, May 2025.

[21] D. Blinder, T. Nishitsuji, and P. Schelkens, “Three-dimensional spline-based computer-generated holography,” *Optics Express*, vol. 31, no. 2, pp. 3072–3082, Jan. 2023.

[22] S. Choi, B. Chao, J. Yang, M. Gopakumar, and G. Wetzstein, “Gaussian Wave Splatting for Computer-Generated Holography,” *ACM Transactions on Graphics*, vol. 44, no. 4, p. 57, Aug. 2025.

[23] T. Yatagai, “Three-dimensional displays using computer-generated holograms,” *Optics Communications*, vol. 12, no. 1, pp. 43–45, Sep. 1974.
[24] D. Blinder and P. Schelkens, “Phase added sub-stereograms for accelerating computer generated holography,” *Optics Express*, vol. 28, no. 11, pp. 16924–16934, May 2020.

[25] H. Zhang, L. Cao, and G. Jin, “Three-dimensional computer-generated hologram with Fourier domain segmentation,” *Optics Express*, vol. 27, no. 8, pp. 11689–11697, Apr. 2019.

[26] D. Min, K. Min, H.-J. Choi, H. Lee, and J.-H. Park, “Non-hogel-based computer generated hologram with occlusion processing between the foreground light field and background hologram,” *Optics Express*, vol. 30, no. 21, pp. 38339–38356, Oct. 2022.

[27] N. Padmanaban, Y. Peng, and G. Wetzstein, “Holographic near-eye displays based on overlap-add stereograms,” *ACM Transactions on Graphics*, vol. 38, no. 6, p. 214, Nov. 2019.

[28] Z. Wang, L. M. Zhu, X. Zhang, P. Dai, G. Q. Lv, Q. B. Feng, A. T. Wang, and H. Ming, “Computer-generated photorealistic hologram using ray-wavefront conversion based on the additive compressive light field approach,” *Optics Letters*, vol. 45, no. 3, pp. 615–618, Feb. 2020.

[29] C. Chang, J. Wu, Y. Qi, C. Yuan, S. Nie, and J. Xia, “Simple calculation of a computer-generated hologram for lensless holographic 3D projection using a nonuniform sampled wavefront recording plane,” *Applied Optics*, vol. 55, no. 28, pp. 7988–7996, Oct. 2016.

[30] A. Lagrange, A. Gilles, K. Heggarty, B. Fracasso, P. Gioia, and V. Savaux, “Real-time calculation of Fresnel holograms using the split-Lohmann method,” *Optics Letters*, vol. 50, no. 8, pp. 2554–2557, Apr. 2025.

[31] J. Zhou, S. Xie, L. Jiang, Y. Luo, Y. Wu, and J. Wang, “Lateral-shift-free Split-Lohmann computer-generated holography,” *Optics Letters*, vol. 50, no. 21, pp. 6558–6561, Nov. 2025.

[32] T. Birnbaum, T. Kozacki, and P. Schelkens, “Providing a Visual Understanding of Holography Through Phase Space Representations,” *Applied Sciences*, vol. 10, no. 14, p. 4766, Jan. 2020.

[33] K. Gröchenig, *Discrete Time-Frequency Representations: Gabor Frames*, 1st ed., ser. Applied and Numerical Harmonic Analysis. Birkhäuser Boston, Dec. 2001, pp. 82–102.

[34] P. L. Søndergaard, “Efficient Algorithms for the Discrete Gabor Transform with a Long Fir Window,” *Journal of Fourier Analysis and Applications*, vol. 18, no. 3, pp. 456–470, Jun. 2012.

[35] V. Guillemin and S. Sternberg, *Introduction*, 1st ed. Cambridge University Press, 1984, pp. 1–150.