Projective Sampling for Dierentiable Rendering of Geometry
ZIYI ZHANG, École Polytechnique Fédérale de Lausanne (EPFL), Switzerland
NICOLAS ROUSSEL, École Polytechnique Fédérale de Lausanne (EPFL), Switzerland
WENZEL JAKOB, École Polytechnique Fédérale de Lausanne (EPFL), Switzerland
(a) Primal rendering and perturbation of lamp rotation
(b) Forward derivative computed using projective sampling
-2.0
2.0
Fig. 1. The visibility function plays a crucial role in dierentiable rendering of geometry. Parameter changes that influence visibility (e.g., a rotation of the
light source in (a) can generate a significant derivative contribution shown in (b) that is diicult to sample using currently existing methods. We project path
segments generated by primal rendering algorithms (e.g., direct illumination sampling) onto nearby silhouees and organize them in a uniform or adaptive
guiding data structure to enhance numerical integration of the challenging boundary term. Compared to prior work [Zhang et al
.
2020], uniform guiding cuts
errors (RMSE) by an average of 8.1×, and adaptive guiding yields an additional 2.7× improvement.
Discontinuous visibility changes at object boundaries remain a persistent
source of diculty in the area of dierentiable rendering. Left untreated, they
bias computed gradients so severely that even basic
optimization tasks fail.
Prior path-space methods addressed this bias by decoupling boundaries
from the interior, allowing each part to be handled using specialized Monte
Carlo sampling strategies. While conceptually powerful, the full potential of
this idea remains unrealized since existing methods often fail to adequately
sample the boundary proportional to its contribution.
This paper presents theoretical and algorithmic contributions. On the
theoretical side, we transform the boundary derivative into a remarkably
simple local integral that invites present and future developments.
Building on this result, we propose a new strategy that projects ordinary
samples produced during forward rendering onto nearby boundaries. The
resulting projections establish a variance-reducing guiding distribution that
accelerates convergence of the subsequent dierential phase.
We demonstrate the superior eciency and versatility of our method
across a variety of shape representations, including triangle meshes, implic-
itly dened surfaces, and cylindrical bers based on Bézier curves.
CCS Concepts: • Computing methodologies → Rendering.
Additional Key Words and Phrases: dierentiable rendering, geometry re-
construction
Authors’ Contact Information: Ziyi Zhang, École Polytechnique Fédérale de Lausanne
(EPFL), Lausanne, Switzerland, ziyi.zhang@ep.ch; Nicolas Roussel, École Polytech-
nique Fédérale de Lausanne (EPFL), Lausanne, Switzerland, nicolas.roussel@ep.ch;
Wenzel Jakob, École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland,
wenzel.jakob@ep.ch.
© 2023 Copyright held by the owner/author(s). Publication rights licensed to ACM.
This is the author’s version of the work. It is posted here for your personal use. Not for
redistribution. The denitive Version of Record was published in ACM Transactions on
Graphics, https://doi.org/10.1145/3618385.
ACM Reference Format:
Ziyi Zhang, Nicolas Roussel, and Wenzel Jakob. 2023. Projective Sampling for
Dierentiable Rendering of Geometry. ACM Trans. Graph. 42, 6, Article 212
(December 2023), 14 pages. https://doi.org/10.1145/3618385
1 Introduction
Dierentiation is a core component of recent work on inverse ren-
dering. Methods in this area combine the derivative of a rendering
algorithm with a variant of gradient descent to optimize a tentative
scene until renderings become consistent with observed data.
However, a frustrating problem that invariably arises when dif-
ferentiating a Monte Carlo renderer is that the computed gradient
is wrong, to the extent that even basic optimizations fail. More
precisely, derivatives with respect to certain kinds of parameters
including vertex positions, curve control points, and transformation
matrices are contaminated by bias. The relevant shared characteris-
tic of these examples is their eect on silhouettes, which refers to
the outlines of objects that occlude the background, other objects,
or even a part of themselves. At silhouettes, small perturbations
to a ray’s origin or direction cause a sudden jump in the closest
ray-surface intersection.
The combination of this property with the ubiquitous Monte
Carlo method leads to a curious mathematical artifact: changing
“problematic”
parameters causes sudden jumps in sampled ray inter-
sections, which should have a noticeable inuence on the rendered
image and its derivative. However, the eect is lost under dierenti-
ation, since the change at individual Monte Carlo samples is discon-
tinuous. This defect is not present in the original continuous theory,
which motivates techniques that seek to remove the discrepancy.
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
212:2 • Zhang et al.
Prior work in this area falls into two broad categories: boundary-
based methods generate additional Monte Carlo samples on silhou-
ettes to account for their eect, while area-based methods compute a
modied integral using ordinary samples covering the full domain.
Our method is a hybrid in this classication: it ingests ordinary
samples and converts them into boundary samples. We feed it with
the output of standard (“primal”) rendering strategies, which follows
Griewank and Walther [2008]’s rule #2 for automatic dierentiation:
“What is go od for function values is good for their derivatives.”. In par-
ticular, BRDF-, emitter-, and multiple importance sampling [Veach
and Guibas 1995] are indispensable tools for ecient primal render-
ing, and we similarly wish to use them to improve the quality of
gradient estimates.
Our method builds on the path-space dierentiable rendering
by Zhang et al
.
[2020] and makes both theoretical and practical
contributions to this framework. We revisit the local parameteriza-
tion of the boundary integral and bring the primary equation into
its ultimate reduced form. We also examine visibility discontinu-
ities arising from the interior of smooth geometry and propose a
generalized local formulation that subsumes all cases.
Contemporary rendering systems support a great variety of shape
representations. Our technique requires users to supply a projection
operator per type that maps ray segments onto a nearby silhouette.
We demonstrate how this modular approach enables dierentiable
rendering of triangle meshes, implicit functions, and cylindrical
bers based on Bézier curves, while producing gradients with an
exceptionally low level of variance compared to prior work. Besides
improving computational eciency, low gradient variance can also
improve the stability of optimizations (Figure 2).
Following a review of prior work, Section 3 will rst introduce
the core idea in a simple 2D setting followed by a mathematical for-
malization of the boundary integral in Section 4. Sections 5 covers
algorithmic and implementation-level details, and Section 6 provides
experimental comparisons to prior area and boundary-based meth-
ods. An open-source implementation of this method is available at
https://github.com/mitsuba-renderer/mitsuba3.
Optimization task With high-quality gradients With low-quality gradients
Converged
Failed
Fig. 2. Gradient quality and robustness. Stochastic optimizers like Adam
support noisy gradients, but excessive variance can be detrimental. In this
example, we infer the X/Y oset of the knot shape from a single reference
view. The middle and right images depict the convergence status following
initialization at the associated point within the parameter domain (white
indicates a failure to converge to the known parameter value). The only
dierence between these experiments is the quality of the gradient, which
was computed using either finite dierences with 128 samples/pixel (middle),
or a reparameterization with 32 auxiliary rays [Bangaru et al
.
2020] (right).
Physically based rendering LanguagesRendering
Edge sampling
Path-space sampling
Reparameterization /
Divergence theorem
TEG (Analytic)
Boundary sampling
for Vector Graphics
Smooth rasterization
SDF Reparam.
Aδ (Finite dierences)
Boundary-based Area-based
Fig. 3. A taxonomy of prior work. Visible discontinuities can be handled
by smoothing boundaries [Liu et al
.
2019] or explicitly sampling them [Li
et al
.
2020]. The physically-based seing features indirectly observed dis-
continuities. Previous boundary-based sampling methods used 6D Hough
trees [Li et al
.
2018] or path space [Zhang et al
.
2020]. Alternatively, an
area-based integral can be reparameterized to freeze discontinuities [Lou-
bet et al
.
2019], which is equivalent to an application of the divergence
theorem [Bangaru et al
.
2020]. SDFs are particularly amenable to reparam-
eterization [Vicini et al
.
2022; Bangaru et al
.
2022]. Recent dierentiable
programming languages dierentiate discontinuous integrals using analytic
methods [Bangaru et al. 2021] and finite dierences [Yang et al. 2022].
This article includes numerous visualizations of forward deriva-
tives (for instance, see Figure 1). These characterize the change of
rendered pixels in response to a perturbation of a single scene pa-
rameter. In reality, the proposed method would more likely be used
to compute backward derivatives of all parameters relative to a given
image loss. We opt to show forward derivatives since they more
intuitively convey information about sources of bias and variance.
2 Related Work
This section reviews relevant prior work for handling discontinuities
in derivatives of rendering algorithms. Figure 3 organizes a subset
of them in a visual taxonomy.
2.1 Rendering
Dierentiable rasterization. In its simplest form, rasterization
converts polygonal meshes into fragments that are shaded and
then depth-tested to form a raster image. Dierentiable rasteriza-
tion [Loper and Black 2014; Kato et al
.
2018; Liu et al
.
2019; Cole
et al
.
2021; Petersen et al
.
2022] replaces specic steps of this pro-
cess with smooth analogs, e.g. by blurring polygon boundaries or
smoothly blending fragments. Dierentiable rasterization can be
made impressively ecient [Laine et al
.
2020] but does not account
for indirect eects like shadows or indirect illumination.
Dierentiable vector graphics. Li et al. [2020] propose a dieren-
tiable antialiasing technique for vector graphics built from cubic
Bézier curves. They apply the Reynolds transport theorem [1903] to
turn derivative contributions from boundaries into a 1D boundary
integral and sample it using the Monte Carlo method.
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
Projective Sampling for Dierentiable Rendering of Geometry • 212:3
Other methods for primary visibility. A large number of prior
works propose specialized methods for dierentiable rendering of
SDFs, volumes, neural elds, and other representations involving
mask or silhouette losses to account for visible discontinuities. To
keep the discussion focused, we skip them along with work on
inverse rendering pipelines, losses, and regularization techniques.
2.2 Physically based rendering
The physically based setting introduces ve additional challenges:
1.
Methods must account for indirectly observed discontinuities
on a higher-dimensional space of light paths. For example, the
amount of shadowing below an object depends on the degree
to which its geometry occludes surrounding sources of direct
and indirect illumination. Figure 4 visually decomposes the
components of the resulting derivative.
2.
The reectance integral at every scene position observes a dif-
ferent set of silhouettes. This means that detecting and storing
these silhouettes ahead of time is not straightforward.
3.
Compared to the vector graphics setting [Li et al
.
2020], there are
too many discontinuities to enumerate. They must be handled
randomly, ideally in some way that guides the computation
towards important ones.
4.
Not all discontinuities are equal: their contribution is propor-
tional to the discontinuous change in the primal integrand (Fig-
ure 5b), which is often highly heterogeneous.
5.
As the scene complexity grows, an increasing fraction of poten-
tial silhouettes becomes irrelevant because they are themselves
occluded by surrounding geometry.
No existing solution fully addresses this challenging set of con-
straints. The following ideas have been proposed in the past:
Boundary sampling. Li et al. [2018] employ a 6D Hough tree to
sample potentially observed boundaries relative to a given scene
position. This pioneering method was the rst to address indirect
visibility, though it comes with various practical limitations: poor
scaling of the high-dimensional data structure with respect to scene
complexity and the diculty of drawing samples proportional to
their contribution while excluding invisible silhouettes.
Zhang et al. [2020] subsequently approached the problem through
Veach’s [1997] path space formalism which reveals that silhouette
paths are more easily generated “from the middle” by sampling
a direction tangential to any edge and expanding it outwards to
form a complete path connecting the sensor to an emitter. The
core sampling problem entails selecting an edge position
𝑡
and a
tangential direction
𝜃, 𝜙
proportional to the expected value of a
Monte Carlo estimator
⟨𝑓 (𝑡, 𝜃, 𝜙)⟩
. This is not easy, as
𝐸[𝑓 ]
can
uctuate by many orders of magnitude. Zhang et al. [2020] tabulate
𝐸[𝑓 ] on a dense 3D grid to sample relevant parts of this domain.
The expense of tabulating at a suciently ne resolution can
become a critical bottleneck, which has motivated subsequent work
on adaptive representations that can more easily accommodate nar-
row peaks in
𝐸[𝑓 ]
. Yan et al
.
[2022] generate a sequence of separate
kd-trees over chains of adjacent edges and subdivide them recur-
sively until random samples in leaf nodes are suciently uniform.
As with any other adaptive integration technique, the method may
fail to detect sharp peaks and stop the subdivision prematurely.
Continuous ValueDiscontinuous
∂ component
1
st
path segment 2
nd
path segment 3
rd
path segment
Dierentiation task
Primal rendering
Fig. 4. Image and derivative decomposition. A physically based renderer
must account for various illumination sources, which further grows when
considering their derivative. The columns below separate the contribution
of visible emiers (le) and shapes subject to direct (middle) and indirect
illumination (right). The continuous derivative (row 3) is easily handled
using automatic dierentiation or adjoint methods [Nimier-David et al
.
2020; Vicini et al
.
2021]. The contribution of discontinuities (boom row)
poses severe challenges and is the primary focus of this article.
Current methods in this category only support polygonal meshes.
This is unfortunate, since continuous shape representations oer
desirable qualities: implicit representations like signed distance func-
tions (SDFs) support topological changes and exhibit improved
stability in non-regularized optimizations compared to discrete
meshes [Vicini et al
.
2022]. Despite these benets, it is still unclear
how the idea of boundary sampling could be adapted to continuous
representations that lack a notion of discrete edges.
Reparameterization. Instead of sampling the boundary, methods
in this category reparameterize the interior integral using coordi-
nates that smoothly deform the evolving domain in such a way
that discontinuous boundaries become stationary. Boundary-related
derivatives then arise from the interior deformation, which is cap-
tured by the parameterization’s Jacobian determinant. The rst work
in this sequence [Loubet et al
.
2019] proposed a heuristic parameter-
ization, which was later improved using a construction with lower
bias [Bangaru et al. 2020].
Visibility changes at object silhouettes are not the only source of
discontinuities: for example, geometric normals can introduce shad-
ing discontinuities at the edges of polygonal meshes, and this can
similarly be handled by a reparameterization of the interior known
as the material form [Zhang et al
.
2020]. We only focus on visibility-
induced discontinuities but note that these ideas are orthogonal and
could be combined to handle both types of discontinuities.
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
212:4 • Zhang et al.
(a) Geometric setup
Boundary
Light source
(b) Discontinuous integral at
Tangential
(c) Projective sampling
Fig. 5. High-level overview. (a) Rendering of a bunny with translation parameter
𝜃
. Increasing the value of
𝜃
brightens the partially shadowed surface
position x
a
. Image (b) shows the spherical integral that determines the reflectance at x
a
, which shows how increasing
𝜃
shis the silhouees (red curves)
towards the le and reveals more of the partially blocked light source. To account for this eect during dierentiation, one can place additional Monte Carlo
samples directly onto the boundary by generating tangential path segments
(
x
a
,
x
b
,
x
c
)
. Our method leverages standard primal sampling techniques to find
relevant parts of this boundary. The example in (c) shows a sample from a direct illumination strategy (blue) that was ultimately unsuccessful due to occlusion.
Our method takes this segment and projects it onto a nearby silhouee.
-0.02
0.02
-8.0
8.0
Continuous
part only
Rotating sphereWavy rectangle
Ground truth
Reparameterized
Scene
Bias
Fig. 6. Reparameterizations [Bangaru et al
.
2020; Loubet et al
.
2019] inject
variance and bias into the derivative computation. (top.) This uniform rotat-
ing sphere should not produce visible derivatives. The reparameterization
fails to identify this and produces biased output. (bot.) The derivative of a
translating wavy rectangle has interior and boundary components. A naïve
Monte Carlo estimator misses the boundary derivative but produces an
otherwise well-converged estimate of the interior part. The reparameterized
version fixes the boundary at the cost of a global variance increase.
Reparameterization-based methods have clear benets: they work
with any geometric representation and are tightly integrated into
the original simulation, which allows them to only focus on visible
scene geometry. On the ipside, they must trace many auxiliary rays
to sense the motion of surrounding geometry, which is necessary to
create appropriate coordinate maps for a given path segment. The
numerical aspects of this process are somewhat arcane and involve
convolutions with singular kernels that inject a substantial amount
of variance into the dierentiable rendering process. This variance
even impacts the interior parts far from silhouettes. We have found
it dicult to congure these methods to strike an acceptable bal-
ance between variance and bias (see Figure 6). In some instances,
reparameterizations can be tailored to the underlying geometric
representation, and two recent works explore the specialization to
SDFs [Vicini et al. 2022; Bangaru et al. 2022].
Analytic visibility. A third way of addressing the issue involves
substituting the Monte Carlo method with analytic solutions that
inherently yield correct derivative when dierentiated [Zhou et al.
2021]. The requirement of closed-form integrability is relatively
stringent and restricts this approach to directly illuminated surfaces
with simple material models.
Dierentiable programming. Other works have explored the ef-
fect of discontinuous decision boundaries (e.g., an
if
) in the con-
text of shaders [Yang et al
.
2022] and domain-specic languages
with an integration primitive [Bangaru et al. 2021].
We believe that some ideas in this article may apply to this more
general setting, e.g., by repeatedly running a program with per-
turbed inputs and adding a root-nding iteration that projects sam-
ples onto the decision boundary to account for its derivative.
3 Overview
A brief comment on terminology: our method computes a boundary
integral to account for the eect of visibility during dierentiation,
but this integral itself has an interior and a boundary term. This is
because both the interior (a 2D domain) and boundary (a 1D curve)
of a surface (e.g., a curved patch) can impact visibility. To reduce
severe terminological overload of the word boundary, we will refer
to the 1D domain of the boundary integral as the perimeter.
We now turn to an example to review the mathematical nature
of the problem and present the high-level idea of our approach.
Figure 5a illustrates the geometric setup: a camera observes an object
casting a smooth shadow due to a nearby area emitter. A solitary
scene parameter 𝜃 controls the shape’s horizontal translation
1
.
1
This single-parameter example is only for illustration. Normal usage of a dierentiable
renderer involves dierentiation with respect to millions of parameters.
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
Projective Sampling for Dierentiable Rendering of Geometry • 212:5
The following hemispherical cosine-weighted integral determines
the diuse outgoing radiance 𝐿
𝑜
(x
a
) at position x
a
:
𝐿
𝑜
(x
a
) = 𝜌
∫
H
2
𝐿
𝑖
(x
a
, 𝜔) d𝝎
⊥
, (1)
where
𝐿
𝑖
denotes the incident radiance that implicitly depends on
𝜃
,
and
𝜌
refers to the object’s diuse albedo. Increasing
𝜃
moves the
shadow away from x
a
, which in turn increases the value of
𝐿
𝑜
(
x
a
)
.
The derivative
𝜕
𝜃
𝐿
𝑜
expresses the precise rate of change and is
given by an integral of a derivative (a.k.a. “dierentiation under
the integral sign”) and an additional boundary term illustrated in
Figure 5b:
𝜕
𝜃
𝐿
𝑜
(x
a
) =𝜌
∫
H
2
𝜕
𝜃
𝐿
𝑖
(x
a
, 𝝎 ) d𝝎
⊥
+ 𝜌
∫
B
(𝜕
𝜃
x(𝑡) · n(𝑡)) Δ𝐿
𝑖
d𝑡. (2)
The rst term is readily computable using existing methods like
automatic dierentiation or Path Replay Backpropagation [Vicini
et al
.
2021], but it only plays a minor role in this example. The dom-
inant second term integrates the discontinuous change in incident
radiance
Δ𝐿
𝑖
across the evolving boundary
B
scaled by its perpen-
dicular velocity. The computation of this boundary term presents
serious diculties, as detailed in Section 2.
Our work on this topic was prompted by the realization that a
boundary projection operation could greatly reduce these diculties.
Such a projection would instantly turn any existing technique for
the interior into a specialized method for boundaries that could then
be used to compute the troublesome integral.
Figure 5c depicts a possible application of this idea: suppose that
a light source sample turns out to be occluded as seem from position
x
a
. This provides a helpful clue for the dierential phase: if the
occlusion is only partial, we may be able to “walk” towards the
associated boundary to account for its derivative contribution.
While such a projection can surely be constructed, this imme-
diately raises another question: how would one actually use it in
the Monte Carlo context? The estimator would need to consider
the probability
𝑝 (𝑡)
of the generated boundary sample
𝑡 ∈ B
per
unit length (or more accurately, its reciprocal). The density
𝑝 (𝑡)
represents an intricate marginal dependent on both the input 2D
density and the specics of the projection. We cannot expect that
this marginal will be available in closed form.
We experimented with unbiased methods that draw auxiliary
samples to estimate the reciprocal density [Booth 2007] but found
the overheads of this process to be unacceptably high. In essence,
the intricacy of the projection is such that attempts to characterize
it precisely negate the benet of having a
projection to begin with.
What, then, if we tried to model
𝑝 (𝑡)
imprecisely? This idea is
the basis of our method, which consists of three steps: the rst is
an ordinary primal rendering step that additionally projects sam-
ples onto nearby boundaries (i.e., it draws samples from
𝑝 (𝑡)
). The
second builds a closed-form guiding distribution
ℎ(𝑡) ≈ 𝑝(𝑡)
from
the projections. The last step computes the boundary integral by
importance sampling the guiding distribution. If
ℎ(𝑡)
covers rele-
vant parts of the domain, such a method will still produce unbiased
estimates despite the approximate nature of
ℎ
. Section 5 provides
algorithmic and implementation-level details needed to turn this
high-level idea into a practical method.
Fig. 7. Boundary sample space. Visualization of the boundary integrand
of the Bunny and Filigree meshes from Figure 10 on boundary sample space.
This 3D domain is parameterized by edge position
𝑡 ∈ 𝜕A
and a spherical
direction
(𝜃, 𝜙 ) ∈ S
2
. Note the extremely sparse and high-frequency nature
of this function (white corresponds to an integrated value of 0 along a ray
viewing the 3D volume). The supplemental video contains an animated
version of this figure.
3.1 Discussion
This high-level summary already provides enough context to discuss
the method’s conceptual advantages and disadvantages.
1.
Modularity. Our method does not rely on any particular repre-
sentation and ts into the modular architecture of contempo-
rary rendering systems. To support a new geometric primitive,
the user must implement two, rarely three operations: the pre-
viously mentioned projection, and a parameterization of the
geometric perimeter (if present) or curved interior (if present).
2.
Heterogeneity. The value of the boundary integrand tends to
uctuate by many orders of magnitude. Figure 5b illustrates this:
the boundary curve has a considerable arclength, but only the
small segment directly adjacent to the
light source contributes.
Our method can leverage an existing ecosystem of sampling
techniques designed to cover high-valued regions of the pri-
mal integrand with high probability, which in turn benets the
boundary integral. Without added eort on our part, the com-
putation stays conned to visible scene geometry, which has
been a notable challenge in previous work.
3.
Guiding. Despite the advantage mentioned above, we nd that
the density
𝑝
resulting from the projection is often suboptimal
and not suciently proportional to the value of the actual bound-
ary integrand. The decoupled nature of the guiding distribution
ℎ ≠ 𝑝
provides the means to address this problem, since it al-
lows us to adjust
ℎ
retroactively using more accurate knowledge
about the boundary integrand.
Figure 7 visualizes the exceedingly sparse nature of the inte-
grand on the boundary sample space covered by our guiding
scheme. The primary role of the projection is to nd the narrow
nonzero regions. We subsequently evaluate the actual integrand
at the projected locations to create the nal guiding distribution.
4.
Sampling eciency. Our method improves statistical e-
ciency compared to prior path-space methods [Zhang et al
.
2020;
Yan et al
.
2022], and it improves handling of the interior com-
ponent compared to reparameterizations [Loubet et al
.
2019;
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
212:6 • Zhang et al.
Edge
Subdomain 1
Boundary integral
Subdomain 2
(d) Local (interior)(c) Local (perimeter)(a) Non-local (perimeter) (b) Non-local (interior)
Fig. 8. Formulations and terms of the boundary integral. Visibility-related derivatives arise from the perimeter (e.g., discrete edges of a triangle mesh)
and the interior of shapes (e.g. the surface of an ellipsoid). Path-space methods compute an integral over tangential path segments to account for them.
Decomposing the integration domain (blue and orange sets) reveals dierent formulations: (a) For the perimeter component, one can integrate over source
points x
a
∈ A
and the “shadow” x
c
∈ B
𝑖
(
x
a
)
cast by a discrete edge
𝑖
. (b) This also generalizes to the interior, but parameterizing and sampling the projected
boundary
B (
x
a
)
is diicult in general. (c) The local formulation instead evaluates a spherical integral at boundaries x
b
∈ 𝜕 A
without explicit consideration
of the neighboring vertices x
a
and x
c
. (d) The interior can be handled analogously but requires a dierent partition into an integral over surface positions
(orange) and tangential directions (blue). We propose a new local boundary integral that accounts for this component.
Bangaru et al
.
2020]. The latter methods must reparameterize
every ray at a signicant cost to trace auxiliary rays, and this
furthermore injects variance aecting the interior derivatives
(Figure 6). Our method can skip the projection in regions that
are far away from the boundary, and it preserves the statistical
quality of interior derivatives.
Our method’s main drawbacks all stem from the central projection
step: conversion of interior to boundary samples only works well
when there is sucient surface area to collect samples, and when
the projection itself is well-behaved. For example, a thin shape like
a blade of grass has a large boundary arclength to surface area
ratio and may not receive a sucient number of samples. Similarly,
a function that maps every interior point onto a xed silhouette
location is technically a projection, but not one that is conductive
to solving our problem. The rst point is a limitation of our method.
We address the second point by presenting high-quality projections
for diverse geometric representations in Section 5.
Before delving into the details of projection and guiding, we must
rst turn to a theoretical aspect of the problem that will enable
generalization to a wider array of shapes.
4 A local boundary integral
Zhang et al.’s [2020] path space formulation of dierentiable render-
ing decouples the eect of boundaries from their interior, enabling
targeted computation of each part using specialized methods. Their
equation (43) provides the starting point of our method.
This equation relates the change in pixel intensity
𝐼
due to visi-
bility changes with respect to a perturbation of a scene parameter
𝜃. Expressed with respect to the path segment (x
a
, x
c
), it states
𝜕𝐼
𝜕𝜃
=
∫
A
∫
B (x
a
)
𝐿
𝑖
(x
a
, x
c
) 𝐺 (x
a
, x
c
)𝑊
𝑖
(x
c
, x
a
)
(𝜕
𝜃
x
c
· n
c
) d𝑙(x
c
) d𝐴(x
a
). (3)
The integral is over segments
(
x
a
,
x
c
)
that make contact with surface
boundary along the way. The domain
B(
x
a
)
describes the “shadow”
of this boundary and is further composed of
B(
x
a
) = ∪
𝑖
B
𝑖
(
x
a
)
(one set
B
𝑖
per edge) when the scene consists of discrete geometry
(Figure 8a). The functions
𝐿
𝑖
and
𝑊
𝑖
refer to the incident radiance
and importance and can be computed using existing primal estima-
tors,
𝐺
is the standard geometric term [Veach 1997], and the inner
product measures the perpendicular speed of the “shadow” at x
c
(see Figure 8a).
In practice, it is simpler and far more ecient to build bound-
ary segments “outwards” starting from the tangential point x
b
, and
Zhang et al. therefore also propose a local formulation of the perime-
ter contribution (Figure 8c).
However, their derivation [Zhang et al
.
2020, Appendix 1] is
incomplete in the sense that the nal equation contains derivatives
arising from a ray tracing operation performed using autodi. The
presence of this complex step obscures simplication opportunities
that can bring this equation into its ultimate reduced form. These
simplications can also have implications on sampling eciency.
Our contribution here is two-fold: we re-derive the local formu-
lation of the perimeter to reveal its inherent simplicity, and we
propose the rst local formulation of the interior. Combining both
sources of derivatives leads to the following complete expression,
which provides the theoretical basis of our method.
𝜕𝐼
𝜕𝜃
=
∫
𝜕A
∫
S
2
𝐿
𝑑
(x
b
, 𝝎 ) 𝑊
𝑖
(x
b
, 𝝎 ) sin 𝜙 (𝜕
𝜃
x
b
·n
b
) d𝝎 d𝑙 (x
b
)
+
∫
A
∫
S
1
𝐿
𝑑
(x
b
, 𝜙)𝑊
𝑖
(x
b
, 𝜙) 𝜅 (𝜙) (𝜕
𝜃
x
b
·n
b
) d𝜙 d𝐴(x
b
). (4)
The derivation of this expression is somewhat technical, and we
refer the reader to the supplemental material for a detailed proof.
The term
𝐿
𝑑
stands for the radiance dierence between fore-
ground and background
𝐿
𝑑
(
x
b
, 𝜔) = 𝐿
𝑜
(
x
b
, 𝜔) − 𝐿
𝑖
(
x
b
, −𝜔)
. In the
perimeter term (top integral),
sin𝜙
refers to the angle between
𝝎
and the boundary tangent t
b
(Figure 8c). In the interior term, the
angle
𝜙 ∈ S
1
parameterizes all relevant quantities over tangential
directions at the surface position x
b
(Figure 8d), and
𝜅(𝜙)
denotes
the normal curvature. The following aspects are noteworthy:
1. Neither term involves the complex non-local domain B(x
a
).
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
Projective Sampling for Dierentiable Rendering of Geometry • 212:7
2.
Contrasting with Equation 3, the geometric term
𝐺
is absent in
both integrals, which means that variance arising from this fac-
tor
2
can be avoided in Monte Carlo methods. This is not obvious
when looking at Equation
(3)
, and the reference implementa-
tion of Zhang et al. [2020], e.g., tabulates
𝐺
within its guiding
distribution despite it canceling in subsequent steps.
3.
When the scene geometry is smooth and closed,
𝜕A = ∅
re-
moves the rst term. Polygonal meshes do not have curved
interior (𝜅 = 0) and hence do not require the second term.
4.
The sine (perimeter) and curvature (interior) terms resemble the
cosine factor in the rendering equation: The inuence of an edge
with tangent t
b
exerted in direction
𝝎
tends to zero as t
b
→ 𝝎
,
since the edge becomes invisible. Likewise, the curvature term
accounts for foreshortening in the mapping between silhouette
positions and scene positions observing this silhouette.
5 Method
With this background in place, we can now delve into the details of
interior and guided sampling. We discuss the projection operation
last, since it requires specialization to various geometric representa-
tions. This section already presents a number of results that examine
the quality of computed gradients to motivate algorithmic design
decisions. These are followed by more comprehensive end-to-end
optimization results in Section 6.
Input sampling strategies. Our method repurposes interior sam-
pling strategies into specialized strategies targeting the boundary.
But which distributions should be used as the input of this pro-
cess? In the setting of primal rendering, it is common to rely on
combinations of multiple strategies like BSDF and emitter sampling
via multiple importance sampling (MIS) [Veach 1997]. Figure 9 ex-
plores this possibility, indicating that the established intuition of
such combinations transfers to the dierential setting.
Embracing imperfection. Section 3 introduced the concept of a
path segment projection that maps an input segment x
a
→
x
b
onto
a nearby segment x
a
→
x
′
b
such that x
′
b
is located on a silhouette
observed by x
a
. It is important to note that the current requirements
for this operation are stricter than necessary, and that there could
be benets to relaxing them.
To see why, remember that we abandoned the idea of solving
the boundary integral with respect to a specic position x
a
, since
the density of projected samples was too expensive to characterize.
The projections merely guide a subsequent integration phase that
accounts for all surface positions x
a
simultaneously. Thus, it suces
to nd an approximate boundary segment x
′
a
→
x
′
b
close to x
a
→
x
b
(i.e., allowing for a slight shift of the origin x
a
as well). This added
exibility can be used in two ways:
First, the projection might fail.
In such cases, we snap x
′
b
to the closest local boundary (e.g. a trian-
gle edge) and select the tangential direction closest to the original
direction of x
a
→
x
b
. Further, if a projection requires root-nding,
2
To see why, consider a rod with 45-degree inclination casting a shadow from an
overhead directional light onto a ground plane. Prior work placed more samples on
the far end of the rod’s silhouette, which is undesirable as the contribution per unit
arclength is uniform.
Scene
Reference
BSDF sampling
Emier sampling
Combined
Fig. 9. Impact of the interior sampling distribution. This experiment
analyzes the derivative of a classic test scene by Veach with respect to a
horizontal translation. Closer investigation reveals familiar failure modes:
the projected BSDF sampling strategy fails to generate a suicient number
of silhouee samples to cover the rough reflection of the smallest sphere,
while emier sampling presents the opposite issue in the smooth reflection
of the largest sphere. Analogous to multiple importance sampling [Veach
1997] for primal rendering, it can be beneficial to project a mixture of both.
we terminate the iteration after a few steps without waiting for con-
vergence to machine precision. The former step reduces variance,
while the latter reduces the cost of creating the guiding distribution.
5.1 Guiding distribution
Given a set of projected samples, the next major step consists of
condensing this information into a representation that supports
ecient sampling and density evaluation.
5.1.1 Grid-based guiding. The simplest option is a 3D density grid,
which can be a good choice when the integrand exhibits sucient
smoothness. Zhang et al.’s [2020] path-space dierentiable rendering
(PSDR) method likewise relies on a grid.
One notable advantage of grids is that the projection can directly
accumulate projections without having to store the samples them-
selves, which can enable high-quality statistics accumulation in
memory-constrained environments. Conversely, a high resolution
3D grid can be very memory-intensive and tends to become a bot-
tleneck in complex scenes requiring a ne discretization to resolve
sparse features.
Figure 10 presents a rst set of results for grid-based guiding.
Each pair of rows shows the gradient with respect to a horizontal
translation in the rst row, followed by error visualizations compar-
ing to a nite dierence reference in the second row. We use two
dierent color scales throughout the paper: a standard red/gray/blue
scale for the error images in the lower rows (with gray representing
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
212:8 • Zhang et al.
Ground truth
Uniform
5x samples
Uniform
50x samples
Projective
1x sample
B
E C
F
F
D
-12
12
-1.5
1.5
-12
12
-1.5
1.5
-12
12
-1.5
1.5
-12
12
-1.5
1.5
-12
12
-1.5
1.5
L
1
= 0
.
237 L
2
= 4
.
635
L
1
= 0
.
062 L
2
= 0
.
400
L
1
= 0
.
052 L
2
= 0
.
179
L
1
= 0
.
108 L
2
= 0
.
931
L
1
= 0
.
034 L
2
= 0
.
202
L
1
= 0
.
028 L
2
= 0
.
111
L
1
= 0.142 L
2
= 2.137
L
1
= 0.044 L
2
= 0.259
L
1
= 0.038 L
2
= 0.133
L
1
= 0.438 L
2
= 1.557
L
1
= 0.202 L
2
= 0.726
L
1
= 0.165 L
2
= 0.481
L
1
= 0
.
186 L
2
= 1
.
434
L
1
= 0
.
064 L
2
= 0
.
308
L
1
= 0
.
056 L
2
= 0
.
197
Fig. 10. alitative and numerical evaluation of grid-based guiding com-
paring uniform [Zhang et al
.
2020] and projective sampling on meshes of
varying geometric complexity. Each pair of rows shows image gradients
followed by a visualization of the associated error. See Table 1 for timings.
zero error), and a more varied color scale for the gradient images due
to their large dynamic range. The
L
1
and
L
2
annotations denote
the gradient MAE and RMSE values compared to the reference. The
primal rendering is illustrative and uses dierent viewpoint and
lighting (Section 6 provides more details on the test scene setup).
All results in Figure 10 use a xed grid with 10000
×
100
×
100
voxels for
(𝑡, 𝜃, 𝜙)
. The only dierence is the number of accumu-
lated samples, and whether or not projective sampling was used.
The second column, labeled “uniform, 5x samples” shows the PSDR
baseline
3
using ve samples per voxel. The last column employs our
projection for triangle meshes and was generated using less time
than the competing methods (see Table 1 for precise timings). We
also use a lower number of samples to populate the grid given the
added cost of the projection. The middle result shows the quality
improvement that could be obtained by PSDR when using many
more samples to detect ne features within voxels, though this mul-
tiplies the computation time by a factor of nearly 10
×
. These results
show that projective sampling consistently outperforms PSDR by a
factor of 2.5-4.5
×
(MAE), which grows to a factor of up to 25
×
in
RMSE due to the presence of outliers.
5.1.2 Hierarchical guiding. Figure 7 (on page 5) previously visu-
alized the intricate structure of the integrand on boundary sample
space. For polygonal meshes, this 3D space consists of a single pa-
rameter
𝑡 ∈ [
0
,
1
]
, which parameterizes the set of mesh edges
𝜕A
,
along with a direction
(𝜃, 𝜙)
represented in spherical coordinates.
Sparse features in this domain indicate complexities of the input
scene, including:
1. BSDFs with narrow peaks.
2. Strongly peaked emitters including directional sources and en-
vironment maps featuring the sun.
3.
Second-order visibility eects, where silhouettes are themselves
occluded by surrounding geometry.
4. Discontinuities in the edge parameterization.
The last two points become more pronounced as the geometric
complexity increases. For example, although a sphere is a very
simple shape, its discretization into a polygonal mesh can generate
an extremely challenging integrand.
Yan et al. [2022] propose two ways to mitigate these diculties.
First, they recommend rearranging the edges into longer connected
chains, which reduces the number of discontinuities in the
𝑡
param-
eter. We have not incorporated this optimization into our method,
though it should be benecial here as well.
Secondly, they replace the grid with a set of kd-trees (one tree
per edge chain) to adapt to ne features of the 3D domain requir-
ing increased resolution. The construction of this tree is top-down
and resembles that of an adaptive quadrature rule: a recursive split-
ting criterion subdivides the current node until the integrand be-
comes suciently smooth. The failure mode of this approach is
also that of adaptive quadrature: a stopping criterion based on lo-
cal evaluations will sometimes fail to detect sharp peaks and stop
the subdivision prematurely.
Our hierarchical guiding uses a set of 100 octrees, each cover-
ing an equal-sized interval of the rst axis to handle heterogeneity
3
We use our PSDR reimplementation, which outperforms the reference code. We also
adapted it to use our improved local boundary integral formulation from Equation 4.
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
Projective Sampling for Dierentiable Rendering of Geometry • 212:9
arising from the parameterization of edges, bers, and the concate-
nation of multiple objects into a shared boundary sample space. The
particular type of hierarchy is not a critical factor (we merely chose
octrees since their construction is easy to parallelize on the GPU).
The main dierence lies in how this hierarchy is constructed. Our
hierarchy construction is bottom-up, starting from a set of projected
samples concentrated at the sparse features in boundary sample
space. We then simply subdivide each node until reaching a max-
imum depth (9), or until the node contains one or fewer samples,
which yields a partition matching the spatio-directional structure
of the integrand. In other words, unlike grid-based guiding, where
a projection was used to accumulate density into the voxels of a 3D
grid, we now employ the projection to cheaply construct a high-
delity adaptive space partition. To estimate the actual value of the
boundary integrand, we additionally draw 32 uniformly distributed
samples within each leaf node, which does not require projections.
In our experiments, this results in an average of
∼
6
.
4
×
10
7
samples
for a base budget of 10
8
projections, which roughly doubles the
initial sample budget.
Figure 11 compares this scheme to the method of Yan et al. [2022]
(“Adaptive quadrature”) at equal time (conservatively, see Table 2
for precise timing values) and a grid-based baseline using projec-
tions. The experiment demonstrates that the octree signicantly
outperforms the grid-based baseline. The method of Yan et al. [2022]
produces spatially non-uniform convergence with low error in some
regions and signicant outliers in others, causing a
∼
190
×
RMSE
dierence in the Neptune experiment.
5.2 Projection
We now discuss the operations needed to support a particular geo-
metric representation:
1.
A projection that maps a ray segment x
a
→
x
b
onto a nearby
segment x
′
a
→ x
′
b
tangential at x
′
b
.
2.
A parameterization of the perimeter (if present), and a parame-
terization of the interior (for curved geometry).
5.2.1 Spheres. This is by far the simplest case, which can be help-
ful to test new implementations of projective sampling. Listing 1
provides pseudocode for the projection. Being smooth and closed,
spheres only need a surface mapping
(we use spherical coordinates).
5.2.2 Triangle meshes. We employ two dierent projection strate-
gies for meshes: Jump and Walk.
The Jump projection is a Newton-style iteration based on a local
linear approximation. We model the neighborhood of a surface posi-
tion p using the parameterization
˜
p(𝑢, 𝑣) = p + 𝑢𝜕
𝑢
p + 𝑣𝜕
𝑣
p
with an
interpolated normal
˜
N(𝑢, 𝑣) = N + 𝑢𝜕
𝑢
N + 𝑣𝜕
𝑣
N
. Under the assump-
tion of a xed viewing direction, the silhouette of this approximation
is a line in the parameter space, allowing us to analytically nd a
solution
(𝑢
′
, 𝑣
′
)
. From the inferred silhouette point, we subsequently
trace a perpendicular ray
(
˜
p(𝑢
′
, 𝑣
′
) + 𝜀
˜
N(𝑢
′
, 𝑣
′
), −
˜
N(𝑢
′
, 𝑣
′
))
to nd
the next intersection, at which point the procedure could be repeated.
The Walk projection is a greedy search strategy. We visit the three
neighboring triangles and compute the angles between their normals
and the viewing direction. Walk aims to gradually increase this
Ground truth
Adaptive quadrature
Projective (grid)
Projective (octree)
B TB
N
-12
12
-1.5
1.5
-12
12
-1.5
1.5
-12
12
-1.5
1.5
L
1
= 0.173 L
2
= 4.391
L
1
= 0.047 L
2
= 0.230
L
1
= 0.030 L
2
= 0.198
L
1
= 0
.
118 L
2
= 3
.
697
L
1
= 0
.
032 L
2
= 0
.
132
L
1
= 0
.
021 L
2
= 0
.
084
L
1
= 0
.
292 L
2
= 21
.
971
L
1
= 0
.
071 L
2
= 0
.
244
L
1
= 0
.
038 L
2
= 0
.
115
Fig. 11. Guiding data structure comparison including adaptive quadra-
ture [Yan et al
.
2022] (with edge sorting, MIS, and parameters set to use
all available GPU memory), a uniform grid with projections, and an octree
initialized with projected samples. See Table 2 for timings.
Algorithm 1. Projection for spheres. This code fragment demonstrates
the interface of the projection operation for a simple geometric primitive.
class Sphere():
def project(self, segment: Segment):
O = segment.a # viewpoint
P = segment.b # point on the sphere
N = (P - self.center) / self.radius # the normal at P
CO = O - self.center # center to viewpoint vector
CO_dir = (O - self.center) / norm(CO)
proj_dir = normalize(N - dot(N, CO_dir) * CO_dir)
# phi: the angle from C to the boundary relative to O
sin_phi = r / norm(CO)
B = CO_dir * sin_phi + proj_dir * sqrt(1 - sin_phi**2)
return Segment(shape=self, a=O, b=B)
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
212:10 • Zhang et al.
angle until reaching a perpendicular (90
◦
) or back-facing triangle.
However, consistently walking to the neighboring triangle with the
largest angle tends to attract too many projections on certain mesh
edges. Instead, we randomly pick between the two neighbors with
the largest angle values, using them as discrete probabilities.
Both strategies only rely on local dierential information. Jump
is more aggressive and can quickly bypass plateaus and highly
tesselated regions, while the smaller steps of Walk robustly detect
nearby silhouettes even on bumpy geometry, where the extrapolated
local model can be deceptive. Walk does not require ray tracing,
which makes individual steps much faster compared to Jump.
In practice, we adopt a hybrid strategy: Walk for 30 steps and
Jump once if no silhouette was found, followed by an additional 30
Walk steps. This strategy only requires a single ray tracing step
and produces high-quality projections on all meshes presented in
this paper. If the projection of the segment x
a
→
x
b
fails to nd
a silhouette x
′
b
, we keep the tentative position x
′
b
and generate a
tangential segment x
′
a
→
x
′
b
closest to the original direction of
x
a
→
x
b
, as outlined at the beginning of Section 5. While more
sophisticated projections could likely be pursued to handle various
corner cases, we nd that even this simple scheme already enables
high-quality guiding.
5.3 Fiber curves
We model bers using a parametric base curve C
(𝑣)
and radius
𝑟 (𝑣)
,
which are both given by cubic B-spline interpolants (Figure 12).
The cross-section of the surface for a xed value of
𝑣
yields a cir-
cle with center C
(𝑣)
, radius
𝑟 (𝑣)
, and normal C
′
(𝑣)
. Assigning an
azimuth angle parameter
𝑢
to this circle yields a
𝐶
1
-continuous
surface M
(𝑢, 𝑣)
. We ignore the curve endpoints (these can, e.g., be
closed using spheres), in which case only the curved interior part
of the local boundary integral in Equation (4) matters.
Projection. Given a viewpoint O and surface position P
=
M
(𝑢
0
, 𝑣
0
)
,
we x
𝑣
0
and nd the value of
𝑢
that yields a silhouette projection,
i.e.,
⟨
M
(𝑢, 𝑣
0
)−
O
,
n
(𝑢, 𝑣
0
)⟩ =
0, where n is the parameterized surface
normal. This equation can be expanded into the form
𝐴 cos
2
𝑢 + 𝐵 cos𝑢 sin𝑢 + 𝐶 cos 𝑢 + 𝐷 sin𝑢 + 𝐸 = 0,
where
𝐴, 𝐵, 𝐶, 𝐷, 𝐸
are
𝑢
-independent constants that can be com-
puted analytically. The equation has an analytic solution when the
ber radius is constant. Otherwise, we run 20 bisection iterations,
which is fast and robust.
5.4 Implicit functions
We also experimented with signed distance functions (SDFs) rep-
resented using a grid-based trilinear interpolant (Figure 12). This
representation has two advantages:
1.
Speed. Intersections can be computed with the help of a mesh
proxy and analytic solutions within voxels, which avoids costly
sphere tracing steps [Söderlund et al. 2022].
2.
Simplicity. The low-order representation simplies the map-
ping from R
2
to the SDF surface.
The trilinear representation also has a clear disadvantage: geometric
normals are discontinuous across voxel boundaries, which means
that both interior and perimeter terms of Equation
(4)
must be
Fiber curve
Signed distance function
Perimeter
Interior
Fig. 12. Smooth geometry. Our new local formulation of the boundary
derivative (Equation 4) enables dierentiable rendering of smooth geometry,
such as cylindrical fibers based on Bézier curves (le) and implicitly defined
surfaces represented using a signed distance function (right). The laer case
involves derivatives arising from the curved interior and potential normal
discontinuities at voxel perimeters.
considered
4
. Note that we do not use the signed distance property of
the SDF representation, so much of the following could in principle
generalize to other types of implicitly dened surfaces.
Parameterization. To dene functions parameterizing the interior
and perimeter, we rst make the following observations:
1.
Interior. Given the trilinear interpolation scheme, a given
𝑥, 𝑦
coordinate within a voxel has at most one root along the
𝑧
axis.
2.
Perimeter. Similarly, the surface intersection with a voxel face
can be uniquely parameterized by the face index and a perpen-
dicular coordinate, which can be used to create a attened 1D
mapping across all possible perimeter curves reminiscent of the
mesh case. Two additional dimensions represent the direction.
We use these properties to create a globally discontinuous per-voxel
mapping. This mapping depends on the choice of a dimension used
to parameterize the curve/surface, which is unstable when the cho-
sen dimension is nearly perpendicular. We assign the most numeri-
cally stable dimension to each voxel based on the SDF gradient.
We did not realize a projection operation for SDFs but believe
that such a step could be added along similar lines (related opera-
tions were also explored in prior work [Bremer and Hughes 1998]).
Presented results for SDFs therefore use a grid data structure with
uniform initialization.
6 Results
We implemented our method on top of Mitsuba 3’s [Jakob et al
.
2022b]
cuda_ad_rgb
backend, using the underlying Dr.Jit [Jakob
et al
.
2022a] framework for forward/reverse-mode AD and GPU
kernel compilation. All experiments ran on an AMD Ryzen 3970X
workstation with an NVIDIA RTX 3090 graphics card.
Test scene. The test scene used to evaluate methods throughout
the paper places the shape in front of a uniform area emitter. The
camera views their reection on a conductive plane with isotropic
roughness (Trowbridge and Reitz [1975],
𝛼 =
0
.
005). This ensures
4
In practice, such a low-order interpolation would be combined with a smooth shading
normal approximation to suppress distracting seams at voxel boundaries. But this only
xes the smooth part of the dierentiation problem–the visibility-related discontinuous
part must consider the true geometric normal.
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
Projective Sampling for Dierentiable Rendering of Geometry • 212:11
(a)
(b)
High-frequency Low-frequency
Fig. 13. Shadows and reflections produce a qualitatively similar
boundary integral. (a) Primal rendering of scenes with shadows or reflec-
tions. (b) The corresponding boundary sample space. This figure demon-
strates their similarity for both high-frequency (hard shadows, sharp re-
flections) and low-frequency (so shadows, rough reflection) cases. Conse-
quently, most of our gradient comparisons focus on the reflection case.
that the derivative image only contains contributions from indirectly
observed discontinuities, which is the main focus of this work.
We prefer sharp reections to investigate the performance of
dierent methods: a soft material or shadow would blur the inte-
grand in boundary sample space, potentially concealing inaccuracies
or aws in the method. As illustrated in Figure 13, shadows and
reections can produce a qualitatively similar boundary sample
space. Instead of a reective plane, one could therefore equivalently
employ a diuse receiver with a directionally peaked emitter.
Initialization time. Tables 1 and 2 list the time needed to initialize
the guiding distributions for the previously discussed experiments
in Figures 10 and 11. The runtime of the dierential rendering phase
is not noticeably impacted by the choice of guiding representation.
In these gradient comparison experiments, derivative images were
rendered using 256 samples per pixel (roughly
∼
0
.
8 second of ren-
dering time), a number typically higher than necessary for inverse
rendering tasks. We use this number to render ground truth quality
derivative images as a validation of our method.
Fiber curve. Figure 14 showcases the inuence of the guiding
distribution for dierentiable rendering of curves. It highlights the
eectiveness of projective sampling, which outperforms uniform
sampling even when the latter uses 50
×
more samples. We nd
that the octree representation can greatly reduce variance in com-
plex cases like the Knitting Yarn experiment, where second-order
visibility eects (occlusion of silhouettes) dominate.
End-to-end optimizations. Figure 15 demonstrates the use of our
method as part of several complete reconstruction tasks. In the
rst two experiments, the optimization must infer the shape from
a single view depicting the object and two shadows on a curved
diuse wall. We use the method of Nicolet et al
.
[2021] to smoothly
evolve the triangle mesh and periodically re-mesh the geometry to
a ner resolution to add more degrees of freedom. Listing 2 shows
the high-level pseudocode of the optimization algorithm.
The nal row of Figure 15 constructs a Magic Lens as in the work
of Papas et al
.
[2012]. The objective of this task is to optimize the
Ground truth
Uniform (grid)
5x samples
K
C
W
K Y
S
Uniform (grid)
50x samples
Projective (grid)
1x sample
Projective (octree)
2x samples
-10
10
-2.0
2.0
-10
10
-2.0
2.0
-10
10
-2.0
2.0
-10
10
-2.0
2.0
-10
10
-4.0
4.0
L
1
= 0.278 L
2
= 1.838
L
1
= 0.102 L
2
= 0.505
L
1
= 0.084 L
2
= 0.407
L
1
= 0.036 L
2
= 0.126
L
1
= 1.515 L
2
= 7.851
L
1
= 0.696 L
2
= 2.180
L
1
= 0.424 L
2
= 1.096
L
1
= 0.133 L
2
= 0.285
L
1
= 0.926 L
2
= 5.498
L
1
= 0.422 L
2
= 1.756
L
1
= 0.277 L
2
= 0.845
L
1
= 0.079 L
2
= 0.187
L
1
= 0.679 L
2
= 4.954
L
1
= 0.297 L
2
= 1.170
L
1
= 0.173 L
2
= 0.567
L
1
= 0.056 L
2
= 0.148
L
1
= 3.725 L
2
= 46.649
L
1
= 2.783 L
2
= 21.269
L
1
= 2.195 L
2
= 5.890
L
1
= 0.519 L
2
= 1.967
Fig. 14. Experimental validation of fiber derivatives analogous to previous
experiments. See Table 3 for timings.
surface of a dielectric plate so that an object located behind the
dielectric appears completely transformed into a dierent object
when observed from a specic viewpoint. We place a hat behind the
lens and optimize the lens so that the refraction resembles a bunny.
In addition to lens’ surface, we also optimize the hat’s position. A
sketch of this setup is shown in the leftmost column.
Inuence of gradient quality. Figure 16 illustrates the inuence of
gradient quality in a challenging reconstruction task. The scene is
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
212:12 • Zhang et al.
Reference geometry
Final geometry
Initial geometry
H
M H
Initial rendering Final rendering Target rendering
E
Scene
Fig. 15. Optimization tasks involving discontinuous derivatives. The first two experiments we reconstruct the shape of an object from a single-view
capture. The last row optimizes the position of a hat and the height field of a refractive plate with a small amount of roughness (Beckmann
𝛼 =
0
.
01) to form
an image of a bunny. (The supplemental video contains animated versions of the first two experiments.)
3.0
-3.0
Projective (grid)
Uniform (grid)
3x samples
20x samples
50x samples120x samples
1x sample
Iteration 40 100 160 220 280 340 400
Reference
Optimization view
Geometry (novel view/relit)
0.5s/iter
0.7s/iter
4.8s/iter
11.1s/iter
26.2s/iter
Forward derivative
Ground truthGround truth Proj 1x
Unif 3x
Unif 20x
Unif 50x Unif 120x
Fig. 16. The eect of gradient quality on geometric reconstruction. We reconstruct a triangle mesh from a single reference view with multiple shadows on
a curved floor (upper le), comparing projective and uniform sampling for dierent iteration and sample counts. Timing values in the rightmost column quantify
discontinuity-related costs while ignoring the shared overheads of primal rendering, preconditioned gradient descent [Nicolet et al. 2021], and re-meshing.
designed so that the shadows reveal a sucient amount of informa-
tion to disambiguate the shape and in principle enable convergence
using gradient descent. We reuse the optimization strategy from
Figure 15 and perform grid-based guiding with an equal resolu-
tion in all ve optimization runs. The only dierence lies in the
initialization of the grid and the resulting gradient quality. We also
visualize forward gradients for a horizontal displacement to explain
the dierences in convergence.
Optimization using low-quality gradients (“uniform, 3x samples”)
stagnates, while estimates with lower variance (“uniform, 20x/50x/120x
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
Projective Sampling for Dierentiable Rendering of Geometry • 212:13
Ground truth Perimeter+Interior Perimeter Interior
T BB TD C
FC
Geometry
-10
10
-10
10
-10
10
-10
10
-10
10
Fig. 17. This visualization shows visibility-induced derivatives of signed
distance functions (SDFs), separating the contributions from perimeter and
interior. These examples all use a 128
3
SDF grid with trilinear interpolation.
samples”) lead to progressively better convergence in an equal-
iteration comparison (though this comes at a greatly increased com-
putational cost). Our method (“projective, 1x sample”) achieves a
good balance between gradient quality and computation time.
Signed distance functions. Figure 17 validates the use of our local
boundary integral formulation (Equation 4) for SDFs. We separate
the derivative contributions of perimeter and interior, which super-
impose to produce a gradient matching the ground truth.
Importance of indirect derivatives. Figure 18 shows the impact of
indirectly observed discontinuities in a single-view mesh reconstruc-
tion. Such indirect eects are an inherent property of essentially
any realistic image. When self-shadowing is not taken into account
during the dierentiation, the optimizer lacks the information that
the shadow near the nose is caused by the eyebrow. It attempts
to conceal the undesired shadow by distorting the nose to cover
it. With self-shadowing derivatives computed by our method, the
optimizer instead subtly adjusts the eyebrow. Other aspects (e.g.
nose height) also improve, since shadows reduce ambiguities in the
challenging single-view setting.
7 Conclusion
A central discovery of path-space dierentiable rendering [Zhang
et al
.
2020] was that the interior and boundary derivatives decouple,
allowing each part to be handled independently. In contrast, our
work identies substantial synergies between these two, indicat-
ing that a full separation is generally inadvisable. Our experiments
demonstrate that the boundary term can greatly benet from infor-
mation gathered during the interior term’s simulation.
With shadow gradients
No shadow gradients
Reference
Optimization view
Geometry (novel view/re-lit)
Gradients
Fig. 18. Single-view reconstruction with and without derivatives due to
self-shadowing. Accounting for them resolves ambiguity in this information-
constrained scenario. The gradient image in the last row (with respect to a
horizontal translation) illustrates the missing derivatives. (The supplemental
video contains an animated version of this figure.)
This idea leads to a modular framework of projections and pa-
rameterizations that tie into the established architecture of contem-
porary rendering systems. Possible constructions include hopping
from triangle to triangle, jumping over larger distances, and the nu-
merical solution of nonlinear equations. We nd this an interesting
design space and believe that future work could also specically
target linear elements to overcome the dicult “blade of grass” case
mentioned earlier.
The discontinuous nature of boundary sample space remains a
signicant source of ineciency. Both prior work [Yan et al
.
2022]
and this article demonstrate that edge permutations and adaptive
discretizations can mitigate this problem to a certain extent, but
these are essentially just band-aids. In general, the representation is
intrinsically smooth, but this property is lost in the mapping onto
a parameter space. Truly smooth global parameterizations of this
challenging domain are an interesting topic for future work.
Acknowledgments
This project has received funding from the European Research Coun-
cil (ERC) under the European Union’s Horizon 2020 research and
innovation program (grant agreement No 948846).
References
Sai Praveen Bangaru, Michaël Gharbi, Fujun Luan, Tzu-Mao Li, Kalyan Sunkavalli,
Milos Hasan, Sai Bi, Zexiang Xu, Gilbert Bernstein, and Fredo Durand. 2022. Dier-
entiable rendering of neural sdfs through reparameterization. In SIGGRAPH Asia
2022 Conference Papers. 1–9.
Sai Praveen Bangaru, Tzu-Mao Li, and Frédo Durand. 2020. Unbiased warped-area
sampling for dierentiable rendering. ACM Trans. Graph. 39, 6 (2020), 1–18.
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
212:14 • Zhang et al.
Algorithm 2. Pseudocode of the optimization loop and a vectorized reverse-
mode implementation of the rendering function.
def optimize()
for l in range(N_iter):
loss = l2(render(scene), reference)
loss.backward()
optimizer.step(scene.grad)
@backward
def render(scene, grad_img):
rays = scene.sensor.sample_rays()
# Propagate derivative of the continuous part with PRB and
# return path segments
img_cont, segments = prb_backward(grad_img)
shape = segments.shape
# Project path segments, returns an array of new segments
# and failure markers.
proj_segments = shape.project(segments)
# Convert into a 3D point in boundary sample space
boundary_sample = shape.map_to_sample(proj_segments)
# Evaluate the boundary integrand without motion
value = eval_integrand(proj_segments)
# Initialize the guiding data structure (the octree may
# call eval_integrand internally)
distr = guiding_octree(boundary_sample, value)
# Draw samples from the guiding distribution
boundary_sample, pdf = distr.sample_pdf()
boundary_segment = shape.map_to_segment(boundary_sample)
# Particle tracer style render starting from boundary segment
img_boundary = render_boundary(boundary_segment, pdf)
img_combined = img_boundary + img_cont
img_combined.backward(grad_img)
Sai Praveen Bangaru, Jesse Michel, Kevin Mu, Gilbert Bernstein, Tzu-Mao Li, and
Jonathan Ragan-Kelley. 2021. Systematically dierentiating parametric discontinu-
ities. ACM Trans. Graph. 40, 4 (2021), 1–18.
Thomas Booth. 2007. Unbiased Monte Carlo Estimation of the Reciprocal of an Integral.
Nuclear Science and Engineering 156 (07 2007), 403–407.
David Bremer and John F. Hughes. 1998. Rapid Approximate Silhouette Rendering of
Implicit Surfaces. In Proceesings of Implicit Surfaces 98.
Forrester Cole, Kyle Genova, Avneesh Sud, Daniel Vlasic, and Zhoutong Zhang. 2021.
Dierentiable surface rendering via non-dierentiable sampling. In Proceedings of
the IEEE/CVF International Conference on Computer Vision. 6088–6097.
Andreas Griewank and Andrea Walther. 2008. Evaluating derivatives: principles and
techniques of algorithmic dierentiation. Vol. 105. SIAM.
Wenzel Jakob, Sébastien Speierer, Nicolas Roussel, Merlin Nimier-David, Delio Vicini,
Tizian Zeltner, Baptiste Nicolet, Miguel Crespo, Vincent Leroy, and Ziyi Zhang.
2022b. Mitsuba 3 renderer. https://mitsuba-renderer.org.
Wenzel Jakob, Sébastien Speierer, Nicolas Roussel, and Delio Vicini. 2022a. Dr.Jit:
A Just-In-Time Compiler for Dierentiable Rendering. Transactions on Graphics
(Proceedings of SIGGRAPH) 41, 4 (July 2022).
Hiroharu Kato, Yoshitaka Ushiku, and Tatsuya Harada. 2018. Neural 3D Mesh Renderer.
In IEEE Conference on Computer Vision and Pattern Recognition (CVPR).
Samuli Laine, Janne Hellsten, Tero Karras, Yeongho Seol, Jaakko Lehtinen, and Timo
Aila. 2020. Modular Primitives for High-Performance Dierentiable Rendering.
Transactions on Graphics (Proceedings of SIGGRAPH) 39, 6 (2020).
Tzu-Mao Li, Miika Aittala, Frédo Durand, and Jaakko Lehtinen. 2018. Dierentiable
monte carlo ray tracing through edge sampling. ACM Trans. Graph. 37, 6 (2018),
1–11.
Tzu-Mao Li, Michal Lukáč, Michaël Gharbi, and Jonathan Ragan-Kelley. 2020. Dieren-
tiable Vector Graphics Rasterization for Editing and Learning. ACM Trans. Graph.
39, 6 (nov 2020).
Table 1. Guiding initialization time (in seconds) for Figure 10.
EmptyCube Bunny Fertility Filigree Dragon
Uniform (5x) 0.52 0.49 0.52 0.46 0.47
Uniform (50x) 4.30 4.05 4.36 3.71 3.80
Projective 0.41 0.29 0.30 0.30 0.29
Table 2. Guiding initialization time (in seconds) for Figure 11.
BumpyTorus Botijo Neptune
Adaptive quadrature 1.12 1.11 1.11
Projective (grid) 0.39 0.26 0.26
Projective (octree) 0.73 0.48 0.46
Table 3. Guiding initialization time (in seconds) for Figure 14.
Knot Coil Winding Sausage KnittingYarn
Uniform (grid, 5x) 0.75 0.68 1.92 0.98 0.59
Uniform (grid, 50x) 7.22 6.52 18.98 9.68 5.48
Projective (grid) 0.13 0.11 0.30 0.12 0.14
Projective (octree) 0.26 0.26 0.64 0.30 0.33
Shichen Liu, Tianye Li, Weikai Chen, and Hao Li. 2019. Soft Rasterizer: A Dierentiable
Renderer for Image-based 3D Reasoning. The IEEE International Conference on
Computer Vision (ICCV) (Oct. 2019).
Matthew M. Loper and Michael J. Black. 2014. OpenDR: An Approximate Dierentiable
Renderer. In European Conference on Computer Vision (ECCV). Springer, 154–169.
Guillaume Loubet, Nicolas Holzschuch, and Wenzel Jakob. 2019. Reparameterizing
discontinuous integrands for dierentiable rendering. ACM Trans. Graph. 38, 6
(2019), 1–14.
Baptiste Nicolet, Alec Jacobson, and Wenzel Jakob. 2021. Large Steps in Inverse Ren-
dering of Geometry. ACM Transactions on Graphics (Proceedings of SIGGRAPH Asia)
40, 6 (Dec. 2021).
Merlin Nimier-David, Sébastien Speierer, Benoît Ruiz, and Wenzel Jakob. 2020. Radiative
Backpropagation: An Adjoint Method for Lightning-Fast Dierentiable Rendering.
Transactions on Graphics (Proceedings of SIGGRAPH) 39, 4 (July 2020).
Marios Papas, Thomas Houit, Derek Nowrouzezahrai, Markus Gross, and Wojciech
Jarosz. 2012. The Magic Lens: Refractive Steganography. ACM Transactions on
Graphics (Proceedings of SIGGRAPH Asia) 31, 6 (Nov. 2012).
Felix Petersen, Bastian Goldluecke, Christian Borgelt, and Oliver Deussen. 2022. GenDR:
A Generalized Dierentiable Renderer. In IEEE Conference on Computer Vision and
Pattern Recognition (CVPR). 4002–4011.
Osborne Reynolds. 1903. The sub-mechanics of the universe. Vol. 3. University Press.
Herman Hansson Söderlund, Alex Evans, and Tomas Akenine-Möller. 2022. Ray Tracing
of Signed Distance Function Grids. Journal of Computer Graphics Techniques Vol 11,
3 (2022).
T. S. Trowbridge and K. P. Reitz. 1975. Average irregularity representation of a rough
surface for ray reection. J. Opt. Soc. Am. 65, 5 (May 1975), 531–536.
Eric Veach. 1997. Robust Monte Carlo Methods for Light Transport Simulation. Ph. D.
Dissertation. Stanford University.
Eric Veach and Leonidas J. Guibas. 1995. Optimally Combining Sampling Techniques for
Monte Carlo Rendering. In Proceedings of the 22nd Annual Conference on Computer
Graphics and Interactive Techniques. Association for Computing Machinery.
Delio Vicini, Sébastien Speierer, and Wenzel Jakob. 2021. Path Replay Backpropagation:
Dierentiating Light Paths using Constant Memory and Linear Time. Transactions
on Graphics (Proceedings of SIGGRAPH) 40, 4 (Aug. 2021).
Delio Vicini, Sébastien Speierer, and Wenzel Jakob. 2022. Dierentiable Signed Distance
Function Rendering. Transactions on Graphics (Proceedings of SIGGRAPH) 41, 4 (July
2022).
Kai Yan, Christoph Lassner, Brian Budge, Zhao Dong, and Shuang Zhao. 2022. Ecient
estimation of boundary integrals for path-space dierentiable rendering. ACM
Trans. Graph. 41, 4 (2022), 1–13.
Yuting Yang, Connelly Barnes, Andrew Adams, and Adam Finkelstein. 2022. A
𝛿
:
autodi for discontinuous programs-applied to shaders. ACM Trans. Graph. 41, 4
(2022), 1–24.
Cheng Zhang, Bailey Miller, Kan Yan, Ioannis Gkioulekas, and Shuang Zhao. 2020.
Path-space dierentiable rendering. ACM Transactions on Graphics 39, 4 (2020).
Yang Zhou, Lifan Wu, Ravi Ramamoorthi, and Ling-Qi Yan. 2021. Vectorization for
Fast, Analytic, and Dierentiable Visibility. ACM Transactions on Graphics 40, 3
(July 2021).
ACM Trans. Graph., Vol. 42, No. 6, Article 212. Publication date: December 2023.
