Computer-generated holography (CGH) represent a class of algorithms for digitally generating holographic interference patterns for many different applications, such as beam shaping [1], optical trapping [2], neural photo-stimulation [3], particle imaging [4], and electro-holographic display systems [5]. Because CGH models diffraction, local features positioned in space will emit outwardly propagating waves that can affect objects everywhere, including all hologram pixels. CGH is therefore calculation intensive, which is why developing computationally efficient and physically accurate CGH algorithms is challenging [6], making it an ongoing research problem.

CGH algorithms come in many kinds, encoding 3D objects and scenes using point clouds [7], polygons [8], layers [9], surface elements [10], or ray tracing [11]. Moreover, many acceleration techniques have been proposed, such as leveraging sparsity [12,13], optimizing caching [14], look-up tables [15], and, more recently, deep learning [16]. We will focus on the characteristics of the point-spread function (PSF) (or Fresnel zone plates) used primarily in point-cloud algorithms.

When computing the PSF at short distances, the finite sampling rate needs to be considered. As the PSF’s instantaneous spatial frequency increases outwardly with its radius, and spatial frequencies are proportional to the incidence angle, evaluating the hologram samples illuminated at too steep angles will cause aliasing, leading to artifacts. A valid evaluation of PSFs can be obtained by determining the maximum PSF radius of influence that will not cause aliasing, leading to a disk-shaped PSF. Aside from eliminating aliasing, this procedure has the important added benefit that it significantly reduces the number of pixels that need to be evaluated, leading to important gains in calculation speed. This principle was initially proposed in Ref. [12], computing the wave field in an intermediate plane close to the point cloud, denoted a wavefront recording plane (WRP), followed by a fast convolutional diffraction operation.

This disk-shaped PSF (or circular zone plate) has been consistently used in many works since [7,17–19], sometimes even using the inscribed square to that disk. However, these widely used shapes are sub-optimal: at least $21.5\%$ and $50\%$ of the available bandwidth is discarded, respectively. In this Letter, we propose an optimal PSF shape, which is a hyperbolic square, or pincushion. We obtain the shape and its bounds analytically, derive a geometric model, show that the new shape reduces the error by a factor of $\approx 3$, and show how this improves the visual quality of holograms rendered with a ray-tracing-based CGH algorithm [11].

The PSF $P$ of a point with coordinates $(x_0,y_0,z_0)$ on the hologram plane $z=0$, centered at the origin, is

Numerically evaluating $P$ requires discretization, which will typically be done using a regular pattern. We assume that the samples are spaced by a sampling pitch $p$, which is also called the pixel pitch, generally matching the actual pixel pitch of a holographic display; we will use both terms interchangeably. The classical disk-shaped PSF radius $w$ can be found by determining the maximum diffraction angle $\theta$, given by the relationship $\sin (\theta ) = {\lambda } / {2p}$. Using the geometrical relationships as described in Refs. [10,12], we get

For a hologram sampled regularly on an axis-aligned square lattice with sampling pitch $p$, signals with frequencies surpassing the Nyquist bound ${1} / {2p}$ along the main axes will cause aliasing. Because the PSF $P$ has a well-defined unique instantaneous frequency $\nu$ in every point, we can establish the non-aliased region by determining for what points this frequency does not surpass the Nyquist bound. The instantaneous frequency components $\nu _x$ and $\nu _y$ are given by

If we take the square of the last inequalities of Eqs. (5) and (6), we get an expression for two bounding hyperbolas, as shown in Fig. 1.

 

Fig. 1. Example of pincushion PSF: the optimal PSF is bounded by two hyperbola, horizontally (in red) and vertically (green). This area is always strictly larger than the conventionally used PSF shown within the disk demarcated by the circle with radius $w$ (shown in blue).

Download Full Size | PDF

The origin of the pincushion shape can also be understood geometrically. We can represent the propagating point emitter wavefront by a sphere, where all points on its surface will correspond to an instantaneous frequency proportional to the incidence angle of the corresponding wave vector at that point with the hologram plane. We can demarcate the valid frequencies by making an orthogonal projection of the bandwidth, represented by a square in Fourier space, onto the sphere, as shown in Fig. 2(a). The relative size of the square will be inversely proportional to $p$. Then, we only consider points on the hologram plane whose rays traversed the projected square from the center of the sphere, also known as a gnomonic projection, since only those rays (or wave vectors) will correspond to frequencies below the Nyquist bound. The resulting shape will be a pincushion, as shown in Fig. 2(b).

 

Fig. 2. Geometric illustrations of PSF shape: (a) orthographic projection of a square on a sphere, representing the hologram bandwidth; (b) gnomonic projection of an outwardly propagating PSF to a pincushion shape on top (in red), representing all pixels whose incidence angle does not surpass the maximum diffraction angle $\theta$.

Download Full Size | PDF

The disk-shaped PSF will only induce frequencies up to magnitude ${1} / {2p}$, which will map to a disk with radius ${1} / {2p}$ in Fourier space; but since the hologram is sampled using a square lattice, its actual bandwidth is a square with side length ${1} / {2p}$. Therefore, only the frequencies within the inscribed circle of that square are used, constituting only of $\pi /4 \approx 78.5\%$ of the total available bandwidth.

When efficiently implementing the PSF calculation in software, it is useful to know the bounding box of the proposed PSF, i.e., the circumscribed axis-aligned square to the pincushion shape. We can find these bounds by finding out where the two hyperbola intersect. Since one intersection will happen on the half-diagonal $x=y>0$, we can find the $x$-coordinate by substituting it in the last equality of Eq. (5):

This also shows that the precise pincushion shape depends on the ratio $p/\lambda$. To illustrate, the PSF shape is plotted for various ratios in Fig. 3. As the ratio gets smaller, the hyperbola recede, making the pincushion shape more pronounced, back to a point when they do not intersect anymore; the PSF is then unbounded, because the maximally allowed frequency will never be reached along some directions. Conversely, as the ratio gets larger, the PSF increasingly approaches the shape of the circumscribed square. This is consistent with the quadratic Fresnel approximation valid for small diffraction angles $\theta \approx \sin (\theta )$, using the approximated PSF

 

Fig. 3. (a)–(e) Example proposed pincushion PSFs sampled at different pixel pitches $p$, for $\lambda =660$ nm, depth $d = {p^2} / {\lambda ^2}\cdot 0.313$ mm, resolution = $1024\times 1024$. Only the real part is shown; unevaluated values beyond the aliasing limit lie outside of the pincushion shapes (shown in light purple). Note how the PSF shape approaches a $w\times w$ square for larger $p$, becoming a circumscribed square to the PSF disk with radius $w$.

Download Full Size | PDF

In the first series of experiments, we evaluated the increased accuracy of the pincushion PSF. To do so, we calculated $4096\times 4096$ pixel test holograms with $\lambda ={660}\mathrm{nm}$ and variable $p$ of a single PSF, from a laterally centered point at a distance ${p^2} / {\lambda ^2} = {0.2}\mathrm{mm}$ from the hologram plane. This variable distance $d$ was chosen to obtain roughly constant PSF sizes along the $x$/$y$ axes, as shown in Fig. 3. The corresponding respective reference holograms were computed each time at the maximum sampling rate of $\lambda /2$, i.e., the resolution of the reference hologram increased by a factor of ${2p} / {\lambda }$ to match the resolution of the respective test holograms. Then the reference holograms were down-sampled to the matching $4096\times 4096$ pixel resolution using a sinc filter (i.e., by cropping the $4096\times 4096$ center frequencies in the Fourier domain and normalizing).

The error between the reference and test holograms was computed with the normalized mean square error (NMSE):

 

Fig. 4. Comparison between ground-truth and proposed pincushion PSF: (a) the reference signal was oversampled at $p=\lambda /2$, followed by down-sampling to match the pixel pitch and resolution of (b) the target PSF, with $p=\lambda$ only non-zero at points satisfying Eqs. (5) and (6).

Download Full Size | PDF

Table 1. Error Comparison between Newly Proposed and Disk-shaped PSFs, for Different Pixel Pitch Sizes^a

View Table

For the visual quality experiments, we use the “photorealistic CGH” algorithm using ray tracing, detailed in Ref. [11]. This algorithm will sample the 3D scene into a point cloud, whose zone plates are spatially modulated using ray tracing, accounting for occlusion and complex lighting effects. The original paper uses a disk-shaped zone plate (i.e., PSF), which we compare with an amended version using our proposed pincushion-shaped PSF instead. We generate a hologram for each PSF shape to compare the effect on the perceived views.

The virtual scene is based on the wooden Cornell box from [11] with both a square area ceiling light and a front-facing light source, several detailed textured models, showing specular highlights, soft shadows, and global illumination. We generated holograms using $2$ million point samples with a wavelength of $\lambda ={660}\mathrm{nm}$, a pixel pitch of $p={2}\mathrm{\mu m}$ and a resolution of $16384\times 16384$ pixels. The algorithm ran on a machine with an Intel Xeon E5-2687W v4 processor @ 3 GHz, 256 GB of RAM, with a NVIDIA TITAN RTX running on Windows 10 OS. The program was implemented in C++17, with CUDA 11.3 and OptiX 7.3. The calculation took $\sim 6$ h per hologram. Both a central view and a top-left corner view were extracted from each hologram using a $2048\times 2048$ pixel aperture in the Fourier domain. The views are shown side-by-side in Fig. 5.

 

Fig. 5. Orthographic views extracted from two holograms generated with different PSFs, all generated with numerical simulations. (a), (b) As expected, no difference is visible in the center views. (c), (d) The corner views are markedly different: using (c) a disk PSF will degrade the visual quality, resulting in a significant loss of detail and ghosting, while (d) the proposed PSF gives a sharp view; (e) and (f) are enlarged views taken from (c) and (d), respectively.

Download Full Size | PDF

The chosen shape of the zone plate or PSF has a significant effect on signal accuracy and visual quality, reducing the error by about three times and better encodes high-frequency signals. This novel PSF shape can easily be integrated in many CGH algorithms using PSFs, such as WRPs, analytical methods, look-up tables, and ray tracing CGH. Moreover, this finding can have uses beyond CGH for holographic displays utilizing numerical diffraction models or circular zone plates, including applications at non-optical wavelengths.