Designing double freeform surfaces for collimated beam shaping with optimal mass transportation and linear assignment problems

We propose a method for designing refractive optical elements for collimated beam shaping. 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.


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][2][3][4][5]. In the general two-dimensional (2D) case, the design of an optical element is essentially more challenging.
In Refs. [6][7][8], the problem of collimated beam shaping was reduced to finding a solution to an elliptic nonlinear partial differential equation (PDE) [6][7][8]. 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 [9][10][11]. 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 [10,11]. 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 [8,12].
In the basic theoretical works [13,14], the problem of collimated beam shaping was formulated as a Monge-Kantorovich mass transportation problem, and the cost functions for the design of a refractive optical element [13] and of a two-reflector system [14] were derived. It is important to note that the cost function for the refractive optical element is non-quadratic [13], 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) [15]. Such an approach is utilized in a recent work by the present authors [12] 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 [12]. 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 [16]. 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 [6][7][8][9][10][11]. 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.

Formulation of the problem
Consider an incident beam with a plane wavefront perpendicular to the z axis and with the irradiance distribution E 0 (u), u ∈ G, where u = (u 1 , u 2 ) are the Cartesian coordinates in the plane z = 0 (Fig. 1a). 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 E(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 E(x) should satisfy the following normalization condition: where the integrals with respect to du (dx) denote double integrals with respect to du 1 du 2 (dx 1 dx 2 ). Besides, since the output beam has a plane wavefront, the irradiance E(x) will remain unchanged in any plane after the optical element (Fig. 1a). Moreover, since both the input and output wavefronts are plane, an 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 through the z u 1 optical element thickness h 0 at a certain point in the form L = (n − 1)h 0 + f 0 . The quantity h 0 can be treated as a parameter of the problem. 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. 1a). According to the energy conservation law along the ray tubes, the mapping T : G → D conserves the light flux. The energy conservation law can be written in the integral form [12,13]: where ω 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 [13] where ) is the inverse ray mapping. According to Eq. (3), 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 [13] 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: where (5) is defined on energy-conserving mappings 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 E 0 (u), u ∈ G and the sought-for irradiance distribution E(x), x ∈ D are used. Let us note that the cost function defined by Eq. (6) differs from the one utilized in [13] 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 h 0 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 the increase in the parameters h 0 and n the achievable value of the shift also increases.
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 |x(u) − u| 2 ≪ h 0 ) C(s) ≈ −γ + |s | 2 2γ , 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 [8][9][10][11].
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 [6][7][8]. 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.

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) [15]. 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 ∫ 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 i of the domain D the cells g 1 , . . . , g N are mapped. In terms of the LAP, the mapping of the cell g i to the cell d j can be interpreted as the assignment of the j-th task to the i-th agent. According to Eq. (6), the cost of the tasks is described by the matrix where u i , x j are the centers of the cells g i and d j . Thus, the mapping x(u i ) =T(u i ) can be found by solving the following assignment problem: where the minimum is sought over all permutations ( j 1 , . . . , j N ). For the solution of the LAP, efficient polynomial algorithms have been proposed, including Hungarian algorithm [15], Jonker-Volgenant algorithm [17], and auction algorithm [18].

Reconstruction of the optical surface
Due to the discrete character of the LAP in Eq. ] 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.

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 [19]. Assume that the domain D, where g(x) is defined, is bounded by a rectangleD 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 the order q, denoting them by B m (x 1 ), m = 1, . . . , N 1 + q − 2 [19]. 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 where p m,n are the expansion coefficients. Note that if the splines 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 x i =T(u i ), i = 1, . . . , N constitute a least-squares approximation of the derivatives ∂ĝ ∂x k (x i ). The latter are derived from Eq. (3) on the basis of a discrete mappinĝ T, i.e. by assuming u(x i ) =T −1 (x i ) in Eq. (3). Thus, the coefficients p m,n can be obtained by solving the following minimization problem: where p is the vector composed of the spline coefficients p m,n . 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 derivatives of an optical surface, allows one to compensate for the errors resulting from computing the mapping by use of the 'discrete' LAPbased approach [20,21]. Note that for a nonrectangular domain D ⊂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 λ m,n p 2 m,n , where λ is a small number, is added to the minimized function S(p).

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) [22][23][24][25]. 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 g . Using the spline representation of the optical surface g(x) [Eq. (10)], the array of points X i can be calculated on a uniform grid with an arbitrarily small step. In the vicinity of each point u ∈ G, the function f (u) corresponds to a lens segment that focuses the incident plane beam into one of the points X i on the upper optical surface.
To derive the lens equation, let us use Fermat's principle (the condition of a constant optical path length from the points on the incident wavefront to the focal point): where the function Φ i (u) describes the lens surface and is the optical path length of the rays propagating from the plane z = 0 to the focal point X i = (x i , g(x i )). Equation (11) defines the lens surface 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 − 1)h 0 /n, and one of the foci located at the point X i [26].
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 [22][23][24][25]: 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 g . Note that Eq. (13) is also convenient for calculating the mapping T (u) at arbitrary values of u. Indeed, according to Eq. (13), we obtain

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 E 0 (u) = 1/(πR 2 ), u ∈ G = (u 1 , u 2 ) | u 2 1 + u 2 2 ≤ R 2 to a beam with rectangular cross-section and uniform irradiance E( . 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 incident beam radius R = 1 mm, the sizes of the output rectangular beam w 1 = 5 mm and w 2 = 2.5 mm, the optical element thickness h 0 = 5 mm at u = (0, 0), and the refractive index of the element material n = 1.5. For these parameters, the optical path length from the plane z = 0 to the plane z = f 0 = 10 mm is L = (n − 1)h 0 + f 0 = 12.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 i with the size 0.0124 × 0.0124 mm 2 . For the solution of the LAP (8), the auction algorithm [18] implemented in [27] 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 [6][7][8] requires the solution of a system of 20449 nonlinear equations and is by several orders of magnitude slower [16].
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 twodimensional 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. 2a. For the sake of illustration, Figs. 3a and 3b depict the mapping T (u) calculated using Eq. (14). Figure 3a shows a square grid in the domain G, whereas Fig. 3b 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 [28]. To do this, the optical surfaces shown in Fig. 2a were exported to a computer-aided design software Rhinoceros [29], in which a 3D-model of the element was created (Fig. 2a). Then this model was exported to TracePro for simulation. Figures 2b and 2c 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.

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 , h 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. 4a and Figs. 3a and 3c, respectively. Figures 4b and 4c show the simulated normalized irradiance distributions generated by the optical element in the plane 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.

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. 5a. 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. 5a. 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. 3d. For instance, the image of a straight blue line shown in Fig. 3a 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. 5b and 5c demonstrate that  Fig. 5a). The break is marked with a red line.
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 piecewisesmooth 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 [6][7][8][9][10][11] are either not applicable or numerically unstable.

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 [6][7][8][9][10][11], 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.