Analytic computation of line-drawn objects in computer generated holography

: Digital holography is a promising display technology that can account for all human visual cues, with many potential applications i.a. in AR and VR. However, one of the main challenges in computer generated holography (CGH) needed for driving these displays are the high computational requirements. In this work, we propose a new CGH technique for the eﬃcient analytical computation of lines and arc primitives. We express the solutions analytically by means of incomplete cylindrical functions, and devise an eﬃciently computable approximation suitable for massively parallel computing architectures. We implement the algorithm on a GPU (with CUDA), provide an error analysis and report real-time frame rates for CGH of complex 3D scenes of line-drawn objects, and validate the algorithm in an optical setup.


Introduction
Computer generated holography (CGH) consists of numerical diffraction algorithms for efficiently calculating interference for various applications in holography, such as display technology [1], beam shaping [2] and interferometry [3], among many others.One of the main computational challenges is that every (virtual) scene element can potentially affect all hologram pixels.This problem is exacerbated for holographic video displays, requiring high-resolution hologram generation at video rates.
Holographic displays are often cited as the ultimate 3D display technology [4].In principle, holograms can display virtual objects which are indistinguishable from real objects since they account for all human visual cues.This makes them promising candidates for augmented and virtual reality applications, since they do not suffer from the drawbacks of head-mounted displays with conventional screens, as they typically cannot account for eye-focusing and exhibit accommodation-vergence conflicts.
In this work, we investigate the problem of diffracting line-drawn segments.We target applications such as head-mounted displays, near-eye displays, navigation systems; but this algorithm could be used as a building block for CGH of more complex virtual scenes too.This may also have broader applications in pure diffraction theory beyond displays, e.g. the efficient computation of straight or curved line apertures.This problem was recently addressed in [21], which we designate as the reference "CG line" (computer graphics) CGH method.It subdivides the scene into depth planes, utilizing a pre-computed one-dimensional signal per plane encoding a converged light distribution.For every curve, the buffer is added on the hologram plane for each of a series of equidistant points along the curve with a period equal to the pixel pitch.This buffer is sweeped along the curve, drawn perpendicularly to the curve's tangent at every point.The algorithm was successful in creating line-drawn object CGH at multiple depths in real-time.
However, the current version of that method has a few drawbacks.It is primarily designed for shorter distances to the hologram plane and operates within a fixed set of depth planes.Furthermore, noticeable errors are present at the curve's edges, especially at lower resolutions.
To address this, we introduce a novel way for computing the hologram of line-drawn objects.We propose an analytical technique for directly computing the value of diffracted line and circle arc apertures in any point (e.g.pixel).This eliminates the error present at edges, and imposes no constraints on the number of depth levels, up to one per line if needed.Moreover, this analytical formulation makes a parallel implementation more straightforward: we achieve video-rate CGH for driving a HD spatial light modulator (SLM) with complex 3D line-drawn content on a GPU.
The paper is structured as follows: in sections 2 and 3 we derive analytical expressions for the diffraction of curve apertures.We begin with the straight line case, and proceed with the derivation of general circular arc segments.We derive an efficiently computable numerical approximation and provide an error analysis.Then, in section 4, we compare the calculation speed and accuracy with multiple CGH algorithms, including the reference point-based CGH and the CG line CGH algorithms.We examine the visual quality in both simulation environment and validate it in a holographic display setup.Finally, we conclude in section 5.

Line apertures
The expression for a single point-spread function P, describing the diffraction pattern induced by a luminous point at coordinates (δ, , ζ), is given by where ζ is the distance to the hologram plane, i is the imaginary unit, λ is the wavelength and a ∈ C is the point amplitude.This expression will be generalized for diffracting line-drawn segments of two kinds: straight line segments and arbitrary circle segment apertures, or arcs.The diffraction pattern for an arc aperture A, cf.Fig. 1(a), can be described by the sum of the set of all P(x, y) points lying on the arc with radius ρ and spanning an angle α, and is given by The integral is multiplied by ρ to get a constant intensity per unit arc length; e.g. when ρ doubles, the integral "sweeps" twice as fast over the arc, requiring doubling the integral to exactly compensate for the energy loss.For notational simplicity, we assume the arc to be centered at the origin and anchored in the x-axis.The solution at the end can be easily adapted for arbitrary arcs later by a change of coordinates.Similarly, we consider a straight line segment parallel lying on the x-axis, cf.Fig. 1(b), beginning at the origin with length : The goal is to solve these integrals and find an efficiently computable analytical expression for the real-time CGH calculation of line-drawn objects.Any curve shape can be approximated to arbitrary precision with lines and arcs, such as arc splines [22].We will begin with the special case of a straight line, whose solution will serve as a basis for the computation of arbitrary arcs.

Straight line apertures
We rewrite (3) as where erf(c) is the error function.Note that the input argument c to erf(c) is complex-valued, so we cannot rely on standard built-in functions and algorithms.Although c is complex, we only have to consider input values along the main diagonals of the complex plane, since c = (1 − i)s, for some s ∈ R.This function is therefore closely related to the Fresnel integrals, for which many numerical solutions exist such as table-based calculations.We will use the same approximation function E(t), t ∈ R as derived in [11]: where E S (t) and E L (t) are valid for small and large values of |t|.E S (t) was derived using a Taylor expansion, in separate real ( ) and imaginary ( ) parts: while E L (t) was obtained by truncating an asymptotic expansion: where sgn(t) is the sign function, equal to 1 when t is positive and equal to −1 when t is negative.Using E(t), we can rewrite (4) and obtain the final expression for straight line segment diffraction, where γ = π λζ .

Arc apertures
We convert expression (2) to polar coordinates (x, y) = (r cos θ, r sin θ): which can be simplified further using the identity cos φ cos θ + sin φ sin θ = cos (φ − θ): The main remaining hurdle is to solve an integral of the form where w and s are arbitrary real-valued parameters of a two-dimensional function we want to be able to efficiently compute.Unfortunately, this integral cannot be expressed in terms of elementary functions.They are closely related to the non-homogeneous Bessel's differential equations.We start from the expressions [23] where J ν (s) is the Bessel function of the first kind of order ν, H ν (s) is the Struve function of order ν and Γ is the gamma function.Using Euler's formula exp(is) = cos(s) + i sin(s) and taking the order ν = 0, we obtain for the fixed integration bounds w = π 2 , π To obtain a solution for all w, we will need to further analyze E(s, w), which we will call an "incomplete cylindrical function", following the naming convention in [24].Although E(s, w) can be evaluated ∀s, w ∈ R, we can reduce the space of values under consideration by leveraging the following symmetries: where the overbar stands for complex conjugation.Thus we will only consider s>0 and 0 ≤ w< π 2 for the numerical model, since values for s and w outside of these bounds can always be reduced using the symmetries of (16).To efficiently numerically approximate E(s, w), we will consider two cases: large and small values of w √ s, which will be expounded on in the subsections below.

Asymptotic approximation for large parameters
Substituting u = cos w in E(s, w) and using (15), we obtain The latter integral can be calculated [25] by using the Fourier series expansion of (1 − u 2 ) − 1 2 : Although this series can be used to compute E(s, w) to arbitrary precision for large s, the convergence is too slow for our intended application, requiring too many calculations.
Instead, we propose a different approach.We use the following equality: showing that the left-hand side of ( 19) is close to an elementary function solution of the sought integral (18), except for the second term in the right-hand side.But because that term has a s in the denominator, its contribution will be small for large s.We thus have Finally, we would like to find a way to compute H 0 (s); unlike Bessel functions, there are no standard C++/CUDA implementations of Struve functions at the time of writing.We therefore use the following asymptotic form for Struve functions: where Y ν (s) is the Neumann function of order ν, also known as a Bessel function of the second kind, which is easier to calculate.Substituting this for ν = 0 gives us: this first error term will exactly compensate the i s term derived in (20), resulting in the final expression valid for larger values of s and sin w, where H 0 (s) = J 0 (s) + iY 0 (s) is the Hankel function of order 0.

Series approximation for small parameters
For small values of w √ s, we need a different method.A natural first approach is to take the Taylor expansion in s of the expression within the integral of E(s, w), which can be integrated term by term to obtain Unfortunately, this expression is only useful for very small values of |s|<10 −2 ; (25) cannot capture the oscillatory behavior of E(s, w) with a small number of terms.
We propose a different approach by taking the Taylor expansion of the exponent instead: valid for small values of w ≈ sin w, after integration Note that this is exactly the same kind of restricted domain error function we approximated with (5) for straight line apertures.Using this function, we get the final expression for E S (s, w), only valid for small values of w √ s:

Error analysis
To know whether to choose E S (s, w) or E L (s, w) for any given (s, w), we have to quantify which one has the smallest error.For this we use the following measure of relative precision: where a is the reference value and â is the approximation.The relative precision of E S and E L w.r.t.E are shown in Figs.2(a)-2(b) and their combined maximum in Fig. 2(c).In Fig. 2(d), we evaluate whether RP (E, E L )> RP (E, E S ) or not for all tested pairs of (s, w).There is a clear separation in two regions, so it suffices to know on what side of the divide a point (s, w) lies, for which we want to find a simple expression.Using linear regression with the curve w √ s = c (where c is the parameter to be fitted), we get that with a goodness-of-fit R 2 = 0.9985 when tested in the region s ∈ ]0, 10 5 ].
Finally, we have to consider the case s ≈ 0, which will cause numerical problems since s appears in the denominator for expressions in both E S and E L .For this we can use the (truncated) Taylor expansion in (24); in conclusion we get otherwise. (31)

Evaluating arc aperture expressions
With the help of the incomplete cylindrical function E(s, w), we can rewrite (11) as where s = −2γ 2 rρ and γ = π λζ as defined for (9).However, we approximated E(s, w) only for s ≥ 0 and 0 ≤ w ≤ π 2 , so we need to extend the domain of E(s, w) using the symmetries in ( 16).Extending for s is more straightforward: the sign of s only depends on γ 2 , whose sign is equal to that of the depth position ζ.Typically, the sign of ζ for all elements in a scene will be the same, so we can fill in |s| instead of s and take the complex conjugate of the final hologram if s is negative.But this can be done on a per-line-segment basis as well if needed.
The case for w is more tricky, as the input can be many multiples of π 2 , including negative ones.Recursively applying the two last symmetries of ( 16) would be computationally sub-optimal, as this would introduce branching divergence across computation threads, which is unfavorable for parallel computation implementations.Instead, we combine all cases for w using the following expression: where • is the floor operator (rounding down), with where mod stands for the modulo operation.Thus, "bitsgn (t)" is a type of sign function returning a value based on the least significant bit of an integer.In practice, this means we compute E(s, |w − mπ|) once using e.g.floating point numbers, and we xor the sign bits of the real and imaginary part with bits 0 and 1 of integer n, respectively.The value of πJ 0 (s) only has to be computed once per evaluated value of A(r, θ) and can be re-used for both evaluations of E(s, w) as well as in the approximated case E L (s, w).

Experiments
The experiments are partitioned in two parts.First, we perform numerical experiments to compare the calculation time and accuracy of the proposed algorithm with other references.Then, we test the algorithm in an optical setup and compare visual quality differences.

Numerical calculation experiments
We utilized two line segment data sets for our experiments: (1) the "Simple shapes" data set (Fig. 3), consisting of 40 straight lines and 6 circular segments, forming simple geometric shapes placed at multiple depth planes; (2) the "Seigaiha" data set (Fig. 4) made out of 144 arc segments, which is a traditional Japanese motif of superposed wave patterns.We implemented multiple algorithms for computing line holograms.
• Point-based: this is the exact reference algorithm, where we decompose all lines into points with a period equal to the pixel pitch.We sum over all point-spread functions according to (1).
• FFT-based: this will subdivide the scene into D depth planes depending on the data set, directly draw the lines in their respective depth planes, followed by an FFT-based • CG line method: this is the method proposed in [21], using a 1D line buffer which are drawn over the hologram plane along the curve.
• Analytical method: the proposed algorithm used in this paper.
The algorithms may be implemented using multi-threading on a CPU, or using massively parallel processing on a GPU.The multi-threaded CPU implementations ran on a machine with an Intel Core-i7 8850H 2.60GHz with 16GB of RAM, using either Visual Studio's "concurrency::parallel_for" or OpenMP with maximal CPU core utilization.The GPU versions of the algorithm were run on a machine with an Intel Xeon E5-2687W v4 processor 3Ghz, 256 GB of RAM, a NVIDIA TITAN RTX GPU with CUDA 11.0, enabling CUDA compute capability 7.5.All versions of the algorithm were implemented in C++17 running on Windows 10 OS, using 32-bit floating point precision.
The calculation times for all tested algorithms applied on both data sets are summarized in Table 1.The reference point-based algorithm is embarrassingly parallel, and thus easily implemented on a GPU; it will have calculation times proportional to the number of points, which will depend in turn on the number of line segments and their lengths.This explains why the "Seigaiha" with more line elements takes about half a second, which is 3 times longer than the "Simple shapes".On the other hand, the FFT-based algorithm's time will mostly depend on the number of required FFT operations, proportional to the number of depth planes.Unlike all other tested algorithms, this one does not scale linearly with the hologram resolution.
The CG method was implemented in a multi-threaded CPU environment.Direct implementation of this algorithm on GPU is not straightforward, as the drawn line buffers are self-intersecting in the hologram plane introducing dependencies, and requiring more complex memory access patterns.A successful GPU implementation may yield further significant speedups.
The proposed analytical method can compute pixel values independently, making it highly suitable for GPU: all threads can execute the same operations with built in trigonometric and Bessel function operations.A functional, but non-optimized implementation on CPU is also provided for reference.The calculation time is proportional to the number of line segments, yet independent of their parameters.We report calculation times of 2ms and 21ms for "Simple shapes" and "Seigaiha" respectively, which ranges from 20x to over 50x faster than the point-based CGH reference algorithm.
Furthermore, we want to quantify and compare the numerical accuracy of the different algorithms for the tested data sets.We utilize two quality metrics: the peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM).The used quality metrics require a reference signal H, measuring the degree by which an approximated signal Ĥ differs from the reference.The general definition of the PSNR for a real-/ complex-valued signal H is PSNR(H, Ĥ) := 10 log 10 max (|H|) where • 2 is the Euclidean norm.The PSNR was directly applied on the complex-valued holographic signal and the results are reported in Table 1.
The SSIM uses the definition found in [26], which has been shown to significantly outperform PSNR for measuring perceptual image quality.Thus we also report both metrics on the numerically backpropagated hologram amplitude at multiple different depth planes where some of the objects are in focus.This is done for the proposed analytical CGH algorithm, the FFT-based and the CG-based CGH algorithms, using the Point-based CGH as a reference.The results for the "Simple shapes" are reported in Table 2 and for "Seigaiha" are reported in Table 3.The algorithm outperforms the CG line algorithm on average by 3.4 dB PSNR (and 0.27 SSIM units) for "Simple shapes" and by 16.1 dB PSNR (and 0.74 SSIM units) for "Seigaiha".This gain difference is likely attributable to the stronger approximation of the CG method at line edges, which are more prominent in "Seigaiha".The FFT-based algorithm performs slightly better than the CG line algorithm overall.However, please note that the current version of the FFT-based algorithm draws 1-pixel thick lines using nearest pixel neighbors.The quality may be improved considerably by using anti-aliasing.Moreover, note that the analytical method, unlike the other methods, requires a decomposition of curves into lines and arcs, so this may affect the resulting curve shape depending on the chosen approximation [22].
Finally, we want to confirm that the proposed analytical algorithm has consistently better quality for any arc radius, not just for specific cases.To do so we measure the approximation error using the the normalized mean-square error (NMSE) testing the values by varying the arc radius from 160 up to 6000 pixels.The holograms were computed at 0.1m and at 0.8 for both the CG line and the analytical CGH versions and graphed in Fig. 5.The NMSE error curves for the analytical CGH are always strictly lower than the CG line CGH counterpart, on average by 0.04 units for the 0.1m case and 0.15 units for the 0.8m case.

Optical experiments
We also verify the proposed algorithms using optically reconstructed image of the tested holograms.We used a phase-modulation-type SLM (Holoeye Photonics AG, "PLUTO") and a green laser with 532 nm wavelength (Showa Optronics, "J150GS", Japan), cf. Figure 6.The numerical and optical reconstructions of the "Simple shapes" and "Seigaiha" data sets are shown in Tables 4 and 5, respectively.For "Simple shapes", we showed cropped zoomed-in images of the holograms brought into focus different depth planes.The numerical reconstructions look largely the same for all three methods, but there is some intensity loss at the line's edges for the CG line method as previously observed in [21].These optical reconstructions mirror the numerical reconstructions quite well.For "Seigaiha", we show reconstructions of the full hologram refocused in the middle of the virtual depth volume at 0.5m.The analytical method's reconstructions resemble the reference Point-based CGH closer than the CG line method, whose edge darkening effect is visible just like in the other data set.
The visible noise in the optical reconstructions is due to several factors: the zero-order diffraction light, the lack of amplitude modulation in the SLM and the partial visibility of objects in other focal planes due to their long depth of fields.This effect could be mitigated by using i.a. a 4f optical system with a spatial filter, using SLMs with higher resolutions and modulation capabilities, time-multiplexing or rapidly changing the relative phase values of the different geometric shapes over time.Table 5. Reconstructions of the "Seigaiha" data set CGH, refocused at 0.50m.

Conclusion
We propose a novel algorithm for the fast computation of CGH for line-drawn objects.We analytically derive the solution with incomplete cylindrical functions and simplify the expressions to an efficiently computable approximation suitable for GPU implementations.We report 1 to 2 orders of magnitude speedups and higher accuracy than previous solutions, and validate the rendered CGH in an holographic display setup.This algorithm enables the real-time generation of line-drawn CGH, opening the door to novel display applications in digital holography.Future work may include adding phase modulation to the line segments, different curve primitives and extensions to 3D curves.

Fig. 1 .
Fig. 1.Diagram of luminous line drawn apertures with its parameters.The example arc (a) is centered at the origin and anchored at the x-axis, with radius ρ and spanning an angle α.The exemplary straight line (b) lies on the x-axis, begins at the origin and has length .

Fig. 2 .
Fig. 2. Error maps showing the relative precision of (a) E S (s, w) and (b) E L (s, w) w.r.t.E(s, w).The color scale effectively shows the number of correct digits of computed values, clipped between 0 and 5 for better contrast.Note how the high and low precision regions are opposed in the respective maps, shown by taking the maximum of both maps in (c), which will be the effective precision of the algorithm.(d) the transition map shows for a given point (s, w) whether E S has higher precision (black) or E L has higher precision (gray); it can be accurately modelled by the curve (w √ s = 2.26), shown by the red dashed line.

Fig. 3 .
Fig. 3. "Simple shapes" data set, seen from the front (a) and the side (b), color-coded according to the depth plane.

Fig. 4 .
Fig. 4. Diagram of the "Seigaiha" data set, seen from the front (a) and the side (b), color-coded according to the depth plane.

Fig. 5 .
Fig. 5. Graphs of the NRMSE for the CGH of a single arc line segment as a function of its radius.This was tested for both the Analytical methods and the CG line method, each placed at 0.1 m and at 0.8 m.

Fig. 6 .
Fig. 6.Images of the used holographic display setup.(a) shows a simplified diagram, (b) is an annotated photograph of the optical setup.

Table 4 .
Zoomed-in numerical and optical reconstructions of the "Simple shapes" data set CGH at various depths.The first three columns show numerical reconstructions, the latter three show the corresponding optical reconstructions.