Proximal Algorithms for Discrete-Level Phase-Shifting Mask Design with Application to Optogenetics

This work studies the problem of designing computer-generated holograms using phaseshifting masks limited to represent only a small number of discrete phase levels. This problem has various applications, notably in the emerging field of optogenetics and lithography. A novel regularized cost function is proposed for the problem at hand that penalizes the unfeasible phase levels. Since the proposed cost function is non-smooth, we derive proper proximal gradient algorithms for its minimization. Simulation results, considering an optogenetics application, demonstrate that the proposed proximal gradient algorithm yields better performance as compared to other algorithms proposed in the literature.


Introduction
The research area of optogenetics studies the targeted stimulation of neurons that have been previously functionalized with light-sensitive compounds, called opsins [1]. In particular, genetically designed opsins are properly placed by physicians to some area of the brain and are able to produce electrical signals when hit by light of a specific wavelength. In the sequel, structured light is used to target specific areas of the brain and thus stimulate the neurons to which the opsins have been attached. Achieving precise spatial specificity is very important in order to minimize unnecessary heating of the tissue that should not be stimulated. Due to the prospects that it offers for specified neuronal stimulation, the field of optogenetics has been embraced by the research community as an efficient tool towards demystifying the operation of the brain [2]. Indeed, the optogenetic method has the potential of stimulating a complex of neurons down to a single neuron. This unprecedented specificity and spatial flexibility cannot be matched by other deep brain stimulation techniques.
Furthermore, in many studies, short-term excitation is required to study the dynamics of the neuronal activity, since the impulse response function of networks is best obtained with spatiotemporally well-defined stimuli. In practice, such short-term stimulation of neurons can be achieved only via the use of ferroelectric Spatial Light Modulators (fSLMs) that are characterized by very fast response times when a new phase-shifting mask is loaded to the device [3,4]. Thus, short-term stimulation is achieved via such fast swapping of different phase-shifting masks that have been designed beforehand for this task. However, typical fSLMs are restricted to using only a small number of discrete phase levels, most commonly two or four, where in the latter case, two binary fSLMs are employed [5].
To this end, it is becoming important to devise methods that design accurate computergenerated holograms, using discrete level phase-shifting masks (PSMs). Previous works on this problem have mainly focused on lithography applications. In particular, the branch and bound algorithm and simulated annealing were employed in [6] for both binary and phase-shifting mask design. Furthermore, an iterative approach to generate binary masks was reported in [7]. The method of projections onto convex sets (POCS) was employed in [8] for the PSM design problem. Furthermore, a genetic optimization algorithm was employed in [9] to jointly optimize for the mask and the light source parameters. In contrast to the previous methods that rely upon optimization approaches for discrete variables, the work in [10] introduced a method suitable for gradient-based optimization algorithms that, in general, are less computationally expensive. The approach in [11] further extended the work in [10] and derived a gradient-based algorithm that was able to generate more general mask patterns. More recently, the emerging field of algorithm unrolling [12] was employed to derive fast versions of the algorithm that was proposed in [11]. In particular, in [13,14], the iterative algorithm of [11] was unfolded to neural networks of suitable structure, which were then trained in an unsupervised manner.
While most previous works for the problem of discrete-level phase-shifting mask design have mainly considered lithography applications, in this work the focus is on the application of optogenetics. When considering an optogenetics application, the problem is similar to that encountered in lithography applications, in the sense that a phase-shifting mask with limited phase levels must be designed. However, it must be emphasized that in lithography applications there are a number of other restrictions regarding manufacturability [15] that are not present in our case.
In this work, we extend our previous work in [4] that considered the problem of discrete-level phase-shifting mask design in optogenetics applications. Similar to [11], we follow a gradient-based regularized optimization approach. In more detail, in our previous work [4], it was shown that performance gains can be obtained by properly modifying the regularization term proposed in [11]. However, these modifications make the regularization term non-differentiable and, thus, proper optimization algorithms must be employed. The subgradient [16] method was used in our previous work in [4] to deal with the non-smooth regularization term. The present is different from our previous work in [4] in two aspects: (i) we propose a new form for the regularization term that is easier to handle analytically, and (ii) we derive a proximal gradient [17] type algorithm for the minimization. Simulation results are provided that demonstrate the effectiveness of the proposed approach.
The rest of this manuscript is organized as follows. Section 2, contains our theoretical work. In more detail, Section 2.1 gives a mathematical formulation of the problem at hand and demonstrates why the globally optimal solution cannot be computed in practice. Section 2.2 provides the main ideas of our proposed approach, and Section 2.3 details the derivation of the proposed optimization algorithm. Numerical results that compare the performance of the proposed approach to other approaches are given in Section 3. Finally, Section 4 discusses the results in this work and outlines our future efforts.

Problem Formulation
Let us consider that the input phase mask is denoted as φ(x, y), where x, y ∈ {1, 2, . . . , N} denote the discrete spatial coordinates (pixels) on the SLM plane. As explained in the introduction, in this work we consider an application in which, due to technological limitations of the fSLM devices, φ(x, y) is not allowed to take any value in R, rather, only q equidistant phases are permitted, as denoted by where q denotes the number of discrete phase levels, typically 2 or 4 for the considered application, and Φ is the corresponding set of permissible phase values. We assume that monochromatic light from a proper laser source is reflected by the SLM device, thus generating the (input) light field where A ∈ R + denotes the intensity of the input light field and j denotes the complex imaginary unit. We also consider that I(x, y) is the input to an optical system that generates the output light field O(x, y) ∈ C, as described by the relation where the function H(·) describes the optical system under study. As an example, for lithography applications, H(·) is usually modeled as a low-pass Gaussian filter followed by an operation that computes the magnitude of the complex light field and the application of a sigmoid function [11]. In this work, for conducting our numerical simulations, we consider the optical path demonstrated in Figure 1, which consists of a 2f setup with focal length f followed by free-space propagation at a distance z. Since our focus is on an optogenetics application, the parameters of our simulation setup were selected so that the field-of-view at the target plane is an area of 266 × 266 µm. Clearly, the output light field O(x, y) is a function of the input light amplitude A and the phases φ(x, y). Given a desired response output intensity field D(x, y) ∈ R + , the optimal values A * and φ * (x, y) for the parameters A and φ(x, y) are given as the solution to the following optimization problem where d(·, ·) denotes some suitable distance/cost function, for example, the mean square error (MSE) or any other, more elaborate, cost function as in [18]. It is easily seen that problem (4) is a combinatorial optimization problem in the sense that the optimal solution can be obtained by examining all q N 2 possible phase-shifting masks and finding the one that yields the minimum cost, while the problem becomes even more complicated due to the continuous variable A. Thus, in practice, the optimal solution of (4) cannot be computed. A simple approach to circumvent this problem, and derive a sub-optimal solution, is to solve (4) by neglecting the constraint φ(x, y) ∈ Φ, letting φ(x, y) ∈ R, and then to replace each of the resulting phase values with the closest phase level from the constraint set Φ.
In the sequel, we refer to this approach as quantization.

The Proposed Approach
Trying to solve (4) while neglecting the constraint φ(x, y) ∈ Φ makes the problem significantly easier; however, the resulting solution (after quantization) exhibits poor performance, in the sense that the resulting output intensity |O(x, y)| 2 will be significantly different from the desired intensity D(x, y). To circumvent this effect, we employ the notion of regularized optimization. In essence, the method utilizes a modified cost function by adding a proper regularization part that penalizes values of the variables that do not satisfy the initial constraint. In more detail, we consider a cost function of the form where the first term is recognized as the mean squared error between the intensity |O(x, y)| 2 of the output light field O(x, y) and the desired intensity D(x, y), whereas the second term is the regularization term, that is, it obtains greater value for phase values away from those in the set Φ. The parameter λ, usually referred to as the regularizer, is some positive scalar that sets the relative importance between the mean squared error and the impact of the regularization. In our previous work [4], the regularization term was defined using where ρ > 0 is a positive scalar that controls the exact shape of the regularization terms in (6). It was experimentally shown in [4] that the value ρ = 0.5 gives good results. However, it can be seen that the regularization term in (6) is a non-smooth function of φ(x, y) when ρ < 1. Thus, the so-called subgradient method [16] was employed in [4] to minimize the resulting cost function that is suitable for such cases. The approach followed in this work is different from our previous work in [4] in two aspects. Firstly, in place of the regularization term in (6) we propose here the term where k(x, y) = φ(x,y) 2pi/q . Secondly, in this work, we derive a proximal gradient [17] type algorithm to deal with the fact that the considered regularization terms are non-smooth, while a subgradient type algorithm was used in [4]. It can be verified that the regularization terms in (6) and (7) give very similar values, when the value ρ = 0.5 is used in (6). However, in this work, we prefer the form in (7) because it can lead to a closed form proximal operator, shown in Equation (14). The regularization term in (7) consists of a series of parts of seconddegree polynomial functions, properly placed one next to the other, so as to have the desired regularization effect. In Figure 2a,b, we demonstrate the regularization terms r 2 (φ(x, y)) and r 4 (φ(x, y)), respectively.

Derivation of the Proximal Gradient Algorithm
The cost function in (5) using the regularization term in (7) can be recognized as a cost function that consists of a smooth term (the mean square error) and a non-smooth term (the regularization part). Thus, the method of proximal gradient [17] can be employed for its minimization. The method is more commonly employed in convex optimization problems, but it can also be employed in non-convex problems [19] as it is the case in this work. Defining the vector φ that contains all the phase variables of our optimization problem, we also define and Using the above definitions, the considered optimization problem becomes The proximal gradient algorithm is an iterative procedure that generates the estimates φ (n) and A (n) at the n-th iteration. In particular, for the problem considered here, the algorithm consists of the following two equations where γ is a properly selected positive step size parameter and Prox γg q (·) is the so-called proximal operator, defined by It can be easily verified that each of the N 2 variables in vector θ contributes two positive terms in the last minimization problem, and these two terms are independent of the terms that correspond to different variables. Thus, the minimization can be performed separately for each of the N 2 variables. After some calculations, it can be shown that the proximal operator for each scalar variable θ i , i = 1, . . . , N 2 , in the vector θ is given by where k i = θ i 2pi/q and β = 2γλq N 2 π . Finally, the condition must be satisfied for the previous expressions to be true. In Figure 2c,d, we demonstrate the form of the proximal operators for q = 2 and q = 4, respectively. It can be observed that the effect of the proximal operator is (i) to quantize phase variables that are close to the discrete levels in the set Φ, as shown by the flat regions in Figure 2c,d, and (ii) to apply a small "force" to the rest of the phase variables towards the closest discrete phase level in the set Φ, as shown by the "ramp" regions in Figure 2c,d. Furthermore, it is very interesting to note that the derived proximal operator has a very similar form to the so-called soft-thresholding operator, employed in the context of the sparse coding problem [20].
To sum up, the proposed optimization approach consists of computing Equations (11) and (12) until convergence, where the proximal operator employed in Equation (12) is shown in Equation (14) and it is applied to each of the N 2 phase variables.

Results
In this section, the results of our study are demonstrated. The simulation setup and parameters used were as follows: A square 512 × 512 fSLM device, measuring 10.24 mm along each dimension, was considered. We used monochromatic light with wavelength w = 532 nm. The light reflected by the SLM travels through a 2f setup with focal length f = 100 mm. In the sequel, for the scopes of reducing the considered pixel size, we consider free space propagation of the light along a distance equal to 10 mm, using a two step Frenel propagation method and a magnification parameter m = 1/10, as detailed in [21]. Thus, we arrive to a pixel size of around dx = 0.52 µm at the target plane (see Figure 1).
Since we are considering an optogenetics application, the desired response intensity image D(x, y) was created as five, randomly placed, disks with a radius equal to 10 µm, as it is demonstrated in Figure 3a, where each disk corresponds to the approximate size and shape of a neuronal soma in the mouse cerebral cortex [18]. In the same figure, we demonstrate the resulting images for various approaches, namely, simple phase quantization in Figure 3b, regularized optimization using the smooth regularizer of [11] in Figure 3c, regularized optimization using the non-smooth regularizer of [4] and the subgradient method in Figure 3d, regularized optimization using the proposed regularizer in Equation (7) and the proposed proximal gradient algorithm in Figure 3e and regularized optimization using a proximal gradient algorithm with an increasing λ, where λ is increased from zero to its final value 0.005 linearly over the 5000 iterations, shown in Figure 3f.
Apart from the optical demonstration of the resulting images shown in Figure 3, the methods were also compared in terms of the Peak Signal-to-Noise Ratio (PSNR), defined using the MSE as where the maximum value of the desired response intensity was max x,y D(x, y) 2 = 1, and the MSE is defined as The results obtained can be seen in Table 1. It is evident from these results that the proposed approach offers significant performance improvements against the simple phase quantization approach as well as against the previously reported method [11]. In particular, for q = 2, the proposed approach offers more than 1.8 dB better performance as compared to the simple phase quantization approach, and 1 dB better performance as compared to the method of [11]. For q = 4, the proposed approach offers more than 1 dB better performance as compared to the simple phase quantization approach and about 0.5 dB as compared to the previously reported approach [11]. It is also interesting to note that our previously proposed algorithm in [4] and the approach proposed here give quite similar results, but with the approach proposed here always having better performance. In addition, Table 1 reveals that the approach in which the proximal gradient algorithm that utilizes a gradually increasing value for the regularizer λ, gave the best results. Finally, it can be seen that the performance gains of the various schemes become smaller as the number q of discrete phase levels increases. Note that for all the considered approaches, several experiments were conducted to yield the best values for the parameter λ, and these values also appear in Table 1. The PSNR values obtained for various algorithms as a function of the parameter λ can be seen in Figure 4, for the case q = 2. . Demonstrationof the resulting images (where the fifth root of the intensities multiplied by 1/2 has been used to amplify the details). All the results correspond to q = 2 phase levels. Furthermore, the square regions have been magnified to better demonstrate the results. (a) Desired response intensity image, (b) non-regularized optimization followed by simple quantization to 2 phase levels, (c) regularized optimization using the regularizer in (6) with λ = 0.006 and ρ = 2 (as used in [11]), (d) regularized optimization using the regularizer in (6) with λ = 0.02 and ρ = 0.5 (as used in [4]), (e) regularized optimization using the regularizer in (7) with λ = 0.02 (proposed approach), and (f) regularized optimization using the regularizer in (7) where λ was linearly increased over 5000 iterations from 0 to 0.005.

Discussion
In this work, the problem of designing few-level phase-shifting masks for computergenerated holograms was studied. This problem has applications in optogenetics but also in lithography, where such few-level phase-shifting masks are needed. For the problem at hand, we followed a gradient-based regularized optimization approach. Since the considered cost function is non-differentiable, we developed a novel proximal gradient type algorithm for its minimization. Simulation results were used to verify the effectiveness of the proposed approach.
Future work will focus on the extension of the proposed approach so as to be able to take into account manufacturability constraints, as used in lithography applications. Furthermore, another line of research will focus on the unrolling of the proposed algorithm into a neural network of suitable structure, trained properly so as to derive a scheme with significantly smaller computational complexity.