definitiontheoremdefinition remarktheoremremark exampletheoremexample
lemmatheoremlemma corollartheoremcorollar propositiontheoremproposition 1295
Research algorithm/technique Nicolas Klenert and Daniel Baum are with Zuse
Institute Berlin. E-mail: {klenert $|$<EMAIL_ADDRESS>Verena Lepper is with
Ägyptisches Museum und Papyrussammlung. E-mail<EMAIL_ADDRESS>Klenert et al.: Extraction of 2D Manifolds from Thin-Layer Structures

# A Local Iterative Approach for the Extraction of 2D Manifolds
from Strongly Curved and Folded Thin-Layer Structures

Nicolas Klenert0009-0006-4443-8620  Verena Lepper  and Daniel
Baum0000-0003-1550-7245

###### Abstract

Ridge surfaces represent important features for the analysis of 3-dimensional
(3D) datasets in diverse applications and are often derived from varying
underlying data including flow fields, geological fault data, and point data,
but they can also be present in the original scalar images acquired using a
plethora of imaging techniques. Our work is motivated by the analysis of image
data acquired using micro-computed tomography (µCT) of ancient, rolled and
folded thin-layer structures such as papyrus, parchment, and paper as well as
silver and lead sheets. From these documents we know that they are
2-dimensional (2D) in nature. Hence, we are particularly interested in
reconstructing 2D manifolds that approximate the document’s structure. The
image data from which we want to reconstruct the 2D manifolds are often very
noisy and represent folded, densely-layered structures with many artifacts,
such as ruptures or layer splitting and merging. Previous ridge-surface
extraction methods fail to extract the desired 2D manifold for such
challenging data. We have therefore developed a novel method to extract 2D
manifolds. The proposed method uses a local fast marching scheme in
combination with a separation of the region covered by fast marching into two
sub-regions. The 2D manifold of interest is then extracted as the surface
separating the two sub-regions. The local scheme can be applied for both
automatic propagation as well as interactive analysis. We demonstrate the
applicability and robustness of our method on both artificial data as well as
real-world data including folded silver and papyrus sheets.

###### keywords:

Ridge surface, crease surface, 2D manifold extraction, fast marching, virtual
unfolding, historical documents.

K.6.1Management of Computing and Information SystemsProject and People
ManagementLife Cycle; K.7.mThe Computing ProfessionMiscellaneousEthics
Extraction of three pieces of a 2000-year-old, folded silver sheet from the
Museum _Varusschlacht im Osnabrücker Land_ , Germany. Left: Volume rendering
of the folded silver sheet. Middle: Three 2D manifolds extracted from the
silver sheet. Right: The corresponding flattened silver pieces rendered with
thin-volume rendering [16]. Colors indicate correspondence.

Introduction

With the advances in imaging techniques in recent decades, a new era in the
analysis of historical written and decorated documents has begun. Instead of
physically investigating the often very fragile documents, which bears the
risk of damaging or even destroying them, the documents are instead imaged and
virtually analyzed [30, 34, 29, 17, 36, 8, 26, 41, 42, 9]. One important step
in this process is the reconstruction of the geometry of the thin-layer
document described by, for example, the medial surface of the writing
material. For rolled or simply folded documents, the geometry reconstruction
can often be done in 2D cross-sections of the image data, resulting in
1-dimensional (1D) contours for each cross-section. The contours of
consecutive cross-sections can then be concatenated to form a 2D manifold
describing the document’s structure [30, 34, 36, 8]. However, many historical
documents show a more complex folding structure, that is, they are folded
along at least two different directions [8, 41, 9]. Such documents are usually
either only reconstructed piece-wise [41], that is, not in their entirety, or
the reconstructed documents might be subject to strong deformations [8]. A
notable exception is the recent work on unfolding letters [9]. But even for
rolled documents, the cross-sectional approach does not always work, for
example in case of deformations due to buckling perpendicular to the folding
direction [29, 28].

In this work, we address the limitations posed by previous analysis methods
for folded historical documents by tackling the reconstruction of the 2D
manifold directly in 3D. Since we know that the historical documents which
motivated this work were open, orientable 2D manifolds before being folded, we
are particularly interested in methods that allow the extraction of 2D
manifolds with such properties. For this, we carefully analyzed previous
methods for ridge-surface extraction for their suitability regarding the
geometric reconstruction of 2D manifolds of historical documents. Among the
previous works, we identified the one by Algarni and Sundaramoorthi [1] as the
most promising one for our purpose. Similarly to their method, we employ fast
marching to compute optimal paths. However, as our modifications to the
algorithm of Algarni and Sundaramoorthi are extensive, we consider our
proposed method a novel method for extracting 2D manifolds. In fact, all steps
of their algorithm have been either modified or replaced:

  * •

The overall global approach was replaced by a local one that besides enabling
the reconstruction of 2D manifolds from strongly curved and folded structures
also improves the performance.

  * •

The algorithm for ridge extraction was improved to compute the optimum instead
of an approximation.

  * •

For the calculation of the ridge surface, a novel idea and implementation
details are presented.

  * •

Our method does not use cubical complexes and thus is easier to implement and
modify to yield subvoxel precision.

  * •

The surface generated by our algorithm is guaranteed to be an orientable 2D
manifold. Their method only guarantees the creation of a complex which is
homotopy-equivalent to a point, that is, it is retractable. More details can
be found in Sect. S1 of the Supplementary Material.

In summary, our main contributions can be described as follows:

  * •

A novel method for the extraction of 2D manifolds from (µ)CT scans of thin and
highly curved documents that is less prone to noise than previous methods
based on differential operators like the Hessian. It also does not make use of
scale-space approaches that may distort the data and may lead to topological
changes in order to be robust against noise.

  * •

The resulting surface is guaranteed to be an orientable 2D manifold without
any constraints on the input data or parameters, and small holes as well as
small surface patches are less likely to occur. Both properties are important
for historical documents to make their full contents accessible.

  * •

Our local, iterative scheme enables straightforward integration of user input,
allowing the extraction of 2D manifolds from even very challenging data for
which a full automation is not yet possible. At the same time, it supports
easy fine-tuning of the parameters.

  * •

Our novel methods are rigorously mathematically defined in a continuous
setting for easy implementation and modification.

The rest of the paper is structured as follows. Previous methods are described
in Sect. 1. The proposed algorithmic pipeline is then presented in Sect. 2,
followed by a short description of the most important implementation details
in Sect. 3. In Sect. 4, we present the results for several artificial and real
data sets and analyze the robustness of the tested methods w.r.t. noise and
other challenging image characteristics. Our findings are then discussed in
Sect. 5 and conclusions as well as ideas for future work are presented in
Sect. 6.

##  1 Related Work

Our work is motivated by the analysis of 3D image data of strongly curved and
folded, written and decorated, historical documents. In the volumetric image
data, the documents are usually represented as thin-layer structures, which
often resemble ridges or valleys that are also called creases. Since they
represent important features for the analysis of 3D data, they have raised a
great interest in many visualization applications [21, 33, 32, 31, 13, 20].

Early definitions of such extremal structures in the context of image analysis
were given in the works by Haralick [14], Koenderink and van Doorn [23], and
Eberly et al. [11]. A good summary of these early definitions is given by
Schultz et al. [35]. In their work, they also showed that when extracting
crease surfaces, which are 2D manifolds in 3D space, it is important to
consider degeneracies of the Hessian. Moreover, they presented an efficient
algorithm for the extraction of crease surfaces that outperformed previous
methods, both in terms of accuracy and speed. Based on the work by Schultz et
al. [35], Barakat et al. first presented an efficient rendering approach for
crease surfaces [7]. Later, they proposed an efficient meshing scheme
intertwining ridge point sampling and meshing to yield accurate approximations
of crease surfaces [6]. For the efficient and more robust sampling of
locations of points on the crease surface, they make use of the scale-space
particle sampling approach presented by Kindlmann et al. [22].

In contrast to the methods described above, which all make use of the Hessian,
Algarni and Sundaramoorthi [1] developed an approach based on optimal paths
[37] in combination with a topological reduction of a volume to a retractable
surface. Their method requires a single seed point, which can often be
computed automatically. In extensive tests, they showed that their approach is
more accurate and produces less holes than the one by Schultz et al. [35].

For completeness, Poisson reconstruction [18] and medial surface extraction
[3] methods should also be mentioned, since both of these approaches allow the
generation of surfaces. However, just like the isosurface computation, the
Poisson reconstruction method yields only closed surfaces, or in the case of
screened Poisson reconstruction [19], also surfaces reaching the boundary of
the domain itself. In our data, the boundary of our surface lies in the domain
and is not known beforehand. Medial surface methods, on the other hand, do not
guarantee to produce 2D manifolds. Since the surfaces we want to extract are
usually open 2D manifolds, neither of the two approaches is applicable to our
use-case and data.

Apart from the more methodologically oriented works described above, recently,
a processing pipeline for the automatic unfolding of folded letters was
presented [9]. Their pipeline combines many processing steps in a
sophisticated way but their surface reconstruction step does not guarantee a
2D manifold. They also make some assumptions, for example, that the thickness
of the writing material should not vary too much, which cannot be guaranteed
in our application, in particular in case of papyrus but also folded silver
sheets.

We closely investigated and assessed the methods mentioned above for their
suitability to extract 2D manifolds from 3D image data of historical, written
documents. As result of this assessment, we develop a new algorithmic pipeline
that is described in the next section and compared to other methods in Sect.
4.

##  2 Algorithmic Pipeline

\begin{overpic}[width=433.62pt]{figures/pipeline_small.jpg}\sf\small\put(1.0,1.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}A}}
\put(8.8,10.5){\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$s$}
\put(9.5,9.0){\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>\put(26.0,1.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}B}}
\put(51.0,1.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}C}}
\put(59.5,7.0){\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$I_{\tilde{t}}$}
\put(57.0,8.5){\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>\put(57.0,5.5){\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>\put(76.0,1.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}D}}
\put(84.5,8.0){\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$S_{\tilde{t},1}$}
\put(81.5,8.7){\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>\put(84.5,5.7){\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$S_{\tilde{t},2}$}
\put(81.5,5.5){\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>\end{overpic} Fig. 1: Illustration of the proposed method for a single seed
point $s$. (A) Image data showing a smooth thin-layer structure with some
noise and the seed point $s$ for fast marching. (B) Time field $t$ resulting
from fast marching starting in $s$. The lines from near the seed point in both
directions illustrate the flooding sets to identify starting points for the
labelling of the isocontour of the time field. (C) The length field of the
computed optimal paths that was used for the labeling and the labeled
isocontour $I_{\tilde{t}}$, which here is equivalent to $S_{\tilde{t}}$. (D)
Propagation of the isocontour labels $S_{\tilde{t},1}$ and $S_{\tilde{t},2}$
to its enclosed region, the interface of which represents part of the desired
manifold $M_{R}$.

In this section, we motivate and describe in detail each processing step of
the proposed algorithm. Note that this description is for a continuous setting
and that we use the term _surface_ for a set of points that are almost
everywhere equivalent to a 2D manifold. Details regarding the implementation
in the discrete setting are provided in Sect. 3.

###  2.1 Problem Statement

Let $\Omega\subseteq\mathbb{R}^{3}$ be a domain containing a 2D orientable
manifold $M_{R}$. Let also $p:\Omega\to[0,1]$ be a probability field, where
$p(x)$ tells us how likely it is for $M_{R}$ to contain point $x$. Given the
probability field, we try to reconstruct $M_{R}$. For a given manifold
$M_{R}$, the best possible probability field would be the characteristic
function

$p_{\textit{opt}}(x):=\begin{cases}1&\text{if }x\in M_{R},\\\
0&\text{otherwise.}\end{cases}$

In this case, $M_{R}$ could be directly extracted from $p_{\textit{opt}}$. In
real-world applications, however, such a probability field is not realistic.
Instead, the probability field $p$ can be quite blurry, with large values near
$M_{R}$, and may also contain a lot of inhomogeneous noise. For all examples
used in this work, the probability field $p$ is derived from the volumetric
input data by clamping the image values to a user-defined range with linear
mapping inside the range. For the real-world data, the volumetric input data
are given by the intensity values of the µCT scans.

###  2.2 Overview of the Algorithmic Pipeline

The algorithmic pipeline can now be summarized as follows. Let $\mathcal{S}$
be a set of seed points for which we assume that they are on or close to the
orientable 2D manifold $M_{R}$ that we are looking for. Also, let $M_{R^{*}}$
be an empty surface. Then, for each $s\in\mathcal{S}$, we do the following:

  1. 1.

Compute optimal paths within the probility field $p$ w.r.t. some cost function
$c_{\textit{fm}}$ (see Sect. 2.3) with a user-defined maximum distance
starting from $s$, generating a connected region $I_{\leq\tilde{t}}$ around
the seed point, using fast marching in $p$. Fast marching determines for each
point inside this region a distance and a time value.

  2. 2.

Compute isosurface $I_{\tilde{t}}$ of the generated time field $t$ with the
smallest time value $\tilde{t}$ such that $I_{\tilde{t}}$ contains a point
with distance to $s$ larger than a user-defined value (Fig. 1B). Select the
connected component $S_{\tilde{t}}$ of $I_{\tilde{t}}$ that fully contains
$I_{\leq\tilde{t}}$.

  3. 3.

Label $S_{\tilde{t}}$ with two labels, $S_{\tilde{t},1}$ and
$S_{\tilde{t},2}$, such that the 1D contour on $S_{\tilde{t}}$ separating the
two labels lies in $M_{R}$. For this, two poles are generated on
$S_{\tilde{t}}$ that are flooded w.r.t. the generated distance mapped to
$S_{\tilde{t}}$ (Fig. 1C).

  4. 4.

Propagate the two labels starting from $S_{\tilde{t}}$ to the inside of region
$I_{\leq\tilde{t}}$, resulting in a labelling of $I_{\leq\tilde{t}}$ with two
distinct labels (Fig. 1D).

  5. 5.

Extract the surface patch $M$ separating the previously generated labels for
$I_{\leq\tilde{t}}$.

  6. 6.

Merge the surface patch $M$ into $M_{R^{*}}$.

In our implementation (Sect. 3), we do not merge the local surface patches but
instead the local, labeled regions, taking into account the time values. Also,
we implemented an automatic, iterative seed point generation step to automate
the extraction of $M_{R}$. Hence, in general, only an initial single seed
point is required.

A suitable approach to reconstruct $M_{R}$ from the probability field $p$ is
the extraction of its 2D ridge structures (ridge surfaces). In the following,
we turn to the problem of extracting ridge surfaces from $p$. To define a
ridge surface, we first define ridge points akin to the height definition of
Eberly et al. [11].

######  . Definition Ridge Point

Let $(M,g)$ be an $n$-dimensional Riemannian manifold with
$M\subset\mathbb{R}^{n}$ denoting the differentiable manifold and $g$ the
Riemanian metric, which induces the local norm on the tangent bundles of the
manifold. Further, let $f:M\to\mathbb{R}$ be a continuous function. A point
$x\in M$ is called a ridge point of type $n-l$ for some $l\in\mathbb{N},l\leq
n$, if there exists an $\epsilon>0$, an exponential map $\exp_{x}:E_{n}\to M$
with $E_{n}:=\\{y\in T_{x}M\mid\lvert y\rvert_{x}<\epsilon\\}$, and a linear
$l$-dimensional subspace $T_{l}$ of $T_{x}M$ such that for all $y\in
E_{d}:=E_{n}\cap T_{l}$ it holds that
$f(x)>f({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\exp_{x}(y)})$.

A ridge point of type $n-l$ is thus a point on the n-dimensional manifold $M$
that is a local maximum in $l$ independent directions. And the set of ridge
points of type $n-l$ forms an $(n-l)$-dimensional manifold. We call a 2D
manifold of ridge points a ridge surface. For the problem stated in Sect. 2.1,
we have $n=3$, $l=1$, and the ridge points we are interested in are of type 2.

Assuming that a surface that fits well to the probability field coincides with
a prominent ridge surface of that field, methods for calculating ridge
surfaces by computing the Hessian for each point can be used. However, since
the Hessian is very sensitive to noise, another approach was utilized by
Algarni and Sundaramoorthi [1] based on the fast marching method to find
optimal paths through the probability field. The idea of their approach is
that the ridge points of type $1$ of some function $d$ on the fast marching
front are points of the surface $M_{R}$ we are looking for. Here, the value of
$d$ for a point $x$ is defined as the length of an optimal path from $x$ to
some start point, where the optimal path was generated by the fast marching
method. Thus, instead of calculating ridge points of type $2$ in
$\mathbb{R}^{3}$, the ridge points of type $1$ on a surface $S$ in
$\mathbb{R}^{3}$ are calculated.

###  2.3 Optimal Path Computation by Fast Marching

We use the extended real numbers
$\overline{\mathbb{R}}:=\mathbb{R}\cup\\{-\infty,+\infty\\}$ and write
$f(\infty)=y$ for any set $Y$ and any function $f:\mathbb{R}\to Y$ if the
limit $\lim_{u\to\infty}f(u)=y$ exists. Furthermore, for some closed volume
$C\subset\mathbb{R}^{3}$, let $n(x)$ be the normal of the surface $\delta C$
at point $x\in\delta C$.

######  . Definition Eikonal Equation [39]

Let $\Omega\subseteq\mathbb{R}^{3}$ be an open set,
$\mathcal{S}\subseteq\Omega$ the set of starting points,
$t:\Omega\to\mathbb{R}$ and $c_{\textit{fm}}:\Omega\to[0,\infty]$. The
classical eikonal equation is known as

$\displaystyle\lvert\nabla t(x)\rvert$ $\displaystyle=c_{\textit{fm}}(x)\quad$
$\displaystyle\forall x\in\Omega,$ $\displaystyle t(x)$ $\displaystyle=0$
$\displaystyle\forall x\in\mathcal{S}.$

The fast marching method [37] finds an approximate solution $t$ to this
equation given a cost function $c_{\textit{fm}}$. As all variations of the
algorithm are using the causality property, better known as the upwind
condition [38], the solution is calculated from the starting points
$s\in\mathcal{S}$ outwards into the domain $\Omega$, creating one or multiple
fronts of points for which a possible value $t(x)$ is calculated.

The time field $t$ can then be seen as the time the fast marching algorithm
took to march from a point in $\mathcal{S}$ to point $x$, which is the minimal
amount of time required, given the cost function. Note that as $t$ can be
interpreted as the time of arrival for an optimal path, which starts from $s$
and is created by fast marching, the maximal integral curves of $\nabla t$
represent the optimal paths themselves (see, e.g., [24] for a definition of
maximal integral curves). We assume the length of an optimal path to be given
by $d(x)$ and we say that the distance between $x$ and $\mathcal{S}$ is
$d(x)$.

Given a probability field $p:\Omega\to[0,1]\subset\mathbb{R}$, we solve the
eikonal equation (Sect. 2.3) with cost function
$c_{\textit{fm}}:\Omega\to[0,\infty]$,

$c_{\textit{fm}}(x):=\begin{cases}\frac{1-p(x)}{p(x)}&\text{if }p(x)\neq 0,\\\
\infty&\text{otherwise.}\end{cases}$

Note that the starting points $\mathcal{S}$ should lie in the manifold
$M_{R}$.

###  2.4 Bundling of Optimal Paths

If $\mathcal{S}$ contains exactly one point $s$, then the region
$I_{\leq{\tilde{t}}}:=\\{x\in\Omega\mid t(x)\leq{\tilde{t}}\\}$ is connected
for any ${\tilde{t}}\in[0,\infty]$. However, the isosurface
$I_{\tilde{t}}:=\\{x\in\Omega\mid t(x)={\tilde{t}}\\}$ is not necessarily
connected. Let $S_{\tilde{t}}$ be the connected component of the isosurface
$I_{\tilde{t}}$ which encloses $I_{\leq{\tilde{t}}}$, and let
$S_{\leq{\tilde{t}}}$ be the closed volume of $S_{\tilde{t}}$. Thus, we have
$I_{\leq{\tilde{t}}}\subseteq S_{\leq{\tilde{t}}}$.

For our case, we are mostly interested in isosurfaces of $t$ and integral
curves of $\nabla t$. An isosurface $I_{\tilde{t}}$ gives us the front of the
fast marching algorithm at time $\tilde{t}$ and an integral curve for some
point in $I_{\tilde{t}}$ follows the fronts $I_{t^{\prime}}$ of the fast
marching algorithm through time $t^{\prime}$. Integral curves of $\nabla t$
can be seen as an “orthogonal” element of the isosurfaces of $t$ in the sense
that the normalized tangent of the integral curve of any point $x$ with
$t(x)=\tilde{t}$ is the same as the normal of the isosurface $I_{\tilde{t}}$
at $x$.

Due to the upwind condition of the fast marching algorithm, $s$ is the only
local minimum of the time field. Together with the fact that for any point $x$
on $I_{\leq\tilde{t}}$ and any $u\in\mathbb{R}$, any integral curve
$\Gamma_{x}:T\to I_{\leq\tilde{t}}$ with $\Gamma_{x}(u)\in I_{\tilde{t}}$ must
have as tangent vector $\Gamma_{x}^{\prime}(u)=\nabla t(x)$ a vector with a
direction pointing to the outside of $I_{\leq\tilde{t}}$ (again because of the
upwind condition), a maximal integral curve can only start at $s$, that is,
$T=[-\infty,u_{\textit{max}}],u_{\textit{max}}\in\overline{\mathbb{R}}$ and
$\Gamma_{x}(-\infty)=s$. If the integral curve ends in a local maximum, the
same reasoning can be applied to show that $u_{\textit{max}}=\infty$ is the
case. Otherwise $\Gamma_{x}(u_{\textit{max}})\in I_{\tilde{t}}$, that is, it
can only end on $I_{\tilde{t}}$, with $u_{\textit{max}}\in\mathbb{R}$.

The other direction also holds: No integral curve on $I_{\leq{\tilde{t}}}$ can
pass through a local maximum or a point $x$ on $I_{\tilde{t}}$, as the
gradient at a local maximum is $0$ and the gradient at $x$ on $I_{\tilde{t}}$
is the normal of $I_{\tilde{t}}$ at $x$.

With these properties in mind, we are able to define a mapping from
$I_{\tilde{t}}$ to $I_{\leq{\tilde{t}}}$, that also defines a non-overlapping
partition. A new partition, which we will call flooding sets, are created by
slight modifications of the mapping induced partition. These flooding sets
will allow us to project labels from $S_{\tilde{t}}$ to $S_{\leq\tilde{t}}$
and vice versa.

######  . Definition Integral Curve End Mapping

Let $\Omega\subseteq\mathbb{R}^{3}$ be open, $C\subset\Omega$ be a smooth,
compact manifold with a boundary, and $t:\Omega\to\mathbb{R}$ a differentiable
function. We define a mapping on $C$:

$\zeta(x):=\Gamma_{x}(u_{\textit{max}})\ \forall x\in C$

with
$\Gamma_{x}:[u_{\textit{min}},u_{\textit{max}}]\subseteq\overline{\mathbb{R}}\to
C$ being the maximal integral curve of $\nabla t$.

The integral curve for each point $x$ is unique and well defined, because the
solution $t$ of the eikonal equation is differentiable and bounded on all
compact sets. The uniqueness and existence is a classical result from the
theory of ordinary differentiable equations [4, p. 305]. In consequence, the
mapping $\zeta$ is also well defined.

######  2.1. Theorem Integral Curve End Mapping Induced Relation

Let $\Omega$, $t$, $C$ and $\zeta$ be the same as in Sect. 2.4. Also, let

$n(x)=\frac{\nabla t(x)}{\lvert\nabla t(x)\rvert}$

for all $x\in\delta C$ with $\nabla t(x)\neq 0$. Then the following statements
hold:

  1. 1.

$\zeta$ induces an equivalence relation,

  2. 2.

the surjective image of $\zeta$ is $\delta C\cup E$, where $E$ are the local
maxima in $C$ regarding $t$, and

  3. 3.

for any $x\in\delta C\cup E$, we have $\zeta(x)=x$.

For the proof, see Sect. S4 of the Supplementary Material.

The partition defined through $\sim_{\zeta}$ depends on $C$ since the local
maxima bundle multiple sets together into one.

This can be seen in the image on the right, in which all integral curves on
the left side end in the local maximum and, thus, belong to one equivalence
class as long as the local maximum is contained in $C$.

One can observe that not only $I_{\leq\tilde{t}}$ but also $S_{\leq\tilde{t}}$
for any $\tilde{t}\in\mathbb{R}$ is a valid choice for $C$. Thus, given a
labeling of $\delta S_{\leq\tilde{t}}\cup E=S_{\tilde{t}}\cup E$, we are able
to label $S_{\leq\tilde{t}}$ by giving all elements in the equivalence class
the same label. However, our aim is to create a partition of
$S_{\leq\tilde{t}}$ with representatives only in $S_{\tilde{t}}$. To achieve
this, we can combine equivalence classes with an element in $E$ with another
class with an element in $S_{\tilde{t}}$.

Our solution is to use a flooding algorithm and map integral curves which end
in a local maximum to another nearby integral curve which ends in
$S_{\tilde{t}}$. We call the partition obtained from the flooding algorithm
starting at point $x$ the flooding set of $x$. If no local maximum exists in
$S_{\leq\tilde{t}}$, then this is equal to the inverse image $\zeta^{-1}(x)$.

###  2.5 Ridge Curve Extraction

Given a surface $M_{R}$ and a discrete probability field

$p_{\textit{disc}}(x):=\begin{cases}y_{1}&\text{if }x\in M_{R},\\\
y_{2}&\text{otherwise,}\end{cases}$

with $y_{1}>y_{2}$ and $y_{1},y_{2}\in[0,1]$, Algarni and Sundaramoorthi [1]
proved that points of $M_{R}$ must be ridge points in $S_{\tilde{t}}$ w.r.t.
the distance field $d$. Furthermore, we assume $M_{R}$ to intersect
$S_{\tilde{t}}$ such that $S_{\tilde{t}}\setminus M_{R}$ are two connected
components, $S_{\tilde{t},1}$ and $S_{\tilde{t},2}$, and $S_{\tilde{t}}\cap
M_{R}$ is a closed contour on $S_{\tilde{t}}$ (see Fig. 1). This condition
allows us to describe a robust labeling algorithm for $S_{\tilde{t}}$, and we
say that the surface $M_{R}$ intersects $S_{\tilde{t}}$ fully. Thus, we search
$S_{\tilde{t}}$ for closed contours containing ridge points w.r.t. the
distance field $d$. Such a closed contour will simply be called a (1D) ridge
in $S_{\tilde{t}}$.

A classic hierarchical watershed algorithm on the distance field of the fast
marching evaluated on $S_{\tilde{t}}$ can be used to find such ridges.
However, we found that in practice, a flooding algorithm creates better and
more stable results. Given a closed surface $S_{\tilde{t}}\subset\Omega$ and
two starting points $u_{1},u_{2}\in S_{\tilde{t}}$ on the surface, we split
the surface into two parts $S_{\tilde{t},1}$ and $S_{\tilde{t},2}$. The set
$S_{\tilde{t},1}$ contains all points $x$ on $S_{\tilde{t}}$ for which the
minimal cost of a path from $x$ to $u_{1}$ is less than a minimal cost path
from $x$ to $u_{2}$. The set $S_{\tilde{t},2}$ is defined similarly. The cost
of a path $\gamma:[0,1]\to S_{\tilde{t}}$ from $x$ to $y$, that is
$\gamma(0)=x$ and $\gamma(1)=y$ holds, is defined by

$c_{\textit{rc}}(\gamma):=\max_{z\in[0,1]}d(z).$

To find the two points $u_{1}$ and $u_{2}$ on the surface, we make use of the
flooding sets of $t$ as well as the structure tensor of the 3D probability
field $p$ at the given seedpoint. Let $A$ be the structure tensor at the
seedpoint $s$ for which the fast marching created the given closed surface
$S_{\tilde{t}}$ as a part of an isosurface. Let $v_{1}$ be the eigenvector of
the structure tensor with the largest eigenvalue. We set
$u_{1}^{\prime}:=s+\lambda\frac{v_{1}}{\lvert v_{1}\rvert}$ and
$u_{2}^{\prime}:=s-\lambda\frac{v_{1}}{\lvert v_{1}\rvert}$ for a user-defined
constant $\lambda$. In practice, $\lambda$ depends on the thickness of the
layers, the distance between layers of the probability field and the curvature
of the manifold. If $u_{1}^{\prime}$ is on or outside of the surface
$S_{\tilde{t}}$, we define $u_{1}$ to be the point on the surface which is
between $s$ and $u_{1}^{\prime}$. Otherwise, if $u_{1}^{\prime}$ is inside of
the surface $S_{\tilde{t}}$, we start a flooding algorithm at
$u_{1}^{\prime}$, which prefers the highest possible time values $t$. The
first point of $S_{\tilde{t}}$ which is flooded by the algorithm is defined as
the point $u_{1}$. In other words, $u_{1}\in S_{\tilde{t}}\cap
F(u_{1}^{\prime})$, where $F(u_{1}^{\prime})$ is the flooding set of
$u_{1}^{\prime}$. We define $u_{2}$ analogously.

###  2.6 Extraction and Merging of Surface Patches

\begin{overpic}[width=433.62pt]{figures/iterativeapproach-
small.jpg}\sf\small\put(12.0,-1.2){A} \put(37.0,-1.2){B} \put(62.0,-1.2){C}
\put(87.0,-1.2){D} \end{overpic} Fig. 2: Comparison of global vs. local,
iterative approach. Top row: Time fields of the fast marching method are shown
using the viridis colormap. Bottom row: Label fields resulting from our
approach. From the interface between the labels, the manifold to be
reconstructed is computed. (A) Illustration of the global approach with a
single seed point, as proposed by Algarni and Sundaramoorthi [1]. The 2D
manifold of the strongly curved and folded structure cannot be fully
reconstructed. (B) A single local step starting from the same seed point as in
A. (C) Two more seed points were added. (D) Full reconstruction of the
manifold after four additional seed points have been added.

Given a closed surface and a partition of the surface into two pieces, one may
also partition the inside of the surface by following the flooding sets. Let
$\Omega\subseteq\mathbb{R}$ be a domain, $S_{\tilde{t}}\subseteq\Omega$ a
closed surface and
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}S_{\tilde{t},1},S_{\tilde{t},2}}$
a partition of the surface into two parts. Then, we have a partition of the
inner volume given by the flooding sets, their representatives and the
labeling of the surface.

Since $S_{\tilde{t},1}$ and $S_{\tilde{t},2}$ are retractable and as union
homeomorphic to a sphere, each of them is homeomorphic to a disk. As the front
of integral curves can only change their topology through time by the
existence of a local maximum – ‘splitting’ the front – and flooding sets
always relate the integral curves ending at the local maximum to a neighboring
integral curve, the topology of the front of flooding sets does not change
until it contracts to the single seed point. Thus, the partitions of the
volume $S_{\tilde{t}}$ are homeomorphic to a solid sphere. The border between
the two components defines the surface patch of interest and is guaranteed to
be a 2D manifold and in fact homeomorphic to a disc.

Our calculations are based on the assumption that the surface $M_{R}$
intersects $S_{\tilde{t}}$ fully. In consequence, the value of $\tilde{t}$
cannot be chosen too high as the front of the fast marching algorithm would
otherwise grow further outside beyond the surface itself. Thus, the surface
patch generated by $S_{\tilde{t}}$ is only a part of the whole 2D manifold
$M_{R}$ we try to reconstruct. Instead of extracting multiple surface patches
and merging them, we merge the two components of $S_{\leq\tilde{t}}$ for all
seed points (Fig. 2). The manifold $M_{R}$ is then given by the interface
between these two components.

Let $\mathcal{S}:=\\{s_{0},s_{1},\ldots,s_{n}\\}\subset\mathbb{R}^{3}$,
$n\in\mathbb{N}$ be a finite list of seed points. For each seed point
$s_{i}\in\mathcal{S}$, we obtain for all points $x\in S_{\leq{\tilde{t}}}$ a
time $t_{i}(x)$ and a label $l_{i}(x)\in\\{2i+1,2i+2\\}$ that encodes to which
seed point and component, $S_{\tilde{t},1}$ or $S_{\tilde{t},2}$, the point
belongs. For any other point $y\notin S_{\leq{\tilde{t}}}$, we define
$t_{i}(y):=\infty$ and $l_{i}(y):=0$. We also define the dominant seed point
of a point $x\in\mathbb{R}^{3}$ to be the seed point with the smallest arrival
time

$\operatorname*{dsp}(x):=\operatorname*{arg\,min}_{i=0,\ldots,n}\ t_{i}(x).$

In the case that there is more than one seed point with smallest arrival time,
any one can be chosen. The resulting time and label fields are then defined to
coincide with the dominating seed point

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}t_{\Sigma}}(x):=t_{\operatorname*{dsp}(x)}(x)\text{
and
}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}l_{\Sigma}}(x):=l_{\operatorname*{dsp}(x)}(x).$

Note that the choice of the labels allows us to find the dominating seed point
through the label value itself, as it holds true that

$\operatorname*{dsp}(x)=\left\lfloor\frac{l_{\Sigma}(x)-1}{2}\right\rfloor.$

Each seed point has its own two labels and we also have a label for the
background. To generate a merged surface, we must cast these $2n+2$ labels to
only two (the background label is excluded here). For this, we want to define
an equivalence relation $\sim_{2}$ between labels. Let $\mu$ be a measure (see
for example [5, p. 41-42]) and $i,j\in\\{0,\ldots,n\\}$ arbitrarily chosen but
fixed. Furthermore let $L_{pq}:=l_{i}^{-1}(2i+p)\cap l_{j}^{-1}(2j+q)$ for any
$p,q\in\\{1,2\\}$ and
$c_{rel}:=\mu(L_{11})+\mu(L_{22})-\mu(L_{12})-\mu(L_{21})$. We set the
relations

$\begin{cases}2i+1\sim_{2}2j+1\text{ and }2i+2\sim_{2}2j+2&\text{if
}c_{rel}>0,\\\ 2i+1\sim_{2}2j+2\text{ and }2i+2\sim_{2}2j+1&\text{if
}c_{rel}<0.\end{cases}$

Note that we do not set any relations if $c_{rel}=0$. The relation $\sim_{2}$
is reflexive as $L_{11}>0$, $L_{22}>0$, $L_{12}=0$ and $L_{21}=0$ holds for
$i=j$, and symmetric by definition. We furthermore enforce the transitivity,
that is, if $k_{1}\sim_{2}k_{2}$ and $k_{2}\sim_{2}k_{3}$, then we also set
$k_{1}\sim_{2}k_{3}$ for all $k_{1},k_{2},k_{3}\in\\{1,\ldots,2n+2\\}$. Thus,
the relation is an equivalence relation.

The seed points $i$ and $j$ have an overlap if and only if
$\mu(L_{11})+\mu(L_{12})+\mu(L_{21})+\mu(L_{22})>0$. In the edge case that the
seed point $i$ and $j$ have an overlap but $c_{rel}$ is $0$ we may assume the
patch generation to contain errors. Otherwise if all seed points are connected
via overlap and the patch generation does not contain an error, the relation
$\sim_{2}$ has at most two equivalence classes. If the generated surface is
non-orientable, the relation has exactly one class. As we are interested in
surfaces which can be unfolded, we may assume the relation to have two
equivalence classes. Thus, we are able to switch the labels $2i+1$ and $2i+2$
for all $i\in\\{0,\ldots,n\\}$ such that $c_{rel}$ is always positive. This
simplifies our relation $\sim_{2}$ as now it is completely described by

$2i+1\sim_{2}2j+1\land 2i+2\sim_{2}2j+2\,\,\forall i,j\in\\{0,\ldots,n\\},\
\text{or alternatively}$ $\forall
k_{1},k_{2}\in\\{1,\ldots,2n+2\\}:k_{1}\sim_{2}k_{2}\Leftrightarrow
k_{1}\equiv k_{2}\pmod{2}\ .$

Together with the background label, we now have $3$ labels but are only
interested in the shared border of the two labels $1$ and $2$, disregarding
the border with the background label.

###  2.7 Seed Point Neighbors

Merging patches by the use of a label field sidesteps many specialized
operations on surfaces, however it creates another problem: the creation of
unwanted surfaces. The image on the right gives an example which showcases a
label field with two labels and a manifold (black curve) which could have
created such a label field. However, one can see that the border between the
labels not only includes the manifold itself (black curve), but that there is
a second border between the same labels.

To solve this, we make use of the fact that the border between the two labels
shall only describe a surface if the seed points responsible for the labels
are neighbors. The intuition is that a neighbor relationship between seed
points describes which seed points are near each other on the manifold such
that labels from another layer do not interfere with the current layer.

Let $\mathcal{S}:=\\{s_{0},s_{1},\ldots,s_{n}\\}$, $n\in\mathbb{N}$ again be a
finite set of seed points. We define a neighborhood graph $G$ with vertices
$\mathcal{S}$. The points $s_{i}\in\mathcal{S}$ and $s_{j}\in\mathcal{S}$ are
adjacent to each other in $G$ if and only if

$\displaystyle\delta L_{s_{i}}\cap\delta L_{s_{j}}\neq\emptyset\text{ with}$
$\displaystyle\delta
L_{s_{i}}:=\\{x\in\mathbb{R}^{3}\mid\exists\epsilon>0,y_{1},y_{2}\in
U_{\epsilon}(x)\text{ such that }$ $\displaystyle l_{s_{i}}(y_{1})=2i+1\text{
and }l_{s_{i}}(y_{2})=2i+2\\}$

and $\delta L_{s_{j}}$ analogously. Then $\delta L_{s_{i}}$ is the set of
border points between the two labels of the seed point $s_{i}$ and two seed
points are connected exactly when they have a border point in common.

##  3 Implementation

Fast Marching.

We implemented the classical fast marching algorithm [37]. If the upwind
condition is broken in 3D, we reduce the dimension of the update step, similar
to the work by Deschamps and Cohen [10] to guarantee that the upwind condition
holds. As the values of the solution $t$ can become huge, we update the values
depending on the mean of all neighboring values. This is numerically more
stable. The Euclidean length $d$ of the optimal paths can be calculated at the
same time as $t$. Fast marching is done on $t$ but the same update steps for
$t$ are also applied to $d$, with the only difference that $d$ has as cost
function the constant function $c_{d}(x)=1$. Our stopping criteria for the
fast marching front is specified by a user-defined upper bound for the path
length $d$. The fast marching is stopped when the first value higher than that
is encountered. The choice of $d$ is related to the average thickness of the
thin-layer structure that we are interested in. As a rule of thumb, it should
be at least 5 to 10 times the thickness. However, the choice of $d$ mainly
influences the performance in terms of run time, not quality. When $d$ becomes
too large, run time might increase drastically.

Flooding.

In Sect. 2.5, we use a standard flooding algorithm over voxels. Given a set of
voxels $V$, the algorithm greedily takes the voxel neighboring a voxel in $V$
with the highest time value and inserts it into $V$. This is repeated until a
voxel is found which is at the border of $S_{\leq\tilde{t}}$. The same
algorithm is used with faces instead of voxels to label $S_{\tilde{t}}$ given
two faces on $S_{\tilde{t}}$ to grow out from. Instead of the time values, the
distance values generated by the fast marching method are used. Flooding could
also be used to find the labels of the voxels in $S_{\leq{\tilde{t}}}$ in
Sect. 2.6, however, it is faster to go the other direction and propagate the
labels from the outside to the inside, that is, starting from $S_{\tilde{t}}$,
using the time values again. To do this, we insert all voxels at the border,
which are already labeled, in a sorted set $V$, take the voxel with the
highest time value out of the sorted set and insert all its neighbors which
are not yet labeled into $V$. This is repeated until $V$ is empty and
therefore all voxels in $S_{\leq{\tilde{t}}}$ are labeled.

Automation.

We implemented a simple, iterative approach to automatically add further seed
points. Given a surface and the seed points which were used to construct the
surface, we want to generate new seed points. This could be done by generating
points on the boundary of the surface, but small inconsistencies of the ridge
curve extraction might occur. However, these inconsistencies tend to vanish
when moving only a slight distance away from the surface boundary. Therefore,
a more stable method is to set seed points some distance away from the
boundary of the surface. Furthermore, if an already generated seed point is
near the boundary of the surface, the fast marching algorithm did not prefer
growing in the direction of the boundary. Therefore, we may assume to have
reached the end of the 2D manifold we want to extract. In summary, we can
automate the generation of seed points by finding points on the surface that
are a specific distance away from the boundary (10% of the marching distance
in our case) and are far away from all seed points (30-50% of the marching
distance is often a good value). If no such point exists, the automation stops
and we assume that the whole surface was extracted.

Iterative Construction.

Note that even though we defined the relation of different labels in Sect. 2.6
explicitly, they can be calculated iteratively by adding one seed point after
another to the set $\mathcal{S}$. Hence, we only need to save the global time
field $t$, the global label field $l$, the neighboring graph and the
constructed surface of the current set of seed points. Given a new seed point,
the time and label fields as well as the neighborhood graph are updated and
the newly created patch is merged with the surface. For this, a modified
marching cubes algorithm [15] can be used to generate triangle meshes based on
Sect. 2.6 and Sect. 2.7. The new triangles are then concatenated to the
surface. Old appendages of the surface that were created by voxels for which
the times were updated, are deleted. The modified marching cubes algorithm
[15] allows the handling of a background label. Alternatively, by changing the
labeling of a voxel from the background label to one of the other two labels,
depending which label appears the most in the neighboring voxels, the
classical marching cubes algorithm [25] can be used instead.

Neighborhood Graph.

The update of the neighborhood graph is done by checking for common border
points of the newly added seed point and the old seed points. As we work with
voxels, our border points are faces between the two labels generated from the
seed points. Given the transitivity rule of the label relation as well as our
assumption of an orientable surface, we make use of the simplified description
of label relation and switch the labels of the newly added seed point such
that the relation holds. Furthermore, as we discretize the volume in voxels,
we work with a finite number of elements. Thus, the sets $L_{11}$, $L_{12}$,
$L_{21}$ and $L_{22}$ are also finite. Also, as all volume elements are of
equal size, we can therefore use the cardinality of the sets themselves as the
measure $\mu$. In the marching cubes step of generating the surface, we only
generate a triangle inside a cube if the vertices of the cube contain both
labels $1$ and $2$ and only if the dominating seed points of the vertices are
connected by the graph $G$.

\begin{overpic}[width=411.93767pt]{figures/simple_small.jpg}\sf\small\put(9.0,82.0){1
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(data)}}
\put(43.0,82.0){2
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\cite[cite]{[\@@bibref{}{Schultz:2009:CS}{}{}]})}}
\put(77.0,82.0){3
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(ours)}}
\put(-3.0,75.7){A} \put(-3.0,63.7){B} \put(-3.0,51.7){C} \put(-3.0,39.7){D}
\put(-3.0,27.7){E} \put(-3.0,15.7){F} \put(-3.0,3.7){G} \end{overpic} Fig. 3:
Reconstruction of a simple plane with different types and degrees of noise
shown in the first column. Comparison between Schultz et al.’s method [35],
shown in the second column, and ours, shown in the third column. (A) Original
scalar image showing a clear ridge. (B) Gaussian noise (10%) added to A. (C)
Gaussian filtered image from B. (D) Simplex noise (30%) added to A. (E)
Simplex noise (50%) added to A. (F) Warping of image A with strength 2. (G)
Warping of image A with strength 4.
\begin{overpic}[width=433.62pt]{figures/folded_new_small.jpg}\sf\small\put(14.0,-1.0){A
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(data)}}
\put(42.0,-1.0){B
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\cite[cite]{[\@@bibref{}{Dambrogio:2021:Unlocking}{}{}]})}}
\put(65.0,-1.0){C
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\cite[cite]{[\@@bibref{}{Schultz:2009:CS}{}{}]})}}
\put(87.0,-1.0){D
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(ours)}}
\end{overpic} Fig. 4: Reconstruction of a simply folded plane that was
slightly warped (A). Comparison between Dambrogio et al.’s method [9], shown
in B, Schultz et al.’s method [35] (C), and ours (D). In B, many holes and one
additional surface patch near the folding edge can be seen. In C, only holes
at the folding edge are present, but multiple additional surface patches near
the folding edge were created. Only our method recreates a topologically
correct surface without any additional surface patches.

##  4 Results

In this section, we present results for both artificial and real data. In
order to compare our method with one of the state-of-the-art methods in ridge
surface extraction [35], we tested both approaches on the same data. Note that
we tried to optimize the results of the method by Schultz et al.[35] by
testing different strength parameter values. Shown are always what we
considered as the best results using their method.

We also present a qualitative comparison with the letter unfolding pipeline
[9]. Note that a quantitative comparison is not possible, since the mesh
generated by the pipeline does not represent a real surface mesh but rather a
graph. Seen as a surface mesh, faces overlap and intersect, and no possible
orientation (even for small areas) can be found. Thus, a topological analysis
is impossible and we are also unable to show the mesh in two different colors.
See Sect S2 in the Supplementary Material for more information. Also, we
remark that their method takes 6 minutes to segment our simple surface
dataset. For more information, we refer to Sect. S2 in the Supplementary
Material.

We also do not include the algorithm of Algarni and Sundaramoorthi [1] in our
quantitative comparison, as their method generates many non-manifold points.
Since both their and our methods make use of fast marching to calculate
optimal paths and extract surfaces containing these optimal paths, the
generated surfaces are comparable if the right conditions for their algorithm
were met and was also done only locally (see Fig. 2). We were also able to
improve the performance of their algorithm to reduce the time needed for a
patch from over 15 minutes to a few seconds in the case of the silver sheet.
For a more detailed description, we again refer to Sect. 1 of the
Supplementary Material.

The probability fields necessary for our method were generated by normalizing
the image data to $[0,1]\subset\mathbb{R}$. The surface meshes shown in this
section should be orientable and are therefore always colored with two colors,
blue for one side of the mesh and green for the other side. This allows to
easily distinguish the two sides and to make artifacts in the reconstruction
more apparent. Note that the orientation of all meshes was made consistent as
a post-processing step after reconstruction. For the visualization of the
surface meshes, we deliberately used triangle normals in Fig. 5 and Fig. 6 to
give the user a visual feedback of the mesh resolution and quality.

All computations were carried out on a desktop PC with the following
configuration: CPU: _Intel Core i9-10980XE_ ; GPU: _Nvidia GeForce RTX 3090_ ,
24 GiB VRAM; Memory: 128 GiB DDR4.

\begin{overpic}[width=433.62pt]{figures/dragon_small.jpg}\sf\small\put(14.0,100.0){{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\cite[cite]{[\@@bibref{}{Schultz:2009:CS}{}{}]}}}
\put(50.0,100.0){{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}ours}}
\end{overpic} Fig. 5: Reconstruction of the Stanford Dragon. Left: Result from
the method by Schultz et al. [35]. Right: Our method. The green ellipses
highlight regions where the method by Schultz et al. had difficulties while
our algorithm did not show any artifacts in the reconstruction.

###  4.1 Artificial Data

Artificial data was used to provide a ground truth for the surface
reconstruction, which is not available for the real-world data.

####  4.1.1 Plane

The images used to compare Schultz et al.’s algorithm and our method are
visualized in Column 1 of Fig. 3. Sect. 3 of the Supplementary Material
contains a description on how they were created. From top to bottom the files
are called: original, g10, g10.gf100, s30, s50, ws200.gf100 and ws400.gf100.
Column 2 of Fig. 3 shows the results by applying Schultz et al.’s method,
while in column 3 our results are shown. Notice the degenerate case of A2.
This is a known behaviour of Schultz et al.’s method for synthetic data [35].
For real world data, or in our case, synthetic data with noise, such
degeneracy is highly unlikely to occur. The results of the quantitative
analysis can also be found in Sect. S3 of the Supplementary Material.

Our algorithm always creates a surface that is topologically equivalent to a
plane, while the method by Schultz et al. shows some artifacts like a triangle
soup (B2) and different degrees of holes (D2-G2), only sometimes generating a
valid 2D manifold. Our topological analysis confirms that. While our method is
robust against all types of noise, the result is for low noise images not as
precise as their method. We assume that this is mainly due to the fact that we
have not implemented sub-voxel accuracy. However for images with strong noise
and especially for warped images (which most closely represents our papyrus
data), our method creates more accurate surfaces. Also, our method usually
creates slightly larger surfaces, reaching out to the very boundary, thereby
sometimes resulting in strongly bent surfaces. Last but not least, for low
noise images Schultz et al.’s method is slightly faster. However, the
performance is highly variable and strongly depends on the strength value. As
this value has to be experimentally found (ranging from 150 to 12,000 for this
dataset), the real time cost is much higher. Our algorithm in contrast shows
similar time cost for all test data.

####  4.1.2 Simply Folded Sheet

In the second artificial test case, we folded a plane once, added Gaussian
blurring to it to simulate a thin-layer structure and slightly warped the
data. The data as well as the results using Dambragio et al.’s pipeline,
Schultz et al.’s method and ours are shown in Fig. 4. In Fig. 4A we show a
cross-section through the data. In B, the result of the pipeline by Dambragio
et al. is shown. It has a lot of holes and a small wrong patch was calculated
near the ridge. This artifact is also visible in C, which shows the result of
Schultz et al.’s method. The patches near the ridge are more pronounced. The
folding edge is also leaky, but in general the result shows a lot fewer holes
than B. Our method in D creates a surface that is topologically equivalent to
a plane, in particular, it does not contain any holes and consists of only one
patch.

####  4.1.3 Dragon

As a third artificial test data set, we used the Stanford Dragon [40].
Similarly to the previous data, the surface was sampled on a uniform grid,
here of dimensions $841\times 599\times 388$, and smoothed with a Gaussian
filter, resulting in a smooth thin-layer structure. No noise was added, but
since the dragon shows different degrees of geometrical intricacies, we found
it a good test case. Both algorithms were applied to the image data and were
able to reconstruct the original surface, as can be seen in Fig. 5, first row.
However, the algorithm by Schultz et al. created artifacts, some of which are
highlighted in the figure. Our algorithm, which used $820$ automatically
placed seed points in addition to the first manually placed one, did not
create any artifacts.

###  4.2 Real-World Data

The real data shown in this section motivated the work. Computation of the 2D
manifold from the set of seed points took 3 minutes for the papyrus data.

####  4.2.1 Complexly Folded Silver Sheet

Here, we analyzed a $\mu$CT scan with dimensions $375\times 271\times 614$ of
a 2000-year-old complexly folded silver sheet (data courtesy by Christiane
Matz and Stefan Burmeister from VARUSSCHLACHT im Osnabrücker Land gGmbH,
Museum und Park Kalkriese). A volume rendering of the folded sheet together
with three pieces that were reconstructed with the proposed method are shown
in A Local Iterative Approach for the Extraction of 2D Manifolds from Strongly
Curved and Folded Thin-Layer Structures. Note that in order to visualize the
flattened pieces, the extracted 2D manifold meshes were drastically simplified
using the algorithm by Garland and Heckbert [12], subsequently flattened into
the 2D plane utilizing quasi-isometric flattening proposed by Ambellan et al.
[2], and finally visualized with the thin-volume rendering by Herter et al.
[16]. Simplification was done for two purposes. First, it made the flattening
more robust. Second, and more importantly, in order to visualize the ornaments
that are impressed into the silver sheet, the approximating surfaces should
not follow the impressed ornaments but be smooth in these regions. Otherwise,
the ornaments would be flattened away. A simple way to create a smooth mesh
that does not follow the ornamental details is achieved using standard mesh
simplification methods.

In Fig. 6, we compare the results of the surface reconstruction using the
method by Schultz et al. and our method. Their global approach, that computes
all ridge surfaces in the image data, worked remarkably well for the silver
sheet. Nevertheless, many floating triangles and triangle patches can be
observed in Fig. 6 when zooming in to the electronic version of the paper.
Many of these triangles could be removed easily. However, as is indicated by
the color change of the surface shown in Fig. 6A,B, wrong connections occurred
resulting in a surface mesh that was not orientable anymore, which is a
requirement for the subsequent unfolding step. Identifying and removing these
connections to end up with an orientable 2D manifold might be very tedious.
With our method, we have full control over the reconstruction process that is
visually supported. One of the three pieces of the sheet was fully automically
reconstructred starting from a single seed point. The other two pieces were
semi-automatically reconstructed by letting the automatic reconstruction only
run for a few steps after which the result was visually checked. This way, we
guided the algorithm when the automatic processing made a mistake by removing
automatically placed seed points and adding manual ones if needed.

####  4.2.2 Doubly Folded Papyrus

Here, we analyze a $\mu$CT data set with dimensions $454\times 107\times 374$
of an ancient papyrus package ([26], Greek papyrus; data courtesy by Eve Menei
and Marc Etienne, Musée du Louvre, Paris). The package is folded twice in
orthogonal directions. Physically opening it will most likely result in
damaging the precious object. Therefore, a virtual unfolding (Fig. 7) is
necessary to access the writing.

\begin{overpic}[width=433.62pt]{figures/silversheet-2_small.jpg}\sf\small\put(22.0,65.0){{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\cite[cite]{[\@@bibref{}{Schultz:2009:CS}{}{}]}}}
\put(72.0,65.0){{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}ours}}
\put(22.0,32.0){A} \put(72.0,32.0){B} \put(22.0,-1.0){C} \put(72.0,-1.0){D}
\end{overpic} Fig. 6: Comparison of reconstruction results for the folded
silver sheet shown in A Local Iterative Approach for the Extraction of 2D
Manifolds from Strongly Curved and Folded Thin-Layer Structures. (A) Complete
reconstruction of the silver sheet using Schultz et al.’s method. (B)
Reconstruction of three pieces of the silver sheet using our method. (C+D)
Zoom-in to the rectangular region of A+B, respectively. Note the color flip in
A+C in the highlighted region. This shows that the extracted surface mesh is
not orientable anymore but instead contains a (wrong) connection changing the
orientation.
\begin{overpic}[width=433.62pt]{figures/papyrus_small.jpg}\sf\small\put(13.0,1.0){A}
\put(37.0,1.0){B} \put(63.0,1.0){C} \put(87.0,1.0){D} \end{overpic} Fig. 7:
Doubly folded papyrus document. (A) Volume rendering of the $\mu$CT image of
the document, also showing the position of the cross-section depicted in B.
(B) Cross-section of the $\mu$CT scan shown in A. (C) Volume rendering and
extracted 2D manifold mesh. (D) Thin-volume rendering [16] of the flattened
papyrus document. Note the writing on the papyrus document.

Since papyrus consists of two layers of papyrus fibers that are oriented
orthogonal to each other and then pressed to form the papyrus document, the
reconstruction of the document structure is particularly challenging, since
these two layers often separate from each other and have a high variability in
thickness and intensity. This is especially true for old papyrus and is also
present in the data analyzed here. Therefore, it comes as no surprise that the
pipeline by Dambragio et al. fails to generate one connected surface, missing
a substantial amount of papyrus (Fig. 8, first row). The method by Schultz et
al. fails to reconstruct a 2D manifold (Fig. 8, second row). Instead, the
separate layers of the papyrus document are locally reconstructed with many
connections between these two layers. As can be seen in Fig. 8, third row, our
algorithm successfully reconstructs an orientable 2D manifold. In order to
achieve this, we had to modify the image data by setting the scalar values in
regions where the folded document came too close to itself to 0. Note that the
same modified image data field was used for all approaches.

\begin{overpic}[width=390.25534pt]{figures/papyruscomparison2_small.jpg}\sf\small\put(0.0,32.0){{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\cite[cite]{[\@@bibref{}{Dambrogio:2021:Unlocking}{}{}]}}}
\put(0.0,19.0){{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\cite[cite]{[\@@bibref{}{Schultz:2009:CS}{}{}]}}}
\put(0.0,6.0){{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}ours}}
\end{overpic} Fig. 8: Comparison of reconstruction results for the folded
papyrus document shown in Fig. 7. First row: Method by Dambrogio et al. [9].
Note that the surface is not orientable at all, hence there is only one color.
Second row: Method by Schultz et al.[35]. Third row: Our method, resulting in
a fully orientable 2D manifold.

##  5 Discussion

In this section, we discuss the results of comparing our novel approach for 2D
manifold extraction from thin-layer structures with three other state-of-the-
art approaches [35, 1, 9].

Regarding methods based on the Hessian, we compared our method with the one by
Schultz et al. [35]. In agreement with previous work [1], we found that the
approach based on the Hessian is not robust to noise (see Paragraph 4.1.1 and
Sect. S2 in the Supplementary Material). Furthermore, we showed that it has
problems with strongly curved, thin-layer structures (see Paragraphs 4.2.1,
4.2.2). In all real test cases, their method, even though not being developed
for such data, did remarkably well. In particular, we were surprised how well
it worked on the silver sheet data. However, we could also clearly show its
limitations w.r.t. the characteristics of the data that motivated our work.
Even though improvements have been made in terms of speed and robustness [6]
compared to the method by Schultz et al., we still consider the latter method
as representative for the Hessian approach.

The optimal-path approach, on the other hand, as used by Algarni and
Sundaramoorthi [1] and in this work, can deal better with strong noise.
Unfortunately, the approach presented by Algarni and Sundaramoorthi [1] does
also not fulfill the requirements for the virtual unfolding of ancient
documents, motivating this work. We observed several problems for the type of
data we were dealing with. First, their approach makes use of one fast
marching solution for the entire data field. This has two implications. (1) As
the domain $\Omega$ is therefore the upper bound for the number of elements
for the fast marching algorithm, the performance is $\mathcal{O}(n\log n)$,
with $n$ being the number of voxels $\Omega$ is discretized in. This is
noteworthy as complexly folded thin-layer structures necessitate high-
resolution image data to allow differentiation between different layers and,
hence, $n$ can be very large. (2) As we illustrated in Fig. 2, a global front
is not able to follow thin-layer structures which are folded tightly back to
themselves, as it will take a shortcut to the next layer instead of following
the ridge. Yet, this is exactly what we often have to deal with in case of
historical, written documents. Second, their approach does not guarantee that
the manifold $M_{R}$ intersects $S_{\tilde{t}}$ as a closed 1D manifold.
Therefore, the identification of the ridge points on $S_{\tilde{t}}$ is
unstable and can result in incorrectly extracted ridge surfaces (see Sect. S1
in the Supplementary Material). We also found the ridge-curve extraction
method using a recursive watershed approach utilized by Algarni and
Sundaramoorthi [1] not stable enough for our data. This is particularly true
for papyrus, which often differs too much from a clear ridge structure. Third,
the cubical complex reduction of the volume to calculate the surface is not
guaranteed to be a 2D manifold. Thus, their approach to extract a surface from
a closed ridge line was not viable in our case (see Sect. S1 in the
Supplementary Material). To render the approach of optimal paths more suitable
for the type of data we have to deal with, we replaced all steps of the
algorithm of Algarni and Sundaramoorthi [1]. The resulting algorithmic
pipeline has the following properties: (1) it is faster; (2) it creates a
surface that is guaranteed to be an orientable 2D manifold; (3) it is viable
for highly curved surfaces; (4) it allows adjusting the local parameters to
the data, that is, to the denseness of the folded layers as well as the
strength of the folding; (5) for very challenging data, it allows the ridge
surface to be extracted partially or even completely user-driven, thus giving
the user full control over the extraction process.

Finally, we compared our method with the one proposed by Dambrogio et al. [9]
for unfolding letters. This very interesting approach seems to be rather
specific to the problem of unfolding letters. Limitations are that is assumes
a rather similar thickness of the document to be unfolded. This assumption
holds true for letters but not for papyrus documents and silver sheets.
Nevertheless, their algorithm also did remarkably well on both real datasets.
But none of the resulting reconstructed meshes was a 2D manifold and large
parts were missing in the initial reconstructed mesh (see Sect. S3 in the
Supplementary Material). In particular, the presence of an abundance of non-
manifold vertices, edges and faces leads to many problems when trying to post-
process the data with state-of-the-art mesh operations.

In terms of our own method, a few things are worth mentioning. As each seed
point only affects the resulting surface locally and our algorithm can be
implemented iteratively, addition, deletion and modification of seed points
can be implemented in such a way that these operations are very fast. In our
test cases, it took approximately 0.2 seconds to add one seed point and update
the surface. Of course, this time depends on the distance fast marching is
allowed to march, and hence, the distance should not be too large. However,
empirically, we observed that as long as the distance is a few times larger
than the thickest layer of the structure (10 times seems a good heuristic),
the parameter does not have a big influence on the quality of the generated
surface. Yet, the time performance gets worse for larger values of the
parameter and the patches themselves are bigger. For all our datasets, we used
the default distance. Finally, we would like to remark that our algorithm is
stable under the location of the starting seed point, which in our case is set
manually.

Regarding locality vs. globality, our method can be considered as being in-
between the methods based on the Hessian (e.g. [35]), which is a very local
one, and a purely global one as proposed by Algarni and Sundaramoorthi [1]. We
can adjust the degree of locality depending on the characteristics of the
thin-layer structures by changing the maximal distance that the fast marching
is allowed to march. The maximal distance parameter as well as a custom cost
function per seed point gives us a lot of flexibility that might be needed for
the challenging data we deal with. Even more flexibility is provided through
the combination of automation and visual interaction allowed by our approach.

##  6 Future Work and Conclusion

The biggest memory and performance bottleneck in our current method are the
label and time fields. For large image data, these fields are also large and
this results in many cache misses. A method to merge surface patches directly
could therefore render these fields obsolete. Furthermore, the consideration
made in Sect. 2.7 would not be necessary, allowing for greater marching
distances for strongly rolled structures. Last but not least, it would also
allow the creation of unorientable 2D manifolds, one limitation our method is
currently facing.

A more sophisticated watershed ridge curve extraction in Sect. 2.5 may allow
to discern between different equally important ridge contours. This
information could be used to decide if the front of the fast marching
algorithm took a shortcut to another layer or if intersecting or tangentially
connected surfaces inside the marched volume exist. A possible generalization
could also allow the construction of a surface and a graph with information of
their non-manifold connections, allowing the use of graph cut algorithms or
similar, to globally find the best possible 2D manifold or to allow
intersecting surfaces. Furthermore, artifacts on the border of the surface
(seen, e.g., in Fig. 3) are results of our method trying to maximize the
length of the optimal paths. Another ridge curve extraction approach could
mitigate this bias.

Automation can also be approved upon. We believe that a hybrid approach based
on both optimal paths and the Hessian could create even better results. One
could use as probability field $p(x)$ the probability of a point to be a ridge
point defined by the Hessian. This might have the advantages of both worlds,
the accuracy of the Hessian approach together with the guarantees of the
optimal path approach. A different option would be to place seed points
dependent on the Hessian.

Last but not least, anisotropic fast marching [27] might also increase
performance and quality of our algorithm, however some questions have to be
answered beforehand. How can one calculate the Euclidean distance efficiently
while calculating the time field, for example. A further problem is that the
anisotropic fast marching algorithm uses stencils reaching over many voxels
wide (dependent on the given tensor). This allows the front to “jump over”
some high cost voxels, for example, in densely folded areas. Hence,
anisotropic fast marching makes it difficult to work with a non-constant
scalar field.

Another remaining limitation of all currently available methods is the
extraction of surfaces for tangentially connected structures. The separation
of such structures will also be subject to future work.

In this paper, we have developed a local optimal path approach for extracting
smooth, simple surfaces in noisy 3D data sets. Our novel method is able to
extract such surfaces given highly complexly folded and rolled thin-layer
structures while taking up to only a single seed point as input. In contrast
to other methods, our approach guarantees the construction of an orientable 2D
manifold, which is a requirement of many virtually unfolding methods.
Furthermore, the locality of our method and, hence, the assignability of local
volumes to a seed point enables user interaction and modification of the
generated surface, resulting in a semi-automatic workflow that allows one to
address even very challenging data.

## Supplemental Materials

All synthetic data for the _Plane_ test case (Paragraph 4.1.1, Sect. 3 of the
Supplementary Material) used to do the quantitative evaluation of the methods,
including the Excel sheet containing the analysis results of Fig. 3, are
available on figshare at https://doi.org/10.6084/m9.figshare.23600115,
released under a CC BY 4.0 license.

###### Acknowledgements.

This work was supported by the German Research Foundation (DFG) through grants
BA 5042/2-1 and LE 1837/2-1 awarded to Daniel Baum and Verena Lepper,
respectively. We would like to thank Christiane Matz and Stefan Burmeister
from VARUSSCHLACHT im Osnabrücker Land gGmbH, Museum und Park Kalkriese, for
providing the folded silver sheet, and Eve Menei and Marc Etienne, Musée du
Louvre, Paris for providing the papyrus document. Image acquisitions of the
silver sheet and the papyrus document were funded by the VolkswagenStiftung
(Az 92029) and the European Research Council (Project ”ELEPHANTINE”, ID
637692), respectively. We would like to express our great thanks to Amanda
Ghassaei, one of the corresponding authors of the letter unfolding paper [9],
for her great support when running out data sets with their algorithm. Many
thanks also to our colleague Ben Schmitt for getting the letter unfolding
approach running on our computers and applying it to the various simple data,
as well as our colleague Finn Schwörer for implementing most scripts to export
our data sets to files compatible with the methods we tested our algorithm
against.

## References

  * [1] M. Algarni and G. Sundaramoorthi.  Surfcut: Surfaces of minimal paths from topological structures.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 41(3):726–739, 2019. doi: 10 . 1109/TPAMI . 2018 . 2811810
  * [2] F. Ambellan, S. Zachow, and C. von Tycowicz.  Rigid motion invariant statistical shape modeling based on discrete fundamental forms: Data from the osteoarthritis initiative and the Alzheimer’s disease neuroimaging initiative.  Medical Image Analysis, 73:article no. 102178, 18 pages, 2021. doi: 10 . 1016/j . media . 2021 . 102178
  * [3] N. Amenta, M. Bern, and M. Kamvysselis.  A new voronoi-based surface reconstruction algorithm.  In Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’98, pp. 415––421. Association for Computing Machinery, New York, NY, USA, 1998. doi: 10 . 1145/280814 . 280947
  * [4] V. I. Arnold.  Ordinary Differential Equations.  Springer Berlin, Heidelberg, 1 ed., 2006.
  * [5] S. Axler.  Measure, Integration & Real Analysis.  Springer Cham, 1 ed., 2019. doi: 10 . 1007/978-3-030-33143-6
  * [6] S. Barakat, N. Andrysco, and X. Tricoche.  Fast extraction of high-quality crease surfaces for visual analysis.  In Computer Graphics Forum, vol. 30, pp. 961–970. Wiley Online Library, 2011. doi: 10 . 1111/j . 1467-8659 . 2011 . 01945 . x
  * [7] S. Barakat and X. Tricoche.  An image-based approach to interactive crease extraction and rendering.  Procedia Computer Science, 1:1709–1718, 2010. doi: 10 . 1016/j . procs . 2010 . 04 . 192
  * [8] D. Baum, N. Lindow, H.-C. Hege, V. Lepper, T. Siopi, F. Kutz, K. Mahlow, and H.-E. Mahnke.  Revealing hidden text in rolled and folded papyri.  Applied Physics A, 123:1–7, 2017. doi: 10 . 1007/s00339-017-0808-6
  * [9] J. Dambrogio, A. Ghassaei, D. S. Smith, H. Jackson, M. L. Demaine, G. Davis, D. Mills, R. Ahrendt, N. Akkerman, D. Van der Linden, et al. Unlocking history through automated virtual unfolding of sealed documents imaged by x-ray microtomography.  Nature communications, 12(1):article no. 1184, 10 pages, 2021. doi: 10 . 1038/s41467-021-21326-w
  * [10] T. Deschamps and L. D. Cohen.  Fast extraction of minimal paths in 3d images and applications to virtual endoscopy.  Medical Image Analysis, 5(4):281–299, 2001. doi: 10 . 1016/S1361-8415(01)00046-9
  * [11] D. Eberly, R. Gardner, B. Morse, S. Pizer, and C. Scharlach.  Ridges for image analysis.  Journal of Mathematical Imaging and Vision, 4:353–373, 1994. doi: 10 . 1007/BF01262402
  * [12] M. Garland and P. S. Heckbert.  Surface simplification using quadric error metrics.  In Proceedings of the 24th annual conference on Computer graphics and interactive techniques, SIGGRAPH ’97, pp. 209–216. ACM Press/Addison-Wesley Publishing Co., 1997. doi: 10 . 1145/258734 . 258849
  * [13] D. Hale.  Fault surfaces and fault throws from 3D seismic images.  In 2012 SEG Annual Meeting, pp. 1–6. OnePetro, 2012. doi: 10 . 1190/segam2012-0734 . 1
  * [14] R. M. Haralick.  Ridges and valleys on digital images.  Computer vision, graphics, and image processing, 22(1):28–38, 1983. doi: 10 . 1016/0734-189X(83)90094-4
  * [15] H.-C. Hege, D. Stalling, M. Seebass, and M. Zöckler.  A generalized marching cubes algorithm based on non-binary classification.  Technical report, Konrad-Zuse Zentrum für Informationstechnik, Berlin, 1997.  SC-97-05.
  * [16] F. Herter, H.-C. Hege, M. Hadwiger, V. Lepper, and D. Baum.  Thin-volume visualization on curved domains.  In Computer Graphics Forum, vol. 40, pp. 147–157. Wiley Online Library, 2021. doi: 10 . 1111/cgf . 14296
  * [17] G. Hoffmann Barfod, J. M. Larsen, A. Lichtenberger, and R. Raja.  Revealing text in a complexly rolled silver scroll from Jerash with computed tomography and advanced imaging software.  Scientific reports, 5(1):1–10, 2015. doi: 10 . 1038/srep17765
  * [18] M. Kazhdan, M. Bolitho, and H. Hoppe.  Poisson surface reconstruction.  In Proceedings of the fourth Eurographics symposium on Geometry processing, vol. 7 of SGP ’06, pp. 61–70. Eurographics Association, 2006. doi: 10 . 2312/SGP/SGP06/061-070
  * [19] M. Kazhdan and H. Hoppe.  Screened poisson surface reconstruction.  ACM Trans. Graph., 32(3), jul 2013. doi: 10 . 1145/2487228 . 2487237
  * [20] G. Kindlmann, C. Chiw, T. Huynh, A. Gyulassy, J. Reppy, and P.-T. Bremer.  Rendering and extracting extremal features in 3D fields.  In Computer Graphics Forum, vol. 37, pp. 525–536. Wiley Online Library, 2018. doi: 10 . 1111/cgf . 13439
  * [21] G. Kindlmann, X. Tricoche, and C.-F. Westin.  Anisotropy creases delineate white matter structure in diffusion tensor mri.  Lecture Notes in Computer Science, 4190:126–133, 2006. doi: 10 . 1007/11866565_16
  * [22] G. L. Kindlmann, R. S. J. Estépar, S. M. Smith, and C.-F. Westin.  Sampling and visualizing creases with scale-space particles.  IEEE Transactions on Visualization and Computer Graphics, 15(6):1415–1424, 2009. doi: 10 . 1109/TVCG . 2009 . 177
  * [23] J. J. Koenderink and A. J. van Doorn.  Local features of smooth shapes: ridges and courses.  In Geometric methods in computer vision II, vol. 2031, pp. 2–13. SPIE, 1993. doi: 10 . 1117/12 . 146617
  * [24] J. M. Lee.  Integral Curves and Flows, pp. 205–248.  Springer New York, New York, NY, 2 ed., 2013. doi: 10 . 1007/978-1-4419-9982-5
  * [25] W. E. Lorensen and H. E. Cline.  Marching cubes: A high resolution 3D surface construction algorithm.  SIGGRAPH Computer Graphics, 21(4):163–169, Aug. 1987. doi: 10 . 1145/37402 . 37422
  * [26] H.-E. Mahnke, T. Arlt, D. Baum, H.-C. Hege, F. Herter, N. Lindow, I. Manke, T. Siopi, E. Menei, M. Etienne, et al. Virtual unfolding of folded papyri.  Journal of cultural heritage, 41:264–269, 2020. doi: 10 . 1016/j . culher . 2019 . 07 . 007
  * [27] J.-M. Mirebeau.  Anisotropic fast-marching on cartesian grids using lattice basis reduction.  SIAM Journal on Numerical Analysis, 52(4):1573–1599, jan 2014. doi: 10 . 1137/120861667
  * [28] V. Mocella.  personal communication.
  * [29] V. Mocella, E. Brun, C. Ferrero, and D. Delattre.  Revealing letters in rolled Herculaneum papyri by x-ray phase-contrast imaging.  Nature communications, 6(1):article no. 5895, 6 pages, 2015. doi: 10 . 1038/ncomms6895
  * [30] D. Neuber and C. Reinhart.  Geheimnisvolle bleischriftrolle der mandäer: Moderne computertomographie erlaubt virtuelle entschlüsselung eines antiken rätsels.  ZfP-Zeitung, 132:39–41, 2012.
  * [31] R. Peikert and F. Sadlo.  Height ridge computation and filtering for visualization.  In 2008 IEEE Pacific Visualization Symposium, pp. 119–126. IEEE, 2008. doi: 10 . 1109/PACIFICVIS . 2008 . 4475467
  * [32] F. Sadlo and R. Peikert.  Efficient visualization of lagrangian coherent structures by filtered AMR ridge extraction.  IEEE Transactions on Visualization and Computer Graphics, 13(6):1456–1463, 2007. doi: 10 . 1109/TVCG . 2007 . 70554
  * [33] J. Sahner, T. Weinkauf, N. Teuber, and H.-C. Hege.  Vortex and strain skeletons in Eulerian and Lagrangian frames.  IEEE Transactions on Visualization and Computer Graphics, 13(5):980–990, 2007. doi: 10 . 1109/TVCG . 2007 . 1053
  * [34] O. Samko, Y.-K. Lai, D. Marshall, and P. L. Rosin.  Virtual unrolling and information recovery from scanned scrolled historical documents.  Pattern Recognition, 47(1):248–259, 2014. doi: 10 . 1016/j . patcog . 2013 . 06 . 015
  * [35] T. Schultz, H. Theisel, and H.-P. Seidel.  Crease surfaces: From theory to extraction and application to diffusion tensor mri.  IEEE Transactions on Visualization and Computer Graphics, 16(1):109–119, 2009. doi: 10 . 1109/TVCG . 2009 . 44
  * [36] W. B. Seales, C. S. Parker, M. Segal, E. Tov, P. Shor, and Y. Porath.  From damage to discovery via virtual unwrapping: Reading the scroll from En-Gedi.  Science advances, 2(9):article no. e1601247, 9 pages, 2016. doi: 10 . 1126/sciadv . 160124
  * [37] J. A. Sethian.  A fast marching level set method for monotonically advancing fronts. Proceedings of the National Academy of Sciences, 93(4):1591–1595, 1996. doi: 10 . 1073/pnas . 93 . 4 . 1591
  * [38] J. A. Sethian.  Fast marching methods.  SIAM Review, 41(2):199–235, 1999. doi: 10 . 1137/S0036144598347059
  * [39] J. A. Sethian and A. M. Popovici.  3-d traveltime computation using the fast marching method.  GEOPHYSICS, 64(2):516–523, 1999. doi: 10 . 1190/1 . 1444558
  * [40] Stanford University Computer Graphics Laboratory.  The Stanford 3D Scanning Repository.  https://graphics.stanford.edu/data/3Dscanrep, 2014.  [Online; accessed 10-March-2023].
  * [41] D. Vavřík, K. Knauber, D. Urbanová, I. Kumpová, K. Blažková, and Z. Šámal.  Unveiling magic from the middle ages: tomographic reading of a folded lead amulet from dřevíč fortress (czech republic).  Archaeological and Anthropological Sciences, 12:1–12, 2020. doi: 10 . 1007/s12520-019-00976-4
  * [42] B. Wilster-Hansen, D. C. Mannes, K. L. Holmqvist, K. Ødeby, and H. Kutzke.  Virtual unwrapping of the BISPEGATA amulet, a multiple folded medieval lead amulet, by using neutron tomography.  Archaeometry, 64(4):969–978, 2022. doi: 10 . 1111/arcm . 12734