Open Access
Add to BibSonomy
Abstract
We consider the problem of calculating the eikonal function defined on a certain curved surface from the condition of generating a prescribed irradiance distribution on a target surface. We show that the calculation of the “ray mapping” corresponding to the eikonal function is reduced to the solution of a linear assignment problem (LAP). We propose an iterative algorithm for calculating a refractive optical surface from the condition of generating a prescribed near-field irradiance distribution in a non-paraxial case. The algorithm is based on sequential calculation of eikonal functions defined on curved surfaces using the LAP-based approach. The proposed algorithm is applied to the calculation of refractive optical elements generating uniform irradiance distributions in a rectangular region and in a region in the form of the letters “IPSI” in the case of a circular incident beam. The presented ray-tracing simulations of the designed optical elements demonstrate high efficiency of the proposed iterative algorithm.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
The problem of design of an optical element generating a prescribed irradiance distribution in a certain region belongs to the inverse problems of non-imaging optics and is extremely challenging. In the most cases, this problem is reduced to finding a solution to a fully nonlinear differential equation (NDE) of Monge–Ampère type. The solution of such an NDE is a complex theoretical and computational problem.
In the recent years, methods for the design of optical surfaces based on numerical finite-difference solution of the NDEs of Monge–Ampère type were proposed [1–8]. In these methods, the derivatives are replaced by their finite-difference approximations, and the solution of the NDE is reduced to the solution of a system of nonlinear equations with respect to the values of the sought-for function describing the optical surface and defined on a certain grid. Depending on the complexity of the problem being solved, this system contains from thousands to tens of thousands nonlinear equations. For the solution of the obtained system, iterative methods, e.g. Newton’s method, are used.
In a number of approaches, the Monge–Ampère equation is solved using an explicit finite-difference method based on the replacement of the original equation by a non-stationary parabolic Monge–Ampère equation [9–11]. The paper [12] introduces an iterative least-squares method for solving the Monge–Ampère equation. In [13, 14], the so-called B-spline collocation method was proposed for the solution of NDEs of Monge–Ampère type. This method represents the surface of the optical element as a system of B-splines. The coefficients of these splines are calculated by solving a system of nonlinear equations. In the opinion of the present authors, this method is among the leading methods for the design of refracting and reflecting optical surfaces, which allow one to generate complex irradiance distributions, e.g. photographs of urban landscapes, animals etc.
The formulation of the inverse problems of non-imaging optics as the problems of the solution of NDEs of Monge–Ampère type assumes that the surface of the designed optical element is smooth. This imposes significant restrictions on the class of the irradiance distributions that can be generated by the element designed using the methods of [1–14]. For example, a smooth optical surface cannot generate an irradiance distribution defined on a “disconnected region” (a set of disjoint sub-regions). A typical example of such a region is a word consisting of several letters. Restrictions arise also in the case of generation of irradiance distributions defined on “connected regions”, which are not simply connected or have complex non-smooth boundaries. For instance, the numerical method for the solution of the Monge–Ampère equation [15] used in [9–11] for the design of refractive optical elements works under the assumption that both the optical element and the illuminated region have a square shape. Using this method, it is possible to design an optical element generating a complex irradiance distribution, e.g. an image of a landscape or a logo [9]. However, these images are generated in a square region on a substantially non-zero background. For the application of the method in the case of a non-square target region, it is necessary to define the boundary conditions and define some initial solution of the MA equation satisfying these boundary conditions. This problem is non-trivial, especially in the case when the illuminated region is not simply connected.
A number of inverse problems of non-imaging optics can be formulated as a Monge–Kantorovich mass transportation problem (MTP) with a special cost function. The solution of the MTP allows one to obtain an integrable ray mapping between the coordinates of the incident rays from a light source and the coordinates of the rays in the target plane after refraction by or reflection from the optical surface. The form of the cost function is determined by the type of the problem being solved. In particular, in the problems described by a standard Monge–Ampère equation, the cost function is quadratic [15]. Let us note that an MTP with a quadratic cost function enables describing the problem of the calculation of an optical surface only in the paraxial approximation [9, 16]. In the non-paraxial case, the cost function is no longer quadratic. In particular, in the problem of calculation of the light field eikonal function from the condition of focusing into a prescribed region in the non-paraxial approximation, the cost function is equal to the distance from a point of the source region, in which the eikonal function is specified, to a point of the target region [16]. In the theoretical works [17–20], the form of the cost function for designing mirrors [17, 18] and refractive surfaces [19, 20] generating prescribed far-field irradiance distributions was derived. In [21, 22], the cost functions for the problems of calculating refractive optical elements transforming a given incident beam with a plane wavefront into an output beam with a prescribed irradiance distribution and a plane wavefront were obtained.
The formulation of an inverse problem of non-imaging optics as a Monge–Kantorovich mass transportation problem enables working with discontinuous ray mappings and with the corresponding piecewice-smooth continuous optical surfaces [16–18,23], in contrast to the smooth surfaces satisfying the NDEs of Monge–Ampère type. Indeed, in [16–18], a concept of the so-called weak (generalized) solution is used, which assumes that the optical surface is differentiable not everywhere, but almost everywhere. Therefore, the MTP-based approach to the solution of non-imaging optics problems allows us to extend the class of generated irradiance distributions to the case of “disconnected regions” or regions with complex non-smooth boundaries.
In the discrete case, the MTP can be formulated as a linear assignment problem (LAP) [16, 24]. For the solution of the LAP, efficient polynomial algorithms have been proposed, including Hungarian algorithm [24] and auction algorithm [25]. Let us note that the approach for the solution of the non-imaging optics problems through the solution of an LAP was for the first time proposed in the recent works of the present authors [16, 23]. The design examples presented in [16, 23] demonstrate the applicability of the LAP-based approach for generating irradiance distributions defined on complex-shaped regions. In particular, the paper [23] considers the design of an optical element with a piecewise-smooth continuous optical surface transforming a plane-wavefront beam with a circular cross-section into a plane-wavefront cross-shaped beam. In this case, the known methods based on the numerical solution of an NDE of Monge–Ampère type are either inapplicable or numerically unstable.
In the present work, we consider the problem of the design of a refracting surface generating a prescribed irradiance distribution in the near field in a general non-paraxial case. This problem cannot be directly formulated as an MTP [26, 27] and therefore the previously considered LAP-based methods [16,23] are not applicable. In this work, we propose a new efficient iterative algorithm for solving this problem. At each iteration, we refine the surface G of the designed optical element by calculating an eikonal function defined on this surface and providing the generation of the prescribed irradiance distribution on the target surface F. We show that the eikonal function can be calculated within the LAP-based approach by solving an MTP with the cost function corresponding to the distance between the points on the surfaces G and F. As examples, we design refractive optical elements generating uniform irradiance distributions in a rectangular region and in a disconnected region corresponding to the letters “IPSI” on a zero background. The presented numerical simulation results confirm high efficiency of the proposed iterative method.
Let us note that the proposed iterative method is quite different from the LAP-based approaches considered in the recently published papers [16, 23]. Comparing to the previous paper [16], where a thin optical element approximation was adopted, the proposed iterative algorithm enables solving the problem in an essentially non-paraxial case. Comparing to the paper [23] dedicated to the design of optical element with two free-form surfaces for collimated beam shaping, the proposed iterative method does not assume that the transmitted optical beam is collimated and solves the problem using a single free-form optical surface.
2. Problem statement
Suppose that there are two surfaces G and F defined by the equations z = g(x, y) and z = f (x, y), respectively, in a Cartesian coordinate system (x, y, z). We assume that the surface F is located above the surface G: f(x, y) > g(x, y) (Fig. 1). In Section 5, we will refer to the surface F as the target surface, whereas the surface G will correspond to the surface of an optical element. In what follows, we describe the coordinates (x, y) of the points on the surfaces G and F by the vectors u = (u1, u2) and v = (v1, v2), respectively. For further reasoning, it will be convenient for us to assume that the points u and v are defined in two separate 2D coordinate systems (Fig. 1).
Let us assume that at u ∈ Ω ⊂ ℝ2, an eikonal function Φ(u) and an irradiance distribution EG(u) are defined on the surface G. We denote by p(u) = (p1(u), p2(u), p3(u)) the unit-length 3-vector of the ray in the medium with the refractive index n = 1, which is defined by the eikonal function Φ(u) defined on the surface G. According to the eikonal equation, the vector p(u) of the ray outgoing from a point on the surface G can be expressed through the derivatives of the eikonal function from the following system of equations [28,29]:
The eikonal function Φ(u) defines a “ray mapping” PΦ : u ↦ v in the following way. The image of the point u ∈ Ω is the point v corresponding to the intersection of the ray outgoing from the point (u, g(u)) ∈ G with the surface F (Fig. 1).
Let us define the integrability property of a mapping. In what follows, we say that a certain mapping P : u ↦ v is integrable if there exists an eikonal function Φ(u), u ∈ Ω such that P = PΦ. Let us formulate the integrability condition in a form convenient for the further analysis. The mapping P transforms a point (u, g(u)) on the surface G into the point (P(u), f(P(u))) on the surface F. Consider the 3-vector t = (P(u) − u, f(P(u)) − g(u)). Let us denote its length by ρ(u, P(u)), divide by it and obtain the following unit vector:
If the mapping P is defined by the eikonal function Φ(u), the integrability condition can be obtained from Eq. (1) in the form
where ∇ = (∂/∂u1, ∂/∂u2) is the gradient operator. Let us note that due to the symmetry of the problem with respect to the surfaces G and F, the integrability condition for the inverse mapping P−1 : v ↦ u has a similar form: where Ψ(v) is the eikonal function on the second surface F.Let us formulate the problem of calculation of the eikonal function Φ(u) from the condition of generating a prescribed irradiance distribution EF(v), v ∈ Θ ⊂ ℝ2 on the second surface F (Fig. 1). This problem consists in finding the eikonal function Φ (u) defining such a mapping PΦ : Ω → Θ, so that for any measurable set ω ⊂ Θ the following light flux conservation law holds:
where and are the surface area differentials, and du and dv stand for du1 du2 and dv1 dv2, respectively. Making the change of variables v = PΦ(u) in the second integral in Eq. (4), we obtain where and are the irradiance distributions “projected” onto the regions Ω and Θ, and JPΦ (u) is the Jacobian of the coordinate transformation (v1, v2) = (PΦ,1(u), PΦ,2(u)): According to Eq. (5), the differential form of the light flux conservation law reads asThus, the problem of calculation of the eikonal function Φ(u) is reduced to finding a solution to Eq. (6), where PΦ(u) is defined by Eq. (1).
3. Mass transportation problem
The considered problem of eikonal function calculation can be formulated as a Monge–Kantorovich mass transportation problem with the cost function
equal to the distance between the corresponding points on the surfaces F and G. In order to formulate this MTP, let us consider the transport plans, i.e. the mappings P: Ω → Θ preserving the measure (the light flux) in the sense of Eq. (6). Let us define the following functional on the transport plans:The solution of the MTP consists in finding the mapping P that minimizes (or maximizes) the functional (8) and satisfies the condition (6).
In the recent work [16], it was shown that an integrable mapping in the eikonal calculation problem [see Eq. (2)] can be obtained from the solution of the MTP of Eq. (8). Let us note that the proof of this statement in [16] is given for a special case when the surfaces G and F are two parallel planes. However, by repeating the reasoning presented in [16], one can easily obtain the proof in the general case of arbitrary surfaces G and F. This proof can also be obtained by finding the extremum of the functional of Eq. (8) under the condition of Eq. (6) using the method of Lagrange multipliers [30]. According to this method, a conditional extremum of the functional of Eq. (8) coincides with an ordinary (unconditional) extremum of the following functional:
where the function λ(u) is the Lagrange multiplier. The calculation of the variation of the functional of Eq. (9) was considered in [31]. By equating the variation to zero, one can obtain the integrability condition of Eq. (2) and the light flux conservation law of Eq. (6) [31]. Let us note that the critical points of the functional of Eq. (9) are not necessarily minima or maxima of the initial functional defined by Eq. (8). Therefore, not only the conditional extrema of the initial functional (8) give integrable ray mappings, but also any local extremum or critical point of this functional [31].Thus, the solution of the considered problem of the eikonal function calculation can be reduced to finding a mapping P(u) through the solution of the MTP of Eq. (8) and to the subsequent reconstruction of the eikonal function using Eq. (2).
4. Reduction to the linear assignment problem
The discrete version of the mass transportation problem of Eq. (8) can be formulated as a linear assignment problem (LAP) [24]. Indeed, let us assume that the domains Ω and Θ are divided into equal-flux cells (or approximated by such cells), with the equality being valid for any pair of cells ωi ⊂ Ω, θj ⊂ Θ. In the simplest case, when E0(u) = E0, u ∈ Ω and E(v) = E, v ∈ Θ, the partition of the domains Ω and Θ into N equal-area cells can be obtained by choosing the steps of rectangular grids superimposed on the said domains [16,23]. When the irradiance distributions E0(u) and E(v) are non-uniform, an equal-flux set of N = M1 M2 cells in the domain Ω (or Θ) can be constructed using the following general approach. The domain Ω can be divided by parallel lines u1 = u1,i, i = 1, ..., M1 − 1 into M1 equal-flux sub-domains Ωi. Afterwards, each sub-domain Ωi is further divided into M2 sub-domains by parallel lines , j = 1, ..., M2 − 1. This approach for the construction of equal-flux cell sets is widely used in the so-called source-target mapping methods based on mapping of equal-flux grids between the source and the target areas [32–34]. When the target irradiance distribution E(v) is defined as a complex raster image, one can use Lloyd’s sampling algorithm based on Voronoi tessellations to approximate E(v) by a set of equal-flux cells. In papers [35,36], the reader can find a review of other advanced methods (referred to as image stippling algorithms) that can be used to obtain equal-flux partitions of the domains Ω and Θ.
For the grids (approximations) introduced in the domains Ω and Θ, all energy-conserving mappings P : Ω → Θ can be described by permutations of N numbers (j1, ..., jn), which determine to which cells θj of the domain Θ the cells ω1, ..., ωn are mapped. In terms of the LAP, the mapping of the cell ωi to the cell θj can be interpreted as the assignment of the j-th task to the i-th agent. According to Eq. (8), the cost of the tasks is described by the matrix
where ui, vj are the centers of the cells ωi and θj. Thus, the mapping v = P(u) can be found by solving the following linear assignment problem: where the minimum is sought over all permutations (j1, ..., jn). Let us remind that the LAP of Eq. (11) is simply a discrete approximation of the corresponding MTP 𝒞(P) → min under the condition of the light flux conservation law [Eq. (6)]. For the solution of the LAP, efficient polynomial algorithms are available, in particular, Hungarian algorithm [24] and auction algorithm [25].5. Iterative algorithm for the calculation of the optical element surface
If the mapping P(u) is calculated using the LAP of Eq. (11), the eikonal function Φ(u) can be obtained by integrating Eq. (2) [15,23]. The details of the reconstruction of the eikonal function are presented in Section 6 describing the designed optical elements.
Having reconstructed the eikonal function Φ(u), one can calculate a refractive optical element located just before the surface G and generating the eikonal function Φ(u) on this surface. For the sake of simplicity, let us consider an incident beam with a plane wavefront propagating in the positive direction of the z axis. The irradiance distribution of the beam is described by the function E0(u), u ∈ Ω. As the lower, “non-working”, surface of the optical element, we will use a plane surface perpendicular to the propagation direction of the incident beam. In this case, the upper, “working”, surface of the element can be calculated in the thin optical element approximation [37,38]:
where n0 is the refractive index of the element material. The optical element defined by Eq. (12) will generate the required eikonal function Φ(u) on the surface z = g(u) and, therefore, will provide the generation of the required irradiance distribution E(v), v ∈ Θ on the surface F. However, the thin element approximation used in Eq. (12) is applicable only in the paraxial case.We propose the following iterative algorithm for the calculation of the surface of the optical element generating a prescribed irradiance distribution E(v), v ∈ Θ on the surface F in the general non-paraxial case.
Algorithm
Figure 2 schematically shows two iterations of this algorithm in the problem of calculation of an optical element generating a uniform irradiance distribution in a rectangular region for a circular incident beam. In this example, the surface F corresponds to a certain plane z = z0 > 0. At the first iteration, the eikonal function defined in the plane z = 0 (g(u) = 0) is calculated from the condition of generating a rectangular irradiance distribution at z = z0. The calculation of the mapping P1(u) (the subscript denotes the iteration number) is carried out through the solution of an LAP of Eq. (11) with the cost function . Then, the eikonal function Φ1(u) is calculated using Eq. (2), and the optical surface zel,1(u) is reconstructed using Eq. (12). This surface is shown in the right part of Fig. 2 depicting the second iteration. Then (at the second iteration), the mapping P2(u) is calculated through the solution of an LAP of Eq. (11) with the cost function , and the eikonal Φ2(u) on the surface zel,1(u) is calculated. Next, the surface of the optical element zel,2(u) (not shown in Fig. 2) is calculated from the function Φ2(u), and the process is repeated. Thus, each iteration of the algorithm requires solving an LAP with the appropriate cost function. In the following Section 6, we solve the LAPs using the auction algorithm [25].
Fig. 2 Schematic representation of two iterations of the algorithm for calculating the working surface of the optical element (CF stands for the cost function).
Let us explain the essence of the presented algorithm. Let zel,n−1(u) and zel,n(u) be the optical surfaces obtained at the (n − 1)-th and n-th iterations, respectively. These surfaces are related by a certain operator equation zel,n(u) = T̂(zel,n−1(u)). According to the proposed algorithm, the operator T̂ includes the calculation of the mapping Pn(u) through the solution of an LAP of Eq. (11) with the cost function , the calculation of the eikonal function Φn(u) using Eq. (2), and the reconstruction of the optical surface zel,n(u) using Eq. (12) with g(u) = zel,n−1(u). The eikonal of the incident beam with a plane wavefront on the obtained surface zel,n(u) equals Φinc,n(u) = n0zel,n(u). If at the n-th iteration of the algorithm Φinc,n(u) = Φn(u), then zel,n(u) = zel,n−1(u) = T̂(zel,n−1(u)). The equality zel,n(u) = zel,n−1(u) means that the iterative algorithm has converged and the problem of calculation of the optical element has been solved. Indeed, in this case the calculated eikonal function Φn(u) providing the generation of the prescribed irradiance distribution E(v) on the surface F coincides with the eikonal of the incident beam on the working surface of the optical element. Thus, the solution of the problem is reduced to solving the operator equation zel(u) = T̂(zel(u)). Therefore, we proposed an algorithm that solves this operator equation using the simple iteration method.
Let us note that the ray-tracing in step 4 of the algorithm is carried out to estimate the algorithm convergence in terms of the error δ =‖Ẽ(v) − E(v)‖ representing the difference between the “current” and the required irradiance distributions. Instead, one could control the convergence state by calculating the value of δz =‖zel,n(u) − zel,n−1(u)‖, which represents the difference between the optical surfaces obtained at two consequent iterations.
In the next section, we use the proposed algorithm to design refractive optical elements generating prescribed irradiance distributions on a plane surface F. The presented design examples demonstrate fast convergence of the algorithm and high performance of the designed optical elements.
Fig. 3 Rectangular mesh on a circle (a) and its image in the target rectangle (b) calculated at the first iteration of the algorithm. The inset shows the magnified fragment along with the discrete mapping obtained by solving the LAP of Eq. (11). Only the images of the points lying on the square mesh (a) are shown in the inset.
6. Design examples
6.1. Optical element generating a uniform irradiance distribution in a rectangular region
As a first example, let us consider the calculation of an optical element focusing a circular incident beam with a constant irradiance distribution E0(u) = (πR2)−1, into a rectangle v ∈ Θ = {(v1, v2)||v1| ≤ w1/2, |v2| ≤ w2/2} with a uniform irradiance distribution E(v) = (w1w2)−1 located in the plane z = z0. The wavefront of the incident beam is plane and is perpendicular to the z axis. The calculation was carried out for the following parameters: radius of the incident beam R = 0.5 cm, rectangle dimensions w1 = 5 cm, w2 = 2.5 cm, distance to the target plane z0 = 4 cm, refractive index of the element material n0 = 1.5.
The optical element was calculated using the algorithm described above in Section 5. The mapping P(u) at each iteration was found through the solution of an LAP of Eq. (11). For the calculation of the mapping, the rectangle Θ was represented as a set of N = 1502 = 22500 rectangular cells θi with the size 0.033×0.017 mm2. The circular domain Ω was also approximated by 22500 square cells ωi with the size 0.0059 × 0.0059 mm2. For the solution of the LAP (11), we used an open-source implementation [39] of the auction algorithm [25]. For the considered example, this implementation allowed us to solve the LAP in about ten minutes. Let us note that the computation time can be essentially reduced (up to an order of magnitude) by using a multi-scale strategy similar to the one discussed in [40, 41].
The eikonal function Φ(u) can be calculated from the mapping P(u) by numerical integration of Eq. (2). However, due to the discrete character of the LAP of Eq. (11), the obtained discrete mapping describes the sought-for continuous mapping P(u) with an error. To compensate this error, the eikonal function Φ(u) at each iteration was calculated using the regularized method presented in [23]. In this method, the eikonal function Φ(u) is approximated by a system of two-dimensional B-splines:
where Bm(u1) and Dn(u2) are the B-splines of a certain order q defined in the domain Ω on a rectangular grid containing N1 · N2 vertices. The coefficients of the B-splines cmn at each iteration were determined from the condition that the derivatives of the splines (13) provide a least-squares approximation of the derivatives of the function Φ(u) given in Eq. (2). In the considered example, B-splines of the 4-th order with N1 = N2 = 16 were used for the representation of the function Φ(u).For the sake of illustration, Fig. 3 shows the continuous mapping obtained by solving Eq. (2) where the eikonal function Φ(u) was represented as a system of B-splines [Eq. (13)]. The presented mapping corresponds to the first iteration of the algorithm. The points in the inset show the discrete solution of the LAP (11) with the cost function ρ1. As it is evident from Fig. 3, the representation in Eq. (13) compensates the errors resulting from the discrete form of the LAP and provides a continuous mapping P(u).
Let us note that the eikonal function Φ(u) is defined up to an additive constant. This constant at each iteration was chosen so that the plane z = 0 was tangent to the surface of the optical element defined by Eq. (12) at the point u = 0 [i.e. g(0) = 0].
Fig. 4 Irradiance distributions generated by the designed optical elements at the iterations 1–4 (a). Optical element calculated at the fifth iteration (b) and the corresponding irradiance distribution (c).
Figure 4(a) shows the irradiance distributions generated in the plane z0 = 4 cm by the optical elements calculated at the iterations 1–4. For the calculation of the irradiance distributions, an in-house implementation of the ray-tracing technique was used. The distributions shown in Fig. 4(a) were calculated by tracing 107 rays. Let us note that ray-tracing computation in such a simple geometry including single refractive surface takes of about 2 seconds which is negligible comparing to the time required for solving the LAP. At the first iteration, the generated irradiance distribution is very different form the required one. This is due to the thin element approximation used for the calculation of the optical surface [Eq. (12)]. Indeed, the chosen geometrical parameters of the problem (R = 0.5 cm, w1 = 5 cm, w2 = 2.5 cm, and z0 = 4 cm) correspond to an essentially non-paraxial case. At these parameters, the maximum deflection angle of the incident rays upon the refraction on the surface of the optical element exceeds 31°. At the next iterations, the quality of the generated rectangular distribution monotonically increases. At the fifth iteration [Fig. 4(c)], the generated distribution almost coincides with a perfect uniform rectangular distribution: the RMS deviation of the obtained distribution from a constant value amounts to 3.8%. The optical element calculated at the fifth iteration is shown in Fig. 4(b). Maximum height of the element equals 0.238 cm and is comparable to the element radius R = 0.5 cm.
6.2. Optical element a generating uniform irradiance distribution in the form of the letters “IPSI”
As a second, more challenging example, we designed an optical element generating uniformly illuminated letters “IPSI” on a zero background (Fig. 5). The letters are bounded by a rectangle with the dimensions of 5 × 2 cm2 and are located in the plane z0 = 4 cm. As in the previous example, the incident beam is circular (R = 0.5 cm) and has a constant irradiance.
Fig. 5 Irradiance distributions generated by the optical elements calculated at the iterations 1–4 (a). Optical element calculated at the fifth iteration (b), and the corresponding irradiance distribution (c). The curves on the element surface show the sharp bends caused by the discontinuous character of the mapping P(u).
In this example, the calculation of the eikonal function Φ(u) from the mapping P(u) (step 2 of the iterative algorithm) has a number of specific features. The generated irradiance distribution “IPSI” is “disconnected” because it consists of several separate letters. Therefore, the mapping P(u) is discontinuous. At the same time, the inverse mapping P−1(v) is continuous, in the sense that the close points from the domain Θ (the letters “IPSI”) are mapped to the close points in the domain Ω (the circle). The reconstruction of the eikonal function Φ(u) from a discontinuous mapping in the general case requires the solution of an LAP having a large dimension. Instead, for the calculation of the eikonal function Φ(u), a special method was used, which supports discontinuous mappings P(u). As a first step of this method, the eikonal function Ψ(v), v ∈ Θ defined on the target surface F is calculated from the mapping P−1(v). Note that the mappings P(u) and P−1(v) are obtained simultaneously from the solution of the LAP of Eq. (11). Then, the function Ψ(v) is approximated by a system of two-dimensional B-splines. Similarly to Eq. (13), the splines were defined on a rectangular grid of N1 · N2 vertices located in the rectangular region containing the letters “IPSI”. The coefficients of the B-splines are determined from the condition that the spline derivatives provide a least-squares approximation of the derivatives of the function Ψ(v) given in Eq. (3). Let us note that if the surface F is plane (z = z0), the second term in Eq. (3) vanishes. Finally, the calculated function Ψ(v) is used to obtain the eikonal function Φ(u) [16]:
Equation (14) implies that the eikonal on the surface G is represented as the envelope of a family of eikonal functions of lenses with the foci lying in the domain Θ [16]. This makes the described method capable of finding a piecewise-smooth continuous eikonal function Φ(u) corresponding to a discontinuous mapping P(u). Let us note that a somewhat similar approach was used in [23] for designing two-surface optical elements for collimated beam shaping.The calculation of the optical element generating the illuminated letters “IPSI” was carried out using the algorithm presented in Section 5 and the special eikonal calculation method based on Eq. (14) described above. In the calculations, both the circular source domain Ω and the target domain Θ were represented by sets of N = 2502 = 625000 square cells. At the chosen N value, the solution of the LAP takes of about an hour. The rest operations associated with the calculation of the eikonal function and of the optical surface take 2–3 minutes.
The irradiance distributions generated in the plane z0 = 4 cm by the optical elements calculated at the iterations 1–4 are shown in Fig. 5(a). As in the first example, the distribution generated at the first iteration is very different from the required one because of the essentially non-paraxial character of the problem being solved. Indeed, the maximum deflection angle of the incident rays upon the refraction on the optical surface exceeds 30°. Similarly to the previous example, at the next iterations the quality of the generated distribution increases, and at the fifth iteration [Fig. 5(c)] the obtained irradiance distribution almost coincides with the prescribed one. The resulting optical element calculated at the fifth iteration is shown in Fig. 5(b). Maximum height of the element amounts to 0.257 cm and is comparable to the element radius R = 0.5 cm.
Figure 6 shows the Gaussian curvature of the designed optical surface. Along the red lines, the curvature diverges, thus indicating the sharp bends on the surface of the calculated optical element. These bends result from the discontinuities of the mapping P(u) and are also shown with black lines in Fig. 5(b).
The presented example demonstrates one of the advantages of the proposed method, namely, the possibility of calculating piecewise-smooth continuous optical surfaces. In this case, the known methods based on numerical solution of an NDE of Monge–Ampère type [1–14] are either not applicable, or are numerically unstable. Indeed, in these methods, the optical surface is assumed to be a smooth, twice differentiable function. This significantly limits the class of the generated distributions, and, in particular, does not allow one to generate “disconnected” distributions, e.g. text on a zero background.
7. Conclusion
We considered the application of an LAP-based approach for the calculation of an eikonal function defined on a certain curved surface from the condition of generating a prescribed irradiance distribution on another surface. A new iterative algorithm for calculating a refractive surface generating a prescribed irradiance distribution was proposed. The algorithm is based on sequential calculation of eikonal functions defined on curved surfaces. We demonstrated that this algorithm corresponds to the solution of an operator equation describing the problem of the optical surface design using the simple iteration method. High performance of the proposed iterative algorithm is confirmed by the design examples of the optical elements generating uniform irradiance distributions in a rectangle and in the region in the form of the letters “IPSI” on a zero background.
The presented iterative algorithm can also be applied for the design of reflecting optical surfaces. In this case, the calculation of the surface is performed using the same Eq. (12) with n0 = −1. It is important to note that the problem of calculation of a reflecting surface generating a given distribution in the near field cannot be formulated as a mass transportation problem [27].
Therefore, the direct (non-iterative) application of the LAP-based approach is not available for the solution of this problem.
Funding
Russian Science Foundation (18-19-00326); Ministry of Science and Higher Education of the Russian Federation (State assignment to the FSRC “Crystallography and Photonics” RAS).
Acknowledgments
The development of the design method and the investigation of the designed elements was supported by Russian Science Foundation, and the implementation of the software for the optical element design was supported by Ministry of Science and Higher Education of the Russian Federation.
References
1. R. Wu, Z. Feng, Z. Zheng, R. Liang, P. Benítez, and J. C. Miñano, “Design of freeform illumination optics,” Laser & Photonics Rev. 12(7), 1700310 (2018). [CrossRef]
2. R. Wu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “A mathematical model of the single freeform surface design for collimated beam shaping,” Opt. Express 21(18), 20974–20989 (2013). [CrossRef] [PubMed]
3. R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Xiu, “Freeform illumination design: a nonlinear boundary problem for the elliptic Monge–Ampère equation,” Opt. Lett. 38(2), 229–231 (2013). [CrossRef] [PubMed]
4. R. Wu, Y. Zhang, M. M. Sulman, Z. Zheng, P. Benítez, and J. C. Miñano, “Initial design with L2 Monge–Kantorovich theory for the Monge–Ampère equation method in freeform surface illumination design,” Opt. Express 22(13), 16161–16177 (2014). [CrossRef] [PubMed]
5. Y. Ma, H. Zhang, Z. Su, Y. He, L. Xu, X. Lui, and H. Li, “Hybrid method of free-form lens design for arbitrary illumination target,” Appl. Opt. 54(14), 4503–4508 (2015). [CrossRef] [PubMed]
6. X. Mao, S. Xu, X. Hu, and Y. Xie, “Design of a smooth freeform illumination system for a point light source based on polar-type optimal transport mapping,” Appl. Opt. 56(22), 6324–6331 (2017). [CrossRef] [PubMed]
7. R. Wu, S. Chang, Z. Zheng, L. Zhao, and X. Liu, “Formulating the design of two freeform lens surfaces for point-like light sources,” Opt. Lett. 43(7), 1619–1622 (2018). [CrossRef] [PubMed]
8. C. R. Prins, J. H. M. ten Thije Boonkkamp, J. van Roosmalen, W. L. IJzerman, and T. W. Tukker, “A Monge–Ampère-solver for free-form reflector design,” SIAM J. Sci. Comput. 36(3), B640–B660 (2014). [CrossRef]
9. C. Bösel and H. Gross, “Ray mapping approach for the efficient design of continuous freeform surfaces,” Opt. Express 24(13), 14271–14282 (2016). [CrossRef] [PubMed]
10. C. Bösel and H. Gross, “Single freeform surface design for prescribed input wavefront and target irradiance,” J. Opt. Soc. Am. A 34(9), 1490–1499 (2017). [CrossRef]
11. C. Bösel and H. Gross, “Double freeform illumination design for prescribed wavefronts and irradiances,” J. Opt. Soc. Am. A 35(2), 236–243 (2018). [CrossRef]
12. C. R. Prins, R. Beltman, J. H. M. ten Thije Boonkkamp, W. L. IJzerman, and T. W. Tukker, “A least-squares method for optimal transport using the Monge–Ampère equation,” SIAM J. Sci. Comput. 37(6), B937–B961 (2015). [CrossRef]
13. K. Brix, Y. Hafizogullari, and A. Platen, “Designing illumination lenses and mirrors by the numerical solution of Monge–Ampère equations,” J. Opt. Soc. Am. A 32(11), 2227–2236 (2015). [CrossRef]
14. K. Brix, Y. Hafizogullari, and A. Platen, “Solving the Monge–Ampère equations for the inverse reflector problem,” Math. Model. Methods Appl. Sci. 25(6), 803–837 (2015). [CrossRef]
15. M. M. Sulman, J. F. Williams, and R. D. Russell, “An efficient approach for the numerical solution of the Monge– Ampère equation,” Appl. Numer. Math. 61(3), 298–307 (2011). [CrossRef]
16. L. L. Doskolovich, A. A. Mingazov, D. A. Bykov, E. S. Andreev, and E. A. Bezus, “Variational approach to calculation of light field eikonal function for illuminating a prescribed region,” Opt. Express 25(22), 26378–26392 (2017). [CrossRef] [PubMed]
17. T. Glimm and V. Oliker, “Optical design of single reflector systems and the Monge–Kantorovich mass transfer problem,” J. Math. Sci. 117(3), 4096–4108 (2003). [CrossRef]
18. X.-J. Wang, “On the design of a reflector antenna II,” Calc. Var. 20(3), 329–341 (2004). [CrossRef]
19. C. E. Gutiérrez, “Refraction problems in geometric optics,” Fully Nonlinear PDEs in Real and Complex Geometry and Optics, Vol. 2087 of the Series Lecture Notes in Mathematics (Springer, 2014), pp. 95–150. [CrossRef]
20. C. E. Gutiérrez and Q. Huang, “The refractor problem in reshaping light beams,” Arch. Ration. Mech. Anal. 193, 423–443 (2009). [CrossRef]
21. J. Rubinstein and G. Wolansky, “Intensity control with a free-form lens,” J. Opt. Soc. Am. A 24(2), 463–469 (2007). [CrossRef]
22. V. Oliker, “Designing freeform lenses for intensity and phase control of coherent light with help from geometry and mass transport,” Arch. for Ration. Mech. Analysis 201(3), 1013–1045 (2011). [CrossRef]
23. L. L. Doskolovich, D. A. Bykov, E. S. Andreev, E. A. Bezus, and V. Oliker, “Designing double freeform surfaces for collimated beam shaping with optimal mass transportation and linear assignment problems,” Opt. Express 26(19), 24602–24613 (2018). [CrossRef]
24. J. Munkres, “Algorithms for the assignment and transportation problems,” Journal of the Society for Industrial and Applied Mathematics 5(1), 32–38 (1957). [CrossRef]
25. D. P. Bertsekas, “The auction algorithm: A distributed relaxation method for the assignment problem,” Ann. Oper. Res. 14(1), 105–123 (1988). [CrossRef]
26. C. E. Gutiérrez and Q. Huang, “The near field refractor,” Ann. Inst. Henri Poincaré C Non Linear Anal. 31(4), 655–684 (2014). [CrossRef]
27. S. A. Kochengin and V. I. Oliker, “Determination of reflector surfaces from near-field scattering data,” Inverse Probl. 13(2), 363–373 (1997). [CrossRef]
28. Yu. A. Kravtsov and Yu. I. Orlov, Geometrical Optics of Inhomogeneous Media (Springer-Verlag, 1990). [CrossRef]
29. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University, 1999). [CrossRef]
30. L. C. Evans, “Partial differential equations and Monge–Kantorovich mass transfer,” in Current Developments in Mathematics, R. Bott, A. Jaffe, D. Jerison, G. Lusztig, I. Singer, and S.-T. Yau, eds. (International Press of Boston, 1999).
31. A. A. Mingazov, D. A. Bykov, L. L. Doskolovich, and N. L. Kazanskiy, “Variational interpretation of the eikonal calculation problem from the condition of generating a prescribed irradiance distribution,” Comput. Opt. 42(4), 567–573 (2018).
32. X. Mao, H. Li, Y. Han, and Y. Luo, “Polar-grids based source-target mapping construction method for designing freeform illumination system for a lighting target with arbitrary shape,” Opt. Express 23(4), 4313–4328 (2015). [CrossRef] [PubMed]
33. Y. Ding, X. Liu, Z.-R. Zheng, and P.-F. Gu, “Freeform LED lens for uniform illumination,” Opt. Express 16(17), 12958–12966 (2008). [CrossRef] [PubMed]
34. L. L. Doskolovich, N. L. Kazansky, S. I. Kharitonov, and V. A. Soifer, “A method of designing diffractive optical elements focusing into plane areas,” J. Mod. Opt. 43(7), 1423–1433 (1996). [CrossRef]
35. M. Balzer, T. Schölmer, and O. Deussen, “Capacity-constrained point distributions: a variant of Lloyd’s method,” ACM Trans. Graph. 28(3), 86 (2009). [CrossRef]
36. F. de Goes, K. Breeden, V. Ostromoukhov, and M. Desbrun, “Blue noise through optimal transport,” ACM Trans. Graph. 31(6), 171 (2012). [CrossRef]
37. J. Turunen and F. Wyrowski, Diffractive Optics for Industrial and Commercial Applications (Wiley, 1997).
38. V. A. Soifer, V. V. Kotlyar, and L. L. Doskolovich, Iterative Methods for Diffractive Optical Elements Computation (Taylor & Francis, 1997).
39. Implementation of Bertsekas’ auction algorithm. http://www.mathworks.com/matlabcentral/fileexchange/48448.
40. A. M. Oberman and Y. Ruan, “An efficient linear programming method for optimal transportation,” https://arxiv.org/abs/1509.03668.
41. B. Schmitzer and C. Schnörr, “A hierarchical approach to optimal transport,” in Scale Space and Variational Methods in Computer Vision, A. Kuijper, K. Bredies, T. Pock, and H. Bischof, eds. (Springer, 2013), pp. 452–464. [CrossRef]
