## Abstract

The efficient design of continuous freeform surfaces, which maps a given light source to an arbitrary target illumination pattern, remains a challenging problem and is considered here for collimated input beams. A common approach are ray-mapping methods, where first a ray mapping between the source and the irradiance distribution on the target plane is calculated and in a subsequent step the surface is constructed. The challenging aspect of this approach is to find an *integrable* mapping ensuring a *continuous* surface. Based on the law of reflection/refraction and an 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 - a mathematical framework for the calculation of a mapping between two positive density functions. We show that the surface can be constructed by solving a linear advection Eq. 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, which shows the ability of the presented approach to handle target illumination patterns with steep irradiance gradients and numerous gray levels.

© 2016 Optical Society of America

## 1. Introduction

In recent years much progress has been made in the field of freeform surface design without the assumption of symmetries for nonimaging [1–21] and imaging applications [22–24]. In nonimaging optics, which is considered here, the goal of these design methods is the solution of an inverse problem: for a given light source of arbitrary shape and a given target illumination pattern, one or more freeform surfaces have to be calculated which map the source and the target onto each other. We will concentrate in the following on the design of *continuous* freeform surfaces.

The first successful method at calculating a continuous freeform surface utilizing a complex target illumination pattern was developed by Ries and Muschaweck [1], but unfortunately the numerical method was not published [1]. Their approach is able to handle the near field and far field design problem for single freeform surfaces illuminated by a point source [1].

Many other methods have since been developed by different research groups. A common approach are the Monge-Ampère methods [2–8]. They are based on the modelling of the design problem by a nonlinear partial differential Eq. of Monge-Ampère type and solving it with complex numerical techniques. These methods are able to handle the single freeform design problem in the far field for point sources [2, 6, 8] and collimated beams [4, 6, 7] as well as ir-radiance and phase control with double freeform surfaces [5], but are in some cases limited to simple boundary shapes of the freeform [2–6].

Another common method 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 point on the target plane to build the required irradiance 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 problems.

Ray-mapping techniques [10, 12, 14–21], which are based on the calculation of a ray mapping between the light source and the target irradiance and a subsequent construction of the freeform surface [10, 12, 14, 15, 18–21], are also often used. 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–21]. These gained some interest in recent years and utilize 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 irradiance control [18,19] as well as double freeform surfaces for irradiance and phase control [20, 21]. The connection between mass transport and freeform design was also discussed in [25].

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 irradiances is calculated via OMT. In the second step the freeform surface is constructed with the help of the law of refraction/reflection and the commonly used integrability condition [1], which ensures the continuity of the surface.

The difference between the OMT methods mentioned above is the second step. In first approach by Bäuerle *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, this point is clarified in our work for a single freeform illuminated by a collimated beam by deriving a condition for an integrable map which is fulfilled approximately by the OMT map. Based on this, we present an efficient numerical freeform surface construction technique differing from the previously published OMT methods [18–21]. To do so this paper is structured as follows.

In section 2 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. Thus it will be shown that the freeform surface can be constructed from a linear advection Eq. with appropriate boundary conditions. In section 3 we then argue that the advection Eq. 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 4 by applying it to two design examples with complex irradiance distributions, followed by a discussion of the results in section 5.

## 2. Design method

#### 2.1. 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

*I*(

_{T}**x**) with

**u**(

**x**). If

*M*is defined as the set of mappings fulfilling Eq. (2), we are searching for a mapping minimizing the transport cost according to the Kantorovich-Wasserstein distance

*inf*

_{u}_{∈}

*denotes the mapping for which the integral is minimal. This mapping defined by the Eqs. (1), (2) and (3) is unique [26] and characterized by a vanishing curl [27]. The latter property will be important for section 2.2.*

_{M}In the particular case of freeform surface design which is considered here, the densities *I _{S}*(

**x**) and

*I*(

_{T}**x**) correspond to the source and target irradiance distributions with the units

*W · m*

^{−2}(see Fig. 1) and the integration in Eq. (1) is done over the corresponding apertures. Therefore Eq. (1) and Eq. (2) describe a local and global energy conservation, respectively.

For the numerical examples in section 4, we have implemented the OMT method developed by Sulman *et al.* [28]. It allows 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.

#### 2.2. Freeform surface construction

### 2.2.1. Ray-mapping condition and surface Eq

In the following, a differential Eq. for the direct calculation of a freeform surface for a given *general* ray mapping **u**(*x,y*) is derived. The derivations will lead to a condition which an *integrable* mapping and its corresponding freeform surface have to fulfill. As it will be demonstrated, this condition can be fulfilled approximately over a wide range of lens-target distances by the OMT mapping defined in the section 2.1.

To do so, two basic Eqs. are considered. On the one hand, for an incoming beam described by the ray direction vector field **s**_{1}(*x,y*) and the refracted vector field **s**_{2}(*x,y*), the law of refraction

*n*

_{1}of the lens and

*n*

_{2}of the surrounding medium, needs to be fulfilled, whereby the ”hat” denotes a unit vector field. On the other hand, the continuity of the surface

*z*(

*x,y*) needs to be ensured by the integrability condition [1]

Since the collimated beam **s**_{1} and **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):

*z*(

*x,y*). Plugging Eq. (4) into Eq. (5) the integrability condition can be written in the form (see Appendix 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

**v**= (

**u**(

*x,y*)−

**Id**)

^{⊥}, the identity vector

**Id**= (

*x,y*)

*and ∇*

^{T}^{⊥}= (−∂

*, ∂*

_{y}*). This Eq. is a semilinear two dimensional advection Eq., whereby the unknown surface*

_{x}*z*(

*x,y*) corresponds to conserved transport quantity and the right-hand side (RHS) to a source term.

In principle one could solve Eq. (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 by considering the condition that the normalized vector field (4) needs to be equal to the gradient of the surface *z*(*x,y*):

Plugging this relation into the LHS of Eq. (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 of whether there is a way to fulfill 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 the OMT mapping defined in section 2.1 is used, it follows from ∇ × **s**_{3} *≡* 0 that ∇**v** = 0. By additionally using the small-angle approximation

**n**(

*x,y*) and the outgoing ray

**s**

_{2}(

*x,y*), the condition (10) can be fulfilled approximately by using an OMT map. Because

**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, Eq. (8) reduces to

*z*(

*x,y*).

### 2.2.2. Boundary conditions

In order to solve the linear advection Eq. (12), the function *z*(*x,y*) needs to be calculated on the inflow part of the boundary where the velocity field points into the integration area Ω* _{S}* [29]. Because of energy conservation, this area is defined as
${\mathrm{\Omega}}_{S}:=\overline{\left\{(x,y)\in {\mathbb{R}}_{2}|{I}_{S}(x,y)\ne 0\right\}}$ and therefore its inflow part by
$\partial {\mathrm{\Omega}}_{S-}:=\left\{(x,y)\in \partial {\mathrm{\Omega}}_{S}|\mathbf{v}\cdot \widehat{\mathbf{r}}<0\right\}$ with the outward boundary normal
$\widehat{\mathbf{r}}$. Together with this boundary condition Eq. (12) |has at most one solution [29].

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 we 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)

Since *z*(*s*) with *s* ∈ *∂*Ω* _{S}* only determines the tangential deflection, it is sufficient to consider the law of refraction (4) in the tangential plane spanned by
$\widehat{\mathbf{t}}(\mathbf{s})$ and

**e**

*. Hence, we can interpret the boundary value calculation as a two dimensional problem which allows us to derive a differential Eq. for the boundary values*

_{z}*z*(

*s*) by the two dimensional equivalent of Eq. (9) projected on the $\widehat{\mathbf{t}}(\mathbf{s})-{\mathbf{e}}_{z}$ - plane

*l*(

*s*) was introduced for dimensional reasons. From this we get:

*z*

_{0}. Since Eq. (12) itself can be interpreted as a far field approximation, as explained above, Eq. (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 Eq. (12) presented in the next section.

At the end of this section, we want to note that the Eqs. (12) and (16) can be derived analogously for freeform mirrors by replacing Eq. (4) by the law of reflection and keeping in mind that *z _{T}* <

*z*(

*x,y*).

## 3. Numerical algorithm

We could solve Eq. (12) by standard computational fluid dynamic approaches, like finite volume methods, which are appropriate for the numerical treatment of linear advection Eqs.. Based on the nature of Eq. (16), a different approach is proposed in the following.

Considering Eq. (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 Ω

*for each of which we can calculate the boundary values by Eq. (16). Therefore, the freeform surface can be calculated on each subarea Ω*

_{S,i}*and the solution on Ω*

_{S,I}*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 the Eqs. (12) and (16) the OMT map is approximately integrable as long as Eq. (11) holds.*

_{S}Hence, the most convenient way to get the solution of Eq. (12) seems to be a line-by-line integration of Eq. (16), which along the *x*- and *y*-direction is equivalent to Eq. (9) in a far field approximation. One possible way of integrating Eq. (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.

Therefore the design approach can be summarized as follows:

- Calculate a mapping
**u**(*x,y*) between the light source and target illumination pattern represented respectively by the densities*I*(_{S}*x,y*) and*I*(_{T}*x,y*) by using an optimal mass transport algorithm (for example from [28]).

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 design examples.

## 4. 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 irradiance on the logo of the Institute of Applied Physics (IAP) in Jena with a resolution of 500 × 500 pixels (see Fig. 4(a)). It shows strong irradiance gradients between the letters and the background. To omit a division by zero within the implemented OMT algorithm [28] we have to use a background irradiance *I* > 0 for the input and output irradiance distributions. For an appropriate speed of convergence the background irradiance is set to 20 per cent of the maximum irradiance, but we want to emphasize that the method by Sulman et al. [28] is also able to handle lower background irradiances. The second example, which shows smooth irradiance variations and a lot of details, is the picture of Lena with a resolution of 500 × 500 pixels (see Fig. 4(b)).

Since the freeforms were calculated by integrating the OMT mapping according to (16) along straight lines (see Fig. 3) every pixel corresponds to one grid point and the specific characteristics of both pictures do not have any influence on the speed of the lens construction step explained in section 3, 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 R2015b 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 from (16) on the upper left side of the integration area at (*x,y*) = (−0.5,−0.5) was chosen to be *z*_{0} = 0 arb. unit. 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 arb. unit. The input and output beam as well as the freeform lens had side length of one arb. unit. Since every spatial value is normalized the results are scalable.

At this point we want to note that according to Eq. (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 irradiance 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 13 to verify the results by a raytracing simulation. The imported lens data were interpolated by ZEMAX with a linear interpolation algorithm, since the bicubic spline interpolation in ZEMAX led to oscillations in the illuminations patterns. The results from the raytracing can be seen in Figs. 4(e) and 4(f) and a comparison between the prescribed irradiance patterns and the simulation results in Figs. 4(g) and 4(h). Especially in the regions with steep gradients, we can observe a reduced quality of the irradiance distribution from the raytracing simulation. This is mainly due to the limited precision of the OMT map algorithm and the linear interpolation of the lens data by ZEMAX. According to Eq. (11) the quality of the simulation results can also be improved by choosing a larger distance between the freeform surface and the target plane. According to ZEMAX the efficiency was about 99.999% for both examples.

We want to point out that the presented approach is not restricted to designs where the target irradiance distribution has the same size and shape as the freeform surface. The target size can be scaled by scaling the mapping **u**(*x,y*) and/or the target position can be shifted by shifting the components of the mapping since the curl of **s**_{3} is invariant under these transformations and therefore remains zero. It is also possible to calculate a second freeform lens for the examples in this section by mirroring the mapping at the point of origin in the *x*-*y*-plane, since the curl of **s**_{3} is again invariant.

## 5. 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 Eq.. Based on finding appropriate boundary conditions for the advection Eq., 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 speed, 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 target irradiance distributions or source and target irradiances 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 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.

## Appendix A

We want to derive Eq. (8) from the law of refraction (4) and the integrability condition (5). Since the curl of the incident field ${\mathbf{s}}_{1}$ vanishs, plugging Eq. (4) into Eq. (5) gives

With

Inserting Eq. (4) and using ${\widehat{\mathbf{s}}}_{1}=(0,0,1)$ and ${\widehat{\mathbf{s}}}_{2}\xb7({\mathbf{s}}_{2}\times \dots )=0$ leads to

**s**

_{2}=

**s**

_{3}

*−*

**s**

_{1}it follows

Using the definition **v** = (*−*(**s**_{2})* _{y}*,(

**s**

_{2})

*) and (6) the terms in Eq. (21) can be written as*

_{x}## Acknowledgments

The authors thank M. Esslinger, and M. Tessmer for valuable discussions, R. Hambach and S. Schmidt for valuable discussions and comments on the manuscript, C. Liu and D. Lokanathan for the help with the ZEMAX implementation and D. Musick for the spelling and grammar check. We also acknowledge the Federal Ministry of Education and Research Germany for financial support through the project KoSimO (FKZ:031PT609X).

## References and links

**1. **H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A **19**(3), 590–595 (2002). [CrossRef]

**2. **R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “Freeform illuminaton design: a nonlinear boundary problem for the elliptic Monge-Ampère equation,” Opt. Lett. **38**(2), 229–231 (2013). [CrossRef] [PubMed]

**3. **R. Wu, K. Li, P. Liu, Z. Zheng, H. Li, and X. Liu, “Conceptual design of dedicated road lighting for city park and housing estate,” Appl. Opt. **52**(21), 5272–5278 (2013). [CrossRef] [PubMed]

**4. **R. Wu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “A mathematical model of the single freeform surface design for collimated beam shaping,” Opt. Express **21**(18), 20974–20989 (2013). [CrossRef] [PubMed]

**5. **Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge-Ampère equation method,” Opt. Comm. **331**, 297–305 (2014). [CrossRef]

**6. **R. Wu, Y. Zhang, M. M. Sulman, Z. Zheng, P. Benítez, and J. C. Miñano, “Initial design with L2 Monge-Kantorovich theory for the Monge–Ampère equation method in freeform surface illumination design,” Opt. Express **22**(13), 16161–16177 (2014). [CrossRef] [PubMed]

**7. **C. R. Prins, J. H. M. ten Thije Boonkkamp, J. Van Roosmalen, W. L. I. Jzerman, and T. W. Tukker, “A Monge-Ampère-Solver for free-form reflector design,” SIAM J. Sci. Comput. **36**(3), B640–B660 (2014). [CrossRef]

**8. **K. Brix, Y. Hafizogullari, and A. Platen, “Designing illumination lenses and mirrors by the numerical solution of Monge–Ampère equations,” J. Opt. Soc. Am. A **32**(11), 2227–2236 (2015). [CrossRef]

**9. **V. I. Oliker, “Mathematical aspects of design of beam shaping surfaces in geometrical optics,” in *Trends in Nonlinear Analysis*, M. Kirkilionis, S. Kromker, R. Rannacher, and F. Tomi, eds. (Springer-Verlag, 2003), pp. 193–222. [CrossRef]

**10. **F. R. Fournier, W. J. Cassarly, and J. P. Rolland, “Fast freeform reflector generation using source-target maps,” Opt. Express **18**(5), 5295–5304 (2010). [CrossRef] [PubMed]

**11. **C. Canavesi, W. J. Cassarly, and J. P. Rolland, “Target flux estimation by calculating intersections between neighboring conic reflector patches,” Opt. Lett. **38**(23), 5012–5015 (2013). [CrossRef] [PubMed]

**12. **D. Ma, Z. Feng, and R. Liang, “Tailoring freeform illumination optics in a double-pole coordinate system,” Appl. Opt. **54**(9), 2395–2399 (2015). [CrossRef] [PubMed]

**13. **D. Michaelis, D. Schreiber, and A. Brauer, “Cartesian oval representation of freeform optics in illumination systems,” Opt. Lett. **36**(6), 918–920 (2011). [CrossRef] [PubMed]

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

**15. **V. Oliker, J. Rubinstein, and G. Wolansky, “Ray mapping and illumination control,” J. Photonics Energy **3**(1), 035599 (2013). [CrossRef]

**16. **V. Oliker, “Differential equations for design of a freeform single lens with prescribed irradiance properties,” Opt. Eng. **53**(3), 031302 (2014). [CrossRef]

**17. **V. Oliker and B. Cherkasskiy, “Controlling light with freeform optics: recent progress in computational methods for optical design of freeform lenses with prescribed irradiance properties,” Proc. SPIE **9191**, 919105 (2014). [CrossRef]

**18. **A. Bauerle, A. Bruneton, P. Loosen, and R. Wester, “Algorithm for irradiance tailoring using multiple freeform optical surfaces,” Opt. Express **20**(13), 14477–14485 (2012). [CrossRef] [PubMed]

**19. **A. Bruneton, A. Bauerle, P. Loosen, and R. Wester, “High resolution irradiance tailoring using multiple freeform surfaces,” Opt. Express **21**(9), 10563–10571 (2013). [CrossRef] [PubMed]

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

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

**22. **J.C. Miñano, P. Benítez, W. Lin, J. Infante, F. Muñoz, and A. Santamaría, “An application of the SMS method for imaging designs,” Opt. Express **17**(26), 24036–24044 (2009). [CrossRef]

**23. **A. Hofmann, J. Unterhinninghofen, H. Ries, and S. Kaiser, “Double tailoring of freeform surfaces for off-axis aplanatic systems,” Proc. SPIE **8550**, 855014 (2012). [CrossRef]

**24. **F. Duerr, P. Benítez, J.C. Miñano, Y. Meuret, and H. Thienpont, “Analytic design method for optimal imaging: coupling three ray sets using two free-form lens profiles,” Opt. Express **20**5, 5576–5585 (2012).

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

**26. **J.-D. Benamou, Y. Brenier, and K. Guittet, “The Monge-Kantorovich mass transfer and its Compuational Fluid Mechanics formulation,” Int. J. Numer. Meth. Fluids **40**(1–2), 21–30 (2002). [CrossRef]

**27. **S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent, “Optimal mass transport for registration and warping,” Int. J. Comp. Vis. **60**(3), 225–240 (2004). [CrossRef]

**28. **M. M. Sulman, J. F. Williams, and R. D. Russel, “An efficient approach for the numerical solution of Monge-Ampère equation,” Appl. Numer. Math. **61**(3), 298–307 (2011). [CrossRef]

**29. **D. Kuzmin, *A Guide to Numerical Methods for Transport Equations* (University Erlangen-Nuremberg, 2010).