1. INTRODUCTION

Computer-generated holography (CGH) can create custom 3D intensity patterns by modulating the wavefront of a coherent light beam. This ability to control light volumetrically is critical to applications such as 3D displays [1], optical trapping [2], and neural photostimulation [3–6]. Typically, 3D CGH is implemented by encoding the phase of a coherent wavefront in the Fourier space (pupil plane) of an imaging system with a spatial light modulator (SLM). The SLM modulates the 2D Fourier phase of the wavefront, allowing one to control the 3D intensity distribution near the image plane.

Computing the 2D Fourier phase mask that generates a specific 3D intensity pattern is usually a non-trivial and ill-posed problem, for several reasons. First, most 3D intensity patterns are not physically feasible. The wave-field must propagate through free space while satisfying Maxwell’s equations. Intensity must be conserved at each depth plane, and resolution is limited to the coherent diffraction limit set by the numerical aperture (NA). As a result, arbitrary 3D intensity patterns cannot be created exactly. Another limitation is the finite degrees-of-freedom (DOF) for patterning. Holograms must be phase-only, and the SLM dynamic range and pixel count impose further limits. Consider, for example, the DOF mismatch of our setup: we have $\sim {10}^9$ resolvable voxels, but only $\sim {10}^6$ SLM pixels. Given all these considerations, 3D CGH algorithms must search for the best approximate solution to the problem, within the constraints imposed by physics and hardware. The resulting optimization task is non-convex with non-unique solutions in a high-dimensional space (millions of pixels) that cannot be explored exhaustively.

Most existing algorithms can be classified into two categories. The first, superposition algorithms [7–9], decompose the desired 3D intensity into discrete 2D planes and digitally back-propagate each to the SLM plane. Each plane is considered independently and the solution is the average of the phases of the resulting complex fields. This class of algorithms is computationally efficient and always yields a solution. However, interference between planes is not accounted for, so the results are sub-optimal, and quality degrades quickly with more than 2–3 depth levels.

The second category of 3D CGH algorithms is iterative projection algorithms [10]. These methods iteratively update the solution by projecting it onto a subspace computed from the set of constraints. Gerchberg and Saxton showed that alternating projections onto the intensity constraints at the hologram plane and one focal plane will not diverge (though it is not guaranteed to converge) [11]. Similar concepts apply in 3D [12,13], and recent variations use multiple constraints in each iteration [2,14,15]. This class of algorithms is easy to implement and gives good results for single-plane intensity patterns. However, the strong enforcement of constraints may cause the algorithm to diverge when attempting to generate physically infeasible patterns. Introducing hyper-parameters can improve results [16], but it is difficult to tune without an explicitly defined cost function.

Here, we present a new algorithm, termed Non-Convex Optimization for VOlumetric CGH (NOVO-CGH) that combines a customizable cost function with an optimization strategy inspired by prior work in phase retrieval [17]. Phase retrieval from through-focus intensity measurements represents a similar mathematical problem to 3D CGH, since both take a 3D intensity distribution as input and solve for the 2D phase that generates it. However, there is a major difference in that the 3D intensity patterns in the phase retrieval case are always physically feasible, whereas they are generally not in the 3D CGH case.

For 3D CGH, where there is no exact solution and many approximate solutions, results depend heavily on the cost function used. Our NOVO-CGH algorithm enables a search for optimal holograms with custom cost functions that can be chosen by the user to best respond to a specific application. For example, by specifying dark and bright regions with intensity thresholds, as well as regions where target intensity is unspecified, we give the algorithm flexibility to find a suitable solution. While our algorithm cannot expand the set of physically feasible intensity patterns, we will demonstrate that, within this same set, it can find better solutions than other algorithms.

In Section 2, we mathematically describe the 3D CGH problem, review previous methods, and introduce our algorithm. In Section 3, we discuss the specific case of binary 3D intensity patterns for which we propose a customized multi-part cost function. In Section 4, we demonstrate improved performance over existing methods, both in simulation and experiment. Finally, we show a specific application in neural photostimulation, achieving enhanced spatial confinement of two-photon absorption when using NOVO-CGH. In Supplement 1, we provide additional details and implementations of NOVO-CGH for applications with non-binary target holograms.

2. THEORY

3D CGH algorithms take the target 3D intensity pattern, $V(x,y,z)$, as an input and return a 2D phase, $\varphi$, to be displayed on the SLM. Given the general setup in Fig. 1(a), a coherent plane wave, $A_0$, is incident on a phase-only SLM, which patterns the wavefront to generate the complex field $P_1(x,y)=A_0\exp (i\varphi (x,y))$. ${\mathcal{T}}_f$ denotes the transfer function from Fourier space at $P_1$ to image space at $P_2$, and ${\mathcal{T}}_z$ denotes propagation between $P_2$ and $P_3$ for the simple case of a two-plane target intensity pattern. Extension to more depth planes is straightforward.

 

Fig. 1. (a) Experimental setup for 3D CGH. A SLM imposes a phase pattern (hologram) on an incoming collimated plane wave. A Lens (focal length, $f$) in $2f$ configuration transforms the resulting wavefront into a 3D intensity distribution in the volume-of-interest. 3D CGH algorithms aim to generate a target 3D intensity distribution, $V(x,y,z)=V_z(x,y)$, within the volume-of-interest by designing, $\varphi$, the phase pattern on the SLM. (b) For our binary NOVO-CGH algorithm, we introduce a custom cost function with two parts, $L_1$ and $L_2$, which are truncated square functions defined in Eq. (9). $L_1$ penalizes intensities above a certain threshold, $l$, for dark voxels of the target pattern, and $L_2$ is the opposite (penalizes intensity below a threshold, $h$, for bright voxels). Total cost is the sum of all pixels’ cost values.

Download Full Size | PDF

Under the paraxial approximation, the complex fields can be written as $P_2={\mathcal{T}}_f{P}_1$ and $P_3={\mathcal{T}}_z{P}_2$:

A. Existing Algorithms

1. Superposition Algorithm

Superposition algorithms aim to simplify the 3D CGH problem by splitting the desired intensity into segments and solving each independently. Our implementation uses the random superposition (RS) algorithm [2], except that instead of decomposing the intensity pattern into points, we decompose it into planes. This accelerates the computation and reduces errors due to intra-plane interference. As with other algorithms, the number of planes along the $z$-axis can be adjusted, and uniform sampling is not required. The phase hologram is computed as

The simplicity of the superposition algorithm is also what limits its accuracy. Since amplitude information is discarded in the last step, it may not provide the optimal phase-only pattern. More importantly, the algorithm achieves 3D multiplexing by averaging the holograms corresponding to each plane. This assumes that each plane is independent, ignoring wave interference effects between planes. This approximation rapidly becomes invalid for dense 3D patterns with multiple depth planes.

2. Iterative Projection Algorithms

We also compare our algorithm with two popular iterative projection algorithms. Both are based on the Gerchberg–Saxton (GS) algorithm: sequential GS [12] and global GS [14] (see Supplement 1), and both use the superposition algorithm for initialization.

The sequential GS algorithm (also known as the ping-pong algorithm) digitally propagates the estimated complex field through each depth plane, one-by-one, while replacing the amplitude of the estimate with the desired amplitude at that plane. Then, the complex field is propagated to the SLM, where only phase is retained. This hologram becomes the estimate for starting the next iteration. This procedure is equivalent to stochastic gradient descent (SGD) with a fixed step size [19]. It diverges easily for infeasible target intensities, so we make a small modification that relaxes constraints slightly to enable convergence in most cases [16]; instead of forcing the dark areas to be zero, we let ${\mathcal{P}}_{V,z}E(x,y)=\mu \exp (j\arg [E(x,y)])$, where $\mu$ is a small number.

The global GS algorithm corresponds to batched SGD. Instead of sequentially propagating through all depth planes, the current phase hologram estimate is propagated to all depth planes at once, and the holograms from each are averaged at each iteration.

Though sequential and global GS algorithms are similar, their performance is very different. Sequential GS tends to diverge when the target 3D intensity is not physically feasible or becomes dense, while global GS is more stable because it simultaneously averages the contributions of all planes at each iteration.

B. Non-Convex Optimization for Volumetric CGH (NOVO-CGH)

Here, we propose a new algorithm, NOVO-CGH, that gives improved performance when it is not possible to exactly satisfy all the intensity constraints (the target pattern is physically infeasible). In this case, the choice of cost function and optimization procedure become critical. Previous algorithms used either fixed or implicit cost functions. Here, to enable more flexible addition of priors and customized cost functions, we start by formulating the 3D CGH problem as an optimization problem:

This problem is non-convex (the Hessian matrix of the cost function is not positive semi-definite), but when the forward model, $I(\varphi )$, and the cost function have a well-defined derivative, we can still solve the problem using first or second-order optimization methods. Here, we use L-BFGS [20], a quasi-Newton gradient descent procedure that always converges to a stationary point [21] and only requires limited amounts of computer memory. Our NOVO-CGH algorithm is summarized in Algorithm 1.

Algorithm 1. NOVO-CGH Algorithm

View Table

To compute the gradient, $\frac{\partial L}{\partial \varphi }$, directly, we consider the 2D complex field at $P_1$ given by $E=\exp (i\varphi )$ and derive

In the following notation, we convert $I$ and $V$ into 1D vectors of size $N\times 1$, where $N$ is the total number of voxels in $I$. $\frac{\partial L}{\partial \varphi }$ has size $1\times M$, and $\frac{\partial I}{\partial E}$ has size $N\times M$, where $M$ is the number of pixels on the SLM. For these vectors, we let ${(·)}^T$, ${(·)}^{*}$, and ${(·)}^H$ denote the transpose, conjugate, and Hermitian operators, respectively. To find $\frac{\partial I}{\partial E}$, we first define $I_z(x,y)=I(x,y,z)$, then compute the gradient for each value of $z$.

In discrete space, we can represent linear operations as matrix multiplication, giving ${\mathcal{T}}_f=\gamma F^H$ and ${\mathcal{T}}_z=F^H{K}_{z}F$, where $F$ is the discrete Fourier transfrom (DFT) matrix, $\gamma$ is a scaling coefficient determined by the number of pixels, and $K_z$ is a diagonal matrix representing Fresnel propagation. Hence, ${\mathcal{T}}_z{\mathcal{T}}_f={\gamma }^{\prime }{F}^H{K}_z$. Since intensity is point-wise amplitude squared, that is, $I_z={\Vert {\mathcal{T}}_z{\mathcal{T}}_{f}E\Vert }^2=\mathrm{diag}{[{\gamma }^{\prime }{F}^H{K}_{z}E]}^H{\gamma }^{\prime }{F}^H{K}_{z}E$, we can write

The only remaining requirement is that the cost function should be defined with respect to the intensity in the volume-of-interest and differentiable. NOVO-CGH can therefore be implemented with any cost function for which both quantities, $L(I(\varphi ),V)$ and $\frac{\partial L}{\partial I}$, can be expressed within the discretized volume. By combining $\frac{\partial L}{\partial I}$ with Eqs. (5)–(3), we deduce the gradient which is required for $L$-BFGS optimization.

In Sections 3 and 4, we propose cost functions for several specific applications. In Supplement 1, we consider alternate implementations of NOVO-CGH for 3D holography with variable-intensity patterns and multi-photon holography. In each case, we derive the cost function and corresponding gradient.

3. BINARY 3D HOLOGRAPHY

We first consider the case of binary intensity targets, which is relevant to applications in optical trapping and neural photostimulation [22,23], where precise control over light-sculpting capabilities is critical. Neurons can be modified optogenetically to be made photosensitive, which in theory enables the triggering [24] or inhibition [25] of neural activity on demand with light. In practice, a neuron’s response to light is highly nonlinear, and activating custom neuron ensembles requires the ability to concentrate light on targeted neurons while keeping the intensity received by non-targeted neurons below some threshold value to prevent unwanted activation. With NOVO-CGH, custom cost functions can be tailored to match the nonlinear physiological response of the neurons.

Although the choice of cost function is flexible, we propose the following cost function:

Calculating the derivative of this cost function gives

Since all matrix products in the derivative are either Fourier transforms, inverse Fourier transforms, or element-wise multiplications, the memory requirement is only $O(N)$, and computing the gradient remains spatially and temporally efficient.

4. RESULTS

We demonstrate our algorithm’s improved performance over previous techniques, in both simulation and experiment, for two applications: multi-plane displays and neural photostimulation.

In simulations, we have ground truth information and use two metrics to numerically evaluate performance—accuracy, $\alpha$, and efficiency, $\eta$. Accuracy measures the resemblance of the generated intensity pattern to the target pattern, defined as

A. Multi-Plane Display

To compare our method to the superposition and iterative projection algorithms, we use each method to generate phase holograms of the same 3D target intensity patterns. We first simulate a relatively dense pattern that will not have a feasible solution [Fig. 2(a)]. The target intensity pattern contains 10 depth planes with different letters at evenly spaced $z$ distances. Each letter extends an area of $1\mathrm{mm}\times 1\mathrm{mm}$, and the axial range spans 45 mm. The system properties for simulations match our experimental setup (Fig. S1a of Supplement 1). The system NA is limited by the size of the SLM ($\mathrm{NA}\approx 12.5\mathrm{\mu m}*1280/(2\times 200\mathrm{mm})=0.04$).

 

Fig. 2. Simulation results comparing NOVO-CGH to existing algorithms. (a) Target 3D intensity pattern containing 10 letters at evenly spaced depth planes (5 mm separation). (b) Quantitative comparison of simulation results between superposition [2], sequential Gerchberg–Saxton (GS) [12], global GS [14], and our binary NOVO-CGH algorithm. Normalized intensity histograms compare intensity values within bright and dark voxels separately, showing that NOVO-CGH yields more uniform intensity. (c) (Left) Computed 2D phase mask to be displayed on the SLM. (Center) Resulting intensity distributions at each patterned depth plane. (Right) Intensity slice in the $x\text{-}z$ plane. (d) Intermediate intensity at planes between the “A” and “B” planes for NOVO-CGH. All intensity images have the same colorbar. Propagation images through intermediate planes are shown in Visualization 1.

Download Full Size | PDF

As expected, the superposition algorithm generates a noisy result, because it ignores coherent interference from wave propagation between planes (see Fig. 2). Sequential GS also performs poorly, since it heavily favors the last plane used as a constraint. Global GS produces a better result, but it still has unwanted stray light in each plane. Our NOVO-CGH algorithm produces the most accurate result, particularly in the dark areas. Figure 2(b) plots the intensity distribution within bright and dark voxels, showing that our algorithm yields less variation. Note that all simulated intensity patterns have the same total energy, but some images appear darker because energy is spread or cropped by the colorbar. Additional slice views at the depth planes where the letters are [Fig. 2(c)] and at intermediate depths [Fig. 2(d)] show enhanced light confinement capabilities with NOVO-CGH. Media 1 shows the intensity at all planes through the 3D volume.

To quantify performance, Fig. 3 shows accuracy and efficiency plots for the simulated results. The three iterative algorithms generally improve with each iteration, eventually converging. Sequential GS exhibits a drop in performance during the first iteration because the algorithm is initialized with the superposition algorithm result, which is reasonably good but does not satisfy any single constraint. Therefore, projecting it onto any one of the constraints drags the solution further away from other constraints. This would not happen if the pattern was physically feasible. In Fig. 3(b), we show the accuracy and efficiency of each algorithm’s converged solution as a function of the number of intensity planes in $V$ (increasing complexity of the 3D pattern). As expected, when more planes are added, performance degrades, since the solution becomes less feasible. NOVO-CGH always produces the most accurate result, though the global GS algorithm has slightly higher efficiency.

 

Fig. 3. Numerical comparison of algorithm performance in simulation. (a) Efficiency and accuracy versus number of iterations for a target intensity with 3 depth planes. Superposition algorithm is not iterative, so is drawn as a dotted line. (b) Efficiency and accuracy degrade as the target intensity pattern becomes more complicated (contains more depth planes).

Download Full Size | PDF

To experimentally verify our algorithm, we conducted experiments (using the same parameters as simulations) over a volume of $1.5\mathrm{mm}\times 1.5\mathrm{mm}\times 35\mathrm{mm}$. First, we show a simple 3D target intensity pattern consisting of only two planes [Fig. 4(a)]. In this case, all algorithms perform well. However, when the target intensity contains more planes with different letters, the solution becomes less feasible and performance varies. In the case of six letter planes [Fig. 4(b)], NOVO-CGH achieves higher efficiency, less speckle, and better fidelity than the other algorithms.

 

Fig. 4. Experimental results demonstrating improved 3D intensity patterning with NOVO-CGH, as compared to previous algorithms. (a) Comparison of simulation and experimental results for a two-plane intensity target. Calculated phase holograms are shown at left. (b) Experimental measurements for a six-plane intensity pattern. All images use the same camera settings.

Download Full Size | PDF

B. Two-Photon NOVO-CGH for Neural Photostimulation

1. Enhanced Confinement Around a Single Target

Holographic two-photon photostimulation [26,27] is an immediate application that directly benefits from custom cost functions. To illustrate the benefits of NOVO-CGH, we first consider the simple case where we wish to illuminate the entire soma of a specific neuron, while avoiding surrounding neurons. The target intensity pattern is shown in Fig. 5(a). One plane with a disc pattern (approximately the size and shape of a neuronal soma in the mouse cerebral cortex) is surrounded by dark planes. By specifying a protected volume around the target, our algorithm seeks a solution with good axial confinement for use with dense neural networks. For this task, the sequential GS algorithm fails to converge with the added dark planes, so we use only the center plane with the disc as input. The global GS and NOVO-CGH algorithms will have their computation scale linearly with number of dark planes added in the target intensity pattern.

 

Fig. 5. Experimental results for neural photostimulation. (a) Target 3D intensity is a bright disc (10 μm radius) at one depth, surrounded by darkness. Planes are evenly spaced by 0.03 mm. (b) $x\text{-}y$ and $y\text{-}z$ cross-sections of two-photon absorption (fluorescence) measured along the center of the volume. (c) Two-photon absorption axial cross-cuts. (d) Sorting all voxels in descending order, we plot two-photon absorption versus volume. Given a neuron of volume $8\times {10}^3{\mathrm{\mu m}}^3$, values in the pink region represent intensity in undesired regions.

Download Full Size | PDF

Figure 5 shows that NOVO-CGH provides improved axial confinement in two-photon experiments. We use a femtosecond infrared laser source (Fidelity 18 W, $\lambda =1040\mathrm{nm}$) and a microscope objective ($20\times$, 1.0 NA water immersion) to de-magnify the volume of operation. The target pattern [Fig. 5(a)] is a disc of radius 10 μm at a single depth plane (dark everywhere else). For characterization, we measure the two-photon absorption in 3D by observing fluorescence induced in a thin film at various depths (see Supplement 1).

Figures 5(b) and 5(c) shows cross-sections of the measured 3D two-photon absorption. All algorithms successfully generate a lateral disc pattern, but NOVO-CGH provides the best axial confinement. To explore the implications of our algorithm’s improvement on the relative threshold values for two-photon photostimulation, Fig. 5(d) shows normalized two-photon absorption as a function of an artificially increasing volume threshold. In any given volume size, the two-photon absorption threshold is consistently lower than NOVO-CGH compared to other algorithms, indicating better volumetric confinement of light.

2. Multiple Targets with Unspecified Regions

A major advantage of NOVO-CGH for 3D light sculpting in optogenetic photostimulation and photolithography is the ability to define spatially dependent cost functions that enforce both bright and dark voxels, while leaving other areas unspecified. Dark voxels enforce a penalty in locations where illumination should be avoided, instead of enforcing it uniformly in all non-bright voxels (see Section 3 of Supplement 1). This strategy relaxes constraints and searches for holograms for which there exists a pathway for photons to reach the bright targets while avoiding the dark targets, better protecting sensitive parts of the sample from receiving inadvertent illumination.

To demonstrate in simulation, we consider a 3D distribution of sphere-like targets within the accessible volume [Fig. 6(a)]. There are two categories of targets: red regions denote bright targets, where we wish to deposit energy, and blue regions denote dark targets, which should be protected from unwanted illumination. The remainder of the volume is unspecified (may have any intensity value). With other algorithms [Figs. 6(b)–6(d)] and binary NOVO-CGH [Fig. 6(e)], the user would only specify $V(x,y,z)=1$ within the red targets, then attempt to enforce darkness everywhere else within the volume. Here, we selectively relax the darkness constraint by applying a strong penalty on dark target voxels (see Supplement 1). We specify $V^{+}(x,y,z)=1$ within the bright targets, $V^{+}=0$ elsewhere, and $V^{-}(x,y,z)=10$ within the dark targets (see Fig. S2) and 0 elsewhere. Outside the bright and dark targets, $V^{+}=V^{-}=0$. The resulting 3D intensity generated by this variant of NOVO-CGH [Fig. 6(f)] finds a suitable optical pathway to illuminate the bright targets while avoiding the dark targets. For all targets within the volume, we compute the average received power and repeat the simulation 10 times with randomized target distributions (no overlap). Statistical data is shown in Fig. 6(g), confirming that dark targets are better protected from unwanted illumination with NOVO-CGH using unspecified regions.

 

Fig. 6. NOVO-CGH with unspecified regions. (a) Simulated 3D target intensity with 70 randomly distributed spheres. Red targets are desired bright regions, and blue targets are desired dark regions, with the remainder of the volume unspecified. (b)–(f) Resulting 3D intensity distributions from (b) superposition algorithm, (c) sequential GS, (d) global GS, (e) binary NOVO-CGH, and (f) NOVO-CGH with unspecified regions. We display in each case threshold envelopes of the volume intensity distribution (3% of maximal intensity, in red) alongside the blue regions. (g) For an equal amount of illumination within the red targets, we compute the amount of unwanted illumination in the blue targets, repeating the simulation 10 times with randomized distributions. NOVO-CGH out-performs other algorithms and further enhances confinement when unspecified regions are allowed.

Download Full Size | PDF

C. Two-Photon and Variable-Intensity NOVO-CGH

In Supplement 1, we provide additional variations for NOVO-CGH. One enables the rendering of grayscale (variable intensity) 3D intensity patterns. For image rendering applications, we propose a threshold-based and Euclidean cost function, then derive their respective gradients. Simulation results (see Fig. S3 of Supplement 1 and Visualization 2) show that our method enhances accuracy. Finally, we look at the case of two-photon imaging and propose a cost function (and its derivative) that directly optimizes for $I^2$ instead of $I$, tailoring holograms to nonlinear light-matter interactions (see Fig. S4 of Supplement 1).

D. Computation

The NOVO-CGH algorithm is implemented in Matlab with the $L$-BFGS method [20]. All algorithms used a graphics processing unit (GPU) for accelerated computation. Runtimes for generating holograms of size $1280\times 1024$ on our GPU (NVIDIA Tesla K20c) are shown in Fig. 7 as a function of the number of planes in the target volume intensity (which scales linearly with the number of voxels). The runtimes for superposition, GS, and NOVO-CGH algorithms all scale linearly with the number of planes. NOVO-CGH has the largest computational burden from computing the cost and the derivative, as well as a significant offset due to large overhead. Hence, the improvement in hologram quality comes at a cost of more computation time.

 

Fig. 7. Computation times on log scale for varying numbers of target depth planes. We use GPU computing to generate holograms of size $1272\times 1024$, with 50 iterations for the GS algorithms and the NOVO-CGH algorithm.

Download Full Size | PDF

5. CONCLUSION

We have demonstrated a new algorithm for improved 3D CGH by optimization with customizable cost functions that can be explicitly defined as a function of light intensity or space. The method can be tailored to particular applications or samples. For the example case of binary intensity targets, we demonstrated improved performance for multi-plane displays and enhanced confinement for light sculpting. We further included results for different cost functions to show that our algorithm can be easily adapted to different applications.

These results show that the customizable cost function of NOVO-CGH can be adapted not only to a specific type of holography but also to sample-specific features. The result is an algorithm that is better suited to accomplish specific tasks. For instance, in neuroscience applications, the cost function can be tuned to match the response of a specific opsin to light and can also take into account the spatial distribution of neurons within the 3D volume-of-interest to enable photostimulation of specific sets of neurons on demand. The same principle may also be applied to holographic photolithography [28] to match the optical response of a specific photoresist and to prevent a material from inadvertently curing in locations that would affect the shape of the structures being created.

We hope that NOVO-CGH can provide a different prospective for CGH algorithms since it is different from existing algorithms. Instead of paying attention to the solver, the user can now focus on designing the cost function to define what a good solution means. Such an approach addresses problems more directly, giving improved performance at the expense of higher computational cost. We want to emphasize that our algorithm does not expand the set of possible solutions, as that is constrained by the physics and hardware. Instead, NOVO-CGH finds a better solution within the set, based on particular tasks. Open-source code for NOVO-CGH can be found in our GitHub repository [29].