A least-squares method for the design of two-reflector optical systems

The purpose of this paper is to present a method for the design of two-reflector optical systems that transfer a given energy density of the source to a desired energy density at the target. It is known that the two-reflector design problem gives rise to a Monge–Ampère (MA) equation with transport boundary condition. We solve this boundary value problem using a recently developed least-squares algorithm (Prins et al 2015 J. Sci. Comput. 37 B937–61). It is one of the few numerical algorithms capable to solve these type of problems efficiently. The least-squares algorithm can provide two solutions of the MA problem, one is concave and the other one is convex. The reflectors are validated for several numerical examples by a ray-tracer based on Monte-Carlo simulation.


Introduction
The optical design problem can be formulated as an inverse problem: determine an optical system consisting of reflectors and/or lenses that converts a given light distribution of the source into a desired target light distribution. Solution methods for inverse problems, in short inverse methods, can significantly speed up the design process, and even provide designs that could realistically never be achieved without these methods [1,2]. On the other hand, we can consider the forward problem: given a light source and an optical system, compute the resulting target light distribution. The prototypical solution method for this problem is Monte-Carlo ray tracing [3]. In this method, a large collection of rays is traced through an optical system until they hit a target receiver. From the distribution of target rays one can compute an estimate of the light output. The procedure is simple but slow, since the error decreases proportional to the reciprocal of the square root of the number of rays. Consequently, optical design by ray tracing is a slow process, all the more since the resulting light output is most likely not equal to the desired output. Therefore ray tracing has to be embedded in an iterative procedure to update the optical system. Alternatively, ray tracing can be used to validate the optical design computed by an inverse method.
In this article, we design an optical system consisting of two reflector surfaces using inverse methods. A partial differential equation for the reflector surfaces can be derived using the law of reflection and conservation of luminous flux. This partial differential equation is a Monge-Ampère (MA) type equation with a transport boundary condition [4,5]. The MA-equation has a unique convex and a unique concave solution [6, p 125].
This so-called second boundary value problem is closely related to the optimal mass transport problem: 'given a pile of sand and a hole, find a plan to transport the sand into the hole while minimizing the total transport cost' [4]. In fact, the optical design problem is equivalent to find an optimal transfer plan, i.e. a mapping m which transfers the given energy distribution of the source to the desired energy distribution at the target, that also minimizes the total transport cost. Glimm and Oliker have shown that the optical design of tworeflector systems for a light source emitting parallel light rays is equivalent to the Monge-Kantorovich mass transfer problem with a quadratic cost function [1].
There are several methods which can be employed to design these optical systems. Benamou et al [7][8][9] presented numerical methods based on the theory of viscosity solutions of MA-type equations. These are alternative methods for the optimal transport problem with quadratic costs, and can also be employed for optical design problems. Furthermore, these numerical methods are also applied to optical design problems by Feng et al [10,11]. In [10] the authors first solve the MA-equation using the optimal transport algorithm (theory of viscosity solutions) described in [8,9], then refine the solution using an iterative Fourier-transform algorithm. A drawback of the method is that it can only be implemented in the paraxial approximation. In [11] the authors formulated a MA-equation in terms of the stereographic coordinates and solved it in a least-squares sense using Newton's method. Brix et al [12,13] solved the optical design problem using a collocation method with a tensorproduct B-spline basis. Glimm and Oliker [1] showed that the illuminance control problem can be solved using an optimization approach instead of solving a MA-type equation. Further, a similar approach was used by Rubinstein and Wolansky [14] to design freeform surfaces of a lens.
We solve the MA boundary value problem using a variant of the least-squares method introduced by Prins et al [15]. The least-squares algorithm is an iterative minimization procedure. At each iteration, three steps are performed: two of these are nonlinear minimization steps, which can be performed point-wise, and the third step requires the solution of two Poisson problems. The least-squares method gives the mapping which transforms the given energy density at the source to the desired energy density at the target and the optical system is obtained via this mapping. The least-squares method has shown good performance for single surface optical systems in the far field approximation [2]. In this article, we challenge the algorithm for two-reflector optical systems for a given emittance profile with parallel light rays to achieve a desired output illuminance profile at the target, also with parallel light rays. We compute both solutions, i.e. convex and concave, of the MA-equation using the least-squares algorithm although Prins et al [15] obtained only the convex solution. We show that the algorithm works smoothly even for complex target distributions.
This article is organized as follows. In section 2 we give a generalized approach for geometrical construction of optical systems. The least-squares algorithm to solve the MA-equation is presented in section 3. In section 4 we provide the outline of the algorithm for two-reflector systems. In section 5 we calculate the solution of the MA-equation from the optimal mapping which is nothing but the surface of the first reflector, the second reflector is calculated via a relation which is established in section 2. The reflectors are validated by a ray-tracing algorithm based on Monte-Carlo simulation, which is explained in appendix. We give several test results in section 6 and we conclude with discussion and conclusions in section 7.

Formulation of the optical system
Let us consider the two-reflector system shown schematically in figure 1 3 denote the Cartesian coordinates, with z the vertical coordinate and 2 the coordinates in the plane z=0, denoted by α 1 , and let  be a bounded domain in the plane α 1 . Consider a beam B 1 of parallel light rays emitted from the source  , which is propagating in the positive z-direction. The emittance, i.e. luminous flux per unit area, of the light at the plane α 1 is given by , where f is a non-negative integrable function on the domain  (for an introduction on photometric quantities see, e.g. [2, pp 7-9]).
The incoming beam B 1 is intercepted by the first reflector  1 , defined as the graph of a function = ( ) . The beam B 1 is reflected off  1 and forms a new beam B 2 . The beam B 2 is intercepted by reflector  2 , which transforms it into the output beam B 3 , that is required to have parallel light rays again. The reflector  2 is defined as the graph of a function y y , 1 2 are the Cartesian coordinates of the target plane α 2 : z = ℓ, and obtained by a mapping    m: . The beam B 3 consists of parallel light rays propagating in the same direction as B 1 . The target at a distance ℓ>0 from the plane α 1 is denoted by  .
The goal is to design a two-reflector system such that after two reflections the outgoing rays again are parallel with a prescribed illuminance lm m 2 at the plane α 2 , where g is a non-negative integrable function on the domain  . It is assumed that both  1 and  2 are perfect reflectors and no energy is lost in the reflections. Conservation of energy gives the constraint The key tool for the design of such an optical system is to find a mapping   =  ( ) y m x : that satisfies the energy conservation constraint (1), i.e. for each   Ì , we have and after a change of variables the constraint becomes where m D is the Jacobian of the mapping    m : , which measures the expansion/contraction of a tube of rays due to the two reflections.
The mapping m can be derived by tracing a typical ray through the optical system. Let us consider a ray emitted from a position  Î x on the source and propagating in the positive z-direction and let =ŝ e z be the unit direction vector of the incident ray. Throughout this article, we use the convention that a hat denotes a unit vector. The ray strikes the first reflector  1 , reflects off in the directiont , strikes the second reflector  2 , and reflects off, again in the directionŝ . The downward unit normal vectorn 1 on  1 is given by According to the law of reflection, the directiont of the reflected ray is given by stating that the normal component of the incident ray is reversed. Note that =ˆ( ) t t x . Let us denote by ( ) x d the distance from reflector  1 to reflector  2 along the ray reflected in the direction ( ) t x . The total optical path length ( ) x L corresponding to the ray associated with a point  Î x , is given by The theorem of Malus and Dupin (the principle of equal optical path) states that the total optical path length between any two orthogonal wave-fronts is the same for all rays [16, p 130], so the total optical path length will be independent of the position vectorxand is given by On the other hand, the vertical distance between source and target planes is given by where the second term is the projection ofˆ( ) t x d on -ŝ , which is counted twice in the sum + ( ) ( ( )) x m x u u 1 2 . Subtracting (7) from (6) and using the relation = ·ŝ t t 3 , where t 3 is the third component oft , we obtain the following expression for the distance between the reflectors The image on the target of the point . This mapping can be obtained by the projection oft on the plane a 1 , i.e., where t 1 , t 2 are the first two components of the vectort . Using the expression fort in equation (5) we can rewrite (9) as follows Combining equations (8) and (10), we find the relation where β=L−ℓ is the 'reduced' optical path length [1].
Next, we will derive the mathematical expressions for the reflector surfaces. The optical distance d between the reflectors is given by Thus, from equations (6) and (12) we obtain which can be rewritten as This is a mathematical expression for the location of the reflectors. The right hand side, say ( ) represents the sum of the reflectors heights. Now the problem is to find such a pair of reflectors that satisfies relation (13) using the energy conservation equation (3) and mapping relation (11). Also, this problem is closely related to the mass transfer problem as Glimm and Oliker have shown in their article [1], with a quadratic cost function which equals the right hand side of equation (13). Thus, once we have found such a mapping m that satisfies the energy conservation relation (3) then the reflector surfaces can be obtained from equations (11) and (13), respectively. Note that relation (11) can be written as where Φ is some potential function on  . There is an important theorem by Brenier [6, p 125] stating that such an optimal mapping is the unique gradient of a convex or a concave function Φ. Substituting the mapping ( ) m x from (14) in equation (3) and by choosing the positive sign for the determinant, we obtain , 15 2 where F D 2 is the Hessian of Φ. This equation is known as the MA equation. The accompanying boundary condition states that the boundary of the source  is mapped to the boundary of the target  and is a consequence of the edge-ray principle [17]. This MA problem has a unique convex and a unique concave solution. Both solutions can be obtained using the least-squares algorithm. In this article, we will give a detailed explanation for the concave solution Φ and a brief summary for the convex solution.
Note that equation (15a) only makes sense if F > ( ( )) x g 0, a condition that is taken care of in our numerical simulations. On the other hand = ( ) x f 0 is possible, implying that the Hessian matrix F D 2 is singular, however, this does not pose a problem for our least-squares algorithm.

Numerical method
In this section, we extend a recently developed numerical method by Prins et al [15]. They introduced a leastsquares method to solve MA type equations and used the method to design single-surface optical systems. Here, we extend the least-squares method to solve boundary value problem (15a)-(15b) for two-reflector optical systems. Moreover, we compute both convex and concave solutions. We require the Jacobi matrix ofm, given by The matrix P approximates D u 2 1 (for simplicity of notation we omit the arguments). The matrix Q, approximating F D 2 , should be negative semi-definite to obtain the concave solution Φ, which in turn gives a concave reflector surface ( ) x u 1 as well, because the sum of two concave functions gives a concave function, see (14). We know that the 2×2 matrix Q is negative semi-definite if and only if Since the source emittance f and the target illuminance g are non-negative, from (17) we conclude that second inequality in (18) is already satisfied. We only have to verify the inequality 1 , which is the condition for m to be a conservative field on an open and simply-connected domain  , and thus the gradient of a function; for more details see [18, p 866].
We enforce the equality = m Q D by minimizing the following functional: The norm   · F used in this functional is the Frobenius norm, which is defined as follows. Let A B : denote the Frobenius inner product of the matrices = ( ) A a i j , and = ( ) the Frobenius norm is then defined as F . Next, we address the boundary condition (15b) by minimizing the functional where |·| denotes the Euclidian of ℓ 2 -norm for vectors. We combine the functionals J I for the interior and J B for the boundary by a weighted average:

B
The parameter α (0<α<1) controls the weight of the first functional compared to the second functional. The functionals J J , I B and J are defined on the following spaces, , 23 respectively. The minimizer gives us the mappingm which is the gradient of the potential Φ. We calculate this minimizer by repeatedly minimizing over the three spaces separately. We start with an initial guess m 0 , which will be specified shortly. Subsequently, we perform the iterations . 24 m n n n 1 1 1 We initialize our minimization procedure by constructing an initial guess m 0 which maps the source area  to a bounding box of the target area  . Let   Note that the corresponding Jacobian matrix = Q m D 0 0 of the initial condition is symmetric (in fact diagonal) negative definite, and this initial condition satisfies our requirement  ( ) Q tr 0. We solve the MA problem (15) by discretizing the source  with a standard rectangular N 1 × N 2 grid for some  Î N N , 1 2 , so the grid points Finally, we obtain the converged mapping m and calculate the reflector surfaces ( ) x u 1 and ( ) y u 2 via relation (11) and equation (13), respectively, see section 5.

Convex solution of the MA equation
Note that for the convex case, we have following condition instead of (18) on matrixQ Taking into account the above conditions we obtain the convex solution Φ of the MA problem (15) but note that the convexity of solution Φ does not guarantee the convexity of reflector surface u 1 which can be seen from relation (14). Also, we have to make the required changes in the initial condition (25) to have a symmetric positive definite matrix = Q m D 0 0 , which can be simply achieved by swapping c min to c max , and d min to d max in initial condition (25).

Minimizing procedures for b Q
, and m We start the minimizing iteration process (24) using initial guess m 0 . Here each iteration consists of the three steps (24a)-(24c), the minimization steps (24a) and (24b) are performed point-wise because no derivatives of b and Q appear in the corresponding functionals. In the third minimization step, we solve two Poisson problems for the components of the mappingm. In this article we give the detailed description of the minimization step (24b), for the details of minimization steps (24a) and (24c) see [15].

Minimizing procedure forQ
The minimizing procedure for Q differs from the minimizing procedure prescribed in [15] because here we minimize the functional J I to obtain a concave solution instead of a convex solution of the problem (15). Thus, we have the concavity condition (18) in the minimization (24b). We assume m fixed and minimize ( ) m Q J , I over the matrices under the condition (17).
Since the integrand of ( ) m Q J , I does not contain derivatives of Q, the minimization procedure can be done point-wise. Thus we need to minimize - , where D is the standard second order central difference approximation of m D . Let us define  for a threshold value e > 0. In our computations we choose ε=10 −8 . Note that the matrix D need not be symmetric and d 12 , d 21 >0 is possible. Next we define the function Also we define the matrix D S , the symmetric part of the matrix D, i.e., 2  + ( ) q q c 0. 33 11 22 Prins et al [15] solved the above problem analytically, and the authors have shown that for given d d d , , 11 22 S and f/g there exist at least one and at most four possible solutions that minimize ( ) H q q q , , S 11 22 12 and satisfy the nonlinear constraint -= q q q f g 11 22 12 2 , which are local minima. From these we have to select the ones that give rise to a negative semi-definite matrix Q, i.e. satisfy q 11 +q 22 0, and we will show that this is always possible. Finally, we compare the values of ( ) H q q q , , S 11 22 12 to find the global minimum. The possible minimizers of problem (33) are obtained using the Lagrangian function, defined as , ; det , 34 11  where λ is the Lagrange multiplier. By setting all partial derivatives of Λ to 0, we find the critical points of Λ, and this gives the following algebraic system l . 37

S
We can solve the quartic equation (37) analytically using Ferrari's method [19, p 32]. This equation has four roots, for the detailed solution see [15]. Furthermore, it can be shown that the quartic equation has at least two real roots. For the real symmetric matrix D S , we can deduce

Computation of reflector surfaces and algorithm summary
The minimization procedure steps (24a)-(24c) are repeated until ( ) m Q b J , , is not decreasing anymore. Since we are interested in the solution of the MA equation which is nothing but the first reflector surface ( ) x u 1 , we compute the reflector ( ) x u 1 from the converged mappingm, i.e. we calculate u 1 such that b +  = x m u 1 . Let us rewrite this relation as 1 , which impliesm is a conservative vector field [18, p 866] and thus there exist a potential u 1 such that  =m u 1 . However, we will most likely not be in this ideal situation, therefore we look for a function that satisfies equation (15a) in the least-squares sense, i.e., We calculate the minimizing function ( ) x u 1 using calculus of variations. The first variation of the functional (43) in a direction v is given by The minimizer is given by Letn denote the unit outward normal at the boundary  ¶ . Using the following identity derived from Gauss's theorem, we conclude from equations (44) and (45) that Applying the fundamental lemma of calculus of variation [20, p 15] we find This is a Neumann problem, and only has a solution if the compatibility condition is satisfied, which reads

Numerical results
Two solutions for the optical system can be obtained using the least-squares algorithm, one of them is concave and the other one is convex. We test the algorithm for two problems for the concave solution: in the first problem we map a square parallel beam of light with uniform emittance into a circular parallel beam with uniform illuminance, and in the second problem we transform a square uniform parallel beam into a picture on a screen. We also show the differences between the concave and convex solutions of the optical system for the second problem. The numerical results are verified by our self-build ray tracer based on the Monte-Carlo method, see appendix.

From a uniform square to a uniform circle
The first test case is the design of an optical system of two reflectors that transforms the uniform emittance of a square into a circle. The source domain is given by 4 1, 1 and the target domain by . The light source emits a parallel beam of light with uniform emittance, i.e., The target plane is at a distance = ℓ 20 from the source plane and we have fixed the reduced optical path length, more specifically b = + ( ) 5 1 2 , which keeps the incident angle π/8 for both reflectors corresponding to the central ray, i.e. β fixes the reflectors in such a way that the tangent of the reflectors at the central point have inclination of π/4. The target  is illuminated by a parallel beam of light rays with uniform illuminance, i.e., Note that the energy conservation condition (1) is satisfied. We solve the boundary value problem (15) for a uniform 200×200 grid. We found from various experiments that the number of grid points N b on the boundary  ¶ to minimize the functional J B , does not influence the convergence of the algorithm if it is chosen large enough. Since a large value of N b does not significantly increase the calculation time, we have chosen = N 1000 b . Another issue is the choice of α. Obviously, α should not be too close to either 0 or 1, since then either the minimization of J I or J B does not have enough weight in the overall iteration. Ideally, we have to choose α such that » J B I B for large iteration numbers. Unfortunately, we do not have a convergence analysis of the iteration (24), therefore, we have to determine α empirically. From previous numerical simulations, we know that convergence of (24) weakly depends on the choice of α, provided α is bounded away from 0 and 1; see e.g. [15,21]. Furthermore, a cautious conclusion from these references is that for targets with a non-smooth boundary, it is good to choose a rather small α, a result which is confirmed by the convergence histories in figure 2.
The resulting mapping after 200 iterations on a 200×200 grid is shown in figure 3(a). We tested the reflectors using the ray trace algorithm. We ran the algorithm for 10 6 uniformly distributed random points on the source, the corresponding illuminance on the target is plotted in figure 3(b). The functions u 1 (x 1 , x 2 ) and u 2 (y 1 , y 2 ) representing the reflector surfaces  1 and  2 , respectively, are shown with flat shading in figure 4 for the concave case and in figure 5 for the convex case. The reflector surfaces look almost flat for the convex solution.

From a uniform square to a picture
The second test case is the design of an optical system of two reflectors that transforms a square uniform parallel beam of light into a target distribution corresponding to a given picture. Here, we challenge our algorithm for two different pictures on the target. The emittance of the light source is again the same as defined in (51) and the parameters ℓ and β are also the same as defined in section 6.1. The desired target illuminance g(y 1 , y 2 ) is given by two gray-scale test images shown in figure 6. Note that the pictures are converted into gray-scale and contain some black regions, which results in g(y 1 , y 2 )=0 for some  Î ( ) y y , 1 2 . This would give division by 0 in the MA equation. Therefore, we have increased the illuminance to 5% of the maximum value if it is less than this threshold value.
The two test pictures pose different challenges for our optical design algorithm. The first test image is a photo of one of the contributing authors 'Jan ten Thije Boonkkamp', the original picture is shown in figure 6(a). The photo has different features, e.g. the dark and sharp edges of the glasses and hair of the beard. In the second test image, Mandrill, see figure 6(b), the challenge is to depict the whiskers of the animal. Our algorithm was able to reproduce all these details well.  Consider the 'Jan ten Thije Boonkkamp' test case. We covered the source  with a 500×500 grid, with 1000 boundary points. The convergence history of the algorithm is shown in figure 7(a) for α=0.2, which from several numerical simulations proved to be the best choice for convergence. The algorithm exhibits the same convergence for the two test images. We stopped the algorithm after 150 iterations, because J I and J B did no longer seem to decrease. The resulting mapping is shown in figure 7(b). The mapping contains the gradient of the reflector surface u 1 (x 1 , x 2 ). Another representation of the mapping can be seen as contour of grids by plotting the gradient of the reflectors, which is shown in figure 8.
The reflector surfaces are validated using our ray tracing algorithm. We ran the algorithm for 10 7 uniformly distributed random points on the source to compute the actual illumination pattern produced on the target. The target illuminance is plotted in figure 6. We observed that 99.99% of the rays hit the target, the remaining 0.01% does not arrive at the target due to numerical errors. Each of the two output images closely resembles the corresponding original version, even the finest details are reproduced, although the images are slightly blurred. To improve the resolution, we have to refine the grid covering  to compute a more accurate representation of the reflector surfaces, and moreover, trace (substantially) more rays.
The least-squares algorithm is very memory-efficient. All our numerical calculations are performed on a laptop with 32 GB of RAM. The calculation times and numerical residuals in the algorithm corresponding to the iteration counts for the two test images are shown in table 1 for a 500×500 grid. Furthermore, the calculation times and residuals corresponding to several grids and a fixed number of 150 iterations is presented in table 2. The calculation times are independent of the two test images. The parameters for the results in tables 1 and 2 are α=0.2 and N b =1000.

Discussion and conclusions
We have presented a geometrical formulation for a two-reflector optical system which is based on the principle of equal optical path length. This formulation is rigorously used to obtain a mathematical expression for the location/shape of the reflectors. This is also suitable for more general optical systems consisting of one or more reflecting/refracting surfaces.
We have extended the recently developed least-squares method to design an optical system with two reflectors described by a MA type equation with transport boundary condition. Compared to other available  methods, we know of [1,5,13,14], none of them is able to provide both convex and concave solutions. Our method is capable to solve complex problems efficiently as is evident from the previous section.
The ability to solve the optical design problem for both convex and concave solutions makes our method a valuable addition to existing methods. The algorithm is very time and memory efficient which makes it very suitable to use the algorithm for these type of problems. The least-squares method has shown very good performance to achieve the desired target, even for complex images. We have also developed a ray tracing algorithm to verify our results obtained by the least-squares method.  In future work we plan to apply this approach to more general optical systems consisting of reflectors/lenses. A partial differential equation governing these optical systems can be derived using similar principles of optics. We will work on systems for a point light source as well as for a parallel beam of light.

Appendix. The ray-tracing algorithm
Here, we describe a ray-tracing algorithm to verify our optical systems obtained by the least-squares method. We use our self-build ray tracing algorithm to have full control and to avoid additional discretization error at the target illuminance. Given a ray emitted at  Î x ij , the algorithm steps are the following.
• Find the quadrilateral on  containingx ij and the interpolated value ( ) x u ij 1 of the reflector  1 .
• Compute the direction of the reflected ray.
• Compute the intersection of the reflected ray with the reflector  2 . First of all, we generate a random point = p x ij on the source domain  using the Monte-Carlo random generator. Secondly, we calculate the interpolated value of ( ) p u 1 using bilinear interpolation. For that we have to Figure 8. Gradients of the reflectors for the 'Jan ten Thije Boonkkamp' test case.  know the quadrilateral that contains the point p. In our calculation the source domain has a regular rectangular grid. We find the nearest grid point to the point p via the minimum distance approach, then we search for the point p in the four surrounding rectangular cells and find the proper cell using MATLAB's inbuilt function inpolygon. Let us, say v 1 is the interpolated z-coordinate of the reflector  1 at the point p, then the corresponding coordinates on the reflector  1 are (p 1 , p 2 , v 1 ).

A.2. Compute the direction of the reflected rays
The direction of the reflected ray = ( ) t t t t , , 1 2 3 is calculated using the law of reflection (5), with = ( ) s 0, 0, 1 the direction of the incident rays, andn 1 the unit downward normal of the reflector  1 which is given by relation (4).
A.3. Compute the intersection of the reflected ray with the reflector  2 Once we obtained the interpolation value ( ) p p v , , 1 of the reflector  1 corresponding to the point p, and the reflected directiont , we find the intersection of the reflected ray with the second reflector  2 . The reflected ray goes through the cylinder generated by the target domain  (see figure A1) but it is hard to locate the intersection point due to the irregularity of the grid on  , see figure 8(a). The target grid depends on the illuminance pattern of the target which is not uniform in most cases. We use the fact that the ray propagating towards the reflector  2 intersects the cylinder generated by the target  at least twice. The reflected ray can be written in the vector form as where 0μ(i)1 for all i=1, 2, K, N b and N b is number of boundary grid points. The reflected ray given by equation (A.1) intersects a boundary segment of the cylinder if where the denominators in above expressions are not 0 in our calculations. We repeat this process for all line segments of the target boundary. There will be an intersection between the reflected ray and the cylinder Figure A1. Sketch of the cylinder generated by the target. generated by the target if and only if 0μ(i)1 for at least one i. If a light ray enters the cylinder then it also has to leave the cylinder so there will be at least two intersection points (may be more due to irregularity of the boundary or no intersection if the ray is going in wrong direction). If there is no intersection (for a ray almost parallel to a boundary segment) between the reflected ray and line segments of the target boundary we will go for the next ray and if there is an intersection we go to the next step of our algorithm. Let us say there are at least two intersection points then we start with the first two s-values and find the interpolated value u 2 of the reflector  2 at these points using bilinear interpolation. Then we find the position of these points by comparing the z-coordinate of these points with the cylinder and the corresponding interpolated value of the reflector  2 . The possibilities are: • If any of these points is approximately on the reflector  2 then go to the next step (find corresponding illuminance on the target).
• If both points are either above or below the reflector  2 then go for the next two intersection points.
• If one of them is above the reflector  2 and other one is below the reflector  2 then go to the bisection process explained below.
Let us say a=s(1) is above the reflector and b=s(2) below the reflector. Now, we have two points to start the bisection process, so we can calculate the desired point through the bisection method. Note that interpolation of the reflector  2 at all iteration steps is calculated using bilinear interpolation.
A.4. Find the quadrilateral on  containing = ( ) y m x ij ij and the interpolated value ( ) y u ij 2 of the reflector  2 Before calculating the interpolated value of ( ) y u 2 at the specified point through bilinear interpolation we have to know the proper quadrilateral which contains the point. The target grid is in general not regular since it depends on the desired target illuminance. We go from a regular grid on the source to an irregular grid on the target via the mapping m which is given by relation (11). To find the grid point closest to the specified point, we minimize the following norm to calculate the nearest grid point - Due to the skewness of the target grid and numerical noise it is not always possible to find the required point in the nearest four surrounding quadrilaterals, thus we extended our search algorithm to sixteen surrounding quadrilaterals. Once we find the proper quadrilateral we can find the corresponding interpolated value as described above. We can verify that the direction of the reflected ray is the same as the direction of the incident ray using the law of reflection (5).
A.5. Calculate the illumination on the target Finally, we calculate the target illuminance by calculating the intersection of the reflected rays with the plane α 2 and collect the rays in bins.