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

Abstract

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.

1 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.

2 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.

2.1 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.

Figure 1

Shape of the two-faceted cup. Each line of the system is labeled with a number. The source $\mathsf{S}= [-2,2]$ (line 1) is located on the x-axis. The target $\mathsf{T}= [-17, 17]$ (line 4) is parallel to the source and is located at a height $z= 40$. The left and right reflectors (lines 2 and 3) connect the source with the target

The light source $\mathsf{S}= [-a, a]$ (line 1) and the target $\mathsf{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 \in \{1, \ldots , 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, \ldots , 4}$ denote the intersection of a ray with line i and, $\boldsymbol{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}\in (-\pi /2, \pi /2)$. Therefore, a ray segment between $(x_{i}, z_{i})$ and $(x_{j}, z_{j})$ with $j\neq i$ is parameterized in real space by:

$$ \boldsymbol{r}(s)= \begin{pmatrix} x_{i}-s\sin (t_{i}) \\ z_{i}+s\cos (t_{i})\end{pmatrix},\quad 0\leq s\leq s_{\mathrm{max}}, $$

(2.1)

where s denotes the arc-length and $s_{\mathrm{max}}$ is the maximum value that it can assume. A Lambertian optical source is considered; hence, the intensity at the source $\mathsf{S} = [-a, a]$ emitted in the direction $t_{1}$ is given by:

$$ I(t_{1}) = 2aL\cos (t_{1}), $$

(2.2)

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:

$$ \boldsymbol{t}_{2} = \boldsymbol{t}_{1}-2( \boldsymbol{t}_{1},\boldsymbol{\nu })\boldsymbol{\nu }, $$

(2.3)

where $(\boldsymbol{t}_{1},\boldsymbol{\nu })$ is the inner product between $\boldsymbol{t}_{1}$ and ν, ν is the unit normal to line j that the ray hits and $\boldsymbol{t}_{1}$ and $\boldsymbol{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.

2.2 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\times P$, where Q is the set of the position coordinates q and P is the set of the direction coordinates $p=n\sin {\tau }$ 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 $\tau \in [-\pi /2, \pi /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\in \{ 2, 3, 4\}$ are indicated with $(q_{\mathrm{t}, i}, p_{\mathrm{t}, i})$ on $T_{i}$. After reflection, the ray leaves line $i\in \{1, 2, 3\}$ at the same position and with a new direction, the new coordinates are indicated with $(q_{\mathrm{s}, i}, p_{\mathrm{s}, i})$ on $S_{i}$. Note that $q_{\mathrm{s}, i}= q_{\mathrm{t}, i}$ while $p_{\mathrm{s}, i}$ is obtained applying the reflection law to the direction coordinate $p_{\mathrm{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\neq i$ is the index of the line that is illuminated by i and $k\neq i$ is the index of the line that illuminates i. Hence, we indicate with $S_{i,j} \subset S_{i}$ the part of $S_{i}$ corresponding to rays that illuminate line j and with $T_{i,k}\subset 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:

$$\begin{aligned}& S_{1} = S_{1,2}\cup S_{1,3} \cup S_{1,4}, \\& S_{2} = S_{2,3} \cup S_{2,4}, \\& S_{3} = S_{3,2} \cup S_{3,4}, \\& T_{2} = T_{2,1} \cup T_{2,3}, \\& T_{3} = T_{3,1}\cup T_{3,2}, \\& T_{4} = T_{4,1}\cup T_{4,2}\cup T_{4,3}. \end{aligned}$$

(2.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 $\partial S_{i,j}$ are mapped into the boundaries $\partial T_{j,i}$ for every $i\in \{1, 2, 3\}$ and $j\in \{2, 3,4\}$ with $j\neq 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.

2.3 Computation of the boundaries of the patches in PS with positive luminance

Given two lines i and k with $i\neq 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, \ell }, z_{i, \ell })$ and $(x_{i, \mathrm{r}},z_{i, \mathrm{r}})$ the coordinates of the points located at the left and the right end point of line i, respectively. Similarly, $(x_{k, \ell }, z_{k, \ell })$ and $(x_{k, \mathrm{r}},z_{k, \mathrm{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\neq k$, $\partial S_{i,k}$ and $\partial 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 $\partial S_{i,k}$ and $\partial T_{k,i}$ are given by

$$ \begin{aligned} &\partial S_{i,k} = \partial S_{i,k}^{1}\cup \partial S_{i,k}^{2} \cup \partial S_{i,k}^{3}\cup \partial S_{i,k}^{4}, \\ &\partial T_{k,i} = \partial T_{k,i}^{1}\cup \partial T_{k,i}^{2}\cup \partial T_{k,i}^{3} \cup \partial T_{k,i}^{4}. \end{aligned} $$

(2.5)

In the following we explain in more details the case of $i=1$ and $k=4$; see Fig. 2. The boundaries $\partial S_{1,4}$ and $\partial T_{4,1}$ are given in Figs. 3 and 4, respectively. $\partial S_{ 1,4}^{1}$ and $\partial T_{4,1}^{1}$ are obtained tracing out line 4 from $q_{4, \mathrm{min}} = -b$ to $q_{4, \mathrm{max}} = b$ by rays leaving $q_{1, \mathrm{min}}= -a$ with varying $p_{1}$, these rays are shown in Fig. 2(a), and the boundary segments $\partial S_{1,4}^{1}$ and $\partial T_{4,1}^{1}$ are the orange line segments labeled with c. $\partial S_{1,4}^{2}$ and $\partial T_{4,1}^{2}$ are given tracing out line 1 from $q_{1, \mathrm{min}}= -a$ to $q_{1, \mathrm{max}}= a$ with varying $p_{1}$, such that all rays hit $q_{4, \mathrm{max}} = b$, these rays are shown in Fig. 2(b), the boundary segments $\partial S_{ 1,4}^{2}$ and $\partial T_{4,1}^{2}$ are depicted in blue (lines segments labeled with d). Likewise, $\partial S_{1,4}^{3}$ and $\partial T_{4,1}^{3}$ are obtained tracing out line 4 from $q_{4, \mathrm{max}}= b$ to $q_{4, \mathrm{min}}= -b$ by rays leaving $q_{1, \mathrm{max}}=x_{1, \mathrm{r}} = a$ with varying $p_{1}$. These rays are shown in Fig. 2(c), $\partial S_{1,4}^{3}$ and $\partial T_{4,1}^{3}$ are the red line segments labeled with e. Finally, $\partial S_{1,4}^{4}$ and $\partial T_{4,1}^{4}$ are given tracing out line 1 from $q_{1, \mathrm{max}} = a$ to $q_{1, \mathrm{min}} = -a$ with varying $p_{1}$, such that all rays hit $q_{4, \mathrm{min}} = -b$, these rays are shown in Fig. 2(d), $\partial S_{1,4}^{4}$ and $\partial T_{4,1}^{4}$ 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, \mathrm{min}}=x_{1, \ell }$, $q_{1, \mathrm{max}}=x_{1, \mathrm{r}}$, $q_{4, \mathrm{min}}=x_{4, \ell }$ and $q_{4, \mathrm{max}}=x_{4, \mathrm{r}}$.

Figure 2

Rays from source (line 1) to target (line 4). $A = (x_{1 \ell }, z_{1, \ell }) =(-2, 0)$ and $D = (x_{1, \mathrm{r}}, z_{1, \mathrm{r}}) = (2, 0)$ are the left and right corner points (or end points) of S (line 1), respectively. $B = (x_{4, \ell }, z_{4, \ell }) = (-17, 40)$ and $C = (x_{4, \mathrm{r}}, z_{4, \mathrm{r}}) = (17 , 40)$, are the left and right corner points (or end points) of T (line 4), respectively

Figure 3

Source phase space of line 1. Boundary of the region $S_{1,4}$

Figure 4

Target phase space of line 4. Boundary of the region $T_{4,1}$

For the two-faceted cup there is an analytic expression for every line segment $\partial S_{i,k}^{j}$ and $\partial T_{k,i}^{j}$ in Eq. (2.5) with $j\in \{1, \ldots , 4\}$. For instance, the rays on the boundaries $\partial S_{i,k}^{1}$ and $\partial T_{k,i}^{1}$ are parameterized in the $(x, z)$-plane by

$$ \boldsymbol{r}_{i, k}(t)= \begin{pmatrix} x_{k, \ell }-x_{i, \ell }+t( x_{k, \mathrm{r}}-x_{k, \ell }) \\ z_{k, \ell }-z_{i, \ell }+t( z_{k, \mathrm{r}}-z_{k,\ell }) \end{pmatrix},\quad 0 \leq t\leq 1. $$

(2.6)

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 $\partial S_{i,k}^{1}$ and $\partial T_{k,i}^{1}$ are

$$\begin{aligned}& \partial S_{i,k}^{1}(t)= \bigl\{ ( q_{i}, p_{i}) = \bigl(q_{ i, \mathrm{min}}, \bigl\vert \boldsymbol{\nu }_{i}\times \hat{\boldsymbol{r}}_{i, k}(t) \bigr\vert \bigr) \bigr\} , \end{aligned}$$

(2.7)

$$\begin{aligned}& \partial T_{k,i}^{1}(t)= \bigl\{ ( q_{k}, p_{k}) = \bigl(q_{ k, \mathrm{max}}-q_{i, \mathrm{min}}+t( q_{k, \mathrm{max}}-q_{k, \mathrm{min}}), \bigl\vert \boldsymbol{\nu }_{k}\times \hat{\boldsymbol{r}}_{i, k}(t) \bigr\vert \bigr) \bigr\} , \end{aligned}$$

(2.8)

where we have indicated with $\hat{\boldsymbol{r}}_{i, k}(t)$ the normalization of the ray in Eq. (2.6) and, $\boldsymbol{\nu }_{i}$ and $\boldsymbol{\nu }_{k}$ are the normalized inward normals to lines i and k, respectively. Note that, $\sin {\tau _{i}} = |\nu _{i}\times \hat{\boldsymbol{r}}_{i, k}(t)|$ and $\sin {\tau _{k}} = |\nu _{k}\times \hat{\boldsymbol{r}}_{i, k}(t)|$. Likewise, the boundaries $\partial S_{i,k}^{ j}$ and $\partial T_{k,i}^{j}$ are calculated for every $j\in \{2,3,4\}$ and $\partial S_{i,k}$ and $\partial T_{k,i}$ are found using Eq. (2.5).

In Figs. 5(a)–5(f), $(\partial S_{i,j})_{i\neq j=2, 3, 4}$ and $(\partial T_{i,k})_{i\neq 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)\rightarrow (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}$.

Figure 5

Source and target phase spaces of every optical line

2.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_{\mathrm{t}, 4}, p_{\mathrm{t}, 4})$. The intensity I along a given direction $p\in [-1,1]$ in target phase space $T_{4}$ is defined as a function of the luminance $L(q, p)$:

$$ I_{\mathrm{PS}}(p) = \int _{-b}^{b} L( q,p) \, \mathrm{d}q. $$

(2.9)

Note that the intensity is a function of $p= \sin (\tau )$ 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:

$$\begin{aligned}& L(q, p) >0\quad \forall (q, p)\in T_{4,1}, \end{aligned}$$

(2.10a)

$$\begin{aligned}& L(q, p) \geq 0\quad \forall (q, p) \in (T_{4,i})_{i=2,3}. \end{aligned}$$

(2.10b)

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 $\mathrm{M}_{1,4}$: $S_{1}\rightarrow T_{4}$ is introduced such that $\mathrm{M}_{1,4}(q_{\mathrm{s}, 1},p_{\mathrm{s}, 1})=(q, p)$. As not all parts of $T_{4}$ are illuminated by the source S, the map $\mathrm{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\neq j$, the propagation map $\mathrm{P}_{i,j}$: $S_{i,j}\mapsto T_{j,i}$ relates $S_{i,j}$ with $T_{j,i}$ and, it is defined as follows:

$$ \mathrm{P}_{i,j}(q_{\mathrm{s}, i},p_{\mathrm{s}, i})=(q_{ \mathrm{t}, j},p_{\mathrm{t}, j}), $$

(2.11)

where $q_{\mathrm{t}, j}$ is given by the x-coordinate of the intersection point between the ray and line j, and $p_{\mathrm{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 $\mathrm{R}_{j,k,h}:T_{j,k} \mapsto S_{j,h}$ relates the regions $T_{j,k}\subset T_{j}$ and $S_{j,h}\subset S_{j}$. To simplify the notation, from now on we omit the dependence of $\mathrm{R}_{j,k,h}$ from k and h, i.e., $\mathrm{R}_{j,k,h} = \mathrm{R}_{j}$. The reflection map is defined as follows:

$$ \mathrm{R}_{j}(q_{\mathrm{t},j}, p_{\mathrm{t},j})=(q_{\mathrm{s},j}, p_{\mathrm{s},j}), $$

(2.12)

where $p_{\mathrm{t}, j}$ changes according to the reflection law and $q_{\mathrm{t}, j}= q_{\mathrm{s}, j}$ as $\mathrm{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 $\mathrm{M}_{1,4}$ is described by the composition of mappings $\mathrm{P}_{i,j}$ and $\mathrm{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 $\mathrm{M}_{1,4}(\Pi )$ the map $\mathrm{M}_{1,4}$ restricted to path Π and with $R(\Pi )\subset 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(\Pi )$ 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 $\Pi =(1, 2, 4)$ and the corresponding map $\mathrm{M}_{1,4}(\Pi ): S_{1}\mapsto R(\Pi )$ that describes the propagation of all rays that follow the path Π is defined by:

$$ \mathrm{M}_{1,4}(\Pi ): S_{1,2}\mapsto T_{2,1}\mapsto S_{2,4}\mapsto T_{4,2}, $$

(2.13)

which can be written as:

$$ \mathrm{M}_{1,4}({\Pi }) = \mathrm{P}_{2,4} \circ \mathrm{R}_{2}\circ \mathrm{P}_{1,2}. $$

(2.14)

In general, to construct the map $\mathrm{M}_{1,4}(\Pi )$ 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\in [-1,1]$ we consider the intersection points between the line $p=\mathrm{const}$ and $(\partial 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 $\mathrm{P}_{i,4}^{-1}: T_{4,i}\mapsto S_{i,4}$ and $\mathrm{R}_{i}^{-1}: S_{i}\mapsto 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 $\mathrm{M}_{1,4}^{-1}(\Pi )$ of the map $\mathrm{M}_{1,4}(\Pi )$ is constructed for every possible path Π. The map $\mathrm{M}_{1,4}^{-1}(\Pi )$ is the composition of the inverses of the propagation and the reflection maps in reverse order according to the path Π. For instance, for path $\Pi = (1,2,4)$, $\mathrm{M}_{1,4}^{-1}(\Pi )$ is given by:

$$ \mathrm{M}_{1,4}^{-1}({\Pi }) = \mathrm{P}_{1,2}^{-1} \circ \mathrm{R}_{2}^{-1} \circ \mathrm{P}_{2,4}^{-1}. $$

(2.15)

The steps of the procedure are shown in the graph in Fig. 6 where the map in Eq. (2.15) is written in red.

Figure 6

Tree that describes how to detect all the possible paths from S to T

Using the procedure explained above, given a ray with coordinates $(q, p)\in T_{4}$ we can establish whether it is located inside one of the regions $R(\Pi )$ with positive luminance or not. In case the ray is inside a region $R(\Pi )$, its corresponding coordinates $(q_{\mathrm{s}, 1},p_{\mathrm{s}, 1})\in S_{1}$ are obtained using $\mathrm{M}_{1,4}^{-1}(\Pi )$, where Π is the path followed by this ray. Equation (2.10a)–(2.10b) becomes:

$$ \begin{aligned} &L(q, p) >0\quad \forall (q, p)\in R(\Pi ), \\ &L(q, p) = 0\quad \mbox{otherwise}, \end{aligned} $$

(2.16)

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(\Pi )$ and it has no contribution on the other parts of $T_{4}$. Indicating with $q^{\mathrm{min}}(\Pi ,p)$ and $q^{\mathrm{max}}(\Pi ,p)$ the minimum and maximum position coordinates of the intersection points between the boundaries $\partial R(\Pi )$ and line $p= \mathrm{const}$, Eq. (2.9) reduces to:

$$ I_{\mathrm{PS}}(p) = \sum_{\Pi } \int _{q^{\mathrm{min}}(\Pi , p)}^{q^{\mathrm{max}}(\Pi ,p)}L(q, p)\, \mathrm{d}q = \sum _{\Pi } \bigl(q^{\mathrm{max}}(\Pi ,p)-q^{\mathrm{min}}(\Pi , p) \bigr), $$

(2.17)

where the sum is over all the possible paths and the second equation holds as we assume $L=1$ in $R(\Pi )$. 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(\Pi )$ do not overlap each other, i.e.

$$ \bigcap_{\Pi }R(\Pi )= \emptyset , $$

(2.18)

where the intersection is over all the possible paths. In the next paragraph the details of the procedure to compute the coordinates $q^{\mathrm{min}}(\Pi , p)$ and $q^{\mathrm{max}}(\Pi , p)$ are explained.

2.5 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=\mathrm{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 = \mathrm{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 = \mathrm{const}$ and with the parallel rays corresponding to the end points $(q^{\mathrm{min}}, p) = (- b, p)$ and $(q^{\mathrm{max}}, p) = ( b, p)$. We set the initial intensity $I(p)=0$ along direction $p = \mathrm{const}$. Considering the intersection between the line $p=\mathrm{const}$ and the boundaries $(\partial T_{4,i})_{i=1,2,3}$ three intervals are found. Each interval corresponds to rays emitted by line i ($i \in \{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_{\mathrm{t}, i}$ and $\partial T_{i,j}$ (with $j\neq i,4$) are calculated, where $p_{\mathrm{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\neq j$ to line j are indicated with $(q_{\mathrm{t}, j}^{1}, p_{\mathrm{t}, j})$ and $(q_{\mathrm{t}, j}^{2}, p_{\mathrm{t}, j})$. The minimum and the maximum position coordinates are $q_{\mathrm{t}, j}^{\mathrm{min}} = \min \{ q_{\mathrm{t}, j}^{1}, q_{\mathrm{t}, j}^{2}\}$ and $q_{\mathrm{t}, j}^{\mathrm{max}} = \max \{ q_{\mathrm{t}, j}^{1}, q_{\mathrm{t}, j}^{2}\}$, respectively. The coordinates of the intersection points of $p= p_{\mathrm{t}, j}$ with boundaries $\partial T_{j,i}$ need to be determined for every $i \in \{1,2,3\}$ and $j \in \{2,3,4\}$ with $j\neq i$. They are indicated with $(u_{j,i}^{\mathrm{min}}, p_{ \mathrm{t}, j})$ and $(u_{j,i}^{\mathrm{max}}, p_{ \mathrm{t}, j})$ where $u_{j,i}^{\mathrm{min}}< u_{ j,i}^{\mathrm{max}}$. Since not all the rays whose corresponding coordinates are located inside the segment $[q_{\mathrm{t}, j}^{\mathrm{min}}, q_{ \mathrm{t}, j}^{\mathrm{max}}]$ with direction $p =p_{\mathrm{t}, j}$ follow the same path, the intersection segment $[v_{j,i}^{\mathrm{min}}, v_{ j,i}^{\mathrm{max}}] = [q_{\mathrm{t}, j}^{\mathrm{min}}, q_{\mathrm{t}, j}^{\mathrm{max}}] \cap [u_{j,i}^{\mathrm{min}}, u_{j,i}^{\mathrm{max}}]$ needs to be calculated. $(v_{j,i}^{\mathrm{min}},p_{ \mathrm{t}, j})$ and $(v_{j,i}^{\mathrm{max}},p_{ \mathrm{t}, j})$ are the coordinates of the rays that need to be traced back from line j to line i.

The method can be outlined as follows.

• 1. Calculate the intersection points $(u_{4,i}^{\mathrm{min}}, p)$ and $(u_{4,i}^{\mathrm{max}}, p)$ between line $p= \mathrm{const}$ and $\partial T_{ 4,i}$ for every $i\in \{1,2,3\}$, where $u_{4,i}^{\mathrm{min}}< u_{4,i}^{\mathrm{max}}$. This can be done analytically because the exact expression of the boundaries $\partial T_{4,j}$ is found as explained in Sect. 2.3.

• 2. Calculate the intersection segment

  $$ \bigl[v_{4, i}^{\mathrm{min}}, v_{4, i}^{\mathrm{max}} \bigr] = \bigl[u_{4, i}^{\mathrm{min}}, u_{4, i}^{\mathrm{max}} \bigr]\cap \bigl[q^{\mathrm{min}}, q^{\mathrm{max}}\bigr]. $$

• 3. If $i= 1$, the coordinates $(v_{4,1}^{\mathrm{min}}, p)$ and $(v_{4,1}^{\mathrm{max}}, p)$ are equal to the coordinates $(q^{\mathrm{min}}(\Pi , p), p)$ and $(q^{\mathrm{max}}(\Pi , p), p)$, respectively, of the rays located on the boundary $\partial R( \Pi )$ with $\Pi = (1,4)$. All the parallel rays with direction coordinate p and q-position coordinate with $u_{4,1}^{\mathrm{min}}\leq q \leq u_{4,1}^{\mathrm{max}}$ are emitted by the source and they directly hit the target.

  Update the intensity using Eq. (2.17)

  $$ I(p)= I(p)+q^{\mathrm{max}}(\Pi , p)-q^{\mathrm{min}}(\Pi , p). $$

• 4. If $i\neq 1$, continue with the following steps.

• 5. Trace back $(v_{4,i}^{\mathrm{min}}, p)$ and $(v_{4,i}^{\mathrm{max}}, p)$ from line 4 to line i to find their corresponding coordinates on $T_{i}$

  $$ \begin{aligned} &\bigl(q_{\mathrm{t}, i}^{1}, p_{\mathrm{t}, i}\bigr) =\mathrm{R}_{i}^{-1}\circ \mathrm{P}_{i,4}^{-1}\bigl(v_{4,i}^{\mathrm{min}}, p\bigr), \\ &\bigl(q_{\mathrm{t}, i}^{2}, p_{\mathrm{t}, i}\bigr) = \mathrm{R}_{i}^{-1}\circ \mathrm{P}_{i,4}^{-1} \bigl(v_{4,i}^{\mathrm{max}}, p\bigr). \end{aligned} $$

• 6. Update the path $\Pi = (i, 4)$.

• 7. Determine $q_{\mathrm{t}, i}^{\mathrm{min}}= \min \{ q_{\mathrm{t}, i}^{1}, q_{\mathrm{t}, i}^{2}\}$ and $q_{\mathrm{t}, i}^{\mathrm{max}}= \max \{ q_{\mathrm{t}, i}^{1}, q_{\mathrm{t}, i}^{2}\}$.

• 8. Calculate the intersection points $(u_{i,j}^{\mathrm{min}}, p)$ and $(u_{i,j}^{\mathrm{max}}, p)$ between line $p_{\mathrm{t}, i}$ and $\partial T_{ i,j}$ for every $j\in \{1,2,3\}$ with $j\neq i$.

• 9. Since not all rays whose corresponding coordinates are located inside the segment $[q_{\mathrm{t}, i}^{\mathrm{min}}, q_{ \mathrm{t}, i}^{\mathrm{max}}]$ follow the same path, compute the intersection segment

  $$ \bigl[v_{i, j}^{\mathrm{min}}, v_{ i, j}^{\mathrm{max}} \bigr] = \bigl[u_{i, j}^{\mathrm{min}}, u_{i, j}^{\mathrm{max}} \bigr]\cap \bigl[q_{\mathrm{t}, i}^{\mathrm{min}}, q_{\mathrm{t}, i}^{\mathrm{max}} \bigr]. $$

• 10. For $j\neq 1$

  • (a) Trace back $(v_{i, j}^{\mathrm{min}}, p_{ \mathrm{t}, i})$ and $(v_{i, j}^{\mathrm{max}}, p_{ \mathrm{t}, i})$ from i to j

    $$ \begin{aligned} &\bigl(q_{\mathrm{t}, j}^{1}, p_{\mathrm{t}, j}\bigr) =\mathrm{R}_{j}^{-1}\circ \mathrm{P}_{j,i}^{-1}\bigl(v_{i, j}^{\mathrm{min}}, p_{\mathrm{t}, i}\bigr), \\ &\bigl(q_{\mathrm{t}, j}^{2}, p_{\mathrm{t}, j}\bigr) = \mathrm{R}_{j}^{-1}\circ \mathrm{P}_{j,i}^{-1} \bigl(v_{i, j}^{\mathrm{max}}, p_{\mathrm{t}, i}\bigr). \end{aligned} $$

  • (b) Update the path $\Pi = (j, \Pi )$.

  • (c) Put $i=j$ and repeat the procedure from point 7.

• 11. If $j=1$, the rays reached the source and a possible path $\Pi = (1, \ldots , 4)$ is found.

  • (a) Apply

    $$ \begin{aligned} &\bigl(q_{\mathrm{s}, 1}^{1}, p_{ \mathrm{s}, 1}\bigr) = \mathrm{P}_{1,i}^{-1} \bigl(v_{i, 1}^{\mathrm{min}}, p_{\mathrm{t}, i}\bigr), \\ &\bigl(q_{\mathrm{s}, 1}^{2}, p_{ \mathrm{s}, 1}\bigr) = \mathrm{P}_{1,i}^{-1}\bigl(v_{i, 1}^{\mathrm{max}}, p_{\mathrm{t}, i}\bigr). \end{aligned} $$

  • (b) Apply the direct map $\mathrm{M}_{1,4}(\Pi )$ restricted to the path Π found:

    $$ \begin{aligned} &\bigl(q^{1}(\Pi , p), p \bigr) = \mathrm{M}_{1,4}(\Pi ) \bigl(q_{\mathrm{s}, 1}^{1}, p_{\mathrm{s}, 1}\bigr), \\ &\bigl(q^{2}(\Pi , p), p\bigr) = \mathrm{M}_{1,4}(\Pi ) \bigl(q_{\mathrm{s}, 1}^{2}, p_{\mathrm{s}, 1}\bigr). \end{aligned} $$

  • (c) Update intensity

    $$ I(p) = I(p)+q^{ \mathrm{max}}(\Pi , p)- q^{\mathrm{min}}(\Pi , p) $$

    with $q^{\mathrm{min}} = \min \{ q^{1}(\Pi , p), q^{2}(\Pi , p)\}$ and $q^{\mathrm{max}} = \max \{q^{1}(\Pi , p), q^{2}(\Pi , 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.

Figure 7

Inverse ray mapping for the two-faceted cup

The procedure starts with the rays with direction $p = -0.2$ on $T_{4}$, where $q^{\mathrm{min}}= -b$ and $q^{\mathrm{max}}= b$ are the left and the right end points of the target T, respectively. The intersection points $(u_{4,i}^{\mathrm{min}}, p)$ and $(u_{4,i}^{\mathrm{max}}, p)$ of the line $p= -0.2$ with boundaries $\partial T_{4,i}$ are computed for every $i\neq 4$.

We start from $i=1$. Therefore the coordinates $(u_{4,1}^{\mathrm{min}}, p)$ and $(u_{4,1}^{\mathrm{max}}, p)$ of the intersection points between line $p = -0.2$ and the boundary $\partial 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_{4,1}^{\mathrm{min}}, p)$ and $(u_{4,1}^{\mathrm{max}}, 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 $\partial R(\Pi _{1})$ with $\Pi _{1} = (1,4)$. Therefore, the contribution to the intensity formed by the rays that follow the path $\Pi _{1} = (1,4)$ is given by $u_{4,1}^{\mathrm{max}}-u_{4,1}^{\mathrm{min}}$.

We continue with $i = 2$. The boundary $\partial T_{4,2}$ is considered in order to find other paths. The intersection points $(u_{4,2}^{\mathrm{min}}, p)$ and $(u_{4,2}^{\mathrm{max}}, p)$ of line $p=-0.2$ with the boundary $\partial T_{4,2}$ are calculated. They are depicted in Fig. 7(b) with the magenta dots, also the intersection segment

$$ \bigl[v_{4,2}^{\mathrm{min}}, v_{4,2}^{\mathrm{max}} \bigr] = \bigl[ u_{4,2}^{\mathrm{min}}, u_{4,2}^{\mathrm{max}} \bigr]\cap \bigl[ q^{\mathrm{min}}, q^{\mathrm{max}}\bigr] $$

(2.19)

is calculated.¹ Their corresponding position coordinates $q_{\mathrm{s}, 2}^{1}$ and $q_{\mathrm{s}, 2}^{2}$ on $S_{2}$ are obtained from:

$$ \begin{aligned} &\mathrm{P}_{2,4}^{-1} \bigl(v_{4,2}^{\mathrm{min}}, p\bigr) = \bigl(q_{\mathrm{s}, 2}^{1}, p_{\mathrm{s}, 2}\bigr), \\ &\mathrm{P}_{2,4}^{-1} \bigl(v_{4,2}^{\mathrm{max}}, p\bigr) = \bigl(q_{\mathrm{s}, 2}^{2}, p_{\mathrm{s}, 2}\bigr). \end{aligned} $$

(2.20)

Note that $p_{\mathrm{s}, 2}=\mathrm{Const}$ because all the lines are straight and their normals do not depend on the location at which they are computed. Then, the corresponding direction $p_{\mathrm{t}, 2}$ on $T_{2}$ is calculated from:

$$ \begin{aligned} &\mathrm{R}_{2}^{-1} \bigl(q_{ \mathrm{s}, 2}^{1}, p_{\mathrm{s}, 2}\bigr) = \bigl(q_{\mathrm{t}, 2}^{1}, p_{\mathrm{t}, 2}\bigr), \\ &\mathrm{R}_{2}^{-1}\bigl(q_{\mathrm{s}, 2}^{2}, p_{\mathrm{s}, 2}\bigr) = \bigl( q_{\mathrm{t}, 2}^{2}, p_{\mathrm{t}, 2}\bigr). \end{aligned} $$

(2.21)

Note that $q_{\mathrm{s}, 2}^{1} = q_{\mathrm{t}, 2}^{1}$ and $q_{\mathrm{s}, 2}^{2}= q_{\mathrm{t}, 2}^{2}$ since the reflection map does not change the position coordinates. Equations (2.20) and (2.21) lead to:

$$ \begin{aligned} &\mathrm{R}_{2}^{-1}\circ \mathrm{P}_{2,4}^{-1}\bigl(v_{4,2}^{\mathrm{min}}, p\bigr) = \bigl(q_{\mathrm{t}, 2}^{1}, p_{\mathrm{t}, 2}\bigr), \\ &\mathrm{R}_{2}^{-1}\circ \mathrm{P}_{2,4}^{-1} \bigl(v_{4,2}^{\mathrm{max}}, p\bigr) = \bigl(q_{\mathrm{t}, 2}^{2}, p_{\mathrm{t}, 2}\bigr). \end{aligned} $$

(2.22)

The map $\mathrm{R}_{2}^{-1}\circ \mathrm{P}_{2,4}^{-1}$ is depicted in red in Fig. 6. The minimum $q_{\mathrm{t}, 2}^{\mathrm{min}} = \min \{ q_{\mathrm{t}, 2}^{1}, q_{\mathrm{t}, 2}^{2}\}$ and the maximum $q_{\mathrm{t}, 2}^{\mathrm{max}} = \max \{ q_{\mathrm{t}, 2}^{1}, q_{\mathrm{t}, 2}^{2}\}$ are calculated. The points with coordinates $(q_{\mathrm{t}, 2}^{\mathrm{min}}, p_{ \mathrm{t}, 2})$ and $(q_{\mathrm{t}, 2}^{\mathrm{max}}, p_{ \mathrm{t}, 2})$ are depicted in Fig. 7(c) where $p_{\mathrm{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_{\mathrm{t}, 2}=0.82$.

Next, the intersection points $(u_{2,i}^{\mathrm{min}}, p_{\mathrm{t}, 2})$ and $(u_{2,i}^{\mathrm{max}}, p_{\mathrm{t}, 2})$ of line $p_{\mathrm{t}, 2}= 0.82$ with boundaries $\partial T_{2,i}$ are computed for every $i\in \{1,3\}$. We start from the boundary $\partial T_{2,1}$ obtaining the points $(u_{2, 1}^{\mathrm{min}}, p_{\mathrm{t}, 2})$ and $(u_{2, 1}^{\mathrm{max}}, p_{\mathrm{t}, 2})$ shown in Fig. 7(c). Now, the position coordinates $v_{2, 1}^{\mathrm{min}} = \max \{q_{\mathrm{t}, 2}^{\mathrm{min}},u_{2,1}^{\mathrm{min}}\}$ and $v_{2, 1}^{\mathrm{max}} = \min \{q_{\mathrm{t}, 2}^{\mathrm{max}},u_{2, 1}^{\mathrm{max}}\}$ need to be considered. All the rays located inside the segment $[v_{2, 1}^{\mathrm{min}},v_{2, 1}^{\mathrm{max}}]$ in $T_{2}$ and with direction $p_{\mathrm{t}, 2}$ follow the path $\Pi _{2} = (1,2,4)$. In particular, the rays corresponding to the coordinates $(v_{2, 1}^{\mathrm{min}},p_{\mathrm{t}, 2} )$ and $(v_{2, 1}^{\mathrm{max}},p_{\mathrm{t}, 2} )$ are located on the boundaries of the region $R( \Pi _{2})$ on $T_{4}$ formed by all the rays that follow path $\Pi _{2}$. Their corresponding coordinates $(q^{1}(\Pi _{2}, p), p)$ and $(q^{2}(\Pi _{2}, p), p)$ on $T_{4}$ are obtained from:²

$$ \begin{aligned} &\mathrm{P}_{2,4} \circ \mathrm{R}_{2} \bigl(v_{2, 1}^{\mathrm{min}}, p_{\mathrm{t}, 2} \bigr) = \bigl(q^{1}, p\bigr), \\ &\mathrm{P}_{2,4} \circ \mathrm{R}_{2} \bigl(v_{2, 1}^{\mathrm{max}}, p_{\mathrm{t}, 2} \bigr) = \bigl(q^{2}, p\bigr). \end{aligned} $$

(2.23)

The rays corresponding to the coordinates $(q^{1}, p)$ and $(q^{2}, p)$ are located on the boundary $\partial R(\Pi _{2})$ along direction $p=-0.2$. Indicating with $q^{\mathrm{min}} = \min \{q^{1}, q^{2}\}$ and $q^{\mathrm{max}} = \max \{q^{1}, q^{2}\}$, the distance $q^{\mathrm{max}}-q^{\mathrm{min}}$ gives the contribution to the intensity $I(p)$ of the rays located in $R(\Pi _{2})$ where $p=-0.2$.

$T_{2}$ can also be illuminated by line 3, therefore the intersection points $(u_{2, 3}^{\mathrm{min}}, p_{\mathrm{t}, 2})$ and $(u_{2, 3}^{\mathrm{max}}, p_{\mathrm{t}, 2})$ of line $p_{\mathrm{t}, 2}=0.82$ and $\partial T_{2,3}$ are calculated, these points are depicted in Fig. 7(d). The coordinates $(v_{2, 3}^{\mathrm{min}}, p_{\mathrm{t}, 2})$ and $(v_{2, 3}^{\mathrm{max}}, p_{\mathrm{t}, 2})$ are shown in the same figure. As the source is not reached yet ($i=3$), the rays corresponding to $(v_{2, 3}^{\mathrm{min}}, p_{\mathrm{t}, 2})$ and $(v_{2, 3}^{\mathrm{max}}, p_{\mathrm{t}, 2})$ are followed back using the inverses of the propagation and the reflection maps. The coordinates on $T_{3}$ are shown in Fig. 7(e) with blue circles and they are obtained from:

$$ \begin{aligned} &\mathrm{R}_{3}^{-1}\circ \mathrm{P}_{3,2}^{-1}\bigl(v_{2, 3}^{\mathrm{min}}, p_{\mathrm{t}, 2}\bigr) = \bigl(q_{ \mathrm{t}, 3}^{1}, p_{\mathrm{t}, 3}\bigr), \\ &\mathrm{R}_{3}^{-1}\circ \mathrm{P}_{3,2}^{-1} \bigl(v_{2, 3}^{\mathrm{max}}, p_{\mathrm{t}, 2}\bigr) = \bigl(q_{ \mathrm{t}, 3}^{2}, p_{\mathrm{t}, 3}\bigr). \end{aligned} $$

(2.24)

The minimum and the maximum position coordinates are $q_{\mathrm{t}, 3}^{\mathrm{min}}= \min \{ q_{\mathrm{t}, 3}^{1}, q_{\mathrm{t}, 3}^{2}\}$ and $q_{\mathrm{t}, 3}^{\mathrm{max}}= \max \{ q_{\mathrm{t}, 3}^{1}, q_{\mathrm{t}, 3}^{2}\}$, respectively. Since $[q_{\mathrm{t}, 3}^{\mathrm{min}}, q_{ \mathrm{t}, 3}^{\mathrm{max}}]\subset [u_{3,2}^{\mathrm{min}},u_{3,2}^{\mathrm{max}}]$, we have that $v_{3,2}^{\mathrm{max}}\neq u_{3,2}^{\mathrm{max}}$, this means that the rays with corresponding position coordinates inside the interval $[q_{\mathrm{t}, 3}^{\mathrm{max}},u_{3,2}^{\mathrm{max}}]$ will follow a different path. The procedure continues recursively. It stops either when the rays encounter the source, i.e., when $i=1$, or when no intersection points between the direction $p=p_{\mathrm{t}, j}$ and the boundaries $\partial T_{j,i}$ are found for any $i={1,2,3}$ with $i\neq j$.

If the source is reached, then a valid path $\Pi = (1,3,2,4)$ is found. Using the inverse of the propagation map, we compute

$$ \begin{aligned} &\mathrm{P}_{1,3}^{-1} \bigl(q_{\mathrm{t}, 3}^{\mathrm{min}}, p_{\mathrm{t}, 3}\bigr) = \bigl( q_{\mathrm{s}, 1}^{1}, p_{\mathrm{s}, 1}\bigr), \\ &\mathrm{P}_{1,3}^{-1}\bigl(q_{\mathrm{t}, 3}^{\mathrm{max}}, p_{\mathrm{t}, 3}\bigr) = \bigl( q_{\mathrm{s}, 1}^{2}, p_{\mathrm{s}, 1}\bigr). \end{aligned} $$

(2.25)

The direct map $\mathrm{M}_{1,4}(\Pi ):S_{ 1}\mapsto R(\Pi )$ restricted to path $\Pi = (1,3,2,4)$, i.e.,

$$ \mathrm{M}_{1,4} = \mathrm{P}_{2,4}\circ \mathrm{R}_{2}\circ \mathrm{P}_{3,2}\circ \mathrm{R}_{3}\circ \mathrm{P}_{1,3} $$

(2.26)

is applied to the coordinates $(q_{\mathrm{s}, 1}^{1}, p_{\mathrm{s}, 1})$ and $(q_{\mathrm{s}, 1}^{2}, p_{\mathrm{s}, 1})$ giving:

$$ \begin{aligned} &\mathrm{M}_{1,4}\bigl(q_{ \mathrm{s}, 1}^{1}, p_{\mathrm{s}, 1}\bigr) = \bigl(q^{1}(\Pi , p), p\bigr), \\ &\mathrm{M}_{1,4}\bigl(q_{\mathrm{s}, 1}^{2}, p_{\mathrm{s}, 1}\bigr) = \bigl( q^{2}(\Pi , p), p\bigr). \end{aligned} $$

(2.27)

The coordinates $( q^{1}, p ) = ( q^{1}(\Pi ,p), p )$ and $( q^{2}, p ) = ( q^{2}(\Pi ,p), p )$ located on $\partial R(\Pi )$ in $T_{4}$ are found. Introducing $q^{\mathrm{min}}( \Pi ,p ) = \min \{q^{1}, q^{2}\}$ and $q^{\mathrm{max}}( \Pi ,p ) = \max \{q^{1}, q^{2}\}$, the contribution to the intensity due to the rays that follow the path Π is given by

$$ I(p) = I(p)+q^{\mathrm{max}}(\Pi , p)-q^{\mathrm{min}}(\Pi , p). $$

(2.28)

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_{2, 3}^{\mathrm{min}}, 0.82)$ and $(v_{2, 3}^{\mathrm{max}}, 0.82)$ on $T_{2}$ in Fig. 7(d). Below we explain this case in detail.

In Fig. 7(e), the coordinates $(q_{\mathrm{t}, 3}^{\mathrm{min}}, p_{ \mathrm{t}, 3})$ and $(q_{\mathrm{t}, 3}^{\mathrm{max}}, p_{ \mathrm{t}, 3})$ in $T_{3}$ with $p_{\mathrm{t}, 3}=-0.29$ are shown. They are obtained from:

$$ \begin{aligned} &\mathrm{R}_{3}^{-1}\circ \mathrm{P}_{3,2}^{-1}\bigl(v_{2, 3}^{\mathrm{min}}, 0.82\bigr) = \bigl(q_{\mathrm{t}, 3}^{1}, p_{\mathrm{t}, 3} \bigr), \\ &\mathrm{R}_{3}^{-1}\circ \mathrm{P}_{3,2}^{-1} \bigl(v_{2, 3}^{\mathrm{max}}, 0.82\bigr) = \bigl(q_{\mathrm{t}, 3}^{2}, p_{\mathrm{t}, 3}\bigr). \end{aligned} $$

(2.29)

From Fig. 7(e) we note that there are no intersection points of line $p_{\mathrm{t}, 3}=-0.29$ with $\partial T_{3,1}$. So, only the coordinates of the intersections $(u_{3, 2}^{\mathrm{min}}, -0.29)$ and $(u_{3, 2}^{\mathrm{max}}, -0.29)$ between line $p_{\mathrm{t}, 3}=-0.29$ and $\partial T_{3,2}$ are calculated. Next, the intersection interval

$$ \bigl[v_{3, 2}^{\mathrm{min}},v_{3, 2}^{\mathrm{max}} \bigr] = \bigl[ u_{3, 2}^{\mathrm{min}}, u_{3, 2}^{\mathrm{max}} \bigr] \cap \bigl[q_{\mathrm{t}, 3}^{\mathrm{min}}, q_{\mathrm{t}, 3}^{\mathrm{max}} \bigr], $$

(2.30)

formed by parallel rays with direction $p_{\mathrm{t}, 3} = -0.29$, is considered. Using:

$$ \begin{aligned} &\mathrm{R}_{2}^{-1}\circ \mathrm{P}_{2,3}^{-1}\bigl(v_{3, 2}^{\mathrm{min}},-0.29 \bigr) = \bigl(q_{\mathrm{t}, 2}^{ \mathrm{min}}, p_{\mathrm{t}, 2}\bigr), \\ &\mathrm{R}_{2}^{-1}\circ \mathrm{P}_{2,3}^{-1} \bigl(v_{3, 2}^{\mathrm{max}},-0.29\bigr) = \bigl(q_{\mathrm{t}, 2}^{ \mathrm{max}}, p_{\mathrm{t}, 2}\bigr), \end{aligned} $$

(2.31)

the corresponding coordinates $(q_{\mathrm{t}, 2}^{\mathrm{max}}, p_{ \mathrm{t}, 2})$ and $(q_{\mathrm{t}, 2}^{\mathrm{min}}, p_{ \mathrm{t}, 2})$ on $T_{2}$ are found (Fig. 7(f)) with $p_{\mathrm{t}, 2}=-0.41$. Now the procedure is repeated again for $T_{2}$ along the direction $p_{\mathrm{t}, 2}=-0.41$. No intersection points between line $p_{\mathrm{t}, 2}=-0.41$ and $\partial T_{2,1}$ occur. Only, the intersection points $(u_{2,3}^{\mathrm{min}},p_{\mathrm{t}, 2})$ and $(u_{2,3}^{\mathrm{max}},p_{\mathrm{t}, 2})$ of line $p_{\mathrm{t}, 2}=-0.41$ and $\partial T_{2,3}$ are found (see Fig. 7(f)). The intersection segment

$$ \bigl[v_{2,3}^{\mathrm{min}},v_{2,3}^{\mathrm{max}} \bigr] = \bigl[ q_{\mathrm{t}, 2}^{\mathrm{min}},q_{ \mathrm{t}, 2}^{\mathrm{max}} \bigr] \cap \bigl[q_{ \mathrm{t}, 2}^{\mathrm{min}},q_{\mathrm{t}, 2}^{\mathrm{max}} \bigr] $$

(2.32)

is calculated. The coordinates on $T_{3}$ corresponding to the end points of the intersection interval are found using:

$$ \begin{aligned} &\mathrm{R}_{3}^{-1}\circ \mathrm{P}_{3,2}^{-1} \bigl(v_{2,3}^{\mathrm{min}},p_{\mathrm{t}, 2} \bigr) = \bigl(q_{ \mathrm{t}, 3}^{\mathrm{min}},p_{\mathrm{t}, 3}\bigr), \\ &\mathrm{R}_{3}^{-1}\circ \mathrm{P}_{3,2}^{-1} \bigl(v_{2,3}^{\mathrm{max}},p_{\mathrm{t}, 2}\bigr) = \bigl(q_{ \mathrm{t}, 3}^{\mathrm{max}},p_{\mathrm{t}, 3}\bigr), \end{aligned} $$

(2.33)

where $p_{\mathrm{t}, 3} = 0.91$ (see Fig. 7(g)).

Considering the PS $T_{3}$ and the direction $p_{\mathrm{t}, 3}=0.91$, we note that there are no intersection points between line $p_{\mathrm{t}, 3}=0.91$ and both $\partial T_{3,1}$ and $\partial T_{ 3,2}$. Indeed, the whole segment $[q_{\mathrm{t}, 3}^{\mathrm{min}},q_{ \mathrm{t}, 3}^{\mathrm{max}}]$ is outside both $T_{3,2}$ and $T_{3,1}$. Because of this, all the rays with q-coordinates inside the interval $[q_{\mathrm{t}, 3}^{\mathrm{min}},q_{ \mathrm{t}, 3}^{\mathrm{max}}]$ and with direction $p= p_{\mathrm{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}\in [-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(\Pi )$ 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 $\partial R(\Pi )$ are determined for a given path Π as well as the coordinates $q^{\mathrm{max}}(\Pi ,p)$ and $q^{\mathrm{min}}(\Pi ,p)$ for every $p\in [-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$.

Algorithm 2.1

Recursive procedure for the intensity calculation

In the next section we provide the numerical results for the two-faceted cup.

3 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} < \cdots < p^{\mathrm{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\in [p^{h}, p^{h+1}]$ is given by the ratio of the number of rays $\mathrm{Nr}[p^{h}, p^{h+1}]$ that arrive at the bin $[p^{h}, p^{h+1}]$ and the total number of rays $\mathrm{Nr}[-1, 1]$ that arrive at the target, that is:

$$ \hat{I}_{\mathrm{MC}}(p)= \frac{\mathrm{Nr}[p^{h}, p^{h+1}]}{\mathrm{Nr}[-1, 1]}, $$

(3.1)

for every $p\in [p^{h}, p^{h+1}]$. Hence, the MC intensity is piecewise constant. Note that $\hat{I}_{\mathrm{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} = \frac{1}{2}(p^{ h+1}+p^{h}) )_{h = 0, \ldots , \mathrm{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⁷ rays and taking $\mathrm{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^{\mathrm{min}}(\Pi , p^{h}), p^{h})$ and $(q^{\mathrm{max}}(\Pi , p^{h}), p^{h})$ of the rays located on $\partial R(\Pi )$ are determined for every $p = (p^{h})_{h=0, \ldots , \mathrm{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:

$$ \hat{I}_{\mathrm{PS}}\bigl(p^{h+1/2}\bigr) = \frac{\int _{p^{h}}^{p^{h+1}}{I}_{\mathrm{PS}}(p)\, \mathrm{d}p}{\int _{-1}^{1}{I}_{\mathrm{PS}}(p)\, \mathrm{d}p} , $$

(3.2)

for $h = 0,1, \ldots , \mathrm{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.

Figure 8

Target phase space for the two-faceted cup divided into 100 bins. Five different paths are found. The rays with coordinates $(q^{\mathrm{min}}, p)$ and $(q^{\mathrm{max}}, p)$ in $T_{4}$ that are located at the boundaries $\partial R(\Pi )$ 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⁷ rays. The dotted green line shows the intensity found with the new method

For the two-faceted cup, the intensity can be computed analytically. This intensity is taken as the reference intensity $\hat{I}_{\mathrm{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 $\hat{I}_{\mathrm{A}}$ ($\mathrm{A} = \mathrm{MC}, \mathrm{PS}$) and the exact intensity $\hat{I}_{\mathrm{exact}} = \hat{I}_{\mathrm{ref}}$. This error is defined as:

$$ \mathrm{error} = \frac{\sum_{h = 1}^{\mathrm{Nb}} \vert \hat{I}_{\mathrm{A}}(p^{h}) - \hat{I}_{\mathrm{ref}}(p^{h}) \vert }{\mathrm{Nb}} . $$

(3.3)

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⁻⁴ is required.

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

4 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 $\mathrm{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 where all the lines are labeled with numbers. The source $\mathsf{S}= [-a, a]$ (line 1) and the target $\mathsf{T}= [-b, b]$ (line 22) are two segments both perpendicular to the optical axis, with $a=2$ and $b=17$. S is located at the height $z=0$ while T has a height $z=40$. Both sides of the system are divided into 10 segments which connect S with T. The 10 adjacent segments at the left of the system (lines $2, \ldots , 11$) connect the left extreme of the source with the left extreme of the target. Similarly, 10 adjacent segments at the right of the system (lines $12, \ldots , 21$) connect the right extreme of the source with the right extreme of the target. These segments are designed as follows. The intervals $[-b, -a]$ and $[a, b]$ are divided into 10 subintervals of the same length $(b-a)/10$. The x-coordinates of the end points of the line segments $12, \ldots , 21$ are equal to the x-coordinates of the subintervals of $[a,b]$, while the x-coordinates of the end points of the line segments $2, \ldots , 11$ are equal to the x-coordinates of the subintervals of $[-a,-b]$. The z-coordinates of every end point of the line segments $2, \ldots , 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 \in \{1, \ldots , 22\}$ are located in air, therefore the refractive index $n_{i}=1$ for every i.

Figure 11

Shape of the 20-faceted cup. The system is formed by 22 different line segments: the source S, the target T, 10 left reflectors and 10 right reflectors. $\mathsf{S}=[-2,2]$ is located at $z=0$. $\mathsf{T} = [-17,17]$ is parallel to the source and it is located at a height $z=40$. The x-coordinates of the extremes of the reflectors are equidistributed. The z-coordinates are on a parabola that passes through the end points of the target and has as symmetry axis the z-axis. All the lines are located in air

Similarly to the two-faceted cup, also for the multi-faceted cup we define the phase spaces of all the lines $i\in \{1, \ldots , \mathrm{Nl}\}$ as in Sect. 2.5, for the 20-faceted cup $\mathrm{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, $2\mathrm{Nl}-2$ phase spaces are considered. For all the systems formed by straight line segments, the boundaries $(\partial S_{i,j})_{i\neq j=2, \ldots , \mathrm{Nl}}$ and $(\partial T_{i,k})_{i\neq k=1, \ldots , \mathrm{Nl}-1}$ of the regions that form every PS can be determined as explained in Sect. 2.3.

The boundaries $(\partial T_{\mathrm{Nl},k})_{k=1, \ldots , \mathrm{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 20-faceted 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(\Pi )$ with positive luminance at target PS $T_{\mathrm{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=\mathrm{const}$ only the position coordinates $q^{\mathrm{min}}(\Pi ,p)$ and $q^{\mathrm{max}}(\Pi ,p)$ of the intersection points between the boundaries $\partial R(\Pi )$ and the line $p= \mathrm{const}$ are needed for every possible path Π. Finally, the target intensity $I_{\mathrm{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.

Figure 12

Target phase space of the 20-faceted cup. The red lines are the boundaries $(\partial T_{22,k})_{k=1, \ldots , 21}$ which are determined analytically. The numbers inside the regions $T_{22,k}$ indicate the value of the index k

5 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 $\hat{I}_{\mathrm{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 $\hat{I}_{\mathrm{ref}}$ (red line) which is computed using MC ray tracing with 10⁸ 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]).

Figure 13

Intensity for the 20-faceted cup. The red line shows the reference intensity computed using MC ray tracing with 10⁸ rays. The blue line depicts the MC intensity with 10⁷ rays. The green and dotted line shows the intensity found with the ray mapping method

In order to show the performance of the new method, we calculate the error between the approximate intensities $\hat{I}_{\mathrm{A}}$ ($\mathrm{A} = \mathrm{MC}, \mathrm{PS}$) and the reference intensity $\hat{I}_{\mathrm{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⁸ 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 $\hat{I}_{\mathrm{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. Nevertheless, we calculate the averaged intensity because we want to compare it with the MC intensity $\hat{I}_{\mathrm{MC}}$.

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⁻⁴, is desired

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⁻⁶. 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 $\hat{I}_{\mathrm{exact}}$ and the normalized approximate intensity $\hat{I}_{\mathrm{A}}$ is given by:

$$ \begin{aligned}[b] &\frac{1}{\mathrm{Nb}}\sum_{h = 1}^{\mathrm{Nb}} \bigl\vert \hat{I}_{ \mathrm{exact}}\bigl(p^{h}\bigr) - \hat{I}_{\mathrm{A}}\bigl( p^{h}\bigr) \bigr\vert \\ &\quad \leq \frac{1}{\mathrm{Nb}} \Biggl(\sum_{ h = 1}^{\mathrm{Nb}} \bigl\vert \hat{I}_{\mathrm{exact}}\bigl(p^{h}\bigr) - \hat{I}_{\mathrm{ref}}\bigl(p^{h}\bigr) \bigr\vert + \sum _{h = 1}^{\mathrm{Nb}} \bigl\vert \hat{I}_{\mathrm{ref}}\bigl( p^{h}\bigr) - \hat{I}_{\mathrm{A}} \bigl(p^{h}\bigr) \bigr\vert \Biggr). \end{aligned} $$

(5.1)

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 with the cyan dot. From numerical simulation we obtain the difference between the extrapolated value and the exact intensity

$$ \sum_{h = 1}^{\mathrm{Nb}} \bigl\vert \hat{I}_{\mathrm{exact}}\bigl( p^{h}\bigr) - \hat{I}_{\mathrm{ref}} \bigl(p^{h}\bigr) \bigr\vert / \mathrm{Nb}\approx 1.68*10^{-6}, $$

where $\mathrm{Nb}=100$ The results show that

$$ \sum_{h = 1}^{\mathrm{Nb}} \bigl\vert \hat{I}_{\mathrm{exact}}\bigl( p^{h}\bigr) - \hat{I}_{\mathrm{ref}} \bigl(p^{h}\bigr) \bigr\vert / \mathrm{Nb} \approx \sum _{h = 1}^{\mathrm{Nb}} \bigl\vert \hat{I}_{ \mathrm{ref}} \bigl(p^{h}\bigr) - \hat{I}_{\mathrm{PS}}\bigl( p^{h} \bigr) \bigr\vert /\mathrm{Nb}. $$

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.

6 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.

7 Nomenclature

$$\begin{array}{rlrl}& \mathit{r}(s)& & \text{Parametrization of a ray in the optical system}\\ & s& & \text{arc-length}\\ & \mathsf{S}& & \text{Light source}\\ & \mathsf{T}& & \text{Target}\\ & \mathrm{Nr}& & \text{Number of rays}\\ & \mathrm{Nb}& & \text{Number of bins}\\ & \mathrm{Nl}& & \text{Number of optical lines}\\ & Q& & \text{Set of the position coordinates of the rays}\\ & P& & \text{Set of the direction coordinates of the rays}\\ & S=Q\times P& & \text{Phase space of the rays with position in }Q\text{ and direction in }P\\ & S_i& & \text{Source phase space of line }i\\ & S_{i,j}& & \text{Part of }{S}_i\text{ that illuminates line }j\text{ with }i\ne j\\ & T_i& & \text{Target phase space of line }i\\ & T_{i,j}& & \text{Part of }{T}_i\text{ that is illuminated by line }j\text{ with }i\ne j\\ & n_i& & \text{Index of refraction of the medium in which the line }i\text{ is located}\\ & {\mathit{\nu }}_i& & \text{Normal to line }i\\ & q_{\mathrm{s},i}& & \text{Position coordinates of the rays on }{S}_i\\ & q_{\mathrm{t},j}& & \text{Position coordinates of the rays on }{T}_j\\ & {\tau }_i& & \text{Angle between the ray that hits line }i\text{ and }{\mathit{\nu }}_i\text{, measured}\\ & & & \text{counterclockwise}\\ & p_{\mathrm{s},i}& & \text{Direction coordinates on }{S}_i\text{, }{p}_{\mathrm{s},i}=n_{i}sin({\tau }_i)\\ & p_{\mathrm{t},i}& & \text{Direction coordinates on }{T}_i\text{, }{p}_{\mathrm{t},i}=n_{i}sin({\tau }_i)\\ & \mathrm{\Pi }& & \text{Path followed by the rays}\\ & R(\mathrm{\Pi })& & \text{Regions in target PS formed by rays that follow path }\mathrm{\Pi }\text{, }R(\mathrm{\Pi })\subset T_{\mathrm{Nl}}\\ & {\mathrm{M}}_{1,\mathrm{Nl}}& & \text{Map from }{S}_1\text{ into }{T}_{\mathrm{Nl}}\\ & {\mathrm{M}}_{1,\mathrm{Nl}}(\mathrm{\Pi })& & \text{Map from }{S}_1\text{ to }R(\mathrm{\Pi })\\ & {\mathrm{P}}_{i,j}& & \text{Propagation map from }{S}_i\text{ to }{T}_j\\ & {\mathrm{R}}_j& & \text{Reflection map from }{T}_j\text{ to }{S}_j\\ & q^{\min }(\mathrm{\Pi },p)& & \text{Minimum position coordinate on }\partial R(\mathrm{\Pi })\text{ along the direction }p\\ & q^{\max }(\mathrm{\Pi },p)& & \text{Maximum position coordinate on }\partial R(\mathrm{\Pi })\text{ along the direction }p\\ & q_{\mathrm{t},j}^{\min }(p)& & \text{Minimum position coordinate on }{T}_j\text{ with momentum }p\\ & q_{\mathrm{t},j}^{\max }(p)& & \text{Maximum position coordinate on }{T}_j\text{ with momentum }p\\ & u_{j,i}^{\min }(p)& & \text{Minimum position coordinate of the intersection points of }p=p_{\mathrm{t},j}\\ & & & \text{with }\partial T_{j,i}\\ & u_{j,i}^{\max }(p)& & \text{Maximum position coordinate of the intersection points of }p=p_{\mathrm{t},j}\\ & & & \text{with }\partial T_{j,i}\\ & v_{j,i}^{\min }& & max\{q_{\mathrm{t},j}^{\min },u_{j,i}^{\min }\}\\ & v_{j,i}^{\max }& & min\{q_{\mathrm{t},j}^{\max },u_{j,i}^{\max }\}\\ & L(q,p)& & \text{Luminance at the target}\\ & {\hat{I}}_{\mathrm{MC}}(p)& & \text{Monte Carlo normalized intensity at the target}\\ & {\hat{I}}_{\mathrm{PS}}(p)& & \text{Phase space averaged and normalized intensity at the target}\\ & {\hat{I}}_{\mathrm{ref}}(p)& & \text{Reference normalized intensity at the target}\end{array}$$