Path-Space Differentiable Rendering of Implicit Surfaces
Siwei Zhou

Youngha Chang

Nobuhiko Mukai

siwei.zhou@ist.osaka-u.ac.jp
Osaka University
Japanychang@tcu.ac.jp
Tokyo City University
Japan

nmukai@tcu.ac.jp
Tokyo City University
Japan

Hiroaki Santo

Fumio Okura

Yasuyuki Matsushitasanto.hiroaki@ist.osaka-u.ac.jp
Osaka UniversityJapanokura@ist.osaka-u.ac.jp
Osaka University
Japanyasumat@ist.osaka-u.ac.jp
Osaka University
Japan

Shuang Zhao
shz@ics.uci.edu
UC Irvine & NVIDIA
USA
Positive

Negative(a) Ordinary(b) Derivative

Figure 1: We apply path-space differentiable rendering to implicit surfaces. Specifically, we derive the necessary mathematical
tools and then propose Monte Carlo estimators of the boundary integrals. In this example, we showcase our method can handle
the boundary integrals even with complex shapes.

ABSTRACT
Physics-based differentiable rendering is a key ingredient for integrating forward rendering into probabilistic inference and machinelearning pipelines. As a state-of-the-art formulation for differentiable rendering, differential path integrals have enabled the development of efficient Monte Carlo estimators for both interior
and boundary integrals. Unfortunately, this formulation has been
designed mostly for explicit geometries like polygonal meshes.
In this paper, we generalize the theory of differential path integrals to support implicit geometries like level sets and signeddistance functions (SDFs). In addition, we introduce new Monte

Carlo estimators for efficiently sampling discontinuity boundaries
that are also implicitly specified. We demonstrate the effectivenessof our theory and algorithms using several differentiable-renderingand inverse-rendering examples.

CCS CONCEPTS• Computing methodologies → Rendering.

KEYWORDS
Differentiable rendering, stochastic sampling

Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than theauthor(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or
republish, to post on servers or to redistribute to lists, requires prior specific permission
and/or a fee. Request permissions from permissions@acm.org.
SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA© 2024 Copyright held by the owner/author(s). Publication rights licensed to ACM.
ACM ISBN 979-8-4007-0525-0/24/07
https://doi.org/10.1145/3641519.3657473

ACM Reference Format:
Siwei Zhou, Youngha Chang, Nobuhiko Mukai, Hiroaki Santo, Fumio Okura,
Yasuyuki Matsushita, and Shuang Zhao. 2024. Path-Space Differentiable
Rendering of Implicit Surfaces. In Special Interest Group on Computer Graphics and Interactive Techniques Conference Conference Papers ’24 (SIGGRAPH
Conference Papers ’24), July 27–August 01, 2024, Denver, CO, USA. ACM, New
York, NY, USA, 11 pages. https://doi.org/10.1145/3641519.3657473

SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA

1

INTRODUCTION

Physics-based rendering is the core technology for high quality
computer-generated imagery such as those seen in feature films.
Physics-based differentiable rendering brings this technology toother fields such as machine learning. The main object being studied
in this field is the boundary integrals. Those boundary integrals
occur when shape is to be differentiated with respect to. Early
research focused on deriving these boundary integrals in different
contexts [Li et al. 2018; Zhang et al. 2019; Loubet et al. 2019; Zhang
et al. 2020, 2021]. Recent works focus more on how to efficiently
estimate those boundary integrals [Yan et al. 2022; Xu et al. 2023;
Zhang et al. 2023].
A common assumption in those researches is the representationof the shape being mesh or parametric surfaces. While mesh is
a popular shape representation, another popular one is implicit
surface. An implicit surface is defined by the zero-level set of somescalar-valued function, i.e., where the function takes its value zero.
In practice, there are two families of Monte Carlo estimators tohandle the boundary integrals, namely, reparameterization [Loubet
et al. 2019; Bangaru et al. 2020; Xu et al. 2023] and edge-sampling [Liet al. 2018; Zhang et al. 2020, 2023]. Reparameterization replaces
the boundary integrals with equivalent interior ones and, thus,
avoids explicitly sampling discontinuities boundaries. In addition,these techniques can be adopted for implicit geometries relatively
easily [Vicini et al. 2022; Bangaru et al. 2022]. On the other hand, thereparameterized interior integrals remain difficult to importancesample, which can cause the reparameterization-based methods tosuffer from higher variance.
On the other hand, edge-sampling-based methods directly estimate the boundary integrals using Monte Carlo integration. When
importance sampled properly [Zhang et al. 2023], they tend to
produce lower-variance estimates than the reparameterizationbased methods. Unfortunately, state-of-the-art techniques in thiscategory—such as the path-space formulation by Zhang et al. [2020]—
have been mostly derived for explicit geometries. For example,
Monte Carlo sampling over silhouettes of implicit surfaces has notbeen addressed by prior methods.
In this paper, we bridge this gap by introducing a path-space
theory for differentiable rendering of implicit surfaces. Concretely,
Our contributions include:• Generalizing the path-space differentiable rendering framework[Zhang et al. 2020] to include implicit surfaces (§4);• Introducing Monte Carlo estimators for the boundary integralswith implicit surfaces (§5).
We demonstrate the effectiveness of our method using several
synthetic differentiable-rendering and inverse-rendering examples.

2

RELATED WORKS

Forward rendering of implicit surfaces. Implicit surface rendering or isosurface rendering can be used to visualized volumetricdata such as those obtained by CT or MRI. One way to render
implicit surfaces is to triangulated the isosurface using Marching
Cubes [Lorensen and Cline 1998] and then render the triangularmesh using regular rendering methods. One can also render theZhou, Chang, Mukai, Santo, Okura, Matsushita, and Zhao

Table 1: Commonly used symbols in this paper. The rightmost column indicates 𝜋-dependency.
Symbol Definition
M
B
X(·, 𝜋 )𝜋 -dep.

The union of all surfaces in the scene
Reference surface
Differentiable one-to-one mappingtransforming B to M (𝜋 )

Yes
No
Yes

ˆ
𝛀
ˆ
𝜕𝛀
𝑓ˆ

material path spacematerial boundary path spacematerial measurement function

No
Yes
Yes

ΔB𝐾
𝒏𝜕

Jump discontinuities of 𝑓ˆ with respect to 𝒑 𝐾
Normal of discontinuity curve

Yes
Yes

𝒗
𝑉𝜕

Velocity field over B defined in Eq. (15)Curve normal velocity defined in Eq. (8)

Yes
Yes

implicit surface from point samples generated by particle repulsion [Witkin and Heckbert 1994]. Yet another way is to use ray tracing [Barr 1986; Levoy 1988]. Ray tracing requires a ray-surface intersection routine, examples of which include sphere tracing [Hart1996] and interval analysis [Flórez et al. 2008]. Besides the above
surface rendering methods, volume rendering [Drebin et al. 1988]
can also be used to render implcit surfaces.
Physics-based differentiable rendering. Physics-based differentiable rendering correctly accounts for derivatives of rendered images with respect to shapes. Early works [Li et al. 2018; Loubet et al.
2019; Zhang et al. 2019; Bangaru et al. 2020; Zhang et al. 2020, 2021]
derive the exact mathematical form of the derivatives in various
contexts. Following works [Xu et al. 2023; Yan et al. 2022; Yu et al.
2022, 2023; Zhang et al. 2023; Vicini et al. 2021; Nimier-David et al.
2020] attempt to estimate the derivatives more efficiently. An common assumption shared by these works is the underlying surface
representation being mesh or parametric surface.
Differentiable rendering of implicit surfaces. Bangaru et al. [2020]’s
method despite being designed for meshes, can be easily applied
to implicit surfaces. Bangaru et al. [2022] and Vicini et al. [2022]
propose more efficient methods to specifically do differentiable rendering on signed distance functions. Our method applies to generalimplicit surfaces defined as the zero-level set of some function.

3

PRELIMINARIES

We now briefly revisit some mathematical and algorithmic preliminaries related to the differential path integral formulation [Zhang
et al. 2020]. Table 1 presents a list of symbols that are commonlyused in this paper and their definitions.
Path integral. In forward rendering, the path integral formulation
introduced by Veach [1997] has allowed the development of many
advanced Monte Carlo techniques (e.g., bidirectional path tracing).
Under this formulation, the response of a radiometric detector isexpressed as a path integral:∫𝐼=𝑓 ( 𝒙¯ ) d𝜇 ( 𝒙¯ ),(1)𝛀Path-Space Differentiable Rendering of Implicit Surfaces

SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA

where 𝛀 := ∪𝑁 ≥1 M 𝑁 +1 is the path space comprising light transport
paths 𝒙¯ = (𝒙 0, 𝒙 1, . . . , 𝒙 𝑁 ) with M being the union of all object
surfaces, 𝑓 is the measurement contribution function, and 𝜇 is thecorresponding area-product measure.
Material-form reparameterization. Differentiating Eq. (1) withrespect to some arbitrary parameter 𝜋 ∈ R can be complicatedwhen the scene geometry M—and hence the path space 𝛀—evolves
with this parameter. To address this problem, Zhang et al. [2020]
have proposed to reparameterize the path integral (1) using somefixed reference surface B coupled with a differentiable one-to-one
mapping X(·, 𝜋) that transforms the fixed reference surface B to
the evolving surface M (𝜋) for all 𝜋.
This mapping allows a change of variable from light paths 𝒙¯ =(𝒙 0, . . . , 𝒙 𝑁 ) to material paths 𝒑¯ = (𝒑 0, . . . , 𝒑 𝑁 ) with 𝒙 𝑛 = X(𝒑𝑛 , 𝜋)
for all 0 ≤ 𝑛 ≤ 𝑁 . Applying this change of variable to Eq. (1) yields
the material-form path integral:∫¯ d𝜇 ( 𝒑),
¯
𝐼=
𝑓ˆ(𝒑)
(2)
ˆ
𝛀

ˆ = ∪𝑁 ≥1 B 𝑁 +1 . Further,which is over the material path space 𝛀
the integrand of this path integral 𝑓ˆ is the material measurement
contribution defined as𝑁Ö¯ = 𝑓 ( 𝒙¯ )𝑓ˆ( 𝒑)
𝐽 (𝒑𝑛 , 𝜋),
(3)
𝑛=0

where
𝐽 (𝒑, 𝜋) := ∥ d𝐴(𝒙 )/d𝐴(𝒑) ∥(4)is the Jacobian resulting from the reparameterization with 𝐴 indicating the surface-area measure.
Differential path integral. Zhang et al. [2020] have shown thatthe derivative of Eq. (2) with respect to an arbitrary parameter 𝜋
can be expressed as material-form differential path integrals of the
form
∫
∫
d ˆ
d𝐼
¯ d𝜇 ( 𝒑)¯ +¯ 𝑉𝜕 (𝒑 𝐾 ) d𝜇¤ ,=𝑓 ( 𝒑)𝑓ˆ( 𝒑)
(5)
ˆ d𝜋
ˆ
d𝜋
𝛀
𝜕𝛀
|
{z
} |
{z
}
interior

boundary

where the integrand of the interior term is the derivative of thematerial measurement contribution 𝑓ˆ defined in Eq. (3).
The boundary integral in Eq. (5) is over the material boundary
ˆ comprising material boundary paths
path space 𝜕 𝛀𝒑¯ 𝐾 = (𝒑 0, . . . , 𝒑 𝑁 ),𝒑 𝐾 ∈ ΔB (𝒑 𝐾 −1 ),(6)containing exactly one vertex 𝒑 𝐾 constrained to a set of visibility
boundaries ΔB (𝒑 𝐾 −1 ). Precisely, 𝒑 𝐾 is a jump discontinuity pointof the mutual visibility V(X(𝒑 𝐾 −1, 𝜋) ↔ X(·, 𝜋)) with 𝒑 𝐾 −1 fixed.
We define boundary segments to be the pair
(𝒑 𝐾 −1, 𝒑 𝐾 ) ∈ B × ΔB (𝒑 𝐾 −1 ),(7)and refer to them as 𝒑 𝐾 −1 𝒑 𝐾 .
Additionally, 𝑉𝜕 captures the rate at which the visibility boundary shifts (with respect to 𝜋) and is defined asd𝒑𝑉𝜕 (𝒑 𝐾 ) = 𝒏 𝜕 (𝒑 𝐾 ) · 𝐾 ,(8)d𝜋where 𝒏 𝜕 (𝒑 𝐾 ) denotes the unit normal of the visibility boundary(that is perpendicular to both the surface normal 𝒏 B (𝒑 𝐾 ) and the

𝒑$%
𝒑"%

𝒑$#

𝒑!

𝒑#
"

𝝎!

Figure 2: Boundary segments: A boundary light path (6) iscomprised of a boundary segment (7) 𝒑 S0 𝒑 D
0 (shown in red) as
well as a source and a detector subpath (illustrated in orange
and blue, respectively). The two endpoints of the boundary
segment (7) are located on the visibility boundary of eachother. A boundary segment (6) 𝒑 S0 𝒑 D
0 can be uniquely determined with a surface point 𝒑 B ∈ B interior to the segment,along with a direction 𝝎 B .

tangent of the curve). We assume without loss of generality that𝒏 𝜕 (𝒑 𝐾 ) points toward the occluded side of the visibility boundary.
Multi-directional form of the boundary path integral. Numericallyestimating the boundary integral in Eq. (5) is known to be challenging due to the need of sampling boundary segments which requires
silhouette detection.
To avoid silhouette detection at every vertex, Zhang et al. [2020]
propose to re-index a material boundary path asD
𝒑¯ = (𝒑𝑠S, . . . , 𝒑 S0 , 𝒑 D
0 , . . . , 𝒑𝑡 )(9)with 𝒑 S0 𝒑 D
0 being the boundary segment that further separates thefull path 𝒑¯ into a source subpath (𝒑𝑠S, . . . , 𝒑 S0 ) and a detector subpathD(𝒑 D
0 , . . . , 𝒑𝑡 ). Then, the boundary component of Eq. (5) can be
rewritten in a multi-directional form as
∫

∫

∫ ∫
S
𝑓ˆS d𝜇 ( 𝒑¯ S ) 𝑓ˆB𝑓ˆD d𝜇 (𝒑¯ D ) dℓ (𝒑 D
0 ) d𝐴(𝒑 0 ),B

ΔB (𝒑 𝐾 −1 )

ˆ
𝛀

ˆ
𝛀

(10)
¯ 𝑉𝜕 (𝒑 D
where 𝑓ˆB , 𝑓ˆS , and 𝑓ˆD are components of the integrand 𝑓ˆ(𝒑)0)capturing the contributions of the boundary segment, the source
subpath, and the detector subpath, respectively.
Using the multi-directional-form integral (10), one can construct
a boundary path by (i) sampling the material boundary segment
𝒑 S0 𝒑 D
0 ; and (ii) using standard methods such as unidirectional andbidirectional path sampling to build the source and detector subpaths. This avoids the need of expensive silhouette detection [Liet al. 2018] and, thus, scales significantly better to scenes withdetailed geometries.
Previously, Zhang et al. [2020] have derived the exact form ofEq. (10) and the corresponding sampling procedures for explicit
geometries (i.e., polygonal meshes). In this paper, we will derive themulti-directional-form boundary path integral for implicit geometries and develop new Monte Carlo sampling methods to efficiently
estimate this integral.

SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA

4Zhou, Chang, Mukai, Santo, Okura, Matsushita, and Zhao
X(𝒑, 𝜋)

DIFFERENTIAL PATH INTEGRALS FORIMPLICIT SURFACES

In what follows, we explain how the formulation of differential
path integrals can be realized for implicit surfaces.
Specifically, we first discuss how the material-form parameterization can be defined for implicit surfaces in §4.1. Then, we derive
the multi-directional form of the boundary path integral—which is
crucial for efficient estimation of this term—in §4.2.
To simplify the notations, we assume without loss of generalitythat all derivatives are to be evaluated at 𝜋 = 0 and define dℎ𝜕𝜋 ℎ :=(11)d𝜋 𝜋 =0
for all ℎ.

4.1

Material-Form Parameterization

In this paper, we focus on differentiable rendering of evolving closed
implicit surfaces defined as zero-level sets of some scalar-valued
function 𝜙 (𝒙; 𝜋):
M (𝜋) := {𝒙 ∈ R3 : 𝜙 (𝒙; 𝜋) = 0}.

(12)

To facilitate the differentiation of the path integral, we follow
Zhang et al. [2020] by reparameterizing the evolving implicit surface M (𝜋) using a closed reference surface B that is independentof the parameter 𝜋 and defined as
B := M (0) = {𝒙 ∈ R3 : 𝜙 (𝒙; 0) = 0}.

(13)

Unfortunately, unlike when using explicit geometries, it is hard—
if not impossible—to define a one-to-one mapping X(·, 𝜋) that transforms the reference surface B to the evolving M (𝜋) for all 𝜋. Instead, we define a mapping illustrated in Figure 3 that is only oneto-one at 𝜋 = 0 via
X(𝒑, 𝜋) := 𝒑 + 𝑡 (𝒑, 𝜋) 𝒏 B (𝒑),(14)

where:• 𝒏 B (𝒑) denotes the unit normal of the reference surface B at 𝒑
and is assumed to be smooth and differentiable;• 𝑡 (𝒑, 𝜋) ∈ R indicates the distance from 𝒑 to the first intersection
between the (fixed) line passing through 𝒑 via the direction𝒏 B (𝒑) and the evolving surface M (𝜋).
X(·, 0) is the identity map and thus one-to-one.
We define the velocity of 𝒑 ∈ B as the derivative of X w.r.t. 𝜋:𝒗 (𝒑) := (𝜕𝜋 𝑡)(𝒑) 𝒏 B (𝒑),(15)

and with our definition of B and 𝜋 = 0 [Stam and Schmidt 2011] [Yariv
et al. 2020],
𝜕𝜋 𝜙 (𝒑).
(16)(𝜕𝜋 𝑡)(𝒑) = −∥∇𝒑 𝜙 (𝒑)∥Recall that 𝐽 (𝒑, 𝜋) := ∥ d𝐴(𝒙 )/d𝐴(𝒑) ∥ is the ratio between surface
elements at 𝒙 = X(𝒑, 𝜋) ∈ M (𝜋) and 𝒑 ∈ B. It holds thatd(𝜕𝜋 𝐽 ) (𝒑) :=𝐽 (𝒑, 𝜋)= 𝜅 (𝒑) 𝑉 (𝒑).
(17)
d𝜋
𝜋 =0
We leave the proof in the supplemental.

𝒑
ℬ
ℳ(𝜋)

Figure 3: X(𝒑, 𝜋) is the intersection of 𝒑’s normal with the
evolving surface M (𝜋).

4.2

Deriving Boundary Segment

Sampling the boundary segment as in Eq. (10) still requires silhouette detection. Zhang et al. [2020] performs a change of integration
variable to avoid the silhouette detection with meshes for most typesof boundary segment except for those with camera as one of its
endpoints. We derive the change of integration variable for higherorder smooth surfaces. Note Zhang et al. [2023] also derive the samechange of integration variable in a concurrent work. We provide ahigh-level overview here and leave the complete derivation in the
supplemental.
For brevity, we name the integrand of Eq. (10) 𝐹 and drop thesubscripts and start from there:∫ ∫𝐹 (𝒑 S, 𝒑 D ) dℓ (𝒑 D ) d𝐴(𝒑 S ).
(18)
B

ΔB

Our goal is to change the integration variable from (𝒑 S, 𝒑 D ) to(𝒑 B, 𝝎 B ) illustrated in Figure 4.
𝒑 D are silhouette projected onto surface, and we first project
them back onto the silhouette 𝒑 B . The projection results in thefollowing change of variable:∫ ∫dℓ (𝒑 D )𝐹 (𝒑 S, 𝒑 D )dℓ (𝒑 B ) d𝐴(𝒑 S ),(19)dℓ (𝒑 B )B ΣB
where ΣB denotes the silhouette of B as seen from the perspectiveof 𝒑 S .
𝒑 S can vary in two directions, namely dℓ𝑛 and dℓ𝑡 , where dℓ𝑛
is the projection of 𝝎 B onto the tangent plane at 𝒑 S and dℓ𝑡 is the
direction orthogonal to dℓ𝑛 . When 𝒑 S varies along dℓ𝑛 , 𝒑 B varies
along 𝝎 B , so we can map the variation of 𝒑 S along dℓ𝑛 onto thatof 𝒑 B along 𝝎 B . The mapping results in the following change ofvariable:∫ ∫dℓ𝑛 (𝒑 S )dℓ (𝒑 D )𝐹 (𝒑 S, 𝒑 D )dℓ𝑡 (𝒑 S ) d𝐴(𝒑 B ),(20)dℓ𝑛 (𝒑 B )dℓ (𝒑 B )B ΓB
where ΓB denotes the intersection of B and the tangent plane at
𝒑B.
𝒑 S is the intersection of 𝒑 B ’s tangent vector 𝝎 B with B, meaning
there is a one-to-one mapping between 𝒑 S and 𝝎 B defined by theintersection. we use this mapping to map 𝒑 S ’s other variation onto𝝎 B ’s. And the mapping results in the following change of variable:∫ ∫dℓ𝑡 (𝒑 S )dℓ𝑛 (𝒑 S )dℓ (𝒑 D )𝐹 (𝒑 S, 𝒑 D )d𝜃 (𝝎 B ) d𝐴(𝒑 B ). (21)d𝜃 (𝝎 B )dℓ𝑛 (𝒑 B )dℓ (𝒑 B )B S
We leave the detailed derivation in the supplemental.

Path-Space Differentiable Rendering of Implicit Surfaces

SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA𝒑"

two vectors equals zero, illustrated in Figure 5-(a). That dot productnaturally defines another implicit surface we call the orientation
function associated with the surface and camera, defined as follows:𝜓𝜙,𝒄 (𝒑) := ⟨∇𝜙 (𝒑), v𝒄 (𝒑)⟩,

𝝎!

𝒑#
dℓ%

𝒑

dℓ$

Figure 4: We change the integration variable from (𝒑 S, 𝒑 D )to (𝒑 B, 𝝎 B ).

5

OUR MONTE CARLO ESTIMATORSWe use Monte Carlo integration to estimate the differential integrals of Eq. (5). In practice, the interior integral in Eq. (5) can beeasily computed through the differentiation of the forward rendering process. However, acquiring the boundary integrals is notstraightforward since the domain of those boundary integrals is
hard to sample efficiently. To address this problem with implicit surfaces, we introduce new Monte Carlo estimators for the boundary
integrals.
There are two kinds of boundary integrals. Aforementioned
Eq. (21) is called the secondary boundary integral whose boundary
segment has surface points as its two endpoints. On the other hand,
Eq. (10) with 𝒑 S0 being a camera is called the primary boundary
integral whose boundary segment has the camera as one of its
endpoints. The change of integration variables in Section 4.2 helps
to avoid silhouette detection, i.e., finding and sampling points fromthe silhouette of the surface from the point of view of some point,
when sampling the boundary segment. However, when a camera is
one of the endpoints, the change of variable no longer holds. And
we have to perform silhouette detection for the camera to estimate
the primary boundary integral.
We estimate the primary boundary integral by sampling points
directly from the silhouette called the primary edges of the implicit
surface. Estimating the primary boundary integral and sampling
primary edges is easy with explicit geometries such as mesh. Withmesh, one only needs to project the edges onto the image planeand sample from them. However, with implicit surfaces, we don’t
have a finite set of edges at our disposal or a parameterizationof the primary edges. We instead use range analysis [Stolfi and
De Figueiredo 1997] to obtain a tight bound and then utilize a proxyedge.
We estimate the secondary boundary integral by first samplinga surface point and then a tangent direction. We do not have aparameterization of the implicit surface, so we instead use a proxysurface bounding the surface. As for the tangent direction, we use
next-event estimation to importance sample it.

5.1

(22)

!

Sampling Primary Edges

Discontinuities in the camera importance function form the primaryedges of the surface. Primary edges is where the surface normal isperpendicular to the view vector, i.e., where the dot product of thesewhere 𝒄 is the camera position and v := 𝒄 − 𝒑 is the view vector.
We define the primary edges S to be the intersection of 𝜙 and 𝜓 ,visualized in Figure 5-(a), as follows:
S := {𝒑 | 𝜙 (𝒑) = 0,𝜓𝜙,𝒄 (𝒑) = 0}.

(23)

We sample points from S by marching rays over some proxysurface as in Figure 5-(b). Let us use the surface’s bounding box as an
example. The projection of S on one of the faces of the boundingbox is some curve on the face. Let us pretend this projection is
one-to-one for now. Then sampling S is equivalent to sampling the
projection. However, sampling the projection is still not tractable. So
let us simplify the problem once more. Let us project the projection
again onto some proxy curve we know how to sample. In thisexample, that is one of the edges of the face. Then we can sample apoint on S by first sampling a point on the edge and then mapping
this point onto S, illustrated in Figure 5-(b).
Specifically, after sampling a point on the edge, we march along
the direction perpendicular to the edge starting from this point.
At each step, we shoot a ray in the direction normal to the face
toward the surface. We evaluate 𝜓 at the intersection and compare
the value with that from the last step. If the signs are different, we
take this intersection as the sample point on S. This procedure is
written in pseudocode in Algorithm 1.
In general, the ray might intersect either surface multiple times.
That means we have to keep track of all the intersections from theprevious step and keep marching all the way to the other edge. We
instead use range analysis [Stolfi and De Figueiredo 1997] to boundS tighter and assume there is only one intersection in this tighter
bound.
Range analysis [Stolfi and De Figueiredo 1997] is a method topropagate bounds through functions. We specifically use affinearithmetic. Given an axis-aligned bounding box (AABB), affinearithmetic tells us a bound of the function in this AABB. If theoutput bound does not contain zero, then this AABB is guaranteed
not to intersect with the surface.
We cut the local space of the implicit surface into AABBs of thesame size. The exact size is a hyperparameter and we use 6.1e −5 ,i.e., cutting 42 times a unit cube in half, in all our experiments. Afterthe cutting, we apply affine arithmetic to each AABB against 𝜙 and𝜓 . We throw away those AABBs guaranteed not to intersect withboth surfaces and keep the others that possibly intersect both of
the surfaces. Those AABBs bound S as in Figure 6.

5.2

Sampling Secondary Edges

We sample the secondary boundary integral by first sampling asurface point and then a tangent direction.
5.2.1 Sampling Implicit Surfaces. The difficulty of sampling implicit surfaces is our surface representation is implicit, meaning we
do not have a global parameterization. Our solution is to use some
simple proxy surface amenable to parameterization.

SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA

Zhou, Chang, Mukai, Santo, Okura, Matsushita, and Zhao

𝒮 = {𝒑𝟏 , 𝒑𝟐 }

𝒑!

ℰ

𝒮
𝜙

𝜓

ALGORITHM 1: Our primary-edge sampling routine as in Figure 5(b)

ℱ

𝜓

𝜙

𝒒#
𝒑

𝒑"

(a)

𝒒!

(b)

Figure 5: (a) Primary edges is the intersection of 𝜙 and 𝜓 . In
2D, 𝜙 and 𝜓 are curves, and S is a finite set of points. (b) We
sample the primary edges by marching over a proxy surface

𝒮

𝜓
𝜙

(a)

Sample_Primary_Edge()
Input: Face F of bounding box and edge E of F, and max numberof steps 𝑖 max and step size 𝜖
Output: 𝒑 ∈ S
2 begin/* Sample a point from the proxy*/3
Sample a point 𝒒 ∈ E;/* 𝜓 prev stores 𝜓 ’s value at the previous step andis initialized to zero*/4
𝜓 prev ← 0;/* March along the normal of E over F*/
5
for 𝑖 = 1, 2, . . . , 𝑖 max do/* Intersect the ray starting from 𝒒 along thenormal 𝒏 F of face F toward implicit surface
𝜙*/6𝒑 ← ray_intersect(𝒒, 𝒏 F );7
if 𝒑 == NOT_HIT then8𝜓 prev ← 0;9else/* Evaluate 𝜓 at the current step*/10
𝜓𝑖 ← 𝜓 (𝒑);/* If sign of 𝜓 at this step is different
from that at the last step*/11
if 𝜓 prev · 𝜓𝑖 < 0 then/* We just intersected 𝜓 , meaning 𝒑
intersects both 𝜙 and 𝜓*/
12
return 𝒑;
13
end14
𝜓 prev ← 𝜓𝑖 ;15
end/* Step ahead along the normal 𝒏 E of edge E */16𝒒 ← 𝒒 + 𝜖 · 𝒏E
17
end
18
return NULL
19 end
1

𝒒!"#

(b)

Figure 6: AABBs classified by affine arithmetic as possiblyintersect with (a) 𝜙 or 𝜓 and (b) both of them. We boundsilhouette S by (b) AABBs that possibly intersect with boththe implicit surface 𝜙 and its associated orientation function
𝜓.

Imagine our implicit surface being a sphere. We assume there isa bounding box surrounding the sphere. Then we can use one of
the faces of the box as the proxy.
Let F be one of the faces of the bounding box of the sphere
B. We consider the following integral defined over B and try tochange the domain to F .
∫ℎ(𝒑) d𝒑(24)
B

We next define the transformations needed. For any point 𝒒 ∈ F ,we can cast a ray 𝑟 := {𝒒 + 𝑡 𝒏 F (𝒒) | 𝑡 ∈ R+ } from 𝒒 toward B.
This ray will intersect with B two times. Let T1 : F → B be thetransformation mapping 𝒒 ∈ F to the first intersection and T2
the second. T1 and T2 cut B into two disjoint parts and so do theintegral:∫∫ℎ(𝒑) d𝒑 +T1 ( F)ℎ(𝒑) d𝒑(25)

T2 ( F)

Next we change the domains to F :∫∫ℎ(T1 (𝒒)) 𝐽 T1 (𝒒) d𝒒 +ℎ(T2 (𝒒))𝐽 T2 (𝒒) d𝒒
F

=

F

∫ ∑︁
2
F 𝑖=1

(26)
ℎ(T𝑖 (𝒒))𝐽 T𝑖 (𝒒) d𝒒

where 𝐽 T𝑖 are the Jacobian determinants of T𝑖 .
In general, B could be any surface, meaning there could be morethan two intersections, and each point in F could have a different
number of intersections. For the former problem, we simply continue
the summation on. The latter means T𝑖 might not be defined forsome points in F since the 𝑖-th intersection might not exist. We let
those points have zero contribution as follows:ℎ(T𝑖 (𝒒))𝐽 T𝑖 (𝒒), T𝑖 (𝒒) defined𝑔𝑖 (𝒒) :=(27)
0,otherwise
Letting 𝑁 be the max number of intersections, we can change thedomain of Eq. (24) from the implicit surface B to the face F of thebounding box as follows:𝑁 ∫∑︁𝑔𝑖 (𝒒) d𝒒.
(28)
𝑖=1

F

Then we can estimate Eq. (28) by directly sampling points from its
domain F .
5.2.2 Next Event Estimation. After sampling a surface point, we
could simply uniformly sample a tangent direction. But we cando it more efficiently using next-event estimation. We do it by intersecting the tangent plane with light and then sampling points
from the intersection as in Figure 7-(b). For environment maps, wePath-Space Differentiable Rendering of Implicit Surfaces

ℱ

ℒ

T𝒑 ℬ

𝒯𝟏 (𝒒)

Ordinary

FD

Ours

(b)

pre-compute some plane-sphere intersections with uniformly distributed normals determining the plane orientation. We then build
a 1D distribution conditioned on the normal direction off of thoseintersections. This way, we do not have to build distributions on the
fly. We assume the number of lights is one in our experiments andleave efficient plane intersection with a large number of polygons
to future work.

Figure 8: We validate our method by comparing derivative
images estimated by our method to that by finite difference(FD) with a large number of samples. Our NEFERTITI andSUZANNE results are computed using our unidirectionaland bidirectional estimators, respectively. Our results closely
match the FD reference and demonstrate the correctness of
our implementation.

EXPERIMENTS

To assess the effectiveness of our method, we conduct experimentson several test scenes.

Settings

We implement our Monte Carlo estimator with both unidirectionaland bidirectional path sampling based on the cuda_ad_rgb backendof Mitsuba 3 [Jakob et al. 2022]. We conduct all our experiments on
an NVIDIA RTX 4090 GPU.
We represent an implicit surface as a grid-based SDF with tricubic interpolation as in diffsdf [Vicini et al. 2022] so that thesurface normal is continuous.
Optimization setup. We use the Adam optimizer [Kingma and Ba
2014] and the L2 loss in all our inverse-rendering experiments. The
number of views is eight, and the batch size is two. The runtime ofone optimization is typically 30 to 90 minutes which usually uses
64 to 256 epochs.

6.2

Ours

ℬ

Figure 7: (a) We sample points from the implicit surface 𝜙by sampling a point 𝒒 from the proxy F and project 𝒒 onto
the implicit surface 𝜙. We then take all the intersections,
T1 (𝒒) and T2 (𝒒) in this example, as the sample points fromthe implicit surface 𝜙. (b) We intersect the tangent plane T𝒑 B
at 𝒑 with light L and only sample from tangent directions
intersecting the light L

6.1

FD

𝒒

(a)

6

Ordinary𝒑

𝜙
𝒯𝟐 (𝒒)SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA

Differentiable Rendering

We validate our estimators in Figure 8 by comparing derivative
estimated by our estimators and that by finite difference using alarge number of samples. In the top row, we differentiate the image
w.r.t. the horizontal translation of the shape. In the bottom row, we
differentiate the image w.r.t. the vertical-axis rotation of the shape.
The area light near the ceiling is flipped upside down, and the sceneis thus dominated by indirect illumination.
In Figure 9, we show that our bidirectional estimator performs
better than our unidirectional one with hard-to-sample lights andspecular materials.

6.3

Inverse Rendering

We quantitatively compare our optimized shape with diffsdf’s using
Chamfer distance (between target meshes and those extracted from
optimized implicit surfaces).
In Figure 12, the light is blocked by the ceiling with only a small
hole letting light in. And We only observe the shadow of the shape,so the gradient only consists of the secondary boundary integral.
Since light is occluded by the ceiling except for a small region, it isvery hard to sample light in this configuration. Our guiding grid
helps us to find the un-occluded region of the light and makes ourderivative estimates less noisy.
In Figure 13, the scene is lit by a small light flipped upside down.
The light is hard to sample in this scene and so our bidirectional
estimator is more efficient than the other unidirectional estimators.
Our bidirectional estimator’s derivative image is cleaner, and the
optimized shape is closer to the target and smoother than the otherestimators.
In Figure 10 and Figure 11, our estimator computes cleaner gradient for all the three components, leading to much better results
than diffsdf.

7

DISCUSSION AND CONCLUSION

Limitations. As is the case with all PSDR methods, samplingboundary paths efficiently is nontrivial and typically requires someform of guiding. Our current implementation uses a simple regulargrid-based scheme (similar to [Zhang et al. 2020]) which is known
to have difficulties handling reflections off of highly glossy surfaces. Fortunately, we expect using concurrent projective sampling [Zhang et al. 2023] to greatly mitigate this problem.

SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA
OrdinaryFD

Zhou, Chang, Mukai, Santo, Okura, Matsushita, and Zhao

are known to be essential for efficient estimation of these integrals.
Based on our generalized formulation, we introduced new Monte
Carlo estimators—including bidirectional ones—capable of rendering implicit surfaces with complex light transport effects.

ACKNOWLEDGMENTS
Ours (unidir.)Ours (bidir.)We thank the anonymous reviewers for their constructive suggestions. We also thank Heng Guo for his advice and insights. This
work was partially supported by JSPS KAKENHI grant JP23H05491
and NSF grant 1900927.

REFERENCES

Figure 9: Differentiable rendering comparison (SUZANNE)
between our bidirectional and unidirectional estimator.
Scenes where lights are difficult to sample and with specularmaterials need to be handled by bidirectional estimators.

Future work. How to efficiently guide the new secondary boundary integral for implicit surfaces is something one can explore. One
thing that is hard about it is that for one sample on the proxy shape,
there are multiple corresponding surface points. We only use a single guiding grid, and that means the tangent direction is fixed overthose possibly multiple surface points. We only did experiments
on a grid-based SDF while neural SDF is a more popular representation. Theoretically, our theory supports arbitrary continuous
level-set functions —including neural SDFs. In practice, we build
our system based on diffsdf’s implementation (a custom branch of
Mitsuba) that uses grid-based SDFs. Integrating simple neural SDFs
(e.g., those using small MLPs) into the system should be possible,
but supporting more complex neural representations using Dr.Jit
efficiently can be challenging (at the system level).
Bias. Theoretically, our formulation is unbiased. In practice,small amounts of bias can be introduced by the ray-surface intersection and the sampling of silhouette points (needed by the
primary boundary paths) as follows. Computing the intersection
between a ray and an implicit surface (using, for example, spherical
tracing) is normally approximated. In other words, the computed
intersection is normally very close to the surface but may not be
exactly on it. We note that existing methods like diffsdf also havethis type of bias. When finding the silhouette points, our method
assumes at most one intersection for each sample point on theproxy curve—which may be violated when the AABBs are notsmall enough. Fortunately, we found the bias to have little impact
on the reconstruction quality in all our experiments.
Conclusion. In this paper, we generalized path-space differentiable rendering to support implicit surfaces. Specifically, we demonstrated how material-form parameterization can be realized without
requiring a global parameterization of the reference surface. Also,
we derived multi-directional form of boundary path integrals which

Sai Bangaru, Michael Gharbi, Tzu-Mao Li, Fujun Luan, Kalyan Sunkavalli, Milos
Hasan, Sai Bi, Zexiang Xu, Gilbert Bernstein, and Fredo Durand. 2022. Differentiable Rendering of Neural SDFs through Reparameterization. In ACM SIGGRAPH
Asia 2022 Conference Proceedings (Daegu, Republic of Korea) (SIGGRAPH Asia ’22).
Association for Computing Machinery, New York, NY, USA, Article 22, 9 pages.
https://doi.org/10.1145/3550469.3555397
Sai Praveen Bangaru, Tzu-Mao Li, and Frédo Durand. 2020. Unbiased Warped-Area
Sampling for Differentiable Rendering. ACM Trans. Graph. 39, 6 (2020), 245:1–
245:18.
Alan H Barr. 1986. Ray tracing deformed surfaces. ACM SIGGRAPH Computer Graphics
20, 4 (1986), 287–296.
Robert A Drebin, Loren Carpenter, and Pat Hanrahan. 1988. Volume rendering. ACM
Siggraph Computer Graphics 22, 4 (1988), 65–74.
Jorge Flórez, Mateu Sbert, Miguel A Sainz, and Josep Vehí. 2008. Efficient ray tracing
using interval analysis. In Parallel Processing and Applied Mathematics: 7th International Conference, PPAM 2007, Gdansk, Poland, September 9-12, 2007 Revised Selected
Papers 7. Springer, 1351–1360.
John C Hart. 1996. Sphere tracing: A geometric method for the antialiased ray tracingof implicit surfaces. The Visual Computer 12, 10 (1996), 527–545.
Wenzel Jakob, Sébastien Speierer, Nicolas Roussel, Merlin Nimier-David, Delio Vicini,
Tizian Zeltner, Baptiste Nicolet, Miguel Crespo, Vincent Leroy, and Ziyi Zhang.
2022. Mitsuba 3 renderer. https://mitsuba-renderer.org.
Diederik P Kingma and Jimmy Ba. 2014. Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980 (2014).
Marc Levoy. 1988. Display of surfaces from volume data. IEEE Computer graphics andApplications 8, 3 (1988), 29–37.
Tzu-Mao Li, Miika Aittala, Frédo Durand, and Jaakko Lehtinen. 2018. Differentiable
Monte Carlo ray tracing through edge sampling. ACM Trans. Graph. 37, 6 (2018),
222:1–222:11.
William E Lorensen and Harvey E Cline. 1998. Marching cubes: A high resolution 3D
surface construction algorithm. In Seminal graphics: pioneering efforts that shaped
the field. 347–353.
Guillaume Loubet, Nicolas Holzschuch, and Wenzel Jakob. 2019. Reparameterizing
discontinuous integrands for differentiable rendering. ACM Trans. Graph. 38, 6(2019), 228:1–228:14.
Merlin Nimier-David, Sébastien Speierer, Benoît Ruiz, and Wenzel Jakob. 2020. Radiative backpropagation: an adjoint method for lightning-fast differentiable rendering.
ACM Trans. Graph. 39, 4 (2020), 146:1–146:15.
Jos Stam and Ryan Schmidt. 2011. On the Velocity of an Implicit Surface. ACM
Trans. Graph. 30, 3, Article 21 (may 2011), 7 pages. https://doi.org/10.1145/1966394.
1966400
Jorge Stolfi and Luiz Henrique De Figueiredo. 1997. Self-validated numerical methods
and applications. In Brazilian mathematics colloquium monograph. IMPA ˆ eRio deJaneiro Rio de Janeiro, 1–116.
Eric Veach. 1997. Robust Monte Carlo methods for light transport simulation. Vol. 1610.
Stanford University PhD thesis.
Delio Vicini, Sébastien Speierer, and Wenzel Jakob. 2021. Path Replay Backpropagation:
Differentiating Light Paths Using Constant Memory and Linear Time. ACM Trans.
Graph. 40, 4, Article 108 (2021), 108:1–108:14 pages.
Delio Vicini, Sébastien Speierer, and Wenzel Jakob. 2022. Differentiable Signed Distance
Function Rendering. ACM Trans. Graph. 41, 4 (2022), 125:1–125:18.
Andrew P Witkin and Paul S Heckbert. 1994. Using particles to sample and control
implicit surfaces. In Proceedings of the 21st annual conference on Computer graphics
and interactive techniques. 269–277.
Peiyu Xu, Sai Bangaru, Tzu-Mao Li, and Shuang Zhao. 2023. Warped-Area Reparameterization of Differential Path Integrals. ACM Trans. Graph. 42, 6 (2023),
213:1–213:18.
Kai Yan, Christoph Lassner, Brian Budge, Zhao Dong, and Shuang Zhao. 2022. Efficientestimation of boundary integrals for path-space differentiable rendering. ACM
Trans. Graph. 41, 4 (2022), 123:1–123:13.

Path-Space Differentiable Rendering of Implicit Surfaces

Lior Yariv, Yoni Kasten, Dror Moran, Meirav Galun, Matan Atzmon, Basri Ronen, andYaron Lipman. 2020. Multiview Neural Surface Reconstruction by Disentangling
Geometry and Appearance. Advances in Neural Information Processing Systems 33(2020).
Z. Yu, C. Zhang, O. Maury, C. Hery, Z. Dong, and S. Zhao. 2023. Efficient Path-Space
Differentiable Volume Rendering With Respect To Shapes. Computer Graphics
Forum 42, 4 (2023).
Zihan Yu, Cheng Zhang, Derek Nowrouzezahrai, Zhao Dong, and Shuang Zhao. 2022.
Efficient Differentiation of Pixel Reconstruction Filters for Path-Space Differentiable
Rendering. ACM Trans. Graph. 41, 6 (2022), 191:1–191:16.
Cheng Zhang, Bailey Miller, Kai Yan, Ioannis Gkioulekas, and Shuang Zhao. 2020.
Path-space differentiable rendering. ACM Trans. Graph. 39, 4 (2020), 143:1–143:19.
Cheng Zhang, Lifan Wu, Changxi Zheng, Ioannis Gkioulekas, Ravi Ramamoorthi, andShuang Zhao. 2019. A differential theory of radiative transfer. ACM Trans. Graph.
38, 6 (2019), 227:1–227:16.
Cheng Zhang, Zihan Yu, and Shuang Zhao. 2021. Path-space differentiable renderingof participating media. ACM Trans. Graph. 40, 4 (2021), 76:1–76:15.
Ziyi Zhang, Nicolas Roussel, and Wenzel Jakob. 2023. Projective Sampling for Differentiable Rendering of Geometry. Transactions on Graphics (Proceedings of SIGGRAPH
Asia) 42, 6 (Dec. 2023). https://doi.org/10.1145/3618385

SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA

SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA

Config.

Initial (Chamfer Dist.)

Zhou, Chang, Mukai, Santo, Okura, Matsushita, and Zhao

Ordinary

FD Reference

Target (0)Ours (7.8e!" )Ours

diffsdf (339.3e!" )diffsdf

Loss

ForwardShape Vis.

Figure 10: Differentiable and inverse rendering comparison (CHAIR) between our unidirecitional estimator and diffsdf [Vicini
et al. 2022]. Our estimator can be used in inverse rendering to recover complex topology and thin features such as the chair inthis scene.
Config.

Initial (Chamfer Dist.)

Ordinary

FD Reference

Target (0)Ours (7.9e!" )Ours

diffsdf (37.0e!" )diffsdf

Loss

ForwardShape Vis.

Figure 11: Differentiable and inverse rendering comparison (BUNNY) between our unidirecitional estimator and diffsdf [Vicini
et al. 2022]. Our estimator estimates all the three components better than diffsdf and leads to the optimized shape being
smoother.

Path-Space Differentiable Rendering of Implicit Surfaces

Config.

Initial (Chamfer Dist.)

SIGGRAPH Conference Papers ’24, July 27–August 01, 2024, Denver, CO, USA

Ordinary

FD Reference

Target (0)Ours (6.24e!" )Ours

diffsdf (6.73e!" )diffsdf

Loss

ForwardShape Vis.

Figure 12: Differentiable and inverse rendering comparison (SHADOW) between our unidirecitional estimator anddiffsdf [Vicini et al. 2022]. We compare the secondary boundary integral of our estimator and diffsdf’s in this scene where we
only observe the shadow of the shape. We achieve less variance in the secondary boundary integral by explicit sampling thiscomponent and guiding.
Ordinary

Initial (Chamfer Dist.)

FD Reference

Target (0)Ours (bidir.)Ours (unidir.)Ours-bidir (2.7e!" )diffsdf

Loss

Ours-unidir (10.2e!" )diffsdf (2000e!" )ForwardShape Vis.

Figure 13: Differentiable and inverse rendering comparison (BOB) between our bidirectional estimator, unidirectional estimatorand diffsdf [Vicini et al. 2022]. With scenes where lights are difficult to sample, our bidirectional estimator can achieve lessvariance than other unidirectional estimators.