Carsten Dachsbacher

Marc Stamminger

George Drettakis

Implicit Visibility and Antiradiance for Interactive Global Illumination

Keywords: Global illumination, visibility, GPU

: Left: An office which receives light through a window on the right wall while the rest of the scene is lit indirectly. Our method avoids explicit visibility computation and our GPU implementation enables the manipulation of the light or objects at 9 frames per second (fps). This includes four iterations of radiance and antiradiance per frame. Center: We can also include an animated character; indirect light is updated at 7.8 fps. Right: Our method naturally handles glossy surfaces; the glossy bust is lit indirectly, since the light points at the ceiling.

Introduction

Interactive global illumination is still an elusive goal, despite an abundance of heuristic solutions in recent years. The main expense in most methods is the cost of visibility; a primary goal of this paper is to enable global illumination without requiring explicit visi-bility computation. We start from the observation that if visibility is ignored, some light is transmitted extraneously through opaque objects and needs to be canceled out. This is why we introduce a new quantity called antiradiance that corresponds to light that needs to be removed. We present a reformulation of the rendering equation that requires the consideration of both radiance and antiradiance but enables the treatment of visibility in an implicit manner.

Our reformulation dramatically simplifies the treatment of pairwise interaction between objects: the propagation of radiance and antiradiance between two objects does not depend on occlusion due to other objects, which greatly facilitates parallelization on graphics hardware. The price to pay is that we now need to simulate two quantities, radiance and antiradiance, and the directional distribution of antiradiance must be considered. However, for indirect illumination (or any low frequency illumination), visibility quickly becomes "averaged out" after a few bounces of light and the directional discretization of antiradiance does not need to be extremely fine. For high-frequency direct lighting, we can use a separate traditional technique such as shadow maps. Furthermore, the treatment of directional quantities is desirable to handle glossy materials, although they incur extra cost for the local shading integral.

Our paper makes the following contributions:

A reformulation of the rendering equation, where visibility is implicit and antiradiance compensates for extraneous transport.

Two iterative solutions, one that provably converges but requires more iterations, and one that converges in practice and is cheaper.

Hierarchical discretization based on a spatial and directional data structure, with appropriate refinement and pre-filtering.

An efficient GPU implementation, which allows real-time global illumination for scenes which are hard to treat with all previous methods (moving objects and lights, varying material properties, mainly indirect illumination).

We show that our method can compute solutions which are close to reference images from a path-tracer [START_REF] Pharr | Physically Based Rendering from Theory to Implementation[END_REF], and illustrate our algorithm on several test examples.

Related Work

Global illumination is a vast field and excellent introductory books exist, e.g. [Dutré et al. 2006;[START_REF] Pharr | Physically Based Rendering from Theory to Implementation[END_REF]. We only discuss work that is most directly related to our approach. "Negative light" was used to allow for incremental radiosity updates, where it compensates for the different visibility configuration between two frames [START_REF] Buckalew | Illumination networks: Fast realistic rendering with general reflectance functions[END_REF][START_REF] Pharr | Physically Based Rendering from Theory to Implementation[END_REF][START_REF] Chen | Incremental radiosity: An extension of progressive radiosity to an interactive image synthesis system[END_REF]]. Shadow photons [START_REF] Jensen | Efficiently Rendering Shadows Using the Photon Map[END_REF] are also related to the idea of negative light. Note that visibility computation is still required, in contrast to our use of antiradiance. In particular, methods for the efficient update of global illumination [START_REF] Shaw | Hierarchical radiosity for dynamic environments[END_REF][START_REF] Drettakis | Interactive of global illumination using a line-space hierarchy[END_REF]] require searching for regions that need recomputation of visibility. Our reformulation avoids this search and does not require any explicit visibility computation.

Bunnel [2005] uses "negative light" to approximate ambient occlusion and simulates the effect of inter-reflection. In contrast to our work, his negative light is not directional, and this heuristic does not respect the rendering equation.

The most related work is by [START_REF] Pellegrini | Rendering equation revisited: How to avoid explicit visibility computations[END_REF] who also derives a new rendering equation where explicit visibility is avoided. We use similar ideas where negative light is transmitted through each surface to compensate for the lack of occlusion treatment. However, his work does not demonstrate an implementation and does not discuss convergence and directional discretization issues.

A number of solutions for dynamic global illumination perform partial computations and use caching or reprojection of results of previous frames. An example is the Shading Cache [START_REF] Tole | Interactive global illumination in dynamic scenes[END_REF]; this approach caches samples from a path tracer and uses graphics hardware for interpolation. However, these approaches still require visibility to compute the sparse samples.

Our approach builds on ideas from clustering methods for hierarchical radiosity [START_REF] Sillion | A unified hierarchical algorithm for global illumination with scattering volumes and object clusters[END_REF][START_REF] Smits | A clustering algorithm for radiosity in complex environments[END_REF]], since we create links between clusters of objects to accelerate computations for more complex scenes. These approaches however require sophisticated data structures and global random-access visibility computations for form-factors. Extensions to this method stored directional radiance using wavelets [START_REF] Stamminger | Efficient glossy global illumination with interactive viewing[END_REF] or spherical harmonics [START_REF] Sillion | A clustering algorithm for radiance calculation in general environments[END_REF]. In contrast, we store a directional information of "negative light", thus avoiding visibility computation.

Precomputed Radiance Transfer is a method to quickly compute the reflection integral for a given environment, e.g. [START_REF] Sloan | Precomputed radiance transfer for real-time rendering in dynamic, lowfrequency lighting environments[END_REF][START_REF] Kautz | Fast, arbitrary BRDF shading for low-frequency lighting using Spherical Harmonics[END_REF][START_REF] Sloan | Local, deformable precomputed radiance transfer[END_REF]. Local inter-reflection between an object and the environment can be approximated, e.g. [START_REF] Sloan | Precomputed radiance transfer for real-time rendering in dynamic, lowfrequency lighting environments[END_REF]], but all these approaches require significant precomputation and do not simulate the balance of light due to global illumination. Our technique is complementary and combining our global simulation with precomputed local interaction is an exciting avenue of future work. [START_REF] Kristensen | Precomputed local radiance transfer for real-time lighting design[END_REF] simulate global illumination but for a fixed geometry and require substantial precomputation. [START_REF] Kontkanen | Wavelet radiance transport for interactive indirect lighting[END_REF] use a 4D discretization and a hierarchy, but they need to compute explicit visibility and require lengthy precomputation to treat static scenes only. Hasan et al. [2006] precompute a direct-to-indirect light transfer matrix for interactive GPU-based relighting, but are restricted to a fixed camera and static scene. [START_REF] Ren | Realtime soft shadows in dynamic scenes using spherical harmonic exponentiation[END_REF] proposed a different reformulation of visibility using the exponentiation of spherical harmonics in the context of soft shadows. They however still need to determine the occluders between an object and the light and do not treat global illumination.

A number of authors have mapped global illumination algorithms on graphics hardware, e.g. [START_REF] Coombe | Radiosity on graphics hardware[END_REF][START_REF] Purcell | Photon mapping on programmable graphics hardware[END_REF]]. However, they need to cope with visibility computations that are not GPU-friendly, raising important implementation challenges and performance limitations. [START_REF] Keller | Instant radiosity[END_REF] places a large number of virtual light sources, requiring visibility treatment, and accumulates their contribution with hardware shadow mapping and blending.

Overview

To achieve rapid and GPU-friendly global illumination, we get rid of explicit visibility computation such as ray casting and reorganize the traditional rendering equation by reformulating occlusion in an implicit manner. This allows us to consider pairwise object interaction as if all other objects were transparent and did not cause occlusion, thereby greatly simplifying computations.

To compensate for missing occlusion, every surface point emits incident light backwards as "negative light". We refer to this "negative light" as antiradiance. Antiradiance is propagated in the same way as normal radiance and cancels out the light that incorrectly traversed surfaces due to omitted occlusion. We thus propagate two quantities: radiance and antiradiance and it is critical to consider the directional distribution of antiradiance. At this price, we can replace the very costly occluded propagation of traditional global illumination by a much cheaper unoccluded propagation step. We illustrate the basic concept in Fig. 2. We formalize this process by reformulating the rendering equation and demonstrating the correctness of this reformulation. We also provide an iterative solution method, and discuss speed and convergence choices.

Since antiradiance is a directional quantity, we use spatial and directional finite elements to solve the system of equations. We use appropriate pre-filtering in space and directions for form-factor computation, and develop a hierarchical solution to deal with complexity. Our solution has been designed to be stream-processor friendly. The actual GPU implementation includes several important design choices, allowing us to achieve interactive global illumination updates, with moving objects and moving lights.

Reformulating the Rendering Equation

Traditional rendering equation The Rendering Equation [ [START_REF] Kajiya | The rendering equation[END_REF]] can be expressed in terms of linear operators that act on radiance L(x, ω), where x is a point in space and ω a direction [START_REF] Arvo | A framework for the analysis of error in global illumination algorithms[END_REF]] (Fig. 3). Please note that we follow the convention of [START_REF] Arvo | A framework for the analysis of error in global illumination algorithms[END_REF] for the directions and signs in the operator formulation. The reflection operator K takes incident radiance L in and performs the shading integral to output outgoing radiance L at each point x,

(KL in )(x, ω) = Ω nx f (x; ω in → ω)L in (x, ω in ) n x | -ω in dω in (1)
where n x is the normal at x, f the BRDF and | the dot product. Incident light L in is obtained through the (occluded) geometry operator G which uses an explicit visibility function ray(x, ω) that determines the first front-facing point starting at x in direction -ω (rays that reach no hit point return a virtual black point p inf ):

L in (x, ω) = (GL)(x, ω) = L(ray(x, ω), ω) (2) E(x, ω)
is the self-emission, and the rendering equation becomes:

L(x, ω) = E(x, ω) + (KGL)(x, ω).
(3)

Reformulation

We seek to replace the transport operator G by an operator U that does not require explicit visibility. That is, we define U using a modified ray function RAY that returns the set of all intersection points (Fig. 3b) in direction -ω:

L unocc in (x, ω) = (UL)(x, ω) = ∑ y∈RAY(x,ω) L(y, ω) (4)
Informally, using U directly would propagate excessive light and we need to compensate for it. This is why we simulate a second quantity, antiradiance, that cancels out this extraneous radiance. Note that each time light reaches a surface, it gets extraneously

L GL JGL GJGL UL (a) (b) (c) (d) (e) UJGL (f) 0 1 2 3
Figure 4: Radiance after unoccluded transport U is the sum of radiance after a series of occluded transport G and "go-through" J.

propagated through it by U. This is why we create antiradiance that corresponds exactly to the amount of light incident to an object. In this manner, the light wrongly transmitted during a propagation gets canceled by antiradiance during the following step.

We now formalize this and relate the unoccluded U and occluded G operators. To make the fact that U ignores occlusion explicit, we define a "go-through" operator J that lets incident light through opaque objects (Fig. 3b lower right). It is essentially the identity:

(JL in )(x, ω) = L in (x, ω). ( 5
)
With the help of J, we can describe the relation between G and U. The idea is shown in Fig. 4 with L a Dirac impulse at the top patch 0. GL is non-zero at patch 1 only. By letting light pass (J) and propagate it again using G, it reaches patch 2. Every further application of GJ propagates light one layer further. U is finally the sum over all layer depths:

UL = GL + GJGL + GJGJGL + ... ( 6
)
where the terms become zero when G(JG) i reaches scene depth. This gives us the fundamental equation that relates G and U:

UL = GL + (G + GJG + ...)JGL = GL + UJGL → GL = UL -UJGL (7)
where the left (positive) term UL is unoccluded propagation of radiance and the right term -UJGL propagates "negative light" to compensate for the lack of occlusion. We call JGL antiradiance, and denote it by

A: A = JGL (8)
Our definition of antiradiance still includes visibility, but we can replace GL by Eq. 7 to get a recursive equation for A without G:

A = J(UL -UJGL) = JU(L -A) (9)
This equation allows us to determine the antiradiance field A for any radiance distribution L. We start with A = 0 and iteratively update A. After n iterations, where n is the maximum scene depth complexity, the computation converges and we obtain the antiradiance distribution A corresponding to L by an iteration using U (instead of a single application of G). With the knowledge of A, we can compute one iteration step of classical light transport:

GL = UL -UJGL = U(L -A) (10)
We finally reformulate the rendering equation ( 3) by replacing G with its expression in terms of U and antiradiance A (Eq. 7, 8)

L = E + KGL = E + KU(L -A),
i.e., we replace G by U and propagate L -A instead of L. Our new rendering equation is

L = E + K U(L -A) (11) A = J U(L -A). ( 12
)
Note the similarities between Eqs. 11 and 12. Both are recursive and propagate L-A using the unoccluded transport U. L is reflected via a BRDF in K, antiradiance A with the simple operator J, which is essentially the BRDF of a transparent surface. The recursion is interdependent, since L and A are defined based on each other.

Iterative Solution

Before presenting our finite-element solution to the new rendering equation, we discuss iterative strategies and convergence. We propose two schemes that differ in the relative number of steps of Eqs. 11 and 12. The first scheme has been described in the previous derivation and is asymmetric. It applies one radiance propagation and then antiradiance steps until convergence to emulate one step of GL. It corresponds to traditional global illumination and hence provably converges. The second method is symmetric and alternates Eqs. 11 and 12. While we cannot formally prove convergence, our experiments have shown that it is stable in practice; we present our intuition on why this is the case.

Asymmetric Iteration

In our first technique, after each step of radiance propagation (Eq. 11), we iterate over Eq. 12 to propagate antiradiance until convergence. So after one step of

L (i+1) = E + K U(L (i) -A ( j) ), (13)
we apply multiple steps of

A ( j+1) = J U(L (i) -A ( j) ), (14)
where L (i) is frozen while j iterates. The antiradiance iteration effectively computes A = JGL; the required number of iterations is bound by the scene's depth complexity.

Symmetric Iteration

In our second technique, we solve the joint equations 11 and 12 by iteratively propagating radiance (Eq. 13) and antiradiance (Eq. 14). Such a scheme converges only if the iteration matrix does not have eigenvalues with magnitude higher than one. In traditional radiosity, this can be easily proven because the matrix is diagonal-dominant. We cannot make such a claim here, because the unoccluded transport operator U can create energy.

In fact, in simplified 2D examples with multiple occlusion and perfect mirror-BRDFs, eigenvalues larger than one occur. However, in practice the symmetric propagation seems to converge faster and we did not observe divergence. We hypothesize that there is much natural dampening due to the shading integral. Highly-specular BRDFs might create more problems and deserve further investigation in future work.

Discussion

We measured convergence for the office scene from apply one radiance propagation and four antiradiance propagation steps, which become visible as steps in the graph.

For both schemes, we did more than 300 iterations. No divergence could be observed, the difference between both final results was less than 10 -4 . The graph shows that the symmetric iteration is better in the beginning, however the asymmetric iteration overtakes after 90 steps. Visually, convergence is achieved after less than 10 iterations, and in this range the symmetric iteration has better behavior. However, the graph does not differentiate between antiradiance and radiance iteration steps. In particular for complex BRDFs antiradiance propagation is generally faster in computation time, because the application of J is very cheap compared to K. For the rest of this paper, we focus on the symmetric scheme.

Discrete Finite Element Solution

To solve our new rendering equation, we use spatial and directional finite elements. Our approach uses a hierarchy to handle complexity. The refinement strategy uses both spatial and directional information, and is well adapted to moving objects and lights. We also show how different direct lighting approaches can be used such as directional lights or environment maps.

Discretization

We discretize the scene into a hierarchy of patches P i with centers x i , and the space of directions into m bins Ω j with main directions ω j . Each bin covers a solid angle of ω bin = 4π/n, in practice bins span between 6 and 20 degrees.

For each patch P i , we need to store a discretization of exitant radiance L(x, ω) and antiradiance A(x, ω). To make hierarchy maintenance simpler, it is convenient to use intensities I L (x i , ω j ) and I A (x i , ω j ), i.e. flux over solid angle. We only treat opaque solids, so that operators K and J update separate hemispheres (Fig. 6(a)). Thus we can store the quantity exitant intensity

I ex (x i , ω j ) = I L (x i , ω j ) -I A (x i , ω j
) in a single data structure. We also store incident light as radiance L in (x i , ω j ). Finally, we store total intensity I total , the final intensity used for display.

The iterative computation consists of two steps: in the global pass exitant intensity I ex = I L -I A is propagated from a sending to a receiving patch (depending on geometric properties of the two samples, but not on visibility) and contributes these to L in (x i , ω). This corresponds to the operator U. In the subsequent local pass, this incident radiance is transformed into exitant intensity and antiradiance for the next iteration, thus computing operators K and J.

Global pass

The global pass is very similar to the link pass of standard radiosity, and in its simplest form it is an O(n 2 ) process, because each of the n patches exchanges light with each other. Consider a particular pair of patches (P i , P k ) with P i being the receiver and P k the sender. We first assume that, thanks to our hierarchical decomposition, patch P k occupies only one directional bin as seen from the center of P i ; we lift this restriction later.

The emitted intensity at P k towards P i is

I i←k = I ex (x k , x i ← x k ),
where x i ← x k is the direction linking the centers of the two patches.

To compute the influence at x i , we position a virtual infinitesimal surface element dA at x i that directly faces x k (see Fig. 6 (b)). The flux from x k emitted towards dA is then:

dΦ = I i←k dA ||x i -x k || 2 . ( 15
)
With our discretization, we first have to find the directional bins, via which sender and receiver interact (see Fig. 6 (c)). Let the directional bins at x i and x k , in which this light arrives, be ω j and ω l .

We have to distribute the incident flux uniformly over this receiving bin, so the incident radiance at x i is:

L i j←kl in (x i , ω j ) = dΦ dAω bin = I ex (x k , ω l ) ||x i -x k || 2 ω bin . ( 16
)
The "form factor" coupling the exitant light at P k and the incident light at P i is thus:

F i j←kl = 1 ||x i -x k || 2 ω bin . ( 17
)
Note the difference to standard radiosity form factors: the cosine at the sender and the sender's area are implicitly stored in the intensity distribution I(x k , ω l ). The cosine at the receiver is computed during the local pass described below.

(a) (b) Interacting pairs of sender and receiver are stored as links. For storage convenience, we only consider solid objects. This implies that we establish bidirectional links between bins which are both in the positive hemisphere, while bins in the negative hemi-sphere only have outgoing links; see Fig. 7.

Up to now we computed the incident light at the center of P i and assumed that the energy of the sender P k is contracted to its center. The first assumption corresponds to point collocation in standard radiosity [START_REF] Sillion | Radiosity and Global Illumination[END_REF]. The latter assumption, however, can lead to artifacts if the sender patch covers a large solid angle and spreads over several bins. Our oracle refines a link if the solid angle of the sender is too large. However, when maximum refinement or patch area thresholds are reached, this is not always the case.

For an accurate CPU implementation, we advocate the precise computation of exchanges between all pairs of bins in the cone subtended by the sender as seen from the receiver. However, for graphics hardware, we need to make computations as local as possible and use the following heuristic: for each link, we compute the solid angle of the sender. If this is larger than that of a bin, we "spread" and filter its energy to the bin's neighbors. That is, we select a circular neighborhood of size σ , the solid angle subtended by the receiver and in a sense filter the form factor with a Gaussian blur similar to prefiltering for antialiasing. Assuming an optimal uniform bin distribution the angle between a bin and its closest neighbor is 4π m ; this is the standard deviation of our Gaussian. This means that the appropriate bins receive energy, but that this energy is taken from the wrong bin: that defined by the direction between the patch centers. However, link refinement ensures that an exchange usually covers only one or very few bins; Fig. 8 demonstrates the effectiveness of this heuristic.

Local pass

In the local pass, we transform incident radiance L in to reflected radiance and antiradiance, represented as intensities I L and I A . The reflection is defined by Eq. 1, however we have to transform exitant radiance to intensity. At this step, the surface area a i represented by a sample i is needed as well as the sample's normal n i . We can then derive:

I L (x i , ω j ) = a i n i |ω j L(x i , ω j ) = a i n i |ω j ω bin ∑ m-1 l=0 f (x i , ω l → ω j ) n i | -ω l L in (x i , ω l ) (18)
For antiradiance, we have to transfer the light incident from bin ω j to the opposite direction and transform it to intensity:

I A (x i , ω j ) = a i n i | -ω j L in (x i , ω j ) (19)

Hierarchy and Refinement

To make computations tractable, we solve the above finite element problem hierarchically in a spirit similar to Hierarchical Radiosity with Clustering [START_REF] Smits | A clustering algorithm for radiosity in complex environments[END_REF][START_REF] Sillion | A unified hierarchical algorithm for global illumination with scattering volumes and object clusters[END_REF]]. We use a hierarchy over the scene and simulate exchanges through links at the appropriate level. The refinement is the same for both quantities.

Our algorithm starts with the creation of links, which is the only step which must happen on the CPU since it requires random global access to the scene. We start by considering exchanges from the root cluster to itself. These self-links are always subdivided, i.e. replaced by links between all pairs of children. For a non-selflink we use a simple oracle for refinement: if the solid angle of the sender with respect to the receiver is larger than the solid angle of a bin of the directional discretization, its light will arrive at the sender within several bins of the receiver's incident light discretization. In this case, we subdivide the link, i.e. we subdivide sender and/or receiver and consider the links between the children recursively. Otherwise, we establish the link by computing its form factor and storing the link. By this, the epsilon threshold of Hierarchical Radiosity that controls the accuracy of the solution is linked to the directional discretization accuracy. Evidently, link refinement is stopped if sender and receiver become smaller than a discretization threshold. This step works very well with clustering: we do not need a normal and thus do not have to use average normals, nor consider self-occlusion.

Note that we do not consider the emitted light of the sender as this is usually done in radiosity (BF-refinement). This would make it necessary to refine links after each iteration, which does not fit to our GPU implementation.

Maintenance of the hierarchy

The hierarchical representation must be maintained using a push-pull method. We push all light gathered in interior nodes to the leaves after the global pass. Then the local pass is performed at leaves. After the light has been reflected and antiradiance generated, the light is pulled up to update the inner nodes' light. Since we pull intensity, an inner node's intensity is simply the sum of its children's intensities.

Direct light and Dynamic Scenes

Direct Light Initialization Our approach is geared towards indirect illumination. For direct illumination, we can initialize the incident radiance of the samples in the scene, using one of the many existing methods. We can for example initialize all the L in using an environment map. For direct shadows, we can use shadow maps or one of the more recent soft shadow methods [Guennebaud et al. 2006]. We evaluate direct lighting on each of the samples in the scene, for example with the shadow map, and initialize the exitant radiance. No antiradiance is generated in this pass, because shadowing of the direct illumination is done by the shadow map. During rendering, both parts of the lighting computation, fast per-pixel direct lighting with shadows and indirect lighting from the coarser finite-element illumination is summed in screen-space.

Moving Objects and Lights

Moving or deforming objects are handled naturally and easily with our approach. Since we do not consider visibility, we only have to update links to and from moving objects, and we do not have to search for all links that are modified by changing visibility. The situation is more involved for hierarchical computations, because for a moving object the optimal link hierarchy changes with the position. However, if the position of moving objects can be reasonably bound, we can also bound the oracle and generate a link mesh based on the worst case. The link hierarchy thus does not need to be updated, but only the unoccluded form factors need to be recomputed, which in turn is a cheap operation. At price of a less stringent oracle, the sender's emission usually also can be limited by reasonable bounds. One powerful feature of our approach is that moving lights can be handled in the same way, albeit with the same limitations.

Implementation on Graphics Hardware

One motivation for our approach is to obtain a GPU friendly global illumination algorithm. The global and local passes described above can be easily and efficiently performed on graphics hardware.

The CPU only needs to compute the hierarchy of links which is cheap, because no visibility has to be computed. We now describe data layout and the mapping of the most important operations.

Data Structures

The element hierarchy is laid out in depth-first order, which allows for efficient push/pull operations as we will see later. We use three tables with n rows (n being the number of elements) each with m (number of bins) columns, storing exitant intensity I ex , incident radiance L in and total intensity I total for every element (Fig. 9). Note that when stored in textures the tables may be split into multiple texture columns as texture sizes are restricted by hardware (Fig. 10). All per-element attributes required for computations (normal, position, area, color, specular exponent) are stored in textures with the same layout but smaller width, so that we can use the same texture coordinates as we use for accessing the element bins.

Links are stored as render primitives in vertex buffers (Fig. 11,left) that includes the indices of the sender, receiver, and bin as well as the form factor value and the cone angle for filtering.

Operations

For each iteration, we perform the following four operations working on the aforementioned tables stored as textures and used as render targets. As GPUs do not allow the use of a render target as texture at the same time we use intermediate render targets or pingpong rendering where necessary.

Global Pass

In the global pass (Fig. 11) energy is transferred from one outgoing bin of a sender element, i.e. one entry in the I ex table, to the bins of a receiver. Sender and receiver positions and sizes are used to find the bins used in our form-factor approximation (Sec. 3.1). The bins affected and appropriate filtering are precomputed for a finite number of cone sizes σ . For each σ we store an m 2 convolution table (Fig. 11 upper right): for each outgoing direction (rows), we store the contribution using the filtered form-factor to all bins of a receiver (columns). Some receiver bins may not receive energy; in this case the corresponding lookup table entry is zero.

Using this method, global energy transfer becomes rendering of the line covering all bins of the receiver element in L in . As a result exitant energy weighted by the filter is transferred to the receiver bins as incident radiance (Fig. 11). For dynamic links, the affected bins and filter parameters are determined on the fly, otherwise they are computed once in the linking step.

Push Operation After the global transfer, we push energy down the hierarchy. Due to depth-first order all child elements are stored consecutively and we can add a node's radiance to several (or all) children at once by rendering quad primitives.

Local Pass For the local pass, we iterate over all bins of every element: we convolve the incident radiance with the BRDF to obtain the reflected energy and transform the incident radiance of the opposite bin to antiradiance. The reflected intensity is also accumulated in the I total table, which is used for display. For diffuse surfaces, the BRDF computation can be accelerated: The reflected intensity in the normal direction is computed, stored in intermediate textures and used to compute reflected energy for all bins.

Pull Operation After the local pass we pull energy from the leaves to their parents. We read the intensity bins of a leaf element (from I ex ) and add it to the bins of all its predecessors in the hierarchy (in I ex ). This is achieved by rendering textured lines. Vertex positions are chosen such that the lines cover all receiver bins, while texture coordinates define the source bins. We store render primitive data in vertex and index buffers. Note that intensity I ex is stored after weighting by element area and thus no further scaling is required.

Interpolation by Splatting

For final rendering, we need to interpolate the sample values to obtain a smooth solution. For shared vertex meshes, we compute the total intensity at the vertex locations and use linear hardware interpolation. For flat surfaces, we use a three pass rendering approach.

Similarly to splatting for point-based rendering [START_REF] Lavignotte | Scalable photon splatting for global illumination[END_REF][START_REF] Ren | Object space EWA surface splatting: A hardware accelerated approach to high quality point rendering[END_REF], we start with a pass which writes surface IDs to a render target: we interpolate across surfaces patches sharing the same ID, which is determined after linking as a preprocess. In a second pass we render triangle fans centered at ele-ments' centers and accumulate radiances weighted by a smooth falloff. The distance to neighbor elements can be used to deform the fans and better preserve local features. A third pass re-normalizes the radiance values by the total weights.

Results

All the tests were run on a Xeon processor running at 3Ghz, with an NVIDIA GeForce 8800GTX graphics card. Reference solutions are computed with PBRT [START_REF] Pharr | Physically Based Rendering from Theory to Implementation[END_REF], using the pathtracer module (faster algorithms exist but they tend to be biased and we focus here on accuracy).

Validation on Simple Scenes

We present a few simple scenes to investigate some sampling issues and overall validation of our method. We first show a geometrically simple maze scene (Fig. 12), in which all lighting is indirect. A directional light enters from above, lighting the wall facing away from the viewer; all light and shadows in this image are indirect. This is a particularly tough case as can be seen with the noisy path-tracing solution, despite the 12K rays/pixel. Our solution shown here, with 1024 bins, spatial subdivision for a total of 6088 elements, and 8 iterations took 3min 45s (including linking) computed on the CPU, while the path tracing solution for the image (350 × 330 pixels) took 13h 39min. As with any finite element basis, sampling rate affects the result.

The effect of directional sampling rate is intuitive, essentially resulting in blurring of shadows. Consider Fig. 13 in which the left image is computed with 128 bins and the right-hand image with 512 bins. We can clearly see the sharpening of shadow detail with higher directional sampling. Example Scenes We show three more involved example scenes. The first is an office scene, which receives a directional light coming through the window, illuminating the right wall. The interior is illuminated almost exclusively by indirect light. The second scene is a simple non-diffuse environment, with a glossy object. The third scene is a oriental living room (Fig. 16). We use shadow maps for direct light in the office and living room. Recordings of interactive sessions for all three sequences are shown in the video.

In Fig. 1 (left), we show the office scene. This scene contains 5,552 input polygons and is subdivided to 18,253 elements and 632,098 links. We use 128 bins for directions. We can move the position of the light, or any object in the scene at 7-10 fps for 4 iterations. For example, we can include an animated character with full interactive update of global illumination in the scene (center). This scene has 19,716 hierarchy elements, 660,783 links and 8,072 input polygons.

To compare the lighting solution against other methods, we chose to remove textures which tend to hide artifacts. A texture-free rendition of the office is shown in Fig. 14(a) (running at about 9.0fps), while Fig. 14(b) shows a solution on the CPU with improved lighting simulation; Fig. 14(c) is a path-traced solution computed in PBRT, which we use as a reference. This solution uses 2048 rays per pixel, a path length of 4, and took 8h 25min to compute and is still noisy. We also show the absolute difference between the GPU solution and the reference (Fig. 14(d)). An Instant Radios-ity [START_REF] Keller | Instant radiosity[END_REF]] solution running at 9 fps with 190 virtual light sources (Fig. 14(e)) suffers from artifacts like bright spots, incorrect lighting levels, shadows which are too sharp, and banding. To generate Fig. 14(f) we used our link structure to compute a Hierarchical Radiosity with Clustering solution using explicit visibility which takes 80min.

The second scene includes a glossy object (3072 triangles) lit by a spot light. We use a Lafortune BRDF as implemented in PBRT.

For low-frequency glossy lighting, our approach can achieve a good approximation (Fig. 15, left) compared to a reference solution computed by PBRT with path-tracing (Fig. 15, right). This example at interactive quality with 4550 elements and 177,000 links running at about 12 fps for 3 iterations is shown in the video. A more complex object (12,940 input polygons for the bust) is shown in Fig. 1 right. The glossy bust is lit indirectly, and we can move the spot light, change view point, and modify material parameters at about 6 fps.

We show statistics for all scenes in Tables 1 and2. Table 2 shows the breakdown of computation time and memory usage. Unsurprisingly, the local pass dominates for glossy scenes, while it is the other way round for diffuse scenes. Overall GPU memory usage (< 140Mb for all scenes, including optional look-up textures for speed-up) is reasonable for modern high-end graphics cards. On the GPU we store intensities and radiances as 16-bit floats and link data as 32-floats; no precision problems were observed.

Discussion

As illustrated in the results and in the accompanying video, our method provides fast interactive indirect lighting for interesting scenes. The main approximation in our approach is the directional discretization. For indirect light, the shadows and the overall lighting effects are generally smooth anyway, and thus this does not create too many visible artifacts. We consider that the comparison in overall lighting levels with our reference path-traced solutions is satisfactory for the interactive application context we target.

There are limitations in our current implementation: Increased memory usage for the directional data structures and excessive subdivision for dynamic scenes. At the numerical level, the convergence of the symmetric scheme cannot be ensured; however, the asymmetric one converges, albeit slower. The excessive subdivision for dynamic links is a design choice. In the context of a content creation workflow, we could imagine a tool which allows the designer to specify regions in which this should happen, based on knowledge of the usage scenarios. Update of the patch subdivisions and links can of course be made at runtime, but would involve CPU computations and upload to the GPU. Finally, in all the examples we have shown, the computation of 3-4 iterations is sufficient both in quality and in speed to avoid artifacts. Clearly, for more complex scenes, this could be a problem: hierarchical or space-subdivision schemes could potentially provide a solution to this issue.

Our algorithm shares some of the limitations of hierarchical radiosity methods in what concerns refinement and meshing. Both are hard problems; we hope that the availability of interactive global illumination using our approach will accelerate research in this area and hopefully open the way to practical solutions.

We could also handle transparent surfaces, however the operators J and K would interfere and could not share a single data structure thus increasing the memory required.

Our method is related to the idea of "negative light" (Sect. 1.1) re-cently used by [START_REF] Bunnell | 2:Programming Techniques for High Performance Graphics and General-Purpose Computation[END_REF] on the GPU. However we present a proper definition of "negative light" (antiradiance) and the appropriate theoretical basis showing that we solve the rendering equation. In addition, our algorithmic solution is completely different, since we store a directional representation for radiance/antiradiance, and we have reorganized light-transfer into global and local computation. Our hierarchical solution also permits the treatment of more complex scenes.

Figure 16: Oriental room scene; indirect light is updated at 6.8 fps.

Conclusion

We have presented a reformulation of the rendering equation which uses implicit visibility. This is achieved by using a new quantity, called antiradiance, which is propagated with radiance, and is stored in directional elements in the scene. This results in a reorganization of the computation and data structures which are local, and much easier to parallelize. As a result, we can achieve interactive updates of global illumination for moderately complex scenes, with moving objects and lights, and glossy reflectors.

In future work, we believe that the reformulation provides a different way of thinking about visibility, and could also lead to results in other applications in computer graphics such as occlusion culling. A more theoretical analysis of the sampling and filtering issues would be interesting and should lead to more efficient and highquality solutions. We are also considering more efficient rendering, extensions to more involved lighting phenomena and resolving the meshing/refinement issues discussed above.

Figure 2 :

 2 Figure 2: Top, left to right: after one iteration there are no shadows, since antiradiance has not been propagated. At two iterations shadows are too dark; after three they are too bright (yellow box). Bottom: with 4 iterations the solution is almost converged, compared to 8 iterations and the path-traced reference.

Figure 3 :

 3 Figure 3: (a) Operator formalization of the rendering equation. (b) Our reformulation uses unoccluded transport but creates antiradiance to compensate for extraneous transport.

Fig. 1 .Figure 5 :

 15 Figure 5: Convergence comparison of asymmetric (red) and symmetric iteration (blue) for our office scene (Fig. 1). (a) linear plot. (b) log plot.

Figure 6 :

 6 Figure 6: (a) We can store I ex , which is the difference I L -I A because the operators J and K affect separate hemispheres. (b) The interaction between x i and x k . (c) Finding the interacting bins.

Figure 7 :

 7 Figure 7: We only treat solids: (a) Bidirectional links are only established between bins both in the positive hemisphere. (b) Bins in the negative hemisphere only have outgoing links.

Figure 8 :

 8 Figure 8: Image (a) without filtering (b) with filtering. We use 4x4 samples on the blocker, accentuating the effect for illustration.

Figure 9 :

 9 Figure 9: The global pass transfers exitant energy between elements. The local pass first reflects light and generates antiradiance and second tracks the element's intensity for display.

Figure 11 :

 11 Figure 11: Energy transfer via a link by rendering a line primitive: the receiver index of p gives the position, sender index q and outgoing bin ω determine the bin and precomputed spreading/filtering.

Figure 12 :

 12 Figure 12: Left, our solution for an entirely indirect lighting situation: a directional light enters the window at the top of the scene; all shadows here are indirect. Right a path-tracing reference solution.

Figure 13 :

 13 Figure 13: Effect of directional sampling. (a): 128 bins, (b) 512 bins; lower directional sampling rates result in blurring of shadows.

Figure 14 :

 14 Figure 14: The untextured office scene for comparison (left to right, top to bottom): (a) The GPU solution at 9 fps. (b) A more refined solution computed on the CPU in 1m45s with 44,982 elements, 208 directions and 4.6 million links. (c) path-tracing reference, 2048 samples/pixel (path length 4). (d) absolute difference between GPU and path-tracing (variations at discontinuities mainly due to different anti-aliasing). (e) Indirect illumination computed with Instant Radiosity running at 9 fps. (f) Using the same link mesh as for Antiradiance, computing form-factors with visibility (141 rays per link) takes 80min, resulting in a solution akin to Hierarchical Radiosity with Clustering. Note that a visibility-based refinement would improve the quality.

Figure 15 :

 15 Figure 15: (a) A scene with a glossy object rendered with high quality settings (12.5K elements, 480K links, 8 iterations, 1.1 fps). Note the bluish directional indirect illumination from the glossy cube onto the floor and ceiling. (b) A PBRT reference solution using path-tracing with 8K rays/pixel computed in more than 40h.

Table 1 :

 1 Scene data and timings. T: number of input triangles, E: number of elements, L: number of links, N i : number of iterations; 128 bins are used. fps s : with splatting and 4x anti-aliasing. Shared vertex mesh elements are created at vertices yielding fewer samples.

	Scene	T	E	L	N i	fps	fps s
	Office	5.6K	18.2K	632K	4	10.1	9.0
	Character	8.1K	19.7K	661K	4	8.3	7.1
	Oriental	11.9k	15.6K	776K	4	7.5	6.8
	Glossy box	3.5k	4.6K	177K	3	14.9	11.9
	Glossy bust	13.4k	11.5K	412K	3	6.4	5.6

Table 2 :

 2 Timing and memory breakdown. IL is the time for initial linking. LP is the time (total, for all iterations) for the local pass, GP for the global pass, PP for push pull, M T is texture memory for I ex /L in tables and ping-pong textures and M L is the memory required for links (vertex buffers). Element attributes require less than 500Kb in all scenes.

	Scene	IL(s)	LP(ms)	GP(ms)	PP(ms)	M T (Mb)	M L (Mb)
	Office	8.7s	7.1	63.4	14.0	71.3	29.5
	Character	13.1s	8.2	72.0	17.0	77.0	31.1
	Oriental	13.0s	10.8	81.4	18.9	60.9	36.0
	Glossy box	2.4s	43.4	2.4	2.8	17.8	3.8
	Glossy bust	13.0s	100.9	5.6	6.6	44.9	8.8

Acknowledgements

The first author initiated this work at the University of Erlangen, and was partially funded by the DFG project "Interaktive Visualisierung Prozeduraler Modelle"; he also acknowledges support of the Marie-Curie Fellowship "Scalable-GlobIllum" (MEIF-CT-2006-041306). F. Durand acknowledges a Microsoft Research New Faculty Fellowship, a Sloan fellowship and an NSF CAREER award 0447561 "Transient Signal Processing for Realistic Imagery." We thank Autodesk for the generous donation of Maya, A. Olivier, F. Firsching and P. Richard for modeling help, and P. Green for suggesting the term Antiradiance.


Keywords: Spatial audio rendering, reverberation, geometrical acoustics, statistical acoustics

Topological Sound Propagation with Reverberation Graphs

Abstract

Reverberation graphs is a novel approach to estimate global soundpressure decay and auralize corresponding reverberation effects in interactive virtual environments. We use a 3D model to represent the geometry of the environment explicitly, and we subdivide it into a series of coupled spaces connected by portals. Off-line geometrical-acoustics techniques are used to precompute transport operators, which encode pressure decay characteristics within each space and between coupling interfaces.

At run-time, during an interactive simulation, we traverse the adjacency graph corresponding to the spatial subdivision of the environment. We combine transport operators along different sound propagation routes to estimate the pressure decay envelopes from sources to the listener. Our approach compares well with off-line geometrical techniques, but computes reverberation decay envelopes at interactive rates, ranging from 12 to 100 Hz.

We propose a scalable artificial reverberator that uses these decay envelopes to auralize reverberation effects, including room coupling. Our complete system can render as many as 30 simultaneous sources in large dynamic virtual environments.

Introduction

Modeling of sound propagation is very important for virtual acoustics and virtual reality applications. In particular, indoor reverberation effects due to sound scattering off wall surfaces carry major cues related to the size of the environment and distance to sound sources [START_REF] Begault | 3D Sound for Virtual Reality and Multimedia[END_REF]]. Therefore, reverberation helps users to establish a better sense of presence in virtual environments and thus is arguably one of the most important audio effects to simulate.

A number of geometrical approaches have been presented to model early specular reflections (and possibly edge diffraction) interactively [START_REF] Lauterbach | Interactive sound rendering in complex and dynamic scenes using frustum tracing[END_REF][START_REF] Lentz | Virtual reality system with integrated sound field simulation and reproduction[END_REF][START_REF] Funkhouser | A beam tracing method for interactive architectural acoustics[END_REF]]. Such approaches can handle complex geometry but are too computationally intensive to be applicable to interactive reverberation effects. Traditionally, interactive reverberation effects in single rooms are implemented by estimating energy decay rates and deriving artificial reverberation filters [START_REF] Jot | Real-time spatial processing of sounds for music, multimedia and interactive human-computer interfaces[END_REF]]. Parameters of the decay envelope can be directly manipulated by the sound designer to achieve the desired effect without requiring any geometrical modeling. Alternatively, the energy decay rates can also be determined from few geometrical parameters using well-known statistical models, such as the Sabine and Eyring formulas [START_REF] Kuttruff | Room Acoustics[END_REF]]. Such approaches are widely used in current audio production environments and consumer audio hardware [START_REF] Eax | Environmental audio extensions 4.0, Creative c[END_REF]], but they cannot model the global propagation of sound throughout a complex multiroom environment.

Consider the typical example of a game-like environment containing several connected rooms such as in Figure 1 (a). The late energy decay in the listener's room will be governed by the trans-mission of energy from the neighboring rooms through the open doors. The corresponding indirect energy contributions will also have strong directional components that need to be rendered correctly from the location of the coupling apertures. In the context of architectural acoustics, coupled-room configurations have been increasingly studied for the design of modern concert halls. While a preference for so-called "double-sloped" decays has not been established [START_REF] Bradley | The effects of simple coupled volume geometry on the objective and subjective results from nonexponential decay[END_REF][START_REF] Ermann | Double sloped decay: Subjective listening test to determine perceptibility and preference[END_REF]], they have been shown to be clearly audible under a variety of experimental conditions [START_REF] Picard | Audibility of non-exponential reverberation decays[END_REF][START_REF] Knight | Audibility of non-exponential energy decay in running reverberation[END_REF][START_REF] Ermann | Double sloped decay: Subjective listening test to determine perceptibility and preference[END_REF]] and thus are necessary to replicate for perceptually accurate reverberation effects.

In this paper, we concentrate on interactive, yet physically-based, modeling of reverberation effects in indoor architectural environments containing arbitrarily coupled rooms. Our approach combines characteristics of geometrical and statistical acoustics to achieve realistic results at interactive rates.

Starting with a 3D input model, an adjacency graph that describes the spatial coupling of acoustic spaces is constructed. Using a divide and conquer strategy, energy decay envelopes are then determined off-line for each space using path tracing simulations. This decouples complex geometrical calculations from global propagation, which is performed on-line by exploring all energy propagation routes along the spatial adjacency graph.

Similar to statistical acoustics, we assume that energy is scattered from space to space through the coupling apertures, which act as diffuse emitters [START_REF] Summers | Statisticalacoustics models of energy decay in systems of coupled rooms and their relation to geometrical acoustics[END_REF]]. We also assume that the fine-grain temporal structure of reverberation filters can be modeled as statistical noise. From the global propagation routes physicallybased artificial reverberation filters can be efficiently constructed and used for real-time spatial audio rendering.

Our contributions are:

• The definition of transport operators for general acoustically coupled rooms.

• An efficient global energy propagation algorithm based on the topology of the room adjacency graph.

• An efficient filter construction algorithm leveraging the coherence of sound propagation routes in the graph.

• A scalable interactive system for auralizing physically-based late reverberation effects for environments with many sources.

• An initial validation against alternative geometrical acoustics solutions.

Our complete pipeline allows for rendering up to 30 sources in complex dynamic environments.

Related Work

A large body of work has been devoted to simulating sound propagation effects in virtual environments. In the following section we review the techniques most directly related to our work.

Statistical acoustics

Statistical acoustics (SA) are widely used to determine energy decay rates in rooms and lead to efficient approaches to model reverberation effects. For instance, the Eyring model describes the energy decay in an enclosure as a decaying exponential:

E(t) = E o e cS 4V t log(1-ᾱ) , (1)
where c is the speed of sound (≈ 340m.s -1 ), V is the volume of the space, S the total surface area of the walls and ᾱ the surfaceaveraged absorption coefficient.

While such approaches can give good results in the case of single rooms, they fail to capture the multi-slope decays created when rooms are connected. In Figure 2, decays of two independent rooms are compared with a decay obtained by coupling the two spaces. Decays were computed using a ray-tracing approach to construct image sources [Vorländer 1989]. One million rays were used to compute a "reference" solution up to 200 reflection orders. Note that the slope of the decay is not linear in dB scale and is a combination of the decay slopes of the individual rooms. SA methods specifically designed to predict reverberation time in coupled volumes also exist. The earliest results for two coupled rooms are due to Davis [START_REF] Davis | Reverberation equations for two adjacent rooms connected by an incompletely sound-proof partition[END_REF]] and Eyring [START_REF] Eyring | Reverberation time measurements in coupled rooms[END_REF]]. Further details for the two-room case are given by Lyle [Lyle 1981a;Lyle 1981b] and Cremer and Müller [Cremer and Müller 1982], and Kuttruff extended the technique to an arbitrary number of coupled volumes [START_REF] Kuttruff | Room Acoustics[END_REF]]. All of these assume that the reverberant decays for each room in the coupled system are exponentially decaying and can be predicted independently. They also assume that modifications to the single-room decays occur due to diffuse energy exchange through portals between spaces. [START_REF] Summers | Statisticalacoustics models of energy decay in systems of coupled rooms and their relation to geometrical acoustics[END_REF]] provide an improved SA model for coupled volumes which uses corrections based on geometrical acoustics to increase the accuracy of single-room decay predictions as well as room-to-room energy exchanges.

Our approach can be seen as a generalization of previous SA approaches but also aims at providing a finer-grain estimation of the short-time integrated impulse response so that it can be directly auralized. A solution to more accurate modeling is to use methods based on geometrical acoustics.

Geometrical acoustics

Geometrical acoustics (GA) is probably the most popular approach for physically-based acoustic modeling. GA is a high-frequency approximation that models sound propagation along ray-paths. Specular reflection paths can be constructed using techniques such as ray or beam tracing [START_REF] Lauterbach | Interactive sound rendering in complex and dynamic scenes using frustum tracing[END_REF][START_REF] Lentz | Virtual reality system with integrated sound field simulation and reproduction[END_REF][START_REF] Funkhouser | A beam tracing method for interactive architectural acoustics[END_REF]]. Monte-Carlo path tracing or radiositylike approaches can be used for diffuse exchanges [START_REF] Heinz | Binaural room simulation based on an image source model with addition of statistical methods to include the diffuse sound scattering of walls and to predict the reverberant tail[END_REF]Dalenbäck 1996;[START_REF] Siltanen | The room acoustic rendering equation[END_REF]].

GA-based techniques have been previously used for the offline prediction of reverberation in coupled volumes. Using ODEON [Christensen ], a commercially available GA-modeling program, Bradley and Wang [START_REF] Bradley | Comparison of measured and computer-modeled objective parameters for an existing coupled volume concert hall[END_REF] evaluated the effects of various architectural parameters on the energy decay in a coupled-volume system. They also have shown good agreement at high frequencies between ODEON simulations and measurements from an existing coupled-volume concert hall [START_REF] Bradley | Comparison of measured and computer-modeled objective parameters for an existing coupled volume concert hall[END_REF]. [START_REF] Summers | Adapting a randomized beam-axis-tracing algorithm to modeling of coupled rooms via late-part ray tracing[END_REF] have shown that reverberant tails in coupled volumes can be predicted accurately with ray tracing in the software CATT [Dalenbäck ].

In a recent paper [START_REF] Lentz | Virtual reality system with integrated sound field simulation and reproduction[END_REF]], Lentz et al. report update times of up to 3 sec. to update full reverberation filters in a Monte-Carlo path tracing context. Early reflection visibility checking takes about 0.7 sec. for up to order 3 in their 105 polygon test model, leading to 300000 potentially visible image sources to test. While acceptable in the context of acoustics design, such update rates are limiting for typical consumer entertainment applications, such as games, where the simulation should typically run between 10 and 60Hz to accommodate fast listener or source motions and geometrical changes in the environments. Another paper [START_REF] Laine | Accelerated beam tracing algorithm[END_REF]], using a beamtracing approach, reports constructing all third order reflections of a sonund source in 1 msec. for a 1190 polygon model. However, even the fastest interactive GA approaches remain limited to low-orders of reflection (typically 8) and few sources. Hence, they cannot be applied to compute full reverberation effects interactively since reverberation filters typically contain audible reflections of order 100 or more.

To avoid the cost of an on-line GA simulation, Pope et al. [START_REF] Pope | Realtime room acoustics using ambisonics[END_REF]] have proposed to precompute and store the reverberation filters for a number of locations in a virtual environment. However, this approach does not scale well, since it would require computing a number of filters quadratic with the number of source/listener locations or rooms. It also does not support any interactive modification of the environment. In this paper however, we will retain the approach of pre-computing costly GA-based decay envelopes in an off-line process.

Artificial Reverberators

Statistical acoustics formulations lead to the development of efficient artificial reverberators, which are currently the most widespread techniques used to auralize late reverberation effects in interactive virtual environments [START_REF] Gardner | Reverberation Algorithms, Applications of Digital Signal Processing to Audio and Acoustics[END_REF][START_REF] Rocchesso | Spatial Effects, DAFX -Digital Audio Effects[END_REF]. Artificial reverberators do not model the fine-grain temporal structure of a reverberation filter but assume that reverberated components can be modeled as a temporal noise process modulated by slowly-varying pressure-amplitude envelopes in different frequency subbands. Following Eq 1, these envelopes are often considered as exponentially decaying, which leads to the design of efficient recursive Feedback Delay Network (FDN) filters [START_REF] Schroeder | Natural sounding artificial reverberation[END_REF][START_REF] Jot | Real-time spatial processing of sounds for music, multimedia and interactive human-computer interfaces[END_REF][START_REF] Gardner | Reverberation Algorithms, Applications of Digital Signal Processing to Audio and Acoustics[END_REF][START_REF] Rocchesso | Spatial Effects, DAFX -Digital Audio Effects[END_REF]. This approach of modeling the fine-grain structure of reverberation filters as a noise process has also been used to auralize Monte-Carlo path tracing results [START_REF] Kuttruff | Auralization of impulse responses modeled on the basis of ray-tracing results[END_REF]Schröder et al. 2007]. In this case, a short-time integrated impulse response is constructed via path-tracing and is used to modulate the noise process. In the general case, this approach cannot be directly implemented using the previous FDN formulation due to non-exponential decays, and results in a high computation cost. In this work, we will also be using this latter solution but show that it can be efficiently implemented within a scalable Fourier-domain framework [START_REF] Tsingos | Scalable perceptual mixing and filtering of audio signals using an augmented spectral representation[END_REF]].

Overview

Our approach requires a 3D model of the environment to simulate. In a first pre-processing step, we decompose the 3D environment into a cell-and-portal structure [START_REF] Teller | Visibility Computations in Densely Occuded Polyhedral Environments[END_REF][START_REF] Lefebvre | Automatic cell-and-portal decomposition[END_REF][START_REF] Haumont | Volumetric cell-and-portal generation[END_REF]]. We build an adjacency graph that encapsulates the physical connectivity of the individual building blocks of the geometry, i.e. the cells and portals. The nodes of the graph correspond to the rooms and the edges to the different coupling apertures (see Figure 1 (b)).

In contrast to previous GA work [START_REF] Lentz | Virtual reality system with integrated sound field simulation and reproduction[END_REF][START_REF] Funkhouser | A beam tracing method for interactive architectural acoustics[END_REF]], we do not use the cell-portal decomposition to speed up visibility checks. Rather, we use the topology of the corresponding adjacency graph to decouple the acoustic spaces and treat them independently. As such, our work shares some similarities with the approach concurrently suggested in [START_REF] Foale | Portal-based sound propagation for first-person computer games[END_REF].

By following this divide-and-conquer strategy, we seek to model reverberation as an energy transfer between rooms, using operators similar in spirit to form-factors [START_REF] Siltanen | The room acoustic rendering equation[END_REF]]. These operators encode not only direct propagation from room to room but also higher order reflections within rooms that encounter the portals and thus lead to additional, delayed energy exchange. For a more accurate modeling of the acoustical behaviour of each subspace, these operators are pre-computed for each individual room using a GA approach. Section 4 describes the operators and their off-line evaluation.

During an on-line simulation, global room-to-room propagation is solved by traversing the adjacency graph and constructing all possible propagation routes from all rooms to the current listening room (Figure 1 (c)). These routes will be re-used to auralize all audible sound sources in the listening room. Once sound-propagation routes have been determined, global short-time integrated pressure decays are computed for each source room. This involves convolving adjacent transport operators as energy radiates from subspace to subspace along each route. We propose to compute the required convolutions numerically in the Fourier domain. In addition, we propose a convolution caching mechanism that leverages the coherence between all propagation routes in the graph. This process is described in Section 5.

A set of matching artificial reverberation filters are then generated for auralization (Figure 1(d)). In Section 6, we show that efficient auralization of the reverberation effects can be achieved using a scalable rendering approach, similar to [START_REF] Tsingos | Scalable perceptual mixing and filtering of audio signals using an augmented spectral representation[END_REF]], which we extend to include spatial rendering of the reverberated components.

Finally in Section 7, we show that our approach compares well with full GA solutions and that our rendering engine can deal with complex environments at interactive rates.

Transport operators for coupled spaces

In this section, we introduce transport operators, or TOPs, describing global energy exchanges inside and in-between spaces. Transport operators model how energy decays locally within a single building block of the virtual environment due to reverberation. We also introduce operators modeling the response of incoming sound at the interface between each building block, i.e. portals. A specific pressure decay function is computed and associated with each transport operator.

Definitions

We define several transport operators for any given room in the environment. To handle general environments, we pre-compute transport operators using Monte-Carlo path tracing, see Figure 3 for an example. The Source-Room-Listener (SRL) operator is the point-to-point response of a source in a room for a listener in the same room (Figure 3 right). It is used for sources that directly share the same space as the listener. This operator typically can be obtained by summing up the contributions of image-sources identified by the path tracer [Vorländer 1989]. In our current implementation, SRL operators do not include the direct contributions to the listener, which we maintain separately for faster update rates (see Section 6).

The Source-Room-Portal (SRP) operator is the response of a point source in a room, integrated on a given portal in the same room. In our current implementation, we use a uniform solid angle sampling from the point source and compute the energy response of the SRP operator as:

SRP(t) = 4π nbrays nbrays ∑ i=1 test i ∏ k (1 -α k )δ (t -r/c),
where nbrays is the total number of rays fired and test i = 1 if the ray interects portal P and test i = 0 otherwise. c is the speed of sound, r the length of the path, δ the Dirac function. The α k 's are the surfaces' energy absorption coefficients along the path, with k being the reflection order. For better integration, russian roulette could alternatively be used to handle reflection coefficients.

The Portal-Room-Portal (PRP) operator is the response of a portal acting as a diffuse emitter integrated for all points on a portal of the same room (Figure 3 left). Note that the same portal can act as both a source and a receiver in a PRP operator, in which case only the indirect paths are sampled. To compute a PRP operator, we sample simultaneously the surface of the source portal P i and the hemisphere of shooting directions and obtain the energy at the destination portal P j as:

P i RP j (t) = 2π nbrays A i A j nbrays ∑ i=1 test i ∏ k (1 -α k )δ (t -r/c),
where A i and A j are the respective areas of the source and destination portals.

The Portal-Room-Listener (PRL) operator is the response of a portal acting as a diffuse emitter for a point receiver in the same room. In our current implementation PRL operators are computed in the same way as SRP operators.

Portals can be fully absorbent (e.g., open door), fully reflective (e.g., closed thick door) or a mixture of both (e.g., a thin wall between rooms). If the portals are fully absorbent, a ray hitting the portal is propagated no further.

We also assume that operators do not depend on the position of the source and listener within a given space and estimate them for a representative source and listener located at the centroid of each room. However, this is not a hard limitation since multiple operators could be precomputed for several source or receiver locations within each space at the expense of more processing and a larger memory footprint.

As a result, for a room with N P portals, N 2 P +2N P +1 operators must be determined, corresponding to N 2 P PRPs, N P PRLs, N P SRPs and one SRL. Note that this number is a function of the local number of portals in each room and is independent of the total number of rooms, portals or the complexity of the underlying geometry present in each room. Example operators are shown in Figure 4.

Pre-computing and storing transport operators

Our current implementation supports specular reflections, which we compute for up to 200 orders, but could be extended to include diffuse scattering as well for improved modeling [Vorländer 1995]. A regular grid was used to optimize ray-triangle intersections. Any alternative GA prediction approach supporting point and area receivers could be used. Each operator is computed for a number of frequencies to account for absorption and propagation effects. In our current implementation we used four chosen as 200, 1000, 8000 and 16000 Hz.

In the case of PRP and SRP operators, energy of all paths hitting the portal is integrated over short time-frames consistent with human hearing resolution (typically ∆T ≈ 10 ms).

For PRL and SRL operators, we use the rays to construct imagesources that we render at the listener's location [Vorländer 1989]. We also integrate the energy of the corresponding image-sources using the same ∆T. To limit aliasing we use a 25cm-radius detection sphere which is roughly consistent with the size of the listener's head. Aliasing-free approaches, such as pyramidal beam tracing, could alternatively be used.

From the obtained energy responses, we take the square-root to obtain a short-time averaged pressure response. Table 1 shows precomputation time for all operators in environments of varying complexity. As can be seen, the timings remain reasonable even for our highly conservative settings and pre-computation time largely depends on the topological complexity of the scene.

Operators are stored as short-time integrated pressure histograms (fixed size and ∆T), sampled at less than 200 Hz. We use a size of 1024 samples and ∆T ≈ 10 ms., which can encode decays up to 10 sec., and requires a storage of 4Kb per histogram assuming 32-bit floating point values. To account for frequency-dependent absorption effects, we store 4 echograms corresponding to different frequency subbands. Table 1 shows memory requirements to store the operators for environments of various complexities. As it can be seen, memory requirements remain reasonably low even for scenes with high polygon counts since transport operators only depend on the topological complexity of the scene.

In comparison to the brute force solution of storing all frequency responses for possible couples of source and listening rooms, this approach is more efficient, since we do not have to propagate through portals and thus can trace fewer rays. In addition, storing short time-integrated pressure decay envelopes, instead of full audio sampling-rate data (e.g., 44.1 KHz) reduces storage and memory footprint.

Combining operators

When sound propagates from room to room, we successively apply the transfer function of each operator along a desired sound propagation route. For example, Figure 1 (c) shows a decay envelope obtained by convolving all operators along a possible route. Combining operators is achieved by convolving their decay envelopes. An analytical solution for convolution of two exponential decay envelopes has been proposed in [START_REF] Svensson | Energy-time relations in a room with an electroacoustic system[END_REF]]. To handle more one convolution and support non-exponential decays, we compute the convolution numerically. Individual transport operators are pre-stored in the frequency domain so that fast convolutions can be achieved using only complex multiplies.

We further assume that all decay envelopes start at t = 0 and account for possible initial time-delay gaps separately. The time-delays of all operators convolved along a route can be factored out and added together leading to the total time-delay for the route. This delay can be subsequently applied to the result of the convolution. This limits aliasing problems, especially for small time-delay gaps.

We also need an estimation of energy contributed by a given propagation route. This energy is used as an oracle to direct graph traversal during the simulation of the global propagation and can also be used to prioritize rendering of propagation routes.

There is no general formula to estimate the energy of a convolution. However, an upper bound is given by:

|| f * g|| 2 ≤ min(|| f || 2 ||g|| 1 , || f || 1 ||g|| 2 ),
(2)

|| f * g|| 1 ≤ || f || 1 ||g|| 1 , (3)
where

|| f || 2 = ∑ i f 2 (i) and || f || 1 = ∑ i | f (i)|.
Given the starting energy of an SRP operator, we can successively apply Eq. 2 for successive PRP operators to estimate the energy carried along the route. Note that Eq. 3 is required to successively estimate the || f || 1 term in Eq. 2.

On-line topological sound propagation

Considering the initial cell-portal connectivity graph is already established and the transport operators for each room are computed, we describe in this section how global coupling is solved on-line by propagating energy through the graph and combining successive operators along all possible propagation routes. The pseudo-code of the complete algorithm can be seen in Algorithm 1.

Graph traversal

To model sound reverberation in the environment, we emit energy from the listener's room and, according to the connectivity graph, we select and apply the energy decays of respective transport operators to the propagated energy. In this way we approximate the exchanges that take place between the rooms of the environment.

The tree formed by selecting the listening room as the root node enables a fast, depth-limited, traversal algorithm to be applied. The criterion driving our traversal algorithm is that the propagated energy remains above a conservative energy threshold (see Eq. 2). In our simulations, we used a 10 -10 threshold (i.e., a 100 dB dynamic range).

The traversal algorithm takes advantage of the connectivity graph between the rooms and portals of the environment. Each time we enter a new room, we identify and recursively follow all possible adjacent nodes. At each step we update the propagating energy by applying the energy decay of the respective transport operator (Eq. 2). Followed routes are typically coded as transport operator sequences, e.g. P n R v P n-1 * . . . * [P i-1 R m-1 P i ] * P i R m L. Note that since operators already include all reflections within rooms, edges in the reverberation graph model coupling between adjacent spaces and always connect different nodes. As a result, two successive operators in a sequence always correspond to different rooms.

This traversal only extracts the longest possible routes that energy can follow from room to room until it decays below our minimal threshold, and is independent of the actual location of the actual sound sources that need to be auralized. To account for the potential presence of sound sources in any room along these routes, we simply consider that a unique virtual source is present in each room.

The initial extracted routes are further subdivided to construct all possible routes pairing these virtual sources with the listener, e.g. routes of the form SR v P n-1 * . . . * P i R m L.

Since the transport operators of a route will later be convolved to construct decay filters, the computational requirement of creating a filter along a route, made up of N T OPS number of TOPs, will be equal to performing N T OPS -1 convolutions. One way to further control the performance in an interactive system is to set a global budget and prune all routes that have a higher computational payload based on the estimated energy they carry.

Convolution caching

Once the significant propagation routes have been determined, we must calculate how energy is spread through time along each route. This involves combining the transport operators along a route, which corresponds to a set of convolutions as described in Sec. 4.3.

To estimate the global decay shape for all source rooms to the listening room, we cluster these routes first by source room, then by incoming portal in the listening room. The sum of all decay envelopes (SR v P * . . . * P i RL) must be computed for each portal P i in the listener's room and a source room R v .

Convolutions are computed in the Fourier domain. Operators' FFT data is pre-cached so that only the complex multiplies required for convolution are performed on-line. Since the number of samples used to encode the operators is small, this is efficient. Typically, 850 convolutions can be performed at 60 Hz on a Intel Xeon 2.33GHz processor using an unoptimized implementation. However, as the environments become more complex, convolving the operators can become the bottleneck of the approach.

A key aspect of formulating the energy propagation as a sequence of operators that can be convolved together is that we can take advantage of the associative and commutative properties of convolution to significantly reduce the number of convolutions actually performed. Recall that energy will sequentially propagate between the rooms of the environment in an iterative manner, thus in many routes the same operators may occur multiple times. The commutativity of convolution enables us to arbitrarily re-order transport operators, by using a pre-assigned unique numeric value that acts as an identification number of each transport operator. Furthermore, the associative property of convolution allows us to apply a simple dictionary-based compression scheme in which transport operators are convolved pair-wise and new transport operators, denoting this convolution, are stored in a cache. Therefore, when filters are created for each route, this convolution cache is queried for precomputed convolutions of operators. The high correlation of transport operators in the different sound propagation routes, emanating from the listener's room, provides a significant decrease in the number of required convolutions, as can be seen in Table 2. We obtained compression ratios of more than 10 to 1 in some of our test cases.

Handling changes in portal state

Our approach allows for modifying the topology of the graph interactively, for instance to model opening or closing doors in the environment. However, closing a door in a room will also affect transport operators which must then be recomputed. Our current solution to smoothly blend between portal states is to pre-compute each operator for all possible portal state (i.e., open/closed) in the room. We smoothly interpolate between states by linear combination of the transport operators. If a portal becomes fully opaque, we can remove the edge from the graph to speed up propagation.

Real-time auralization

Once our pressure-response envelopes have been computed, they can be used to auralize sound sources with consistent reverberation effects. For auralization, the decay envelopes are transformed back to time-domain and used as the envelope of a noise process, similar to artifical reverberators. To optimize the reverberation process, we use a progressive Fourier-domain approach, similar to [START_REF] Tsingos | Scalable perceptual mixing and filtering of audio signals using an augmented spectral representation[END_REF]], to perform block-based convolution of the input audio signals with the reverberation filters.

The next sections detail our spatial audio rendering model and our implementation of a global pipeline, shown in Figure 5, where energy propagation and spatial audio processing are performed asynchronously. A key feature of this pipeline is to decouple the propagation performed from room to room from the actual sources that need to be processed, which can be much larger.

Spatial audio rendering model

Consider a sound source in room R S that is audible in the listener's room R L . We auralize three components: 1) the direct sound to the listener, 2) the directional incoming contributions through each portal in the listener's room, and 3) the reverberated contributions including the reverberation effect of the listening room. Our current implementation uses binaural rendering.

The direct sound to the listener is rendered using positional audio processing.

Our current framework uses Head Related Transfer Functions (HRTF) from the LISTEN database (http://www.ircam.fr/listen). If the sound source is not in the listener's room, we use the shortest route from the graph propagation to derive an approximate diffracted component, which we render from the centroid of the portal, using a simplified Geometrical Theory of Diffraction formalism [START_REF] Tsingos | Modeling acoustics in virtual environments using the uniform theory of diffraction[END_REF]].

For each portal P i in the listening room R L , we render an additional directional source located at the centroid of P i . The source signal, in this case, is first convolved by a monophonic reverberation filter and then positioned using the corresponding HRTF. The decay envelope of the reverberation filter corresponds to the sum of all SR S P * PRP * . . . * PRP i routes ending at this portal. The obtained decay envelope is then converted back to time-domain to be used by our scalable artificial reverberation engine. This involves N P inverse Fourier transforms for the N P portals in the current listener's room.

Finally, we render spatially diffuse reverberated components for both the direct sound and the portal contributions. We convolve the signal of the source with a 2-channel random noise process, whose decay envelope includes the reverberation effect of the listener's room R L . For sources located in R L , this is directly a pre- For multiple sources in the same room the same decay envelopes are used, thus making our approach scalable. Another advantage of separating portal and diffuse contributions is that the relative level of the various components (portal, diffuse, room) can be individually adjusted by a sound designer which gives more control on the desired effect in the context of entertainment applications.

Implementation of the rendering pipeline

We implemented our real-time audio processing pipeline using two separate processes running asynchronously (see Figure 5).

The first process handles the energy propagation and construction of reverberation decay envelopes. If the listener changes room or the portal state is modified, we re-propagate energy in the environment from the room containing the listener (see Section 5).

In the second process, the artificial reverberation engine retrieves time-domain decay envelopes for each source and applies blockbased FIR convolution to the input sound signals. FIR convolution is done directly in the Fourier domain following the approach introduced in [START_REF] Tsingos | Scalable perceptual mixing and filtering of audio signals using an augmented spectral representation[END_REF]]. This approach performs lossy blockbased convolution, accounting for masking between reverberation blocks of all sources. Additionally, the engine dynamically allocates processing to blocks contributing more energy to the final mix according to a user-specified computational budget. Figure 7 illustrates how masking and budget allocation allows for controlling the computational load.

We process 1024 sample blocks at 44.1 KHz sampling-rate and use 50% overlap-add blending to avoid artefacts in the reconstruction. Hence, our audio processing thread runs at 86Hz outputing frames of 512 reconstructed samples.

To optimize processing, all our input signals are pre-stored in the Fourier frequency domain so that only a single inverse FFT is required to reconstruct an audio frame. To synthesize the reverberation filters, we use a circular array of successive white-noise blocks, which we also pre-transform into the Fourier domain. To render spatially diffuse reverberation components, we generate uncorre- lated noise sequences by picking a random starting block in the array for each output channel. Each block of noise is equalized in the Fourier domain using the values stored in the time-domain decay envelopes obtained after propagation. Hence, the synthesized impulse response corresponds to a noise sequence multiplied by a frequency-dependent staircase envelope with a staircase-length equal to our frame size (11 msec. in this case). To avoid audible artefacts, we use a 15-band equalization by interpolating the 4 frequency-dependent values stored in the envelopes at each timeframe (see Section 4.2). The 15 obtained gains are directly used to multiply the complex Fourier coefficients of the noise data within each band prior to the complex multiplication with the source signal. Figure 6 illustrates our reverberation process. Although we do not use zero-padding on the source signal (thus performing circular convolution), we found that the use of a Hann window limits audible aliasing. Overlap-add reconstruction also introduces timedomain smoothing and eliminates possible artefacts due to our discontinuous staircase envelope.

Results

We integrated our approach in the Penta G game engine (http://www.gebauz.com/index.php?page=penta-g) developed in the context of the EU IST Gametools project (http://www.gametools.org/).

A snapshot of the game is shown in Figure 1 (d).

Please visit http://wwwsop.inria.fr/reves/projects/revGraphs/ for additional material including demo sounds and videos.

Performance

As can be seen in Table 2, our topological sound propagation approach can compute global reverberation filters at interactive rates for environments exhibiting complex coupling structure. Using our approach, we were able to compute the contributions of 256 propagation routes, resulting in 1800 convolutions at 17 Hz. Although the compression ratio of our current caching algorithm is very good, performance does not increase sub-linearly due to the overhead of our current algorithm. However, we believe this could be significantly improved . Clustering of the obtained decay filters according to source room and destination portal and final time-domain reconstruction were achieved at 30 Hz in this environment. Inverse Fast Fourier Transforms (FFT) are only performed on-demand for rooms that contain the sources to auralize. One inverse FFT is required per active room and per portal in the current listener's room.

Our scalable reverberation engine runs at 86 Hz and is able to render about a 100 individual reverberation filters at each time-frame.

Figure 8 shows timings for each component of our scalable reverberation approach: masking, budget allocation and processing.

On a dual core Intel Xeon 2.33GHz workstation, our complete pipeline allows for rendering in real-time all reverberated components obtained from our propagation algorithm for 30 individual sources.

Comparison to path-tracing solutions

We compared our interactive propagation algorithm to an off-line image-source simulation implemented using a ray-casting process in different test environments. To generate these solutions 3 million rays were fired from the source and traced up to 200 orders of reflection. Valid image-sources were collected at the listening point using a 25 cm-radius detection sphere. simulation and our approach. As can be seen, both approaches lead to similar decay envelopes, although our results slightly overestimates the slope of the late decay.

Discussion and limitations

The fundamental limitation of our approach is that portals act as diffuse emitters. As a result, our approach does not accurately model early reflections. However, it would be possible to precompute a room-to-room matrix of possible early reflection paths or image-sources [START_REF] Lentz | Virtual reality system with integrated sound field simulation and reproduction[END_REF], and utilize our approach to render higher order reflections. Our current implementation only includes specular reflection and does not model diffuse reflection on surfaces, nor scattering/diffraction at portals which have been shown to significantly improve acoustical modeling [Vorländer 1995;[START_REF] Funkhouser | A beam tracing method for interactive architectural acoustics[END_REF]]. Our pre-computing step could be extended to include these effects. Another expected benefit of introducing diffuse reflection would be to reduce the error introduced by our "diffuse-portal" assumption by generating more diffuse exchanges between spaces. We also currently pre-computed operators for a single representative source and listener position in each space.

If the geometry of the space is complex, for instance if a room and a corridor have been merged into the same space, this might lead to unrealistic results. A solution is to define cell and portals so that the geometry of the space boundaries is convex which will limit this problem.

The performance of our approach and quality of the results is also directly influenced by the temporal sampling resolution of the decay envelopes. Sampling resolution is related to the modulation spectrum of the modeled transfer functions [START_REF] Schroeder | Modulation transfer function: Definition and measurement[END_REF]]. In particular, properly taking into account early reflections would require a finer resolution or preferably a separate modeling. Our choice of a 10 msec. integration window is also directly related to the processing frame size for our auralization engine (currently running at 83 Hz). For improved efficiency, we currently use the same DT (11.6 msec.) for both the length of an output audio frame and the sampling rate of the decay envelopes. Hence, we can directly obtain the decay value at each audio frame without re-sampling the envelopes.

Currently we render our reverberation effects as decorrelated noise signals that do not capture the fine-grain spatial structure of the reverberation within each room. We believe an approach simi- lar to the Spatial Impulse Response Rendering [START_REF] Merimaa | Spatial impulse response rendering[END_REF]] could be used to introduce finer grain modeling of the spatial aspects of the reverberation effect by pre-computing timedependent average incidence directions with SRL and PRL operators.

Our current work focuses on indoor architectural environments but we believe possible application to more general situations, such as outdoor scenes (cities, etc.). These raises several issues, for instance how to define cells/portals and their connectivity graph in this context. Our assumption of diffuse energy exchanges might also not be valid in this particular case.

Conclusion

We presented a novel graph-based topological sound propagation algorithm that can compute interactive reverberation effects in complex coupled environments. In our auralization system we use auditory masking and scalable Fourier domain processing to render a large number of reverberated components.

Our approach supports multiple moving sources and listener. It also allows on-line topology changes and some geometrical modifications and would be suitable for procedurally-generated environment, e.g. using pre-defined building blocks as in [START_REF] Merrell | Example-based model synthesis[END_REF]].

Results obtained using our interactive approach compare well to off-line geometrical acoustics solutions. In the future, the results of our technique could be improved by incorporating diffuse scattering, instead of pure specular reflections. While further numerical and perceptual validation is necessary, we believe the proposed method opens interesting perspectives for real-time auralization that could find applications both for acoustical design and virtual environments.

Figure 1 :

 1 Figure 1: (a) Our approach uses a 3D scene, decomposed into cells and portals, to compute reverberation effects interactively. (b) We estimate global pressure decays along possible sound propagation routes in the cell adjacency graph describing the 3D geometry. (c) The pressure decay for a route (highlighted in red in (b)) is obtained by convolving room-to-room transport operators computed off-line. (d) Our graph-based reverberation engine has been integrated into a 3D game with scalable reverberation processing and supports up to 30 sources interactively. Inset shows the spectrogram of a reconstructed reverberation filter.

Figure 2 :

 2 Figure 2: Integrated energy decays for two uncoupled and coupled rooms. The volume of room #1 is 29 m 3 while room #2 is much larger at 578 m 3 . Energy absorption is uniform on all wall surfaces and was set to 5%.The double decay slope due to coupling is clearly visible on the blue curve.

Figure 3 :

 3 Figure 3: Illustration of the path tracing process used to compute a PRP and SRL transport operator.

Figure 4 :

 4 Figure 4: Example transport operators pre-computed for the highlighted room.

Figure 5 :

 5 Figure 5: Overview of the rendering pipeline. Each box runs as an asynchronous process.

Figure 8 :Figure 9 :

 89 Figure8: Breakout times for the various components of our scalable reverberation engine for the same sequence as Figure7. Note that the masking and budget allocation introduce very little overhead.

Figure 10 :

 10 Figure 10: Comparison of our approach to an image-source simulation for different configurations of sources and listeners. Stars and spheres indicate source and listener positions, respectively. Graph colors correspond to <source,listener> pairs for which the simulations were computed. Darker colors correspond to the reference path tracing simulation while light colors are the result of our topological graph-based approach.

Table 1 :

 1 Precomputing

		Environment			Time	Memory
	3D Model	#cells(tris)	#portals	#TOPs	(mins)	(Kb)
		3(39)	2	17	6	272
		2(12729)	1	8	22	128
		14(195)	15	140	42	2240
		29(1614)	33	339	151	5424

transport operators. Statistics collected on an Intel Xeon 2.33GHz by tracing 300000 rays/operator, with 200 the maximum number of reflections per ray. Memory statistics are provided for operators with 4 frequency bands with a length of 1024 samples. We show color-coded cells, and portals are highlighted in yellow.

Table 2 :

 2 Topological sound propagation statistics, averaged over 600 runs on an Intel Xeon 2.33Ghz. computed SRL envelope. For portal contributions, we generate the corresponding decay envelope as the sum of all SR S P * PRP * . . . * PRP i * P i R L L routes. This only involves a single additional convolution of the previously computed directional source contribution with the final P i RL operator. The obtained envelope is also conback to time-domain before being transfered to the artificial reverberation engine.

	Traversal Time (msecs)	0.2	0.07		7.25	1.71
	Route Construction Time (msecs)	0.36	0.05		12.93	9.32
	#Routes	23		9			256	59
	Convolution Cache	off	on	off	on	off	on	off	on
	#Convolutions	123	22	43	10	1819	420	852	62
	Convolutions Time (msecs)	2.82	1.79	1.19	0.83	40.23	37.65	15.69	6.73
	Compression ratio	-	5.6:1	-	4.3:1	-	4.3:1	-	13.7:1
	Cache size (Kb)	-	432	-	208	-	7200	-	1152
	Total Time (msecs)	3.38	2.35	1.31	0.95	60.41	57.83	26.72	17.76

  Figures 9 and 10 show comparisons of pressure decay envelopes obtained using the image-source

	700 700
			# sources
	600 600		# frames
			# audible frames
		# FFT coeffs * 100
	500 500	# allocated coeffs * 100 budget (30000)
	400 400
	300 300
	200 200
	100 100
	0 0
	0 0	200 200	400 400	600 600	800	1000 1000	1200 1200
				frames
	Figure 7: for our scalable reverberation engine on a se-
	quence of 1200 audio processing frames where up to 20 sources
	are active. The plot displays the number of sources, number of
	frames of input audio data mixed for block-based convolution, and
	corresponding number of FFT coefficients before and after budget
	allocation (we set a global budget of 30000 coefficients in this ex-
	ample).
	seconds
	0.012
				audio frame length
	0.01					update masking
						alloc
					processing
	0.008				reconstruct total
	0.006
	0.004
	0.002
	0
	0	200	400	600	800	1000	1200
				frames

Acknowledgments

The authors would like to thank Fanouris Moraitis for modeling the test 3D environments, Michael Wimmer and Chris Chiu for making the Penta-G game engine available to us and providing support, David Alloza for initial discussions on graph-based reverber-ation and Xavier Granier. We also would like to thank the anonymous reviewers for their helpful suggestions. This research was funded by the EU FET Open project IST-014891-2 CROSSMOD (http://www.crossmod.org).