Differentiable Geometric Acoustic Path Tracing

using Time-Resolved Path Replay Backpropagation

Ugo Finnendahl $^{\dagger}$ , Markus Worchel $^{\dagger}$ , Tobias Jüterbock, Daniel Wujecki,
Fabian Brinkmann, Stefan Weinzierl, Marc Alexa

Computer Graphics Group

Audio Communication Group

$^{\dagger}$ : Equal contribution

Light

(acoustic) scene parameter: $\pi$

$\text{Loss}$

$\[[\]$

$\[]\]$

differentiable rendering

inverse rendering

image: $$I$$

gradient step

Input

Output

objective function

Sound

differentiable simulation

spectrogram: $$S$$

inverse simulation

gradient step

$\text{Loss}$

$\frac{\partial}{\partial \pi}\text{Loss}(S)$

$\[[\]$

$\[]\]$

Because similar to light rendering, where, given scene parameters, the renderer produces an image,
There is sound rendering, that, given an acoustic scene, renders a so called Spectrogram.

And this is heavily used in architectural acoustics where engineers try to find the suitable input parameters to archive a certain output, for example to reduce noise or archive better listening experiences in theaters. Which is analog to inverse rendering.

But this inverse simulation is largely a manual task so far, where engineers tweak the parameters by hand and rerun the simulation.

So our goal of this paper was to make this process automatic like we do in graphics
by making the acoustic simulation differentiable with respect to the scene parameters so we can optimize an objective function by applying gradient steps on the input.

So how does an acoustic simulation work?

There are different ways on how to simulate sound but one common way,

wavelength

time

energy

wavelength

energy

time

energy

$\text{Histogram}(t)$

(R,G,B)

per ray

independent

dimensions

✔ spatial (2D) ✘

✘ time (1D) ✔

✔ wavelength (1D) ✔

pixels

no pixels

There are different ways on how to simulate sound but one common way,
used by a lot of commercial software,
is to ray trace sound energy.

This might be new to you but we can use the same gerneralized model for sound and light.

We have emitters that induce waves based energy
and the energy propagation can be simulated using rays that bounce around in the scene and get reflected according to a BSDF
until they eventually get measured by a sensor.

If we take a continuous perspective, the sensor receives energy over different wave lengths over time.
And this for every point on its surface.

And here is the first difference, the sensor in sound rendering is usually a sphere instead of a plane segment.

But in both cases the sensor receives energy over a continuous four dimensional space.
Two spacial dimensions, one time dimension and one wavelength dimension.

In light rendering we usually assume that light travels instantanously so we kind of integrate over time
and we also discretize the spatial dimensions by using pixels on the sensor.

whereas in sound rendering, time is very important, as sound waves travel much slower than light waves.
But on the otherhand we usually integrate over the spatial dimension so we do not have pixels on the sensor.

For the wavelength dimension we reduce the complexity in both cases, by representing it in a lower dimensional space.

Which results in RGB values for light and the already shown Spectrogram for sound. Usually displayed as heat map.

As the basis elements of the spectrum are considered independent, we can simulate them independently, (click)
revealing the major difference between light rendering and acoustic rendering.
In light rendering each ray collects scalar intensities whereas in sound rendering each ray
collects a histogram that is energy over time, so a single row of a spectrogram.

So if we formalize this we get two rendering equations that differ in two aspects.

$E_{\out}$

$E_{\inp}\frag{40}{ = E_{\out}^{\text{next}}}$

$$f$$

$E_{\emit}$

$\htmlId{leq}{\fragclass{1}{custom gray}{\fragclass{65}{bluebox box custom}{E_{\out}} = } \fragclass{3}{custom gray}{\fragclass{60}{redbox box custom}{E_{\emit}} +} \fragclass{7}{custom gray}{\int_{S}}\fragclass{6}{custom gray}{\fragclass{61}{redbox box custom}{f}}\,\, \fragclass{5}{custom gray}{\fragclass{65}{bluebox box custom}{E_{\del{45}{\inp}\add{45}{\out}}^{\add{45}{\text{next}}}}} \, \fragclass{7}{custom gray}{d\solid'^\bot}}$

light rendering equation

$E_{\out}$

$E_{\out}^{\text{next}}$

$$f$$

$E_{\emit}$

time

$E_{\inp}\frag{40}{\neq E_{\out}^{\text{next}}}$

$\flip{15}{leq}{\fragclass{65}{bluebox box custom}{E_{\out}\add{20}{(t)}}\, = \fragclass{60}{redbox box custom}{E_{\emit}\add{20}{(t)}} +\int_{S}\fragclass{61}{redbox box custom}{f}\,\, \fragclass{65}{bluebox box custom}{E_{\del{45}{\inp}\add{45}{\out}}^{\add{45}{\text{next}}}\add{20}{(t\del{45}{)}\add{45}{\,-\tau)}}} \, d\solid'^\bot}$

acoustic rendering equation

Light

Sound

local

recursive

time between

two points

(R,G,B)

$\in$

$\R^n$

$\fragclass{35}{fade-out}{\mathcal{F}}$

$\in$

So if we formalize this we get two rendering equations that differ in two aspects.

The light rendering equation should be familiar. click click click click click

Where as the acoustic rendering equation needs to consider time,
which results in an additional parameter t.

So each ray contributes a function to the total energy, instead of a RGB vector.
This function is usually discretized into time bins so it can be stored in a vector as well, but in a much higher dimension.

And the second difference is that incoming energy is not equal to the outgoing energy like it was in light rendering.
We need to consider the time the energy takes from one point to the next intersection point.
So the incoming energy curve needs to be translated by the time it takes for the sound to reach the next intersection.

Or in other words if an emitter has a certain emitter curve the path length defines when this emitter curve is received at the sensor.

But the structure shows that this can still be simulated in a single for loop
as in both cases we have
a local term plus the integral of a local term times a recursive term.

We just need to track the path length in the acoustic setting.

So the Histogram measured at the sensor is just the integral of the incoming energies. Or if we use Monte Carlo Simulation the sum of the energies collected by the sampled rays.

$\htmlId{histt}{\text{Histogram}(t)}=\int_S E_{\out}(t) \add{3}{\approx}\add{3}{\sum E_{\out}(t)}$

AD( $\text{Loss}(\text{Histogram}(t))$ ) $\rightarrow$ $\htmlData{fragment-index=4,class=custom added delayed}{\frac{\partial}{\partial\pi}\text{Loss}(}\flip{4}{histt}{\text{Histogram}(t)}\htmlData{fragment-index=4,class=custom added delayed}{)}$

AD Graph

too much memory

So the Histogram measured at the sensor is just the integral of the incoming energies. Or if we use Monte Carlo Simulation the sum of the energies collected by the sampled rays.

So now we are interessted in the gradient of a Loss function on this histogram, with respect to an arbitrary scene parameter $\pi$ , for example an absorption factor.

The straight forward way would be to use automatic differentiation.
And although automatic differentiation of a Monte Carlo simulation does not always results in the correct gradients, we assume for this presentation that it does. Check the paper for details.

And automatic differention works, but the AD graph that records all the operation grows with the depth of the simulation path
and it quickly overgrowed the storage space of our GPU.
While the path depth in light rendering is usually low, for acoustic scenes more bounces are relevant for the characteristics of a room.

For that reason we need to compute the gradient directly instead of relying on AD.

In the past it was shown that this is surprisingly easy in light rendering using so called "Path Replay Backpropagation" or PRB.

Light

Sound

Primal Simulation

$\htmlId{PRBEout}{\fragclass{60}{bluebox box custom}{E_{\out}}} \htmlId{PRBEouteqsign}{=} \htmlId{PRBEemit}{\fragclass{40}{redbox box custom}{E_{\emit}}} \htmlId{PRBplusint}{+ \int_S} \htmlId{PRBf}{\fragclass{50}{redbox box custom}{f}}\, \htmlId{PRBEoutnext}{\fragclass{60}{bluebox box custom}{E_{\out}^{\text{next}}}} \,\htmlId{PRBdsolid}{\text{d}\solid'}$

$\bluebox{110}{\add{230}{\deltaH}E_{\out}(t)} = \redbox{110}{\add{230}{\deltaH}E_{\emit}(t)} + \int_S \redbox{110}{f(t)}\, \bluebox{110}{\add{230}{\deltaH}E_{\out}^{\text{next}}(t-\tau)} \,\text{d}\solid'$

$\greenbox{140}{\add{230}{\deltaH}\frac{\partial}{\partial t}E_{\out}(t)} = \redbox{140}{\add{230}{\deltaH}\frac{\partial}{\partial t}E_{\emit}(t)} + \int_S \redbox{140}{f(t)}\, \greenbox{140}{\add{230}{\deltaH}\frac{\partial}{\partial t}E_{\out}^{\text{next}}(t-\tau)} \,\text{d}\solid'$

Derivative

$\begin{alignat*}{6}\fragclass{60}{orangebox box custom}{\add{30}{\frac{\partial}{\partial \pi}} \flip{20}{PRBEout}{E_{\out}}} \flip{20}{PRBEouteqsign}{=} \fragclass{40}{redbox box custom}{\add{30}{\frac{\partial}{\partial \pi}}\flip{20}{PRBEemit}{E_{\emit}}} \flip{20}{PRBplusint}{+ \int_S}&&\flip{20}{PRBf}{\fragclass{50}{redbox box custom}{\htmlId{PRBf2}{f}}}\, &&\fragclass{60}{orangebox box custom}{\add{30}{\frac{\partial}{\partial \pi}}\flip{20}{PRBEoutnext}{\htmlId{PRBEoutnext2}{E_{\out}^{\text{next}}}}}&\add{30}{+}\,\fragclass{30}{fade-out instant}{\flip{20}{PRBdsolid}{\htmlId{PRBdsolid2}{\text{d}\solid'}}}\\&&\fragclass{70}{redbox box custom}{\add{30}{\frac{\partial}{\partial \pi}\flip{30}{PRBf2}{f}}}\, &&\flip{30}{PRBEoutnext2}{\fragclass{80}{bluebox box custom}{E_{\out}^{\text{next}}}}&\,\flip{30}{PRBdsolid2}{\text{d}\solid'}\end{alignat*}$

$\begin{alignat*}{6}\orangebox{110}{\frac{\partial}{\partial \pi} E_{\out}(t)} = \redbox{110}{\frac{\partial}{\partial \pi}E_{\emit}(t)} + \int_S&&\redbox{110}{f(t)}\, &&\orangebox{110}{\frac{\partial}{\partial \pi}E_{\out}^{\text{next}}(t-\tau)}&+\\&&\redbox{110}{\frac{\partial}{\partial \pi}f(t)}\, &&\bluebox{110}{E_{\out}^{\text{next}}(t-\tau)}&-\\&&\redbox{135}{f(t)}\, &&\greenbox{137}{\frac{\partial}{\partial t}E_{\out}^{\text{next}}(t-\tau)}& \redbox{135}{\frac{\partial}{\partial \pi}\tau}\,\text{d}\solid'\end{alignat*}$

local

$\R^n$

$\in$

$\R^{\del{235}{n}}$

$\in$

$\R^{\del{235}{n}}$

$\in$

$\frac{\partial}{\partial \pi}\text{Loss}(\text{H}(t))\add{225}{= [\deltaH\text{Loss}]\frac{\partial}{\partial \pi}\text{H}(t)}$

scalar

linear time &
constant memory

In the past it was shown that this is surprisingly easy in light rendering using so called "Path Replay Backpropagation" or PRB.

Basically they took the derivative of the rendering equation and saw that the gradient can be computed using a similar simulation to normal rendering.

beacause the derivative reveals a familiar pattern
A local term plus the integral of a local term times a recursive term.
There is just one additional term in the integral that is new due to the product rule.
But this is just a local term times the result of the primal evaluation.

Path Replay Backpropagation showed that this incoming energy from the primal simulation, can be precomputed using a forward simulation and then retrieved in a second simulation for the derivative using constant memory.

Therefore the derivative is computable in linear time and constant memory with respect to the path depth.

So we did the same with the acoustic rendering equation
and you can see the structure of the derivative is pretty similar,
but we have a third term in the integral introduced by time.

Remmeber, moving geometry might change the length of a path and if the path length changes, all incoming energy is shifted,
therefore we have a contribution to the gradient that captures how much pi influences the path length and how much the path length influences the incoming energy.

But this term consists of two terms that can be evaluated locally and another term that follows a familiar pattern once again.

I hope that this is convincing enough that you believe me that this aditional term can also be precomputed like the incoming radiance in a primal simulation and then retrieved in a second simulation for the derivative.

So temporal path replay backpropagation seems straight forward.

Unfortunately there is one slight problem. Remember, we work on function spaces or high dimensional vectors.
That is a problem as we need to precompute and store the green and blue highlighted function for each ray
which would take a lot of storage per ray.

I will not explain how we solved that in detail
but as we want to backpropagate the gradient of a loss function, we need to evaluate a dot product due to the chain rule.
And luckily, as the derivative comes down to a bunch of sums, we can build the dot product wherever we need.
So we do it as early as possible and therefore reduced the necessary storage to two scalars.

I know this is a bit handwavy and I skipped a lot of things, but its the intuition and you can find the full story in the paper.

We implemented our prototype on top of Mitsuba3 and our approach does in fact only need constant memory, is faster than AD and returns the same results.

	AD		Temporal PRB
depth	time	memory	time	memory
4	0.730s	13.5MiB	0.213s	6.75MiB
8	1.873s	27MiB	0.343s	6.75MiB
16	7.111s	54MiB	0.616s	6.75MiB
32	17.178s	108MiB	1.246s	6.75MiB
64	DNF	DNF	2.699s	6.75MiB

faster & constant memory

Material Properties

absorption factors (visualized)

optimization

Emitter Position

Geometry

lecture hall

side view

No Wave effects

Diffraction
Interference

Biased Gradients

(if geometry moves)

1. First Differentiable Geometric Acoustics Path Tracer

2. Temporal extension to Path Replay Backpropagation

Differentiable Time-Gated Rendering [Wu et al. 2021]

Differentiable Transient Rendering [Yi et al. 2021]

Jakob, W., Speierer, S., Roussel, N., et al. 2022. Mitsuba 3 renderer. Software.

Worchel, M., Finnendahl, U., and Alexa, M. 2025. Radiative Backpropagation with Non-Static Geometry. Eurographics Symposium on Rendering, The Eurographics Association.

Wu, L., Cai, G., Ramamoorthi, R., and Zhao, S. 2021. Differentiable time-gated rendering. ACM Trans. Graph. 40, 6.

Yi, S., Kim, D., Choi, K., Jarabo, A., Gutierrez, D., and Kim, M.H. 2021. Differentiable transient rendering. ACM Trans. Graph. 40, 6.

Differentiable Geometric Acoustic Path Tracing

using Time-Resolved Path Replay Backpropagation

Ugo Finnendahl\(^{\dagger}\), Markus Worchel\(^{\dagger}\), Tobias Jüterbock, Daniel Wujecki,Fabian Brinkmann, Stefan Weinzierl, Marc Alexa

\(^{\dagger}\): Equal contribution

Differentiable (Geometric Acoustic) Rendering

Light

Input

Output

Sound

Energy Path Tracing

\(\text{Histogram}(t)\)

(R,G,B)

dimensions

Rendering Equations

Light

Sound

Automatic Differentiation (AD)

too much memory

Temporal Path Replay Backpropagation (TPRB) [Vicini et al. 2021]

Light

Sound

Primal Simulation

Derivative

local

linear time & constant memory

Implementation

TPRB (our)

AD

spectrogram

scene

gradients

faster & constant memory

Results

Material Properties

Emitter Position

Geometry

Limitations

No Wave effects

Biased Gradients

Summary

Thank You!

References

Test

Test

Ugo Finnendahl $^{\dagger}$ , Markus Worchel $^{\dagger}$ , Tobias Jüterbock, Daniel Wujecki,
Fabian Brinkmann, Stefan Weinzierl, Marc Alexa

$^{\dagger}$ : Equal contribution

$\text{Histogram}(t)$

linear time &
constant memory