1. Introduction

The problem of design of an optical element generating a required illuminance distribution belongs to the class of extremely challenging inverse problems of nonimaging optics [1]. In many cases such problems can be reduced to finding solutions to nonlinear elliptic partial differential equations of second order. In optical design, direct numerical solutions of such equations in special cases have been considered recently in [2, 3]. In general, such inverse problems of optical design are approached using various iterative techniques [4–13]. One of the universal iterative methods widely used for the design of nonimaging optical elements is the so-called supporting quadric method (SQM) [10–13]. Originally, SQM was developed for solving problems in design of mirrors and lenses transforming incoming beam of light so that the outgoing beam illuminates a given set of points with prescribed intensity distribution. Depending on the type of the problem being solved, the required optical surface is defined as an envelope of a family of quadrics, such as paraboloids, ellipsoids, or hyperboloids. In particular, when designing a mirror for illuminating a target set in the far-field, the required mirror is represented as a set of segments each of which is a piece of a paraboloid. In the near-field case the required mirror is a set of segments of ellipsoids [10, 11]. In [12], one of the present authors used SQM to give a mathematical treatment of the problem of design of a freeform segmented refracting lens transforming an incident collimated beam with given intensity into an output light beam illuminating a given target set with a prescribed intensity distribution. In either case, each segment of the appropriate quadric is defined by certain parameters determined by a special iterative procedure whose convergence is established in a mathematically rigorous manner. The same ideas were applied in [13] to design problems in which quadrics should be replaced by Cartesian ovals.

In this work, a modification of SQM is proposed enabling the method to be applied to the problem of calculating the eikonal function of a light field (specified in a certain plane) providing the focusing into a given discrete set of points. This problem is relevant in the design of diffractive optical elements (DOEs), especially of multifocal diffractive lenses, beam splitters, and diffusers [14–18]. Indeed, the height of the DOE microrelief is in direct proportion to the eikonal function [14]. In addition, the eikonal function can be used to reconstruct refractive (or reflective) optical element generating this eikonal distribution at a preset illuminating beam [19]. This makes the proposed approach well-suited for design of refractive (or reflective) optical elements under the assumption that illuminance distribution of the incident beam is changed insignificantly upon transmission through the optical element [19]. In our setting, we determine the eikonal function providing the focusing into a given set of points as an envelope of a two-parameter set of eikonals of convergent spherical beams with different foci. The proposed method is remarkable for its computational simplicity. As an example, we calculate the eikonal function providing the focusing into an array of points corresponding to a calculator keyboard image. The obtained eikonal function is then used to construct a refractive optical element and DOEs. We demonstrate that the proposed method can be used for designing DOEs producing prescribed continuous intensity distributions within the scalar theory of diffraction.

2. The eikonal function for focusing into a two-dimensional region and discrete point set

The problem of calculating the eikonal for generating a required discrete illuminance distribution can be formulated as follows. Let $E_0\left(u\right),u\in G$ be the illuminance distribution in a region (aperture) $G$ in the plane $z=0$, where $u=\left(u,v\right)$ are the Cartesian coordinates in that plane. We wish to determine the eikonal distribution $\Psi \left(u\right),u\in G$ of a light field providing the focusing into a given set of points $x_i=\left(x_i,y_i\right),\text{}\text{}i=1,\mathrm{...},N$ in the plane $z=f,f>0$ with prescribed energy distribution $I_i,\text{}\text{}i=1,\mathrm{...},N$ (Fig. 1). Note that in the framework of geometrical optics the light ray directions are determined by the partial derivatives of the eikonal function: $\nabla \Psi \left(u\right)=\left(\partial \Psi \left(u\right)/\partial u,\partial \Psi \left(u\right)/\partial v\right)$. Assuming that the surrounding medium has unity refractive index, we can write the unit vector of a ray originating from a point $u$ as

 

Fig. 1 Geometry of the problem of focusing into a set of focal points.

Download Full Size | PDF

On the basis of the function $\Psi \left(u\right)$, a refractive or diffractive optical element focusing into a given set of points can be designed. Several such examples are given below in sections 3 and 4.

This discrete problem of focusing into a set of points may be considered as a particular case of focusing into a desired one or two-dimensional (2D) region. Actually, assume that a mesh $x_i=\left(x_i,y_i\right),\text{}\text{}i=1,\mathrm{...},N$ approximates a 2D region $D$. Let us calculate the eikonal function in the plane$z=0$, requiring that the rays originating in G arrive in the region $D$. For that, let us write the eikonal function providing the focusing at point $x$ in the plane $z=f$:

It can be easily shown that Eq. (2) corresponds to the eikonal of a spherical beam converging to point $x$, written at $z=0$. In this case the ray directions (Eq. (1)) are given by unit vectors $p\left(u\right)=\left(x-u,f\right)/\sqrt{{\left(u-x\right)}^2+f^2}$ directed at the focus. Since such a beam can be formed by an ideal lens, in what follows we refer to Eq. (2) as the “lens eikonal function”. In Eq. (2), the coordinates $x=\left(x,y\right)$ are interpreted as parameters, with ${\Psi }_f\left(x\right)$ defining a constant equal to the eikonal value at the focal point.

It is proposed that the eikonal function for focusing into a desired region $D$ should be defined as the envelope of a two-parameter family of functions in Eq. (2) with respect to the parameters $\left(x,y\right)\in D$. By definition, this envelope is a function tangent to each eikonal function from Eq. (2) at a certain point. Accordingly, the partial derivatives of the function defining the envelope are equal to the partial derivatives of the eikonal of the corresponding lens at the common point of tangency. In geometrical optics, the ray direction is defined by partial derivatives of the eikonal function, Eq. (1). Hence, the rays defined by the envelope will be directed to the points of region D where the foci of spherical beams defined by Eq. (2) are located.

Thus, the calculation of the eikonal for focusing into a desired region $D$ is reduced to constructing an envelope of a two-parameter family of lens eikonals defined by Eq. (2). In Cartesian coordinates, the envelope is given by [20]

The general representation of the eikonal function in Eqs. (4)–(6) depends on the ${\Psi }_f\left(x\right)$ function which is to be determined from the requirement that the energy distribution in the region of focusing $D$ is equal to a prescribed distribution.

In the discrete version of the problem, light is focused into a set of focal points $x_i,\text{}\text{}i=1,\mathrm{...},N$ of the region $D$. So, Eqs. (4)–(6) take the form

3. Solution of the discrete problem of focusing into a set of points

Now, we explain how to determine the values ${\Psi }_{f,i},\text{}\text{}i=1,\mathrm{...},N$ using the requirement that the desired energy distribution $I_i,\text{}\text{}i=1,\mathrm{...},N$ is to be generated at points $x_i=\left(x_i,y_i\right),\text{}\text{}i=1,\mathrm{...},N$. This problem can be addressed using an analog of the supporting quadric method described in Refs [10, 11]. for design of mirrors. In this case, the required mirror is represented as a set of segments of ellipsoids, with the number of segments equal to that of the prescribed focal points. For each ellipsoid, one focus is located at the light source while the other is at one of the focal points $x_i,\text{}\text{}i=1,\mathrm{...},N$. The large axes of the ellipsoids are calculated by an iterative algorithm using the condition that the required energy distribution must be achieved at the focal points [10, 11]. Note that the large axes of the ellipsoids are equal to the eikonal at the focal points (i.e. the path length).

The problem of determination of the eikonal function from the condition of focusing into points with prescribed energy distribution $I_i,\text{}\text{}i=1,\mathrm{...},N$ is nearly identical to the design problem solved by the supporting quadric method of Refs [10, 11]. The functions describing ellipsoids are replaced by the lens eikonal functions:

The iterative algorithm for the calculation of the values ${\Psi }_{f,i}$ includes multiple solution of a direct problem that involves the calculation of energy at the focal points based on a ray tracing technique. When designing mirrors composed of segments of ellipsoids, on each iteration, the rays from the source are traced and the distance from the source to each ellipsoid along the current ray is determined. Two kinds of mirrors can be constructed: the first one is built from ellipsoidal segments nearest to the source (along the ray) and the second one is built from segments which are the farthest from the source. The minimum condition is used to design elliptic mirrors (with the rays intersecting the optical axis), whereas the maximum condition is used for hyperbolic mirrors (with the rays not intersecting the optical axis).

For our discrete problem, it follows from Eqs. (7)–(9) that in order to trace a ray from the point $u$ it suffices to determine the minimal or maximal element of the array ${\Psi }_{l,i}\left(u\right),i=1,\mathrm{...},N$.

In the general case, the energy distribution between the focal points for specified values of ${\Psi }_{f,i},i=1,\mathrm{...},N$ is calculated as follows. Let $u_j=\left(u_j,v_j\right),\text{}\text{}j=1,\mathrm{...},N_G$ denote the nodes of a rectangular mesh with cell size ${\delta }_u\times {\delta }_v$ that approximates the domain $G$ in which the eikonal function is defined. Each ray is assumed to carry the energy $E_j=E_0\left(u_j\right){\delta }_u{\delta }_v,\text{}\text{}j=1,\mathrm{...},N_G$. For every starting point $u_j$ of the ray, the destination point is derived by finding the minimal (maximal) element of the array ${\Psi }_{l,i}\left(u_j\right),i=1,\mathrm{...},N$. At the destination point (denoted by index $m=\underset{i}{\arg \min }{\Psi }_{l,i}\left(u_j\right)$), the energy is incremented by the magnitude $E_j$.

In the case of mirrors, a detailed description of the monotonically converging iterative algorithm for calculating the ellipsoid parameters can be found in Refs [10, 11]. The said algorithm can be extended without any alterations to calculate the values of ${\Psi }_{f,i},i=1,\mathrm{...},N$ in Eq. (10). At each step of the algorithm, the energy distribution between the focal points is calculated, followed by the correction of the values of ${\Psi }_{f,i},i=1,\mathrm{...},N$. If the energy at a point $x_i$ is lower than the specified value $I_i$ then, based on the minimum condition in Eq. (8), the value of ${\Psi }_{f,i}$ needs to be decreased by a magnitude $\Delta$. With ellipsoids, a decrease in the large axis makes the ellipsoid closer to the source, in this way increasing the number of rays directed to the destination point (second focus of the ellipsoid).

Thus, the eikonal function to focus the beam into a set of focal points can be derived using the following iterative algorithm:

Note that the software implementation of the algorithm involves a single computation of the matrix ${\Psi }_{l,ij}$ at step 3. Afterwards, the solution process consists in the iteration of steps 4–6. In that case, the energy distribution between the focal points in step 4 is deduced by calculating the minimal (maximal) elements in the columns of the matrix ${\Psi }_{l,ij}$, whereas the matrix ${\Psi }_{l,ij}$ is recalculated in step 6 via altering the row elements by values ${\Delta }_i$.

For a detailed description of the method for correcting the values ${\Psi }_{f,i},i=1,\mathrm{...},N$ in step 5, which enables the error to decrease monotonically during the iteration process, see Refs [10, 11]. At the same time, the method is characterized by a relatively slow convergence rate, especially for a large number of focal points. The analysis conducted by the authors has shown the following simple heuristic approach to the correction of the values of ${\Psi }_{f,i}$ to be efficient:

4. Design example

To demonstrate the capability of the proposed method to generate complex images, we calculated the eikonal function for the generation of a calculator image defined on a 2D equidistant mesh of $56\times 82$ pixels (Fig. 2). Note that Fig. 2 contains $N_{\text{im}}=1414$ white pixels that form the required image. The eikonal was calculated from the condition of focusing into a discrete set of points corresponding to white pixels in Fig. 2. The points are located in the plane $z=100$mm on a mesh with step ${\delta }_{\text{xy}}=0.\text{2586}\text{mm}$, the energy distribution between the points is assumed to be constant $\left(I_i=I_0,\text{}\text{}i=1,\mathrm{...},N\right)$. The eikonal function was specified in the plane $z=0$ in a square-shaped domain $G$ with the side $a=8.192\text{mm}$ on an equidistant mesh of $256\times 256$ pixels with step $d_{\text{uv}}=32\mu m$. The illuminance $E_0\left(u\right)$ in the domain $G$ was assumed to be constant.

 

Fig. 2 An image of a calculator on an equidistant mesh of $56\times 82$ pixels.

Download Full Size | PDF

The calculated eikonal function is depicted in Fig. 3 and provides the generation of an almost uniform illuminance distribution at the desired points, with the root-mean-square deviation from the specified uniform distribution being less than 6% (Fig. 4). Let us note that the resolution of Fig. 3 does not allow to see the details of the constructed eikonal function. In the next section, a central fragment of this function is shown which demonstrates that it consists of multiple segments connected in a continuous manner.

 

Fig. 3 The eikonal function to focus into a set of points in Fig. 2.

Download Full Size | PDF

 

Fig. 4 Illuminance distribution generated in the target plane $z=100\text{mm}$ using the eikonal function in Fig. 3.

Download Full Size | PDF

On the basis of the derived eikonal function, a refractive surface was reconstructed. Exact analytical formulas to calculate a refractive surface that would generate a required eikonal distribution were presented in Ref [19]. At the same time, for the chosen “paraxial” parameters, the height of the refractive surface can be obtained using the thin element approximation [14]:

To verify the validity of the design technique given above, the optical element defined by Eq. (12) and the eikonal in Fig. 3 was simulated using the commercial software TracePro implementing a ray tracing technique [21]. For the simulation, the surface defined by Eq. (12) was approximated by a system of non-uniform rational Bezier splines (NURBS) using the computer-aided design software Rhinoceros [22]. The input surface of the optical element was assumed to be the plane $z=0$ perpendicular to the propagation direction of a uniform incident beam. Figure 5 shows the refractive optical element and the simulated illuminance pattern generated by the element. The illuminance distribution was calculated using 1,000,000 rays and confirms that the incident beam is focused into a desired pattern. A certain non-uniformity in the energy distribution between the focal points is due to the spline-aided approximation of the element surface and the inherent random noise of the ray tracing technique.

 

Fig. 5 A refractive optical element to focus into a set of points in Fig. 2 and the generated illuminance pattern obtained using the TracePro software.

Download Full Size | PDF

5. Analysis of the diffraction effects

The problem of designing phase diffractive optical elements (DOEs) is often formulated as that of calculating the eikonal function (or a related phase function $\phi \left(u\right)=k\Psi \left(u\right)$, where $k=2\pi /\lambda$ is the wavenumber and $\lambda$ is the free-space wavelength) from the condition of generating a desired illuminance pattern [14]. The DOE microrelief height is given by

The eikonal function in Fig. 3 to focus into a set of points in Fig. 2 was calculated in the framework of geometrical optics. In this context, it is worth studying the performance of a DOE defined by Eq. (13) within the scalar theory of diffraction. Below, the illuminance distribution generated by the eikonal function of Fig. 3 is analyzed using the Fresnel-Kirchhoff approximation. The complex amplitude of the light field in the DOE plane $z=0$ is defined as $w_0\left(u\right)=\sqrt{E_0\left(u\right)}\exp \left[ik\Psi \left(u\right)\right],u\in G$. In the Fresnel-Kirchhoff approximation, the illuminance distribution in the plane $z=f$ takes the form [14]

Fig. 6 (a) depicts the central fragment of the function${\Psi }_{\text{add}}\left(u\right)$ modulo $\lambda =532\text{nm}$. For the calculation of the integral in Eq. (16), the said eikonal function was calculated using Eqs. (7) and (9) on an equidistant mesh of $4096\times 4096$ pixels with the step $d_{\text{uv}}=2\text{µm}$. The resulting illuminance pattern obtained in the Fresnel approximation and corresponding to the calculator image is shown in Fig. 6(d). The illuminance $E\left(x\right)$ in Fig. 6(d) was calculated on an equidistant mesh composed of $4096\times 4096$ pixels, with the step $d_{\text{xy}}=\lambda f/a$. Note that in the Fresnel approximation, the magnitude $2\Delta =2\lambda f/a$ corresponds to a focal spot size generated by a lens of focus f and square aperture G with the side length $a$ [14]. For the given parameters ($\lambda =532\text{nm}$, $f=100\text{mm}$, $a=8.192\text{mm}$), $d_{\text{xy}}=\Delta \approx 6.5\text{µm}$. In contrast with Fig. 5, individual points forming the keyboard image are significantly less evident in Fig. 6(d). The reason is that due to diffraction the points turned into finite-size spots, which have merged to form an almost continuous image of the calculator. Let us derive a qualitative estimate of an average size of the diffraction spots for the aforesaid example. In accordance with the employed calculation method, the eikonal function in Figs. 3 and 6(a) corresponds to a set of lens eikonal segments (subapertures) of equal-area that are connected in a continuous manner. As an estimate of the focal-spot average size, we will use the quantity

 

Fig. 6 (a)–(c) Fragments of the eikonal function ${\Psi }_{\text{add}}\left(u\right)$, $0.5{\Psi }_{\text{add}}\left(u\right)$, $0.25{\Psi }_{\text{add}}\left(u\right)$ in the central regions of size $3\times 3{\text{mm}}^2$; (d)–(f) Illuminance patterns generated in the Fresnel-Kirchhoff approximation by the eikonal functions in Fig. 6 (a)–(c).

Download Full Size | PDF

These arguments show that the proposed method is suitable for calculating the eikonal function when designing DOEs intended to form given continuous illuminance patterns. Note that the continuity of the distribution is achieved by properly choosing the distance between points while calculating the eikonal function. Let us study the magnitude of diffraction effects with decreasing the mesh step ${\delta }_{\text{xy}}$ relatively to the value of ${\Delta }_{\text{sub}}$. It can be easily shown [14] that in the paraxial approximation, the eikonal function

Within the scalar theory of diffraction, iterative and gradient methods have been widely used to design DOEs [14, 23]. The most popular methods are the Gerchberg-Saxton (GS) algorithm and its modifications [14, 23]. Both the convergence of the GS algorithm and the quality of the obtained solution essentially depend on the initial guess for the DOE phase function. Let us analyze the performance of the GS algorithm using the geometrical optics solution $\phi \left(u\right)=k{\Psi }_{\text{add}}\left(u\right)$ as an initial guess. As the required distribution, the continuous image of the calculator in Fig. 2 is used with the points being replaced by square pixels of size ${\delta }_{\text{xy}}$.

Figs. 7 (а)-(c) depict the eikonal functions derived after 30 iterations of the GS algorithm. The eikonal functions in Figs. 6 (а)-(c) were used as initial guesses. The illuminance distributions generated by the eikonal functions in Figs. 7 (а)-(c) are presented in Figs. 7 (d)–(f), showing the required calculator images at different scales. In all three cases, the root-mean-square deviation values of the numerically calculated illuminance patterns from the desired ones do not exceed 5%. Thus, the use of the geometrical optics solution as the initial guess enables the fast convergence of the iterative algorithm. Note, for comparison, that with a random initial guess of the phase, the root-mean-square deviation of the numerically obtained illuminance pattern from the required one stays in the range of 10–15% after 50–100 iterations. The authors believe that the problem of choosing the initial guess is of special significance when generating required illuminance patterns on non-planar or inclined planar surfaces. In that case, the relation between the field in the aperture region and the target region cannot be expressed by the Fourier integral (16), which can be calculated using known fast algorithms. At the same time, as we indicated at the end of Section 1, the proposed algorithm for calculating the eikonal function can be extended almost without change to the case when the target points are located on an arbitrary surface.

 

Fig. 7 (a)–(c) Fragments of the eikonal functions after 30 iterations of the GS algorithm; (d)–(f) Illuminance patterns generated using the Fresnel-Kirchhoff approximation for the eikonal functions in Fig. 7 (a)–(c).

Download Full Size | PDF

6. Conclusion

An iterative method for calculating the eikonal function of a light field from the condition of focusing into a set of points was proposed. The developed method is a modification of the supporting quadric method developed for designing mirrors and refractive optical elements generating discrete intensity or illuminance patterns on a given set of points. The method has been developed as a particular solution of the problem of calculating the eikonal to focus light into a desired region, with the eikonal function being expressed as the envelope of a two-parameter family of eikonals of lenses. High efficiency of the proposed approach is demonstrated by calculating the eikonal function to form a complex illuminance pattern corresponding to a calculator keyboard image. The method is suited for the design of DOEs to generate desired continuous illuminance patterns using the scalar theory of diffraction. The continuity of the resulting patterns has been achieved by choosing a proper distance between the focal points while calculating the eikonal function. The eikonal function offers a good initial guess when designing DOEs using the iterative Gerchberg-Saxton algorithm.