Freeform lens design for a point source and far-field target

The ﬁeld of freeform illumination design has surged since the introduction of new fabrication techniques that allow for the production of non-axially symmetric surfaces. Freeform surfaces aim to efﬁciently control the redistribution of light from a particular source distribution to a target irradiance, but designing such surfaces is a challenging problem in the ﬁeld of nonimaging optics. Optical design strategies have been developed in both academia and industry. In this paper, we consider the design of a single freeform lens that converts the light from an ideal (zero-étendue) point source into a far-ﬁeld target. We present a mathematical approach and numerically solve the corresponding generalized Monge–Ampère equation of the optical system. We derive this equation using optimal transport theory and energy conservation. We use a generalized least-squares algorithm that can handle a non-quadratic cost function in the corresponding optimal transport problem. The algorithm ﬁrst computes the optical map and subsequently constructs the optical surface. We demonstrate that the algorithm can generate a peanut-shaped lens for roadlighting purposes and a highly detailed lens that produces an image on a projection screeninthefarﬁeld.


INTRODUCTION
LED lighting systems have become increasingly popular in the last decade due to their high-energy efficacy and long lifetime. A major advantage is that an LED light operates at lower temperatures than conventional light sources, which allows for the use of optical components composed of plastic materials. The function of these plastics is to transform the light from the LED source into a required light output pattern of the lighting system. Optical diamond turning techniques are capable of producing plastics via injection molding in arbitrary shapes at high precision. These techniques have pushed the field of illumination optics to develop sophisticated and highly precise methods to compute freeform (i.e., non-axially symmetric) shapes that convert the energy of LED light sources to a desired energy (intensity) distribution in the far field.
In this paper, we use a point light source approximation in the freeform lens design for LEDs. We compute the shape of the freeform lenses using inverse methods, rather than forward methods such as ray tracing. This involves finding numerical solutions to a fully nonlinear second-order elliptic partial differential equation known as a generalized Monge-Ampère equation. Solutions to this partial differential equation for a lighting system are the optical surfaces we are interested in.
The generalized Monge-Ampère equation of an optical system can be derived using concepts and techniques in geometrical optics, optimal transport theory, and energy conservation. The calculation of the freeform shape is equivalent to solving an optimal mass transport problem, i.e., all the light from the source needs to arrive at the target such that all light rays traverse the least optical path length. For many optical systems, we can derive a cost function that describes a necessary relation between the direction vectors of incoming and outgoing light rays (p. 55, [1]). This cost function can be derived using Hamilton's characteristic functions, which provide a measure for the optical path length traversed by the rays. Existence and uniqueness of the solution to the optimal transport problem have been studied extensively, under the condition of convexity of the domain and c -convex or c -concave Kantorovich potential [2][3][4].
Currently, there exists a wide variety of methods to solve the standard Monge-Ampère equation for quadratic cost functions [5][6][7][8]. In this paper, we consider a lens that transforms the energy distribution of a point source into a required far-field intensity. This system has a logarithmic cost function in the optimal transport problem. Generalized Monge-Ampère equations with conditions for existence, uniqueness, and smoothness which also have a logarithmic cost function, were derived in [9][10][11][12].
The generalized Monge-Ampère equation for a point source and freeform lens is increasingly treated from an optimal transport point of view. In Gutiérrez et al. [13,14] a logarithmic cost function is derived, and the existence and uniqueness of weak solutions are established, based on a similar approach in [11,15] dedicated to reflective surfaces.
The solution strategies to compute the optical surface can be roughly divided into three categories: (1) direct methods to numerically solve the generalized Monge-Ampère equation [16][17][18][19], (2) optimization methods that consider the associated Monge-Kantorovich mass transportation problem [11,[13][14][15][20][21][22][23], and (3) indirect methods to compute the surface employing ray-mapping techniques [24][25][26][27][28][29]. (1) Wu et al. [16] derive the Monge-Ampère equation and solve the corresponding nonlinear boundary value problem using Newton's method. Brix et al. [17,18] derive the Monge-Ampère equation for a point source with a near-field target and use a collocation method with a tensor-product B-spline basis to calculate reflectors and lenses capable of producing a detailed image on a near-field projection screen. Bosel and Gross [19] present a design algorithm for a single freeform lens for prescribed ideal sources. The generalized Monge-Ampère equation is solved using a root-finding algorithm. The mapping and surface of the lens are initialized using a quadratic cost function and adaptive mesh method. (2) The generalized Monge-Ampère equation in [11,[13][14][15] is reduced to finding a minimizer or a maximizer of a linear functional subject to a linear constraint, which allows for the use of linear programming algorithms. Moreover, Gutierrez [13] shows that the far-field case can also be treated using a nearfield method, taking into account loss of energy via internal reflection and recovering the Fresnel equations. Oliker et al. [20] use the supporting quadric method that involves a pixelation of the target domain and an iterative adjustment of the parameters of tangent quadrics to the optical surface. In Doskolovich et al. [21][22][23], both the far field and near field are considered. In its discrete version, the optimal transport problem is reduced to a linear assignment problem (LAP) for the mapping, based on the construction of an equal-flux grid in the source and target domains. Methods such as the Hungarian algorithm and auction algorithm can be used to solve the LAP. The surface is computed from the mapping using B-splines. Using the LAPbased approach, Doskolovich et al. [21] consider the design of two refractive surfaces for generating irradiance distributions with small angular dimensions as well as off-axis irradiance distributions, designing piecewise-smooth continuous optical surfaces. (3) In ray-mapping techniques, first the standard Monge-Ampère equation is solved (det(D 2 u) = f (x , u, ∇u), with D 2 u the Hessian matrix and f a positive function) and an optical map is constructed as the gradient of the solution m = ∇u. Subsequently, the surface is computed from the mapping using Snell's law and an integrability condition to ensure continuity of the surface. Ray-mapping techniques work well under the paraxial approximation and thin lens approximation, and allow for extensions to multiple freeform surfaces [28]. Feng et al. [30] derive a Monge-Ampère equation of a parameterized outgoing wavefront. The ray mapping is obtained with a Newton-Krylov solver. The freeform surface is calculated in a least-squares sense and is iteratively revised by reconstructing the outgoing wavefront. The method is capable of constructing a double freeform lens with two freeform surfaces (the first surface is pre-defined by an analytical formula) and producing complicated images as target irradiance.
Numerical methods based on optimal transport have also been developed for other non-quadratic cost functions corresponding to different optical systems. Notable examples are [8,31,32], which consider two freeform reflector/lens surfaces for collimated beam shaping.
In this paper, we present a numerical algorithm for a point source with corresponding logarithmic cost function in the optimal transport problem. We use a novel formulation in (Cartesian and polar) stereographic coordinates to represent the source and target domains. In particular, the representation of the source domain in polar stereographic coordinates is convenient for light sources emitting a cone-shaped bundle. The corresponding generalized Monge-Ampère equation in stereographic coordinates can be formulated using the cost function and solved numerically using a generalization of the leastsquares approach presented in [7,33]. The original least-squares approach worked only for a quadratic cost function, extended by Beltman et al. [34] to arbitrary orthogonal coordinate systems. Moreover, the method was extended to non-quadratic cost functions by Yadav et al. [32] and Romijn et al. [35].
In this paper, we further elaborate on this approach and apply it to freeform lens design for a point source and far-field target. The method works by first computing the optical map in an iterative procedure that minimizes the defect in the energy balance. It also imposes a transport boundary condition by minimizing the deviation of the boundary of the optical map to the boundary of the target. Upon convergence of the iterative procedure, the location of the optical surface is calculated from the mapping also in a least-squares sense. We show that the algorithm described in Ref. [32] can also be applied to point light sources using a logarithmic cost function, and we introduce a convenient representation of the source domain in polar stereographic coordinates. First, we present the optimal transport formulation of computing a freeform lens surface, since this is not generally known in the field of applied optics. Subsequently, we will present a broad outline of the numerical method and present the results of two complicated test problems.

MATHEMATICAL FORMULATION
In this section, we derive the generalized Monge-Ampère equation for a freeform lens that redistributes the light from a point source into a far-field target. We first derive the corresponding cost function in optimal transport theory, which we will combine with energy conservation to derive the generalized Monge-Ampère equation. Figure 1(a) schematically illustrates a point source, the lens with refractive index n and far-field target. A beam of light emanates from the point source located at O of the Cartesian coordinate system with (x , y , z) ∈ R 3 . The point source emits beams of light radially outward in the directionŝ =ê r , wherê e r is the radial basis vector in the spherical coordinate system. Unit vectors are denoted by hats. The first surface is spherical, and the freeform lens surface L is described by the parametric equation L : r(φ, θ ) = u(φ, θ )ê r , where u(φ, θ ) > 0 is the radial parameter that describes the location of the lens surface, 0 ≤ φ ≤ π is the zenith, and 0 ≤ θ < 2π is the azimuth in the spherical coordinate system. The intensity of the source is given by f (φ, θ ) [lm/sr], and the required target intensity in the far field is denoted by g (ψ, χ ) [lm/sr], where (ψ, χ ) represents a different set of spherical coordinates, with zenith 0 ≤ ψ ≤ π and azimuth 0 ≤ χ < 2π. The origin of the coordinate system describing the target is the lens surface approximated as point in space. This approximation is called the far-field approximation.
The direction of the incident rayŝ =ê r is unaltered by the first spherical surface of the lens and subsequently refracted by the second freeform surface L in the directiont.
We transform coordinates on the source and target domains from spherical to stereographic. This is convenient, since the vectorsŝ = (s 1 , s 2 , s 3 ) T andt = (t 1 , t 2 , t 3 ) T are defined on the unit sphere S 2 . Hence, |ŝ| = |t| = 1. We define with corresponding inverse projectionŝ Schematic representation of the stereographic projections of the unit sphere S 2 (the points P are projected to P and the points Q to Q ).
We represent the incoming raysŝ by using a stereographic projection from the south pole (0, 0, −1) onto the plane z = 0, as drawn schematically in Fig. 2(a). The stereographic projection in Eq. (1a) is undefined at the south pole, and we consider s 3 = −1 and 0 ≤ φ < π. Likewise, for the outgoing rays, we use a stereographic projection from the south pole (0, 0, −1) with the lens surface as origin in the far-field approximation, as shown in Fig. 2(b). The stereographic projection in Eq. (1b) is undefined at the south pole, and we consider t 3 = −1 and 0 ≤ ψ < π. Transforming to stereographic coordinates, we obtain a bounded source domain that is circular for a coneshaped incoming beam. Likewise, we choose the south pole for the outgoing rays to ensure that the stereographic projection is always defined, assuming the lens surface does not refract rays downwards parallel to the z axis.
We define our source domain X as the supporting domain off (x ) = f (φ(x ), θ (x )), and our target domain Y as the supporting domain ofg ( y) = g (ψ( y), χ ( y)). We refer to m : X → Y as the optical map y = m(x ) from the source set of stereographic coordinates X to the target set of stereographic coordinates Y.
In this paper, we will derive the optical mapping m and corresponding surface u using the cost function in optimal transport theory. We note that we can also derive the mapping m by writing down an expression for the unit surface normaln in terms of the gradient of u and using the vectorial law of refraction t = nŝ − (n(n ·ŝ) + 1 − n(1 − (n ·ŝ) 2 ))n.

A. Cost Function Approach
We derive the cost function using Hamilton's characteristic functions, which are measures of the optical path length between specified source and target planes [36,37]. The characteristics are classified as the point characteristic V (equal to the optical path length between two points), the two mixed characteristics W and W * , and the angle characteristic T. Figure 3 displays a point light source located at the origin O. We consider a lens of refractive index n with a spherical surface closest to O and a freeform surface L.
Let us consider an incident ray propagating in the directionŝ that intercepts the freeform lens surface L at the point P and refracts in the directiont. We assume a target plane at z = L, where the ray intersects at point Q. The freeform lens surface is given by the radial coordinate r = u(ŝ)ê r , where we let u(ŝ) = u(φ, θ ).
The space and direction coordinate vectors on the source plane are given by the two-vectors q s = 0 and p s = (ns 1 , ns 2 ) T , respectively. The space and direction coordinates on the target plane are given by q t and p t = (t 1 , t 2 ) T , respectively. We define where nu(ŝ) is the optical path length from O to P , where we can assume the spherical surface infinitesimally close to O for simplicity, since it does not change the direction of the light rays. d (P , Q) denotes the Euclidean distance between P and Q.
Hamilton's angle characteristic, which depends on the direction of the ray at the source plane and the direction at the target plane, is given by It can be interpreted geometrically as the optical path length between the intersection of the source ray with the plane perpendicular to it going through the origin of the source plane, and the intersection of the target ray with the plane perpendicular to it going through the origin of the target plane. The spatial coordinate q s at the source plane and spatial coordinate q t at the target plane can be derived as (p. 105, [36]) By the former equation, we can thus conclude that the angle characteristic T is independent of the direction coordinate p s , and the expression for T reads We can derive and Hence, we arrive at gives the relatioñ Note: n −ŝ ·t > 0. Transforming Eq. (10) to the stereographic coordinates in Eqs. (1a) and (1b), we arrive at where u 1 (x ) =ũ 1 (ŝ(x )), and u 2 ( y) =ũ 2 (t( y)). In summary, we have derived a relation of the form u 1 (x ) + u 2 ( y) = c (x , y) for the location of the optical surface u, where u 1 (x ) = log(u), and c (x , y) is a non-quadratic cost function in optimal transport theory. Equation (11) has many solutions for u 1 (x ), u 2 ( y), and consequently for u(x ). We can find a unique solution by assuming that u 1 and u 2 are either c-convex or c -concave functions (p. 58, [1]). The surfaces u 1 and u 2 are c -convex if which we call the maximum solution, or c -concave if which we call the minimum solution.
For a continuously differentiable function, the c-convex/ concave functions u 1 , u 2 are Lipschitz continuous (p. 60, [1]), and the expression for the optical map y = m(x ) is implicitly given by the critical point of Eqs. (12b) or (13b), i.e., where ∇ x c is the gradient of c with respect to x , under the condition that the Jacobi matrix C, defined by

Research Article
is invertible. A sufficient condition for a maximum/minimum solution requires to be positive/negative semi-definite (SPD/SND), respectively, where D 2 u 1 is the Hessian matrix of u 1 , and D x x c is the Hessian matrix of c with respect to x . Hence, for a c-convex pair, we require tr( P) ≥ 0 and det( P) ≥ 0. On the other hand, for a c -concave pair, we need tr( P) ≤ 0 and det( P) ≥ 0. Note that P is a symmetric matrix. Differentiating Eq. (14) again with respect to x gives where Dm(x ) is the 2 × 2 Jacobi matrix of m with respect to x . Combining Eqs. (16) and (17) gives the matrix equation Note: in the cost function approach of this paper, we assume that each rayŝ refracts at the optical surface and always reaches the target plane. Hence, we circumvent the occurrence of total internal reflection (TIR), and the logarithmic cost function in Eq. (10) is never complex-valued. However, deriving the mapping via the vectorial law of refraction requires that the term under the square root 1 − n(1 − (n ·ŝ) 2 ) ≥ 0. Deriving an expression forn gives the TIR condition

B. Energy Conservation
By transferring the light from source to target, we require that all light from the source ends up at the target and energy is conserved, i.e., for an arbitrary set A ⊂ S 2 and image sett(A) ⊂ S 2 . Note that this image set corresponds to the far-field approximation. If we substituteŝ =ŝ(x ) andt =t( y) from Eq. (2), we can write Eq. (20) as wheref (x ) = f (φ, θ ), andg ( y) = g (ψ, χ ). We can derive that Substituting Eq. (22) and the mapping y = m(x ) into the energy conservation relation Eq. (21) gives where we omit the absolute value sign of the determinant and restrict ourselves to a positive Jacobian of the mapping. Using Eq. (18), we can rewrite Eq. (23) to the generalized Monge-Ampère equation , m(x ))) .

(24)
Our goal is to find a mapping y = m(x ) : X → Y for a particular source domain X and target domain Y. We also require global energy conservation, i.e., Eq. (21), with x (A) = X and y(t(A)) = Y under any mapping y = m(x ) : X → Y, implying that the light from the entire source is mapped to the target. For this reason, we cannot choose any combination of X , Y, f (x ), andg ( y). We eliminate this dependency by normalizing the source intensityf (x ) by the total intensity over X : and normalizing the target intensityg ( y) by the total intensity over Y: (25b) Note that we do not need to perform this normalization step if we have chosenf (x ) andg ( y) such that I (X ) = I (Y). Using Eq. (25), we rewrite the generalized Monge-Ampère equation in Eq. (24) to where we introduce F (x , m(x )) to denote the total righthand side. We define the corresponding transport boundary condition to Eq. (26a) as stating that all light from the boundary of the source X is mapped to the boundary of the target Y [7,33].

C. Polar Stereographic Coordinates
For an incoming cone-shaped beam, the source domain X is circular. Hence, we perform a change of coordinates for the stereographic variables x ∈ X on the source domain to polar stereographic coordinates ω = (ρ, ζ ) with the transformation where ρ ≥ 0 is the radial coordinate and 0 ≤ ζ < 2π the azimuth (angle with respect to positive x 1 axis). We maintain Cartesian stereographic coordinates for the target domain. We let u * (ω) = u(ŝ) and u * 1 (ω) = u 1 (x ) = log(u) and denote the optical map by y * = m * (ω) from source domain X * to target domain Y * . For ease of notation, we drop all asterisks, and in the following, we use the notation u 1 (ω) = log(u(ω)) for the surface and y = m(ω) for the mapping m : X → Y from the polar stereographic source domain to the Cartesian stereographic target domain, etc.
The mapping is implicitly given by , m(ω)), (28) where ∇ ω c is the gradient of c (ω, m(ω)) = c (x , m(x )) [the cost function in Eq. (11) rewritten in polar stereographic coordinates] with respect to ω, under the condition that the Jacobi matrix C, defined by is invertible. Following the same arguments as in the previous section, we derive the matrix equation in Eq. (18) as where and is SPD/SND with Note that D 2 u 1 (ω) and D ωω c (ω, m(ω)) are not Hessian matrices.

NUMERICAL METHOD
We will explain the broad steps of the least-squares algorithm using polar (stereographic) coordinates for the source domain X . The numerical method has previously been explained in detail using Cartesian coordinates in [7,32,33] and Cartesian stereographic coordinates in [35]. It works by first computing the mapping m in an iterative procedure whereafter the surface u is calculated also in a least-squares sense. The Monge-Ampère equation Eq. (37a) can be written as the matrix equation C Dm(ω) = P, where P(ω) is a symmetric, positive definite (SPD) matrix satisfying det ( P(ω)) =F (ω, m(ω))det(C(ω, m(ω)). We write m = m(ω) and minimize the functional under the constraint det( P) =F det(C). The norm used is the Frobenius norm. To impose the transport boundary condition, we minimize the functional where | · | denotes the l 2 -norm. Minimizing both functionals J I and J B together as a weighted average gives with 0 < α < 1 as weighting parameter for the interior domain and the boundary. We iterate starting from an initial guess m 0 and cost function matrix C(·, m 0 ): where the minimization steps are performed over the spaces m)det(C(·, m))}, (42b) i.e., spaces of (once and twice) continuously differentiable vector fields. After each iteration, we update the matrix C(·, m n ).
As initial guess m 0 , we map the circular source domain X centered around O enclosing X to a disc over Y. We assume the source X has radius ρ max and the bounding box of the target Y has rectangular shape [a min , a max ] × [b min , b max ]. In order to find a c -convex u 1 , we specify the initial guess m 0 = (m 0 1 , m 0 2 ) T for (ρ, ζ ) ∈ X , with 0 ≤ ρ ≤ ρ max and 0 ≤ ζ < 2π , as i.e., the initial guess is a disc of radius ρ max R centered around ((a min + a max )/2, (b max − b min )/2) T . The corresponding Jacobi matrix Dm 0 is SPD. After setting the initial guess m 0 , we perform the minimization steps in Eq. (41) and subsequently update C in every iteration. The minimization steps Eqs. (41a) and (41b) are performed point-wise and described in detail in [35]. The operations in these point-wise minimization steps are completely analogous to the Cartesian case. In contrast, the minimization step Eq. (41c) and the subsequent calculation of the lens surface cannot be performed point-wise and require more alterations when using polar stereographic coordinates. We will describe this in Sections 3.A and 3.B.
The minimizer is given by δ J [m, P, b](η) = 0 for all η ∈ M. Using Gauss' law and the fundamental lemma of calculus of variations [38], we obtain the coupled elliptic boundary value problem where the divergence operator works on the matrix A = C T C Dm as follows: Hence, we obtain two coupled elliptic equations for the components m 1 and m 2 of m. We solve for m using the finite volume method on a polar grid and use a closed control volume to incorporate the origin of the coordinate system, explained in detail in Appendix A.

B. Computation of the Lens Surface
Upon convergence of the iterative procedure for the mapping m, we can calculate the location of the optical surface from Eq. (28). In the ideal situation, ∇ ω c (ω, m(ω)) is a conservative vector field, so that there exists a u 1 that satisfies Eq. (28). Then we also have that C Dm = P in Eq. (30) (symmetric), which is most likely not satisfied after running the iterative procedure. In the numerical algorithm ∇ ω c (ω, m(ω)) is not conservative due to numerical errors. Using direct integration methods to determine u 1 is not unique and depends on the order of integration. Hence, we compute the generalized least-squares solution by minimizing the functional where u 1 (ω) = u 1 (x ), and we use short-hand notation for ∇ ω c (·, m) = ∇ ω c (ω, m(ω)) = ∇c (ω, y)| y=m(ω) .
We cannot perform this step point-wise, and analogous to the minimization procedure for m, we compute the first variation δ I [u 1 ](v) with respect to u 1 for v ∈ C 2 (X ) as The minimizer is given by δ I [u 1 ](v) = 0 for all v ∈ C 2 (X ). Once more, using Gauss' law and the fundamental lemma of calculus of variations [38], we obtain the boundary value problem This is a Neumann problem that has a unique solution up to an additive constant, and consequently a corresponding finite difference matrix with incomplete rank. We calculate a unique least-squares solution by using the Q R-decomposition of the finite difference matrix. The compatibility condition is satisfied for this Neumann problem up to discretization errors. Finally, we compute the location of the optical surface u by calculating u(ω) = e u 1 (ω) , since u 1 = log(u). We transform our polar stereographic source coordinates defined in Eqs. (1a) and (27) to Cartesian coordinates denoting (x , y , z) T = u(ω)ŝ and plot the reflector surface using cf. Eq. (2a), for all ω ∈ X .

NUMERICAL RESULTS
In this section, we apply the generalized least-squares algorithm to two example problems. First, we show that we can calculate the shape of a peanut lens typically used in roadlighting applications. Second, we compute the freeform surface that transforms light into an image on a projection screen in the far field. For both test cases, we choose an α = 0.1 as weighting parameter for the functional J , cf. Eq. (40). The laptop used for the calculations has an Intel Core i7-7700HQ CPU 2.80 GHz with 32.0 GB of RAM.

A. Peanut Lens
We compute a lens surface that transforms the light from a point source into a far-field target intensity distribution corresponding to a typical roadlighting profile. We consider a circular  source domain X = {(ρ, ζ ) ∈ R 2 | 0 ≤ ρ ≤ 1, 0 ≤ ζ < 2π } with a Gaussian light distribution: We determine I (X ) from Eq. (36) numerically using MATLAB's inbuilt quad2d . The target intensity g (ψ, χ ) is shown in Fig. 5 on the left, with zenith 0 ≤ ψ ≤ π with respect to the negative z axis and azimuth 0 ≤ χ < 2π extending into the far field on the street, as illustrated schematically in Fig. 4. The opening angle of the cone-shaped bundle is π rad = 180 • .
Changing to stereographic coordinates using Eq. (1b) transforms the left of Fig. 5 into the right figure. We determine I (Y) from Eq. (25b) by dividing the target in spherical coordinates into quadrants and approximating the integral of sin(ψ)g (ψ, χ ), since an infinitesimal area element of the unit sphere has size sin(ψ)dψdχ , by using MATLAB's inbuilt 2D integration method quad2d .
We use the least-squares algorithm to compute the optical map m and the lens surface. We discretize the source domain  using a 100 × 100 grid. The optical map, plotted using a coarsened version of the source grid, the lens surface, and convergence results are shown in Fig. 6. The total computation time of performing 200 iterations to calculate m is 24.9 s. The subsequent computation time of u 1 is 0.7 s. The left figure in Fig. 7 shows the target intensity converted to the local Cartesian coordinates (x , y ) on the street for a lamp a distance 6 m above the street. The conversion from the target intensity g (ψ, χ ) to local Cartesian coordinates is explained in detail in (p. 78, [33]). The right figure shows the ray trace results of our peanut lens, using a self-programmed ray-tracing algorithm in MATLAB performed with 1 million rays from a polar stereographic grid with quasi-random positions (quasi-Monte Carlo). For each ray, we calculate the interpolated value of the surface u and the surface normaln. Subsequently, we compute the direction of the light ray using the vectorial law of refraction and determine the corresponding grid cell on the target domain. The target on the street, cutoff at 15 m from the lamp in the x and y directions, is divided into a uniformly spaced grid. The white dot indicates the position of the lamp, and the black line marks the separation between the curb (left) and the street (right). The ray-trace result matches the desired profile on the street, with more light directed onto the street. As a metric for measuring the deviation between the original image and the raytrace result, we take the difference between the target intensity on the street (interpolated bilinearly onto the ray-tracing grid) and the ray-tracing irradiance. The mean absolute difference in intensity is 11%, and the maximum absolute difference is 20% (with the maximum value of the target intensity interpreted as the 100% value). The correlation coefficient between the intensities is 0.56, and the energy efficiency of the ray tracer is 34%, due to the cutoff at 15 m. TIR does not occur, and we can verify that condition Eq. (19) is not satisfied for all rays traced.

B. Picture on a Projection Screen in the Far Field
We challenge our numerical algorithm to compute a lens surface that converts the light from a point source into a far-field target intensity distribution corresponding to a picture. The source domain in polar stereographic coordinates is given by the circle X = {(ρ, ζ ) ∈ R 2 | 0 ≤ ρ < 0.5, 0 ≤ ζ < 2π } and has a uniform light distributionf (ω) = 1. Using Eq. (36), we get I (X ) = 4 5 π . The refracted rays are projected on a screen in the far field, parallel to the x y plane at 1 m distance above the lens. The required illumination L(ξ, η) [lm/m 2 ], with (ξ, η) the local Cartesian coordinates on the projection screen, is derived from the gray scale values of a photograph of the first author's dog (Fig. 8). The target distributiong ( y) is a deformation of the illuminance L(ξ, η) and shown in Fig. 9; the conversion from L(ξ, η) tog ( y) is explained in detail in (p. 78, [33]).
The gray scale values of the picture prescribe the illuminance. However, the conversion from the colored image to gray scale values creates black regions in the target distribution for which g ( y) = 0. To avoid division by 0 in the right-hand side of Eq. (26a), we increase values ofg ( y), which are below a threshold of 5% of its maximum value to this threshold. We determine I (Y) by dividing the region on the projection screen into quadrants with the pixels of the picture as corners and approximating the integral of L(ξ, η) using MATLAB's inbuilt 2D integration method quad2d .
We use the least-squares algorithm to compute the optical map m and the lens surface. We discretize the source domain using a 500 × 500 grid. The optical map m, plotted using a coarsened version of the source grid, the lens surface, and convergence results are shown in Fig. 9. The error J I for the interior converges to a larger value than the error J B for the boundary, as the finite grid is unable to capture the finest details of the dog. The total computation time of performing 1000 iterations to calculate m is 4347.8 s. The subsequent computation time of u 1 is 247.1 s.
Subsequently, we validated the resulting lens image using ray tracing. We traced 3000 × 3000 rays with positions taken from a regular dense polar stereographic grid. The resulting target illuminance L(ξ, η) is plotted in Fig. 8. The ray-trace image closely resembles the original picture, showing details such as the hairs on the dog's coat. The mean absolute difference in intensity is 22% and the maximum absolute difference is 100% with the maximum difference occurring at the top and bottom boundaries. The correlation coefficient between the intensities is 0.87, and the energy efficiency of the ray tracer is 99.8%. TIR does not occur, and we can verify that condition Eq. (19) is not satisfied for all rays traced. The ray-trace image is not perfectly rectangular since the mapping does not perfectly align to the top and bottom boundaries of the target in Fig. 9, due to the iterative nature of the numerical algorithm. This creates small errors at the boundaries of the lens surface.

CONCLUSION
In this paper, we presented a method to compute freeform lens surfaces that convert the light from a point source into a far-field target. We used a least-squares approach to solve the generalized Monge-Ampère equation of the system, which we derived by combining the cost function in optimal transport theory and energy conservation. We introduced a coordinate transformation for the source domain to polar stereographic coordinates. This transformation facilitates the use of circular source domains.
We tested the numerical algorithm for two lighting applications. We could compute a peanut-shaped lens surface typically used in roadlighting. Moreover, we applied the numerical algorithm to a complicated test problem involving a picture on a screen in the far field. Other methods currently available that produce this level of detail with a freeform lens are [17,30], although [17] computes a near-field example, and [30] uses an iterative wavefront tailoring approach.
In future work, we would like to consider extended light sources and double freeform surfaces. Also, we would like to take into account scattering phenomena and abberation theory, and explore the applicability of the numerical algorithm in other fields of science and engineering relating to optimal transport theory.

APPENDIX A: FINITE VOLUME METHOD IN POLAR STEREOGRAPHIC COORDINATES
For the minimization procedure of m, we derived the coupled elliptic boundary value problem in Eq. (45). This is a system of two coupled elliptic equations with Robin boundary conditions for the components of m 1 and m 2 of m.
Using the cost function matrix C(ω, m(ω)) in Eq. (29) and In a polar basis, we can rewrite Eq. (45) as with Research Article r p = r 11 r 12 r 21 r 22 . (A3b)

B. Inner Boundary
We consider the grid adjacent to the origin, i.e., grid points (ρ 1 , ζ j ), with 1 ≤ j ≤ N ζ . We denote the values of m 1 (ρ 0 , ζ j ) = m 1 (ρ 0 ) and m 2 (ρ 0 , ζ j ) = m 2 (ρ 0 ) at the origin by m 1,O and m 2,O , as the vector m O . By considering a control volume as a disc around the origin, we will derive expressions for m 1,O and m 2,O and remove them from the main equations. We consider the disc D h ρ /2 , centered at the origin with radius h ρ /2, as the control volume, as shown in Fig. 12.

(A18)
We can rewrite Eq. (A17) by approximating the integrals to second-order accuracy to and Eq. (A18) to Using Eqs. (A19) and (A20), we can solve for m 1,O and m 2,O . For the grid points (ρ 1 , ζ j ), with 1 ≤ j ≤ N ζ , we can also elaborate on the main discretization Eq. (A13) for both Eqs. (A4a) and (A4b) and substitute the expressions found for m 1,O and m 2,O into the discretization. Note that after every iteration, we update the matrix C and also calculate c 1,O and c 2,O at the origin. We use this to be able to calculate values at the interface points via linear interpolation, such as |c 1 | 2 1 2 , j .