[Open-source platforms for fast room acoustic simulations in complex
structures]
Matthieu AUSSAL(1) ,Robin GUEGUEN(2)(1) CMAP - École polytechnique, Palaiseau, France, matthieu.aussal@polytechnique.edu(2) ISCD - Sorbonne Université, Paris, France, robin.gueguen@gmail.com

Abstract
This article presents new numerical simulation tools, respectively developed in the Matlab and Blender environments. Available in open-source under the GPL 3.0 license, it uses a ray-tracing/image-sources hybrid method
to calculate the room acoustics for large meshes. Performances are optimized to solve problems of significant
sizes (typically more than 100,000 surface elements and about a million of rays). For this purpose, a Divide
and Conquer approach with a recursive binary tree structure has been implemented to reduce the quadratic complexity of the computation of ray/element interactions to near-linear. Thus, execution times are less sensitive
to the mesh density, which allows simulations of complex geometries. After ray propagation, a hybrid method
leads to image-sources, which can be visually analyzed to localize sound map. Finally, impulse responses areconstructed from the image-sources and FIR filters are proposed natively over 8 octave bands, taking into account material absorption properties and propagation medium. This algorithm is validated by comparisons withtheoretical test cases. Furthermore, an example on a quite complex case, namely the ancient theater of Orange
is presented.
Keywords: Room acoustics, ray-tracing, fast methods, open-source

1

INTRODUCTION

Today, digital technologies allow research to explore previously inaccessible areas, as virtual reality for archaeology. In this domain, many works focus on the visual restitution, but acoustic studies can reinforce researches
to improve the understanding of the ancient world. For example, during the Roman Empire, architects havedesigned buildings using acoustic rules [3]. In this study, we focus on the ancient theater of Orange which has
a significant size (100m wide), a complex geometry (ornaments, bleachers, columns, arches, etc.) and which is
open-air. In a previous work, a complete mesh was designed using Blender CAD software [12]. To be representative, this mesh, shown in the figure 1, possesses 436 000 elements (triangular faces) and is accurate withrespect to the actual archeological knowledge [5]. As the mesh size makes difficult the use of precise methods(FEM, BEM, etc.), ray-tracing approximation was performed, in order to compute fastly the full-band room impulse response (50 to 15000Hz). Under this assumption, open-source platforms were developed following two
steps. We first build a prototype using Gypsilab, an open source M ATLAB library for fast prototyping [1]. This
preliminary work was useful to construct and validate ideas and algorithms. It leads to the creation of a new
toolbox, openRay, now appended to the master branch of Gypsilab and freely downloadable [14]. In a second
step, all algorithms were retranscrypted in C++ using Qt Creator, leading to an autonomous tool Just4RIR. A
python interface was added, in order to use this library as a Blender plug’in [15]. At the end, this plug’inallows archeologists to only work on Blender, modifying easily meshes and materials, run acoustic simulationand visualize results.
After reminders on acoustical energy propagation represented by ray-tracing, this paper gives implementation
details of the method that was used to obtain a fast computation for large meshes. At the end, application on a
virtual model of the ancient theater of Orange is performed.

6059

Figure 1. Left : Mesh of restituted theater of Orange modelized on Blender (436 000 triangles). Right :Representation of a r-radius measurement sphere centered in xm , receiving energy from a sound source in xs .

2

ACOUSTICAL ENERGY MODELIZATION

2.1 Continuous domain equation
Considering a sound source as a point source centered in R3 [6] and neglecting the losses due to the propagation
medium, the normalized acoustical energy E is given by:ZE=

I · ds = 1,

(1)

Sd

where Sd is a centered d-radius sphere and I the acoustical intensity, such as:
I=

d
.
4πd 3

(2)

The energy obtained by integration on a portion σd of Sd satisfies:
d·n
Ω
ds =
3
4π
σd 4πd

Z

Eσ =

(3)

where Ω define a solid angle, characterizing a flow through the oriented surfaces σd . Thus, Eσ is constant for
all radius d and corresponds to a portion of the total energy E. Finally, subdividing S in N portions σi , E canbe decomposed as a sum of elementary energies:
Z

N
1 N
d·n
1 N
E = ∑ Eσi =
ds
=
(4)
∑
∑ Ωi .
4π i=1 σi d 34π i=1i=1
2.2 Discrete model
To numerically represent the energy propagation of an omnidirectional sound source located in xs , we have tobuild a solid angle basis (Ωi )i∈[1,N] in equation (4). For this purpose, we introduce rays objects with:• Source coordinates xi , with xi = xs ∀i ∈ [1, N],• Unitary direction ui , with an uniform sampling on the unitary sphere (e.g. [8]),1• Energy carried Ei , with Ei = 4π
N ∀i ∈ [1, N].
1 For directional sources, energies can be function of the direction.

6060

As the acoustic source is discretized, we have to define a discrete measure of energy propagation. To this end,
we consider a r-radius measurement sphere S(xm , r), centered on xm (fig. 1). We can then add the contributionsof a n-ray beam that intersects this sphere to calculate the acoustic energy Em at the point xm . In the particularcase of an omnidirectionnal source, we have:
Em ≈

n
1 n
∑ Ωi ≈ N ,4π i=1(5)which means that the measured energy Em is statistically and naturally represented by the ratio between thenumber of rays forming a beam to the total number of rays. Furthermore, considering Ωm the solid angle ofthe measurement sphere in figure 1:
!
r
r2
(6)
Ωm = 2π(1 − cos α) = 2π 1 − 1 − 2 ,dwe observe for r  d that:

r2
n
πr2
⇒ Em ≈ ≈
.
(7)
2
d
N
4πd 2
For far measurements, approximation (7) expresses that Em is well approximated by the ratio of the areas of
the r-radius disk to the d-radius sphere. Otherwise, to ensure the existence of this last approximation, beams
formed by n rays has to be measurable and count at least one ray (n ≥ 1). This assumption is crucial to ensure
the validity of the concept. Thus, fixing a measurement radius r, approximation (7) gives a maximum range ofthe discrete model:
r
r N
d≤
.
(8)
2 n
In addition, figure 2 shows how this statistic modelization represent the distance between source and measures.
The accuracy of the measurement depends strongly on the number of rays counted, then, the more n increases,
the more accurate will be the measurement. Nevertheless, in practice, values for a short distance between thesource and the measurement sphere represent direct sound and first reflections, whereas long distances describe
the diffuse field. Under this assumption, we can consider this model acceptable for all beam such as n ≥ 1.
Ωm ≈ π

2.3 Presence of an obstacle
For the case of acoustic propagation in the presence of an obstacle, we choose to consider only specular reflections (Snell-Descartes laws). Indeed, this approximation is suitable when surfaces are larges in comparison to
wavelengths and diffraction effects can be neglected [6]. For a room, this condition is reached if:
ka  1,(9)with k the wave number and a the characteristic diameter of the room [9]. This approach is currently usedby room acoustic softwares (e.g. Odeon [16], Grasshopper [13], etc.) regarding to audible frequency range(62,5 to 15000Hz). In particular, for the case of the theater of Orange considered in this paper which has acharacteristical diameter of about 50 meters, fixing sound celerity to 340 m/s, the high frequency approximation
is valid for f  1 Hz.
When an incident ray intersects a flat surface, a reflected ray is generated from the collision point. Noting ui
the direction vector of the incident ray, the reflected direction vector ur is defined by:
ur = (ui · T)T − (ui · n)n,(10)

with T the tangent basis and n the normal vector to the surface. Moreover, we modify the energy of thereflected ray by:
Er ( f ) = Ei ( f )(1 − α( f )),
(11)

6061

Figure 2. Left : Measured energy (dB) in function of distance between xs and xm , for r = 0.36 meter and
N = 106 rays. Blue crosses stand for the statistical measure f (r) = n(r)N and red line the analytic function2πr
f (r) = 4πd2 . Right : Sketch of the creation of an image-source by successive reflections of rays on the wallsof a room.

with α( f ) the absorption coefficient of the surface that depends on of the frequency f . In practice, the absorption coefficients are often given per octave bands and can be found in various databases. Both Gypsilab andJust4RIR use the open access Odeon database [16] defined on eight octave bands.
Finally, considering wall absorption, energy measured statistically (eq. 7) is extended by:n
Em ( f ) ≈ (1 − α( f )),(12)N
that we generalize to the case of m reflexions:
Em ( f ) ≈n m(1 − α j ( f )).
N∏
j=1

(13)

2.4 Image-sources
Although the generalized formulation (13) may be sufficient to generate room acoustic data, we also construct
images-sources from the path of rays. To this end, when rays intersect the measurement sphere and following
the reverse return principle, they are retro-propagated along the last direction vector. Thus, rays focus onpunctual images-sources (fig. 2). Each image-source is then located relatively to the listener and carries anenergy according to formulation (13).
By noting (xs )s∈[1,Ns ] the relative position of the Ns image-sources and (Es )s∈[1,Ns ] the associated energy, couples(xs ; Es ( f ))s∈[1,Ns ] contain all useful informations for room acoustic analysis and auralization. First of all, relative
distance of each image-source (ds )s∈[1,Ns ] can be computed. These distances are also used to take into accountthe air absorption, by modifying equation (13) into:
Es ( f ) ≈n −β ( f )ds me∏ (1 − α j ( f )),
N
j=1

(14)

with β ( f ) a frequency dependent absorbing coefficient [11]. Furthermore, fixing the sound celerity c, room
impulse response can be generated, converting each distance ds in time of arrival. Paying attention to the6062

fact that the energy is proportional to the square of the pressure, the finite impulse response can be generated
and analyzed using standard metrics (e.g. T30 , C80 , D50 , etc.). For auralization, this room impulse response
is convolved with an audio signal in order to listen the acoustical rendering. In particular, this convolution
can involve relative position of predominant images sources, in order to realize a spatialized auralization withmultichannel or binaural renderers. Finally, to complete acoustic studies with visual analysis, images sources
can be projected on the room used for computation in order to see where are located listened reflections (seelast impact on figure 2).

3

IMPLEMENTATION

3.1 Standard algorithm
As common principles are introduced, we focus now on the numerical implementation of an acoustic renderer
by ray-tracing. Before any acoustic computation, a numerical room has to be modelized with surfaces andmaterials. In our case, we use the classical representation with a mesh composed of flat triangles.
Using a discrete source point to initalize rays, geometrical intersections are computed between rays (L) andmesh elements (P) using parametric equations:(L) : a + δ u,

δ ∈ R,(P) : b + λ v + µw,

λ , µ ∈ R.

(15)
(16)

Considering v and w driven by two edges of each triangle, the following conditions give pairs (rays;elements)
with uniqueness:• (0 ≤ λ ≤ 1), (0 ≤ µ ≤ 1) and (λ + µ ≤ 1) to ensure that the intersection is inside the triangle,• δ > 0 to respect the propagation direction,• δ minimum not to go through the whole mesh.
Practically, to find these pairs, we can solve directly the underlying linear system or use the Moller-Trumber
algorithm [7]. For N rays and M triangular elements, this process has a quadratic numerical cost (proportional
to NM), which is critical if both N and M are large (see section 3.2).
Then, once all pairs are found, discrete measure has to be done in order to build images-sources (see section
2.4). To this end, the rays are intersected to the measurement sphere S(xm , r), which leads to a linear numericalcost proportional to N.
Finally, a ray is reflected according to equation (10) and propagated while its travelled distance verify condition(8). This iterative strategy ensure the energy propagation by the elimination of all rays that would be in nonmeasurable beams. In the particular case of an open-air room, rays which don’t intersect any surface of themesh are also eliminated. Once all rays are eliminated, images-sources can be built and post-treated (room
impulse response, auralization, etc.).
3.2 Tree-base acceleration
As we have seen, the most critical stage of the standard algorithm is the research of intersections between
rays and triangular elements, leading a priori to a quadratic complexity O(NM). Indeed, each ray has to betested with each face, for each iteration of the ray-tracing algorithm. For a large number of mesh elements andrays which are needed to ensure reasonable accuracy (e.g. M > 105 and N > 106 for the Orange theater), thecalculation time may be prohibitive. To solve this problem, a "Divide and Conquer" approach using a binary
tree is performed.
The general principle consists in creating a mother-box, containing all the mesh elements. This mother-box isthen subdivided along the largest dimension to create two daughter-boxes, each one possesing the same number

6063

Figure 3. Left : Computation time for one iteration of ray-tracing in function of the number of face and rays,
such as N = M (log scale). An omnidirectional source is located at the center of a mesh of a unitary tetrahedral.
Right : All images-sources projected on the mesh (up to RT60).
of elements (median spatial subdivision). This process is then applied recursively, until a stopping criterion isreached. In our case, we stop when the leaves contain only one element. This hierarchical tree is completely
mesh dependent, computed in O(M log M) operations, and gives a structure which permits to quickly navigate
inside the mesh.
Then, we initialize the ray sorting process starting from the mother-box, containing all rays and elements.
Using the first tree-subdivision, we distribute rays inside the two daughter-boxes. This stage is done in O(N)operations by the algorithm proposed by W. Amy et al. [10]. Indeed, each box has N1 and N2 rays, such as
N = N1 + N2 . Assuming this subdivision is performed recursively to the level p, the ith box contains (Ni ) rays
with:
p
2

N = ∑ Ni .

(17)

i=1

Then, ray sorting at the (p + 1)-level also conducts to O(N) operations. To reach the level of the leaves, wehave to perform O(N log M) operations, where log M is close to the depth of the binary tree. At the end, as
we have only one element per box, the ray-element intersection only needs O(N) operations. Finally, instead of
O(NM) operations, we compute all the intersections in:O(M log M) + O(N log M) + O(N),(18)

which operations is a near-linear complexity. Moreover, if the binary tree is precomputed, each iteration of the
ray-tracing algorithm becomes:O(N log M) + O(N).
(19)
To evaluate numerically complexities with or without binary tree acceleration, we measure the computation timeof one iteration by increasing the number of rays and the number of faces in the mesh (N = M). As we can see
in figure 3, the complexity of the algorithm is therefore quasi linear by using tree-based method. This allows totreat large meshes with millions of rays, by keeping a reasonable computation time on a standard laptop (single
core at 2.7 GHz and 8 Go ram).

4

APPLICATION TO ORANGE THEATER

In this paper, we focus only on the Orange Theater acoustics by ray-tracing. However, to evaluate and validate
methods and algorithms, several non-regression tests have been detailled in R. Gueguen’s PhD manuscript [5],associated to open-source implementation on Gypsilab [14] and just4rir [15] websites.

6064

Figure 4. Impulse response in dB (early reflections, up to RT20).
First of all, a restituted version of the Orange theater has been realized on the software Blender (fig. 1),
according to the archeological surveys performed by l’Institut de Recherche sur l’Architecture Antique [2]. As
this mesh just needs to fill the geometry, there is no need of Delaunay properties, but the architecture complexity
however leads to a triangulation with 436 000 faces. In particular, this virtual resitution is mainly composedby the Postscaenium (stage wall) partially ornamented with two basilicas on either side, the Pulpitum (stage),
the Orchestra, the Cavea (bleachers) with Porticus at the top (column gallery) and various covers (stage roof,
velum, etc.).
Several materials are assigned to each part of the theater, in order to define specific absorption coefficientstaken from Odeon database [16]. The ray-tracing solver is then used to compute the spatial impulse response.
To reach the reverberation time close to -60dB (RT60 ), we fix one million rays and a 2m-radius measurement
sphere. An omnidirectional source is located at the front stage, 1.60 m above the floor (this correspond to theposition of the mouth of an average actor). The listener is on the same axis, in the bleachers.
The full ray-tracing propagation was done in few minutes on a standard laptop (2.7 GHz core and 8 Go ram).
At the end, all sources-images projections are generated (fig. 3), illustrating a spatial diffusion of the soundsource. Indeed, as some areas carry a lot of images-sources, reflections seem to surround the listener. Moreprecisely, on figure 4, we see the early reflections of the multi-band echogramm of the theater (until RT20 ). We
can notice a high contribution of the orchestra, the wall-stage, the stage and the roof, as F. Canac demonstrated
in the 60’ [4].
To go further, many others results are available in R. Gueguen’s PhD manuscript [5]. All results correspondto a room adapted for musical playback, more than a speech transmission. These results are confirmed by anequivalent simulation in the Odeon software (commercial license), using the same mesh and parameters.

5

CONCLUSION

In this paper, a full-chain engineering process is given, leading by an acoustical study of an imposing ancientmonument. From archaeological needs, a fine mesh of the monument was constructed, associated to a complete
room acoustic application suite, both for M ATLAB and Blender. The high complexity provided by this typeof architecture and its ornaments leads to approximate calculation methods. Indeed, by only simulating specular reflections and wall absorption, energy propagation and measurement can be simulated by beams, carriedby ray-tracing. From this representation basis, it is possible to generate a multi-band impulse response, while
respecting the laws of high-frequency acoustics. Moreover, a fast algorithm with a near-linear complexity has
been implemented, allowing users to quickly evaluate architectural assumptions, modifying their meshes regardless of the number of elements. At the end, various post-treatments have been added, as the Room Impulse
Response generation, the source-image visualization, the classical perceptive factors and an auralization process.

6065

Even if current versions of proposed softwares (Gypsilab [14] and Just4RIR [15]) are complete enough to be
used for various studies, there are many possibilities of improvement. First of all, as image-source positions
are known, a spatial audio renderer could be added to improve current auralization tool (e.g. binaural or multichannel, eventually with trackers). Secondly, as virtual reality is becoming more and more important in today’s
applications, we could consider moving the listener in real time and thus, allow a complete virtual tour of thebuilding. Finally, as ray-tracing modelization is an high-frequency approximation of waves phenomena, diffraction effects should be added in order to get a better fit with the physical phenomena.

6

ACKNOWLEDGMENTS

The authors particularly wish to thank François Alouges, Titien Bartette, Pascal Frey and Emmanuelle Rosso
for the help they all provided at the different stages of this project. Thanks also to Jean-Dominique Polack
for advices on architectural acoustics and Martin Lesellier for various contributions. This work is part of R.
Gueguen PhD thesis founded by Sorbonne Université.

REFERENCES[1] Alouges, F., & Aussal, M. (2018). FEM and BEM simulations with the Gypsilab framework. SMAI-JCM, vol.
4, p. 297-318.
[2] Badie, A., Fincker, M., Moretti, J. C., Rabatel, L., Rosso, E. & Tardy, D. (2013). Le théâtre d’Orange Rapport final d’opération - Texte. PACA Vaucluse, Orange Théâtre antique 84 087 0031 - Patriarche 9827
n2012-203.
[3] Callebat, L. (1991). Vitruve. De l’architecture. Livre I.
[4] Canac, F., Lejeune, M., & Coulomb, J. (1967). L’acoustique des théâtres antiques: ses enseignements. CNRS.
[5] Gueguen, R. (2018). Virtualisation architecturale visuelle et auditive du théâtre antique d’orange. PhD Thesis.
[6] Jouhaneau, J. (1997). Acoustique des salles et sonorisation. Conservatoire national des arts et métiers Acoustique appliquée, vol 3.
[7] Möller, T., & Trumbore, B. (2005). Fast, minimum storage ray/triangle intersection. In ACM SIGGRAPH
2005 Courses (p. 7). ACM.
[8] Swinbank, R., & Purser, R. J. (2006). Fibonacci grids: A novel approach to global modelling. Quarterly
Journal of the Royal Meteorological Society, 132(619), 1769-1793.
[9] Terrasse, I., & Abboud, T. (2007). Modélisation des phénomenes de propagation d’ondes. Ecole Polytechnique.
[10] Williams, A., Barrus, S., Morley, R. K., & Shirley, P. (2005). An efficient and robust ray-box intersection
algorithm. Journal of graphics tools, 10(1), 49-54.
[11] Iso-9613-1 (1993). Acoustics - Attenuation of sound during propagation outdoors.
[12] https://www.blender.org/[13] https://www.grasshopper3d.com/[14] https://github.com/matthieuaussal/gypsilab/[15] https://github.com/RobinGueguen/Just4RIR/[16] https://odeon.dk/

6066