Variational approach to calculation of light field eikonal function for illuminating a prescribed region

Leonid L. Doskolovich; Albert A. Mingazov; Dmitry A. Bykov; Evgeniy S. Andreev; Evgeni A. Bezus

doi:10.1364/OE.25.026378

1. Introduction

The problem of designing an optical element generating a required illuminance distribution in a certain plane belongs to the inverse problems of non-imaging optics and is extremely challenging. In the most cases, this problem can be reduced to finding a solution to a nonlinear partial differential equation (PDE) of elliptic type. The design methods based on ‘direct’ numerical solution of an equation of this type appeared only in the recent years [1–4]. In these methods, the derivatives are replaced by their finite-difference approximations, and the solution of the PDE is reduced to the solution of a system of nonlinear equations with respect to the values of the function describing the optical surface and defined on a certain grid. Depending on the complexity of the problem being solved, this system typically consists of several thousands to several tens of thousands equations. For the solution of the obtained system, iterative methods, e.g. Newton’s method, are used. A common drawback of this approach is high computational complexity of solving a system of nonlinear equations, as well as the choice of the initial approximation to the solution of the system. Therefore, many other methods are used to solve inverse problems of non-imaging optics.

One of the commonly used iterative methods is the supporting quadric method (SQM) [5–11]. The method is intended for the design of mirrors or refractive surfaces generating discrete illuminance distributions consisting of a set of points. In this method, an optical surface is represented as a set of quadric segments with certain parameters. Depending on the problem being solved, paraboloids [5, 6], ellipsoids [8, 9], or hyperboloids [7, 10] are used as quadrics. The number of segments is equal to the number of points in the required discrete distribution. The calculation of the quadric parameters is carried out using an iterative method, the convergence of which has been rigorously proved [6]. In Ref. [11], a modification of SQM for the calculation of the eikonal function for focusing into a set of points was proposed. This problem is important for the design of diffractive and Fresnel optical elements, in particular, beam shapers [12]. In addition, a refractive optical surface can be reconstructed from the eikonal function, which makes it possible to apply the method for the calculation of refractive optical elements [13].

Ray-mapping methods, based on the calculation of an equal-flux source-target mapping between the light source and the target region, constitute another class of widely used methods [8, 14–18]. The main problem here is to find an integrable mapping, which enables the calculation of a continuous optical surface. In many of these methods, the mapping is found either with use of the SQM [8] or by solving a Monge–Kantorovich mass transportation problem (MTP) with a quadratic cost function (CF) [15–18]. For a quadratic CF, the solution of an MTP (determination of the ray mapping) is reduced to the numerical solution of the Monge–Ampère equation. For the solution of this equation, the methods of Refs. [1–4], which are based on numerical solution of a system of nonlinear equations, and the explicit finite-difference method based on the replacement of the original equation by the non-stationary Monge–Ampère equation [15, 18, 19] are used. These MTP-based mapping methods can handle complex design problems including the design of double freeform surfaces for illuminance control [18]. At the same time, they lack sufficient theoretical justification. In fact, the integrability condition of the MTP-based mapping calculated for a quadratic CF is strictly fulfilled only for a very specific class of problems, e.g. for the problem of calculating the eikonal function providing focusing into a given region in paraxial approximation (see Section 5). The ‘correct’ CF providing the fulfillment of the integrability condition depends on the type of the problem being solved. The determination of the CFs for different inverse problems of non-imaging optics is of great theoretical and practical importance.

In the basic theoretical works [20, 21], it was shown that the problem of calculating a mirror generating a prescribed far-field intensity distribution can be formulated as a variational problem and as an MTP with a special CF. In this work, the approach of [20, 21] is for the first time applied to the solution of the inverse problem of calculating the eikonal function. We show that the problem of calculating the eikonal function providing focusing into a given region in general non-paraxial case can also be formulated as a variational problem and as an MTP with a non-quadratic CF. In the MTP, the masses are understood as the initial and target illuminance distributions, and the CF corresponds to the distance between the points of the initial and the target regions, respectively. This result demonstrates that the sought-for eikonal function corresponds to the mapping minimizing the total distance between the points of the initial and the target regions. In our opinion, this result is important for further development of the MTP-based mapping methods currently using paraxial quadratic CFs [15–18].

The proposed formalism makes it possible to reduce the calculation of the eikonal function to a linear programming problem (LPP). We also show that the calculation of the corresponding mapping can be further reduced to the solution of a linear assignment problem (LAP) [22]. As examples, we calculated the eikonal functions and designed refractive optical elements for focusing a circular beam into a rectangle and a beam with a square cross-section into a ring. The latter example demonstrates the efficiency of the proposed LAP-based approach: the computation time for square grids of 100 × 100 points on a standard PC is less than 1 minute. In contrast, the calculation of such a mapping using finite-difference methods for the solution of PDE of Monge–Ampère type requires solving a system of 10000 nonlinear equations and is by several orders of magnitude slower [16, 17]. Based on our experience with the SQM [9, 11], we also expect the SQM-based mapping computation [8] to be two orders of magnitude slower.

The paper is organized as follows. Section 2 is dedicated to the problem statement. In Section 3, the so-called weak (generalized) solution of the problem is defined. In Section 4, the corresponding variational problem is formulated. Section 5 describes the connection of the variational problem to the Monge–Kantorovich mass transportation problem. In Section 6, the MTP is reduced to a linear assignment problem. In our opinion, this section contains the most important practical result of this work. Section 7 is dedicated to the examples. For reader’s convenience, the mathematical proofs are given in the Appendices.

2. Problem statement

Let a light field with the illuminance distribution E₀(u) and the eikonal function Φ(u) be defined on a region G in the plane z = 0, where u = (u₁ u₂) are the Cartesian coordinates. We will refer to the region G as the aperture. Let us denote by $r (u) = (Φ_{u_{1}}^{'}, Φ_{u_{2}}^{'} \sqrt{1 - {(\nabla Φ)}^{2}})$ the unit vector of a the unit vector of a ray in the medium with unity refractive index, where $Φ_{u_{i}}^{'} = \partial Φ / \partial u_{i}$ , i = 1, 2, $\nabla Φ = (Φ_{u_{1}}^{'}, Φ_{u_{2}}^{'})$ . Let us consider a “ray mapping” γ_Φ of the region G into a region in the plane z = f > 0 defined by the eikonal function Φ(u). In this mapping, the point x = (x₁, x₂) = γ_Φ(u) is the point of intersection of the ray and the plane z = f (Fig. 1):

x = u + \nabla Φ f / \sqrt{1 - {(\nabla Φ)}^{2}} .

Let us note that Eq. (1) is written in a general non-paraxial form. The illuminance distribution in the target region can be found from the energy conservation law along the ray tubes:

E (x) = E (γ_{Φ} (u)) = E_{0} (u) / J (γ_{Φ} (u)),

where J(γ_Φ(u)) = ∂x₁/∂u₁ · ∂x₂/∂u₂ − ∂x₂/∂u₁ · ∂x₁/∂u₂ is the Jacobian of the coordinate transformation defined by Eq. (1). Equation (2a) can be written in the integral form

\int_{ω} E_{0} (u) d u = \int_{γ_{Φ} (ω)} E (x) d x,

where ω ⊆ G is an arbitrary region. The integrals with respect to du (dx) in Eq. (2b) denote double integrals with respect to du₁du₂ (dx₁dx₂). The energy conservation law can also be represented in the following general form (see Lemma 3.3 in [20]):

\int_{G} h (γ_{Φ} (u)) E_{0} (u) d u = \int_{D} h (x) E (x) d x,

where h(x) is an arbitrary continuous function.

Fig. 1 Geometry of the problem of eikonal function calculation.

Download Full Size | PDF

Let us formulate the inverse problem of calculating the eikonal function of the light field from the condition of generating a prescribed illuminance distribution: it is required to calculate the eikonal function Φ(u), u ∈ G from the condition of generating a given illuminance distribution E(x), x ∈ D in a region D in the plane z = f (Fig. 1).

3. Weak solution

In the papers [20, 21] dedicated to the design of mirrors generating prescribed intensity distributions in the far field, the concept of the so-called weak (generalized) solution was introduced. In this section, we present a similar concept for the problem of eikonal function calculation.

Substituting Eq. (1) into Eq. (2a), one can see that the calculation of the function Φ(u) is reduced to finding a solution to a nonlinear elliptic PDE. The requirement of differentiability of the eikonal function Φ(u) and the mapping γ_Φ(u) significantly limits the class of problems, in which this approach is applicable (i.e. the class of generated illuminance distributions). In this regard, let us introduce the concept of the weak (generalized) solution that makes it possible to extend the class of admissible functions Φ(u). To define the weak solution, let us consider the construction of the eikonal function Φ(u) from eikonal functions of lenses with the foci at the points x ∈ D [11]. The eikonal function of a lens has the form

Φ_{l} (u; x) = Ψ (x) - ρ (u, x)

where

ρ (u, x) = \sqrt{f^{2} + {(u - x)}^{2}}

is the distance between the points with coordinates (u₁, u₂, 0) and (x₁, x₂, f). It is easy to see that Eq. (3) at z = 0 corresponds to the eikonal function of a converging spherical beam with the focus at the point x ∈ D. The constant Ψ(x) in Eq. (3) corresponds to the eikonal value at the focus. The eikonal function for focusing into a prescribed region D can be represented as the envelope of a family of functions of Eq. (3) with respect to the parameters x = (x₁, x₂) ∈ D [14]. Indeed, by definition, the envelope surface is tangent to each of the eikonal functions of the family defined by Eq. (3) at a certain point. At the point of tangency, the partial derivatives of the envelope surface are equal to the partial derivatives of the lens eikonal. Since the direction of the ray is determined by these partial derivatives, the rays corresponding to the eikonal in the form of the envelope surface will be directed to the points of the region D, which were chosen as the foci of the lenses in Eq. (3). The envelope surface of the family of lens eikonals defined by Eq. (3) has the form [11]

{\begin{cases} Φ_{l} (u; x) = Ψ (x) - ρ (u, x), \\ \frac{\partial}{\partial x_{1}} [Ψ (x) - ρ (u, x)] = 0, \\ \frac{\partial}{\partial x_{2}} [Ψ (x) - ρ (u, x)] = 0 . \end{cases}

The system (5) contains the lens eikonal (first equation) and the derivatives of this eikonal with respect to the parameters (x₁, x₂) (second and third equations). Let us note that the last two equations in the system (5) describe the necessary condition for an extremum of the function (3) with respect to the variables (x₁, x₂) at a fixed value of u = (u₁, u₂). This condition corresponds to the Fermat’s principle and allows one to represent the eikonal function with respect to the variables (u₁, u₂) in one of the following forms:

Φ (u) = \max_{x \in D} [Ψ (x) - ρ (u, x)],

Φ (u) = \min_{x \in D} [Ψ (x) - ρ (u, x)] .

The representations in Eqs. (6) depend on the function Ψ(x), which defines the energy distribution in the focusing region. By construction, the function Ψ(x) corresponds to the eikonal in the region D, so it can also be represented as an envelope of a family of eikonals of lenses defined on the region D and focusing light into the points u ∊ G. The eikonal function of such a lens has the following form:

Ψ_{l} (x; u) = Φ (u) + ρ (u, x),

where the constant Φ(u) corresponds to the eikonal value at the focal point. The envelope of the family defined by Eq. (7) with respect to the parameters u = (u₁, u₂) has the form

{\begin{cases} Ψ_{l} (x; u) = Φ (u) + ρ (u, x), \\ \frac{\partial}{\partial u_{1}} [Φ (u) + ρ (u, x)] = 0, \\ \frac{\partial}{\partial u_{2}} [Φ (u) + ρ (u, x)] = 0 . \end{cases}

The last two equations of the system (8) describe the necessary condition for an extremum of function (7) with respect to the variables (u₁, u₂) at a fixed value of x = (x₁, x₂). This enables the representation of the function Ψ(x) in one of the following forms:

Ψ (x) = \max_{u \in G} [Φ (u) + ρ (u, x)],

Ψ (x) = \min_{u \in G} [Φ (u) + ρ (u, x)] .

Let us now show that if the direct mapping γ_Φ : G → D is defined by Eq. (6a), then the inverse mapping

γ_{Φ}^{- 1} : D \to G

is defined by Eq. (9b). To do this, we additionally assume that the mapping γ_Φ is a bijection. Indeed, let x₀ = γ_Φ (u₀), where x₀ = arg max_x_∈_D [Ψ(x) − ρ (u₀, x)]. Let us prove by contradiction that in this case u₀ = γ_Ψ(x₀). Assume that u₁ = γ_Ψ(x₀), where u₁ ≠ u₀. In this case, according to Eq. (9b), we obtain the inequalities Ψ(x₀) = Φ(u₁) + ρ(u₁, x₀) < Φ(u₀) + ρ(u₀, x₀) and Φ(u₀) > Ψ(x₀) − ρ(u₀, x₀). The latter inequality contradicts the initial condition x₀ = γ_Φ(u₀). according to which Φ(u₀) = Ψ(x₀) − ρ(u₀, x₀). Thus,

γ_{Ψ} = γ_{Φ}^{- 1}

. Similarly, one can show that if the direct mapping γ_Φ : G → D is defined by Eq. (6b), the inverse mapping

γ_{Φ}^{- 1} : D \to G

is defined by Eq. (9a). In what follows, we consider the case when the eikonal functions at the aperture G and the target region D are related by the following equations:

Φ (u) = \max_{x \in D} [Ψ (x) - ρ (u, x)], Ψ (x) = \min_{u \in G} [Φ (u) + ρ (u, x)] .

Both extrema in Eq. (10) are reached when u and x satisfy the equality x = γ_Φ(u)

(i . e . u = γ_{Φ}^{- 1} (x) = γ_{Ψ} (x))

. From the second equation of (10), it follows that the eikonal function Φ(u) is the “upper” envelope, since for every lens eikonal function Φ_l (u; x) from the family defined by Eq. (5), the inequality Φ(u) ≥ Φ_l (u; x) holds.

Let us note that the functions Φ(u) and Ψ(x) are related by the inequality Ψ(x) − Φ(u) ≤ ρ(u,x), which turns into an equality if u and x satisfy the equality γ_Φ(u) = x. Taking this into account, we can rewrite the mapping γ_Φ in the following form:

γ_{Φ} (u) = {x \in D | Ψ (x) - Φ (u) = ρ (u, x)} .

Let us generalize the concept of the solution of the considered problem without assuming the differentiability of the eikonal functions Φ(u) and Ψ(x).

Definition 1. Let us call a pair of functions (α(u), β(x)) a “proper pair” if α(u) is continuous on G, β(x) is continuous on D, and the functions α(u) and β(x) are related by the equalities

α (u) = \max_{x \in D} [β (x) - ρ (u, x)], β (x) = \min_{u \in G} [α (u) + ρ (u, x)] .

Let us note that the functions α(u) and β(x) do not have to be differentiable. By definition, the equality ∀_u_∈_G∀_x_∈_D β(x)−α(u) ≤ ρ(u, x) holds. It is shown in Appendix A that the functions α(u) and β(x) are Lipschitz functions in the regions G and D, respectively. According to Rademacher’s theorem, a Lipschitz function is differentiable almost everywhere in the corresponding region that is, the points at which the function is not differentiable form a set of Lebesgue measure zero). Thus, the functions α(u) and β(x) forming a proper pair are differentiable almost everywhere and the corresponding mapping γ_Φ [see Eq. (11)] is single-valued at the points of differentiability.

Finally, let us define the term “weak (generalized) solution”.

Definition 2. A weak solution is a proper pair of functions (Φ(u), Ψ(x), which satisfy the energy conservation law of Eq. (2b) for the mapping defined by the following expression:

γ_{Φ} (u) = {x \in D | Ψ (x) - Φ (u) = ρ (u, x)} .

Note that in the weak solution the derivatives of the function Φ(u) and, consequently, the mapping (1) are defined almost everywhere. Since the set of points where the derivatives do not exist has measure zero, these points do not affect the values of the integrals in Eq. (2b).

4. Variational formulation of the problem

The problem of finding a weak solution defined in the previous section can be formulated as a variational problem. Let us consider the set of pairs of continuous functions

Ω = {(α, β) \in C (G) \times C (D) | \forall_{u \in G}, \forall_{x \in D} β (x) - α (u) \leq ρ (u, x)} .

Let us define the following functional on Ω:

I (α, β) = \int_{D} β (x) E (x) d x - \int_{G} α (u) E_{0} (u) d u .

It is obvious that I(α, β) is continuous with respect to the norm ║(α, β)║_Ω = max {║α║_∞, ║β║_∞}, where the norm ║·║_∞ is defined as the maximum of the absolute value of the function over the corresponding region.

Theorem 1. The following conditions are equivalent:

The pair (Φ, Ψ) is a weak solution;
The pair (Φ, Ψ) maximizes the functional of Eq. (15).

The proof of Theorem 1 is completely analogous to the proof of Theorem 3.4 in the paper [20], where the variational formulation of the problem of calculating a mirror for generating a prescribed intensity distribution is considered. The only difference consists in the form of the function ρ(u, x), which does not alter the proof.

In discrete form, the problem of maximization of the functional of Eq. (15) can be represented as a linear programming problem (LPP). Indeed, let us consider two sets of points ${u_{i} \in G}_{i = \bar{1, N}}$ and ${x_{i} \in D}_{j = \bar{1, M}}$ , which are nodes of rectangular grids Λ and Θ on G and D. Let us denote by $e_{i}^{0}$ and e_j the average values of the functions E₀(u) and E(x) in the cells of the grids. Besides, we denote ρ_ij = ρ(u_i, x_j). The values $e_{i}^{0}$ , e_j, and ρ_ij are known, and α_i = α(u_i), β_j = β(x_j) are the unknown variables. In discrete form, the problem of maximization of the functional of Eq. (15) on the set Ω is reduced to the following LPP:

\begin{array}{l} I^{Λ, Θ} (α_{1}, \dots, α_{N}; β_{1}, \dots, β_{M}) = \sum_{j = 1}^{M} β_{j} e_{j} | τ_{j} | - \sum_{i = 1}^{N} α_{i} e_{i}^{0} | τ_{i} | \to \max, \\ β_{j} - α_{i} \leq ρ_{i j}, i = 1, \dots, N j = 1, \dots, M, \end{array}

where |τ_j| and |τ_i| are the areas of the corresponding grid cells. Let us note that the function I^Λ,Θ in Eq. (16) corresponds to the Riemann sum of the continuous functional of Eq. (15). It is natural to expect that when the grid step in the discrete problem (16) is decreased, its solution converges to the solution of the original continuous problem [21]. Note that the constraint matrix of the LPP of Eq. (16) has a special form: the matrix elements take the values 0 or ±1. Such LPPs can be solved using a polynomial complexity algorithm [23].

5. Connection with the mass transportation problem

Let us consider the connection of the variational problem of maximizing the functional of Eq. (15) with a mass transportation problem (MTP) from the region G to the region D with the cost function ρ(u, x). The mass distributions in the regions are described by the functions E₀(u), u ∈ G and E(x), x ∈ D. We call a transport map the mapping P : G → D, which conserves the energy (measure), i.e. for any measurable subset the following equality holds:

\int_{ω} E_{0} (u) d u = \int_{P (ω)} E (x) d x .

The mapping P(u) can be defined almost everywhere. Let us introduce the following functional on the set of transport maps:

C (P) = \int_{G} ρ (u, P (u)) E_{0} (u) d u .

Theorem 2. Let (Φ, Ψ) be a weak solution of the problem. Then the mapping γ_Φ minimizes the functional of Eq. (18). Furthermore, if P is another optimal map, it coincides with γ_Φ almost everywhere. The proof of this theorem is given in Appendix B.

Theorem 2 yields an interesting result, which can be interpreted as a generalization of the Fermat’s principle. Indeed, among all mappings G → D which conserve the energy, the weak solution of the inverse problem of calculating the eikonal function Φ(u) corresponds to a mapping, for which the total distance between the points of the regions G and D with the weight function E₀(u) is minimized.

Let us consider the calculation of the eikonal function Φ(u) from the mapping x = γ_Φ(u). For the implementation of the mapping γ_Φ(u), partial derivatives of the function Φ(u), which determine the directions of the rays, must have the form

\nabla Φ = (Φ_{u_{1}}^{'}, Φ_{u_{2}}^{'}) = (x - u) / ρ (u, x) .

The eikonal function can be recovered from its partial derivatives using the following well-known formula:

Φ (u_{1}, u_{2}) = Φ (0, 0) + \int_{0}^{u_{1}} Φ_{u_{1}}^{'} (t, 0) d t + \int_{0}^{u_{2}} Φ_{u_{2}}^{'} (u_{1}, s) d s .

Indeed, according to the definition of the weak solution (see Section 2) the function Φ(u), u ∈ G is a Lipschitz function. It is easy to show that the function Φ(u) will also be a Lipschitz function along any curve u(τ) with natural parametrization. Thus, the function F(τ) = Φ (u(τ)) is an absolutely continuous function of a natural parameter τ. Due to the properties of absolutely continuous functions, the function F(τ) can be uniquely recovered from its derivative, which exists almost everywhere, by integrating it along the curve. In a special case when the curve u(τ) is a polygonal chain with the segments parallel to the coordinate axes, we arrive to Eq. (20).

Let us consider the paraxial approximation, when |x − u|/f ≪ 1 (|∇Φ| ≪ 1). In this case ρ(u, x) ≈ f + (x − u)²/2f, and the problem of finding the mapping x = γ_Φ(u) is reduced to the minimization of the functional

C (γ_{Φ} (u)) = \int_{G} {[γ_{G} (u) - u]}^{2} E_{0} (u) d u \underset{γ_{Φ} (u)}{\to} \min

under the condition of Eq. (2a) or Eq. (2b), corresponding to the energy conservation law. Equations (21) and (2a) [or (2b)] define an MTP with a quadratic cost function. For such a cost function, the solution of the MTPis unique and has the form [19]

x (u) = γ_{Φ} (u) = \nabla Φ_{add} (u),

where Φ_add(u) is a certain convex function. For a smooth mapping x(u), the function Φ_add(u) satisfies the Monge–Ampère equation [19]:

\det [\nabla^{2} Φ_{add} (u)] = E_{0} (u) / E (\nabla Φ_{add} (u)),

where ∇²Φ_add(u) is the Hessian matrix. Let us show that the function Φ_add(u) provides the generation of a prescribed illuminance distribution in paraxial approximation. Under this approximation, Eq. (19) takes the form

\nabla Φ = [x (u) - u] / f .

Let us represent the eikonal function in the form Φ(u) = −u²/(2f) + Φ_add(u)/f−f/2. Using this representation, we obtain from Eq. (24):

x (u) = \nabla Φ_{add} (u) .

In this case, Eq. (2a) becomes Eq. (23). In addition, the function Φ_add(u) is convex. Indeed, in paraxial approximation one can obtain from Eq. (6a) that

Φ_{add} (u) = \max_{x \in D} [x \cdot u - Ψ_{add} (x)],

where Ψ_add(x) = x²/2 − fΨ(x) + f²/2. According to Eq. (26), the function Φ_add(u) is the Legendre transform of the function Ψ_add(x). Since the Legendre transform is a convex function, the resulting function Φ_add(u) is convex. Moreover, in paraxial approximation one can easily obtain from Eq. (9a) that the function Ψ_add(x) is the Legendre transform of the function Φ_add(u) = u²/2 + fΦ(u) + f²/2 and, consequently, Φ_add(x) is also convex.

Thus, we have shown that, in paraxial approximation, Theorem 2 yields a known result for the mass transportation problem with a quadratic cost function [Eqs. (21)–(23)].

6. Reduction to the linear assignment problem

In the discrete case, the mass transportation problem can be formulated as a linear assignment problem (LAP) [22]. Indeed, assume that the regions G and D are divided into N cells each (or approximated by N cells) with equal energies, i.e. for any pair of cells g_i ⊂ G, d_j ⊂ D the equality $\int_{g_{i}} E_{0} (u) d u = \int_{d_{j}} E (x) d x = e$ holds. In this case, all mappings G → D satisfying the energy conservation law can be described as permutations of the set of N numbers (i₁,…, i_N), which determine to which cells d_i of the region D the cells g₁,…, g_N are mapped. In terms of the LAP, the mapping of the cell g_i to the cell d_j can be interpreted as the assignment of the j-th task to the i-th agent. According to Eq. (18), the cost of the tasks is described by the matrix C_i_,_j = ρ (u_i, x_j), i, j = 1, …, N, where u_i, x_j are the centers of the cells g_i and d_j, and ρ(u, x) is defined by Eq. (4). Thus, for the construction of the mapping γ_Φ corresponding to the weak solution, the following assignment problem can be solved:

C (j_{1}, \dots, j_{N}) = \sum_{i = 1}^{N} ρ (u_{i}, x_{j_{i}}) \to \min,

where the minimum is sought over all permutations (j₁,…, j_N) For the solution of the LAP, efficient polynomial algorithms have been proposed [22, 24, 25]. Substituting the obtained mapping γ_Φ into Eqs. (19) and (20) and integrating them, we arrive at the eikonal function.

7. Design examples

7.1. Focusing a circular beam into a rectangle: comparison of LPP- and LAP-based approaches

As a first example, let us consider the calculation of the eikonal function from the condition of focusing a circular uniform-intensity beam (E₀(u) = E₀, u ∈ G) into a uniformly illuminated rectangle (E(x) = E, x ∈ D). This example is purely illustrative and serves to demonstrate the different computational complexity of the two considered approaches based on the solution of an LPP and an LAP, respectively. In the calculations, we used the following parameters: beam radius (the radius of the region G) R = 3 mm, focal distance (distance to the target plane) f = 50 mm, the size of the rectangle (the size of the region D) $d_{x}_{_{1}} = 12 mm$ , $d_{x}_{_{2}} = 4 mm$ (the center of the rectangle is located at the origin). Due to the symmetry of the problem, it is sufficient to calculate the eikonal function in the first quadrant of the coordinate plane (in the region $G_{1} = {u | u_{1} > 0, u_{2} > 0, u_{1}^{2} + u_{2}^{2} \leq R^{2}}$ ) from the condition of focusing into the rectangle D₁ with the size of 6 × 2 mm², which is a quarter of the original rectangle D located in the first quadrant. In order to formulate the LPP of Eq. (16), let us impose rectangular grids over G₁ and D₁ with the discretization steps Δu₁ = Δu₂ = R/N_u = 3/13 ≈ 0.231 mm and $Δ x_{1} = d_{x_{1}} / N_{x_{1}} = 6 / 16 = 0.375 mm$ , $Δ x_{2} = d_{x_{2}} / N_{x_{2}} = 2 / 9 \approx 0.222 mm$ , respectively, and take the nodes of the grids that lie inside these regions. In this case, the function Φ(u) in the region G₁ and the function Ψ(x) in the rectangle D₁ at these nodes are both represented by sets of N = M = 144 values. The obtained values of N and M are consistent with the performance of the PC used for the solution of the constructed LPP containing 288 variables and N² = 20736 inequality constraints. For the solution of this LPP, the standard MATLAB^® routine linprog was used. The time of solving the LPP on a standard modern PC (Intel^® Core™ i7-2600 CPU @3.4 GHz) was of about 5 hours. The found eikonal function Φ(u) is shown in Fig. 2.

Fig. 2 The eikonal function in the first quadrant calculated using an LPP.

Download Full Size | PDF

Using the eikonal function Φ(u), one can reconstruct a refractive optical element located immediately before the plane z = 0, which generates this eikonal distribution in the region G. Let us consider a plane illuminating beam propagating in the positive direction of the z axis. We assume that the first surface of the optical element is the plane perpendicular to the propagation direction of the incident beam. The second surface can be calculated from the condition of generating the given eikonal in the region G of the plane z = 0 using the Fermat’s principle [13]. At the same time, the considered example has “paraxial” parameters, which makes it possible to calculate the height of the refracting surface in the approximation of a thin optical element [12]:

z (u) = Φ (u) / (n - 1),

where n is the refractive index of the material of the optical element.

In order to verify the theoretical results obtained in the present work, the eikonal function Φ(u) shown in Fig. 2 was used in Eq. (28) for the calculation of a refracting surface. Then, the designed optical element was simulated in the commercial software TracePro^® implementing a ray tracing technique [26]. For this, the surface defined by Eq. (28) was approximated by non-uniform rational basis splines (NURBS). Figure 3(a) shows the designed refractive optical element with the refractive index n = 1.5 (inset) and the simulated illuminance distribution which is generated upon illumination of the optical element by a collimated uniform-intensity beam. For all examples in the present work, 500 000 rays were traced in the illuminance distribution simulations. Figure 3(b) confirms the focusing into the rectangle of the specified size. The normalized root-mean-square deviation (NRMSD) of the obtained illuminance distribution from the required uniform distribution amounts to 16.2%. The nonuniformity of the distribution is due to the insufficient number of points N = 144 used for the representation of the eikonal function (Fig. 2). At the same time, the increase in the number of points N leads to a substantial increase in the computational complexity of the LPP, since the size of the constraint matrix determining the complexity equals to N². As mentioned above, at N = 144 the time of solving the corresponding LPP with N² = 20736 inequality constraints is of about 5 hours. Let us note that at N > 150 the MATLAB^® routine linprog produces an internal error due to the high dimensionality of the problem.

Fig. 3 (insets) Optical elements for the focusing of a circular beam with 3 mm radius into a rectangle of the size 12 × 4 mm² calculated using an LPP (a) and an LAP (c). Generated illuminance distributions (a) and (c) and its cross sections along the coordinate axes (b) and (d), respectively.

Download Full Size | PDF

A significantly more efficient approach to the solution of the problem consists in the calculation of the mapping x = γ_Φ(u) using the linear assignment problem (LAP) of Eq. (27). For the calculation of the mapping γ_Φ, the aperture G was approximated by a set of N = 1060 square cells g_i with the size 0.162 × 0.162 mm². The rectangle D was also represented as 1060 rectangular cells d_i with the size 0.231 × 0.21 mm². For the solution of the LAP, the auction algorithm [24] implemented in [27] was used. For the considered example, the time of solving the LAP was less than a second. Using Eq. (19), the derivatives of the eikonal function were calculated from the found mapping, and then the eikonal was reconstructed by numerical integration of Eq. (20). Similarly, to the previous case, the surface of the optical element focusing into the required rectangle was calculated using Eq. (28). The obtained optical element and the generated illuminance distribution calculated in TracePro^® are shown in Fig. 3(c). In contrast to Figs. 3(a) and 3(b), the illuminance distribution in Figs. 3(c) and 3(d) is significantly more uniform, and the NRMSD of the obtained distribution from the uniform distribution amounts to 5.6% only. The presented examples demonstrate that the calculation of the eikonal function using the LAP of Eq. (27) is optimal in terms of the computation time.

7.2. Focusing a square beam into a ring: LAP-based approach

The utilization of the proposed approach for the eikonal calculation on the basis of the solution of an LAP is of great interest in the cases when the methods based on the direct solution of the equations of the Monge–Ampère type [1–4, 15, 18] are too complicated to be applied. In these methods, one of the stages of the solution consists in the formulation of the boundary conditions. In the case of an arbitrary shape of the boundaries of the optical element and the illuminated region, the formulation of these conditions is a non-trivial problem. For instance, the numerical method for the solution of the Monge–Ampère equation [19] used in [15] 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 in paraxial approximation a sophisticated illuminance distribution, e.g. corresponding to a logo or an image of a human face [4]. These images are generated in a square region on some substantially non-zero background. At the same time, this approach has limited applicability for generating non-square illuminated regions (triangle, circle, ring) on a zero background. For the application of the method in this case, it is necessary to define some initial differentiable mapping of the boundary of a square to the boundary of the focusing region, and to obtain the corresponding solution of the Monge–Ampère equation. This problem is non-trivial even for relatively simple regions such as a triangle or a ring. In particular, a square is simply-connected (1-connected) while a ring is not. Due to this reason, the formulation of the boundary conditions describing the mapping of the boundary of a square into the boundary of a ring is a difficult problem. Moreover, the boundary of a square is not smooth, and, therefore, the construction of a differentiable mapping is impossible.

In this connection, as a second example we consider the application of the eikonal calculation method based on the solution of an LAP to the design of an optical element focusing an incident uniform-intensity beam with a square section into a uniformly illuminated ring. The calculation was carried out for the following parameters: the size of the incident beam (the size of the aperture G) is 1 × 1 mm², the distance to the focusing plane f = 5 mm, the inner and outer radii of the ring (of the region D) R₁ = 1 mm, R₂ = 2.5 mm. Let us note that these parameters correspond to a non-paraxial case. Indeed, the maximum value of the angle between the z axis and the rays connecting the boundary of the region G and the outer boundary of the region D exceeds 21°.

For calculating the mapping γ_Φ, the square aperture G and the ring D were represented by sets of N = 100² = 10000 square cells with the dimensions 0.01 × 0.01 mm² and 0.04 × 0.04 mm², respectively. In this case, the time of solving the LAP using the auction algorithm [27] was less than one minute. The obtained mapping is shown in Fig. 4. For comparison, let us note that the calculation of this mapping using finite-difference methods for solving PDE of Monge–Ampère type [1–4, 16, 17] requires the solution of a system of 10000 nonlinear equations and is by several orders of magnitude slower [16].

Fig. 4 Rectangular grid on the aperture (a) and its mapping to the focusing plane (b).

Download Full Size | PDF

Using the found mapping, the eikonal function and the surface of the optical element were reconstructed using Eqs. (20) and (28), respectively. The designed element is shown in the inset of Fig. 5(a). Figures 5(a) and 5(b) show the generated illuminance distribution calculated using the TracePro^® software, and its cross sections along the coordinate axes, respectively. Figures 5(a) and 5(b) demonstrate generation of a highly uniform illuminance distribution in the required ring. The NRMSD of the obtained distribution from the uniform distribution amounts to 5.8%. It is interesting to note that the surface of the optical element is non-differentiable at u = 0. This is due to the fact that the neighborhood of this point is mapped to a neighborhood of the inner boundary of the ring D. Indeed, the focusing into a thin ring is described by the eikonal function $Φ (| u |) ~ - \sqrt{f^{2} + {(R_{1} - | u |)}^{2} + C}$ [28], which is non-differentiable at u = 0.

Fig. 5 (inset) Optical element for the focusing of a beam with a square section (1 × 1 mm²) into a ring (inner radius 1 mm, outer radius 2.5 mm) calculated by solving an LAP. Generated illuminance distribution (a) and its cross sections along the coordinate axes (b).

Download Full Size | PDF

8. Conclusion

In the present work, it was shown that the problem of calculation of the eikonal function from the condition of focusing into a given region can be formulated as a variational problem of Eq. (15) and as a Monge–Kantorovich mass transportation problem. It was obtained that the cost function in the MTP corresponds to the distance between a point in the region in which the eikonal function is defined and a point in the focusing region. Such a cost function yields a mapping, which minimizes the total distance between the points of the initial plane and the focusing region. The proposed formalism makes it possible to reduce the calculation of the eikonal function to finding a solution to an LPP of Eq. (16). In this case, the calculation of the corresponding mapping is reduced to the solution of an LAP (27).

The presented design examples demonstrate that the eikonal function calculation using the LAP of Eq. (27) is significantly less computationally expensive than the solution of an LPP of Eq. (16). In addition, in contrast to the methods based on the numerical solution of the Monge–Ampère equation, the proposed LAP-based method can be applied for illuminating regions with a complex boundary on a zero background. Let us note that this approach can also be used for the design of mirrors generating prescribed illuminance distributions. The only difference consists in the used cost function [20, 21]. In terms of computational complexity, the LAP-based approach is expected to be significantly more efficient than the approach based on the solution of an LPP considered in [29, 30].

Appendix A. Proof of the Lemma

Lemma. Let (α(u), β(x)) be a proper pair of functions. Then the function α(u) [the function β(x)] is a Lipschitz function on G (on D), i.e. there exists a constant L, for which |α(u″) −α(u′)| < L · |u″ − u′| for any u″, u′ ∈ G.

Proof. Let us choose arbitrary u″, u′ ∈ G. Let us assume that α(u″) ≥ α(u′) [the case α(u″) < α(u′) can be treated similarly]. Let x″ ∈ D be a point, for which the equality β(x″) − α(u″) = ρ(u″, x″) holds. Then

α (u^{″}) - α (u^{'}) = β (x^{″}) - ρ (u^{″}, x^{″}) - α (u^{'}) \leq ρ (u^{'}, x^{″}) - ρ (u^{″}, x^{″}) .

According to the mean value theorem, the right-hand side of the last inequality can be written as the following scalar product:

ρ (u^{'}, x^{″}) - ρ (u^{″}, x^{″}) = \nabla ρ (\tilde{u}, x^{″}) \cdot (u^{'} - x^{″}),

where

\tilde{u}

∈ G is a certain point from a segment connecting the points u′ and u″, and the gradient of the function ρ(u, x) is calculated with respect to the variables (u₁, u₂) Since the left-hand side of the inequality (29) is non-negative, we can estimate its absolute value using Eq. (30) and the Cauchy-Schwarz inequality as

| α (u^{″}) - α (u^{'}) | \leq | \nabla ρ (\tilde{u}, x^{″}) | \cdot | u^{'} - u^{″} | \leq L \cdot | u^{'}, u^{″} |,

where L = max_u_∈_G_,_x_∈_D |∇ρ(u, x)|. Thus, the Lipschitz condition is fulfilled with the constant L. The proof for the function β(x) is completely analogous.

Appendix B. Proof of Theorem 2

Let P be an arbitrary transport map. According to Eq. (9b), the inequality ρ(u, P(u)) ≥ Ψ (P(u)) − Φ(u) holds, and it turns into equality if and only if P(u) = γ_Φ(u). Let us multiply this inequality by E₀(u) and integrate it with respect to du:

C (P) = \int_{G} ρ (u, P (u)) E_{0} (u) d u \geq \int_{G} Ψ (P (u)) E_{0} (u) d u - \int_{G} Φ (u) E_{0} (u) d u .

According to Eq. (2c), we have:

\int_{G} Ψ (P (u)) E_{0} (u) d u = \int_{D} Ψ (x) E (x) d x = \int_{G} Ψ (γ_{Φ} (u)) E_{0} (u) d u .

Therefore,

\begin{array}{l} \int_{G} Ψ (P (u)) E_{0} (u) d u - \int_{G} Φ (u) E_{0} (u) d u \\ = \int_{G} Ψ (γ_{^{Φ}} (u)) E_{0} (u) d u - \int_{G} Φ (u) E_{0} (u) d u \\ = \int_{G} [Ψ (γ_{^{Φ}} (u)) - Φ (u)] E_{0} (u) d u \\ = \int_{G} ρ (u, γ_{Φ} (u)) E_{0} (u) d u \\ = C (γ_{Φ}) . \end{array}

By comparing Eqs. (31) and (32), we obtain that C(P) ≥ C(γ_Φ), which proves the first statement of the theorem. Furthermore, if P(u) ≠ γ_Φ(u) on a set of nonzero measure, this inequality becomes strict.

Funding

Ministry of Education and Science of Russian Federation; Russian Foundation for Basic Research (RFBR) grant 17-47-630164.

References and links

1. R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “Freeform illumination design: a nonlinear boundary problem for the Monge–Ampère equation,” Opt. Lett. 38(2), 229–231 (2013). [CrossRef] [PubMed]

2. R. Wu, P. Benítez, Y. Zhang, and J. C. Miñano, “Influence of the characteristics of a light source and target on the Monge–Ampère equation method in freeform optics design,” Opt. Lett. 39(3), 634–637 (2014). [CrossRef] [PubMed]

3. 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]

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

5. V. I. Oliker, “Mathematical aspects of design of beam shaping surfaces in geometrical optics,” in Trends in Nonlinear Analysis, M. Kirkilionis, S. Krömker, R. Rannacher, and F. Tomi, eds. (Springer, 2003). [CrossRef]

6. S. A. Kochengin and V. I. Oliker, “Computational algorithms for constructing reflectors,” Comput. Vis. Sci. 6(1), 15–21 (2003). [CrossRef]

7. V. I. Oliker, “Controlling light with freeform multifocal lens designed with supporting quadric method (SQM),” Opt. Express 25(4), A58–A72 (2017). [CrossRef] [PubMed]

8. F. R. Fournier, W. J. Cassarly, and J. P. Rolland, “Fast freeform reflector generation using source-target maps,” Opt. Express 18(5), 5295–5304 (2010). [CrossRef] [PubMed]

9. L. L. Doskolovich, K. V. Borisova, M. A. Moiseev, and N. L. Kazanskiy, “Design of mirrors for generating prescribed continuous illuminance distributions on the basis of the supporting quadric method,” Appl. Opt. 55(4), 687–695 (2016). [CrossRef] [PubMed]

10. D. Michaelis, P. Schreiber, and A. Bräuer, “Cartesian oval representation of freeform optics in illumination systems,” Opt. Lett. 36(6), 918–920 (2011). [CrossRef] [PubMed]

11. L. L. Doskolovich, M. A. Moiseev, E. A. Bezus, and V. Oliker, “On the use of the supporting quadric method in the problem of the light field eikonal calculation,” Opt. Express 23(15), 19605–19617 (2015). [CrossRef] [PubMed]

12. V. A. Soifer, V. V. Kotlyar, and L. L. Doskolovich, Iterative Methods for Diffractive Optical Elements Computation (Taylor & Francis Ltd., 1997), p. 245.

13. L. L. Doskolovich, A. Y. Dmitriev, and S. I. Kharitonov, “Analytic design of optical elements generating a line focus,” Opt. Eng. 52(9), 091707 (2013). [CrossRef]

14. 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]

15. 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]

16. 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]

17. R. Wu, Y. Zhang, M. M. Sulman, Z. Zheng, P. Benítez, and J. C. Miñano, “Initial design with L² Monge–Kantorovich theory for the Monge–Ampère equation method in freeform surface illumination design,” Opt. Express 22(13), 16161–16177 (2014). [CrossRef] [PubMed]

18. A. Bäuerle, A. Bruneton, R. Wester, J. Stollenwerk, and P. Loosen, “Algorithm for irradiance tailoring using multiple freeform optical surfaces,” Opt. Express 20(13), 14477–14485 (2012). [CrossRef] [PubMed]

19. 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]

20. 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]

21. X.-J. Wang, “On the design of a reflector antenna II,” Calc. Var. 20(3), 329–341 (2004). [CrossRef]

22. J. Munkres, “Algorithms for the assignment and transportation problems,” SIAM J. Appl. Math. 5(1), 32–38 (1957). [CrossRef]

23. É. Tardos, “A strongly polynomial algorithm to solve combinatorial linear programs,” Oper. Res. 34(2), 250–256 (1986). [CrossRef]

24. D. P. Bertsekas, “The auction algorithm: A distributed relaxation method for the assignment problem,” Ann. Oper. Res. 14(1), 105–123 (1988). [CrossRef]

25. R. Jonker and A. Volgenant, “A shortest augmenting path algorithm for dense and sparse linear assignment problems,” Computing 38(4), 325–340 (1987). [CrossRef]

26. Opto-mechanical software TracePro. https://www.lambdares.com/tracepro/

27. Implementation of Bertsekas’ auction algorithm. http://www.mathworks.com/matlabcentral/fileexchange/48448

28. L. L. Doskolovich, S. N. Khonina, V. V. Kotlyar, I. V. Nikolsky, V. A. Soifer, and G. V. Uspleniev, “Focusators into a ring,” Opt. Quant. Electron. 25(11), 801–814 (1993). [CrossRef]

29. C. Canavesi, W. J. Cassarly, and J. P. Rolland, “Direct calculation algorithm for two-dimensional reflector design,” Opt. Lett. 37(18), 3852–3854 (2012). [CrossRef] [PubMed]

30. C. Canavesi, W. J. Cassarly, and J. P. Rolland, “Observations on the linear programming formulation of the single reflector design problem,” Opt. Express 20(4), 4050–4055 (2012). [CrossRef] [PubMed]

Variational approach to calculation of light field eikonal function for illuminating a prescribed region

Abstract

1. Introduction

2. Problem statement

3. Weak solution

4. Variational formulation of the problem

5. Connection with the mass transportation problem

6. Reduction to the linear assignment problem

7. Design examples

7.1. Focusing a circular beam into a rectangle: comparison of LPP- and LAP-based approaches

7.2. Focusing a square beam into a ring: LAP-based approach

8. Conclusion

Appendix A. Proof of the Lemma

Appendix B. Proof of Theorem 2

Funding

References and links

Cited By

Figures (5)

Equations (38)

Optics Express