## Abstract

We propose a method for designing refractive optical elements for collimated beam shaping in the geometrical optics approximation. In this method, the problem of finding a ray mapping is formulated as a linear assignment problem, which is a discrete version of the corresponding mass transportation problem. A method for reconstructing optical surfaces from a computed discrete ray mapping is proposed. The method is suitable for designing continuous piecewise-smooth optical surfaces. The design of refractive optical elements transforming beams with circular cross-section to variously shaped (rectangular, triangular, and cross-shaped) beams with plane wavefront is discussed. The presented numerical simulation results confirm high efficiency of the designed optical elements.

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

## 1. Introduction

The problem of collimated beam shaping consists in the design of a refractive optical element transforming a given incident beam with a plane wavefront into an output beam with a prescribed irradiance distribution and a plane wavefront. In the general case, an optical element that performs this transformation has two working surfaces. As a rule, the problem under study is addressed within geometrical optics and has a wide range of applications, such as illumination, laser beam processing and optical lithography. The problem is easy to solve only in the axisymmetric case. In this case, the problem is reduced to a single spatial dimension, enabling the optical surfaces to be calculated by solving ordinary differential equations [1–6]. In the general two-dimensional case, the design of an optical element is essentially more challenging.

In [7–9], the problem of collimated beam shaping was reduced to finding a solution to an elliptic nonlinear partial differential equation (PDE). The PDE is solved using finite-difference methods, in which the derivatives are replaced by their finite-difference approximations. Thus, the solution of the PDE is reduced to the solution of a system of nonlinear equations with respect to the values of the sought-for function, which describes the first optical surface and is defined on a certain grid. The resulting system of equations is solved using Newton’s method. The second optical surface is then reconstructed analytically from the first surface by applying the requirement of a constant optical path length. General shortcomings of this approach include high computational complexity of solving a system of nonlinear equations as well as the choice of the initial approximation to the solution of the system and the formulation of the boundary conditions.

The considered problem can also be addressed by using another, mathematically less rigorous approach based on the so-called ray-mapping methods [10–12]. In this case, the ray mapping, i.e. the mapping between the coordinates of the rays on the wavefronts of the incident and the output beams is found by solving a standard Monge–Ampère (MA) equation. The MA equation is solved either using the above-mentioned finite-difference method based on reducing the task to the solution of a system of nonlinear equations, or using an explicit finite-difference method based on the replacement of the original equation by the non-stationary parabolic Monge–Ampère equation [11,12]. The optical surfaces are calculated on the basis of the computed ray mapping using numerical integration. In addition to having high computational complexity, the described ray-mapping methods can offer only an approximate solution because the standard MA equation is valid only in the paraxial approximation [9,13].

In the basic theoretical works [14,15], the problem of collimated beam shaping was formulated as a Monge–Kantorovich mass transportation problem, and the corresponding cost functions were derived. It is important to note that these cost functions are non-quadratic, meaning that in this case the ray mapping is not described by a standard MA equation. In a discrete version, the mass transportation problem can be formulated as a linear assignment problem (LAP) [13,16]. Such an approach is utilized in a recent work by the present authors [13] for solving an inverse problem of calculating the eikonal of a light field.

In this work, we for the first time apply this LAP-based approach to a problem of collimated beam shaping. Due to the fact that the element under study has two working optical surfaces, the application of the LAP-based method for the solution of this problem has a number of specific features and differs from [13]. As examples, we design refractive optical elements transforming a plane circular beam into plane beams with uniform irradiance distribution in a rectangular and a triangular domains in a nonparaxial case. These examples demonstrate high efficiency of the proposed LAP-based approach: the computation time for square grids of 143 × 143 points on a standard PC is of about 3 minutes. In contrast, the calculation of such a mapping using finite-difference methods for the solution of a PDE of Monge–Ampère type requires solving a system of 143^{2} = 20449 nonlinear equations and is by several orders of magnitude slower [17]. Besides, the computation of the ray mapping by solving a mass transportation problem makes it possible to design elements with piecewise-smooth continuous optical surfaces, as distinct from smooth optical surfaces resulting from the solution of an elliptic nonlinear PDE [7–12]. It widens the range of the problems that can be solved, in particular, enabling the generation of uniform irradiance distributions in complex domains with non-smooth boundaries. As a demonstration of such a capability, we design an element with a piecewise-smooth optical surface transforming a circular beam into a cross-shaped beam.

## 2. Formulation of the problem

In this section, we formulate the problem within the geometrical optics framework using the associated notions of wavefronts, rays, and irradiances. The geometrical optics approximation implies that the characteristic dimensions of the optical surfaces and generated light distributions are much larger than the light wavelength.

Consider an incident beam with a plane wavefront perpendicular to the *z* axis and with the irradiance distribution *I*_{0}(**u**), **u** ∈ *G*, where **u** = (*u*_{1}, *u*_{2}) are the Cartesian coordinates in the plane *z* = 0 [Fig. 1(a)]. The beam is being transformed by a refractive optical element with the refractive index *n* > 1, which is surrounded by a medium with the refractive index *n*_{0} = 1. The element has two optical surfaces *f* and *g*. The problem to be solved consists in the design of the optical surfaces generating an output beam with a plane wavefront perpendicular to the *z* axis and a prescribed irradiance distribution *I*(**x**), **x** ∈ *D*, where **x** = (*x*_{1}, *x*_{2}) are the Cartesian coordinates in a certain plane *z* = *f*_{0} located after the optical element. Note that the function *I*(**x**) should satisfy the following normalization condition:

**u**(d

**x**) denote double integrals with respect to d

*u*

_{1}d

*u*

_{2}(d

*x*

_{1}d

*x*

_{2}). Besides, since the output beam has a plane wavefront, the irradiance

*I*(

**x)**will remain unchanged in any plane after the optical element [Fig. 1(a)]. Moreover, since both the input and output wavefronts are plane, the optical path

*L*will be the same for any ray propagating from the plane

*z*= 0 to the plane

*z*=

*f*

_{0}. It is convenient to define the value of

*L*in the following form:

*L*=

*f*

_{0}+

*δ*, where

*δ*is a change of the optical path length introduced by the optical element. The quantity

*δ*can be treated as a parameter of the problem. In particular, when the irradiance distributions

*I*(

**x**) and

*I*

_{0}(

**u**) are centrally symmetric, the central ray will propagate along the optical axis

*z*without refraction. In this case,

*δ*is easily expressed as

*δ*= (

*n*− 1)

*h*

_{0}, where

*h*

_{0}is the thickness of the optical element at

**u**= 0.

Let **x** = *T*(**u**) = (*x*_{1}(*u*_{1}, *u*_{2}), *x*_{2}(*u*_{1}, *u*_{2})) denote a ray mapping that represents the Cartesian coordinates **x** of a ray in the output plane *z* = *f*_{0} through the coordinates of the ray **u** in the plane *z* = 0 [Fig. 1(a)]. In the geometrical optics approximation, we can use the energy conservation law along the ray tubes. According to this law, the mapping *T*: *G* → *D* conserves the light flux. The light flux conservation law can be written in the following integral form [13,14]:

*ω*is any measurable subset of the domain

*G*.

Note that for a ray to remain parallel to the *z* axis after the refraction by the optical surfaces described by the functions *z* = *f*(**u**) and *z* = *g*(**x**), the normal vectors to the surfaces calculated at the points **u** and **x**(**u**) = *T*(**u**) have to be parallel. In this case, the optical surfaces can be locally considered as a plane-parallel plate, which shifts a ray without changing its propagation direction. Using the vector form of the refraction law, the partial derivatives of the functions *z* = *f*(**u**) and *z* = *g*(**x**) can be easily derived in the form [14]

**u**(

**x**) =

*T*

^{−}(

**x**) = (

*u*

_{1}(

*x*

_{1},

*x*

_{2}),

*u*

_{2}(

*x*

_{1},

*x*

_{2})) is the inverse ray mapping. According to Eq. (3), $\frac{\partial f(\mathbf{u})}{\partial {u}_{i}}=\frac{\partial g(\mathbf{x}(\mathbf{u}))}{\partial {x}_{i}}$,

*i*= 1, 2. Besides, Eqs. (3) suggest that not any mapping

**x**(

**u**) [

**u**(

**x**)] that conserves the energy in the sense of relationship (2) is a solution of the initial problem. Indeed, for the reconstruction of the functions

*z*=

*f*(

**u**) and

*z*=

*g*(

**x**) from Eqs. (3), the partial derivatives have to satisfy the following conditions:

It has been shown [14] that a mapping **x** = *T*(**u**) that corresponds to the solution of the problem in question satisfies the integrability conditions in Eq. (4) and provides an extremum of the following functional:

*T*:

*G*→

*D*. According to Eqs. (5) and (6), the mappings can be computed by solving a Monge–Kantorovich mass transportation problem. As the masses, the initial irradiance distribution

*I*

_{0}(

**u**),

**u**∈

*G*and the sought-for irradiance distribution

*I*(

**x**),

**x**∈

*D*are used. Let us note that the cost function defined by Eq. (6) differs from the one utilized in [14] by an additive constant that does not affect the extremum of the functional (5).

As mentioned above, the optical surfaces in the vicinity of the points **u**_{0} and **x**(**u**_{0}) locally correspond to a plane-parallel plate that shifts the ray going from the point **u**_{0} to the point **x**(**u**_{0}) in the output plane. It can be easily shown that the cost function in Eq. (6) turns to infinity when, at given parameters *δ* and *n* it is impossible to shift the ray by the magnitude |**s**_{0}| = |**u**_{0} − *T*(**u**_{0})|. If the cost function turns to infinity on a set with nonzero measure, there may be no solution to the problem. Note that with an increase in the parameter *γ* the achievable value of the shift also increases. In particular, we can increase the achievable shift by increasing *δ* and, hence, the thickness of the optical element.

One can show that a mapping **x** = *T*(**u**) that minimizes the functional (5) corresponds to a Galilean configuration of rays. The maximum of the functional corresponds to a mapping with a Keplerian configuration of rays.

In the paraxial approximation (i.e. under the condition $|\mathbf{x}(\mathbf{u})-\mathbf{u}{|}^{2}\ll {\gamma}^{2})C(\mathbf{s})\approx -\gamma +\frac{|\mathbf{s}{|}^{2}}{2\gamma}$, and the mapping **x**(**u**) = *T*(**u**) is computed by solving the mass transportation problem with a quadratic cost function. In this case, should there exist a differentiable mapping **x**(**u**), its computation is reduced to finding a solution to a standard Monge–Ampère equation [9–12].

It is important to note that the direct determination of the ray mapping from the solution of a mass transportation problem in Eqs. (5) and (6) makes it possible to operate with mappings that are differentiable not everywhere, but almost everywhere. In this case, the optical surfaces [functions *z* = *f* (**u**) and *z* = *g*(**x**)] are continuous and piecewise-smooth, as distinct from smooth surfaces obtained as the solution of elliptic PDEs [7–9]. This makes it possible to extend the range of problems that can be solved, in particular, to generate prescribed irradiance distributions on a zero-intensity background in domains with complex and non-smooth boundaries.

## 3. Reduction to the linear assignment problem

The discrete version of the mass transportation problem in Eqs. (5) and (6) can be formulated as a linear assignment problem (LAP) [16]. Indeed, let us assume that the domains *G* and *D* are divided into *N* equal-flux cells (or approximated by *N* cells), with the equality ${\int}_{{g}_{i}}{I}_{0}(\mathbf{u})\mathrm{d}\mathbf{u}}={\displaystyle {\int}_{{d}_{j}}I(\mathbf{x})\mathrm{d}\mathbf{x}=e$ being valid for any pair of cells *g* ⊂ *G*, *d _{j}* ⊂

*D*. For the sake of illustration, Fig. 1(b) shows approximations of a circular domain

*G*and a triangular domain

*D*by grids of

*N*= 196 square cells. While constructing these approximations, the grid steps were chosen so that

*N*points fall within the domains

*G*and

*D*. This approach based on the choice of the rectangular grid steps is suitable for generating uniform irradiance distributions in the domains of interest. In the case of nonuniform distributions

*I*

_{0}(

**u**) and

*I*(

**x**), an equal-flux set of

*N*cells in the domain

*G*(or

*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 [18–20]. In particular, an equal-flux set of

*N*=

*M*

_{1}

*M*

_{2}cells can be constructed in the following way. The domain

*G*can be divided by parallel lines

*u*

_{1}=

*u*

_{1,}

*= 1,*

_{i}, i*…, M*

_{1}into

*M*

_{1}equal-flux sub-domains. Afterwards, each sub-domain

*G*is further divided into

_{i}*M*

_{2}sub-domains by parallel lines

*u*

_{2}=

*u*

_{2,}

*= 1, …,*

_{j}, j*M*

_{2}.

For the grids (approximations) introduced in the domains *G* and *D*, all energy-conserving mappings *G* → *D* can be described by permutations (*i*_{1}, …, *i _{N}*) of

*N*numbers, which determine to which cells

*d*of the domain

_{i}*D*the cells

*g*

_{1}, …,

*g*are mapped. In terms of the LAP, the mapping of the cell

_{N}*g*to the cell

_{i}*d*can be interpreted as the assignment of the

_{j}*j*-th task to the

*i*-th agent. According to Eq. (6), the cost of the tasks is described by the matrix

**u**

*,*

_{i}**x**

*are the centers of the cells*

_{j}*g*and

_{i}*d*. Thus, the mapping $\mathbf{x}({\mathbf{u}}_{i})=\widehat{T}({\mathbf{u}}_{i})$ can be found by solving the following assignment problem:

_{j}*j*

_{1},

*…, j*). For the solution of the LAP, efficient polynomial algorithms have been proposed, including Hungarian algorithm [16], Jonker–Volgenant algorithm [21], and auction algorithm [22].

_{N}## 4. Reconstruction of the optical surface

Due to the discrete character of the LAP in Eq. (8), the resulting discrete mapping ${\mathbf{x}}_{i}=\widehat{T}({\mathbf{u}}_{i})$, *i* = 1, …, *N* describes the sought-for continuous mapping **x**(**u**) = *T*(**u**) with an error. Besides, if the domain *D* has a non-smooth boundary, the mapping *T*(**u**) may have discontinuities. In this case, the reconstruction of the optical surfaces by direct numerical integration of Eqs. (3) may result in a significant error. Because of this, in this section we propose a method for reconstructing optical surfaces from a computed discrete mapping. Within the proposed method, the second optical surface [function *g*(**x**)] is assumed to be smooth, while the first surface [function *f*(**u**)] is assumed to be continuous and piecewise-smooth. This case corresponds to the most challenging example of an optical element described in Subsection 5.3 and designed to generate a cross-shaped beam.

#### 4.1. Reconstruction of the upper optical surface

Let us first consider the calculation of the function *g*(**x**), which describes the second (upper) optical surface of the element. Because the function *g*(**x**) is assumed to be smooth, it can be approximated by two-dimensional B-splines [23].

Assume that the domain *D*, where *g*(**x**) is defined, is bounded by a rectangle $\overline{D}$ with the sides *w*_{1} and *w*_{2}. Let *x*_{1} = *w*_{1}· [*m*/(*N*_{1} − 1) − 1/2, *m* = 0, …, *N*_{1} − 1 be a set of *N*_{1} equidistant knots on the interval [−*w*_{1}/2, *w*_{1}/2] of the *x*_{1} axis. Using these knots, let us introduce basis splines (B-splines) of order *q*, denoting them by *B _{m}*(

*x*

_{1}),

*m*= 1, ….,

*N*

_{1}+

*q*− 2 [23]. In a similar way, let us introduce B-splines

*P*(

_{n}*x*

_{2}) for the second coordinate

*x*

_{2}∈ [−

*w*

_{2}/2,

*w*

_{2}/2]. The function

*g*(

**x**) =

*g*(

*x*

_{1},

*x*

_{2}) can be represented using the introduced basis functions as

*p*are the expansion coefficients. Note that if the splines

_{m,n}*B*(

_{m}*x*

_{1}) and

*P*(

_{n}*x*

_{2}) have the same order and the same number of knots, Eq. (9) contains

*N*= (

*N*

_{1}+

*q*− 2)

^{2}unknown coefficients

*p*.

_{m,n}The coefficients *p _{m,n}* are determined from the condition that the values of the derivatives of the spline

*g*(

**x**) at the points ${\mathbf{x}}_{i}=\widehat{T}({\mathbf{u}}_{i})$,

*i*= 1, …,

*N*constitute a least-squares approximation of the derivatives $\frac{\partial \widehat{g}}{\partial {x}_{k}}({\mathbf{x}}_{i})$. The latter are derived from Eq. (3) on the basis of a discrete mapping $\widehat{T}$, i.e. by assuming $\mathbf{u}({\mathbf{x}}_{i})={\widehat{T}}^{-1}({\mathbf{x}}_{i})$ in Eq. (3). Thus, the coefficients

*p*can be obtained by solving the following minimization problem:

_{m,n}**p**is the vector composed of the spline coefficients

*p*. It is easy to show that the solution of the problem (10) is reduced to a standard least-squares method. The proposed approach, which consists in approximating the derivatives of an optical surface, allows one to compensate for the errors resulting from computing the mapping by use of the ‘discrete’ LAP-based approach [24,25].

_{m,n}Note that for a nonrectangular domain $D\subset \overline{D}$, the solution of the problem (10) may be numerically unstable. In this case, one should apply a regularized least-squares method, in which the term $\lambda {\displaystyle {\sum}_{m,n}{p}_{m,n}^{2}}$, where *λ* is a small number, is added to the minimized function *S*(**p**).

#### 4.2. Reconstruction of the lower optical surface

Next, let us consider the calculation of the function *f*(**u**), which describes the first (lower) optical surface. As we mentioned earlier, the function *f*(**u**) is assumed to be piecewise-smooth, which does not allow us to approximate it by splines. We propose to reconstruct the function *f*(**u**) through the function *g*(**x**) using the representation of the optical surface from the supporting quadric method (SQM) [26–29]. The SQM is intended for designing a piecewise-smooth (refracting or reflecting) optical surface to focus the incident beam into a given array of points. In this case, the optical surface is composed of continuously stitched lens segments.

To apply the SQM to our problem, we assume that the second optical surface is defined on a sufficiently fine grid and is described by a set of points **X*** _{i}* = (

**x**

*(*

_{i}, g**x**

*)),*

_{i}*i*= 1, …,

*N*. Using the spline representation of the optical surface

_{g}*g*(

**x**) [Eq. (10)], the array of points

**X**

*can be calculated on a uniform grid with an arbitrarily small step. In the vicinity of each point*

_{i}**u**∈

*G*, the function

*f*(

**u**) corresponds to a lens segment that focuses the incident plane beam into one of the points

**X**

*on the upper optical surface.*

_{i}To derive the lens equation, let us use the condition of a constant optical path length from the points on the incident wavefront to the focal point [30,31]:

*(*

_{i}**u**) describes the lens surface and

*L*=

_{i}*L*− [

*f*

_{0}−

*g*(

**x**

*)] =*

_{i}*δ*+

*g*(

**x**

*) is the optical path length of the rays propagating from the plane*

_{i}*z*= 0 to the focal point

**X**

*= (*

_{i}**x**

*(*

_{i}, g**x**

*)). Equation (11) defines the lens surface*

_{i}*z*= Φ

*(*

_{i}**u**) in an implicit form. Through simple transformations, the function Φ

*(*

_{i}**u**) can be derived from Eq. (11) in the following explicit form:

It is easy to show that the lens surface is described by an ellipsoid of rotation with the rotation axis parallel to the *z* -axis, eccentricity *ε* = 1/*n*, focal parameter *p* = *δ*/*n*, and one of the foci located at the point **X*** _{i}* [30, 31].

According to the SQM, the surface that focuses into the set of points **X*** _{i}* consists of a set of ellipsoid segments defined by Eq. (12). At each point

**u**∈

*G*, the surface is calculated from the following expression [26–29]:

Thus, we propose to calculate the function *f*(**u**) using Eqs. (12)–(13) through the second optical surface defined on a set of points **X*** _{i}* = (

**x**

*(*

_{i}, g**x**

*)),*

_{i}*i*= 1, …,

*N*. Note that Eq. (13) is also convenient for calculating the mapping

_{g}*T*(

**u**) at arbitrary values of

**u**. Indeed, according to Eq. (13), we obtain

## 5. Design examples

#### 5.1. Transformation of a circular beam to a rectangular beam

As a first example, we consider the design of an optical element transforming a beam with circular cross-section and uniform irradiance ${I}_{0}(\mathbf{u})=1/(\pi {R}^{2})$, $\mathbf{u}\in G=\{({u}_{1},{u}_{2})|{u}_{1}^{2}+{u}_{2}^{2}\le {R}^{2}\}$ to a beam with rectangular cross-section and uniform irradiance *I*(**x**) = 1/(*w*_{1}*w*_{2}), **x** ∈ *D* = {(*x*_{1}, *x*_{2})| |*x*_{1}| ≤ *w*_{1}/2, |*x*_{2}| ≤ *w*_{2}/2}. The wavefronts of the incident and output beams are plane and perpendicular to the *z* axis. The element was designed for the following parameters: the refractive index of the element material *n* = 1.5, the incident beam radius *R* = 1 mm, the sizes of the output rectangular beam *w*_{1} = 5 mm and *w*_{2} = 2.5 mm, the output plane *z* = *f*_{0} = 10 mm, and the parameter *δ* = 2.5 mm. For the considered parameters, the optical path length from the plane *z* = 0 to the plane *z* = *f*_{0} is *L* = *f*_{0} + *δ* = 12.5 mm. In this example, the irradiances *I*_{0}(**u**) and *I*(**x**) are centrally symmetric. Therefore, the central ray originating from the point **u** = (0, 0) will propagate through the optical element without refraction. This allows us to calculate the element thickness at **u** = (0, 0) as *h*_{0} = *δ*/(*n* − 1) = 5 mm.

The mapping **x**(**u**) = *T*(**u**) was found by solving the LAP of Eq. (8). For the calculation of the mapping, the rectangle *D* was represented as a set of *N* = 143^{2} = 20449 rectangular cells *d _{i}* with the size 0.0282 × 0.0141 mm

^{2}. The circular domain

*G*was also approximated by 20449 square cells

*g*with the size 0.0124 × 0.0124 mm

_{i}^{2}. For the solution of the LAP (8), the auction algorithm [22] implemented in [32] was used. For the considered example, the time of solving the LAP was less than 3 minutes. For comparison, let us note that the calculation of this mapping using finite-difference methods for solving PDE of Monge–Ampère type [7–9] requires the solution of a system of 20449 nonlinear equations and is by several orders of magnitude slower [17].

The optical surfaces of the element were reconstructed from the computed discrete mapping using the technique described above in Section 4. The function *g*(**x**) was represented as a two-dimensional B-spline at *q* = 4 and *N*_{1} = *N*_{2} = 10. The coefficients of the spline were found using the least squares method from the condition of minimizing expression (10). Then, the second optical surface defined as spline in Eq. (9) was recalculated on a uniform 400 × 400 rectangular grid (*N _{g}* = 400

^{2}). Then, using Eqs. (12) and (13), the first optical surface was calculated on a uniform 143 × 143 square grid. The resulting optical surfaces are shown in Fig. 2(a). For the sake of illustration, Figs. 3(a) and 3(b) depict the mapping

*T*(

**u**) calculated using Eq. (14). Figure 3(a) shows a square grid in the domain

*G*, whereas Fig. 3(b) shows its image by the mapping

*T*(

**u**). It should be noted that in the considered example, the angles of ray deflection by the optical surfaces (which reach maximum at the corners of the rectangle

*D*) are as large as 25°. Hence, the optical element is designed in a non-paraxial case.

In order to verify the proposed design method, the performance of the optical element was numerically simulated using the commercial ray-tracing software TracePro [33]. To do this, the optical surfaces shown in Fig. 2(a) were exported to a computer-aided design software Rhinoceros [34], in which a 3D-model of the element was created [Fig. 2(a)]. Then this model was exported to TracePro for simulation. Figures 2(b) and 2(c) show the simulated normalized irradiance distributions generated by the optical element in the planes *z* = 10 mm and *z* = 30 mm. The simulation results demonstrate the formation of a rectangle with the required size and near-uniform irradiance. The normalized root-mean-square deviations (NRMSD) of the numerically simulated irradiance distributions from a constant value amount to 5.3% and 6.1% at *z* = 10 mm and *z* = 30 mm, respectively. The fact that the size of the rectangle remains unchanged at different distances from the element shows that the wavefront of the beam generated by the element is nearly plane. Actually, the root-mean-square deviation of the optical path length from a constant value at *z* = 10 mm does not exceed 1 nm.

#### 5.2. Transformation of a circular beam to a triangular beam

As a second example, let us consider the design of an optical element transforming a uniform circular incident beam to a uniform triangular beam. The equilateral triangle (the domain *D*) has a 3-mm side, with the rest geometric and simulation parameters (*R*, *f*_{0}, δ, *n*, *N*, *q*, *N*_{1}, *N*_{2}, *N _{g}*) remaining the same. The designed optical surfaces and the mapping

*T*(

**u**) are depicted in Fig. 4(a) and Figs. 3(a) and 3(c), respectively. Figures 4(b) and 4(c) show the simulated normalized irradiance distributions generated by the optical element in the planes

*z*= 10 mm and

*z*= 30 mm obtained using TracePro. The simulation results show that a triangle-shaped irradiance domain with the required size is generated. The NRMSD of the calculated irradiance distributions from a uniform value amount to 5.6% and 6.4% at

*z*= 10 mm and

*z*= 30 mm, respectively. As in the previous example, the root-mean-square deviation of the optical path length from a constant value at

*z*= 10 mm does not exceed 1 nm.

#### 5.3. Transformation of a circular beam to a cross-shaped beam

The final and the most challenging example involves the design of an optical element transforming a uniform incident beam with circular cross-section into a beam with cross-shaped uniform irradiance distribution. In this case, the domain *D* corresponds to a cross obtained by a superposition of two mutually perpendicular rectangles with the sizes 2.8 × 1 mm^{2} and 1 × 2.8 mm^{2}. The rest geometric and calculation parameters remain the same. The resulting optical surfaces are shown in Fig. 5(a). Note that the lower surface of the optical element has ‘breaks’ (sharp bends) along the bisector in each quadrant of the coordinate system. Figure 6 depicts a magnified fragment of the lower optical surface with a break, which is highlighted by a red rectangle in Fig. 5(a). These breaks are caused by the discontinuities in the mapping **x**(**u**)= *T*(**u**). Indeed, the rays that strike the lower optical surface on the opposite sides of the bisector are mapped to different parts of the cross. This is confirmed by the form of the mapping *T*(**u**) shown in Fig. 3(d). For instance, the image of a straight blue line shown in Fig. 3(a) has discontinuities near the inner corners of the cross and is composed of several curvilinear segments. It is worth noting that since the inverse mapping *T*^{−1}: *D* → *G* is continuous, the obtained upper optical surface is smooth.

The simulation results obtained using TracePro and shown in Figs. 5(b) and 5(c) demonstrate that a high-quality cross-shaped irradiance distribution is generated. The NRMSD of the calculated irradiance distributions from a constant value are 7.1% and 7.9% at *z* = 10 mm and *z* = 30 mm, respectively. The root-mean-square deviation of the optical path length from a constant value is 1.1 nm at *z* = 10 mm. The authors believe that the somewhat larger NRMSD of the generated irradiance distributions from a constant value is due to the approximation error of the piecewise-smooth lower optical surface in the used CAD software Rhinoceros. This software adopts non-uniform rational basis splines (NURBS) for representing a surface from a point cloud. As a result, the sharp bends of the lower optical surface are somewhat smoothed.

The third example demonstrates one of the advantages of the proposed method, which is the capability of operating with discontinuous mappings. In this case, the known methods based on the numerical solution of an elliptic PDE [7–12] are either not applicable or numerically unstable.

## 6. Conclusion

We proposed a method for designing refractive optical elements for collimated beam shaping. The computation of a ray mapping is reduced to a mass transportation problem with a non-quadratic cost function. In the proposed method, the problem of computing a ray mapping is formulated as a linear assignment problem, which constitutes a discrete version of the corresponding mass transportation problem. A technique for reconstructing an optical surface from the computed ray mapping has also been proposed. As distinct from the methods based on the numerical solution of an elliptic PDE [7–12], the method proposed here is suitable for the design of optical elements with continuous piecewise-smooth optical surfaces.

We have designed refractive optical elements transforming a beam with circular cross-section into variously shaped (rectangular, triangular, and cross-shaped) beams with plane wavefronts. The proposed method is computationally efficient. It requires less than 3 minutes on a standard computer to compute a ray mapping on a 143 × 143 grid. High performance of the designed optical elements is confirmed by the numerical simulation results.

## Funding

Russian Foundation for Basic Research (18-29-03067, 17-47-630164, 18-07-00982); Ministry of Science and Higher Education of the Russian Federation (007-GZ/Ch3363/26); US Air Force Office of Scientific Research (FA9550-18-1-0189).

## Acknowledgments

L. L. Doskolovich, D. A. Bykov, E. S. Andreev, and E. A. Bezus acknowledge support from the Russian Foundation for Basic Research (the development of the design method and the investigation of the designed elements), and from Ministry of Science and Higher Education of the Russian Federation (the implementation of the software for the optical element design). V. Oliker acknowledges support from the US Air Force Office of Scientific Research (application of the SQM for the reconstruction of the optical surfaces).

## References

**1. **B. R. Frieden, “Lossless conversion of a plane laser wave to a plane wave of uniform irradiance,” Appl. Opt. **4**(11), 1400–1403 (1965). [CrossRef]

**2. **J. L. Kreuzer, “Coherent light optical system yielding an output beam of desired intensity distribution at a desired equiphase surface,” U. S. Patent No. 3,476,463 (4 Nov. , 1969).

**3. **P. W. Rhodes and D. L. Shealy, “Refractive optical systems for irradiance redistribution of collimated radiation: their design and analysis,” Appl. Opt. **19**(20), 3545–3553 (1980). [CrossRef] [PubMed]

**4. **J. A. Hoffnagle and C. M. Jefferson, “Design and performance of a refractive optical system that converts a Gaussian to a flattop beam,” Appl. Opt. **39**(30), 5488–5499 (2000). [CrossRef]

**5. **H. Ma, Z. Liu, P. Jiang, X. Xu, and S. Du, “Improvement of Galilean refractive beam shaping system for accurately generating near diffraction-limited flattop beam with arbitrary beam size,” Opt. Express **19**(14), 13105–13117 (2011). [CrossRef] [PubMed]

**6. **X. Hui, J. Liu, Y. Wan, and H. Lin, “Realization of uniform and collimated light distribution in a single freeform-Fresnel double surface LED lens,” Appl. Opt. **56**(15), 4561–4565 (2017). [CrossRef] [PubMed]

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

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

**9. **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] [PubMed]

**10. **Z. Feng, L. Huang, M. Gong, and G. Jin, “Beam shaping system design using double freeform optical surfaces,” Opt. Express **21**(12), 14728–14735 (2013). [CrossRef] [PubMed]

**11. **Z. Feng, L. Huang, G. Jin, and M. Gong, “Designing double freeform optical surfaces for controlling both irradiance and wavefront,” Opt. Express **21**(23), 28693–28701 (2013). [CrossRef]

**12. **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] [PubMed]

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

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

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

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

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

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

**19. **Y. Ding, X. Liu, Z.-R. Zheng, and P.-F. Gu, “Freeform LED lens for uniform illumination,” Opt. Express **16**(17), 12958–12966 (2008). [CrossRef] [PubMed]

**20. **L. L. Doskolovich, N. L. Kazansky, S. I. Kharitonov, and V. A. Soifer, “A method of designing diffractive optical elements focusing into plane areas,” J. Mod. Opt. **43**(7), 1423–1433 (1996). [CrossRef]

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

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

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

**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. **J. Rubinstein and G. Wolansky, “Reconstruction of optical surfaces from ray data,” Opt. Rev. **8**(4), 281–283 (2001). [CrossRef]

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

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

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

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

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

**31. **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] [PubMed]

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

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

**34. **Computer-aided design software Rhinoceros. http://www.rhino3d.com