Concatenated Backward Ray Mapping on the Compound Parabolic Concentrator

Concatenated backward ray mapping is an alternative for ray tracing in 2D. It is based on the phase-space description of an optical system. Phase space is the set of position and direction coordinates of light rays intersecting a surface. The original algorithm is limited to optical systems consisting of only straight surfaces; we generalize it to accommodate curved surfaces. The algorithm is applied to a standard optical system, the compound parabolic concentrator. We compare the accuracy and speed of the generalized algorithm, the original algorithm and Monte Carlo ray tracing. The results show that the generalized algorithm outperforms both other methods.


Introduction
The illumination optics industry deals with the design of optical systems.An optical system consists of a light source, optical components such as lenses and reflectors and a target which is either a receiver or an aperture.Light emitted at the source of the optical system propagates through the system and forms an intensity distribution at the target.The shape of the target distribution depends on the optical system and the intensity distribution at the source.The goal in illumination optics is to obtain a desired intensity distribution at the target.Designing an optical system is an iterative process during which the target distribution is computed many times.Therefore, there is a need for fast and accurate simulation methods.
The target distribution is typically computed using Monte Carlo (MC) or Quasi-Monte Carlo (QMC) ray tracing.MC ray tracing is based on a probabilistic interpretation of the source distribution [2].Many randomly distributed rays are traced from source to target; the intensity distribution is found by dividing the target into equal cells and counting the number of rays in each cell.QMC ray tracing was introduced as a faster alternative to MC ray tracing.The difference between them is that QMC ray tracing distributes the rays along low discrepancy sequences [3].MC ray tracing is a slow and expensive procedure with a convergence rate of O(1/ √ Nr) [4], where Nr is the number of rays traced.QMC ray tracing performs better with a convergence rate of O(1/Nr) [4], but it is still an expensive procedure.
Concatenated backward ray mapping (CBRM) is an alternative to (Q)MC ray tracing in 2D that uses the phase space (PS) [1].The PS of a surface is defined by the position and direction coordinates of all rays that interact with the surface.An optical surface in 2D is a line segment or a curved segment, but we will still refer to it as a surface.The algorithm determines which light rays emitted by the source reach the target at a certain angle by tracing backward; only those rays are then traced from source to target, resulting in a significant reduction in the number of rays needed.Doing this for all angles at the target gives the intensity distribution.Numerical results show that CBRM computes the intensity distribution more accurately and with less computation time than (Q)MC ray tracing [1].However, CBRM requires the optical system to consist of straight surfaces.We generalize CBRM to accommodate curved optical surfaces.
The purpose of this paper is to introduce the generalized CBRM algorithm and apply it to the compound parabolic concentrator (CPC) [5].The CPC (Fig. 1) is a standard optical system which collects light from a Lambertian source and reshapes it to a focused beam.Before we introduce the generalized CBRM algorithm we first explain CBRM and apply it to the two-faceted cup (Fig. 2).The performance of CBRM and generalized CBRM is compared on the CPC.Since CBRM requires the optical system to consist of straight surfaces, it is applied to a discretized CPC.
The structure of this paper is as follows.Section 2 defines the phase space of an optical surface.In Section 3 we describe the concatenated backward ray mapping algorithm.Section 4 explains the generalized algorithm.Section 5 describes the numerical experiments that compare the algorithms on the (discretized) CPC.Section 6 shows the results of the experiments.In Section 7 we draw conclusions.Fig. 1: The CPC with source (1), reflectors (2,3), and target (4).Fig. 2: The two-faceted cup with source (1), reflectors (2,3), and target (4).

Phase Space
Concatenated backward ray mapping is based on the phase-space (PS) [4] description of an optical system.Phase space is the set of position and direction coordinates of rays intersecting a surface.It is a two-dimensional space for 2D surfaces.The position coordinate q ∈ Q is the x-coordinate of the intersection of the ray with the surface.The direction coordinate p ∈ P is given by p = n sin θ, where θ ∈ [−π/2, π/2] is the angle between the ray and the inward facing unit surface normal ν and n is the index of refraction.We consider only optical systems formed by reflective surfaces, therefore refraction is not taken into account and n = 1.The index n will be omitted from now on.PS is indicated with S = Q × P .Every optical surface has a source PS and a target PS.Source PS describes light emitted by the surface and target PS describes light reaching the surface.The light source only emits light, so it only has a source PS; the target only receives light, so it only has a target PS.The optical surfaces of the system are numbered 1, . . ., N where 1 is the light source and N is the target.The source and target PS of surface j are indicated with S j and T j , respectively.The coordinates of a ray reaching surface j is indicated with (q t,j , p t,j ) ∈ T j .After reflection, the ray is emitted by the surface from the same position but with a new direction.The ray now has the coordinates (q s,j , p s,j ) ∈ S j .Note that q t,j = q s,j while p s,j is obtained by applying the law of reflection to the ray.
The phase spaces S j and T j are divided into regions S j,k and T j,l .S j,k is the region of S j containing all light rays emitted by j, illuminating surface k ∈ {2, . . ., N} assuming that j acts as a light source.T j,l is the region of T j containing all light rays emitted by surface l ∈ {1, . . ., N − 1}, illuminating surface j assuming that l acts as a light source.Note that S j and T j may contain empty regions.The boundaries ∂S j,k are connected to the boundaries ∂T k,j for every j ∈ {1, . . ., N − 1} and k ∈ {2, . . ., N} by the edge-ray principle [4].These boundaries can be determined analytically and are formed by four curves: Given two surfaces j and k, ∂S j,k and ∂T k,j are determined as follows.∂S 1 j,k and ∂T 1 k,j are formed by the set of rays originating at the left endpoint of j tracing out surface k. ∂S 2 j,k and ∂T 2 k,j are formed by the set of rays tracing out surface j reaching the right endpoint of k. ∂S 3 j,k and ∂T 3 k,j are formed by the set of rays originating at the right endpoint of j tracing out surface k. ∂S 4 j,k and ∂T 4 k,j are formed by the set of rays tracing out surface j reaching the left endpoint of k.As an example Fig. 3 shows these sets of rays in the two-faceted cup where j = 1 and k = 4, i.e., the source and target.The segments formed by these rays in S 1 are shown in Fig. 4. The segments formed in T 4 are shown in Fig. 5.
Previously [1,4] we only had analytic expressions of the boundaries of Eq. ( 1) when j and k were straight surfaces.Here, we introduce a general analytic expression for each boundary of Eq. ( 1) when surfaces j and k are described by parametric equations.Let surface j be described by the parameterization P j (γ) = x j (γ), z j (γ) (γ min ≤ γ ≤ γ max ) and surface k by the parameterization P k (λ) = x k (λ), z k (λ) (λ min ≤ λ ≤ λ max ), then the rays that form the boundaries ∂S 1 j,k and ∂T 1 k,j are parameterized by The rays form a vertical segment in S j as only the direction coordinate changes; they form a curved segment in T k because both the position and direction coordinates change.The analytic expressions for ∂S 1 j,k and ∂T 1 k,j are We indicate with r1 j,k (λ) the normalization of the ray in Eq. ( 2) and with τ j (γ) and τ k (λ) the normalized tangent vectors to surfaces j and k respectively.The tangent vectors are obtained by rotating the inward facing surface normals νj and νk by an angle of π/2 counterclockwise.Note that the parameter in Eq. (3) corresponds to the surface that is traced out.The expressions for the other boundary segments are similar.

Phase Space of the Two-Faceted Cup
The two-faceted cup depicted in Fig. 2 is a simple optical system consisting of four surfaces.A Lambertian light source (surface 1), two reflectors formed by straight line segments (surfaces 2, 3) and a target (surface 4).A ray leaving the source of the cup can reflect many times between the reflectors before reaching the target.Light that reflects on one of the reflectors always propagates to another surface.The phase spaces of the two-faceted cup can be seen in Fig. 6 and are given by the following expressions: Fig. 3: Rays on the boundaries of the regions S 1,4 and T 4,1 in the two-faceted cup.  of the two-faceted cup consisting of segments ∂T 1 4,1 (green), ∂T 2 4,1 (blue), ∂T 3  4,1 (purple) and ∂T 4 4,1 (orange).

Phase Space of the Compound Parabolic Concentrator
The CPC [5] depicted in Fig. 1 is a standard optical system that collects light from a Lambertian source and reshapes it to a focused beam.The system consists of an aperture (surface 4), a receiver (surface 1) and two reflectors that are segments of parabolas (surfaces 2, 3).We take the receiver as the light source and the aperture as the target, thus considering light traveling in the opposite direction.A light ray reflects on at most one reflector of the CPC.It can however reflect infinitely many times on a reflector.The parameterization of the right reflector follows from the polar equation of a parabola [5] and is given by: where the axis of the parabola is rotated θ radians around the origin and the parabola is translated horizontally along a distance a > 0. Note that the parabola has focus (−a, 0) and intersects the horizontal axis at the point (a, 0).The parameterization of the left reflector is similar but it is rotated −θ radians, it is translated in the opposite direction and it has different bounds for φ.The value of h(φ) is given by: A ray leaving the source of the CPC can reflect many times on a single reflector before reaching the target.The phase spaces of the CPC can be seen in Fig. 7 and are given by the following expressions: 3 Concatenated Backward Ray Mapping For each optical system there exists an optical map M 1,N : S 1 → T N such that M 1,N (q s,1 , p s,1 ) = (q t,N , p t,N ) for every (q s,1 , p s,1 ) ∈ S 1 .All rays that follow the same path Π from source to target form a region R s (Π) ⊂ S 1 and R t (Π) ⊂ T n .A path is the sequence of surfaces encountered by a ray traveling from source to target.The map M 1,N (Π) is the map M 1,N restricted to the path Π and relates R s (Π) to R t (Π).The areas covered by light rays in S 1 and T n are equal because of étendue conservation [5].Empty regions occur where no light rays exist that travel from source to target.The light intensity I(p) at the target for a given p = const is computed by integrating the target luminance over the position coordinate q and is defined by: The light intensity in T n depends on the luminance, which is positive in non-empty regions of PS.Assuming positive luminance on the source gives the following relation: It is partitioned into regions T 2,j where j ∈ {1, 3} containing rays reaching surface 2 and emitted by surface j.It is partitioned into regions T 2,j where j ∈ {1, 2, 3} containing rays reaching surface 2 and emitted by surface j.L(q, p) > 0 The phase spaces of an optical system are connected through maps that relate the coordinates on every PS.Propagation maps P j,k : S j,k → T k,j describe light that travels from a surface j to another surface k; they relate coordinates of S j to T k such that P j,k (q s,j , p s,j ) = (q t,k , p t,k ).Reflection maps R k : T k → S k describe light that reflects on a surface k; they relate coordinates of T k to S k such that R k (q t,k , p t,k ) = (q s,k , p s,k ).Note that q t,k = q s,k .Every map M 1,N (Π) can be described by a composition of propagation and reflection maps.Considering all paths Π from source to target, the positive luminance regions R t (Π) ⊂ T n can be determined.From Eq. ( 9) follows that the luminance for some path Π connecting the source and target is: A Lambertian source emits light with a constant luminance [5].Assuming a Lambertian source with luminance equal to 1, computing I(p) reduces to computing the boundaries ∂R t (Π) of the regions of all possible paths Π.The intensity is given by the sum of the interval lengths formed by the intersections of the support of the luminance and the line p = const.The intersection points between p = const and the boundary ∂R t (Π) have position coordinates q min (Π, p) and q max (Π, p), where q min (Π, p) < q max (Π, p).Using Eq. ( 10), it follows that Eq. ( 8) reduces to: CBRM computes the light intensity I(p) using the phase spaces of all surfaces of an optical system.A parallel light beam is represented by a straight line segment on the line p = const in the source/target PS of a straight surface since the direction coordinate of all rays is the same.The intersections between p = const and the boundaries in PS are computed several times by the algorithm; this is done analytically since p = const is a horizontal line, and we have analytical descriptions of all boundaries of all phase spaces.Therefore, the light beam is required to always stay parallel which in turn requires the optical system to consist of only straight surfaces.
The algorithm uses the map M 1,N (Π) : R s (Π) → R t (Π) for all possible paths Π to compute I(p) for a given direction p = const in T n .To construct the map M 1,N (Π), the corresponding path Π should be known.The algorithm computes all possible paths Π by considering rays in T n along a given direction p ∈ [−1, 1] and tracing them backward recursively in PS [1].The endpoints of the light beam in T n are (q min t,N , p t,N ) and (q max t,N , p t,N ).The line p = const intersects various regions in PS.The intersection points with boundaries ∂T N,j are (u min N,j , p t,N ) and (u max N,j , p t,N ).The intersection segment with region T N,j is given by [v min N,j , v max N,j ] = [q min t,N , q max t,N ] ∩ [u min N,j , u max N,j ] and corresponds to rays emitted by another surface j = N.The endpoints of the intervals are transformed to coordinates (q min t,j , p t,j ) and (q max t,j , p t,j ) in T j by sequentially applying the maps P −1 j,N : T N,j → S j,N and R −1 j : S j → T j .The procedure is repeated in T j .The recursion ends when an intersection segment [v min k,j , v max k,j ] is traced back to S 1 or when an intersection segment is empty.If S 1 is found, the endpoints (q min s,1 , p s,1 ) and (q max s,1 , p s,1 ) of the light beam are traced to T n along Π by applying M 1,N (Π).This gives two points (q min (Π, p), p) and (q max (Π, p), p) at the boundary of a region R t (Π) ⊂ T n with positive luminance.The main steps to calculate I(p) are given in Algorithm 1.
The range of angular coordinates at the target is divided into Ni equidistant intervals with endpoints p m and p m+1 where m ∈ {0, . . ., Ni − 1}.The averaged and normalized intensity Î is given for every interval p m+1/2 = 1 2 (p m + p m+1 ) by: and is computed using the trapezoidal rule.U t denotes the total étendue at the target, which in PS corresponds to the area of the region covered by the light rays [4].The intensity distribution is obtained by plotting p m+1/2 on the horizontal axis and Î(p m+1/2 ) on the vertical axis for each interval.For more details on the CBRM algorithm see [1,4].Fig. 10 shows the first steps of the algorithm for light on p = −0.2 in T 4 on the twofaceted cup.In step 1 (Fig. 10a) we find light traveling directly from source to target.The algorithm updates the intensity and continues to the next region of PS in T 4 .In step 2 (Fig. 10b) we find light traveling from surface 2 (left reflector) to the target.The light is traced to T 2 where it lies on p = −0.82 and the procedure is repeated.In step 3 (Fig. 10c) we find light traveling from the source to surface 2. The algorithm traces this light back to the target and updates the intensity, then it continues to the next region of T 2 .In step 4 (Fig. 10d) we find light traveling from surface 3 (right reflector) to surface 2. The light is traced to T 3 where it lies on p = 0.29 and the procedure is repeated.In step 5 (Fig. 10e) we find light traveling from surface 2 to surface 3. The light is traced to T 2 where it lies on p = 0.41 and the procedure is repeated.In step 6 (Fig. 10f) we find light traveling from surface 3 to surface 2. The light is traced to T 3 where it does not intersect any region of PS meaning it was not emitted by the source.Therefore, the computation for this part of the light beam stops; the algorithm is finished for the light of step 2 reaching the target from surface 2. The process is repeated for the light reaching the target from surface 3.

Intensity Distribution of the Two-Faceted Cup
The cup is a simple system for which we can compute the intensity distribution exactly [6].The target is divided into 100 bins for (Q)MC ray tracing.The intensity in each bin is also computed with CBRM using Eq. ( 12).The intensity distribution found with MC ray tracing, shown in Fig. 8, is noisy and not close to the exact solution.The intensity distribution computed with CBRM, shown in Fig. 8, matches the exact solution precisely.MC ray tracing required 10 4 rays, but CBRM required only 10 3 rays.CBRM is more accurate than MC ray tracing and requires fewer rays to compute the intensity.We can use CBRM to find the boundaries of the positive luminance regions in T 4 ; they are shown in Fig. 9. (q min t,4 , p) and (q max t,4 , p) are the endpoints of the light beam.The intersection points between p = −0.(q min t,2 , p) and (q max t,2 , p) are the endpoints of the light beam.The intersection points between p = −0.82 and ∂T  (e) Light in T 3 on the line p = 0.29.(q min t,3 , p) and (q max t,3 , p) are the endpoints of the light beam.The intersection points between p = 0.29 and ∂T 3,2 are (u min 3,2 , p) and (u max 3,2 , p).The intersection segment has endpoints (v min 3,2 , p) and (v max 3,2 , p) with v min 3,2 = max{q min t,3 , u min 3,2 } and v max 3,2 = min{q max t,3 , u max 3,2 }.
-16 -14 -12 -10 - T 2,1 u min =v min q min q max =u max =v max (f) Light in T 2 on the line p = 0.41.(q min t,2 , p) and (q max t,2 , p) are the endpoints of the light beam.The intersection points between p = 0.41 and ∂T  Input: k = N, q min t,k is the x-coordinate of the left endpoint of surface N, q max t,k is the x-coordinate of the right endpoint of surface N, p t,k = p = const, I(p) = 0, Π = (N).
1: procedure INTENSITY COMPUTATION(k, q min t,k , q max t,k , p, I(p), Π) 2: for j ← 1, N do 3: Compute (u min k,j , p t,k ) and (u max k,j , p t,k ) 5: if [v min k,j , v max k,j ] is not empty then 7: Update the path Π 8: if j = 1 then 9: Compute (q 1 t,j , p t,j 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 } 12: INTENSITY COMPUTATION(j, q min t,j , q max t,j , p t,j , I(p), Π) 13: else 14: if k = N then 15: Compute (q 1 s,1 , p s,1 Compute (q 2 s,1 , p s,1 Compute (q 1 (Π, p), p) = M 1,N (Π)(q 1 s,1 , p s,1 ) 18: Compute (q 2 (Π, p), p) = M 1,N (Π)(q 2 s,1 , p s,1 ) 19: q min (Π, p) = min{q 1 , q 2 } and q max (Π, p) = max{q 1 , q 2 } 20: where q 1 = q 1 (Π, p) and q 2 = q 2 (Π, p) 21: else 23: end for 29: end procedure 4 Generalized Algorithm CBRM only handles parallel beams of light, limiting the algorithm to optical systems consisting of only straight surfaces.We generalize the algorithm to accommodate curved surfaces.Eq. (11) defines the intensity in T N as the sum of the interval lengths formed by the intersections of the support of the luminance and the line p = const.Therefore, we must still consider a parallel beam of light at the target.The direction of all rays of the beam will vary as it reflects on a curved surface.As a result, the light beam is not necessarily parallel at the other surfaces of the system.In PS this means that position and direction coordinates of all rays in the beam can vary.The beam is no longer represented by a straight line segment on p = const, but by a curved segment for which we generally have no analytical expression.
Intersections in PS cannot be computed analytically since there is no expression for the light beam.We instead discretize the line p = const in T n at the start of the procedure by taking equidistant points and connect them with straight line segments.This discretization is the PS representation of the light beam; when traced to other surfaces it forms a discretized curve C.
In addition to the light beam, the phase spaces of the optical system are also discretized.Recall from Section 2 that there are analytic descriptions of all boundaries of all phase spaces.We discretize each boundary in PS by taking points on these boundaries equidistant in the q-direction, and connect them with straight line segments.Fig. 11 shows T 4 of the CPC with light at p = −0.1 and the discretization used by the algorithm.Each discretized PS is stored in a doubly connected edge list (DCEL) [7].The DCEL stores each region of the PS as a face, each point as a vertex and each line segment as a pair of half-edges.All edges are incident to two faces of the DCEL; therefore, they are split into two half-edges such that each half-edge has exactly one incident face.Each face is bound by the vertices and half edges that form the boundary of the PS region.The half-edges connect pairs of vertices and are ordered counterclockwise around the face they bound.The DCEL is a nice data structure to store geometric information.It makes it easy to perform operations such as traversing the boundary of a given face, accessing a face from an adjacent one if a common edge is given or visiting all edges around a given vertex.
Computing an intersection between C and the boundaries in PS reduces to computing the intersection between two straight line segments; one segment of C and one segment of the discretized boundaries.However, since the PS and the light beam are discretized with many segments we have to check for many pairs of segments if they intersect.To solve this we build a KD-tree [8] for each PS.The KD-tree places a bounding box around the PS and subdivides it into increasingly smaller regions storing the boundary segments in its leafs.Non-overlapping regions of PS are stored in the internal nodes and leafs of the KD-tree, i.e., the data structure is a space partition.A boundary segment that intersects different regions of the partition is stored in multiple leafs of the KD-tree since the regions of the KD-tree do not overlap, and they all contain the segment.The tree is a binary tree and every node is split in half along an axis-aligned split plane.We use the surface area heuristic (SAH) [8] to select the best splitting plane for every potential split.With the SAH we compute a cost for all possible split planes of a region of the KD-tree.The split plane with the lowest cost is considered to be the best split plane.A region is split when the cost of the best split plane is less than the cost of not splitting the region.The region is not split in half when the cost of the best split plane is greater than the cost of not splitting the region.
Given a segment of C we create a half line by extending the origin of the segment beyond the outer boundaries of PS.The origin of the segment is the endpoint with the smallest q-coordinate; it is extended along the direction of the segment in increasing q direction.The segments of the boundaries in PS that can be intersected by the half line are in the cells (leafs) of the KD-tree intersected by the half line; the segments in all cells not intersected by the half line can be ignored.We are interested in the intersection point closest to the origin of the half line.Therefore, we first check for intersections in the cell closest to the origin and continue through the other cells along the half line.When an intersection is found between the half line and a PS segment all subsequent cells of the tree along the ray can be ignored.If the intersection point lies on the segment of C then an intersection between C and a discretized boundary is found; if the intersection point is not on the segment of C then the segment is contained inside a PS region.
The luminance for all possible paths Π for a given direction p ∈ [−1, 1] in T n is computed recursively, similarly to the original algorithm.The endpoints of the light beam in T n are C min t,N and C max t,N .The discretized curve C t,N intersects various regions in T n .C N,j is the subset of C t,N intersecting the region T N,j .The intersection points C min N,j and C max N,j with ∂T N,j are computed analytically.Each subset C N,j corresponds to rays emitted by another surface j = N.All segments of C N,j are transformed to a new discretized curve C t,j with endpoints C min t,j and C max t,j in T j .This is done by sequentially applying the maps P −1 j,N : T N,j → S j,N and R −1 j : S j → T j .The procedure is repeated in T j .The recursion ends when a subset C k,j is traced back to S 1 or when a subset is empty, like in the original algorithm.If S 1 is found, all points of C s,1 are traced to T n along Π by applying M 1,N (Π).The endpoints C min (Π, p) and C max (Π, p) of C(Π, p) are on the boundary of a region R t (Π) ⊂ T n with positive luminance.The steps of the generalized CBRM algorithm are given in Algorithm 2.

Intensity Distribution of the CPC
Recall from Section 2.2 that light in the CPC can reflect an infinite number of times, generalized CBRM can run indefinitely as a result.To prevent this we do not consider light that reflects more than 10 times.The reference solution of the CPC is computed with QMC ray tracing using 10 9 rays that reflect no more than 10 times.The target is divided into 110 bins for (Q)MC ray tracing.The intensity in each bin is also computed with generalized CBRM using Eq. ( 12).The intensity distribution computed with MC ray tracing, shown in Fig. 12, is again noisy like for the two-faceted cup.Generalized CBRM is able to match the reference solution closely as can be seen in Fig. 12. MC ray tracing required 10 5 rays reflecting no more than 10 times, but generalized CBRM required only 1.6 • 10 4 rays.The light beam and the boundaries in PS were discretized with 10 3 segments each for the generalized CBRM method.The generalized CBRM algorithm is much more accurate than MC ray tracing and requires far fewer rays to compute the intensity.We can also use generalized CBRM to find the boundaries of the positive luminance regions in T 4 ; they are shown in Fig. 13.
Algorithm 2 Recursive procedure to compute the intensity for optical systems containing curved surfaces C k,j has endpoints C min k,j and C max k,j 4: Update the path Π C t,j has endpoints C min t,j and C max t,j 11: then inverse the ordering end if INTENSITY COMPUTATION(j, C t,j , I(p), Π) 13: else 14: C s,1 has endpoints C min s,1 and C max s,1 17:

Numerical Experiments
We compare generalized CBRM to CBRM and MC ray tracing on the CPC.Recall from Section 2.2 that light in the CPC can reflect an infinite number of times, generalized CBRM can therefore run indefinitely.To prevent this we do not consider light that reflects more than 10 times.This restriction is also applied to CBRM and MC ray tracing for correct comparison.Since CBRM only handles optical systems consisting of straight surfaces, it is applied to a discretized CPC.We discretize each reflector by taking points on the parabola equidistant in the x-direction, and connect them with straight line segments.The maximum number of reflections that can occur in the discretized CPC is limited by the number of discrete segments.Recall from 2.2 that a light ray reflects on at most one reflector of the CPC.Therefore, a ray in the discretized CPC can reflect on at most all segments that discretize a reflector.MC ray tracing and generalized CBRM are applied to the regular CPC.We first compare intensity distributions of the CPC computed with generalized CBRM, CBRM and MC ray tracing to a reference solution.We apply CBRM to a discretized CPC that has 10 segments discretizing each reflector, such that the maximum number of reflections that can occur is equal to the maximum number of reflections we consider.The boundaries in PS and the light beam used by generalized CBRM are all discretized with 10 3 segments.We compute the reference solution by QMC ray tracing 10 9 rays that reflect at most 10 times.The target is divided into 110 bins for MC ray tracing.The intensity in each bin is also computed with (generalized) CBRM using Eq. ( 12) where the endpoints of the bin are p m and p m+1 .
Next, we compare the performance of generalized CBRM for various parameter settings of the algorithm.This is done twice, using 110 and 1010 intervals/bins at the target.We define the performance as the error of the solution compared to the computation time.The error is calculated with: where Ni is the number of bins (intervals) at the target and Îref denotes the reference intensity.We compare five different discretizations of the boundaries in PS using 100•2 i segments, where i ∈ {0, . . ., 4}.We compute the intensity distribution for each PS discretization 10 times using different discretizations of the light beam consisting of 100 • 2 i segments, where i ∈ {0, . . ., 9}.The reference solution is again computed by QMC ray tracing a number of rays that reflect at most 10 times.We use a reference solution of 10 9 rays when the target is divided into 110 bins and of 10 10 rays when the target is divided into 1010 bins.The best performing discretization of the boundaries in PS is the one that reaches the smallest error.If the smallest errors are similar then the performance is determined by the time it takes to compute the smallest error.In this case, the best performing discretization of the boundaries in PS is the one with the least computation time.
Finally, we compare the discretizations (of the boundaries in PS) that give the best performance of generalized CBRM for 110 and 1010 bins to the performance of CBRM.The intensity distribution of CBRM is computed using discretizations of 2 i segments on the reflectors, where i ∈ {0, . . ., 9}.The performance of CBRM is also defined as the computation time compared to the error of the solution given in Eq. ( 13).We compute the error of CBRM using the reference solutions from the previous experiment.

Results
The aim of the research was to introduce a generalization to the concatenated backward ray mapping algorithm.We discussed the original backward method and our generalized method.Both algorithms were used to compute the intensity distribution on the compound parabolic concentrator.
Fig. 14a and Fig. 14b show three intensity distributions of the CPC computed by generalized CBRM, CBRM and MC ray tracing compared to a reference solution.The figures clearly show the differences between the intensity distributions found by the algorithms.MC ray tracing computed a noisy intensity distribution with a profile similar to that of the reference solution.CBRM on the other hand found an intensity distribution with a profile that differs from the profile of reference solution.Finally, generalized CBRM found an intensity distribution closely matching the profile of the reference solution and without any noise.It took 3.6 • 10 4 rays to compute the CBRM solution, 1.6 • 10 4 rays to compute the generalized CBRM solution and 10 6 rays to compute the MC ray tracing solution.MC ray tracing gives an intensity distribution with steep sides that is similar to the reference solution because it was applied to the regular CPC.However, the intensity distribution is noisy due to the random selection of rays at the source.This is in agreement with the noisy distributions of the cup and the CPC that were previously discussed in section 3.1 (Fig. 8) and section 4.1 (Fig. 12).CBRM [1] on the other hand was applied to a discretization of the CPC where each reflector is replaced by 10 straight line segments.This optical system is different from the CPC, so it also behaves differently.As a result CBRM gave an intensity distribution without steep sides that differs from the reference solution.A finer discretization of the CPC makes the systems more alike, so the distributions should also be more alike.Generalized CBRM did not require the optical system to be discretized, and it did not use a random set of light rays.It computed an intensity distribution with steep sides that is similar to the reference solution without any noise.What stands out in these results is that generalized CBRM is more accurate than MC ray tracing and CBRM but required fewer rays to compute the intensity distribution.Fig. 15a compares the performance of generalized CBRM for different parameter settings of the algorithm using 110 intervals/bins at the target surface.Every curve was computed using a different discretization of the boundaries in PS.What is interesting about Fig. 15a that all curves have a similar shape, which implies that using a finer discretizations of the boundaries has only a small effect on the performance of generalized CBRM.The points of each curve were computed using different discretizations of the light beam.The error of all curves became much smaller as the number of light beam segments increased.We can also see in Fig. 15a that all curves stopped improving at the same discretization of the light beam.The discretization of the light beam had a larger effect on performance than the discretization of the boundaries in PS because of the recursion that happens in generalized CBRM.During every iteration, generalized CBRM computes the intersection between the light beam and the regions in PS; the algorithm continues for each intersection segment.As a result, generalized CBRM computes using a smaller subset of the discretized light beam after each reflection.Therefore, it is important to discretize the light beam with enough segments to ensure that it still closely resembles the light beam after many reflections.Still there is a point after which increasing the number of discrete segments on the light beam no longer improves the accuracy of the solution.
The test of Fig. 15a was repeated in Fig. 15b this time using 1010 bins at the target.The results are similar to those of the previous test, but there are two notable differences.The algorithm was more accurate because more bins were used at the target which allowed it to use more PS information.Furthermore, using only 100 segments to discretize the boundaries in PS significantly reduced the performance of the algorithm.This happened because generalized CBRM used more PS information to compute the intensity profile.As a result it also required a better discretization of the boundaries in PS.This shows that there is a correlation between the number of bins at the target and the minimum discretization of the boundaries in PS required by generalized CBRM.We compared generalized CBRM, using the PS discretization with the best performance, to the performance of CBRM.In Fig. 15a, all discretizations of the boundaries in PS reach a similar error value.So, the best result for the case of 110 bins at the target is the discretization with 100 boundary segments since it takes the least computation time.In Fig. 15b, the discretizations also reach a similar error except for the discretization of 100 boundary segments.So, the best result for the case of 1010 bins at the target is the discretization with 200 boundary segments.The comparison of CBRM and generalized CBRM for the case of 110 bins at the target is shown in Fig. 16a.The most striking observation in this figure is the sudden decrease of the CBRM error for a discretization of 64 reflector segments.The figure also shows that the generalized CBRM error is smaller than the CBRM error in most cases but that CBRM reaches maximum precision slightly faster than generalized CBRM.The sudden decrease of the CBRM error can be explained by the target PS of the target of the discretized CPC.The discretized CPC is different from the CPC meaning that the non-empty regions in the target PS of the target are also different from those of the applied to the compound parabolic concentrator (CPC), a standard optical system that reshapes light from a Lambertian source into a focused beam.
MC ray tracing is limited by the random selection of rays at the source because of which the intensity distribution is noisy.The number of bins at the target has a large effect on CBRM.The CPC must be discretized by many reflection segments if the number of bins is large.Otherwise, the CBRM is not accurate enough.However, if the CPC is discretized with many segments CBRM is slow.Generalized CBRM on the other hand is effected by the discretization of the light beam and the discretization of the boundaries in PS.The number of bins at the target has no effect on generalized CBRM.Of the two discretizations the discretization of the light beam has the biggest effect on performance.The number of light beam segments required for the discretization depends on the maximum number of reflections for which is the algorithm is computed; more segments are required for more reflections.Generalized CBRM computes equally good or better compared to CBRM depending on the number of bins at the target.
We have shown that there are general expressions of the boundaries in PS and that they can be used to compute the phase spaces of an optical system with curved surfaces.We introduced a generalized backward ray mapping algorithm that uses this PS information to compute the intensity distribution of an optical system.Our results showed that generalized CBRM computes the intensity distribution of the CPC faster and more accurately than the original CBRM method and Monte Carlo ray tracing.
Future research may involve exploring splines instead of straight segments to discretize the light beam and the boundaries in PS which could have an interesting effect on performance.It is also interesting to add refraction to the algorithm.This would allow us to compute the intensity distribution of many more 2D optical systems.Another research direction is to extend the method to 3D optical systems which will make it more applicable; a first step could be to find a suitable representation of the boundaries in PS for 3D optical systems.
Rays that leave the left endpoint of the source and trace out the target forming ∂S 1 Rays that trace out the source and hit the right endpoint of the target forming ∂S 2 Rays that leave the right endpoint of the source and trace out the target form-Rays that trace out the source and hit the left endpoint of the target forming ∂S4  1,4 and ∂T 4 4,1 .
Source PS of surface 1 (light source).It is partitioned into regions S 1,k where k ∈ {2, 3, 4} containing rays emitted from surface 1 and reaching surface k.Target PS of surface 4 (target).It is partitioned into regions T 4,j where j ∈ {1, 2, 3} containing rays reaching surface 4 and emitted by surface j.Source PS of surface 2 (left reflector).It is partitioned into regions S 2,k where k ∈ {3, 4} containing rays emitted from surface 2 and reaching surface k.Target PS of surface 2 (left reflector).
Source PS of surface 3 (right reflector).It is partitioned into regions S 3,k where k ∈ {2, 4} containing rays emitted from surface 3 and reaching surface k.
Target PS of surface 3 (right reflector).It is partitioned into regions T 3,j where j ∈ {1, 2} containing rays reaching surface 3 and emitted by surface j.

Fig. 6 :
Fig. 6: Source and target phase spaces of the two-faceted cup.
Source PS of surface 1 (light source).It is partitioned into regions S 1,k where k ∈ {2, 3, 4} containing rays emitted from surface 1 and reaching surface k.
Target PS of surface 4 (target).It is partitioned into regions T 4,j where j ∈ {1, 2, 3} containing rays reaching surface 4 and emitted by surface j.Source PS of surface 2 (left reflector).It is partitioned into regions S 2,k where k ∈ {2, 3, 4} containing rays emitted from surface 2 and reaching surface k.
Target PS of surface 2 (left reflector).
Source PS of surface 3 (right reflector).It is partitioned into regions S 3,k where k ∈ {2, 3, 4} containing rays emitted from surface 3 and reaching surface k.
Target PS of surface 3 (right reflector).It is partitioned into regions T 3,j where j ∈ {1, 2, 3} containing rays reaching surface 3 and emitted by surface j.

Fig. 7 :
Fig. 7: Source and target phase spaces of the CPC.

Fig. 8 :Fig. 9 : 1 } and v max 4
Fig. 8: Intensity distributions of the two-faceted cup computed with MC ray tracing (green) and CBRM (red) compared to the exact solution.

Fig. 10 :
Fig. 10: Concatenated backward ray mapping on the two-faceted cup for p = −0.2 in T 4 .The intersection segment computed at each step is colored green.
Discretization of the PS and the line p = −0.1 used by generalized CBRM.

Fig. 12 :
Fig. 12: Intensity distributions of the CPC computed with MC ray tracing (green) and generalized CBRM (red) compared to a reference solution.

Fig. 13 :
Fig. 13: Points on the boundaries of the positive luminance regions in T 4 of the CPC for all paths Π with a maximum of 10 reflections.Each color represents a different path.
Zoom of the indicated area of the intensity distributions.

Fig. 14 :
Fig. 14: Comparison of intensity distributions of the CPC.We compare MC ray tracing (green), CBRM (blue) and generalized CBRM (red) to a reference solution computed with QMC ray tracing.
CBRM performance for 110 bins at the target.
CBRM performance for 1010 bins at the target.

Fig. 15 :
Fig.15: Performance of generalized CBRM on the CPC.The error is the difference between the intensity distribution computed with generalized CBRM and a reference solution computed with QMC ray tracing.We compare PS discretizations of 100 (blue), 200 (orange), 400 (yellow), 800 (purple) and 1.6 • 10 3 (green) segments.