Analytic computation of line-drawn objects in computer generated holography

David Blinder; David Blinder; Takashi Nishitsuji; Takashi Kakue; Tomoyoshi Shimobaba; Tomoyoshi Ito; Peter Schelkens; Peter Schelkens

doi:10.1364/OE.405179

1. 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.

There are many kinds of CGH algorithms [4,5], whose relative advantages and drawbacks are reflected in attributes such as computational speed, supported graphical effects and scene geometry, numerical precision and visual realism. One major way of classifying CGH algorithms is by their used elementary constituents: point cloud methods [6–11], layer-based methods [12–14], polygon methods [15–17] and ray-based methods [18–20].

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.

2. Line apertures

The expression for a single point-spread function $P$, describing the diffraction pattern induced by a luminous point at coordinates $(\delta , \epsilon , \zeta )$, is given by

(1)$$P(x,y) = a\cdot\exp\Big(\frac{\pi i}{\lambda \zeta}\big[(x-\delta)^{2} + (y-\epsilon)^{2} \big]\Big)$$

where $\zeta$ is the distance to the hologram plane, $i$ is the imaginary unit, $\lambda$ is the wavelength and $a\in \mathbb {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 $\rho$ and spanning an angle $\alpha$, and is given by

(2)$$A(x,y) = a\rho\int_0^{\alpha} \exp\Big(\frac{\pi i}{\lambda \zeta}\big[(x-\rho\cos{\phi})^{2} + (y-\rho\sin{\phi})^{2} \big]\Big) d\phi.$$

The integral is multiplied by $\rho$ to get a constant intensity per unit arc length; e.g. when $\rho$ 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 $\ell$:

(3)$$L(x,y) = a\int_0^{\ell} \exp\Big(\frac{\pi i}{\lambda \zeta}\big[(x-u)^{2} + y^{2}\big]\Big) du.$$

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.

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 $\rho$ and spanning an angle $\alpha$. The exemplary straight line (b) lies on the x-axis, begins at the origin and has length $\ell$.

Download Full Size | PDF

2.1 Straight line apertures

We rewrite (3) as

(4)$$\begin{aligned} L(x,y) &= a \exp\left(\frac{\pi i}{\lambda \zeta}y^{2}\right) \int_0^{\ell} \exp\left(\frac{\pi i}{\lambda \zeta}(x-u)^{2}\right)du \\ &= \tfrac{1}{4}a(1+i)\sqrt{2\lambda\zeta} \exp\left(\frac{\pi i}{\lambda \zeta}y^{2}\right) \left[\operatorname{erf}{\left( \sqrt{\frac{\pi}{2\lambda\zeta}}(1-i)x \right)} + \operatorname{erf}{\left( \sqrt{\frac{\pi}{2\lambda\zeta}}(1-i)(\ell-x) \right)}\right] \end{aligned}$$

where $\operatorname {erf}(c)$ is the error function. Note that the input argument $c$ to $\operatorname {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\in \mathbb {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 $\mathcal {E}(t)$, $t\in \mathbb {R}$ as derived in [11]:

(5)$$(i-1)\sqrt{\tfrac{\pi}{2}} \operatorname{erf}(\tfrac{1-i}{\sqrt{2}}t) = \mathcal{E}(t) \approx \begin{cases} \mathcal{E}_S(t), & \textrm{if}\ \lvert t \rvert < 1.609 \\ \mathcal{E}_L(t), & \textrm{otherwise}. \end{cases}$$

where $\mathcal {E}_S(t)$ and $\mathcal {E}_L(t)$ are valid for small and large values of $\lvert t \rvert$. $\mathcal {E}_S(t)$ was derived using a Taylor expansion, in separate real $(\Re )$ and imaginary $(\Im )$ parts:

(6)$$\Re\{\mathcal{E}_S(t)\} = 2t - \tfrac{1}{5}t^{5} + \tfrac{1}{108}t^{9} $$

(7)$$\Im\{\mathcal{E}_S(t)\} = -\tfrac{2}{3}t^{3} + \tfrac{1}{21}t^{7} - \tfrac{1}{660}t^{11} $$

while $\mathcal {E}_L(t)$ was obtained by truncating an asymptotic expansion:

(8)$$\mathcal{E}_L(t) = (1-i)\sqrt{\tfrac{\pi}{2}}\cdot\mathrm{sgn}(t) + \frac{i\exp\big(-it^{2}\big)}{t}$$

where $\mathrm {sgn}(t)$ is the sign function, equal to $1$ when $t$ is positive and equal to $-1$ when $t$ is negative.

Using $\mathcal {E}(t)$, we can rewrite (4) and obtain the final expression

(9)$$L(x,y) = -\frac{ia}{2\gamma}\exp\left(i\gamma^{2} y^{2}\right) \big[\mathcal{E}\left( \gamma x \right) + \mathcal{E}\left( \gamma(\ell-x) \right)\big]$$

for straight line segment diffraction, where $\gamma =\sqrt {\frac {\pi }{\lambda \zeta }}.$

3. Arc apertures

We convert expression (2) to polar coordinates $(x,y)=(r\cos {\theta }, r\sin {\theta })$:

(10)$$\begin{aligned} A(r,\theta) &= a\rho\int_0^{\alpha} \exp\Big(\frac{\pi i}{\lambda \zeta}\big[(r\cos{\theta}-\rho\cos{\phi})^{2} + (r\sin{\theta}-\rho\sin{\phi})^{2} \big]\Big) d\phi \\ &= a\rho\int_0^{\alpha} \exp\Big(\frac{\pi i}{\lambda \zeta}\big[ r^{2} + \rho^{2} -2r\rho\cdot(\cos{\phi}\cos{\theta} + \sin{\phi}\sin{\theta}) \big]\Big) d\phi \\ &= a\rho\exp\Big(\frac{\pi i}{\lambda \zeta}(r^{2} + \rho^{2})\Big)\int_0^{\alpha} \exp\Big(-i\frac{2\pi r\rho}{\lambda \zeta}\big(\cos{\phi}\cos{\theta} + \sin{\phi}\sin{\theta}\big)\Big) d\phi \end{aligned}$$

which can be simplified further using the identity $\cos {\phi }\cos {\theta } + \sin {\phi }\sin {\theta } = \cos {(\phi -\theta )}$:

(11)$$A(r,\theta) = a\rho\exp\Big(\frac{\pi i}{\lambda \zeta}(r^{2} + \rho^{2})\Big)\int_{\theta}^{\theta-\alpha} \exp\Big(-i\frac{2\pi r\rho}{\lambda \zeta}\cos{\phi}\Big) d\phi.$$

The main remaining hurdle is to solve an integral of the form

(12)$$E(s,w) = \int_0^{w} \exp(is\cos{t})dt$$

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]

(13)$$ J_\nu(s) = \frac{s^{\nu}}{2^{\nu-1}\Gamma(\nu+\tfrac{1}{2})\sqrt{\pi}}\int_0^{\frac{\pi}{2}}{\cos{(s\cos{t})}\sin^{2\nu}{(t)}dt} $$

(14)$$H_\nu(s) = \frac{s^{\nu}}{2^{\nu-1}\Gamma(\nu+\tfrac{1}{2})\sqrt{\pi}}\int_0^{\frac{\pi}{2}}{\sin{(s\cos{t})}\sin^{2\nu}{(t)}dt} $$

where $J_\nu (s)$ is the Bessel function of the first kind of order $\nu$, $H_\nu (s)$ is the Struve function of order $\nu$ and $\Gamma$ is the gamma function. Using Euler’s formula $\exp (is)=\cos (s)+i\sin (s)$ and taking the order $\nu =0$, we obtain for the fixed integration bounds $w=\frac {\pi }{2},\pi$

(15)$$E(s,\pi) = \pi\cdot J_0(s); \qquad E\left(s,\frac{\pi}{2}\right) = \frac{\pi}{2}\big(J_0(s) + i H_0(s)\big)$$

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 $\forall s,w \in \mathbb {R}$, we can reduce the space of values under consideration by leveraging the following symmetries:

(16)$$E(-s,w) = \overline{E(s,w)}; \qquad E(s,w+\pi) = \frac{\pi}{2}J_0(s) + \overline{E(s,w)}; \qquad E(s,-w) = -E(s,w).$$

where the overbar stands for complex conjugation. Thus we will only consider $s>0$ and $0\le w < \tfrac {\pi }{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\sqrt {s}$, which will be expounded on in the subsections below.

3.1 Asymptotic approximation for large parameters

Substituting $u=\cos {w}$ in $E(s,w)$ and using (15), we obtain

(17)$$E(s,w) = \int_{\cos{w}}^{1} \frac{\exp(isu)}{\sqrt{1-u^{2}}}du = \frac{\pi}{2}\big(J_0(s) + i H_0(s)\big) - \int_0^{\cos{w}} \frac{\exp(isu)}{\sqrt{1-u^{2}}}du.$$

The latter integral can be calculated [25] by using the Fourier series expansion of $(1-u^{2})^{-\frac {1}{2}}$:

(18)$$\begin{aligned} & \int_0^{\cos{w}} \frac{\exp(isu)}{\sqrt{1-u^2}}du = \int_0^{\cos{w}}\exp(isu)\left(\frac{\pi}{2} + \pi\sum_{m=1}^\infty{J_0(m\pi)\cos{(m\pi u)}}\right)du = \\ & \frac{i\pi}{2s}(1 - e^{i\cos{w}}) + \sum_{m=1}^\infty{\frac{\pi J_0(m\pi)}{m^2\pi^2-s^2}\Big( e^{i\cos{w}}\big[is\cos{(m\pi\cos{w})+m\pi\sin{(m\pi\cos{w})}} \big] - is\Big)}. \end{aligned}$$

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:

(19)$$\frac{\partial}{\partial u}\frac{-i\exp(isu)}{s\sqrt{1-u^2}} = \frac{\exp(isu)}{\sqrt{1-u^2}} - \frac{iu\exp(isu)}{s(1-u^2)^\frac{3}{2}}$$

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

(20)$$\int_0^{\cos{w}} \frac{\exp(isu)}{\sqrt{1-u^{2}}}du \approx \frac{-i\exp(is\cos{w})}{s\sin{w}} + \frac{i}{s}.$$

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:

(21)$$H_\nu(s) - Y_\nu(s) = \frac{\left(\frac{s}{2}\right)^{\nu-1}}{\sqrt{\pi} \Gamma \left (\nu+\frac{1}{2} \right )} + \mathcal{O}\left(\left (\tfrac{s}{2}\right)^{\nu-3}\right)$$

where $Y_\nu (s)$ is the Neumann function of order $\nu$, also known as a Bessel function of the second kind, which is easier to calculate. Substituting this for $\nu =0$ gives us:

(22)$$H_0(s) = Y_0(s) + \frac{2}{\pi s} + \mathcal{O}\left(\frac{s^{-3}}{8}\right)$$

this first error term will exactly compensate the $\frac {i}{s}$ term derived in (20), resulting in the final expression

(23)$$E(s,w) \approx E_L(s,w) = \frac{\pi}{2}\mathcal{H}_0(s) - \frac{i\exp(is\cos{w})}{s\sin{w}}$$

valid for larger values of $s$ and $\sin {w}$, where $\mathcal {H}_0(s) = J_0(s) + i Y_0(s)$ is the Hankel function of order 0.

3.2 Series approximation for small parameters

For small values of $w\sqrt {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)$,

(24)$$\exp(is\cos{t}) = 1 + is\cos{t} - \tfrac{1}{2}s^{2}\cos{t}^{2} - \tfrac{1}{6}is^{3}\cos{t}^{3} + \mathcal{O}(s^{4})$$

which can be integrated term by term to obtain

(25)$$\int_0^{w} \exp(is\cos{t})dt = w + is\sin{w} - \tfrac{1}{4}s^{2}(w+\cos{w}\sin{w}) - \tfrac{i}{18}s^{3}\sin{w}(2+\cos{w}^{2}) + \mathcal{O}(s^{4}).$$

Unfortunately, this expression is only useful for very small values of $\lvert s \rvert < 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:

(26)$$\exp(is\cos{t}) = \exp\Big( is\left(1 - \tfrac{1}{2}t^{2} + \mathcal{O}(t^{4})\right)\Big)$$

valid for small values of $w\approx \sin {w}$, after integration

(27)$$\int_0^{w} \exp(is\cos{t})dt \approx \int_0^{w} \exp\Big(is\left(1-\tfrac{1}{2}t^{2}\right)\Big)dt = \frac{\pi}{2} \exp{\left(i\left(s-\tfrac{\pi}{4}\right)\right)}\operatorname{erf}{\left(\tfrac{1+i}{2}w\sqrt{s}\right)}.$$

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\sqrt {s}$:

(28)$$E(s,w) \approx E_S(s,w) = \frac{1}{\sqrt{2s}} \exp{(is)}\mathcal{E}{\left(\frac{w\sqrt{2s}}{2}\right)}.$$

3.3 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:

(29)$$\operatorname{RP}{(a,\hat{a})} = -\log_{10}{\left\lvert\frac{a - \hat{a}}{a}\right\rvert} = \log_{10}{\left\lvert a\right\rvert} - \log_{10}{\left\lvert a - \hat{a}\right\rvert}$$

where $a$ is the reference value and $\hat {a}$ 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 $\operatorname {RP}{(E, E_L)} > \operatorname {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\sqrt {s} = \mathit {c}$ (where $\mathit {c}$ is the parameter to be fitted), we get that

(30)$$\big[\operatorname{RP}{(E(s,w), E_L(s,w))} = \operatorname{RP}{(E(s,w), E_S(s,w))}\big] \approx \left(w\sqrt{s} = 2.26\right)$$

with a goodness-of-fit $R^{2} = 0.9985$ when tested in the region $s \in \, ]0, 10^{5}].$

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\sqrt{s} = 2.26)$, shown by the red dashed line.

Download Full Size | PDF

Finally, we have to consider the case $s\approx 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

(31)$$E(s,w) \approx \begin{cases} w + is\sin(w), & \textrm{if}\ s < 10^{-2} \\ E_S(s,w), & \textrm{else if}\ w\sqrt{s} < 2.26 \\ E_L(s,w), & \textrm{otherwise}. \end{cases}$$

3.4 Evaluating arc aperture expressions

With the help of the incomplete cylindrical function $E(s,w)$, we can rewrite (11) as

(32)$$A(r,\theta) = a\rho\exp\left(i\gamma^{2}(r^{2} + \rho^{2})\right) \big[ E(s,\alpha + \theta) - E(s,\theta) \big]$$

where $s = -2\gamma ^{2} r\rho$ and $\gamma =\sqrt {\frac {\pi }{\lambda \zeta }}$ as defined for (9). However, we approximated $E(s,w)$ only for $s\ge 0$ and $0\le w \le \tfrac {\pi }{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 $\gamma ^{2}$, whose sign is equal to that of the depth position $\zeta$. Typically, the sign of $\zeta$ for all elements in a scene will be the same, so we can fill in $\lvert s \rvert$ 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 $\tfrac {\pi }{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:

(33)$$E(s,w) = m\pi J_0(s) + \operatorname{bitsgn}{(n)}\Re{\{E(s,\lvert w - m\pi \rvert)\}} + i\operatorname{bitsgn}{\left(\lfloor\tfrac{n}{2}\rfloor\right)}\Im{\{E(s,\lvert w - m\pi \rvert)\}}$$

where $\lfloor \cdot \rfloor$ is the floor operator (rounding down), with

(34)$$n = \left\lfloor\frac{2w}{\pi}\right\rfloor, \qquad m = \left\lfloor\frac{n+1}{2}\right\rfloor, \qquad \operatorname{bitsgn}{(t)} = \begin{cases} +1, & \textrm{if}\:t \equiv 0\:(\textbf{mod}\:2) \\ -1, & \textrm{otherwise}. \end{cases}$$

where $\textbf {mod}$ stands for the modulo operation. Thus, “$\operatorname {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,\lvert w - m\pi \rvert )$ 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 $\pi J_0(s)$ only has to be computed once per evaluated value of $A(r,\theta )$ and can be re-used for both evaluations of $E(s,w)$ as well as in the approximated case $E_L(s,w).$

4. 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.

4.1 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.

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

Download Full Size | PDF

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

Download Full Size | PDF

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 convolution for computing the Fresnel diffraction. This requires $(D+1)$ FFT calculations: $D$ for each plane and a final extra (I)FFT to obtain the hologram after adding all the Fourier transforms together.
• 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”.

Table 1. Run times and numerical accuracy of the different algorithms for both tested data sets.

View Table | View all tables in this article

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 $\hat {H}$ differs from the reference. The general definition of the PSNR for a real- / complex-valued signal $H$ is

(35)$$\mathrm{PSNR}(H, \hat{H}) := 10\log_{10}{\left(\frac{\max{(\lvert H \rvert)}}{\left\lVert{H-\hat{H}}\right\rVert_2}\right)}$$

where $\left \lVert {\cdot }\right \rVert _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.

Table 2. Error measurements for the CGH rendered images from the “Simple shapes” data set, for different algorithms w.r.t. the reference Point-based CGH algorithm, numerically backpropagated at different depths.

View Table | View all tables in this article

Table 3. Error measurements for the CGH rendered images from the “Seigaiha” data set, for different algorithms w.r.t. the reference Point-based CGH algorithm, numerically backpropagated at different depths.

View Table | View all tables in this article

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.

Download Full Size | PDF

The proposed 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)

(36)$$\mathrm{NMSE}(H, \hat{H}) := \sqrt{\frac{\left\lVert{H-\hat{H}}\right\rVert_2}{\left\lVert{H}\right\rVert_2}}$$

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.

4.2 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.

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

Download Full Size | PDF

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.

View Table | View all tables in this article

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 at 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.

Table 5. Reconstructions of the “Seigaiha” data set CGH, refocused at 0.50m.

View Table | View all tables in this article

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.

5. 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.

Funding

Fonds Wetenschappelijk Onderzoek (12ZQ220N, VS07820N); Kenjiro Takayanagi Foundation; Inoue Foundation for Science; Japan Society for the Promotion of Science (20K19810).

Disclosures

The authors declare no conflicts of interest.

References

1. F. Yaras, H. Kang, and L. Onural, “State of the art in holographic displays: A survey,” J. Disp. Technol. 6(10), 443–454 (2010). [CrossRef]

2. T. Dresel, M. Beyerlein, and J. Schwider, “Design of computer-generated beam-shaping holograms by iterative finite-element mesh adaption,” Appl. Opt. 35(35), 6865–6874 (1996). [CrossRef]

3. J. H. Burge, “Applications of computer-generated holograms for interferometric measurement of large aspheric optics,” in International Conference on Optical Fabrication and Testing, vol. 2576T. Kasai, ed., International Society for Optics and Photonics (SPIE, 1995), pp. 258–269.

4. 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 Commun. 70, 114–130 (2019). [CrossRef]

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

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

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

8. 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(16), 15205–15211 (2011). [CrossRef]

9. 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(1), 112–123 (2017). [CrossRef]

10. 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(26), 34259–34265 (2018). [CrossRef]

11. D. Blinder, “Direct calculation of computer-generated holograms in sparse bases,” Opt. Express 27(16), 23124–23137 (2019). [CrossRef]

12. Y. Zhao, L. Cao, H. Zhang, D. Kong, and G. Jin, “Accurate calculation of computer-generated holograms using angular-spectrum layer-oriented method,” Opt. Express 23(20), 25440–25449 (2015). [CrossRef]

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

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

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

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

17. H. Nishi and K. Matsushima, “Rendering of specular curved objects in polygon-based computer holography,” Appl. Opt. 56(13), F37–F44 (2017). [CrossRef]

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

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

20. 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(4), 3901–3913 (2015). [CrossRef]

21. T. Nishitsuji, T. Shimobaba, T. Kakue, and T. Ito, “Fast calculation of computer-generated hologram of line-drawn objects without fft,” Opt. Express 28(11), 15907–15924 (2020). [CrossRef]

22. D. Meek and D. Walton, “Approximating smooth planar curves by arc splines,” J. Comput. Appl. Math. 59(2), 221–231 (1995). [CrossRef]

23. G. N. Watson, A treatise on the theory of Bessel functions (MacMillan, 1945).

24. M. M. Agrest, M. S. Maksimov, H. E. Fettis, J. Goresh, and D. Lee, Theory of incomplete cylindrical functions and their applications, vol. 160 (Springer, 1971).

25. H. Nagaoka, “Diffraction phenomena produced by an aperture on a curved surface,” The journal Coll. Sci. Imp. Univ. Jpn. pp. 301–322 (1891).

26. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. on Image Process. 13(4), 600–612 (2004). [CrossRef]

	Simple shapes		Seigaiha
Algorithm version	time (ms)	PSNR (dB)	time (ms)	PSNR (dB)
Point-based (GPU)	157	(reference)	484	(reference)
FFT-based (GPU)	277	38.1	285	35.8
CG method (CPU)	106	39.4	159	36.0
Analytical method (CPU)	1144	73.7	14048	51.5
Analytical method (GPU)	2	73.7	21	51.5

	Analytical		FFT-based		CG-based
depth (m)	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
0.24	23.31	0.61	22.07	0.38	21.81	0.46
0.28	24.22	0.65	23.25	0.41	22.22	0.47
0.32	25.24	0.68	24.19	0.45	23.71	0.52
0.36	24.36	0.69	23.21	0.42	22.49	0.49
0.40	25.21	0.73	23.77	0.44	23.38	0.51
0.44	26.99	0.78	25.15	0.50	24.15	0.54
0.48	28.00	0.80	25.44	0.55	22.49	0.47
0.52	25.97	0.76	23.09	0.46	23.31	0.49
0.56	27.13	0.77	23.01	0.43	23.62	0.47
0.60	29.70	0.81	25.39	0.53	21.93	0.37
0.64	30.70	0.80	25.02	0.50	25.34	0.47
0.68	28.17	0.74	23.47	0.41	23.79	0.42
average	26.58	0.74	23.92	0.46	23.19	0.47

	Analytical		FFT-based		CG-based
depth (m)	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
0.10	28.22	0.86	19.16	0.20	18.09	0.11
0.20	30.40	0.90	20.42	0.27	18.30	0.15
0.30	34.97	0.93	21.38	0.31	18.53	0.22
0.40	36.72	0.94	21.43	0.31	18.47	0.20
0.50	37.13	0.94	21.39	0.30	18.36	0.21
0.60	36.96	0.94	21.43	0.29	18.08	0.20
0.70	36.66	0.93	20.84	0.26	18.14	0.18
0.80	34.30	0.90	20.35	0.23	18.49	0.17
average	34.42	0.92	20.80	0.27	18.31	0.18

	Simple shapes		Seigaiha
Algorithm version	time (ms)	PSNR (dB)	time (ms)	PSNR (dB)
Point-based (GPU)	157	(reference)	484	(reference)
FFT-based (GPU)	277	38.1	285	35.8
CG method (CPU)	106	39.4	159	36.0
Analytical method (CPU)	1144	73.7	14048	51.5
Analytical method (GPU)	2	73.7	21	51.5

	Analytical		FFT-based		CG-based
depth (m)	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
0.24	23.31	0.61	22.07	0.38	21.81	0.46
0.28	24.22	0.65	23.25	0.41	22.22	0.47
0.32	25.24	0.68	24.19	0.45	23.71	0.52
0.36	24.36	0.69	23.21	0.42	22.49	0.49
0.40	25.21	0.73	23.77	0.44	23.38	0.51
0.44	26.99	0.78	25.15	0.50	24.15	0.54
0.48	28.00	0.80	25.44	0.55	22.49	0.47
0.52	25.97	0.76	23.09	0.46	23.31	0.49
0.56	27.13	0.77	23.01	0.43	23.62	0.47
0.60	29.70	0.81	25.39	0.53	21.93	0.37
0.64	30.70	0.80	25.02	0.50	25.34	0.47
0.68	28.17	0.74	23.47	0.41	23.79	0.42
average	26.58	0.74	23.92	0.46	23.19	0.47

Analytic computation of line-drawn objects in computer generated holography

Abstract

1. Introduction

2. Line apertures

2.1 Straight line apertures

3. Arc apertures

3.1 Asymptotic approximation for large parameters

3.2 Series approximation for small parameters

3.3 Error analysis

3.4 Evaluating arc aperture expressions

4. Experiments

4.1 Numerical calculation experiments

4.2 Optical experiments

5. Conclusion

Funding

Disclosures

References

Cited By

Figures (6)

Tables (5)

Equations (36)

Optics Express