Ray mapping approach for the efficient design of continuous freeform surfaces

The efficient design of continuous freeform surfaces, which transform a given source into an arbitrary target intensity, remains a challenging problem. A popular approach are ray-mapping methods, where first a ray mapping between the intensities is calculated and in a subsequent step the surface is constructed. The challenging part hereby is the to find an integrable mapping ensuring a continuous surface. Based on the law of reflection/refraction and the well-known integrability condition, we derive a general condition for the surface and ray mapping for a collimated input beam. It is shown that in a small-angle approximation a proper mapping can be calculated via optimal mass transport. We show that the surface can be constructed by solving a linear advection equation with appropriate boundary conditions. The results imply that the optimal mass transport mapping is approximately integrable over a wide range of distances between the freeform and the target plane and offer an efficient way to construct the surface by solving standard integrals. The efficiency is demonstrated by applying it to two challenging design examples.

The efficient design of continuous freeform surfaces, which transform a given source into an arbitrary target intensity, remains a challenging problem. A popular approach are ray-mapping methods, where first a ray mapping between the intensities is calculated and in a subsequent step the surface is constructed. The challenging part hereby is the to find an integrable mapping ensuring a continuous surface. Based on the law of reflection/refraction and the well-known integrability condition, we derive a general condition for the surface and ray mapping for a collimated input beam. It is shown that in a small-angle approximation a proper mapping can be calculated via optimal mass transport. We show that the surface can be constructed by solving a linear advection equation with appropriate boundary conditions. The results imply that the optimal mass transport mapping is approximately integrable over a wide range of distances between the freeform and the target plane and offer an efficient way to construct the surface by solving standard integrals. The efficiency is demonstrated by applying it to two challenging design examples.

I. INTRODUCTION
In recent years much progress has been made in the field of freeform surface design without the assumption of symmetries . The goal of these design methods is the solution of the so called inverse problem of nonimaging optics. This means that for given arbitrary source and target intensities I S (x, y) and I T (x, y) one or more freeform surfaces have to be calculated, which map the intensities onto each other. Especially the design of continuous freeform surfaces, on which we will concentrate on in the following, is a challenging problem and of great interest for practical applications. The first successful method at calculating a continuous freeform surface utilizing a complex target intensity was developed by Ries and Muschaweck [1], but unfortunately the numerical method was not published [1]. Their approach is able to handle the far field design problem for single freeform surfaces illuminated by a point source [1]. Nowadays, many other methods have been developed by different research groups. A quite popular approach are the so-called Monge-Ampre methods [2][3][4][5][6][7][8]. They are based on the modelling of the design problem by a nonlinear partial differential equation of Monge-Ampre type and solving it with sophisticated numerical techniques. These methods are able to handle the design problem in the far field for intensity control of point sources [2,6,8] and collimated beams [4,6,7] as well as intensity and phase control with double freeform surfaces [5]. Another popular approach for the single freeform surface design with point sources is the supporting ellipsoids method developed by Oliker [9]. With this method a freeform mirror is constructed by putting a point source in the focal point of a unification of ellipoids, whereby every ellipsoid has a different position of the second focal * christoph.boesel@uni-jena.de point on the target plane to build the required intensity pattern. The challenge of this method is the calculation of a smooth freeform surface by the unification of the ellipsoids. Therefore it was further developed by other research groups [10,11] and generalized to calculate freeform lenses [13]. It can handle the far field as well as the near field design problem. Also quite often used are ray-mapping techniques [10,12,[14][15][16][17][18][19][20][21], which are frequently based on the calculation of a ray mapping between the source and the target intensity and a subsequent construction of the freeform surface [10,12,14,15,[18][19][20][21]. The aim and challenging part of these methods is to find an integrable ray mapping, which allows the calculation of a continuous surface. The approach we will concentrate on is a subgroup of the ray-mapping techniques, which are called optimal mass transport (OMT) methods [18][19][20][21]. These gained some interest in recent years and are partly based on the mathematical concept of optimal mass transport as explained in the following paragraph. They can handle the design problem of a single and double freeform surfaces for intensity control [18,19] as well as double freeform surfaces for intensity and phase control [20,21]. The connection between mass transport and freeform design was also discussed in [22]. This approach for the freeform surface design consists of two separate steps. In the first step a ray mapping between the given input and output intensities is calculated via OMT. In the second step the freeform surface is constructed with the help of the law of refraction/reflection and the well-known integrability condition, which ensures the continuity of the surface. The difference between the OMT methods mentioned above is the second step. In first approach by Buerle et al. the freeform is constructed by an optimization procedure [18,19], while the second approach by Feng et al. uses a simultaneous point-by-point construction method to design a double freeform surface [20,21].
Since these attempts seem to be quite successful but do not give theoretical insights about the integrability of the OMT map, we want to clarify this point in our work for a single freeform illuminated by a collimated beam. This will be done by deriving a condition for an integrable map and showing that it can be fulfilled (approximately) by the OMT map. Based on our findings, we present an efficent and easy-to-implement numerical freeform surface construction technique differing from the previously published OMT methods [18][19][20][21]. To do so this paper is structured as follows.
In section II after a short introduction to the OMT and a presentation of its basic properties, we will derive from the law of reflection/refraction and the integrability condition a general condition for an integrable ray mapping and its corresponding surface for collimated input beams. It will be shown that in a small-angle approximation this condition can be fulfilled by using an OMT mapping and therefore the freeform surface design process indeed decouples into the two steps described above. Thereby it will be shown that the freeform surface can be constructed from a linear advection equation with appropriate boundary conditions. In section III we then argue that the advection equation for the freeform construction can be solved by simple integrations, which is different to the OMT freeform design methods mentioned above and implies the approximate integrability over a wide range of freeform-target plane distances. The efficiency of this approach is then demonstrated in section IV by applying it to two challenging design examples, followed by a discussion of our results in section V.

A. Optimal mass transport
The problem statement of OMT, also called the Monge-Kantorovich problem, is as follows: two positive density functions I S (x) and have to be mapped onto each other according to the Jacobi equation with a smooth, bijective mapping u(x). If M is defined as the set of mappings fulfilling equation (2), we are searching for a mapping minimizing the transport cost according to the Kantorovich-Wasserstein distance whereby inf u∈M denotes the mapping for which the integral has its minimal value. This mapping, which is defined by (1), (2) and (3), has the useful property that it is unique [23] and it is characterized by a vanishing curl [24]. The latter property will be important for our findings in the next subsection.
In the special case of freeform surface design considered here, the densities I S (x) and I T (x) correspond to the source and target intensities with the units W · m −2 (see Fig. 1). Therefore equation (1) describes a global and equation (2) a local energy conservation. For the numerical examples in section IV, we have implemented the OMT method developed by Sulman et al. [25]. It provides a good compromise between an easy implementation and an efficient mapping calculation and is thus sufficient for our test purposes. The result of this design process step is therefore the mapping u(x, y), but we have to keep in mind, that it is not obvious at this point, whether the OMT mapping is integrable and if it is, for which lens-target distances (see Fig. 1) this is the case. This point will be clarified in the next subsection.
B. Freeform surface construction

Ray-mapping condition and surface equation
In the following, we want to derive a differential equation for the direct calculation of a freeform surface for a given general ray mapping u(x, y). Our derivations will lead to a condition an integrable mapping and its corresponding freeform surface have to fulfill. As we will show, this condition can be fulfilled approximately over a wide range of lens-target distances by the OMT mapping defined in the previous subsection. To do so, two basic equations are considered. On the one hand, for an incoming beam described by the ray direction vector field s 1 and the refracted vector field s 2 , the law of refraction with the refractive indices n 1 of the lens and n 2 of the surrounding medium, has to be fulfilled. On the other hand, we want to ensure the continuity of the surface z(x, y) by the well-known integrability condition [1] n · (∇ × n) = 0.
Since the collimated beam s 1 as well as s 2 can be expressed in terms of the unknown freeform surface z(x, y) and the given ray mapping u(x, y) (see Fig. 1): the equations (4) and (5) represent a differential equation for z(x, y). Plugging (4) into equation (5) the integrability condition can be written in the form (see Appendix 1. a) For the given input intensity IS(x, y) and output intensity IT (x, y) the OMT ray mapping u(x, y) is calculated. The mapping defines a vector field s3 between the source plane z = 0 and the target plane z = zT . ∂ΩS and ∂ΩT are the source and target intensity boundaries, respectively. b) In the second process step the freeform surface z(x, y) is constructed in a way that it is continuous and redirects the incoming collimated beam, described by the vector field s1, according to the given OMT map. Because of energy conservation, the boundary of the freeform corresponds to the shape of the source intensity ∂ΩS. A) which holds for a general ray mapping u(x, y). Equation (7) is organized in a way that only the left-hand side (LHS) depends on the derivatives of z(x, y). Equation (7) takes a more familiar form by inserting the vector fields (6), which leads to with the velocity field v = (u(x, y) − Id) ⊥ , the identity vector Id = (x, y) T and ∇ ⊥ = (−∂ y , ∂ x ). This equation is a semilinear two dimensional advection equation, whereby the unknown surface z(x, y) corresponds to conserved transport quantity and the right-hand side (RHS) to a source term. In principle one could try to solve equation (8) after applying suitable boundary conditions, but as it will be demonstrated, it is not appropriate for finding a continuous freeform surface. This can be seen easily by considering the condition that the normalized vector field (4) has to be equal to the gradient of the surface z(x, y): (9) Plugging this relation into the LHS of (8), we get v∇z(x, y) ≡ 0, which can only be fulfilled if the RHS vanishs: The importance of this condition is due to the fact that is has to hold for every integrable ray mapping. Since v ⊥ (u(x, y) − Id) and therefore (u(x, y) − Id)||∇z(x, y) it reflects the nature of the law of refraction, that according to the definitions (6) and (9) the vectors s 1 (x, y), s 2 (x, y) and n(x, y) have to lie in the same plane. We now know, that the source term of (8) has to vanish, but we are still left with the question, if we can find a way to fullfil condition (10). The main task is obviously to find a ray mapping for which relation (10) holds, which is nontrivial, since it couples the mapping with the unknown function z(x, y). But if we use the OMT mapping, it follows from its vanishing curl that ∇v = 0. If we use in addition to that, the small-angle approximation between the surface normal n(x, y) and the outgoing ray s 2 (x, y), we see that the condition (10) can be fulfilled approximately by using an OMT map. Because of the fact that n · s 2 is locally proportional to z T − z(x, y) in contrast to the RHS, (11) can be interpreted as a far field approximation. This also implies, that the integrability of OMT map is only asymptotically exact. Hence for an OMT mapping and the small-angle approximation, we get our final equation v∇z(x, y) = ∇(v · z(x, y)) = 0, v = (u(x, y) − Id) ⊥ (12) which has to be solved to get the required freeform surface z(x, y).

Boundary conditions
If we want to solve a linear advection equation like (12), we have to know the function z(x, y) on the inflow part of the boundary, where the velocity field points into the integration area Ω S [26]. Because of energy conservation, this area is defined as Ω S := {(x, y) ∈ R 2 | I S (x, y) = 0} and therefore its inflow part by ∂Ω S− := {(x, y) ∈ ∂Ω S |v· r < 0} with the outward boundary normalr. Together with this boundary condition equation (12) has at most one solution [26]. In our case the boundary conditions can be deduced in the following way. First, we have to realize, that for an incoming collimated beam the boundary of the freeform surface can only determine the tangential deflection of a ray refracted at the boundary. The normal deflection is determined by the inner parts of the surface z(x, y). Therefore it makes sense to parameterize the boundary ∂Ω S by a parameter s and define the local coordinate system at each point of the boundary by the vectors(see Fig. 2 (13) Since z(s), s ∈ ∂Ω S only determines the tangential deflection, it is sufficient to consider the law of refraction (4) in the tangential plane spanned byt(s) and e z . Hence, we can interpret the boundary value calculation as a two dimensional problem which allows us to derive a differential equation for the boundary values z(s) by the two dimensional equivalent of equation (9) where path length l(s) was introduced for dimensional reasons. From this we get: which reduces in the far field to whereby the position of the surface in space compared to the target plane is fixed by integration constant. Since (12) itself can be interpreted as a far field approximation, as explained above, equation (16) seems more suitable for our purposes. It provides us with a simple way of calculating the boundary values. The only degree of freedom left is the integration constant. Equation (16) will build the basis of the numerical algorithm for solving equation (12) presented in the next section. At the end of this section, we want to note that (12) and (16) can be derived analogously for freeform mirrors by replacing (4) by the law of reflection and keeping in mind that z T < z(x, y).

III. NUMERICAL ALGORITHM
We could solve equation (12) by standard computational fluid dynamic approaches, like finite volume methods, which are appropriate for the numerical treatment of linear advection equations. Based on the nature of equation (16), a different approach is proposed in the following. Considering (16), we recognize that the boundary values are calculated by the velocity field itself. This is in contrast to the usual fluid dynamical framework and allows us to separate Ω S into an arbitrary number of subareas Ω S,i for each of which we can calculate the boundary values by (16). Therefore, the freeform surface can be calculated on each subarea Ω S,i and the solution on Ω S by their unification. This implies that the freeform surface can be constructed by an integration of the OMT map along arbitrary paths on Ω S , which characterizes the integrability of ray mappings [10,12]. Thus according to (12) and (16) the OMT map is approximately integrable as long as (11) holds. Hence, the most convenient way to get the solution of (12) seems to be a line-by-line integration of (16), which along the x-and y-direction is equivalent to equation (9) in a far field approximation. One possible way of integrating (16) is shown as an example in Fig. 3. Thus, only the integration constant of one integral has to be fixed from which the others follow automatically. The proposed approach has the useful feature that we do not need to parameterize the boundary ∂Ω S , which allows the calculation of freeform surfaces with complex boundary shapes. The efficiency of the line-by-line integration approach is shown in the next chapter for two challenging design examples.

IV. EXAMPLES
To show the efficiency of the algorithm, we want to apply it to two design examples. In the first one, we will calculate a freeform lens that maps a collimated beam of uniform intensity on the logo of the Institute of Applied Physics (IAP) in Jena with a resolution of 500 x 500 pixels (see Fig. 4). It shows strong intensity gradients between the letters and the background. To omit a division by zero within the implemented OMT algorithm [25] we have to use a background intensity I > 0 for the input and output intensities. For an appropriate speed of convergence the background intensity is set to 20 per cent of the maximum intensity. The second example, which shows smooth intensity variations and a lot of details, is the well-known picture of Lena with a resolution of 500 x 500 pixels (see Fig. 4). Since the freeforms were calculated by integrating equation (16), the specific characteristics of both pictures do not have any influence on the speed of the lens construction step explained in section III, but they increase the mapping calculation time. The pictures both have a square format. Hence, it is convenient to integrate first along the upper side of the square region and then line-by-line in the orthogonal direction with the starting values given by the first integration. Therefore we have to solve 501 integrals, which took in both cases less than one second in MATLAB on an Intel Core i3 at 2.4Ghz with 16GB RAM. This time has to be added up with the mapping calculation time, which strongly depends on the implemented method and the specific features of the picture.
The integration constant on the upper left side of the integration area was chosen to be z 0 = 1. Values of n 1 = 1.5 and n 2 = 1 were used for the refractive indices, and the source-target distance was chosen to be z T = 5. The input and output beam as well as the freeform lens had side length of one. Since every spatial value is normalized the results are scalable. At this point we want to note that according to (16) the validity of the approximation (11) can simply be checked by scaling the numerical results with 1/z T , which of course has to be done for each intensity and configuration individualy. For our examples the quality of the illumination patterns produced by the raytracing did not change significantly even for distances between the lens and the target plane smaller than the side length of the lens. In both cases the calculated freeform lens was imported as a grid sag surface into ZEMAX to verify our results by a raytracing simulation. The imported lens data are interpolated by ZEMAX automatically. The results can be seen in Fig. 4.

V. CONCLUSION
We presented an efficient numerical method for the single freeform surface design for the shaping of collimated beams. It is based on the derivation of the condition (10) for integrable mappings by combining the law of reflaction/refraction and the well-known integrability condition in a suitable way. We showed that the condition can be fulfilled in the small-angle approximation (11) by using a ray mapping calculated by optimal mass transport. Equation (11) therefore represents a quantitative estimate for the applicability of the OMT map. This serves as a theoretical basis for the decoupling of the design process into two separate steps: the calculation of the OMT mapping and the construction of the freeform surface with a linear advection equation. On the basis of the finding of appropriate boundary conditions for the advection equation, we presented a simple numerical algorithm for the surface construction by solving the standard integrals (16), which differs from the previously published OMT freeform design methods. Besides its simplicity, accuracy and quickness, a useful feature of the construction technique is the independence of the freeform boundary shape (see Fig. 3). By using a proper method for the calculation of the OMT mapping, this allows for example the calculation of a freeform for disconnected intensities or source and target intensties with two concave boundaries, which is in general a nontrivial problem, but important for applications. The results imply the approximate integrability of the OMT map over a wide range of freeform-target plane distances and gives the OMT freeform design methods a theoretical basis for collimated input beams. Especially the condition (10) is interesting, also for other ray mapping techniques than OMT methods, since it has to hold  (16) at the upper left side was z0 = 1 and the source-target distance zT = 5. The imported numerical data was interpolated automatically by ZEMAX. e), f): intensity pattern from a ZEMAX raytracing with 4 · 10 7 rays.
for every integrable mapping. Hence, it opens up new possibilities for nearfield calculations as well as generalizations to it for point sources and double freeforms, which are currently under work.