Iterative Scheme for Solving Optimal Transportation Problems Arising in Reflector Design
We consider the geometric optics problem of finding a system of two reflectors that transform a spherical wavefront into a beam of parallel rays with prescribed intensity distribution. Using techniques from optimal transportation theory, it has been shown previously that this problem is equivalent to an infinite-dimensional linear programming (LP) problem. Here we investigate techniques for constructing the two reflectors numerically by considering the finite-dimensional LP problems which arise as approximations to the infinite-dimensional problem. A straightforward discretization has the disadvantage that the number of constraints increases rapidly with the mesh size, so only very coarse meshes are practical. To address this well-known issue we propose an iterative solution scheme. In each step, an LP problem is solved, where information from the previous iteration step is used to reduce the number of necessary constraints. As an illustration, we apply our proposed scheme to solve a problem with synthetic data, demonstrating that the method allows for much finer meshes than a simple discretization. We also give evidence that the scheme converges. There exists a growing literature for the application of optimal transportation theory to other beam shaping problems, and our proposed scheme is easy to adapt for these problems as well.
The following beam shaping problem from geometric optics was described in ; see Figure 1. Suppose a point source emits a spherical wavefront with a given intensity distribution. The problem we are concerned with consists of transforming this input beam into an output beam of parallel light rays with a prescribed intensity distribution. This transformation is to be achieved with a system of two reflectors. The problem has some practical importance in engineering; see further literature cited in .
The first rigorous mathematical solution to the problem was provided in , using an approach based on the theory of optimal transportation [3–6]. See also the references [7–11] which deal with other beam shaping problems using related techniques.
The result of  is summarized in Section 2 below, with Theorem A.5 stating the main result. The central feature is that the original reflector design problem is reformulated as an infinite-dimensional constrained optimization problem, namely the problem of minimizing a certain linear functional on a function space. It is the dual problem for the problem of finding a map from the input aperture to the output aperture which minimizes a certain cost functional (see Theorem 3 in Section 2 for the exact statement).
This reformulation of the problem not only is of theoretical value for questions of existence and uniqueness of solutions, but it also translates into a practical method for finding the solution. In fact, the discretization of the infinite-dimensional constrained optimization problem is a standard linear programming problem and can be solved numerically. This was already observed for similar beam shaping problems by Glimm and Oliker [7, 8] and independently by Wang .
There has been increased interest in numerical methods for optimal transportation problems in the last 10 years. Most work has concentrated on the Monge-Ampère equation, which arises in the special case of optimal transport in with quadratic cost (also called optimal transport) [12–20], although some authors have treated costs proportional to the distance ( optimal transport) . These methods are based on fluid mechanical approaches  or various finite difference approaches . They are generally faster and allow for much larger mesh sizes than methods based on a discretization of the linear programming problem, but since they use the special structure of the Monge-Ampère equation, they are not directly applicable to more general cost functions or more general manifolds.
A linear programming approach similar to the one we investigate here for more general situations for optimal approach was investigated by Rüschendorf and Uckelmann in 2000  and more recently in two papers by Canavesi et al. [23, 24] who also propose a new algorithm inspired by, but distinct from, a linear programming approach for the single reflector problem treated in [8, 9].
As noted in  and other works, an immediate obstacle for the linear programming approach is that the number of constraints becomes very large even for relatively coarse meshes—if there are sample points in the input aperture (denoted by in Figure 1) and points in the output aperture ( in Figure 1), then the number of constraints is . For instance, Rüschendorf and Uckelmann note in  that the LP solver they used, SOPLEX, was designed to handle up to 2 million variables. This corresponds to mesh sizes of approximately 1400 points on the input and output aperture in our problem if one employs a straightforward discretization scheme. We found that, in our implementation, the solver MOSEK (Version 7) could handle about 3 million constraints, corresponding to approximately 1750 sample points on the domains and (see Table 1 for details). Rüschendorf and Uckelmann noted in  that, for better results, one needs carefully designed programs. This is what we supply here for our problem: We devised a more elaborate iterative method, where discretized systems are solved on finer and finer grids, in each step, using information from the previous solution to choose only a subset of all possible constraints. This drastically slows the growth of the number of constraints. Details are described in Section 3. With this step-wise mesh refinement scheme, as an illustration, we were able to obtain solutions on meshes with about 10,300 points on each aperture using MATLAB and MOSEK 7, a commercially available LP solver. The number of constraints needed for this was only a fraction of the size of the “full” system. For instance, the run displayed in Table 2 used only 0.33% of all possible constraints.
We note that our proposed algorithm does not make any a priori symmetry assumptions on the form the reflectors. It can also be adapted for the numerical solution of other beam shaping problems for which a formulation using optimal transportation theory has been found, for example those in [7, 8, 25].
This article is organized as follows: in Section 2, we recall the reformulation of the reflector construction problem as a linear programming problem as given in  and fix some notation. We then describe a basic discretization scheme for the numerical solution and propose an improved “mesh refinement” scheme in Section 3. The following Section 4 is devoted to numerical illustrations and tests of this scheme. We first derive an explicit analytical solution to be used as synthetic data for our numerical work. Then we compute the solution numerically and analyze the error of approximation. We conclude with a brief summary of the advantages and drawbacks of the method and propose future work.
2. Formulation as Linear Programming Problem
We now briefly review the notation for the problem posed in , as well as the result of , which reformulates the reflector design problem as a linear programming problem. More details are contained in Appendix A.
Consider the configuration shown in Figure 1. A point source at the origin generates a spherical wavefront over a given input aperture contained in the unit sphere . We use the bar in the notation for to stress that is a closed set. The input beam has a given intensity distribution. By means of two reflectors, this wavefront is to be transformed into a beam of parallel rays propagating in the direction of the negative -axis. This output beam is required to have a prescribed intensity distribution. The cross section of the output beam in a plane perpendicular to the direction of propagation is called the output aperture and is denoted by . Again, we use the notation to stress that is closed. Certain regularity conditions apply to and ; see  for the technical details.
Denote points in space by pairs , where is the position vector in a plane perpendicular to the direction of propagation and is the coordinate in the (negative) direction of propagation. See again Figure 1 for our convention on the direction of the -axis. Points on the unit sphere will typically be denoted by ; their components are also written as with .
We fix the output aperture in the plane . We will seek to represent the first reflector as the graph of its polar radius (for and the second reflector as the graph of a function (for ). See Figure 1. That is
The geometrical optics approximation is assumed. It follows from general principles of geometric optics that all rays will have equal length from to the plane ; this length is called the optical path length and will be denoted by . We define the reduced optical path length as .
Oliker  showed that local energy conservation translates into a complicated partial differential equation of Monge-Ampère type for . As noted in , the resulting equation is quite involved, and a rigorous analysis of this equation seems very difficult. See equation (59) in .
To amend this, the problem was reformulated in  as an infinite dimensional linear programming problem, which makes a complete analysis possible, both concerning theoretical results on existence and uniqueness, and gives a method for practical computations.
For this, the following function , called the cost function in analogy with the theory of optimal transportation, plays an important role:
In further preparation, the following two transformations are needed.
Definition 1. Let be a continuous function defined on . Then define the function
Definition 2. Let be a continuous function defined on with . Then define the function
The following characterization of the solution of the reflector problem described above was given in .
Theorem 3. If the pair is a solution of the above reflector design problem, then the transformed pair minimizes the functional on the function space for all.
One can solve the reflector construction problem by solving the infinite-dimensional LP problem in the above theorem and then applying the inverse transformations of (3) and (4) to recover and . In the remainder of this paper, we will concentrate on numerical solutions for this problem.
It is also possible to characterize the ray-tracing map analytically.
Definition 4. Let be a solution of the reflector problem. Define its ray-tracing map as a set-valued map via
In , it is shown that is in fact single-valued for almost all . Therefore, may be regarded as a transformation . Indeed, one can show that a ray emitted from the origin in the direction will be reflected to a ray in the negative -direction labeled by ; see .
More details and rigorous statements are found in Appendix A.
3. Numerical Schemes
In the following, we describe two schemes for solving the minimization problem in Theorem 3 numerically. First, we list the given data on which the solution depends.(i)Domains , (The rigorous result in  required the additional technical assumption . For practical purposes, these assumptions can often be dropped.)(ii)Nonnegative integrable functions , , and , , with .(iii)A reduced optical length . (iv)The value for some fixed , where is the radial function of the first reflector. (Note that without this constraint, the reflectors are not uniquely determined. This can be seen in the expression for the functional in Theorem 3. Indeed, any transformation of the form , for constant will leave the objective functional unchanged.)
As described in the introduction, this scheme is memory-intensive because every pair of points from and gives rise to a constraint. To address this issue, we propose a more sophisticated scheme, Scheme 6, which is an iterative scheme using finer and finer meshes, where, in each iteration, information from the previous solution is utilized to reduce the number of constraints.
3.1. “Simple Discretization” Scheme
The characterization of solutions to the reflector problem by means of the minimization of a linear functional in Theorem 3 immediately suggests a solution algorithm by means of a discretization. Essentially, the same algorithm was used for different problems in [22–24].
Scheme 5 (simple discretization). (i) Create a mesh in the input domain by choosing sets , where the interiors of any and are disjoint (for ) and . For instance, may be a triangulation of .
(ii) Similarly, create a mesh in the output domain by choosing sets , where the interiors of any and are disjoint (for ) and
(iii) Choose sample points (for ) and (for ). Here, we may assume that is the same point as given in the data in (iv) above.
(iv) Find the solution of the following LP problem: (v) Find the numbers and such that and . This is straightforward by taking the inverse of the transformations given in Definition 1 and Scheme 5.
Then is an approximation for the true value of the radial function of the first reflector, evaluated at the sample point for . Similarly, is an approximation for the true value of the function describing the second reflector, evaluated at the sample point for . We are using the superscript “(0)” here to distinguish this solution from additional iterative approximations we will obtain with the iterative scheme described below.
Note that (7) is a discretization of the sum of integrals in the functional . We have chosen the simplest discretization here. Depending on the geometry of the problem, other discretization schemes  may yield better results. In practice, we have chosen partitions and of approximately uniform size and diameter, but depending on the form of the intensity distributions, other partitions may be more appropriate.
It is important to point out that this discretization also yields a discretized version of the ray-tracing map (see Definition 4). Namely, for each index , there is at least one corresponding index , where the constraint corresponding to the pair is active; that is, where we have This means that the point is approximately the image of the point under the ray-tracing map .
3.2. “Iterative Refinement” Scheme
As noted in the introduction, one of the drawbacks of Scheme 5 is that the constraint set in the linear programming problem in the penultimate step becomes large very fast with finer meshes, and the corresponding problem becomes too large to handle for standard LP solvers. We were not able to solve problems with more than about 1750 mesh points on each aperture on a standard PC with 4 GB RAM ().
We addressed this problem by developing an iterative scheme. First, the problem is solved for a mesh with relatively few sample points (say ). Then, a finer mesh is chosen, and the previous solution is used to reduce the number of constraints.
Specifically, we have the following scheme, depending on a number , which we call the “quasi-activity threshold,” and an integer , the nearest neighbor search depth, to be explained in detail below.
Scheme 6 (iterative refinement). Use Scheme 5 to find an initial solution for sample points on and sample points on , respectively.
Definition 7. Call a pair with and “quasi-active” if Note that the pair is “active” in the standard sense of linear programming theory if the above term is equal to zero.(i)Create a second, finer mesh on by choosing sets , where the interiors of any and are disjoint (for ) and Here, is chosen larger than , meaning that we have a finer mesh. (ii) Similarly, create a second mesh on by choosing sets , where the interiors of any and are disjoint (for ) and (iii) Choose sample points (for ) and (for ). Here, we may assume that is the same point as given in the data in (iv) above.(iv) For each pair with , and , find the nearest neighbors of from the previous mesh on and the nearest neighbors of from the previous mesh on .
Definition 8. Say that the pair is potentially active if at least one of the pairs consisting of a nearest neighbor of and a nearest neighbor of is quasi-active in the above sense. (i) Find the solution of the following LP problem: (ii) Find the numbers and such that and
The idea of Scheme 6 is to solve the discretized LP on a coarse mesh first and then use this information to reduce the number of constraints needed for the LP on a finer mesh. A key difference between Schemes 5 and 6 is that in the discretized LP problem, in Scheme 6, not all pairs of sample points from and , represented by pairs of indices , are included in the list of constraints. In fact, we would in principle only need to include those where the constraint is active, that is, where the corresponding inequality holds with equality. Of course, this information is not available a priori. Instead, we use the following heuristic. A constraint should only be included in the LP problem if the corresponding point on the second reflector is close to the image of the corresponding point on the first reflector under the ray-tracing map. In terms of the LP problem formulation, the corresponding pair of points is located close to an active pair from the previous iteration. We call such a constraint “potentially active.” Our heuristic is that a constraint should be considered potentially active if there is a pair of its respective nearest neighbors which corresponds to a quasi-active constraint form the previous iteration. Here, “quasi-active” constraints include active constraints (where the constraint is attained, i.e., the term is zero) and those which are “almost active” in the sense that the term is positive, but small. “Smallness” is encoded in the quasi-activity threshold . Clearly, increasing means more constraints are included, but also potentially the solution of the LP represents a better approximation of the exact reflector pair.
In our numerical tests, we found it advantageous to scale the functions and so that the approximations for the integrals and yield exactly the same result.
Note that (7) is a discretization of the sum of integrals in the functional . We note that we have chosen again the simplest discretization here for the sum of integrals . As in Scheme 6, depending on the geometry of the problem, other discretization schemes  may yield better results.
3.3. Choice of Thresholds for Iterative Refinement Scheme 6
For the iteration of Scheme 6 described above, we have to choose a corresponding sequence of mesh sizes . We also need to choose a corresponding sequence of thresholds .
This last question is of a certain practical importance. Indeed, if the thresholds are “high” (e.g. if they decreases slowly), then “many” constraints are included in each iteration, meaning that the size of the LP problems increases quickly. If the sequence is “low”, then “few” constraints are included in each iteration. This is of course in principle desirable, but if the thresholds are too low, this could cause that some index (with ) may not even be included in any of the constraints with any of the indices (with or vice versa. This would cause the LP problem to be unbounded, and hence, there would be no solution.
We can make a very rough estimate for a good choice of the sequence assuming for simplicity that . Let us assume that each sample point in the input aperture is paired via the ray-tracing map with a unique point and vice versa. Thus, there are such pairs. Let and denote the solutions of the reflector problem. Then the set of all points in -space close to a given pair that satisfy is approximately an ellipsoid for small . This can be seen by using the Taylor expansion around the pair in -space. The volume of this ellipsoid is proportional to . So the proportion of pairs of points that satisfy the above constraint among all constraints is roughly proportional to . Since each pair that satisfies this inequality corresponds to a constraint included in the LP of the iteration, the total number of constraints should be roughly proportional to . So to keep the number of constraints growing linearly throughout the iterations, based on these heuristics, one may choose where is a constant. These heuristics are very rough. In practice, we have found that when using the formula , there is a range of values of which all produce similar results. See Section 4.2.
4. Numerical Tests
To test the validity of the numerical scheme described above, we used it on a case where the solution is known in analytic form. In the next subsection, we first describe this special analytic solution and then discuss our results.
4.1. A Data Set
We have constructed a particular data set using an explicit analytical solution described in more detail in Appendix B. The data set of this section is obtained from the more general solution of Appendix B by choosing and , , and . We consider the input aperture which is a spherical cap centered at the point . The data do not represent a physically possible set of reflectors since there would be blockage. However, the example is useful as an illustration, and the numerical scheme itself is completely independent of the shape of the apertures or any a priori symmetries of the problem. Using the results from Appendix B, we can now find the output aperture in a straightforward manner as follows: We choose a constant output intensity Using the relation , this gives the input intensity The two reflectors are given as follows: The corresponding reduced optical path length is The two reflectors are shown in Figure 2. is part of a spheroid, and is part of a paraboloid; see Appendix B for details.
Using the data set from Section 4.1, we conducted a series of numerical tests to establish the validity of the numerical Scheme 6. Since there are explicit analytical expressions for the two solutions and describing the two reflectors as given in (18) and (19), respectively, we were able to compute the error of approximation directly. See also Figure 2 for plots of these reflectors.
The Schemes 5 and 6 were implemented in MATLAB with the package distmesh  for mesh generation in each iteration step and the commercial linear programming solver MOSEK (version 7)  for solving the resulting LP problems. Computations were performed on an Intel Core i7-2620M/2.70 GHz with 4 GB RAM, running Windows 7 (64 bit).
The mesh generation algorithm is based on a relaxation scheme of forces in a truss structure . For this, the input is the desired edge length , not the number of mesh points (for the domain ) or (for the domain ), respectively. However, it can be seen from geometric considerations that the relation between the average edge length and the number of mesh points on obeys In the tables below, we show both the number of mesh points and the edge length , as we believe both are informative.
Our primary aim is to demonstrate the practicability of the proposed Scheme 6. For comparison, we first tested the performance of the “Simple Discretization” Scheme 5. The results are summarized in Table 1. The maximum mesh sizes attainable are about 1750 points, corresponding to approximately 3 million constraints. We found that this was the maximum size for the LP solver.
To test Scheme 6, we chose a sequence of desired edge lengths by reducing by in each step; that is, . Then, based on the heuristics given above and in Section 3.3, we determined the corresponding constraint thresholds with the following formula: where and are constants. The heuristic arguments in Section 3.3 yield , but we also tested other values. As indicated in Section 3.3, “large” values for mean that “many” constraints are included in the LP in each step, which potentially means good approximations, but also high memory usage. In contrast, “small” values for improve memory usage but may mean that the LP problem may become unbounded. This is the case if constraints that would be active in the LP problem comprising of all constraints are left out in an iterated LP problem. In our tests, the only time this happened was when we chose , that is, if the quasi-activity threshold is set to zero and only active constraints are considered “quasi-active.” We conclude that our results are very robust with respect to and , meaning for a broad range of values, we reached the terminal iteration with ; see Table 2. This is because the memory limiting issue in our tests was not the LP solver but the maximum size of arrays in MATLAB when computing the cost matrix; see paragraph below. We also conducted some informal tests of the algorithm dependence on the output and input intensities and , specifically to know whether the algorithm is stable under small perturbations of and . For this, we added small perturbation to and to , then renormalized so that the two integrals were the same. We found that the algorithm was indeed stable in the sense that the same maximum number of points was reached for different values of .
Results of a typical run are summarized in Table 2. Our results show that the proposed scheme is easy to implement and makes it possible to use much finer meshes than a simple discretization scheme. The largest number of points were about 10,000 on each aperture. It is notable that what prevented our tests to include more plot points was not the LP solver but memory problems related to the maximum size of matrices in MATLAB during the set-up phase of the LP problem when computing the matrix of cost coefficients , which is of size . For the configuration we used, the maximum size of arrays in MATLAB is about 280 million variables, see .
These results give indication that the scheme converges and that the maximum error for the first reflector decreases proportional to , or equivalently proportional to , where , data from Table 2. The same formulation for scheme 5 yields an exponent of , unsurprisingly a slightly more rapid decrease, data from Table 1. The geometric distribution of errors is indicated in Figure 3.
To investigate the role of the nearest neighbor search depth , we performed several tests with different values. The results are summarized in Table 3. The value of has very little influence on the error of approximation. This makes intuitive sense as the crucial issue in each iteration is that those constraints which are active in the “full” LP problem (i.e., the one that includes all possible constraints) are included. As long as all active constraints are included, one expects to get the same solution. Note, however, that the number of constraints changes quite drastically with , ranging from 96,256 ( of all possible constraints) for to 48,2058 ( of all possible constraints) for in this example. The time to solve the LP problems, thus, also varies from 2.82 seconds for to seconds for , although those gains are dwarfed by the overall run times of the algorithm which include assembling the LP problem.
5. Conclusion and Future Work
We investigated two numerical schemes for solving an infinite-dimensional optimal transportation problem arising in reflector design, a straightforward discretization and a novel iterative scheme, which uses knowledge of the previous solution in each step to reduce the number of constraints. The latter scheme is easily adapted to similar transportation problems arising in beam shaping problems, for example see [7–9].
As a proof of concept, we implemented both schemes using MATLAB and the commercial LP solver MOSEK and showed that the new scheme is easy to implement and makes it possible to solve the problem on much finer meshes. In particular, we demonstrated that the number of constraints are greatly reduced by the new scheme, to as little as 0.09% of the full number of constraints. Our aim was to show the validity of the scheme, and we contend that performance could be further improved, for instance, by using a compiled computer language or a more specialized LP solver for transportation problems.
One drawback of the scheme is that it is quite slow. See Table 2. The most time-consuming part of the algorithm was not the solution of the LP problem but assembling the problem, that is, deciding which constraints to be included. We believe that this step, in particular, could be improved by using a compiled computer language instead of MATLAB.
There are a number of possible directions for future research. We did not rigorously prove that the scheme converges, although we strongly expect that it does. In fact, the decreasing error shown in Table 2 gives evidence for this. It would be valuable to have a rigorous proof of convergence. We also expect that PDE-based schemes for solving the transportation cost with quadratic costs [13, 17–19] could be adapted to the problem at hand, leading to faster algorithms.
A. Details on Theorem 3
In , the following notion of a reflector pair and its associated ray-tracing map is used.
Definition A.2. Let be a reflector pair. Define its reflector map, or ray-tracing map, as a set-valued map via In , it is shown that is, in fact, single-valued for almost all . Therefore, may be regarded as a transformation .
If is a “reflector pair” in the above sense, the following can be shown , justifying the choice of nomenclature. If physical copies of the corresponding surfaces (1) are made from a reflective material, then a ray emitted from the origin in the direction will be reflected to a ray in the negative direction labeled by .
With this, the reflector problem can be formulated rigorously.
Problem A.3 (reflector problem). For given input and output intensities , and , respectively, satisfying global energy conservation , find a pair that satisfies the following conditions: (i) is a reflector pair in the sense of Definition A.1; (ii)the ray-tracing map satisfies for any Borel set .
Here, denotes the standard area element on the sphere . The second condition is local energy conservation.
As we indicate below, the reflector problem can be reformulated in the following form.
Problem A.4 (minimize the functional). on the set.
Problem A.4 is the dual problem to the problem of finding a measure-preserving map from to which minimizes the transportation functional . Problems A.3 and A.4 are equivalent, as expressed in the following theorem.
Thus, solving the reflector construction problem in Problem A.3 is equivalent to solving the infinite-dimensional LP problem in Problem A.4. Indeed, it can be shown that a solution as in Theorem A.5 exists. See Corollary 6.5 in .
B. An Explicit Analytic Solution to the Reflector Problem
In order to obtain a data set where an explicit analytic solution is known, we consider a special given configuration of reflectors and then solve the “forward” problem of determining the output intensity this system produces for a given input intensity. This general family of solutions we obtained was used to generate the explicit synthetic data set in Section 4.1.
B.1. Construction of a Pair of Reflectors ,
Consider the following set of two reflectors and , sketched in Figure 4. Let be a given point in , and let and be two positive numbers with . Now, let the first reflector be the boundary of a prolate spheroid (ellipsoid of revolution) with foci at the origin and at and with major diameter . Thus, for each point on , the sum of the distances from each of the two foci equals . Let the second reflector be the boundary of a paraboloid whose main axis is the negative axis and whose focus is at with focal parameter . Thus, for each point on , the sum of the distance to the focus and the shifted -component equals .
Note that by definition of the two reflectors, any cone of light rays emitted from the origin O will be transformed into a beam of parallel rays traveling in the direction of the negative -axis. Indeed, a ray emitted from O will be reflected off toward the focus . This ray will then be reflected in the direction of the negative -axis by . See again Figure 4.
It is not hard to find explicit expressions for the two reflectors. If we write and , then we have the following expressions: where we used the shifted coordinates .
B.2. Ray-Tracing Map
One can now determine the ray-tracing map corresponding to the reflector pair , . See again Figure 4.
Let be a given direction. Thus, with and . To determine the image , consider the ray emitted from the origin in the direction . The ray will encounter the first reflector at the point and will then be reflected towards the second focus . Thus, the reflected ray can be parameterized by , where is a parameter. With (B.2), this yields that the point where the reflected ray hits the second reflector is , where The projection of this point to the plane is the value of the ray-tracing map. Thus, we have the following explicit formula for the ray-tracing map:
B.3. Solution of the Forward Problem
We can now solve the forward problem. Given the reflector pair as defined above, an input aperture and an intensity distribution , , find the output aperture and the output intensity , .
The output aperture is simply the image of under the ray-tracing map : . To find the induced intensity , use the defining property for all Borel sets . This is an energy balance equation.
The above integral equation allows us to find an explicit expression for given or vice versa. We consider, for simplicity, the case that is contained in the left hemisphere . We can then use coordinates: In these coordinates, the standard measure on is given by , where . Using this, (B.5) is then equivalent to the equation Here, denotes the Jacobian of the map . With the help of a computer algebra system like Mathematica, this can be evaluated explicitly as where again .
Tilmann Glimm gratefully acknowledges the contribution of V. Oliker for the unpublished collaborative work on an algorithm similar to the one described in this paper for a different reflector design problem. This work was done during 2002 to 2004 at Emory University.
T. Graf, On the near-field reector problem and optimal transport [Ph.D. thesis], Emory University, 2010.
G. L. Delzanno, L. Chacón, J. M. Finn, Y. Chung, and G. Lapenta, “An optimal robust equidistribution method for two-dimensional grid adaptation based on Monge-Kantorovich optimization,” Journal of Computational Physics, vol. 227, no. 23, pp. 9841–9864, 2008.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
G. Awanou, “Spline element method for the Monge-Ampere equation,” December 2010.View at: Google Scholar
C. Canavesi, W. J. Cassarly, and J. P. Rolland, “Implementation of the linear programming algorithm for freeform reflector design,” p. 84850E-84850E-6, 2012.View at: Google Scholar
I. The MathWorks. Support pages for matlab, 2012, http://www.mathworks.com/support/solutions/en/data/1-IHYHFZ/index.html.