Phase space ray tracing for a two-dimensional parabolic reflector

Ray tracing is a technique used in geometric optics for calculating the light distribution at the target of an optical system. Monte Carlo (MC) ray tracing is very common in non-imaging optics. We propose a new ray tracing method that employs the phase space of the source and the target of the system. The new method gives a more accurate target distribution than classical MC ray tracing and requires less computational time. It is tested for two-dimensional optical systems. The results for the paraboloid reflector are provided.


Introduction
Given a light source and an optical system, ray tracing is a forward method to describe how light propagates through the system. Each ray is traced from the source to the target of the system. Every time that a ray hits a surface of the system its intersection point with the surface is computed. At that point the ray changes the direction, which needs to be calculated in order to compute the next surface that the ray hits. This procedure is repeated for a large number of rays. Finally, considering their distribution at the target the output photometric variables are found.
Currently, a very common ray tracing procedure used in optics is Monte Carlo (MC) ray tracing [9]. Using this technique the rays are traced randomly within the system and the target photometric variables are calculated using a probabilistic approach. The method is very robust, easy both to understand and implement. Despite this, it is very computationally expensive and it suffers from MC noise due to random fluctuations at the target. The error decreases at the order of the inverse of the square root of the number of rays traced.
We propose a new ray tracing method which computes the intensity accurately, tracing far less rays compared to MC ray tracing. To test our ideas, we restrict ourselves to two-dimensional optical systems, therefore we work with optical lines instead of surfaces. Note that, to implement the ray tracing procedure, only the coordinates of the intersection point between the ray and the optical lines and the direction of the ray are needed. The space where both these variables are described simultaneously is called phase space (PS). Therefore, in two dimensions the PS of a line is a two-dimensional space where every ray is identified with its position and direction on that line, see [10]. In our method both the source and the target PS are considered.
In general, not the entire target PS is reached by the rays emitted by the source. Hence, the output luminance has jumps from zero to a positive value. Defining a ray path as the sequence of lines that the ray hits, the target PS shows that rays that follow the same path are located in the same region in PS. There exists a map from the source to the target PS that maps rays inside the regions with positive luminance at the source PS into rays located inside the regions with positive luminance at the target PS. The edge-ray principle [7] states that also the boundaries of the regions at the source PS are mapped into the boundaries at the target PS. The boundaries of these regions give location of the discontinuities of the luminance at the target PS. The idea of the method is that, assuming constant luminance, only the rays located on the boundaries of these regions need to be computed. A grid refinement on the source PS is constructed in such a way that only the rays close to the boundaries of the regions in PS are traced within the system. Employing the edge-ray principle the regions with positive luminance on target PS are obtained. We demonstrate that using the source and the target PS representation leads to a significant reduction of the number of rays. Also, the computational time is reduced without loss of accuracy.
In this paper we show the method for a twodimensional parabolic reflector where only the reflection law is considered. Generally, the method can also be implemented for systems where refraction plays an important role, [5]. The intensity at the target of the system is calculated and the result is compared with a reference intensity. Finally, we compare the computational time needed for ray tracing on PS and for MC ray tracing. Numerical results show that our method converges proportionally to the inverse of the number of rays traced. The new method outperforms MC ray tracing, which converges proportionally to the inverse of the square root of the number of rays traced, obtaining a significant reduction in computational time.
The paper is organized as follows. In Section 2, the geometry of the optical system is described and MC ray tracing is briefly introduced. Next, the new ray tracing method is explained in detail. Section 3 shows a comparison of the numerical results obtained with MC ray tracing and ray tracing on PS. The last section of the article provides conclusions and insights about future possible prospectives to extend the method.
2 Ray tracing for a 2D system: the parabolic reflector In this section we present ray tracing for twodimensional optical systems. We focus on the twodimensional system depicted in Figure 1.

Monte Carlo ray tracing
Given a Cartesian coordinate system (x, z), the twodimensional system in Figure 1 is defined in the (x, z)plane. It consists of a source S (line 1), a target T (line 4) parallel to S and two reflectors (lines 2 and 3) which are arcs of the same parabola. The minimum of the parabola is located at the point with x-coordinate equal to 0. S = [-a, a] (with a = 1.9) and T = [−b, b] (with b = 17) are lines perpendicular to the optical axis and are located at z = 0 and z = 40, respectively. All the optical lines are located in air, therefore the index of refraction n i = 1 for every line i = 1, · · · , 4. The optical axis of the system in Figure 1 corresponds to the z -axis. The position coordinates of a ray intersecting line i are indicated with (x i , z i ) and the direction s i is given by the vector s i = (− sin t i , cos t i ), where t i is the angle that the ray forms with the optical axis, measured counterclockwise. Note that t i ∈ (−π/2, π/2) as we consider only forward rays. A ray segment between (x i , z i ) and (x j , z j ) for i = j is parameterized by: where s denotes the arc-length. Assuming a Lambertian source, the input intensity at S emitted in the direction t 1 is given by: where L is the luminance and a is the half width of S. In order to compute the target intensity, we need to find a relation between the intensities at S and T. Hence, we need to know how the optical system influences the direction of the rays when they hit an optical line. To this purpose, the ray tracing procedure is often used in optics. Ray tracing relates the position coordinates (x 1 , z 1 ) and the direction vector s 1 of every ray at the source S with the corresponding position (x 4 , z 4 ) and direction s 4 at the target T. As in the following we will use often the target coordinates of the rays, from now on, to simplify the notation, we write t instead of t 4 and (x, z) instead of (x 4 , z 4 ) for the target coordinates. The ray tracing algorithm can be schematized as follows. For every ray that leaves S with initial position (x 1 , z 1 ) and initial angle t 1 , its ray parametrization is implemented according to Eq. (1). Then, the coordinates (x i , z i ) of the intersection point between the ray and the line i that it hits are computed. The unit normal ν i to the line i at the point (x i , z i ) is calculated to compute the change of direction of the ray. Since all the lines of the system are located in air, only the reflection law plays a role, [4]. Therefore, denoting with t 1 the direction of the incident ray, the direction t 2 of the reflected ray is given by: where the vectors t 1 and t 2 are unit vectors. The procedure explained above is repeated for every line that the ray encounters until it reaches the target and for every ray traced through the system. There are different ways to implement the ray tracing procedure. An often used method is MC ray tracing which calculates the target photometric variables considering a sample of many rays that are traced randomly from S to T. The output intensity is computed as a function of the angular coordinate t and is calculated dividing the target into intervals of the same length, the so-called bins. A partitioning P 1 : −π/2 = t 0 < t 1 < · · · < t Nb = π/2 of the interval [−π/2, π/2] is defined where Nb is the number of bins in P 1 . We remark that, with a slight abuse of notation, we indicated the angular coordinates of the rays at the target with t j instead of t 4,j for every j ∈ {0, · · · , Nb}. The normalized approximated intensity g MC (t) is a piecewise constant function and its value over the j-th bin is the ratio between the number of rays that fall into that bin Nr[t j−1 , t j ) and the total number of rays traced Nr[−π/2, π/2]. Hence, g MC is defined by: Furthermore, the output intensity is computed from the value of the intensity g MC (t j−1/2 ) along the direction t j−1/2 = (t j−1 +t j )/2 for every bin [t j−1 , t j ) j=1,··· ,Nb . The intensity g MC (t j−1/2 ) gives an estimate of the probability that a ray reaches the target with an angle in the j-th interval [t j−1 , t j ) of the partitioning P 1 . This probability P j,∆t is given by: where G(t) is the output intensity (not normalized) and it is measured in lumen per radian [lm/rad]. Note that Nb j=1 P j,∆t = 1. Using the mean value theorem, the integral at the numerator of the previous equation can be written as: Hence, P j,∆t is proportional to the size ∆t = (t Nb − t 0 )/Nb of the intervals, i.e., inversely proportional to the number of bins Nb of the partitioning P 1 . Indicating with Φ = π/2 −π/2 G(t)dt the total flux (measured in lumen [lm]), the error between the intensity G(t j−1/2 ) and the averaged MC intensity Φg MC (t j−1/2 )/∆t is given by: The first term of the right hand side of inequality (7) gives an estimate of how much the averaged intensity 1 ∆t tj tj−1 G(t)dt differs from the exact intensity G(t j−1/2 ). This term is due to the discretization of the target and therefore it depends on the number of bins Nb considered. Substituting G(t) with its Taylor expansion around the point t j−1/2 we obtain that this term is proportional to the square of the size of the bins, therefore the following equality holds: The second part of the right hand side of inequality (7) gives an estimate of the MC error and therefore it depends also on the number of rays traced. In order to show how this term decreases as a function of the number of rays traced, we define the random variable X j (t) as the variable that is equal to 1 if the ray with angular coordinate t is inside the interval [t j−1 , t j ) and equal to 0 otherwise, The Bernoulli trial X j follows a binomial distribution B(1, P j,∆t ). Considering a sample of Nr rays, the variable Y j = Nr k=1 X j (t k ) follows a binomial distribution B(Nr, P j,∆t ), where t k is the angle that the kth ray forms with the optical axis. Then, using the de Moivre-Laplace theorem, we conclude that the variable Y j is approximated by a normal distribution with mean value E[Y j ] = NrP j,∆t and variance σ 2 [Y j ] = NrP j,∆t (1 − P j,∆t ) when a large number of rays is considered, see [8,11]. Thus, the normalized intensity along the direction t j−1/2 is given by: The mean value E[g MC (t j−1/2 )] = P j,∆t and the variance for some C 2 > 0. σ j can be used to give an estimate of the difference between the intensity g MC (t j−1/2 ) and its mean value P j,∆t . Therefore, the second term of the right hand side of relation (7) becomes: for some C 3 > 0, where the approximation holds because σ j gives a measure for the error between g MC (t j−1/2 ) and the probability P j,∆t , [3]. The second equality follows from Eq. (11). To conclude, the MC error over the j-th bin is estimated by: for C 4 > 0. Considering a fixed number of rays, we obtain that the minimal error is reached when Nb ≈ Nr 1/5 . Hence, if 10 10 rays are considered the target has to be divided into 10 2 bins to minimize the MC error. This leads to computational efforts resulting in a very slow procedure.

Ray tracing on phase space
In two dimensions the optical phase space of a given line i is a set containing ordered pairs (q i , p i ) describing the possible position coordinates q i and the direction p i = n i sin(t i ) of a ray intersecting that line, with n i the index of refraction of the medium in which the line i is located and t i the angle that the ray forms with respect to the optical axis (z -axis) measured counterclockwise, [10]. Furthermore, we discard rays parallel to the source S. Hence, we consider angles t 1 in the interval [− π 2 + δ, π 2 − δ] with δ > 0 a small number. Since all the lines i are located in air (n i = 1 for every i ∈ {1, · · · , 4}), from now on we omit the refractive index in the expression for direction p i . We indicate with S = [−a, a] × [−1 + ε, 1 − ε] the source PS and with T= [−b, b] × [−1, 1] the target PS where ε > 0 is a small number. In particular we take ε = 10 −2 . This is physically reasonable, the rays that leave the source with an angle almost equal to π/2 and −π/2 have multiple reflections with the left and the right reflectors, respectively. These rays hit the extremes of the target giving a very small contribution to the output intensity.
The distribution of the rays in PS gives information about the path Π that they follow, where we refer to a path as the sequence of the lines that the rays encounter during their propagation from S to T. The regions in PS formed by the rays that follow the path Π are indicated with R 1 (Π) and R 4 (Π) in S and T, respectively. For the system in Figure 1 multiple reflections of the rays along the parabolic mirrors 2 and 3 can occur. Hence, many different paths are found. The number of paths depends on the geometry of the system and cannot be determined a priori.
In the following we refer to the coordinates on the target PS with (q, p) instead of (q 4 , p 4 ). Ray tracing in PS takes advantage of the fact that there exists an optical map M: S → T such that M(q 1 , p 1 ) = (q, p) for every (q 1 , p 1 ) ∈ S. For the parabolic reflector in Figure  1 and for others optical systems, the map M is not continuous. Despite this, given a path Π, the map M(Π): R 1 (Π) → R 4 (Π) which maps the regions R 1 (Π) in S into the regions R 4 (Π) in T is a continuous and bijective map. The edge ray principle guarantees that M(Π) maps R 1 (Π) into R 4 (Π) preserving topological features. In particular, the boundaries ∂R 1 (Π) are mapped into the boundaries ∂R 4 (Π), see [7,2,6]. Employing the maps M(Π) for all the possible paths Π, the output photometric variables can be calculated. For instance, the luminance L(q, p) at the target PS is given by: for some path Π. Since the luminance is conserved along a ray and a Lambertian source is considered, the luminance is constant inside the regions R 4 (Π). The target intensity along a given direction p = const is computed through an integration of the target luminance L(q, p) and it is defined in T by: Eq. (15) implies that, assuming a Lambertian source, the problem of computing the target intensity is reduced to the problem of calculating the boundaries ∂R 4 (Π) for all the possible paths Π. To this end, rather than choosing randomly the rays that leave the source, as MC ray tracing does, a triangulation on S is defined according to the paths followed by the rays. The procedure starts tracing four rays located exactly at the corners of S. Then, joining two opposite vertices of the grid two equal right-angled triangles are defined. The paths followed by the rays corresponding to the vertices of these triangles are considered. If there are at least two different paths, then one or more boundaries ∂R 1 (Π) are expected to cross the triangles. In order to determine the location of these boundaries, three more rays corresponding to the middle points of every side of the triangle are added to the sample of the rays traced. Moreover, the paths of these three rays are stored. The procedure explained above is repeated for every new triangle. Every refinement step leads to four new right-angled triangles, see  When a boundary is detected, smaller and smaller triangles are defined around it. The rays at the vertices of the triangles close to the boundaries have always a different path. To stop the refinement once it is fine enough, two parameters ε q 1 ,min and ε p 1 ,min are introduced for the q and p-axis, respectively. If all the vertices of the triangle correspond to rays that follow the same path, there is no need to refine new triangles unless its size is too big. There is indeed the possibility that a region formed by rays that follow a given path Π 1 is located completely inside a triangle whose vertices are all related to another path Π 2 . To avoid this, two other parameters ε q 1 ,max > ε q 1 ,min and ε p 1 ,max > ε p 1 ,min are defined for the q-axis and the p-axis, respectively. When the length of the sides of the triangle are greater than both ε q 1 ,max and ε p 1 ,max , a new triangle is defined also in case the rays corresponding to the vertices of the triangle follow the same path. The refinement procedure terminates either when the paths of the rays corresponding to the vertices of the triangle are equal and the triangle's sides are smaller than the parameters ε q 1 ,max and ε p 1 ,max or when those paths are different and the triangle's sides are smaller than the parameters ε q 1 ,min and ε p 1 ,min . Hence, the values of the parameters ε q 1 ,min , ε p 1 ,min , ε q 1 ,max and ε p 1 ,max establish the size of the smallest and largest triangle involved in the triangulation and the number of rays traced. Decreasing the values of these parameters, the number of rays traced increases. Also, the values of the parameters determine the number of possible paths that we are able to detect. Considering a very fine grid, also rays which are close to the edge of the source and form an angle almost equal to −π/2 or π/2 with the optical axis are traced. These rays have multiple reflections with the reflectors. Therefore, the finer the triangulation is the more paths are detected. To select the correct parameters, conservation of the energy is employed, see [5] for more details.
With the refinement procedure, the rays closest to the boundaries ∂R 1 (Π) are selected from an initial ray set and more rays in their vicinity are created to get progressively better estimates of them. In Figure 3, the triangulation refinement of S with the values of parameters ε q 1 ,min = 0.01, ε q 1 ,max = 0.8, ε p 1 ,min = ε q 1 ,min /2 and ε p 1 ,max = ε q 1 ,max /2 is shown. With these values, approximately 1.5 · 10 4 rays are traced and 11 different paths are found. The edge-ray principle is employed in order to obtain the refinement at the target. Note that, by construction, all the triangles that overlap the boundaries have two vertices corresponding to the same ray path and one vertex corresponding to another path. Given a path Π, all the rays of the triangulation that belong to the region R 4 (Π) are considered. Then, to approximate the boundary ∂R 4 (Π) the coordinates of the vertices of the smallest triangles corresponding to the coordinates of the rays that follow the path Π are joined obtaining the approximated boundary ∂R 4 (Π) at target PS, T.
Using the triangulation refinement the distribution of the rays both on source and target PS of the system is found. The phase spaces S and T are represented at the top and bottom of Figure 4, respectively. Every dot in PS is a ray and the dots depicted with the same colour refers to rays that follow the same path. Figure 4 shows that S and T are partitioned into regions. All the rays that follow the same path are located in the same region in PS. S is almost completely covered by those regions (as we assume that the light is uniformly distributed at the light source), while T has some parts not covered by any ray. Nevertheless, because ofétendue conservation, the total area covered by the regions in S and T is always the same. Note that the scale of the q-axis for the source and target phase spaces is different.
To conclude, the target intensity along a direction p = constant is computed. Denoting with q min (p, Π) < q max (p, Π) the position coordinates in T of the rays located on the intersection points between the boundaries ∂R 4 (Π) and the line p = constant, Eq. (15) becomes: where the second equality holds assuming a Lambertian source with L = 1.
With a slight abuse of notation, from now on we do not specify the dependence of the variables q min and q max on the path Π and the direction p. Note that Eq. (16) is valid when only two intersection points q min and q max between the line p = constant and the boundaries rays are traced through the system. The triangulation refinement is implemented using the parameters ε q 1 ,min = 0.01, εq 1 ,max = 0.8, ε p 1 ,min = ε q 1 ,min /2 and εp 1 ,max = εp 1 ,max/2. 11 different paths are found. All rays that follow the same path are depicted with the same color. The boundaries are depicted with black lines and they are determined connecting all the rays that follow the same path. ∂R 4 (Π) are found (which is the case for the system in Figure. 1). In case ∂R 4 (Π) has more than two intersection points with p = constant, that equation should be modified accordingly. This situation occurs when more complicated optical systems are considered and the shape of the regions R 4 (Π) is not very regular (see as an example the case of the TIR collimator described in [5]).
To compute the intensity for all possible directions, a uniform partitioning P 2 : − 1 = p 0 < p 1 < · · · < p Nb = 1 of the interval J = [−1, 1] is considered. The intensity for all the directions p h defined in P 2 with h = 0, · · · , Nb is obtained using Eq. (16).
To validate the PS method, the MC intensity is calculated using the same partitioning P 2 . The number of rays that fall into each bin [p h−1 , p h ] h=1,··· ,Nb is calculated for all h ∈ {1, · · · , Nb}. The averaged and normalized intensity in the direction p h−1/2 ∈ [p h−1 , p h ) is then given by:ḡ number of rays. Since MC ray tracing computes the averaged intensity over every bin of the partitioning, the PS intensity G PS has to be averaged in order to compare the two methods. The averaged PS intensity g PS (p) along a given direction p ∈ P 2 is given by the average of G PS over the bin to which p belongs. Then, the averaged and normalized intensityḡ PS is calculated dividing g PS by the total flux, see [5] for more details.
The numerical results of the PS method are given in the next section.

Numerical results
Ray tracing on PS is compared with MC ray tracing and numerical results are shown in this section. We compute the target intensities for both methods and we compare them with a reference intensity. For very few and simple systems the target intensity can be computed analytically but in the case of the system in Figure 1 that is not possible. Therefore, we take as a reference intensity the result of a MC simulation where 10 8 rays are traced. The intensities are calculated as a function of the angular coordinate p = sin(t). Figure 5 shows that the PS intensity (dotted red line) and the MC intensity (dotted blue line) approximate the intensity very well. To obtain the same accuracy far less rays need to be traced with the PS procedure. The intensities shown in Figure 5 are obtained tracing around 1.5 · 10 4 rays for PS ray tracing while approximately 10 7 rays are used to compute MC intensity. In addition, we calculate the error between the approximate intensityḡ A (A = PS, MC) and the reference intensity using the following relation: Eq. (18) is implemented several times both for the MC and PS method. Every time the number of sampling rays is increased in order to reduce the error. In Table  1 and Table 2 the errors values are given for MC and PS ray tracing, respectively. These values are obtained considering a fixed number of bins, in particular we chose Nb = 100. In Figure 6, a comparison between the speed of convergence between MC ray tracing and PS ray tracing is shown. In these pictures, we depict with a blue line the MC-error and with a red line the PS-error. We provide errors both as a function of the number of rays traced and as a function of the CPU-time.
In the picture at the top of Figure 6, the errors are depicted as a function of the number of rays traced in a logarithmic scale. The slope of statistical error (the black line) calculated computing the standard deviation as explained in Section 2.1 is drawn to demonstrate that MC ray tracing converges proportionally to 1/ √ Nr. Then, drawing a line with slope −1 (the green line) we show that ray tracing in PS converges proportionally to 1/Nr. Note that to obtain an accuracy given by an error equal to, say 1.6 · 10 −5 , around 10 7 ray need to be traced for MC and only around 6.8 · 10 4 rays for the PS method resulting in a reduction of more than 100. In the graph at the bottom of Figure 6, the error as a function of CPU-time is shown both for MC (blue line) and PS (red line) method. With our method we need to take into account the path followed by any ray and compute the boundaries of all the regions formed by the rays that follow the same path. Despite this overhead, we observe big advantages also in terms of CPU-time using our method.

Conclusions
A new ray tracing method which employs the PS representation of source and target of an optical system was presented. The method was tested for a two-dimensional optical system with two parabolic mirrors. The idea of PS ray tracing is to partition the phase spaces in several regions according to the different paths followed by the rays. This allows tracing only the rays close to the boundaries of these regions. Eventually, applying the edge-ray principle, we concluded that the rays traced are located close to the boundaries of those regions also in target PS. These boundaries determine where the luminance has a jump from 0 to a positive value. Assuming a Lambertian source, the output intensity can be computed from the location of the intersection points of the regions with the lines with p = const. As a result, the output intensity is calculated using far fewer rays, and the computational time is reduced. The efficiency of the method was demonstrated using a highly accurate intensity profile as a reference and a comparison with the classical ray tracing approach was illustrated. Numerical results show that the PS method outperforms MC ray tracing. For both methods, a computation of the error between a reference intensity and the approximated intensities is provided. For the PS method the error as a function of the number of rays traced decreases proportionally to the inverse of the number of rays traced versus a speed of convergence proportional to the inverse of the square root of the number of rays traced obtained with MC ray tracing. A computation of the error as function of the CPU-time shows that ray tracing in PS is significantly faster than MC ray tracing. For the optical system considered in this work, only the reflection law was involved. The method is also valid for systems with refraction, see [5]. Furthermore, ray tracing in PS can be generalized to systems with a non-constant luminance.
In the future, we plan to consider systems where also Fresnel reflection is of importance. Our expectation is that the method is suitable also for three-dimensional systems. Future investigations should go in this direction as well.