# Solid NURBS Conforming Scaffolding for Isogeometric Analysis

Stefano Moriconi  King’s College LondonStrandLondonUnited Kingdom
<EMAIL_ADDRESS>,  Parashkev Nachev  University College LondonGower
St.LondonUnited Kingdom<EMAIL_ADDRESS>,  Sébastien Ourselin  King’s
College LondonStrandLondonUnited Kingdom<EMAIL_ADDRESS>and  M.
Jorge Cardoso  King’s College LondonStrandLondonUnited Kingdom
<EMAIL_ADDRESS>
###### Abstract.

This work introduces a scaffolding framework to compactly parametrise solid
structures with conforming NURBS elements for isogeometric analysis. A novel
formulation introduces a topological, geometrical and parametric subdivision
of the space in a minimal plurality of conforming vectorial elements. These
determine a multi-compartmental scaffolding for arbitrary branching patterns.
A solid smoothing paradigm is devised for the conforming scaffolding achieving
higher than positional geometrical and parametric continuity. Results are
shown for synthetic shapes of varying complexity, for modular CAD geometries,
for branching structures from tessellated meshes and for organic biological
structures from imaging data. Representative simulations demonstrate the
validity of the introduced scaffolding framework with scalable performance and
groundbreaking applications for isogeometric analysis.

NURBS, Solid, Scaffolding, Conforming Lattice, Branching, Organic,
Isogeometric Analysis

††booktitle: Woodstock ’18: ACM Symposium on Neural Gaze Detection, June
03–05, 2018, Woodstock, NY††price: 15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs:
Theory of computation Computational geometry††ccs: Mathematics of computing
Geometric topology††ccs: Computing methodologies Parametric curve and surface
models††ccs: Computing methodologies Volumetric models††ccs: Computing
methodologies Physical simulation††ccs: Applied computing Computer-aided
design††ccs: Applied computing Imaging Graphical Abstract

##  1\. Introduction

Devising an end-to-end framework to trace parametric geometries fully
compatible with continuous-domain computational simulations is an extremely
complex task, and currently stands as an open challenge. The design of
geometrical structures traditionally employs computer-graphics techniques
based on modelling raster surfaces, which often require high polygonal density
to accurately represent smooth and irregular shapes (Mortenson, 1997; Russo,
2006). In computer-aided design (CAD) vectorial and piece-wise continuous
entities, i.e. Non-Uniform Rational B-Spline (NURBS) patches, leverage similar
polygonal subdivision schemes in a more compact and resolution-free
parametrisation by representing the underlying geometry exactly (Piegl, 1991;
Piegl and Tiller, 1996). Nonetheless, most of the commercially available CAD
tools exclusively focus on surface geometries and do not integrate parametric
simulations or a solid finite-element analysis framework in the same vectorial
domain. Solid and volumetric meshing is predominantly addressed in a raster
fashion, where discrete polyhedral subdivision schemes fill with an
unstructured plurality of solid elements the hollow space delimited by the
superficial (closed) boundary (Ho-Le, 1988). In this case, conventional
finite-elements methods ultimately approximate partial differential equations
(PDEs) solutions over such discrete solid tessellations (Rao, 2017).

Isogeometric analysis (IGA) (Hughes et al., 2005; Cottrell et al., 2009) is a
relatively recent computational approach that directly integrates finite-
element analysis with vectorial NURBS-based geometries. By means of highly
compact shape parametrisations, IGA has the advantage of high-performance and
high-throughput continuous simulations while reducing geometrical
approximation errors and providing higher numerical stability for the
simulation solution profiles. Raster polygonal or polyhedral meshes, however,
are not compatible with IGA. Conversely, a consistent NURBS-based geometrical
domain is required to satisfy stringent conforming parametric conditions. The
lack of a parametric framework to compactly trace solid shapes in IGA-
compatible domains currently limited the development and application of
Isogeometric Analysis to simple 2D/3D toy examples and few case studies.

With the aim of bridging the technological gap, this work first proposes a
novel solid scaffolding framework jointly leveraging a topological,
geometrical and parametric formulation. A solid scaffolding of minimal
parameters and vectorial elements is sought to model the shape of the
underlying structure by enforcing a conforming condition and preserving
smoothness properties of a solid continuum, while optimising for the
computational complexity of downstream simulations. The introduced scaffolding
framework paves the way towards novel approaches to rationalising space for
industrial design, CAD modelling, organic structures representation, and
compact digital shape tracing with a particular focus on scalable, accurate
and high-performance finite-elements analysis using IGA.

###  1.1. Related Work

##### NURBS Geometries

The pioneering work of (Piegl, 1991; Piegl and Tiller, 1996) first introduced
a parametrisation framework to exactly represent 3D geometries with vectorial
elements as opposed to raster finite-resolution meshes employed in polygonal
modelling. Curvilinear and surface NURBS elements have de-facto become the
standard for computed-aided design (CAD), manufacturing and engineering in a
number of commercial and industrial applications. Beyond the established
standards, the NURBS formulation generalises also for the parametrisation of
higher-dimensional entities (Chen et al., 2016; Xu et al., 2014). Among those,
solid NURBS build on extending and combining linear NURBS operators on
univariate primitives, thus constituting a free-form cuboid in the most simple
form. Recent research toolboxes (De Falco et al., 2011; Vázquez, 2016; Bingol
and Krishnamurthy, 2019; Dalcin et al., 2016; Jüttler et al., 2014) have
developed a multi-platform formulation for high-dimensional vectorial
representations, including hierarchical NURBS (Garau and Vázquez, 2018; Bracco
et al., 2018) and IGA-compatible frameworks (Hesch et al., 2016; Bazilevs et
al., 2010) for the solution of PDEs based on splines. However, the design and
construction of higher-dimensional geometrical domains has only been addressed
for simple toy examples or individual case studies (Hughes et al., 2005;
Cottrell et al., 2009).

##### Scaffolding Construction

The design and construction of geometrical partitioning domains was initially
addressed in (Suárez and Hubert, 2017, 2018; Panotopoulou et al., 2018) with a
scaffolding approach. In particular, a quadrilateral mesh is first introduced
as the coarsest shape representation around a one-dimensional skeleton of high
genus. Proof is provided for the existence of such coarse representation based
on Voronoi diagrams and the problem of constructing an optimal scaffold is
formalised as an integer linear program (Suárez and Hubert, 2018). In these
cases, affinities between the scaffolding and tensor product splines are
mentioned, yet a structural subdivision and a parametric formulation specific
for NURBS elements and their arrangement are not addressed. Also, strategies
to recover the geometrical embedding of the underlying solid structure from
the scaffolding are limited to coarse polygonal refinement schemes.

##### Shape Estimation

In (Livesu et al., 2016) a solid meshing approach builds from a skeleton
topology and integrates an adaptive refinement strategy to locally recover the
shape of a structure delimited by a surface mesh. In a similar fashion,
(Livesu et al., 2020) proposed a hexahedral meshing method which optimises for
cuboid elements leveraging a cascade of loop-wise cuts based on a surface-
field aware block decomposition of a generic mesh. Alternative solid meshing
techniques (Tarini et al., 2004; Livesu et al., 2013; Gao et al., 2017;
Gregson et al., 2011; Si, 2015; Wang et al., 2008) approximate the geometry of
the tessellation by locally sampling the boundary of the structure and further
subdivide the solid mesh in smaller polyhedra or finer inner partitions with
adaptive or pre-defined resolution. Shape and curvature features as well as
geometrical cues are often retrieved from the boundary mesh with directional
fields (Panozzo et al., 2014; Pietroni et al., 2016) and geodesic (Kimmel and
Sethian, 1998; Martinez et al., 2004; Surazhsky et al., 2005; Crane et al.,
2020) techniques, by first partitioning and segmenting superficial patches of
the geometry and ultimately fitting the 3D coordinates of the solid polyhedral
mesh. Dense solid meshes recover accurately the underlying geometry of the
structure, often in detriment to a compact representation, where the resulting
unorganised plurality of polyhedra exhibits a piece-wise linear and faceted
lattice with pure positional continuity among neighbouring elements.

##### Smoothness and Continuity

A number of studies (Catmull and Clark, 1978; Doo and Sabin, 1978; Loop, 1987;
Stam, 1998a) have proposed methods and strategies to recover a smoother
lattice and higher degrees of continuity for associated limit surfaces
focusing on polygonal meshes. In general, they build on refining first the
mesh with a denser lattice and subsequently adjusting the coordinates of the
geometry based on the local neighbourhood. In (Stam, 2001) a new class of
subdivision surfaces is introduced by bridging the gap between polygonal
meshes and uniform B-spline surfaces of arbitrary degree. In (Stam, 1998b), an
exact solution is provided for similar spline-based surfaces, and a closed
form is devised for patches joining at extraordinary vertices of arbitrary
valence. Irregular surface structures exactly recover the typically organic
smoothness and higher degree of parametric continuity in their differential
forms, which provided accurate simulation profiles for surface geometries
using IGA (Pan et al., 2016). Rectification routines and extended subdivisions
are devised for convex polyhedral meshes (Livesu et al., 2015; Burkhart et
al., 2010), where the coordinates of the inner lattice are constrained by the
local connected neighbourhood. Yet, a solid smoothing strategy for higher
degrees of continuity in a conforming solid NURBS scaffolding for IGA has not
been investigated, and currently remains an open challenge.

##### IGA Applications

An increasing number of studies (Bazilevs et al., 2006; Bazilevs et al., 2008;
Cottrell et al., 2009; Urick et al., 2019; Carraturo et al., 2019; Bucelli et
al., 2021) are investigating and employing finite-element analysis striving
for a unified framework embedding vectorial geometries with continuous-domain
accurate and fast simulations. Among those, the pioneering work by (Hughes et
al., 2005; Bazilevs et al., 2006) first introduced a CAD-compatible
integration and a monolithic formulation of computational analyses with IGA,
where spline-basis functions solve for linear elastic problems and
(bio)-mechanical characterisations, as opposed to conventional polynomial
finite-element methods (Rao, 2017). In early medical applications, (Zhang et
al., 2007; Bazilevs et al., 2008) presented a subject-specific approach for
tracing vessels with a set of conforming solid NURBS geometries. In a similar
fashion, (De Falco et al., 2011; Vázquez, 2016; Bingol and Krishnamurthy,
2019; Jüttler et al., 2014) developed vectorial tools for modelling complex
interactions between composite and multi-compartmental structures in a
compact, general-purpose and light-weight manner. Quantitative IGA analyses
proved higher accuracy and numerical stability (Bazilevs et al., 2006)
compared with conventional approaches, as the solution profile of PDEs
eliminates the propagation of geometrical approximation errors. Most works,
however, have focused on a limited set of geometries, addressing feasibility
studies and proof-of-concept applications. This is mainly due to the lack of a
principled scaffolding framework able to capture the underlying structure of
more complex geometries and generalise for compatible real case scenarios on
large scale (Bazilevs et al., 2010; Urick et al., 2019).

###  1.2. Outline

Aiming to tackle the aforementioned open problems in a unified manner, a
conforming solid scaffolding for a generic 3D structure is described in the
following sections. In section 2, a graph configuration is first considered as
the underlying skeleton of the solid structure. In section 3, the generic form
of a solid scaffolding is defined as a structured set of adjacent vectorial
elements, and the conforming condition is introduced. In section 3.2 the
luminal solid scaffolding construction is introduced, where the minimal number
of vectorial elements, the associated spatial arrangement and the conforming
parametric configurations are detailed. A straightforward extension, i.e. the
wall scaffolding, is addressed in appendix A and alternative formulations are
mentioned in appendix B. An interfacing configuration is introduced in section
3.2 for vectorial elements joining at a generic junction with arbitrary
branching pattern. In section 4 the geometrical embedding of the scaffolding
elements is described, including a data-driven fitting scheme. Lastly, in
section 5, a smoothing paradigm is devised for the scaffolding, achieving
higher degrees of geometrical and parametric continuity, as in an organic
medium in the form of a continuum.

##  2\. Graph Configuration of the Structure

A solid structure is often characterised by one or more branches, as elongated
regions protruding from its main body. Complex structures may exhibit holes or
concave regions, others have the form of laminar thin or thick shells.
Convoluted and irregular organic structures, as well as composite CAD
geometries, may show several of the above components, which can be locally
decomposed in sub-blocks and modules (Chuang et al., 2000; Lu et al., 2017;
Wang et al., 2017; Takayama, 2019). In any case, a graph can be determined as
the skeleton underlying a branching portion of the structure by means of a
number of (semi-)automatic computational routines and methods for different
data types (Saha et al., 2016, 2017). Such graph embeds both spatial and
topological configurations, where generally the nodes correspond to junctions
or terminal end-points, and the edges correspond to the elongated connecting
regions underlying each branch fig. 1.

\begin{overpic}[width=429.28616pt]{fig1R.png} \put(0.0,19.0){\small{(a)}}
\put(31.0,19.0){\small{(b)}} \put(63.0,19.0){\small{(c)}} \end{overpic} Figure
1. Graph configuration of the structure: (a) Skeleton underlying a 3D surface
mesh: Armadillo; (b) Connecting graph of a biological structure from imaging
data: Phantom; (c) Skeleton underlying the modular (red) CAD geometry: Gear by
construction.

##  3\. Conforming Solid Scaffolding

The scaffolding builds on the graph configuration of the branching structure
in section 2, as it subdivides each connecting branch in a pre-defined minimal
number of vectorial elements. A rational partitioning of the structure’s space
is introduced, and the scaffolding configures a set of conforming elements
being jointly defined over high-dimensional, geometrical and parametric
domains. In particular, the scaffolding comprises an organised and structured
lattice of vectorial elements, each conforming to the neighbouring ones, and
constitutes a conforming solid NURBS multi-patch fig. 2. The conforming
condition is necessary for an IGA-compatible domain, and it represents the key
parametric constraint and the structural criterion at the base of the solid
scaffolding.

\begin{overpic}[width=429.28616pt]{fig2aR.png} \put(1.0,20.0){\small{(a)}}
\put(35.0,20.0){\small{(b)}} \put(68.0,20.0){\small{(c)}} \end{overpic}
---
\begin{overpic}[width=429.28616pt]{fig2bR.png} \put(3.0,30.0){\small{(d)}}
\put(43.0,30.0){\small{(e)}} \put(68.0,30.0){\small{(f)}}
\put(39.0,10.0){\small{$u$}} \put(38.0,3.5){\small{$v$}}
\put(41.0,7.0){\small{$w$}} \end{overpic}
Figure 2. Solid Scaffolding of a Torus: (a) Minimal element, (b) Equivalent
multi-patch (4 elements) with reflection lines; (c) Multi-patch explosion
view; (d) Individual element: control points lattice and inner solid
sectioning; (e) Parametric embedding: individual element isomorphic to a cube;
(f) Multi-patch sides (edges) facing (along) directions $(u,v,w)$,
respectively: dark/light red for $u=\\{0,1\\}$, dark/light green for
$v=\\{0,1\\}$, and dark/light blue for $w=\\{0,1\\}$ \- conforming interfaces
facing $w$; boundaries facing $u$ and $v$.

###  3.1. Vectorial Element

A vectorial element is the atomic component of the scaffolding and consists of
a 3D free-form solid patch determined by NURBS primitives fig. 2. As
introduced in (Piegl and Tiller, 1996), the vectorial element
$\textrm{H}(u,v,w)$ is defined as a trivariate tensor product of NURBS curves
along three independent parametric directions $(u,v,w)$, and has the form

(1)
$\begin{matrix}\textrm{H}(u,v,w)=\sum_{i=1}^{l_{u}}\sum_{j=1}^{l_{v}}\sum_{k=1}^{l_{w}}R_{i,j,k}(u,v,w)\mathbf{P}_{i,j,k},&\textrm{with}\end{matrix}$
(2)
$R_{i,j,k}(u,v,w)=\frac{N_{i,d_{u}}(u)N_{j,d_{v}}(v)N_{k,d_{w}}(w)\textrm{w}_{i,j,k}}{\sum_{\hat{i}}^{l_{u}}\sum_{\hat{j}}^{l_{v}}\sum_{\hat{k}}^{l_{w}}N_{\hat{i},d_{u}}(u)N_{\hat{j},d_{v}}(v)N_{\hat{k},d_{w}}(w)\textrm{w}_{\hat{i},\hat{j},\hat{k}}},$

where $R_{i,j,k}(u,v,w)$ are the rational basis functions;
$\mathbf{P}_{i,j,k}\in\mathbb{R}^{3}$ are the control points of cardinality
$l_{u}$, $l_{v}$, and $l_{w}$ in each independent parametric direction;
$N_{i,d_{u}}(u)$, $N_{j,d_{v}}(v)$, and $N_{k,d_{w}}(w)$ are the univariate
B-spline basis functions of degree $d_{u}$, $d_{v}$, and $d_{w}$ respectively
defined on the three independent parametric directions; and
$\textrm{w}_{i,j,k}$ are the strictly positive weights of the rational basis
functions. The span of the basis functions’ parametric domain is sampled in
each independent direction by a set of intervals defined by the knot vectors
$\mathbf{k}_{u}$, $\mathbf{k}_{v}$ and $\mathbf{k}_{w}$ respectively. The
generic knot vector $\mathbf{k}$ comprises a set of real-valued knots
${k_{1},...,k_{r}}\in[0,1]$ and is defined as

(3)
$\begin{matrix}\mathbf{k}=\\{~{}\underbrace{k_{1},\dots,k_{1}}_{\textrm{m}}~{},~{}k_{2},\dots,k_{r-1}~{},~{}\underbrace{k_{r},...,k_{r}}_{\textrm{m}}~{}\\},&\textrm{with}&r=l+d+1.\end{matrix}$

The knot vector is assumed open, i.e. $k_{1}=0$, $k_{r}=1$ and the
multiplicity $\textrm{m}=d+1$, with $r$ being the total cardinality of the
knots, $l$ being the number of control points in the generic parametric
direction, and $d$ being the arbitrary degree of the associated univariate
B-spline basis function.

The parametrisation in eq. 1 maps a solid tensor grid defined on the
parametric domain to a physical-space, where each vectorial element represents
a solid free-form cuboid. The free-form shape is modulated by the 3D
coordinates of the control points $\mathbf{P}_{i,j,k}$ and by the weights
$\textrm{w}_{i,j,k}$ of the rational basis functions. The boundary of the
element comprises 6 sides of quadrilateral profile, which are oriented, in
opposite pairs, to face the directions of the parametric domain.

####  3.1.1. Conforming Condition

A pair of vectorial elements is said to be adjacent when there is at least one
side in common for each element. Such shared (or adjacent) sides are referred
as interfaces, whereas the free-end unshared sides are referred as boundaries.
The same pair of adjacent elements is also said to be conforming when all the
interfaces are conforming, i.e. the parametrisation of each pair of adjacent
sides is matching. In other words, for each interface, all the control points,
as well as, all the knot vectors and all the weights of the rational basis
functions and the degrees of the associated univariate B-spline basis
functions defining the sides of the adjacent sides must coincide, up to a
different orientation of the elements in the parametric domain fig. 2.

A conforming solid scaffolding requires all adjacent elements in the lattice
jointly meet the conforming condition.

###  3.2. Luminal Scaffolding Construction

A luminal scaffolding builds on the branching structure graph in section 2,
where higher density of branches and greater complexity of the branching
pattern require an increasing number of elements to minimally parametrise the
underlying spatial regions. The luminal scaffolding comprises the organised
set of vectorial elements filling the innermost space of a branching
structure, e.g. the fluid region as in a set of tubes or connected chambers,
or the volumetric space occupied by a solid object. Leveraging the free-form
shape of each element, the luminal scaffolding models the branching structure
as a set of connecting portions that elongate for each structure protrusion
encoded in the graph (fig. 3).

\begin{overpic}[width=429.28616pt]{fig3R.png} \put(1.0,18.0){\small{(a)}}
\put(22.0,18.0){\small{(b)}} \put(44.0,18.0){\small{(c)}}
\put(63.0,18.0){\small{(d)}} \end{overpic} Figure 3. Luminal Scaffolding
Construction from Graph: (a) Bifurcation; (b) 5-way-junction: 4 elements
required per branch; (c) Cross-sectional quadrants arrangement and elements
adjacency per branch - note the conforming interfaces of adjacent elements
couple pairs of sides facing $u$ and $v$; (d) Scaffolding construction for a
3D mesh structure (Fertility) with multiple junctions and elongated portions.
For simplicity, all scaffolding elements are shown as piece-wise linear
segments.

##### Minimal Elements

The scaffolding subdivides the branching structure in a pre-defined number of
vectorial elements based on the topology of the graph. In general, the luminal
scaffolding requires 4 vectorial elements per branch, assuming branches
arbitrarily connect with generic $n$-way-junctions (fig. 3).

##### Arrangement and Adjacency

For a generic graph with arbitrary branching pattern, the luminal scaffolding
arranges each branch in 4 organised cross-sectional quadrants, each accounting
for an individual element. This results in 4 adjacent vectorial elements
elongated along the longitudinal direction of the branch, either joining a
junction with a terminal endpoint, alternatively connecting two junctions.
Assuming the parametric direction $w$ mapping the longitudinal direction along
the branch, each element is adjacent to other two elements of the same branch
along one of its cross-sectional sides. In particular, two interfaces are
determined for each vectorial element in the correspondence of the interior
part of the quadrant: one interface facing $u$ and the other interface facing
$v$ respectively (dark red for $u$ and light green for $v$ in fig. 3). The
remaining pair of sides, i.e. the boundaries, determines the exterior part of
the quadrant: respectively, one boundary side facing $u$ (light red side) and
another boundary side facing $v$ (dark green side), opposite to the
interfaces. The boundaries define a pseudo-circumferential profile relative to
the elongation of the branch.

At any terminal branch, each element of the scaffolding terminates with a
boundary side facing the longitudinal direction $w$. At junctions, each
element determines an interface with another adjacent element of a
neighbouring incident branch. The scaffolding configuration at the junction
interfaces is defined by a quadrilateral junction simplex.

##### Quadrilateral Junction Simplex

\begin{overpic}[width=429.28616pt]{fig4R.png} \put(2.0,21.0){\small{(a)}}
\put(27.0,21.0){\small{(b)}} \put(52.0,21.0){\small{(c)}}
\put(77.0,21.0){\small{(d)}}
\put(4.0,15.0){\small{$\hat{\mathbf{n}}^{(5)}_{1}$}}
\put(10.0,17.0){\small{$\hat{\mathbf{n}}^{(5)}_{2}$}}
\put(16.0,15.0){\small{$\hat{\mathbf{n}}^{(5)}_{3}$}}
\put(12.0,2.0){\small{$\hat{\mathbf{n}}^{(5)}_{4}$}}
\put(8.0,8.0){\small{$\hat{\mathbf{n}}^{(5)}_{5}$}}
\put(28.0,11.0){\small{$\hat{\mathbf{n}}_{1}$}}
\put(43.0,18.0){\small{$\hat{\mathbf{n}}_{2}$}}
\put(37.0,2.0){\small{$\hat{\mathbf{n}}_{3}$}}
\put(45.0,4.0){\small{$\hat{\mathbf{v}}$}}
\put(34.0,20.0){\small{$\hat{\mathbf{n}}_{1,2}^{\perp}$}}
\put(47.0,11.0){\small{$\hat{\mathbf{n}}_{2,3}^{\perp}$}}
\put(30.0,3.0){\small{$\hat{\mathbf{n}}_{3,1}^{\perp}$}}
\put(70.0,4.0){\small{$\hat{\mathbf{v}}$}} \end{overpic} Figure 4.
Quadrilateral junction simplex: (a) Set of incident directions
$\hat{\mathbf{n}}^{(n)}$ of a generic $n$-way-junction, here $n$ = 5. (b)
Associated atomic simplex $\mathcal{Q}^{\textrm{atom}}$ from three clustered
principal directions: white nodes (original extraordinary vertices in black),
colour-coded quadrant arcs surrounding the initial directional partition with
a quadrilateral profile. (c) Resulting simplex $\mathcal{Q}$ after $n$-3
iterative bisections: new extraordinary vertices in white. (d) Quadrilateral
and quadrant-based interfaces all facing $w$ constructed on $\mathcal{Q}$.

In a generic junction, all the elements of each incident branch adjacently
meet in the neighbourhood of the junction node. This determines a structured
set of interfaces, all facing the $w$ parametric direction, where each
interface (as any other side) of the elements are configured with a
quadrilateral profile.

These structural requirements are met by introducing the quadrilateral
junction simplex $\mathcal{Q}$ (fig. 4). The simplex is a rational
partitioning of the junction’s space and determines an interfacing scaffolding
configuration, whose atomic form underlies a bifurcation.

The topology of the atomic simplex $\mathcal{Q}^{\textrm{atom}}$ accounts for
a total of 5 nodes and 6 edges being configured in 3 loops, each comprising 4
nodes and 4 edges, which are jointly connected by 2 common nodes, i.e. the
initial extraordinary vertices of the atomic simplex. The associated spatial
embedding of the atomic simplex delimits the interfacing cross-sectional areas
of three incident branches, where each loop of $\mathcal{Q}^{\textrm{atom}}$
surrounds the incidental direction of the branch with 4 connecting arcs. In
this form, each loop defines a quadrilateral profile, and each arc of each
loop underlies a quadrant.

The geometrical embedding of an atomic simplex is formulated by considering
the directions of the incident branches as three non-coincident unit vectors
$\hat{\mathbf{n}}_{i}\in\mathbb{R}^{3}$, with $i=\\{1,2,3\\}$, centred at the
origin of a local reference system. The initial extraordinary vertices lie
along a common axis in opposite directions, where the common axis
$\hat{\mathbf{v}}$ is defined as

(4)
$\hat{\mathbf{v}}=\frac{\sum_{i=1}^{3}\hat{\mathbf{n}}_{i}\times\hat{\mathbf{n}}_{j}}{|\sum_{i=1}^{3}\hat{\mathbf{n}}_{i}\times\hat{\mathbf{n}}_{j}|}.$

The unit vector $\hat{\mathbf{v}}$ is determined as the average of the cross-
products of the incidental directions considered in pairs, in circular shift,
i.e. $j=i~{}(\textrm{mod}~{}3)+1$. The remaining three independent nodes of
the atomic simplex $\mathcal{Q}^{\textrm{atom}}$ lie along the intermediate
directions $\hat{\mathbf{n}}_{i,j}^{\perp}$, which are orthogonal to the
common axis $\hat{\mathbf{v}}$. The intermediate directions are defined as

(5)
$\begin{matrix}\hat{\mathbf{n}}_{i,j}^{\perp}=f_{i,j}~{}\frac{\hat{\mathbf{n}}_{i,j}-(\hat{\mathbf{n}}_{i,j}\cdot\hat{\mathbf{v}})\hat{\mathbf{v}}}{|\hat{\mathbf{n}}_{i,j}-(\hat{\mathbf{n}}_{i,j}\cdot\hat{\mathbf{v}})\hat{\mathbf{v}}|},&\textrm{with}&\hat{\mathbf{n}}_{i,j}=\frac{\hat{\mathbf{n}}_{i}+\hat{\mathbf{n}}_{j}}{|\hat{\mathbf{n}}_{i}+\hat{\mathbf{n}}_{j}|},\end{matrix}$

being $\hat{\mathbf{n}}_{i,j}$ the bisecting unit vector between each pair of
incident directions. The flip factor $f_{i,j}$ determines whether the
associated intermediate direction are considered in opposite sign. This avoids
inconsistent partitions of the space for particularly close and narrow sets of
incidental directions. Each flip factor $f_{i,j}$ is defined as

(6)
$\begin{matrix}f_{i,j}=\bigg{\\{}\begin{matrix}-1~{}~{}&\textrm{when}~{}~{}~{}p_{i,j}>0~{}\land~{}p_{j,k}\leq
0~{}\land~{}p_{k,i}\leq 0\\\
1~{}~{}&\textrm{otherwise}\end{matrix},&\textrm{where}\end{matrix}$ (7)
$\begin{matrix}p_{i,j}=\hat{\mathbf{n}}_{i,j}\cdot\hat{\mathbf{n}}_{k}^{\perp},&\textrm{and}&\hat{\mathbf{n}}_{i}^{\perp}=\frac{\hat{\mathbf{n}}_{i}-(\hat{\mathbf{n}}_{i}\cdot\hat{\mathbf{v}})\hat{\mathbf{v}}}{|\hat{\mathbf{n}}_{i}-(\hat{\mathbf{n}}_{i}\cdot\hat{\mathbf{v}})\hat{\mathbf{v}}|}.\end{matrix}$

The projection $p_{i,j}$ is computed as scalar-product between the bisecting
unit vector $\hat{\mathbf{n}}_{i,j}$ and $\hat{\mathbf{n}}_{k}^{\perp}$, being
the latter the unit vector associated to the $k$-th incident direction and
orthogonal to the common axis $\hat{\mathbf{v}}$, with
$k=(i+1~{}(\textrm{mod}~{}3))+1$.

In the general form, the quadrilateral junction simplex $\mathcal{Q}$ of a
junction with $n$ incident branches topologically consists in an adjacent
stack of atomic simplexes, where each pair mutually shares a common loop.

In a physical-space embedding, such topological composition is equivalent to a
cascade of spatial bisections of an initial atomic simplex, which first
considers three principal incidental directions clustered from the $n$-way-
junction (fig. 4). In this case, the total number of required bisections
equals to $n-3$, for $n$ incident branches. At each bisection, the initial
parent loop is split into a pair of child loops following a binary spatial
partition by connecting any alternating pair of vertices of the parent loop
with an arc, and by splitting the connecting arc at the midpoint. This binary
spatial partition is iterated for the newly generated child loops until
convergence, where each child loop surrounds every incidental direction of the
$n$-way-junction. At convergence, each loop of the simplex inherits a
quadrilateral profile.

The resulting spatial configuration of $\mathcal{Q}$ is data-driven and is
determined upon the incidental directions and sizes of the branches at the
junction. Each bisection is based on maximising the volumetric partition for
each incident branch.

A set of interfacing quadrants is obtained for each loop in $\mathcal{Q}$, by
first connecting all the vertices of the simplex with its centroid, then by
splitting each arc of each loop at the midpoint (fig. 4). This determines a
consistent quadrant-based configuration and an initial geometrical embedding
of the incident elements. In particular, for each interface facing $w$, the
subset of vertices belonging to the quadrant-based simplex coincides with the
a subset of control points in eq. 1.

##### Conforming Constraints

The conforming condition of the considered scaffolding introduces a certain
dependency among the parameters of each vectorial element in the lattice. In
line with the construction, the parametric components along $u$ and $v$
directions are equivalent, i.e. the the univariate basis functions degrees
$d_{u}=d_{v}$ in eq. 2, as well as, the knot vectors
$\mathbf{k}_{u}=\mathbf{k}_{v}$, and the cardinality of control points
$l_{u}=l_{v}$ as in eq. 1 and eq. 2, respectively for each vectorial element.
Geometrically, the control points coordinates associated to each adjacent side
of any pair of conforming vectorial elements coincide as well. Each element of
the luminal scaffolding has independent parametric components along $w$, i.e.
arbitrary degree $d_{w}$, arbitrary spacing of the knot vector
$\mathbf{k}_{w}$ and arbitrary cardinality of control points $l_{w}$.

##  4\. Geometrical Embedding

The conforming scaffolding associates a 3D geometrical embedding to the
lattice of control points in eq. 1. Each free-form element of the scaffolding
recovers both shape and size of the underlying portion of the structure,
similarly to (Krishnamurthy and Levoy, 1996). The directional formulation of
the interfacing quadrilateral junction simplex in section 3.2 provides an
initial directional embedding, whereas the geometrical extent and
configuration of each branch of the scaffolding is detailed in the following
sections, by viable fitting strategies for different input data.

###  4.1. Curvilinear Primitives and Lattice

The geometrical embedding of each vectorial element is determined by a solid
profiling operator. This is controlled by a cross-sectional quadrilateral loop
of varying shape sweeping along the branch-wise skeletal curvilinear
primitive. In particular, the quadrant cross-sectional quadrilateral loop
comprises the curvilinear profiles at the interfacing simplexes and the
outlines estimated from the superficial boundary of the structure and the
underlying skeleton, at the extremities and along each branch, respectively.
Internal control points of the solid lattice are derived with a composition of
blended interpolations from such curvilinear primitives (Piegl and Tiller,
1996).

####  4.1.1. Skeletal Primitives

\begin{overpic}[width=429.28616pt]{fig5R.png} \put(2.0,21.0){\small{(a)}}
\put(32.0,21.0){\small{(b)}} \put(57.0,21.0){\small{(c)}}
\put(81.0,21.0){\small{(d)}}
\put(12.0,8.0){\small{$\textrm{C}_{\textrm{skel}}(w)$}}
\put(1.0,8.0){\small{$\hat{\mathbf{e}}_{i}^{(a)}$}}
\put(3.0,0.0){\small{$\hat{\mathbf{f}}_{i}^{(a)}$}}
\put(7.0,2.0){\small{$\hat{\mathbf{n}}_{i}^{(a)}$}}
\put(22.0,13.0){\small{$\hat{\mathbf{e}}_{i}^{(b)}$}}
\put(26.5,21.5){\small{$\hat{\mathbf{f}}_{i}^{(b)}$}}
\put(27.0,17.5){\small{$\hat{\mathbf{n}}_{i}^{(b)}$}} \end{overpic} Figure 5.
Curvilinear primitive: (a) Skeletal curve $\textrm{C}_{\textrm{skel}}(w)$:
curvature-preserving (magenta) and uniformly spaced (cyan) samples, with
terminal orthogonal bases; (b) Normal $\hat{\mathbf{e}}$ (orange) and binormal
$\hat{\mathbf{f}}$ (purple) vectors matching the terminal quadrant bases for
$\theta$; (c) Minimal torsion scaffolding and control points; (d) Scaffolding
between two junctions (Phantom) matching the quadrant bases for $\theta$.

The underlying skeletal centerline of the $i$-th branch is parametrised as a
NURBS curve $\textrm{C}_{\textrm{skel}}(w)$, conforming to the univariate
formulation in eq. 1 along the longitudinal direction $w$. The control points
and the associated knot vector are estimated following (Piegl and Tiller,
1996; Ma and Kruth, 1998; Pagani and Scott, 2018), considering an adaptive
sampling scheme accounting for uniform-, or non-uniform and curvature-
preserving samples, up to a desired spatial resolution (fig. 5). The
curvilinear primitive $\textrm{C}_{\textrm{skel}}(w)$ links either a pair of
junctions, or a junction with a terminal node, and each point of
$\textrm{C}_{\textrm{skel}}(w)$ is associated with a set of orthogonal bases,
comprising a tangential unit vector $\hat{\mathbf{n}}_{i}$ and a pair of
normal $\hat{\mathbf{e}}_{i}$ and binormal $\hat{\mathbf{f}}_{i}$ unit
vectors, as in a Frenet–Serret apparatus (Hanson, 1994). Both normal and
binormal bases along the skeletal primitive determine the directional
embedding of the branch-wise scaffolding quadrants.

At a junction simplex $\mathcal{Q}$, the normal $\hat{\mathbf{e}}_{i}$ and
binormal $\hat{\mathbf{f}}_{i}$ unit vectors are defined for the $i$-th branch
as

(8)
$\begin{matrix}\hat{\mathbf{e}}_{i}=\frac{\hat{\mathbf{x}}_{i}+\hat{\mathbf{x}}_{i}^{\perp}}{|\hat{\mathbf{x}}_{i}+\hat{\mathbf{x}}_{i}^{\perp}|}&\textrm{and}&\hat{\mathbf{f}}_{i}=\frac{\hat{\mathbf{y}}_{i}+\hat{\mathbf{y}}_{i}^{\perp}}{|\hat{\mathbf{y}}_{i}+\hat{\mathbf{y}}_{i}^{\perp}|},&\textrm{with}\end{matrix}$
(9)
$\begin{matrix}\hat{\mathbf{x}}_{i}^{\perp}=\frac{\hat{\mathbf{x}}_{i}-(\hat{\mathbf{x}}_{i}\cdot\hat{\mathbf{y}}_{i})\hat{\mathbf{y}}_{i}}{|\hat{\mathbf{x}}_{i}-(\hat{\mathbf{x}}_{i}\cdot\hat{\mathbf{y}}_{i})\hat{\mathbf{y}}_{i}|}&\textrm{and}&\hat{\mathbf{y}}_{i}^{\perp}=\frac{\hat{\mathbf{y}}_{i}-(\hat{\mathbf{x}}_{i}\cdot\hat{\mathbf{y}}_{i})\hat{\mathbf{x}}_{i}}{|\hat{\mathbf{y}}_{i}-(\hat{\mathbf{x}}_{i}\cdot\hat{\mathbf{y}}_{i})\hat{\mathbf{x}}_{i}|},\end{matrix}$

where, $\hat{\mathbf{x}}_{i}$ and $\hat{\mathbf{y}}_{i}$ are the unit vectors
orthogonal to the tangential $\hat{\mathbf{n}}_{i}$, respectively joining the
local pair of extraordinary vertices of the quadrilateral loop and the
opposite pair of loop’s midpoints. The normal and binormal bases of a terminal
node follow the Frenet–Serret formulas (Hanson, 1994) starting from the bases
at the junction simplex as in eq. 8, and osculate with minimal intrinsic
torsion.

By directly matching the the normal and binormal bases of two connected
junctions, the parametrisation of the $i$-th branch results in an additive
structural twist. To reduce such structural twist, the bases are quadrant-wise
matched by minimising a torsional angle and by diluting it along the skeletal
primitive (fig. 5).

Assuming $\hat{\mathbf{n}}_{i}^{(a)}$, $\hat{\mathbf{e}}_{i}^{(a)}$,
$\hat{\mathbf{f}}_{i}^{(a)}$ and $\hat{\mathbf{n}}_{i}^{(b)}$,
$\hat{\mathbf{e}}_{i}^{(b)}$, $\hat{\mathbf{f}}_{i}^{(b)}$ being the
orthogonal bases at the endpoints of $\textrm{C}_{\textrm{skel}}(w)$
underlying the $i$-th connecting branch, at the corresponding junctions $a$
and $b$, the torsional angle $\theta$ is given by

(10)
$\begin{matrix}\theta=\min\Bigg{(}\arccos{\Bigg{(}\hat{\mathbf{e}}_{i}^{(a)\perp}\cdot\underbrace{\left(\mathrm{R}_{\hat{\mathbf{n}}_{i}^{(b)}}\left(\frac{j\pi}{2}\right)\hat{\mathbf{e}}_{i}^{(b)}\right)}_{\big{\\{}\hat{\mathbf{e}}_{i}^{(b)},\hat{\mathbf{f}}_{i}^{(b)},-\hat{\mathbf{e}}_{i}^{(b)},-\hat{\mathbf{f}}_{i}^{(b)}\big{\\}}}\Bigg{)}}\Bigg{)},&\textrm{for}&j=\\{0,1,2,3\\},\end{matrix}$

where $\hat{\mathbf{e}}_{i}^{(a)\perp}$ is the orthogonal projection of
$\hat{\mathbf{e}}_{i}^{(a)}$ on $\hat{\mathbf{n}}_{i}^{(b)}$, and the rotation
matrix $\mathrm{R}_{\hat{\mathbf{n}}_{i}^{(b)}}(\frac{j\pi}{2})$ accounts for
quadrant rotations of the bases $\hat{\mathbf{e}}_{i}^{(b)}$,
$\hat{\mathbf{f}}_{i}^{(b)}$ around $\hat{\mathbf{n}}_{i}^{(b)}$.

Leveraging the symmetry of the branch-wise quadrants, the maximal structural
torsion along the skeletal primitive corresponds to
$\theta_{\max}=\frac{\pi}{4}$.

####  4.1.2. Cross-Sectional Quadrilateral Loops

\begin{overpic}[width=429.28616pt]{fig6R.png} \put(1.0,16.0){\small{(a)}}
\put(18.0,16.0){\small{(b)}} \put(45.0,16.0){\small{(c)}}
\put(81.0,16.0){\small{(d)}} \put(24.5,4.0){\small{$\textrm{C}_{0}(u)$}}
\put(33.0,14.0){\small{$\textrm{C}_{0}(v)$}}
\put(39.0,6.5){\small{$\textrm{C}_{1}(u)$}}
\put(34.0,0.5){\small{$\textrm{C}_{1}(v)$}} \end{overpic} Figure 6. Cross-
Sectional Quadrilateral Loops: (a) Projected simplex(es) - gold arcs, white
nodes on the surface of Genus3; (b) Qruadrant-wise oriented parametrisation of
boundary curves at the junction: $\textrm{C}_{0}(u)$ and $\textrm{C}_{0}(v)$
in a pseudo-radial direction, and $\textrm{C}_{1}(u)$ and $\textrm{C}_{1}(v)$
in a pseudo-circumferential direction, respectively; (c) Qruadrant-wise
parametrisation of boundary curves (lilac) for the elongation of the branch
delimited by longitudinal seam-cuts (cerulean). (d) Bilinearly blended Coons
patch delimited by the boundary curves: the interpolated cross-sectional
control points determine the internal lattice - splines with cubic degree(s).

The projected image of a junction simplex $\mathcal{Q}$ on the superficial
boundary of the structure determines the curvilinear profiles of the cross-
sectional quadrilateral loops in each respective quadrant. In particular, each
pair of extraordinary vertices and loops’ midpoints of $\mathcal{Q}$ is first
mapped following the directional formulation in section 3.2. Then, each arc of
the simplex loops is delineated with geodesics (Kimmel and Sethian, 1998;
Martinez et al., 2004; Surazhsky et al., 2005; Crane et al., 2020).

The parametrisation of each geodesic arc is constrained by the conforming
conditions along the parametric directions $u$ and $v$. For a luminal
scaffolding, each geodesic arc of a quadrant is first split in two segments at
the midpoint, and each segment is parametrised along $u$ and $v$ respectively,
in an alternating and opposite scheme (fig. 6). This determines a pair of
associated curvilinear primitives $\textrm{C}_{1}(u)$ and $\textrm{C}_{1}(v)$
for each arc of the quadrant, following the pseudo-circumferential arrangement
as in section 3.2. The inner part comprises pairs of curvilinear primitives
along $u$ and $v$ fitting the segments connecting each pair of extraordinary
vertices and each pair of loop’s midpoints to the junction node itself. This
determines another pair of curvilinear primitives, namely $\textrm{C}_{0}(u)$
and $\textrm{C}_{0}(v)$, originating from the junction node, each underlying a
pseudo-radial direction. For each quadrant, the 4-tuple of curvilinear
primitives $\textrm{C}_{0}(u)$, $\textrm{C}_{0}(v)$, $\textrm{C}_{1}(u)$ and
$\textrm{C}_{1}(v)$ constitute the sides of a quadrilateral loop on $(u,v)$ as
in fig. 6.

In a similar manner, cross-sectional quadrilateral loops are determined for
the elongation of each branch (fig. 6), where consecutive 4-tuple of
curvilinear primitives are sampled orthogonal to the longitudinal direction
$w$. Seam-cuts conforming to the skeletal primitive
$\textrm{C}_{\textrm{skel}}(w)$ connect the projected images of the simplexes
along the branch as in the underlying graph (fig. 6). At the outer part of the
quadrant along the branch, similar curvilinear geodesic arcs are determined by
leveraging a linear combination of the normal $\hat{\mathbf{e}}_{i}$ and
binormal $\hat{\mathbf{f}}_{i}$ directions along
$\textrm{C}_{\textrm{skel}}(w)$. The orthogonal bases are (non-)linearly
mapped onto the superficial boundary of the structure, and the associated
projections constrain the local geodesics outlines delineation.

####  4.1.3. Internal lattice

The internal lattice of each quadrant, and of each vectorial element, is
determined with a composition of blended interpolations from the cross-
sectional quadrilateral loops along the skeletal primitive. In particular,
each 4-tuple of curvilinear primitives $\textrm{C}_{0}(u)$,
$\textrm{C}_{0}(v)$, $\textrm{C}_{1}(u)$ and $\textrm{C}_{1}(v)$ encloses a
bilinearly blended Coons patch (Farin and Hansford, 1999), whose conforming
parametrisation on $(u,v)$ accounts for the interpolated inner control points
coordinates of the lattice (fig. 6). The conforming parametrisation of the
internal lattice is given by orderly stacking and further interpolating the
cross-sectional Coons patches along $\textrm{C}_{\textrm{skel}}(w)$, with a
profiling operator as in (Piegl, 1991; Piegl and Tiller, 1996).

###  4.2. Geodesic Fitting Strategies

Curvilinear arcs are determined using geodesics by connecting two end-points
on the superficial boundary of the branching structure. Fitting and
parametrising the curvilinear primitives follows (Piegl and Tiller, 1996; Ma
and Kruth, 1998; Pagani and Scott, 2018), in keeping with the conforming
constraints.

##### Polygonal Surface Meshes

Figure 7. Geodesic Fitting – Gallery of Polygonal Surface Meshes: mapping the
geometrical embedding of the scaffolding in quad-based segmentation seam-cuts.
Junction simplexes and terminal nodes: (gold), elongated branches (cerulean)
and sampling outlines (lilac).

Superficial outlines underlie pseudo-circumferential and longitudinal profiles
along each branch. Geodesic curves are extracted with state-of-the-art methods
as in (Surazhsky et al., 2005; Crane et al., 2020) leveraging exact or
approximated level-sets on a tessellated polygonal mesh. The directional
embedding of the scaffolding is first projected on the polygonal surface mesh
by means of a non-linear inverse skteleton-to-mesh mapping as in (Au et al.,
2008), or by semi-automatic user adjustments. The projected images of each
junction simplex nodes constitute the geodesic end-points. Along each branch,
the orthogonal bases of the Frenet–Serret apparatus elicit preferential
diffusion directions on the mesh, which constrain the geodesics along
directional fields (Panozzo et al., 2014; Pietroni et al., 2016) in quad-
regularised segmentation seam-cuts and outlines (fig. 7).

##### Volumetric Imaging

\begin{overpic}[width=390.25534pt]{fig8R.png} \put(-4.0,32.0){\small{(a)}}
\put(31.5,32.0){\small{(b)}} \put(55.5,32.0){\small{(c)}}
\put(77.0,32.0){\small{(d)}} \end{overpic} Figure 8. Geodesic Fitting –
Anatomical Structures (Phantom): mapping the geometrical embedding of the
scaffolding on snakes. (a) Volumetric image segmentation with snakes (lilac);
(b) Cross-sectional curve-to-planar projections along each branch and contour
delineation; (c) Scaffolding branch fitting on cross-sectional snakes; (d)
Junction simplex fitting on overlapping snakes point-cloud.

Geodesic active contours (Cheng et al., 2015), i.e. snakes, delineate the
superficial boundary of a structure from 3D volumetric imaging. In particular,
snakes contour pseudo-circumferential profiles along the longitudinal axis of
each branch (fig. 8). The volumetric image is first cross-sectionally
resampled with a series of curve-to-planar projections (Kanitsar et al., 2002)
forming a branch-wise stack of slices orthogonal to the centerline. Snakes
evolve on each slice minimising a level-set potential (Sethian, 1996, 1999) in
the correspondence of the boundary of the structure. For the elongation of the
branch, the directional embedding of the scaffolding is linearly projected
onto the snakes coordinates, directly fitting the curvilinear arcs (fig. 8).
At the junctions, the directional embedding of the simplexes fits a 3D point-
cloud determined by the overlapping snakes with least-squares (fig. 8).

##  5\. Smoothing Paradigm

The parametric and geometrical continuity of each scaffolding element depends
on the degrees of the basis functions and on the multiplicity of each knot in
the associated knot vector. In particular, referring to eq. 3 and considering
a generic parametric direction of a vectorial element, the number of
continuous derivatives at the knot $k$ equals to $d-\textrm{m}_{(k)}$, with
$d$ and $\textrm{m}_{(k)}$ being the univariate B-spline degree and the
multiplicity of the knot $k$, respectively.

Considering open knot vectors as in eq. 3, the parametric and geometrical
continuity of any pair of conforming elements of the scaffolding is only
positional, i.e. the conforming elements are mutually $\textrm{G}^{0}$.

This may result in discontinuous curvature, geometrical creases and visible
sharp edges, as well as in differential discontinuities of computational
simulation profiles at the contiguous boundaries and across interfaces of the
scaffolding. Similarly to isoparametric Hermite elements (Petera and Pittman,
1994), a solid smoothing paradigm is devised for the scaffolding, achieving
higher smoothness and continuity, typical of an organic structure as a
continuum.

The solid smoothing paradigm first builds on a topological subdivision of the
scaffolding control points into a set of nested sub-lattices, i.e. medial
shells. Then, the geometrical embedding of each sub-lattice is adjusted with a
smoothing scheme leveraging concepts from state-of-the-art smoothing
algorithms employed in subdivision surface modelling.

###  5.1. Medial Shells Subdivision

\begin{overpic}[width=429.28616pt]{fig9R.png} \put(10.0,29.0){\small{(a)}}
\put(63.5,29.0){\small{(b)}} \put(63.5,13.5){\small{(c)}}
\put(70.0,28.5){\small{$\mathcal{M}_{2}$}}
\put(81.5,28.5){\small{$\mathcal{M}_{1}$}}
\put(91.0,28.5){\small{$\mathcal{M}_{0}=\mathcal{M}_{\tilde{\tau}}$}}
\put(70.0,13.5){\small{$\mathcal{F}_{(\mathcal{M}_{2})}$}}
\put(81.0,13.5){\small{$\mathcal{F}_{(\mathcal{M}_{1})}$}}
\put(93.0,13.5){\small{$\mathcal{F}_{(\mathcal{M}_{\tilde{\tau}})}$}}
\end{overpic} Figure 9. Medial Shells Subdivision – Single Cuboid Element: (a)
Depth-wise colour-coded lattice of control points for each closed
$\mathcal{M}_{\tau}$. (b) Individual medial shell decomposition down to the
medial locus $\mathcal{M}_{\tilde{\tau}}$. (c) Minimal fringe
$\mathcal{F}_{(\mathcal{M}_{\tau})}$ of each medial shell: corner-lattice
vertices ($\Diamond$) and edge-lattice vertices ($\bigtriangleup$).

A topological subdivision of the scaffolding control points determines a set
of nested sub-lattices at different depths, i.e. a set of medial shells (figs.
9 and 10). A medial shell $\mathcal{M}$ is a lattice graph of regular tiling,
comprising a subset of control points $P$ and a subset of edges $E$ connecting
points of the same shell in quadrilateral cells. Such medial shell can also be
regarded as a quadrilateral (open or closed) mesh. A set of medial shells
(fig. 9) can be defined over the structured lattice of control points in an
individual element as

(11)
$\begin{matrix}\mathcal{M}_{\tau}=\\{(P,E)|P=\\{\mathbf{P}_{i,j,k}|\mathrm{t}_{i,j,k}=\tau\\},E=\mathcal{N}_{6}(P)\\},&\mathrm{with}&\tau=0,1,\dots,\bar{\tau}.\end{matrix}$

The scalar value $\mathrm{t}_{i,j,k}\in\mathbb{N}$ is the depth tag of each
control point $\mathbf{P}_{i,j,k}$ in the structured lattice. All control
points with same depth tag $\mathrm{t}_{i,j,k}=\tau$ constitute the vertices
of the associated medial shell $\mathcal{M}_{\tau}$, and the respective edges
are defined by the 6-connected neighbourhood $\mathcal{N}_{6}$ pattern, as in
the original structured lattice of the element. The depth-wise medial shell
subdivision follows a recursive definition of the tag $\mathrm{t}_{i,j,k}$.
All tags are first initialised as infinite values for the entirety of the
structured lattice, whereas $0$ values are imposed for the boundary sides. A
recursive scheme assigns a finite value to each tag in the lattice as

(12)
$\textrm{t}_{i,j,k}=\min\left(~{}\textrm{t}_{i,j,k}~{},~{}\min\left(~{}\textbf{t}_{\mathcal{N}_{26}(i,j,k)}~{}\right)+1~{}\right),$

where $\textbf{t}_{\mathcal{N}_{26}(i,j,k)}$ is the array of tag values
associated to the 26-connected neighbourhood $\mathcal{N}_{26}$ relative to
the index tuple $(i,j,k)$ in the structured lattice of each (multi-patch)
element.

Multiple nested shells are generated from the superficial boundary of the
lattice, i.e at null depth, following an iterative integer-step shrinking of
the exterior shell towards the innermost region of the structure. In
particular, each inner shell represents a discrete contraction of the relative
outer shell, up to an ultimate irreducible shell $\mathcal{M}_{\tilde{\tau}}$,
i.e. the medial locus, which lies in the neighbourhood of the structure
skeleton and can be a degenerate lattice graph.

\begin{overpic}[width=429.28616pt]{fig10R.png} \put(20.0,30.0){\small{(a)}}
\put(51.0,30.0){\small{(b)}} \put(51.0,14.0){\small{(c)}}
\put(58.5,30.0){\small{$\mathcal{M}_{3}$}}
\put(70.0,30.0){\small{$\mathcal{M}_{2}$}}
\put(81.5,30.0){\small{$\mathcal{M}_{1}$}}
\put(90.5,30.0){\small{$\mathcal{M}_{0}=\mathcal{M}_{\tilde{\tau}}$}}
\put(57.5,13.5){\small{$\mathcal{F}_{(\mathcal{M}_{3})}$}}
\put(69.5,13.5){\small{$\mathcal{F}_{(\mathcal{M}_{2})}$}}
\put(81.0,13.5){\small{$\mathcal{F}_{(\mathcal{M}_{1})}$}}
\put(92.5,13.5){\small{$\mathcal{F}_{(\mathcal{M}_{\tilde{\tau}})}$}}
\end{overpic} Figure 10. Medial Shells Subdivision – Multi-patch 3-way-
junction: (a) Depth-wise colour-coded lattice of control points for a section
of closed $\mathcal{M}_{\tau}$ (exploded view). (b) Individual medial shell
decomposition down to the medial locus $\mathcal{M}_{\tilde{\tau}}$. (c)
Composite fringe $\mathcal{F}_{(\mathcal{M}_{\tau})}$ of each medial shell:
corner-lattice vertices ($\Diamond$), edge-lattice vertices
($\bigtriangleup$), and face-lattice vertices ($\square$).

###  5.2. Medial Shells Smoothing

Given the analogy between medial shells and quadrilateral meshes, the
smoothing paradigm leverages concepts introduced by state-of-the-art
subdivision surface algorithms (Catmull and Clark, 1978; Stam, 2001).
Reference literature generally employs first a local refinement of the mesh
with local subdivisions. Then, vertices coordinates are adjusted according to
a smoothing scheme based on the vertex valence $\nu$, i.e. the number of edges
that meet at each vertex. In a similar manner, a lattice refinement by knot-
insertion can be alternated to the following smoothing paradigm at each
smoothing step, reaching the underlying limit geometry (fig. 13).

####  5.2.1. Medial Shell Fringe

\begin{overpic}[width=429.28616pt]{fig11R.png} \put(19.0,39.0){\small{(a)}}
\put(38.0,39.0){\small{(b)}} \put(38.0,25.0){\small{(c)}}
\put(38.0,11.0){\small{(d)}} \end{overpic} Figure 11. Closed Medial Shell
Fringe – Multi-patch 3-way-junction: (a) Parametric embedding and exterior
shell fringe - selected vertices (magenta outlines); (b) Neighbourhoods of
corner-lattice vertices ($\Diamond$): edge-connected neighbours (green), face-
connected neighbours (blue); (c) Neighbourhoods of face-lattice vertices
($\square$): edge-connected neighbours (green) and face-connected neighbours
(blue). (d) Neighbourhoods of edge-lattice vertices ($\bigtriangleup$): edge-
connected neighbours as in eq. 15 (green).

The smoothing paradigm operates on each medial shell fringe $\mathcal{F}$ of
the (sub-)lattice(s). The fringe is defined as the union of the boundary rims
of a lattice within a certain outskirt of influence. By way of example, the
outermost and minimal-influence fringe of a solid scaffolding coincides with
the set of edges belonging to the boundary sides of each vectorial element.

Considering the topological subdivision in section 5.1, the definition of a
shell-wise fringe $\mathcal{F}_{(\mathcal{M}_{\tau})}$ is applied to each
contracted medial shell, in keeping with the parametric arrangement of the
scaffolding elements. In detail, the medial shell fringe
$\mathcal{F}_{(\mathcal{M}_{\tau})}$ includes the vertices belonging to the
boundaries rims of the contracted sub-lattice. These are labelled either as
corner-, or edge-lattice vertices (or face-lattice vertices, with increasing
outskirt of influence). Additional vertices at the rims of the scaffolding
interfaces are also included in the shell-wise fringe
$\mathcal{F}_{(\mathcal{M}_{\tau})}$, these all being labelled as edge-lattice
vertices for a minimal outskirt of influence, or labelled as face-lattice
vertices when they are corners of the respective individual adjacent elements.
Note that all the extraordinary vertices of a simplex $\mathcal{Q}$ in section
3.2 coincide with face-lattice vertices of the fringe.

##### Smoothing Scheme

The smoothing scheme determines a geometric dual configuration for the
vertices in each fringe $\mathcal{F}_{(\mathcal{M}_{\tau})}$. Similarly to
(Stam, 2001), in the simplest form, the smoothing scheme adjusts the
coordinates of each vertex by adopting an iterative weighting scheme based on
the vertex valence $\nu$ and on the label of the vertex.

For each corner- or face-lattice vertex (fig. 11), the geometrical coordinates
are iteratively adjusted as

(13)
$\mathbf{P}_{i+1}=\alpha(\nu)\mathbf{P}_{i}+\sum_{\nu}~{}\beta(\nu)\mathbf{P}_{i}^{(\mathrm{e}_{\nu})}+\sum_{\nu}~{}\gamma(\nu)\mathbf{P}_{i}^{(\mathrm{f}_{\nu})},$

where $\mathbf{P}_{i}$ represents the coordinates of the considered corner- or
face-lattice vertex of the fringe $\mathcal{F}_{(\mathcal{M}_{\tau})}$,
$\mathbf{P}_{i}^{(\mathrm{e}_{\nu})}$ represents the coordinates of the
respective $\nu$-th edge-wise connected neighbour in $\mathcal{M}_{\tau}$, and
$\mathbf{P}_{i}^{(\mathrm{f}_{\nu})}$ represents the coordinates of the
respective $\nu$-th face-wise connected neighbour in $\mathcal{M}_{\tau}$,
respectively at the $i$-th iteration. The weights are defined by the valence
$\nu$ of the fringe vertex as

(14)
$\displaystyle\begin{matrix}~{}~{}~{}\alpha(\nu)=\frac{\nu-3}{\nu}~{}~{}~{}&~{}~{}~{}\beta(\nu)=\frac{2}{\nu^{2}}~{}~{}~{}&~{}~{}~{}\gamma(\nu)=\frac{1}{\nu^{2}},~{}~{}~{}\end{matrix}$

and satisfy an affine invariant scheme, i.e.
$\alpha(\nu)+\nu~{}\beta(\nu)+\nu~{}\gamma(\nu)=1$.

In this case, corner-lattice vertices of the fringe usually have valence
$\nu=3$, whereas face-lattice vertices of the fringe are associated either to
regular vertices ($\nu=4$), or to extraordinary vertices with arbitrary
valence $\nu>4$.

Conversely, for each edge-lattice vertex (fig. 11), the geometrical
coordinates are iteratively adjusted as

(15)
$\mathbf{P}_{i+1}=\tilde{\alpha}\mathbf{P}_{i}+\sum_{\nu}~{}\tilde{\beta}\mathbf{P}_{i}^{(\mathrm{e}_{\nu}^{\perp})},$

where, at the $i$-th iteration, $\mathbf{P}_{i}$ represents the coordinates of
the considered edge-lattice vertex of the fringe
$\mathcal{F}_{(\mathcal{M}_{\tau})}$, and
$\mathbf{P}_{i}^{(\mathrm{e}_{\nu}^{\perp})}$ represents the coordinates of
the $\nu$-th edge-wise connected neighbour in $\mathcal{M}_{\tau}$, whose
parametric direction in the (sub-)lattice lies perpendicularly to the one of
the considered fringe vertex. In other words, each edge-wise connected
neighbour in $\mathcal{M}_{\tau}$ belonging to the same rim in the
(sub-)lattice is excluded. The latter condition imposes each edge-lattice
vertex of the fringe has valence $\nu=2$. The weights, in this case, are fixed
as $\tilde{\alpha}=\frac{1}{2}$ and $\tilde{\beta}=\frac{1}{4}$.

The iterative smoothing scheme ultimately converges for a machine error
evaluated on the norm $|\mathbf{P}_{i+1}-\mathbf{P}_{i}|$.

##### Smoothing Scheme of Open Medial Shells

\begin{overpic}[width=429.28616pt]{fig12R.png} \put(0.0,35.0){\small{(a)}}
\put(30.0,35.0){\small{(b)}} \put(60.0,14.0){\small{(c)}}
\put(81.0,14.0){\small{(d)}} \put(15.0,33.0){\small{$\mathcal{M}_{\tau}$}}
\put(44.0,33.0){\small{$\mathcal{F}_{(\mathcal{M}_{3})}$}}
\put(65.0,33.0){\small{$\mathcal{F}_{(\mathcal{M}_{2})}$}}
\put(79.0,33.0){\small{$\mathcal{F}_{(\mathcal{M}_{1})}$}}
\put(92.5,33.0){\small{$\mathcal{F}_{(\mathcal{M}_{\bar{\tau}})}$}}
\end{overpic} Figure 12. Open Medial Shells Fringe – Multi-patch 3-way-
Junction: (a) Depth-wise colour-coded lattice for an open $\mathcal{M}_{\tau}$
(section) - open-fan (red); (b) Medial shell fringe decomposition; (c)
Neighbourhood of open-fan corner-lattice vertices ($\Diamond$), and (d) face-
lattice vertices ($\square$). Respective edge-connected neighbours as in eq.
16 (green). Edge-lattice vertices of the open-fan ($\bigtriangleup$) remain
unaltered.

Certain configurations of the medial shells result in vertices of the fringe
$\mathcal{F}_{(\mathcal{M}_{\tau})}$ having an open-fan in the quadrilateral
tiling. In this particular case, the smoothing scheme of the corner-, edge-
and face-lattice vertices belonging to the open-fan changes to account for the
open boundary of the lattice.

Each edge-lattice vertex of the open-fan is kept unaltered, whereas each
corner- and face-lattice vertex of the open-fan is treated as a particular
edge-lattice vertex of $\mathcal{M}_{\tau}$ (fig. 12), and the geometrical
coordinates are iteratively adjusted as

(16)
$\mathbf{P}_{i+1}=\tilde{\alpha}\mathbf{P}_{i}+\sum_{\nu}~{}\tilde{\beta}\mathbf{P}_{i}^{(\mathrm{e}_{\nu}^{\parallel})}.$

Here, $\mathbf{P}_{i}^{(\mathrm{e}_{\nu}^{\parallel})}$ represents the
coordinates of the $\nu$-th edge-wise connected neighbour in the open-fan of
$\mathcal{M}_{\tau}$. As in eq. 15, the resulting valence is $\nu=2$, and the
weights coincide to $\tilde{\alpha}=\frac{1}{2}$ and
$\tilde{\beta}=\frac{1}{4}$.

##### Smoothing Scheme of the Medial Locus

The medial locus $\mathcal{M}_{\tilde{\tau}}$ makes exception in the smoothing
scheme, as the fringe of a degenerate lattice graph may not be defined. For
each vertex of the medial locus, the original edge-wise 6-connected
neighbourhood ($\mathcal{N}_{6}$) of the conforming lattice is considered,
regardless of the medial shell of belonging. The associated geometrical
coordinates are iteratively adjusted as

(17)
$\mathbf{P}_{i+1}=\bar{\alpha}\mathbf{P}_{i}+\sum_{\nu}~{}\bar{\beta}(\nu)\mathbf{P}_{i}^{(\mathrm{e}_{\nu}^{\mathcal{N}_{6}})},$

where, at the $i$-th iteration, $\mathbf{P}_{i}$ represents the coordinates of
the considered vertex of the medial locus $\mathcal{M}_{\tilde{\tau}}$, and
$\mathbf{P}_{i}^{(\mathrm{e}_{\nu}^{\mathcal{N}_{6}})}$ represents the
coordinates of the $\nu$-th edge-wise connected neighbour in $\mathcal{N}_{6}$
of the conforming lattice. The weights are defined as
$\bar{\alpha}=\frac{1}{2}$ and $\bar{\beta}(\nu)=\frac{1}{2\nu}$, satisfying
an affine invariant scheme.

Open-fan vertices of the medial locus ($\nu<6$) are adjusted by considering
only the neighbours in $\mathcal{N}_{6}$ belonging to the same boundary side
of the conforming scaffolding.

\begin{overpic}[width=429.28616pt]{fig13R.png} \put(0.0,49.0){\small{(a)}}
\put(0.0,23.0){\small{(b)}} \end{overpic} Figure 13. Solid Smoothing Paradigm
at Convergence: (a) Single cuboid; (b) Multi-patch 4-way-junction with locally
mixed open and closed medial shells. Resulting geometry, depth-wise medial
shells of control points, reflection lines and inner solid trabecular view.

##  6\. Results

##### Input Data

A set of synthetic and CAD geometries with varying shape, complexity,
composite graph structure and branching patterns were obtained by
construction. Popular benchmarking geometries used in computer graphics
applications have been considered form (Zhou and Jacobson, 2016; Attene et
al., 2016) as reference surface meshes for the scaffolding parametrisation.
Also, representative multi-modal scans, i.e. computed tomography (CT) and
magnetic resonance (MR), from available medical repositories (Piccinelli et
al., 2009; Bullitt et al., 2005) were considered for tracing anatomical
structures from raster volumetric imaging.

##### Parametric and Implementation Details

All parametrisations are considered for solid vectorial elements of cubic
degree, i.e. $d_{u}=d_{v}=d_{w}=3$, with unitary weights of the rational basis
functions. Minimal number of control points $l_{u},l_{v},l_{w}\geq 4$ are
adopted by enforcing an adaptive spatial resolution, based on (Pagani and
Scott, 2018). For simplicity, and in keeping with the conforming constraints,
the parametrisations along $u$ and $v$ directions in the luminal scaffolding
employed either uniformly-spaced, or non-uniformly symmetrically-spaced knot
vectors ($\mathbf{k}_{u}$, $\mathbf{k}_{v}$). The technical implementation of
the introduced scaffolding framework builds on available tools from (De Falco
et al., 2011; Vázquez, 2016; Bingol and Krishnamurthy, 2019).

##### Parametrisation Evaluation

Features of the parametrised structures account for the total number of
scaffolding elements, the density of control points lattice, the computational
time, and the accuracy for fitted geometries.

For synthetic and CAD structures, the compactness of the proposed approach is
first evaluated by comparing our parametrisation with state-of-the-art solid
meshing techniques employed in computational simulations. For each structure,
an optimal tetrahedral mesh is generated using (Si, 2015), reporting the
number of quadratic-element tetrahedra, total number of vertices and
computational time.

Similarly, the scaffolding parametrisation is evaluated for polygonal meshes.
The fitting accuracy of the geometric embedding is reported considering the
input surface mesh as reference. High accuracy corresponds to minimal
deviation from the reference mesh, in particular, the accuracy is evaluated as
the complement of the percentage parametrisation error, i.e.
$\text{Acc.}=(1-\varepsilon)\cdot 100\%$, and the error $\varepsilon$ is
computed as the Euclidean (or Hausdorff: $\varepsilon_{95}$) distance between
the closest points on the boundary surface of the parametrised structure and
on the reference surface using (Aspert et al., 2002), normalised by the local
radius of the underlying inscribed sphere (Stolpner et al., 2011). For
representative structures, the proposed parametrisation is also compared
against the state-of-the-art linear-elements hexahedral meshing method (Livesu
et al., 2020).

The parametric accuracy of fitted anatomical structures is evaluated on the
normalised distance between the closest points on the parametrised boundary
surface and the sampled geodesic active contours as reference. A qualitative
comparison with seminal work by (Zhang et al., 2007; Urick et al., 2019) is
also provided.

Table 1. Compactness of synthetic shapes and CAD geometries: Solid NURBS Scaffolding: $|\mathrm{H_{3}}|$ total cubic-degree elements; $|\mathbf{P}|$ total control points in the multi-patch lattice; the computational time (s). Equivalent Solid Tetrahedral Meshes (Si, 2015): $|\mathrm{T_{2}}|$ total quadratic tetrahedra, $|\mathbf{V}|$ total vertices; the computational time (s). Coupled luminal ($L$) and wall ($W$) multi-compartments in subscript. | Solid NURBS Scaffolding | Solid Tetrahedral Mesh
---|---|---
Structure | $|\mathrm{H_{3}}|$ | $|\mathbf{P}|$ | Time | $|\mathrm{T_{2}}|$ | $|\mathbf{V}|$ | Time
3-Junct.L | 12 | 1.3k | 5.3 | 26.7k | 46.9k | 0.4
3-Junct.W | 24 | 2.7k | 5.3 | 40.0k | 80.4k | 0.7
5-Junct.L | 20 | 2.2k | 10.2 | 41.2k | 74.9k | 0.7
5-Junct.W | 40 | 4.5k | 10.2 | 68.2k | 135.5k | 1.3
7-Junct.L | 28 | 3.1k | 20.3 | 62k | 109.2k | 0.9
7-Junct.W | 56 | 6.3k | 20.3 | 97.6k | 191.8k | 1.9
9-Junct.L | 36 | 4.0k | 29.1 | 40.4k | 137.1k | 1.3
9-Junct.W | 72 | 8.0k | 29.1 | 128.0k | 249.2k | 2.8
Cross | 7 | 875 | 2.4 | 38.0k | 74.0k | 0.9
EggL | 4 | 256 | 1.4 | 6.2k | 12.0k | 0.1
EggW | 16 | 1.0k | 1.4 | 18.8k | 28.5k | 0.2
Frame | 8 | 1.0k | 2.7 | 40.0k | 78.5k | 1.0
Net 1 | 40 | 5.2k | 21.7 | 224.4k | 345.3k | 8.8
Net 2 | 72 | 8.4k | 58.5 | 317.6k | 569.9k | 21.9
Net 3 | 80 | 8.0k | 54.6 | 408.9k | 604.1k | 15.2
Pinocchio | 24 | 1.5k | 5.6 | 115k | 196.7k | 1.9
Quadball | 1 | 64 | 0.2 | 5.8k | 10.9k | 0.1
Torus | 4 | 448 | 0.3 | 10.5k | 21.1k | 0.2
Twist | 1 | 625 | 1.6 | 16.1k | 29.7k | 0.3
Air | 36 | 2.3k | 3.3 | 61.8k | 114.1k | 1.4
Blades | 48 | 3.1k | 3.3 | 64.3k | 125.3k | 1.3
Gear | 54 | 3.5k | 6.5 | 73.2k | 145.1k | 1.7
Hooks | 14 | 0.9k | 2.9 | 54.1k | 104.4k | 0.9
Plate | 72 | 4.6k | 5.1 | 73.1k | 145.7k | 1.6
Socket | 32 | 2.0k | 7.3 | 36.3k | 68.3k | 0.6
Stent | 120 | 11.5k | 29.7 | 580.8k | 789.9k | 10.1
T-pipe | 12 | 2.4k | 5.5 | 39.9k | 80.3k | 0.8

###  6.1. Synthetic and CAD Geometries

Figure 14. Gallery of Solid Synthetic and CAD Geometries – Minimal elements
scaffolding and modular decomposition in the explosion view; smooth and sharp
features are preserved for different medial shells configurations. Control
points lattice, parametric embedding, reflection lines and inner solid
trabecular and sectioning views.

Simple to complex shapes are modelled using the introduced scaffolding. More
complex geometries combine a composition of interfacing luminal and wall
scaffoldings to model holes and complementary regions, leveraging modular and
local blocks decomposition. Individual elements are assessed in the explosion
view (fig. 14), where minimal multi-patches exactly recover the geometrical
shape of large and complex geometries. Smooth and sharp features of the
structures are preserved by enforcing arbitrary smoothing configurations to
the medial shells. The continuum structure of the parametrised solids is shown
in the trabecular view (fig. 14), where free-form cuboid cells are sampled
from the vectorial domain. The compactness of the proposed solid
parametrisation accounts for dozens to hundreds of cubic-degree elements,
comprising few thousands of control points in the solid lattice for increasing
structure complexity (table 1). The compactness is compared against adaptive
quadratic-element tetrahedral meshing methods employed in computational
simulations (Si, 2015). To accurately represent the same structure,
tetrahedral meshes required a dense tessellation of elements and vertices,
which differs by several orders of magnitude compared to the proposed
parametrisation. Coarser tessellations are feasible (fig. 15), yet higher
deviation from the boundary structure is observed, introducing larger
geometrical errors and approximations for computational simulations. The
structures parametrised with our approach exhibit higher degree of smoothness
and parametric continuity. This is shown for the superficial boundary in the
reflection lines (Balzer and Werling, 2010), where nearly G2 smoothness is
found for organic and smoothly blended portions of the structures (fig. 14).
Also, structural parametric continuity is observed in the smooth profiles of
the trabecular view, where inner sections and structural cells of the
vectorial solid domains are arbitrarily sampled (fig. 14). Different smoothing
configurations can be defined by altering the medial shells subdivision for
mixed sharp and organically blended sub-portions, as in conventional smoothing
groups. The construction of fully-coupled geometries, e.g. $n$-way-junctions,
accounting for both luminal and wall multi-compartments are simultaneously
generated, determining a conforming multi-compartmental domain fig. 14.

\begin{overpic}[width=429.28616pt]{fig15R.png} \put(2.0,94.0){\small{(a)}}
\put(32.0,94.0){\small{(b)}} \put(66.0,94.0){\small{(c)}}
\put(20.5,92.0){\small{$|\mathrm{H_{3}}|:4$}}
\put(20.5,90.0){\small{$|\mathbf{P}|:448$}}
\put(20.5,54.0){\small{$|\mathrm{H_{3}}|:32$}}
\put(20.5,52.0){\small{$|\mathbf{P}|:2.0$k}}
\put(20.5,16.0){\small{$|\mathrm{H_{3}}|:12$}}
\put(20.5,14.0){\small{$|\mathbf{P}|:2.4$k}}
\put(53.0,92.0){\small{$|\mathrm{T_{2}}|:10.5$k}}
\put(53.0,90.0){\small{$|\mathbf{V}|:21.1$k}}
\put(53.0,88.0){\small{$\text{Acc.}\approx 100\%$}}
\put(53.0,54.0){\small{$|\mathrm{T_{2}}|:36.3$k}}
\put(53.0,52.0){\small{$|\mathbf{V}|:68.3$k}}
\put(53.0,50.0){\small{$\text{Acc.}\approx 100\%$}}
\put(53.0,16.0){\small{$|\mathrm{T_{2}}|:39.9$k}}
\put(53.0,14.0){\small{$|\mathbf{V}|:80.3$k}}
\put(53.0,12.0){\small{$\text{Acc.}\approx 100\%$}}
\put(86.0,92.0){\small{$|\mathrm{T_{2}}|:204$}}
\put(86.0,90.0){\small{$|\mathbf{V}|_{\mathrm{coarse}}:427$}}
\put(86.0,88.0){\small{$\text{Acc.}\,74.2\pm 17.4\%$}}
\put(86.0,86.0){\small{$\text{Acc.}_{\varepsilon_{95}}\,42.3\%$}}
\put(86.0,54.0){\small{$|\mathrm{T_{2}}|:1.5$k}}
\put(86.0,52.0){\small{$|\mathbf{V}|_{\mathrm{coarse}}:2.9$k}}
\put(86.0,50.0){\small{$\text{Acc.}\,87.6\pm 10.3\%$}}
\put(86.0,48.0){\small{$\text{Acc.}_{\varepsilon_{95}}\,67.2\%$}}
\put(86.0,16.0){\small{$|\mathrm{T_{2}}|:1.1$k}}
\put(86.0,14.0){\small{$|\mathbf{V}|_{\mathrm{coarse}}:2.2$k}}
\put(86.0,12.0){\small{$\text{Acc.}\,88.4\pm 10.3\%$}}
\put(86.0,10.0){\small{$\text{Acc.}_{\varepsilon_{95}}\,68.9\%$}}
\end{overpic} Figure 15. Compactness of the Parametrisation – representative
examples: (a) Solid NURBS scaffolding: exact geometrical domain - sectioning
inner trabecular view. (b) High-resolution equivalent tetrahedral mesh using
(Si, 2015) with adaptive sampling and quadratic-elements. (c) Coarse
quadratic-tetrahedral mesh optimised for
$|\mathbf{V}|_{\textrm{coarse}}\approx|\mathbf{P}|$. Superficial degradation
for increasing coarsening.

###  6.2. Polygonal Surface Meshes

Table 2. Fitting polygonal surface meshes: features of parametrised structures. Input surface geometry: $|\mathrm{T}|$ total triangles, $|\mathbf{V}|$ total vertices. Solid NURBS Scaffolding: $|\mathrm{H_{3}}|$ total cubic-degree elements; $|\mathbf{P}|$ total control points in the multi-patch lattice; the computational time (s); geometrical embedding fitting accuracy ($\%$) is reported mean$\pm$sd and relative to Hausdorff distance ($\varepsilon_{95}$). | Surface Mesh | Solid NURBS Scaffolding
---|---|---
Structure | $|\mathrm{T}|$ | $|\mathbf{V}|$ | $|\mathrm{H_{3}}|$ | $|\mathbf{P}|$ | Time | Acc. % | $\text{Acc.}_{\varepsilon_{95}}~{}\%$
Armadillo | 83.4k | 41.7k | 96 | 20.5k | 463.9 | 95.37$\pm$6.44 | 82.76
Botijo | 82.3k | 41.2k | 52 | 22.0k | 426.8 | 97.18$\pm$2.61 | 92.54
Bunny | 104.3k | 52.2k | 40 | 13.8k | 362.1 | 97.70$\pm$3.19 | 92.40
Elephant | 78.3k | 39.2 | 76 | 18.2k | 324.3 | 97.74$\pm$2.79 | 93.08
Elk | 108.3k | 54.1 | 48 | 17.2k | 401.4 | 98.05$\pm$2.53 | 93.19
Fertility | 483.2k | 241.6k | 32 | 11.4k | 862.1 | 98.67$\pm$1.64 | 95.68
Genus3 | 13.3k | 6.7k | 24 | 4.5k | 15.9 | 96.87$\pm$3.29 | 90.63
Hand | 73.2k | 36.6k | 28 | 10.8k | 206.8 | 97.83$\pm$3.78 | 93.08
Horse | 61.3k | 30.6k | 36 | 14.2k | 239.8 | 97.21$\pm$4.36 | 91.42
Kitten | 110.1k | 55.4k | 28 | 10.0k | 249.3 | 98.67$\pm$1.80 | 95.67
Figure 16. Gallery of Solid Parametrised Geometries from Polygonal Surface
Meshes – Explosion view, reflection lines and inner solid trabecular view.
Nearly G2 smoothness and parametric continuity overall for boundaries and
interfaces with organic profiles.

A set of benchmark polygonal meshes are parametrised as reference (fig. 16).
The luminal scaffolding recovers the underlying shape by fitting the mesh
geometries with high accuracy (table 2) with respect to the reference boundary
surface. Complex structures with multiple busy branching patterns are
compactly parametrised with minimal element saffoldings. Free-form cuboids of
the conforming scaffolding recover the geometrical embedding of the underlying
structure with localised deviations and minor distortions for high-curvature
regions. The spatial partitions of junctions are automatically determined for
busy skeletal branch-points. The automatic mapping of the scaffolding
initially produced asymmetric quadrants and projected irregular quad-based
segmentations on the boundary surface, especially for busy junctions with
severe changes in size and shape. The resulting smoothness and continuity of
the parametrised solid structures is shown in the reflection lines for the
outer boundary and in the trabecular view for the innermost region in fig. 16.
Nearly G2 smoothness is achieved globally, and high structural regularity of
the lattice is obtained in the form of an organic continuum. In the
neighbourhood of extraordinary vertices, the solid scaffolding exhibits G1
smoothness and continuity, whereas pure positional continuity (G0) is observed
exactly for the extraordinary vertices. Representative structures are
parametrised using state-of-the-art hexahedral meshing (Livesu et al., 2020)
in fig. 17. Higher compactness and higher accuracy is observed for the
proposed approach with respect to linear-element hexahedral meshes. Note that
a further NURBS-based parametrisation is required in (Livesu et al., 2020)
prior to employing the hexahedral mesh in IGA.

\begin{overpic}[width=429.28616pt]{fig17R.png} \put(0.0,90.0){\small{(a)}}
\put(31.0,90.0){\small{(b)}} \put(65.0,90.0){\small{(c)}}
\put(2.5,88.0){\small{$|\mathrm{T}|:110.1$k}}
\put(2.5,86.0){\small{$|\mathbf{V}|:55.4$k}}
\put(33.5,88.0){\small{$|\mathrm{H_{3}}|:28$}}
\put(33.5,86.0){\small{$|\mathbf{P}|:10$k}}
\put(33.5,84.0){\small{$\text{Acc.}\,98.67\pm 1.8\%$}}
\put(33.5,82.0){\small{$\text{Acc.}_{\varepsilon_{95}}\,95.67\%$}}
\put(67.5,88.0){\small{$|\mathrm{H_{1}}|:1.5$k}}
\put(67.5,86.0){\small{$|\mathbf{P}|:1.9$k}}
\put(67.5,84.0){\small{$\text{Acc.}\,96.7\pm 4.74\%$}}
\put(67.5,82.0){\small{$\text{Acc.}_{\varepsilon_{95}}\,88.54\%$}}
\put(23.0,36.0){\small{$|\mathrm{T}|:104.3$k}}
\put(23.0,34.0){\small{$|\mathbf{V}|:52.2$k}}
\put(55.0,36.0){\small{$|\mathrm{H_{3}}|:40$}}
\put(55.0,34.0){\small{$|\mathbf{P}|:13.8$k}}
\put(55.0,32.0){\small{$\text{Acc.}\,97.7\pm 3.19\%$}}
\put(55.0,30.0){\small{$\text{Acc.}_{\varepsilon_{95}}\,92.4\%$}}
\put(86.0,36.0){\small{$|\mathrm{H_{1}}|:2.2$k}}
\put(86.0,34.0){\small{$|\mathbf{P}|:2.8$k}}
\put(86.0,32.0){\small{$\text{Acc.}\,97.53\pm 4.12\%$}}
\put(86.0,30.0){\small{$\text{Acc.}_{\varepsilon_{95}}\,90.14\%$}}
\end{overpic} Figure 17. Comparison with Hexahedral Meshing: (a) Reference
triangle surface mesh; (b) Solid NURBS scaffolding; (c) Solid hexahedral mesh
using (Livesu et al., 2020). Quantitative attributes reporting the compactness
and the accuracy at the boundary. $|\mathrm{H}_{3}|$ and $|\mathrm{H}_{1}|$
respectively for cubic-degree and linear-degree total elements.

###  6.3. Anatomical Structures from Clinical Imaging

Figure 18. Gallery of Anatomical Structures Parametrised from Volumetric Scans. Contours, slice intersections, solid NURBS scaffolding with reflection lines and inner solid trabecular views. Scalable parametrisation to large anatomical structures. Table 3. Tracing anatomy from volumetric imaging: features of parametrised structures. Input scan: Imaging modality, $\mathbf{v}$ voxel size (mm). Solid NURBS Scaffolding: $|\mathrm{H_{3}}|$ total cubic-degree elements; $|\mathbf{P}|$ total control points in the multi-patch lattice; the computational time (s); geometrical parametrisation error (mm) is reported mean$\pm$sd. Hausdorff distance ($\varepsilon_{95}$). | Input Scan | Solid NURBS Scaffolding
---|---|---
Structure | Modality |  $\mathbf{v}$ (mm) | $|\mathrm{H_{3}}|$ | $|\mathbf{P}|$ | Time |  $\varepsilon$ (mm) |  $\varepsilon_{95}$ (mm)
Angio 1 | CT | 0.26 isotropic | 52 | 10.1k | 42.7 | 0.26$\pm$0.15 | 0.55
Angio 2 | CT | 0.26 isotropic | 36 | 10.7k | 54.3 | 0.23$\pm$0.13 | 0.48
Angio 3 | CT | 0.35 isotropic | 52 | 11.5k | 70.1 | 0.26$\pm$0.14 | 0.49
Angio 4 | CT | 0.35 isotropic | 44 | 10.6k | 62.2 | 0.34$\pm$0.21 | 0.71
Angio 5 | CT | 0.35 isotropic | 52 | 11.4k | 70.7 | 0.25$\pm$0.15 | 0.54
Angio 6 | CT | 0.35 isotropic | 36 | 11.4k | 73.5 | 0.28$\pm$0.19 | 0.61
Brain Net 1 | MR | [0.36$\times$0.36$\times$0.5] | 100 | 19.6k | 122.2 | 0.34$\pm$0.21 | 0.72
Brain Net 2 | CT | [0.5$\times$0.5$\times$0.7] | 76 | 29.5k | 597.8 | 0.41$\pm$0.25 | 0.87
Phantom | MR | [0.36$\times$0.36$\times$0.5] | 20 | 3.5k | 12.5 | 0.27$\pm$0.15 | 0.53

Anatomical structures from volumetric scans are parametrised using the
introduced scaffolding (fig. 18). The compact parametrisation produced a
minimal luminal scaffolding comprising exactly 4 elements per branch. Tubular
structures with arbitrary branching and connectivity patterns are accurately
recovered (table 3) for the elongation of the branches. Also, irregular and
convoluted tubular profiles as well as size-varying junctions are captured by
the solid scaffolding with an organic outline using an optimal smoothing and
fitting strategy. Nearly G2 smoothness and continuity is achieved globally, by
preserving local sharp features at the anatomical inlets and outlets allowing
for boundary conditions in bio-mechanic computational simulations. Similar
tubular structures are parametrised in previous work (Zhang et al., 2007;
Urick et al., 2019), where elongated portions employ hexahedral meshes with
dense cross-sectional subdivisions. In these studies, neither data-driven, nor
generic parametric form is given for bifurcations and busy junctions, which
are conversely modelled with a dictionary of user-defined branching templates.
In (Zhang et al., 2007), the parametrised anatomical structures reach
circumferential G1 smoothness and continuity, whereas junctions locally
exhibits sharp creases and longitudinal discontinuities.

###  6.4. Isogeometric Analysis Simulations

Representative computational problems are considered as proof-of-concept IGA
applications to validate the parametrised scaffolding. A linear elasticity
problem (Auricchio et al., 2007) is solved considering the Plate (CAD) and
Stent (CAD) structures, as well as the Bunny (Mesh) geometry. In this cases,
boundary conditions and simulated solid deformations are shown in fig. 19. A
steady-state fluid problem (Gomez et al., 2010) is solved for representative
luminal synthetic junctions, as well as for BrainNet 1 anatomical structure.
Boundary conditions and simulation solution profiles of pressure and flow
distributions are shown in fig. 19. Lastly, the Maxwell eigenvalue problem
(Vázquez and Buffa, 2010) is solved for Genus3 and Kitten (Mesh), as well as
for the EggW (Synthetic) structure. In this case a set of eigenfunctions and
the associated distributions are shown for both luminal and hollow wall
geometries in fig. 19. All simulations are performed using (De Falco et al.,
2011), and the computational time as well as the memory footprint are provided
in fig. 20. The scalability of IGA simulations and associated performance is
also evaluated on equivalent denser lattices. An iterative cascade of
conforming domain subdivisions (fig. 20) is progressively performed for the
considered geometries, achieving comparable density as in hexahedral meshes
(Livesu et al., 2020; Gao et al., 2017, 2019; Takayama, 2019).

Successful simulations converged for all our scaffolding parametrisations
(fig. 20), with smooth, continuous and consistent solution profiles. The
introduced framework exhibits low computational time and minimal memory
footprint. As in fig. 20, performance benchmarks shows a quadratic trend for
increasing density of the lattice, suggesting improved scalability for
structures parametrised with the introduced scaffolding. Equivalent hexahedral
meshes determined a progressive increase of elements per conforming
subdivision, resulting in a parametric lattice of substantially higher
density. In this case, simulations only converged for meshes with less than
1.5k elements, as the memory footprint determined the limiting factor for the
adopted IGA engine and solver. All simulations were computed on a single-core
CPU 3.60 GHz with 64GB RAM machine.

\begin{overpic}[width=429.28616pt]{fig19AR.png} \put(0.0,60.0){\small{(a)}}
\put(0.0,44.0){\small{(b)}} \put(35.0,18.0){\small{(c)}} \end{overpic} Figure
19. IGA Simulations using (De Falco et al., 2011): (a) Linear-elastic
deformation; (b) Fluid pressure and flow; (c) Set of Maxwell eigenfunctions.
\begin{overpic}[width=260.17464pt]{fig19BR.png} \put(2.0,98.0){\small{(a)}}
\put(2.0,67.0){\small{(b)}} \put(13.0,75.0){\small{$|\mathrm{H}|:1$}}
\put(13.0,72.0){\small{$|\mathbf{P}|:125$}}
\put(37.0,75.0){\small{$|\mathrm{H}|:8$}}
\put(37.0,72.0){\small{$|\mathbf{P}|:512$}}
\put(61.0,75.0){\small{$|\mathrm{H}|:64$}}
\put(61.0,72.0){\small{$|\mathbf{P}|:4$k}}
\put(85.0,75.0){\small{$|\mathrm{H}|:512$}}
\put(85.0,72.0){\small{$|\mathbf{P}|:32.8$k}} \end{overpic} Figure 20. IGA
Simulations using (De Falco et al., 2011) \- Benchmarking performance: (a)
Cascade of conforming subdivisions; (b) Computational time and memory
footprint as function of elements and control points density in the lattice;
our scaffolding (orange dots) vs equivalent hexahedral meshes (blue dots).
OoM: Out of Memory threshold (red dashed line).

##  7\. Discussions and Conclusion

This work introduced a novel scaffolding framework to trace 3D structures with
solid NURBS elements for computational simulations using isogeometric
analysis.

Results obtained form different input data showed the scaffolding framework is
able to compactly generalise for the parametrisation of arbitrarily complex
structures with minimal elements and low parametric density. A novel
interfacing configuration was introduced, i.e. the quadrilateral junction
simplex, modelling the spatial partitioning of generic junctions for an
underlying network of arbitrary connectivity patterns. Conforming and data-
driven fitting strategies were presented to determine the geometrical
embedding of the scaffolding with (semi-)automatic pipelines for polygonal
surface meshes and volumetric image-based data, respectively. A novel solid
smoothing paradigm is devised for the scaffolding lattice to obtain higher
than positional geometric and parametric continuity at the interfaces and at
the external boundaries of the structures, globally recovering an organic
shape in the form of a continuum. State-of-the-art tracing techniques employ
raster (un-)structured polyhedral meshes to define the innermost lattice of 3D
solid objects. Despite these intrinsically differ from parametric NURBS
geometries, both qualitative and quantitative evaluations are proposed for the
introduced framework based on analogies with the resulting scaffolding
elements, lattice density, computational time and accuracy of fitted
geometries.

To recover negligible geometrical deviations, tetrahedral meshes for
conventional finite-element analyses (Si, 2015) required a dense tessellation
of elements and vertices, which differs by several orders of magnitude to the
proposed parametrisation. Geometrical degradation is observed for iso-
parametric optimised meshes, in particular for organic synthetic shapes and
CAD designs with irregular and organic shapes.

The compactness of the introduced scaffolding was validated also on irregular
geometries form 3D tessellated surfaces, where the resulting density of
control points in the solid lattice was considerably lower than the input
surface vertices. By leveraging geodesic strategies and a conforming sampling,
the proposed scaffolding parametrisation reported high geometrical accuracy
with respect to the reference geometries, with minor deviations for sharp and
high-curvature features. Smoothness of the parametrised structures achieved
nearly G2 continuity globally.

Similarly, convoluted anatomical structures were compactly parametrised from
volumetric imaging data with a fully automatic pipeline. Conforming fitting
strategies based on geodesics and computer vision algorithms reported minor
geometric deviations comparable to the input data resolution. Smoothness and
continuity of the parametrised biological structures achieved an organic
profile underlying a homogeneous medium in the form of a continuum.

Several IGA computational simulations validated the compatibility of the
parametrised scaffolding with benchmarking performance. The simulations
scalability was evaluated on equivalent hexahedral meshes obtained with
cascades of conforming domain subdivisions to achieve similar element density
as in state-of-the-art techniques(Livesu et al., 2020). IGA simulations
successfully converged for the proposed scaffolding framework in all cases,
with smooth, continuous and consistent solution profiles. The memory
footprint, however, determined the limiting factor for simulations on
equivalent hexahedral tessellations, with a quadratic trend for increasing
density of the meshes.

Previous studies on hexahedral meshes (Livesu et al., 2020; Gao et al., 2017;
Takayama, 2019), sought a solid re-parametrisation from raster surface
geometries as a viable option towards an IGA-compatible domain. In similar CAD
applications, the input tessellated geometry employed in finite-element
analyses is often even obtained from rasterising an original vectorial design.
For the sake of finite-element analyses, such re-parametrisation results in a
convoluted and design-wise inefficient approach, specifically when IGA is
used. Although mesh degradation can be mitigated with denser lattices, the
migration from a parametric continuous domain to finite-resolution meshes
introduces geometrical approximations and downgrades higher degrees of
continuity to a pure positional one. Moreover, generic tetrahedral or
hexahedral-dominant meshes cannot be directly fed to IGA, since all mesh
elements require a conforming cuboid topology. Any other cell type in the
mesh, e.g. pyramidal or wedge-like elements, invalidate the computational
domain when a Galerkin discretisation (Cottrell et al., 2009) is considered
for IGA. Given the affinity of linear-element cuboid hexahedral meshes with
tensors of splines, a parametric equivalent can be determined. However,
fitting NURBS cuboids a-posteriori results in a dramatic increase of control
points density for higher-degree spline elements. By elevating a linear-degree
cuboid (minimal lattice of 8 control points) to an equivalent of cubic-degree
(minimal lattice of 64 control points) the lattice density increases by
700$\%$ for the individual cell. Such increase of lattice density does not
introduce, per se, higher degrees of continuity at the boundaries and
interfaces of the solid mesh, instead the latter would simply remain
positional (G0) in the absence of a solid smoothing strategy. From an IGA-
conforming perspective, applying the required parametric constraints to an
unstructured solid mesh may result in a complete dependency of all univariate
components, i.e. along $u$, $v$, and $w$ in all elements, with limited
parametrisation flexibility with respect to the introduced scaffolding. Also,
without an organised aggregation strategy, the scalability and simulation
performance advantage observed on compact representations is dramatically
diluted against dense tessellations, ultimately resulting in prohibitive or
infeasible finite-element computations.

Based on the introduced framework, we argue that a streamlined solid design
for finite-element analysis with IGA should rather begin with a consistent
construction of a conforming scaffolding. In the general form, this could
leverage similar graphic user-interfaces, modelling routines and virtual
environments employed for the construction of more complex composite CAD
geometries. This process would make use of guided and semi-automatic
interactions for industrial and engineering designs, as in commercially
available software, where modelling geometries rely on interactive user
experiences. Nonetheless, simpler geometries may benefit from a more
streamlined and fully-automatic scaffolding pipeline (e.g. anatomical
structures), which can dramatically impact on in-silico medicine and digital
twins formulations on large scale, addressing complex systems, diverse
physical scenarios and continuous solutions of PDEs.

##### Limitations and Future Work

Current limitations of the introduced scaffolding framework stand in
automatically mapping on reference surface meshes complex junction partitions
for particularly busy branching patterns. Although the existence of the
junction simplex is guaranteed (Suárez and Hubert, 2018), the overly compact
geometrical embedding may result counter-intuitive and difficult to visualise,
especially for the parametrisation of shape-varying structures as in complex
polygonal tessellations. As described in appendices A and B, alternative IGA-
compatible scaffolding configurations are feasible, and may locally
accommodate for geometrical simplifications and topological adjustments. A
simplified junction partition, e.g. the grafting scaffolding in appendix B,
can be developed for a non-minimal yet highly compact scaffolding, which
locally decomposes the busy junctions with extra cuts and interleaving
elements, diluting therefore the cascade of bisections in a more regular and
distributed fashion. Also, a more flexible formulation of the quadrilateral
junction simplexes may locally relax the introduced geometrical and
directional embedding to better model irregular and size-varying branching
patterns. On the one hand, this would improve the accuracy of the in-house
fitting paradigm for polygonal structures with a minimal distortion
parametrisation of the boundary surface, e.g. adopting standard UV-mapping and
seam-cuts unwrapping techniques (Lévy et al., 2002). On the other hand, a more
regular and distributed branching pattern of the junction simplex would reduce
the valence of extraordinary vertices, therefore improving the geometrical
smoothness and parametric continuity close to the singularities. Although the
introduced solid smoothing paradigm achieves nearly G2 geometrical smoothness
and parametric continuity for the entirety of the solid scaffolding, a minor
degradation is observed at the medial locus and in the neighbourhood of
extraordinary vertices. This is however in line with the considerations of
(Stam, 2001) for topological singularities. The continuity observed in our
parametrised geometries replicates the same experimental outcomes in (Stam,
2001; Bajaj et al., 2002) for vertices of varying valence. The devised
smoothing strategy builds on the same averaging scheme, here with ad-hoc
topological variations for multi-patch solid lattices. Based on experimental
observations, the introduced smoothing paradigm can be similarly applied also
to solid multi-patch meshes which exhibit a different scaffolding
configuration, comprising (un-)structured hexahedral cuboid elements. As shown
in appendix B, both Cross and Frame geometries, exhibit a face-wise
6-connected scaffolding, similarly to hexahedral meshes, where each branch is
extruded, without employing a quadratic junction simplex.

An organised structured and conformal aggregation strategy would potentially
convert a cuboid hexahedral mesh into a compact scaffolding by means of an
adaptive trade-off between splitting, refining and merging of cuboid elements.
This would represent a viable strategy for composite geometries requiring
extensive sub-block decomposition (Chuang et al., 2000; Lu et al., 2017; Wang
et al., 2017; Takayama, 2019), without further geometric re-parametrisation.
An adaptive refinement may also address individual vectorial elements; recent
IGA developments support hierarchical and T-NURBS elements (Hesch et al.,
2016; Garau and Vázquez, 2018; Bracco et al., 2018; He et al., 2006; Bazilevs
et al., 2010). With this view, a local lattice refinement may be selectively
employed in the neighbourhood of discontinuous regions and singular points
(i.e. extraordinary vertices) with the aim of improving the local geometrical
support, the smoothness and parametric continuity and towards adaptive degrees
of freedom for computational simulations. Adaptive and hierarchical local
subdivisions of the vectorial elements (Garau and Vázquez, 2018; Bracco et
al., 2018) could optimally circumvent the conforming propagation of the
lattice refinement, as in a continuous Galerking method (Cockburn et al.,
2009; Chen et al., 2016), without globally affecting the topological
formulation and the parametric arrangement of the entire scaffolding.
Prospectively, the extension to hierarchical NURBS scaffolding elements could
be further coupled with higher degrees splines and an optimal solid smoothing
strategy along the lines of (Cashman et al., 2009; Stam, 1998b) where
localised hierarchical refinements follow geometrical adjustments based on
non-uniform knot vectors and on the diagonalisation of the lattice subdivision
matrices.

##### Concluding Remarks

The introduced scaffolding framework constitutes a fundamental advance towards
bridging the gap between compactly tracing vectorial solid structures and
continuous-domain computational simulations. We believe our contributions
could further stimulate and inspire novel developments with revolutionising
applications on multi-disciplinary grounds for next-generation computational
simulations.

## References

  * (1)
  * Aspert et al. (2002) N. Aspert, D. Santa-Cruz, and T. Ebrahimi. 2002.  MESH: measuring errors between surfaces using the Hausdorff distance. In _Proceedings. IEEE International Conference on Multimedia and Expo_ , Vol. 1. 705–708 vol.1.  https://doi.org/10.1109/ICME.2002.1035879
  * Attene et al. (2016) Marco Attene, Daniela Cabiddu, Stefano Gagliardo, Franca Giannini, and Marina Monti. 2016\.  A web repository to describe and execute shape oriented workflows.  _Computer-Aided Design and Applications_ 13, 5 (2016), 637–646.
  * Au et al. (2008) Oscar Kin-Chung Au, Chiew-Lan Tai, Hung-Kuo Chu, Daniel Cohen-Or, and Tong-Yee Lee. 2008\.  Skeleton extraction by mesh contraction.  _ACM transactions on graphics (TOG)_ 27, 3 (2008), 1–10.  https://doi.org/10.1145/1399504.1360643
  * Auricchio et al. (2007) F Auricchio, L Beirao da Veiga, Annalisa Buffa, C Lovadina, A Reali, and G Sangalli. 2007\.  A fully “locking-free” isogeometric approach for plane linear elasticity problems: A stream function formulation.  _Computer methods in applied mechanics and engineering_ 197, 1-4 (2007), 160–172.  https://doi.org/10.1016/j.cma.2007.07.005
  * Bajaj et al. (2002) Chandrajit Bajaj, Scott Schaefer, Joe Warren, and Guoliang Xu. 2002.  A subdivision scheme for hexahedral meshes.  _The visual computer_ 18, 5-6 (2002), 343–356.  https://doi.org/10.1007/s003710100150
  * Balzer and Werling (2010) Jonathan Balzer and Stefan Werling. 2010.  Principles of shape from specular reflection.  _Measurement_ 43, 10 (2010), 1305–1317.  https://doi.org/10.1016/j.measurement.2010.07.013
  * Bazilevs et al. (2006) Yuri Bazilevs, L Beirao da Veiga, J Austin Cottrell, Thomas JR Hughes, and Giancarlo Sangalli. 2006\.  Isogeometric analysis: approximation, stability and error estimates for h-refined meshes.  _Mathematical Models and Methods in Applied Sciences_ 16, 07 (2006), 1031–1090.  https://doi.org/10.1142/s0218202506001455
  * Bazilevs et al. (2010) Yuri Bazilevs, Victor M Calo, John A Cottrell, John A Evans, Thomas Jr R Hughes, S Lipton, Michael A Scott, and Thomas W Sederberg. 2010.  Isogeometric analysis using T-splines.  _Computer Methods in Applied Mechanics and Engineering_ 199, 5-8 (2010), 229–263.  https://doi.org/10.1016/j.cma.2009.02.036
  * Bazilevs et al. (2008) Yuri Bazilevs, Victor M Calo, Thomas JR Hughes, and Yongjie Zhang. 2008.  Isogeometric fluid-structure interaction: theory, algorithms, and computations.  _Computational mechanics_ 43, 1 (2008), 3–37.  https://doi.org/10.1007/s00466-008-0315-x
  * Bingol and Krishnamurthy (2019) Onur Rauf Bingol and Adarsh Krishnamurthy. 2019.  NURBS-Python: An open-source object-oriented NURBS modeling framework in Python.  _SoftwareX_ 9 (2019), 85–94.  https://doi.org/10.1016/j.softx.2018.12.005
  * Bracco et al. (2018) Cesare Bracco, Carlotta Giannelli, and Rafael Vázquez. 2018\.  Refinement algorithms for adaptive isogeometric methods with hierarchical splines.  _axioms_ 7, 3 (2018), 43.  https://doi.org/10.3390/axioms7030043
  * Bucelli et al. (2021) Michele Bucelli, Matteo Salvador, Alfio Quarteroni, et al. 2021\.  Multipatch Isogeometric Analysis for electrophysiology: Simulation in a human heart.  _Computer Methods in Applied Mechanics and Engineering_ 376 (2021), 113666\.  https://doi.org/10.1016/j.cma.2021.113666
  * Bullitt et al. (2005) Elizabeth Bullitt, Donglin Zeng, Guido Gerig, Stephen Aylward, Sarang Joshi, J. Keith Smith, Weili Lin, and Matthew G. Ewend. 2005\.  Vessel Tortuosity and Brain Tumor Malignancy: A Blinded Study.  _Academic Radiology_ 12, 10 (2005), 1232–1240.  https://doi.org/10.1016/j.acra.2005.05.027
  * Burkhart et al. (2010) Daniel Burkhart, Bernd Hamann, and Georg Umlauf. 2010\.  Iso-geometric Finite Element Analysis Based on Catmull-Clark: ubdivision Solids. In _Computer Graphics Forum_ , Vol. 29. Wiley Online Library, 1575–1584.  https://doi.org/10.1111/j.1467-8659.2010.01766.x
  * Carraturo et al. (2019) Massimo Carraturo, Carlotta Giannelli, Alessandro Reali, and Rafael Vázquez. 2019.  Suitably graded THB-spline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes.  _Computer Methods in Applied Mechanics and Engineering_ 348 (2019), 660–679.  https://doi.org/10.1016/j.cma.2019.01.044
  * Cashman et al. (2009) Thomas J Cashman, Ursula H Augsdörfer, Neil A Dodgson, and Malcolm A Sabin. 2009.  NURBS with extraordinary points: high-degree, non-uniform, rational subdivision schemes.  In _ACM SIGGRAPH 2009 papers_. 1–9.  https://doi.org/10.1145/1576246.1531352
  * Catmull and Clark (1978) Edwin Catmull and James Clark. 1978.  Recursively generated B-spline surfaces on arbitrary topological meshes.  _Computer-aided design_ 10, 6 (1978), 350–355.  https://doi.org/10.1145/280811.280992
  * Chen et al. (2016) L Chen, B Simeon, and S Klinkel. 2016.  A NURBS based Galerkin approach for the analysis of solids in boundary representation.  _Computer Methods in Applied Mechanics and Engineering_ 305 (2016), 777–805.  https://doi.org/10.1016/j.cma.2016.03.019
  * Cheng et al. (2015) Yuanzhi Cheng, Xin Hu, Ji Wang, Yadong Wang, and Shinichi Tamura. 2015.  Accurate vessel segmentation with constrained B-snake.  _IEEE Transactions on Image Processing_ 24, 8 (2015), 2440–2455.  https://doi.org/10.1109/tip.2015.2417683
  * Chuang et al. (2000) Jen-Hui Chuang, Chi-Hao Tsai, and Min-Chi Ko. 2000.  Skeletonisation of three-dimensional object using generalized potential field.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 22, 11 (2000), 1241–1251.  https://doi.org/10.1109/34.888709
  * Cockburn et al. (2009) Bernardo Cockburn, Jayadeep Gopalakrishnan, and Raytcho Lazarov. 2009.  Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems.  _SIAM J. Numer. Anal._ 47, 2 (2009), 1319–1365.  https://doi.org/10.1137/070706616
  * Cottrell et al. (2009) J Austin Cottrell, Thomas JR Hughes, and Yuri Bazilevs. 2009\.  _Isogeometric analysis: toward integration of CAD and FEA_.  John Wiley & Sons.
  * Crane et al. (2020) Keenan Crane, Marco Livesu, Enrico Puppo, and Yipeng Qin. 2020\.  A Survey of Algorithms for Geodesic Paths and Distances.  _arXiv preprint arXiv:2007.10430_ (2020).
  * Dalcin et al. (2016) L. Dalcin, N. Collier, P. Vignal, A.M.A. Côrtes, and V.M. Calo. 2016.  PetIGA: A framework for high-performance isogeometric analysis.  _Computer Methods in Applied Mechanics and Engineering_ 308 (2016), 151–181.  https://doi.org/10.1016/j.cma.2016.05.011
  * De Falco et al. (2011) Carlo De Falco, Alessandro Reali, and R Vázquez. 2011\.  GeoPDEs: a research tool for isogeometric analysis of PDEs.  _Advances in Engineering Software_ 42, 12 (2011), 1020–1034.  https://doi.org/10.1016/j.advengsoft.2011.06.010
  * Doo and Sabin (1978) Daniel Doo and Malcolm Sabin. 1978.  Behaviour of recursive division surfaces near extraordinary points.  _Computer-Aided Design_ 10, 6 (1978), 356–360.  https://doi.org/10.1145/280811.280991
  * Farin and Hansford (1999) Gerald Farin and Dianne Hansford. 1999.  Discrete coons patches.  _Computer aided geometric design_ 16, 7 (1999), 691–700.  https://doi.org/10.1016/s0167-8396(99)00031-x
  * Gao et al. (2017) Xifeng Gao, Wenzel Jakob, Marco Tarini, and Daniele Panozzo. 2017\.  Robust hex-dominant mesh generation using field-guided polyhedral agglomeration.  _ACM Transactions on Graphics (TOG)_ 36, 4 (2017), 1–13.  https://doi.org/10.1145/3072959.3073676
  * Gao et al. (2019) Xifeng Gao, Hanxiao Shen, and Daniele Panozzo. 2019.  Feature Preserving Octree-Based Hexahedral Meshing. In _Computer Graphics Forum_ , Vol. 38. Wiley Online Library, 135–149.  https://doi.org/10.1111/cgf.13795
  * Garau and Vázquez (2018) Eduardo M Garau and Rafael Vázquez. 2018.  Algorithms for the implementation of adaptive isogeometric methods using hierarchical B-splines.  _Applied Numerical Mathematics_ 123 (2018), 58–87.  https://doi.org/10.1016/j.apnum.2017.08.006
  * Gomez et al. (2010) Hector Gomez, Thomas JR Hughes, Xesús Nogueira, and Victor M Calo. 2010.  Isogeometric analysis of the isothermal Navier–Stokes–Korteweg equations.  _Computer Methods in Applied Mechanics and Engineering_ 199, 25-28 (2010), 1828–1840.  https://doi.org/10.1016/j.cma.2010.02.010
  * Gregson et al. (2011) James Gregson, Alla Sheffer, and Eugene Zhang. 2011\.  All-hex mesh generation via volumetric polycube deformation. In _Computer graphics forum_ , Vol. 30. Wiley Online Library, 1407–1416.  https://doi.org/10.1111/j.1467-8659.2011.02015.x
  * Hanson (1994) Andrew J Hanson. 1994\.  Quaternion Frenet frames: Making optimal tubes and ribbons from curves.  _Computer Science Department, Indiana University Bloomington, In_ 47405 (1994).
  * He et al. (2006) Ying He, Kexiang Wang, Hongyu Wang, Xianfeng Gu, and Hong Qin. 2006.  Manifold T-spline. In _International conference on geometric modeling and processing_. Springer, 409–422.  https://doi.org/10.1007/11802914_29
  * Hesch et al. (2016) C Hesch, M Franke, M Dittmann, and I Temizer. 2016\.  Hierarchical NURBS and a higher-order phase-field approach to fracture for finite-deformation contact problems.  _Computer Methods in Applied Mechanics and Engineering_ 301 (2016), 242–258.  https://doi.org/10.1016/j.cma.2015.12.011
  * Ho-Le (1988) K Ho-Le. 1988.  Finite element mesh generation methods: a review and classification.  _Computer-aided design_ 20, 1 (1988), 27–38.  https://doi.org/10.1016/0010-4485(88)90138-8
  * Hughes et al. (2005) Thomas JR Hughes, John A Cottrell, and Yuri Bazilevs. 2005\.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  _Computer methods in applied mechanics and engineering_ 194, 39-41 (2005), 4135–4195.  https://doi.org/10.1016/j.cma.2004.10.008
  * Jüttler et al. (2014) Bert Jüttler, Ulrich Langer, Angelos Mantzaflaris, Stephen E Moore, and Walter Zulehner. 2014\.  Geometry+ simulation modules: Implementing isogeometric analysis.  _PAMM_ 14, 1 (2014), 961–962.  https://doi.org/10.1002/pamm.201410461
  * Kanitsar et al. (2002) A. Kanitsar, D. Fleischmann, R. Wegenkittl, P. Felkel, and E. Groller. 2002.  CPR - curved planar reformation. In _IEEE Visualization, 2002. VIS 2002._ 37–44.  https://doi.org/10.1109/VISUAL.2002.1183754
  * Kimmel and Sethian (1998) Ron Kimmel and James A Sethian. 1998.  Computing geodesic paths on manifolds.  _Proceedings of the national academy of Sciences_ 95, 15 (1998), 8431–8435.  https://doi.org/10.1073/pnas.95.15.8431
  * Krishnamurthy and Levoy (1996) Venkat Krishnamurthy and Marc Levoy. 1996.  Fitting smooth surfaces to dense polygon meshes. In _Proceedings of the 23rd annual conference on Computer graphics and interactive techniques_. 313–324.  https://doi.org/10.1145/237170.237270
  * Lévy et al. (2002) Bruno Lévy, Sylvain Petitjean, Nicolas Ray, and Jérome Maillot. 2002.  Least squares conformal maps for automatic texture atlas generation.  _ACM transactions on graphics (TOG)_ 21, 3 (2002), 362–371.  https://doi.org/10.1145/566654.566590
  * Livesu et al. (2016) Marco Livesu, Alessandro Muntoni, Enrico Puppo, and Riccardo Scateni. 2016.  Skeleton-driven adaptive hexahedral meshing of tubular shapes. In _Computer Graphics Forum_ , Vol. 35. Wiley Online Library, 237–246.  https://doi.org/10.1111/cgf.13021
  * Livesu et al. (2020) Marco Livesu, Nico Pietroni, Enrico Puppo, Alla Sheffer, and Paolo Cignoni. 2020.  LoopyCuts: Practical Feature-Preserving Block Decomposition for Strongly Hex-Dominant Meshing.  _ACM Transactions on Graphics (SIGGRAPH)_ 39, 4 (2020).  https://doi.org/10.1145/3386569.3392472
  * Livesu et al. (2015) Marco Livesu, Alla Sheffer, Nicholas Vining, and Marco Tarini. 2015.  Practical hex-mesh optimization via edge-cone rectification.  _ACM Transactions on Graphics (TOG)_ 34, 4 (2015), 1–11.  https://doi.org/10.1145/2766905
  * Livesu et al. (2013) Marco Livesu, Nicholas Vining, Alla Sheffer, James Gregson, and Riccardo Scateni. 2013.  Polycut: Monotone graph-cuts for polycube base-complex construction.  _ACM Transactions on Graphics (TOG)_ 32, 6 (2013), 1–12.  https://doi.org/10.1145/2508363.2508388
  * Loop (1987) Charles Loop. 1987\.  Smooth subdivision surfaces based on triangles.  _Master’s thesis, University of Utah, Department of Mathematics_ (1987).
  * Lu et al. (2017) Jean Hsiang-Chun Lu, William Roshan Quadros, and Kenji Shimada. 2017.  Evaluation of user-guided semi-automatic decomposition tool for hexahedral mesh generation.  _Journal of Computational Design and Engineering_ 4, 4 (2017), 330–338.  https://doi.org/10.1016/j.jcde.2017.05.001
  * Ma and Kruth (1998) Weiyin Ma and J-P Kruth. 1998.  NURBS curve and surface fitting for reverse engineering.  _The International Journal of Advanced Manufacturing Technology_ 14, 12 (1998), 918–927.  https://doi.org/10.1007/bf01179082
  * Martinez et al. (2004) Dimas Martinez, Luiz Velho, and P Cezar Carvalho. 2004\.  Geodesic paths on triangular meshes. In _Proceedings. 17th Brazilian Symposium on Computer Graphics and Image Processing_. IEEE, 210–217.
  * Mortenson (1997) Michael E Mortenson. 1997\.  _Geometric modeling_.  John Wiley & Sons, Inc.
  * Pagani and Scott (2018) Luca Pagani and Paul J Scott. 2018.  Curvature based sampling of curves and surfaces.  _Computer Aided Geometric Design_ 59 (2018), 32–48.  https://doi.org/10.1016/j.cagd.2017.11.004
  * Pan et al. (2016) Qing Pan, Guoliang Xu, Gang Xu, and Yongjie Zhang. 2016\.  Isogeometric analysis based on extended Catmull–Clark subdivision.  _Computers & Mathematics with Applications_ 71, 1 (2016), 105–119.  https://doi.org/10.1016/j.camwa.2015.11.012
  * Panotopoulou et al. (2018) Athina Panotopoulou, Elissa Ross, Kathrin Welker, Evelyne Hubert, and Géraldine Morin. 2018.  Scaffolding a skeleton.  In _Research in Shape Analysis_. Springer, 17–35.  https://doi.org/10.1007/978-3-319-77066-6_2
  * Panozzo et al. (2014) Daniele Panozzo, Enrico Puppo, Marco Tarini, and Olga Sorkine-Hornung. 2014.  Frame fields: Anisotropic and non-orthogonal cross fields.  _ACM Transactions on Graphics (TOG)_ 33, 4 (2014), 1–11.  https://doi.org/10.1145/2601097.2601179
  * Petera and Pittman (1994) J. Petera and J. F. T. Pittman. 1994.  Isoparametric Hermite Elements.  _Internat. J. Numer. Methods Engrg._ 37, 20 (1994), 3489–3519.  https://doi.org/10.1002/nme.1620372006
  * Piccinelli et al. (2009) Marina Piccinelli, Alessandro Veneziani, David A Steinman, Andrea Remuzzi, and Luca Antiga. 2009\.  A framework for geometric analysis of vascular structures: application to cerebral aneurysms.  _IEEE transactions on medical imaging_ 28, 8 (2009), 1141–1155.  https://doi.org/10.1109/tmi.2009.2021652
  * Piegl (1991) L. Piegl. 1991.  On NURBS: a Survey.  _IEEE Computer Graphics and Applications_ 11, 01 (jan 1991), 55–71.  https://doi.org/10.1109/38.67702
  * Piegl and Tiller (1996) Les Piegl and Wayne Tiller. 1996.  _The NURBS book_.  Springer Science & Business Media.  https://doi.org/10.1007/978-3-642-59223-2
  * Pietroni et al. (2016) Nico Pietroni, Enrico Puppo, Giorgio Marcias, Roberto Scopigno, and Paolo Cignoni. 2016\.  Tracing field-coherent quad layouts. In _Computer Graphics Forum_ , Vol. 35. Wiley Online Library, 485–496.  https://doi.org/10.1111/cgf.13045
  * Rao (2017) Singiresu S Rao. 2017\.  _The finite element method in engineering_.  Butterworth-heinemann.
  * Russo (2006) Mario Russo. 2006\.  _Polygonal modeling: basic and advanced techniques_.  Jones & Bartlett Learning.
  * Saha et al. (2016) Punam K Saha, Gunilla Borgefors, and Gabriella Sanniti di Baja. 2016\.  A survey on skeletonization algorithms and their applications.  _Pattern recognition letters_ 76 (2016), 3–12.  https://doi.org/10.1016/j.patrec.2015.04.006
  * Saha et al. (2017) Punam K Saha, Gunilla Borgefors, and Gabriella Sanniti di Baja. 2017\.  _Skeletonization: Theory, methods and applications_.  Academic Press.
  * Sethian (1996) James A Sethian. 1996\.  A fast marching level set method for monotonically advancing fronts.  _Proceedings of the National Academy of Sciences_ 93, 4 (1996), 1591–1595.  https://doi.org/10.1073/pnas.93.4.1591
  * Sethian (1999) James Albert Sethian. 1999\.  _Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science_. Vol. 3.  Cambridge university press.
  * Si (2015) Hang Si. 2015.  TetGen, a Delaunay-based quality tetrahedral mesh generator.  _ACM Transactions on Mathematical Software (TOMS)_ 41, 2 (2015), 1–36.  https://doi.org/10.1145/2629697
  * Stam (1998a) Jos Stam. 1998a.  Evaluation of loop subdivision surfaces. In _SIGGRAPH’98 CDROM Proceedings_. Citeseer.
  * Stam (1998b) Jos Stam. 1998b.  Exact evaluation of Catmull-Clark subdivision surfaces at arbitrary parameter values. In _Proceedings of the 25th annual conference on Computer graphics and interactive techniques_. 395–404.  https://doi.org/10.1145/280814.280945
  * Stam (2001) Jos Stam. 2001.  On subdivision schemes generalizing uniform B-spline surfaces of arbitrary degree.  _Computer Aided Geometric Design_ 18, 5 (2001), 383–396.  https://doi.org/10.1016/s0167-8396(01)00038-3
  * Stolpner et al. (2011) Svetlana Stolpner, Paul Kry, and Kaleem Siddiqi. 2011\.  Medial spheres for shape approximation.  _IEEE transactions on pattern analysis and machine intelligence_ 34, 6 (2011), 1234–1240.  https://doi.org/10.1109/tpami.2011.254
  * Suárez and Hubert (2017) AJ Fuentes Suárez and Evelyne Hubert. 2017.  Scaffolding skeletons using spherical Voronoi diagrams.  _Electronic Notes in Discrete Mathematics_ 62 (2017), 45–50.  https://doi.org/10.1016/j.endm.2017.10.009
  * Suárez and Hubert (2018) AJ Fuentes Suárez and Evelyne Hubert. 2018.  Scaffolding skeletons using spherical Voronoi diagrams: Feasibility, regularity and symmetry.  _Computer-Aided Design_ 102 (2018), 83–93.  https://doi.org/10.1016/j.cad.2018.04.016
  * Surazhsky et al. (2005) Vitaly Surazhsky, Tatiana Surazhsky, Danil Kirsanov, Steven J Gortler, and Hugues Hoppe. 2005\.  Fast exact and approximate geodesics on meshes.  _ACM transactions on graphics (TOG)_ 24, 3 (2005), 553–560.  https://doi.org/10.1145/1186822.1073228
  * Takayama (2019) Kenshi Takayama. 2019\.  Dual sheet meshing: An interactive approach to robust hexahedralization. In _Computer Graphics Forum_ , Vol. 38. Wiley Online Library, 37–48.  https://doi.org/10.1111/cgf.13617
  * Tarini et al. (2004) Marco Tarini, Kai Hormann, Paolo Cignoni, and Claudio Montani. 2004.  Polycube-maps.  _ACM transactions on graphics (TOG)_ 23, 3 (2004), 853–860.  https://doi.org/10.1145/1186562.1015810
  * Urick et al. (2019) Benjamin Urick, Travis M Sanders, Shaolie S Hossain, Yongjie J Zhang, and Thomas JR Hughes. 2019\.  Review of patient-specific vascular modeling: template-based isogeometric framework and the case for CAD.  _Archives of Computational Methods in Engineering_ 26, 2 (2019), 381–404.  https://doi.org/10.1007/s11831-017-9246-z
  * Vázquez (2016) Rafael Vázquez. 2016\.  A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0.  _Computers & Mathematics with Applications_ 72, 3 (2016), 523–554.  https://doi.org/10.1016/j.camwa.2016.05.010
  * Vázquez and Buffa (2010) Rafael Vázquez and Annalisa Buffa. 2010.  Isogeometric analysis for electromagnetic problems.  _IEEE Transactions on Magnetics_ 46, 8 (2010), 3305–3308.  https://doi.org/10.1109/tmag.2010.2044563
  * Wang et al. (2008) Hongyu Wang, Ying He, Xin Li, Xianfeng Gu, and Hong Qin. 2008.  Polycube splines.  _Computer-Aided Design_ 40, 6 (2008), 721–733.  https://doi.org/10.1016/j.cad.2008.01.012
  * Wang et al. (2017) Rui Wang, Chun Shen, Jinming Chen, Haiyan Wu, and Shuming Gao. 2017.  Sheet operation based block decomposition of solid models for hex meshing.  _Computer-Aided Design_ 85 (2017), 123–137.  https://doi.org/10.1016/j.cad.2016.07.016
  * Xu et al. (2014) Gang Xu, Bernard Mourrain, André Galligo, and Timon Rabczuk. 2014.  High-quality construction of analysis-suitable trivariate NURBS solids by reparameterization methods.  _Computational Mechanics_ 54, 5 (2014), 1303–1313.  https://doi.org/10.1007/s00466-014-1060-y
  * Zhang et al. (2007) Yongjie Zhang, Yuri Bazilevs, Samrat Goswami, Chandrajit L Bajaj, and Thomas JR Hughes. 2007\.  Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow.  _Computer methods in applied mechanics and engineering_ 196, 29-30 (2007), 2943–2959.  https://doi.org/10.1016/j.cma.2007.02.009
  * Zhou and Jacobson (2016) Qingnan Zhou and Alec Jacobson. 2016.  Thingi10k: A dataset of 10,000 3d-printing models.  _arXiv preprint arXiv:1605.04797_ (2016).

##  Appendix A Wall Scaffolding Construction

\begin{overpic}[width=429.28616pt]{fig20R.png} \put(1.0,17.0){\small{(a)}}
\put(19.0,17.0){\small{(b)}} \put(40.0,17.0){\small{(c)}}
\put(50.0,17.0){\small{(d)}} \put(69.0,17.0){\small{(e)}} \end{overpic} Figure
21. Wall Scaffolding Construction: (a) 4-way-junction: 8 elements per
elongated branch forming an octagonal ring; (b) Additional capping elements
per terminal branch. (c) Boundaries and conforming interfaces of adjacent
elements facing the parametric directions; (d) Smooth limit geometry; (e) EggW
thin wall as closed hollow chamber: scaffolding comprising two sets of capping
elements.

The wall scaffolding is modelled as an extension of the luminal one (section
3.2). It comprises the organised set of vectorial elements occupying the
hollow volumetric space enclosing a luminal region, e.g. a set of connecting
hollow pipes or hollow chambers with thin or thick walls, in the form of a
solid shell (fig. 21). The wall scaffolding can be considered either
separately, or additionally coupled with a luminal scaffolding.

##### Minimal Elements

The minimal number of elements depends on the complexity of the underlying
structure. Building on concepts introduced in section 3.2, the wall
scaffolding requires 8 elements longitudinally surrounding each branch.
Additional 4 elements are required when one end of any branch is closed, i.e.
the wall covers and closes the terminal portion of the hollow structure. In
general, the elements of the wall scaffolding are proportional to the boundary
sides of a virtually nested luminal scaffolding, and the construction of the
wall resembles a conforming extrusion along the luminal interfaces.

##### Arrangement and Adjacency

In general the wall scaffolding cross-sectionally configures the elements of
each elongated branch in an octagonal ring (fig. 21). As in the luminal
region, these elements are elongated along the parametric direction $w$,
mapping the longitudinal direction of a branch. The parametric directions $u$
and $v$ of the wall scaffolding independently map a pseudo-circumferential
direction, and a pseudo-radial direction, respectively, relative to the
elongation of the branch. Each adjacent element of the wall scaffolding has
interfaces facing the pseudo-circumferential direction $u$. In the absence of
a coupled luminal scaffolding, two sets of boundaries are determined, both
facing the pseudo-radial direction $v$; one facing the interior hollow luminal
region, another facing the exterior of the wall structure. However, when the
wall is coupled with a luminal scaffolding, both scaffoldings share a mutual
set of interfaces, constituting a complementary and multi-compartmental nested
structure.

For a closed terminal branch the wall scaffolding additionally introduces 4
capping elements (fig. 21). These form a quadrant-like configuration and their
parametric arrangement is consistent with a luminal scaffolding element, i.e.
$u$ and $v$ no longer map independently the pseudo-circumferential and the
pseudo-radial directions, respectively. Each capping element determines at
least 4 interfaces. Two contiguous sides (one facing $u$, and the other facing
$v$) are shared with 2 adjacent capping elements; whereas other two contiguous
sides (one facing $u$, and the other facing $v$, on opposite directions) are
shared with 2 adjacent elongated elements.

Without a coupled luminal region, both sides of the capping element facing $w$
are referred as boundaries, each belonging to the interior and to the exterior
wall boundaries, respectively. Conversely, with a coupled luminal scaffolding,
each capping element shares an extra interface between the nested
scaffoldings.

The formulation of a wall scaffolding configuration at the interface of a
$n$-way-junction is a straightforward extension of the luminal quadrilateral
junction simplex $\mathcal{Q}$ in section 3.2. Both interior
($\mathcal{Q}_{\textrm{int}}$) and exterior ($\mathcal{Q}_{\textrm{ext}}$)
simplexes of the of the wall scaffolding replicate the same configuration and
differ in size, accounting for the pseudo-radial thickness of the wall. In
this case, the resulting set of interfaces facing $w$ is determined by
connecting one-to-one each vertex of $\mathcal{Q}_{\textrm{int}}$ with the
respective vertex of $\mathcal{Q}_{\textrm{ext}}$, and the cross-sectional
octagonal configuration is kept for all the branches.

##### Conforming Constraints

The conforming conditions depend on coupling the wall with a nested luminal
scaffolding.

For an uncoupled open wall scaffolding, no mutual dependency is determined for
the parametric components of each element along the different parametric
directions. Conversely, for an uncoupled closed wall scaffolding, i.e. in
presence of capping elements, the parametric components of each element along
$u_{(\textrm{wall~{}branch})}$ and $v_{(\textrm{wall~{}branch})}$ are
equivalent, whereas only the parametric component along
$w_{(\textrm{wall~{}branch})}$ is independent. The capping elements of each
terminal branch, however, make exception: all the parametric components are
dependent and mutually equivalent along $(u,v,w)_{(\textrm{wall~{}cap})}$ to
$(u,v)_{(\textrm{wall~{}branch})}$.

When a wall scaffolding is coupled with a nested luminal one, the conforming
constraints are mutually related between scaffoldings. In particular, for a
coupled closed wall scaffolding, the parametric components of each element
along $(u,v)_{(\textrm{wall~{}branch})}$ coincide with those of
$(u,v)_{(\textrm{luminal})}$. Also, the parametric components along
$w_{(\textrm{wall~{}brach})}$ coincide with those along
$w_{(\textrm{luminal})}$. All the capping elements of the closed wall
scaffolding are mutually dependent, i.e. $(u,v,w)_{(\textrm{wall~{}cap})}$
coincide to $(u,v)_{(\textrm{wall~{}branch})}$ as well as to
$(u,v)_{(\textrm{luminal})}$. Lastly, for a coupled open wall scaffolding, the
parametric components of each element along $u_{(\textrm{wall~{}branch})}$
coincide with those of $(u,v)_{(\textrm{luminal})}$, whereas
$v_{(\textrm{wall~{}branch})}$ and $w_{(\textrm{wall~{}branch})}$ are
independent.

##  Appendix B Other Conforming Scaffolding Configurations

\begin{overpic}[width=429.28616pt]{fig21R.png} \put(-1.0,36.0){\small{(a)}}
\put(33.0,36.0){\small{(b)}} \put(67.0,36.0){\small{(c)}}
\put(0.0,17.0){\small{(d)}} \put(31.0,17.0){\small{(e)}}
\put(66.0,17.0){\small{(f)}} \end{overpic} Figure 22. Other Scaffolding
Construction: (a) Cross and (b) Frame geometries with 6-face-connectivity
scaffolding. (c) Exceptionally reduced planar scaffolding for bifurcating
structures; (d) Associated configuration of the wall scaffolding. (e) Grafting
scaffolding: Pinocchio geometry with primary and secondary branches; (f)
Pinocchio topological variation, arrangement and adjacency.

For completeness, alternative conforming scaffolding configurations are simply
mentioned below. These may result in local deviations from the introduced
framework with minor topological and parametric variations affecting the
density of vectorial elements, their arrangement and the associated conforming
constraints.

By way of example, a 6-face-connectivity scaffolding is shown in fig. 22, for
the synthetic Cross and Frame geometries. Starting from an innermost cube,
such configuration consist of an independent conformal extrusion along each
boundary face into a set of adjacent vectorial elements. Note that the
extension of such configuration includes structured (e.g. voxel-based) or
unstructured hexahedral meshing approaches.

Based on section 3.2, a planar scaffolding exceptionally halves the number of
vectorial elements per branch. This is suitable for isolated 3-way-junctions
as in fig. 22. Alternatively, in case of a network, this configuration is
applicable only if all the junctions are either bifurcations, or they exhibit
a planar arrangement. In such case, note that the structural torsional twist
along connecting branches of the luminal scaffolding may increase up to
$\theta_{\text{max}}=\frac{\pi}{2}$ by adopting the same geometrical embedding
in section 4. Note that such exceptional reduction can affect the
configuration of a coupled wall scaffolding, or symmetrical and modular
structures (e.g. Gear).

Lastly, a grafting scaffolding defines a hierarchical branching relationship
and labels branches of the skeleton as primary, secondary, and so on (e.g.
based on size), these being interleaved by cross-sectional cuts in fig. 22.
This configuration is suitable for particularly busy or spatially close
junctions by locally simplifying the interfacing simplexes. Note that, the
grafting configuration reduces the scaffolding elements of the secondary
branch relative to the primary one, with potential limitations for downstream
branching and connectivity patterns.