Optimizing Ray Tracing Techniques for Generating Large-Scale 3D Radio Frequency Maps

Accurate estimation of the signal power received at a given location can be obtained through Ray Tracing (RT) propagation model. Nevertheless, processes involved in RT are extremely complex and not adapted for Radio Frequency (RF) mapping purposes. Particularly, the reception test process, responsible for validating received rays at given locations, contributes to over 98% of the overall RT complexity in large-scale scenarios. In this paper, we introduce a novel RT algorithm able to divide by a factor of 50 the RF map rendering time while limiting the difference in bitrate estimation to less than 1Mbps on average as compared to existing approaches. Our algorithm integrates hardware acceleration techniques for RT, proposes a new modeling of the terrain in the form of a triangular mesh, and introduces a new reception test able to efficiently calculate the signal power received per each triangle of the mesh. Our simulations show how these processes together can reduce the complexity of RF map generation, hence allowing the algorithm to scale to large-scale and precise scenarios.


I. INTRODUCTION
Generating Radio Frequency (RF) maps at large-scale with the objective to give the signal strength level and bitrate capacity estimation everywhere in a city or a country is gaining a lot of momentum nowadays [1], [2].These maps can help end users to better understand the Quality of Service (QoS) offered by their Mobile Network Operators (MNOs) depending on their living place.They can also help MNOs themselves to troubleshoot their own network and improve their QoS when needed, by for example appropriately managing the antennas' placement and configuration [1].These RF maps are also useful to telecommunication regulatory authorities to check if mobile network operators' specifications are met [2].
Nevertheless, generating RF maps is quite challenging given the large complexity of the environment (topography, buildings, obstacles, etc.).Indeed, to meet the above objectives and be useful, RF maps should cover a large area and reach a sufficient level of accuracy to be able to estimate the signal strength across different locations in a city or even a country.Due to this complexity, telecommunication regulators and most MNOs provide end-users with coverage maps illustrating areas covered by MNO's antennas [3].These maps often use circles to indicate the extent of the coverage area.While existing maps are valuable for identifying areas lacking coverage, they do not provide precise and accurate information about the bitrate received at those specific locations.
To mitigate this challenge, Ray Tracing (RT) is a wellknown propagation model known to be accurate for calculating the paths followed by the radio waves and for estimating the signal power received in different locations of the area of interest.Indeed, in RT, different multi-path propagation events such as reflection, refraction and diffraction [4], are captured in a deterministic way across the different obstacles of the environment.This level of detail contributes to improving the quality of the produced RF maps.However, this comes at the cost of high computational load and long execution time due to the complexity of the processes involved in it.For example, the icosahedron technique used to generate rays in RT is proven to be computationally heavy [5] [6].Also, a naive ray/object intersection test leads to a high execution time in complex scenarios [7] [8].Further, the reception test that consists of checking whether a ray is received or not by a set of receivers also suffers from a heavy matrix computation and a greedy memory consumption [9].
Despite its appealing features, RT was not primarily designed for RF mapping purposes [9].As any propagation model, it is designed to compute the path loss between a given transmitter Tx and a given receiver Rx.The intuitive way to generate RF maps using RT is to discretize the area of interest using a limited number of reception points as outlined in previous studies [10] [11] [12].These reception points are distributed throughout the area according to some statistical distribution to ensure comprehensive coverage.The traditional RT algorithm is then performed for each pair of transmitting antenna and reception point.This approach is used for example in [10] and [11] to generate a multi-cell coverage map in the case of millimeter-wave frequency band and smallcell networks.Although intuitive, this approach has several drawbacks.First, the map produced is not always precise as its precision is related to the number of reception points.To the best of our knowledge, no work exists about the appropriate distribution of reception points that gives the best precision over the entire map.Second, depending on the number of reception points, some rays may simply not be considered, since only rays that are received by at least one reception point are considered.This can lead to serious accuracy and precision problems when some areas are considered not to receive rays, while in reality the issue stems from not having placed enough reception points.To solve this problem, a large number of reception points must be used to ensure capturing all the rays.However, this large number of points leads to excessive computational load and potential memory limitations, making it challenging to generate the RF maps effectively.
To overcome these problems of RT complexity, different sets of acceleration techniques were developed [13] [14] [15].However, none of these techniques is designed for RF mapping purposes.The visibility algorithm [16] for instance was developed to reduce the overhead incurred by the ray/object intersection test.It is meant to reduce the number of candidate tests between Tx and Rx.Indeed, a naive intersection test consists of testing for each pair of ray and building whether there is an intersection or not.By reducing the number of candidate tests, the visibility algorithm reduces the number of intersection tests performed, hence decreasing the complexity of RT.Nevertheless, when it gets to RF mapping, many receivers are involved, making this algorithm suffer from its Tx/Rx dependency [8].Indeed, the algorithm must be repeated N × M times, with N being the number of transmitting antennas and M the number of reception points.This results in many unnecessary computations performed without leveraging the inherent redundancy, thus having a high execution time.
In this paper, we present a novel ray tracing algorithm able to accurately and rapidly generate RF maps at largescale whilst taking into account the 3D topography of the area of interest.We introduce new schemes based on 3D triangulation to model the area and we leverage efficient hardware acceleration techniques to generate RF maps at large-scale.Instead of calculating the received signal power over a finite set of reception points and inferring the map by interpolation, we propose a new reception test able to directly calculate the received signal power over the triangles composing the area, which allows for the consideration of all the rays hitting the surface and the direct rendering of the 3D map.Our contributions can be summarized as follows: • Introduction of a new acceleration scheme efficient for RF mapping by adapting Embree [17], an Intel ray tracing engine used in computer graphics.• Introduction of a surface topography modeling based on 3D triangulation to account for all variations in the surface's elevation and all obstacles composing the area.• Design of a new reception test between 3D triangles and rays based on their geometries.• Proposition of a new RF map rendering approach based on the signal power received by the triangular mesh.
Simulations of our RF mapping solution show that our approach outperforms existing ones based on discrete reception points both in terms of execution time and memory consumption.Indeed, our approach is 50 times more efficient than traditional RT approaches to RF mapping as the one implemented in the Matlab official distribution [18].Due to memory limitations, existing solutions are unable to produce the RF map in urban areas involving multiple antennas and requiring a precise map.On the other hand, our solution generates a precise and accurate RF map in all types of terrains with a difference as small as 1Mbps on average compared to existing solutions.We also show that our approach preserves the accuracy of the RF map when changing the resolution step, i.e., with a 50 meters resolution step we achieve almost the same accuracy as when the resolution step gets near zero.Furthermore, we demonstrate in this paper that the reception test accounts for more than 98% of the complexity of traditional RT approaches in urban complex areas.We then highlight how our reception test is able to tackle this issue with only a slight increase of the execution time with regards to the resolution step.We equally show how the ray/object intersection test in our solution which is based on Embree behaves independently from the number of triangles, i.e., one can choose a fine-grained size for triangles to improve the RF map precision without degrading the performance of the intersection test.
The rest of the paper is organized as follows.In Section II, we overview the state-of-the-art for related works on the generation of RF maps.Section III is an overview of our solution with a main focus on its acceleration part.In Section IV we present our RF map generation method with its novel reception test process, while simulation results and discussions are shown in Section V. Conclusions are given in Section VI.

II. RELATED WORKS
To mitigate Ray Tracing (RT) complexity in complex environments, different sets of techniques were proposed.The visibility algorithm explained in Section I suffers from its Tx/Rx dependency, as it is optimized for each pair of transmitter and receiver.Another technique based on triangular grids approach [8] was proposed to overcome this dependency.In the latter, the vertices of the buildings serve as point clouds that are used to generate a triangular mesh.The obtained triangles are used afterwards to track the rays and check in which triangles the intersections with rays happen.Proceeding this way, the intersection test does not depend anymore on Tx and Rx but on the number of buildings itself, hence giving better results and accelerating the ray/object intersection test for large-scale scenarios.However, the technique in [8] has several drawbacks.First, it is designed for a 2D scenario with a heuristic proposed to generalize it to 3D scenarios.According to this heuristic, when rays are received by a receiver in 2D, a reverse work is performed to check whether these rays hit buildings or not, taking into consideration their respective heights.This approach provides a rough approximation of rays received.In fact, in 3D rays are usually launched by an antenna placed at the center of a unit sphere in such a way to cover all azimuths and elevations; while in 2D, rays are launched in all azimuths of a single chosen elevation.It follows that all the other elevations are not accounted for.Second, due to this 2D nature of the technique, propagation mechanisms from the ground are not considered.Hilly or mountainous terrains are approximated as being flat without their elevation being considered.Finally, the number of triangles used to speed up the intersection test is a function of the number of buildings' vertices.In complex environments, triangles' size may be close to zero, leading to inefficiencies and increased complexity of the ray/object intersection process.
The second aspect that is widely covered in the literature is the sampling of the area of interest to produce the RF maps as discussed in Section I. Sampling is used in [19] to analyze the propagation characteristics of drone transmission at different frequencies.The simulations of the drone were performed through RT and the RF map was displayed on a 450m x 500m terrain with a 5m resolution step between the reference reception points.In [20], it is used to show the RT computation of the received power distribution over different antenna frequency settings in an indoor scenario with 10cm resolution step.Finally, in [21], the RF map is constructed for the received signal strength computed via RT at predefined reference points for indoor positioning purposes.
Sampling at high rates comes with an inherent complexity caused by the reception test, which is needed to determine if rays have been received or not by the reception points.Indeed, for each pair of ray and reception point, one must check whether the ray is received or not.In scenarios requiring high precision as in [19] and [20], large matrix computations need to be performed to determine which rays have been received by which reception point.Moreover, in large-scale scenarios, namely at city-level, millions of rays are launched and thousands or even millions of reception points are considered depending on the level of precision sought.Along with the large number of reflections that may be needed in some scenarios, the reception test can quickly get out of control [9].One must give up on the accuracy and/or the precision to be able to produce the RF map.
In [22], the discretization approach is used to show the RF map over a university campus.Authors designed a new technique to optimize the reception test and the field computation process.They built a triangular mesh of their area of interest using high resolution LIDAR data and ran 3D RT on it.Their approach requires for each pair of Tx and Rx an exhaustive search amongst all the triangular surfaces to determine triangles involved in the propagation mechanisms (reflection and diffraction) between the transmitter and the receiver.In complex scenarios including large number of reception points and transmitters, this approach leads to an impractical execution time.To enhance their model's computational efficiency, a subset of surfaces is a priori identified for each pair of Tx/Rx based on the distance between them.This optimization, although useful in reducing the computational complexity of their model, compromises its accuracy.Indeed, in case of high order reflections and diffraction, many triangles may be missed hence altering the real signal power received at a given geographical location.
To the best of our knowledge, none of the existing RT techniques are designed for RF mapping purposes and/or are not sufficiently optimized to meet the requirements of largescale RF map generation with account for the topography of the surface.To mitigate this, we propose in this paper a fast RT solution fully designed to generate RF maps at large-scale.The acceleration technique leveraged in this paper as well as the cartography solution that we propose are both tuned to meet the requirements of generating RF maps at largescale while considering the 3D topography of the terrains and the information on the buildings.For this, we speed up the ray/object intersection test by adapting a fully optimized RT engine called Embree [17].As shown later, Embree's complexity is only ray-dependent, i.e., the cost of the ray/object intersection test remains almost the same while changing the number of triangles (aka the resolution of the map).This property helps us to generate high quality RF maps with a large number of triangles without increasing the complexity of the ray/object intersection process.Moreover, we leverage the geometrical properties of the triangular mesh used during the ray/object intersection test to propose a new reception test scheme that is suited to RF maps generation and does not require the consideration of reception points.Our scheme consists of directly performing a spatial intersection between the triangles and the circular cross-sections of the rays, instead of going through the complex matrix computations in the current reception test approach.Thanks to these underlying RT processes being introduced, we can move directly from the triangular mesh of the terrain to the RF map production after proposing signal power estimators for triangles and coloring them accordingly.This forms an optimized pipeline fully designed for RF mapping purposes and capable of reaching the goal of RF map generation at large-scale, regardless of the type of the terrain.Our simulation results show that this optimized pipeline can divide by a factor of 50 the RF map generation time as opposed to a traditional RT approach.

A. Real-life application
Our technique is meant to generate RF maps in real-life scenarios.All our work is made on real-life open data sets obtained from governmental open-source databases in France.We collected antennas information from the French National Frequency Agency (ANFR) [23].ANFR's antennas data set is publicly available, updated monthly and contains information about all the antennas in France.Buildings' data sets were obtained from the French National Institute of Geographic and Forest Information (IGN) [24].

B. Topography and buildings modeling
To consider the variations of the 3D topography of the terrains and to account for the different buildings and obstacles over it, we proceed by triangulation.Triangulating the terrain consists of uniformly dividing the terrain into triangles of the same size determined by the resolution step.With the latter, we represent a 3D terrain using a uniform mesh of triangles able to accurately account for altitude variations in the terrain (hills, valleys, etc).Based on these 3D triangles, we can comprehensively capture and model fine-grained changes in the terrain needed for our study.This setup is very interesting, since we will be able to generate precise RF maps showing the signal strength and interference without any altitude approximation.Additionally, we model buildings in the same fashion as the topography using triangulation.We carry out buildings' triangulation in two folds: 1) Since the facets of a building are vertical, we simply divide each facet into two triangles.With this division, we could capture all the information on the building's facets without any approximation.The vertices of the obtained triangles have 3D coordinates each representing the longitude, the latitude, and the elevation of the vertex.2) Buildings' rooftops are horizontal with a polygonal shape.The form and size of the polygon is not known a priori as it depends on each building.There may be as many rooftop shapes as the number of buildings in the area of interest.We then simply perform a Delaunay Triangulation [25], a well-known triangulation technique, on the rooftops' vertices to obtain their triangular mesh.

C. Ray/Object intersection test
To the best of our knowledge, none of the acceleration techniques used in the literature were designed to account for the variations of the ground when performing the ray/object intersection test.Most of these techniques focus on optimizing the ray/building intersection test without any regard to ground reflections.In reality, almost all the rays launched below antennas' horizon hit the ground and are reflected back in the area of interest, hence the importance to take this reflected energy into account when performing RT at large-scale.On the other hand, the large number of ray reflections over the buildings' facets and rooftops is also very important mainly in urban areas.There is thus a need for a technique fast enough to account for intersections due to both the ground and the buildings.This is why we chose Embree, a high-performance RT library developed by Intel for x86 Central Processing Units (CPUs).Embree is primarily designed to accelerate RT computations in 3D rendering applications.The low-level kernel of Embree is designed to leverage the computational power of modern CPUs to optimize RT operations.In particular, the SIMD (Single Instruction, Multiple Data) feature of Embree allows parallel processing of multiple data elements in a single instruction and leads to significant performance gains.Embree also uses Ray Sorting, which consists of rearranging rays based on their proximity to one another, ensuring that nearby rays are processed together.It also provides a collection of functions and data structures that enable efficient ray traversal and intersection calculations.Furthermore, the Bounding Volume Hierarchy (BVH) technique used by Embree allows to organize the geometries in the area of interest into a hierarchical tree structure.Ray/Objects intersections are then performed by traversing the tree, hence reducing further the computational overhead incurred by RT.
We build upon Embree and develop a tool for a fast ray/object intersection test.When the rays and the geometries in the area of interest are given to Embree as input, we get as output the intersection points between the rays and the objects.
To identify these intersection points, Embree takes as input a triangular mesh of the objects (in our case, the ground, and the buildings), together with the origins and directions of rays.This property of Embree is of great interest to us because it is fully compliant with the terrain triangulation scheme explained in Section III-B.Since we model the terrain surface by 3D triangles, our model can be automatically considered by Embree when performing the intersection tests.Once the intersection points between rays and objects are found, our tool calculates the reflected rays and proceeds in a recursive fashion until the maximum number of reflections is reached.

IV. RF MAP GENERATION
Above we presented an overview of our approach with a focus on the models it embeds and the acceleration technique it uses to speed up RT.Since our final goal is to produce large-scale RF maps, we focus here on the rendering part.
As said in Section II, the state-of-the-art approach to produce RF maps has performance, precision and accuracy issues.To overcome these issues, we leverage the triangular mesh obtained from the terrain modeling scheme introduced in Section III-B to directly produce the RF maps.The signal power corresponding to each triangle is computed based on the rays it receives.With this power computed for every triangle, we can produce the RF map as a heat map of triangles, thus giving the level of signal in every single location of the area of interest at the granularity of considered triangles.This approach for rendering has the following advantages: • Completeness: All the rays hitting the ground and buildings (i.e., the triangular mesh) will be considered.None of them will be discarded as opposed to the traditional approach where some rays are not considered because they are not received by the chosen reception points.• Precision: Since the triangles cover the entire 3D space, we can capture the power received at every single location without the need for further interpolation.Whatever is the zoom level of the RF map, our approach can always provide an estimate of the signal power received at any location in the area, which is not the case of the state-ofthe-art approach where power is calculated over a set of reception points.We recall here that the precision of our approach is still determined by the size of triangles we consider.Later, we will show sensitivity results on the impact of this size on the quality of the RF map.• Applicability in real life: The triangles used are in 3D and representative of the terrain's topography.Since we are leveraging those triangles to compute the signal power received, we can produce accurate RF maps without any altitude approximation.• Scalability: Our approach is applicable to any city or country regardless of its size.Triangles being independent from each other, we can treat them separately by splitting a big area of interest into smaller ones and merging them afterwards to produce the final RF map.To achieve all these goals, a new reception test is needed to account for the geometrical interaction between the triangles and the rays.We then build upon all the rays received by triangles to calculate both the power of the main signal and the interference for each triangle and use them to color the triangles and produce the RF maps.

A. Reception test
The reception test used in the literature to verify whether rays are received by the reception points is a cumbersome process.The large matrix computations performed in this process lead to high memory consumption that prevents largescale applicability of RT.It also leads to a high execution time, thus preventing RT from producing maps in many scenarios.Since our goal is to produce large-scale RF maps in a reasonable time, we remove this cumbersome process and propose a new one more suitable for our use case.Our proposed test consists of leveraging the geometrical properties of the rays and the triangles they intersect with and performing spatial intersection between them.Indeed, in reality rays are cone-shaped and the radius of a ray cross-section is given by R = αd √ 3 [26], with α being the maximum angular separation of the ray with its neighbors and d is its unfolded distance from the antenna to its reception by the triangle.Given this, we start by computing the radius of all the rays hitting the ground at their points of intersection.Since the cross-sections of the rays are circular, this corresponds to directly performing a spatial intersection between the two geometrical figures: circles and triangles.Instead of going through all the complex matrix computations of the traditional RT reception test, we propose the following test of reception for a certain ray received by a certain triangle: a ray is said to be received by a triangle if its cross-section overlaps with that triangle.Leveraging the geometrical properties of rays and 3D triangles defined by their geographical locations and shapes, we avoid performing the reception test of each pair of ray and Rx as is the case in the literature.Instead, the reception test now consists of an intelligent spatial intersection operation that accounts for the locations of rays and triangles.Due to these, we are even able to perform these operations at once, by just joining the two geometries (rays and triangles).As result we obtain which triangle overlaps with which ray in a rapid and efficient fashion.In this paper, we use Shapely [27], a geometrical tool to perform this operation.
Furthermore, our new test can easily account for the fact that a ray can be received by more than one triangle.Indeed, due to the new spatial intersection test that does not only consider the center of the ray but its whole geometry, the contribution of a ray is considered in all its overlapping triangles.This property is of great interest since all the geometries of rays are accounted for, which leads to the production of accurate RF maps.This final information on the rays that are received by given triangles is needed for the computation of the downlink signal strength and the level of interference, which can then be used to calculate the Signal to Noise plus Interference Ratio (SINR) and so to estimate the channel capacity (bitrate) as well as to produce the RF map for either of these QoS metrics.

B. Power computation
Once the rays received by the different triangles are known, we compute the SINR for each of them.The intuition behind computing the power received by a triangle is to assume each triangle to be the equivalent of a reception point.The only difference is that a point corresponds to a given location and that a triangle can cover several locations.To calculate the power across a triangle, we apply a similar methodology to the one used to calculate the power at a reception point while getting rid of the phase shift between received rays.Indeed, we consider the phase shift between rays to be a uniform random variable caused by the small value of the wavelength at high frequencies being far less than triangles' size.By considering the phase of a ray as a uniform random variable in the range [0; 2π] independent of the phase shifts of the other rays, we can easily derive an expression for the expectation of the received power.This expectation corresponds to the average power received by a triangle.
Consider a set of received signals in the form s i = A i e j(wt+φi) , with A i , w and φ i being respectively the amplitude, the frequency and the phase of the i th signal.All these signals are supposed to be in the same frequency band.The power of the total received signal can be expressed as: Hence, It follows from Equation (1) that the average power received by a triangle is the sum of the powers of all the rays it receives for the same frequency band.This power P rx is given by Equation ( 2) with P i being the received power of the i th ray, P tx is the transmitted power, G tx and G rx are the antenna gains, λ is the wavelength, d is the unfolded distance crossed by the ray and R n is the n th reflection coefficient of the i th ray at the n th obstacle along its path [6].
This power expression helps to derive the SINR, which can be computed using Equation (3).In this latter equation, the signal source (S) is the antenna with the highest signal power, whereas other antennas at the same frequency act as sources  of interference (I), and N denotes the total noise power at the receiver side.
Once the SINR of all the triangles is computed, the RF map can be generated in a heat map fashion showing the SINR at a every single location in the area of interest.Note that with the SINR in hand, one can compute the bitrate depending on the information about the technologies used by the antennas: 3G, 4G or 5G.From the SINR, the downlink Channel Quality Indicator (CQI) can be derived for instance, then mapped into the download bitrate at the given location using the specifications related to each wireless technology.For our simulation purposes, we opt for the Shannon capacity which gives the maximum channel capacity given the SINR value and the channel frequency bandwidth.

V. RESULTS
To compare our solution against existing ones, we implemented it in a Ray Tracer that we built from scratch using the Python programming language.This Ray Tracer includes all the aforementioned optimizations brought by our solution.Using our tool, we conducted several simulations to evaluate our solution's performance and assess its capacity to generate RF maps at large-scale while considering the terrain's topography.These simulations have three main objectives: (i) validate our solution against a state-of-the-art approach, (ii) evaluate its execution time as compared to existing RT tools and, (iii) show the effect of the triangles' sizes on the overall accuracy of the obtained map.We conducted the simulations over three different real terrains, all having a surface area around 79km 2 .Characteristics of these three terrains are given in Table I.From the number of antennas and buildings, we can see that the two first ones are urban areas, and the 3rd terrain is a rural area.Figure 1 shows sample of buildings of Terrain 1.We only show 20000 buildings in this figure for better visualization.The figure is not completely circular as its down part corresponds to the sea, with no buildings.For our simulations, we use a Linux server featuring 32 processors and 94 Gigabytes of memory.General parameters values used in our simulations are given in Table II.

A. Validation
We first start by validating our model and checking the correctness of the new concepts involved in its processes.We use for this a state-of-the-art traditional implementation of Ray Tracing (RT), namely the RT tool developed by Matlab [18].Indeed, Matlab offers a toolbox that helps to generate RF maps using the space sampling approach explained in Section II.
According to Matlab, one can discretize the area of interest following a uniform distribution with a given resolution step between the reception points.Since our approach is trianglebased, the received signal power of a given location is the signal power of the triangle to which it belongs, contrary to Matlab which directly considers discrete reception points.To be as close as possible to Matlab's estimation, we set the resolution step in our model to the smallest value that we could take constrained by memory limitations, which is 2 meters.Indeed, the smaller the size of the triangle is, the closer we are to the power received by a single reception point.Based on this, we run Matlab's RT on 24804 reception points and obtain their received signal power.Then, we check to which triangles these points belong in our model and assign to them the signal power of these triangles.This results in signal power estimations for both Matlab and our model for the same locations.
To assess the robustness of our approach, we perform the simulations over different settings.For all the three terrains, we vary the number of reflections and the angular distance between rays and compare our model's estimation to the one of Matlab.We choose three different ray angles known as angular separation values corresponding to Matlab's predetermined values: large: 1.07°, medium: 0.55°, small: 0.28°.It is important to note that small angular separation means that the distance between neighboring rays is small due to the large number of rays launched.On the other hand, three different values of maximum number of possible reflections are chosen: 0, 2 and 5.The former corresponds to Line-Of-Sight (LOS) scenarios while the latter are Non-LOS (NLOS) ones.
Figures 2 and 4 show the cumulative distribution functions (CDFs) of our model's bitrate estimation against the one of Matlab in respectively LOS scenario and with 5 maximum number of reflections.Figures 3 and 5 show how these  estimations are correlated together in these two scenarios.In the LOS case, the two distributions overlap exactly with each other, i.e., the two estimations have the same distributions, with their bitrate values ranging from 15 to 32Mbps.The percentiles of the bitrate are the same for the two models, since 25% of Matlab's and our model's estimations have a bitrate lower than 21Mbps, the two models have the same median of 22.5Mbps and so on.This result shows that our model completely aligns with the one of Matlab even though our model uses a different approach to generate RF maps.Furthermore, the Root Mean Square Error (RMSE) which gives the mean of the point-to-point deviation between our model's estimation and Matlab's is 0.63Mbps.This value is explained by Figure 3, where we see a linear correlation of almost 96% between the two models.We can observe from the figure that the slope of the direct line is 1, i.e., that in almost 96% of cases, our model and Matlab estimations are the same.This shows the ability of our model to accurately estimate the signal power received at a given location.The figure also shows the outliers on both sides of the main line.This slight difference is because in LOS scenarios, Matlab takes the exact path between Tx and Rx, while in our model we consider the rays intersecting the triangle in which Rx lies in.This results in a slight difference in the distance followed by the ray, and therefore a difference in the power estimate, where the power is a function of 1/d 2 .These outliers explain this RMSE value, which sounds reasonable in view of the range of the bitrate.
Similarly in Figure 4, we see that in the case of 5 maximum number of reflections, the two models still maintain the same distribution, i.e., they have the same overall statistics in terms of bitrate estimation as in the LOS case.This shows that the maximum number of reflections does not affect the accuracy of our model and that it is able to accurately generate RF maps including high order of reflections despite the new schemes introduced.Indeed, the RMSE of 1.09Mbps remains acceptable in this context.Figure 5 provides the explanation for this small error.From this Figure, we can see that in almost 89% of cases, our model gives the same estimation as Matlab.The other 11% correspond to reception points where the two models estimate the bitrate differently.These latter points seem to be distributed in a symmetric way around the diagonal.We interpret this phenomenon by a difference in the set of rays received between the two models on these pointsIndeed, as the number of reflections increases, the likelihood of receiving the same number of rays for our model and Matlab is reduced for some locations.This phenomenon leads to slight outliers on both sides of the diagonal.
Overall, our model maintains the accuracy of RT by minimizing the differences with traditional RF mapping ap- We can observe that the larger the number of rays launched is, the lower is the error.This is because when the angular separation between rays becomes small, the rays' cross-sections decrease as explained in Section IV-A, and so the less will be the difference between a reception point and a reception triangle.The same thing happens with the number of reflections; when the number of reflections gets higher, we notice a slight increase in the error, as more reflections means longer distances for rays and larger rays' cross-sections.Nevertheless, we can see from the table that despite the terrain, the error is almost constant.This shows the robustness of our approach and its applicability over different terrains without altering the error.The case of Terrain 3 is interesting, as the error does not seem to change with the number of reflections.Since this terrain is a rural area, a large number of reflections is not necessary, due to the limited number of buildings.Overall, our model maintains a very good accuracy with an average RMSE that is less than 1Mbps, which is a small error regarding the bitrate values that is between 0 and 35Mbps in the scenarios we considered.

B. Execution time
In this section, we conduct a comparative analysis of the execution time between our RF mapping solution and the traditional approach presented in Section II.Since Matlab allows the generation of RF maps based on this approach, we compare our model against it.To ensure fairness along this comparison, we run the two solutions with comparable configurations.Indeed, in Matlab the terrain is discretized using a uniform distribution with a specific distance between the reception points.In our model, this distance corresponds to the size of the square used for triangulating the environment.To ensure that the number of reception points in Matlab is equivalent to the number of triangles in our model, we set the two resolution steps in Matlab and in our solution in a way to satisfy the following equality Res M atlab = Res M odel / √ 2. We then conduct the simulations over the terrains described in Table I.For all these three terrains, we assess the execution time to generate their RF maps in both Matlab and our model.Since the resolution step and the number of rays launched determine respectively the precision and the accuracy of the map, we assess how the execution time evolves with this spatial step and the number of rays.The distance is chosen in the range of 50m to 15m with a 5m step, and the number of rays is set to the three values predetermined by Matlab: large, medium, and high.Figures 6, 7 and 8 show the results of our simulation over the three terrains respectively.In each of these figures are given our model and Matlab execution time in hours for the three sets of angular separations as a function of the resolution step.These figures clearly show that Matlab is more complex in terms of execution time than our model for the three terrains, the three angular separations and for all resolution steps.Additionally, Matlab shows an exponential increase in time as the resolution step becomes smaller, thus showing its inability to produce accurate maps in large-scale complex scenarios.On the other hand, our model's execution time remains almost constant for all the resolution steps.This finding is of great interest, since there is only a small overhead incurred by our model when the number of triangles increases, thus allowing the generation of highresolution RF maps capable of precisely giving the bitrate at any location.Moreover, our model shows little sensitivity in its execution time regarding the number of rays, since the difference of execution times between the large and small angular separations is small.We are therefore able to produce accurate and precise maps with an execution time that is considerably smaller than the one of Matlab.Indeed from Figure 6, to obtain an accurate map with Matlab for a 35m precision, one must wait more than 50 times the execution time of our model.Moreover, in Figures 6 and 7, with respectively 133 and 61 number of antennas, Matlab's plots stop at steps 35m and 25m respectively as the tool is unable to produce the RF map for all the three angular separations.This is due to the high complexity of the traditional approach used for RF mapping.Indeed, the large number of reception points leads to intensive matrix calculations exceeding the memory capacity of the machine and therefore to the impossibility of producing the RF map.This demonstrates that with the current RF mapping approach it is impossible to obtain precise maps in large-scale urban scenarios involving hundreds of antennas.Whereas our model easily generates accurate and precise maps in these terrains with a quasi-constant time, thus confirming its out-performance over the traditional approach and its optimized design for RF mapping.
From the above figures, we can also see the impact of the number of antennas on the execution time.The maximum ratio between our model's execution time and Matlab's execution time is 50, 43 and 5 respectively for the three terrains.This ratio decreases with the number of antennas.Fewer antennas mean fewer rays launched, thus leading to less complex matrix calculation to perform by the traditional approach.However, even in those cases, our model still outperforms the traditional  approach and exhibits a slight increase in the execution time as opposed to the exponential increase in Matlab.
Figure 9 shows the contribution of Embree, namely the ray/object intersection test on the overall performance in Terrain 1 (The results hold for the other terrains, as they exhibit the same pattern).From this figure, we can see that Embree's complexity grows with the number of rays launched but is not sensitive to the resolution step.The larger the number of rays launched is, the more ray/object intersection tests Embree will have to perform.However, we see that as the resolution step decreases, Embree's execution time remains constant.This property brought by Embree's underlying technologies explained in Section III-C helps us to choose a fine resolution step without exploding the cost of the intersection test.The combination of Figures 6 and 9 shows the efficiency of our reception test.Indeed, as the resolution step decreases, the execution time of Matlab is exponential due to the large matrix calculation necessary to carry out the reception test.On the other hand, the efficiency of our reception test prevents the explosion of our execution time.This efficiency is further emphasized when Matlab is launched with almost no reception points to decouple the cost of reception test from other processes.Since Matlab requires at least one reception point, we run the same simulation with only 1 reception point.With this setting applied to Terrain 1, Matlab took 67, 59 and 57 seconds respectively for small, medium, and large angular separations (results here hold for Terrains 2 and 3 as well).These values correspond to the cost of Matlab RT process without RF mapping, hence without the reception test.This is typically the cost of the ray launching and the ray/object intersection test processes, which remains unchanged with the number of reception points.Hence, all the other costs incurred by Matlab come from the reception test process, showing how complex this process is.From this finding, we can see that more than 98% of Matlab's complexity comes from its reception process.This highlights the fact that traditional RT approaches are not designed for RF mapping purposes due to the complexity of their reception test process.On the other hand, this shows the effectiveness of our reception test scheme and that it represents most of the gain in terms of execution time.Overall, by combining the efficiency of Embree in performing the intersection test regardless of the number of triangles and that of our new reception test, we can outperform existing models both in terms of execution time and map precision.

C. Sensitivity study
Our technique being based on triangles, we evaluate the impact of the size of the triangles, namely the resolution step, on the accuracy of the bitrate estimation.A first observation from this figure is that all the resolution steps have close distribution to the one of Matlab, in particular for the median and percentiles.In addition, our model is stable to the variations of the resolution step.This result is of great interest since we can choose a resolution step of 50m in our model without any loss in terms of accuracy while minimizing the execution time.However, we can notice from the figure that the smaller the resolution step is, the more outliers there are.Indeed, as the resolution step decreases, the points we consider in this experiment tend to be more likely located in different triangles, therefore leading to this high variance and increased spread in the outliers.This is also the case for Matlab estimates, as they consider individual points.As one would expect, the smaller the resolution step is, the closer our distribution gets to the one of Matlab.Still, our approach is almost not sensitive to the resolution step in terms of the percentiles of bitrate estimation.This helps to produce accurate maps with larger resolution steps in a smaller execution time.

VI. CONCLUSION
This paper presents a novel RF mapping approach able to generate accurate and precise large-scale RF maps.Our approach outperforms existing ones in both execution time and RF map precision.We prove that traditional RT approaches are not adapted for RF mapping purposes.Our experiments show that the reception test process accounts for more than 98% of the total RT time in traditional approaches.
To overcome this limitation, we propose a new reception test that exploits the geometric properties of triangles and rays to quickly determine which rays are received by which triangle.Moreover, thanks to the triangulation schemes introduced in our approach, we consider the 3D topography of the terrain and propose a new ray/object intersection testing scheme that outperforms existing solutions for this process.Simulations confirm that our model can divide the ray tracing end-to-end execution time by 50 for RF mapping purposes.Finally, we show that our approach is not sensitive to the resolution step when generating RF maps, thus allowing the choice of a less precise map without any loss in terms of accuracy.
We are planning in the future to compare our model against real life measurements to assess its compliance with reality both in terms of signal power received and bitrate estimation.Afterwards, this model will be integrated into a tool meant to give the world-wide up-to-date signal and bitrate coverage.
Figure 10 is obtained from Terrain 2, with the low angular separation and the maximum number of reflections set to 5.This figure is a box plot showing the distribution of our model's bitrate estimates regarding the resolution step.The simulation is carried out by choosing 28400 locations in Terrain 2. The bitrate of each of these locations corresponds to the bitrate of the triangle in which they are located.For each resolution step, the bitrates corresponding to these points are obtained and used to plot the figure.We also add to this figure Matlab's distribution of bitrate for the same 28400 locations.

TABLE I :
Characteristics of terrains.

TABLE II :
Simulation parameters.

TABLE III :
RMSE (Mbps) on the 3 terrains.TableIIIsummarizes our findings.From this table we can see that all our simulations show the same trend as the results explained earlier.