DOI: 10.1111/cgf.

70206

Eurographics Symposium on Geometry Processing 2025 COMPUTER GRAPHICS forum

M. Attene and S. Sellán Volume 44 (2025), Number 5
(Guest Editors)

Arrange and Traverse Algorithm for Computation of Reeb Spaces

of Piecewise Linear Maps

Petar Hristov 1 , Daisuke Sakurai2 , Hamish Carr3 , Ingrid Hotz1 and Talha Bin Masood1
1 Scientific Visualization Group, Department of Science and Technology (ITN), Linköping University, Norrköping, Sweden
2 Fujitsu Research Institute, Fukuoka, Japan
3 School of Computer Science, University of Leeds, Leeds, United Kingdom

1 4 5 6
7 8 9
2
10

Figure 1: Computing a Reeb space and using it to extract features from a torus data set with two scalar functions. The first scalar function
f is the implicit equation of a torus and the second g is a horizontal plane. We show isosurfaces of f at different levels and contours of g
restricted to those isosurfaces. As f increases, the torus becomes a topological sphere and then splits into multiple components that disappear
at the boundary. The Reeb space (on the right) allows us to detect changes in the topology of the level sets of the combined map ( f , g). Each
two-dimensional sheet of the Reeb space corresponds to a feature. Extracting multiple fibers (level sets of ( f , g)) for several sheets of the
Reeb space yields the main topological features of the dataset — the two handles of the torus as well as its expansion into a topological
sphere. The colors of the sheets of the Reeb space match the colors of the fibers. Where sheets overlap, the colors blend together.
Abstract
We present the first combinatorial algorithm for efficiently computing the Reeb space in all dimensions. The Reeb space is a
higher-dimensional generalization of the Reeb graph, which is standard practice in the analysis of scalar fields, along with
other computational topology tools such as persistent homology and the Morse-Smale complex. One significant limitation
of topological tools for scalar fields is that data often involves multiple variables, where joint analysis is more insightful.
Generalizing topological data structures to multivariate data has proven challenging and the Reeb space is one of the few
available options. However, none of the existing algorithms can efficiently compute the Reeb space in arbitrary dimensions
and there are no available implementations which are robust with respect to numerical errors. We propose a new algorithm
for computing the Reeb space of a generic piecewise linear map over a simplicial mesh of any dimension called arrange and
traverse. We implement a robust specialization of our algorithm for tetrahedral meshes and evaluate it on real-life data.

1. Introduction climate simulation involving multiple interrelated variables such as

Topological data analysis (TDA) offers robust and efficient tools for humidity, pressure, and vertical velocity. To analyze and visualize
studying scalar fields, such as Reeb graphs, persistent homology, such data effectively, it is not sufficient to examine each scalar field
and Morse-Smale complexes [HLH*16; YMS*21]. However, real-
world application data is often multivariate. Take, for instance, a
© 2025 The Author(s). Computer Graphics Forum published by Eurographics - The European Asso-
ciation for Computer Graphics and John Wiley & Sons Ltd.
This is an open access article under the terms of the Creative Commons Attribution License, which
CGF 44-5 | e70206
permits use, distribution and reproduction in any medium, provided the original work is properly
cited.
2 of 17 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps

in isolation. Instead, it is crucial to understand how these scalar plane in section 5, and extend this to arbitrary dimensions in sec-
fields interact, revealing the system’s joint features and patterns. tion 6. Implementation details are discussed in section 7, and an
Generalizing existing topological tools to the multivariate set- evaluation is presented in section 8. We conclude with a summary
ting remains one of the biggest open challenges in computational and directions for future work in section 9 and section 10.
topology. While some generalizations of basic concepts exist, they
remain limited in practical applications. For example, persistent ho- 2. Background
mology can be extended to the multiparameter setting, but the lack
of a complete invariant has hindered its practical use. Similarly, In this section, we introduce relevant mathematical background
the Reeb graph generalizes to the Reeb space, which completely from general topology, PL topology, and computational topology.
describes the connectivity of level sets in higher dimensions. How-
2.1. Topological Spaces
ever, the structural complexity of the Reeb space poses significant
challenges for both computation and interpretation. Even in the Topology offers a fundamental mathematical framework for de-
simplest case, such as a bivariate map defined over a tetrahedral scribing shapes and spaces. Among the most important topologi-
mesh, the Reeb space forms a two-dimensional cell complex, con- cal spaces are d-manifolds, which locally resemble d-dimensional
sisting of vertices and edges, connected by two-dimensional sur- Euclidean space. Although manifolds possess a simple local struc-
faces known as sheets. ture, their global topological and geometric properties can be quite
State-of-the-art algorithms to compute the Reeb space either complex. For this reason, it is often beneficial to work with more
generalize earlier algorithms for computing Reeb graphs through structured objects, such as cell complexes, which offer a more man-
contouring [EHP08; TC17; CD14] or represent the Reeb space as a ageable framework for analysis and computation.
collection of nested Reeb graphs called the multi-dimensional Reeb A cell complex is a combinatorial structure formed by gluing
graph (MDRG) [CRS24]. Contouring-based methods are compu- together basic building blocks called cells via continuous maps.
tationally efficient for three-dimensional data but can suffer from Many manifolds, especially in three dimensions, can be represented
robustness issues, often resulting in approximate representations of as cell complexes, allowing the use of combinatorial methods to
the Reeb space. Methods based on the MDRG, on the other hand, study their properties. A particular type of cell complex is the sim-
have at worst cubic running time and rely on mathematical assump- plicial complex, which consists of piecewise linear (PL) elements
tions, which are impractical for many real-life data sets. None of the and is especially well-suited for practical implementations.
methods proposed so far can compute a practical representation of
the Reeb space in any dimension. 2.2. Piecewise Linear Topology
We propose the first combinatorial algorithm for computing and Let {v0 , . . . , vd } be d + 1 affinely independent points in Rd . The
efficiently representing the Reeb space of piecewise linear maps sum ∑di=0 λi vi , where ∑di=0 λi = 1 and λi > 0 for all i is called a
in arbitrary dimensions. We call it the arrange and traverse algo- convex combination. The coefficients λi are known as barycentric
rithm. Our approach operates under very general assumptions on coordinates. The set of all convex combinations defines a d-simplex
the input. We assume the input K is a simplicial mesh of dimen- σ, where the points {v0 , . . . , vd } are called the vertices of the sim-
sion d along with a multivariate generic PL map f : |K| → Rm . We plex. The number of vertices defines its dimension dim(σ) = d. A
build on the work of Hristov [Hri22], which characterizes combina- simplex defined by a subset of the vertices of σ is called a face of σ.
torial transitions of fibers. We leverage these insights to generalize A 0-simplex is a point, a 1-simplex is a line segment, a 2-simplex
state-of-the-art Reeb graph algorithms by Parsa [Par12] and Do- is a triangle, a 3-simplex is a tetrahedron, and a d-simplex is their
raiswamy et al. [DN09] to efficiently keep track of level sets. Our higher-dimensional generalization. Given two simplices σ1 and σ2 ,
algorithm runs in time O((m + 1)(Nm−1 )m−1 Nm log2 Nm ), where their join, denoted σ1 σ2 , is defined as the simplex spanned by the
Nm and Nm−1 are the number of simplices of dimension m and union of their vertex sets.
m − 1 respectively. A collection of simplices K is called a simplicial complex if ev-
For tetrahedral meshes with two input functions, a common ery face of each simplex σ1 ∈ K is also an element of K and the
practical case, we present the first numerically robust implemen- intersection of any two simplices σ1 , σ2 ∈ K is either empty or a
tation of a Reeb space algorithm. Its worst-case running time is face of both. The dimension d of a simplicial complex is the highest
O(Nt Ne log Ne ), where Ne and Nt are the number of edges and tri- dimension of any simplex in the complex. We will use the notation
angles. We eliminate a logarithmic factor from the general case Kk for all simplices in K with dimension k. A simplicial complex
by efficiently tracking one-dimensional fiber connectivity. In prac- where all simplices of dimension less than d are the face of a sim-
tice, this bound is not tight, because our implementation is output- plex of dimensions d is called a simplicial mesh. Three dimensional
sensitive with with respect to the size of the arrangement of the simplicial meshes are called tetrahedral meshes.
images of the edges in the range space. We also prove correctness A subcomplex of K is a simplicial complex that consists of a
and validate on synthetic and real-world datasets. subset of the simplices of K. We will make use of a particular type
In section 2, we provide the background from piecewise lin- of simplicial complex called the i-skeleton K(i) , which is the sub-
ear (PL) topology, followed by a review of related work on Reeb complex of K that includes all simplices of dimension i or lower.
graph and Reeb space algorithms in section 3. Our main contri- Topological neighborhoods in a simplicial complex are defined by
butions are summarized in section 4. We introduce our algorithm the star and link of a simplex. The star of a simplex, Star(σ), is the
for bivariate maps from three-dimensional simplicial meshes to the set of simplices that have σ as a face. Generally, it is not a simplicial

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps 3 of 17

complex because it does not contain the faces of all of its simplices. the Reeb graph correspond to equivalence classes of regular level
We refer to the set of faces of simplices in the star, which do not sets components and the vertices correspond to singular level set
have σ as a face as the link of σ, denoted as Link(σ). components where connectivity changes. See subsection 3.1 for an
example of Reeb graph computation.
A simplicial complex allows us to model a complex geomet-
ric space from simple primitives. The underlying space of a sim-
plicial complex, denoted |K|, is the union of all its simplices 2.4. Multifield Topology
|K| = σ∈K σ. When the underlying space of a simplex matches
S

that of a manifold, we say that K triangulates that manifold. To When multiple scalar PL functions f1 , . . . , fm are defined on a d-
guarantee the regular structure of the stars and links, we can im- dimensional simplicial complex K, they can be combined into a
pose the stronger condition that K is a combinatorial d-manifold. single multivariate PL map f = ( f1 , . . . , fm ) : K → Rm . We re-
In a combinatorial d-manifold, the link of every vertex triangulates strict our attention to the case where d > m. A PL map with two
the (d − 1) sphere and is itself a combinatorial (d − 1) manifold. scalar functions is referred to bivariate. A PL map is generic if it is
injective on the simplices of dimension m or less, and the images
of those simplices are in general position in the range [EHP08].
2.3. Scalar Field Topology A bivariate PL map is generic if the images of any three vertices
Functions that preserve PL structure are called PL functions. We are non-collinear and the images of the interiors of any three edges
can define a PL function on a simplicial complex by only using val- are non-concurrent [Hri22]. To avoid ambiguity, we will refer to
ues at the vertices and interpolating to the rest of the points within the images of edges in the range space as segments throughout the
remainder of the paper.
simplices using barycentric interpolation. Let f¯ : K(0) → R be a
function defined on the vertex set of K and let x ∈ |K| be such that
x ∈ σ, where σ is a simplex in K. We can then define a PL function
f such that f (x) = ∑ki=0 λi f (vi ), where λi are the barycentric coor-
dinates of x. PL functions are also called scalar since their range
is the real line. A scalar PL function is generic when all vertices
are mapped to distinct values. This ensures no degenerate cases
and can be guaranteed efficiently in practice with symbolic per-
turbation techniques such as the simulation of simplicity [EM90].
Any scalar PL function can be approximated arbitrarily well by a
generic scalar PL function.
Critical points of generic PL functions occur only at vertices. Figure 2: Example simplicial complex (left) and bivariate generic
They are characterized using the lower and upper link of a vertex. PL map (right). The three tetrahedra abv1 v2 , abv2 v3 and abv3 v4
The lower link of a vertex u is a subcomplex induced by the vertices form the star of the edge ab. The upper link of ab is the subcomplex
of the link v with a lower value f (v) < f (u) than u. Conversely, the {v1 , v4 }, induced by the vertices in the star that are mapped above
upper link consists of vertices in the link with higher function value. the segment f (ab). The vertices in the star mapped below the seg-
The topology of these links (lower or upper), captured via their ment f (ab) induce the lower link {v2 , v3 , v2 v3 }.
reduced Betti numbers [EH09], provides a complete classification
of the vertices of K. For a triangulated 2-manifold K, vertices are
Level sets of multivariate PL maps are defined as the preimage
either regular, local minima, local maxima, or saddles.
of a point in the range f −1 (p) = {x ∈ |K| : f (x) = p ∈ Rm }. We
We can study the structure of scalar generic PL functions through refer to such level sets as fibers, following terminology from fiber
the topology of their level sets. A level set is the preimage of a topology [Sae04], and the corresponding point p in the range as a
point in the range f −1 (c) = {x ∈ |K| : f (x) = c}. When K is two- fiber point. A simplex is said to be active with respect to a fiber
dimensional, a level set is generically a set of contours, and when if it is intersected by that fiber. We refer to the connected compo-
K is three-dimensional, a level set is generically a set of surfaces nents of a fiber as its fiber components. A fiber f −1 (p) is generi-
called an isosurface. The topology of a level set changes at critical cally a simplicial complex of dimensions d − m [EHP08]. Fibers
vertices — at local minima and maxima, new connected compo- can be computed exactly for tetrahedral meshes [KTCG17] or re-
nents appear or disappear, and at saddles, multiple components can constructed otherwise [KLAK23]. In the bivariate case, taking the
merge or split. Saddles where more than two components meet or preimage of a polygon produces a surface swept by the fiber points
split can be broken up into multiple simple saddles, by editing the along the polygon called a fiber surface [CGT*15]. The preimages
input K and adding additional simplices. Scalar generic PL maps of higher dimensional geometry in multivariate maps are known as
without such saddles are called simple. While they are preferred in feature level sets [JH19].
theoretical work, they can complicate algorithmic development, so
The higher dimensional analog of the set of critical points of a
they are usually allowed in practical algorithms [CWSA16].
scalar function is the Jacobi set [EH02]. In the smooth case, the Ja-
To capture the topology of level sets of a scalar generic PL cobi set is the set of points where the gradients of the m component
function, we can use a topological construction called the Reeb functions are parallel — it is a smoothly embedded m − 1 manifold
graph [BDF*08]. The Reeb graph is a contraction of each con- in the domain. In the PL case, the Jacobi set is a subset of the m − 1
nected component of each level set to a single point. The edges in simplices. It is defined using the topology of the upper and lower

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
4 of 17 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps

lated [HS13] into a simplicial complex. However, that triangulation

can be too fine for practical use, so the Reeb space can be also rep-
resented as an m-dimensional cell complex.
In order to represent the Reeb space as a cell complex we project
the Jacobi set to the Reeb space via the quotien map that contracts
fibers. The image of the Jacobi set in the Reeb space is called the
Jacobi structure [CCDG14]. We create a cell complex where the
highest-dimensional cells, called sheets, are maximal areas that do
not intersect the Jacobi structure. They correspond to equivalence
classes of regular fiber components. The lower-dimensional cells
are the Jacobi structure. They correspond to equivalence classes of
singular fiber components. See Figure 3 for examples.
The Reeb space has been studied in the smooth case in singu-
larity theory and fiber topology, but less so in the PL case. Only
two works are known to the authors related to the study of the
structure of the Reeb space for simple generic PL maps [EHP08]
and more generally for generic PL maps [Hri22]. We describe the
state of the art in algorithms for computing the Reeb space in sub-
section 3.2. Examples of constructing simple Reeb spaces can be
Figure 3: Changes in fiber connectivity across a Jacobi edge in the found on pages 12 and 13 in the following paper [CTW20] and on
domain. The first column shows the location of the fiber point in page in 3 the following paper [CD14]. We show an example of how
the range relative to a segment (the image of an edge). The remain- the Reeb space relates to the Reeb graph in Figure 7.
ing columns show the corresponding fibers for different types of
Jacobi edges. Bold lines highlight the fiber components in a small
neighborhood around the Jacobi edge. The bottom row depicts the 3. Related Work
associated Reeb space, with one sheet per fiber component.
We begin with a review of background work on Reeb graph and
contour tree algorithms, then discuss methods for computing Reeb
spaces, and conclude with geometric arrangement techniques.
link of a simplex. The upper link of a simplex σ of dimension m − 1
is the subcomplex induced by the vertices in the link of σ, which
are mapped above the hyperplane defined by f (σ), where above is a 3.1. Reeb Graph Algorithms
matter of orientation and can be resolved with a geometric orienta-
tion test. The number of fiber components in a small neighborhood The Reeb graph was first introduced in the field of Morse theory
of σ for a fiber point slightly above f (σ) is equal to the number by Georges Reeb in 1946 [Ree46]. It was initially used to study the
of connected components of the upper link [EHP08; Hri22] The structure of manifolds equipped with Morse functions. The first al-
lower link is analogously defined using the vertices in the link of σ gorithm for computing the Reeb graph was proposed by Shinagawa
mapped below f (σ) and the number of fiber components in a small and Kunii [SK91]. This algorithm constructs the Reeb graph from
neighborhood of σ for a fiber point slightly below f (σ) is equal to cross-sectional contours of a surface in Θ(n2 ) time. The running
the number of connected components of the lower link. time was later improved to O(n log n) through more efficient con-
tour tracking [CEH*03]. Contours are modeled as graphs, whose
Simplices whose upper and lower links are simply connected nodes correspond to edges of the input and arcs to connecting tri-
are called regular, as they do not change the topology of passing angles. Representing these graphs as cyclic or linear linked lists
fibers [TC17]. Fibers that only intersect regular simplices are called enables efficient updates, such as edge insertions, removals, and
regular fibers. In contrast, Jacobi simplices can alter fiber topology. connectivity checks, using binary search trees in O(log n) time.
They are called definite when a fiber appears or disappears, simi-
lar to minima or maxima in scalar fields, and indefinite when fibers Follow-up algorithms [DN09] extend the computation of the
merge and split, analogous to saddles. If exactly two components Reeb graph to 3-manifolds, where contours become surfaces. The
merge and split, the Jacobi simplex is simple; otherwise, it is non- key idea is that the connectivity of a PL surface can be entirely
simple, analogous to simple and monkey saddles in scalar fields, derived from its 1-skeleton, which can be treated as a graph. Up-
respectively. In tetrahedral meshes, the Jacobi set is a subset of the date operations on graphs of surfaces can be processed efficiently
edges. Examples illustrating how regular and Jacobi edges affect with the tree-cotree data structure. This yields an algorithm with
fiber connectivity are shown in Appendix C and Figure 3. a running time of O(n log n + n log g(log log g)3 ), where g is the
maximum genus among all contours. Building on the idea of using
The topology of fibers can be described by a generalization of the
the 1-skeleton of preimages, the authors apply dynamic graph con-
Reeb graph called the Reeb space [EHP08], where all fiber compo-
nectivity techniques [HDT01] to track contour connectivity over
nents are contracted to a single point. Since there are m degrees of
n-manifolds, achieving a running time of O(n log n(log log n)3 ).
freedom in the range, the Reeb space locally has the structure of
an m manifold [EHP08]. Globally, the Reeb space can be triangu- The first deterministic optimal algorithm for computing Reeb

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps 5 of 17

difference is that instead of vertices, consider m − 1 dimensional

simplices for a multivariate generic PL map f : |K| → Rm , where
K is a combinatorial manifold of dimension d. The first step of
the algorithm is to compute the arrangement of the hyperplanes
induced by the image of the m − 1 skeleton of |K| in the range.
The connected components of the preimage of each face in the ar-
rangement correspond to regions of uniform fiber connectivity. In
the second step, each face is taken with multiplicity equal to the
number of connected components of its preimage and all faces are
glued those together along their shared boundaries, based on preim-
age connectivity. This produces a triangulation of the Reeb space,
that is significantly finer than is necessary and needs to be coars-
ened for practical use [EHP08].
The coarsest representation of a Reeb space is its canonical
stratification, which decomposes the space into manifold pieces
with uniform local topology. The algorithm uses a PL version
Figure 4: Computing a Reeb graph with contour track- of Goresky and MacPherson’s approach from intersection homol-
ing [CEH*03; DN09; Par12].
The preimage graph is computed for ogy [GM80]. A key step is the DoesBlend subroutine, which deter-
each interval f (vi ), f (vi+1 ) . The Reeb graph is induced by the mines if two strata can be merged. This relies on checking isomor-
change of connectivity in the preimage graphs as we sweep past all phism between triangulated spaces which, however, is an undecid-
vertices. able problem in dimensions greater than four [EHP08]
The Joint Contour Net (JCN) [CD14] captures the topology of
graphs in O(n log n) time was introduced by Parsa [Par12], leverag- multivariate data in scientific visualization by approximating the
ing offline dynamic graph connectivity [Epp94]. Building on earlier Reeb space. The range is quantized using a regular grid with a
work [DN09], contours are reduced to preimage graphs. Knowing fixed resolution. The connected components of the preimage of
the update order in advance allows for efficient reordering, enabling each cell in this quantization are then computed. These compo-
O(log n) updates. Figure 4 illustrates an example of using contour nents are connected through shared boundaries and contracted into
tracking to compute a Reeb graph. Reeb graph algorithms that are points to form the graph structure called the JCN. The algorithm’s
based on preimage graphs can be computed robustly using floating- running time is O(mNQ Nα(NQ N)), where m is the range dimen-
point arithmetic, as they rely only on floating-point comparisons. sion, N is the number of simplices, NQ is the number of fragments
Degenerate cases are handled using a symbolic perturbation tech- produced by the quantization, and α is the extremely slow-growing
nique called simulation of simplicity [EM90]. inverse Ackermann function. Empirical results show that NQ in-
creases significantly with the resolution of the quantization. Finally,
Other algorithms for computing Reeb graphs decompose the do- Mapper [MW16] is a related concept from topological data analysis
main into simply connected regions where the Reeb graph has that also approximates the Reeb space, but that is more applicable
no loops [TGS09; DN13]. In regions where the Reeb graph con- to point cloud data.
tains no loops, it simplifies to a contour tree, which can be com-
puted efficiently in any dimension [CSA03]. Finally, there are algo- The first practical Reeb space algorithm is an improved ver-
rithms for online computation [PSBM07], output-sensitive compu- sion of the original [EHP08] that is specialized to the bivariate
tation [DN12a] as well as randomized [HWW10] and parallel and case [TC17]. This algorithm generalizes the critical contouring
distributed computation [GFJT19; CWSA16; CRW22; CSLA15]. strategy [DN12b]. Instead of mapping all edges of the input to the
range, only the Jacobi edges are mapped, as they are the only ones
where fiber connectivity changes. For each Jacobi edge, a fiber sur-
3.2. Reeb Space Algorithms face (see section 2) of its image is extracted, referred to as a Jacobi
There are currently two main approaches for computing the Reeb fiber surface. The 3D volumes resulting from the segmentation of
space. The first approach is based on contouring, and the sec- the domain by the Jacobi fiber surfaces correspond to the preimages
ond is based on a hierarchical decomposition into multiple Reeb of sheets in the Reeb space.
graphs. In the contouring approach, the domain is segmented into
Computing all Jacobi fiber surfaces is embarrassingly parallel
regions of uniform fiber connectivity, which are glued together
with a work complexity of O(N j NT ), where N j is the number of
to obtain the Reeb space. Contouring can be based on the Ja-
Jacobi edges and NT is the number of tetrahedra in the input. Ja-
cobi edges [TC17], or more generally on higher dimensional sim-
cobi fiber surfaces can be used to efficiently approximate the sheets
plices [EHP08], or on quantizing the range [CD14]. In the decom-
in the Reeb space by growing connected components of the input
position approach [CRS24], the Reeb space is represented as a
mesh vertices until they encounter an edge intersected by a Jacobi
nested hierarchy of Reeb graphs. Refer to Table 1 for an overview
fiber surface. However, to compute the exact Reeb space, we must
of the key properties of each algorithm.
accurately track the intersections of Jacobi fiber surfaces [CRS24].
The first Reeb space algorithm [EHP08] uses a similar contour- This can be done by using algorithms for three-dimensional surface
ing approach to the first Reeb graph algorithm [SK91]. The key arrangements or, as noted by the authors in [TC17, Section 3.3], via

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
6 of 17 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps

Algorithm Edelsbrunner et al. Carr et al. Tierny et al. Chattopadhyay et al. Ours
[EHP08] [CD14] [TC17] [CRS24]
Domain Combinatorial d-dimensional PL 3-manifold Triangulation of a closed d-dimensional
d-manifold simplicial mesh 3-manifold simplicial mesh
PL Map Generic Rm Rm Generic R2 Simple and generic R2 Generic Rm
Output Reeb space Quantized Reeb Reeb space Reeb space Reeb space
space
Complexity Polynomial O(NQ Nα(NQ N)) O(N j NT ) + O(?) O(N 2 + NK j log N + O(Ne Nt log Nt )
(m = 2) Ncmax )
Limitations Impractical for Unknown speed Practical issues Simple Jacobi set, no Computing the
m>4 of convergence with robustness boundary arrangement

Table 1: Overview of algorithms for computing the Reeb space and its approximations. The notation is as follows — d is the dimension
of the input; m is the dimension of the range; Nv , Ne , Nt , NT are the number of vertices, edges, triangles and tetrahedra in the input;
N = Nv + Ne + Nt + NT is the total number of simplices; N j is the number is the number of Jacobi edges in the input and K j is the number of
intersections of their images, the Jacobi segments in the range; cmax is an upper bound on the number of simplices in the link of any edge;
NQ is the number of fragments in the quantization Q of the JCN algorithm [CD14]; α is the slowly growing inverse Ackermann function. The
running time of the Jacobi fiber surfaces algorithm [TC17] includes an additional term O(?), which is not specified in the original paper.
This term accounts for the final step of computing the intersections of Jacobi fiber surfaces to obtain the Reeb space sheets.

constrained triangulations. Both methods however increase time which cannot be embedded in three dimensions, limiting applica-
complexity beyond O(N j NT ) (the overall bound is not specified) bility for visualizing spatial data. Finally, the Jacobi set is required
and reduce available parallelism. to be simple, which requires the addition of new simplices to sim-
While the algorithm runs efficiently in practice, its robustness plify non-simple indefinite edges.
is not guaranteed. Numerical precision issues can lead to incor-
rect combinatorial structures in the Reeb space. Unlike scalar Reeb 3.3. Geometric Arrangements
graphs, Reeb space computation requires 2D orientation tests (e.g.,
A finite set of geometric objects S, such as lines, curves, surfaces,
for Jacobi sets), which demand exact arithmetic and filtering for
or hyperplanes, decomposes Euclidean space into connected open
performance [GOT18]. The problem is exacerbated when Jacobi
cells [GOT18]. This decomposition is called the arrangement A(S)
fiber surfaces and their intersections introduce many new vertices
of S. The combinatorial structure of an arrangement is given by the
with exact coordinates in 3D. Robust handling of these challenges
incidence graph of the arrangement [GOT18], where the nodes are
has yet to be addressed in the literature.
the open cells and the arcs are their adjacency. Algorithms for com-
The final approach for computing the Reeb space is based on puting arrangements are some of the most well-studied in compu-
computing a multi-dimensional Reeb graph (MDRG) [CRS24]. tational geometry, especially in lower dimensions [GOT18].
The MDRG is a hierarchical decomposition of the Reeb space into
multiple nested Reeb graphs across different dimensions. In the The arrangement of line segments in the plane can be repre-
bivariate case, the Reeb graph of the first scalar function is com- sented with a doubly-connected edge list (DCEL) that defines each
puted initially. Subsequently, the points at which the topology of the face (two-dimensional cell) through its boundary edges. The in-
nested Reeb graphs changes are identified. For each such point, the tersection points of n segments can be computed in O(n log n + k)
Reeb graph of the second scalar function, restricted to the preim- time [Bal95; CE92] and the full arrangement can be computed
age of that point under the first function, is then computed. The re- in O((n + k) log n) time [DCVO08a], where k is the number of
sulting MDRG structure is homeomorphic to the Reeb space. The intersections of the line segments [DCVO08a]. The number of
running time is O(N 2 + NK j log N + Ncmax ), where N is the total intersections is at worst quadratic in the number of line seg-
number of simplices, K j is the number of intersections of images of ments [DCVO08a]. The arrangement of n hyperplanes can be com-
Jacobi edges in the range and cmax is an upper bound of the number puted in time O(nm ) [GOT18].
of simplices in the link of an edge.
While this algorithm provides a nice mathematical decomposi- 4. Contributions
tion, it has a few practical limitations. First, in the expression for In this paper, we present a novel combinatorial algorithm for Reeb
the worst case running time, the quadratic term N 2 can be large space computation. While all previous algorithms [EHP08; TC17;
even for small datasets and K j is at worst cubic in Ne . This is be- CRS24] assume a type of triangulated manifold, we relax this as-
cause the number of Jacobi edges N j is at worst on the order of sumption to allow for a more general type of input — a simpli-
the number of edges Ne and K j is at worst N 2j . Second, the authors cial mesh. Furthermore, the only existing general Reeb space algo-
assume the domain is a closed PL 3-manifold (without boundary), rithm [EHP08], cannot represent the Reeb space with its canonical

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps 7 of 17

(a) Arrangement (b) Preimage graphs (c) Correspondence graph (d) Reeb space.

Figure 5: The algorithm proceeds in three main steps. (b) Preimage graphs Gi for each face in the arrangement are computed using a
breadth-first search traversal; (c) A union-find data structure is then used to determine correspondences between fiber components across
adjacent faces; (d) These correspondences yield the sheets of the Reeb space along with their connectivity.

stratification (with the least number of elements) for m > 4, be- c2 ⊆ f −1 (p2 ) be two regular fiber components in K for two fiber
cause of an undecidable subroutine. In high-dimensional cases that points p1 , p2 ∈ R2 . We say that c1 corresponds to c2 if there is a
algorithm can only produce a triangulation with significantly more path P ⊂ R2 between p1 and p2 such that both c1 and c2 belong
elements than needed [EHP08]. Our algorithm also works for any to the same connected component C of f −1 (P) and C does not
dimension of the domain d and the range m such that d > m (oth- intersect the Jacobi set. The surface C is a continuous deformation
erwise the Reeb space is trivial). We avoid the use of any undecid- between the two fiber components.
able procedures by considering preimages of points as graphs, for The correspondence between regular fiber components across all
which there are efficient connectivity algorithms. This allows us to fibers defines an equivalence relation whose equivalence classes
represent the Reeb space with a cell complex that potentially has have a one-to-one matching with the sheets of the Reeb space. In-
more elements than the canonical stratification but allows for more deed, any two corresponding fiber components belong to the same
practical implementations and applications in any dimension. sheet of the Reeb space, and conversely, any two fiber components
We show that our algorithm lends itself to a robust implementa- that lie on the same sheet correspond to each other. Under this for-
tion. We present the first implementation of a Reeb space algorithm mulation, computing the sheets of the Reeb space reduces to identi-
in the case of tetrahedral meshes that is robust with respect to nu- fying all regular fiber components and determining their correspon-
merical errors. Unlike previous implementations, which require an dences to obtain the equivalence classes. Once these classes are es-
extension to handle this robustness in 3D, we reduce this problem tablished, computing adjacency between sheets becomes straight-
to the robustness of computing an arrangement of line segments in forward: two sheets are adjacent if and only if their corresponding
2D, which is well-studied and already has robust implementations equivalence classes share a boundary in the domain K.
[The24]. The algorithmic complexity of our algorithm is on par
To compute the fiber components of fibers, we will use a combi-
with competing approaches, and it is output-sensitive with respect
natorial representation known as a preimage graph [Par12, Section
to the size of the arrangement of the image of the edge set of the
4.1]. The preimage graph G p of a point p ∈ R2 is a graph whose ver-
input mesh. Finally, our approach works for generic maps and does
tices are the active simplices of K for the fiber f −1 (p) and whose
not require special treatment of non-simple Jacobi edges [CRS24].
edges are the adjacency of the active simplices in K. In the general
case where the fiber f −1 (p) does not intersect any edge or vertex,
5. Bivariate Reeb Space Algorithm its intersection with a tetrahedron is a line segment connecting two
For clarity, we will begin with a description of the algorithm in points located on distinct boundary triangles of the tetrahedron. In
the base case, where a generic PL map with two components is de- this case, we can simplify the preimage graph by only using the ac-
fined over a tetrahedral mesh. The general case is a direct extension, tive triangles and connecting the ones that are the faces of the same
which we will discuss in section 6. In order to compute the Reeb tetrahedron, see Figure 6 left column.
space, we will reframe its sheets as equivalence classes of corre- The active simplices of a fiber f −1 (p) can only change when
sponding fibers with the same topology. First, we use a geometric the fiber point p crosses a segment in the range (i.e., the image of
arrangement algorithm to compute regions in the range where cor- an edge under the map f ). Consequently, all fiber points not sepa-
responding fibers are mapped. Then, we traverse the arrangement rated by such segments have isomorphic preimage graphs. Maximal
to combine the fibers of those regions into the equivalence classes regions where the preimage graphs G p remain isomorphic are the
that represent the sheets of the Reeb space and their connectivity. faces of the geometric arrangement A = Arrangement( f (K(1) )) in-
duced by the image of the 1-skeleton of K. We use the notation GF
5.1. Algorithm Description to denote the preimage graph of a face F.
Let K be a tetrahedral mesh, let f = ( f1 , f2 ) : |K| → R2 be a generic Next, we describe how preimage graphs change across incident
PL map and let R be the Reeb space of f . Let c1 ⊆ f −1 (p1 ) and faces in the arrangement [Hri22]. Let F1 and F2 be two adjacent

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
8 of 17 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps

faces in A, and let e ⊆ f (ab) be the segment of the arrangement

separating them, where ab is an edge in K. When a fiber point
crosses e, the resulting change in the fiber is localized to the star
of the edge ab (see section 2). Specifically, the triangles abu, where
u belongs to the lower link of ab, cease to be active, while the tri-
angles abv, where v belongs to the upper link of ab, become active.
Conversely, when crossing the segment in the opposite direction,
the triangles abu with u in the upper link become inactive, and those
with v in the lower link become active. Refer to Appendix C – Fig-
ure 9 for illustrations of the various cases.

This leads to a combinatorial algorithm for computing the preim-

age graphs {GF }F∈A of all faces in the arrangement A. The local
changes in fibers allow for an efficient computation of the preimage
graph of a face, given the preimage graph of one of its neighboring
faces. The algorithm begins at the outside face, which has an empty
preimage graph, and then traverses the faces of A using a graph
search, updating the preimage graphs by adding and removing tri-
angles as necessary. This can be done efficiently by representing
preimage graphs as linked lists and tracking changes in their con-
nectivity with binary search trees in logarithmic time [Epp94].

Before proceeding with the remainder of the algorithm, we il-

lustrate the computation of preimage graphs with a small example.
Consider the three tetrahedra sharing the edge ab and the bivari- Figure 6: An example of how preimage graphs evolve along a path
ate generic PL map f in Figure 6. The left column shows three through the arrangement. The step IDs correspond to those in Fig-
preimage graphs, one per row, corresponding to consecutive fiber ure 5. Each row shows a representative fiber for the preimage graph
points along a traversal. We begin with G0 = ∅, representing the of a face in the arrangement. On the left, the upper and lower stars
preimage graph for a fiber point located in the outside face of the of the segment being crossed are shown. Due to the direction of
arrangement. This graph is empty because the fiber of a point out- traversal, triangles in the upper star are removed (striped fill pat-
side the image f (|K|) of the domain K is empty. In the first row, tern), while triangles in the lower star are added (solid fill).
the fiber point is moved to a neighboring face by crossing the seg-
ment f (bv4 ). The triangles that become active are highlighted in
light gray. Since bv4 is a definite edge, a new fiber component is (or vice versa) along e. We establish a one-to-one correspondence
born, and no triangles become inactive. With slight abuse of nota- between all fiber components that have not been affected because
tion, we update the preimage graph as G3 = G0 + {abv4 , bv3 v4 }. they have not changed. If the edge ab is regular, we also establish
The edges of G3 are implicitly defined by the adjacency of trian- correspondence between the two affected fiber components. If the
gles in the mesh K. In the second row, the fiber point crosses the edge ab is Jacobi, we do not establish any further correspondence.
segment f (av1 ), the image of another definite edge, where a sec- Finally, once we have computed all preimage graphs {GF }F∈A
ond fiber component appears. The preimage graph is updated to and the correspondence graph H, we are ready to compute the
G4 = G3 + {av1 v2 , abv1 }, having two fiber components. Finally, Reeb space R. To compute the sheets of R, we need the connected
in the third row, the fiber point crosses the image of ab. The edge components of H. Every connected component of H represents an
ab is indefinite because its upper link has two connected compo- equivalence class of corresponding fibers, and every equivalence
nents and the lower link has one, which causes the two fiber com- class of corresponding fibers represents a sheet in the Reeb space.
ponents to merge into one. The updated graph is given by G8 = To obtain the geometry of the sheets, we take the union of the faces
G4 + {abv3 , abv2 } − {abv1 , abv4 }. The two triangles {abv1 , abv4 }, of the preimage graphs that have a fiber component in that equiv-
which become inactive, are marked with a striped pattern. alence class (see Theorem B.3). If we keep track of where definite
and indefinite edges prevent correspondence in the construction of
Having computed the preimage graphs for all faces, the next step
H we can determine which sheets of R are adjacent.
is to establish correspondences between their fiber components.
To this end, we define the correspondence graph H, whose ver- We present pseudocode for our Reeb space algorithm in Ap-
tices represent all connected components of the preimage graphs pendix A and proof of correctness in Appendix B. We show an
{GF }F∈A . Initially, this graph contains no edges. Let F1 and F2 be example of the steps of our algorithm for our running example in
two adjacent faces in the arrangement A, separated by a segment Figure 5. Next, we discuss the algorithmic complexity.
of the arrangement e ⊆ f (ab), where ab is an edge in K. Let GF1
and GF2 denote the corresponding preimage graphs for these faces,
5.2. Complexity
with fiber components {α1 , . . . αn } and {β1 , . . . βm } respectively. A
component αi or β j is said to be affected if it contains any triangle First, we compute the arrangement A of the image of the 1-skeleton
that is either added or removed during the transition from F1 to F2 of K. This can be done with a deterministic output-sensitive al-

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps 9 of 17

gorithm [DCVO08b] in time O((Ne + k) log Ne ), where Ne is the In the general case, the Jacobi set of f is a subset of the simplices
number of edges in K and k is the number of intersections of their of dimension m − 1 [EH02]. Therefore, the connectivity of a fiber
images in R2 . can only change when we cross the image of a simplex σ of dimen-
sion m − 1 or one of its faces. How a fiber changes as we cross the
In the second step, we compute all preimage graphs {GF }F∈A ,
image of σ is analogous to the bivariate case. All m-dimensional
and, in the third step we compute the correspondence graph H. Both
simplices in the star of σ with a vertex in the upper link become
steps consider each half-edge e in the arrangement at most once
inactive, while those with a vertex in the lower link become active
and perform at most deg(ab) dynamic graph operations, where ab
– or vice versa, depending on the direction of travel. The change
is the originating edge for the segment such that e ⊆ f (ab) and
of the fiber is, again, local to the simplex, whose image we cross.
deg is the number of adjacent triangles in K. Representing fibers
as linked lists allows us to track their connectivity with balanced The first step of our general Reeb space algorithm is to discard
binary search trees [CEH*03] for log Nt time per update or query, the simplices of K of dimension higher than m + 1 and to compute
where Nt is the number of triangles in K. the arrangement A of the image of the m − 1 skeleton of K. Since
If we define kab as the number of times the image of the edge the image of the m − 1 skeleton of K is a collection of (m − 1)-
ab is intersected in the range then each triangle abc ∈ K is con- dimensional hyperplane segments, the arrangement will be a cell
sidered at most kab + kbc + kac times. Summing all over triangles, complex of dimension m. All preimage graphs in the interior of a
we get a total running time of ∑abc∈K3 (kab + kbc + kca ) log Nt for face (cell of maximum dimension) of that arrangement are isomor-
constructing the preimage graphs and the correspondence graph, phic because none of their corresponding fibers intersect any Jacobi
where K3 is the set of triangles of K. Every term kab in this simplices. This allows us to extend the definition of a preimage
sum is present deg(ab) times, so we can reformulate the sum as graph to a face of the arrangement, like in the bivariate case.
∑ab∈K2 kab deg(ab), where K2 is the set of edges of K. By def- In the second step, we compute all preimage graphs {GF }F∈A
inition kab < Ne , since every edge can be crossed by at most by traversing the faces of the arrangement by using their incidence
Ne − 1 other edges and with a counting argument we can show via the (m − 1)-dimensional cells. In the third step, we compute the
that ∑ab∈K2 deg(ab) = 3Nt , so then ∑ab∈K2 kab deg(ab) < 3Nt Ne , correspondence graph H, by examining the correspondence of fiber
which yields a worst-case running time of O(Ne Nt log Nt ). components of adjacent faces. Finally, we represent the sheets of
Finally, we can compute the sheets of the Reeb space in linear the Reeb space as a union of faces that all correspond to connected
time in the size of the correspondence graph H, since this only in- components of H. The incidence relations of the sheets of the Reeb
volves computing the connected components of H. We can identify space are similarly obtained from the incidence relations of their
the adjacency of the sheets by keeping track of the Jacobi edges in corresponding faces in the range.
the previous step with no extra cost.
Overall, the worst-case time complexity of our Reeb space al- 6.1. Complexity
gorithm for a bivariate generic PL map defined over a tetrahedral The algorithmic complexity of the arrangement algorithm of n hy-
mesh is O(Ne Nt log(Nt )). An important note here, is that this worst perplanes in Rm is O(nm ) (see subsection 3.3). For hyperplane seg-
case bound is not tight and the practical running time of our algo- ments, we can compute the arrangement of the hyperplanes they
rithm depends on size of the arrangement. induce and then merge adjacent faces if their shared boundary is
not in the original input. This can be done in linear time in the size
6. General Algorithm of the arrangement.

In this section, we discuss the generalization of our algorithm to a In the construction of the preimage graphs and the correspon-
simplicial mesh K of dimension d and a multivariate generic PL dence graph, every hyperplane segment s of dimension m − 1 in
map f = ( f1 , f2 , ..., fm ) : K → Rm where d > m. This generaliza- the arrangement is considered at most once with deg(σ) updates or
tion relies on the Skeleton lemma [EHP08], which states that the queries, where σ ∈ K is a simplex of dimension m − 1 such that
Reeb space of f is homeomorphic to the Reeb space of the restric- s ⊆ f (σ), and deg counts the number of m-dimensional simplices
tion of f to the m + 1 skeleton of K. A special case of this lemma that have σ as a face. We define kσ to be the number of m − 1 di-
has been used to compute Reeb graphs in any dimension by dis- mensional cells in the arrangement A that are contained in f (σ).
carding all simplices of dimension higher than two [DN09]. At most Nm−1 other simplices can intersect σ, where Nm−1 is the
number of simplices of dimension m − 1, so kσ ≤ (Nm−1 )m−1 .
Using the Skeleton lemma, we can simplify the computation of
the Reeb space by discarding all simplices of dimension higher than Next, we determine how many times a simplex τ of dimension
m + 1. If K is a simplicial mesh of dimension d then, K(m+1) is a m is considered by our algorithm. We denoted the boundary of τ as
simplicial mesh of dimension m + 1. Let f¯ : |K(m+1) | → Rm be ∂τ = {σ1 , . . . , σm+1 }, where each σi is of dimension m − 1. With
the restriction of f . Note that f¯ is generic because f is generic. this notation τ is considered ∑σ∈∂τ kσ times. Summing over all
Furthermore, we know from previous work [EHP08] that the fiber m-dimensional simplices in K we obtain ∑τ∈Km ∑σ∈∂τ kσ , where
of the restricted map f¯−1 (p) is the 1-skeleton of the fiber f −1 (p) Km is the set of m-dimensional simplices in K. Every term kσ ap-
for p ∈ Rm . Therefore, we can represent f¯−1 (p) with a preimage pears exactly deg(σ)
times. Therefore, we can rewrite the sum as
graph and track changes in its connectivity with dynamic graph ∑σ∈Km−1 kσ deg(σ) . Since kσ ≤ (Nm−1 )m−1 and by a counting


connectivity algorithms [HDT01]. Next, we describe how this con- argument we can show that ∑σ∈Km−1 deg(σ) = (m + 1)Nm we
have ∑σ∈Km−1 kσ deg(σ) ≤ (m + 1)(Nm−1 )m−1 Nm .


nectivity changes in a high-dimensional range.

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
10 of 17 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps

Update and query operations for computing the preimage graphs in Figure 1, where we have shown the evolution of the isosurfaces
and the correspondence graph can be performed in time log2 (Nm ) of f . The contours of g, restricted to isosurfaces of f , correspond
for general graphs, where Nm is the number of simplices of di- to the fibers of the bivariate map ( f , g).
mension m [HDT01]. The total running time of our general Reeb
The Reeb space in this example has 456 sheets in total, but only
space algorithm is therefore O((m + 1)(Nm−1 )m−1 Nm log2 Nm ). In
ten of them are large enough to be clearly visible in Figure 1. The
higher dimensions, of course, the size of the arrangement suffers
rest of the 445 sheets correspond to noise from the resolution and
from a combinatorial explosion. Even so, our algorithm does not
artifacts from the tetrahedralization. In order to determine whether
contain any undecidable subroutines, which is in contrast to the
sheets of the Reeb space correspond to features in the data, we have
only other method for computing the Reeb space [EHP08] in di-
shown some of the fibers of the first five sheets. For each of these
mensions higher than four.
sheets, we take a vertical line, fixing an isovalue of f , and sample
a range of values for g, thereby generating multiple fibers to vi-
7. Implementation sualize the volume they sweep through the domain. Sheets 1 and
2, which overlap in our particular visualization of the Reeb space,
We implemented a proof of concept version of the algorithm in correspond to the two handles of the torus. Sheet 3 corresponds to
C++ for the bivariate case described in section 5. The input to the the topological sphere we get by expanding the torus until the hole
algorithm is a tetrahedral mesh K, along with two scalar values per in the center pinches and disappears. The rest of the ten sheets cor-
vertex representing the PL map f : |K| → R2 . To reduce the number respond to the fibers breaking apart and disappearing as the isosur-
of degenerate cases, we perturb the values of f numerically. faces of f expands to intersect the boundary. We show an example
For the arrangement of line segments, we used the CGAL library, of how the Reeb graphs of each isosurface in Figure 1 correspond
version 6.0.1 [The24]. We used CGAL’s kernel for exact predicates to the Reeb space in Figure 7.
and exact constructions (EPEC) to ensure a topologically correct
arrangement [WBF*24]. Using faster kernels with floating point
operations is feasible, as long as the relative order of intersections
in the arrangement is handled robustly.
We track preimage graph connectivity in close to linear time in
their size with union-find. This suffices for examples we present in
section 8, where the size of fibers was significantly less than the size
of the input. Since fibers are intersections of isosurfaces this is also
supported by previous work that shows estimates on the number of
triangles produced by marching cubes of n0.82 [HDD06].
If we combine the computation of the preimage graphs and the
correspondence graph in a single BFS search, we can only keep the (a) Reeb space. (b) Embeded Reeb space.
preimage graphs in the current level or higher of the BFS search
and discard the rest to save memory. Finally, computing fibers and Figure 7: Embedding of the top left quadrant of the Reeb space
assigning them to sheets in the Reeb space can be done efficiently from Figure 1 in three dimensions. The Reeb space captures the
by taking all the active triangles from the preimage graph of the evolution of the Reeb graphs (bottom right (b)) of the second scalar
face of the arrangement that contains the fiber point and computing function of the bivariate map, restricted to the contours of the first.
barycentric coordinates.

For comparison, we also computed the Reeb space with the only
8. Experiments
other available Reeb space implementation [TC17] known to us
In this section, we demonstrate the application of our algorithm to from the Topology ToolKit (TTK) library [TFL*18]. While the
two example datasets. The first example is a bivariate field of an TTK implementation was considerably faster, almost none of the
analytical function and the other dataset comes from chemistry. We sheets matched the ones we computed nor the main features in the
conclude this section with a performance evaluation. data set, which can be verified by hand. Out of a total of 44 sheets
computed with the TTK implementation, one sheet covered the en-
tire area in the range where the whole of the Reeb space projects.
8.1. A Synthetic Example Most of the points in the domain correspond to this sheet, which is
A torus equipped with a height function is one of the first exam- incorrect. Most other sheets are substantially smaller and the han-
ples one encounters when learning about Reeb graphs [EH09]. We dles of the torus, arguably the most important features of this data
can construct a similar example for the Reeb space. We will use set, were not identified.
a tetrahedralization of a regular cubic grid with 20480 tetrahedra.
For each vertex in the tetrahedral mesh, we sample two analytic
8.2. A Case Study from Chemistry
functions
p f and g. The first function f is the implicit function of a
torus ( x2 + y2 − R)2 + z2 = r2 and the second function g = z is a Next, we consider a dataset from Chemistry related to electronic
horizontal plane cutting through the torus. You can see an example bond identification where bivariate analysis has been shown to play

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps 11 of 17

(a) Input bivariate field (b) Reeb space (c) Extracted features

Figure 8: (a) The input is a bivariate field—electron density and reduced gradient—computed for ethanediol. Translucent grey and orange
isosurfaces represent electron density and reduced gradient, respectively. Electron density maxima indicate atom positions, marked by red
(Oxygen), grey (Carbon), and white (Hydrogen) spheres. While reduced gradient reveals some chemical bonds, identifying all interactions
and distinguishing covalent from non-covalent bonds requires bivariate analysis. (b, c) Reeb space analysis automatically identifies key
features like atoms and bonds. The non-covalent interaction between the two OH groups appears as a distinct sheet in the Reeb space.

a crucial role [JKM*10; BCTa16]. The two scalar fields of interest were unable to find any chemically relevant features corresponding
are electronic density and reduced gradient computed for a small to the extracted volumes representing sheets of the Reeb space.
molecule called ethanediol. Examining each of these scalar fields
in isolation can provide some hints about important features such as
the location of atoms and bonds, see Figure 8a. But a more detailed 8.3. Performance evaluation
analysis can only be done when both fields are considered simulta- We report the runtime performance of our prototype implementa-
neously. An initial study on this dataset using interactive bivariate tion in Table 2. To evaluate scalability we generate down-sampled
analysis techniques such as continuous scatter plots and fiber sur- versions of the Torus and Ethanediol datasets introduced in subsec-
faces was successful in distinguishing different types of atoms and tion 8.1 and subsection 8.2. For both datasets, runtime and mem-
the type of molecular interactions between atoms [CGT*15]. How- ory usage scale super-linearly with input size, consistent with the
ever, the authors left an interesting note in their discussion: “We worst-case complexity of O(Ne Nt log Nt ). We make an important
leave automatic detection of chemical features in range space for observation here that the performance also depends on the com-
the future, as it is likely to depend on topological analysis of fiber plexity of the Reeb space. For instance, the Torus dataset, with
surfaces” [CGT*15, Section 8]. With the Reeb space, we can now fewer Reeb sheets, requires significantly less time and memory than
take first steps toward addressing this challenge of automatic detec- the Ethanediol dataset of similar mesh size. This is further reflected
tion of chemical features. in smaller intermediate structures – such as the arrangement and
correspondence graph – for the Torus dataset. Overall, while the
Using the segmentation of the spatial domain induced by sheets
evaluation suggests near-quadratic scaling, it also highlights the al-
in Reeb space, we can not only distinguish between the type of
gorithm’s output-sensitive behavior. The second key observation
atoms and bonds, we can extract each atom and bond individually
from Table 2 is that the bottleneck in the performance of this al-
and study their local topological environment, see Figure 8. The
gorithm is not the computation of the arrangement itself, but the
Reeb space shown in Figure 8b consists of ~17000 sheets, how-
traversal of the arrangement to construct the correspondence graph.
ever we limited our focus to the top 50 sheets based on their area
in the range space. These top sheets captured all the key features
such as the atoms and covalent and non-covalent interactions in 9. Discussion and Future Work
the molecule. A selected subset of these sheets are highlighted in
the spatial domain using fibers extracted automatically using a few To the best of our knowledge, this is the first implementation of a
sample points as shown in Figure 8c. We highlight one particular Reeb space algorithm that offers robustness guarantees for generic
bond, the non covalent interaction between the two OH groups in PL maps. This foundational work opens up a lot of future possi-
ethanediol which is not easy to detect using traditional scalar field bilities, for example, towards improving algorithmic efficiency and
analysis methods. In our analysis, this bond automatically comes towards applications of Reeb space-guided exploration and analy-
out as a single sheet in Reeb space. The atoms, as well as other co- sis of multivariate fields. However, we note some limitations of the
valent bonds in this molecule, are also identifiable, but we note that current algorithm and its proof-of-concept implementation. One of
the corresponding Reeb sheets suffer from over segmentation. This the main limitations of the current implementation is the practi-
can be addressed via Reeb space simplification which is an open cal performance and memory usage as evident from Table 2. The
research challenge in itself and not the focus of our current work. key bottleneck is the computation of the preimage graphs and the
correspondence graph since it has to be done for every face in the
Lastly, we note that similar to the synthetic data example, the arrangement. In the future, we plan to address this by investigating
Reeb space computed by TTK had some inconsistencies and we techniques that utilize the arrangement of the image of the Jacobi

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
12 of 17 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps

#Tets #Faces Size of H #Sheets Time (s) Memory (Mb)

Arrangement Reeb space Arrangement Reeb space
625 15,845 35,721 56 0.06 0.48 30 83
1080 32,468 73,455 95 0.12 1.25 54 113
1715 55,909 109,981 105 0.22 2.40 91 157
Torus

2560 102,993 224,526 128 0.43 5.16 164 244

5000 242,259 527,694 201 1.12 15.34 390 498
10985 670,588 1,375,264 322 3.72 57.41 1,100 1,300
20480 1,524,240 3,094,903 456 9.93 170.88 2,500 2,900
400 22,471 171,421 274 0.09 2.32 39 80
875 75,245 685,615 548 0.30 9.85 114 212
Ethanediol

1440 146,185 1,329,737 847 0.65 19.07 215 376

2800 386,890 4,058,648 1,240 2.18 71.72 563 976
5800 1,081,479 11,800,052 2,257 8.06 298.49 1,600 2,600
13260 3,076,395 32,106,056 3,913 26.45 1401.26 4,600 8,500
25200 7,472,595 131,086,345 8,436 46.98 3577.06 11,000 29,000

Table 2: Performance metrics for the Torus and Ethanediol datasets. Each dataset is downsampled multiple times to evaluate the algorithm’s
scalability. We report input/output sizes (#Tets, #Sheets), intermediate structures (#Faces in the arrangement, size of correspondence graph
H), computation times and memory usage for both the arrangement and for subsequently computing the Reeb space.

set, which is much smaller in practice, and with parallelization of Acknowledgments

the arrangement traversals with parallel BFS algorithms [BM06;
This work was partially supported by the Wallenberg Autonomous
BM11]. Another challenge that must be addressed before wider
Systems and Software Program (WASP) funded by the Knut and
adoption for the analysis of practical datasets is multi-scale sim-
Alice Wallenberg Foundation; the Swedish e-Science Research
plification for automatic feature identification.
Center (SeRC); the Swedish Research Council (VR) grants 2019-
In the examples we tested, the numerical perturbation was 05487 and 2023-04806; and, JSPS KAKENHI Grant Numbers
enough to avoid most degenerate cases. In particular, there were 22K18267 and 20K19809. The chemistry dataset used in subsec-
no issues with degeneracies for the torus example and only two tion 8.2 is taken from the TTK [TFL*18] dataset repository with
out of ~8500 sheets of the Reeb space were degenerate and could acknowledgments to Roberto Alvarez Boto. Finally, we would like
not be computed correctly for the chemistry dataset. Ensuring that to thank the anonymous reviewers for their valuable comments.
f is truly generic, however, and that the algorithm remains robust
for any PL map is an open problem that requires a more system-
References
atic solution. In the future, we plan to address this with symbolic
perturbation techniques such as simulation of simplicity [EM90]. [Bal95] BALABAN, I VAN J. “An Optimal Algorithm for Finding Seg-
ments Intersections”. Proceedings of the Eleventh Annual Symposium
on Computational Geometry - SCG ’95. Vancouver, British Columbia,
Canada: ACM Press, 1995, 211–219. ISBN: 978-0-89791-724-7. DOI:
10. Conclusion 10.1145/220279.220302. (Visited on 03/24/2025) 6.
[BCTa16] B OTO, ROBERTO A., C ONTRERAS -G ARCÍA, J ULIA, T IERNY,
In this paper, we presented a combinatorial algorithm for comput- J ULIEN, and AND, J EAN -P HILIP P IQUEMAL. “Interpretation of the re-
ing the Reeb space in all dimensions called arrange and traverse. duced density gradient”. Molecular Physics 114.7-8 (2016), 1406–1414.
The first step of this algorithm is to compute the geometric arrange- DOI : 10.1080/00268976.2015.1123777 11.
ment of the images of the m − 1 dimensional simplices in the range. [BDF*08] B IASOTTI, S., D E F LORIANI, L., FALCIDIENO, B., et al. “De-
In the second step, we track the changes in the connectivity of fibers scribing Shapes by Geometrical-Topological Properties of Real Func-
combinatorially during a traversal of the faces of the arrangement. tions”. ACM Computing Surveys 40.4 (Oct. 2008), 1–87. ISSN: 0360-
0300, 1557-7341. DOI: 10.1145/1391729.1391731. (Visited on
Finally, we combine fibers with the same connectivity into equiva-
04/08/2025) 3.
lence classes that represent the sheets of the Reeb space. We derive
[BM06] BADER, DAVID A. and M ADDURI, K AMESH. “Designing Mul-
the adjacency of the sheets through the adjacency of their images
tithreaded Algorithms for Breadth-First Search and st-connectivity on
in the range. the Cray MTA-2”. 2006 International Conference on Parallel Process-
ing (ICPP’06). 2006, 523–530. DOI: 10.1109/ICPP.2006.34 12.
We implement a specialization of this algorithm to the practi-
cal case for a generic PL map from a tetrahedral mesh to the plane. [BM11] B ULUÇ, AYDIN and M ADDURI, K AMESH. “Parallel breadth-first
search on distributed memory systems”. Proceedings of 2011 Interna-
Finally, we discuss implementational details and present two exam- tional Conference for High Performance Computing, Networking, Stor-
ples of how to use the Reeb space for analyzing bivariate data. We age and Analysis. SC ’11. Seattle, Washington: Association for Com-
show that the Reeb space does allow us to detect significant fea- puting Machinery, 2011. ISBN: 9781450307710. DOI: 10 . 1145 /
tures in both analytic and real-world data sets. In both examples, 2063384.2063471 12.
our application is more robust than competing implementations.

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps 13 of 17

[CCDG14] C HATTOPADHYAY, A MIT, C ARR, H AMISH, D UKE, DAVID, [DN12a] D ORAISWAMY, H. and NATARAJAN, V. “Output-Sensitive Con-
and G ENG, Z HAO. Extracting Jacobi Structures in Reeb Spaces. struction of Reeb Graphs”. IEEE Transactions on Visualization and
2014. DOI: 10 . 2312 / EUROVISSHORT . 20141156. (Visited on Computer Graphics 18.1 (Jan. 2012), 146–159. ISSN: 1077-2626. DOI:
07/10/2024) 4. 10.1109/TVCG.2011.37. (Visited on 02/17/2025) 5.
[CD14] C ARR, H AMISH and D UKE, DAVID. “Joint Contour Nets”. IEEE [DN12b] D ORAISWAMY, H. and NATARAJAN, V. “Output-Sensitive Con-
Transactions on Visualization and Computer Graphics 20.8 (Aug. struction of Reeb Graphs”. IEEE Transactions on Visualization and
2014), 1100–1113. ISSN: 1941-0506. DOI: 10.1109/TVCG.2013. Computer Graphics 18.1 (Jan. 2012), 146–159. ISSN: 1077-2626. DOI:
269. (Visited on 10/22/2024) 2, 4–6. 10.1109/TVCG.2011.37. (Visited on 03/23/2025) 5.
[CE92] C HAZELLE, B ERNARD and E DELSBRUNNER, H ERBERT. “An [DN13] D ORAISWAMY, H. and NATARAJAN, V. “Computing Reeb
Optimal Algorithm for Intersecting Line Segments in the Plane”. Jour- Graphs as a Union of Contour Trees”. IEEE Transactions on Visualiza-
nal of the ACM 39.1 (Jan. 2, 1992), 1–54. ISSN: 0004-5411, 1557-735X. tion and Computer Graphics 19.2 (Feb. 2013), 249–262. ISSN: 1077-
DOI : 10.1145/147508.147511. URL : https://dl.acm.org/ 2626. DOI: 10.1109/TVCG.2012.115. (Visited on 03/23/2025) 5.
doi/10.1145/147508.147511 (visited on 05/28/2025) 6.
[EH02] E DELSBRUNNER, H ERBERT and H ARER, J OHN. “Jacobi Sets of
[CEH*03] C OLE -M C L AUGHLIN, K REE, E DELSBRUNNER, H ERBERT, Multiple Morse Functions”. Foundations of Computational Mathemat-
H ARER, J OHN, et al. “Loops in Reeb Graphs of 2-Manifolds”. Pro- ics, Minneapolis (2002), 37–57 3, 9.
ceedings of the Nineteenth Annual Symposium on Computational Ge-
[EH09] E DELSBRUNNER, H ERBERT and H ARER, J OHN. Computational
ometry. San Diego California USA: ACM, June 2003, 344–350. ISBN:
Topology. Providence, Rhode Island: American Mathematical Society,
978-1-58113-663-0. DOI: 10.1145/777792.777844. (Visited on
Dec. 2009. ISBN: 978-0-8218-4925-5 978-1-4704-1208-1. DOI: 10 .
03/23/2025) 4, 5, 9.
1090/mbk/069. (Visited on 03/04/2025) 3, 10.
[CGT*15] C ARR, H AMISH, G ENG, Z HAO, T IERNY, J ULIEN, et al. “Fiber
[EHP08] E DELSBRUNNER, H ERBERT, H ARER, J OHN, and PATEL, A MIT
Surfaces: Generalizing Isosurfaces to Bivariate Data”. Computer Graph-
K. “Reeb Spaces of Piecewise Linear Mappings”. Proceedings of the
ics Forum 34.3 (2015), 241–250. DOI: 10.1111/cgf.12636 3, 11.
Twenty-Fourth Annual Symposium on Computational Geometry. College
[CRS24] C HATTOPADHYAY, A MIT, R AMAMURTHI, YASHWANTH, and Park MD USA: ACM, June 2008, 242–250. ISBN: 978-1-60558-071-5.
S AEKI, O SAMU. An Algorithm for Fast and Correct Computation of DOI : 10.1145/1377676.1377720. (Visited on 03/14/2025) 2–7, 9,
Reeb Spaces for PL Bivariate Fields. 2024. arXiv: 2403 . 06564 10.
[cs.CG]. URL: https : / / arxiv . org / abs / 2403 . 06564 2,
[EM90] E DELSBRUNNER, H ERBERT and M ÜCKE, E RNST P ETER. “Sim-
5–7.
ulation of Simplicity: A Technique to Cope with Degenerate Cases
[CRW22] C ARR, H AMISH A., RÜBEL, O LIVER, and W EBER, G UNTHER in Geometric Algorithms”. ACM Transactions on Graphics 9.1 (Jan.
H. “Distributed Hierarchical Contour Trees”. 2022 IEEE 12th Sympo- 1990), 66–104. ISSN: 0730-0301, 1557-7368. DOI: 10.1145/77635.
sium on Large Data Analysis and Visualization (LDAV). Oklahoma City, 77639. (Visited on 04/18/2021) 3, 5, 12.
OK, USA: IEEE, Oct. 2022, 1–10. ISBN: 978-1-6654-9156-3. DOI: 10.
[Epp94] E PPSTEIN, D. “Offline Algorithms for Dynamic Minimum Span-
1109/LDAV57265.2022.9966394. (Visited on 03/23/2025) 5.
ning Tree Problems”. Journal of Algorithms 17.2 (Sept. 1994), 237–250.
[CSA03] C ARR, H AMISH, S NOEYINK, JACK, and A XEN, U LRIKE. ISSN : 01966774. DOI : 10 . 1006 / jagm . 1994 . 1033. (Visited on
“Computing Contour Trees in All Dimensions”. Computational Geome- 03/24/2025) 5, 8.
try 24.2 (Feb. 2003), 75–94. ISSN: 09257721. DOI: 10.1016/S0925-
[GFJT19] G UEUNET, C HARLES, F ORTIN, P IERRE, J OMIER, J ULIEN, and
7721(02)00093-7. (Visited on 03/23/2025) 5.
T IERNY, J ULIEN. “Task-Based Augmented Reeb Graphs with Dynamic
[CSLA15] C ARR, H AMISH, S EWELL, C HRISTOPHER, L O, L I -TA, and ST-Trees”. Eurographics Symposium on Parallel Graphics and Visual-
A HRENS, JAMES. Hybrid Data-Parallel Contour Tree Computation. ization (2019), 11 pages. ISSN: 1727-348X. DOI: 10 . 2312 / PGV .
Tech. rep. LA-UR-15-24579. Los Alamos National Laboratory, 2015 5. 20191107. (Visited on 02/17/2025) 5.
[CTW20] C ARR, H AMISH, T IERNY, J ULIEN, and W EBER, G UNTHER H. [GM80] G ORESKY, M ARK and M AC P HERSON, ROBERT. “Intersection
“Pathological and Test Cases for Reeb Analysis”. Topological Meth- Homology Theory”. Topology 19.2 (1980), 135–162. ISSN: 0040-9383.
ods in Data Analysis and Visualization V. Ed. by C ARR, H AMISH, DOI : 10.1016/0040-9383(80)90003-8 5.
F UJISHIRO, I SSEI, S ADLO, F ILIP, and TAKAHASHI, S HIGEO. Cham:
[GOT18] G OODMAN, JACOB E., O’ROURKE, J OSEPH, and T ÓTH,
Springer International Publishing, 2020, 103–120. ISBN: 978-3-030-
C SABA D., eds. Handbook of Discrete and Computational Geometry.
43036-8 4.
Third edition. Discrete Mathematics and Its Applications. Boca Raton
[CWSA16] C ARR, H AMISH A., W EBER, G UNTHER H., S EWELL, London New York: CRC Press, 2018. ISBN: 978-1-4987-1139-5 978-1-
C HRISTOPHER M., and A HRENS, JAMES P. “Parallel Peak Pruning for 351-64591-1 6.
Scalable SMP Contour Tree Computation”. 2016 IEEE 6th Symposium
[HDD06] H AMISH, C ARR, D UFFY, B RIAN, and D ENBY, B RAIN. “On
on Large Data Analysis and Visualization (LDAV). Baltimore, MD, USA:
Histograms and Isosurface Statistics”. IEEE Transactions on Visualiza-
IEEE, Oct. 2016, 75–84. ISBN: 978-1-5090-5659-0. DOI: 10 . 1109 /
tion and Computer Graphics 12.5 (2006), 1259–1266. DOI: 10.1109/
LDAV.2016.7874312. (Visited on 03/23/2025) 3, 5.
TVCG.2006.168 10.
[DCVO08a] D E B ERG, M ARK, C HEONG, OTFRIED, VAN K REVELD,
[HDT01] H OLM, JACOB, D E L ICHTENBERG, K RISTIAN, and T HORUP,
M ARC, and OVERMARS, M ARK. Computational Geometry: Algorithms
M IKKEL. “Poly-Logarithmic Deterministic Fully-Dynamic Algorithms
and Applications. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008.
for Connectivity, Minimum Spanning Tree, 2-Edge, and Biconnectivity”.
ISBN : 978-3-540-77973-5 978-3-540-77974-2. DOI : 10.1007/978-
Journal of the ACM 48.4 (July 2001), 723–760. ISSN: 0004-5411, 1557-
3-540-77974-2. (Visited on 03/24/2025) 6.
735X. DOI: 10.1145/502090.502095. (Visited on 03/24/2025) 4,
[DCVO08b] D E B ERG, M ARK, C HEONG, OTFRIED, VAN K REVELD, 9, 10.
M ARC, and OVERMARS, M ARK. Computational Geometry: Algorithms
[HLH*16] H EINE, C., L EITTE, H., H LAWITSCHKA, M., et al. “A Survey
and Applications. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008.
of Topology-based Methods in Visualization”. Computer Graphics Fo-
ISBN : 978-3-540-77973-5 978-3-540-77974-2. DOI : 10.1007/978-
rum 35.3 (2016), 643–667. DOI: https://doi.org/10.1111/
3-540-77974-2. (Visited on 03/18/2025) 9.
cgf.12933 1.
[DN09] D ORAISWAMY, H ARISH and NATARAJAN, V IJAY. “Efficient Al-
[Hri22] H RISTOV, P ETAR G EORGIEV. “Hypersweeps, Convective Clouds
gorithms for Computing Reeb Graphs”. Computational Geometry 42.6-7
(Aug. 2009), 606–616. ISSN: 09257721. DOI: 10.1016/j.comgeo. and Reeb Spaces”. PhD thesis. University of Leeds, June 2022. (Visited
2008.12.003. (Visited on 02/17/2025) 2, 4, 5, 9. on 07/08/2024) 2–4, 7.

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
14 of 17 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps

[HS13] H IRATUKA, J ORGE T and S AEKI, O SAMU. “Triangulating Stein [TFL*18] T IERNY, J ULIEN, FAVELIER, G UILLAUME, L EVINE, J OSHUA
Factorizations of Generic Maps and Euler Characteristic Formulas: A., et al. “The Topology ToolKit”. IEEE Transactions on Visualization
Dedicated to Professor Masahiko Suzuki on the Occation of His Six- and Computer Graphics 24.1 (2018), 832–842. DOI: 10.1109/TVCG.
tieth Birthday (Singularity Theory, Geometry and Topology)”. RIMS 2017.2743938 10, 12.
Kôkyûroku Bessatsu B38 (2013) 4. [TGS09] T IERNY, J., G YULASSY, A., and S IMON, E. “Loop Surgery
[HWW10] H ARVEY, W ILLIAM, WANG, Y USU, and W ENGER, for Volumetric Meshes: Reeb Graphs Reduced to Contour Trees”.
R EPHAEL. “A Randomized O ( m Log m ) Time Algorithm for IEEE Transactions on Visualization and Computer Graphics 15.6 (Nov.
Computing Reeb Graphs of Arbitrary Simplicial Complexes”. Pro- 2009), 1177–1184. ISSN: 1077-2626. DOI: 10.1109/TVCG.2009.
ceedings of the Twenty-Sixth Annual Symposium on Computational 163. (Visited on 02/17/2025) 5.
Geometry. Snowbird Utah USA: ACM, June 2010, 267–276. ISBN: [The24] T HE CGAL P ROJECT. CGAL User and Reference Manual. 6.0.1.
978-1-4503-0016-2. DOI: 10 . 1145 / 1810959 . 1811005. (Visited
CGAL Editorial Board, 2024. URL: https://doc.cgal.org/6.
on 02/17/2025) 5.
0.1/Manual/packages.html 7, 10.
[JH19] JANKOWAI, J OCHEN and H OTZ, I NGRID. “Feature Level-Sets:
[WBF*24] W EIN, RON, B ERBERICH, E RIC, F OGEL, E FI, et al. “2D Ar-
Generalizing Iso-Surfaces to Multi-Variate Data”. IEEE Transactions on
rangements”. CGAL User and Reference Manual. 6.0.1. CGAL Editorial
Visualization and Computer Graphics 26.2 (2019), 1308–1319. ISSN: Board, 2024. URL: https://doc.cgal.org/6.0.1/Manual/
1077-2626, 1941-0506, 2160-9306. DOI: 10 . 1109 / TVCG . 2018 .
packages.html#PkgArrangementOnSurface2 10.
2867488 3.
[YMS*21] YAN, L IN, M ASOOD, TALHA B IN, S RIDHARAMURTHY,
[JKM*10] J OHNSON, E RIN R., K EINAN, S HAHAR, M ORI -S ÁNCHEZ,
R AGHAVENDRA, et al. “Scalar Field Comparison with Topological De-
PAULA, et al. “Revealing Noncovalent Interactions”. Journal of the
scriptors: Properties and Applications for Scientific Visualization”. Com-
American Chemical Society 132.18 (2010). PMID: 20394428, 6498–
puter Graphics Forum 40.3 (2021), 599–633. DOI: 10 . 1111 / cgf .
6506. DOI: 10.1021/ja100936w 11.
14331 1.
[KLAK23] KOHLBRENNER, M., L EE, S., A LEXA, M., and K AZHDAN,
M. “Poisson Manifold Reconstruction — Beyond Co-dimension One”.
Computer Graphics Forum 42.5 (Aug. 2023), e14907. ISSN: 0167-7055,
1467-8659. DOI: 10.1111/cgf.14907 3.
[KTCG17] K LACANSKY, P., T IERNY, J., C ARR, H., and G ENG, Z. “Fast
and Exact Fiber Surfaces for Tetrahedral Meshes”. IEEE Transactions
on Visualization and Computer Graphics 23.7 (July 2017), 1782–1795.
ISSN : 1077-2626. DOI : 10.1109/TVCG.2016.2570215 3.

[MW16] M UNCH, E LIZABETH and WANG, B EI. Convergence between

Categorical Representations of Reeb Space and Mapper. Apr. 12, 2016.
DOI : 10.48550/arXiv.1512.04108. arXiv: 1512.04108. (Vis-
ited on 10/22/2024). Pre-published 5.
[Par12] PARSA, S ALMAN. “A Deterministic O(m log m) Time Algorithm
for the Reeb Graph”. Proceedings of the Twenty-Eighth Annual Sympo-
sium on Computational Geometry. SoCG ’12: Symposium on Computa-
tional Geometry 2012. Chapel Hill North Carolina USA: ACM, June 17,
2012, 269–276. ISBN: 978-1-4503-1299-8. DOI: 10.1145/2261250.
2261289. (Visited on 02/17/2025) 2, 5, 7.
[PSBM07] PASCUCCI, VALERIO, S CORZELLI, G IORGIO, B REMER,
P EER -T IMO, and M ASCARENHAS, A JITH. “Robust On-Line Compu-
tation of Reeb Graphs: Simplicity and Speed”. ACM SIGGRAPH 2007
Papers. San Diego California: ACM, July 2007, 58. ISBN: 978-1-
4503-7836-9. DOI: 10 . 1145 / 1275808 . 1276449. (Visited on
02/17/2025) 5.
[Ree46] R EEB, G. “Sur les points singuliers d’une forme de Pfaff com-
pletement integrable ou d’une fonction numerique [On the Singular
Points of a Completely Integrable Pfaff Form or of a Numerical Func-
tion]”. Comptes Rendus Acad. Sciences Paris 222 (1946), 847–849. URL:
https://cir.nii.ac.jp/crid/1571417125676878592 4.
[Sae04] S AEKI, O SAMU. Topology of Singular Fibers of Differentiable
Maps. Vol. 1854. Lecture Notes in Mathematics. Berlin, Heidelberg:
Springer Berlin Heidelberg, 2004. ISBN: 978-3-540-23021-2 978-3-540-
44648-4. DOI: 10.1007/b100393. (Visited on 03/04/2025) 3.
[SK91] S HINAGAWA, Y. and K UNII, T.L. “Constructing a Reeb Graph
Automatically from Cross Sections”. IEEE Computer Graphics and Ap-
plications 11.6 (Nov. 1991), 44–51. ISSN: 0272-1716. DOI: 10.1109/
38.103393. (Visited on 02/17/2025) 4, 5.
[TC17] T IERNY, J ULIEN and C ARR, H AMISH. “Jacobi Fiber Surfaces for
Bivariate Reeb Space Computation”. IEEE Transactions on Visualization
and Computer Graphics 23.1 (Jan. 2017), 960–969. ISSN: 1077-2626.
DOI : 10.1109/TVCG.2016.2599017. (Visited on 03/05/2025) 2,
4–6, 10.

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps 15 of 17

Appendix A: Pseudocode for the Reeb Space Algorithm Algorithm 3 Compute the correspondence graph.
1: function C OMPUTE C ORRESPONDENCE G RAPH(K, A, f , G)
Algorithm 1 Compute the Reeb Space of f : K → R2 2: H ← DisjointSet()
1: function C OMPUTE R EEB S PACE(K, f ) 3: for F ∈ A.faces do ▷ Initialize the vertices of H
2: K(1) ← 1-skeleton of K (vertices and edges). 4: for c ∈ G[F].connectedComponents do
3: A ← C OMPUTE 2DA RRANGEMENT( f (K(1) )) 5: H.insert(c)
4: G ← C OMPUTE P REIMAGE G RAPHS(K, f , A) 6: end for
5: H ← C OMPUTE C ORRESPONDENCE G RAPH(K, f , A, G) 7: end for
6: R ← C OMPUTE R EEB S PACE(K, f , A, G, H) 8: for e ∈ A.halfedges do ▷ Establish correspondence
7: return R 9: K ← e.face
8: end function 10: L ← e.twin.face
11: CK ← affected fibers of G(K).connectedComponents
12: CL ← affected fibers of G(L).connectedComponents
Algorithm 2 Compute preimage graphs of all faces in the ar- 13: if e is NOT Jacobi then
rangement 14: H.join(CK [1],CL [1])
1: function C OMPUTE P REIMAGE G RAPHS(K, f , A) 15: end if
2: G←[] 16: for c1 ∈ G(K).connectedComponents, c ∈ / CK do
3: F0 ← the unbounded cell of A 17: t ← any triangle in c1
4: Q = queue() 18: c2 ← G(L).findComponent(t)
5: Q.push(F0 ) ▷ All cells adjacent to ∂A. 19: H.join(c1 , c2 )
6: F0 .visited = true 20: end for
7: while Q = ̸ ∅ do 21: end for
8: Fc = Q.pop() 22: return H
9: for Fn ∈ Fc .neighbours do 23: end function
10: if Fn .visited is false then
11: continue ▷ Skip if the face has been visited Algorithm 4 Compute the Reeb space.
12: end if
1: function C OMPUTE R EEB S PACE(K, A, f , G, H)
13: e ← half-edge between Fc and Fn
14: c ← originating curve of e in A 2: Sheets ← [ ]
15: ab ← edge of K such that f (ab) = c 3: for F ∈ A.faces do
16: if e.orientation = ab.orientation then 4: for c ∈ G(K).connectedComponents do
17: G[Fn ] ← G[Fc ] + upperStar(e) − lowerStar(e) 5: Sheet[H.find(c)] ← F
18: else 6: end for
19: G[Fn ] ← G[Fc ] + lowerStar(e) − upperStar(e) 7: end for
20: end if ▷ If two adjacent faces are parts of different sheets, those
21: Q.push(Fn ) sheets are connected because they share a boundary
22: Fn .visited = true 8: Connections ← [ ]
23: end for 9: for F ∈ A.faces do
24: end while 10: for Fn ∈ F.neighbours do
25: return G 11: CF ← G(K).connectedComponents
26: end function
12: CFn ← G(L).connectedComponents
13: RF ← H.find(CF ) ▷ Component-wise
14: RFn ← H.find(CFn ) ▷ Component-wise
15: SF ← RF \ RFn
16: SFn ← RFn \ RF
17: for s1 ∈ SF do
18: for s2 ∈ SFn do
19: Connections[s1 ].insert(s2 )
20: Connections[s2 ].insert(s1 )
21: end for
22: end for
23: end for
24: end for
25: return {Sheets, Connections}
26: end function

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
16 of 17 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps

Appendix B: Algorithm Correctness and we only establish correspondence between fiber components
that are not separated by a Jacobi edge.
Let K be a three-dimensional simplicial complex, let f = ( f1 , f2 ) :
|K| → R2 be a generic PL map. Let {GF }F∈A be the set of all Conversely, suppose that d1 ⊆ f −1 (p1 ) and d2 ⊆ f −1 (p2 ) are two
preimage graphs, let H be the correspondence graph, and let R be corresponding fiber components for two points p1 and p2 outside
the Reeb space, all computed with the algorithm we have shown in in the image of the 1-skeleton of K. Let c1 and c2 be the two con-
section 5. We show that all those computations are correct. nected components of preimage graphs of the faces containing p1
and p2 that have the same active triangles as d1 and d2 . Our goal is
Lemma B.1 The preimage graphs {GF }F∈A are computed cor- to show that c1 and c2 are in the same connected component of H.
rectly. Since there is a correspondence between d1 and d2 there is a path P
path between p1 and p2 in R2 such both d1 and d2 are in the same
Proof Let F be a face in the arrangement with preimage graph GF .
connected component C of f −1 (P) and C does not intersect the
Let Tp = {t1 , . . . ,tn } be set of active triangles of a fiber f −1 (p)
Jacobi set. The sequence of faces intersected by P gives a path in
defined for a point p in the interior of F. Let VF be the vertex set of
the incidence graph of the arrangement. Since C has not intersected
GF . We wish to show that Tp = VF .
the Jacobi set, in each pair of incident intermediate faces on the
Let t ∈ Tp be an active triangle. Since t is active, its image f (t) path we have established a correspondence between two connected
must contain the point p. Since both F and f (t) contain p then components which in turn correspond to c1 and c2 . Therefore, c1
F ∩ f (t) ̸= ∅. Since all sides of f (t) are part of the arrangement, and c2 are in the same connected component of H.
then F ⊆ f (t). Consider the spanning tree of our graph traversal and
go up the parent list of the face F until we reach a face that shares
the boundary of f (t). Suppose that boundary is a half-edge e of the Theorem B.3 Reeb space R is computed correctly.
arrangement that originates from the image of the edge ab. Then ab Proof Let q : |K| → R be the quotient function that maps the points
is an edge of t, so let v be the third vertex of t. Since p ∈ f (t), the of every fiber to a point on the sheet that the fiber belongs to in the
points p and v are in the same half-plane of the line defined by ab. Reeb space. Let h : R → R2 be the map from Reeb space to the
So, our algorithm must have added the triangle t when crossing e. range such that f = h ◦ q. Let U be a sheet in the Reeb space. First,
Since e is the closest half-edge that shares a boundary with f (t) in we will show that we compute the geometry of U correctly. By the
the parent list by our assumption, our algorithm has not removed t. geometry of U we mean its shape h(U) ⊂ R2 in the range. Note that
Therefore t ∈ VF . we consider sheets to be closed, i.e. they include their boundary.
Conversely, suppose that t ∈ VF . Suppose that t is not active so that Let V be the equivalence class in terms of regular fiber component
p∈ / f (t). Again, consider the spanning tree of our graph traversal. correspondence that matches U. Since the union of the fiber compo-
Since t ∈ VF we must have added it in one of the parent graphs nents of V is an open set in K we will consider its closure |V| ⊂ K.
by crossing one of its boundary edges and must not have removed Let C be the connected component of H that matches V. Let
it after the last time it was added. Since we are now at F and we {α1 , . . . , αn } be the fiber components that correspond to vertices
have not removed f (t), then F ⊆ f (t). This is a contradiction since of C. Let {G1 , . . . , Gn } be the preimage graphs that contain those
p ∈ F ⊆ f (t). fiber components, such that αi ⊂ Gi . Finally, let {F1 , . . . , Fn } be the
Lemma B.2 The correspondence graph H is computed correctly. faces of the arrangement for each of the preimage graphs Gi respec-
tively. By our construction we have that f (|V|) = ∪ni=0 Fi . Since V
Proof The correspondence graph H is computed correctly when its matches U we have that q(|V|) = U. Therefore, h(q(|V|)) = h(U),
connected components match the classes of the equivalence rela- but since f = h ◦ q we have that f (|V|) = h(U). Finally, using that
tion defibed by fiber component correspondence. By Lemma B.1 f (|V|) = ∪ni=0 Fi , we obtain h(U) = ∪ni=0 Fi . This means that we can
the preimage graphs are computed correctly. represent sheets as the union of their faces in the arrangement.
Let c1 , c2 be two connected components of the preimage graphs Next, we discuss the adjacency of the sheets of the Reeb space.
G1 and G2 for the faces F1 and F2 . Suppose that c1 , c2 are in the Two sheets U and U ′ in the Reeb space are adjacent if and only
same connected component of H. Let p1 and p2 be two points in if their boundaries have a non-empty intersection. This can also
the interiors of F1 and F2 and let f −1 (p1 ) and f −1 (p2 ) be their be decided in R2 by whether ∂h(U) ∩ ∂h(U ′ ) is the empty set. Let
fibers. Finally, suppose that d1 ⊆ f −1 (p1 ) is the fiber component h(U) = ∪ni=0 Fi and h(U ′ ) = ∪ni=0 Fi′ be the faces in the arrangement
that has the same active triangles as c1 and d2 ⊆ f −1 (p2 ) is the for each sheet. The boundary of the sheet ∂U is equal to ∂ ∪ni=0
fiber component with the same active triangles as c2 We wish to Fi , because we assume that each sheet is simply connected. Two
show that d1 corresponds to d2 . sheets are adjacent if and only if ∂ ∪ni=0 Fi ∩ ∂ ∪ni=0 Fi′ ̸= ∅, which
Since c1 and c2 are in the same connected component of H, there completes this proof.
is a path P in the incidence graph of the arrangement between their
corresponding faces. Over this path, we have established a corre-
spondence between c1 and c2 via the connected components of the Appendix C: Local Change of Fiber Topology for Bivariate Maps
preimage graphs of the intermediate faces. If we take a point in each
face in that path and connect the points in the incident faces we ob-
tain a path P in R2 with starting point p1 and endpoint p2 . This
path satisfies the definition of correspondence between d1 and d2
since all fiber components within face interiors are homeomorphic

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.
 Petar Hristov et al. / Arrange and Traverse Algorithm for Computation of Reeb Spaces of Piecewise Linear Maps 17 of 17

Figure 9: Case table of how fiber connectivity changes in the local neighborhoods of different types of Jacobi edges ab. The first row lists
the type of Jacobi edges, and the second row shows the image of the star of the edge ab in the range. Triangles that contain ab and a vertex
from the lower link are shown in white and triangles that contain ab and a vertex from the upper link are shown in gray. The next three rows
show the fiber in the star of ab for a fiber point above the segment f (ab), at and below that segment. In the fiber above the segment all gray
triangles are active, and in the fiber below the segment, all white triangles are active. Note that there are infinitely many types on non-simple
indefinite Jacobi edges with n number of fiber components where n > 2.

© 2025 The Author(s).

Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.