Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

On the use of the supporting quadric method in the problem of designing double freeform surfaces for collimated beam shaping

Open Access Open Access

Abstract

We propose a version of the supporting quadric method for calculating a refractive optical element with two working surfaces for collimated beam shaping. Using optimal mass transportation theory and generalized Voronoi cells, we show that the proposed method can be regarded as a gradient method of maximizing a concave function, which is a discrete analogue of the Lagrange functional in the corresponding mass transportation problem. It is demonstrated that any maximum of this function provides a solution to the problem of collimated beam shaping. Therefore, the proposed method does not suffer from “trapping” at a local extremum, which is typical for gradient methods. We present design examples of refractive optical elements illustrating high performance of the method.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

The problem of designing freeform optical elements generating light beams with prescribed irradiance distribution and wavefront has a wide range of applications including the design of lighting and backlight systems, laser beam processing and shaping, and optical lithography, among others. In the present work, we consider the problem of collimated beam shaping in the geometrical optics approximation. This problem consists in calculating a refracting (or reflecting) optical element, which transforms a given incident beam with a plane wavefront into an output beam with a prescribed irradiance distribution and a plane wavefront. Simultaneous control of both the irradiance distribution and the wavefront of the output beam requires the use of an optical element with two “working” optical surfaces.

In the geometrical optics approximation, the problem of collimated beam shaping can be reduced to finding a solution to a complex, fully nonlinear differential equation (NDE) of Monge–Ampère type. For the solution of this NDE, various finite-difference methods have been proposed [17]. The formulation of the problem of designing an optical element as an NDE assumes that the calculated surfaces of the optical element are smooth. This imposes a number of restrictions on the class of irradiance distributions that can be generated. In particular, an optical element with smooth surfaces cannot generate an irradiance distribution defined on a disconnected or a multiply connected region, or on a region with complex non-smooth boundary [8,9].

In the basic theoretical works [1013], the problem of collimated beam shaping (PCBS) was formulated as a Monge–Kantorovich mass transportation problem (MTP) for the cases of a refractive optical element [10,12,13] and a system of two mirrors [11]. This MTP describes a problem of calculating an integrable ray mapping (i. e. a mapping relating the coordinates of the rays incident on an optical element to the coordinates of the refracted (or reflected) rays in the target plane), which provides the formation of an output beam with the required irradiance distribution. The calculation of an optical element in this MTP-based approach can be reduced to finding a solution to a linear programming problem [9] or a linear assignment problem [8]. In contrast to the NDE-based formulation, the MTP-based formulation of the PCBS enables calculating optical elements with continuous piecewise-smooth optical surfaces generating prescribed irradiance distributions defined on disconnected regions and on regions with complex non-smooth boundaries [8,9].

One of the widely used methods of designing refracting and reflecting optical surfaces is the supporting quadric method (SQM) [1420]. It was initially proposed by V. Oliker and co-workers for solving the problem of calculating mirrors generating prescribed irradiance distributions [1416]. The method was then generalized to the case of refracting optical surfaces [1820]. The main advantages of the SQM are its versatility and simplicity. In the method, the prescribed irradiance distribution is approximated by a discrete distribution defined on a finite set of $N$ points. Then, the optical surface is represented as a set of $N$ segments of quadrics, each focusing the incident beam to one of these points. Depending on the problem being solved, paraboloids, ellipsoids, hyperboloids or more complex surfaces (e. g., Cartesian ovals [18]) are used as quadrics. In particular, in the problem of calculating a mirror generating a given discrete irradiance distribution in the far field, the mirror surface is represented as a set of segments of paraboloids. The calculation of the parameters of the quadrics is carried out iteratively. Note that in the above-mentioned problem of calculating a mirror, the convergence of this iterative method has been strictly proved [14].

In the present work, we consider a version of the SQM for designing a refractive optical element with two working surfaces for collimated beam shaping. The proposed method is simple to implement and is based on just a few main equations. We demonstrate that despite its simplicity, the method has firm theoretical foundations. Using the MTP-based formulation of the problem, we show that the proposed SQM is, in its essence, a gradient descent method of maximizing a concave function being a discrete analogue of the Lagrange functional in the MTP. It is important to note that any maximum (any critical point) of this function provides a solution to the PCBS. Therefore, the proposed version of the SQM does not suffer from the “false” convergence to a local maximum, which is typical for gradient methods. We illustrate the capabilities of the proposed method by designing optical elements confirming its high performance. The presented examples demonstrate the possibility of calculating optical elements with continuous piecewise-smooth surfaces generating prescribed irradiance distributions in doubly connected regions.

The paper is organized as follows. In Section 2, we formulate the problem of nonimaging optics being solved. In Section 3, we represent the surfaces of the designed element as envelopes of quadrics. Section 4 is dedicated to the description of the proposed version of the SQM. In Sections 5, we establish the connection between the SQM and the gradient method of maximizing a certain concave function associated with the MTP. In Section 6, we present two design examples illustrating the efficiency of the method. Section 7 concludes the paper.

2. Problem statement

The geometry of the considered problem is presented in Fig. 1(a). A light beam, which has a plane wavefront parallel to the plane $z=0$, impinges on an optical element with two refracting surfaces [Fig. 1(a)]. The irradiance distribution of the incident beam is defined by $I(\mathbf {x})$, $\mathbf {x} = (x_1, x_2)\in G$, where $G$ is a certain domain in the plane $z=0$. The first surface $R_1$ of the optical element is defined as $z=u_1 (\mathbf {x})$, $\mathbf {x}\in G$. It is convenient to define the second surface $R_2$ with respect to a certain “output” plane $z=f>0$, in which the irradiance distribution of the transmitted beam will be defined. Let us denote by $\mathbf {y} = (y_1, y_2)$ the Cartesian coordinates in this plane and define the second surface through the function $u_2 (\mathbf {y})$ as $z=f-u_2 (\mathbf {y})$. We assume that the refractive index of the medium between the surfaces [at $u_1 (\mathbf {x})\le z\le f-u_2 (\mathbf {x})$] equals $n_1$, and the refractive index of the surrounding medium [at $z < u_1 (\mathbf {x})$ and at $f-u_2 (\mathbf {x}) < z$] equals $n_0$.

 figure: Fig. 1.

Fig. 1. (a) Geometry of the problem of collimated beam shaping. (b) Representation of the first surface as the envelope of a family of surfaces focusing light to the points of the second surface (two-dimensional case). The envelope is shown with a dashed black curve. Solid colored curves show the cross-sections of different surfaces of the family.

Download Full Size | PPT Slide | PDF

An incident ray from a point $\mathbf {x}\in G$ of the plane $z=0$ is successively refracted by the surfaces $R_1$ and $R_2$ and arrives at some point $\mathbf {y}$ of the “target” domain $D$ in the output plane $z=f$ [Fig. 1(a)]. In this sense, the surfaces $R_1$ and $R_2$ define a ray mapping $\gamma : G\to D$, which relates the coordinates of the refracted rays $\mathbf {y}=\gamma (\mathbf {x})$ to the coordinates of the incident rays. This mapping determines the irradiance distribution $L(\mathbf {y})$, $\mathbf {y}\in D$ generated in the plane $z=f$ according to the light flux conservation law:

$$I(\mathbf{x})=L(\gamma (\mathbf{x}))\cdot J_{\gamma} (\mathbf{x}),$$
where $J_{\gamma } (\mathbf {x})$ is the Jacobian of the mapping $\gamma$. The light flux conservation law can also be written in the integral form
$$\int_{\gamma^{{-}1} (B)} I(\mathbf{x})\mathrm{d}{\mathbf{x}} = \int_{B} L(\mathbf{y})\mathrm{d}{\mathbf{y}},$$
where $B$ is an arbitrary subset of $D$.

The problem of collimated beam shaping, which we consider in the present work, is formulated in the following way. Assuming that the irradiance distribution of the incident beam $I(\mathbf {x})$, $\mathbf {x}\in G$ and the required irradiance distribution $L(\mathbf {y})$, $\mathbf {y}\in D$ in the output plane $z=f$ [Fig. 1(a)] are given, it is necessary to find such functions $u_1(\mathbf {x})$ and $u_2 (\mathbf {y})$ so that the transmitted light beam would have the given irradiance distribution $L(\mathbf {y})$ and a plane wavefront parallel to the plane $z=f$.

3. Representation of the optical surfaces

Let us consider a ray from a point $\mathbf {x}$ of the plane $z=0$, which, after being refracted by the surfaces $R_1$ and $R_2$, arrives at a point $\mathbf {y}=\gamma (\mathbf {x})$ of the plane $z=f$. Since the wavefronts of both the incident and the output beams are plane, the optical path length $F$ between these points is constant and does not depend on $\mathbf {x}$. This condition can be written as

$$n_0 u_1 (\mathbf{x})+n_1 \sqrt{\left\| \mathbf{y}-\mathbf{x}\right\|^{2} +\left[f-u_1 (\mathbf{x})-u_2 (\mathbf{y})\right]^{2}} +n_0 u_2 (\mathbf{y})=F,$$
where $\left \|\cdot \right \|$ is the Euclidean norm. By solving Eq. (3) with respect to the sum $u_1 (\mathbf {x})+u_2 (\mathbf {y})$, we obtain [10]
$$u_1 (\mathbf{x})+u_2 (\mathbf{y})={\mathcal{K}}(\mathbf{x}, \mathbf{y}),$$
where
$${\mathcal{K}}(\mathbf{x}, \mathbf{y})=\frac{n_1^{2} f-n_0 F}{n_1^{2} -n_0^{2}} -\frac{n_1}{\left|n_1^{2} -n_0^{2} \right|} \sqrt{\left(F-n_0 f\right)^{2} -\left(n_1^{2} -n_0^{2} \right)\left\| \mathbf{y}-\mathbf{x}\right\|^{2}}.$$
In what follows, we will assume that the square root argument is non-negative at $\mathbf {x}\in G$, $\mathbf {y}\in D$.

Let us show that the surfaces $R_1$ and $R_2$ of the optical element can be represented as envelope surfaces of a special kind. Indeed, for differentiable surfaces $R_1$ and $R_2$, the ray mapping $\mathbf {y}=\gamma (\mathbf {x})$ can be described in the following way. The ray coming from a point $\mathbf {x}\in G$, after being refracted at the point $\left (\mathbf {x}, u_1(\mathbf {x})\right )$ of the surface $R_1$, goes to the point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ of the surface $R_2$. Therefore, the surface $R_1$ can be represented as the envelope of a family of surfaces, each focusing the incident plane beam to one of the points $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ of the surface $R_2$ [Fig. 1(b)] [9].

The equation of the refracting surface $z=u_1^{\mathbf {y}} (\mathbf {x})$ focusing the incident beam to a certain point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ can be obtained by requiring the optical path length from the points $\mathbf {x}$ of the plane $z=0$ to the point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ to be constant and equal to $F-n_0 u_2 (\mathbf {y})$. Note that Eqs. (4) and (5) describe exactly this condition, and therefore the function $u_1^{\mathbf {y}} (\mathbf {x})$ has the form

$$u_1^{\mathbf{y}} (\mathbf{x})={\mathcal{K}}(\mathbf{x}, \mathbf{y})-u_2 (\mathbf{y}).$$
Substituting Eq. (5) into Eq. (6), one can easily obtain that the surface $z=u_1^{\mathbf {y}} (\mathbf {x})$ is an ellipsoid of revolution if $n_0 <n_1$, and a hyperboloid of revolution if $n_0 >n_1$ [21]. The axis of this ellipsoid or hyperboloid is parallel to the $z$ axis, and one of its foci coincides with the point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$.

Considering $\mathbf {y} = (y_1, y_2)$ as parameters, we can treat the equation $z=u_1^{\mathbf {y}} (\mathbf {x})$ as a two-parameter family of surfaces. The first surface of the designed optical element, which is described by the function $u_1 (\mathbf {x})$, is the envelope of this family [9] and is defined by the following equations [22]:

$$u_1 (\mathbf{x})={\mathcal{K}}(\mathbf{x}, \mathbf{y})-u_2 (\mathbf{y}),$$
$$\frac{\partial}{\partial y_i} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y})-u_2 (\mathbf{y})\right]=0,\; \; i=1,2.$$

Similarly, the second surface $R_2$ can be considered as the envelope of a family of surfaces

$$u_2^{\mathbf{x}} (\mathbf{y})={\mathcal{K}}(\mathbf{x}, \mathbf{y})-u_1 (\mathbf{x}),$$
each focusing a plane beam, which propagates “backwards” from the plane $z=f$, to one of the points $\left (\mathbf {x}, u_1 (\mathbf {x})\right )$ of the first surface $R_1$. In this case, the function $u_2 (\mathbf {y})$ is the envelope of the family of surfaces of Eq. (8) with respect to the parameters $\mathbf {x} = \left (x_1, x_2 \right )\in G$ and is defined by the following equations:
$$u_2 (\mathbf{y})={\mathcal{K}}(\mathbf{x}, \mathbf{y})-u_1 (\mathbf{x}),$$
$$\frac{\partial}{\partial x_i} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y})-u_1 (\mathbf{x})\right]=0, \; \; i=1,2.$$

In what follows, we will consider a particular case of the envelope of Eqs. (7a) and (7b), when the condition of Eq. (7b) corresponds not just to some critical point of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})-u_2 (\mathbf {y})$ with respect to $\mathbf {y}$, but to its minimum. Such a choice provides the so-called Galilean configuration of rays [10]. In this case, the envelope surface defined by Eqs. (7a) and (7b) takes the form

$$u_1 (\mathbf{x})=\mathop{\min}_{\mathbf{y}\in D} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y})-u_2 (\mathbf{y})\right].$$
The envelope surface corresponding to a maximum point of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})-u_2 (\mathbf {y})$ can be considered in a similar way and will provide the Keplerian configuration of rays. Equation (10) enables defining the ray mapping through the function $u_2 (\mathbf {y})$ in the form
$$\gamma (\mathbf{x})=\mathop{\arg \min}_{\mathbf{y}\in D} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y})-u_2 (\mathbf{y})\right].$$
For a bijective mapping $\gamma :G\to D$, the function $u_2 (\mathbf {y})$ can be expressed through $u_1 (\mathbf {x})$ as
$$u_2 (\mathbf{y})=\mathop{\min}_{\mathbf{x}\in G} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y})-u_1 (\mathbf{x})\right].$$
This can be easily proved by contradiction (see, e. g., a similar proof in Refs. [23,24]). Note that Eq. (12) is also a particular case of the envelope surface defined by Eqs. (9a) and (9b).

Equations (10) and (11) make it possible to formulate the problem of collimated beam shaping as a problem of calculating such a function $u_2 (\mathbf {y})$, so that the corresponding ray mapping (11) satisfies the light flux conservation law of Eq. (1) or Eq. (2). Let us note that not every mapping satisfying the light flux conservation law is integrable, i. e. can be implemented by a pair of optical surfaces. In order for the mapping $\gamma$ to be integrable, it is necessary for the function $u_2 (\mathbf {y})$ defining it to exist and to satisfy Eq. (7b). Therefore, below we will refer to Eq. (7b) as the integrability condition.

4. Supporting quadric method

The supporting quadric method (SQM) [1420] is widely used for designing mirrors and refracting surfaces generating prescribed irradiance distributions. In the method, the required continuous irradiance distribution is approximated by a discrete distribution consisting of a finite set of points. The designed optical surface is assumed to be piecewise-smooth and consisting of segments of quadrics. The parameters of the quadrics are found iteratively from the condition of obtaining the required light flux values at the points of the discrete distribution.

In this section, we consider a version of the SQM for solving the problem of collimated beam shaping [Fig. 1(a)]. Strictly speaking, a rigorous discussion of the SQM requires the introduction of the so-called “weak solution” concept, which makes it possible to correctly define the solution in the cases, in which the ray mapping is discontinuous on a set of measure zero [9]. However, for the sake of simplicity, we will concentrate on the practical side of the method. Let the prescribed continuous irradiance distribution $L(\mathbf {y})\in D$ be approximated by a discrete distribution $L_{\textrm {d}} (\mathbf {y})=\sum _{i=1}^{N} L_i \delta (\mathbf {y}-\mathbf {y}_i)$ (a sum of Dirac delta functions) concentrated at $N$ points $\mathbf {y}_i$, $i=1, \ldots , N$. In this case, the function $u_2 (\mathbf {y})$ in Eqs. (10) and (11) is replaced by a set of $N$ numbers $u_{2, i}$, $i=1, \ldots , N$. It is these values that are varied in the method. In this case, refractive surface $R_1$ is a piecewise-smooth surface consisting of segments of refractive lenses of Eq. (6) (ellipsoids or hyperboloids) focusing the incident plane beam to the points $(\mathbf {y}_i, f-u_{2, i})$, $i=1, \ldots , N$. This surface is schematically shown in Fig. 1(b) with bold segments of the colored curves showing the ellipsoids focusing the incident light to the individual points of the second optical surface. According to Eq. (10), the function $u_1 (\mathbf {x})$ defining the surface $R_1$ has the form

$$u_1 (\mathbf{x})=\mathop{\min}_i \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y}_i)-u_{2, i} \right].$$
For the fixed values $u_{2, i}$, $i=1, \ldots , N$, it is possible to calculate the discrete energy distribution generated at the points $\left (\mathbf {y}_i, f-u_{2, i} \right )$ by the surface defined by the function $u_1 (\mathbf {x})$ of Eq. (13). We do this by calculating the domains $C(\mathbf {y}_i)\subset G$, from which the light is focused to the points $u_{2, i}$. According to Eqs. (11) and (13), these domains are defined by the following condition:
$$C(\mathbf{y}_i)=\left\{\mathbf{x}\in G\; |\; {\mathcal{K}}(\mathbf{x}, \mathbf{y}_i)-u_{2, i} \le {\mathcal{K}}(\mathbf{x}, \mathbf{y}_j)-u_{2, j} ,\; \; \; \forall j \in \{1,\ldots,N\}\right\}.$$
By integrating the incident irradiance distribution over this domain, we find the light flux focused to the point $u_{2, i}$:
$$L_{\textrm{est}, i} =\int_{C(\mathbf{y}_i)} I(\mathbf{x})\mathrm{d}{\mathbf{x}}.$$

In practical calculations, the integral in the right-hand side of Eq. (15) can be calculated numerically using the ray tracing method, which is schematically shown in Fig. 1(b). In this approach, the incident beam is approximated by a set of $N_\textrm {tr}$ rays going from the nodes of a certain grid $\mathbf {x}_j$, $j=1, \ldots , N_\textrm {tr}$ in the plane $z=0$. These rays “carry” the light flux values $I_j$, $j=1, \ldots , N_\textrm {tr}$, which are proportional to the values of the irradiance of the incident beam at the grid nodes. According to Eq. (13), for each ray coming from a certain point $\mathbf {x}_j$, we find an index $i(j)=\mathop {\arg \min }\nolimits _i \left [{\mathcal {K}}(\mathbf {x}_j, \mathbf {y}_i)-u_{2,i} \right ]$ corresponding to the lens, which the ray intersects first [Fig. 1(b)]. This ray is directed to the corresponding point $\left (\mathbf {y}_{i(j)}, f-u_{2,i(j)} \right )$. Then, at the points of the second optical surface, the light flux values are calculated, which are equal to the sums of the light flux values of the rays focused to these points. In Fig. 1(b), these values are shown by numbers near each of the points of the “discretized” second optical surface. It is easy to see that with an increase in the $u_{2, i}$ value (in the case of fixed $u_{2, j}, j\ne i$), the domain $C(\mathbf {y}_i)$ increases, so that the light flux $L_{\textrm {est}, i}$ increases as well. Thus, if at some point $u_{2, i}$, the required light flux is greater than the obtained one [$L_{\textrm {d}} (\mathbf {y}_i)=L_i >L_{\textrm {est}, i}$], it is necessary to increase the $u_{2i}$ value (and to decrease it in the opposite case). It is natural to choose the magnitude of this increase or decrease proportional to the difference between the required and calculated light flux values $\left (L_i -L_{\textrm {est}, i} \right )$ [20]. In this way, the values $u_{2, i}$ can be found iteratively. At each iteration, the generated distribution $L_{\textrm {est}, i}$, $i=1, \ldots , N$ is calculated, and the $u_{2, i}$ values are updated according to the rule

$$u_{2, i} \Rightarrow u_{2, i} + \alpha \cdot \left(L_i -L_{\textrm{est}, i} \right)=u_{2, i} + \alpha \cdot \left(L_i -\int_{C(\mathbf{y}_i)} I(\mathbf{x})\mathrm{d}{\mathbf{x}}\right), i=1, \ldots, N,$$
where $\alpha > 0$ is the algorithm step size. The described method is a version of the SQM and enables obtaining an approximate solution of the problem. The accuracy of this solution is determined by the number of points $N$ used for approximating the prescribed continuous irradiance distribution $L(\mathbf {y})\in D$ and for representing the surface $R_2$. In the following section, we show that the SQM of Eq. (16) is equivalent to the gradient descent method applied to a certain function associated with the mass transportation problem.

It is important to note that according to the light beam reversibility principle, the considered PCBS can be formulated in a “reverse” form. In this case, the propagation direction is reversed and the input and required irradiance distributions are interchanged. In such a reverse formulation, one can discretize the irradiance distribution $I(\mathbf {x})$ [instead of $L(\mathbf {y})$] and then use an equation similar to Eq. (16) for the calculation of the values $u_{1,i}$ representing the first surface.

5. Connection of the SQM with a gradient method

5.1 Mass transportation problem

In Ref. [10], it was shown that in the considered collimated beam shaping problem, the calculation of a ray mapping satisfying the light flux conservation law (1) can be reduced to finding a solution to a Monge–Kantorovich mass transportation problem (MTP) with the cost function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$ of Eq. (5). The strict formulation of this fact is given by the following theorem:

Theorem 1. Let the mapping $\gamma : G\to D$ provide a minimum to the functional

$${\mathcal{F}}(T)=\int_G {\mathcal{K}}(\mathbf{x}, T(\mathbf{x}))I(\mathbf{x})\mathrm{d}{\mathbf{x}},$$
defined on the space of bijective mappings $T:G\to D$ satisfying the light flux conservation law
$$I(\mathbf{x})=L(T(\mathbf{x}))\cdot J_{T} (\mathbf{x}).$$
In this case, $\gamma$ is the sought-for mapping for the system of two optical surfaces.

Proof. Let us consider the Lagrange functional for the considered variational problem [25]:

$${\mathcal{H}}(T, \mu) = \int_G \left\{{\mathcal{K}}\left(\mathbf{x}, T(\mathbf{x})\right)I(\mathbf{x})+\mu (\mathbf{x})\left[L(T(\mathbf{x}))\cdot J_{T} (\mathbf{x})-I(\mathbf{x})\right]\right\}\mathrm{d}{\mathbf{x}},$$
where $\mu (\mathbf {x})$ is an arbitrary continuous function (Lagrange multiplier). Due to the bijectivity of the mapping $T$, instead of the function $\mu (\mathbf {x})$ in Eq. (19), we can write $\lambda \left (T(\mathbf {x})\right )$, where $\lambda (\mathbf {y})$ is an arbitrary continuous function [26]. In this case, we rewrite functional (19) as
$${\mathcal{H}}(T, \lambda) = \int_G \left\{{\mathcal{K}}\left(\mathbf{x}, T(\mathbf{x})\right)I(\mathbf{x})+\lambda (T(\mathbf{x}))\left[L\left(T(\mathbf{x})\right)\cdot J_{T} (\mathbf{x})-I(\mathbf{x})\right]\right\}\mathrm{d}{\mathbf{x}}.$$
The conditional maximum of the functional (17) under the constraint (18) is a critical point of the functional (20).

The variations of the functional ${\mathcal {H}}(T, \lambda )$ can be easily calculated. From Eq. (20), one can easily show that the condition $\delta _{\lambda } {\mathcal {H}}=0$ of zero variation with respect to $\lambda$ is equivalent to the light flux conservation law [26]. To calculate the variation $\delta _{T} {\mathcal {H}}$ with respect to $T$, let us split the integral in the right-hand side of Eq. (20) into two integrals and make the change of variables $\mathbf {y}=T(\mathbf {x})$ in the second integral:

$$\begin{aligned} {\mathcal{H}}(T, \lambda) & = \int_G \left[{\mathcal{K}}\left(\mathbf{x}, T(\mathbf{x})\right)-\lambda \left(T(\mathbf{x})\right)\right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} + \int_G \lambda \left(T(\mathbf{x})\right)L\left(T(\mathbf{x})\right)\cdot J_{T} (\mathbf{x})\mathrm{d}{\mathbf{x}} \\ & = \int_G \left[{\mathcal{K}}\left(\mathbf{x}, T(\mathbf{x})\right)-\lambda \left(T(\mathbf{x})\right)\right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} + \int_D \lambda(\mathbf{y})L(\mathbf{y})\mathrm{d}{\mathbf{y}}. \end{aligned}$$
From Eq. (21), one can obtain that the condition $\delta _{T} {\mathcal {H}}=0$ of zero variation with respect to $T$ is equivalent to the following pair of equations:
$$\left. \frac{\partial}{\partial y_i} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y})-\lambda (\mathbf{y})\right]\right|_{\mathbf{y}=T(\mathbf{x})} =0, \; \; \; \; \; i=1, 2.$$

The obtained equations are equivalent to the integrability condition (7b). According to the conditional extremum theorem [25], if the mapping $\gamma$ provides a conditional extremum to the problem of Eqs. (17) and (18), then there exists a function $\lambda (\mathbf {y})$, which satisfies Eqs. (22). This means that the mapping $\gamma$ is integrable, i. e. can be implemented by a system of two refracting surfaces described by the functions $u_1 (\mathbf {x})$ and $u_2 (\mathbf {y})$, where $u_2 (\mathbf {y})\equiv \lambda (\mathbf {y})$, and the function $u_1 (\mathbf {x})$ is defined through $u_2 (\mathbf {y})$ using Eq. (10). This concludes the proof of the theorem.

From the presented proof, it is evident that the solution of the collimated beam shaping problem can be obtained from any critical point of the Lagrange functional (20), and not only from the point corresponding to the minimum of the functional (17). In particular, the mapping $\gamma$, which provides a maximum to the functional (17), also gives a solution of the problem. From the form of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$ [see Eq. (5)], it is clear that the conditions of minimization and maximization of the functional (17) correspond to the Galilean and Keplerian configurations of rays, respectively [8,9].

Let us note that the ray mapping appears under the integral in the Lagrange functional (20). This makes it possible to use the MTP approach for calculating refracting surfaces implementing not only smooth mappings $\gamma$, but also mappings, which are discontinuous on a set of measure zero. In the case of a discontinuous mapping $\gamma$, the optical surface $u_1 (\mathbf {x})$ calculated using Eq. (10) will be piecewise-smooth.

5.2 Modification of the Lagrange functional

From the proof of Theorem 1, it follows that if a pair $(T, \lambda )$ is a critical point of the Lagrange functional ${\mathcal {H}} (T, \lambda )$, then the mapping $\mathbf {y}=T(\mathbf {x})$ satisfies the light flux conservation law and is defined by a system of refracting surfaces $u_1 (\mathbf {x})$ and $u_2 (\mathbf {y})$, where $u_2 (\mathbf {y})=\lambda (\mathbf {y})$, and $u_1 (\mathbf {x})$ is found from Eq. (10). At a critical point, the mapping $\mathbf {y}=T(\mathbf {x})$ and the function $\lambda (\mathbf {y})$ are related by the integrability condition (22). Therefore, the mapping can be expressed through the function $\lambda (\mathbf {y})$ using Eq. (22). Similarly to the discussion above, where we restricted our consideration to the envelopes defined by Eq. (10), in what follows, we consider the critical points of a special kind. Namely, we assume that the condition of the equality of the derivatives to zero in Eq. (22) corresponds not just to some critical point of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})-\lambda (\mathbf {y})$, but to a minimum point. In this case, the mapping $T$ can be written in a form similar to the expression for the ray mapping (11):

$$T^{\lambda} (\mathbf{x})=\mathop{\arg \min}_{\mathbf{y}\in D} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y})-\lambda (\mathbf{y})\right].$$
As a result, instead of the Lagrange functional ${\mathcal {H}}(T, \lambda )$ [Eq. (21)], we can consider the following functional depending only on the function $\lambda (\mathbf {y})$:
$${\mathcal{H}}(\lambda)=\int_G \left[{\mathcal{K}}\left(\mathbf{x}, T^{\lambda} (\mathbf{x})\right)-\lambda \left(T^{\lambda} (\mathbf{x})\right)\right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} + \int_{D} L(\mathbf{y})\lambda (\mathbf{y})\mathrm{d}{\mathbf{y}}.$$

5.3 Generalized Voronoi cells and gradient descent method

In this section, we establish the connection between the functional (24) and the supporting quadric method (16). For this, let us consider the form of the functional (24) when the continuous irradiance distribution $L(\mathbf {y})$ is approximated by a discrete distribution $L_{\textrm {d}} (\mathbf {y})=\sum _{i=1}^{N} L_i \delta \left (\mathbf {y}-\mathbf {y}_i \right )$ defined on $N$ points $\mathbf {y}_i$, $i=1, \ldots , N$. In this case, the function $\lambda (\mathbf {y})$ is replaced by a vector $\boldsymbol \lambda = (\lambda _1, \ldots , \lambda _N)$, where $\lambda _i =\lambda (\mathbf {y}_i)$, and the mapping $T^{\lambda }$ becomes a mapping of the domain $G$ to a discrete set of points $\mathbf {y}_i$, $i=1, \ldots , N$.

Definition. Let us define the generalized weighted Voronoi cells as the following domains $C^{\boldsymbol \lambda } (\mathbf {y}_i)$: a point $\mathbf {x}\in G$ belongs to the domain $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ if and only if

$${\mathcal{K}}(\mathbf{x}, \mathbf{y}_i )-\lambda_i \le {\mathcal{K}}(\mathbf{x}, \mathbf{y}_j)-\lambda_j, \; \; \; \; \; \forall j\in\{1, \ldots, N\}.$$
We refer to the domains $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ as generalized weighted Voronoi cells, since in the case of a quadratic cost function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})=\left \| \mathbf {y}-\mathbf {x}\right \|^{2}$, these domains would become conventional weighted Voronoi cells [27]. Note that the cell $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ coincides with the inverse image $\left (T^{\lambda } \right )^{-1} (\mathbf {y}_i)$. Accordingly, the partition of the domain $G$ into Voronoi cells is the partition of this domain into the inverse images of the points $\mathbf {y}_i$, $i=1, \ldots , N$ associated with the mapping $T^{\lambda }$. As an example, Fig. 2 shows the generalized Voronoi cells $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ associated with the points $\mathbf {y}_i$ lying on a square grid. The values $\lambda _i$ used for obtaining Fig. 2 were chosen to provide equal areas of the cells (and, hence, equal flux values $L_i$).

 figure: Fig. 2.

Fig. 2. Partition (a) of the domain $G$ into equal-area generalized Voronoi cells $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ associated with the points $\mathbf {y}_i$ (b).

Download Full Size | PPT Slide | PDF

In the case of a discrete distribution $L_{\textrm {d}} (\mathbf {y})$, the functional ${\mathcal {H}}(\lambda )$ becomes a function defined on ${\mathbb {R}}^{N}$:

$$h(\boldsymbol\lambda)=\sum_{i=1}^{N}\int_{C^{\boldsymbol\lambda} (\mathbf{y}_i)} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y}_i)-\lambda_i \right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} +\sum_{i=1}^{N}\lambda_i L_i,$$
where $\boldsymbol \lambda$ is the vector of weights of the generalized Voronoi cells. It is worth noting that a similar function was considered in [27,28] in the context of the optimal mass transportation problem in the case, in which the initial mass distribution was continuous, and the target mass distribution was given by a weighted sum of Dirac masses. In the mentioned papers, the calculation of the gradient of the function similar to Eq. (26) in the case of a quadratic function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$ was considered.

In the Appendix, we consider the calculation of the partial derivatives of the function (26) for an arbitrary function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$. According to the results presented in the Appendix, these derivatives have the form

$$\frac{\partial h(\boldsymbol\lambda)}{\partial \lambda_i} =L_i -\int_{C^{\boldsymbol\lambda} (\mathbf{y}_i)} I(\mathbf{x})\mathrm{d}{\mathbf{x}}, \, \, \, i=1, \ldots, N.$$
Moreover, it is easy to show that the function $h(\boldsymbol \lambda )$ is concave [27,29]. The proof of this fact is completely analogous to the proof in Ref. [27] carried out for the case of a quadratic function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$.

Note that since the function $\lambda (\mathbf {y})$ coincides with the function $u_2 (\mathbf {y})$ defining the surface $R_2$, then $u_{2,i}=\lambda _i$ and the generalized Voronoi cells $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ coincide with the domains $C (\mathbf {y}_i)$ [Eq. (14)], on which the SQM quadrics are defined, which focus the incident beam to the points $\mathbf {y}_i$, $i=1, \ldots , N$. The values $L_{\textrm {est}, i} = \int _{C^{\boldsymbol \lambda } (\mathbf {y}_i)} I(\mathbf {x})\mathrm {d}{\mathbf {x}}$ are equal to the light flux focused to the points $\mathbf {y}_i$, $i=1, \ldots , N$. According to Eq. (27), at a critical point of the function $h(\boldsymbol \lambda )$, the following conditions hold:

$$L_{\textrm{est}, i} =\int_{C^{\boldsymbol\lambda} (\mathbf{y}_i)} I(\mathbf{x})\mathrm{d}{\mathbf{x}} = L_i ,\, \, \, i=1, \ldots, N.$$
Thus, the problem of calculating the second optical surface is (in the discrete case) reduced to finding a critical point of the function $h(\boldsymbol \lambda )$ [Eq. (26)]. Due to the concavity of $h(\boldsymbol \lambda )$, this critical point is the maximum point. It is easy to see that Eq. (16) describing the rule for correcting the $z$ coordinates of the points of the second surface in the SQM corresponds to the gradient descent method of maximizing the function (26). It is important to note that since the function $h(\boldsymbol \lambda )$ is concave, any local maximum will also be the global one. Therefore, for the SQM of Eq. (16), there is no trapping near local extrema, which is typical for “conventional” gradient methods.

It is important to note that the above theoretical analysis is different from the theoretical results of papers [10,12], which, in the discrete case, allow one to reduce the PCBS to a linear programming problem (LPP). In contrast, in the present work we reduce the PCBS to an unconstrained optimization problem, which consists in maximizing the concave nonlinear function $h(\boldsymbol \lambda )$ being a discrete analogue of the modified Lagrange functional (24). This approach, in comparison with the LPP-based formulations in Refs. [10,12], turns out to be much simpler in terms of the computational complexity. Indeed, the LPP representing the PCBS has large dimensions and, in most cases, cannot be solved using the commonly available LPP solvers (e.g. SOPLEX, MOSEK, or linprog MATLAB routine). For example, if both surfaces of the optical element are defined on grids with the dimensions of $100\times 100$ containing $N=10000$ points each (in this case, the incident beam and the generated collimated beam are also defined on grids with the same dimensions), the number of constraints in the mentioned LPP will be equal to $N^{2}=10000^{2}=10^{8}$ [10]. The LPP with such a large number of constrains is extremely complex for the existing LPP solvers. Thus, the “direct” application of the LPP-based technique of Refs. [10,12] is possible only for rather coarse meshes (with the sizes of about $50\times 50$), which is insufficient for many practical problems. The design examples presented in the next section demonstrate that the proposed SQM version of Eq. (16) [the gradient descent method of maximizing the function (26)] enables designing high-performance optical elements generating collimated beams with required irradiance distributions defined on the grids with sizes $100\times 100$ and $150\times 150$.

6. Design examples

6.1 Generation of a “square-frame” beam with uniform irradiance

As a first example, let us consider the calculation of an optical element transforming a beam with a circular cross-section and uniform irradiance

$$I(\mathbf{x})=\frac{1}{\pi R^{2}},\;\;\mathbf{x}\in G=\left\{ (x_1, x_2) \,|\, x_1^{2} +x_2^{2} \le R^{2} \right\}$$
into a “square-frame” beam. This beam has uniform irradiance in a square-frame region
$$L(\mathbf{y})=\frac{1}{w_1^{2} -w_2^{2}},\;\;\mathbf{y}\in D=S_1 \backslash S_2,$$
where $S_i =\left \{ (x_1, x_2) \, |\, |x_1|\le w_i /2,\, |x_2|\le w_i /2\right \}$, $i=1,2$ are two squares with the sides $w_1$ and $w_2$ ($w_1 > w_2$). The calculation was carried out with the following parameters: radius of the incident beam $R=1\textrm {~mm}$; dimensions of the generated square-frame beam $w_1 =2\textrm {~mm}$ and $w_2 =1\textrm {~mm}$; optical path length $F=11\textrm {~mm}$ in the output plane $z=f=10\textrm {~mm}$; refractive indices of the element material and the surrounding medium $n_1 =1.5$ and $n_0 = 1$, respectively.

For the calculation of the optical surfaces, the required irradiance distribution of the output beam (30) was approximated by a discrete distribution defined on a square grid with spacing $\Delta _y = 0.0173\textrm {~mm}$. In this case, the discrete distribution contains $N=100^{2}$ points. For the calculation of the second surface (of the values $u_{2,i}$, $i=1,\ldots ,N$), we used the SQM version described in Section 4. Under this approach, at each iteration the light flux at the points $\left (\mathbf {y}_i, f-u_{2,i} \right )$, $i=1,\ldots ,N$ is calculated using the ray tracing technique, and then the values $u_{2,i}$ are updated using Eq. (16). In the ray tracing method, the incident beam was approximated by a set of $N_\textrm {tr} =2N=20000$ rays defined on a square grid in the plane $z=0$. In the considered example, 25000 iterations were performed. The total calculation time on a desktop computer (Intel Core i9-7940X CPU) was of about 6 minutes.

After the calculations, the second surface defined by a discrete set of values $u_{2,i}$, $i=1,\ldots ,N$ was approximated using the least-squares approach by a system of two-dimensional B-splines [30]. Then, from the calculated second surface, the first surface was reconstructed using Eq. (10). The obtained optical surfaces are shown in Fig. 3(a). It is interesting to note that in the vicinity of the origin the first surface has sharp bends (derivative discontinuities) and resembles a tetrahedral pyramid. This is due to the fact that the central part of the first surface maps the central part of the incident beam (the neighborhood of the point $\mathbf {x}=0$) to the inner contour of the square-frame. This mapping is discontinuous. The possibility to design piecewise-smooth surfaces implementing discontinuous ray mappings is an important advantage of the SQM. In this case, the existing methods based on finding a numerical solution to a partial differential equation of elliptic type [17] are either inapplicable or numerically unstable.

 figure: Fig. 3.

Fig. 3. Optical surfaces of the designed element (a) and generated normalized irradiance distributions in the planes $z=10\textrm {~mm}$ (b) and $z=30\textrm {~mm}$ (c). Red and blue lines show the central irradiance cross-sections.

Download Full Size | PPT Slide | PDF

Figures 3(b) and 3(c) show the calculated normalized irradiance distributions generated by the designed optical element in the output plane $z=f=10\textrm {~mm}$ and in the plane $z=30\textrm {~mm}$. Note that these distributions were calculated on a denser grid with the spacing $\Delta _{y} /4=0.0043\textrm {~mm}$ by tracing $10^{7}$ rays. The presented simulation results demonstrate the formation of a square-frame beam with the required dimensions and almost uniform irradiance. The normalized root-mean-square deviations of the obtained irradiance distributions from the uniform one amount to 3.3% and 3.8% at $z=10\textrm {~mm}$ and $z=30\textrm {~mm}$, respectively. Let us note that the dimensions of the square-frame beam at different distances from the optical element remain the same, therefore, the wavefront of the generated output beam is plane. Indeed, the calculated root-mean-square (RMS) deviation of the optical path length from the constant design value $F=11\textrm {~mm}$ amounts to $9\textrm {~nm}$ at $z=10\textrm {~mm}$.

6.2 Generation of an annular beam with non-uniform irradiance

As the second example, we designed an optical element transforming the incident beam (29) with the radius $R=1\textrm {~mm}$ into an annular beam with a non-uniform irradiance distribution

$$L(\rho, \phi)=1.5+\cos 4\phi ,\, \, \, 0.5\,\textrm{mm}\le \rho \le 1\,\textrm{mm},\, \, 0\le \phi <2\pi,$$
where $(\rho , \phi )$ are the polar coordinates in the output plane $z=f$. The design was performed with the parameters $F=11\textrm {~mm}$, $f=10\textrm {~mm}$, $n_1 =1.5$, and $n_0 = 1$, coinciding with the parameters of the previous example. The process of calculating the optical surfaces was completely similar to the previous example. For the calculations, the irradiance distribution (31) was approximated by a discrete distribution consisting of $N=150^{2}$ points. Note that in contrast to the previous example, the required flux values $L_i$, $i=1,\ldots ,N$ [see Eq. (16)] are no longer constant but calculated from Eq. (31). In the ray tracing procedure, $N_\textrm {tr} =25\cdot N$ rays were used at each iteration. The utilization of greater values of $N$ and $N_\textrm {tr}$, as compared to the previous example, is due to the more complex form of the generated distribution. As in the previous example, 25000 iterations were performed. The total computation time amounted to approximately 2.5 hours. The significant increase in the computation time, as compared to the first example, is mainly due to approximately 28-fold increase in the number of rays traced at each iteration.

The calculated surfaces of the optical element are shown in Fig. 4(a). As in the first example, the first surface has sharp bends near the origin of coordinates, which is due to the discontinuity of the ray mapping, which maps the neighborhood of the point $\mathbf {x}=0$ to the inner boundary of the ring. Figures 4(b) and 4(c) show the simulated normalized irradiance distributions generated by the designed element in the output plane $z=f=10\textrm {~mm}$ and in the plane $z=30\textrm {~mm}$. As in the previous example, the distributions were calculated on a grid with $0.0043\textrm {~mm}$ spacing by tracing $10^{7}$ rays. The presented results demonstrate the formation of an annular beam with the required dimensions. The normalized root-mean-square deviations of the calculated irradiance distributions from the required distribution defined by Eq. (31) amount to 6.2% and 6.7% at $z=10\textrm {~mm}$ and $z=30\textrm {~mm}$, respectively. The wavefront of the generated beam is almost plane and the calculated root-mean-square deviation of the optical path length from the constant design value $F=11\textrm {~mm}$ amounts to $13\textrm {~nm}$ at $z=10\textrm {~mm}$.

 figure: Fig. 4.

Fig. 4. Optical surfaces of the designed element (a) and generated normalized irradiance distributions in the planes $z=10\textrm {~mm}$ (b) and $z=30\textrm {~mm}$ (c). Red and blue lines show the central irradiance cross-sections.

Download Full Size | PPT Slide | PDF

6.3 Discussion of the convergence criteria

The SQM convergence is usually controlled by using criteria that represent the difference between the calculated and prescribed discrete distributions (e. g., the RMS deviation) [15,20]. However, such a choice is not optimal in the problem under consideration. Indeed, according to the results presented in Section 5.3, the SQM of Eq. (16) is a gradient method of maximizing a concave function $h(\boldsymbol \lambda )$ [Eq. (26)], which is a discrete counterpart of the Lagrange functional. To illustrate the concavity of this function, we chose two random vectors, $\boldsymbol \lambda _1$ and $\boldsymbol \lambda _2$, and plotted the function along a line segment connecting them: $h(t \boldsymbol \lambda _1 + (1-t)\boldsymbol \lambda _2)$, $t\in [0, 1]$ [see the upper panel in Fig. 5(a)]. For comparison, the lower panel of Fig. 5(a) shows the RMS deviation between the vectors consisting of the prescribed, $\mathbf {L}$, and the calculated, $\mathbf {L}_\textrm {est}(t \boldsymbol \lambda _1 + (1-t)\boldsymbol \lambda _2)$, light flux values. The presence of several local minima and noisy behavior of the RMS deviation make it a poorer criterion for controlling the convergence of the iterative process than $h(\boldsymbol \lambda )$ or a certain monotonic function of it.

 figure: Fig. 5.

Fig. 5. (a) Convergence criteria calculated at $\boldsymbol \lambda = t \boldsymbol \lambda _1 + (1-t)\boldsymbol \lambda _2$, $t\in [0, 1]$: $h(\boldsymbol \lambda )$ (upper panel) and RMS of the difference between the prescribed and calculated distributions (lower panel). (b) Convergence rate of the method: the difference $h_\textrm {opt} - h(\boldsymbol \lambda )$ (upper panel) and the RMS of the difference between the prescribed and calculated distributions (lower panel) vs. the iteration number.

Download Full Size | PPT Slide | PDF

Let us illustrate this fact using as an example the optical element generating a square-frame beam (Fig. 3). The upper panel of Fig. 5(b) shows the plot of the quantity $h_\textrm {opt} - h(\boldsymbol \lambda )$, where $h(\boldsymbol \lambda )$ is calculated using Eq. (26) at each SQM iteration at $\boldsymbol \lambda =(u_{2,1}, \ldots ,u_{2,N})$, and $h_\textrm {opt}$ is the final value of this function obtained at the last iteration. From Fig. 5(b), it is evident that the difference $h_\textrm {opt} - h(\boldsymbol \lambda )$ monotonically decreases. An abrupt decrease in this value at iteration 5001 is caused by a decrease in the step $\alpha$ in Eq. (16): the first 5000 iterations were performed with the step $\alpha =1$, and then the step was decreased to $\alpha = 0.1$. A similar convergence plot for the RMS criterion is shown in the lower panel of Fig. 5(b). In contrast to $h(\boldsymbol \lambda )$, the RMS deviation graph exhibits a counterintuitive increase starting from the 5001st iteration, and is much more noisy. Therefore, the RMS criterion is not suitable for controlling the convergence, and it is the $h(\boldsymbol \lambda )$ function that should be used for this purpose.

7. Conclusion

We proposed a version of the supporting quadric method for designing refractive optical elements for collimated beam shaping. We demonstrated that this method is, in its essence, the gradient method of maximizing a concave function, which is a discrete analogue of the Lagrange functional in a corresponding mass transportation problem. We established that any maximum of this function provides a solution to the problem of collimated beam shaping. Therefore, the proposed method does not suffer from the convergence to a local extremum, which is typical for the “conventional” gradient methods. Another important advantage of the method is simplicity, since its implementation is based on only two simple equations: Eqs. (13) and (16). It is important to note that instead of the simple gradient descent method of Eq. (16), one can use more sophisticated convex optimization algorithms (e. g. Nesterov’s accelerated gradient method) to obtain a higher convergence rate.

As examples, we designed refractive optical elements, which transform a circular incident beam into a square-frame beam with uniform illuminance and to an annular beam with non-uniform illuminance. The presented examples demonstrate the possibility of designing optical elements with continuous piecewise-smooth surfaces and show high performance of the proposed method.

The authors believe that the connection between the SQM and the gradient method of maximizing a function related to a certain MTP can also be obtained for other problems of nonimaging optics, which can be formulated as an MTP. In particular, similar versions of the SQM can be formulated for the problems of calculating mirrors and refracting surfaces generating prescribed irradiance distributions in the far field. The differences will only appear when defining the generalized weighted Voronoi cells, the form of which is determined by the cost function, depending on the problem of nonimaging optics being solved.

Appendix

Let us consider the calculation of the derivatives of the function $h(\boldsymbol \lambda )$ in Eq. (26). For the sake of convenience, let us introduce the function

$$M(\mathbf{x}, \boldsymbol\lambda) = \mathop{\min}_{i\in \{1,\ldots,N\}} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y}_i)-\lambda_i \right].$$
Then, using the definition of the Voronoi cells, the first sum in Eq. (26) can be rewritten as
$$\Psi (\boldsymbol\lambda)=\sum_{i=1}^{N}\int_{C^{\boldsymbol\lambda} (\mathbf{y}_i)} \left[{\mathcal{K}}(\mathbf{x}, \mathbf{y}_i)-\lambda_i \right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} =\int_G M(\mathbf{x}, \boldsymbol\lambda)I(\mathbf{x})\mathrm{d}{\mathbf{x}}.$$

Let $\boldsymbol \delta _m =(0,\ldots ,0,\delta ,0,\ldots ,0)$ be a vector of dimension $N$ having a single ($m$-th) nonzero component. Let us consider the difference

$$\Psi \left(\boldsymbol\lambda +\boldsymbol\delta_m \right)-\Psi (\boldsymbol\lambda ) = \int_G \left[M\left(\mathbf{x},\boldsymbol\lambda + \boldsymbol\delta_m \right)-M(\mathbf{x}, \boldsymbol\lambda)\right]I(\mathbf{x})\mathrm{d}{\mathbf{x}}.$$
To be specific, let us assume that $\delta <0$ (the case of $\delta >0$ can be treated in a similar way). From Eq. (25), it follows that at $\delta <0$, the generalized Voronoi cell $C^{\boldsymbol \lambda } (\mathbf {y}_m )$ becomes smaller, and therefore, the following inclusion chain holds:
$$C^{\boldsymbol\lambda + \boldsymbol\delta_m} \left(\mathbf{y}_m \right)\subset C^{\boldsymbol\lambda} \left(\mathbf{y}_m \right)\subset G.$$

Let us divide the last integral in (34) into three integrals over the regions (35):

$$\begin{aligned} \Psi (\boldsymbol\lambda + \boldsymbol\delta_m )-\Psi (\boldsymbol\lambda) = & \int_{C^{\boldsymbol\lambda + \boldsymbol\delta_m} \left(\mathbf{y}_m \right)} & \left[M(\mathbf{x},\boldsymbol\lambda + \boldsymbol\delta_m)-M(\mathbf{x}, \boldsymbol\lambda)\right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} & \\ + & \int_{C^{\boldsymbol\lambda} (\mathbf{y}_m) \backslash C^{\boldsymbol\lambda + \boldsymbol\delta_m} (\mathbf{y}_m)} & \left[M(\mathbf{x},\boldsymbol\lambda + \boldsymbol\delta_m )-M(\mathbf{x}, \boldsymbol\lambda)\right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} & \\ + & \int_{G\backslash C^{\boldsymbol\lambda} (\mathbf{y}_m)} & \left[M( \mathbf{x}, \boldsymbol\lambda + \boldsymbol\delta_m ) - M(\mathbf{x}, \boldsymbol\lambda)\right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} & . \end{aligned}$$
Now, we consider each integral separately, starting with the first one.

Since $C^{\boldsymbol \lambda + \boldsymbol \delta } (\mathbf {y}_m ) \subset C^{\boldsymbol \lambda } (\mathbf {y}_m)$, both minima defining the values $M(\mathbf {x},\boldsymbol \lambda + \boldsymbol \delta _m )$ and $M(\mathbf {x}, \boldsymbol \lambda )$ [see Eq. (32)] are reached at the point $\mathbf {y}_m$:

$$\begin{aligned} M(\mathbf{x}, \boldsymbol\lambda + \boldsymbol\delta_m ) & = {\mathcal{K}}(\mathbf{x},\mathbf{y}_m )-\left(\lambda_m +\delta \right), \\ M(\mathbf{x}, \boldsymbol\lambda) & = {\mathcal{K}}\left(\mathbf{x},\mathbf{y}_m \right)-\lambda_m. \end{aligned}$$
Then, we get
$$\int_{C^{\boldsymbol\lambda + \boldsymbol\delta_m} \left(\mathbf{y}_m \right)}\left[M(\mathbf{x}, \boldsymbol\lambda + \boldsymbol\delta_m )-M(\mathbf{x}, \boldsymbol\lambda)\right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} ={-}\delta \cdot \int_{C^{\boldsymbol\lambda + \boldsymbol\delta_m} (\mathbf{y}_m)} I(\mathbf{x})\mathrm{d}{\mathbf{x}}.$$

For the second integral in (36), the minima defining the values $M(\mathbf {x},\boldsymbol \lambda + \boldsymbol \delta _m)$ and $M(\mathbf {x},\boldsymbol \lambda )$ are reached at different points $\mathbf {y}_i$. Nevertheless, for every point $\mathbf {x} \in C^{\boldsymbol \lambda } (\mathbf {y}_m ) \backslash C^{\boldsymbol \lambda + \boldsymbol \delta _m} (\mathbf {y}_m )$, the weighted distance ${\mathcal {K}}(\mathbf {x}, \mathbf {y}_i)-\lambda _i$ to an arbitrary point $\mathbf {y}_i$ either did not change, or changed by $|\delta |$. Accordingly, the minimum of these distances $M(\mathbf {x},\boldsymbol \lambda + \boldsymbol \delta _m )$ changed by no more than $|\delta |$. Taking into account the sign of $\delta$, we can write

$$0\le M(\mathbf{x}, \boldsymbol\lambda + \boldsymbol\delta_m )-M(\mathbf{x}, \boldsymbol\lambda )\le -\delta.$$
As a result, for the second integral we obtain the estimate
$$0\le \int_{C^{\boldsymbol\lambda} (\mathbf{y}_m ) \backslash C^{\boldsymbol\lambda + \boldsymbol\delta_m} (\mathbf{y}_m )}\left[M\left(\mathbf{x},\boldsymbol\lambda + \boldsymbol\delta_m \right)-M(\mathbf{x}, \boldsymbol\lambda)\right]I(\mathbf{x})\mathrm{d}{\mathbf{x}} \le - \delta \cdot \int_{C^{\boldsymbol\lambda} (\mathbf{y}_m ) \backslash C^{\boldsymbol\lambda + \boldsymbol\delta_m} (\mathbf{y}_m )} I(\mathbf{x})\mathrm{d}{\mathbf{x}}.$$
At $\delta \to 0$, the difference $C^{\boldsymbol \lambda } (\mathbf {y}_m ) \backslash C^{\boldsymbol \lambda + \boldsymbol \delta _m} (\mathbf {y}_m )$ tends to an empty set, and hence the integral in the right-hand side of Eq. (38) tends to zero. That is, we established that
$$\mathop{\lim}_{\delta \to 0} \left\{ \frac{1}{\delta} \cdot \int_{C^{\boldsymbol\lambda} (\mathbf{y}_m ) \backslash C^{\boldsymbol\lambda + \boldsymbol\delta_m} (\mathbf{y}_m )}\left[M(\mathbf{x}, \boldsymbol\lambda + \boldsymbol\delta_m )-M(\mathbf{x},\boldsymbol\lambda )\right]I(\mathbf{x}) \mathrm{d}{\mathbf{x}} \right\}=0.$$

The third integral equals zero, because for any point $\mathbf {x}$ from the integration domain both minima $M(\mathbf {x},\boldsymbol \lambda + \boldsymbol \delta _m)$ and $M(\mathbf {x},\boldsymbol \lambda )$ are reached at the same $\mathbf {y}_i ,\, \, i\ne m$ and, consequently, their difference is zero.

Finally, using Eqs. (36), (37) and (39), one can easily obtain the derivatives of the function $h(\boldsymbol \lambda )$ [Eq. (26)] in the form of Eq. (27).

Funding

Ministry of Science and Higher Education of the Russian Federation (State assignment to FSRC "Crystallography and Photonics" RAS); Russian Foundation for Basic Research (18-07-00982, 18-29-03067).

Acknowledgments

The development of the design method and the investigation of the designed elements was supported by Russian Foundation for Basic Research, and the implementation of the software for the optical element design was supported by Ministry of Science and Higher Education of the Russian Federation.

Disclosures

The authors declare no conflicts of interest.

References

1. Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014). [CrossRef]  

2. N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Computation of double freeform optical surfaces using a Monge–Ampère solver: Application to beam shaping,” Opt. Commun. 439, 251–259 (2019). [CrossRef]  

3. S. Chang, R. Wu, A. Li, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge–Ampère type equation,” J. Opt. 18(12), 125602 (2016). [CrossRef]  

4. C. Bösel, N. G. Worku, and H. Gross, “Ray-mapping approach in double freeform surface design for collimated beam shaping beyond the paraxial approximation,” Appl. Opt. 56(13), 3679–3688 (2017). [CrossRef]  

5. X. Mao, J. Li, F. Wang, R. Gao, X. Li, and Y. Xie, “Fast design method of smooth freeform lens with an arbitrary aperture for collimated beam shaping,” Appl. Opt. 58(10), 2512–2521 (2019). [CrossRef]  

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

7. Z. Feng, B. D. Froese, C.-Y. Huang, D. Ma, and R. Liang, “Creating unconventional geometric beams with large depth of field using double freeform-surface optics,” Appl. Opt. 54(20), 6277–6281 (2015). [CrossRef]  

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

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

10. J. Rubinstein and G. Wolansky, “Intensity control with a free-form lens,” J. Opt. Soc. Am. A 24(2), 463–469 (2007). [CrossRef]  

11. T. Glimm and V. Oliker, “Optical design of two-reflector systems, the Monge–Kantorovich mass transfer problem and Fermat’s principle,” J. Math. Mech. 53(5), 1255–1278 (2004). [CrossRef]  

12. V. Oliker, “Designing freeform lenses for intensity and phase control of coherent light with help from geometry and mass transport,” Arch. Ration. Mech. Anal. 201(3), 1013–1045 (2011). [CrossRef]  

13. J. Rubinstein, G. Wolansky, and Y. Weinberg, “Ray mappings and the weighted least action principle,” J. Math. Ind. 8(1), 6 (2018). [CrossRef]  

14. L. A. Caffarelli, S. A. Kochengin, and V. I. Oliker, “On the numerical solution of the problem of reflector design with given far-field scattering data,” Ser. Contemp. Appl. Math. 226, 13–32 (1999). [CrossRef]  

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

16. V. 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).

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

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

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

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

21. L. L. Doskolovich, A. Y. Dmitriev, E. A. Bezus, and M. A. Moiseev, “Analytical design of freeform optical elements generating an arbitrary-shape curve,” Appl. Opt. 52(12), 2521–2526 (2013). [CrossRef]  

22. L. P. Eisenhart, A Treatise on the Differential Geometry of Curves and Surfaces (Schwarz, 2008).

23. D. A. Bykov, L. L. Doskolovich, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation problem in the design of freeform optical elements generating far-field irradiance distributions for plane incident beam,” Appl. Opt. 58(33), 9131–9140 (2019). [CrossRef]  

24. L. L. Doskolovich, A. A. Mingazov, D. A. Bykov, and E. A. Bezus, “Formulation of the inverse problem of calculating the optical surface for an illuminating beam with a plane wavefront as the Monge–Kantorovich problem,” Comput. Opt. 43(5), 705–713 (2019). [CrossRef]  

25. 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, 1999).

26. L. L. Doskolovich, D. A. Bykov, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions,” Opt. Express 27(9), 13083–13097 (2019). [CrossRef]  

27. Q. Mérigot, “A multiscale approach to optimal transport,” Comput. Graph. Forum. 30(5), 1583–1592 (2011). [CrossRef]  

28. F. Aurenhammer, F. Hoffmann, and B. Aronov, “Minkowski-type theorems and least-squares clustering,” Algorithmica 20(1), 61–76 (1998). [CrossRef]  

29. G. Wolansky, “On semi-discrete Monge–Kantorovich and generalized partitions,” J. Optim. Theory App. 165(2), 359–384 (2015). [CrossRef]  

30. C. de Boor, A Practical Guide to Splines (Springer-Verlag, 2001).

References

  • View by:

  1. Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014).
    [Crossref]
  2. N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Computation of double freeform optical surfaces using a Monge–Ampère solver: Application to beam shaping,” Opt. Commun. 439, 251–259 (2019).
    [Crossref]
  3. S. Chang, R. Wu, A. Li, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge–Ampère type equation,” J. Opt. 18(12), 125602 (2016).
    [Crossref]
  4. C. Bösel, N. G. Worku, and H. Gross, “Ray-mapping approach in double freeform surface design for collimated beam shaping beyond the paraxial approximation,” Appl. Opt. 56(13), 3679–3688 (2017).
    [Crossref]
  5. X. Mao, J. Li, F. Wang, R. Gao, X. Li, and Y. Xie, “Fast design method of smooth freeform lens with an arbitrary aperture for collimated beam shaping,” Appl. Opt. 58(10), 2512–2521 (2019).
    [Crossref]
  6. 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]
  7. Z. Feng, B. D. Froese, C.-Y. Huang, D. Ma, and R. Liang, “Creating unconventional geometric beams with large depth of field using double freeform-surface optics,” Appl. Opt. 54(20), 6277–6281 (2015).
    [Crossref]
  8. 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]
  9. 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]
  10. J. Rubinstein and G. Wolansky, “Intensity control with a free-form lens,” J. Opt. Soc. Am. A 24(2), 463–469 (2007).
    [Crossref]
  11. T. Glimm and V. Oliker, “Optical design of two-reflector systems, the Monge–Kantorovich mass transfer problem and Fermat’s principle,” J. Math. Mech. 53(5), 1255–1278 (2004).
    [Crossref]
  12. V. Oliker, “Designing freeform lenses for intensity and phase control of coherent light with help from geometry and mass transport,” Arch. Ration. Mech. Anal. 201(3), 1013–1045 (2011).
    [Crossref]
  13. J. Rubinstein, G. Wolansky, and Y. Weinberg, “Ray mappings and the weighted least action principle,” J. Math. Ind. 8(1), 6 (2018).
    [Crossref]
  14. L. A. Caffarelli, S. A. Kochengin, and V. I. Oliker, “On the numerical solution of the problem of reflector design with given far-field scattering data,” Ser. Contemp. Appl. Math. 226, 13–32 (1999).
    [Crossref]
  15. S. Kochengin and V. Oliker, “Computational algorithms for constructing reflectors,” Comput. Visualization Sci. 6(1), 15–21 (2003).
    [Crossref]
  16. V. 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).
  17. 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]
  18. 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]
  19. V. Oliker, “Controlling light with freeform multifocal lens designed with supporting quadric method (SQM),” Opt. Express 25(4), A58–A72 (2017).
    [Crossref]
  20. 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]
  21. L. L. Doskolovich, A. Y. Dmitriev, E. A. Bezus, and M. A. Moiseev, “Analytical design of freeform optical elements generating an arbitrary-shape curve,” Appl. Opt. 52(12), 2521–2526 (2013).
    [Crossref]
  22. L. P. Eisenhart, A Treatise on the Differential Geometry of Curves and Surfaces (Schwarz, 2008).
  23. D. A. Bykov, L. L. Doskolovich, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation problem in the design of freeform optical elements generating far-field irradiance distributions for plane incident beam,” Appl. Opt. 58(33), 9131–9140 (2019).
    [Crossref]
  24. L. L. Doskolovich, A. A. Mingazov, D. A. Bykov, and E. A. Bezus, “Formulation of the inverse problem of calculating the optical surface for an illuminating beam with a plane wavefront as the Monge–Kantorovich problem,” Comput. Opt. 43(5), 705–713 (2019).
    [Crossref]
  25. 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, 1999).
  26. L. L. Doskolovich, D. A. Bykov, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions,” Opt. Express 27(9), 13083–13097 (2019).
    [Crossref]
  27. Q. Mérigot, “A multiscale approach to optimal transport,” Comput. Graph. Forum. 30(5), 1583–1592 (2011).
    [Crossref]
  28. F. Aurenhammer, F. Hoffmann, and B. Aronov, “Minkowski-type theorems and least-squares clustering,” Algorithmica 20(1), 61–76 (1998).
    [Crossref]
  29. G. Wolansky, “On semi-discrete Monge–Kantorovich and generalized partitions,” J. Optim. Theory App. 165(2), 359–384 (2015).
    [Crossref]
  30. C. de Boor, A Practical Guide to Splines (Springer-Verlag, 2001).

2019 (5)

2018 (4)

2017 (2)

2016 (1)

S. Chang, R. Wu, A. Li, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge–Ampère type equation,” J. Opt. 18(12), 125602 (2016).
[Crossref]

2015 (3)

2014 (1)

Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014).
[Crossref]

2013 (1)

2011 (3)

Q. Mérigot, “A multiscale approach to optimal transport,” Comput. Graph. Forum. 30(5), 1583–1592 (2011).
[Crossref]

V. Oliker, “Designing freeform lenses for intensity and phase control of coherent light with help from geometry and mass transport,” Arch. Ration. Mech. Anal. 201(3), 1013–1045 (2011).
[Crossref]

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]

2010 (1)

2007 (1)

2004 (1)

T. Glimm and V. Oliker, “Optical design of two-reflector systems, the Monge–Kantorovich mass transfer problem and Fermat’s principle,” J. Math. Mech. 53(5), 1255–1278 (2004).
[Crossref]

2003 (1)

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

1999 (1)

L. A. Caffarelli, S. A. Kochengin, and V. I. Oliker, “On the numerical solution of the problem of reflector design with given far-field scattering data,” Ser. Contemp. Appl. Math. 226, 13–32 (1999).
[Crossref]

1998 (1)

F. Aurenhammer, F. Hoffmann, and B. Aronov, “Minkowski-type theorems and least-squares clustering,” Algorithmica 20(1), 61–76 (1998).
[Crossref]

Andreev, E. S.

Aronov, B.

F. Aurenhammer, F. Hoffmann, and B. Aronov, “Minkowski-type theorems and least-squares clustering,” Algorithmica 20(1), 61–76 (1998).
[Crossref]

Aurenhammer, F.

F. Aurenhammer, F. Hoffmann, and B. Aronov, “Minkowski-type theorems and least-squares clustering,” Algorithmica 20(1), 61–76 (1998).
[Crossref]

Bezus, E. A.

L. L. Doskolovich, A. A. Mingazov, D. A. Bykov, and E. A. Bezus, “Formulation of the inverse problem of calculating the optical surface for an illuminating beam with a plane wavefront as the Monge–Kantorovich problem,” Comput. Opt. 43(5), 705–713 (2019).
[Crossref]

L. L. Doskolovich, D. A. Bykov, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions,” Opt. Express 27(9), 13083–13097 (2019).
[Crossref]

D. A. Bykov, L. L. Doskolovich, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation problem in the design of freeform optical elements generating far-field irradiance distributions for plane incident beam,” Appl. Opt. 58(33), 9131–9140 (2019).
[Crossref]

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]

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]

L. L. Doskolovich, A. Y. Dmitriev, E. A. Bezus, and M. A. Moiseev, “Analytical design of freeform optical elements generating an arbitrary-shape curve,” Appl. Opt. 52(12), 2521–2526 (2013).
[Crossref]

Bösel, C.

Bräuer, A.

Bykov, D. A.

Caffarelli, L. A.

L. A. Caffarelli, S. A. Kochengin, and V. I. Oliker, “On the numerical solution of the problem of reflector design with given far-field scattering data,” Ser. Contemp. Appl. Math. 226, 13–32 (1999).
[Crossref]

Cassarly, W. J.

Chang, S.

S. Chang, R. Wu, A. Li, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge–Ampère type equation,” J. Opt. 18(12), 125602 (2016).
[Crossref]

de Boor, C.

C. de Boor, A Practical Guide to Splines (Springer-Verlag, 2001).

Dmitriev, A. Y.

Doskolovich, L. L.

D. A. Bykov, L. L. Doskolovich, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation problem in the design of freeform optical elements generating far-field irradiance distributions for plane incident beam,” Appl. Opt. 58(33), 9131–9140 (2019).
[Crossref]

L. L. Doskolovich, A. A. Mingazov, D. A. Bykov, and E. A. Bezus, “Formulation of the inverse problem of calculating the optical surface for an illuminating beam with a plane wavefront as the Monge–Kantorovich problem,” Comput. Opt. 43(5), 705–713 (2019).
[Crossref]

L. L. Doskolovich, D. A. Bykov, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions,” Opt. Express 27(9), 13083–13097 (2019).
[Crossref]

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]

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]

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]

L. L. Doskolovich, A. Y. Dmitriev, E. A. Bezus, and M. A. Moiseev, “Analytical design of freeform optical elements generating an arbitrary-shape curve,” Appl. Opt. 52(12), 2521–2526 (2013).
[Crossref]

Eisenhart, L. P.

L. P. Eisenhart, A Treatise on the Differential Geometry of Curves and Surfaces (Schwarz, 2008).

Evans, L. C.

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, 1999).

Feng, Z.

Fournier, F. R.

Froese, B. D.

Gao, R.

Glimm, T.

T. Glimm and V. Oliker, “Optical design of two-reflector systems, the Monge–Kantorovich mass transfer problem and Fermat’s principle,” J. Math. Mech. 53(5), 1255–1278 (2004).
[Crossref]

Gross, H.

Hoffmann, F.

F. Aurenhammer, F. Hoffmann, and B. Aronov, “Minkowski-type theorems and least-squares clustering,” Algorithmica 20(1), 61–76 (1998).
[Crossref]

Huang, C.-Y.

IJzerman, W. L.

N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Computation of double freeform optical surfaces using a Monge–Ampère solver: Application to beam shaping,” Opt. Commun. 439, 251–259 (2019).
[Crossref]

Kochengin, S.

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

Kochengin, S. A.

L. A. Caffarelli, S. A. Kochengin, and V. I. Oliker, “On the numerical solution of the problem of reflector design with given far-field scattering data,” Ser. Contemp. Appl. Math. 226, 13–32 (1999).
[Crossref]

Li, A.

S. Chang, R. Wu, A. Li, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge–Ampère type equation,” J. Opt. 18(12), 125602 (2016).
[Crossref]

Li, H.

Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014).
[Crossref]

Li, J.

Li, X.

Liang, R.

Liu, P.

Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014).
[Crossref]

Liu, X.

Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014).
[Crossref]

Ma, D.

Mao, X.

Mérigot, Q.

Q. Mérigot, “A multiscale approach to optimal transport,” Comput. Graph. Forum. 30(5), 1583–1592 (2011).
[Crossref]

Michaelis, D.

Mingazov, A. A.

Moiseev, M. A.

Oliker, V.

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]

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]

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

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]

V. Oliker, “Designing freeform lenses for intensity and phase control of coherent light with help from geometry and mass transport,” Arch. Ration. Mech. Anal. 201(3), 1013–1045 (2011).
[Crossref]

T. Glimm and V. Oliker, “Optical design of two-reflector systems, the Monge–Kantorovich mass transfer problem and Fermat’s principle,” J. Math. Mech. 53(5), 1255–1278 (2004).
[Crossref]

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

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

Oliker, V. I.

L. A. Caffarelli, S. A. Kochengin, and V. I. Oliker, “On the numerical solution of the problem of reflector design with given far-field scattering data,” Ser. Contemp. Appl. Math. 226, 13–32 (1999).
[Crossref]

Rolland, J. P.

Rubinstein, J.

J. Rubinstein, G. Wolansky, and Y. Weinberg, “Ray mappings and the weighted least action principle,” J. Math. Ind. 8(1), 6 (2018).
[Crossref]

J. Rubinstein and G. Wolansky, “Intensity control with a free-form lens,” J. Opt. Soc. Am. A 24(2), 463–469 (2007).
[Crossref]

Schreiber, P.

ten Thije Boonkkamp, J. H. M.

N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Computation of double freeform optical surfaces using a Monge–Ampère solver: Application to beam shaping,” Opt. Commun. 439, 251–259 (2019).
[Crossref]

Wang, F.

Weinberg, Y.

J. Rubinstein, G. Wolansky, and Y. Weinberg, “Ray mappings and the weighted least action principle,” J. Math. Ind. 8(1), 6 (2018).
[Crossref]

Wolansky, G.

J. Rubinstein, G. Wolansky, and Y. Weinberg, “Ray mappings and the weighted least action principle,” J. Math. Ind. 8(1), 6 (2018).
[Crossref]

G. Wolansky, “On semi-discrete Monge–Kantorovich and generalized partitions,” J. Optim. Theory App. 165(2), 359–384 (2015).
[Crossref]

J. Rubinstein and G. Wolansky, “Intensity control with a free-form lens,” J. Opt. Soc. Am. A 24(2), 463–469 (2007).
[Crossref]

Worku, N. G.

Wu, R.

S. Chang, R. Wu, A. Li, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge–Ampère type equation,” J. Opt. 18(12), 125602 (2016).
[Crossref]

Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014).
[Crossref]

Xie, Y.

Yadav, N. K.

N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Computation of double freeform optical surfaces using a Monge–Ampère solver: Application to beam shaping,” Opt. Commun. 439, 251–259 (2019).
[Crossref]

Zhang, Y.

Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014).
[Crossref]

Zheng, Z.

S. Chang, R. Wu, A. Li, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge–Ampère type equation,” J. Opt. 18(12), 125602 (2016).
[Crossref]

Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014).
[Crossref]

Algorithmica (1)

F. Aurenhammer, F. Hoffmann, and B. Aronov, “Minkowski-type theorems and least-squares clustering,” Algorithmica 20(1), 61–76 (1998).
[Crossref]

Appl. Opt. (5)

Arch. Ration. Mech. Anal. (1)

V. Oliker, “Designing freeform lenses for intensity and phase control of coherent light with help from geometry and mass transport,” Arch. Ration. Mech. Anal. 201(3), 1013–1045 (2011).
[Crossref]

Comput. Graph. Forum. (1)

Q. Mérigot, “A multiscale approach to optimal transport,” Comput. Graph. Forum. 30(5), 1583–1592 (2011).
[Crossref]

Comput. Opt. (1)

L. L. Doskolovich, A. A. Mingazov, D. A. Bykov, and E. A. Bezus, “Formulation of the inverse problem of calculating the optical surface for an illuminating beam with a plane wavefront as the Monge–Kantorovich problem,” Comput. Opt. 43(5), 705–713 (2019).
[Crossref]

Comput. Visualization Sci. (1)

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

J. Math. Ind. (1)

J. Rubinstein, G. Wolansky, and Y. Weinberg, “Ray mappings and the weighted least action principle,” J. Math. Ind. 8(1), 6 (2018).
[Crossref]

J. Math. Mech. (1)

T. Glimm and V. Oliker, “Optical design of two-reflector systems, the Monge–Kantorovich mass transfer problem and Fermat’s principle,” J. Math. Mech. 53(5), 1255–1278 (2004).
[Crossref]

J. Opt. (1)

S. Chang, R. Wu, A. Li, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge–Ampère type equation,” J. Opt. 18(12), 125602 (2016).
[Crossref]

J. Opt. Soc. Am. A (2)

J. Optim. Theory App. (1)

G. Wolansky, “On semi-discrete Monge–Kantorovich and generalized partitions,” J. Optim. Theory App. 165(2), 359–384 (2015).
[Crossref]

Opt. Commun. (2)

Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. 331, 297–305 (2014).
[Crossref]

N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Computation of double freeform optical surfaces using a Monge–Ampère solver: Application to beam shaping,” Opt. Commun. 439, 251–259 (2019).
[Crossref]

Opt. Express (6)

Opt. Lett. (1)

Ser. Contemp. Appl. Math. (1)

L. A. Caffarelli, S. A. Kochengin, and V. I. Oliker, “On the numerical solution of the problem of reflector design with given far-field scattering data,” Ser. Contemp. Appl. Math. 226, 13–32 (1999).
[Crossref]

Other (4)

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

C. de Boor, A Practical Guide to Splines (Springer-Verlag, 2001).

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, 1999).

L. P. Eisenhart, A Treatise on the Differential Geometry of Curves and Surfaces (Schwarz, 2008).

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (5)

Fig. 1.
Fig. 1. (a) Geometry of the problem of collimated beam shaping. (b) Representation of the first surface as the envelope of a family of surfaces focusing light to the points of the second surface (two-dimensional case). The envelope is shown with a dashed black curve. Solid colored curves show the cross-sections of different surfaces of the family.
Fig. 2.
Fig. 2. Partition (a) of the domain $G$ into equal-area generalized Voronoi cells $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ associated with the points $\mathbf {y}_i$ (b).
Fig. 3.
Fig. 3. Optical surfaces of the designed element (a) and generated normalized irradiance distributions in the planes $z=10\textrm {~mm}$ (b) and $z=30\textrm {~mm}$ (c). Red and blue lines show the central irradiance cross-sections.
Fig. 4.
Fig. 4. Optical surfaces of the designed element (a) and generated normalized irradiance distributions in the planes $z=10\textrm {~mm}$ (b) and $z=30\textrm {~mm}$ (c). Red and blue lines show the central irradiance cross-sections.
Fig. 5.
Fig. 5. (a) Convergence criteria calculated at $\boldsymbol \lambda = t \boldsymbol \lambda _1 + (1-t)\boldsymbol \lambda _2$ , $t\in [0, 1]$ : $h(\boldsymbol \lambda )$ (upper panel) and RMS of the difference between the prescribed and calculated distributions (lower panel). (b) Convergence rate of the method: the difference $h_\textrm {opt} - h(\boldsymbol \lambda )$ (upper panel) and the RMS of the difference between the prescribed and calculated distributions (lower panel) vs. the iteration number.

Equations (43)

Equations on this page are rendered with MathJax. Learn more.

I ( x ) = L ( γ ( x ) ) J γ ( x ) ,
γ 1 ( B ) I ( x ) d x = B L ( y ) d y ,
n 0 u 1 ( x ) + n 1 y x 2 + [ f u 1 ( x ) u 2 ( y ) ] 2 + n 0 u 2 ( y ) = F ,
u 1 ( x ) + u 2 ( y ) = K ( x , y ) ,
K ( x , y ) = n 1 2 f n 0 F n 1 2 n 0 2 n 1 | n 1 2 n 0 2 | ( F n 0 f ) 2 ( n 1 2 n 0 2 ) y x 2 .
u 1 y ( x ) = K ( x , y ) u 2 ( y ) .
u 1 ( x ) = K ( x , y ) u 2 ( y ) ,
y i [ K ( x , y ) u 2 ( y ) ] = 0 , i = 1 , 2.
u 2 x ( y ) = K ( x , y ) u 1 ( x ) ,
u 2 ( y ) = K ( x , y ) u 1 ( x ) ,
x i [ K ( x , y ) u 1 ( x ) ] = 0 , i = 1 , 2.
u 1 ( x ) = min y D [ K ( x , y ) u 2 ( y ) ] .
γ ( x ) = arg min y D [ K ( x , y ) u 2 ( y ) ] .
u 2 ( y ) = min x G [ K ( x , y ) u 1 ( x ) ] .
u 1 ( x ) = min i [ K ( x , y i ) u 2 , i ] .
C ( y i ) = { x G | K ( x , y i ) u 2 , i K ( x , y j ) u 2 , j , j { 1 , , N } } .
L est , i = C ( y i ) I ( x ) d x .
u 2 , i u 2 , i + α ( L i L est , i ) = u 2 , i + α ( L i C ( y i ) I ( x ) d x ) , i = 1 , , N ,
F ( T ) = G K ( x , T ( x ) ) I ( x ) d x ,
I ( x ) = L ( T ( x ) ) J T ( x ) .
H ( T , μ ) = G { K ( x , T ( x ) ) I ( x ) + μ ( x ) [ L ( T ( x ) ) J T ( x ) I ( x ) ] } d x ,
H ( T , λ ) = G { K ( x , T ( x ) ) I ( x ) + λ ( T ( x ) ) [ L ( T ( x ) ) J T ( x ) I ( x ) ] } d x .
H ( T , λ ) = G [ K ( x , T ( x ) ) λ ( T ( x ) ) ] I ( x ) d x + G λ ( T ( x ) ) L ( T ( x ) ) J T ( x ) d x = G [ K ( x , T ( x ) ) λ ( T ( x ) ) ] I ( x ) d x + D λ ( y ) L ( y ) d y .
y i [ K ( x , y ) λ ( y ) ] | y = T ( x ) = 0 , i = 1 , 2.
T λ ( x ) = arg min y D [ K ( x , y ) λ ( y ) ] .
H ( λ ) = G [ K ( x , T λ ( x ) ) λ ( T λ ( x ) ) ] I ( x ) d x + D L ( y ) λ ( y ) d y .
K ( x , y i ) λ i K ( x , y j ) λ j , j { 1 , , N } .
h ( λ ) = i = 1 N C λ ( y i ) [ K ( x , y i ) λ i ] I ( x ) d x + i = 1 N λ i L i ,
h ( λ ) λ i = L i C λ ( y i ) I ( x ) d x , i = 1 , , N .
L est , i = C λ ( y i ) I ( x ) d x = L i , i = 1 , , N .
I ( x ) = 1 π R 2 , x G = { ( x 1 , x 2 ) | x 1 2 + x 2 2 R 2 }
L ( y ) = 1 w 1 2 w 2 2 , y D = S 1 S 2 ,
L ( ρ , ϕ ) = 1.5 + cos 4 ϕ , 0.5 mm ρ 1 mm , 0 ϕ < 2 π ,
M ( x , λ ) = min i { 1 , , N } [ K ( x , y i ) λ i ] .
Ψ ( λ ) = i = 1 N C λ ( y i ) [ K ( x , y i ) λ i ] I ( x ) d x = G M ( x , λ ) I ( x ) d x .
Ψ ( λ + δ m ) Ψ ( λ ) = G [ M ( x , λ + δ m ) M ( x , λ ) ] I ( x ) d x .
C λ + δ m ( y m ) C λ ( y m ) G .
Ψ ( λ + δ m ) Ψ ( λ ) = C λ + δ m ( y m ) [ M ( x , λ + δ m ) M ( x , λ ) ] I ( x ) d x + C λ ( y m ) C λ + δ m ( y m ) [ M ( x , λ + δ m ) M ( x , λ ) ] I ( x ) d x + G C λ ( y m ) [ M ( x , λ + δ m ) M ( x , λ ) ] I ( x ) d x .
M ( x , λ + δ m ) = K ( x , y m ) ( λ m + δ ) , M ( x , λ ) = K ( x , y m ) λ m .
C λ + δ m ( y m ) [ M ( x , λ + δ m ) M ( x , λ ) ] I ( x ) d x = δ C λ + δ m ( y m ) I ( x ) d x .
0 M ( x , λ + δ m ) M ( x , λ ) δ .
0 C λ ( y m ) C λ + δ m ( y m ) [ M ( x , λ + δ m ) M ( x , λ ) ] I ( x ) d x δ C λ ( y m ) C λ + δ m ( y m ) I ( x ) d x .
lim δ 0 { 1 δ C λ ( y m ) C λ + δ m ( y m ) [ M ( x , λ + δ m ) M ( x , λ ) ] I ( x ) d x } = 0.

Metrics

Select as filters


Select Topics Cancel
© Copyright 2022 | Optica Publishing Group. All Rights Reserved