A Heat Method for Generalized Signed Distance

NICOLE FENG and KEENAN CRANE, Carnegie Mellon University

We introduce a method for approximating the signed distance function (SDF)
of geometry corrupted by holes, noise, or self-intersections. The method
implicitly de�nes a completed version of the shape, rather than explicitly
repairing the given input. Our starting point is a modi�ed version of the
heat method for geodesic distance, which di�uses normal vectors rather
than a scalar distribution. This formulation provides robustness akin to
generalized winding numbers (GWN), but provides distance function rather
than just an inside/outside classi�cation. Our formulation also o�ers several
features not common to classic distance algorithms, such as the ability to
simultaneously �t multiple level sets, a notion of distance for geometry that
does not topologically bound any region, and the ability to mix and match
signed and unsigned distance. The method can be applied in any dimension
and to any spatial discretization, including triangle meshes, tet meshes,
point clouds, polygonal meshes, voxelized surfaces, and regular grids. We negative 0 positive input curve
evaluate the method on several challenging examples, implementing normal
o�sets and other morphological operations directly on imperfect curve and Fig. 1. The signed heat method computes a well-behaved signed distance
surface data. In many cases we also obtain an inside/outside classi�cation function (SDF), even for imperfect or incomplete input geometry. Here we
dramatically more robust than the one obtained provided by GWN. compute geodesic signed distance to the broken curve shown in magenta.
CCS Concepts: • Computing methodologies ! Shape analysis; • Math-
ematics of computing ! Partial di�erential equations.
Additional Key Words and Phrases: geometry processing, signed distance More precisely, we aim to approximate the signed distance to
the boundary of a shape within a �xed domain "—for example, a
ACM Reference Format: curve on a surface, or a surface in R3 . When " is itself a curved
Nicole Feng and Keenan Crane. 2024. A Heat Method for Generalized Signed
manifold, we consider the signed geodesic distance, i.e., the minimal
Distance. ACM Trans. Graph. 43, 4, Article 92 (July 2024), 20 pages. https:
//[Link]/10.1145/3658220
length among all paths along ". We assume we are given as input
a defective description ⌦ ⇢ " of the shape boundary: a “broken”
1 INTRODUCTION curve or surface which can have holes, self-intersections, and noise.
We then seek a reasonable approximation of the SDF for the true
A signed distance function (SDF) describes the distance to the clos- (unknown) geometry.
est point on the boundary of a given shape, using sign to indicate Our solution to this problem, the signed heat method (SHM), pro-
whether points are inside or outside. SDFs are a basic component of vides a “least squares” approach: we �nd a signed distance approxi-
numerous algorithms from geometric modeling [Museth et al. 2002], mation that matches the input geometry as well as possible, even
physical simulation [Osher et al. 2004], rendering [Quilez 2008], path when it does not perfectly bound any region. There are three basic
planning [Oleynikova et al. 2016], geometric learning [Yariv et al. steps (Figure 2):
2023] and computer vision [Vicini et al. 2022]. In many problems,
geometry exhibits severe defects due to errors in modeling or acqui- I. Di�use the normals of ⌦ for a small time C > 0, by solving a
sition. Unfortunately, traditional SDF algorithms assume a “perfect” vector-valued heat equation.
watertight description of the boundary, relying on a well-de�ned II. Normalize the di�used vectors, yielding a unit vector �eld
inside and outside to determine the sign [Sethian 1999]. At the same that approximates the gradient of the (unknown) SDF.
time, robust inside/outside tests do not provide a useful signed dis- III. Find the function whose gradient best matches the unit vector
tance approximation [Jacobson et al. 2013]. Our method provides �eld, by solving a scalar Poisson equation.
the best of both worlds: a generalized signed distance (GSD) function A di�usion-based approach is valuable because it extends in a robust
computed directly from broken geometry, which provides not only way to imperfect geometry. For instance, if a curve has gaps, or a
an inside/outside classi�cation, but also distance information. surface has holes, di�usion averages together normals at nearby
points—providing smooth interpolation of the observed data. Nor-
Authors’ address: Nicole Feng; Keenan Crane, Carnegie Mellon University, 5000 Forbes
Ave, Pittsburgh, PA, 15213. malization of gradients then ensure that we recover a distance ap-
proximation, rather than just a smooth function.
Permission to make digital or hard copies of part or all of this work for personal or Why not simply “�x” the broken geometry before computing
classroom use is granted without fee provided that copies are not made or distributed
for pro�t or commercial advantage and that copies bear this notice and the full citation distance? For example, one could �t an implicit function, extract
on the �rst page. Copyrights for third-party components of this work must be honored. an explicit mesh, then compute signed distance using any standard
For all other uses, contact the owner/author(s). algorithm. Here, however, errors in reconstruction can propagate
© 2024 Copyright held by the owner/author(s).
0730-0301/2024/7-ART92 into distance computation (Figure 3, left). Moreover, each step (say,
[Link] computing signed distance) already costs as much as our entire

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
92:2 • Feng and Crane

diffuse normals normalize vectors integrate offset surfaces

zero
level
set

regular generalized
input generalized grid signed
signed distance distance
Fig. 2. The three basic steps of the signed heat method.

Fig. 4. Our method can be applied on virtually any spatial data structure, in
any dimension. Here for instance we compute generalized signed distance to
method. We avoid intermediate representations altogether, provid- a badly broken surface (le�), by solving on a regular grid in R3 . Contouring
ing more accurate distance for about an order of magnitude lower this function yields well-behaved and evenly-spaced o�set surfaces (right).
cost (Section 9.5.3).
By formulating our method from the perspective of short-time
2.1 Inside/Outside Classification
heat di�usion, it also inherits many bene�ts of the original heat
method: it can compute (geodesic) distance on curved domains, it Given a corrupted region boundary m , there are two basic strategies
applies directly to most discretizations (triangle/tet meshes, gen- for estimating whether a given point G 2 " is inside or outside:
eral polygon meshes, point clouds, etc.), is robust to noise in both • Pseudonormal test. Let =(G) be the
the input geometry and the underlying domain (since it is based normal at the closest point G 2 m . The sign
on solving “nice” elliptic problems), and is easily accelerated and of the function q? (G) := h=(G), G Gi indicates
parallelized using existing systems for sparse linear algebra (since whether G is inside or outside [Bærentzen
the main cost is solving two standard linear systems). More broadly, and Aanæs 2005]. This test can be applied as-
a variational, PDE-based approach to SDF computation enables ex- is to broken geometry, e�ectively using the
tensions not possible even with other geodesic distance algorithms— tangent plane of the closest boundary point as
such as �nding the “best �t” distance function for a collection of a proxy for missing geometry. However, since
partially-observed level sets (Section 7.1), or computing a notion of tangent planes are not globally consistent, neither is the pseudonor-
generalized signed distance for input shapes that do not topologi- mal distance q? , which in general is not even ⇠ 0 continuous (inset).
cally bound regions (Section 7.2). • Winding number. The winding number captures the number
of times m wraps around the point G, and is hence nonzero only
when G is inside . It is equivalent (up to a constant) to signed
2 RELATED WORK
solid angle—known in graphics as the generalized winding number
Our algorithm complements a large set of existing methods: (GWN) [Jacobson et al. 2013], which is well-de�ned even for broken
• Methods for robust inside/outside classi�cation, which do geometry. GWN supplies robust inside/outside queries in several
not directly provide signed distance (Section 2.1). algorithms [Zhou et al. 2016; Hu et al. 2018], but de�nes a harmonic
• Methods for signed distance computation, which generally function rather than an SDF—and hence cannot be used for a broad
do not work well for broken geometry (Section 2.2). variety of geometric tasks (e.g., Figure 26, top).
input pseudo-
Each method implicitly applies a prior
It also of course builds on earlier heat methods (Section 2.3). Un- normal GWN (linear vs. harmonic extension), suitable
fortunately, simply signing the unsigned distance yields a function
in di�erent scenarios (inset). GSD e�ec-
completely di�erent from the true SDF—see Figure 6. In general,
tively interpolates between these options:
there are only a few methods that explore robust SDF computation—
as C ! 0, points inherit the normal at the
and none suitable for curved geometry.
closest point; as C ! 1, di�used vectors
become componentwise harmonic (Fig-
ure 5). Yet unlike the pseudonormal, GSD
generalized generalized signed distance is based on global integration, making it
winding number
less sensitive to small perturbations. Unlike GWN, which must (by
original de�nition) interpolate the input, GSD is robust to noise in positions
shape (Figure 3), and yields better surface completions due to smooth
extrapolation of normal information (Figure 27).
corrupted (constrained) (unconstrained)
2.1.1 Curved Domains. Inside/outside tests are also not easily gen-
Fig. 3. Both GWN and GSD implicitly define a completed surface—but GSD eralized to curved domains. Here, closest point queries entail com-
also provides distance information. Moreover, GWN must interpolate the puting the geodesic distance at every point—at which point one
input geometry, leading to noisy output (le�). GSD can either interpolate may as well compute signed distance rather than the pseudonormal.
or approximate the input, yielding a more faithful completion (far right). Likewise, the recent method of Feng et al. [2023] generalizes GWN

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
 A Heat Method for Generalized Signed Distance • 92:3

pseudonormal generalized signed distance winding signing unsigned distance our solution
number

Fig. 6. For broken geometry, simply signing unsigned distance does not
Fig. 5. For di�erent di�usion times, our method e�ectively interpolates
work, yielding completely di�erent level sets (le�) from the distance to the
between the “linear” prior used by the pseudonormal test (as C ! 0) and
unbroken curve (right).
the “harmonic” prior used by generalized winding numbers (as C ! 1).

ours fast marching

to curved surfaces, but again provides no notion of distance (and

has a similar cost to GSD).

2.2 Signed Distance Computation

For watertight geometry in R= , unsigned distance can be computed
via fast, exact closest point queries [Sawhney et al. 2020] and signed
via basic inside-outside tests like ray shooting [Haines 1994]. Alter- Fig. 7. Our elliptic formulation regularizes the available information, grace-
natively, one can sample geometry onto a grid and use methods like fully handling broken curves (le�). In contrast, hyperbolic methods like
fast sweeping [Osher et al. 2004]. fast marching consider an expanding wavefront that propagates sign and
For broken geometry, it is tempting to again just sign the unsigned distance errors (right).
distance—but unless gaps are very small, the resulting function will
be quite di�erent from the true SDF (Figure 6). Likewise, wavefront-
based methods like fast marching [Kimmel and Sethian 1998] and 2.2.2 Geodesic Distance. With few exceptions, work over the past
learning-based variants [Lichtenstein et al. 2019; Huberman et al. 40 years on geodesic distance computation has aimed largely at im-
2023] can be executed on broken geometry, but propagation of proving speed and accuracy—Crane et al. [2020] provides a detailed
orientation errors yields signi�cant artifacts (Figure 7). survey. We focus not on this performance race, but rather on a more
Some “signed distance” methods do not actually compute a true robust notion of signed distance that can be applied in a broader
distance function, but rather just a smooth implicit function that range of real-world scenarios. However, robustness does not come
vanishes near the input geometry. Like GWN, such functions are at a major cost: our method remains competitive with widely-used
unreliable for tasks like morphological modeling or sphere trac- methods such as the heat method and fast marching (Section 9.5).
ing [Hart 1996]. For instance, smooth signed distance (SSD) [Calakli Most past work considers distance to isolated points, which is
and Taubin 2011] is actually a biharmonic function with unit-norm always unsigned, though a few methods consider unsigned distance
gradient only at the zero set. Likewise, many so-called neural SDFs to curves embedded in surfaces [Bommes and Kobbelt 2007; Trettner
encourage the signed distance property only near the zero set, or et al. 2021]; as in R= , signing unsigned distance will fail to yield
weakly enforce the eikonal condition |rq | = 1 [Gropp et al. 2020], satisfactory results for broken curves. Overall, we are not aware of
which is insu�cient to characterize SDFs [Marschner et al. 2023]. any work that computes robust signed distance to broken curves
on surfaces, or more generally, submanifolds of a curved domain.
2.2.1 Robust Signed Distance. Some past work considers regular-
ized signed distance for broken geometry. E.g., Xu and Barbič [2014] 2.3 Heat Methods
e�ectively apply morphological operations on unsigned distance The signed heat method (SHM) builds on the unsigned heat method
to “heal” holes of a user-speci�ed size—at the cost of losing small (UHM) of Crane et al. [2013b], and the vector heat method (VHM) of
details. Brunton and Rmaileh [2021] use smoothing to mitigate arti- Sharp et al. [2019c]. UHM di�uses a scalar distribution concentrated
facts arising from pseudonormal-like sign estimation (Section 2.1); on source points for a small time C > 0; Varadhan’s formula then
Bærentzen [2005] likewise smooth a pseudonormal distance ob- implies that, as C ! 0, the gradient of the resulting heat distribu-
tained from voxelization. Mullen et al. [2010] sign the unsigned dis- tion becomes parallel with the gradient of the distance function
tance by performing ray intersection tests on a small band around q [Varadhan 1967]. Since the true distance function satis�es the
the zero set. As noted above, neural networks have also been used eikonal property |rq | = 1, normalizing and integrating the initial
to �t distance-like functions to noisy/incomplete input [Park et al. gradient �eld (via a Poisson equation) recovers an accurate distance
2019; Atzmon and Lipman 2019; Gropp et al. 2020], but do not pro- approximation. VHM instead starts by di�using tangent vectors at
duce true SDFs even for watertight geometry. Unlike GSD, none points; a similar normalization yields the parallel transport of these
of the methods in this section apply to curved domains, nor gen- vectors along minimal (shortest) geodesics.
eralize to alternative spatial discretizations, nor handle the richer Relative to UHM we change only the �rst step, di�using normal
constraints furnished by our variational formulation. vectors concentrated along curves or surfaces, rather than a scalar

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
92:4 • Feng and Crane

measure concentrated at points. Relative to VHM, we change its 3.1 Signed Heat Method
last three steps, integrating (normalized) di�used vectors to obtain The �rst step of our algorithm is to di�use the normals # from ⌦
distance. Our basic insight, not considered in prior work, is that to the rest of the domain " for a short time C > 0, which extends
di�used normals will be parallel to the gradient of signed distance, information about surface orientation to the rest of the domain. In
since parallel transport ensures they are tangent to oriented min- particular we solve a vector-valued di�usion equation
imal geodesics (Section 3.1). Hence, normalizing and integrating
d r- C > 0,
these vectors recovers an accurate SDF approximation. For broken dC -C = C,
(1)
geometry the story remains largely the same, except that di�usion -0 = # `⌦ .
e�ectively interpolates normals from nearby points to obtain well-
On R= this equation simply di�uses each scalar component of the
behaved gradients.
vector �eld; more generally it accounts for how vectors di�use across
Overall, our method inherits many bene�ts of past heat methods:
a domain with curvature (see Sharp et al. [2019c, Section 4.3]). As
for instance, it is not tied to a particular spatial discretization (Sec-
with scalar di�usion, the magnitude of -C decays exponentially
tion 8), and is robust not only to broken source geometry, but also
with distance from ⌦, but its direction remains well-de�ned almost
poor discretization of the domain itself (Section 9.4). Other work im-
everywhere since the support of the vector heat kernel is all of ".
proves or applies heat methods in various ways—for instance, Sharp
Using di�usion to extrapolate orienta-
et al. [2019b] and Gillespie et al. [2021] show how the accuracy and
tion is not a heuristic, but is rather moti-
robustness of the heat method can be signi�cantly improved via use
vated by a key observation from di�eren-
of intrinsic triangulations—a strategy we also employ here. Several
tial geometry: as C ! 0 the di�used vector
works improve e�ciency and scalability via parallel or iterative
-C (G) at each point G 2 " becomes paral-
solvers [Tao et al. 2019; Rawat and Biswas 2022a,b], Belyaev et al.
lel to the normal # (Ḡ) at the closest point
[2013] improve accuracy by iterating the integration step, Sun and
Ḡ 2 ⌦. More precisely, it aligns with the
Liu [2022] extend the method to inhomogeneous and anisotropic
vector obtained via parallel transport of # (Ḡ) along a minimal geo-
metrics, and Litman and Bronstein [2016] approximate all-pairs
desic W [Berline et al. 1992, Theorem 2.30]. Since parallel transport
distance using spectral methods. These extensions o�er rich oppor-
along geodesics preserves tangency, this vector will be tangent to W
tunities for future improvement and generalization of the signed
itself—and since traveling along W is the quickest way back to ⌦, it
heat method.
must be parallel to the unsigned distance gradient rq. Moreover,
since we transport oriented normals, we get the correct sign. We
3 SMOOTH THEORY can hence normalize -C to obtain an approximation .C : -C /k-C k
Fundamentally, our algorithm is de�ned in terms of operations in of the signed distance gradient.
the continuous setting, as described here. It can then be discretized The vector �eld .C will not describe exact gradients for any SDF,
in many di�erent ways, as explored in Sections 5 and 8. due to both the di�usion approximation—and more signi�cantly—
Throughout we consider an =-dimensional errors in the input. However, we can still look for the function q
Riemannian manifold " with metric 6, and whose gradient is as close as possible, in a least-squares sense, to .C .
want to compute the signed distance function In particular, we seek a minimizer for the problem
q for a codimension-1 submanifold ⌦ ⇢ " π
(e.g., curves within a surface, or surfaces min krq .C k 2 . (2)
q: "!R "
within a volume). In general we assume that
this data might represent a corrupted version of ideal input, mean- Using integration by parts, one can show that a minimizer satis�es
ing " and/or ⌦ may have holes, self-intersections, noise, and may the Poisson equation
not be consistently oriented. (We also consider the more general
q = r · .C on "
case where ⌦ can include isolated points—see Appendix A.1.) We mq (3)
use # to denote the unit normals of ⌦, and = for the unit normals m= = = · .C on m",
of the domain boundary m". We use ` ⌦ to denote a measure con- whose solution is determined up to a constant shift. To exactly in-
centrated on ⌦, similar in spirit to an indicator function1 ; # ` ⌦ terpolate the input (and get a unique solution) we could also require
is likewise a vector �eld (or vector measure) equal to zero away that q = 0 along ⌦—though if the input is corrupted, interpolation
from ⌦, and determined by # for points in ⌦. Finally, we use may be ill-advised (Figure 3); see Section 7.1 for further discussion.
for the negative-de�nite Laplace-Beltrami operator on ", and r In summary, we arrive at the following algorithm:
for the negative-de�nite connection Laplacian [Gallier et al. 2020];
intuitively, these operators measure the deviation of a scalar func- (1) Solve a vector di�usion equation 3C
3 - = r - (Equation 1).
C C
tion and tangent vector �eld (resp.) from their average in a local (2) Evaluate the vector �eld .C = -C /k-C k.
neighborhood. (3) Solve a Poisson equation q = r · .C (Equation 3).

4 TIME DISCRETIZATION
Ø
1 More formally, for any Borel measurable set * ⇢ " , ` ⌦ (* ) := ⌦\*
3+ , where As in past heat methods we discretize Equation 1 in time via one
3+ is the usual volume measure on ⌦. step of backward Euler [Crane et al. 2013b; Sharp et al. 2019c], and

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
 A Heat Method for Generalized Signed Distance • 92:5

solve a linear equation A corresponding basis for vector �elds is

expressed by identifying tangent vectors with
(id C r
)-C = - 0 (4) complex numbers. In particular, at each edge
for a single, �xed time C > 0 (where id is the identity). To get a feel 89 we choose a coordinate system such that 1
for the behavior of this equation, we can consider the Euclidean and the imaginary unit y correspond to 4̂89 and
case (" = R= ), where applying r to a vector �eld is equivalent 4̂89? , resp. The function k89 := i89 + 0y then de-
to applying the scalar Laplacian to each component. Then as �nes a basis vector �eld parallel to the edge, and
C ! 1 we approach a standard Poisson equation; moreover, since Ik89 describes a locally supported vector �eld
- 0 = # ` ⌦ is zero almost everywhere, the solution will look nearly parallel to I 2 C (right inset).
harmonic away from ⌦.
5.3 Edge-Based Laplacian and Mass Matrix
5 SPATIAL DISCRETIZATION If we use the CR bases to discretize the scalar Laplacian (via a
To implement the signed heat method on any geometric data struc- standard Ritz-Galerkin approach), we get a real-valued positive
ture, one need only provide sparse matrices representing the Lapla- semide�nite matrix L 2 R |⇢ | ⇥ |⇢ | with nonzero entries
cian, connection Laplacian, and divergence operator. Section 8 ex- L89,9: = 2 cot \ :8 889 2 ⇢, 89: >89,
9 ,
plores several possibilities (regular grid, polygon mesh, point cloud, 9: (6)
2 cot \8 ,
Õ
L89,:8 = 889 2 ⇢, 89: >89,
etc.); in this section we focus on triangle meshes.
L89,89 = 89:>89 L89,9: + L89,:8 , 889 2 ⇢,
Here we use edge-based operators for the vector di�usion step
(Step I), which makes it straightforward to discretize curve sources; where either the length/area or angle-based expressions can be used.
we then use standard vertex-based operators for the Poisson equa- The mass matrix M 2 R |⇢ | ⇥ |⇢ | is diagonal, with nonzeros
tion (Step III), so that our �nal SDF is stored at vertices. As in Õ
Djerbetian and Chen [2016] we adopt a complex encoding of the M89,89 = 13 89:>89 |89: |, 889 2 ⇢. (7)
connection Laplacian), which uses half as much storage/bandwidth
as the real-valued version used by Stein et al. [2020]: each 2x2 real 5.4 Edge-Based Connection Laplacian
block (four �oats) is replaced by a single complex value (two �oats). Intuitively, the connection Laplacian mea-
sures the deviation of a vector �eld from
5.1 Notation its average value in a small local neighbor-
Our domain is a triangle mesh " = (+ , ⇢, ) with hood. To compute this average on a mesh,
arbitrary topology (e.g., nonmanifold, nonorientable, we must parallel transport vectors into a
and/or with boundary); we use ⇠ for the set of all common coordinate system at each edge
triangle corners. We denote :-simplices by (: + 1)- (see inset). Since edges are intrinsically �at,
tuples of vertex indices, i.e., vertices 8 2 + , edges parallel transport is simply a translation—
89 2 ⇢, and faces 89: 2 . Likewise, we denote the but we must also account for a change of
corner of triangle 89: at vertex 8 as 9: coordinates. In particular, for each triangle 89:, we encode the paral-
8 2 ⇠. These indices are also used lel transport of a vector from edge 89 to the next edge 9: in counter-
to express quantities stored on mesh elements—for instance, corner
9:
angles are denoted by \8 . We use < and > to indicate summation clockwise order as a rotation
over all elements contained by or containing another element (resp.).
Õ
:8 :8
A89!9: := B89!9: 4 y (c \ 9 ) = B89!9: 4y\ 9 , (8)
For instance 89:>89 sums over all triangles 89: containing edge 89.
We use | · | to denote the volume of a simplex—e.g., |89 | is an edge where the sign B89!9: 2 {+1, 1} is equal to 4y0 = +1 if 89 and 9:
length and |89: | is a triangle area. To avoid trigonometric functions, have the same relative orientation with respect to 89:, and 4yc = 1
we evaluate all angle cotangents via the formula otherwise. The positive-semide�nite connection Laplacian matrix
Lr 2 C |⇢ | ⇥ |⇢ | is then obtained by multiplying o�-diagonal entries
cot \8 = (|89 | 2 | 9: | 2 + |:8 | 2 )/(4|89: |).
9:
(5) of L by edge rotations A89,9: , yielding nonzero entries
For each edge 89, we let 489 be a vector parallel to the edge with Lr 89,9: = 1
2A89!9: cot \ :8 889 2 ⇢, 89: >89,
9 ,
arbitrary (but �xed) orientation, and magnitude equal to the edge
Lr = 2A:8!89 cot \8
9:
889 2 ⇢, 89: >89, (9)
length. We let 489? be the 90 rotation of 489 in the counter-clockwise 89,:8
Lr 89,89 = L89,89 , 889 2 ⇢,
direction, and use 4̂89 , 4̂89? for the corresponding unit vectors.
with cotans evaluated as in Equation 5. Since A89!9: = A 9:!89 1 =
5.2 Edge Basis Functions A 89!9: , the matrix L is Hermitian. It is positive semide�nite because
r
The Crouzeix-Raviart (CR) basis functions it is the Hessian of a convex energy [Stein et al. 2020, Equation 5].
associate each edge 89 2 ⇢ with a face-
wise linear function i8 9 : " ! R inter- 5.5 Vertex-Based Laplacian and Divergence
polating the value 1 at the midpoint of 89, Our vertex-based operators are identical to those used in Step III
and 0 at all other edge midpoints (inset). of the unsigned heat method [Crane et al. 2013b, Section 3.2.1]. In

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
92:6 • Feng and Crane

Fig. 8. The unsigned heat method exhibits bias near the domain boundary
for large di�usion times (figure reproduced from Crane et al. [2013b], Figure
11 using C = 100⌘ 2 ). Using their proposed boundary condition heuristic
only slightly improves results. In contrast, our method has correct boundary
conditions.

Fig. 9. We can mix and match signed and unsigned distance, selectively
particular, we use the cotan Laplacian C 2 R |+ | ⇥ |+ | [Crane et al. treating open curves as either broken region boundaries—or as literal open
2013a, Chapter 6], which has nonzero entries curves. We can also incorporate distance to isolated points.
1Õ 89
C8,9 = 2Õ 89: >89 cot \ : , 889 2 ⇢ (10)
C8,8 = 89 >8 C8,9 , 88 2 + . than vectors in R3 , by writing all vectors as di�erences of points in
barycentric coordinates (see Sharp et al. [2021, Section 3.1]).
At interior vertices 8 2 + , the divergence of a per-
The formula in Equation 11 is obtained by projecting the contin-
face vector �eld .89: is
’ uous measure # `W onto the vector CR basis. A nice feature of CR
ô
(r · . )8 := 12
89
(cot \: 4 89ô +cot \ :8 bases is that they are orthogonal. Hence, for a given segment W we
9 4 :8 ) · .89: ,
ô

8 9:>8 need only integrate the function k89 with respect to the measure, or
where 4 89ô denotes the vector from vertex 8 to vertex 9 (see Polthier equivalently, take a Hermitian inner product along the segment:
π π π
and Preuß [2003, Section 4]). For boundary vertices, the discrete
divergence includes (due to integration by parts) an additional = · . k 89 #W d`W = k 89 #W dB = #W k89 dB = #W |W |i89 (<W ).
" W W
term—which we omit because it exactly cancels with our desired
These equalities hold since #W is constant and k89 = i89 + 0y is linear
boundary conditions in Equation 3.
along W. Equation 11 expresses the �nal quantity in the edge basis.
5.6 Source Discretization 6 ALGORITHM
We discretize the initial conditions from Equation 1 as values (X0 )89 2
Our discrete algorithm amounts to solving two sparse linear systems
C at edges 89, where the source set ⌦ can be a mix of oriented and
using the discrete operators de�ned in Section 5. First we solve the
unoriented curves, as well as isolated points (Figure 9). Here we
discrete vector heat equation,
treat oriented curves ; Appendix A treats other source types.
In principle our formulation works for arbi- (M + CLr )X = X0, (12)
trary curves, though to establish explicit formu- obtaining a di�used vector �eld X. Following Sharp et al. [2019c,
las we will assume is comprised of straight seg- Section 7.3], we let C = ⌘ 2 , where ⌘ is the mean distance between
ments, each of which is contained entirely inside nodes—in our case, edge midpoints, yielding half the mean edge
one face or edge. Note that when the domain " length.
is orientable, the normal # to is uniquely de- Next, we average the di�used vectors X to each
termined by a 90 counter-clockwise rotation; on face 89: 2 via -89: := (-89 + - 9: + -:8 )/3 (taking
nonorientable domains, one must explicitly specify normals. care to express all vectors in the same basis, as
For each triangle 89: 2 , a segment W contained in its interior in Section 5.4), and compute unit vectors .89: :=
contributes to the initial values at all three of its three edges. For X89: /kX89: k which represent the gradient of our
instance, the contribution to (- 0 )89 is given by the complex value (generalized) SDF. Finally, to obtain the SDF q 2 R |+ | at vertices,
⇣ ⇣ ⌘ ⌘ we solve a sparse linear system
|W | #W · 4̂89 + #W · 4̂89? y i89 (<W ) (11)
Cq = b (13)
where |W |, #W , and <W are the length, normal, and midpoint of
where C is the cotan Laplacian (Section 5.5), b 2 R |+ | is the vector
the segment, resp. (and similarly for (- 0 ) 9: , (- 0 ):8 ). Empirically,
of discrete divergences given in Section 5.5.
however, results are more accurate if we omit the factor i89 (<W ). If
W runs along an edge 89, then it makes the same contribution, but Boundary Behavior. Unlike the unsigned heat method, our signed
only to (- 0 )89 . These contributions are summed over all segments to heat method exhibits the correct behavior at the boundary (Figure 8),
obtain �nal values for - 0 . For intrinsic retriangulation (Figure 22), without any special boundary treatment (as in Edelstein et al. [2023,
Equation 11 can be expressed using purely intrinsic data, rather Section 4.2]). The reason is that UHM obtains the vector �eld - as

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
 A Heat Method for Generalized Signed Distance • 92:7

no constraints constrain all to zero constant per component

Fig. 11. So� level set constraints can easily be incorporated into our PDE-
based framework. Here so� constraints encourage the zero set to lie near
the input curves, without distorting the SDF near self-intersections as much
as a hard constraint.

saddle-point problem
Fig. 10. We can also fit a signed distance function to several partial level sets.   
Le�: Without constraints, isolines deviate slightly from source geometry. C AT q r · YC
= , (14)
Center: Blindly constraining to zero along all curves grossly violates the A 0 ` 0
distance property. Right: Constraining values to be constant along each
where ` 2 R< are Lagrange multipliers. Since q | ⌦ is now constant
component nicely matches the input geometry.
by construction, we can shift it to be exactly zero on all of ⌦—so
long as the constraints are compatible. We apply this procedure
the gradient of a scalar heat distribution D with either zero-Dirichlet in all �gures unless otherwise noted. More generally, suppose the
or zero-Neumann boundary conditions—in either case, the gradient input ⌦ represents multiple, distinct level sets q 1 (2 1 ), q 1 (2 2 ), . . .
of D cannot point in the right direction (either purely normal or for values 28 which are not known a priori. Here we can apply an
purely tangential, resp.); Crane et al. [2013b, Section 3.4] suggests identical set of constraints per connected component, ensuring that
to simply take a �xed linear combination. In contrast, even at the the value along each component is constant (Figure 10, right).
boundary our vector di�usion step directly provides the normal at Finally, similar to Kazhdan and Hoppe [2013, SectionØ4], we can
the closest point, which agrees with the gradient of the true SDF. The replace the hard constraint with a soft penalty term ⌦ kq (G)
basic reason is that the discrete connection Laplacian (Equation 9) 2 (G)k 2 dG, where _ > 0 is the penalty strength (Figure 11). The
encodes zero-Neumann boundary conditions on the vector �eld minimizer is obtained by solving the system (C _A> A)q = b.
itself [Gelfand et al. 2000, I.6], rather than a scalar potential D.
7.2 Nonbounding Loops and Discontinuous Distance
7 OPTIONAL EXTENSIONS The gradient of unsigned distance
Here we discuss some optional extensions of our basic method. can be discontinuous on a lower-
These extensions take advantage of the global variational nature dimensional set (analogous to the cut lo- nonbounding
of our method to provide capabilities not possible with iterative cus of a point), but is still tangent contin-
region-growing methods based on MMP [Mitchell et al. 1987] or uous relative to this set, meaning vectors
fast marching [Kimmel and Sethian 1998]. on either side project to the same direc-
tion (Figure 13, top). In contrast, when
7.1 Preserving Level Sets ⌦ does not bound a region of ", the vec-
tor �eld .C from Step II of our method
The �nal Poisson solve (Equation 13) recovers the signed distance
can fail to be tangent continuous—and
function q only up to a constant shift. To make the zero level set
hence fail to be integrable via a globally
of q approximate the source geometry ⌦, a common strategy is
continuous function (Figure 13, bottom).
to shift q by its average over ⌦ [Kazhdan et al. 2006; Calakli and
One possibility is to simply �lter out
Taubin 2011; Crane et al. 2013b]. An alternative in our setting is to Fig. 12. GSD extends to
nonbounding components [Feng et al. curves that do not even
add an explicit linear constraint to Equation 2 that ensures q takes
2023]; we instead consider a generalized bound a region.
the same value at all points of ⌦.
notion of signed distance that remains
On a surface mesh, for instance, suppose ⌦ has
meaningful. In particular, we now seek a piecewise continuous solu-
at most one segment per triangle, with endpoints
tion by replacing the ! 2 problem in Equation 2 with an ! 1 problem
W 0, . . . , W< on edges. Then we can impose linear con-
that ignores neighborhoods where .C is highly nonintegrable. On
straints of the form
triangle meshes, we quantify nonintegrability via an edge-based
(1 C? )q8? + C? q 9? = (1 C 0 )q8 0 + C 0q 90 , ? = 1, . . . , <, curl inspired by Polthier and Preuß [2003, Section 4]:
where C? 2 [0, 1] encodes the location of G? along edge 8? , 9? (see (r ⇥ Y)89 := (4̂89 + 4̂ 98 ) · Y89: . (15)
inset). We encode these constraints by a matrix A 2 R<⇥ |+ | . Mini- This quantity directly captures tangent discontinuity along edges:
mizing krq .C k 22 subject to Aq = 0 then corresponds to solving a taking edge orientation into account, it is simply the di�erence

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
92:8 • Feng and Crane

tangent inconsistently-oriented curves

continuous standard consistently-
integration (L2) oriented curves
piecewise continuous
integration (L1)
ideal
solution

tangent
discontinuous
Fig. 13. Unlike the gradient of an unsigned distance function (top), the unit integration: standard piecewise
continuous
vector field we compute for signed distance might not be integrated by any
continuous function (bo�om). In such scenarios, we can instead compute a Fig. 14. Our algorithm works out-of-the-box for orientable input curves,
piecewise continuous SDF (right). even on nonorientable domains (top). By adopting piecewise continuous
integration (Section 7.2) we can also handle inconsistently-oriented curves,
whether or not the domain itself is orientable (bo�om).
of projected vectors. Discontinuous functions q are represented
as piecewise linear functions interpolating values at corners. To
integrate Y in a piecewise continuous sense, we then minimize a
weighted ! 1 norm that penalizes discontinuous jumps across edges: Note that robust sign information is available only thanks to our
’
min |89 |4 | r⇥Y|89 |q8
9: ;9
q8 | method—Equation 18 cannot be applied directly to broken curves.
q 2R |⇠ | 89 2⇢ In practice ADMM is quite slow, requiring repeated linear solves.
(16)
9: Since we have a good initializer, it might be better to use a �rst-order
s.t. q :8
9 q8 = Y89 · 4 89ô, 8 9:
8 2 ⇠, method like the primal-dual hybrid gradient method, which is often
Exponential weights incur a lower penalty when the �eld is not faster than ADMM [Chambolle and Pock 2011, Figure 5] and easily
integrable; simultaneously, the ! 1 norm tries to minimize the length implemented in parallel. We leave such exploration to future work.
of this discontinuity as much as possible. The constraints ensure
Y is integrated exactly within each triangle. In practice, we apply 8 OTHER SPATIAL DISCRETIZATIONS
the same transformations as Feng et al. [2023, Section 3.4] to obtain 8.1 Tetrahedral Meshes
a linear program with only | | degrees of freedom. Examples are
The discretization from Section 5 generalizes directly to tet meshes—
shown in Figure 13 (right) and Figures 14 and 12. In Figure 12 we also
extending our method from curve processing to broader surface
use a small amount of heat �ow to smooth out minor discontinuities
processing tasks. We again use a Crouzeix-Raviart �nite element
from !1 optimization.
discretization, detailed in Appendix B. Note that we do not need
7.3 Distance Sharpening a separate discretization for the connection Laplacian, since on a
�at domain we can just apply the scalar Laplacian componentwise.
Unsigned geodesic distance can also be expressed as the solution to For simplicity, we generate a mesh that conforms to input triangles
a convex optimization problem, akin to the convex formulation of (via TetGen [Hang 2015]), so that the source term is just - 0 (G) :=
graph distance [Dantzig 1963, Ch. 17], [Erickson 2019, Ch. H]: Õ
89: 2⌦ |89: |#89: k89: (G). Examples are shown in Figures 19, 26.
π
max q (G) dG
q "
(17)
s.t. |rq |  1 convex formulation signed heat method
q = 0 on ⌦. solve
vex
con
Belyaev and Fayolle [2020] solve Equation 17 via ADMM to compute
the distance to point sources. This formulation tends to be more
accurate than our method—but is an order of magnitude slower
(Section 9.5), and more importantly, can compute only unsigned
distance. If desired, however, one can “sharpen” our results using a
generalization of Equation 17. We simply replace the objective in solve time: 0.70s before sharpening after sharpening
Equation 17 with solve time: 0.51s additional time: 0.66s
π Fig. 15. Le�: Methods based on convex optimization yield more accurate
max sign(q 0 (x))q (x) dx (18) distances, but compute only unsigned distance. Right: Using our method
q " as a warm start, we can “sharpen” distance while preserving the inside-
where q 0 is the distance computed by the signed heat method, and outside classification. Here we start with a large di�usion time (C = 100⌘ 2 )
use q 0 as an initial guess for q. An example is shown in Figure 15. to visually emphasize the e�ect.

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
 A Heat Method for Generalized Signed Distance • 92:9

Fig. 16. Our method extends to polygon meshes, point clouds, and digital surfaces. Digital surface meshes are from Coeurjolly and Levallois [2015].

8.2 Regular Grids 8.4 Point Clouds

To avoid the cost of tet mesh generation, we can also use a regular For point clouds, we use the Laplacian of
grid on R3 . Here, rather than solve a linear system in Step I, we Sharp and Crane [2020a, Section 5.7], which
can directly approximate -C via the fundamental solution of the also de�nes per-point tangent planes, and the
operator id C , given in R3 by the Yukawa potential mass matrix from their tufted cover [Sharp and
p Crane 2020a, Section 3.3]. We again augment
⌧C (G, ~) = 4c kG1 ~ k 4 1/C kG ~ k . this Laplacian with o�-diagonal rotations to
get a connection Laplacian, as described by Sharp et al. [2019c, Sec-
(In practice, we �nd this approach works better than convolving tion 6.2.2] and implemented in geometry-central [Sharp et al. 2019a].
with the heat kernel directly.) The di�used vector �eld at each Each component of is a sequence of points W8 from the cloud,
grid node G 2 R3 is then approximated by midpoint quadrature at where normal #8 can be computed by projecting W8 1, W8 , W8+1 onto
triangle barycenters 189: : the tangent plane at point 8 and taking the in-plane normal (inset).
’ π ’ As in Section 8.3, we compute - 0 by accumulating length-weighted
-C (G) = N89: ⌧C (G, ~) 3~ ⇡ |89: |N89: ⌧C (G, 189: ), normals at each point. To compute divergence, we average di�used
89: 2 89: 89: 2 vectors onto faces of the tufted cover and apply the formula from
Section 5.5. An example is shown in Figure 16, center.
where N89: are triangle normals; this step could be accelerated from
$ (| |) to $ (log | |) via the method of Barnes and Hut [1986]. Step
III uses a standard seven-point Laplacian, and evaluates r · YC using
8.5 Voxelizations
forward di�erences. Figure 4 shows that we obtain high-quality Our method also extends to digital surfaces, which are the bound-
distances and o�set surfaces even for badly corrupted geometry aries of voxelized volumes. We use the vertex-based Laplacian, diver-
(with contouring via marching cubes [Lorensen and Cline 1998]). gence, and connection Laplacian de�ned by Coeurjolly and Lachaud
[2022], using the implementation in DGtal [Coeurjolly et al. 2010].
8.3 Polygon Meshes Analogous to Section 8.3, we encode the source term via extrinsic
projection of curve normals onto tangent bases de�ned at vertices,
On surface meshes made of general polygons
where vertex normals are normalized averages of per-face corrected
(which may be nonconvex or even nonplanar),
normals [Coeurjolly et al. 2014]. We assume that is piecewise
we adopt the Laplacian, divergence, and mass ma-
linear on faces of the digital surface, and obtain initial conditions - 0
trices de�ned by Bunge et al. [2020] in terms of
by projecting segment normals onto tangent bases at vertices, and
a carefully-chosen triangular re�nement. Since
multiplying by half the edge length of incident segments. Examples
Bunge et al. do not de�ne a connection Laplacian, we de�ne tan-
are shown in Figure 16, right.
gent planes at vertices using the area-weighted sum of �ne triangle
normals (picking an arbitrary basis in each plane). As in Equation 9,
a connection Laplacian is then obtained by multiplying o�-diagonal 9 EVALUATION
entries with complex numbers that encode the smallest rotation be- We next evaluate our method via a broad range of numerical exper-
tween tangent bases at adjacent vertices. For simplicity we assume iments. As noted in Section 1, our primary focus is not on speed or
the curve runs along edges; the initial conditions - 0 at each vertex accuracy, but rather on improving the robustness and generality of
8 of are then given by the average edge normal # times half the signed distance computation. Nonetheless, our signed heat method
sum of incident edge lengths ✓1, ✓2 (inset). To compute divergence, (SHM) is quite competitive in speed and accuracy, since it boils down
we prolong the di�used vector �eld X to the �ne vertices, average to the same kinds of sparse linear systems as previous methods. In
them onto faces, then apply the operator D from Bunge et al. [2020, particular, the method takes only a few hundred milliseconds for
Section 4.5]. Examples are shown in Figure 16, left. meshes with hundreds of thousands of triangles, and exhibits the

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
92:10 • Feng and Crane

ground truth topological errors

geometric errors

orientation errors red = flipped orientation

Fig. 17. Our method provides robust and reliable signed distance approximation, failing gracefully in the presence of significant topological, geometric, or
orientation errors. Errors n in geodesic distance are displayed relative to the exact polyhedral SDF of a finely sampled version of the original curve.

same linear convergence with respect to mesh re�nement as earlier Zhang 2010] using marching triangles, and remesh to a constrained
schemes. More importantly, for broken geometry it still provides an intrinsic Delaunay triangulation containing these curves [Sharp and
accurate approximation of the distance to the original, uncorrupted Crane 2020b, Section 6.3].
geometry—whereas past methods can fail to provide a reasonable
distance approximation (Figure 7, right). 9.3 Examples
9.3.1 Morphological Operations. One natural use case for general-
9.1 Implementation ized signed distance is to contour broken geometry, and generate
All methods were implemented in C++, in double precision, using accurate �xed-distance o�sets (Figures 26, 19, 4). Generating such
geometry-central [Sharp et al. 2019a] for mesh data structures and o�sets from imperfect geometry is useful, for example, for 3D print-
intrinsic retriangulation, and Eigen [Guennebaud et al. 2010] for ing, or for downstream mesh processing tasks that require closed
sparse linear algebra. Timings were taken on an Apple M1 with 8GB or manifold surfaces. One can also “in�ate” or “shrink” shapes by
RAM. The error n in any distance approximation q is quanti�ed via taking positive or negative o�sets—and combining these two opera-
!2 error normalized by the range of the true distance q 0 , i.e., tions in sequence can be used to simplify high-frequency features
✓ ⇣ ⌘ 2 ◆ 1/2 of broken shapes as if they were whole (Figure 20). Note that in
1 Õ 0 Figure 26 we sample both functions onto the same tet mesh, using
n (q) := (max(q 0 ) min(q 0 ) ) A 1/2 8 2+ 8 q8 q8 ,
libigl to evaluate GWN [Jacobson et al. 2018].
Õ
where 8 = 13 89:>8 |89: | is the area associated with vertex 8, and
Õ 9.3.2 Illustration on Surfaces.
A := 8 2+ 8 is the total surface area. Geodesic distance is also
used to design curves on
9.2 Dataset surfaces—where prior work
For evaluation, we consider a largely considers perfect closed
dataset of closed, region-bounding curves [Nazzaro et al. 2022].
loops (since past methods do Our algorithm can be used to
not handle open or nonbounding pre-process curves into the smaller t larger t
curves). This dataset is derived requisite closed format, e.g.,
from all 44 genus-zero meshes those hastily sketched on a Fig. 18. Adjusting di�usion time fills in
without boundary from Myles surface (Figure 21). To con- broken le�ers with either round or sharp
corners, yielding e�ects similar to di�er-
et al. [2014]. In order to perform trol the behavior of this com-
ent line joins for 2D strokes.
a study of convergence under re�nement (Figure 25) we �rst remesh pletion operation, one can
each model to approximately 1.25k vertices using quadric error sim- use di�usion time to control completion behavior, providing an
pli�cation [Garland and Heckbert 1997], then compute four levels analog of line join options from 2D vector graphics (Figure 18).
of loop subdivision [Loop 1987]. At each level, we then extract level
sets of the same low-frequency Laplacian eigenfunctions [Lévy and

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
 A Heat Method for Generalized Signed Distance • 92:11

Fig. 19. By extracting level sets of generalized signed distance, we can convert broken, noisy, and nonmanifold input geometry (far le�) into closed, regular,
manifold surfaces and evenly-spaced o�set surfaces.

Fig. 21. Broken curves easily arise from a�empting to draw curves on sur-
faces of high genus, with overhangs, and with holes and scanner noise.
Our method yields signed distance functions robust to these challenges.
(Scanned bench from from Choi et al. [2016].)

Fig. 20. One can simplify high-frequency features of broken curves and
surfaces by taking consecutive positive/negative o�sets of a generalized 9.5 Accuracy and Performance
signed distance function, akin to dilation/erosion.
We compare against the unsigned heat method (UHM) of Crane
et al. [2013b], which provides a useful reference point since nearly
all work on geodesic distance algorithms from the past decade com-
9.4 Robustness pares against this method. As a more recent reference point, we also
compare with the convex formulation of Belyaev and Fayolle [2020]
Since our method is built on (labeled BF). Since neither method directly handles curve sources,
well-behaved elliptic PDEs, it per- we either integrate the initial scalar heat distribution against hat
forms well on not only broken functions (for UHM) or simply use the set of curve vertices as the
input with corrupted topology, source set (for BF). (Methods that directly handle curve sources do
geometry, and orientation (Fig- not have an open source implementation [Bommes and Kobbelt
ure 17) but also challenging surface domains (Figure 23). Since our 2007], or do not include curve sources in their public release [Tret-
method is purely intrinsic, it applies to meshes that are only im- tner et al. 2021].) Finally, since BF must constrain the zero set, we
mersed (rather than embedded); it also applies out of the box to impose the same constraints on UHM/SHM (Section 7.1), and do
nonmanifold and nonorientable meshes (inset), since all our di�er- not pre-factor any matrices. Note, however, that for multiple source
ential operators are local and de�ned per-face, and hence oblivious
to any nonmanifold features. Our method is robust across varying
degrees of nonmanifoldness and missing data (Figure 23, top and
bottom left). As with all methods that rely on discretizing PDEs, the
quality of the solution can degrade with poor tessellations of the
geometry, though we can easily apply intrinsic Delaunay triangula-
tion to get good-quality solutions (Figure 22). For surfaces in R3 our
method remains robust, even for extremely corrupted input surfaces
(Figures 4,19).
As seen in Figure 27 we also obtain more natural surface com-
pletion than GWN, which for general surfaces is hard to contour
with any single level set value. Similar to spline interpolation, the
zero set of GSD nicely matches both positions and normals along Fig. 22. The quality of our method depends on the underlying mesh quality;
hole boundaries. Here meshes for both methods are extracted using but since our formulation is purely intrinsic, we can trivially improve accu-
the marching tetrahedra implementation in libigl [Jacobson and racy and robustness by invoking intrinsic Delaunay refinement [Gillespie
Panozzo 2017]. et al. 2021], without changing anything else about our implementation.

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
92:12 • Feng and Crane

more nonmanifold nonmanifold & self-intersecting

GSD

ground
truth

large holes

Fig. 23. Our method robust not only to errors in the source geometry, but also in the domain mesh itself. Here we obtain well-behaved SDFs even for meshes
found “in the wild,” such as amateur-created 3D scans [Choi et al. 2016]. Even in cases where a notion of inside and outside is meaningless (such as the
rightmost mesh), our method fails gracefully—still producing a good signed distance approximation near the input curve.

terms, heat methods can achieve about two orders of magnitude truth” distance as the exact polyhedral distance to a �nely-sampled
speedup by omitting factorization [Crane et al. 2013b, Table 1]. version of the input curve (100 samples per curve edge) using MMP
[Mitchell et al. 1987], which itself has $ (⌘ 2 ) error. We plot conver-
9.5.1 Planar Domains. As noted by Crane et al. [2013b, Figure 21],
gence and solve times in Figure 25. We observe that our method has
even exact polyhedral distance (including MMP) provides only a
approximately linear convergence in mean edge length, with better
2nd-order accurate estimate of true (smooth) geodesic distance, due
consistency on curve sources compared to other methods. We �nd
to errors in the approximation of the domain itself. To avoid con-
the same trend in solve times as in Section 9.5.
�ating these two sources of error, we �rst consider closed, planar
curves—where the exact SDF is easily computed via closest-point 9.5.3 Broken Curves. Finally, we compare against the end-to-end
queries [Sawhney et al. 2020], and sign can unambiguously be de- pipeline of repairing broken geometry, computing unsigned dis-
termined via standard inside/outside tests [Haines 1994]. As seen in tance to the �xed geometry, and signing the unsigned distance. In
Figure 24, our method is slightly slower but slightly more accurate particular, we use surface winding numbers (SWN) to contour broken
than UHM. Without sharpening, it is not as accurate as BF—but is an curves [Feng et al. 2023], and compute exact polyhedral distance us-
order of magnitude faster. Moreover, BF must trade o� between bias ing MMP [Mitchell et al. 1987] using the curve vertices as the source
near the boundary [Edelstein et al. 2023, Figure 3, left], or distortion set. As input surface domains, we use the meshes with ⇠5k vertices
in the presence of curve sources, depending on whether Hessian from the same dataset as Section 9.5.2. As input curves, we take
regularization is omitted or included (resp.). level sets of �ve di�erent low-frequency Laplacian eigenfunctions,
9.5.2 Surface Domains. We next consider closed curves on sur- and add geometric and topological errors by taking the union of the
face meshes. Here we can no longer obtain the true distance on an curves with their o�sets (found by taking boundaries of triangle
unknown underlying smooth surface; we hence compute “ground strips), and deleting about 50% of the curve at random intervals.
The repair-distance-sign pipeline is particularly sensitive to er-
rors in the input, since any errors made during contouring are
permanent and will destroy the quality of the �nal SDF no matter
how accurate the subsequent distance computation. In particular,
contouring the winding number is notoriously di�cult, and often
0.37s 0.14% 0.17s 0.15% 1.89s 0.13%
SHM (ours) UHM BF

0.27s 0.38% 0.11s 0.57% 1.04s 0.16%

Ours (0.83) FMM (0.23) UHM (0.59) ADMM-BF (0.70)

Fig. 24. Distribution of error in distance approximation for a perfect, unbro-
ken curve on a high- and low-quality mesh (top/bo�om). Both meshes have Fig. 25. We observe approximately linear convergence in distance accuracy
about 100k faces. Inset numbers on SDF and error plots indicate compute on a benchmark of unbroken (closed) curves on 44 di�erent meshes. The
time and mean error (resp.). Overall our method is quite comparable to legend shows median orders of accuracy. Note that if we omit the zero set
the original heat method, and less accurate than BF on the low-quality constraint, enabling us to re-use both factorizations, our method and UHM
mesh—but about 4–5x faster. become 1–2 orders of magnitude faster.

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
 A Heat Method for Generalized Signed Distance • 92:13

ours winding number winding number

Fig. 26. Top: Generalized winding number (GWN) cannot be used for o�set (alternate contouring)
surfaces, since it provides only a smooth indicator function—and not a
signed distance. Bo�om: Our generalized signed distance (GSD) provides Fig. 28. We compare, on 220 examples, the accuracy of GSD versus a hybrid
much nicer o�sets on the same broken mesh. scheme that repairs broken geometry using winding numbers, computes
unsigned distance to the repaired geometry, then signs the unsigned dis-
tance. Winding numbers are worse at classifying inside/outside, regardless
of contouring method. As a result, the hybrid scheme yields about 3x lower
generalized signed distance distance accuracy on average (0.11 vs. 0.04 !2 error, resp.), even though it
offset benefits from exact distance (via MMP).
input zero set surfaces

10 LIMITATIONS AND FUTURE WORK

As with past heat methods, di�usion time cannot be made arbitrarily
small: similar to Crane et al. [2013b, Appendix A], we observe em-
pirically that GSD behaves like graph distance as C ! 0 (see inset).
generalized winding number As with past heat methods, however, SHM
consistently provides accurate distance ap-
proximation for C = ⌘ 2 . Our strategy for topo-
logically nonbounding and nonorientable
curves (Section 7.2) incurs a more expensive ours

bulges bulges ! 1 problem—but is already a substantial gen-

bulges in
out in eralization over anything handled by previ-
ous geodesic or SDF algorithms. A more care-
UHM
Fig. 27. GWN completes surfaces with saddle-shaped harmonic patches ful study of this situation (including mitiga-
that exhibit poor normal continuity with the observed geometry (across tion of small per-edge discontinuities) is an
many contour values). Our method directly incorporates normal information, interesting direction for future work. Like-
providing more plausible reconstruction even for large holes. wise, since broken input sometimes comes
with quanti�able measurement error, it might be interesting to ex-
plore a notion of “con�dence” in the computed distance values,
along the lines of work by Sellán and Jacobson [2022].
leads to misclassi�ed regions (Figure 28). Though Feng et al. [2023] Finally, like other robust inside-outside and reconstruction algo-
suggest a rounding procedure, it remains unclear which half-integer rithms [Bærentzen and Aanæs 2005; Kazhdan et al. 2006; Jacobson
level set to take as the boundary between inside and outside; we et al. 2013; Feng et al. 2023], we assume input has mostly consistent
use the values of the rounded winding number function that appear orientation. Though we are fairly robust to orientation errors (Sec-
most often along the input curves, though this results in 12% of tion 9.4), it may be interesting to explore di�usion of more general
examples with >50% of surface area misclassi�ed. We also try con- line �elds to factor out normal orientation, à la Alliez et al. [2007].
touring winding number according to the average of the winding
number function along the input curves, but we �nd similar results. ACKNOWLEDGMENTS
In contrast, our method achieves greater robustness by averaging This work was funded by an NSF CAREER Award (IIS 1943123),
normal vectors over the whole domain. The repair-distance-sign NSF Award IIS 221229, a Packard Fellowship, and gifts from Adobe
pipeline is also about 10x more expensive due to its intermediate Research. The authors also thank Alec Jacobson for helpful conver-
steps; SWN alone has cost asymptotically equivalent to our method. sations about winding numbers and the pseudonormal.

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
92:14 • Feng and Crane

REFERENCES 10.1145/3592401
Marc Alexa, Philipp Herholz, Max Kohlbrenner, and Olga Sorkine-Hornung. 2020. Jean Gallier, Jocelyn Quaintance, Jean Gallier, and Jocelyn Quaintance. 2020. Op-
Properties of Laplace Operators for Tetrahedral Meshes. Computer Graphics Forum erators on Riemannian Manifolds: Hodge Laplacian, Laplace-Beltrami Laplacian,
(2020). [Link] the Bochner Laplacian, and Weitzenböck Formulae. Di�erential Geometry and Lie
P. Alliez, D. Cohen-Steiner, Y. Tong, and M. Desbrun. 2007. Voronoi-Based Variational Groups: A Second Course (2020), 361–401.
Reconstruction of Unoriented Point Sets. In Proceedings of the Fifth Eurograph- Michael Garland and Paul S. Heckbert. 1997. Surface simpli�cation using quadric error
ics Symposium on Geometry Processing (Barcelona, Spain) (SGP ’07). Eurographics metrics. In Proceedings of the 24th Annual Conference on Computer Graphics and
Association, Goslar, DEU, 39–48. Interactive Techniques (SIGGRAPH ’97). ACM Press/Addison-Wesley Publishing Co.,
Matan Atzmon and Yaron Lipman. 2019. SAL: Sign Agnostic Learning of Shapes from USA, 209–216. [Link]
Raw Data. CoRR abs/1911.10414 (2019). arXiv:1911.10414 [Link] Izrail Moiseevitch Gelfand, Richard A Silverman, et al. 2000. Calculus of variations.
10414 Courier Corporation.
J.A. Bærentzen and H. Aanæs. 2005. Signed distance computation using the angle Mark Gillespie, Nicholas Sharp, and Keenan Crane. 2021. Integer Coordinates for
weighted pseudonormal. IEEE Transactions on Visualization and Computer Graphics Intrinsic Geometry Processing. ACM Trans. Graph. 40, 6 (2021).
11, 3 (2005), 243–253. [Link] Amos Gropp, Lior Yariv, Niv Haim, Matan Atzmon, and Yaron Lipman. 2020. Im-
J Andreas Bærentzen. 2005. Robust generation of signed distance �elds from triangle plicit Geometric Regularization for Learning Shapes. CoRR abs/2002.10099 (2020).
meshes. In Fourth International Workshop on Volume Graphics, 2005. IEEE, 167–239. arXiv:2002.10099 [Link]
Josh Barnes and Piet Hut. 1986. A hierarchical O (N log N) force-calculation algorithm. Gaël Guennebaud, Benoît Jacob, et al. 2010. Eigen v3. [Link]
nature 324, 6096 (1986), 446–449. Eric Haines. 1994. Point in Polygon Strategies. Graphics Gems 4, 1 (1994), 24–46.
Alexander Belyaev and Pierre-Alain Fayolle. 2020. An ADMM-based scheme for distance Si Hang. 2015. TetGen, a Delaunay-based quality tetrahedral mesh generator. ACM
function approximation. Numerical Algorithms 84 (2020), 14. [Link] Trans. Math. Softw 41, 2 (2015), 11.
1007/s11075-019-00789-5 John C Hart. 1996. Sphere tracing: A geometric method for the antialiased ray tracing
Alexander Belyaev, Pierre-Alain Fayolle, and Alexander Pasko. 2013. Signed Lp-distance of implicit surfaces. The Visual Computer 12, 10 (1996), 527–545.
�elds. Computer-Aided Design 45, 2 (2013), 523–528. Yixin Hu, Qingnan Zhou, Xifeng Gao, Alec Jacobson, Denis Zorin, and Daniele Panozzo.
Nicole Berline, Ezra Getzler, and Michèle Vergne. 1992. Heat Kernels and Dirac Operators 2018. Tetrahedral meshing in the wild. ACM Trans. Graph. 37, 4 (2018), 60–1.
(1 ed.). Springer Berlin, Heidelberg. Saar Huberman, Amit Bracha, and Ron Kimmel. 2023. Deep Accurate Solver for the
David Bommes and Leif Kobbelt. 2007. Accurate Computation of Geodesic Distance Geodesic Problem. In International Conference on Scale Space and Variational Methods
Fields for Polygonal Curves on Triangle Meshes. VMV, 151–160. in Computer Vision. Springer, 288–300.
Alan Brunton and Lubna Abu Rmaileh. 2021. Displaced Signed Distance Fields for Alec Jacobson, Ladislav Kavan, and Olga Sorkine. 2013. Robust Inside-Outside Segmen-
Additive Manufacturing. ACM Trans. Graph. 40, 4, Article 179 (July 2021), 13 pages. tation using Generalized Winding Numbers. ACM Trans. Graph. 32, 4 (2013).
[Link] Alec Jacobson and Daniele Panozzo. 2017. Libigl: Prototyping geometry processing
Astrid Bunge, Philipp Herholz, Misha Kazhdan, and Mario Botsch. 2020. Polygon Lapla- research in c++. In SIGGRAPH Asia 2017 courses. 1–172.
cian Made Simple. Computer Graphics Forum 39, 2 (2020), 303–313. [Link] Alec Jacobson, Daniele Panozzo, et al. 2018. libigl: A simple C++ geometry processing
10.1111/cgf.13931 arXiv:[Link] library. [Link]
Fatih Calakli and Gabriel Taubin. 2011. SSD: Smooth signed distance surface recon- Michael Kazhdan, Matthew Bolitho, and Hugues Hoppe. 2006. Poisson Surface Re-
struction. In Computer Graphics Forum, Vol. 30. Wiley Online Library, 1993–2002. construction. In Proceedings of the Fourth Eurographics Symposium on Geometry
Antonin Chambolle and Thomas Pock. 2011. A First-Order Primal-Dual Algorithm for Processing (Cagliari, Sardinia, Italy) (SGP ’06). Eurographics Association, Goslar,
Convex Problems with Applications to Imaging. Journal of Mathematical Imaging DEU, 61–70.
and Vision 40 (2011). [Link] Michael Kazhdan and Hugues Hoppe. 2013. Screened Poisson Surface Reconstruction.
Sungjoon Choi, Qian-Yi Zhou, Stephen Miller, and Vladlen Koltun. 2016. A Large ACM Trans. Graph. 32, 3, Article 29 (jul 2013), 13 pages. [Link]
Dataset of Object Scans. arXiv:1602.02481 (2016). 2487228.2487237
David Coeurjolly and Jacques-Olivier Lachaud. 2022. A Simple Discrete Calculus for Ron Kimmel and James A Sethian. 1998. Computing geodesic paths on manifolds.
Digital Surfaces. In IAPR Second International Conference on Discrete Geometry and Proceedings of the national academy of Sciences 95, 15 (1998), 8431–8435.
Mathematical Morphology, Étienne Baudrier, Benoît Naegel, Adrien Krähenbühl, Bruno Lévy and Hao Zhang. 2010. Spectral mesh processing. In ACM SIGGRAPH 2010
and Mohamed Tajine (Eds.). Springer, LNCS. Courses. 1–312.
David Coeurjolly, Jacques-Olivier Lachaud, et al. 2010. DGtal: Digital Geometry tools Moshe Lichtenstein, Gautam Pai, and Ron Kimmel. 2019. Deep eikonal solvers. In Scale
and algorithms library. [Link] Space and Variational Methods in Computer Vision: 7th International Conference, SSVM
David Coeurjolly, Jacques-Olivier Lachaud, and Jérémy Levallois. 2014. Multigrid 2019, Hofgeismar, Germany, June 30–July 4, 2019, Proceedings 7. Springer, 38–50.
convergent principal curvature estimators in digital geometry. Computer Vision and Roee Litman and Alex M Bronstein. 2016. Spectrometer: Amortized sublinear spectral
Image Understanding 129 (2014), 27–41. approximation of distance on graphs. In 2016 Fourth International Conference on 3D
David Coeurjolly and Jérémy Levallois. 2015. VolGallery. [Link] Vision (3DV). IEEE, 499–508.
VolGallery. Charles Loop. 1987. Smooth Subdivision Surfaces Based on Triangles. Ph.D. Dissertation.
Keenan Crane, Fernando De Goes, Mathieu Desbrun, and Peter Schröder. 2013a. Digital William E Lorensen and Harvey E Cline. 1998. Marching cubes: A high resolution 3D
geometry processing with discrete exterior calculus. In ACM SIGGRAPH 2013 surface construction algorithm. In Seminal graphics: pioneering e�orts that shaped
Courses. 1–126. the �eld. 347–353.
Keenan Crane, Marco Livesu, Enrico Puppo, and Yipeng Qin. 2020. A Survey of Algo- Zoë Marschner, Silvia Sellán, Hsueh-Ti Derek Liu, and Alec Jacobson. 2023. Constructive
rithms for Geodesic Paths and Distances. arXiv e-prints, Article arXiv:2007.10430 Solid Geometry on Neural Signed Distance Fields. In SIGGRAPH Asia 2023 Conference
(July 2020), arXiv:2007.10430 pages. arXiv:[Link]/2007.10430 Papers. 1–12.
Keenan Crane, Clarisse Weischedel, and Max Wardetzky. 2013b. Geodesics in heat: Joseph S. B. Mitchell, D. Mount, and C. Papadimitriou. 1987. The Discrete Geodesic
A new approach to computing distance based on heat �ow. ACM Transactions on Problem. SIAM J. Comput. 16, 4 (1987), 647–668.
Graphics (TOG) 32, 5 (2013), 1–11. Patrick Mullen, Fernando De Goes, Mathieu Desbrun, David Cohen-Steiner,
George Dantzig. 1963. Linear programming and extensions. Princeton university press. and Pierre Alliez. 2010. Signing the Unsigned: Robust Surface Re-
Alexandre Djerbetian and Mirela Ben Chen. 2016. Tangent Vector Fields on Triangulated construction from Raw Pointsets. Computer Graphics Forum 29, 5
Surfaces-An Edge-Based Approach. Ph.D. Dissertation. Computer Science Department, (2010), 1733–1741. [Link]
Technion. arXiv:[Link]
Michal Edelstein, Nestor Guillen, Justin Solomon, and Mirela Ben-Chen. 2023. A Ken Museth, David E Breen, Ross T Whitaker, and Alan H Barr. 2002. Level set surface
Convex Optimization Framework for Regularized Geodesic Distances. In ACM editing operators. In Proceedings of the 29th annual conference on Computer graphics
SIGGRAPH 2023 Conference Proceedings (Los Angeles, CA, USA) (SIGGRAPH ’23). and interactive techniques. 330–338.
Association for Computing Machinery, New York, NY, USA, Article 2, 11 pages. Ashish Myles, Nico Pietroni, and Denis Zorin. 2014. Robust �eld-aligned global
[Link] parametrization. ACM Trans. Graph. 33, 4, Article 135 (jul 2014), 14 pages. https:
Je� Erickson. 2019. Algorithms. Je� Erickson. [Link] //[Link]/10.1145/2601097.2601154
algorithms/ Giacomo Nazzaro, Enrico Puppo, and Fabio Pellacini. 2022. geoTangle: Interactive
Alexandre Ern and Jean-Luc Guermond. 2004. Theory and Practice of Finite Elements Design of Geodesic Tangle Patterns on Surfaces. ACM Transactions on Graphics 41,
(1 ed.). Applied Mathematical Sciences, Vol. 159. Springer, New York, NY. https: 2 (2022), 12:1–12:17. [Link]
//[Link]/10.1007/978-1-4757-4355-5 Helen Oleynikova, Alexander Millane, Zachary Taylor, Enric Galceran, Juan Nieto, and
Nicole Feng, Mark Gillespie, and Keenan Crane. 2023. Winding Numbers on Discrete Roland Siegwart. 2016. Signed distance �elds: A natural representation for both
Surfaces. ACM Trans. Graph. 42, 4, Article 36 (jul 2023), 17 pages. [Link] mapping and planning. In RSS 2016 workshop: geometry and beyond-representations,

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
 A Heat Method for Generalized Signed Distance • 92:15

physics, and scene understanding for robotics. University of Michigan. geodesic circle ⇠Y of radius Y > 0 centered on a point source at vertex
Stanley Osher, Ronald Fedkiw, and K Piechor. 2004. Level set methods and dynamic 8, and take the limit as Y ! 0. We let `Y := (⇥8 Y) 1 H⇠1 be a measure
implicit surfaces. Appl. Mech. Rev. 57, 3 (2004), B15–B15. Y
Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Love- of unit mass supported on ⇠Y , where H⇠1 is the Hausdor� measure
grove. 2019. DeepSDF: Learning Continuous Signed Distance Functions for Shape Õ Y
on ⇠Y , and ⇥8 := 89:>8 \8 be the discrete angle sum around vertex
9:
Representation. In Proceedings of the IEEE/CVF Conference on Computer Vision and
Pattern Recognition (CVPR). 8. We let # denote the outward unit normals to ⇠Y . We integrate
Konrad Polthier and Eike Preuß. 2003. Identifying Vector Field Singularities Using a the vector-valued measure # `Y against CR basis functions for each
Discrete Hodge Decomposition. In Visualization and Mathematics III, Hans-Christian
Hege and Konrad Polthier (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, edge 4 2 ⇢ (Section 5.6).
113–134. We restrict our attention to a single triangle 89:
Inigo Quilez. 2008. Raymarching Signed Distance Fields. [Link]
raymarchingdf/.
containing 8 and 4; the �nal contribution to the
Sudhanshu Rawat and Mantosh Biswas. 2022a. Enhanced Heat Method for Computation 89 entry is the sum of contributions from all 89:
of Geodesic Distance on Triangular Meshes. In 2022 IEEE 9th Uttar Pradesh Section incident on 4. W.l.o.g. we let 4 = 89. We express
International Conference on Electrical, Electronics and Computer Engineering (UPCON).
IEEE, 1–5. all quantities in complex numbers with respect to
Sudhanshu Rawat and Mantosh Biswas. 2022b. Modi�ed Heat Method using GMRES the polar coordinate system with origin at vertex 8 and real axis
for Geodesic Distance. In 2022 IEEE Conference on Interdisciplinary Approaches in ô
along 89 , parameterizing ⇠Y by the angle \ between G 2 ⇠Y and 4 89ô .
Technology and Management for Social Innovation (IATMSI). IEEE, 1–5.
Rohan Sawhney, Ruihao Ye, Johann Korndoerfer, and Keenan Crane. 2020. FCPW: Expressed in this basis, the unit normal # (\ ) at G = Y4y\ is simply
4y\ . We let U := \8 , and evaluate
Fastest Closest Points in the West. 9:
Silvia Sellán and Alec Jacobson. 2022. Stochastic Poisson Surface Reconstruction. ACM π π U
Trans. Graph. 41, 6, Article 227 (nov 2022), 12 pages. [Link] 1 (1 4yU )y
3555441 lim h=,k89 i d`Y = lim # (\ )Y d\ .
Y!0 89: Y!0 ⇥8 Y 0 ⇥8
James A Sethian. 1999. Fast marching methods. SIAM review 41, 2 (1999), 199–235.
Nicholas Sharp and Keenan Crane. 2020a. A Laplacian for Nonmanifold Triangle The contribution to edge 89 multiplies this quantity by a sign B 89ô 2
Meshes. Computer Graphics Forum (SGP) 39, 5 (2020). ô
Nicholas Sharp and Keenan Crane. 2020b. You Can Find Geodesic Paths in Triangle {+1, 1} equal to 4y0 = +1 if 89 agrees with the global orientation
Meshes by Just Flipping Edges. ACM Trans. Graph. 39, 6 (2020). of 89, and 4 = 1 otherwise.
yc
Nicholas Sharp, Keenan Crane, et al. 2019a. GeometryCentral: A modern C++ library of
data structures and algorithms for geometry processing. [Link]
To compute the contributions to edge 9: (resp.
net/. (2019). :8), the only di�erence is that G 2 ⇠Y and =(\ )
Nicholas Sharp, Mark Gillespie, and Keenan Crane. 2021. Geometry Processing with must be expressed relative to the tangent basis at
Intrinsic Triangulations. (2021).
Nicholas Sharp, Yousuf Soliman, and Keenan Crane. 2019b. Navigating Intrinsic Trian- edge 9: (resp. :8) instead of 89. This amounts to
gulations. ACM Trans. Graph. 38, 4 (2019). rotating the coordinate system by A89!9: (Equa-
Nicholas Sharp, Yousuf Soliman, and Keenan Crane. 2019c. The Vector Heat Method. tion 8). We arrive at the per-face contributions
ACM Trans. Graph. 38, 3 (2019).
Oded Stein, Max Wardetzky, Alec Jacobson, and Eitan Grinspun. 2020. A Simple yU
Discretization of the Vector Dirichlet Energy. Computer Graphics Forum 39, 5 (2020).
(X0 )89 = B89 (1 ⇥48 )y
yU
[Link] (X0 ) 9: = B89 A89!9: (1 ⇥48 )y (19)
Kaiyue Sun and Xiangyang Liu. 2022. Heat method of non-uniform di�usion for com-
puting geodesic distance on images and surfaces. Multimedia Tools and Applications 1
(X0 ):8 = B89 A:8!89 (1 4yU )y
.
⇥8
81, 25 (2022), 36293–36308.
Jiong Tao, Juyong Zhang, Bailin Deng, Zheng Fang, Yue Peng, and Ying He. 2019. Parallel If we were to use basis functions at vertices, there would not be
and scalable heat methods for geodesic distance computation. IEEE transactions on
pattern analysis and machine intelligence 43, 2 (2019), 579–594. enough degrees of freedom to encode radially symmetric vector
Philip Trettner, David Bommes, and Leif Kobbelt. 2021. Geodesic distance computation �elds at vertices; Sharp et al. [2019c, App. A] ameliorate this short-
via virtual source propagation. In Computer graphics forum, Vol. 40. Wiley Online coming by arbitrarily taking the limit as Y ! 1 instead of zero, but
Library, 247–260.
Sathamangalam R Srinivasa Varadhan. 1967. On the behavior of the fundamental this limit is not scale-invariant like ours.
solution of the heat equation with variable coe�cients. Communications on Pure
and Applied Mathematics 20, 2 (1967), 431–455. A.2 Unoriented Curves
Delio Vicini, Sébastien Speierer, and Wenzel Jakob. 2022. Di�erentiable Signed Distance
Function Rendering. Transactions on Graphics (Proceedings of SIGGRAPH) 41, 4 (July Here we consider open curves to which we
2022), 125:1–125:18. [Link] compute unsigned distance. Similar to Ap-
Hongyi Xu and Jernej Barbič. 2014. Signed Distance Fields for Polygon Soup Meshes.
Graphics Interface 2014 (2014). pendix A.1, we consider a geodesic Y-o�set Y
Lior Yariv, Omri Puny, Natalia Neverova, Oran Gafni, and Yaron Lipman. 2023. Mosaic- of the (piecewise linear) curve , with a mea-
SDF for 3D Generative Models. arXiv preprint arXiv:2312.09222 (2023).
Qingnan Zhou, Eitan Grinspun, Denis Zorin, and Alec Jacobson. 2016. Mesh arrange-
sure of unit density concentrated on Y , and
ments for solid geometry. ACM Transactions on Graphics (TOG) 35, 4 (2016), 1–15. take Y ! 0. For simplicity, we consider unori-
ented curves that lie along edges of the mesh.
Then Y can generically be decomposed into
A ADDITIONAL SOURCE GEOMETRY
four types of curves: ( ) linear segments that
Here we derive discretizations for isolated point sources and unori- intersect faces incident on edges of the curve
ented curves. ; (⌫) segments that intersect faces incident
on interior vertices of ; (⇠) circular arcs that
A.1 Point Sources intersect faces incident on endpoints of ; and (⇡) linear segments
To compute unsigned distance to point sources, we encode a radially that intersect faces incident on endpoints of (inset).
symmetric vector �eld centered at each point source. We adopt As Y ! 0, the contributions of type-⌫ and ⇡ segments go to zero.
the approach of Sharp et al. [2019c, App. A] and consider a small A type- segment lying within face 89:, where edge 89 lies on the

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.
92:16 • Feng and Crane

curve , converges to an oriented curve along edge 89 as Y ! 0; For both = = 2, 3, the =-dimensional CR Laplacian has entries = 2
its contributions to each edge within 89: are given by Equation 11. times the corresponding entries of the so-called primal vertex Lapla-
For type-⇠ circular arcs, we obtain the same per-face formulas as cian Alexa et al. [2020], since CR basis functions are equal to the
Equation 19, except we normalize not by the total angle sum around hat functions on medial =-simplices.
the endpoint vertex but by the sum of corner angles of faces that
intersect the circular portion of Y , and weight by the length of the B.2 Divergence Operator
adjacent edge. As in 2D, we average vectors X from triangles to tetrahedra via
1⇣ ⌘
B CROUZEIX-RAVIART IN 3D X89:; = -89: + - 98; + -8:; + - 9;: ,
4
Here we derive the Crouzeix-Raviart Laplace and mass matrices which corresponds to evaluating the CR-interpolated �eld at the
for tetrahedral meshes. In general, Crouzeix-Raviart basis func- barycenter. Letting Y89:; := X89:; /kX89:; k in each tet, the discrete
tions are piecewise-linear and associated with the barycenters of divergence is then
codimension-one simplices in an =-dimensional simplex. The basis ’
function \8 associated with the (= 1)-dimensional face f8 oppo- (r · Y)8 = | 9;: |= 9;: · Y89:; . (22)
site vertex 8 has value 1 on f8 , and value 1 = at 8. Within each 89:; >8
(= 1)-dimensional simplex incident on f8 , the function \8 is linear B.3 Mass Matrix
and is de�ned as \8 (G) = 1 =_8 (G) where _8 (G) is the value of the
hat function associated with vertex 8 at point G [Ern and Guermond The | | ⇥ | | Crouzeix-Raviart mass matrix M has nonzero entries
given by
2004, §1.2]. We use i89: to denote the basis function associated with π
triangle face 89: in a tet mesh. M89:,98; = hi89: (G), i 98; (G)i dG
π 1 π 1 Dπ 1 C D
"
B.1 Scalar Laplacian ’
= 6|89:; | (1 3_; (G)) (1 3_: (G)) dB dC dD
The | | ⇥ | | Crouzeix-Raviart Laplace matrix L of a tet mesh " = 0 0 0
89:; >89:,98;
(+ , ⇢, ,) ) has entries given by
π where G is parameterized using barycentric coordinates of tet 89:;
L 5 ,5 0 = hri 5 (G), ri 5 0 (G)i dG with vertex positions E8 , E 9 , E: , E; , as G = BE8 +CE 9 +DE: +(1 B C D)E; .
" We obtain entries of the symmetric matrix as
which will be nonzero only if faces 5 , 5 0 are adjacent, so w.l.o.g. we 1 |89:; |,
M89:,98; = M89:,8:; = M89:,9;: = 20
Õ
889:; 2 )
let 5 := 89:, 5 0 := 98; and compute 2
π M89:,89: = 5 89:; >89: |89:; |, 889: 2 .
L89:,98; = hri89: (G), ri 98; (G)i dG (23)
π"
= hr(1 3_; (G)), r(1 3_: (G))i dG
"
’ π
= 9 hr_; (G), r_: (G)i dG .
89:; >89:,98; 89:;
89: 89: 89:
The gradient r_; = =; /⌘; , where =; de-
notes the unit normal to face 89: opposite
89:
vertex ;, and ⌘; the height of the tetrahe-
dron with apex ; and base 89: (and similarly
for r_: ). Since the gradients r_; , r_: are
constant per tet, the inner integral is equal
to
1 |89: | | 98; |
|89:; | cos(c \89:; ) = |89:; | cos \89:; (20)
89: 98;
⌘ ⌘ 3|89:; | 3|89:; |
; :
where \89:; denotes the dihedral angle at edge 89 opposite edge :;.
Since the volume of the tetrahedron 89:; can be expressed |89:; | =
2
3|89 | |89: || 98; | sin \ 89 , we obtain
:;

L89:,98; = 3 |89 | cot \ :; , 889: 2 , 89:; >89:

2 89
L89:,8:; = 3 |:8 | cot \ 9; , 889: 2 , 89:; >89:
2 :8
L89:,9;: = 3 | 9: | cot \ 8; , 889: 2 , 89:; >89:
Õ
2 9:
L89:,89: = 89:; >89: L89:,98; + L89:,8:; + L89:,9;: , 889: 2
(21)

ACM Trans. Graph., Vol. 43, No. 4, Article 92. Publication date: July 2024.