# Topology-preserving Scan-based Immersed Isogeometric Analysis

Sai C Divi<EMAIL_ADDRESS>Clemens V Verhoosel<EMAIL_ADDRESS>Ferdinando Auricchio<EMAIL_ADDRESS>Alessandro Reali
<EMAIL_ADDRESS>E Harald van Brummelen<EMAIL_ADDRESS>Department of Mechanical Engineering, Eindhoven University of Technology,
5600MB Eindhoven, The Netherlands  Department of Civil Engineering and
Architecture, University of Pavia, 27100 Pavia, Italy

###### Abstract

To exploit the advantageous properties of isogeometric analysis (IGA) in a
scan-based setting, it is important to extract a smooth geometric domain from
the scan data (_e.g._ , voxel data). IGA-suitable domains can be constructed
by convoluting the grayscale data using B-splines. A negative side-effect of
this convolution technique is, however, that it can induce topological changes
in the process of smoothing when features with a size similar to that of the
voxels are encountered. This manuscript presents an enhanced B-spline-based
segmentation procedure using a refinement strategy based on truncated
hierarchical (TH)B-splines. A Fourier analysis is presented to explain the
effectiveness of local grayscale function refinement in repairing topological
anomalies. A moving-window-based topological anomaly detection algorithm is
proposed to identify regions in which the grayscale function refinements must
be performed. The criterion to identify topological anomalies is based on the
Euler characteristic, giving it the capability to distinguish between
topological and shape changes. The proposed topology-preserving THB-spline
image segmentation strategy is studied using a range of test cases. These
tests pertain to both the segmentation procedure itself, and its application
in an immersed IGA setting.

###### Abstract

To exploit the advantageous properties of isogeometric analysis (IGA) in a
scan-based setting, it is important to extract a smooth geometric domain from
the scan data (_e.g._ , voxel data). IGA-suitable domains can be constructed
by convoluting the grayscale data using B-splines. A negative side-effect of
this convolution technique is, however, that it can induce topological changes
in the process of smoothing when features with a size similar to that of the
voxels are encountered. This manuscript presents an enhanced B-spline-based
segmentation procedure using a refinement strategy based on truncated
hierarchical (TH)B-splines. A Fourier analysis is presented to explain the
effectiveness of local grayscale function refinement in repairing topological
anomalies. A moving-window-based topological anomaly detection algorithm is
proposed to identify regions in which the grayscale function refinements must
be performed. The criterion to identify topological anomalies is based on the
Euler characteristic, giving it the capability to distinguish between
topological and shape changes. The proposed topology-preserving THB-spline
image segmentation strategy is studied using a range of test cases. These
tests pertain to both the segmentation procedure itself, and its application
in an immersed IGA setting.

###### keywords:

Scan-based analysis, Isogeometric analysis (IGA), Finite Cell Method (FCM),
Image filtering, Topology detection

††journal: CMAME

##  1 Introduction

Computational analyses based on volumetric scan data are of interest in many
fields of research, such as biomechanics, geomechanics, material science,
microstructural analysis, and many more. Scan-based simulations are inherently
three-dimensional and, frequently, the computational domains are complex, both
in terms of geometry and in terms of topology. In addition, the data sets
obtained from, _e.g._ , tomography or photogrammetry techniques are large in
size and represented in data formats which are not directly suitable for
analysis (_e.g._ , DICOM111Digital Imaging and Communications in Medicine,
NIfTI222Neuroimaging Informatics Technology Initiative). For these reasons,
performing high-fidelity simulations at practical computational costs is still
very challenging.

Over the past decade, the use of Isogeometric Analysis (IGA) [1] to perform
accurate scan-based simulations at acceptable computational costs has been
explored. Originally, IGA was proposed as an analysis paradigm to better
integrate analysis and design by employing the spline basis functions from
Computer Aided Design (CAD) directly for the analysis, without intermediate
geometry clean-up and meshing operations [2], as required in traditional
finite element analyses. The advantageous properties of spline basis
functions, in particular their higher-order regularity, have made IGA also
attractive for simulations where the geometry is not represented by a CAD
model [3, 4, 5, 6]. In the context of scan-based isogeometric analysis, the
usage of standard discretization techniques such as the voxel method [7] or
(unstructured) conforming mesh approaches (_e.g._ , generated using a marching
cube algorithm [8, 9]) is not possible, as this deteriorates the favorable
properties of IGA. Therefore, various enhanced isogeometric techniques for
scan-based analysis have been developed. These can be summarized as follows:

  * 1.

Template-fitting techniques construct a spline-based template geometry
(typically consisting of multiple patches) that captures the essential
features of the scanned object, and subsequently apply fitting procedures to
reposition the control points of the template to match the scan data. An
advantage of such techniques is that an explicit parametrization of the scan
object is retrieved, which is beneficial from an analysis point of view. The
downside of template fitting techniques is that they are not useful in cases
where the topology of the object is not known _a priori_. For cases where the
topology of the object is known but complex, template fitting techniques are
typically laborsome. Template fitting techniques using IGA have, for example,
been applied to fluid-structure interaction analyses of arterial blood flow
[10, 11], patient-specific blood-flow analyses in arteries [3, 12] (see [13]
for a review), the mechanical behavior of a femur (consisting of both hard
outer cortical bone and inner trabecular bone) [14], the mechanical behavior
of a carotid artery stent [15], damage in composite laminates [16], the human
heart [17], and fluid-structure interactions of a heart valve [18].

  * 2.

Immersed methods construct a structured mesh for the entire scan domain,
typically a box, and then represent the scanned object by an inside/outside
relation to indicate whether a specific point in the box is part of the
object. In the context of scan-based analysis, the inside/outside relation can
directly be retrieved from the scan data (intensity) in the form of a level
set function [19], but the immersed concept can also be applied to obtain
volumetric representations of objects identified by boundary representations
(_e.g._ , BREP [20, 21, 22, 23, 24, 25], STL333Standard Triangle Language or
Standard Tessellation Language [26, 27, 28, 29, 30]) and to represent trimming
operations in CAD models (_e.g._ , [31, 32, 33]). The advantage of immersed
techniques is their versatility with respect to the geometry and topology of
the scanned objects, in the sense that the analysis procedure is not
essentially affected by increasing the complexity of the objects. From an
analysis perspective, the immersed approach poses additional challenges
compared to mesh-fitting analyses. Most notably, numerical integration of
trimmed elements is challenging from an efficiency point of view, application
of boundary conditions can be non-standard, and the system of equations is
generally ill-conditioned without dedicated treatment. Immersed scan-based IGA
has, for example, been applied for the analysis of trabecular bone [22, 19,
34], coated metal foams [35], porous media [36], metal castings [37] and in
additive manufacturing [38].

The choice for either a template-fitting technique or an immersed technique is
to a large extent dictated by the topological complexity of the scanned
object. If a template with a reasonably low number of control points can be
constructed for the object of interest, template-fitting is generally
favorable. If the creation of a template is impractical, immersed methods are
preferred. It should be noted that it can be favorable to combine the two
techniques for particular scan-based analyses, to exploit the advantages of
both of them. A noteworthy example in this regard is the heart-valve problem
by Kamensky _et al._ [39, 40, 41], where the heart chambers are modeled
through a template-fitting approach, and the moving valves through an immersed
approach.

In this work we build on the immersed scan-based analysis framework originally
proposed in Ref. [19]. In recent years, significant progress has been made to
tackle the above-mentioned challenges associated with immersed methods. A
myriad of advanced integration techniques has been developed to reduce the
computational burden associated with the integration of trimmed elements; see,
_e.g._ , Refs. [42, 43, 44, 45]. Nitsche’s method [46] has been demonstrated
to be a reliable technique to impose essential boundary conditions along
immersed boundaries, _e.g._ , [47, 48, 23, 49], and various techniques have
been developed to construct a parametrization of immersed boundaries to impose
boundary conditions [50, 51, 49, 52]. With respect to the stability and
conditioning, fundamental understanding was obtained in Refs. [23, 53], and
various types of remedies have been demonstrated to be effective, see, _e.g._
, [54, 50, 51, 55, 56, 57, 5, 36, 58, 37]. The various advancements of
immersed techniques have made it into a versatile technique for scan-based
analyses, being able to model problems in different physical domains (solid
mechanics, fluid dynamics), considering high-performance computing problems
[26, 59, 42, 19, 60, 36, 34, 37], and being applicable in the context of
highly nonlinear problems [61, 62, 63, 64].

Although immersed isogeometric analysis is very flexible with respect to
capturing topologically complex volumetric objects, topological anomalies
associated with the image segmentation procedure can occur [65, 66, 67, 68,
69]. This is specifically the case when the resolution of the scan data is
only nearly sufficient to capture the smallest features in the scanned object.
Image smoothing operations, in the case of the procedure of Ref. [19]
associated with the order of the spline level set construction, can trigger
undesirable topology alterations (_e.g._ , closure of a channel,
(dis)connection of structures). The occurrence of topological anomalies
associated with smoothing operations is well understood in the field of image
segmentation (_e.g._ , [70, 71, 72, 73]). Various enhanced image-segmentation
techniques have been proposed in order to ameliorate such problems, like
homology-based preservation techniques [74, 75, 76, 77, 78, 79, 80], topology-
derivative-based techniques [81, 82, 83, 84, 85], and adaptive refinement
techniques [72, 86, 87, 88, 89, 90, 91].

In this contribution we propose an enhancement of the immersed isogeometric
analysis framework of Ref. [19] that rigorously avoids the occurrence of
topological anomalies associated with the B-spline-based image segmentation
procedure. The filtering analysis in Ref. [19] is extended to include the
effect of the level set mesh size. Based on this analysis, it is proposed to
employ truncated hierarchical B-splines (THB-splines) [92, 93, 94] to
discretize the grayscale intensity of the scan data. Inspired by topology
characteristics (_e.g._ , Betti numbers, Euler characteristic) to detect local
topology anomalies [95, 96, 97, 77, 78], this local refinement capability for
the level set representation is complemented with an Euler characteristic
evaluation and moving-window technique to detect local topology anomalies. The
resulting topology-preserving immersed isogeometric analysis is demonstrated
using prototypical test cases in two and three dimensions, considering both
the image-segmentation problem and the complete (stabilized) immersed
isogeometric analysis problem [36]. A scan-based analysis on a representative
data set is considered to demonstrate the practical applicability of the
proposed technique. In this contribution, we restrict ourselves to user-
specified locally-refined meshes, keeping in mind the extension to an adaptive
analysis framework as a possible further development.

This paper is organized as follows. In Section 2 we first extend the filtering
analysis of Ref. [19]. Based on this analysis, in Section 3 we introduce the
topology-preserving extension of the spline-based image segmentation
technique. This technique comprises the moving-window technique to detect
topological anomalies (Section 3.1) and a THB-spline-based discretization of
the level set function (Section 3.2). The adopted immersed isogeometric
analysis framework, including stabilization techniques, is introduced in
Section 4. Numerical examples are then presented in Section 5 to test the
proposed technique. Finally, conclusions and recommendations are drawn in
Section 6.

##  2 Spline-based image segmentation

In this section we review the spline-based image segmentation procedure of
Ref. [19]. The occurrence of topological anomalies on relatively coarse voxel
grids is illustrated and explained using a Fourier analysis. Based on this
analysis, a solution strategy to avoid topological anomalies is proposed.

###  2.1 B-spline level set construction

Consider an $n_{d}$-dimensional image domain, $\Omega_{\rm
img}=[0,L_{1}]\times\ldots\times[0,L_{n_{d}}]$, which is partitioned by a set
of voxels, as illustrated in Figure 1. We denote the voxel mesh by

$\mathcal{V}_{\rm vox}=\\{\Omega_{\rm vox}\subset\Omega_{\rm
img}\mid\exists\mathbf{i}\in\mathbb{Z}_{\geq 0}^{n_{d}},\mbox{ s.t.
}\Omega_{\rm vox}=T_{\mathbf{i}}\circ[0,1]^{n_{d}}\\},$ (1)

where the transformation $T_{\mathbf{i}}:[0,1]^{n_{d}}\rightarrow\Omega_{\rm
vox}$ is defined as

$T_{\mathbf{i}}(\boldsymbol{\xi})=\mbox{diag}(\boldsymbol{\Delta})(\boldsymbol{\xi}+\mathbf{i}),$
(2)

with $\boldsymbol{\Delta}$ the voxel size in each direction and
$\mathbf{i}=(i_{1},\ldots,i_{n_{d}})$ the $n_{d}$-dimensional voxel index. We
denote the number of voxels in the voxel image by $m_{\rm
vox}=\\#\mathcal{V}_{\rm vox}$. The grayscale intensity function, illustrated
in Figure 1(b), is then defined as

$g:\mathcal{V}_{\rm vox}\rightarrow\mathscr{G},$ (3)

with $\mathscr{G}$ the range of the grayscale data (_e.g._ , from 0 to 255 for
8 bit unsigned integers). An approximation of the object $\Omega$ can be
obtained by thresholding the gray scale data,

$\Omega\approx\\{\boldsymbol{x}\in\Omega_{\rm img}|g(\boldsymbol{x})>g_{\rm
crit}\\}\subset\Omega_{\rm img},$ (4)

where $g_{\rm crit}$ is the threshold value (see Figure 1(c)). As a
consequence of the piecewise definition of the grayscale data in equation (3),
the boundary of the segmented object is non-smooth when the grayscale data is
segmented directly. In the context of analysis, the non-smoothness of the
boundary can be problematic, as irregularities in the surface may lead to non-
physical singularities in the solution to the problem.

(a)

(b)

(c)

Figure 1: Representative two-dimensional geometry to illustrate the spline-
based segmentation procedure. _(a)_ The assumed geometry, which in a real-life
scan-based analysis setting is evidently not available. _(b)_ The grayscale
data on a $75\times 75$ voxel grid, which we derived from the exact geometry
by computing the volume fraction of each voxel. In the scan-based analysis
setting this voxel image is obtained directly from the scan data. _(c)_ The
segmented geometry for a grayscale threshold corresponding to 50% of the
volume fraction.

The B-spline segmentation procedure in Ref. [19] – the behavior of which in
the case of linear basis functions closely corresponds with that of marching
volume algorithms [8] – enables the construction of a smooth boundary
approximation based on voxel data. The key idea of this B-spline segmentation
technique is to smoothen the grayscale function (3) by convoluting it using a
B-spline basis of size $n$, $\\{N_{i,k}(\mathcal{V}^{h})\\}_{i=1}^{n}$,
defined over a regular mesh, $\mathcal{V}^{h}$, with fixed element size,
$h_{d}$, per direction. Note that this mesh size can be different from the
voxel size. The B-spline basis can be constructed using the Cox-de Boor
algorithm [98]. We consider full-regularity ($C^{k-1}$-continuous) B-splines
of order $k$, with the order assumed to be constant and isotropic. By
convoluting the grayscale function (3), a smooth level set approximation is
obtained as

$\displaystyle f(\boldsymbol{x})$
$\displaystyle=\sum\limits_{i=1}^{n}N_{i,k}(\boldsymbol{x})a_{i},$
$\displaystyle a_{i}$ $\displaystyle=\cfrac{\int_{\Omega_{\rm
img}}N_{i,k}(\boldsymbol{x})g(\boldsymbol{x}){\rm
d}\boldsymbol{x}}{\int_{\Omega_{\rm img}}N_{i,k}(\boldsymbol{x}){\rm
d}\boldsymbol{x}},$ (5)

where the coefficients $\\{a_{i}\\}_{i=1}^{n}$ are the control point level set
values.

The B-spline level set function corresponding to the voxel data in Figure 1(b)
is illustrated in Figure 2 for the case where the regular mesh,
$\mathcal{V}^{h}$, coincides with the voxel grid (_i.e._ , $\Delta_{d}=h_{d}$
for $d=1,\ldots,n_{d}$) and second order ($k=2$) B-splines. As can be seen,
the object retrieved from the convoluted level set function more closely
resembles the original geometry in Figure 1(a) compared to the voxel
segmentation in Figure 1(c). Also, the boundaries of the domain are smooth as
a consequence of the higher-order continuity of the B-spline basis (see Figure
2(c)). The segmented geometry in Figure 2(c) is constructed using the midpoint
tessellation procedure proposed in Refs. [19, 45], which constructs a
partitioning of the elements that intersect the domain boundary and results in
an accurate parametrization of the interior volume. The reader is referred to
Ref. [19] for a detailed discussion of the properties of the employed B-spline
segmentation technique, and, additionally, to Ref. [45] for details regarding
the midpoint tessellation procedure.

(a) Grayscale data $g(\boldsymbol{x})$

(b) B-spline level set function $f(\boldsymbol{x})$

(c) Segmented geometry

Figure 2: Illustration of the B-spline based segmentation procedure on a
$70\times 70$ voxel grid using second-order B-splines constructed on a regular
mesh that coincides with the voxels. _(a)_ Grayscale data,
$g(\boldsymbol{x})$, of Figure 1(b) shown on the 8 bit signed integer range.
_(b)_ Level set function, $f(\boldsymbol{x})$, computed by equation (5). _(c)_
Segmented domain extracted using the midpoint tessellation procedure outlined
in Refs. [19, 45].

###  2.2 The occurrence of topological anomalies

The spline-based segmentation technique reviewed above has been demonstrated
to yield computational domains that are very well suited for isogeometric
analysis (see, _e.g._ , [36, 58, 72]). However, although the smoothing
characteristic of the technique is frequently beneficial, it may lead to the
occurrence of topological anomalies when the features of the object to be
described are not significantly larger in size than the voxels (_i.e._ , the
Nyquist criterion is not satisfied [19]).

For the object in Figure 1(a) this occurs when a voxel grid of $35\times 35$
is considered, as illustrated in Figure 3 for a second-order ($k=2$) B-spline
basis defined on the $35\times 35$ voxels mesh. As can be seen, small and
slender features, such as the regions highlighted in green in Figure 3(c), are
detectable in the original grayscale data, albeit with a very coarse
representation. The corresponding level set function is smoother, but
introduces topological anomalies in the form of the disappearance of some of
the small and slender features. When we consider the same voxel data, but now
define the B-spline basis for the convolution of the level set function on a
twice as fine mesh (see Figure 3(d)), _i.e._ , $h_{d}=\Delta_{d}/2$ for
$d=1,\ldots,n_{d}$, a topologically correct segmented domain is again
obtained, but still with smoothed boundaries (see Figure 3(e)) compared to the
direct segmentation.

(a) $35\times 35$ grayscale data

(b) Level set for $h_{d}=\Delta_{d}$

(c) Geometry for $h_{d}=\Delta_{d}$

(d) Level set for $h_{d}=\Delta_{d}/2$

(e) Geometry for $h_{d}=\Delta_{d}/2$

Figure 3: Illustration of the B-spline-based segmentation procedure for
$35\times 35$ voxel grayscale data (panel _a_). Using second-order B-splines
constructed on a regular mesh that coincides with the voxels (panel _b_) leads
to topological changes in the segmented domain (panel _c_) compared to the
$70\times 70$ voxels case in Figure 2. When the B-spline basis is constructed
on a mesh that is twice as fine as the voxel grid, a smooth level set function
is obtained (panel _d_) which, after segmentation, correctly represents the
topology of the object (panel _e_).

To elucidate the smoothing behavior of the B-spline segmentation technique, we
generalize the filtering analysis for a univariate B-spline ($n_{d}=1$)
presented in Ref. [19] to the case of non-coinciding voxel and B-spline grid
sizes, _i.e._ , $h\neq\Delta$. Note that, in the univariate setting considered
here, we drop the index for the geometric direction to simplify our notation.
The goal of our analysis is to provide insight into the filtering properties
by considering the level set approximation in the frequency domain.
Additionally, we will obtain an analytical expression for the smoothed level
set function in the spatial domain with a parametrization that provides
insight into the smoothing properties of the level set construction.

(a) $h=1$

(b) $h=0.5$

(c) $h=1$

(d) $h=0.5$

Figure 4: The convolution kernel (7) in the spatial domain (panels _(a)_ and
_(b)_), and in the frequency domain (panels _(c)_ and _(d)_) for $k=2,3,4$
using different mesh sizes. The integration kernel (6) is plotted in the
spatial domain for reference (dashed lines).

We commence our analysis with rewriting the operation (5) as an integral
transform

$\displaystyle f({x})$ $\displaystyle=\int_{\Omega_{\rm
img}}\mathcal{K}(x,y)g({y})\>{\rm d}{y},$ $\displaystyle\mathcal{K}({x},{y})$
$\displaystyle=\sum\limits_{i=1}^{n}\cfrac{N_{i,k}({x})N_{i,k}({y})}{V_{i}},$
(6)

where $\mathcal{K}(x,y)$ is the kernel of the transformation and where
$V_{i}=\int_{\Omega_{\rm img}}N_{i,k}(x)\>{\rm d}x$ is the integral of the
basis function $N_{i,k}$. Following the derivation in Ref. [19] – in which the
essential step is to approximate the B-spline basis functions by rescaled
Gaussians [99] – the integration kernel (6) can be approximated by

$\mathcal{K}(x,y)\approx\kappa(x-y)=\frac{1}{\sigma\sqrt{2\pi}}\exp{\left({-\cfrac{(x-y)^{2}}{2\sigma^{2}}}\right)},$
(7)

where the width of the smoothing kernel is given by

$\sigma=h\sqrt{\frac{k+1}{6}}.$ (8)

This result is very similar to that obtained in Ref. [19], with the important
difference that the width of the kernel depends on the mesh size, $h$, on
which the B-spline is defined, and not on the voxel size, $\Delta$. Note that
the approximate kernel $\kappa$ depends on the difference between the
coordinates $x$ and $y$ only, this in contrast to the kernel $\mathcal{K}$.
Consequently, when the approximate kernel $\kappa$ is considered, the integral
transform (6) becomes a convolution operation.

The dependence of the convolution kernel on the mesh size and on the order of
the B-spline interpolation is illustrated in Figures 4(a) and 4(b). The
B-spline integral transform around $x=0$ computed on a domain of size 10 with
a voxel size 1 is shown for reference, indicating that the Gaussian
approximation improves with an increasing B-spline order. Figure 4 conveys
that, following equation (8), increasing the mesh size and the B-spline order
both increase the zone over which the grayscale data is averaged. The
smoothing width scales linearly with the mesh size, and, for sufficiently
large B-spline orders, with the square root of $k$.

To provide detailed insights into how the filter properties lead to
topological anomalies, we express the convolution operation (6) with the
approximate kernel (7) in the frequency domain as

$\displaystyle F(\xi)=K(\xi)G(\xi),$ (9)

where $F(\xi)$ and $G(\xi)$ are the Fourier transforms of the original
grayscale data and the convoluted level set function, respectively, and where
the Fourier transform of the convolution kernel (7) is given by

$K(\xi)=\exp{\left(-2\pi^{2}\xi^{2}\sigma^{2}\right)}.$ (10)

This frequency-domain form of the kernel is shown in Figures 4(c) and 4(d).
The Fourier-form of the convolution operation in equation (9) conveys that
features corresponding to frequencies for which $K(\xi)$ is close to unity are
preserved in the smoothing operations, whereas features corresponding to
frequencies for which $0<K(\xi)\ll 1$ are filtered.

To further clarify the preservation of features, we consider a one-dimensional
object of size $\ell$, represented by the grayscale function (in both the
spatial and in the frequency domain)

$\displaystyle g(x)$ $\displaystyle=\begin{cases}1&|x|<\ell/2\\\
0&\mbox{otherwise}\end{cases},$ $\displaystyle G(\xi)$
$\displaystyle=\ell\>{\rm
sinc}{\left(\ell\xi\right)}=\frac{\sin{\left(\pi\ell\xi\right)}}{\pi\xi}.$
(11)

The smooth approximation of this feature in the frequency domain follows from
equation (9) as

$F(\omega)=\ell\>{\rm
sinc}{\left(\ell\xi\right)}\exp{\left(-2\pi^{2}\xi^{2}\sigma^{2}\right)}.$
(12)

The corresponding function in the spatial domain can be determined by
expressing the ${\rm sinc}$ function as [100]

${\rm
sinc}(\ell\xi)=\lim\limits_{m\rightarrow\infty}\frac{1}{m}\sum\limits_{n=1}^{m}\cos{\left(2\mu(m,n)\pi\xi\right)},$
(13)

with $\mu(m,n)=\frac{(2n-1)\ell}{4m}$, such that the final form of the
approximated convolution operation (6) can be written as

$f(x)=\frac{\ell}{\sigma\sqrt{2\pi}}\lim\limits_{m\rightarrow\infty}\frac{1}{2m}\sum\limits_{n=1}^{m}\left[\exp{\left(-\frac{\left(x-\mu(m,n)\right)^{2}}{2\sigma^{2}}\right)}+\exp{\left(-\frac{\left(x+\mu(m,n)\right)^{2}}{2\sigma^{2}}\right)}\right].$
(14)

(a) $\hat{\ell}=\frac{\ell}{h}=2$

(b) $\hat{\ell}=\frac{\ell}{h}=2$

(c) $\hat{\ell}=\frac{\ell}{h}=1$

(d) $\hat{\ell}=\frac{\ell}{h}=1$

(e) $\hat{\ell}=\frac{\ell}{h}=\frac{1}{2}$

(f) $\hat{\ell}=\frac{\ell}{h}=\frac{1}{2}$

Figure 5: Smoothed level set approximation (14) of a geometric feature in the
spatial domain (left) and in the frequency domain (right) for various feature-
size-to-mesh ratios, $\hat{\ell}=\ell/h$, and B-spline degrees. The B-spline
level set function (6) is shown in the spatial domain for reference (dashed
lines).

The approximate level set function (14) is illustrated in Figure 5 for various
feature-size-to-mesh ratios, $\hat{\ell}=\ell/h=2,1,\frac{1}{2}$, and B-spline
degrees, $k=2,3,4$. For the considered range of feature-size-to-mesh ratios,
the limit in equation (14) is already approximated well with $m=1$, such that

$f(h\hat{x})\approx\hat{f}_{1}(\hat{x})=\sqrt{\frac{3\hat{\ell}^{2}}{4\pi(k+1)}}\left[\exp{\left(-\frac{3\left(4\hat{x}-\hat{\ell}\right)^{2}}{16(k+1)}\right)}+\exp{\left(-\frac{3\left(4\hat{x}+\hat{\ell}\right)^{2}}{16(k+1)}\right)}\right],$
(15)

where $\hat{x}=x/h$. The value of the smoothed level set function at
$\hat{x}=0$ follows as

$\hat{f}_{1}(0)=\sqrt{\frac{3\hat{\ell}^{2}}{\pi(k+1)}}\exp{\left(-\frac{3\hat{\ell}^{2}}{16(k+1)}\right)}\approx\hat{\ell}\sqrt{\frac{3}{\pi(k+1)}},$
(16)

which conveys that the maximum value of the smoothed level set depends
linearly on the relative feature size $\hat{\ell}$, and decreases with
increasing B-spline order.

The top row in Figure 5 shows the case for which the considered feature is
twice as large as the mesh size, _i.e._ , $\hat{\ell}=2$. Figure 5(a)
illustrates that the sharp boundaries of the original grayscale function are
significantly smoothed, which is also apparent from the frequency domain plot
in Figure 5(b), which shows that the high frequency modes required to
represent the sharp boundary are filtered out by the smoothing operation. In
Figure 5(a), the decrease in the maximum level set value as given by equation
(16) is observed. When the level set function is segmented by a threshold of
$g_{\rm crit}=0.5$, a geometric feature that closely resembles the original
one is recovered.

The middle and bottom rows of Figure 5 illustrate cases where the feature
width is not significantly larger than the mesh size. For the case where the
feature size is equal to the size of the mesh, the maximum of the level set
function drops significantly compared to the the case of $\hat{\ell}=2$. When
considering second-order B-splines, the maximum is still marginally above
$g_{\rm crit}=0.5$. Although the recovered feature is considerably smaller
than the original one, it is still detected in the segmentation procedure.
When increasing the B-spline order, the maximum value of the level set drops
below the segmentation threshold, however, indicating that the feature will no
longer be detected. As a consequence, the B-spline segmentation procedure then
induces a topological alteration. When decreasing the feature length further,
as illustrated in the bottom row of Figure 5, topological changes are
encountered at lower segmentation thresholds. For the case of $g_{\rm
crit}=0.5$, also the quadratic B-splines would lead to a loss of the feature.

In summary, the above analysis shows that topological anomalies occur when the
relative feature length scale, $\hat{\ell}=\ell/h$, becomes too small. The
smallest features in a voxel data file are of the size of a single voxel,
_i.e._ , $\ell=\Delta$. Hence, topological features are lost when the mesh
size on which the B-spline level set is constructed is relatively large
compared to the voxel size. For moderate B-spline orders, this practically
means that for features with the size of a single voxel, the mesh of the
B-spline level set should be a uniform refinement of the voxel mesh, such that
the relative feature size $\hat{\ell}$ is sufficiently large. With this mesh
refinement, the resulting level set function will still be
$C^{k-1}$-continuous, but higher-frequent modes are present in the refined
level set function. In places where the geometric features are sufficiently
large, this is in principle not desirable. Ideally, one only should refine the
B-spline level set function in places where this is needed to preserve the
topology, _i.e._ , around small features. In the next section, we propose a
fully-automated topology-preserving B-spline segmentation strategy that
refines the B-spline level set function in such a manner that topological
anomalies are avoided.

##  3 Topology-preserving image segmentation using THB-splines

In this section we present a topology-preserving B-spline-based image
segmentation strategy relying on the technique proposed in Ref. [19]. The
proposed strategy consists of two steps (schematically illustrated in Figure
6). In the first step, which we will discuss in detail in Section 3.1, a
moving-window strategy is applied to locally compare the topology of the
original voxel data and its smoothly segmented counterpart. The result of this
first step is a function that indicates regions where topological anomalies
occur. In the second step, which is discussed in Section 3.2, truncated
hierarchical (TH)B-splines are employed to locally repair the topology by mesh
refinement, following the analysis presented in Section 2. To demonstrate the
proposed image segmentation technique, various test cases, including the
problematic scenario considered in Section 2, will be discussed in Section
3.3.

Figure 6: Illustration of the topology-preserving image segmentation
procedure. The original grayscale image (panel _a_) is segmented in two ways,
_viz._ , directly by thresholding the voxel data (panel _b_), and through the
B-spline smoothing strategy (panels _f_ and _g_). A moving-window technique
then locally compares the topology between the two segmented images, which
results in an indicator function (panel _d_) to mark topological differences.
THB-spline-based refinements are then introduced to locally increase the
resolution of the smooth level set function (panel _i_), thereby preserving
the topology of the original image.

###  3.1 Moving-window topological anomaly detection

In this section we detail the moving-window strategy to detect topological
anomalies. In Section 3.1.1 we commence with the definition of the moving
window concept, as is commonly used in image processing techniques (_e.g._ ,
[101, 102, 103]). Subsequently, in Section 3.1.2 we introduce the window-
comparison operator used to identify topological changes. A masking operation
is finally introduced in Section 3.1.3 to distinguish between boundary changes
and topological changes.

####  3.1.1 The moving window

To detect local topological changes, local views on the original voxel data
and its smoothed counterpart will be compared. The windows are created by
considering the $r$-neighborhood of each voxel in the original image, _i.e._ ,

$\displaystyle\Omega_{\rm win}$ $\displaystyle=\mathcal{N}_{r}(\Omega_{\rm
vox})$ $\displaystyle\forall\Omega_{\rm vox}\in\mathcal{V}_{\rm vox}.$ (17)

The 0-neighborhood of a voxel is defined as the voxel itself, _i.e._ ,
$\mathcal{N}_{0}(\Omega_{\rm vox}):=\Omega_{\rm vox}$, and the
$r$-neighborhood for $r\geq 1$ is defined recursively as the union of all
voxels that share a vertex with the $(r-1)$-neighborhood. This window
definition is illustrated in Figure 7 for a $10\times 10$ voxel grid.

As discussed in Section 2, it is desirable to keep the mesh refinement to
repair topological anomalies as local as possible. This means that the moving
window should be as small as possible, but large enough to detect topological
anomalies. Hence, the window should be larger than the size of the geometric
features that are subject to topological changes. From the analysis in Section
2.2 it follows that features with a size similar to that of the voxels, _i.e._
, $\ell\approx\Delta$, can be subject to topological changes if the smoothing
is performed on a mesh with a size similar to that of the voxels, _i.e._ ,
$h\approx\Delta$. Practically this means that a window with only a few rings,
_i.e._ , $\ell\approx\Delta<(2r+1)\Delta\ll L$, is expected to be an adequate
choice. The influence of the window size will be studied for the example
discussed in Section 3.3.

(a)

(b)

Figure 7: Illustration of the moving window $\Omega_{\rm win}$ centered at a
voxel $\Omega_{\rm vox}$, with _(a)_ a one-neighborhood ($r=1$) and _(b)_ a
two-neighborhood ($r=2$).

####  3.1.2 Window image comparison

The moving-window technique indicates whether, for the window focused at every
voxel, a topological difference is observed between the directly segmented
image, $V$, and the geometry related to the smoothed image, $S$ (see Figure
6). We herein employ the midpoint tessellation strategy proposed in Ref. [45]
to construct a geometry parametrization, as illustrated in Figure 8(a).
Details of this procedure will be discussed in the context of the immersed
isogeometric analysis framework discussed in Section 4. To enable the window
topology comparison, it is convenient to employ the same representation
concept for both images, $V$ and $S$. Therefore, the tessellated geometry is
first voxelized by segmentation of the B-spline level set function on a grid
that is $n_{\rm sub}$-times uniformly refined with respect to the original
voxel grid, as shown in Figure 8(b). This refinement level should be chosen
such that the voxelization of the smoothed image does not induce topological
changes compared to the geometry on which the simulation is performed. Since
the employed midpoint tessellation strategy also relies on a recursive
subdivision approach (of the computational cells), the number of voxel
subdivisions, $n_{\rm sub}$, can be related directly to the number of
recursive refinements used for the midpoint tessellation (see Figure 8).

(a) Midpoint tessellation with a recursion depth of $\varrho_{\rm max}=2$

(b) Voxelization on a subdivided grid with $n_{\rm sub}=2$

Figure 8: Schematic of _(a)_ the midpoint tessellation procedure of Refs. [19,
45] and _(b)_ the voxelization on a grid that is refined to match the
recursion depth of the tessellation procedure.

To characterize topological similarity, we define the window topology-
comparison operator:

$\mathscr{C}(V,S)=\begin{cases}\mbox{true}&\mbox{if $V$ and $S$ are
topologically equivalent,}\\\ \mbox{false}&\mbox{if $V$ and $S$ are
topologically different.}\end{cases}$ (18)

To compute this boolean operator, the filled voxels in the images $V$ and $S$
are divided into connected regions. For the image $V$, for example, the
regions are denoted by $R\in\mathcal{R}_{V}$, where $\mathcal{R}_{V}$ is the
set of all connected regions. We herein employ vertex-connectivity, meaning
that voxels sharing a vertex are considered to be connected, but the same
procedure can be applied using face-connectivity. For each connected region,
$R$, we determine the Euler characteristic, $\chi(R)$, which is defined as one
minus the number of connected holes, $\\#\mathcal{V}_{R}$, in that region.
Following the definition in Ref. [104], the Euler characteristic of the window
image is given by

$\chi(V)=\sum\limits_{R\in\mathcal{R}_{V}}\chi(R)=\sum\limits_{R\in\mathcal{R}_{V}}\left(1-\\#\mathcal{V}_{R}\right)=\\#\mathcal{R}_{V}-\sum\limits_{R\in\mathcal{R}_{V}}\\#\mathcal{V}_{R}.$
(19)

We note that various methods of evaluating the Euler characteristic have been
rigorously studied in the field of homology, see, _e.g._ , Refs. [105, 106,
107, 108]. In our work we consider the method of Ref. [104] which has been
implemented in scikit-image [109], an open-source image processing library for
Python.

To construct the boolean operator (18), we define sets of regions with a
specific Euler characteristic,
$\bar{\chi}\in\boldsymbol{\chi}_{V}=\\{\chi(R)\mid
R\in\mathcal{R}_{V}\\}\subset\mathbb{Z}$, in an image as

$\displaystyle\mathcal{R}_{V}^{\bar{\chi}}=\\{R\in\mathcal{R}_{V}\mid\chi(R)=\bar{\chi}\\}.$
(20)

Figure 9(c) shows illustrative examples of these region sets. Using the region
sets (20), the following indicator function for topological anomalies is
proposed:

$\mathscr{C}(V,S)=\left[\underset{\bar{\chi}\in\boldsymbol{\chi}_{V}\cup\boldsymbol{\chi}_{S}}{\bigwedge}\left(\\#\mathcal{R}_{V}^{\bar{\chi}}\equiv\\#\mathcal{R}_{S}^{\bar{\chi}}\right)\right]\bigwedge\left[\underset{\bar{\chi}\in\boldsymbol{\chi}_{V^{\prime}}\cup\boldsymbol{\chi}_{S^{\prime}}}{\bigwedge}\left(\\#\mathcal{R}_{V^{\prime}}^{\bar{\chi}}\equiv\\#\mathcal{R}_{S^{\prime}}^{\bar{\chi}}\right)\right]$
(21)

This indicator states that the number of regions with a particular Euler
characteristic is equal in the images $V$ and $S$, as well as in their
complements $V^{\prime}$ and $S^{\prime}$. Note that, by construction, the
operator (21) is invariant to the complement operation, _i.e._ ,
$\mathscr{C}(V,S)=\mathscr{C}(V^{\prime},S^{\prime})$, which makes the
comparison operator objective with respect to the definition of the solid and
void regions. One may note that if $\chi_{V}\neq\chi_{S}$ or
$\chi_{V^{\prime}}\neq\chi_{S^{\prime}}$, then $\mathscr{C}(V,S)$
automatically reverts to false.

(a) $\mathcal{R}_{V}^{1}=\\{{\color[rgb]{0,0,255}R_{1}}\\}$,
$\mathcal{R}_{V}^{0}=\\{{\color[rgb]{0,128,0}R_{2}}\\}$

(b) Trimmed topology

(c) $\mathcal{R}_{S}^{1}=\\{{\color[rgb]{0,0,255}R_{1}}\\}$,
$\mathcal{R}_{S}^{0}=\\{{\color[rgb]{0,128,0}R_{2}}\\}$ (d)
$\mathcal{R}_{V^{\prime}}^{1}=\\{{\color[rgb]{0,0,255}R_{1}},{\color[rgb]{0,128,0}R_{2}}\\}$
(e) Trimmed topology (f)
$\mathcal{R}_{S^{\prime}}^{1}=\\{{\color[rgb]{0,0,255}R_{1}},{\color[rgb]{0,128,0}R_{2}}\\}$
(g) Trimmed topology (h)
$\mathcal{R}_{S}^{1}=\\{{\color[rgb]{0,0,255}R_{1}},{\color[rgb]{0,128,0}R_{2}}\\}$
(i) Trimmed topology (j)
$\mathcal{R}_{S^{\prime}}^{1}=\\{{\color[rgb]{0,0,255}R_{1}}\\}$ (k) Trimmed
topology (l)
$\mathcal{R}_{S}^{1}=\\{{\color[rgb]{0,0,255}R_{1}},{\color[rgb]{200,0,0}R_{3}}\\}$,
$\mathcal{R}_{S}^{0}=\\{{\color[rgb]{0,128,0}R_{2}}\\}$ (m) Trimmed topology
(n)
$\mathcal{R}_{S^{\prime}}^{1}=\\{{\color[rgb]{0,0,255}R_{1}},{\color[rgb]{0,128,0}R_{2}}\\}$
Figure 9: Illustration of the comparison operator (21) for the images $V$ and
$V^{\prime}$ (original gray scale data) in the first and fourth column and the
image $S$ and $S^{\prime}$ (B-spline-based segmented domain extracted using
the midpoint tessellation procedure with $\varrho_{\rm max}=3$) in the second
and fifth column. The third and sixth column are the $n_{\rm sub}=3$
voxelizations of the segmented domain. The images in the first row result in
$\mathscr{C}(V,S)=\mathscr{C}(V^{\prime},S^{\prime})=1$ as $V$ and $S$ are
topologically equivalent. The images in the second row results
$\mathscr{C}(V,S)=\mathscr{C}(V^{\prime},S^{\prime})=0$ as $V$ and $S$ are
topologically different. The images in the third row results in
$\mathscr{C}(V,S)=\mathscr{C}(V^{\prime},S^{\prime})=0$ as $V$ and $S$ are
topologically different due to a boundary spill over.

Figure 9(c) illustrates the comparison operator (21), where use has been made
of a $n_{\rm sub}=3$ sampling of the B-spline-based segmented domain extracted
using the midpoint tessellation procedure with $\varrho_{\rm max}=3$. The
first row in Figure 9(c) shows a case where the topology matches, despite the
substantial changes in geometry. The second row in Figure 9(c) shows a case
where the elliptical region is missing, which is a typical case of a
topological anomaly. The third row in Figure 9(c) shows a case in which the
comparison operator returns false, which is caused by the appearance of a
boundary spillover due to the smoothing procedure at the right bottom border.
From the perspective of the window, the observed change indeed classifies as a
topological change. However, when considering the change from the perspective
of the complete image, the observed difference comes from a boundary that
moves into the view of the window under consideration. This boundary movement
classifies as a shape (geometry) change, and not as a topological change. To
correctly account for this type of changes, in the next section we propose an
image masking operation.

####  3.1.3 Window image masking

To mask the window changes associated with moving boundaries, as a
preprocessing step to the comparison operation (21), the smoothed image $S$ is
masked using

$\displaystyle F$ $\displaystyle=\mathcal{F}(V,S)=(M\cap V)\cup(M^{\prime}\cap
S),$ (22)

where $F$ is the masked image and the mask $M$ depends on arguments $V$ and
$S$, _i.e.,_ $M=\mathcal{M}(V,S)$. This mask corresponds to a set-indication
function which is one in places where a boundary change is detected and zero
everywhere else. The mask is illustrated in Figure 10 for the case considered
in the third row of Figure 9(c). The key idea behind this masking operation is
that in regions where shape changes occur, the masked image $F$ is replaced by
the original voxel image $V$ so that the changes associated with boundary
movement are effectively reverted.

Figure 10: The symmetric difference (24) between the images $V$ and $S$ in
Figure 9(c), with the boundary mask $M=\mathcal{M}(V,S)$ shown in red.

To identify the locations of the changes between the images $V$ and $S$, we
consider the mask to be a subset of the symmetric difference between the
original and the smoothed image, _i.e._ , $M\subseteq V\Delta S=\left(V\cap
S^{\prime}\right)\cup\left(S\cap V^{\prime}\right)$ (see Figure 10). We now
only mask the regions in the symmetric difference which reside completely in
the outer ring of the window _and_ have an Euler characteristic of one (there
are no voids), that is,

$\mathcal{M}(V,S)=\left\\{R\in\mathcal{R}_{V\Delta S}^{1}\mid R\cap
W_{r-1}=\emptyset\right\\},$ (23)

where $W_{r-1}$ is equal to zero in the outer voxel ring of the window and one
everywhere else. Note that since the symmetric difference is unchanged when
the complement images are considered, _i.e._ ,

$(V^{\prime})\Delta(S^{\prime})=\left(V^{\prime}\cap
S\right)\cup\left(S^{\prime}\cap V\right)=V\Delta S,$ (24)

for the mask (23) it holds that

$\mathcal{M}(V^{\prime},S^{\prime})=\mathcal{M}(V,S)=M.$ (25)

From this property of the mask it then follows that the complement of the
masked image $F$ is equal to the masked complement images $V^{\prime}$ and
$S^{\prime}$:

$\displaystyle F^{\prime}=\mathcal{F}(V,S)^{\prime}$ $\displaystyle=(M\cap
V)^{\prime}\cap(M^{\prime}\cap S)^{\prime}=(M^{\prime}\cup
V^{\prime})\cap(M\cup S^{\prime})$ (26)
$\displaystyle=\left[M^{\prime}\cap(M\cup
S^{\prime})\right]\cup\left[V^{\prime}\cap(M\cup
S^{\prime})\right]=(M^{\prime}\cap S^{\prime})\cup(V^{\prime}\cap
M)\cup(V^{\prime}\cap S^{\prime})$ $\displaystyle=(M^{\prime}\cap
S^{\prime})\cup(V^{\prime}\cap M)\cup(V^{\prime}\cap S^{\prime}\cap
M)\cup(V^{\prime}\cap S^{\prime}\cap M^{\prime})=(M\cap
V^{\prime})\cup(M^{\prime}\cap S^{\prime})$
$\displaystyle=\mathcal{F}(V^{\prime},S^{\prime})$

The choice to identify shape changes associated with moving boundaries by
definition (23) is motivated by the idea that if one considers a boundary in
the global image with small curvature, the image smoothing operation discussed
in Section 2 will typically retain the movement of the boundary within one
voxel spacing. When the curvature of the boundary is locally high, the
boundary movement can be larger than a single voxel, which would lead to
falsely identifying a shape change as a topological change. In our algorithm
this would lead to a refinement of the level set function which would not be
strictly required to preserve topology. Avoiding such auxiliary refinements
would, however, severely complicate the algorithm. Moreover, since they
generally result in an improved geometry representation at high-curvature
boundaries, it is considered desirable not to avoid them.

(a) $F=\mathcal{F}(V,S)$

(b) $M\cap V$

(c) $M^{\prime}\cap S$

(d) $F^{\prime}=\mathcal{F}(V^{\prime},S^{\prime})$

(e) $M\cap V^{\prime}$

(f) $M^{\prime}\cap S^{\prime}$

Figure 11: An illustration of the filtering operation (22) for an image $V$
and the corresponding smoothed image $S$ (see Figure 12) and for their
complements.

(a) $V$

(b) $S$

(c) $V\Delta S$

(d) $V^{\prime}$

(e) $S^{\prime}$

(f) $V^{\prime}\Delta S^{\prime}$

Figure 12: An example of _(a)_ the original segmented image $V$, _(b)_ the
spline-based segmented image $S$, and _(c)_ the symmetric difference. Panels
_(d)_ and _(e)_ are the complements of $V$ and $S$, respectively. Panel _(f)_
is the symmetric difference of the complements, which is identical to that in
panel _(c)_.

The masking operation (22) is illustrated in Figure 11 for an exemplifying
image $V$ and its smoothed version $S$ (as defined in Figure 12), as well as
for the complements of these images. The symmetric difference between the
images $V$ and $S$ is shown in Figure 12(c), which, in agreement with equation
(24), is identical to the symmetric difference between $V^{\prime}$ and
$S^{\prime}$ in Figure 12(f). The masked regions are color coded in red,
indicating that only completely solid regions in the outer ring are considered
in the mask. The masked images $F=\mathcal{F}(V,S)$ and
$F^{\prime}=\mathcal{F}(V^{\prime},S^{\prime})$ are shown in Figures 11(a) and
11(d), from which it is observed that, in accordance with equation (26), the
complement of the masked image $F$ is equal to the masked complement image
$\mathcal{F}(V^{\prime},S^{\prime})$.

In Figure 13 we consider the case of the third row of Figure 9(c). As
discussed above, when the images $V$ and $S$ are compared directly, the
comparison operator (21) marks these images to be topologically different on
account of the shape change associated with the boundary movement. When
comparing the image $V$ with the masked image $F=\mathcal{F}(V,S)$, as shown
in Figure 13, the images are considered to be topologically equivalent, which,
in this case, is considered as the correct indicator result.

(a) $V$

(b) Trimmed topology

(c) $F=\mathcal{F}(V,S)$

Figure 13: An illustration of the boundary masking procedure applied to the
case in the third row of Figure 9(c). The image $V$ in panel _(a)_ is compared
with the masked image $F$ in panel _(c)_ to obtain the indicator function
$\mathscr{C}(V,S)$. In contrast to the result in Figure 9(c), the masking
operation results in a $\mathscr{C}(V,S)=1$, indicating that the images $V$
and $S$ are topologically equivalent.

###  3.2 THB-spline-based refinement strategy

After identification of the windows in which topological changes occur, the
mesh on which the level set function is constructed is locally refined to
resolve these anomalies (see Figure 6). More specifically, a support-based
refinement procedure is employed, in which all basis functions that are
supported on a voxel that displays a topological anomaly are replaced by basis
functions from the next higher hierarchical level.

We herein employ truncated hierarchical B-splines [92, 93, 94] to construct a
spline basis over locally refined meshes. The THB-spline construction is
illustrated in Figure 14, which, for the sake of generality, considers the
case of multiple refinement levels, with the level $\ell=0$ corresponding to
the coarsest elements, and the level $\ell=\ell_{\rm max}$ to the most refined
elements. We denote the region covered by elements that are at least $\ell$
times refined by $\Omega^{\ell}$ (note that the refinement regions are nested,
_i.e._ ,
$\Omega^{\ell}\subseteq\Omega^{\ell-1}\subseteq\ldots\subseteq\Omega^{0}=\Omega_{\rm
vox}$). The locally refined mesh corresponding to these refinement regions is
denoted by $\mathcal{V}_{\rm vox}$.

To construct the truncated hierarchical B-spline basis, we consider $\ell_{\rm
max}$ uniform refinements, $\mathcal{V}_{\rm vox}^{\ell}$, of the original
voxel mesh $\mathcal{V}_{\rm vox}^{0}$. The mesh size of the original voxel
mesh is denoted by $h$ and that of its refinements by $2^{-\ell}h$. With each
level we associate a mesh

$\displaystyle\mathcal{V}^{\ell}=\left\\{v\in\mathcal{V}_{\rm vox}^{\ell}\mid
v\cap\Omega^{\ell}\neq\emptyset\right\\}.$ (27)

We define a B-spline basis of degree $k$ and regularity $\alpha$ over each of
these meshes as

$\mathcal{B}(\mathcal{V}^{\ell})=\\{N\in\mathcal{B}(\mathcal{V}_{\rm
vox}^{\ell})\,\mid\,\operatorname{supp}{(N)}\cap\Omega^{\ell}\neq\emptyset\\},$
(28)

where $\mathcal{B}(\mathcal{V}_{\rm vox}^{\ell})$ is the B-spline basis on
$\mathcal{V}_{\rm vox}^{\ell}$.

To construct the truncated hierarchical B-spline basis, splines from the bases
over the uniform meshes (28) are selected and truncated. On the most refined
level, _i.e._ , at $\ell=\ell_{\rm max}$, all basis functions that are
completely inside $\Omega_{\rm vox}^{\ell_{\rm max}}$ are selected:

$\mathcal{H}(\mathcal{V}^{\ell_{\rm
max}})=\\{N\in\mathcal{B}(\mathcal{V}^{\ell_{\rm
max}})\mid\operatorname{supp}{(N)}\subseteq\Omega^{\ell_{\rm max}}\\}.$ (29)

At coarser levels, _i.e._ , $0\leq\ell<\ell_{\rm max}$, the functions that are
completely inside the domain $\Omega^{\ell}$, but not completely inside the
refined domain $\Omega^{\ell+1}$, are selected and truncated:

$\mathcal{H}(\mathcal{V}^{\ell})=\\{\text{trunc}{(N)}\mid
N\in\mathcal{B}(\mathcal{V}^{\ell}),\,\operatorname{supp}{(N)}\subseteq\Omega^{\ell},\,\operatorname{supp}{(N)}\nsubseteq\Omega^{\ell+1}\\}$
(30)

The truncation operation reduces the support of the B-spline functions by
projecting away basis functions retained from the refined levels (see Ref.
[92, 110] for details). The THB-spline basis then follows as

$\mathcal{H}(\mathcal{V})=\bigcup\limits_{\ell=0}^{\ell_{\rm
max}}\mathcal{H}(\mathcal{V}^{\ell}).$ (31)

This selection and truncation procedure is illustrated in Figure 14. It is
noted that THB-splines satisfy the partition of unity property (in contrast to
non-truncated hierarchical B-splines), which is an important property from the
perspective of the image smoothing procedure (5) as it guarantees that the
average grayscale intensity is preserved [19]. In this work we employ the THB-
spline implementation in the Python-based open source numerical library Nutils
[111], which is based on the element-wise construction discussed in Ref.
[112].

Figure 14: An illustration of truncated hierarchical B-splines. The left
column shows the hierarchical levels of a mesh $\mathcal{V}_{\rm vox}$, while
the right column illustrates the concept for a one-dimensional voxel domain
$\Omega_{\rm vox}$.

###  3.3 Examples

To illustrate the topology preserving segmentation strategy presented above we
consider the example presented in the workflow in Figure 6. The voxelized
smooth image with $n_{sub}=2$ is shown in Figure 6g, which closely resembles
the tessellated image with a bi-sectioning depth of $\varrho_{\rm max}=3$ in
Figure 3(c). By comparison with the voxel segmentation in Figure 6b, it is
clearly observed that topological changes occur at two locations. These
locations are highlighted in green in Figure 15c.

Figure 6d presents the comparison indicator for a window size of $3\times 3$
($r=1$). It is observed that the topological changes are indeed detected. Note
that for this example, the mask operation is also active. As an example of the
mask operation, Figure 6h shows a region with a boundary change before the
mask operation. It is observed that the indicator function for topological
anomalies (21) results in zero (false) when compared to the voxel region in
Figure 6c. Figure 6e shows the region after the mask operation. After
application of the mask, the comparison operator results in one (true).

Figure 15: Example of the moving window strategy with a $3\times 3$ window.
Topological changes are correctly detected in the green boxes. In the red
window, an L-shaped exclusion which extends beyond the outer ring of the
window is falsely detected as a topological change.

In Figure 6d it is seen that an additional region, highlighted in red in
Figure 15c, is also marked as a topological change. Strictly speaking, this
would not be necessary, as in both images a hole is present in that region.
Figure 15b shows the window causing this behavior. In the voxel image, the
L-shaped exclusion is not detected as a hole in a region, but as a void region
splitting two solid regions. On the level of the window, this is a
topologically ambiguous case in the sense that one needs to look outside of
the window to see whether the exclusion extends beyond the window. As shown in
Figure 6i and in the zoom in Figure 15e, the effect of refining the level set
in this region is that the voxel geometry is better captured.

(a)

(b)

(c)

Figure 16: The comparison indicator function for the case considered in Figure
6 with a window size of _(a)_ $3\times 3$ ($r=1$), _(b)_ $5\times 5$ ($r=2$),
and _(c)_ $7\times 7$ ($r=3$).

In Figure 16 the influence of the window size is examined. As can be seen, the
refinement regions increase in size with increasing window size. Since the
$3\times 3$ window already adequately corrects the topological anomalies, the
growth of the refinement regions is unnecessary. As can be seen, the effect of
a larger window size on the corrected images is minimal, as elaborated in
Section 2.2. In Figure 16(b) and 16(c) we also observe a case where a shape
change is marked as a topological change, which, as discussed in Section
3.1.3, is caused by the high curvature of the boundary. A zoom of a typical
window in which this occurs is shown in Figure 17. It is also observed in
Figure 16(c) that for the incorrectly identified topology change discussed
above, increasing the window size results in an indicator function in the form
of a ring. This situation is not further explored in this work, but would
require tailoring of the refinement marking strategy to ensure that the
interior of the ring is refined.

(a) $V$

(b) $S$

(c) $\mathcal{F}(V,S)$

Figure 17: A $5\times 5$ window showing an example of a shape change
associated with a high curvature region. Since this shape change extends
beyond the outer ring of the window, it is not masked and hence it is
(incorrectly) detected as a topological change.

##  4 The isogeometric finite cell method

(a)

(b)

Figure 18: Schematic representation of _(a)_ the physical domain $\Omega$
(gray) with boundary $\partial\Omega$ which is embedded in the ambient domain
$\mathcal{A}$, and _(b)_ the ambient domain mesh
$\mathcal{T}_{\mathcal{A}}^{h}$ and the background mesh $\mathcal{T}^{h}$
(yellow), with mesh size parameter $h$.

To provide a basis for the boundary value problems considered in Section 5, in
this section the abstract formulation for the isogeometric finite cell method
[61] is introduced. We consider a physical domain, $\Omega$, formed by the
topology-preserving segmentation procedure outlined above. The domain and its
boundary, $\partial\Omega$, are immersed into an ambient domain
$\mathcal{A}\supset\Omega$ as shown in Figure 18. In the remainder we consider
the ambient domain to coincide with the image (scan) domain, _i.e._ ,
$\mathcal{A}=\Omega_{\rm img}$.

We suppose that the problem under consideration is described by a field
variable $u$ – which can be scalar-valued or vector-valued – and the weak
formulation

$\left\\{\begin{aligned} &\text{Find}\>u\in W\>\text{such that:}\\\
&a(u,v)=f(v)\quad\forall v\in V,\end{aligned}\right.$ (32)

where $W$ is the trial (solution) space, $V$ is the test space, $a:W\times
V\rightarrow\mathbb{R}$ is a continuous bilinear form and
$b:V\rightarrow\mathbb{R}$ is a continuous linear form. The (isogeometric)
finite cell method provides a general framework for constructing the finite
dimensional subspaces $W^{h}\subset W$ and $V^{h}\subset V$, where the
superscript $h$ refers to the mesh parameter associated with the ambient
domain mesh, $\mathcal{T}^{h}_{\mathcal{A}}$, on which the approximation to
the field variable $u$ is computed. Note that the mesh
$\mathcal{T}^{h}_{\mathcal{A}}$ can be different from the level set mesh
$\mathcal{V}^{h}$ discussed in Section 2, as the mesh resolution requirements
following from the approximation of the field variable $u$ generally differ
from those for the level set function.

We define the background mesh as all elements in the ambient domain that touch
the physical domain, _i.e._ ,

$\mathcal{T}^{h}:=\\{K\in\mathcal{T}^{h}_{\mathcal{A}}:K\cap\Omega\neq\emptyset\\},$
(33)

and the interior mesh of the domain $\Omega$ as

$\mathcal{T}^{h}_{\Omega}:=\\{K\cap\Omega:K\in\mathcal{T}^{h}_{\mathcal{A}}\\},$
(34)

where the elements in the background mesh are trimmed to the physical domain
(see Ref. [45] for details). The physical domain boundary mesh is defined as

$\mathcal{T}^{h}_{\partial\Omega}:=\\{K\cap\partial\Omega:K\in\mathcal{T}^{h}_{\mathcal{A}}\\}.$
(35)

The finite dimensional subspaces $W^{h}$ and $V^{h}$ are then constructed
using THB-splines (as elaborated in Section 3.2). We denote the THB-spline
space of degree $k$ and regularity $\alpha$, constructed over the locally-
refined ambient domain $\mathcal{T}^{h}_{\mathcal{A}}$, by

$S^{k}_{\alpha}(\mathcal{A})=\\{N\in C^{\alpha}(\mathcal{A})\,:\,N|_{K}\in
P^{k}(K),\forall K\in\mathcal{T}^{h}_{\mathcal{A}}\\},$ (36)

where $P^{k}(K)$ is the collection of $n_{d}$-variate polynomials on the
element $K\subset\mathbb{R}^{n_{d}}$. The approximation spaces are then
obtained by restricting the THB-splines in $S^{k}_{\alpha}(\mathcal{A})$ to
the physical domain $\Omega$:

$W^{h}=V^{h}=\\{N|_{\Omega}:N\in{}\mathcal{S}^{k}_{\alpha}(\mathcal{A})\\}.$
(37)

Due to the non-mesh-conforming character of the (isogeometric) finite cell
method, it is infeasible to impose Dirichlet boundary conditions by (strongly)
constraining functions in the spaces (37). Instead, Dirichlet conditions are
imposed weakly through Nitsche’s method [46, 48]. By employing a mesh-
dependent consistent stabilization term, a well-posed Galerkin problem is
obtained:

$\left\\{\begin{aligned} &\text{Find}\>u^{h}\in W^{h}\>\text{such that:}\\\
&a^{h}(u^{h},v^{h})=b^{h}(v^{h})\quad\forall v^{h}\in
V^{h}\end{aligned}\right.$ (38)

In this problem the bilinear form $a^{h}:W^{h}\times V^{h}\to\mathbb{R}$ and
linear form $b^{h}:V^{h}\to\mathbb{R}$ are the finite dimensional versions of
the operators in (32), augmented with the above-mentioned Nitsche terms (which
will be specified in Section 5).

Although the Galerkin problem (38) closely resembles that of mesh-conforming
finite element formulations, the immersed setting requires a dedicated
consideration of various computational aspects. In the context of this work,
the following aspects are particularly noteworthy:

##### Ghost-penalty and skeleton-penalty stabilization

To avoid ill-conditioning associated with small volume-fraction trimmed cells,
we apply ghost-penalty stabilization [54]. The idea behind this stabilization
technique is to penalize the jump in the (higher-order) normal gradients of
the solution along all edges in the ghost mesh

$\mathcal{F}^{h}_{\rm ghost}=\\{\partial K\cap\partial
K^{\prime}|K,K^{\prime}\in\mathcal{T}^{h},K\cap\partial\Omega\neq\emptyset,K\neq
K^{\prime}\\},$ (39)

by augmenting the bilinear form with an additional ghost-penalty term. This
ghost-penalty term, which will be detailed for the specific problems
considered in Section 5, also enables scaling of the Nitsche penalty term by
the reciprocal mesh size parameter of the background mesh (independent of the
trimmed-element configurations) [50].

For the flow problems considered in this work, we employ equal-order
discretizations. To make the considered mixed velocity/pressure-
discretizations inf-sup stable, we apply skeleton-stabilization [36] to the
pressure space along all edges in the skeleton mesh

$\mathcal{F}^{h}_{\rm skeleton}=\\{\partial K\cap\partial
K^{\prime}|K,K^{\prime}\in\mathcal{T}^{h},K\neq K^{\prime}\\}.$ (40)

The skeleton-penalty term with which the bilinear form is augmented will be
specified in Section 5.

##### Numerical integration on trimmed elements

To evaluate integrals over the trimmed elements, we consider a recursive
octree bisectioning strategy (see, _e.g._ , Ref. [20, 19]), with the maximum
number of bisections equal to $\varrho_{\rm max}$. On the lowest level of
bisectioning, _i.e._ , $\varrho=\varrho_{\rm max}$, the midpoint tessellation
procedure detailed in Ref. [45] is employed to construct an explicit
parametrization of the trimmed boundary. An illustration of the octree
bisectioning procedure with the midpoint tessellation is shown in Figure 8(a).

Considering equal-order (Gauss) integration schemes on all sub-cells in the
tessellated trimmed elements leads to high computational costs, in particular
when three-dimensional simulations are considered [45]. Various methods have
been proposed to reduce the computational cost, _e.g._ , smart octree methods
[43], moment fitting techniques [44], and error-estimate-based adaptive
integration [45]. We herein consider the error-informed manual selection
strategy proposed in Ref. [45]. At the coarsest level ($\varrho=1$) in the
octree tessellation we set the integration order to $k_{\rm max}$. We then
decrease the order between two levels in such a way that the degree is zero (a
single integration point) at the finest octree levels $\varrho\geq\varrho_{\rm
max}$.

##  5 Immersed isogeometric analysis simulations

In this section we consider three applications of the topology-preserving
image-based immersed isogeometric analysis technique presented above. In
Section 5.1 we start with the case of a single-field problem in two dimensions
by considering an elasticity problem. Subsequently, in Section 5.2 we consider
a multi-field problem in the form of a Stokes flow through a carotid artery
geometry. For this flow case we first consider a representative two-
dimensional test case. The third application pertains to the extension of the
Stokes flow case to a three-dimensional patient-specific geometry based on
scan data. Let us note in advance that the carotid artery test case
realistically pertains to moderate Reynolds number flows [113], but that a
Stokes flow setting is here considered to focus on the topology-preserving
analysis scheme developed in this work.

###  5.1 Uniaxial extension of a two-dimensional structure: a linear
elasticity problem

We consider the two-dimensional specimen shown in Figure 19, which is
represented by $32\times 32$ grayscale voxels. The physical domain, $\Omega$,
with boundary $\partial\Omega$, is immersed into an ambient domain
$\mathcal{A}$ of (dimensionless) size $L\times L$ (with $L=1$), with boundary
$\partial\mathcal{A}$; see Figure 19(c) (with ${\varrho_{\rm max}}=3$). We
consider a linear elasticity problem for which the displacement field,
$\boldsymbol{u}$, is prescribed on the exterior (top and bottom) boundary,
while the interior (immersed) boundary is traction free. In the absence of
inertia effects and body forces, the boundary value problem reads as:

$\left\\{\begin{aligned} \text{Find $\boldsymbol{u}$ such that:}\\\ {\rm
div}(\,\boldsymbol{\sigma}(\boldsymbol{u})\,)&=\boldsymbol{0}&&\mbox{in}\>\Omega\\\
\boldsymbol{u}&=\boldsymbol{0}&&\mbox{on}\>\partial\mathcal{A}_{0}\\\
\boldsymbol{u}&=\bar{u}\boldsymbol{n}&&\mbox{on}\>\partial\mathcal{A}_{\bar{u}}\\\
\boldsymbol{u}\cdot\boldsymbol{n}&=0&&\mbox{on}\>\partial\mathcal{A}\setminus(\partial\mathcal{A}_{0}\cup\partial\mathcal{A}_{\bar{u}})\\\
\left[\boldsymbol{I}-\boldsymbol{n}\otimes\boldsymbol{n}\right]\boldsymbol{\sigma}\boldsymbol{n}&=\boldsymbol{0}&&\mbox{on}\>\partial\mathcal{A}\setminus(\partial\mathcal{A}_{0}\cup\partial\mathcal{A}_{\bar{u}})\\\
\boldsymbol{\sigma}\boldsymbol{n}&=\boldsymbol{0}&&\mbox{on}\>\partial\Omega\setminus\partial\mathcal{A}\end{aligned}\right.$
(41)

The various boundaries are specified in Figure 19(b). The top boundary is
displaced in normal direction by $\bar{u}=0.2$ (20%). In the above problem
definition, the stress is related to the strain by Hooke’s law, _i.e._ ,
$\boldsymbol{\sigma}(\boldsymbol{u})=\lambda{\rm
div}(\boldsymbol{u})\boldsymbol{I}+2\mu\nabla^{s}\boldsymbol{u}$, where
$\nabla^{s}$ denotes the symmetric gradient operator. Throughout this section,
the non-dimensionalized Lamé parameters are set to $\lambda=\frac{1}{2}$ and
$\mu=\frac{1}{2}$. In our analyses, as quantities of interest we consider the
stress state in the specimen and the effective elastic modulus

$\mathcal{Q}=\frac{L}{\bar{u}}\cfrac{1}{V_{\rm
img}}\int\limits_{\Omega}\sigma_{22}\,{\rm d}{V}.$ (42)

In our simulations, the Dirichlet conditions on the external boundary are
applied strongly, _i.e._ , by constraining the degrees of freedom related to
the boundary displacements. This is enabled by the fact that the top and
bottom boundaries are mesh conforming. The Galerkin problem corresponding to
(41) then follows as

$\left\\{\begin{aligned} &\text{Find $\boldsymbol{u}^{h}\in W^{h}(\Omega)$
such that for all $\boldsymbol{v}^{h}\in W^{h}_{0}(\Omega)$:}\\\
&\int_{\Omega}\nabla^{s}\boldsymbol{v}^{h}:\boldsymbol{\sigma}(\boldsymbol{u}^{h})\,{\rm
d}V=0\end{aligned}\right.$ (43)

with the discrete spaces being subsets of $H^{1}(\Omega)$ satisfying the
Dirichlet boundary conditions. The spaces are discretized using second-order
($k=2$) B-spline basis functions defined on a background mesh with uniform
element size $h$. For this setting, immersed analysis results can be obtained
without additional stabilization terms.

In Figure 20 we present the results using a mesh size of $h=L/32$, which is
equal to the voxel size. As can be seen in Figure 20, without application of
the topology-correction algorithm (first column), the left connection (marked
in green) is not reconstructed by the image segmentation procedure. As a
result, the left side of the specimen carries only a small portion of the
load, in the sense that (the vertical component of) the stress is equal to
zero in the left part connected to the top boundary, and relatively small in
the left section that is connected to the bottom boundary. When the topology-
preservation algorithm is applied (second column in Figure 20), the left part
of the structure remains connected, and the left side of the specimen is
appropriately loaded.

From Figure 20 it is observed that the topological anomaly drastically affects
the simulation results. When considering the effective elastic modulus (42),
as shown in Figure 21 for various mesh sizes, it is observed that this
quantity of interest shows fundamentally different behavior between the two
considered cases. Since the topological anomaly occurs independently of the
background-mesh element size, mesh refinement for the determination of the
approximate solution does not repair this problem. Both solutions with and
without the topological anomaly converge under mesh refinement, but the
problem with the anomaly converges to an erroneous result on account of the
incorrect geometry representation.

(a)

(b)

(c)

Figure 19: Illustration of _(a)_ the original grayscale image,
$g(\boldsymbol{x})$, _(b)_ the voxel-segmentation of the image, i.e,
$g(\boldsymbol{x})>0$, and _(c)_ the computational domain (using
${\varrho_{\rm max}}=3$) with the boundary conditions.

(a)

(b)

(c) $\sigma_{22}$

(d) $\sigma_{22}$

Figure 20: Comparison of the segmented geometry obtained from the grayscale
image without (_a_) and with (_b_) topology preservation. The vertical
(dimensionless) stress component for the two cases, computed using $h=L/32$,
is shown in the panels _(c)_ and _(d)_. Figure 21: The effective elasticity
modulus (42) computed with and without application of the topology-
preservation technique for different mesh sizes $h$.

###  5.2 Flow through a carotid artery: a Stokes flow problem

We now consider Stokes flow through a domain $\Omega\subset\mathbb{R}^{n_{d}}$
($n_{d}=2,3$), representative of a carotid artery. This domain is constructed
using the topology-preserving segmentation procedure presented in Section 3,
and is immersed in an ambient domain $\mathcal{A}$. We consider a presssure-
driven incompressible flow of a Newtonian fluid, with viscosity $\mu$, through
the carotid artery. The fluid velocity, $\boldsymbol{u}$, and pressure, $p$,
satisfy the strong formulation

$\left\\{\begin{aligned} \text{Find $\boldsymbol{u}$ and $p$ such that:}\\\
-\nabla\cdot\left(2\mu\nabla^{s}\boldsymbol{u}\right)+\nabla
p&=\boldsymbol{0}&&\mbox{in}\>\Omega\\\
\nabla\cdot\boldsymbol{u}&=0&&\mbox{in}\>\Omega\\\
\boldsymbol{u}&=\boldsymbol{0}&&\mbox{on}\>\Gamma_{D}=\partial\Omega\setminus\partial\mathcal{A}\\\
\left(2\mu\nabla^{s}\boldsymbol{u}-p\boldsymbol{I}\right)\boldsymbol{n}&=\boldsymbol{0}&&\mbox{on}\>\partial\mathcal{A}_{0}\\\
\left(2\mu\nabla^{s}\boldsymbol{u}-p\boldsymbol{I}\right)\boldsymbol{n}&=-\bar{p}\boldsymbol{n}&&\mbox{on}\>\partial\mathcal{A}_{\bar{p}}\end{aligned}\right.$
(44)

where $\bar{p}$ denotes the pressure applied at the inflow (bottom) boundary.

The no-slip boundary condition on the immersed boundary $\Gamma_{D}$ is
imposed weakly through Nitsche’s method. Ghost- and skeleton-stabilizations
are used to avoid oscillations in the velocity and pressure fields (see
Section 4). The mixed Galerkin form is given by

$\displaystyle\left\\{\begin{aligned}
\text{Find}\>\boldsymbol{u}^{h}\in\boldsymbol{V}^{h}\>\text{and}\>p^{h}\in
Q^{h}\>\text{such that:}\>\\\
a(\boldsymbol{u}^{h},\boldsymbol{v}^{h})+b(p^{h},\boldsymbol{v}^{h})+s_{\rm
ghost}(\boldsymbol{u}^{h},\boldsymbol{v}^{h})&=l(\boldsymbol{v}^{h})&\forall\boldsymbol{v}^{h}\in\boldsymbol{V}^{h}\\\
b(q^{h},\boldsymbol{u}^{h})-s_{\rm skeleton}(p^{h},q^{h})&=0&\forall q^{h}\in
Q^{h},\end{aligned}\right.$ (45)

where the bilinear and linear operators are defined as [36]

$\displaystyle a(\boldsymbol{u}^{h},\boldsymbol{v}^{h})$
$\displaystyle:=2\mu(\nabla^{s}\boldsymbol{u}^{h},\nabla^{s}\boldsymbol{v}^{h})-2\mu\left[\langle\nabla^{s}\boldsymbol{u}^{h}\cdot\boldsymbol{n},\boldsymbol{v}^{h}\rangle_{\Gamma_{D}}+\langle\nabla^{s}\boldsymbol{v}^{h}\cdot\boldsymbol{n},\boldsymbol{u}^{h}\rangle_{\Gamma_{D}}\right]+\langle\mu\beta
h^{-1}\boldsymbol{u}^{h},\boldsymbol{v}^{h}\rangle_{\Gamma_{D}}$ (46a)
$\displaystyle b(p^{h},\boldsymbol{v}^{h})$ $\displaystyle:=-(p^{h},{\rm
div}\,\boldsymbol{v}^{h})$ (46b) $\displaystyle l(\boldsymbol{v}^{h})$
$\displaystyle:=-\langle\bar{p},\boldsymbol{v}^{h}\cdot\boldsymbol{n}\rangle_{\partial\mathcal{A}_{\bar{p}}}$
(46c) $\displaystyle s_{\rm skeleton}(p^{h},q^{h})$
$\displaystyle:=\sum_{F\in\mathcal{F}_{\rm
skeleton}}\int_{F}\gamma\mu^{-1}h^{2k+1}\llbracket\partial_{n}^{k}p^{h}\rrbracket\llbracket\partial_{n}^{k}q^{h}\rrbracket\>{\rm
d}S$ (46d) $\displaystyle s_{\rm
ghost}(\boldsymbol{u}^{h},\boldsymbol{v}^{h})$
$\displaystyle:=\sum\limits_{F\in\mathcal{F}_{\rm
ghost}}\int_{F}\tilde{\gamma}\mu
h^{2k-1}\llbracket\partial_{n}^{k}\boldsymbol{u}^{h}\rrbracket\cdot\llbracket\partial_{n}^{k}\boldsymbol{v}^{h}\rrbracket\>{\rm
d}S,$ (46e)

and $(\cdot,\cdot)$ denotes the inner product in $L^{2}(\Omega)$,
$\langle\cdot,\cdot\rangle_{\Gamma_{D}}$ denotes the inner product in
$L^{2}(\Gamma_{D})$, and $\llbracket\cdot\rrbracket$ is the jump operator. The
parameters $\beta$, $\gamma$, and $\tilde{\gamma}$ denote the penalty
constants for the Nitsche term, the Skeleton-stabilization term, and the
Ghost-stabilization term, respectively. We consider second-order ($k=2$)
B-splines constructed on a variety of uniform background meshes with element
size $h$. Let us note that we use equal-order approximations for the velocity
and pressure fields, which is admissible by virtue of the skeleton-
penalization acting on the pressure field.

####  5.2.1 Two-dimensional test case

To demonstrate the developed methodology in the Stokes flow case, we first
consider the idealized two-dimensional geometry shown in Figure 22
(constructed from $32\times 32$ grayscale voxels). In this two-dimensional
case, we consider the ambient domain to be a unit square (with $L=1$), and we
set the non-dimensionalized parameters to $\mu=1$ and $\bar{p}=1$; see Figure
22(c) (constructed with ${\varrho_{\rm max}}=3$). For the penalty parameters
we take, $\beta=100$, $\gamma=0.05$, and $\tilde{\gamma}=0.0005$, which have
been determined empirically.

Figure 23 shows the pressure and velocity contour plots for the case where a
topological anomaly occurs in the form of a pinched-off channel (top row), and
in the case where the topology-preservation algorithm is applied (bottom row).
The topological anomaly evidently obstructs fluid from flowing through the
left branch, resulting in a zero pressure and zero fluid velocity solution in
the top left disconnected domain. The topology-correction strategy proposed in
this work avoids the left channel from being closed and results in a different
flow profile. It is noteworthy that, if Dirichlet conditions were imposed on
the top boundary, the pathological case without the topology-correction would
have become singular.

The influence of the mesh size is studied in Figure 24, which, similar to the
elasticity problem considered above, conveys that both simulation cases
converge under mesh refinement. However, without application of the topology-
correction strategy, the solution converges to an erroneous result.

(a)

(b)

(c)

Figure 22: Illustration of _(a)_ the original grayscale image,
$g(\boldsymbol{x})$, _(b)_ the voxel-segmentation of the image, _i.e_ ,
$g(\boldsymbol{x})>0$, and _(c)_ the computational domain (using
${\varrho_{\rm max}}=3$) with the boundary conditions.

(a) Pressure

(b) Velocity magnitude

(c) Pressure

(d) Velocity magnitude

Figure 23: Comparison of the pressure, $p$, and velocity magnitude,
$|\boldsymbol{u}|$, for the segmented domain constructed without (top row) and
with (bottom row) topology-preservation using $h=L/64$. Figure 24: Total
outflow from the left branch of the carotid artery, computed using different
mesh sizes $h$, with and without the use of the topology-correction strategy.

####  5.2.2 Three-dimensional scan-based simulation

To demonstrate the developed methodology in a real scan-based setting, we
consider a Stokes flow through a carotid artery. The geometry of the carotid
artery is obtained from a CT-scan (see Figure 25). The scan data consists of
$80$ slices of $85\times 70$ voxels. The size of the voxels is $300\times
300\,{\rm\mu m}^{2}$ and the distance between the slices is $400\,{\rm\mu m}$.
Hence, the total size of the scan domain is $25.6\times 21.1\times 32.0\,{\rm
mm}^{3}$. The original grayscale data, represented in DICOM format [113], is
pre-processed using ITK-SNAP (open-source medical image processing tool [114])
from which binary voxel data is exported and read into our Python-based
implementation of the developed topology-correction strategy.

(a)

(b)

(c)

Figure 25: Illustration of _(a)_ the directly segmented scan-data,
$g(\boldsymbol{x})$. The computational domain of the stenotic carotid artery
extracted by the B-spline-based segmentation procedure _(b)_ before and _(c)_
after application of the topology preservation algorithm. The smooth level set
function, $f(\boldsymbol{x})$, is segmented using the midpoint tessellation
procedure with a subdivision level of $\varrho_{\rm max}=2$.

The direct, non-smooth, segmentation of the image data is visualized in Figure
25(a). When applying the spline-based segmentation procedure using second-
order B-splines, the stenotic part of the artery in the original voxel image
is lost (Figure 25(b)). As discussed in Section 2.2, this topological anomaly
resulting from the spline-based segmentation procedure is expected on account
of the feature-to-mesh-size ratio in that particular area.

In Figure 26 a zoom of the stenotic part of the artery is shown. Figure 26(c)
shows the indicator function (21) as determined by the topology-correction
strategy. This image conveys that the topological anomaly in the form of the
missing stenotic part of the artery is appropriately detected. Figure 26(d)
shows the smoothly segmented geometry after THB-spline refinement. From this
figure it is observed that the topological anomaly is corrected by the
proposed strategy. In the case of the complete topology (see Figure 25(c)), it
is observed that additional boundary regions are tagged for refinement on
account of the high-curvature of the boundary surface in these regions.

(a) Original segmentation

(b) Smooth segmentation

(c) Indicator function

(d) Corrected smooth segmentation

Figure 26: Illustration of the topology-preserving procedure focused on the
stenotic part of the carotid artery. Panel _(a)_ shows the original
segmentation obtained directly by thresholding the voxel data, and panel _(b)_
shows the segmentation through the B-spline-based smoothing strategy. A
moving-window technique then locally compares the topology between the two
segmentations, which results in the indicator function _(c)_ that marks
topological differences. THB-spline-based refinements are then introduced to
locally increase the resolution of the smooth level set function, thereby
preserving the topology of the original scan-data _(d)_.

To simulate the flow through the stenotic artery we consider a viscosity of
$\mu=4\,{\rm mPa\,s}$ and a pressure drop of $17.3\,{\rm kPa}$ (130 mm of Hg).
These parameters are selected based on Ref. [113]. The penalty parameters
associated with the weak formulation (45) are set to $\beta=100$,
$\gamma=0.05$, and $\tilde{\gamma}=0.0005$, which have been determined
empirically to yield a stable formulation without adversely affecting the
accuracy of the approximation.

The results for the velocity and pressure fields computed on a uniform mesh
with $h=1.75$ mm in the directions perpendicular to the pressure gradient and
$h=2$ mm in the direction of the pressure gradient are shown in Figure 27.
Similar to the two-dimensional case, the topological anomaly evidently
obstructs fluid from flowing through the stenotic part of the artery,
resulting in a zero pressure and zero fluid velocity solution in the right
artery (see Figure 27). The topology-correction strategy avoids the stenotic
part from being closed and results in a meaningful flow profile. Although the
mesh considered here is relatively coarse, the flux through the stenotic
region of approximately $200\,{\rm mm^{3}/s}$ (corresponding to an average
velocity of approximately $1.3$ m/s, see Figure 27(f)) corresponds reasonably
well with the flow through a circular tube (Hagen–Poiseuille flow) of length
$8$ mm and diameter $0.35$ mm subject to the above-mentioned pressure drop,
which corroborates that the computed speed is meaningful.

In Figure 28 the outflow through the stenotic branch of the artery is depicted
for various uniform meshes, ranging from a very coarse mesh with 8,760
(active) degrees of freedom to a refined mesh with 96,196 degrees of freedom.
Similar to the problems studied above, this mesh verification study conveys
that both simulation cases converge under mesh refinement, but that an
erroneous result is obtained in the case that the topology is not corrected.
Note that, due to the employed uniform meshes, the number of degrees of
freedom increases rapidly under mesh refinement, which will limit the size of
the domain that can be considered by this type of analysis in practice. It is
important to note, however, that mesh refinements throughout most of the
domain do not substantially improve the accuracy of the simulation (in
particular for the considered quantity of interest). Therefore, to properly
leverage the property of immersed techniques that the mesh resolution can be
controlled independently of the geometry (and topology) representation,
significant improvements in computational efficiency can be obtained by means
of adaptive local mesh refinement. The combination of the proposed technique
with an error-estimation-and-adaptivity strategy is an important topic of
further study.

(a) Pressure

(b) Velocity magnitude

(c) Velocity magnitude

(d) Pressure

(e) Velocity magnitude

(f) Velocity magnitude

Figure 27: Comparison of the pressure, $p$, and velocity magnitude,
$|\boldsymbol{u}|$, for the segmented domain constructed without (top row) and
with (bottom row) topology-preservation using a mesh size of $h=2$ mm in the
vertical direction, and $h=1.75$ mm in the other directions. Figure 28: Total
outflow from the right branch of the carotid artery, computed using different
degrees of freedom (DOFs), with and without the use of the topology-correction
strategy.

##  6 Concluding remarks

To leverage the advantageous properties of isogeometric analysis in a scan-
based setting, a smooth representation of the computational domain must be
obtained. This can be achieved by applying a smoothing operation on the voxel-
based gray-scale data and subsequently applying an octree-based tessellation
procedure. A negative side-effect of this smoothing procedure is that it can
induce topological changes when the scan data contains features with a
characteristic length scale similar to the voxel size.

Based on a Fourier analysis of the B-spline-based smoothing operation, it is
proposed to repair smoothing-induced topological anomalies by locally refining
the smoothed gray-scale function using THB-splines. In combination with a
moving-window strategy to detect topological changes, the local refinement
technique is used to develop a topology-preserving image segmentation
technique. Based on a comparison of the Euler characteristic between the
window view on the original voxel data and that on the smoothed
representation, the proposed technique systematically distinguishes shape
changes from topological changes. The algorithm is fail-safe in that it
detects and repairs topological changes, and does not essentially change the
geometry in the (rare) case that a shape change is accidentally marked for
refinement. The proposed location-based masking strategy to detect shape
changes is very effective for the considered test cases, but it is envisioned
that further robustness improvements can be made by the development of a more
advanced masking procedure.

The developed algorithm works for two- and three-dimensional scan data.
Numerical simulations demonstrate the effectivity of the algorithm in both
settings. For all considered test cases, a topologically consistent smoothed
image is obtained after a single topology-correction step. Based on the
presented Fourier analysis this is to be expected, as refining the mesh for
the B-spline grayscale function has a strong effect on the filtering
properties. It can, in principle, occur that topologial changes are not
repaired after a single correction step. Although not considered in this work,
the presented algorithm has the potential to be extended so that it can be
applied recursively in such scenarios.

In this work we have restricted ourselves to immersed isogeometric analyses
based on uniform meshes. In order to optimally benefit from the fact that the
computaional mesh is decoupled from the segmented geometry in the immersed
setting, use should be made of (adaptive) local refinements for the analysis
mesh. Combination of the proposed topology-preservation technique with an
adaptive meshing strategy is therefore an important topic of further study.

## Acknowledgement

We acknowledge the support from the European Commission EACEA Agency,
Framework Partnership Agreement 2013-0043 Erasmus Mundus Action 1b, as a part
of the EM Joint Doctorate Simulation in Engineering and Entrepreneurship
Development (SEED). A.R. acknowledges the partial support of the MIUR-PRIN
project XFAST-SIMS (no. 20173C478N). All the simulations in this work were
performed based on the open source software package Nutils (www.nutils.org)
[111]. We acknowledge the support of the Nutils team. We acknowledge Michele
Conti for providing the CT-scan data of the stenotic carotid and Alice
Finotello for the help in understanding the data.

## References

  * [1] J. A. Cottrell, T. J. R. Hughes, Y. Bazilevs, Isogeometric analysis: Toward integration of CAD and FEA, John Wiley & Sons, 2009.
  * [2] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (39-41) (2005) 4135–4195.
  * [3] Y. Zhang, Y. Bazilevs, S. Goswami, C. L. Bajaj, T. J. R. Hughes, Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow, Computer Methods in Applied Mechanics and Engineering 196 (29-30) (2007) 2943–2959.
  * [4] S. Morganti, F. Auricchio, D. Benson, F. Gambarin, S. Hartmann, T. J. R. Hughes, A. Reali, Patient-specific isogeometric structural analysis of aortic valve closure, Computer Methods in Applied Mechanics and Engineering 284 (2015) 508–520.
  * [5] T. Hoang, C. V. Verhoosel, F. Auricchio, E. H. van Brummelen, A. Reali, Skeleton-stabilized isogeometric analysis: High-regularity interior-penalty methods for incompressible viscous flow problems, Computer Methods in Applied Mechanics and Engineering 337 (2018) 324–351.
  * [6] S. Marschke, L. Wunderlich, W. Ring, K. Achterhold, F. Pfeiffer, An approach to construct a three-dimensional isogeometric model from $\mu$-CT scan data with an application to the bridge of a violin, Computer Aided Geometric Design 78 (2020) 101815.
  * [7] B. van Rietbergen, H. Weinans, R. Huiskes, B. Polman, Computational strategies for iterative solutions of large FEM applications employing voxel data, International Journal for Numerical Methods in Engineering 39 (16) (1996) 2743–2767.
  * [8] D. Rajon, W. E. Bolch, Marching cube algorithm: review and trilinear interpolation adaptation for image-based dosimetric models, Computerized Medical Imaging and Graphics 27 (5) (2003) 411–435.
  * [9] Z. Wang, J. C. M. Teo, C. K. Chui, S. H. Ong, C. H. Yan, S. Wang, H. K. Wong, S. H. Teoh, Computational biomechanical modelling of the lumbar spine using marching-cubes surface smoothened finite element voxel meshing, Computer Methods and Programs in Biomedicine 80 (1) (2005) 25–35.
  * [10] Y. Bazilevs, V. M. Calo, Y. Zhang, T. J. R. Hughes, Isogeometric fluid–structure interaction analysis with applications to arterial blood flow, Computational Mechanics 38 (4) (2006) 310–322.
  * [11] Y. Bazilevs, J. Gohean, T. J. R. Hughes, R. Moser, Y. Zhang, Patient-specific isogeometric fluid–structure interaction analysis of thoracic aortic blood flow due to implantation of the Jarvik 2000 left ventricular assist device, Computer Methods in Applied Mechanics and Engineering 198 (45-46) (2009) 3534–3550.
  * [12] Y. Yu, Y. J. Zhang, K. Takizawa, T. E. Tezduyar, T. Sasaki, Anatomically realistic lumen motion representation in patient-specific space–time isogeometric flow analysis of coronary arteries with time-dependent medical-image data, Computational Mechanics 65 (2) (2020) 395–404.
  * [13] B. Urick, T. M. Sanders, S. S. Hossain, Y. J. Zhang, T. J. R. Hughes, 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.
  * [14] T. Martin, E. Cohen, R. M. Kirby, Volumetric parameterization and trivariate B-spline fitting using harmonic functions, Computer Aided Geometric Design 26 (6) (2009) 648–664.
  * [15] F. Auricchio, M. Conti, M. Ferraro, S. Morganti, A. Reali, R. Taylor, Innovative and efficient stent flexibility simulations based on isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 295 (2015) 347–361.
  * [16] M. S. Pigazzini, Y. Bazilevs, A. Ellison, H. Kim, Isogeometric analysis for simulation of progressive damage in composite laminates, Journal of Composite Materials 52 (25) (2018) 3471–3489.
  * [17] M. Bucelli, M. Salvador, L. Dede, A. Quarteroni, Multipatch isogeometric analysis for electrophysiology: Simulation in a human heart, Computer Methods in Applied Mechanics and Engineering 376 (2021) 113666.
  * [18] T. Terahara, K. Takizawa, T. E. Tezduyar, Y. Bazilevs, M. C. Hsu, Heart valve isogeometric sequentially-coupled FSI analysis with the space–time topology change method, Computational Mechanics (2020) 1–21.
  * [19] C. V. Verhoosel, G. van Zwieten, B. van Rietbergen, R. de Borst, Image-based goal-oriented adaptive isogeometric analysis with application to the micro-mechanical modeling of trabecular bone, Computer Methods in Applied Mechanics and Engineering 284 (2015) 138–164.
  * [20] E. Rank, M. Ruess, S. Kollmannsberger, D. Schillinger, A. Düster, Geometric modeling, isogeometric analysis and the finite cell method, Computer Methods in Applied Mechanics and Engineering 249 (2012) 104–115.
  * [21] D. Schillinger, L. Dede, M. A. Scott, J. A. Evans, M. J. Borden, E. Rank, T. J. R. Hughes, An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Computer Methods in Applied Mechanics and Engineering 249 (2012) 116–150.
  * [22] M. Ruess, D. Tal, N. Trabelsi, Z. Yosibash, E. Rank, The finite cell method for bone simulations: verification and validation, Biomechanics and Modeling in Mechanobiology 11 (3) (2012) 425–437.
  * [23] M. Ruess, D. Schillinger, Y. Bazilevs, V. Varduhn, E. Rank, Weakly enforced essential boundary conditions for NURBS-embedded and trimmed NURBS geometries on the basis of the finite cell method, International Journal for Numerical Methods in Engineering 95 (10) (2013) 811–846.
  * [24] D. Schillinger, M. Ruess, The finite cell method: A review in the context of higher-order structural analysis of CAD and image-based geometric models, Archives of Computational Methods in Engineering 22 (3) (2015) 391–455.
  * [25] M. C. Hsu, C. Wang, F. Xu, A. J. Herrema, A. Krishnamurthy, Direct immersogeometric fluid flow analysis using B-rep CAD models, Computer Aided Geometric Design 43 (2016) 143–158.
  * [26] A. Düster, J. Parvizian, Z. Yang, E. Rank, The finite cell method for three-dimensional problems of solid mechanics, Computer Methods in Applied Mechanics and Engineering 197 (45-48) (2008) 3768–3782.
  * [27] D. Schillinger, M. Ruess, N. Zander, Y. Bazilevs, A. Düster, E. Rank, Small and large deformation analysis with the p\- and B-spline versions of the finite cell method, Computational Mechanics 50 (4) (2012) 445–478.
  * [28] M. Elhaddad, N. Zander, T. Bog, L. Kudela, S. Kollmannsberger, J. Kirschke, T. Baum, M. Ruess, E. Rank, Multi-level hp-finite cell method for embedded interface problems with application in biomechanics, International Journal for Numerical Methods in Biomedical Engineering 34 (4) (2018) e2951.
  * [29] M. Würkner, S. Duczek, H. Berger, H. Köppe, U. Gabbert, A software platform for the analysis of porous die-cast parts using the finite cell method, in: Analysis and Modelling of Advanced Structures and Smart Systems, Springer, 2018, pp. 327–341.
  * [30] S. Duczek, H. Berger, U. Gabbert, The finite pore method: a new approach to evaluate gas pores in cast parts by combining computed tomography and the finite cell method, International Journal of Cast Metals Research 28 (4) (2015) 221–228.
  * [31] H. J. Kim, Y. D. Seo, S. K. Youn, Isogeometric analysis for trimmed CAD surfaces, Computer Methods in Applied Mechanics and Engineering 198 (37-40) (2009) 2982–2995.
  * [32] H. J. Kim, Y. D. Seo, S. K. Youn, Isogeometric analysis with trimming technique for problems of arbitrary complex topology, Computer Methods in Applied Mechanics and Engineering 199 (45-48) (2010) 2796–2812.
  * [33] R. Schmidt, R. Wüchner, K. U. Bletzinger, Isogeometric analysis of trimmed NURBS geometries, Computer Methods in Applied Mechanics and Engineering 241 (2012) 93–111.
  * [34] F. de Prenter, C. V. Verhoosel, E. H. van Brummelen, J. A. Evans, C. Messe, J. Benzaken, K. Maute, Multigrid solvers for immersed finite element methods and immersed isogeometric analysis, Computational Mechanics 65 (3) (2020) 807–838.
  * [35] A. Düster, H. G. Sehlhorst, E. Rank, Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method, Computational Mechanics 50 (4) (2012) 413–431.
  * [36] T. Hoang, C. V. Verhoosel, C. Z. Qin, F. Auricchio, A. Reali, E. H. van Brummelen, Skeleton-stabilized immersogeometric analysis for incompressible viscous flow problems, Computer Methods in Applied Mechanics and Engineering 344 (2019) 421–450.
  * [37] J. N. Jomo, F. de Prenter, M. Elhaddad, D. D’Angella, C. V. Verhoosel, S. Kollmannsberger, J. S. Kirschke, V. Nübel, E. H. van Brummelen, E. Rank, Robust and parallel scalable iterative solutions for large-scale finite cell analyses, Finite Elements in Analysis and Design 163 (2019) 14–30.
  * [38] M. Carraturo, J. N. Jomo, S. Kollmannsberger, A. Reali, F. Auricchio, E. Rank, Modeling and experimental validation of an immersed thermo-mechanical part-scale analysis for laser powder bed fusion processes, Additive Manufacturing 36 (2020) 101498.
  * [39] M. C. Hsu, D. Kamensky, Y. Bazilevs, M. S. Sacks, T. J. R. Hughes, Fluid–structure interaction analysis of bioprosthetic heart valves: significance of arterial wall deformation, Computational Mechanics 54 (4) (2014) 1055–1071.
  * [40] D. Kamensky, M. C. Hsu, D. Schillinger, J. A. Evans, A. Aggarwal, Y. Bazilevs, M. S. Sacks, T. J. R. Hughes, An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves, Computer Methods in Applied Mechanics and Engineering 284 (2015) 1005–1053.
  * [41] D. Kamensky, M. C. Hsu, Y. Yu, J. A. Evans, M. S. Sacks, T. J. R. Hughes, Immersogeometric cardiovascular fluid–structure interaction analysis with divergence-conforming B-splines, Computer Methods in Applied Mechanics and Engineering 314 (2017) 408–472.
  * [42] L. Kudela, N. Zander, T. Bog, S. Kollmannsberger, E. Rank, Efficient and accurate numerical quadrature for immersed boundary methods, Advanced Modeling and Simulation in Engineering Sciences 2 (1) (2015) 1–22.
  * [43] L. Kudela, N. Zander, S. Kollmannsberger, E. Rank, Smart octrees: Accurately integrating discontinuous functions in 3D, Computer Methods in Applied Mechanics and Engineering 306 (2016) 406–426.
  * [44] M. Joulaian, S. Hubrich, A. Düster, Numerical integration of discontinuities on arbitrary domains based on moment fitting, Computational Mechanics 57 (6) (2016) 979–999.
  * [45] S. C. Divi, C. V. Verhoosel, F. Auricchio, A. Reali, E. H. van Brummelen, Error-estimate-based adaptive integration for immersed isogeometric analysis, Computers & Mathematics with Applications 80 (11) (2020) 2481–2516.
  * [46] J. Nitsche, Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind, in: Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg, Vol. 36, Springer, 1971, pp. 9–15.
  * [47] A. Hansbo, P. Hansbo, An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems, Computer Methods in Applied Mechanics and Engineering 191 (47-48) (2002) 5537–5552.
  * [48] A. Embar, J. Dolbow, I. Harari, Imposing dirichlet boundary conditions with Nitsche’s method and spline-based finite elements, International Journal for Numerical Methods in Engineering 83 (7) (2010) 877–898.
  * [49] T. Boiveau, E. Burman, A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity, IMA Journal of Numerical Analysis 36 (2) (2016) 770–795.
  * [50] E. Burman, P. Hansbo, Fictitious domain finite element methods using cut elements: II. a stabilized Nitsche method, Applied Numerical Mathematics 62 (4) (2012) 328–341.
  * [51] E. Burman, P. Hansbo, Fictitious domain methods using cut elements: III. a stabilized Nitsche method for stokes’ problem, ESAIM: Mathematical Modelling and Numerical Analysis 48 (3) (2014) 859–874.
  * [52] F. de Prenter, C. Lehrenfeld, A. Massing, A note on the stability parameter in Nitsche’s method for unfitted boundary value problems, Computers & Mathematics with Applications 75 (12) (2018) 4322–4336.
  * [53] F. de Prenter, C. V. Verhoosel, G. van Zwieten, E. H. van Brummelen, Condition number analysis and preconditioning of the finite cell method, Computer Methods in Applied Mechanics and Engineering 316 (2017) 297–327.
  * [54] E. Burman, Ghost penalty, Comptes Rendus Mathematique 348 (21-22) (2010) 1217–1220.
  * [55] E. Burman, M. A. Fernández, An unfitted Nitsche method for incompressible fluid–structure interaction using overlapping meshes, Computer Methods in Applied Mechanics and Engineering 279 (2014) 497–514.
  * [56] W. Dettmer, C. Kadapa, D. Perić, A stabilised immersed boundary method on hierarchical B-spline grids, Computer Methods in Applied Mechanics and Engineering 311 (2016) 415–437.
  * [57] S. Badia, F. Verdugo, A. F. Martín, The aggregated unfitted finite element method for elliptic problems, Computer Methods in Applied Mechanics and Engineering 336 (2018) 533–553.
  * [58] F. de Prenter, C. V. Verhoosel, E. H. van Brummelen, Preconditioning immersed isogeometric finite element methods with application to flow problems, Computer Methods in Applied Mechanics and Engineering 348 (2019) 604–631.
  * [59] M. Ruess, D. Schillinger, A. I. Oezcan, E. Rank, Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries, Computer Methods in Applied Mechanics and Engineering 269 (2014) 46–71.
  * [60] A. Düster, E. Rank, B. Szabó, The p-version of the finite element and finite cell methods, Encyclopedia of Computational Mechanics Second Edition (2017) 1–35.
  * [61] D. Schillinger, E. Rank, An unfitted hp-adaptive finite element method based on hierarchical B-splines for interface problems of complex geometry, Computer Methods in Applied Mechanics and Engineering 200 (47-48) (2011) 3358–3380.
  * [62] D. Schillinger, A. Düster, E. Rank, The hp-d-adaptive finite cell method for geometrically nonlinear problems of solid mechanics, International Journal for Numerical Methods in Engineering 89 (9) (2012) 1171–1202.
  * [63] S. Heinze, M. Joulaian, A. Düster, Numerical homogenization of hybrid metal foams using the finite cell method, Computers & Mathematics with Applications 70 (7) (2015) 1501–1517.
  * [64] W. Garhuom, S. Hubrich, L. Radtke, A. Düster, A remeshing strategy for large deformations in the finite cell method, Computers & Mathematics with Applications 80 (11) (2020) 2379–2398.
  * [65] J. F. Mangin, V. Frouin, I. Bloch, J. Régis, J. López-Krahe, From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations, Journal of Mathematical Imaging and Vision 5 (4) (1995) 297–318.
  * [66] A. M. Dale, B. Fischl, M. I. Sereno, Cortical surface-based analysis: I. segmentation and surface reconstruction, Neuroimage 9 (2) (1999) 179–194.
  * [67] C. Xu, D. L. Pham, M. E. Rettmann, D. N. Yu, J. L. Prince, Reconstruction of the human cerebral cortex from magnetic resonance images, IEEE Transactions on Medical Imaging 18 (6) (1999) 467–480.
  * [68] B. Fischl, A. Liu, A. M. Dale, Automated manifold surgery: constructing geometrically accurate and topologically correct models of the human cerebral cortex, IEEE Transactions on Medical Imaging 20 (1) (2001) 70–80.
  * [69] F. Ségonne, J. Pacheco, B. Fischl, Geometrically accurate topology-correction of cortical surfaces using nonseparating loops, IEEE Transactions on Medical Imaging 26 (4) (2007) 518–529.
  * [70] H. Lester, S. R. Arridge, A survey of hierarchical non-linear medical image registration, Pattern Recognition 32 (1) (1999) 129–149.
  * [71] Z. Xie, G. E. Farin, Image registration using hierarchical B-splines, IEEE Transactions on Visualization and Computer Graphics 10 (1) (2004) 85–94.
  * [72] A. Pawar, Y. J. Zhang, C. Anitescu, T. Rabczuk, Joint image segmentation and registration based on a dynamic level set approach using truncated hierarchical B-splines, Computers & Mathematics with Applications 78 (10) (2019) 3250–3267.
  * [73] G. Falqui, C. Reina, BRS cohomology and topological anomalies, Communications in Mathematical Physics 102 (3) (1985) 503–515.
  * [74] T. Qaiser, K. Sirinukunwattana, K. Nakane, Y. W. Tsang, D. Epstein, N. Rajpoot, Persistent homology for fast tumor segmentation in whole slide histology images, Procedia Computer Science 90 (2016) 119–124.
  * [75] W. Bae, J. Yoo, J. Chul Ye, Beyond deep residual learning for image restoration: Persistent homology-guided manifold simplification, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops, 2017, pp. 145–153.
  * [76] R. Assaf, A. Goupil, V. Vrabie, T. Boudier, M. Kacim, Persistent homology for object segmentation in multidimensional grayscale images, Pattern Recognition Letters 112 (2018) 277–284.
  * [77] X. Hu, L. Fuxin, D. Samaras, C. Chen, Topology-preserving deep image segmentation, arXiv preprint arXiv:1906.05404.
  * [78] J. Clough, N. Byrne, I. Oksuz, V. A. Zimmer, J. A. Schnabel, A. King, A topological loss function for deep-learning based image segmentation using persistent homology, IEEE Transactions on Pattern Analysis and Machine Intelligence.
  * [79] P. de Dumast, H. Kebiri, C. Atat, V. Dunet, M. Koob, M. B. Cuadra, Segmentation of the cortical plate in fetal brain MRI with a topological loss, arXiv preprint arXiv:2010.12391.
  * [80] N. Byrne, J. R. Clough, G. Montana, A. P. King, A persistent homology-based topological loss function for multi-class CNN segmentation of cardiac mri, arXiv preprint arXiv:2008.09585.
  * [81] A. A. Novotny, R. A. Feijóo, E. Taroco, C. Padra, Topological sensitivity analysis, Computer Methods in Applied Mechanics and Engineering 192 (7-8) (2003) 803–829.
  * [82] M. Burger, B. Hackl, W. Ring, Incorporating topological derivatives into level set methods, Journal of Computational Physics 194 (1) (2004) 344–362.
  * [83] L. He, C. Y. Kao, S. Osher, Incorporating topological derivatives into shape derivatives based level set methods, Journal of Computational Physics 225 (1) (2007) 891–909.
  * [84] I. Larrabide, R. Feijóo, A. Novotny, E. Taroco, Topological derivative: a tool for image processing, Computers & Structures 86 (13-14) (2008) 1386–1403.
  * [85] M. Hintermüller, A. Laurain, Multiphase image segmentation and modulation recovery based on shape and topological sensitivity, Journal of Mathematical Imaging and Vision 35 (1) (2009) 1–22.
  * [86] M. Droske, B. Meyer, M. Rumpf, C. Schaller, An adaptive level set method for medical image segmentation, in: Biennial International Conference on Information Processing in Medical Imaging, Springer, 2001, pp. 416–422.
  * [87] V. Sochnikov, S. Efrima, Level set calculations of the evolution of boundaries on a dynamically adaptive grid, International Journal for Numerical Methods in Engineering 56 (13) (2003) 1913–1929.
  * [88] M. Xu, P. M. Thompson, A. W. Toga, An adaptive level set segmentation on a triangulated mesh, IEEE Transactions on Medical Imaging 23 (2) (2004) 191–201.
  * [89] G. Varadhan, S. Krishnan, T. Sriram, D. Manocha, Topology preserving surface extraction using adaptive subdivision, in: Proceedings of the 2004 Eurographics/ACM SIGGRAPH Symposium on Geometry Processing, 2004, pp. 235–244.
  * [90] Y. Bai, X. Han, J. L. Prince, Octree-based topology-preserving isosurface simplification, in: 2006 Conference on Computer Vision and Pattern Recognition Workshop (CVPRW’06), IEEE, 2006, pp. 81–81.
  * [91] Y. Bai, X. Han, J. L. Prince, Octree grid topology preserving geometric deformable model for three-dimensional medical image segmentation, in: Biennial International Conference on Information Processing in Medical Imaging, Springer, 2007, pp. 556–568.
  * [92] C. Giannelli, B. Jüttler, H. Speleers, THB-splines: The truncated basis for hierarchical splines, Computer Aided Geometric Design 29 (7) (2012) 485–498.
  * [93] D. D’Angella, S. Kollmannsberger, E. Rank, A. Reali, Multi-level Bézier extraction for hierarchical local refinement of isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 328 (2018) 147–174.
  * [94] D. D’Angella, A. Reali, Efficient extraction of hierarchical B-splines for local refinement and coarsening of isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 367 (2020) 113131.
  * [95] X. Han, C. Xu, J. L. Prince, A topology preserving level set method for geometric deformable models, IEEE Transactions on Pattern Analysis and Machine Intelligence 25 (6) (2003) 755–768.
  * [96] H. Lee, H. Kang, M. K. Chung, B. N. Kim, D. S. Lee, Persistent brain network homology from the perspective of dendrogram, IEEE Transactions on Medical Imaging 31 (12) (2012) 2267–2277.
  * [97] D. Lee, J. Yoo, J. C. Ye, Deep residual learning for compressed sensing MRI, in: 2017 IEEE 14th International Symposium on Biomedical Imaging (ISBI 2017), IEEE, 2017, pp. 15–18.
  * [98] L. Piegl, W. Tiller, Curve and surface constructions using rational B-splines, Computer-Aided Design 19 (9) (1987) 485–498.
  * [99] L. Shapiro, Computer vision and image processing, Academic Press, 1992.
  * [100] L. Ortiz-Gracia, C. W. Oosterlee, A highly efficient shannon wavelet inverse fourier technique for pricing european options, SIAM Journal on Scientific Computing 38 (1) (2016) B118–B143.
  * [101] J. Long, E. Shelhamer, T. Darrell, Fully convolutional networks for semantic segmentation, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2015, pp. 3431–3440.
  * [102] M. H. Wilkinson, T. Wijbenga, G. De Vries, M. A. Westenberg, Blood vessel segmentation using moving-window robust automatic threshold selection, in: Proceedings 2003 International Conference on Image Processing (Cat. No. 03CH37429), Vol. 2, IEEE, 2003, pp. II–1093.
  * [103] T. Häikiö, R. Bertram, J. Hyönä, P. Niemi, Development of the letter identity span in reading: Evidence from the eye movement moving window paradigm, Journal of Experimental Child Psychology 102 (2) (2009) 167–181.
  * [104] J. Ohser, W. Nagel, K. Schladitz, The Euler number of discretized sets - on the choice of adjacency in homogeneous lattices, in: Morphology of Condensed Matter, Springer, 2002, pp. 275–298.
  * [105] S. Lichtenbaum, The Weil-étale topology on schemes over finite fields, Compositio Mathematica 141 (3) (2005) 689–702.
  * [106] T. Leinster, The Euler characteristic of a category, arXiv preprint math/0610260.
  * [107] C. Berger, T. Leinster, The Euler characteristic of a category as the sum of a divergent series, Homology, Homotopy and Applications 10 (1) (2008) 41–51.
  * [108] T. M. Fiore, W. Lück, R. Sauer, Finiteness obstructions and Euler characteristics of categories, Advances in Mathematics 226 (3) (2011) 2371–2469.
  * [109] S. van der Walt, J. L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J. D. Warner, N. Yager, E. Gouillart, T. Yu, Scikit-image: image processing in Python, PeerJ 2 (2014) e453.  doi:10.7717/peerj.453.
URL https://doi.org/10.7717/peerj.453

  * [110] E. H. van Brummelen, T. Demont, G. van Zwieten, An adaptive isogeometric analysis approach to elasto-capillary fluid-solid interaction, International Journal for Numerical Methods in Engineering.
  * [111] G. van Zwieten, J. van Zwieten, C. V. Verhoosel, E. Fonn, T. van Opstal, W. Hoitinga, Nutils (version 5.0) (2019).
URL http://doi.org/10.5281/zenodo.3243447

  * [112] E. M. Garau, R. Vázquez, Algorithms for the implementation of adaptive isogeometric methods using hierarchical B-splines, Applied Numerical Mathematics 123 (2018) 58–87.
  * [113] M. Conti, C. Long, M. Marconi, R. Berchiolli, Y. Bazilevs, A. Reali, Carotid artery hemodynamics before and after stenting: A patient specific CFD study, Computers & Fluids 141 (2016) 62–74.
  * [114] P. A. Yushkevich, J. Piven, H. Cody Hazlett, R. Gimpel Smith, S. Ho, J. C. Gee, G. Gerig, User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability, Neuroimage 31 (3) (2006) 1116–1128.