Surface-Filling Curve Flows via Implicit Medial Axes
YUTA NOMA, The University of Tokyo, Japan and University of Toronto, Canada
SILVIA SELLÁN, University of Toronto, Canada
NICHOLAS SHARP, NVIDIA, USA
KARAN SINGH, University of Toronto, Canada
ALEC JACOBSON, University of Toronto and Adobe Research, Canada

Fig. 1. Snakes on a Plane. Given an intricate surface with complex topology, our algorithm generates a curve (63,000 nodes) that completely covers the surface
while avoiding intersections in under one minute.
We introduce a fast, robust, and user-controllable algorithm to generatesurface-filling curves. We compute these curves through the gradient flowof a simple sparse energy, making our method several orders of magnitude
faster than previous works. Our algorithm makes minimal assumptions
on the topology and resolution of the input surface, achieving improved
robustness. Our framework provides tuneable parameters that guide theshape of the output curve, making it ideal for interactive design applications.
CCS Concepts: • Computing methodologies → Mesh models; Meshgeometry models; • Human-centered computing → Interactive systems
and tools.
Additional Key Words and Phrases: geometric flow, space-filling curve, geometry processing

Authors’ addresses: Yuta Noma, noma@akg.t.u-tokyo.ac.jp, The University of Tokyo,
Tokyo, Japan and University of Toronto, Toronto, Canada; Silvia Sellán, sgsellan@cs.
toronto.edu, University of Toronto, Toronto, Canada; Nicholas Sharp, nmwsharp@
gmail.com, NVIDIA, USA; Karan Singh, karan@dgp.toronto.edu, University of Toronto,
Toronto, Canada; Alec Jacobson, jacobson@cs.toronto.edu, University of Toronto and
Adobe Research, Toronto, Canada.

Permission to make digital or hard copies of part or all 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 third-party components of this work must be honored.
For all other uses, contact the owner/author(s).
© 2024 Copyright held by the owner/author(s).
0730-0301/2024/7-ART
https://doi.org/10.1145/3658158

ACM Reference Format:
Yuta Noma, Silvia Sellán, Nicholas Sharp, Karan Singh, and Alec Jacobson.
2024. Surface-Filling Curve Flows via Implicit Medial Axes. ACM Trans.
Graph. 43, 4 (July 2024), 12 pages. https://doi.org/10.1145/3658158

1

INTRODUCTION

A surface-filling curve is a continuous map from a one-dimensional
space onto a two-dimensional manifold, which guarantees that every point on the surface lies within a given maximum distance ofthe curve (see Figure 1). In computer graphics, these curves haveprimarily been leveraged for digital fabrication and artistic creation;yet, their applicability extends to a variety of fields, including robotic
path planning, sensor placement, circuit design, and more. Nonetheless, despite their significant potential, existing methods suffer fromslow computation time, robustness issues, and limited user control.

input

Fig. 2. We compute surface-filling curves through a geometric flow thatevolves to balance total length with distance to the curve’s medial axis.

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

2 •Noma et al.

[Sharp and Crane 2018]#V = 10K#V = 642#V = 42

Ours

Fig. 3. Surface-filling curves can potentially be used in applications likerobotic path planning for tasks including painting or polishing surfaces.

We introduce a fast, robust, and controllable algorithm to generatesurface-filling curves. We model the problem through a two-term
energy that balances curve length and distance to the curve’s medial axis. Instead of computing the medial axis explicitly, we only
implicitly compute the distance to it, making it ideal for speed androbustness. The gradient of this energy yields a geometric flow thatevolves any input curve into a surface-filling one (see Figure 2),
which is efficiently evaluated by solving a sparse linear system withsize proportional to the number of nodes in the curve.
Recently, the Repulsive Curves algorithm [Yu et al. 2021] showed
applications to this problem, generating curves constrained to simple surfaces with a physically-inspired energy. While that method
must deal with a dense linear system solve at each iteration, ourmethod produces qualitatively similar results via sparse ones, leading to a speed-up of several orders of magnitude. Combined withthis massive performance improvement, our algorithm’s superior
robustness enables curve generation on a broader set of surfaces
with diverse styles, resolutions, and topologies. Together with its
speed and robustness, our high degree of user control makes it ideal
for interactive applications in which a designer may wish to change
the spacing of the curve, align it to a given density map or vector
field, or even edit the underlying surface.
We showcase the performance of our method by testing it on adiverse set of inputs inspired by a broad range of contexts. We extensively evaluate our algorithmic choices through ablation studies
and include experiments showing the effect of the different parameters in our energy formulation. Our flow’s significant improved
performance enables a whole new set of applications, which we
exemplify through results in which a user can interactively specify
the curve’s global shape, local density, and alignment.

2

RELATED WORK

Surface-filling curves are a special case of space-filling curves that filla given 𝑑-dimensional domain with one continuous stroke without
self-intersections (e.g., [Hilbert 1891; Peano 1890]). Over the past
decades, space-filling curves have received renewed interest fromcomputer graphics researchers due to their applications ranging
from infill patterns for 3D printing [Chermain et al. 2023] to toolpath
planning for CNC milling [Zhao et al. 2018] and procedural art[Pedersen and Singh 2006]. To address these applications, a key
challenge lies in filling out intricate 2D shapes or 3D surfaces insteadof simpler geometries (e.g., spheres [Gerlach and von der Mosel
ACM Trans. Graph., Vol. 43, No. 4, Article . Publication date: July 2024.

finemesh resolution

coarseFig. 4. Unlike the distortion-minimizing flow by Sharp and Crane [2018],
our algorithm is robust to changes in underlying mesh resolution.

2011]). We classify such methods into two subclasses: initializationbased methods and geometric flows. We also show a feature tableof methods that can compute surface-filling curves in Table 1.

2.1

Initialization-based methods for space-filling curves

Initialization-based methods begin by searching for an approximation near the optimal state, often followed by a local relaxation of thecurve to determine its final form. These methods can be categorized
as either discrete or continuous.
Discrete approaches begin by expressing the domain as a finite
set of points [Giannatsis et al. 2015], clusters [Yan and Mostofi
2016], grid cells [Cheng et al. 2020; Joshi et al. 2019] or more generalpolygonal elements [Akleman et al. 2013; Bedel et al. 2022; Linet al. 2019; Xing et al. 2012] with some notion of adjacency. This
representation is then treated as a graph on which a Hamiltonian
cycle is computed or approximated, often followed by some local,
geometric post-processing in the form of smoothing [Xing et al.
2012] or alignment to a given vector field [Bedel et al. 2022].
Among discrete approaches, the most relevant for our work is
the surface-filling curve generation method by Xing et al. [2012]
and the follow-up work by Akleman et al. [2013]. They computea Hamiltonian cycle on the dual graph of any input surface mesh,followed by optional subdivision and smoothing. By construction,these algorithms are heavily dependent on the input triangulationand provide little to no user control in terms of curve spacing,alignment, and overall shape.
Continuous approaches to space-filling curve generation directlymodel the curve as a smooth, geometric object on the input domain. A subclass of continuous methods extracts the iso-contours
of a continuous function, such as a distance transform [Zhao et al.
2016, 2018], periodic function [Fang et al. 2020; Tricard et al. 2021],
cylindrical parametrization [Schüller et al. 2018], or simply the zcoordinates [Zhong et al. 2023]. However, they must carefully stitch
disconnected components (see, e.g., [Chermain et al. 2023]) or manually segment the mesh into cylindrical pieces (see [Schüller et al.
2018]) to ensure a true single-stroke filling curve.
To escape the inherent problems in initialization-based methods,
we instead opt for growing curves from arbitrary initializations via
geometric flows.

Surface-Filling Curve Flows via Implicit Medial Axes • 3

Table 1. Previous methods fall short in one way or another.

[Xing et al. 2012][Akleman et al. 2013][Schüller et al. 2018][Zhao et al. 2018][Sharp and Crane 2018][Yu et al. 2021]
Ours

2.2

Custom penalties
✗
✗
✗
✗
✓
✓
✓

Low-res meshes
✗
✗
✓
✓
✗
✗
✓

Geometric flows for space-filling curves

Geometric flows evolve a geometry by following the variational
gradient of some energy functional [Brakke 1992]. The most basic
example is curve-shortening flows, in which the curve evolves tominimize its own length. Often combined with area-preservation
terms [Crane et al. 2013] or modified to avoid singularities [Kazhdan
et al. 2012], variants of this flow have become common in geometry
processing tasks like smoothing [Taubin 1995] and even morphological operations [Sellán et al. 2020].
To the best of our knowledge, only two prior works have suggested using geometric flows to compute surface-filling curves. First,
Sharp and Crane [2018] observed that their conformality-inducing
flow serendipitously produced surface-filling curves if the weighton the length penalization was low enough. While an impressive
offshoot result, this method suffers from robustness issues at lowmesh resolutions (see Figure 4) and the inability to handle open
curves.
More recently, Yu et al. [2021] proposed following the gradient
flow of a physics-inspired repulsive energy between every disconnected pair of points on the curve. The authors followed a rigorous
mathematical derivation and introduced a modification of gradientdescent combined with a multi-grid scheme to evolve the curve
according to the proposed energy. While this mathematical andcomputational machinery makes their algorithm tractable, by considering an energy defined between every pair of points, their energy is fundamentally dense. Although they provide a hierarchicalmultigrid solver to circumvent dense linear solves, the linear solves
associated with each level of the hierarchy still remain a bottleneck.

[Yu et al. 2021]

Ours

runtime per iteration [ms][Yu et al. 2021]
w/o multigrid[Yu et al. 2021]
w/ multigrid
Ours

number of nodes

Fig. 5. We qualitatively match the results produced by Yu et al. [2021] (topleft), but instead of a dense linear system solve, we only need a sparseone whose entries are banded near the diagonals (bottom left), providing
several digits of speed-up (right). See Appendix C for details on runtime
measurement.

Spacing control
✗
✗
✓
✓
✗
✗
✓

Unicursality
✓
✓
✓
✗
✗
✓
✓

Mesh-surgery-free
✓
✓
✗
✗
✗
✓
✓

In this paper, we build on the work by Yu et al. [2021] and showthat this repulsive energy can be substituted by a simple geometricproxy based on the distance from each point on the curve to its
medial axis. As we show in Figure 5 and empirically throughout the
paper, this substitution produces qualitatively similar results whilecritically converting the dense solve into a sparse one, obtaining a
massive speed-up and enabling a broader set of applications.

3

METHODIn what follows, let us assume that we are given a surface Ω ∈ R3 anda curve 𝛾 : [0, 𝑙] ↦→ Ω ∈ R3 on it. We define a consistent frame on 𝛾
as follows: at any point 𝛾 (𝑠) on 𝛾, 𝑇 (𝑠) is the (unit) tangent vector
to the curve, 𝑁 (𝑠) is the (orthogonal) normal vector to the surface
at 𝛾 (𝑠) and 𝐵(𝑠) is the (unit) bitangent vector, which is defined as
𝐵(𝑠) = 𝑇 (𝑠) × 𝑁 (𝑠). We will use the notation ∥𝑎 − 𝑏 ∥ Ω to denotethe geodesic distance between two points 𝑎, 𝑏 ∈ Ω, and ∥𝑎 − 𝑏 ∥ todenote the Euclidean distance between two points 𝑎, 𝑏 ∈ R3 .
Our goal will be to find another curve 𝛾 ★ that fills Ω with someradius 𝑟 ; that is, for any point 𝑝 ∈ Ω, there exists a point 𝑞 ∈𝛾 ★ such that ∥𝑝 − 𝑞∥ Ω ≤ 𝑟 . We will do this by constructing anenergy functional 𝐸 (𝛾) and iteratively minimizing it with respectto 𝛾 starting from some initial curve 𝛾 0 .

3.1

Energy construction

To motivate our choice of energy, we will start from the very definition of a surface-filling curve stated above: for any point 𝑝 ∈ Ω, there
must exist a point 𝑞 ∈ 𝛾 ★ such that ∥𝑝 − 𝑞∥ Ω ≤ 𝑟 . Let us now briefly
consider the minimum distance function 𝑑 (𝑝, 𝛾) = min𝑞 ∈𝛾 ∥𝑝 − 𝑞∥ Ω
as a function of 𝑝. We can then rephrase the surface-filling condition as: for any point 𝑝 ∈ Ω, its minimum distance to the curve isbelow 𝑟 . 𝑑 (𝑝, 𝛾) being a continuous function of 𝑝, it is sufficient toimpose this condition on the points 𝑝 ∈ Ω at which 𝑑 (𝑝, 𝛾) is locallymaximized, and they always lie on the medial axis of 𝛾 on Ω, which
we denote M𝛾,Ω . In other words, 𝛾 is a surface-filling curve if forevery point 𝑝 on the medial axis M𝛾,Ω , its distance to the curve isbelow 𝑟 .
We will now define a mappingfrom [0, 𝑙] to M𝛾,Ω , which will allow us to restate the surface-filling
condition in terms of the curve 𝛾
itself. The key observation here isthat we only need the distance tothe medial axis rather than its explicit graph. For every 𝑠 ∈ [0, 𝑙],
ACM Trans. Graph., Vol. 43, No. 4, Article . Publication date: July 2024.

4•

Noma et al.

iterate(1) energy assembly using
implicit medial axes

input(2) descent direction
computation(3) retraction &self-intersection avoidance(4) dynamic refinement

final result

Fig. 6. We obtain a surface-filling curve by iteratively growing an input curve through a geometric flow. At each iteration, we efficiently assemble a quadraticenergy, compute its retracted descent direction on the surface, and evolve the curve along it, dynamically refining it to maintain a homogeneous discretization.

we will grow a (geodesic) sphereof radius 𝑟 tangent at 𝛾 (𝑠) until it touches a different point on 𝛾 (seeinset). The center of this maximal sphere is, by definition, an elementof the medial axis 𝑚(𝑠). In fact, one can grow two such spheres, one
in the direction of 𝐵(𝑠) and one in the opposite direction −𝐵(𝑠), leading to two (not necessarily distinct) points 𝑚 + (𝑠), 𝑚 − (𝑠) ∈ M𝛾,Ω ,which we call the implicit medial axis. When one of the spheres
grows infinitely (e.g., a completely convex curve), we clamp the
distances to a parameter 𝑟𝑚𝑎𝑥 .
These mappings 𝑚 +, 𝑚 − : [0, 𝑙] ↦→ M𝛾,Ω are not necessarily
injective, but they are surjective: as a direct consequence of thedefinition of medial axis, every element in it is the image throughone (or both) of 𝑚 + or 𝑚 − of at least one value of 𝑠. Thus, we can
now rephrase the surface-filling condition as a condition on the
curve: 𝛾 is surface-filling if for every 𝑠 ∈ [0, 𝑙], 𝛾 (𝑠) − 𝑚 + (𝑠) Ω ≤ 𝑟and ∥𝛾 (𝑠) − 𝑚 − (𝑠)∥ Ω ≤ 𝑟 . This reformulation is critical, as it allows
us to define the curve’s medial axis energy as
∫ 

2
𝐸𝑀 =
𝛾 (𝑠) − 𝑚 + (𝑠) Ω + ∥𝛾 (𝑠) − 𝑚 − (𝑠)∥ 2Ω 𝑑𝑠.
(1)𝛾

Since it accumulates the distance values in an 𝐿 2 sense, we can
expect that progressively minimizing this energy will concentrate
the medial axis distances around an increasingly smaller value. Of
course, there is an infinite number of curves that satisfy the surfacefilling constraint as defined above. A reasonable choice between allof them is the shortest one, which we can enforce by adding a term𝐿Ω (𝛾) measuring the total (geodesic) length to the energy
𝐸 = 𝐿Ω (𝛾) + 𝛼𝐸𝑀 ,(2)where 𝛼 is a balancing coefficient, which we will later show to berelated to the desired curve width 𝑟 . Eq. (2) is our final analytical energy, which we will now discretize in the form of a sparse quadratic
form that one can efficiently assemble.

3.2

Discretization

piecewise
linearpiecewise
geodesic

Let us assume that Ω is given asa triangular mesh with vertices
V and faces F . Instead of existing representations (e.g., [Bischoff
et al. 2005]),we represent 𝛾 as apiecewise geodesic curve by storing it as a set of 𝑛 nodes 𝛾 =ACM Trans. Graph., Vol. 43, No. 4, Article . Publication date: July 2024.

#nodes=180K

Fig. 7. Our method is fast enough to obtain surface-filling curves withhundreds of thousands of nodes (180K) in a reasonable time (229 s).

{𝛾 1, 𝛾 2, . . . , 𝛾𝑛 } that are connected by geodesic segments on Ω (see inset). Beyond manifoldness, we make no assumption on the topology
or connectivity of 𝛾 (see Figure 10) or Ω (see Figure 9 or Figure 12).
In what follows, we refer to elements in the underlying surface asvertices and faces whereas elements in the curve as nodes and segments. Each node 𝛾𝑖 is equipped with the vertex-blended normal 𝑁𝑖 ,the normalized bisection tangent 𝑇𝑖 projected onto its tangent plane,
and their cross product 𝐵𝑖 . As noted by Sharp and Crane [2020],
these geodesic segments are intrinsic; i.e., they are independent of
how the surface patch containing the segment is embedded in R3 .
We will begin by making use of this property to approximate thecurve length via a quadratic form.
3.2.1 Curve length. Consider the geodesic segment between 𝛾𝑖 and𝛾𝑖+1 , which we denote 𝛾𝑖,𝑖+1 , with geodesic length 𝑙𝑖 . Let 𝑅𝑖 (𝛾) be the(unique) rigid transformation that maps 𝛾𝑖 to the origin 0 and aligns
the tangent vector 𝑇𝑖 with the 𝑥-axis and the normal vector 𝑁𝑖 withthe 𝑧-axis, and let 𝑄𝑖 (𝛾) be the rigid
transformation that maps 𝛾𝑖,𝑖+1 to thepoint (𝑙𝑖 , 0, 0), similarly aligning the tangent vector 𝑇𝑖+1 with the 𝑥-axis and thenormal vector 𝑁𝑖+1 with the 𝑧-axis. Inrigid
tuitively, these are the rigid transformatransformation
tions one would obtain from unrolling all
the triangles in the surface patch containing 𝛾𝑖,𝑖+1 onto the plane (see inset).
Trivially, the geodesic distance between 𝛾𝑖 and 𝛾𝑖+1 is the Euclidean distance between 𝑅𝑖 (𝛾)𝛾𝑖 andSurface-Filling Curve Flows via Implicit Medial Axes • 5

descent direction

retraction

Euclidean

curve

Euclidean

geodesic

geodesic

Fig. 8. The descent direction of Eq. (5) stays on the tangent plane of thesurface (right), allowing larger steps than merely evaluating the Euclidean
length between each node (left).

Fig. 11. The implicit medial axis in Euclidean space is affected by curves on
nearby regions of the surface without local connectivity.

3.2.2 Medial axis energy. Recall that the medial axis energy is defined as
∫ 

2
𝐸𝑀 =
𝛾 (𝑠) − 𝑚 + (𝑠) Ω + ∥𝛾 (𝑠) − 𝑚 − (𝑠)∥ 2Ω 𝑑𝑠 ,𝛾

where 𝑚 + (𝑠) and 𝑚 − (𝑠) are the implicit medial axis, i.e., the two

Fig. 9. Our method does not assume the topology or connectivity of theunderlying surface, which may be open and have boundaries.
input

result (front)

result (back)Fig. 10. Our method does not assume the topology or connectivity of the
input curve and preserves it throughout the flow.

𝑄𝑖 (𝛾)𝛾𝑖+1 on the XY plane:𝑙𝑖 = ∥𝛾𝑖 − 𝛾𝑖+1 ∥ Ω = ∥𝑅𝑖 (𝛾)𝛾𝑖 − 𝑄𝑖 (𝛾)𝛾𝑖+1 ∥ 2 .

(3)

By assuming only small changes in 𝛾, we can linearize the products
in Eq. (3) and multiply and divide by 𝑙𝑖 to write the length as aquadratic form in 𝛾𝑖 , 𝛾𝑖+1 :𝑙𝑖 ≈1∥𝑅𝑖 𝛾𝑖 − 𝑄𝑖 𝛾𝑖+1 ∥ 22 .
𝑙𝑖

(4)

Accumulating the length of every segment, we can now approximate
the curve length through a quadratic form on 𝜸 :𝐿Ω (𝜸 ) ≈

𝑛−1
∑︁

1
1
∥𝑅𝑖 𝛾𝑖 − 𝑄𝑖 𝛾𝑖+1 ∥ 22 = 𝜸 𝑇 A𝜸 + ℓ𝜸 ,
𝑙
2
𝑖
𝑖=1

(5)

where 𝜸 is the vector of stacked coordinates of the curve nodes, A isa sparse symmetric positive-semidefinite matrix of size 3𝑛 × 3𝑛, andℓ is a vector of size 3𝑛 that accounts for the translation component
in the rigid transformations. We observed that this approximation
allows larger steps than merely adding up the Euclidean distance
between 𝛾𝑖 and 𝛾𝑖+1 (see Figure 8).
Now that we have managed to approximate the curve length as aquadratic form, we can proceed to do the same for the second halfof Eq. (2), the medial axis energy.

points on the medial axis M𝛾,Ω that are closest to 𝛾 (𝑠).
The most immediate question one faces when attempting to discretize 𝐸𝑀 is how to compute 𝑚 ± for a given curve node 𝛾𝑖 . It is
known that computing the explicit medial axis (i.e., the cut locus)
on curved surfaces often becomes tedious (see, e.g., [Mancinelli et al.
2021] Section 7), so we opt for computing the implicit medial axis
for each node. One could loop over every 𝛾 𝑗 , computing the largest
geodesic sphere that is tangent to both 𝛾𝑖 and 𝛾 𝑗 to then choose theminimum radius among all of them. However, this approach would
lead to a quadratic complexity in the number of nodes, which isprohibitive for large curves we wish to model (see Figure 7).
Instead, we once build a closest
point query data structure on the
curve nodes (and the boundary verbinary searchtices of the underlying mesh if they
exist) and then compute 𝑚𝑖+ by carrying out a binary search on the radius of a sphere grown from 𝛾𝑖 in the
direction of 𝐵𝑖 until it lays within atolerance of a given node 𝛾 𝑗 and contains no other node (see inset). We then repeat the same process for𝑚 − , growing a sphere in the direction of −𝐵𝑖 instead.
The geodesic closest point query can be built upon any stateof-the-art geodesic computation algorithm [Crane et al. 2020], but
substituting geodesic spheres with Euclidean spheres can make thealgorithm even faster. This substitution is a common operation ingeometry processing (e.g., [Yuksel 2015]). The Euclidean closest
point query can be built upon a bounding volume hierarchy, leading
to a time complexity of 𝑂 (log 𝑛) per node and 𝑂 (𝑛 log 𝑛) in total.
We note that all of the quantitative measurements in this paper
are based on the Euclidean medial axis. The tradeoff is that 𝑚𝑖± is
not guaranteed to lie on the surface (see Figure 6 (1)), leading tointerference from nearby surfaces without local connectivity (seeFigures 11 or 12) or larger errors on high curvature regions (seeFigure 13). We note that the Euclidean medial axis is more suitablein some situations regardless of the error, e.g., path planning for3D printing where the filament has a constant width in Euclidean
space. We leave this choice to the user, depending on their desired
applications.
ACM Trans. Graph., Vol. 43, No. 4, Article . Publication date: July 2024.

6 •Noma et al.

coarsegeodesic

EuclideanFig. 12. Our method does not assume the topology or connectivity of theunderlying surface, which may be self-intersecting. We also observed thatthe selection of Euclidean/geodesic space for the implicit medial axis affects
the result near the intersecting regions.
input

approximation
error

converged

Fig. 14. Our method preserves the inside/outside of the curve throughout
the flow.

To discretize Eq. (1), we adopt the same idea as the curve length
to unroll the geodesic path between 𝛾𝑖 and 𝑚𝑖± onto the tangent
space of 𝛾𝑖 by:
𝛾𝑖 − 𝑚𝑖± Ω = 𝛾𝑖 − 𝛾𝑖 ± 𝑟𝑖± 𝐵𝑖 2 ,(6)where 𝑟𝑖+ and 𝑟𝑖− are the sphere radii obtained by the binary search.
By treating the unrolled implicit medial axis 𝑚𝑖′± = 𝛾𝑖 ± 𝑟𝑖± 𝐵𝑖
as momentarily fixed for an optimization iteration and summing
contributions from each curve vertex, we approximate the medial
axis energy as:
𝐸𝑀 (𝜸 ) ≈𝑛∑︁𝑙𝑖 −1 + 𝑙𝑖 

2Fig. 15. Our flow dynamically splits and collapses edges to maintain a userspecified edge-length ℎ. In practice, we find our flow to be robust to thisparameter choice as long as it is below 𝑟 /2.5.

3.3

Fig. 13. Our medial axis energy can be defined via Euclidean or geodesic
sphere radii, with the difference between these two being higher in highcurvature regions.

2𝛾𝑖 − 𝑚𝑖′+ 2 + 𝛾𝑖 − 𝑚𝑖′− 2 ,

(7)

𝑖=1

which is yet another quadratic form1
𝐸𝑀 (𝜸 ) ≈ 𝜸 𝑇 M𝜸 + n𝜸 ,2(8)where M is a symmetric matrix (see Appendix B for discussions on
this approximation).
By combining Eq. (5) and Eq. (8), we can now write the approximation of the energy from Eq. (2) as a quadratic form
1
1
𝐸 (𝜸 ) = 𝜸 𝑇 A𝜸 + ℓ𝜸 + 𝛼 𝜸 𝑇 M𝜸 + 𝛼n𝜸 = 𝜸 𝑇 H𝜸 + v𝜸 ,22where, notably, H is a sparse symmetric matrix.
ACM Trans. Graph., Vol. 43, No. 4, Article . Publication date: July 2024.

(9)

dense

Curve evolution

Given a discrete curve 𝜸 𝑘 at the 𝑘-th timestep of the flow, we cancompute its approximated discretized energy 𝐸 (𝜸 𝑘 ) by evaluating
Eq. (9). We can then compute the gradient of this energy with respectto 𝜸 𝑘 as g𝑘 = ∇𝐸 (𝜸 𝑘 ) = H𝑘 𝜸 𝑘 + v𝑘 , and its Hessian as ∇2 𝐸 (𝜸 𝑘 ) =H𝑘 . We can then use both these quantities to get the preconditioned
descent direction via a single sparse linear system solved𝑘 = −(H𝑘 ) −1 g𝑘 .

(10)

Of course, this coordinate-wise descent direction is not guaranteed
to keep the curve on the surface Ω, so we will need to instead compute the retraction (see, e.g.,[Boumal 2023]) of the descent direction
onto Ω. Given the descent direction 𝑑𝑖𝑘 for each node 𝛾𝑖𝑘 , we compute the retraction by finding the projection of 𝑑𝑖𝑘 onto the tangent
plane of Ω at 𝛾𝑖𝑘 , which we denote 𝑑ˆ𝑖 𝑘 , and explicitly tracing thegeodesic path from 𝛾𝑖𝑘 in the direction 𝑑ˆ𝑖 𝑘 of length ∥𝑑ˆ𝑖 𝑘 ∥. Finally,
we reconnect the piecewise geodesic segments using a geodesic
computation algorithm. In practice, we find that the magnitude ofthe descent is larger in the first iterations and then plateaus (seeFigure 16).
3.3.1 Self-intersection avoidance. To guarantee an intersection-free
output, we add a backtracking line search during our curve evolution.
After following the descent direction in its entirety, we loop overevery triangle and check that no two geodesic segments lying on
the same triangle intersect (this is easy enough to do, as geodesics
have the form of straight segments within each triangle). If they
do intersect, we multiply the descent direction by a coefficient 𝑐 ∈[0, 1), evolve the curve, and start the check again. We continue
this procedure until the curve never causes self-intersection. This
procedure was necessary for all of the models we tested.
3.3.2 Dynamic refinement. Our computation of the medial axis
energy only considers distances between pairs of nodes; as such,it relies heavily on the curve being finely discretized. Further, for
our curve evolution to find an optimal surface-filling disposition,
it needs to be provided with enough degrees of freedom during
the optimization. For these reasons, we follow the lead of otherflow-based approaches [Sellán et al. 2023, 2020] and use a dynamicrefinement scheme that adaptively subdivides the curve during theflow: given a user-defined value ℎ, we collapse edges with lengths
smaller than ℎ and subdivide those larger than 2ℎ. We discuss theeffect of ℎ in Section 4.2.

Surface-Filling Curve Flows via Implicit Medial Axes • 7

Table 2. Runtimes reported for several meshes. We kept iterating until theaverage length of the descent direction per node went below 3 × 10 −4𝑟 .

Model
Boot
Bunny
Fish
Hand
Spot
SpringerStuffedtoy

Vertices

Nodes (when converged)

Iterations

Runtime

19.9K
3.3K
10.8K
8.0K
11.5K
9.7K
6.7K

7.5K
11.2K
4.7K
10.2K
3.0K
16.1K
58.1K

41
24
26
25
25
27
32

8.8 s
3.6 s
2.3 s
2.5 s
1.5 s
4.8 s
17.3 saverage length of the descent direction

highest bin in the radius histogram

Fig. 17. We find an empirical relation between the energy weight 𝛼 andthe curve thickness (here computed as the highest bin in the medial axis
distance histogram) that is consistent across input surfaces and curves.

iteration

radiusradiusradiusradiusradiusFig. 16. As the flow progresses, the magnitude of each descent diminishes
and the distance to the medial axis concentrates around a single value.

4

EXPERIMENTS AND APPLICATIONS

Implementation details. We implemented our algorithm in C++
using GeometryCentral [Sharp et al. 2019] and libigl [Jacobson
et al. 2018] for common geometry processing subroutines, TinyAD[Schmidt et al. 2022] for computing the Hessian in Eq. (10), knn-cpp[Meyer 2019] for closest point queries, and the MMP algorithm[Mitchell et al. 1987] to reconnect the geodesic segments. In Table 2,
we report the runtimes of several representative models in thispaper, which we have confirmed scales with O (𝑛 1+𝑓 ), where 𝑓 > 0
accounts for the Laplacian-like sparse system solve. All our reported
runtimes were carried out on a 2022 M2 MacBook Air with 24GB of
RAM.

4.1

Comparisons

We compare the performance of our algorithm to both the first work
to suggest computing surface-filling curves through geometric flows[Sharp and Crane 2018] (Figure 4) as well as to the most recent stateof-the-art algorithm for this task [Yu et al. 2021] (Figure 5).
As we show in Figure 4, our algorithm presents a robustness andmesh independence that is not present in the work by Sharp andCrane [2018], allowing its use in a broader set of surfaces. Further,as seen in Figure 5, by substituting the dense linear solve proposedby Yu et al. [2021] for a sparse one, our algorithm provides massive
speed-ups in several orders of magnitude, enabling us to compute intricate surface-filling curves with hundreds of thousands of vertices(see, e.g., Figures 1 and 7).
This performance improvement comes at a theoretical cost in theform of the approximations carried out in Section 3.2.2; however, as

radiusradiusradiusFig. 18. Given the observed relation in Figure 17, a user can specify thedesired surface-filling thickness through the parameter 𝛼.

we exhaustively show in the rest of this section, these approximations do not seem to have a significant impact on the quality of thecurves produced by our algorithm.

4.2

Experiments and ablations

As we illustrate in Figure 6, our algorithm works by iteratively
following the surface-retracted gradient of the energy in Eq. (2).
In general, as shown qualitatively in Figure 2 and numerically in
Figure 16, we find that the first iterations of the flow tend to produce
the largest changes in the curve followed by a gradual convergence
to a surface-filling curve whose medial axis distances are clusteredaround a value 𝑟 .
A key parameter in this energy is the balancing coefficient 𝛼,
which we find serves as a proxy for the thickness 𝑟 of the final√
surface-filling curve. Experimentally, we find the relation 𝑟 ≈ 1.68/ 𝛼,which appears to be consistent across different curve shapes, thicknesses, and underlying meshes (Figure 18). This matches the dimensional units of Eq. 2, where 𝛼 is presumed to have a dimensionof inverse squared length. While a theoretical derivation of thisrelation escapes us, as shown in Figure 17, one can still utilize it tospecify the width of the curve by simply changing 𝛼.
Our algorithm also requires the specification of two other parameters; namely, the target segment length ℎ and the maximum
sphere size 𝑟𝑚𝑎𝑥 . As shown in Figure 15, we find that our algorithm
is invariant to the selection of ℎ as long as ℎ < 𝑟 /2.5. Similarly, as
shown in Figure 19, we find that 𝑟𝑚𝑎𝑥 can be set to any value larger
than 10𝑟 without too significantly affecting the shape of the final
curve, with values between 2𝑟 and 10𝑟 producing curves with lesspronounced wriggles.
ACM Trans. Graph., Vol. 43, No. 4, Article . Publication date: July 2024.

8 •Noma et al.

energyinput

iteration(linear scale)Fig. 19. 𝑟𝑚𝑎𝑥 controls the numerical behavior of our energy by clamping
the medial axis distance, critical at early flow stages. This parameter can be
set to any value larger than 10𝑟 without much impact on the output, withvalues between 2𝑟 and 10𝑟 producing curves with less pronounced wriggles.

Fig. 21. By adding a field-alignment energy, we can encourage the curve toalign itself to a specified tangent vector field.

w/o bending energy

w/ bending energy

deforming curves computed on frame 0

Fig. 22. By adding a biharmonic bending energy, we can encourage thecurve to be straighter.
with running flow (4 fps)
mapped toGaussian curvature

Fig. 20. Our algorithm is fast enough to run the flow for one step per frame
while allowing real-time deformation of the underlying mesh. Even with asingle step per frame, the spacing remains constant as opposed to merely
deforming the curve.

4.3

mapped to
texture intensity

Fig. 23. An artist may choose to use a spatially-varying 𝛼 to generate curvesof varying density, by mapping it to a property like curvature for procedural
generation (left) or to any texture for more manual control (right).

User control

Our method’s robustness makes it ideal for interactive applications
where a user or designer may want to explore different curves,
surfaces, or parameters. This application guides our algorithmic
design even from a fundamental theoretical level; we make only
minimal (manifoldness) assumptions on the topology of the input
surface and curve (see Figures 9, 12, and 10). Also, as shown in Figure
20, its performance can accommodate even realtime deformation of
the underlying surface while running the flow.
For its use in artistic applications, a designer must be able tomodify the generated curves in an easily controlled manner. Our
framework allows for several modifications, such as aligning the
curve to a given vector field (see Figure 21) or eliminating pronounced wriggles (see Figures 22) by adding simple sparse energies(see Appendix A for details).
Our algorithm also provides a powerful yet intuitive way to control the shape of the curve by specifying a spatially varying density
through the 𝛼 parameter. As shown in Figure 23, an artist may procedurally generate aesthetic patterns by mapping 𝛼 to geometric
ACM Trans. Graph., Vol. 43, No. 4, Article . Publication date: July 2024.

quantities like the Gaussian curvature of the underlying surface
or, for more manual control, by simply painting the desired density on the surface through a texture. Also, as shown in Figure 25,
we can obtain a curve with a minimum number of sharp corners by
progressively decreasing 𝑟 combined with the bending energy.
Finally, as shown in Figure 26,
our method can obtain a unicursally connected curve, thanks to
its topology-preserving feature(see Figure 14). This can be favorable for deploying a robotic toolFig. 24. Surface-filling curves
path (Figure 3), interactive sensor
could be used for fashion design,
network (Figure 24), or heating/- such as an interactive body sensor
cooling system (Figure 24 and Fig- network or sportswear that cools
ure 27) on intricate surfaces.
down body temperature.

Surface-Filling Curve Flows via Implicit Medial Axes • 9

Fig. 25. By progressively decreasing 𝑟 combined with the bending term, we
can produce a surface-filling curve with a small number of sharp corners.

situations where the underlying geometry has noisy normal directions (see inset). This is because our approach depends significantly on the normal direction, and we currently lack an alternative
method to bypass this issue.
Finally, we hope that by releasing ourfast, robust, and user-controllable algorithm to the geometry processing community and beyond, we can enable a novel array ofpreviously impracticable applications of surface-filling curves forboth interactive artistic design and industrial applications at scale.

ACKNOWLEDGMENTS

input

final

Fig. 26. Given two user-specified endpoints as an input curve, our method
can compute a unicursally connected surface-filling path.

Fig. 27. One potential application of our method is a cooling or heating
system deployed on complex surfaces, such as an artificial satellite, a spaceship, or a spacesuit.

5

DISCUSSION AND CONCLUSIONS

We have proposed a geometric flow that evolves a given curve on
an input surface into a surface-filling curve with a user-specified
thickness 𝑟 . We have achieved this by defining a medial axis energy,
which we have discretized through several approximations into aneasily assemblable sparse quadratic energy whose descent direction
we follow via a Newton scheme.
Our proposed algorithm achieves results qualitatively similar to
the most recent work Yu et al. [2021] at a massive speed-up. However, it should be noted that surface-filling curves are only one
example of the more general repulsive energy defined by Yu et al.
[2021], for most of which said work remains state of the art. Specifically, a promising avenue for future work involves generalizing ourapproach to space-filling curves contained on volumetric domains,or even codimension-one space-filling surfaces.
Additionally, a fundamental theoretical inquiry persists regarding
the empirical relationship in Figure 17 between curve thickness 𝑟
and the energy weight 𝛼.
Also, while our method showed excellent robustness on a diverse
set of input surfaces, we observed failure in very few challenging

We thank Chris Yu for early discussions, Abhishek Madan, Chenxi
Liu, Kinjal Parikh, Victor Rong, and Zhecheng Wang for help in
proofreading and rendering figures, and Esther Lin for video narration.
This work was supported in part by JST ACT-X (JPMJAX200L),
JSPS KAKENHI (JP23KJ0582), NSERC Discovery (RGPIN–2022–04680),
the Ontario Early Research Award program, the Canada Research
Chairs Program, a Sloan Research Fellowship, an NSERC Vanier
graduate scholarship, an Adobe Research fellowship, the DSI Catalyst Grant program, and gifts by Adobe Inc.
We acknowledge the authors of the 3D models used throughout this paper and thank them for making them available for academic use. Figures in this work contain the springer [TimEdwards2023], plane [YahooJAPAN 2019], car [kotakotakota 2008], robot
arm [jobojobojobo 2021], Spot [Crane 2019b], boot [inciprocal 2023],
beetle [Sutherland 1972], tower [jansentee3d 2021], Inopos [CosmoWenman 2013], ouroboros [FarnkFella 2012], jellyfish [YahooJAPAN 2013a], stuffed toy [Baran and blobs 2020], cow [3dcs_ch2023], hand [Max Planck Society e.V via SMPL 2021], bunny [TheStanford 3D Scanning Repository 2021], fish [Crane 2019a], Max
Planck [MPI 2020], scorpion [YahooJAPAN 2021], fandisk [Hoppe
2020], Brucewick [Holinaty 2023], falcon statue [colinfizgig 2021],
Demosthenes [ryanfb 2021], Bob [Crane [n. d.]], snail [YahooJAPAN
2013b], armadillo [The Stanford 3D Scanning Repository 2023],
avatar [Adobe [n. d.]], artificial star [tato_713 2023], astronaut [MaxGrueter 2014], fertility [AIM@Shape [n. d.]], Space Shuttle [Microsoft 2014], and octopus [DiabaseEngineer 2014] meshes.

REFERENCES
3dcs_ch. 2023. Cow. Downloaded modified version from odedstein-meshes github.com/
odedstein/meshes/tree/master/objects/cow, originally from thingiverse.com/thing:2431906. Asset licensed under CC BY 4.0..
Adobe. [n. d.]. Mixamo. mixamo.com
AIM@Shape. [n. d.]. Fertility.
Ergun Akleman, Qing Xing, Pradeep Garigipati, Gabriel Taubin, Jianer Chen, and Shiyu
Hu. 2013. Hamiltonian cycle art: Surface covering wire sculptures and duotone
surfaces. Computers & graphics 37, 5 (2013), 316–332.
E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A.
Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. 1999. LAPACK Users’
Guide (third ed.). Society for Industrial and Applied Mathematics, Philadelphia, PA.
Ilya Baran and Cosmic blobs. 2020. Stuffed toy.
Downloaded modified version from odedstein-meshes github.com/odedstein/meshes/tree/master/objects/
stuffedtoy, originally from mit.edu/~ibaran/autorig..
A. Bedel, Y. Coudert-Osmont, J. Martínez, R. I. Nishat, S. Whitesides, and S. Lefebvre. 2022. Closed space-filling curves with controlled orientation for 3D printing.
Computer Graphics Forum 41, 2 (2022), 473–492. https://doi.org/10.1111/cgf.14488
arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.14488

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

10•

Noma et al.

Stephan Bischoff, Tobias Weyand, and Leif Kobbelt. 2005. Snakes on Triangle Meshes. In
Bildverarbeitung für die Medizin 2005, Hans-Peter Meinzer, Heinz Handels, Alexander
Horsch, and Thomas Tolxdorff (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg,
208–212.
Nicolas Boumal. 2023. An introduction to optimization on smooth manifolds. Cambridge
University Press. https://doi.org/10.1017/9781009166164
Kenneth A Brakke. 1992. The surface evolver. Experimental mathematics 1, 2 (1992),
141–165.
Puikei Cheng, Wing Kam Liu, Kornel Ehmann, and Jian Cao. 2020. Enumeration of
additive manufacturing toolpaths using Hamiltonian paths. Manufacturing Letters
26 (2020), 29–32.
Xavier Chermain, Cédric Zanni, Jonàs Martínez, Pierre-Alexandre Hugron, and Sylvain Lefebvre. 2023. Orientable Dense Cyclic Infill for Anisotropic Appearance
Fabrication. ACM Trans. Graph. 42, 4, Article 68 (jul 2023), 13 pages. https://doi.org/10.1145/3592412
colinfizgig. 2021. Falcon statue. Downloaded modified version from odedsteinmeshes github.com/odedstein/meshes/tree/master/objects/falconstatue, originallyfrom thingiverse.com/thing:46631. Asset licensed under CC BY-NC-SA 4.0..
CosmoWenman. 2013. Inopos. Downloaded from thingiverse.com/thing:49168. Assetlicensed under CC BY 4.0..
Keenan Crane. [n. d.]. Bob.
Downloaded from cs.cmu.edu/~kmcrane/Projects/
ModelRepository..
Keenan Crane. 2019a. Fish. Downloaded modified version from odedstein-meshes
github.com/odedstein/meshes/tree/master/objects/fish, originally from cs.cmu.edu/~kmcrane/Projects/ModelRepository..
Keenan Crane. 2019b. Spot. Downloaded modified version from odedstein-meshes
github.com/odedstein/meshes/tree/master/objects/spot, originally from cs.cmu.edu/~kmcrane/Projects/ModelRepository..
Keenan Crane, Marco Livesu, Enrico Puppo, and Yipeng Qin. 2020. A Surveyof Algorithms for Geodesic Paths and Distances. CoRR abs/2007.10430 (2020).
arXiv:2007.10430 https://arxiv.org/abs/2007.10430
Keenan Crane, Ulrich Pinkall, and Peter Schröder. 2013. Robust fairing via conformal
curvature flow. ACM Transactions on Graphics (TOG) 32, 4 (2013), 1–10.
DiabaseEngineer. 2014. Flexible Octopus. Downloaded from thingiverse.com/thing:
159217. Asset licensed under CC BY 3.0..
Guoxin Fang, Tianyu Zhang, Sikai Zhong, Xiangjia Chen, Zichun Zhong, and Charlie CL Wang. 2020. Reinforced FDM: Multi-axis filament alignment with controlled
anisotropic strength. ACM Transactions on Graphics (TOG) 39, 6 (2020), 1–15.
FarnkFella. 2012. Auryn Medallion from The NeverEnding Story. Downloaded fromthingiverse.com/thing:49168. Asset licensed under CC BY 4.0..
Henryk Gerlach and Heiko von der Mosel. 2011. On Sphere-Filling Ropes. Am. Math.
Mon. 118, 10 (2011), 863–876. https://doi.org/10.4169/AMER.MATH.MONTHLY.118.
10.863
John Giannatsis, Aggelos Vassilakos, Vassilios Canellidis, and Vassilis Dedoussis. 2015.
Fabrication of graded structures by extrusion 3D Printing. In 2015 IEEE International
Conference on Industrial Engineering and Engineering Management (IEEM). IEEE,
175–179.
David Hilbert. 1891. Ueber die stetige Abbildung einer Line auf ein Flächenstück. Math.
Ann. (1891).
Josh Holinaty. 2023. Brucewick. Downloaded from odedstein-meshes github.com/
odedstein/meshes/tree/master/objects/brucewick. Asset licensed under CC BY-NC
4.0..
Hughes Hoppe. 2020. Fandisk. Downloaded modified version from common-3d-testmodels https://github.com/alecjacobson/common-3d-test-models..
inciprocal. 2023. Boot. Downloaded modified version from odedstein-meshes github.
com/odedstein/meshes/tree/master/objects/boot, originally from sketchfab.com/3dmodels/harrison-work-boot-d9506d9a318e4c56ac3d8ab273bd2540. Asset licensed
under CC BY 4.0..
Alec Jacobson, Daniele Panozzo, et al. 2018. libigl: A simple C++ geometry processing
library. https://libigl.github.io/.
jansentee3d. 2021. Tower. Downloaded modified version from odedstein-meshes github.
com/odedstein/meshes/tree/master/objects/tower, originally from thingiverse.com/
thing:3155868. Asset licensed under CC BY 3.0..
jobojobojobo. 2021. Industrial robot KUKA KR210 + Flange + Linear + Gripper+
Armature (Bones). Downloaded from blendswap.com/blend/29112. Asset licensed
under CC-0.
Anant A Joshi, Maulik C Bhatt, and Arpita Sinha. 2019. Modification of Hilbert’s
space-filling curve to avoid obstacles: a robotic path-planning strategy. In 2019 Sixth
Indian Control Conference (ICC). IEEE, 338–343.
Michael Kazhdan, Jake Solomon, and Mirela Ben-Chen. 2012. Can mean-curvature flow
be modified to be non-singular?. In Computer Graphics Forum, Vol. 31. Wiley Online
Library, 1745–1754.
kotakotakota. 2008. Body of Nissan 350Z Low-Poly. Downloaded from turbosquid.
com/3d-models/425765..
Sen Lin, Lingwei Xia, Guowei Ma, Shiwei Zhou, and Yi Min Xie. 2019. A maze-like
path generation scheme for fused deposition modeling. The International Journal of

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

Advanced Manufacturing Technology 104 (2019), 1509–1519.
C. Mancinelli, M. Livesu, and E. Puppo. 2021. Practical Computation of the Cut Locus on
Discrete Surfaces. Computer Graphics Forum 40, 5 (2021), 261–273. https://doi.org/
10.1111/cgf.14372 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.14372
Max Planck Society e.V via SMPL. 2021. Hand.
Downloaded modified version
from odedstein-meshes github.com/odedstein/meshes/tree/master/objects/hand,
adapted from M. Loper, N. Mahmood, J. Romero, G. Pons-Moll, M. Black ""SMPL""(SIGGRAPH Asia 2015). Asset licensed under CC BY 4.0..
MaxGrueter. 2014. Apollo Astronaut (The Original. Downloaded from thingiverse.
com/thing:230305. Asset licensed under CC BY 4.0..
Fabian Meyer. 2019. knn-cpp. https://github.com/Rookfighter/knn-cpp.
Microsoft. 2014. Space Shuttle. Downloaded from sketchfab.com/3d-models/spaceshuttle-ff4b00b7ebb24fdd98fb96b08f2c43c9. Asset licensed under CC BY 4.0..
Joseph S. B. Mitchell, David M. Mount, and Christos H. Papadimitriou. 1987. The
Discrete Geodesic Problem. SIAM J. Comput. 16, 4 (1987), 647–668. https://doi.org/
10.1137/0216045 arXiv:https://doi.org/10.1137/0216045
MPI. 2020. Max Planck. Downloaded modified version from common-3d-test-models
https://github.com/alecjacobson/common-3d-test-models, originally from https://gfx.cs.princeton.edu/proj/sugcon/models/..
G Peano. 1890. Sur une courde qui remplit touteune aire plane. Math. Ann. 36 (1890).
Hans Pedersen and Karan Singh. 2006. Organic Labyrinths and Mazes. In Proceedingsof the 4th International Symposium on Non-Photorealistic Animation and Rendering(Annecy, France) (NPAR ’06). Association for Computing Machinery, New York, NY,
USA, 79–86. https://doi.org/10.1145/1124728.1124742
ryanfb. 2021. Demosthenes.
Downloaded modified version from odedsteinmeshes
github.com/odedstein/meshes/tree/master/objects/demosthenes,
originally from sketchfab.com/3d-models/louvre-demosthenes-photoscancf5d9dd9bf4a457c95dc76e17848727b. Asset licensed under CC BY 4.0..
P. Schmidt, J. Born, D. Bommes, M. Campen, and L. Kobbelt. 2022. TinyAD:
Automatic Differentiation in Geometry Processing Made Simple. Computer
Graphics Forum 41, 5 (2022), 113–124.
https://doi.org/10.1111/cgf.14607
arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.14607
Christian Schüller, Roi Poranne, and Olga Sorkine-Hornung. 2018. Shape representation
by zippables. ACM Trans. Graph. 37, 4, Article 78 (jul 2018), 13 pages. https://doi.org/10.1145/3197517.3201347
Silvia Sellán, Christopher Batty, and Oded Stein. 2023. Reach For the Spheres: Tangencyaware surface reconstruction of SDFs. , Article 73 (2023), 11 pages. https://doi.org/
10.1145/3610548.3618196
Silvia Sellán, Jacob Kesten, Ang Yan Sheng, and Alec Jacobson. 2020. Opening andClosing Surfaces. ACM Trans. Graph. 39, 6, Article 198 (nov 2020), 13 pages. https://doi.org/10.1145/3414685.3417778
Nicholas Sharp and Keenan Crane. 2018. Variational Surface Cutting. ACM Trans. Graph.
37, 4, Article 156 (jul 2018), 13 pages. https://doi.org/10.1145/3197517.3201356
Nicholas Sharp and Keenan Crane. 2020. You Can Find Geodesic Paths in Triangle
Meshes by Just Flipping Edges. ACM Trans. Graph. 39, 6, Article 249 (nov 2020),
15 pages. https://doi.org/10.1145/3414685.3417839
Nicholas Sharp, Keenan Crane, et al. 2019. GeometryCentral: A modern C++ library ofdata structures and algorithms for geometry processing. https://geometry-central.
net/. (2019).
Ivan Sutherland. 1972. Sutherland’s Volkswagen.
tato_713. 2023. Death Star scaled one in two million. Downloaded from printables.
com/model/57266. Asset licensed under CC BY 4.0..
Gabriel Taubin. 1995. A signal processing approach to fair surface design. In Proceedingsof the 22nd annual conference on Computer graphics and interactive techniques. 351–
358.
The Stanford 3D Scanning Repository. 2021. Bunny. Downloaded modified version
from odedstein-meshes github.com/odedstein/meshes/tree/master/objects/bunny.
Original mesh at graphics.stanford.edu/data/3Dscanrep..
The Stanford 3D Scanning Repository. 2023. Armadillo. Downloaded modified
version from odedstein-meshes github.com/odedstein/meshes/tree/master/objects/
armadillo. Original mesh by V. Krishnamurthy and M. Levoy ""Fitting SmoothSurfaces to Dense Polygon Meshes"" (SIGGRAPH 1996).
TimEdwards. 2023. Springer.
Downloaded modified version from odedsteinmeshes github.com/odedstein/meshes/tree/master/objects/springer, originally fromthingiverse.com/thing:335658. Asset licensed under GPLv2.
Thibault Tricard, Jimmy Etienne, Cédric Zanni, and Sylvain Lefebvre. 2021. A brick in
the wall: Staggered orientable infills for additive manufacturing. In Proceedings of
the 6th Annual ACM Symposium on Computational Fabrication. 1–8.
Qing Xing, Ergun Akleman, Gabriel Taubin, and Jianer Chen. 2012. Surface covering
curves. In Proceedings of the Eighth Annual Symposium on Computational Aesthetics
in Graphics, Visualization, and Imaging. 107–114.
YahooJAPAN. 2013a. Jellyfish. Downloaded from https://www.thingiverse.com/thing:182222. Asset licensed under CC BY 4.0..
YahooJAPAN. 2013b. Snail. Downloaded from thingiverse.com/thing:182275. Assetlicensed under CC BY 4.0..

Surface-Filling Curve Flows via Implicit Medial Axes •YahooJAPAN. 2019. Plane.
Downloaded modified version from odedsteinmeshes github.com/odedstein/meshes/tree/master/objects/plane, originally fromthingiverse.com/thing:182252. Asset licensed under CC BY 3.0..
YahooJAPAN. 2021. Scorpion.
Downloaded modified version from odedsteinmeshes github.com/odedstein/meshes/tree/master/objects/scorpion, originally fromthingiverse.com/thing:182363. Asset licensed under CC BY 3.0..
Yuan Yan and Yasamin Mostofi. 2016. Efficient clustering and path planning strategies
for robotic data collection using space-filling curves. IEEE Transactions on Control
of Network Systems 4, 4 (2016), 838–849.
Chris Yu, Henrik Schumacher, and Keenan Crane. 2021. Repulsive Curves. ACM Trans.
Graph. 40, 2, Article 10 (may 2021), 21 pages. https://doi.org/10.1145/3439429
Cem Yuksel. 2015. Sample Elimination for Generating Poisson Disk Sample Sets.
Computer Graphics Forum 34, 2 (2015), 25–32. https://doi.org/10.1111/cgf.12538
arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.12538
Haisen Zhao, Fanglin Gu, Qi-Xing Huang, Jorge Garcia, Yong Chen, Changhe Tu,
Bedrich Benes, Hao Zhang, Daniel Cohen-Or, and Baoquan Chen. 2016. Connected
Fermat Spirals for Layered Fabrication. ACM Trans. Graph. 35, 4, Article 100 (jul
2016), 10 pages. https://doi.org/10.1145/2897824.2925958
Haisen Zhao, Hao Zhang, Shiqing Xin, Yuanmin Deng, Changhe Tu, Wenping Wang,
Daniel Cohen-Or, and Baoquan Chen. 2018. DSCarver: Decompose-and-SpiralCarve for Subtractive Manufacturing. ACM Trans. Graph. 37, 4, Article 137 (jul 2018),
14 pages. https://doi.org/10.1145/3197517.3201338
Fanchao Zhong, Yonglai Xu, Haisen Zhao, and Lin Lu. 2023. As-Continuous-As-Possible
Extrusion-Based Fabrication of Surface Models. ACM Trans. Graph. 42, 3, Article 26(mar 2023), 16 pages. https://doi.org/10.1145/3575859

A ENERGY MODIFICATIONS
A.1 Field-alignment energy

radiusradiusradiusradiusradiusradiusradiusradiusradiusradiusFig. 28. We can obtain curves with varying global shapes by changing thepower orders 𝑝 and 𝑞 of the two energy terms (see Eq. (13)).

in which case the Hessian and gradient required for the descentdirection can be computed via automatic differentiation. As shownin Figure 28, we observed that smaller 𝑝 and larger 𝑞 encourageround corners and sharpen the radius histogram.

B

A user may wish to align the curve to a vector field 𝑋 on the input
surface by adding a simple energy

DISCUSSION ON THE GRADIENT

energy vs. wall clock time
Eq.14 + Newton

𝑛−1
∑︁

1
∥ (𝑅𝑖 𝛾𝑖 − 𝑄𝑖 𝛾𝑖+1 ) × (𝑅𝑖 𝑋𝑖 + 𝑄𝑖 𝑋𝑖+1 ) ∥ 22 ,
𝛽𝐹
𝑙
𝑖=1 𝑖

11

fixed medial axis + gradient descent(11)

where 𝑋𝑖,𝑖+1 is the vertex-blended vector field on 𝛾𝑖,𝑖+1 and 𝛽 𝐹 is
the weight which we set to 0.5 to 1 in our experiments. Intuitively,
it encourages the segment 𝛾𝑖 𝛾𝑖+1 to be parallel to 𝑋𝑖 and 𝑋𝑖+1 in therigidly transformed plane.

A.2

Biharmonic bending energy

A user may also wish to encourage straighter segments by adding
a biharmonic bending energy. Let 𝜃𝑖 be the angle of ∠𝛾𝑖 −1𝛾𝑖 𝛾𝑖+1
projected onto the tangent space of 𝛾𝑖 , and 𝑃𝑖 be the rigid transformation matrix that maps 𝛾𝑖 −1 to the point (𝑙𝑖 −1 cos 𝜃𝑖 , 𝑙𝑖 −1 sin 𝜃𝑖 , 0)
while aligning the normal vector 𝑁𝑖 −1 with the 𝑧-axis. Intuitively,
𝑃𝑖 maps 𝛾𝑖 −1 onto the same plane as we evaluate the geodesic curve
length of 𝛾𝑖 𝛾𝑖+1 . We formulate our biharmonic bending energy as𝑛−1∑︁𝑙𝑖 −1 + 𝑙𝑖 𝑃𝑖 𝛾𝑖 −1 − 𝑅𝑖 𝛾𝑖 𝑅𝑖 𝛾𝑖 − 𝑄𝑖 𝛾𝑖+1 2
𝛽𝐵
−
,
2
𝑙𝑖 −1
𝑙𝑖
2
𝑖=1

(12)

where 𝛽𝐵 is the energy weight which we set to 10 to 1000 in our
experiments. We note that both the field-alignment and biharmonic
bending energies are sparse and quadratic, which does not have alarge effect on the runtime.

Eq.14 + gradient descent

fixed medial axis + preconditioning

Fig. 29. Using our approximated Hessian as a preconditioner leads to robust
and faster convergence, compared to treating the medial axis position as afunction as in Eq. 14 or conducting the mere gradient descent as in Eq. 15.

In our gradient computation, we apply two modifications to themere gradient descent of Eq. 2; we (1) fix the unrolled medial axis
𝑚𝑖′± for each step and (2) use the inverse of the (approximated)
Hessian H𝑘 as a preconditioner.
If we do not fix 𝑚𝑖′± , Eq. 7 can be rewritten as
𝐸𝑀 (𝜸 ) ≈𝑛∑︁𝑙𝑖 −1 + 𝑙𝑖 2

A.3

A user can change the global shape of the curve by changing thesquare order of the energies into more general 𝑝 and 𝑞 powers,

𝑖

𝑙𝑖

𝑝

∥𝑅𝑖 𝛾𝑖 − 𝑄𝑖 𝛾𝑖+1 ∥ 2 +

𝑛 
∑︁
𝑖=1

(14)

𝑖=1

Power of the energy

∑︁ 1


2
2
𝑟𝑖+ (𝛾𝑖 ) 2 + 𝑟𝑖− (𝛾𝑖 ) 2 ,𝑞𝑞
𝛾𝑖 − 𝑚𝑖′+ 2 + 𝛾𝑖 − 𝑚𝑖′− 2 , (13)∥𝛾 −𝛾 ∥ 22
where 𝑟𝑖± (𝛾𝑖 ) = 𝛾 −𝛾𝑗 · 𝑖(±𝐵
is the radius of the Euclidean medial( 𝑗 𝑖)𝑖)sphere grown from 𝛾𝑖 that touches 𝛾 𝑗 . However, we observed thatthe gradient directions become so noisy that even the smallest stepwould cause self-intersection, which completely stops the flow due
to our intersection avoidance (Figure 29).
ACM Trans. Graph., Vol. 43, No. 4, Article . Publication date: July 2024.

12•

Noma et al.

Under fixing 𝑚𝑖′± , one can also conduct the mere gradient descentd𝑘 = −g𝑘 ,(15)equipped with an Armijo line search. While this gives qualitatively
similar results, our preconditioned gradient evolves faster, especiallyin regions with thin features (see Figure 29).
Additionally, our approximated Hessian is banded near the diagonals, especially when the curve is singly connected and has an open
end. In this case, a user can take advantage of a Cholesky solver
tailored for symmetric banded matrices (e.g., the dpbtrf functionin LAPACK [Anderson et al. 1999]). In practice, we observed an

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

approximately 30% speed up in the linear solve compared to using
a normal sparse Cholesky solver.

C

COMPARISON WITH [Yu et al. 2021]

The algorithm of Yu et al. [2021] achieved a significant speed-up
thanks to their multigrid method. However, they neither explained
nor implemented how to apply it to the surface constraint. To approximate the effect of the multigrid method in Figure 5, we firstmeasured the runtime with and without the multigrid method on an
example that is supported, and then multiplied the speed-up factor
by the runtime with the surface constraint.