Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Direct calculation of computer-generated holograms in sparse bases

Open Access Open Access

Abstract

Computer-generated holography is computationally intensive, making it especially challenging for holographic displays where high-resolutions and video rates are needed. To this end, we propose a technique for directly calculating short-time Fourier transform coefficients without the need for a look-up table. Because point spread functions are sparse in this transform domain, only a small fraction of the coefficients need to be updated, enabling significant speed gains. Twenty-fold accelerations are reported over the reference implementation. This approach generalizes the notion of the phase-added stereogram, allowing for the calculatiion of an arbitrary number of Fourier coefficients per block, enabling high calculation speed with holograms of good visual quality, targeting minimal memory requirements.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Numerical diffraction simulates how coherent light propagates through free space. Because of the properties of waves, these calculations are computationally costly: every point in the 3D scene can potentially affect every other point in space.

This problem is especially important for holographic displays. They are considered to be the ultimate 3D display technology [1], because holograms can account for all human visual cues, including stereopsis, eye-focusing, parallax and no exhibit accommodation-vergence conflict. These displays require holograms to be calculated at video rates and at high resolutions, which require highly efficient Computer Generated Holography (CGH) algorithms.

CGH algorithms come in many forms and have various trade-offs [1,2], such as: calculation time, what graphical effects are supported (occlusion, shadows, reflectivity, etc.) and what type of 3D scenes and models can be accounted for. They can be classified as point cloud methods [3–8], layer-based methods [9–12], polygon methods [13–16], RGB+depth methods [17, 18], and ray-based methods [19–23]. These categories are not sharply defined or mutually exclusive, as hybrid algorithms exist [24].

This paper will focus on a subset of the point cloud algorithms, which we call “sparse CGH algorithms”. Sparsity is a concept from signal processing literature; it is a measure of how few coefficients are needed to accurately describe a signal in particular basis. For examples, images such as natural photographs typically are sparse in wavelet space, making the wavelet transform highly suitable for compressing images [25], since most coefficients will be (almost) zero.

Sparse CGH algorithms thus calculate the hologram in a well-chosen transform domain, so that most of the signal’s energy is concentrated in a small number of coefficients, thereby requiring much fewer updates than the total hologram pixel count; this can speed up hologram calculation considerably. The challenge is to choose the right transform, balancing sparsity (i.e. minimizing the number of needed coefficients) with the complexity of updating individual coefficients and the calculation cost of the inverse transform, maximizing performance and quality.

Presently, there are only a few examples of sparse CGH in literature [1]. Although initially not described as such, it can be argued that the first instance of sparse CGH is the well-known Wavefront Recording Plane (WRP) [5], which calculates Point Spread Functions (PSFs) in an intermediate plane, which effectively corresponds to performing CGH in a convolved space. Because the propagation distance of a point to the intermediary plane is relatively short, the spatial PSF support is small, rendering them sparser in the spatial domain. After accumulating all the PSFs, the intermediary plane can be propagated efficiently using a convolution, such as the Fresnel transform or the Angular Spectrum Method (ASM) [26]. This notion can be extended further by using multiple WRPs, allocating points to their respective nearest WRP [9,11].

Recently, another type of approach has been proposed that leverages the redundancies within a PSF in the wavefront plane instead of only relying on optimizing WRP depth placement. This was first proposed by Shimobaba et al. [27], by precomputing PSF coefficients in wavelet space using a Look-Up Table (LUT). Because neighboring PSF samples are highly correlated, much like in conventional imagery, fewer coefficients are needed to accurately represent PSFs. Later, sparse CGH was proposed using the Short-Time Fourier Transform (STFT), which was shown to yield even higher sparsity and better preservation of views at large angles [28]. Both these approaches can be combined with WRPs to even further increase sparsity [28,29].

With this outlook, the phase-added stereogram [20] using the Fast Fourier Transform (FFT) [30] can be viewed as a special case of the sparse STFT CGH: for every PSF, a single STFT coefficient is updated in every block. Albeit fast, this approximation is relatively coarse, because the number of available discrete frequencies is limited. This will reduce the visual quality and hence limiting the application domain. Several extensions have been proposed [31] to improve quality, but these may significantly increase memory requirements by oversampling the FFT space by using blocks which are much larger than those of the basic method.

However, typically sparse CGH algorithm rely on a LUT for their operation, which introduces limitations. To this end, we propose the first sparse CGH algorithm without a LUT. This makes the system much more flexible: (1) PSFs positions are not discretized anymore for a fixed set of 3D locations; (2) no need for precomputing and storing sizeable LUT or using overcomplete representations, enabling a wider range of applications such as FPGAs (Field Programmable Gate Arrays) and embedded systems; (3) it reduces memory-bound performance limitations, making the algorithm more cache-friendly. Hence, the proposed algorithm can be viewed as a generalization of the phase-added stereogram.

The main novel contributions are the following:

  • A derivation of an analytical model for directly calculating STFT coefficients, and proposing an efficiently computable and accurate approximation.
  • An algorithm to selectively update a small subset of the highest-magnitude STFT coefficients, with a configurable sparsity level trading-off speed and quality.
  • A validation reporting a 20-fold speedup over the conventional algorithm using a C++ implementation tested on 3D point cloud models, and an error analysis measuring how much the proposed algorithm deviates from a reference PSF for various sparsity levels.

The paper is organized as follows: In section 2, we first derive the analytical expression of the STFT coefficients for PSFs; then an accurate and efficiently computable approximation is proposed, followed by a technique to select which subset of coefficients to update. Then, in section 3, the experiments are divided in two parts: (1) measurment of the accuracy of the proposed algorithm for various sparsity levels, and (2) measuring the gain in calculation speed w.r.t. a reference method, including a comparison in visual quality. Finally, we conclude in section 4 and discuss potential avenues for future work.

2. Methodology

The Fresnel approximation of a PSF P, laterally displaced with a translation (δ, ) and at a depth z from the hologram plane (ignoring the scaling factor) is given by

P(x,y)=exp(iγ2[(xδ)2+(y)2])
where i is the imaginary unit, γ2=πλz with wavelength λ.

This expression can be used to calculate a hologram H by summing over a collection of M points Pj, j ∈ {1, ..., M}, each with their own values (γj, δj, j):

H(x,y)=j=1MPj(x,y)=j=1Mexp(iγj2[(xδj)2+(yj)2])

Sparse CGH will compute the values of Pj indirectly using a linear transform 𝒯. If 𝒯 is well-chosen, the PSF Pj will be sparse and can be computed using only a few coefficients. Because of its linearity, we get that

H(x,y)=𝒯1{𝒯{H(x,y)}}=𝒯1{j=1M𝒯{Pj(x,y)}}
meaning we can first accumulate all Pj in the transform domain, and only need to apply the inverse transform 𝒯−1 once on the summed PSF coefficients. If 𝒯 can be computed quickly compared to the summation, we can obtain a large speedup. Here, we choose 𝒯 to be the STFT.

The STFT divides the hologram into evenly-sized blocks, followed by a local Discrete Fourier Transform (DFT) on each block. The goal is to get an expression for directly calculating these coefficients for any given P(x, y) without relying on a LUT. To obtain the local Fourier coefficients of a hologram block of dimensions 2B × 2B, we have to evaluate the following integral f :

f(ω,η)=B+BB+BP(x,y)exp(i(xω+yη))dxdy
Since both the PSF P as well as the Fourier transform are separable, so is f. We thus need to compute f(ω, η) = fx(ω) · fy(η), where
fx(ω)=B+Bexp(i[γ2(xδ)2+xω])dxandfy(η)=B+Bexp(i[γ2(y)2+yη])dy
which are identical up to constants.

The function fx cannot be written using only elementary functions, but it can be expressed in terms of the error function

erf(x)=2π0xet2dt
with t as the variable of integration. We can now rewrite fx(ω) =
1+i4γ2πexp(i(ωδ+ω24γ2))[erf(1i2((B+δ)γ+ω2γ))+erf(1i2((Bδ)γω2γ))].
This expression can be simplified further. Firstly, we are not interested in the leading complex constant, since this constant phase delay and scaling will not affect the appearance of the hologram. Secondly, we can group the constants to obtain:
fx(ω)1γexp(i(ωδ+α2))[E(Bγ+β)+E(Bγβ)]
where α=ω2γ and β = δγ + α. E(x)=cerf(1i2x) is proportional to the error function evaluated on a diagonal of the complex plane, where c ∈ C can be chosen freely. Directly evaluating the integral from Eq. (6) to calculate E(x) would be too computationally costly, so this will be tackled in the next subsection.

2.1. Efficiently approximating E(x)

The aim is to find an efficiently computable, yet accurate approximation of E(x). This is needed if we want the sparse STFT coefficient computation to be faster than the conventional PSF computation in the spatial domain. To achieve this, an approximation was chosen with two use cases: one for small absolute values of x, and one for large values of |x|. Furthermore, c=(i1)π2 was selected to simplify the subsequent expressions. For large values, we model the limit behavior of E using its asymptotic expansion:

E(x)=c+exp(ix2)xiexp(ix2)2x33exp(ix2)4x5+𝒪(1x7).(forx+)
The series is truncated and simplified to a function EL for x → ±∞:
EL(x)=exp(ix2)x+csgn(x)
where sgn(x) is the sign function, equal to 1 when x is positive and equal to −1 when x is negative. Given that E(x) is an odd function, namely erf(x) = −erf(−x), the same expansion can be reused for very small values when x → −∞.

However, this approximation cannot be used for small values of |x|; as x approaches 0, the singularity in the denominator of the first term of EL(x) will dominate, leading to large deviations from E(x). Instead, we can use the Taylor series around x = 0. We get the following Taylor approximation ES(x) valid for small values of |x|:

{E(x)}{ES(x)}=23x3+121x71660x11
{E(x)}{ES(x)}=2x15x5+1108x9
for the real and imaginary parts, respectively. Note that the terms are separated by powers of x4, meaning that x4 needs only to be computed once and can be reused for evaluating the polynomials more efficiently.

What remains to be done is to find out when EL or ES should be selected for any given value of x. This can be determined by looking at their errors w.r.t. to E, and finding the intersection point. On Fig. 1 we plot the Mean Squared Error (MSE) using a logarithmic plot scale. Numerically, the intersection point is found to be x ≈ 1.609. We thus define the approximation of E(x) to be

E˜(x)={ES(x),if|x|<1.609EL(x),otherwise.
The values of E(x) and (x) are shown (for the imaginary part) on Fig. 2.

 figure: Fig. 1

Fig. 1 MSE of the approximations for small values (ES(x)) and large values (EL(x)) compared to the reference E(x), on a logarithmic scale. The crossing point of both curves (at x ≈ 1.609) is chosen to select the approximation with the smallest error depending on x. Given the symmetry of E(x), the curves look exactly the same (but mirrored) for negative values of x.

Download Full Size | PDF

 figure: Fig. 2

Fig. 2 Comparison of the (imaginary parts of) E(x) and (x). The curves show good agreement (see Fig. 1 for the error analysis).

Download Full Size | PDF

2.2. Accurate approximations of the PSF

Now that we have an expression f(ω, η) for directly computing the PSF in the STFT domain, a new question arises: how do we determine what subset of the coefficients to compute? The location of the highest-magnitude coefficients should be known a priori, since calculating all of them and finding the maximum would defeat the purpose of using a sparse representation.

This problem can be addressed by analyzing the simultaneous time-frequency representation of signals (a.k.a. Wigner space or phase space). Most of the energy of a single STFT coefficient is concentrated within a bounded region in space and bandwidth, which can be represented graphically by Heisenberg boxes (cf. Fig. 3). Every column represents a single STFT block, and every row are all coefficients of the same frequency.

 figure: Fig. 3

Fig. 3 Time-frequency charts of PSF curves superimposed on STFT coefficients. The red curves represent the instantaneous frequency of the two PSF instances, whose slope will depend on z and lateral displacement on δ. The (Heisenberg) boxes correspond to individual coefficients, indicating where their energy is concentrated in space (or time) and frequency. The blue boxes are the coefficients that will be computed for this PSF. Both (a) and (b) have the same sparsity rates s, but (b) has a larger block size N trading off more frequency resolution for less spatial resolution, which is better for points with large |z| values.

Download Full Size | PDF

The energy distribution of the STFT coefficients can be found by looking at the instantaneous frequency of PSFs. The instantaneous frequency ν(x) of a (one-dimensional) PSF P(x) is given by

ν(x)=12πxP(x)=xδλz
where ∠ is the complex argument (or angle). The curve ν(x) is a line with a slope inversely proportional to z and displaced by a lateral translation δ. Because of the non-stationary behavior of P(x) and the discreteness of the transforms, the energy won’t be concentrated at a single coefficient, but rather leak and spread out to neighboring ones. Therefore, a group of STFT coefficients centered around the line ν(x) for every block are selected for calculation to maximize the energy capture, as shown in blue on Fig. 3.

The leftmost STFT coefficient index jx (or equivalently jy) is taken so that the subblock is centered as closely around the expected peak as possible:

jx=max(0,min(Nn,Nn+12Np(bxδ)λz))
clamping the value to ensure valid indices for any point position, and using the floor operator ⌊·⌋ (rounding down). The value bx is the x-coordinate of the STFT block, which is a multiple of N · p (for pixel pitch p). An example is shown on Fig. 4, which also confirms the contiguousness of the coefficients: the energy is localized in a small group of coefficients, leaving all others to be close to zero.

 figure: Fig. 4

Fig. 4 Absolute value of STFT coefficients for a single PSF. The reference PSF is on the left, the proposed sparse PSF is on the right. A single block is magnified for both algorithms in the middle. Note how the highest-magnitude coefficients are concentrated in a single spot for every N × N block. Only the coefficients located in the n × n subblock at position (jx, jy) are calculated.

Download Full Size | PDF

Note that the proposed algorithm will also naturally suppress aliasing when the PSF frequency surpasses the Nyquist rate because of potentially too steep incidence angles of light, since it won’t compute the aliased coefficients. Furthermore, changing the STFT block size N will trade-off time and frequency resolution. The best block size to maximize sparsity will largely depend on the distribution of point distances to the hologram plane, as analyzed more in detail in [28].

Finally, the block-wise inverse FFT should be performed after accumulating all PSF contributions to obtain the final hologram. This has a computational complexity of 𝒪(t log N) for t hologram pixels and block size N. Because N is a relatively small constant, the inverse STFT can be executed very fast. Moreover, individual blocks will generally completely fit in the cache, contrary to global transforms spanning over the whole hologram, further improving speed.

3. Experiments

This section consists of two parts. Firstly, we evaluate the accuracy of the proposed algorithm w.r.t. the reference approach, i.e. calculating the hologram in the spatial domain by directly evaluating Eq. (2). Secondly, we measure the calculation times and report the speedup and compare rendered views from the generated holograms.

3.1. Measuring the accuracy of the sparse STFT algorithm

For this experiment, a PSF is computed with dimensions of 512 × 512 pixels, p = 4 μm, λ = 633 nm and z = 3 cm. The proposed algorithm uses blocks of size N = 32, and is tested for subblock sizes ranging from n = 2 up to n = 32.

In Fig. 5, several renditions are shown for different values of n, besides a reference PSF. Moreover, Fig. 6 depicts a few examples of sparse PSFs for n = 4 at distances z = 2 cm up to z = 8 cm. At higher n, the PSFs are visually indistinguishale from the reference, validating the adequacy of the method. At lower subblock size, block edge defects become increasingly visible, due to the properties of the block-wise DFT. These are reminiscent of JPEG artifacts at high compression rates. Nonetheless, the circular shape of the PSF is still well-preserved even at low n values. To quantify the magnitude of the error, the Normalized Mean Square Error (NMSE) is reported as well, namely

NMSE(R,S)=jRjSj2jRj2
where R is the original reference PSF and S is the sparsified PSF. The subscript j denotes the pixel index, and the norm is Euclidean. On Fig. 7 the NMSE values are plotted for different values of n using the same settings. As expected, the error decreases progressively as n increases. Note that the error for n = 32 is not zero (≈ 3.0 · 10−3), because of the approximation introduced by (x) in section 2.1, and also because the DFT is an approximation of the analytical (continuous) Fourier transform.

 figure: Fig. 5

Fig. 5 Various renderings of a PSF placed at distance z = 3 cm using different algorithm settings. The subfigures denoted with “reference” are rendered using the standard PSF expression P(x, y) in the spatial domain, while the others are generated with the proposed approximate STFT algorithm. n denotes the subblock size and s the corresponding sparsity. ℜ means the real part is shown, φ indicates that the argument or (wrapped) phase is shown. Note how the PSF progressively resembles a phase-added stereogram more as n decreases.

Download Full Size | PDF

 figure: Fig. 6

Fig. 6 Example sparse PSF renderings at various distances z ranging from 2 to 8 cm, showing the wrapped phase for a constant sparsity of s ≈ 1.56% (n = 4).

Download Full Size | PDF

 figure: Fig. 7

Fig. 7 NMSE for various subblock dimensions n compared to a reference PSF (N = 32).

Download Full Size | PDF

Because conventional point-cloud CGH algorithms are a linear combination of many individual PSFs, the distortion will be an accumulation of all PSF errors as well. Therefore these results largely generalize for multiple points or even complete 3D point cloud models.

The gain in calculation time is reported for a collection of M = 105 points in Fig. 8, for n going from 2 up to 32. This is compared to the reference approach, taking up 557.9 seconds. As expected, the calculation time increases progressively as n grows larger. At n = 2, the calculation time is 9.8 s (a 56.8-fold speedup).

 figure: Fig. 8

Fig. 8 Calculation time for computing the hologram of the “Ship” point cloud, for various subblock dimensions n (for N = 32) compared to the reference method, shown by the full horizontal line on the graph.

Download Full Size | PDF

When n ≥ 30, the reference method is faster, mainly because the expression for calculating STFT coefficients (Eq. (7)) is more complex than the standard expression (Eq. (1)).

3.2. Generating holograms of 3D objects

In this subsection, the aim is to compare how directly calculating the hologram of 3D objects in the spatial domain (the reference approach) compares to the proposed sparse algorithm. The reference implementation utilized the separability of the Fresnel transform instead of directly calculating P(x, y): first, the PSF was evaluated to compute the horizontal and vertical components, followed by an outer product to accumulate all hologram pixel values.

The holograms were computed on a machine with an Intel Core i7 6700K CPU, 32GB RAM and ran using a C++17 implementation compiled with Microsoft Visual C++ 2017 on the OS Windows 10. The FFT was calculated using the FFTW3 library. STFT blocks of size N = 32 were used for the holograms of both 3D scenes.

The first 3D scene was configured as shown on Fig. 9. To generate a point cloud, the “Ship” model surface was randomly sampled to create a point cloud consisting of M = 105 points. The hologram dimensions were 2048 × 2048 pixels, with p = 4 μm and λ = 633 nm. The subblock size was n = 4, leading to a sparsity of s = (n/N)2 ≈ 1.56%.

 figure: Fig. 9

Fig. 9 Diagram of the “Ship” scene geometry. (a) front view, showing the lateral dimensions of the model. (b) top view, showing the distance to the hologram plane and the depth extent.

Download Full Size | PDF

The reference implementation took 557.9 seconds, while the proposed implementation only took 26.4 seconds. We thus observe a 21-fold speedup over the default implementation. The STFT only took 13 ms (averaged over 10 runs), which is negligible compared to the total calculation time. Examples of rendered views can be seen on Fig. 10, showing that the overall quality is preserved from multiple viewpoints, even though we calculated only 1.56% if the total amount of coefficients.

 figure: Fig. 10

Fig. 10 Three views taken from the holograms generated by the reference algorithm (top row) and the proposed algorithm (bottom row). The views were rendered by taking a vertically centered 1024 × 1024 subhologram crop at the left side, middle and right side of the hologram, followed by a backpropagation with z = −10.5 cm (using the ASM) and taking the absolute value.

Download Full Size | PDF

For the second 3D scene, a point cloud of the “Train” model was used (cf. Fig. 11), consisting of M = 2 · 105 points assigned with random phases. The hologram dimensions were also 2048 × 2048 pixels, with p = 3 μm and λ = 633 nm. In this case, several different values for n were used to compare differences in visual quality.

 figure: Fig. 11

Fig. 11 Diagram of the “Train” scene geometry. (a) front view of the scene, including lateral dimensions. (b) top view, showing how the train is tilted w.r.t. the z-axis.

Download Full Size | PDF

The image quality can be measured using the Peak Signal-to-Noise Ratio (PSNR), defined by

PSNR(R,S)=10log10(CRmax2jC(RjSj)2)
where R is the reference image and S is the reconstruction from the sparsified representation; C is the total pixel count, j is the pixel index and Rmax is the maximum pixel value, which is 255 for 8-bit images.

The results are reported in Table 1; the calculation times for the reference method is 1099.1 s. The rendered images are displayed in Fig. 12: almost no deterioration can be seen for n ≥ 4, but some noise and loss of detail is visible for n = 2. Taking n = 4 as the baseline, we obtain a speedup factor of 20.6.

Tables Icon

Table 1. “Train” hologram calculation times for different values of n, and PSNR w.r.t. the unsparsified representation for views focused at the front and back of the scene.

 figure: Fig. 12

Fig. 12 Reconstructed images of the “Train” hologram, for various values of subblock size n. The front focus is shown by taking the absolute value after backpropagating at z = −10.3 cm with the ASM. The back focus uses z = −11.2 cm.

Download Full Size | PDF

4. Conclusion

This paper presents an algorithm to analytically calculate the STFT coefficients of PSFs, thereby enabling the sparse CGH needing only a small subset of the total amount of coefficients, without requiring a LUT. A speedup factor of 20 can be achieved w.r.t. the reference implementation with limited quality degradation.

Because of the generality of Fresnel diffraction, this algorithm has use beyond holography for other applications requiring near-field numerical diffraction of 3D point clouds. It can be viewed as a generalization of the phase-added stereogram, trading off accuracy and calculation speed.

Furthermore, this approach is orthogonal to many other CGH acceleration algorithms. Future work will involve combining this algorithm with WRPs, other input elements besides PSFs, and further analysis on tuning of the algorithm’s parameters for optimizing quality and speed for any given 3D scene. Also, because of the low memory requirements, it is a good candidate for a parallelized implementation on specialized hardware such as FPGAs.

Funding

European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013) (617779); Research Foundation - Flanders (FWO) (G024715N).

Acknowledgments

The 3D models are courtesy of TurboSquid (http://www.turbosquid.com): Low Poly Ship (by user “amingin”) and T2-71 Steam Locomotive (by user “Locomotive_works”).

References

1. D. Blinder, A. Ahar, S. Bettens, T. Birnbaum, A. Symeonidou, H. Ottevaere, C. Schretter, and P. Schelkens, “Signal processing challenges for digital holographic video display systems,” Signal Process. Image 70, 114–130 (2019). [CrossRef]  

2. J.-H. Park, “Recent progress in computer-generated holography for three-dimensional scenes,” J. Inf. Displ. 18, 1–12 (2017). [CrossRef]  

3. M. E. Lucente, “Interactive computation of holograms using a look-up table,” J. Electron. Imaging 2, 28–34 (1993). [CrossRef]  

4. S.-C. Kim and E.-S. Kim, “Effective generation of digital holograms of three-dimensional objects using a novel look-up table method,” Appl. Opt. 47, D55–D62 (2008). [CrossRef]   [PubMed]  

5. T. Shimobaba, N. Masuda, and T. Ito, “Simple and fast calculation algorithm for computer-generated hologram with wavefront recording plane,” Opt. Lett. 34, 3133–3135 (2009). [CrossRef]   [PubMed]  

6. P. Tsang, W.-K. Cheung, T.-C. Poon, and C. Zhou, “Holographic video at 40 frames per second for 4-million object points,” Opt. Express 19, 15205–15211 (2011). [CrossRef]   [PubMed]  

7. S. Jiao, Z. Zhuang, and W. Zou, “Fast computer generated hologram calculation with a mini look-up table incorporated with radial symmetric interpolation,” Opt. Express 25, 112–123 (2017). [CrossRef]   [PubMed]  

8. Y. Yamamoto, H. Nakayama, N. Takada, T. Nishitsuji, T. Sugie, T. Kakue, T. Shimobaba, and T. Ito, “Large-scale electroholography by horn-8 from a point-cloud model with 400,000 points,” Opt. Express 26, 34259–34265 (2018). [CrossRef]  

9. A. Symeonidou, D. Blinder, A. Munteanu, and P. Schelkens, “Computer-generated holograms by multiple wavefront recording plane method with occlusion culling,” Opt. Express 23, 22149–22161 (2015). [CrossRef]   [PubMed]  

10. 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, 25440–25449 (2015). [CrossRef]   [PubMed]  

11. A. Symeonidou, D. Blinder, and P. Schelkens, “Colour computer-generated holography for point clouds utilizing the phong illumination model,” Opt. Express 26, 10282–10298 (2018). [CrossRef]   [PubMed]  

12. A. Gilles and P. Gioia, “Real-time layer-based computer-generated hologram calculation for the fourier transform optical system,” Appl. Opt. 57, 8508–8517 (2018). [CrossRef]   [PubMed]  

13. H. Kim, J. Hahn, and B. Lee, “Mathematical modeling of triangle-mesh-modeled three-dimensional surface objects for digital holography,” Appl. Opt. 47, D117–D127 (2008). [CrossRef]   [PubMed]  

14. K. Matsushima and S. Nakahara, “Extremely high-definition full-parallax computer-generated hologram created by the polygon-based method,” Appl. Opt. 48, H54–H63 (2009). [CrossRef]   [PubMed]  

15. 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, 24450–24465 (2014). [CrossRef]   [PubMed]  

16. H. Nishi and K. Matsushima, “Rendering of specular curved objects in polygon-based computer holography,” Appl. Opt. 56, F37–F44 (2017). [CrossRef]   [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, 9192–9197 (2013). [CrossRef]   [PubMed]  

18. T. Shimobaba, T. Kakue, and T. Ito, “Acceleration of color computer-generated hologram from three-dimensional scenes with texture and depth information,” Proc.SPIE 9117, 91170B (2014).

19. T. Yatagai, “Stereoscopic approach to 3-d display using computer-generated holograms,” Appl. Opt. 15, 2722–2729 (1976). [CrossRef]   [PubMed]  

20. M. Yamaguchi, H. Hoshino, T. Honda, and N. Ohyama, “Phase-added stereogram: calculation of hologram using computer graphics technique,” Proc. SPIE 1914, 25–31 (1993). [CrossRef]  

21. 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, 4246–4255 (2009). [CrossRef]   [PubMed]  

22. K. Wakunami, H. Yamashita, and M. Yamaguchi, “Occlusion culling for computer generated hologram based on ray-wavefront conversion,” Opt. Express 21, 21811–21822 (2013). [CrossRef]   [PubMed]  

23. H. Zhang, Y. Zhao, L. Cao, and G. Jin, “Fully computed holographic stereogram based algorithm for computer-generated holograms with accurate depth cues,” Opt. Express 23, 3901–3913 (2015). [CrossRef]   [PubMed]  

24. A. Gilles, P. Gioia, R. Cozot, and L. Morin, “Hybrid approach for fast occlusion processing in computer-generated hologram calculation,” Appl. Opt. 55, 5459–5470 (2016). [CrossRef]   [PubMed]  

25. P. Schelkens, A. Skodras, and T. Ebrahimi, The JPEG 2000 Suite (Wiley Publishing, 2009). [CrossRef]  

26. J. W. Goodman, Introduction to Fourier Optics (W. H. Freeman and Company, 2017).

27. T. Shimobaba and T. Ito, “Fast generation of computer-generated holograms using wavelet shrinkage,” Opt. Express 25, 77–87 (2017). [CrossRef]   [PubMed]  

28. D. Blinder and P. Schelkens, “Accelerated computer generated holography using sparse bases in the STFT domain,” Opt. Express 26, 1461–1473 (2018). [CrossRef]   [PubMed]  

29. D. Arai, T. Shimobaba, T. Nishitsuji, T. Kakue, N. Masuda, and T. Ito, “An accelerated hologram calculation using the wavefront recording plane method and wavelet transform,” Opt. Commun. 393, 107–112 (2017). [CrossRef]  

30. H. Kang, T. Fujii, T. Yamaguchi, and H. Yoshikawa, “Compensated phase-added stereogram for real-time holographic display,” Opt. Eng. 46, 095802 (2007). [CrossRef]  

31. H. Kang, E. Stoykova, and H. Yoshikawa, “Fast phase-added stereogram algorithm for generation of photorealistic 3d content,” Appl. Opt. 55, A135–A143 (2016). [CrossRef]   [PubMed]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (12)

Fig. 1
Fig. 1 MSE of the approximations for small values (ES(x)) and large values (EL(x)) compared to the reference E(x), on a logarithmic scale. The crossing point of both curves (at x ≈ 1.609) is chosen to select the approximation with the smallest error depending on x. Given the symmetry of E(x), the curves look exactly the same (but mirrored) for negative values of x.
Fig. 2
Fig. 2 Comparison of the (imaginary parts of) E(x) and (x). The curves show good agreement (see Fig. 1 for the error analysis).
Fig. 3
Fig. 3 Time-frequency charts of PSF curves superimposed on STFT coefficients. The red curves represent the instantaneous frequency of the two PSF instances, whose slope will depend on z and lateral displacement on δ. The (Heisenberg) boxes correspond to individual coefficients, indicating where their energy is concentrated in space (or time) and frequency. The blue boxes are the coefficients that will be computed for this PSF. Both (a) and (b) have the same sparsity rates s, but (b) has a larger block size N trading off more frequency resolution for less spatial resolution, which is better for points with large |z| values.
Fig. 4
Fig. 4 Absolute value of STFT coefficients for a single PSF. The reference PSF is on the left, the proposed sparse PSF is on the right. A single block is magnified for both algorithms in the middle. Note how the highest-magnitude coefficients are concentrated in a single spot for every N × N block. Only the coefficients located in the n × n subblock at position (jx, jy) are calculated.
Fig. 5
Fig. 5 Various renderings of a PSF placed at distance z = 3 cm using different algorithm settings. The subfigures denoted with “reference” are rendered using the standard PSF expression P(x, y) in the spatial domain, while the others are generated with the proposed approximate STFT algorithm. n denotes the subblock size and s the corresponding sparsity. ℜ means the real part is shown, φ indicates that the argument or (wrapped) phase is shown. Note how the PSF progressively resembles a phase-added stereogram more as n decreases.
Fig. 6
Fig. 6 Example sparse PSF renderings at various distances z ranging from 2 to 8 cm, showing the wrapped phase for a constant sparsity of s ≈ 1.56% (n = 4).
Fig. 7
Fig. 7 NMSE for various subblock dimensions n compared to a reference PSF (N = 32).
Fig. 8
Fig. 8 Calculation time for computing the hologram of the “Ship” point cloud, for various subblock dimensions n (for N = 32) compared to the reference method, shown by the full horizontal line on the graph.
Fig. 9
Fig. 9 Diagram of the “Ship” scene geometry. (a) front view, showing the lateral dimensions of the model. (b) top view, showing the distance to the hologram plane and the depth extent.
Fig. 10
Fig. 10 Three views taken from the holograms generated by the reference algorithm (top row) and the proposed algorithm (bottom row). The views were rendered by taking a vertically centered 1024 × 1024 subhologram crop at the left side, middle and right side of the hologram, followed by a backpropagation with z = −10.5 cm (using the ASM) and taking the absolute value.
Fig. 11
Fig. 11 Diagram of the “Train” scene geometry. (a) front view of the scene, including lateral dimensions. (b) top view, showing how the train is tilted w.r.t. the z-axis.
Fig. 12
Fig. 12 Reconstructed images of the “Train” hologram, for various values of subblock size n. The front focus is shown by taking the absolute value after backpropagating at z = −10.3 cm with the ASM. The back focus uses z = −11.2 cm.

Tables (1)

Tables Icon

Table 1 “Train” hologram calculation times for different values of n, and PSNR w.r.t. the unsparsified representation for views focused at the front and back of the scene.

Equations (17)

Equations on this page are rendered with MathJax. Learn more.

P ( x , y ) = exp ( i γ 2 [ ( x δ ) 2 + ( y ) 2 ] )
H ( x , y ) = j = 1 M P j ( x , y ) = j = 1 M exp ( i γ j 2 [ ( x δ j ) 2 + ( y j ) 2 ] )
H ( x , y ) = 𝒯 1 { 𝒯 { H ( x , y ) } } = 𝒯 1 { j = 1 M 𝒯 { P j ( x , y ) } }
f ( ω , η ) = B + B B + B P ( x , y ) exp ( i ( x ω + y η ) ) d x d y
f x ( ω ) = B + B exp ( i [ γ 2 ( x δ ) 2 + x ω ] ) d x and f y ( η ) = B + B exp ( i [ γ 2 ( y ) 2 + y η ] ) d y
erf ( x ) = 2 π 0 x e t 2 d t
1 + i 4 γ 2 π exp ( i ( ω δ + ω 2 4 γ 2 ) ) [ erf ( 1 i 2 ( ( B + δ ) γ + ω 2 γ ) ) + erf ( 1 i 2 ( ( B δ ) γ ω 2 γ ) ) ] .
f x ( ω ) 1 γ exp ( i ( ω δ + α 2 ) ) [ E ( B γ + β ) + E ( B γ β ) ]
E ( x ) = c + exp ( i x 2 ) x i exp ( i x 2 ) 2 x 3 3 exp ( i x 2 ) 4 x 5 + 𝒪 ( 1 x 7 ) . ( for x + )
E L ( x ) = exp ( i x 2 ) x + c sgn ( x )
{ E ( x ) } { E S ( x ) } = 2 3 x 3 + 1 21 x 7 1 660 x 11
{ E ( x ) } { E S ( x ) } = 2 x 1 5 x 5 + 1 108 x 9
E ˜ ( x ) = { E S ( x ) , if | x | < 1.609 E L ( x ) , otherwise .
ν ( x ) = 1 2 π x P ( x ) = x δ λ z
j x = max ( 0 , min ( N n , N n + 1 2 N p ( b x δ ) λ z ) )
NMSE ( R , S ) = j R j S j 2 j R j 2
PSNR ( R , S ) = 10 log 10 ( C R max 2 j C ( R j S j ) 2 )
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.