Comparison of Monte Carlo ray-tracing and photon-tracing methods for calculation of the impulse response on indoor wireless optical channels

: We present a comparison between the modiﬁed Monte Carlo algorithm (MMCA) and a recently proposed ray-tracing algorithm named as photon-tracing algorithm. Both methods are compared exhaustively according to error analysis and computational costs. We show that the new photon-tracing method offers a solution with a slightly greater error but requiring from considerable less computing time. Moreover, from a practical point of view, the solutions obtained with both algorithms are approximately equivalent, demonstrating the goodness of the new photon-tracing method.


Introduction
In the last years, indoor wireless optical communications have drawn the interest of researchers due to the need of high-speed and inexpensive data links, overall in environments where radio frequency links are not always viable. However, the communication capacity depends on the characteristics of the room where this communication is established [1]. The multipath penalty, which limits the maximum baud rate, can be characterized by the impulse response. In order to evaluate the impulse response on indoor wireless optical channels, several deterministic methods were firstly proposed [2]. These methods can only be implemented to determine the impulse response until the third reflection due to their computational complexity. Later on, a modified Monte Carlo-based ray-tracing algorithm (MMCA) was introduced, which presents a lower computational cost and without limit in the number of reflections to consider [3,4]. The algorithm allows easily for the introduction of new models for the devices and materials, such as the Phong's reflection model [5,6], or the increase of the complexity of the simulation environment (including obstacles and rooms with complex geometries), and it also allows to estimate the error committed by the method during the computation of the impulse response [7].
Recently, a new ray-tracing method, which has been denominated photon-tracing algorithm (PTA), has been proposed [8]. This method is based on the MMCA but now the rays are not always propagated after reflection. Instead, they are only propagated in a certain proportion according to the reflection coefficient, which makes the number of considered rays decrease rapidly after each reflection. Zhang et al. proved the reduction in computational cost of the new re-designed algorithm, but they did not comment on the accuracy of the results. In this paper, we compare in more detail the differences between the two algorithms, MMCA and PTA, according to error analysis and computational cost.
The paper is organized as follows: Section 2 describes the Monte Carlo ray-tracing algorithm, highlighting the differences between MMCA and PTA. The computational cost is evaluated in section 3. The equations that allow us to determine the error in computing the impulse response are presented in section 4. Finally, several computer-simulated results are reported in section 5 in order to compare both methods. The conclusions are summarized in section 6.

Algorithms description
In the modified Monte Carlo ray-tracing algorithm, ray directions are randomly generated according to the radiation pattern from the emitter. The contribution of each ray from the source or after a bounce to the receiver is computed deterministically. Consequently, the discretization error is due to the number of random rays. The line-of-sight (LOS) and multiple-bounce impulse responses are considered when calculating the total impulse.

LOS impulse response
Given an emitter E and receiver R in an environment without reflectors, with a large distance d between both, the received power is approximately where the emitter is modeled using a generalized Lambertian radiation pattern R E (φ , n). A eff (ϕ) represents the effective signal collection area of the receiver.
Here n is the mode number of the radiation lobe which specifies the directionality of the emitter, P E the power radiated by the emitter, A r the physical area of the receiver, and FOV the receiver field of view (semi-angle from the surface normal).

Multiple-bounce impulse response
We consider now an emitter E and receiver R in a room with reflectors. The radiation from the emitter can reach the receiver after any number of reflections (see Fig. 1). In the algorithm, many rays are generated at the emitter position with a probability distribution equal to its normalized radiation pattern R E (φ , n) /P E . The power of each generated ray is initially P E /N, where N is the number of rays used to discretize the source. In MMCA, when a ray impinges on a surface, the reflection point is converted into a new optical source, thus, a new ray is generated with a probability distribution provided by the reflection pattern of that surface, R S (θ , θ ). The process continues during the simulation time. After each reflection, the power is reduced by the reflection coefficient of the surface, and the reflected power reaching the receiver is computed by where d is the distance between the reflection point and receiver, and R s (θ , θ ) is the Phong's model, used to describe the reflection pattern of the surface [6]. This model is able to approximate the behavior of those surfaces that present a strong specular component. It considers the reflection pattern as a sum of both diffuse and specular components. In this way, surface characteristics are defined by two new parameters, the percentage of incident signal that is reflected diffusely r d and the directivity of the specular component of the reflection m. This model is described by where ρ is the surface reflection coefficient, P i represents the optical power of the incident ray, θ is the observation angle, and θ represents the incidence angle.
In PTA, when a ray impinges on a surface, the probability that a new ray is generated is given by the reflection coefficient ρ. Moreover, when a new ray is generated, its direction of propagation is determined according to the probability distribution of the reflection pattern as in MMCA, but now the power of the new ray keeps unalterable and equal to P E /N. In the same way, after each reflection, the reflected power reaching the receiver is computed by using Eq. (5), but with P i = P E /N in Eq. (6). Therefore, both algorithms perform similarly because in PTA only ρN k rays of power P ray = P E /N are reflected (being N k the number of rays remained after the previous kth bounce) whereas in MMCA all the rays are reflected but with output power P ray = ρP i . As we will see, the fact that, in PTA, not all the incident rays are reflected leads the number of computational operations to decrease rapidly with each new bounce with respect to MMCA.

Computational complexity
The main advantage of Monte Carlo ray-tracing algorithms is that they allow for the evaluation of the impulse response for rooms of complex geometries [7], in contrast to other methods [2,5], without a meaningful increase in the computational cost. This can be explained by the number of elementary calculations that are performed: N MMCA op = KNN S , where N is the number of rays, K is the number of reflections that are considered, and N S is the number of surfaces that define the room. An elementary calculation is defined as the calculation of power contribution and delay from a point source (emitter or reflection point of a ray) to the receiver, as in Eq. (5), or the assessment of the propagation of the new generated ray to determine a new point source. The previous value of N MMCA op is an upper bound, because not always, after a reflection, the N S − 1 remaining surfaces have to be considered as possible future reflecting surfaces (in rooms with complex geometries, some surfaces can be placed on other greater surfaces, e.g. windows or doors over walls; therefore, when a ray sets up from a surface with others superimposed to it, all of them can be discarded during the searching of the new collision point). Moreover, not always the ray contributes in the received power, because sometimes the emitting point is out of the FOV of the receiver and its contribution does not have to be computed. Similarly, not all the rays reach the maximum number of reflections K during the simulation time t max either. However, this value represents a good approximation to the required number of elementary calculations of MMCA.
In PTA, if we consider that N rays (photons) are initially launched, after the kth bounce, onlỹ ρ k N k−1 rays (photons) continue their path, where N k−1 is the number of photons remained after the (k − 1)-th bounce (with N 0 = N), andρ k is an average parameter of the reflection coefficient at the kth bounce which depends on the reflection coefficients of the surfaces, but also on the radiation and reflection patterns, and the position and other characteristics of the emitter. The number of elementary calculations can be computed as: It must be noticed that the last rays considered to compute their power contributions are those after the (K − 1)-th reflection. These last rays are propagated to determine the reflecting surfaces (and the reflecting points) and then their power contributions are computed. This is the reason why the previous equation has been truncated in the N K−1 term.
Evidently,ρ k depends on which bounce is being considered. However, we can assume an average reflection coefficientρ which allows us to represent any bounce obtained as an average over all theρ k . Then, the previous series is reduced to N S × N +ρN +ρ 2 N + . . . +ρ K−1 N . This numeric series converges to: As an example, we are going to consider the room studied in Section 5. The average reflection coefficient isρ = 0.69 (only taking into account the areas of the surfaces), the number of rays is N = 500 000, the number of surfaces is N S = 6 and the number of considered reflections K = 10. With MMCA, the number of elementary calculations is upper bounded by N MMCA op = 30 × 10 6 , whereas in PTA, according to Eq. (8), we will only need N PTA op = 9.3 × 10 6 elementary calculations, that is, the 31.1% of those required by MMCA. It must be noticed that the previous results are only crude approximations, but they allow us to state that, from a computational point of view, PTA is much more efficient than MMCA.

Error estimation
Both PTA and MMCA are Monte Carlo based ray-tracing algorithms. Therefore the error analysis for MMCA given in [7] is directly applicable to PTA. There, it was demonstrated that the error in computing the power P j reaching the receiver during a small time interval Δt can be estimated from the variance of P j , var P j . The biggest admissible Δt is defined as the largest interval which ensures that the same ray does not contribute twice to a receiver near the walls, being var P j given by [7]: In the previous equation, N f , j is the total number of flights during the jth time interval (rays flying during that interval), N j the number of rays that contribute in the power reaching the receiver during Δt, and p i, j is the power contribution of the ith ray (i = 1, 2, . . . , N j ) arriving during that interval obtained by using Eq. (5). Therefore, P j is the power value in the jth histogram interval computed with the Monte Carlo ray-tracing algorithm. In MMCA, N f , j coincides with N, because during a certain interval j there are exactly N rays flying, those generated by the emitter, which are never discarded. However, in PTA, after each reflection a certain number of rays are removed according to the reflection coefficient of surfaces, then N f , j decreases exponentially along the time, being different for each time interval (which has been indicated by the sub-index j in N f , j : The relative error can be computed as the square root of the variance (one standard deviation) of P j given by Eq. (9) divided by the computed contribution power defined as Then The above equation allows us to estimate the relative error in computing the received power in the jth time interval.

Simulation results
In order to compare PTA and MMCA methods, the simulation room described in [4,8] has been evaluated. The main characteristics of this room are indicated in Table 1. In this example, the reflecting surfaces are purely diffuse Lambertian (r d = 1), but rooms with materials characterized by the Phong's model (such as the examples described in [7]) have also been evaluated, obtaining similar results to those presented below. The emitter is placed at the center of the room rested on the ceiling and pointing straight down. The receiver is located on the floor near a corner pointing straight up. The maximum number of considered reflections is K = 10 and the simulation time t max = 120 ns.  In Fig. 2 we show the impulses responses (received power) obtained by using both methods when N = 500 000 rays (photons) are generated from the emitter: in Fig. 2(a) we have the response computed by the PTA method, whereas that obtained by MMCA is depicted in Fig. 2(b). We can observe how very similar impulse responses are obtained with both methods, only distinguished by a slight greater ripple in that provided by PTA.
The accuracy of both methods can be compared quantitatively by means of the relative error committed by both algorithms during the calculation of the impulse response. Figure 3 shows the relative error for PTA and MMCA calculated by using Eq. (11). It can be observed how MMCA offers a solution (impulse response) with a relative error roughly constant along the time. However, in PTA the relative error tends to increase with time, what is logical because the number of rays (photons) that contribute in the calculated received power decreases along the time, since many of them begin to be discarded when they collide against the room surfaces. This is understood clearly if we observe the number of rays that contribute in the received power along the time or at a certain bounce k (see Fig. 4). One can see how in MMCA the number of contributions always presents a relative high value (> 1/2 peak value) except at extreme time instants (t > 80%t max ) when the rays are moved apart from each other and their contributions become more and more spread. In addition, no more rays are generated after the kth reflection. On the contrary, in PTA the number of contributions decreases rapidly (in an exponential way) with time or with the bounce index. Therefore, the received power (impulse response) is computed by using a lower and lower number of information samples, leading the algorithm to present a greater relative error at longer time instants.   Table 2 we compare PTA and MMCA according to mean relative error and number of elementary calculations. The mean relative error has been weighted by the value of P j : where j = 1, 2, . . . , J and J = t max /Δt. We can see how PTA presents only an approximately 58% higher mean relative error than MMCA, since its greater values for the relative error are found at time instants where the impulse response exhibits quite low values. We have checked that by using N = 1 000 000 and N = 1 500 000 rays, PTA presents a mean relative error of +12.5% and −8.1% with respect to MMCA with N = 500 000. However, MMCA continues displaying a lower relative error for t > 50 ns. The number of elementary calculations required for both algorithms (see Table 2) demonstrates that PTA is much more efficient than MMCA (∼28% out of the total of operations required by the latter), as we had stated in section 3. We can observe how the predicted values in that section are very close to those obtained during the simulations. Therefore these expressions can be used to determine in advance the number of required calculations approximately. In comparison with MMCA using N = 500 000, the simulations performed with the PTA for N = 1 000 000 and N = 1 500 000 needed the 56% and 84%, respectively, of elementary calculations. However, for the latter we obtained a relative error inferior to that of MMCA when t < 50 ns in spite of requiring 16% less computation time.
Finally, we show in Table 2 the power distribution according to the bounce index. We can observe how both algorithms present very similar values, demonstrating a good behavior of the PTA during the distribution of power in the remained rays after each reflection.

Conclusions
In this paper, we have compared the conventional modified Monte Carlo ray-tracing algorithm with a recently proposed one, which has been called photon-tracing algorithm. We have established quantitative parameters to carry out this comparison according to computational cost and accuracy of the provided solution. We have stated analytically and by means of simulation results that PTA presents a lower computational cost than conventional MMCA. However, regarding the error committed by both algorithms, MMCA is more reliable than PTA, although the mean relative error of the latter can be considered acceptable taking into account the reduction in computing time. In addition, more rays can be used by the new method still requiring lower simulation run-time in order to obtain more accurate results. Therefore, we can conclude that PTA begins to appear as a good substitute to MMCA with superior performance.