Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions

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

doi:10.1364/OE.27.013083

1. Introduction

The problem of calculating a refractive or reflective optical surface from the condition of generating a prescribed irradiance distribution in a certain region belongs to the class of inverse problems of nonimaging optics and is extremely challenging. In the geometrical optics approximation, this problem is reduced to finding a solution to a fully nonlinear differential equation (NDE) of Monge–Ampère type. The solution of this NDE is a complex theoretical and computational problem.

In recent years, various methods for the design of optical elements have been proposed, based on numerical solution of NDE of Monge–Ampère type [1–7]. In [1–7], finite-difference methods are used for the solution of this NDE. 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 unknown function decribing the surface of the optical element. For the solution of the obtained system, iterative methods are used, e.g. Newton’s method. In a number of methods for designing optical elements [8–10], the solution of the Monge–Ampère equation is based on the method proposed in [11]. In this method, a steady-state solution of a parabolic Monge–Ampère equation is found using an explicit Euler scheme. In [12, 13], 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 refractive and reflective 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 nonimaging 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–13]. For example, a smooth optical surface cannot generate an irradiance distribution defined on a not connected 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.

A number of inverse problems of nonimaging 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 the 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 [14–19]. In particular, in the problems described by a standard Monge–Ampère equation, the cost function is quadratic [11]. Let us note that an MTP with a quadratic cost function corresponds to the problem of the calculation of an optical surface only in the paraxial approximation [9,18].

The formulation of an inverse problem of nonimaging optics as an MTP enables working with discontinuous ray mappings and with the corresponding piecewise-smooth continuous optical surfaces, in contrast to the smooth surfaces satisfying the NDEs of Monge–Ampère type [14–19]. Although piecewise-smooth optical surfaces are more challenging to manufacture, they allow one to generate much more complex irradiance distributions. In its discrete version, the MTP can be formulated as a linear assignment problem (LAP) [18, 19]. For the solution of LAPs, efficient polynomial algorithms have been proposed [20, 21]. Let us note that the approach consisting in the calculation of a ray mapping through the solution of an LAP was proposed in the recent works of the present authors [18,19]. In particular, in [18], an LAP-based approach for the calculation of the light field eikonal function providing the generation of a required irradiance distribution was proposed. In [19], we presented an iterative method, which uses the LAP-based approach and is intended for the design of a refractive surface generating a required irradiance distribution in the near field in the case of an incident beam with a plane wavefront. The design examples presented in [19] demonstrate the applicability of the LAP-based approach for generating irradiance distributions defined on complex-shaped and not connected regions.

In the present work, we for the first time consider the application of the LAP-based approach for the solution of the problem of calculating a refractive surface, which generates a prescribed irradiance distribution in the far field in the case of a point light source (e.g. in the case when the incident beam has a spherical wavefront). In nonimaging optics, the far field condition assumes that the dimensions of the refractive surface as compared with the distance to the target plane can be neglected. This condition is fulfilled for a wide class of LED lighting systems including searchlights, road, industrial, residential and architectural lighting systems.

The problem of generating a prescribed irradiance distribution in the far field, in contrast to a similar near-field problem [19], can be formulated as an MTP with a special, non-quadratic cost function [16,17]. As the masses in this MTP, the irradiance distribution formed by the light source on a unit sphere and the required irradiance distribution are used. Let us note that the formulation of the considered problem of calculating a refractive surface as an MTP [16,17] is known only to a very limited circle of researchers in optics. This is caused by “purely mathematical” character of the works [16,17]. Therefore, in this article, we give another proof of the equivalence of the problem of calculating a refractive surface and the corresponding MTP, which is oriented at the applied optics community. The proposed proof is based on the direct calculation of the extremum of the functional describing the transportation cost and requires only a basic knowledge of the calculus of variations, which is widely used in geometrical optics. A well-known example is the variational Fermat’s principle.

The formulation of the problem of calculating a refractive surface as an MTP enables applying an LAP-based approach. Within this approach, the calculation of an integrable ray mapping is reduced to the solution of an LAP. Then, the refractive surface is reconstructed from the obtained mapping. For the solution of this problem, a special method is proposed, which makes it possible to work with discontinuous mappings. The present paper also discusses the use of the LAP-based method for the design of optical elements with two “working” refractive surfaces. This problem is relevant for generating irradiance distributions with small angular dimensions as well as off-axis irradiance distributions [22]. As examples, we designed refractive optical elements with two working surfaces generating uniform irradiance distributions in a square and in a region in the form of the letters “AB” on a zero background. The results of the simulation of the calculated optical elements confirm high performance of the proposed method.

Let us emphasize that the mass transportation problems and LAP-based approaches considered in the previously published papers [18,19,23] are different from the results of the present work. In particular, the papers [18,19] consider the calculation of the eikonal function and refractive surfaces generating a required irradiance distribution in the near field in the case of a plane incident beam. The paper [23] considers the problem of calculating an optical element for collimated beam shaping, in which both the incident and the outgoing beams are assumed to have plane wavefronts. These problems are different from the problem of the current paper, where we consider the optical elements generating required far-field irradiance distribution from a spherical incident beam.

2. Formulation of the problem

Let us consider the problem of designing a refractive surface R generating a prescribed irradiance distribution in the far field for a point light source (Fig. 1). The refractive surface separates two media with the refractive indices n₁ and n₂, n₁ > n₂. The light source is located at the origin of coordinates in the medium with the refractive index n₁. Let us define the unit vector of a ray from the source as

\vec{e} (m) = (m_{1}, m_{2}, \sqrt{1 - m_{1}^{2} - m_{2}^{2}}),

where m = (m₁, m₂) are the Cartesian coordinates of the projection of this vector onto the plane z = 0 (Fig. 1). Here and in the rest of the paper, the 2-vectors are typeset in bold, while the 3-vectors are denoted by an arrow. We assume that the light source is described by the intensity function I(m), m ∈ Ω. The intensity is defined as the light flux emitted by the source into an element of the solid angle.

Fig. 1 Geometry of the problem.

Download Full Size | PDF

Let us define the refractive surface in the form

\vec{R} (m) = ρ (m) \vec{e} (m),

where ρ(m) is the length of the radius vector of a point on the surface. The problem being solved consists in the calculation of a refractive surface R⃗(m), which generates a prescribed irradiance distribution E(x), x ∈ D in a target plane z = f, where x = (x₁, x₂) are the Cartesian coordinates in this plane (Fig. 1). We will assume that the target plane lies in the far field, so that the dimensions of the refractive surface can be neglected as compared with the distance f [i. e. that f ≫ ρ(m)]. Using the far-field condition, we describe the direction of the ray refracted by the surface R⃗(m), which then intersects the plane z = f at the point x = (x₁, x₂), by the following unit vector:

\vec{p} (x) = (x_{1}, x_{2}, f) / \sqrt{x_{1}^{2} + x_{2}^{2} + f^{2}} .

Let us note that the problem of generating a required irradiance distribution in the far field can also be considered as a problem of generating the intensity distribution

L (x) = E (x) \sqrt{x_{1}^{2} + x_{2}^{2} + f^{2}}

, x ∈ D.

In [24], it was shown that the refractive surface that forms a prescribed irradiance distribution in the far field can be represented as an envelope of a two-parameter family of ellipsoids of revolution. An ellipsoid, in one of the foci of which a point light source is located, converts a spherical beam from the source into a collimated beam with the propagation direction p⃗(x), which coincides with the direction of the axis of the ellipsoid. The radius vector of such an ellipsoid has the form [24,25]

{\vec{R}}_{el} (m; x) = \vec{e} (m) ρ_{el} (m; x) = \vec{e} (m) \frac{τ (x)}{1 - ν \vec{e} (m) \cdot \vec{p} (x)},

where e⃗(m) is defined in Eq. (1), ν = n₂/n₁ is the ellipsoid eccentricity, τ(x) is the focal parameter, and the dot denotes the scalar product of two vectors. Here, the coordinates x = (x₁, x₂) ∈ D are considered as parameters. In this case, Eq. (4) defines a two-parameter family of ellipsoids with the axis directions p⃗(x) and focal parameters τ(x). The envelope of this family is defined by the equation R⃗(m) = R⃗_el(m; x(m)), where the function x(m) is found from the following pair of equations [24]:

\frac{\partial ρ_{el} (m; x)}{\partial x_{i}} = 0, i = 1, 2 .

By expressing ρ_el(m; x) from Eq. (4) and substituting it into Eq. (5), after some simple transformations we obtain:

\frac{\partial log τ (x)}{\partial x_{i}} = - \frac{\vec{e} (m) \cdot \frac{\partial \vec{p} (x)}{\partial x_{i}}}{ν^{- 1} - \vec{e} (m) \cdot \vec{p} (x)}, i = 1, 2 .

Taking into account Eqs. (1) and (3), which define the vectors e⃗(m) and p⃗(x), we get

\frac{\partial log τ (x)}{\partial x_{i}} = - \frac{m_{i} - p_{i} (x) \vec{e} (m) \cdot \vec{p} (x)}{ν^{- 1} - \vec{e} (m) \cdot \vec{p} (x)} \cdot \frac{p_{3} (x)}{f}, i = 1, 2 .

Note that Eqs. (5)–(7) describe a critical point of the function ρ_el(m; x) of Eq. (4) with respect to the variables (x₁, x₂) at a fixed value of m = (m₁, m₂). In the present work, we consider the case when this critical point is either a maximum or a minimum of the function ρ_el(m; x). In this case, we can represent the length of the radius vector of the surface as

ρ (m) = max_{x \in D} \frac{τ (x)}{1 - ν \vec{e} (m) \cdot \vec{p} (x)},

or

ρ (m) = min_{x \in D} \frac{τ (x)}{1 - ν \vec{e} (m) \cdot \vec{p} (x)} .

Equations (8) and (9) allow one to find the radius vector of the refractive surface [see Eq. (2)] from the function τ(x), x ∈ D. In what follows, we will refer to τ(x) as the focal function.

Since ρ(m) in Eq. (8) is an “upper” envelope of a family of ellipsoids, for each m there exists an ellipsoid R(m; x) that is tangent to R⃗(m). Therefore, the surface defined by Eq. (8) induces the ray mapping γ in the following way:

γ : m \mapsto \underset{x \in D}{arg max} \frac{τ (x)}{1 - ν \vec{e} (m) \cdot \vec{p} (x)},

according to which the point m ∈ Ω (coordinates of the incident beam) is mapped to the point γ(m) = (x₁(m), x₂(m)) ∈ D (coordinates of the intersection of the refracted beam with the plane z = f). By replacing arg max in Eq. (10) with arg min, we will obtain the ray mapping defined by Eq. (9).

To design a refractive surface generating the prescribed irradiance distribution E(x), x ∈ D, one needs to find a focal function τ(x), x ∈ D, which obeys the light flux conservation law. This law can be formulated as follows: for any measurable subset S ⊂ D, the equality

\int_{γ^{- 1} (S)} I (m) d S_{Ω} = \int_{S} E (x) d x,

must hold, where

d S_{Ω} = d m / \sqrt{1 - m_{1}^{2} - m_{2}^{2}}

is the solid angle differential. The light flux conservation law of Eq. (11) can be rewritten in the differential form:

I_{pr} (m) = E (γ (m)) \cdot | J_{γ} (m) |,

where

I_{pr} (m) = I (m) / \sqrt{1 - m_{1}^{2} - m_{2}^{2}}

and |J_γ(m)| is the Jacobian of the coordinate transformation x(m).

Let us note that not every mapping γ(m) satisfying the light flux conservation law of Eq. (12) is integrable, i.e. can be implemented by a refractive surface. However, γ(m) is integrable once it is defined by Eq. (10) through a focal function τ(x) satisfying Eq. (6) at x = γ(m). Thus, the problem of calculating the refractive surface is reduced to finding such an integrable mapping γ.

3. Mass transportation problem

The problem of calculating the refractive surface considered in Section 2 can be formulated as the Monge–Kantorovich mass transportation problem (MTP) with the cost function [16,17]

𝒦 (m, x) = - log {1 - ν \vec{e} (m) \cdot \vec{p} (x)} .

In order to formulate the MTP, let us introduce the concept of a transport plan. Such a plan is a certain differentiable bijection γ : m ↦ x, which conserves the measure (the light flux) in the sense of Eq. (12). Let us define the following functional on the set of plans:

ℐ (γ) = \int_{Ω} 𝒦 (m, γ (m)) I_{pr} (m) d m .

Finding a solution to the MTP consists in determining a mapping γ(m), which minimizes or maximizes the functional of Eq. (14) provided that Eq. (12) holds.

Let us note that the formulation of the problem of calculating a refractive surface as an MTP with the cost function of Eq. (13) [16,17] is known only to a limited circle of researchers working in the field of mathematical theory of optics. This is due to the fact that the papers [16, 17] have purely mathematical character and were published in mathematical journals. Therefore, in what follows, we give another proof of the equivalence of the problem of calculating a refractive surface and the corresponding mass transportation problem, which is oriented at the applied optics community. The presented proof is based on the direct calculation of the extremum of the functional of Eq. (14). We further show that the equations describing the extremum of the functional of Eq. (14) coincide with Eqs. (6) and (12), which determine the solution of the problem of calculating the refractive surface.

In order to find the extremum of the functional of Eq. (14), let us use the method of Lagrange multipliers [26]. According to this method, the extremum of the functional of Eq. (14) under the condition of Eq. (12) (i.e. with an equality constraint) is a stationary point of the functional

ℱ (γ, μ) = \int_{Ω} {𝒦 (m, γ (m)) I_{pr} (m) + μ (m) [I_{pr} (m) - E (γ (m)) \cdot | J_{γ} (m) |]} d m,

where the function μ(m) is the Lagrange multiplier. To find a stationary point, we equate to zero the variations of the functional with respect to the function μ and to the mapping γ.

The variation with respect to μ can be easily calculated in the following form:

δ_{μ} ℱ (γ, μ) = ε \int_{Ω} η (m) [I_{pr} (m) - E (γ (m)) \cdot | J_{γ} (m) |] d m,

where εη(m) is the variation of the function μ(m). According to the fundamental lemma of the calculus of variations, the equality of the variation δ_μℱ(γ, μ) to zero implies that the light flux conservation law of Eq. (12) holds.

For the derivation of the variation of the functional ℱ(γ, μ) with respect to γ, it is convenient to represent the function μ(m) in the form μ(m) = log λ (γ(m)), where λ(x) is a certain function. Using this representation, we obtain the following functional:

ℱ (γ, λ) = \int_{Ω} [𝒦 (m, γ (m)) + log λ (γ (m))] I_{pr} (m) d m - \int_{Ω} log λ (γ (m)) E (γ (m)) \cdot | J_{γ} (m) | d m .

By making the change of variables m = γ⁻¹(x) in the second integral, we obtain

ℱ (γ, λ) = \int_{Ω} [𝒦 (m, γ (m)) + log λ (γ (m))] I_{pr} (m) d m - \int_{D} log λ (x) E (x) d x .

Next, let us calculate the variation of the functional of Eq. (18) with respect to the mapping γ. Let γ_ε(m) = γ(m) + ε · ω(m) denote the varied mapping, where γ(m) = (x₁(m), x₂(m)), ω(m) = (ω₁(m), ω₂(m)). First, we calculate the variation δ_γ e⃗(m) · p⃗(γ(m)):

\begin{array}{l} δ_{γ} \vec{e} (m) \cdot \vec{p} (γ (m)) & = δ_{γ} \sum_{j = 1}^{3} e_{j} (m) p_{j} (γ (m)) = \sum_{j = 1}^{3} e_{j} (m) δ_{γ} p_{j} (γ (m)) \\ = ε \sum_{j = 1}^{3} e_{j} (m) [\frac{\partial p_{j}}{\partial x_{1}} (γ (m)) ω_{1} (m) + \frac{\partial p_{j}}{\partial x_{2}} (γ (m)) ω_{2} (m)] \\ = ε \sum_{i = 1}^{2} \vec{e} (m) \cdot \frac{\partial \vec{p}}{\partial x_{i}} (γ (m)) ω_{i} (m) . \end{array}

Using this expression, one can easily find the variation δ_γ 𝒦(m, γ(m)) in the form

δ_{γ} 𝒦 (m, γ (m)) = \frac{ε}{ν^{- 1} - \vec{e} (m) \cdot \vec{p} (γ (m))} \sum_{i = 1}^{2} \vec{e} (m) \cdot \frac{\partial \vec{p}}{\partial x_{i}} (γ (m)) ω_{i} (m) .

The last equation enables writing the variation of the functional of Eq. (18) in the following form:

δ_{γ} ℱ (γ, λ) = ε \sum_{i = 1}^{2} \int_{Ω} [\frac{\vec{e} (m) \cdot \frac{\partial \vec{p}}{\partial x_{i}} (γ (m))}{ν^{- 1} - \vec{e} (m) \cdot \vec{p} (γ (m))} + \frac{\partial log λ}{\partial x_{i}} (γ (m))] ω_{i} (m) I_{pr} (m) d m .

By equating the variation δ_γℱ(γ, λ) to zero, we arrive at Eq. (6) with x = γ(m) and τ(x) = λ(γ(m)). This equation describes the integrability condition of the mapping γ(m) in the considered problem. Thus, the focal function τ(x), which defines the focal parameters of the ellipsoids in Eq. (4), coincides with the function λ(x) determining the Lagrange multiplier in the functional of Eq. (17).

Thus, we have shown that the mapping γ(m), which is found by solving an MTP with the cost function of Eq. (13), is integrable. After calculating the inverse mapping m = γ⁻¹(x), we can use Eq. (6) to retrieve the focal function τ(x). Finally, the refractive surface is reconstructed using either Eq. (8) or Eq. (9). Following the results of [16,17], one can show that the problem of maximizing the functional of Eq. (14) corresponds to a refractive surface defined by Eq. (8), while the refractive surface of Eq. (9) corresponds to the problem of minimizing the functional of Eq. (14).

According to Eq. (13), the cost function reaches its maximum at e⃗(m) · p⃗(x) = 1, i.e. when the angle δ between the incident and refracted rays equals zero. With an increase in the absolute value of the angle |δ|, the cost function decreases. Therefore, from the form of the cost function, it follows that the solution of the problem of maximizing the functional of Eq. (14) corresponds to a mapping, in which the total ray deflection is minimized. This case corresponds to the so-called Galilean configuration of rays [27]. Accordingly, the minimization of the functional of Eq. (14) leads to the Keplerian configuration of rays, in which the total ray deflection is maximized.

In the derivations given above, we assumed that the mapping γ(m) is differentiable. However, one can show that the MTP of Eqs. (13) and (14) can be used for calculating refractive surfaces implementing not only smooth mappings γ(m), but also mappings that have discontinuities on a set of measure zero. It follows from the fact that the mapping γ(m) is included only in the integrands of the Lagrange functional of Eq. (17). In the case of a discontinuous mapping γ(m), refractive surfaces calculated by Eqs. (8) and (9) turn out to be piecewise-smooth.

Let us also note that an integrable ray mapping γ(m) is obtained not only at a global minimum or maximum of the initial functional of Eq. (14), but also at any stationary point of the Lagrange functional of Eq. (17). In the latter case, the general Eq. (4) with x = γ(m) has to be used for the calculation of the refractive surface.

4. Reduction to the linear assignment problem

In this section, we show that the discrete version of the mass transportation problem, which consists in minimizing or maximizing the functional in Eq. (12), can be formulated as a linear assignment problem (LAP) [18–20].

Let us divide or approximate the domains Ω and D (on which the functions I(m) and E(x) are defined) into N equal-flux cells, with the equality $\int_{ω_{i}} I_{pr} (m) d m = \int_{d_{j}} E (x) d x = σ$ being valid for any cells ω_i ⊂ Ω, d_j ⊂ D. In the simplest case, when the functions I_pr(m) and E(x) are constant (I_pr(m) = I_pr, m ∈ Ω; E(x) = E, x ∈ D), the partition of the domains Ω and D into N equal-area cells can be obtained by choosing the steps of rectangular grids superimposed on the said domains [18, 19]. An example of the approximation of a circular domain Ω and a square domain D with equal-flux cells is shown in Figs. 2(a) and 2(c).

Fig. 2 Equal-flux-cell approximations of the circular domain Ω corresponding to the Lambertian source O (a), of the domain Ω corresponding to the virtual source O′ (b), and of the square domain D (c). Dots correspond to the centers of the cells. (d)–(f) Ray mapping of a square mesh from the source O (d) to the virtual source O′ (e) and to the target plane (f). The presented mappings correspond to the optical element designed in Subsection 5.2.

Download Full Size | PDF

In the case of non-uniform distributions I_pr(m) and E(x), an equal-flux set of N cells in the domain Ω and D can be obtained using one of the approaches utilized in the so-called source-target mapping methods based on the construction of an equal-flux grid in the source and target domains [28].

For the grids (approximations) introduced in the domains Ω and D, all energy-conserving discrete mappings γ_d : m ↦ x can be described by permutations of N numbers (j₁, . . ., j_N), which determine, to which cells d_j of the domain D the cells ω₁, . . ., ω_N are mapped. In terms of the LAP, the mapping of the cell ω_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. (13), the cost of the tasks is described by the matrix

{(C)}_{i, j} = 𝒦 (m_{i}, x_{j}), i, j = 1, \dots, N,

where m_i and x_j are the centers of the cells ω_i and d_j. We can write the discrete version of the functional of Eq. (14) in the form

ℐ_{d} (j_{1}, \dots, j_{N}) = \sum_{i = 1}^{N} 𝒦 (m_{i}, x_{j_{i}}),

which is defined on the permutations (j₁, . . ., j_N). In Eq. (23), we omitted the single-cell flux σ since it is not important for the further analysis. Thus, the discrete approximation γ_d of the mapping γ can be found by solving one of the following linear assignment problems:

ℐ_{d} (j_{1}, \dots, j_{N}) \to max,

ℐ_{d} (j_{1}, \dots, j_{N}) \to min,

where the maximum or minimum is sought over all permutations (j₁, . . ., j_N). The obtained discrete mapping γ_d(m_i) = x_{j_i} is defined on the cell centers m_i, i = 1, . . ., N. Note that the LAP of Eq. (24) defines the mapping corresponding to the Galilean configuration of rays, while the LAP of Eq. (25) leads to the mapping corresponding to the Keplerian configuration of rays.

For the solution of the LAP, efficient polynomial algorithms are available, in particular, the Hungarian algorithm [20] and the auction algorithm [21].

5. Design examples

5.1. Refractive surface generating uniform irradiance distribution in a square region

As a first example, let us consider the calculation of a refractive surface generating a uniform irradiance distribution in the square D = {(x₁, x₂) ||x₁| ≤ w/2, |x₂| ≤ w/2} located in the plane z = f. We assume a Lambertian light source radiating into a solid angle corresponding to a cone with the apex at the origin and the apex angle (the radiation angle) 2θ₀ (Fig. 1). The axis of the cone coincides with the z axis. In this case, the domain $Ω = {(m_{1}, m_{2}) | m_{1}^{2} + m_{2}^{2} \leq {sin}^{2} θ_{0}}$ is a circle and the intensity of the Lambertian source has the form $I (m) = \sqrt{1 - m_{1}^{2} - m_{2}^{2}}$ , m ∈ Ω. Thus, the function $I_{pr} (m) = I (m) / \sqrt{1 - m_{1}^{2} - m_{2}^{2}}$ used in the light flux conservation law of Eq. (12) turns out to be constant: I_pr(m) = 1.

The calculation was carried out for the following parameters: the source radiation angle 2θ₀ = 90°, the square size w = 1200 mm, distance to the target plane f = 1050 mm, refractive indices of the media n₁ = 1.5 and n₂ = 1.

In the general case, the calculation of the refractive surface includes the following steps:

Calculation of the discrete mapping γ_d by solving the LAP of Eq. (24) or (25).
Calculation of the focal function τ(x) from the discrete mapping γ_d using Eq. (6).
Calculation of the refractive surface using Eq. (8) [when using Eq. (24) in step 1] or Eq. (9) [when using Eq. (25) in step 1] from the calculated focal function τ(x).

Let us describe these steps in more detail for the example being considered.

For the calculation of the mapping x = γ(m), the square D was represented as a set of N = 250² = 62500 square cells d_i with the size 4.8 × 4.8 mm². The circular domain Ω was also approximated by 62500 square cells ω_i with the size 0.005 × 0.005. The discrete mapping γ_d(m_i) = x_{j_i} was found from the solution of the LAP of Eq. (24) using the auction algorithm [21] implemented in [29]. For the considered example, the solution of the LAP on a desktop PC (Intel Core i7-3930K, 16 GB RAM) took several minutes. As noted above, the solution to the LAP of Eq. (24) is a mapping, for which the total ray deflection is minimized. This enables reducing Fresnel losses that are not taken into account in the cost function defined by Eq. (13) [27]. As we demonstrate below, the chosen value of N (the dimension of the LAP being solved) provides a high-quality irradiance distribution.

The focal function τ(x) can be calculated from the mapping γ(m) by numerical integration of Eq. (6). However, due to the discrete character of the LAP of Eq. (24), the obtained discrete mapping γ_d describes the sought-for continuous mapping γ with an error [19,23]. In this regard, the focal function τ(x) was calculated using a regularized least-squares method [19,23]. In this method, the function log τ(x) is represented by a system of two-dimensional B-splines:

log τ (x) = \sum_{m, n} p_{m n} B_{m} (x_{1}) P_{n} (x_{2}),

where B_m(x₁) and P_n(x₂) are the B-splines of a certain order q defined in the domain D on a rectangular grid containing N₁ ·N₂ vertices. The coefficients of the B-splines p_mn were determined from the condition that the derivatives of the splines (26) provide a least-squares approximation of the derivatives of the function log τ(x) given by Eq. (7) at x = x_{j_i},

m = γ_{d}^{- 1} (x_{j_{i}})

[23]. In the considered example, B-splines of the 7-th order with N₁ = N₂ = 10 were used for the representation of the function log τ(x).

From the obtained focal function τ(x), a refractive surface was calculated using Eq. (8). Equations (8) and (26) allow one to calculate the refractive surface at any m ∈ Ω and, thus, at an arbitrarily fine mesh required for the simulation or fabrication of the optical element. The designed surface is shown in Fig. 3(a). Let us note that the radius vector length ρ(m) of the refractive surface is defined up to a multiplicative constant, therefore, the designed optical element is defined up to a scale factor. This allows one to set the length of the radius vector of the surface at a certain point, e.g. at m = (0, 0). The surface shown in Fig. 3(a) was calculated at ρ(0, 0) = 3 mm. In this case, the dimensions of the refractive surface amount to 3.7×3.7×1.1 mm³. The designed surface is smooth; its Gaussian curvature varies from 0.15 mm⁻² to 0.33 mm⁻².

Fig. 3 Refractive surface generating a uniform irradiance distribution in a square region (a) and the simulated irradiance distribution (b).

Download Full Size | PDF

Figure 3(b) shows the simulated irradiance distribution generated in the target plane z = f = 1050 mm by the designed optical surface. For the calculation of the irradiance distribution, an in-house implementation of the ray tracing technique was used, which takes into account the Fresnel losses. The distribution shown in Fig. 3(b) was calculated by tracing 10⁷ rays and almost coincides with a perfect uniform square distribution. The simulations of this element using a commercial ray-tracing software TracePro (not presented in this article) confirm the presented results.

5.2. Refractive optical element with two surfaces generating uniform irradiance distribution in a square region

In the example considered above, the source emits into a cone with the apex angle 2θ₀ = 90°, and the angular dimension of the side of the illuminated square domain amounts to α_sq = 2 arctan w/(2f) = 60°. For the practical applications, it is worth considering the case when the source emits to the hemisphere z > 0. In this case, the source radiation angle is 2θ₀ = 180°. At 2θ₀ = 180°, generation of the required irradiance distribution in a domain with the angular size α_sq = 60° with the use of a single refractive surface is impossible. Indeed, in order to illuminate a domain with the angular size α_sq = 60°, it is necessary to deflect the boundary rays from the source (the rays perpendicular to the z axis) by the angle 90° − 60°/2 = 60°. However, the maximum deflection angle upon refraction at the interface of two media with the relative index of refraction 1.5 amounts to α_rot = 90° − arcsin(1/1.5) = 48.2°. Thus, efficient generation of illuminated regions with the angular size less than 2 · (90° − α_rot) = 83.6° is impossible, if a single refractive surface is used. In order to generate irradiance distributions with a small angular size, it is necessary to use optical elements with two “working” surfaces [22].

Let us consider the utilization of the LAP-based approach for the calculation of optical elements with two working surfaces. As the first (lower) surface, let us use a refractive surface generating a virtual image O′ = (0, 0, −Δ), Δ > 0 of the source O = (0, 0, 0) (Fig. 4). We assume that both the source is located in a medium with unit refractive index. The lower surface forms a virtual source with a decreased radiation angle 2θ_0,v < 180°. The angle θ_0,v is determined by the source offset Δ and the z-coordinate of the central point of the lower surface r₀ [Fig. 4(b)]. Then, the second surface of the optical element is calculated as in the above example with the difference that the virtual source O′ is used in the calculations.

Fig. 4 (a) Profile of an optical element with two working surfaces. The lower surface forms a virtual image O′ of the source O. (b) Geometry of refraction of a ray at the lower surface.

Download Full Size | PDF

The lower surface of the element is a surface of revolution with the profile being a Cartesian oval [25]:

\sqrt{x^{2} + z^{2}} - n_{1} \sqrt{x^{2} + {(z + Δ)}^{2}} = c_{0},

where c₀ = r₀ − n₁(r₀ + Δ). In the polar coordinates x = r sin θ, z = r cos θ, where the angle θ is measured from the z axis, one can easily obtain the following expression for the length of the radius vector of the lower surface:

r (θ) = \frac{- c_{0} - n_{1}^{2} Δ cos θ + n_{1} \sqrt{{(c_{0} + Δ cos θ)}^{2} - (n_{1}^{2} - 1) Δ^{2} {sin}^{2} θ}}{n_{1}^{2} - 1}, θ \in [0, \frac{π}{2}] .

The surface defined by Eq. (28) forms a virtual source O′ with a decreased emission angle 2θ_0,v = 2 arctan(r(π/2)/Δ). Let us note that the use of the Cartesian oval as the first surface may be non-optimal when generating irradiance distributions with a high aspect ratio [6].

Let us consider the design of an optical element with two working surfaces for generating a uniform irradiance distribution in a square for a Lambertian source with 2θ₀ = 180°. All other parameters remain the same as in Subsection 5.1. As it was mentioned above, the angular size of the side of the illuminated square amounts to α_sq = 60°, therefore, it cannot be generated by an element with a single working refractive surface.

The lower surface was calculated using Eq. (28) at Δ = 0.7 mm and r₀ = r(0) = 0.5 mm. The calculated profile is shown with a blue curve in Fig. 4. In the considered example, the lower surface decreases the radiation angle from 2θ₀ = 180° to 2θ_0,v ≈ 146°. The calculation of the second (upper) surface is similar to the calculation of the refractive surface in the previous example with the difference that the virtual source O′ is used as the light source. Let us note that the virtual source O′ is not Lambertian, and therefore, in contrast to the previous example, the function I_pr(m) is not constant. In this regard, an equal-flux grid in the domain Ω, which is used in the LAP of Eq. (24), has to be constructed. In Appendix A, it is shown that an equal-flux grid in the domain Ω corresponding to the virtual source O′ can be obtained by a simple transformation [Eq. (30)] of the equal-flux grid constructed for the Lambertian source O.

As an example, Figs. 2(a) and 2(c) show the approximations of the domain Ω corresponding to the actual source O, and of the domain D by equal-flux grids of N = 24² = 576 square cells. Fig. 2(b) shows the equal-flux grid calculated using Eq. (30) for the domain Ω corresponding to the virtual source O′. From Fig. 2(b), it is evident that the function I_pr(m) corresponding to the virtual source O′ is not constant but has a pronounced maximum in the vicinity of the point m = (0, 0). In the case of a Lambertian source, it is easy to show that sin²(60°/2) · 100% = 25% of the light flux is emitted into a cone with the apex angle 60°. At the same time, according to Eq. (29), the virtual source emits 60% of the light flux into the same cone. Thus, the virtual source O′ turns out to be significantly “more collimated” than the actual Lambertian source O.

The equal-flux grid with N = 576 points is shown in Figs. 2(a)–2(c) for illustrative purposes only and is too coarse for the calculation of the outer surface of the optical element. The solution of the LAP of Eq. (24) and the calculation of the focal function τ(x) for the upper surface were carried out for the values N = 62500, N₁ = N₂ = 10, q = 7 coinciding with the parameters of the example discussed in Subsection 5.1. Figures 2(d)–2(f) illustrate the ray mappings from the source O to the virtual source O′ and to the target plane. The latter mapping was reconstructed from the solution of the LAP of Eq. (24).

Figure 5(a) shows the calculated optical element with two working surfaces. Let us note that it is the xz-cross-section of the designed element that is shown in Fig. 4(a). The dimensions of the optical element amount to 5.5 × 5.5 × 3.2 mm³. The designed optical element has smooth surfaces: the Gaussian curvature of the outer surface ranges from 0 to 0.25 mm⁻², whereas the inner surface has curvature ranging from −0.15 mm⁻² to 2.24 mm⁻².

Fig. 5 Refractive optical element with two surfaces generating a uniform irradiance distribution in a square region (a) and the simulated irradiance distribution (b).

Download Full Size | PDF

Figure 5(b) shows the irradiance distribution generated in the plane z = f = 1050 mm by the designed optical element. As in the previous example, the irradiance distribution was calculated by tracing 10⁷ rays and almost coincides with a perfectly uniform square distribution. The lighting efficiency (the percentage of the emitted light flux which is directed to the target domain D) amounts to 89.8%. Let us note that due to Fresnel losses, at least 4% of the energy of the incident beam is lost at each of the surfaces of the element. Thus, the obtained lighting efficiency is close to the theoretical limit of 92%.

5.3. Refractive optical element with two surfaces generating uniform irradiance distribution in a complex-shaped area

As the third, the most complex example, we designed a refractive optical element generating uniformly illuminated letters “AB” on a zero background (Fig. 6). The letters are bounded by a rectangle with the dimensions of 1200 × 650 mm² and are located in the plane f = 1050 mm. As in the previous example, the Lambertian source emits into a hemisphere (2θ₀ = 180°). Note that the angular size of the domain containing the letters amounts to 60° × 35°. Therefore, it is necessary to use an optical element with two working surfaces.

Fig. 6 Refractive optical element with two working surfaces generating uniformly illuminated letters “AB” (a) and the simulated irradiance distribution (b). The inset shows the sharp bends on the surface of the optical element.

Download Full Size | PDF

The lower surface of the optical element was calculated using Eq. (28) with the parameters Δ = 0.6 mm and ρ₀ = 0.5 mm, which coincide with the parameters of the previous example. For the calculation of the upper surface, the mapping x = γ(m) was obtained from the solution of the LAP of Eq. (24). The domains Ω and D were approximated by N = 300² = 90000 square cells. For the calculation of the focal function τ(x), the function log τ(x) was approximated by a system of two-dimensional B-splines of Eq. (26) of the 20-th order with N₁ = N₂ = 50.

Figure 6 shows the designed optical element [Fig. 6(a)] and the generated irradiance distribution [Fig. 6(b)]. The irradiance distribution was calculated by tracing 10⁷ rays and with high accuracy corresponds to the prescribed distribution in the form of the letters “AB”. The lighting efficiency of the calculated element equals 87.3%. The lower lighting efficiency as compared with the previous example is due to the smaller angular size of the illuminated domain. A decrease in the angular size of the generated irradiance distribution leads to an increase in the refraction angles of the rays at the outer optical surface, and, consequently, to increased Fresnel losses [27]. Let us note that in this example, the mapping γ(m) implemented by the upper surface of the optical element is discontinuous due to the disjoint form of the generated irradiance distribution. Therefore, the upper surface of the designed optical element has several sharp bends [19], which are depicted in Fig. 6. At these bends, the curvature of the surface is negative infinity, whereas the maximal Gaussian curvature reaches 0.40 mm⁻². The total dimensions of the optical element amount to 6.1 × 5.3 × 4.0 mm³.

The presented example demonstrates one of the advantages of the developed method, which consists in the possibility of designing piecewise-smooth continuous optical surfaces. In this case, the known methods based on numerical solution of an elliptic partial differential equation [1–13] are either inapplicable or are numerically unstable. Indeed, in these methods, the optical surface is assumed to be represented by a smooth twice differentiable function. It significantly limits the class of the generated irradiance distributions and, in contrast to the proposed method, does not allow the formation of not connected irradiance distributions such as letters on a zero background.

6. Conclusion

In the present work, we demonstrated that the problem of calculating a refractive surface generating a prescribed irradiance distribution in the far field can be formulated as a mass transportation problem with a non-quadratic cost function. This formulation allowed us to use an LAP-based method for designing refractive optical elements. In this method, the calculation of the ray mapping between the rays emitted by the source and the points in the target plane is reduced to finding a solution to a linear assignment problem. We considered the application of the method for the calculation of optical elements with two refractive surfaces. High performance of the proposed method was illustrated by the designed examples of optical elements generating uniform irradiance distributions in a square and in the region in the form of letters “AB”.

Appendix

Let us consider the construction of an equal-flux grid for the domain Ω, which corresponds to the virtual light source O′ (Fig. 4). This virtual source O′ = (0, 0, −Δ), Δ > 0 is generated by the surface of Eq. (28) from the source O = (0, 0, 0). It follows from the geometry of the problem in Fig. 4(b) that the ray initially emitted by the source O at an angle θ will, after refraction by the surface defined by Eq. (28), correspond to a ray emitted from the virtual source O′ at the angle

θ_{v} (θ) = arcsin (sin θ \cdot r (θ) / t (θ)),

where

t (θ) = \sqrt{r^{2} (θ) + Δ^{2} + 2 r (θ) Δ cos θ}

, and r(θ) is the length of the radius vector of the profile of the lower optical surface defined by Eq. (28). The ray mapping θ_v = θ_v(θ) allows one to obtain an equal-flux grid for the virtual source O′ from the equal-flux grid in the domain Ω corresponding to the actual source O. Indeed, let m_i = (m_1,i, m_2,i), i = 1, . . ., N be the centers of the equal-flux cells in the domain Ω. Since

m_{1, i}^{2} + m_{2, i}^{2} = {sin}^{2} θ_{i}

, where θ_i is the angle between the ray e⃗(m_i) and the z axis, the cell centers of an equal-flux grid m̂_i, i = 1, . . ., N for the virtual source O′ can be obtained from m_i by the following simple transformation:

{\hat{m}}_{i} = m_{i} r (θ_{i}) / t (θ_{i}), i = 1, \dots, N .

One can easily show that from Eq. (30), it follows that ${\hat{m}}_{1, i}^{2} + {\hat{m}}_{2, i}^{2} = {sin}^{2} θ_{v} (θ_{i})$ , where the function θ_v = θ_v(θ) is defined by Eq. (29). Thus, for the construction of an equal-flux grid in the domain Ω for the virtual source O′, we first construct an equal-flux grid for the actual source O. In the case of a Lambertian source O emitting into a hemisphere, the initial domain Ω is a circle with unit radius, and the function I_pr(m) is constant. Then, the obtained centers of the equal-flux cells are transformed using Eq. (30).

Funding

Russian Foundation for Basic Research (18-29-03067, 18-07-00982, 18-07-00514); 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 optical surfaces was supported by Russian Foundation for Basic Research, the application of LAP-based approach to the design of double-surface elements was supported by Russian Science Foundation, and the implementation of the ray-tracing simulation software 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 Photon. 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. 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]

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

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

7. C. R. Prins, J. H. M. ten Thije Boonkkamp, J. van Roosmalen, W. L. IJzerman, and T. W. Tukker, “A Monge–Ampèresolver for free-form reflector design,” SIAM J. Sci. Comput. 36(3), B640–B660 (2014). [CrossRef]

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

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

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

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

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

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

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

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

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

17. C. E. Gutiérrez and Q. Huang, “The refractor problem in reshaping light beams,” Arch. Ration. Mech. Anal. 193, 423–443 (2009). [CrossRef]

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

19. D. A. Bykov, L. L. Doskolovich, A. A. Mingazov, E. A. Bezus, and N. L. Kazanskiy, “Linear assignment problem in the design of freeform refractive optical elements generating prescribed irradiance distributions,” Opt. Express 26(21), 27812–27825 (2018). [CrossRef] [PubMed]

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

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

22. M. A. Moiseev, S. V. Kravchenko, and L. L. Doskolovich, “Design of efficient LED optics with two free-form surfaces,” Opt. Express 22(S7), A1926–A1935 (2014). [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] [PubMed]

24. L. L. Doskolovich, E. S. Andreev, S. I. Kharitonov, and N. L. Kazansky, “Reconstruction of an optical surface from a given source-target map,” J. Opt. Soc. Am. A 33(8), 1504–1508 (2016). [CrossRef]

25. R. K. Luneburg, Mathematical Theory of Optics (University of California, 1964).

26. 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).

27. V. Oliker, L. L. Doskolovich, and D. A. Bykov, “Beam shaping with a plano-freeform lens pair,” Opt. Express 26(15), 19406–19419 (2018). [CrossRef] [PubMed]

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

29. Fast linear assignment problem using auction algorithm (mex). http://www.mathworks.com/matlabcentral/fileexchange/48448

Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions

Abstract

1. Introduction

2. Formulation of the problem

3. Mass transportation problem

4. Reduction to the linear assignment problem

5. Design examples

5.1. Refractive surface generating uniform irradiance distribution in a square region

5.2. Refractive optical element with two surfaces generating uniform irradiance distribution in a square region

5.3. Refractive optical element with two surfaces generating uniform irradiance distribution in a complex-shaped area

6. Conclusion

Appendix

Funding

Acknowledgments

References

Cited By

Figures (6)

Equations (30)

Optics Express