Inverse ray mapping in phase space for two-dimensional reflective optical systems

A new method to compute the target photometric variables of non-imaging optical systems is presented. The method is based on the phase space representation of each surface that forms the optical system. All surfaces can be modeled as detectors of the incident light and emitters of the reflected light. Moreover, we assume that the source can only emit light and the target can only receive light. Therefore, one phase space is taken into account for the source and one for the target. For the other surfaces both the source and target phase spaces are considered. The output intensity is computed from the rays that leave the source and hit the target. We implement the method for two-dimensional optical systems, and we compare the new method with Monte Carlo (MC) ray tracing. This paper is a proof of principle. Therefore, we present the results for systems formed by straight lines which are all located in the same medium. Numerical results show that the intensity found with the ray mapping method equals the exact intensity. Accuracy and speed advantages of several orders are observed with the new method.


Introduction
The goal in non-imaging optics is to compute the light distribution at the target of the system. To this purpose, the Monte Carlo (MC) ray tracing procedure is often used, [4]. Rays are randomly traced from the source to the target and each time that a ray hits an optical surface the coordinates of the intersection point with that surface and the new direction are calculated. The output variables are computed dividing the target into intervals, the so-called bins, and counting the rays that arrive at each bin. To obtain the desired accuracy, millions of rays are required, therefore the method is extremely computationally expensive and it converges as the inverse of the square root of the number of rays traced (see [9]).
Recently we introduced a method to speed up MC ray tracing using the phase space (PS) approach (see [2]). In this method we consider the PS of the source and the PS of the target of the system. Given an optical surface in three dimensions, in PS every ray is described by two position coordinates and two direction coordinates. The position coordinates are given by two of the coordinates of the intersection point of the ray with the surface, while the direction coordinates are the momentum coordinates of the (tangent line to the) ray projected on the optical surface (see [6,7,10]). The output intensity is computed through an integration of the luminance. The target PS consists of regions where the luminance is positive, outside these regions the luminance is equal to zero. Assuming the luminance to be a positive constant whenever it is different from 0, only the PS coordinates of the rays corresponding to the boundaries of the regions with positive luminance need to be computed. Ray tracing on PS was tested for two-dimensional systems. Compared to MC ray tracing, far less rays need to be traced to obtain the same accuracy and a convergence proportional to the inverse of the number of rays traced is obtained (details are explained in [2]).
In this paper a new method to compute the photometric variables at the target of a non-imaging optical system is presented. We restrict ourselves to two-dimensional systems therefore, from now on, we refer to the optical surfaces as lines. We assume that the systems are formed by reflective lines, therefore refraction is not taken into account. The method employs not only the source and the target PS, as PS ray tracing does, but also the PS of all the other lines that constitute the optical system. Furthermore, instead of starting from the source, the new method starts tracing back rays from target PS.
Every line of the system (except for the source S and the target T) constitutes the target for incident rays and the source for reflected rays. Therefore, two different phase spaces are considered for the reflectors and one PS for S and T. All the PS considered are divided into regions, the boundaries of which can be determined exactly for systems formed by straight lines. We make the assumption of a Lambertian source; hence, the luminance is a positive constant when different from 0. As a consequence, the output intensity along a given direction is given by the total width of all the patches with positive luminance, measured along that direction.
Compared to MC ray tracing, the new method takes advantage of the fact that only the rays located exactly on the jump discontinuities of the luminance, from zero to a positive value, are needed. In order to determine the coordinates of these rays, an inverse map from the target to the source PS is constructed as a concatenation of several maps. Employing this inverse map we are able to detect the rays that will hit the target and that are located on the boundaries of regions with positive luminance in target PS. These rays are then traced from S to T to compute their coordinates at the target. The ray mapping method has two main advantages. First, the accuracy of the intensity is improved, in particular, for systems formed by straight lines the intensity is found analytically. Second, even less rays than for PS ray tracing are needed. Compared to PS and MC ray tracing the new method improves both the accuracy and the computational time. The only disadvantage of the new method compared to existing procedures is that it is more complicated to implement. This paper is a proof of principle and investigates the method for two different optical systems in 2D: the two-faceted cup and the multi-faceted cup. This work it is organized as follows. In Sect. 2, the method is implemented for the two-faceted cup. The numerical results of this system are presented in Sect. 3. Section 4 describes a generalization of the method to the more complex multi-faceted cup, which is a system formed by many straight lines as reflectors. The results for the 20-faceted cup are shown in Sect. 5. Conclusions are given in the last section.

Description of the phase space method for the two-faceted cup
In this section a description of the method is presented. First, it is explained for a simple optical system: the two-faceted cup, where only the law of reflection plays a role. In the following we describe the geometry of this system, then we introduce the notion of phase space for all the lines that constitute the system. Finally, we explain how to calculate the output intensity.

Rays propagation through the two-faceted cup
Given a Cartesian coordinate system (x, z), an optical system symmetric with respect to the z-axis is defined. The simplest optical system we can image is the so-called two-faceted cup, the profile of which is depicted in Fig. 1.
The light source S = [-a, a] (line 1) and the target T = [-b, b] (line 4) are two segments normal to the z-axis, where a = 2 and b = 17. The left and right reflectors (line 2 and 3) are oblique segments that connect the source and the target. All the optical lines i with i ∈ {1, . . . , 4} are located in air, therefore the refractive index n i = 1 for every i. The optical axis coincides with z-axis. From now on, the coordinates (x i , z i ) i=1,...,4 denote the intersection of a ray with line i and, s i = (-sin t i , cos t i ) is the direction vector of the ray that leaves i, with t i the angle that the ray forms with respect to the optical axis measured counterclockwise. As we consider only forward rays, the angles t i ∈ (-π/2, π/2). Therefore, a ray segment between (x i , z i ) and (x j , z j ) with j = i is parameterized in real space by: where s denotes the arc-length and s max is the maximum value that it can assume. A Lambertian optical source is considered; hence, the intensity at the source S = [-a, a] emitted in the direction t 1 is given by: where L is the luminance, a = 2 is the half width of the source S, and t 1 is the angle that the ray forms with respect to the z-axis, measured counterclockwise. A ray that leaves the source S (line 1) can hit the reflectors (lines 2 and 3) many times before reaching the target T (line 4). When a ray travels from a line i to another line j its new position is given by the intersection point between the ray and line j and the new direction is computed according to the law of reflection: where (t 1 , ν) is the inner product between t 1 and ν, ν is the unit normal to line j that the ray hits and t 1 and t 2 are unit vectors describing the directions of the incident and the reflected ray, respectively (see [1], Chap. 12, pp. 403-409). We assume that the source S cannot receive light and the target T does not emit light. The reflectors 2 and 3 are designed such that they emit and receive light, hence they constitute the target for the incident rays and the source for the reflected rays. The output intensity can be computed constructing a map from the source to the target of the system. This map can be written as the concatenation of many which can be classified as two different kinds of maps, i.e. the map that connects the source and the target PS of two different lines and the map that connects the target and the source PS of the same line. Employing the inverses of these maps we are able to detect the parts on target PS illuminated by the source. In the next paragraph the method is presented for the two-faceted cup depicted in Fig. 1.

Phase spaces of the optical system
Phase spaces of every line are two-dimensional spaces. The position coordinate in the phase space (PS) of line i is the x-coordinate of the intersection point between the ray and the line i. The direction coordinate is the sine of the angle that the ray forms with respect to the normal of line i multiplied by the index of refraction of the medium in which the ray is located. Let's now introduce some notation before explaining the details of the method. We indicate the PS with S = Q × P, where Q is the set of the position coordinates q and P is the set of the direction coordinates p = n sin τ with n the index of refraction of the medium in which the line is located and τ the angle between the ray segment inside the cup and the normal ν of the line which we choose directed inwards of the optical system. The angle τ is measured conterclockwise and τ ∈ [-π/2, π/2]. In this paper we analyze only systems located in air (n = 1), therefore, from now on, we do not write the index n anymore. The source and the target PS of a line i are indicated with S i and T i , respectively. The coordinates of every ray that reaches the line i ∈ {2, 3, 4} are indicated with (q t,i , p t,i ) on T i . After reflection, the ray leaves line i ∈ {1, 2, 3} at the same position and with a new direction, the new coordinates are indicated with (q s,i , p s,i ) on S i . Note that q s,i = q t,i while p s,i is obtained applying the reflection law to the direction coordinate p t,i of the incident ray. The phase spaces S i and T i of each line i are partitioned into different regions, (S i,j ) j=2, 3,4 and (T i,k ) k=1,2,3 , respectively, where j = i is the index of the line that is illuminated by i and k = i is the index of the line that illuminates i. Hence, we indicate with S i,j ⊂ S i the part of S i corresponding to rays that illuminate line j and with T i,k ⊂ T i the part of T i corresponding to rays originating from the line k. For the two-faceted cup, six different phase spaces need to be considered which are given by the following expressions:

4)
We note that, as the source cannot receive light and the target cannot emit light, the regions (S i,1 ) i=2, 3 and (T i,4 ) i=2, 3 are not considered. The boundaries ∂S i,j are mapped into the boundaries ∂T j,i for every i ∈ {1, 2, 3} and j ∈ {2, 3, 4} with j = i (edge-ray principle [8]). For the two-faceted cup and for all systems that are formed by straight lines, they are determined analytically, see next section.

Computation of the boundaries of the patches in PS with positive luminance
Given two lines i and k with i = k, we show how to compute the boundaries of the region formed by the rays that leave line i and hit line k. We do that both for S i and for T k . We indicate with (x i, , z i, ) and (x i,r , z i,r ) the coordinates of the points located at the left and the right end point of line i, respectively. Similarly, (x k, , z k, ) and (x k,r , z k,r ) are the coordinates of the points located at the left and the right end point of line k, respectively. Given two lines i and k with i = k, ∂S i,k and ∂T k,i are formed by four different curves, two of them are given by all the rays that leave the end points of line i and trace out line k and, the others two are given by the rays that trace out line i and hit the end points of line k. The boundaries ∂S i,k and ∂T k,i are given by In the following we explain in more details the case of i = 1 and k = 4; see Fig. 2. The boundaries ∂S 1,4 and ∂T 4,1 are given in Figs. 3 and 4, respectively. ∂S 1 1,4 and ∂T 1 4,1 are obtained tracing out line 4 from q 4,min = -b to q 4,max = b by rays leaving q 1,min = -a with varying p 1 , these rays are shown in Fig. 2(a), and the boundary segments ∂S 1 1,4 and ∂T 1 4,1 are the orange line segments labeled with c. ∂S 2 1,4 and ∂T 2 4,1 are given tracing out line 1 from q 1,min = -a to q 1,max = a with varying p 1 , such that all rays hit q 4,max = b, these rays are shown in Fig. 2(b), the boundary segments ∂S 2 1,4 and ∂T 2 4,1 are depicted in blue (lines segments labeled with d). Likewise, ∂S 3 1,4 and ∂T 3 4,1 are obtained tracing out line 4 from q 4,max = b to q 4,min = -b by rays leaving q 1,max = x 1,r = a with varying p 1 . These rays are shown in Fig. 2(c), ∂S 3 1,4 and ∂T 3 4,1 are the red line segments labeled with e. Finally, ∂S 4 1,4 and ∂T 4 4,1 are given tracing out line 1 from q 1,max = a to q 1,min = -a with varying p 1 , such that all rays hit q 4,min = -b, these rays are shown in Fig. 2(d), ∂S 4 1,4 and ∂T 4 4,1 are the green lines segments labeled with f. We remind that we use the notation (x, z) for the Cartesian coordinates system of real space, while phase space has (q, p) coordinates. It is worth noting that q 1,min = x 1, , q 1,max = x 1,r , q 4,min = x 4, and q 4,max = x 4,r .
For the two-faceted cup there is an analytic expression for every line segment ∂S  These rays are located on a vertical line segment in S i as only the p i -coordinate changes and on a curved line in T k as both the target position and direction vary. The analytic expressions for ∂S 1 i,k and ∂T 1 k,i are where we have indicated withr i,k (t) the normalization of the ray in Eq. (2.6) and, ν i and ν k are the normalized inward normals to lines i and k, respectively. Note that, sin τ i = |ν i ×r i,k (t)| and sin τ k = |ν k ×r i,k (t)|. Likewise, the boundaries ∂S j i,k and ∂T j k,i are calculated for every j ∈ {2, 3, 4} and ∂S i,k and ∂T k,i are found using Eq. (2.5).
In 3,4 and (∂T i,k ) i =k=1,2,3 are depicted in blue and red, respectively. The source and target PS of lines 2 and 3 have some empty regions. These parts correspond to the regions formed by the rays that either go back to the source or are emitted from the target. As explained in Sect. 2.2, we do not consider these regions, see Eq. (2.5). We observe that, because of the symmetry of the optical system, S 3 is the mirror image of S 2 after reflection in the central point (q, p) = (-9.5, 0) and translation from (q, p) → (q + 19, p). Likewise T 3 is the mirror image of S 2 after the same reflection and translation. In the next section, we show how the phase spaces are related to each other and we define the target photometric variables on T 4 .

Target photometric variables
In this section we explain how to compute the target photometric variables in PS. In the following, to simplify the notation, we indicate the target coordinates of the rays on T 4 with (q, p) instead of (q t,4 , p t,4 ). The intensity I along a given direction p ∈ [-1, 1] in target phase space T 4 is defined as a function of the luminance L(q, p): Note that the intensity is a function of p = sin(τ ) instead of τ . The parts of T 4 that are illuminated by S 1 correspond to parts with positive luminance, for the other parts the luminance is equal to zero. Assuming positive luminance on S, the following relations hold: Once a ray leaves the source S it can hit the reflectors several times before hitting the target T. To relate S and T, a map M 1,4 : S 1 → T 4 is introduced such that M 1,4 (q s,1 , p s,1 ) = (q, p).
As not all parts of T 4 are illuminated by the source S, the map M 1,4 is not surjective. Therefore, we need to determine the subsets of T 4 illuminated by S corresponding to the regions where the luminance is positive. To this purpose, we consider two different kinds of maps. The first map relates the coordinates of the source and the target PS of two different lines, we call it the propagation map. The second map relates the coordinates of the target and the source PS of the same line, we call it the reflection map. In particular, given two lines i and j with i = j, the propagation map P i,j : S i,j → T j,i relates S i,j with T j,i and, it is defined as follows: where q t,j is given by the x-coordinate of the intersection point between the ray and line j, and p t,j is computed considering the direction of the incident ray with respect to the normal of line j. For one single line j, the reflection map R j,k,h : T j,k → S j,h relates the regions T j,k ⊂ T j and S j,h ⊂ S j . To simplify the notation, from now on we omit the dependence of R j,k,h from k and h, i.e., R j,k,h = R j . The reflection map is defined as follows: where p t,j changes according to the reflection law and q t,j = q s,j as R j maps the target PS into the source PS of the same line j. Using a procedure similar to the ray transport matrices approach (see [5], Chap. 6), the map M 1,4 is described by the composition of mappings P i,j and R j defined in Eqs. (2.11) and (2.12), respectively. This composition depends on the path followed by the rays where we refer to a path as the sequence of lines that a ray hits during its propagation from S to T. We indicate with M 1,4 ( ) the map M 1,4 restricted to path and with R( ) ⊂ T 4 the regions on T 4 formed by the rays that follow path . Considering all the possible paths from S to T, all the regions R( ) with positive luminance on T 4 can be determined.
To clarify this concept, we provide the following example. Consider a ray that is emitted from the source (line 1), first hits the left reflector (line 2) and finally reaches the target (line 4). The path followed by this ray is defined as = (1, 2, 4) and the corresponding map M 1,4 ( ) : S 1 → R( ) that describes the propagation of all rays that follow the path is defined by: which can be written as: (2.14) In general, to construct the map M 1,4 ( ) we need to know its corresponding path . To determine all possible paths , instead of tracing the rays from S to T, we start considering the rays in T 4 . In particular, along a given direction p ∈ [-1, 1] we consider the intersection points between the line p = const and (∂T 4,i ) i=1,2,3 . These points are traced back to line i from which they are emitted and their corresponding coordinates on S i and T i are computed. This is done applying sequentially the maps P -1 i,4 : T 4,i → S i,4 and R -1 i : S i → T i . Then the same procedure is repeated considering these new coordinates on T i . The computation stops either when the points found are emitted from the source, that is when they are located on S 1 , or when they reach again the target, that is when they are located on T 4 . If a ray reaches S 1 , then a path from S to T is found. If a ray reaches again the target T 4 , then we conclude that it is not emitted by S and therefore, it is located inside the parts of T 4 with luminance equal to zero.
Finally, the inverse M -1 1,4 ( ) of the map M 1,4 ( ) is constructed for every possible path . The map M -1 1,4 ( ) is the composition of the inverses of the propagation and the reflection maps in reverse order according to the path . For instance, for path = (1, 2, 4), M -1 1,4 ( ) is given by: Tree that describes how to detect all the possible paths from S to T The steps of the procedure are shown in the graph in Fig. 6 where the map in Eq. (2.15) is written in red.
Using the procedure explained above, given a ray with coordinates (q, p) ∈ T 4 we can establish whether it is located inside one of the regions R( ) with positive luminance or not. In case the ray is inside a region R( ), its corresponding coordinates (q s,1 , p s,1 ) ∈ S 1 are obtained using M -1 1,4 ( ), where is the path followed by this ray. Equation (2.10a)-(2.10b) becomes: for some path connecting S and T. Assuming a Lambertian source and employing conservation of luminance along a ray (see [1], Chap. 16), we have that L is a positive constant inside R( ) and it has no contribution on the other parts of T 4 . Indicating with q min ( , p) and q max ( , p) the minimum and maximum position coordinates of the intersection points between the boundaries ∂R( ) and line p = const, Eq. (2.9) reduces to: where the sum is over all the possible paths and the second equation holds as we assume L = 1 in R( ). Note that for a given ray with corresponding coordinates (q, p) on T 4 , only one path is possible as we are assuming that all lines are reflective lines. Because of this, the regions R( ) do not overlap each other, i.e.
where the intersection is over all the possible paths. In the next paragraph the details of the procedure to compute the coordinates q min ( , p) and q max ( , p) are explained.

The inverse ray mapping method and the structure of the algorithm
The goal is to determine the intensity of the light that reaches the target with a given direction p = const. Since we assume a Lambertian source, this is equal to the sum of the lengths of the line segments given by the intersection of the line p = const and the support of L (see Eq. (2.17)). To determine these line segments, a recursive procedure is developed. The procedure starts on T 4 with a given direction p = const and with the parallel rays corresponding to the end points (q min , p) = (-b, p) and (q max , p) = (b, p). We set the initial intensity I(p) = 0 along direction p = const. Considering the intersection between the line p = const and the boundaries (∂T 4,i ) i=1,2,3 three intervals are found. Each interval corresponds to rays emitted by line i (i ∈ {1, 2, 3}). The rays corresponding to the end points of these intervals are traced back from T 4 to S i and subsequently to T i where i is the line from which the rays are emitted. Then, another interval of parallel rays along the corresponding direction in T i has to be considered and the intersection points between the line p = p t,i and ∂T i,j (with j = i, 4) are calculated, where p t,i is the new direction of the rays traced back. The procedure continues recursively until the source is found. Before explaining the details, let us introduce some notation. The role of the variables we introduce will become clear later on in this section. The coordinates in T j of the rays traced back from line i = j to line j are indicated with (q 1 t,j , p t,j ) and (q 2 t,j , p t,j ). The minimum and the maximum position coordinates are q min t,j = min{q 1 t,j , q 2 t,j } and q max t,j = max{q 1 t,j , q 2 t,j }, respectively. The coordinates of the intersection points of p = p t,j with boundaries ∂T j,i need to be determined for every i ∈ {1, 2, 3} and j ∈ {2, 3, 4} with j = i. They are indicated with (u min j,i , p t,j ) and (u max j,i , p t,j ) where u min j,i < u max j,i . Since not all the rays whose corresponding coordinates are located inside the segment [q min t,j , q max t,j ] with direction p = p t,j follow the same path, the intersection segment [ 3. If i = 1, the coordinates (v min 4,1 , p) and (v max 4,1 , p) are equal to the coordinates (q min ( , p), p) and (q max ( , p), p), respectively, of the rays located on the boundary ∂R( ) with = (1, 4). All the parallel rays with direction coordinate p and q-position coordinate with u min 4,1 ≤ q ≤ u max 4,1 are emitted by the source and they directly hit the target.
Update the intensity using Eq. (2.17) 4. If i = 1, continue with the following steps.
6. Update the path = (i, 4). 7. Determine q min t,i = min{q 1 t,i , q 2 t,i } and q max t,i = max{q 1 t,i , q 2 t,i }. 8. Calculate the intersection points (u min i,j , p) and (u max i,j , p) between line p t,i and ∂T i,j for every j ∈ {1, 2, 3} with j = i. 9. Since not all rays whose corresponding coordinates are located inside the segment [q min t,i , q max t,i ] follow the same path, compute the intersection segment (c) Update intensity with q min = min{q 1 ( , p), q 2 ( , p)} and q max = max{q 1 ( , p), q 2 ( , p)}. To clarify the technique, we make an example that describes how the target intensity along direction p = -0.2 is calculated. From Fig. 7(a) to Fig. 7(h) the steps used in this example are shown. A detailed description of those figures is given in the following.
The procedure starts with the rays with direction p = -0.2 on T 4 , where q min = -b and q max = b are the left and the right end points of the target T, respectively. The intersection We start from i = 1. Therefore the coordinates (u min 4,1 , p) and (u max 4,1 , p) of the intersection points between line p = -0.2 and the boundary ∂T 4,1 are computed and these points are depicted in Fig. 7(a). The source is now reached because i = 1 and, one possible path is found. The points (u min 4,1 , p) and (u max 4,1 , p) are located on the boundaries of the region formed by the rays that leave the source and directly hit the target, that is the rays located on  (1, 4). Therefore, the contribution to the intensity formed by the rays that follow the path 1 = (1, 4) is given by u max 4,1u min 4,1 . We continue with i = 2. The boundary ∂T 4,2 is considered in order to find other paths. The intersection points (u min 4,2 , p) and (u max 4,2 , p) of line p = -0.2 with the boundary ∂T 4,2 are calculated. They are depicted in Fig. 7(b)   The map R -1 2 • P -1 2,4 is depicted in red in Fig. 6. The minimum q min t,2 = min{q 1 t,2 , q 2 t,2 } and the maximum q max t,2 = max{q 1 t,2 , q 2 t,2 } are calculated. The points with coordinates (q min t,2 , p t,2 ) and (q max t,2 , p t,2 ) are depicted in Fig. 7(c) where p t,2 = 0.82. To understand whether the corresponding rays are illuminated or not by the source, the preceding procedure used for T 4 is now applied to T 2 along direction p t,2 = 0.82.
If the source is reached, then a valid path = (1, 3, 2, 4) is found. Using the inverse of the propagation map, we compute P -1 1,3 q min t,3 , p t,3 = q 1 s,1 , p s,1 , P -1 1,3 q max t,3 , p t,3 = q 2 s,1 , p s,1 . (2.25) The The coordinates (q 1 , p) = (q 1 ( , p), p) and (q 2 , p) = (q 2 ( , p), p) located on ∂R( ) in T 4 are found. Introducing q min ( , p) = min{q 1 , q 2 } and q max ( , p) = max{q 1 , q 2 }, the contribution to the intensity due to the rays that follow the path is given by If no intersection points are found, then the rays traced are not emitted by the source, therefore no contribution to the intensity needs to be added. This is, for instance, the case of rays with coordinates (v min 2,3 , 0.82) and (v max 2,3 , 0.82) on T 2 in Fig. 7(d). Below we explain this case in detail.
Considering the PS T 3 and the direction p t,3 = 0.91, we note that there are no intersection points between line p t,3 = 0.91 and both ∂T 3,1 and ∂T 3,2 . Indeed, the whole segment [q min t,3 , q max t,3 ] is outside both T 3,2 and T 3,1 . Because of this, all the rays with q-coordinates inside the interval [q min t,3 , q max t,3 ] and with direction p = p t,3 are not illuminated by the source and no new real path is found.
Finally, the recursive procedure is applied to T 4,3 . The first step is depicted in Fig. 7(h). We decided not to show all the steps for T 4,3 as they are similar to those used for T 4,2 and explained above.
Finally, to compute the intensity along another direction p k ∈ [-1, 1] on T 4 , the procedure explained for p = -0.2 is repeated for p = p k . In this way we find all the possible paths and the regions R( ) with positive luminance on T 4 . Furthermore, considering every time the coordinates located on the boundaries of the regions T i,j for every j, also the boundaries ∂R( ) are determined for a given path as well as the coordinates q max ( , p) and q min ( , p) for every p ∈ [-1, 1]. In Algorithm 2.1 they are the main steps to calculate the intensity I(p) along a given direction p = p k in T 4 , where for the first step we take j = 4.
In the next section we provide the numerical results for the two-faceted cup.

Numerical results for the two-faceted cup
To demonstrate the accuracy of the method, a comparison with the ray tracing approach is provided. In particular, we compare our method with MC ray tracing. The MC intensity is computed tracing randomly a large number of rays from the source to the target of the system. Then, a partitioning P : -1 = p 0 < · · · < p Nb = 1 of the interval [-1, 1] is considered, where Nb indicates the number of bins in the partitioning. The MC intensity along the direction p ∈ [p h , p h+1 ] is given by the ratio of the number of rays Nr[p h , p h+1 ] that arrive at the bin [p h , p h+1 ] and the total number of rays Nr[-1, 1] that arrive at the target, that is: for every p ∈ [p h , p h+1 ]. Hence, the MC intensity is piecewise constant. Note thatÎ MC is normalized. In order to compute the intensity distribution for all the directions, the partitioning P of the target is considered. Then, the procedure explained above is repeated for (p h+1/2 = 1 2 (p h+1 + p h )) h=0,...,Nb-1 and the MC intensity is calculated over every bin. The profile of the MC intensity is depicted in Fig. 9 with a blue line. There the intensity is calculated tracing 10 7 rays and taking Nb = 100.
Next, we compute the intensity at the target employing the PS ray mapping method. Using the procedure explained in Sect. 2.5, we are able to detect all the possible paths that a ray can follow during the propagation through the system. For the two-faceted cup 5 different paths are found. Given a path , the coordinates (q min ( , p h ), p h ) and (q max ( , p h ), p h )

Algorithm 2.1 Recursive procedure for the intensity calculation
Initialize j = 4, q min t,4 = q min = -b, q max t,4 = q max = b, p t,4 = p = const, = (4). 1: procedure INTENSITY COMPUTATION(j, q min t,j , q max t,j , p t,j , ) 2: for i = 1, 2, 3 do 3: if i = j then 4: Compute the intersection points (u min j,i , p t,j ) and (u max j,i , p t,j ) between p t,j and ∂T j,i 5: if (i = 1) & (i = 4) then 8: Apply  Target phase space for the two-faceted cup divided into 100 bins. Five different paths are found. The rays with coordinates (q min , p) and (q max , p) in T 4 that are located at the boundaries ∂R( ) are depicted with dots, the color of the dots depends on the path followed by the rays. Using the ray mapping method, only these rays need to be traced from S to T for the intensity computation Figure 9 Intensities for the two-faceted cup with three different approaches. The red line shows the exact intensity. The blue line depicts the intensity computed with MC ray tracing with 10 7 rays. The dotted green line shows the intensity found with the new method of the rays located on ∂R( ) are determined for every p = (p h ) h=0,...,Nb where the values p h are given from the partitioning P used for MC ray tracing. These rays are depicted in Fig. 8, where all the rays that follow the same path are shown with the same color. The PS intensity is obtained from Eq. (2.17). Note from Fig. 8 that only 2Nb rays need to be traced through the system for the intensity computation. The averaged normalized PS intensity is given by: for h = 0, 1, . . . , Nb -1, where the integrals in the previous equation are calculated using the trapezoidal rule. The profile of the PS intensity is depicted in Fig. 9 with the dotted green line. For the two-faceted cup, the intensity can be computed analytically. This intensity is taken as the reference intensityÎ ref and, it is depicted in Fig. 9 with a red line. The results in Fig. 9 show that all the intensity's profiles are all similar to each other. Therefore, we can claim that our method computes the intensity correctly.
In order to compare the speed of convergence of the two methods, we consider the error between the approximate intensitiesÎ A (A = MC, PS) and the exact intensityÎ exact =Î ref .

Figure 10
Error between the approximated intensities and the reference intensity as a function of the CPU time (in seconds). The convergence of the ray tracing approach is depicted with the blue line. The error decreases increasing the number of rays traced. The green dot shows the error obtained with the PS method This error is defined as: As the accuracy of the ray tracing method depends on the number of rays traced, the MC intensity is calculated for an increasing numbers of rays traced through the system and, the error between the approximate intensity and the reference intensity is computed using Eq. (3.3). In Fig. 10 the behavior of the MC error as a function of the CPU-time is depicted with the blue line. Increasing the number of rays the MC error decreases proportionally to the inverse of the square root of the number of rays traced. Next, the PS error is computed using Eq. (3.3). It is depicted in Fig. 10 with the green dot. From the numerical results shown in Fig. 10 we can conclude that the PS ray mapping method is able to compute the output intensity of the two-faceted cup exactly. Also, it is much faster than the classical ray tracing approach when an error smaller than 10 -4 is required.

Generalization of the method to the multi-faceted cup
The method can be generalized to more complicated optical systems. In particular, it can be used for all systems formed by straight line segments. The goal of this section is to show the generalization of the method to the multi-faceted cup which is a system with many left and right segments as reflectors. The design of this system is explained below. A multi-faceted cup is an optical system formed by a source, a target and Nl-2 reflectors, where Nl is the number of optical line segments that form the system. Defining a Cartesian coordinate system (x, z), the multi-faceted cup is symmetric with respect to the optical axis (z-axis). An example of this system is depicted in Fig. 11   The z-coordinates of every end point of the line segments 2, . . . , 21 are given substituting their x-coordinates into the equation of the parabola whose symmetric axis is equal to the z-axis and that passes through the points (-17, 40) and (17, 40). The 20-faceted cup is now well defined and can be seen as an approximation of a parabolic reflector. All the optical lines i with i ∈ {1, . . . , 22} are located in air, therefore the refractive index n i = 1 for every i.
Similarly to the two-faceted cup, also for the multi-faceted cup we define the phase spaces of all the lines i ∈ {1, . . . , Nl} as in Sect. 2.5, for the 20-faceted cup Nl = 22. Note that we always choose the index of the target equal to the index of the number of lines that form the system Nl. For the system in Fig. 11, 42 different phase spaces need to be considered. In general, for a system formed by Nl straight line segments, 2Nl -2 phase spaces are considered. For all the systems formed by straight line segments, the boundaries (∂S i,j ) i =j=2,...,Nl and (∂T i,k ) i =k=1,...,Nl-1 of the regions that form every PS can be determined as explained in Sect. 2.3.
The boundaries (∂T Nl,k ) k=1,...,Nl-1 for the 20-faceted cup are depicted in Fig. 12 with red lines. All the possible paths that the rays can follow when propagating within the 20faceted cup are determined using the same algorithm developed for the two-faceted cup and explained in Sect. 2.5. As the number of optical lines increases, the number of possible paths increases as well. Therefore, we have to construct a more complicated tree than the one in Fig. 6. Despite this, the algorithm explained in the previous section still works fine and, also for the multi-faceted cup we are able to determine all the possible paths and all the regions R( ) with positive luminance at target PS T Nl . Assuming a Lambertian source, only the rays located at the boundaries of these regions need to be computed. Therefore, for a given direction p = const only the position coordinates q min ( , p) and q max ( , p) of the intersection points between the boundaries ∂R( ) and the line p = const are needed for every possible path . Finally, the target intensity I PS (p) along the direction p is obtained employing Eq. (3.2). Numerical results for a 20-faceted cup are given in the next section.

Numerical results for the 20-faceted cup
In this section the results for the 20-faceted cup are presented. We compute the target intensity both with the inverse ray mapping method and MC ray tracing. The same partitioning P of the interval [-1, 1] used for the two-faceted cup is considered. The normalized MC intensity and the normalized PS intensityÎ PS are computed. Their profiles are depicted in Fig. 13 with a blue line and a green dotted line, respectively. They are compared with a reference intensityÎ ref (red line) which is computed using MC ray tracing with 10 8 rays. Note that the intensity profile in Fig. 13 is more concentrated around the direction p = 0 than the intensity of the two-faceted cup (see Fig. 9). In particular, increasing the number of left and right reflectors the intensity profile becomes more and more peaked around the center approaching the profile of a parabolic reflector (see [3]).
In order to show the performance of the new method, we calculate the error between the approximate intensitiesÎ A (A = MC, PS) and the reference intensityÎ ref from Eq. (3.3). In Fig. 14 the speed of convergence for MC ray tracing is shown in blue. The better accuracy (the right most red point in Fig. 14) is obtained tracing 10 8 rays. Also, the PS intensity is computed several times increasing every time the number of bins used in the trapezoidal rule to approximate the integrals in Eq. (3.2). This has to be done in order to find the averaged PS intensityÎ PS over every bin. The error for the ray mapping method is calculated for all the approximated intensities. Increasing the number of bins in the trapezoidal rule, the PS error decreases. We remark that the PS method gives the value of the intensity pointwise, therefore we can compute the PS intensity without numeric integration. Nev- Figure 13 Intensity for the 20-faceted cup. The red line shows the reference intensity computed using MC ray tracing with 10 8 rays. The blue line depicts the MC intensity with 10 7 rays. The green and dotted line shows the intensity found with the ray mapping method Figure 14 Error of the approximated intensities as a function of the CPU time. The convergence of MC ray tracing is depicted with the blue line. The PS error is shown with the green dots. The ray mapping method is more accurate than MC ray tracing and it is faster in case an error smaller than, say 10 -4 , is desired ertheless, we calculate the averaged intensity because we want to compare it with the MC intensityÎ MC .
The error convergence is depicted in Fig. 14 with the red line. Since all the boundaries of the regions in PS are calculated exactly, the PS intensity is analytic. From Fig. 14 we observe that the minimum difference between the reference intensity and the approximated intensity found with the ray mapping method has an order of magnitude of 10 -6 . This is due to the fact that for the 20-faceted cup the intensity cannot be computed exactly. Therefore, we took as reference intensity an intensity computed with MC ray tracing using 7.5 * 10 8 rays which is not the exact intensity. The error between the normalized exact intensityÎ exact and the normalized approximate intensityÎ A is given by: Extrapolating the MC error we obtain an approximation of the difference between the reference solution (MC ray tracing with 7.5 * 10 8 rays) and the exact intensity, this error is depicted in Fig. 14  Therefore, we claim that the error found with the inverse ray mapping method is also due to the MC error. We can conclude that the inverse ray mapping method performs very well also for more complicated systems. Compared to MC ray tracing the new method is not only faster but also much more accurate.

Conclusions
In this paper, we presented an inverse method to compute the intensity of light emitted by the source and received by the target of a given optical system. We tested our method for two-dimensional optical systems formed by straight line segments. The method employs the phase spaces of all the lines that form the system. All these phase spaces are related to each other through two different kinds of maps. A concatenation of these two maps gives a map that connects the coordinates of the rays at the source with those at the target. Employing the inverse of the concatenated map, all the possible paths that rays can follow during their propagation are found. Only the rays located on the boundaries of the regions with positive luminance are traced, where every region is formed by rays that follow the same path during their propagation. Assuming constant luminance, only these rays are needed to calculate the output intensity. We presented numerical results for a simple system: the two-faceted cup. We employ the PS of all the lines that form the system, each of them is divided into regions the boundaries of which are determined exactly. Numerical results show that the exact output intensity is found. We compared our method with MC ray tracing showing significant advantages in terms of the accuracy and the computational time.
Then, we explained how the method can be extended to more complicated optical systems. We took as an example a system with 10 left and 10 right segments as reflectors, the so-called 20-faceted cup. Also for the generalized system, the boundaries of the regions that form every PS are determined exactly. This is true for all the systems formed by straight lines. This allows finding the analytic output intensity. To validate our method, the intensity found using the new technique is compared with the MC intensity. To demonstrate both the accuracy and speed advantages, numerical results are presented also for the 20-faceted cup.
In this paper we presented the method for systems where only reflection plays a role. Currently we are extending the method to more complicated systems formed by curved lines and to systems where also refraction occurs. Future work might investigate systems with Fresnel reflections.

r(s)
Parametrization of a ray in the optical system Part of T i that is illuminated by line j with i = j