A Monge–Ampère Problem with Non-quadratic Cost Function to Compute Freeform Lens Surfaces

In this article, we present a least-squares method to compute freeform surfaces of a lens with parallel incoming and outgoing light rays, which is a transport problem corresponding to a non-quadratic cost function. The lens can transfer a given emittance of the source into a desired illuminance at the target. The freeform lens design problem can be formulated as a Monge–Ampère type differential equation with transport boundary condition, expressing conservation of energy combined with the law of refraction. Our least-squares algorithm is capable to handle a non-quadratic cost function, and provides two solutions corresponding to either convex or concave lens surfaces.


Introduction
The optical design problem involving freeform surfaces is a challenging problem, even for a single mirror/lens surface which transfers a given intensity/emittance distribution of the source into a desired intensity/illuminance distribution at the target [1][2][3]. More specifically, the freeform design problem is an inverse problem: "Find an optical system containing freeform refractive/reflective surfaces that provides the desired target light distribution for a given source distribution". Inverse optical design has a wide range of applications from LED based optical products for street lighting and car headlights to applications in medical science, image processing and lithography [1,4].
To convert a given emittance profile with parallel light rays into a desired illuminance profile with parallel light rays, one requires at least two freeform lens/mirror surfaces [2,5]. This freeform problem can be formulated as a second order partial differential equation of Monge-Ampère (MA) type, with transport boundary conditions, applying the laws of geometrical optics and energy conservation [2,3,6,7]. For the two-reflector problem [2,8], one can obtain the following mathematical formulation using properties of geometrical optics, i.e., where a 1 , a 2 are constants and u 1 (x), u 2 (y) represent the location of the optical surfaces, and |.| denotes the 2-norm for vectors. The right hand side function c(x, y) in the above expression is a quadratic (cost) function. Assuming convex/concave reflective surfaces the ray-trace map can be uniquely expressed as the gradient of u 1 , i.e., Furthermore, using conservation of energy, one can derive a second order partial differential equation (PDE) of MA-type. We refer to this equation as the standard MA-equation, representing optical systems characterized by a quadratic cost function.
In this article, we show that a similar mathematical expression can be obtained for the freeform surfaces of a lens with parallel ingoing and outgoing light rays applying the laws of geometrical optics: where b 1 , b 2 , b 3 are constants, and u 1 (x), u 2 (y) represent the first and second refractive surface of the lens, respectively. For the freeform lens the cost function c(x, y) is no longer a quadratic function, and the ray-trace map can not be expressed as the gradient of some function, we provide more details in Sect. 2. Energy conservation results in a complicated MA-type equation. In this article, we will present a numerical algorithm to compute the freeform surfaces of a lens characterized by a non-quadratic cost function. A rigorous analysis of the existence and uniqueness of weak solutions of similar lens design problems is presented in Oliker's work [9,10]. There are several numerical methods which can be employed to compute freeform surfaces of optical systems characterized by a quadratic cost function. However, to the best of our knowledge, this paper is the first to describe a numerical method for the MA-equation with non-quadratic cost function. Froese et al. [11][12][13] solve the standard MA-equation within the framework of optimal mass transport (OMT). Applying the theory of viscosity solutions, they refine the solution using an iterative Fourier-transform algorithm with overcompensation. In recent publications [5,14,15], the authors obtain freeform optical surfaces by solving the standard MA-equation using Newton iteration. These numerical methods require an initial guess which is obtained through the OMT problem. Brix et al. [3,16] solve the standard inverse design problem using a collocation method with a tensor-product B-spline basis. Glimm and Oliker [8,17] show that the illuminance control problem can be solved using an optimization approach instead of solving a MA-type differential equation. Further, a similar approach to design freeform surfaces of a lens is developed by Rubinstein and Wolansky [18].
A least-squares (LS) method [2,7] has been presented to solve the standard MA-equation to compute single reflector/lens or double reflector freeform surfaces optical systems. The method provides the optical mapping which transfers the given emittance of the source into the desired illuminance at the target, and the freeform surfaces are obtained via this mapping.
However, the coupled freeform lens surfaces design problem corresponds to a nonquadratic cost function. The goal of this paper is to present a numerical method which is applicable to design an optical system corresponding to a non-quadratic cost function. Here, we present a fast and effective extended least-squares (ELS) method to construct the freeform surfaces of the lens. The ELS-method is a two-stage procedure like the LS-method: first we determine an optimal mapping by minimizing three functionals iteratively, next, we compute the freeform surfaces from the converged mapping. In the first stage, there are two nonlinear minimization steps, which can be performed point-wise, like in the LS-method. In the third step two elliptic partial differential equations have to be solved. For the LS-method, these are decoupled Poisson equations. However, in the ELS-method these are coupled elliptic equations.
Our least-squares method is quite generally applicable since it can handle arbitrary twice differentiable cost functions c(x, y), also in other fields of science and engineering such as optimal transport theory, shape optimization, compression modeling, relativistic theory, incompressible fluid flow, economics, astrophysics, atmospheric sciences etc. For the interested reader, we refer to the following: Evans' survey notes [19], articles of Bouchittè-Buttazzo [20,21], Gangbo's lecture notes [22], and paper of Benamou-Brenier [23]. However, we restrict ourselves to the computation of freeform optical systems. This paper is structured as follows. In Sect. 2 we explain the geometrical structure of the optical system and formulate the mathematical model. The detailed procedure of the proposed least-squares method is shown in Sect. 3. We apply the numerical method to four test problems in Sect. 4 and verify the solutions using a ray tracing algorithm [2]. Finally, a brief discussion and concluding remarks are given in Sect. 5.

Formulation of the Problem
The geometrical structure of a lens optical system is shown schematically in Fig. 1. Let (x 1 , x 2 , z) ∈ R 3 denote the Cartesian coordinates with z the horizontal coordinate and x = (x 1 , x 2 ) ∈ R 2 the coordinates in the plane z = 0, denoted by α 1 , and let S be a bounded source domain in the plane α 1 . The source S emits parallel light rays which propagate in the positive z-direction. The emittance, i.e., luminous flux per unit area (for an introduction to photometry quantities see e.g. [24, p. 7-9]), of the source is given by where f is a non-negative integrable function on the domain S. The target at a distance > 0 from the plane α 1 is denoted by T .
The incoming light rays are refracted at the first lens surface L 1 , propagate through the lens and are refracted again at the second lens surface L 2 , to create a parallel bundle of light rays in the positive z-direction. The index of refraction of the lens n > 1 and the surrounding medium is air with refractive index of unity. The lens surfaces are defined as z ≡ u 1 (x), x ∈ S and w ≡ − z = u 2 (y), y ∈ T , respectively, where y = (y 1 , y 2 ) ∈ R 2 are the Cartesian coordinates of the target plane α 2 .
The goal is to design a lens system such that after two refractions the refracted rays must form a parallel beam, propagating in the positive z-direction, and provide a prescribed illuminance g(y) [lm/m 2 ] at the plane α 2 : z = , where g > 0 is a positive integrable function on the domain T . It is assumed that both L 1 and L 2 are perfect lens surfaces and no energy is lost in the refraction.

Geometrical Formulation of the Freeform Lens
In this section, we first give an expression for the ray-trace map, and secondly we derive a mathematical formulation for the location of the freeform surfaces using the laws of geometrical optics.
The mapping m can be derived by tracing a typical ray through the optical system. Let us consider a ray emitted from a position x ∈ S on the source and propagating in the positive z-direction, letŝ be the unit direction of the incident ray. The ray strikes the first lens surface L 1 , refracts off in directiont, strikes the second lens surface L 2 , and reflects off, again in the directionŝ. The unit surface normal of the first lens surface L 1 , directed towards the light source, is given byn Throughout this article, we use the convention that a hat denotes a unit vector. According to Snell's law [24,25], the directiont =t(x) of the refracted ray can be expressed aŝ where η = 1/n < 1 with n the refractive index of the lens and If we writet = (t 1 , t 2 , t 3 ) T then the first two components of the vectort, can be written as a function of the third component of the vectort as The image on the target of the point x ∈ S is the point y ∈ T under the ray trace mapping m, i.e., y = m(x), x ∈ S. This mapping can be obtained by the projection oft on the plane α 1 , i.e., where d(x) is the distance between surfaces L 1 and L 2 along the ray refracted in the directiont(x). The distance d(x) between the lens surfaces can be obtained using properties of geometrical optics: The total optical path length L(x) corresponding to the ray associated with a point x ∈ S, is given by The theorem of Malus and Dupin (the principle of equal optical path lengths) states that the total optical path length between any two orthogonal wavefronts is the same for all rays [26, p. 130]. As we deal with two parallel beams of light rays, the wavefront coincides with planes α 1 and α 2 . Therefore, the total optical path length will be independent of the position vector x, i.e., L(x) = L. The horizontal distance between the source and the target plane is given by Subtracting Eqs. (9) and (10), and using Eq. (5), we obtain the following expression where where β = L − is the "reduced" optical path length. Substituting (7) and (11) in (8), Now, substituting t 3 in the above equation from the law of refraction (5), the mapping m is given by the relation Next, we derive a mathematical expressions for the location of the lens surfaces. An alternative expression for the distance d reads Thus, from Eqs. (9) and (14), we obtain which can be rewritten as and after elementary algebraic derivations, we obtain This is a mathematical expression for the location of the lens surfaces but the sign in front of the square root is unknown yet. To determine this we proceed as follows. Using Eqs. (9) (with L(x) = L) and (14), we can show that Substituting, this expression in Eq. (15), we obtain First, we check the sign of the expression nβ − d(n 2 − 1). Substituting d from Eq. (11), the expression becomes Since β > 0 and n − t 3 > 0, it remains to check the sign of 1 − nt 3 . Using the vectorial form of the law of refraction (5) and expression (6), we can write as n > 1. Thus we have to choose the negative sign in front of the absolute value in Eq. (16). Hence, we obtain Substituting d from relation (9), the above expression becomes In the above equation, the right hand side equals the left hand side for the minus sign, therefore we have to choose minus sign in (15). Thus the mathematical expression for the lens surfaces becomes These kind of freeform optical design problems are closely related to the mass transport problem [10,27]. The right hand side function c(x, y) is known as the cost function in OMT theory.
To conclude, we have derived a mathematical formulation representing the freeform lens optical system which is given in (19). Also, we obtained the expression (13) for the ray-trace mapping m. Next, we formulate a second order partial differential equation for the freeform lens.

Energy Conservation for the Freeform Lens
Recall, that f ≥ 0 and g > 0 are integrable functions and no energy is lost in the light transfer process. Thus energy conservation is given by The key tool for the design of such an optical system is to find a mapping y = m(x) : S → T that satisfies the energy conservation constraint (20) for each measurable set A ⊂ S, i.e., and after a change of variables the constraint becomes where Dm is the Jacobian of the mapping m, which measures the expansion/contraction of a tube of rays due to the two refractions. The accompanying boundary condition is derived from the condition that all the light from the source domain S must be transferred into the target domain T , i.e., stating that the boundary of the source S is mapped to the boundary of the target T . This is a consequence of the edge ray principle [28]. Next, we derive a MA-type equation for the freeform lens using the energy conservation constraint (22) and the mathematical formulation (19) for the location of the lens surfaces. We assume that both lens surfaces u 1 and u 2 are either c-convex or c-concave functions. According to the following definition, the lens surfaces u 1 and u 2 are c-convex if alternatively, these are c-concave if For a continuously differentiable function c ∈ C 1 (S × T ), the c-convex/concave functions u 1 and u 2 are Lipschitz continuous [27,29], and the mapping y = m(x) is implicitly given by the relation which is a necessary condition for (24b) and (25b), and holds under the condition that the Jacobi matrix C = D xy c defined by is invertible. For our optical problem the mapping m given by relation (13) satisfies relation (26) indeed. The matrix C is symmetric negative semi-definite which is a consequence of the fact that the function c depends on |x − y|. This can be verified as follows: let us rewrite the cost function (19) as By differentiating (28) with respect to x and y, we obtain Differentiating one more time with respect to x, we conclude that Evaluating all derivatives, we obtain the following expression We can rewrite the above expression as follows: where γ = n 2 /(n 2 − 1) > 0. Sincec < 0, we conclude that det(C) > 0 and tr(C) ≤ 0 hence the matrix C is symmetric negative semi-definite.
Since the function c(x, y) defined in (19) is continuously differentiable, from relation (26), we deduce where D 2 u 1 is the Hessian of u 1 . The matrix P = D 2 u 1 (x) − D xx c is negative semi-definite for a c-concave pair (u 1 , u 2 ) and positive semi-definite for a c-convex pair (u 1 , u 2 ). In the following, we discuss the convex case, thus we require the matrix P to be positive semidefinite. Substituting Dm from (33) into the energy conservation condition (22), we obtain We know that the 2 × 2 matrix P is positive semi-definite if and only if tr( P) ≥ 0 and det( P) ≥ 0.
Because det(C) > 0 and the right hand side functions f ≥ 0, g > 0 in Eq. (34), it is obvious that det( P) ≥ 0. So, the only requirement left is tr( P) ≥ 0 for convex optical surfaces. In the following section, we give a detailed description of the ELS-algorithm to solve the MA-equation (34) with the boundary condition (23) and constraints (35). The method presented here is based on [7]. Compared to [7] we deal with a non-quadratic cost function that results in the presence of the matrix C in (34).
optical system [2], which is also a quadratic cost problem. Our version of the least-squares method was inspired by publications by Caboussat et al. [30,31], who developed a leastsquares method for the Monge-Ampère-Dirichlet problem. An extension of their method to the three-dimensional equation is presented in [32].
In this section, we extend the least-squares method to compute the freeform surfaces of a lens characterized by a non-quadratic cost function. The ELS-method is a two-stage procedure. In the first stage we calculate the optimal mapping by minimizing three functionals iteratively, and in the second stage we compute the freeform surfaces from the mapping in the least squares sense.

First Stage: Calculation of the Mapping
First, we calculate the mapping m using the least-squares method for the lens optical system as follows: we enforce the equality CDm = P by minimizing the following functional The norm used in this functional is the Frobenius norm, which is defined as follows. Let A : B denote the Frobenius inner product of the matrices A = (a i j ) and the Frobenius norm is then defined as || A|| = √ A : A. Next, we address the boundary by minimizing the functional We combine the functionals J I for the interior and J B for the boundary domain by a weighted average: The parameter α (0 < α < 1) controls the weight of the first functional compared to the second functional. The variables b, P, and m are elements of the following spaces respectively. The minimizer gives us the mapping m which is implicitly related to the surface function u 1 . We calculate this minimizer by repeatedly minimizing over the three spaces separately. We start with an initial guess m 0 , which will be specified shortly, and we calculate the matrix C(x, m 0 ) at the initial guess m 0 . Subsequently, we perform the iteration Next, we compute the matrix C(x, m n ) before going to the next iteration. We initialize our minimization procedure by constructing an initial guess m 0 which maps a bounding box of the source area S to a bounding box of the target area T . Without loss of generality we assume the smallest bounding box of the source and the target are rectangular and denote these by [a min , a max ] × [b min , b max ] and [c min , c max ] × [d min , d max ], respectively. Then the initial guess reads: Note that the corresponding Jacobi matrix Dm 0 of the initial condition is symmetric (in fact diagonal) negative definite. The matrix C is also negative definite, moreover from relation (32) we conclude that c 11 , c 22 < 0 which implies that the matrix P = C(x, m 0 )Dm 0 is positive definite. Thus this initialization satisfies our requirement tr( P) ≥ 0.
Obviously, the minimization steps in (41) as well as the computation of the optical surfaces is done numerically. To that purpose we discretize the source S with a standard rectangular N 1 × N 2 grid for some N 1 , N 2 ∈ N, so the grid points x i j = (x 1,i , x 2, j ) are defined as We start the iteration process (41) using initial guess m 0 . Here each iteration consists of four steps: we perform the three minimization steps (41a)-(41c), and fourthly we update the matrix C at every iteration. In this article, we give a detailed description of the minimization steps (41b) and (41c). The minimization step first (41a) is simple and direct, and performed point-wise because no derivative of b with respect to x appears in the functional, more details can be found in [7]. Finally, from the converged mapping m, we compute the first lens surface u 1 via relation (26) in a least-squares sense, and the second lens surface u 2 from relation (19), see Sect. 3.2.
Minimizing procedure for P We assume m fixed and minimize J I (m, P) over the matrices under the condition (34). Since the integrand of J I (m, P) does not contain derivatives of P, the minimization procedure can be done pointwise. So, we need to minimize ||C D − P|| for each grid point x i j ∈ S, where D is the central difference approximation of Dm. Let's define with where δ x 1 and δ x 2 are the central difference approximations of ∂/∂ x 1 and ∂/∂ x 2 , respectively. Note that, the matrices C, D, P and Q all depend on x i j . For the sake of brevity we omit x i j . Moreover, we want to avoid crossing grid lines, i.e., intersection of images of grid lines in S, and for that reason we require d 11 , d 22 > 0. This can be achieved by imposing for a threshold value ε > 0. This implies that m 1 (x i+1, j ) < m 1 (x i, j ) and m 2 (x i, j+1 ) < m 2 (x i, j ) for all x i, j ∈ S, which assures that there is no crossing of grid lines. In our computations we choose ε = 10 −8 . Note that the matrix P is symmetric but the matrix D need not be symmetric and d 12 , d 21 < 0 is possible. Next we define the function Also, we define the matrix Q S as the symmetric part of the matrix Q, i.e., with q S = 1 2 (q 12 + q 21 ). The function H S corresponding to the symmetric matrix Q S is defined as Since (q 12 − q 21 ) 2 is independent of p 11 , p 22  This problem can be solved analytically, and we will show that for given q 11 , q 22 , q S and f /g there exist at least one and at most four real solutions, see "Appendix A". From these we have to select the ones that give rise to a negative semi-definite matrix P, and we will also show that this is always possible. Finally, we compare the values of H S ( p 11 , p 22 , p 12 ) to find the global minimum. The possible minimizers of (49) are obtained introducing the Lagrangian function Λ, defined as where μ is the Lagrange multiplier. By setting all partial derivatives of Λ to 0 we find the critical points of Λ and this gives the following algebraic system where λ = μ/ det(C). The system (51a)-(51c) is linear in p 11 , p 22 and p 12 , and is regular if λ = ±1 (we discuss the singular cases in the "Appendix A"). In the case when λ = ±1, we calculate the critical points by inverting the system, i.e., we express p 11 , p 22 and p 12 in terms of λ as Substituting these expressions in Eq. (51d) gives the following quartic equation with coefficients given by Furthermore, from Eqs. (51a)-(51b) the condition (49c) becomes and consequently, we need to select Lagrange multipliers that satisfy the above condition. It can be shown that the quartic Eq. (53) has at least two real roots, and one of them is less than − 1 and other one is greater than − 1 (see "Appendix A"). The convexity condition (54) can be satisfied by choosing the appropriate values of λ and tr( Q S ), and the minimizers of H S are given by (52).

Minimizing procedure for m
In this section, we describe the minimization step (41c). The minimizing procedure for m differs from the procedure given in [7] because we have an extra matrix C in the function J I which results in two coupled elliptic equations for the components of the mapping m instead of decoupled Poisson equations. We assume P and b are fixed, and minimize J (m, P, b) over the functions m ∈ M using calculus of variations, i.e., P and b are given in all grid points x i j ∈ S. We want to compute m on the grid covering S. Here, we drop the indices n and n + 1, for ease of notation. In the calculations that follow, we use the identity for the Frobenius norm of matrices, i.e., The first variation of the functional J with respect to m in the direction η ∈ [C 2 (S)] 2 is given by δ J (m, P, b)[η] = lim The minimizer is obtained by setting the variation equal to 0, i.e., Let us define the column vectors p 1 , p 2 , c 1 and c 2 as We can split the first integrand of the final expression in (56) as follows where the vectors v 1 and v 2 are column vectors of the matrix and by defining W = V T = [w 1 , w 2 ], we can rewrite the first integral of the final expression in (56) as Letn denote the unit outward normal at the boundary ∂S. Using the vector-scalar product rule [33, p. 576] and the identity derived from the Gauss's theorem, the integrals in the rhs of (56) become Substituting the integral in the final expression of (56), the minimizer can be obtained from the following relation where m k and b k (k = 1, 2) are the components of the vectors m and b, respectively. Choosing η 2 = 0 and applying the fundamental lemma of calculus of variations [34, p. 15] for η 1 , we find that m 1 and m 2 satisfies, almost everywhere, the equation Similarly, choosing η 1 = 0 and applying the fundamental lemma of calculus of variations for η 2 , we obtain We can rewrite these equations as follows in matrix-vector form These are two coupled elliptic equations with Robin boundary conditions for the two components m 1 and m 2 of the mapping m [35, p. 160]. The above equations are in divergence form which motivates us to discretize the equations using the finite volume method [35, p. 84-88], for more details see "Appendix B".

Second Stage: Calculation of the Freeform Surfaces
We compute the lens surfaces assuming that a numerical approximation of m on the grid is available. We compute the first lens surface u 1 (x) from the converged mapping m using relation (26) in the least-squares senses, i.e., (67) We calculate the minimizing function u 1 (x) using calculus of variations. The first variation of the functional I in (67) in a direction v is given by The minimizer is given by Using the Gauss's identity (61), we conclude from (69) that Applying the fundamental lemma of calculus of variations [34, p. 15], we find This is a Neumann problem, and only has a solution if the compatibility condition is satisfied, which reads By Gauss's theorem, this is satisfied automatically. The solution of the Poisson equation with Neumann boundary condition is unique up to an additive constant. To make the solution unique, we have added the constraint u 1 (x 1 , x 2 ) = 1 at the first discretized left most corner point. We solve this problem using standard finite differences, and the discretized system is solved in Matlab using LU decomposition. The second lens surface is calculated from the relation (19), by substituting the converged mapping m(x) and the first lens surface u 1 (x), thus we have The numerical algorithm is summarized as follows. We start the minimization procedure using the initial guess m 0 given by (42) for the discretized source domain S. Subsequently, we iteratively perform the steps given by (41a), (41b) and (41c). The first and second steps are minimization procedure for b and P, respectively, and both are performed pointwise. The third step is a minimization procedure for the mapping m and is performed by solving two coupled elliptic boundary value problems given by (64) and (65). Next, we update the matrix C given by (27). Finally, after convergence of the iteration (41), the first lens surface is computed from the mapping m by solving Poisson problem (71), and the second lens surface is computed from relation (73).

Numerical Results
We apply the algorithm to four test problems to compute c-convex lens surfaces: first, we map a square with uniform emittance into a circle with uniform illuminance, second, we map an ellipse with uniform emittance into another ellipse with uniform illuminance, third, we map a square with uniform emittance into a non-convex (flower) target distribution, and finally, we challenge our algorithm to map the same distribution into a light pattern given by a picture on the target screen. The numerical results are verified by our self-built ray tracer based on the Monte-Carlo method [2].

From a Square to a Circle
In the first test problem, we design an optical system of lens surfaces that transforms the uniform emittance of a square into a circle with uniform illuminance. The source domain is given by S = [−1, 1] × [−1, 1] and the target domain by T = {y ∈ R 2 ||y|| 2 ≤ 1}. The light source emits a parallel beam of light rays with uniform emittance, i.e., The target plane is at a distance = 40 from the source plane, and we have fixed the refractive index n = 1.5 and the reduced optical path length β = 3π for all numerical problems. The target T is illuminated by a parallel beam of light rays with uniform illuminance, i.e., Note that the energy conservation condition (20) is satisfied. We discretize the source domain S uniformly with 200 × 200 grid points. We have a different grid for the boundary, and found from various experiments that the number of boundary grid points N b does not influence the convergence of the algorithm if it is chosen large enough. Since a large value of N b does not significantly increase the calculation time, we have chosen N b = 1000. We also observed from various experiments that α = 0.65 is a good choice for α to have residuals in J I and J B close together. Choosing α too large or too small slows down convergence. We stopped the algorithm after 200 iterations because J I and J B stall. The resulting mapping after 200 iterations is shown in Fig. 2a, and the convergence history of the algorithm is shown in Fig. 2b. The algorithm performed efficiently, the boundary and interior functionals for the circular target have converged well with residuals of approximately 2.35 × 10 −7 .

From an Ellipse to Another Ellipse
In the second test case, we apply the algorithm to map a uniform emittance of an ellipse into another ellipse with uniform illuminance. The source domain is given by S = {(x 1 , x 2 ) ∈ R 2 4x 2 1 + x 2 2 ≤ 4}, see Fig. 3a, and the target domain by T = {(y 1 , y 2 ) ∈ R 2 y 2 1 +4y 2 2 ≤ 4}, i.e., rotate an ellipse distribution to another ellipse over π/2. We have f (x 1 , x 2 ) = g(y 1 , y 2 ) = 1/4π. We use a 200 × 200 grid with 1000 points on the boundary, the reduced optical path length β = 3π, and the weight parameter α = 0.65. The mapping after 200 iterations is shown in Fig. 3b, and the source distribution in grid format shown in Fig. 3a. The algorithm exhibits almost the same convergence as shown in Fig. 2b.

From a Square to a Non-convex (Flower) Target
In the third test case, we test the algorithm for non-convex flower shaped targets. We apply the algorithm to map a uniform emittance of a square into uniformly illuminated non-convex targets. The source domain is given by the square [−1, 1] × [−1, 1] with f (x 1 , x 2 ) = 1/4, and the target domain is defined in polar coordinates as where ρ(θ) is the distance to the origin and θ is the counterclockwise angle with respect to the y 1 -axis in the target plane. We test the algorithm for the four values e ∈ {0.1, 0.15, 0.20, 0.25} which represent the deviation of the target domain from a convex shape. We use 200 × 200 grid with 1000 points on the boundary, the reduced optical path length β = 3π, and the weight parameter α = 0.50. The mappings after 250 iterations are shown in Fig. 4. The residual J after 250 iterations is 2.73 × 10 −7 , 7.05 × 10 −6 , 3.93 × 10 −5 and 9.99 × 10 −5 , respectively. The convergence problem arises for target domains which strongly deviate from a convex shape, but if the shape deviates mildly from convex, the algorithm performs satisfactorily, see Fig. 4a-d.

From a Square to a Picture
The fourth test problem is to design an optical system of lens surfaces which transforms a square uniform bundle of parallel light rays into a target distribution corresponding to a given picture. Here, we challenge our algorithm for a picture showing costumes of the Indian classical dance Bharatanatyam, see Fig. 5. The emittance of the light source is again the same as defined in (74) and the parameters of the optical system are also the same as defined in Sect. 4.1. The desired target illumination g(y 1 , y 2 ) is given by the grayscale test image shown in Fig. 5. Because the target distribution contains many details, e.g., the pattern of costumes and jewellery, it provides a challenging test for our algorithm. Note that the picture is converted into grayscale and contains some black regions, which results in g(y 1 , y 2 ) = 0 for some (y 1 , y 2 ) ∈ T . This would give division by 0 in the least-  squares algorithm. Therefore, we have increased the illuminance to 5% of the maximum value if it is less than this threshold value. We discretized the source S on a 500 × 500 grid, with 1000 boundary points. The convergence history of the algorithm is shown in Fig. 6b for α = 0.70. We stopped the algorithm after 150 iterations, because J I and J B did no longer seem to decrease. The resulting mapping is shown in Fig. 6a, the image details can be recognized in the grid. The optical system is verified using the ray tracing algorithm [2]. We ran our ray tracing algorithm for 10 million uniformly distributed random points on the source to compute the actual illumination pattern produced on the target. The target illuminance for 10 million rays is plotted in the Fig. 5. The output images is very close to the corresponding original image, although the image is slightly blurred, but even complex details can be identified. The functions u 1 (x) and u 2 (y) representing the freeform surfaces L 1 and L 2 of the lens, respectively, are shown in Fig. 7. The lens surfaces are convex on their respective domains and an alternative representation of the mapping can be seen as contour of grids on the second lens surface.

Discussion and Conclusion
We introduced a least-squares method to compute freeform surfaces of an optical system corresponding to a non-quadratic cost function. The method is an extended version of the least-squares method, earlier introduced in [7]. Furthermore, we presented a new generic (in term of cost function) minimization procedure of P for the functional J I . Moreover, we have shown that the minimization procedure of the mapping m for the total functional J consists of coupled elliptic PDEs.
We presented the extended least-squares method to compute coupled freeform surfaces of a lens. Our method can compute freeform surfaces of any optical system corresponding to a twice continuously differential cost function, which demonstrates the wider applicability of the method. The ELS-method also shows good performance for a non-convex target domain: as long as the domain does not deviate too much from a convex shape.
The algorithm is very time and memory efficient, and provides both convex and concave optical surfaces which makes it very suitable to use for these type of problems. Furthermore, we have applied the method to a very challenging problem containing the details of the costumes of the Indian classical dance Bharatanatyam and obtained a high resolution, preserving details of the original picture.
In future work we would like to apply the algorithm to more complex cost functions, e.g., point light sources and far-field problems. Also, we would to explore the applicability of the Monge-Ampère solver in other fields of science and engineering.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A Solution of the Quartic Equation
We obtain four possible solutions of the quartic equation (53) using Ferrari's method [36, p. 32]. The key idea is to rewrite the quartic equation as two quadratic equations, and by solving both we get solutions of the quartic equation. For detailed solution of the quadratic equations we refer to [7].
Solution of (53) when f > 0 f > 0 f > 0. It can be shown that the problem (53) has at least two real roots. For the real symmetric matrix Q S , we can deduce and using the above relation, we conclude that F(−1) = − tr( Q S ) 2 < 0 and F(1) = tr( Q S ) 2 − 4 det( Q S ) ≥ 0, and the coefficient of λ 4 in the quartic equation (53) is positive. Which imply that (53) has at least two real roots, more precisely one of them is less than −1 and other one is greater than −1. The solution of (53) are given by where y is the solution of a cubic equation in the Ferrari's method. The real roots satisfying the convexity condition (54) are substituted in (52) and (49), yielding the possible minimizers of H S ( p 11 , p 22 , p 12 ). Note that we have division by zero in (78) if y = 0. We find that this happens only when a 1 = 0, i.e., q 11 = q 22 = q S = 0. This is a special case which corresponds to the possibility λ = 1, which we discuss later.
Solution of (53) when λ = 1 λ = 1 λ = 1. If λ = 1, i.e., q 11 = q 22 and q S = 0, and in this case, we cannot invert the system (51a)-(51c) for ( p 11 , p 22 , p 12 ). Therefore we determine the minimum of H S ( p 11 , p 22 , p 12 ) as follows. Using q S = 0 and q 11 = q 22 , the minimization (49) simplifies to wheren is the unit outward normal on the boundary ∂A of A. Now, we create a set of non-overlapping control volumes for the computational grid of the domain S and apply the cell-centred finite volume method, i.e., grid points are located in the centre of the control volume. Let us consider the control volume A ≡ Ω C = [x 1,w , x 1,e ] × [x 2,s , x 2,n ] as shown in Fig. 8, where x 1,w is the x 1 -value at centre of the western cell face Γ w , i.e., x 1,w = x 1,i−1/2 and approximated as x 1,w = x 1 (W ) + x 1 (C) /2, etc, and x 1 (C) = x 1,i , x 1 (W ) = x 1,i−1 , etc. The finite volume method is used to transform equation (89) to a system of discrete equations for the centre point C of the control volume Ω C . First, Eq. (89) is applied for the control volume Ω C . This reduces the equation to one involving only first derivatives. These first order derivatives are replaced with central difference approximations, for more details see [35]. Finally, the integral equation (89) can be descretized as follows a E m 1 (E) + a W m 1 (W ) + a N m 1 (N ) + a S m 1 (S) + a C m 1 (C) where a E = |c 1 | r 1 (C) = 1 h 1 (c 1 · p 1 ) e − (c 1 · p 1 ) w + 1 h 2 (c 1 · p 2 ) n − (c 1 · p 2 ) s .