Massively Parallel Path Space Filtering

Nikolaus Binder, Sascha Fricke, and Alexander Keller

arXiv:1902.05942v2 [[Link]] 3 Feb 2021

Abstract Restricting path tracing to a small number of paths per pixel for per-
formance reasons rarely achieves a satisfactory image quality for scenes of interest.
However, path space filtering may dramatically improve the visual quality by sharing
information across vertices of paths classified as proximate. Unlike screen space-
based approaches, these paths neither need to be present on the screen, nor is filtering
restricted to the first intersection with the scene. While searching proximate vertices
had been more expensive than filtering in screen space, we greatly improve over this
performance penalty by storing, updating, and looking up the required information
in a hash table. The keys are constructed from jittered and quantized information,
such that only a single query very likely replaces costly neighborhood searches. A
massively parallel implementation of the algorithm is demonstrated on a graphics
processing unit (GPU).

1 Introduction

Realistic image synthesis consists of high-dimensional numerical integration of

functions with potentially high variance. Restricting the number of samples therefore
often results in visible noise, which efficiently can be reduced by path space filtering
[KDB16] as shown in Figure 1.
We improve the performance of path space filtering by replacing costly neighbor-
hood search with averages of clusters in voxels resulting from quantization. Our new
algorithm is suitable for interactive and even real-time rendering and it enables many
applications trading a controllable bias for a dramatic speedup and noise reduction.

Nikolaus Binder (NVIDIA), Sascha Fricke (University of Braunschweig), and Alexander Keller
(NVIDIA)

1
2 Nikolaus Binder, Sascha Fricke, and Alexander Keller

Fig. 1 Path tracing at one path per pixel (top) in combination with hashed path space filtering
(middle) very closely approximates the reference solution using 1024 paths per pixel (bottom) and
does so with an overhead of about 1.5 ms in HD resolution. Scene courtesy of Epic Games.

2 Light Transport Simulation

Light transport is simulated by tracing rays that connect the surfaces of the light
sources and the camera sensor through a three-dimensional scene, represented by
surfaces, materials, and scattering properties of the materials and in volume. Upon
interaction with the surfaces or participating media in the volume, a new segment
of the path is generated. The connection between adjacent path segments is called
vertex, and all information of its path and the interaction is referred to as its descriptor.
While in reality paths of photons start on emissive surfaces, and the camera sensor
measures the ones terminating on its surface, simulations most often construct paths
backwards. Somewhat counter-intuitively, a simulation in reverse photon direction
is called “forward path tracing” – forward in view direction.
Massively Parallel Path Space Filtering 3

𝐿𝑒 𝐿𝑒
𝑦 𝑦

𝜃𝑥 𝜃𝑦 𝜃 𝑥 0𝜃
𝑦
𝑃 𝑃
𝑥 𝑥𝑥 0

(a) forward path tracing, subpath (b) density estimation and

connections, and next event estimation path space filtering

Fig. 2 Geometry of light transport simulation by path tracing. Starting from the eye on the left
through the image plane 𝑃, a light transport path segment ends in 𝑥. (a) Forward path tracing
continues the path until the path terminates on the surface of a light source. A path can also be
completed by directly connecting to a vertex 𝑦 on the surface of a light source (next event estimation)
or to a vertex 𝑦 of the same or a different path (subpath connection). (b) Radiance in point 𝑥 can
directly be evaluated by accumulating radiance in vertices 𝑥 0 in a local neighborhood either for
density estimation or path space filtering.

At each point 𝑥 and direction 𝜔, the incident radiance 𝐿 𝑖 is the sum of the emitted
radiance 𝐿 𝑒 and the reflected radiance 𝐿 𝑟 in the point 𝑥 and in direction 𝜔:

𝐿 𝑖 (𝑥, 𝜔) = 𝐿 𝑒 (𝑥, 𝜔) + 𝐿 𝑟 (𝑥, 𝜔) (1)

Given the distance 𝑑 := k𝑐 − 𝑥k 2 between the center of the ball 𝑐 and a point 𝑥, the
characteristic function of the ball 𝐵 with radius 𝑟 is defined as
(
1 𝑑2 < 𝑟 2
𝜒𝐵 (𝑑, 𝑟) := (2)
0 otherwise.

The following eqs. (3) to (6) show four different ways to formulate the reflected
radiance 𝐿 𝑟 in a point 𝑥 in direction 𝜔𝑟 and fig. 2 illustrates the resulting techniques:

𝐿 𝑟 (𝑥, 𝜔𝑟 )
∫
= 𝐿 𝑖 (𝑥, 𝜔) 𝑓𝑟 (𝜔𝑟 , 𝑥, 𝜔) cos 𝜃 𝑥 𝑑𝜔 (3)
S−2 ( 𝑥)
cos 𝜃 𝑦
∫
= 𝑉 (𝑥, 𝑦)𝐿 𝑖 (𝑥, 𝜔) 𝑓𝑟 (𝜔𝑟 , 𝑥, 𝜔) cos 𝜃 𝑥 𝑑𝑦 (4)
𝜕𝑉 |𝑥 − 𝑦| 2


𝜒𝐵 𝑥 − ℎ(𝑦, 𝜔), 𝑟 (𝑥)
∫ ∫


lim 2
𝐿 𝑖 ℎ(𝑦, 𝜔), 𝜔 ·
= 𝑟 ( 𝑥)→0 𝜕𝑉 S−2 ( 𝑦) 𝜋𝑟 (𝑥) (5)


· 𝑓𝑟 𝜔𝑟 , ℎ(𝑦, 𝜔), 𝜔 cos 𝜃 𝑦 𝑑𝜔𝑑𝑦
∫
𝜒𝐵 𝑥 − 𝑥 0, 𝑟 (𝑥) 𝑤(𝑥, 𝑥 0)𝐿 𝑖 (𝑥 0, 𝜔) 𝑓𝑟 (𝜔𝑟 , 𝑥, 𝜔) cos 𝜃 𝑥0 𝑑𝑥 0
∫

𝜕𝑉
= lim ∫
𝑑𝜔
𝑟 ( 𝑥)→0 S−2 ( 𝑥) 𝜕𝑉 𝐵
𝜒 𝑥 − 𝑥 0, 𝑟 (𝑥) 𝑤(𝑥, 𝑥 0)𝑑𝑥 0
(6)
4 Nikolaus Binder, Sascha Fricke, and Alexander Keller

Forward Path Tracing: Equation (3) integrates radiance over the upper hemi-
sphere 𝑆−2 by multiplying the incident radiance 𝐿 𝑖 from the angle 𝜔 in the point 𝑥
with the spatio-directional reflectivity 𝑓𝑟 for the two angles in the point 𝑥 and the
cosine of the incident ray to account for the change of area of the projected solid
angle. Inserting the equation into eq. (1) and the result back into eq. (3) allows for
subsequently prolonging paths and is known as (recursive) forward path tracing.
Next Event Estimation and Subpath Connection: Equation (4) changes the in-
tegration domain to the scene surface 𝜕𝑉 so that a path can be constructed in
which the point 𝑥 connects to any point 𝑦 on the scene surfaces. The integrand is
then extended by the mutual visibility 𝑉 of the two points 𝑥 and 𝑦. The fraction of
the cosine of the second angle and the squared distance of the two points accounts
for the change of measure. One often refers to this fraction as the geometric term.
Equation (4) is especially useful since it allows to directly connect to the surface
of a light source (next event estimation) or any other vertex of any path (subpath
connection).
Density Estimation: Equation (5) again integrates over the scene surface. How-
ever, it realizes density estimation by restricting to a local neighborhood in a
sphere with radius 𝑟 (𝑥) using the characteristic function of the ball 𝜒𝐵 . The den-
sity is then obtained by dividing by the area of the circle that stems from the
intersection of the ball and the flat surface. For a radius going to zero, the formu-
lation is equivalent. Density estimation tracing photons from the light sources is
usually referred to as photon mapping.
Path Space Filtering: Equation (6) similarly integrates over a local neighborhood
of 𝑥. In contrast to density estimation, a local weighted average is calculated. This
local average is normalized by the integral of all weights in the neighborhood
instead of the area of a circle. While for a sphere with vanishing radius this
formulation is again equivalent, the method trades a certain bias for a reduction
of variance for any non-zero radius. Due to the filtering of the local average this
technique has been introduced as path space filtering [KDB16].
While implementations of eqs. (3) to (6) each individually come with their own
strengths and weaknesses and are therefore often combined for robustness in offline
rendering applications [Vea97, LW98, GKDS12], time constraints of real-time image
synthesis as well as advances in path guiding [DK18, MGN17, MMR+ 19] lead to
the vast majority of implementations only employing eqs. (3) and (4). Our work aims
at a substantial acceleration of path space filtering so that it can be used in real-time
applications, too, resulting in a considerable variance reduction.

2.1 Previous Work

Filtering results of light transport simulation is gaining more and more attention
in real-time, interactive, and even offline rendering. The surveys by Zwicker et al.
[ZJL+ 15] and Sen et al. [SZR+ 15] present an overview of recent developments.
The fastest approaches use only information available at primary intersections and
Massively Parallel Path Space Filtering 5

perform filtering in screen space. Further recent work is based on deep neural
networks [BVM+ 17, CKS+ 17], hierarchical filtering with weights based on estimated
variance in screen space [SKW+ 17], or on improving performance by simplifying
the overall procedure [MMBJ17].
Path space filtering [KDB16], on the other hand, averages contributions of light
transport paths in path space, which allows for filtering at non-primary intersections
and for a more efficient handling of dis-occlusions during temporal filtering. How-
ever, querying the contributions in path space so far had been significantly more
expensive than filtering the contributions of neighboring pixels in screen space.
Neglecting the fact that locations that are close in path space are not necessarily
adjacent in screen space enables interactive filtering in screen space [GDW+ 14]. As
a consequence, filtering is almost only a good approximation for primary rays or re-
flections from sufficiently smooth and flat surfaces. In fact, such filtering algorithms
are a variant of a bilateral filtering using path space proximity to determine weights.
Sharing information across pixels according to a similarity measure dates back as
early as the 1990s [Kel98]. Since then, several variants have been introduced, for
example by re-using paths in nearby pixels [BSH02], for filtering by anisotropic
diffusion [McC02] or using edge-avoiding Á-Trous wavelets [Dam11].
Fast filtering is also possible in texture space [MHC+ 16], which at least requires
a bijection between the scene surface and texture space. While this may be tricky
already, issues may arise along discontinuities of a parametrization in addition.
Furthermore, only filtering is restricted to locations on surfaces when operating in
texture space, and thus volumetric effects must be filtered separately.
Hachisuka et al. also use a hash table in a light transport simulation on the
GPU [HJ10]. The approach is fundamentally different in two aspects: First, it traces
photons and stores them in the hash table for density estimation, while our method
averages radiance in vertices from arbitrary light paths. Second, their method imple-
ments simple sampling with replacement in voxels, culling all but one photon per
voxel. Our method does the exact opposite: It collects radiance from all paths whose
vertices coincide in a voxel.
Mara et al. summarize and evaluate a number of methods for photon mapping on
the GPU [MLM13]. Their evaluation also includes work from Ma and McCool using
hash tables with lists of photons in per voxel [MM02]. While all examined methods
may be used for path space filtering instead of photon mapping, their performance
is at least limited by the maintenance of lists.
Havran et al. use two trees for final gathering with photon mapping [HHS05].
The overhead of tree construction and traversal as well as iterating through lists of
vertices severely limit the performance in our intended real-time use case.
Multiple Importance Sampling weights, for example those for path space filtering,
can be further optimized [WGGH20].
Kontkanen et al. explore irradiance filtering, a subset of path space filter-
ing [KRK04]. Spatial caching of shading results in a hash table for walkthroughs of
static scenes [DS07] uses similar methods to the ones presented in this work for the
lookup of these results. Again, the method can be seen as a subset of Path Space
6 Nikolaus Binder, Sascha Fricke, and Alexander Keller

Filtering: It is restricted to caching diffuse illumination and neither includes filtering

nor spatial and temporal integration in an arbitrary number of vertices of a light path.

3 Algorithm

Similar to path space filtering [KDB16], the input of the method is a set of vertices in
which radiance should be filtered. In addition, the method receives information from
their paths, such as attenuation and throughput. Therefore, the overall simulation must
run in at least two phases: The first phase generates paths and stores all necessary
information, similar to a path tracing simulation without filtering.
After a set of light transport paths has been generated, the second phase of the
algorithm filters radiance in voxels, as described in section 3.1, and stores as well
as looks up the averages in a hash table, see section 3.2. Techniques described in
section 3.3 improve the robustness of the algorithm in voxels with a small number
of vertices. Section 3.4 discusses additional filtering opportunities and challenges in
interactive light transport simulations.
Finally, after filtering, for each selected vertex its associated average is multiplied
by its throughput and accumulated in its respective pixel. The throughput is the
attenuation from the camera along the light transport path up to the selected vertex.

3.1 Averaging in Voxels

The new algorithm builds upon a different characteristic function to be used in

eq. (6). Unlike the characteristic function of a three-dimensional ball in eq. (2), it
uses the key 𝑘 constructed from the descriptor of a vertex instead of only depending
on its position 𝑥.
Section 3.1.1 defines this characteristic function, and section 3.1.2 describes the
construction of a key 𝑘 from the descriptor of a vertex in detail. Section 3.1.3 and
section 3.1.4 present extensions to the method that reduce the perceived bias.

3.1.1 Characteristic Function of a Voxel

Given a resolution selection function 𝑠(𝑘), the characteristic function of the voxel is
defined as
(
0 1 b𝑠(𝑘)𝑘c = b𝑠(𝑘 0)𝑘 0c ∧ 𝑠(𝑘) = 𝑠(𝑘 0)
𝜒𝑉 (𝑘, 𝑘 ) := (7)
0 otherwise.

In the simplest case, 𝑠(𝑘) is a constant which defines the size of all voxels. In
practice, other heuristics improve the quality of the approximation. For example, it
Massively Parallel Path Space Filtering 7

is often advisable to increase the world space extent of a voxel with its distance to
the camera sensor along the path. Then, one can filter more aggressively in distant
voxels, and increasing the size counteracts the decrease of the density of vertices
from paths directly coming from the camera sensor with increasing distance. Our
implementation parameterizes 𝑠(𝑘) by defining an area on the screen, and then
calculates the projected size of the area on the screen using the projection theorem.
The transitivity of the characteristic function of the voxel, i.e 𝜒𝑉 (𝑘, 𝑘 0) =
𝜒𝑉 (𝑘, 𝑘 00) ⇔ 𝜒𝑉 (𝑘, 𝑘 0) = 𝜒𝑉 (𝑘 0, 𝑘 00), implies that the set of path vertices can
be partitioned into disjoint sets: 𝜒𝑉 (𝑘, 𝑘 0) equals one for any two vertices with keys
𝑘, 𝑘 0 of the same set, and is always zero otherwise. Figure 3 shows examples for
sets of two-dimensional keys defined by this characteristic function. In practice, keys
are at least three-dimensional, i.e. defined by the world space position of the vertex.
Section 3.1.2 explains the selection of components of the key in detail.
Replacing 𝜒𝐵 in eq. (6) with 𝜒𝑉 , setting 𝑤(𝑥, 𝑥 0) := 1, and using the result as an
approximation instead of only considering its limit, we get
∫
∫ 𝜒 (𝑘, 𝑘 0)𝐿 𝑖 (𝑥 0, 𝜔) 𝑓𝑟 (𝜔𝑟 , 𝑥, 𝜔) cos 𝜃 𝑥0 𝑑𝑥 0
𝜕𝑉 𝑉
𝐿 𝑟 (𝑥, 𝜔𝑟 ) ≈ ∫ 𝑑𝜔. (8)
S−2 ( 𝑥) 𝜒 (𝑘, 𝑘 0 )𝑑𝑥 0
𝜕𝑉 𝑉

If 𝑓𝑟 (𝜔𝑟 , 𝑥, 𝜔) is separable into 𝑓𝑟 (𝜔𝑟 , 𝑥) · 𝑓𝑖 (𝑥, 𝜔), and 𝑓𝑖 (𝑥, 𝜔) is – at least approx-
imately – constant within the voxel, we will be able to rewrite and rearrange eq. (8)
to become
∫
∫ 𝜒 (𝑘, 𝑘 0)𝐿 𝑖 (𝑥 0, 𝜔) 𝑓𝑖 (𝑥 0, 𝜔) cos 𝜃 𝑥0 𝑑𝑥 0
𝜕𝑉 𝑉
𝐿 𝑟 (𝑥, 𝜔𝑟 ) ≈ 𝑓𝑟 (𝜔𝑟 , 𝑥) ∫ 𝑑𝜔. (9)
S−2 ( 𝑥) 𝜕𝑉 𝑉
𝜒 (𝑘, 𝑘 0)𝑑𝑥 0

In the following, we call the product 𝐿 𝑖 (𝑥 0, 𝜔) · 𝑓𝑖 (𝜔𝑟 , 𝑥 0, 𝜔) · cos 𝜃 𝑥0 the contribution

of the vertex in 𝑥 0. All terms of the integrand except for 𝜒𝑉 (𝑘, 𝑘 0) are independent
of 𝑥 and 𝜔𝑟 , and 𝜒𝑉 (𝑘, 𝑘 0) is identical for all vertices with key 𝑘 in a set. Therefore
it is now possible to calculate the integral in eq. (9) only once for the set of vertices
in each voxel.

3.1.2 Construction of Keys

The key 𝑘 of a vertex contains a subset of the components of its descriptor. Thereby,
the integral calculated per voxel is independent of those components excluded from
the key. The selection of components of the key is critically important for defining
the tradeoff between bias and variance reduction: Including additional components
of the descriptor in the key may reduce the bias, while excluding components allows
for the inclusion of more vertices in the integral, therefore reducing variance. In the
following we will give a short overview of components of the descriptor that one
would typically include in the key.
The quality of the approximation in eq. (9) highly depends on the deviation of 𝐿 𝑖
in 𝑥 0 from the one in 𝑥. First and foremost, it is therefore recommended to restrict the
8 Nikolaus Binder, Sascha Fricke, and Alexander Keller

(a) 𝑠 (𝑘) ∝ 1 (b) 𝑠 (𝑘) ∝ 2 blog2 𝑥0 c

Fig. 3 The characteristic function of a voxel 𝜒𝑉 (𝑘, 𝑘 0 ) is 1 for all 𝑘, 𝑘 0 inside the same voxel,
here depicted with the same color for (a) constant resolution selection function 𝑠 (𝑘) resulting in
uniform voxel size or (b) increasing 𝑠 (𝑘) from left to right, which subsequently shrinks the voxels.

world space extent of the voxel by including the position 𝑥 of the vertex in the key 𝑘.
While 𝐿 𝑖 is not continuous in practice – for example in edges of sharp shadows – the
perceived error of approximation decreases with the world space extent of a voxel.
In practice, one can furthermore not guarantee that lim 𝑥0 →𝑥 cos 𝜃 𝑥0 = cos 𝜃 𝑥 due
to different surface orientations in the two locations, for example in sharp edges
of objects. Including the normal of the surface in the point 𝑥 in the key avoids the
resulting “smearing” across edges and “flattening” of surfaces. On the other hand,
the quantization of the normal by eq. (7) causes discontinuities in smooth normals.
Section 3.1.4 details how one can make the artifacts less perceptually pronounced
and more amenable to additional filtering.
Splitting 𝑓𝑟 (𝜔𝑟 , 𝑥, 𝜔) into 𝑓𝑟 (𝜔𝑟 , 𝑥) · 𝑓𝑖 (𝑥, 𝜔) is not always possible. On highly
reflective surfaces, 𝑓𝑟 is defined as a Dirac delta function, and filtering is pointless.
Therefore, vertices on such surfaces should not be selected in the first place. On
glossy surfaces, however, filtering may reduce variance efficiently, again at the cost
of a certain bias. While 𝑓𝑟 can not be split on these surfaces without unpleasantly
and undesirably changing the visual appearance, splitting the domain of incident
angles and computing separate averages for each interval may be a viable tradeoff.
Therefore, for vertices on such surfaces one can append the incident angle to the key
so that the quantization of eq. (7) splits the domain of the incident angle, too.
Materials are often composed of different layers with different properties. Then,
filtering the layers independently offers the opportunity to use different resolutions
𝑠(𝑘) as well as constructing keys with different components for the different layers.
For example, a material consisting of a glossy layer on top of a diffuse one could only
include the angle 𝜔𝑟 in the key used for filtering the glossy layer since the attenuation
of the diffuse layer is independent of it. Thus, the integral used for the diffuse layer
Massively Parallel Path Space Filtering 9

benefits from including more samples. Appending an identifier of the layer to the
key splits the average into several individual ones, which can be combined later.
Note that depending on the number of selected components for the key, voxels are
not necessarily three-dimensional, and their extent may vary between components.

3.1.3 Adaptive Resolution

The choice of the resolution selection function 𝑠(𝑘) is crucial for finding a good
tradeoff between perceived bias and reduction of variance. In theory, its value should
be large in areas with a lot of high frequency detail in 𝐿 𝑖 . Unfortunately, those areas
are almost always unknown in practice since 𝐿 𝑖 is unknown. Figure 4 shows how
a sharp shadow is blurred due to averaging radiance in a large voxel. One would
therefore like to adaptively chose a finer resolution along its boundary.
Finite spatial differences may be used in heuristics for adaptation. While their
computation either introduces a certain overhead or reduces the number of indepen-
dent samples, cost may be amortized over frames in environments changing only
slowly over time. Note that finite differences only estimate spatial variations of the
averages, and one must therefore carefully both choose and adjust such heuristics as
well as determine the number samples used for finite spatial differences.

Fig. 4 Increasing the filter size by lowering the resolution 𝑠 (𝑘) (left to right) increasingly blurs
shadows and also increases the amount of light and shadow leaking.

3.1.4 Filter Kernel Approximation by Jittering

The discontinuities of quantization are removed by jittering components of the key,

which in fact amounts to approximating a filter kernel by sampling. Jittering depends
on the kind of component of the key, for example, positions are jittered in the tangent
plane of an intersection, see algorithm 1. The resulting noise is clearly preferable
over the visible discretization artifacts resulting from quantization, as illustrated and
shown in figs. 5 and 6. In contrast to discretization artifacts, noise from jitter is
simple to filter.
10 Nikolaus Binder, Sascha Fricke, and Alexander Keller

Fig. 5 Jittering trades quantization artifacts for noise. Left: Note that the resolution 𝑠 (𝑘) at the
jittered location (red) may differ from the one of the original location (green). Spatial jittering hides
otherwise visible quantization artifacts (middle): The resulting noise (right) is more amenable to
the eye and much simpler to remove by a secondary filter.

(a) (b) (c) (d)

Fig. 6 Two-dimensional example for averages in voxels: The noisy input (a) is filtered in each
vertex by path space filtering (b). Instead, the new method filters in each voxel, resulting in block
artifacts (c). Additionally jittering before accumulation and lookup resolves the artifacts in noise.

3.2 Accumulation and Lookup in a Hash Table

The averages in each voxel can be calculated in two different ways: First, each voxel
can gather radiance of all included vertices. This process may run in parallel over all
voxels and does not require any synchronization. On the other hand, a list of voxels as
well as a list of vertices per voxel must be maintained. The second way to calculate the
average radiance in a voxel runs in parallel over all vertices: Each vertex atomically
adds its contribution 𝐿 𝑖 (𝑥 0, 𝜔) 𝑓𝑖 (𝑥 0, 𝜔) cos 𝜃 𝑥0 to a running sum of the voxel and
increments the counter of the voxel. The average is finally calculated by dividing the
sum by the counter. While the latter approach requires atomic operations, it does not
involve maintenance of any lists. Furthermore, summation can be parallelized over
the paths or over the vertices, matching the parallelization scheme of typical light
transport simulations. Finally, parallelization per path or per vertex exposes more
parallelism, and therefore the second approach outperforms the first one on modern
graphics processing units (GPUs) significantly.
Accumulation with the latter approach needs a mapping from the key of a vertex
to the voxel with its running sum and counter. Typically, the set of voxels is sparse
since vertices are mostly on two-dimensional surfaces in three-dimensional space.
Additional components of the key increase sparsity even further.
Hash tables provide such a mapping in constant time for typical sets of keys: First,
a hash of the key is calculated using a fast hash function. A modulo operator then
wraps this hash into the index range of the table cells. Since both the hash function
as well as the modulo operator are not bijective, different keys may be mapped to the
same index. Therefore, an additional check for equality of keys is required, and keys
Massively Parallel Path Space Filtering 11

must also be stored in the table. Section 3.2.1 details a cheaper alternative for long
keys.
Upon index collision with a different key, linear probing subsequently increments
the index, checking if the table cell at the updated index is empty or occupied with an
entry with same key. Then, the collision has been resolved. There exist various other
collision resolution methods that improve upon several aspects of linear probing and
have proven to be more efficient in certain use cases, especially for hash tables with
high occupancy. On the other hand, we do not primarily aim to minimize the size of
the hash table, and our experiments show that linear probing comes with a negligible
overhead if the table is sufficiently large. We restrict the number of steps taken for
linear probing to avoid performance penalties of extreme outliers, and resort to the
unfiltered contribution of the vertex if the number of steps exceeds this limit. So far,
our choices for the table size and number of steps so extremely rarely resulted in
such failures that further improvement has been deemed unnecessary. Section 3.2.2
broadens the application of linear probing from only collision resolution to an
additional search for similar voxels.

3.2.1 Fingerprinting

Instead of storing and comparing the rather long keys for checking equality of keys,
we calculate a shorter fingerprint [Rab81, SC08] from a second, different hash
function of the same key and use it for this purpose, see algorithm 1. Excluding a
sentinel value from the fingerprint, we can furthermore use this sentinel to mark
empty cells.
Using fingerprints instead of the full keys is a tradeoff between correctness and
performance: In theory, fingerprints of different keys may coincide. In practice, our
choice of 32 bit fingerprints never caused any collision in our evaluation of several
test scenes and numerous simulations. Still, there is a certain probability of failure,
and we deliberately favor the tiny probability of a failure over the performance
penalty of storing and comparing full keys.

Algorithm 1: Computation of the two hashes used for lookup. Note that the
arguments of a hash function, which form the key, may be extended to refine
clustering.
Input: Location 𝑥 of the vertex, the normal 𝑛, the position of the camera 𝑝cam , and the
scale 𝑠.
Output: Hash 𝑖 to determine the position in the hash table and hash 𝑓 for fingerprinting.
𝑙 ← level_of_detail( | 𝑝cam − 𝑥 |)
𝑥 0 ← 𝑥+ jitter(𝑛) · 𝑠 · 2𝑙
𝑙0 ← level_of_detail( | 𝑝cam − 𝑥 0 |)
j 0 k
𝑥
𝑥˜ ← 𝑠·2 𝑙0

𝑖 ← hash( 𝑥,
˜ . . .)
𝑓 ← hash2( 𝑥,
˜ 𝑛, . . .)
12 Nikolaus Binder, Sascha Fricke, and Alexander Keller

3.2.2 Searching by Linear Probing

In addition, linear probing may be used to differentiate attributes of the light transport
path at a finer resolution as shown in Figure 7: For example, normal information may
be included in the key handed to the fingerprinting hash function instead of already
including it in the main key. This allows one to search for similar normals by linear
probing.
Note that due to other voxels also possibly occupying neighboring cells in the
hash table, searching with linear probing must go beyond mismatching fingerprints.
Therefore, the method works best if both the number of additional contributions as
well as the occupancy of the hash table are low.

a)

b)

c)

Fig. 7 Instead of including normals in the key a) to differentiate contributions whose vertices fall
into the same voxel b), the fingerprint of a key may include normal information. This allows one to
differentiate normal information by linear probing as shown in c).

3.3 Handling Voxels with a low Number of Vertices

Often, there exists a tiny number of voxels that only contain very few vertices,
and therefore variance reduction in those voxels is suboptimal. Examples of such
voxels include those that are only slightly overlapped by objects. Section 3.3.1 and
section 3.3.2 present two approaches that avoid high variance in voxels with only
very few vertices at an almost negligible overhead.
Massively Parallel Path Space Filtering 13

3.3.1 Neighborhood Search

Accumulation in voxels by using quantized keys and a hash table requires a single
atomicAdd operation per component of the radiance of each vertex, and the final
average is computed with one additional non-atomic read operations per component
for the sum and one for the counter. So, as long as access to the hash table happens
in constant time, calculation of the average also takes constant time.
This is in sharp contrast to existing methods computing sums or averages in
a spatial neighborhood which for each vertex take linear time for a search within
a given radius or (typically) logarithmic time for a fixed number of neighboring
vertices.
Even if primarily only one average per voxel is computed and looked up, searching
for neighboring voxels can still be valuable: In theory, increasing the resolution 𝑠(𝑘)
and additionally searching for neighboring voxels may result in variance reduction
similar to the one at a lower resolution, however then with a lower bias. Yet, the
number of neighbors grows exponentially with the number of components of the
key. Thus, such an approach is currently ruled out by the time constraints for real-
time applications. Figure 8 shows a comparison between neighborhood search and
clustering by selecting a coarser resolution.
As a fallback, neighborhood search is still very valuable: If the number of vertices
in a voxel falls below a certain threshold, we allow for an additional search. We
observe that given an appropriate threshold, the number of such voxels is so low that
the overhead is completely negligible while the perceptual improvement is clearly
visible.

Fig. 8 Instead of searching the ns neighborhood in 𝑠 dimensions (left), we utilize clustering

resulting from quantization at a lower resolution to accumulate contributions, which allows for a
single look up (right). We only resort to an additional neighborhood search in the rare case that the
number of vertices in a voxel falls below a certain threshold.

3.3.2 Multi-level Accumulation

Special treatment of voxels with averages from only very few vertices is important
for visual fidelity: Even if the number of voxels with a high variance is very low, it
may be very visible, especially since their appearance is so different from the rest.
Besides searching in a local neighborhood to reduce variance in these cases (see
14 Nikolaus Binder, Sascha Fricke, and Alexander Keller

section 3.3.1), selecting a coarser resolution also effectively increases the number
of vertices in the local average – at the price of an increased bias. Using more than
one resolution at a time avoids the chicken-and-egg problem that arises from first
selecting an appropriate resolution, and then, after accumulation according to this
resolution, determining that it has been set too high or too low.
For simulations in interactive scenarios, one may also select the resolution based
on information from previous frames, see section 3.4.1.

3.4 Accumulation over Time

Reusing contributions across frames dramatically increases efficiency. However,

attention should be paid to arising pitfalls and efficiency aspects: Section 3.4.1 details
the differences and similarities between filtering and integration across frames,
section 3.4.2 explains handling of resolution changes across frames, and section 3.4.3
outlines how the amount of information kept over frames can be limited to avoid
running out of available memory in the hash table.

3.4.1 Temporal Filtering and Temporal Integration

For static scenes, the averages will converge. For dynamic environments, maintaining
two sets of averaged contributions and combining them with an exponential moving
average 𝑐 = 𝛼 · 𝑐 old + (1 − 𝛼) · 𝑐 new is a common tradeoff between convergence and
temporal adaptivity.
However, combining the averages 𝑐 old and 𝑐 new by an exponential moving average
is not equivalent to temporal integration. Especially averages in voxels with relatively
few samples do not converge. In fact, setting 𝛼 := 𝑁old𝑁+𝑁
old
new
correctly integrates across
frames. On the other hand, temporal integration is only possible if the underlying
setting, including lighting conditions and object positions, remains unchanged across
frames.
A first, simple heuristic is to accumulate samples over time up to a certain degree.
This may be implemented using a fixed threshold for the number of samples and
accumulating samples across frames until reaching it. Note that the heuristic is
completely unaware of changes in the scene.
A second, more expensive heuristic builds upon temporal finite differences: A
number of paths is re-evaluated with the same parameters, and the difference of
their contribution to the original ones allows to detect changes that affect the current
voxel. Similar to the spatial finite differences in section 3.1.3, the additional cost
may be amortized over frames, and the number of samples used for finite differences
as well as their influence on the balance between temporal adaptation and temporal
integration must be carefully optimized. Note that averaging in voxels can be used
for the samples used for finite differences, too. Figure 9 shows a comparison between
temporal filtering, temporal integration and a hybrid blending both based on temporal
Massively Parallel Path Space Filtering 15

finite differences. A similar approach for screen space filtering has been explored in
detail by Schied et al. [SPD18].

(a) temporal filtering (b) temporal integration (c) hybrid

Fig. 9 In a scene with a moving light source, temporal filtering using an exponential moving
average blurs shadow boundaries (a), and temporal integration averages out the entire shadow (b).
Adaptively blending 𝛼 between zero (for large temporal differences) and 𝛼 := 𝑁old𝑁+𝑁 old
new
(for no
temporal differences) combines both and preserves sharp shadow boundaries (c). Note that jittering
has deliberately been disabled for all to emphasize sharp boundaries further.

3.4.2 Changes in Resolution across Frames

In many simulations, the camera is dynamic, and therefore the resolution of a voxel
may change across frames if it depends on the position of the camera. Then, already
collected contributions in a voxel at one resolution must be copied to a voxel at either
a higher or lower resolution.
If the resolution in the new frame decreases, one can simply add up the contri-
butions in voxels of higher resolution. Since we store sums and counters, both only
need to be added individually.
If the resolution in the new frame increases, the contributions in voxels of lower
resolution must be distributed to voxels at a higher resolution. Due to the lack of
resolution, this case is much nuanced: On the one hand, using already collected
contributions lowers variance, but on the other hand, the coarser resolution may
become unpleasantly visible. One therefore needs to find a good compromise between
the two, and set 𝛼 accordingly.
Finite spatial or temporal differences can also be filtered across frames, requiring
attention in similar aspects.

3.4.3 Eviction Strategy

Evicting contributions of voxels which have not been queried for a certain period
of time is necessary for larger scenes and changing camera. Besides the least re-
cently used (LRU) eviction strategy, heuristics based on longer term observations
are efficient.
16 Nikolaus Binder, Sascha Fricke, and Alexander Keller

A very simple implementation relies on replacing the most significant bits of the
fingerprinting hash by a priority composed of for example the number of vertices in
the voxel and last access time during temporal filtering. Thus the pseudo-randomly
hashed least significant bits guarantee eviction to be uniformly distributed across the
scene, while the most significant bits ensure that contributions are evicted according
to priority. This allows collision handling and eviction to be realized by a single
atomicMin operation.

4 Results and Discussion

While filtering contributions at primary intersections with the proposed algorithm is

quite fast, it only removes some artifacts of filtering in screen space. However, hashed
path space filtering has been designed to target real-time light transport simulation:
It is the only efficient fallback when screen space filtering fails or is not available,
for example, for specular or transparent objects.
Filtering on non-diffuse surfaces requires to include additional parameters in the
key and heuristics such as increasing the quantization in areas with non-diffuse
materials to minimize the visible artifacts.
Filtering, and especially accumulating contributions, is always prone to light
and shadow leaking (see fig. 4), which is the price we pay for performance. Some
artifacts may be ameliorated by employing suitable heuristics as reviewed in [KDB16,
Sec.2.1] and in section 3.1.3.
The new algorithm filters incoherent intersections at HD resolution (1920 × 1080
pixels) in about 3 ms on an NVIDIA Titan V GPU. Filtering primary intersections
doubles the performance due to the more coherent memory access patterns.
The image quality is determined by the filter size, which balances noise versus blur
as shown in Figure 4. Both the number of collisions in the hash table and hence the
performance of filtering also depend on the size of the voxels. We found specifying
the voxel size by 𝑠-times the projected size of a pixel most convenient. Note that
maximum performance does not necessarily coincide with best image quality. The
hash table size is chosen proportional to the number of pixels at target resolution such
that potentially one vertex could be stored per pixel. In practice, filtering requires
multiple vertices to coincide in a voxel, and therefore the occupancy of the hash
table is rather low. Such a small occupancy improves the performance as it lowers
the number of collisions and time spent for collision resolution.
While path space filtering dramatically reduces the noise at low sampling rates
(see fig. 1), some noise is added back by spatial jittering. Instead of selecting the
first sufficiently diffuse vertex along a path from the camera, path space filtering
can be applied at any vertex. For example, filtering at the second sufficiently diffuse
vertex as shown in fig. 10 resembles final gathering or local passes [Jen01]. Further,
it is obviously also possible to filter in several vertices along the path at the same
time. In fact, path space filtering trades variance reduction for controlled bias and
is orthogonal to other filtering techniques. We therefore abstain from comparisons
Massively Parallel Path Space Filtering 17

with these: Temporal anti-aliasing and complimentary noise filters in screen space are
appropriate to further reduce noise [MMBJ17]. A local smoothing filter [SKW+ 17]
can even help reduce the error inherent in the approximation.

Fig. 10 Indirect illumination by hashed path space filtering only at the second bounce: At 16 paths
per pixel (left), the variance of the integrand is dramatically reduced (right).

5 Conclusion

Relying on only a few synchronizations during accumulation, path space filtering

based on hashing scales on massively parallel hardware. Both accumulation as well as
queries run in constant time per vertex, and neither the traversal nor the construction
of a hierarchical spatial acceleration data structure is required. Ray tracing hardware,
accelerating ray intersection with the scene, is orthogonal to the presented method.
In fact, reasonable visual fidelity can often only be achieved if ray tracing hardware
is combined with such a powerful variance reduction technique due to the time
constraints for real-time simulations.
The simplistic algorithm overcomes many restrictions of screen space filtering,
does not require motion vectors, and enables noise removal beyond the first intersec-
tion including specular and transparent surfaces.
The hashing scheme still bears potential for improvement. For example, important
hashes could be excluded from eviction by reducing the resolution, i.e. accumulating
their contributions at a coarser level. Other than selecting the resolution by the length
of the path, path differentials and variance may be used to determine the appropriate
resolution.
18 Nikolaus Binder, Sascha Fricke, and Alexander Keller

Besides the classic applications of path space filtering [KDB16, Sec.3] like
multi-view rendering, spectral rendering, participating media, and decoupling anti-
aliasing from shading, the adaptive hashing scheme can be used for photon mapping
[Jen01, HJ09] and irradiance probes in reinforcement learned importance sampling
[DK18] in combination with final gathering to store radiance probes. Since the first
publication of this work as a technical report variants of the presented method has
been applied to improving the efficiency of ambient occlusion [Gau20] and for an
efficient implementation of reinforcement learning [Pan20].

Acknowledgements

The authors would like to thank Petrik Clarberg for profound discussions and com-
ments.

References

BSH02. P. Bekaert, M. Sbert, and J. Halton. Accelerating Path Tracing by Re-Using Paths.
In P. Debevec and S. Gibson, editors, Eurographics Workshop on Rendering. The
Eurographics Association, 2002.
BVM+ 17. S. Bako, T. Vogels, B. McWilliams, M. Meyer, J. Novák, A. Harvill, P. Sen, T. Derose,
and F. Rousselle. Kernel-predicting convolutional networks for denoising Monte Carlo
renderings. ACM Trans. Graph., 36(4):97:1–97:14, July 2017.
CKS+ 17. C. Chaitanya, A. Kaplanyan, C. Schied, M. Salvi, A. Lefohn, D. Nowrouzezahrai, and
T. Aila. Interactive reconstruction of Monte Carlo image sequences using a recurrent
denoising autoencoder. ACM Trans. Graph., 36(4):98:1–98:12, July 2017.
Dam11. H. Dammertz. Acceleration methods for ray tracing based global illumination. PhD
thesis, Universität Ulm, 2011.
DK18. K. Dahm and A. Keller. Learning light transport the reinforced way. In A. Owen
and P. Glynn, editors, Monte Carlo and Quasi-Monte Carlo Methods. MCQMC 2016.
Proceedings in Mathematics & Statistics, volume 241, pages 181–195. Springer, 2018.
DS07. A. Dietrich and P. Slusallek. Adaptive spatial sample caching. In 2007 IEEE Symposium
on Interactive Ray Tracing, pages 141–147, 2007.
Gau20. P. Gautron. Real-time ray-traced ambient occlusion of complex scenes using spatial
hashing. In Special Interest Group on Computer Graphics and Interactive Techniques
Conference Talks, SIGGRAPH ’20, New York, NY, USA, 2020. Association for Com-
puting Machinery.
GDW+ 14. P. Gautron, M. Droske, C. Wächter, L. Kettner, A. Keller, N. Binder, and K. Dahm. Path
space similarity determined by Fourier histogram descriptors. In ACM SIGGRAPH
2014 Talks, SIGGRAPH ’14, pages 39:1–39:1, New York, NY, USA, 2014. ACM.
GKDS12. I. Georgiev, J. Křivánek, T. Davidovič, and P. Slusallek. Light transport simulation with
vertex connection and merging. ACM Trans. Graph., 31(6):192:1–192:10, November
2012.
HHS05. V. Havran, R. Herzog, and H.-P. Seidel. Fast final gathering via reverse photon mapping.
Computer Graphics Forum, 24(3):323–332, 2005.
HJ09. T. Hachisuka and H. Jensen. Stochastic progressive photon mapping. In SIGGRAPH
Asia ’09: ACM SIGGRAPH Asia 2009 papers, pages 1–8. ACM, 2009.
Massively Parallel Path Space Filtering 19

HJ10. Toshiya Hachisuka and Henrik Jensen. Parallel progressive photon mapping on GPUs.
SIGGRAPH Sketches, 01 2010.
Jen01. H. Jensen. Realistic Image Synthesis Using Photon Mapping. AK Peters, 2001.
KDB16. A. Keller, K. Dahm, and N. Binder. Path space filtering. In R. Cools and D. Nuyens,
editors, Monte Carlo and Quasi-Monte Carlo Methods 2014, pages 423–436. Springer,
2016.
Kel98. A. Keller. Quasi-Monte Carlo Methods for Photorealistic Image Synthesis. PhD thesis,
University of Kaiserslautern, Germany, 1998.
KRK04. J. Kontkanen, J. Räsänen, and A. Keller. Irradiance filtering for Monte Carlo ray tracing.
In D. Talay and H. Niederreiter, editors, Monte Carlo and Quasi-Monte Carlo Methods
2004, pages 259–272. Springer, 2004.
LW98. E. Lafortune and Y. Willems. Bi-directional path tracing. Proceedings of Third Inter-
national Conference on Computational Graphics and Visualization Techniques (Com-
pugraphics’, 93, 01 1998.
McC02. M. McCool. Anisotropic diffusion for Monte Carlo noise reduction. ACM Transactions
on Graphics, 18, 08 2002.
MGN17. T. Müller, M. Gross, and J. Novák. Practical path guiding for efficient light-transport
simulation. Computer Graphics Forum, 36(4):91–100, June 2017.
MHC+ 16. J. Munkberg, J. Hasselgren, P. Clarberg, M. Andersson, and T. Akenine-Möller. Texture
space caching and reconstruction for ray tracing. ACM Trans. Graph., 35(6):249:1–
249:13, November 2016.
MLM13. M. Mara, D. Luebke, and M. McGuire. Toward practical real-time photon mapping:
Efficient GPU density estimation. In Proceedings of the ACM SIGGRAPH Symposium
on Interactive 3D Graphics and Games (I3D’13), March 2013. Interactive 3D Graphics
and Games 2013.
MM02. V. Ma and M. McCool. Low latency photon mapping using block hashing. In T. Ertl,
W. Heidrich, and M. Doggett, editors, SIGGRAPH/Eurographics Workshop on Graphics
Hardware. The Eurographics Association, 2002.
MMBJ17. M. Mara, M. McGuire, B. Bitterli, and W. Jarosz. An efficient denoising algorithm for
global illumination. In ACM SIGGRAPH / Eurographics High Performance Graphics,
page 7, July 2017.
MMR+ 19. R. Müller, B. McWilliams, F. Rousselle, M. Gross, and J. Novák. Neural importance
sampling. ACM Trans. Graph., 38(5):145:1–145:19, October 2019.
Pan20. J. Pantaleoni. Online path sampling control with progressive spatio-temporal filtering,
2020.
Rab81. M.O. Rabin. Fingerprinting by random polynomials. Technical report, Center for
Research in Computing Technology, Harvard University, 1981.
SC08. M. Slaney and M. Casey. Locality-sensitive hashing for finding nearest neighbors
[lecture notes]. IEEE Signal Processing Magazine, 25(2):128–131, 2008.
SKW+ 17. C. Schied, A. Kaplanyan, C. Wyman, A. Patney, C. Chaitanya, J. Burgess, S. Liu,
C. Dachsbacher, A. Lefohn, and M. Salvi. Spatiotemporal variance-guided filtering:
Real-time reconstruction for path-traced global illumination. In Proceedings of High
Performance Graphics, HPG ’17, pages 2:1–2:12, New York, NY, USA, 2017. ACM.
SPD18. C. Schied, C. Peters, and C. Dachsbacher. Gradient estimation for real-time adaptive
temporal filtering. Proc. ACM Comput. Graph. Interact. Tech., 1(2), Aug 2018.
SZR+ 15. P. Sen, M. Zwicker, F. Rousselle, S.-E. Yoon, and N. Kalantari. Denoising your Monte
Carlo renders: Recent advances in image-space adaptive sampling and reconstruction.
In ACM SIGGRAPH 2015 Courses, SIGGRAPH ’15, pages 11:1–11:255, New York,
NY, USA, 2015. ACM.
Vea97. E. Veach. Robust Monte Carlo Methods for Light Transport Simulation. PhD thesis,
Stanford University, 1997.
WGGH20. R. West, I. Georgiev, A. Gruson, and T. Hachisuka. Continuous multiple importance
sampling. ACM Transactions on Graphics (Proceedings of SIGGRAPH), 39(4), Jul
2020.
20 Nikolaus Binder, Sascha Fricke, and Alexander Keller

ZJL+ 15. M. Zwicker, W. Jarosz, J. Lehtinen, B. Moon, R. Ramamoorthi, F. Rousselle, P. Sen,

C. Soler, and S.-E. Yoon. Recent advances in adaptive sampling and reconstruction for
Monte Carlo rendering. Comput. Graph. Forum, 34(2):667–681, May 2015.