# Edge-Selective Feature Weaving for Point Cloud Matching

Rintaro Yanagi1,2, Atsushi Hashimoto2, Shusaku Sone2,3, Naoya Chiba2,4, Jiaxin
Ma2, and Yoshitaka Ushiku2
Hokkaido Univ.1, OMRON SINIC X Corp.2, OMRON3, Waseda Univ.4
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper tackles the problem of accurately matching the points of two 3D
point clouds. Most conventional methods improve their performance by
extracting representative features from each point via deep-learning-based
algorithms. On the other hand, the correspondence calculation between the
extracted features has not been examined in depth, and non-trainable
algorithms (e.g. the Sinkhorn algorithm) are frequently applied. As a result,
the extracted features may be forcibly fitted to a non-trainable algorithm.
Furthermore, the extracted features frequently contain stochastically
unavoidable errors, which degrades the matching accuracy. In this paper,
instead of using a non-trainable algorithm, we propose a differentiable
matching network that can be jointly optimized with the feature extraction
procedure. Our network first constructs graphs with edges connecting the
points of each point cloud and then extracts discriminative edge features by
using two main components: a shared set-encoder and an edge-selective cross-
concatenation. These components enable us to symmetrically consider two point
clouds and to extract discriminative edge features, respectively. By using the
extracted discriminative edge features, our network can accurately calculate
the correspondence between points. Our experimental results show that the
proposed network can significantly improve the performance of point cloud
matching. Our code is available at https://github.com/yanarin/ESFW

##  1 Introduction

Matching of two sample elements is a fundamental research topic in the field
of computer vision [17, 45, 10]. Among various matching topics, 3D point cloud
registration has been an essential, challenging task with the recent spread of
3D sensing technologies and 3D point cloud applications [45, 10] (e.g.,
autonomous driving [2] and immersive communications [6]). The goal of most 3D
point cloud registration methods is to find the correspondences between the
points of two point clouds. Accurate 3D point cloud registration enables
further development in downstream tasks such as motion transfer [31], shape
editing [22], and object localization for industrial robots [4].

Figure 1: Concept of the proposed network for 3D point cloud registration. Our
proposed network can be jointly optimized with the feature extractor network
to achieve robustness against errors in the extracted features and to maximize
the performance of the feature extractor.

The main process of 3D point cloud registration can be roughly divided into
two steps: feature extraction from each point and correspondence matrix
construction. Most traditional methods [35, 41, 3, 44] focus mainly on how to
extract representative features from each point by using non-trainable
algorithms such as keypoint detectors and local geometric descriptors. With
the recent increase in data availability, deep neural networks (DNNs) have
been widely adopted for feature extraction [27, 36], and they improve the
matching performance significantly. Specifically, an unsupervised 3D point
cloud registration method called CorrNet3D was proposed to reduce the cost of
collecting annotated data, and it has achieved state-of-the-art performance
[44]. Although the rough matching of two point clouds has become possible with
the DNN-based registration methods focusing on feature extraction, the finer
matching is still a challenge because the correspondence matrix is adversely
affected by stochastically unavoidable errors.

When the point clouds do not include any missing points (i.e. occlusions), it
is desirable in 3D point cloud registration for the correspondence matrix to
be binary and orthogonal [35, 41, 3, 44]. The DNN-based registration methods
generally calculate the pairwise similarities between the extracted features
of each point, and the resulting similarity matrix is then converted to the
correspondence matrix via the Sinkhorn algorithm [24], which ideally
constructs a binary, orthogonal matrix. In practice, however, stochastically
unavoidable errors in the extracted features causes errors in the pairwise
similarities [15, 14], which are detrimental to accurate construction of the
correspondence matrix. Moreover, such forcible adjustments can limit the
performance of the feature extractor [26]. Accordingly, for advancing the 3D
point cloud registration methods by maximizing the feature extraction
performance, a trainable matching architecture that can be jointly optimized
with the feature extractor is required.

In this paper, we propose a differentiable matching network that can be
jointly optimized with the feature extractor. Specifically, by adopting a
graph neural network [37], we construct an edge-selective feature weaving
(ESFW) module that can discriminatively analyze the relationships between the
points across the point clouds. This concept is illustrated in Fig. 1. First,
our ESFW module constructs graphs with edges connecting the points of each
point cloud. Then, our ESFW module extracts discriminative edge features by
using two graph neural network-based components: a shared set-encoder and an
edge-selective cross-concatenation. The shared set-encoder enables us to
deeply stack graph neural networks for discriminative extraction of edge
features, while the edge-selective cross-concatenation enables us to handle
mutual point cloud information bidirectionally. Finally, by constructing a
correspondence matrix based on the extracted discriminative edge features, our
network accurately performs the finer matching.

The contributions of this paper are summarized as follows.

Joint optimization of feature extraction and matching

To achieve robustness against errors in the extracted features and to maximize
the performance of the feature extractor, we introduce a differentiable graph
neural network that is specialized for point cloud matching.

Designing network architecture for the finer matching

To achieve the finer matching, we carefully design our ESFW module to
discriminatively analyze the relationships between the points of each point
cloud.

Experimental verification

We experimentally confirm that the performance of the state-of-the-art method
is significantly improved on both rigid and non-rigid datasets in both
supervised and unsupervised settings.

##  2 Related works

###  2.1 Rule-based registration algorithm

For 3D point cloud registration, the pioneering works sought to extract hand-
crafted features and search for the nearest neighbor point [12, 46, 46, 32,
39]. Most such algorithms including Spin Image [12], ISS [46], SHOT [32], RCS
[39], and TOLDI [40], first calculate a local reference frame (LRF). The LRF
enables encoding of the 3D spatial information and makes the hand-crafted
feature extractor rotation invariant. Then, via the LRF, these algorithms
extract the features of each key-point and match the nearest neighbor points.
For accurate matching of the extracted features, the nearest-neighbor distance
ratio (NNDR) [13] is a widely known algorithm. In the NNDR, the nearest and
second-nearest neighbors of each extracted feature are calculated. If the
ratio of the distances between the nearest and second-nearest neighbors is
below a threshold, the nearest neighbor is treated as a correspondence point.
Although the matching performance has been improved with these hand-crafted
matching algorithms, point cloud resolutions, noise, clutter, and occlusions
suppress the performance [7].

Figure 2: Overview of the proposed method. At first, the proposed method
constructs a directional bipartite graph. Then our method extracts the
discriminative edge features based on two main components specialized for
matching: a shared set-encoder and an edge-selective cross-concatenation.

###  2.2 Deep learning-based registration algorithm

The recent growth in deep learning methods has led to various approaches for
application of DNNs to 3D point cloud registration [45, 1, 35, 30, 41, 44].
These approaches can be roughly divided into two types: those that handle the
entire scene simultaneously, and those that are based on local feature-
matching results. While the former approach enables global registration, it is
difficult to apply to large point clouds or scenes with occlusions. Our work
thus follows the latter approach, which can extract local features and perform
end-to-end 3D point cloud registration by introducing a trainable matching
network.

Within the latter approach, many methods extract representative features from
each point and search for corresponding points by calculating the distances
between point features [8, 9]. In the early years of research on this
approach, most methods converted point clouds to other representation styles
such as RGB-D images and then extracted local features with a well-established
2D or 3D convolutional neural network (CNN) [5, 23, 43]. For example, [5]
projected a 3D point cloud to RGB-D images and then extracted features by
using a 2D CNN. Among more recent methods, PointNet [27] and DeepSets [42] can
directly process raw unstructured point clouds. These methods demonstrate that
permutation invariance of the neural network can be introduced by using
symmetric functions, and the resulting feature extractors are used for point
cloud matching. Although these feature extraction methods have greatly
advanced point cloud processing [28, 36, 38], they can only be applied for
matching between individual points [30], which makes it difficult to find the
best combination between points.

To achieve point matching on a many-to-many basis, recent methods construct a
similarity matrix for extracted point feature and then convert the similarity
matrix to a correspondence matrix via the Sinkhorn algorithm [41, 44]. Here,
the Sinkhorn algorithm helps find the best combination between points by
iteratively and alternatively performing the softmax operation column-wise and
row-wise [24]. Among these many-to-many methods, the CorrNet3D model achieves
state-of-the-art performance in both supervised and unsupervised settings by
using an improved version of the Sinkhorn algorithm called the DeSmooth module
[44]. Although these methods can handle point matching on a many-to-many
basis, stochastically unavoidable errors in the extracted features causes
errors in the pairwise similarities, which make it difficult to accurately
construct the correspondence matrix [15, 14]. To overcome this difficulty, and
provide a new idea for many-to-many point matching, our proposed network
learns feature extraction and point matching simultaneously via a
differentiable matching network.

##  3 Edge-selective feature weaving

Let ${\bm{X}^{A}}$ and ${\bm{X}^{B}}$ be point clouds with sizes $N$ and $M$,
respectively. We refer to the $n$-th element in ${\bm{X}^{A}}$ as ${\bf
x}^{A}_{n}\in\mathcal{R}^{{D_{x}}}$ ($n=1,\ldots,N)$ and the $m$-th element in
${\bm{X}^{B}}$ as ${\bf x}^{B}_{m}\in\mathcal{R}^{{D_{x}}}$ ($m=1,\ldots,M)$,
where ${D_{x}}=3$ when we represent a point by its position in 3D space, for
example. Point cloud matching is defined as the task of predicting a
correspondence matrix $\bm{P}=[p_{n,m}]_{{N\times M}}$ from the inputs
${\bm{X}^{A}}$ and ${\bm{X}^{B}}$.

Figure 2 shows an overview of our approach. The goal is to make the
correspondence calculation trainable with a certain network $\mathcal{M}$.
Hence, we assume a feature extractor $\mathcal{E}:{\bf x}\rightarrow{\bf f}$
(${\bf f}\in\mathcal{R}^{D_{f}}$), where the features ${\bf
f}^{A}=\mathcal{E}({\bf x}^{A})$ and ${\bf f}^{B}=\mathcal{E}({\bf x}^{B})$
are fed to $\mathcal{M}$ to estimate $\bm{P}$. By jointly optimizing
$\mathcal{E}$ and $\mathcal{M}$, we seek a better solution than those obtained
with a non-trainable correspondence calculation. In the following subsections,
we describe the details of the preprocessing for $\mathcal{M}$ (Subsec. 3.1)
and the network architecture of $\mathcal{M}$ (Subsec. 3.2).

###  3.1 Preprocessing

For each ${\bf x}^{A}_{n}$, $\mathcal{M}$ works to find a correct match in
${\bm{X}^{B}}$ (and similarly, in ${\bm{X}^{A}}$ for each ${\bf x}^{B}_{m}$).
Let us consider a directional bipartite graph
$G(\bm{V}^{A},\bm{V}^{B},\bm{E}^{A\rightarrow\\!B},\bm{E}^{B\rightarrow\\!A})$
where $\bm{V}^{A}=\\{{\bf f}^{A}_{n}|n=1,\ldots,N\\}$ and $\bm{V}^{B}=\\{{\bf
f}^{B}_{m}|m=1,\dots,M\\}$ are two sides of vertices,
$\bm{E}^{A\rightarrow\\!B}$ is the set of edges from side A to B, and
$\bm{E}^{B\rightarrow\\!A}$ is the set of edges from B to A. For simplicity,
we only give the formulae for A here; the same formulae are defined
symmetrically for B. To prune matching candidates that have no chance to be
matched , we propose to simply select candidate edges from
$\bm{E}^{A\rightarrow\\!B}$ by the $K$-nearest neighbors algorithm.

Our concrete implementation of this calculation is as follows. Since we
implement our method on CorrNet3D, we borrow its similarity definition. For
each pair $\\{{\bf f}^{A}_{n},{\bf f}^{B}_{m}\\}$, we calculate the similarity
$s^{A}_{n,m}$ as

$\displaystyle s^{\text{A}}_{n,m}$ $\displaystyle=\frac{1}{\|{\bf
f}^{A}_{n}-{\bf f}^{B}_{m}\|_{2}}.$ (1)

In terms of this similarity, we find the top $K$ most similar neighbors,
denoted as $\mathcal{N}(n)$, for each $n$-th point and we add
$\\{\\{n,m\\}|m\in\mathcal{N}(n)\\}$ to $\bm{E}^{A\rightarrow\\!B}$.

Finally, using $\bm{E}^{A\rightarrow\\!B}$ and $\bm{E}^{B\rightarrow\\!A}$, we
generate ${{\bf z}^{A}_{0,n,m}}$, an input to the first (0-th) layer of
$\mathcal{M}$ that represents an edge $\\{n,m\\}\in\bm{E}^{A\rightarrow\\!B}$.
Then, with a concatenating operation $\rm{cat}(\cdot)$, we define the feature
as

$\displaystyle{{\bf z}^{A}_{0,n,m}}$ $\displaystyle={\rm cat}({\bf
f}^{A}_{n},(s^{A}_{n,m}),{\bf f}^{B}_{m}).$ (2)

###  3.2 Network architecture

Let ${\bf Z}^{A}_{0}$ denote the set of input features $\\{{{\bf
z}^{A}_{0,n,m}}|\\{n,m\\}\in\bm{E}^{A\rightarrow\\!B}\\}$. We design
$\mathcal{M}$ as a network whose input is $\\{{\bf Z}^{A}_{0},{\bf
Z}^{B}_{0}\\}$ and whose output is an estimate of the correspondence matrix
$\hat{\bm{P}}$. When $\mathcal{M}$ is trained ideally, row-wise and column-
wise argmax operations should yield an identical result, which is the optimal
match.

To implement this matching network architecture, we propose an edge-selective
feature weaving layer inspired from PointNet [27] and DeepSets [42]. Let
$\text{FW}_{l}:\\{{\bf Z}^{A}_{l},{\bf Z}^{B}_{l}\\}\rightarrow\\{{\bf
Z}^{A}_{l+1},{\bf Z}^{B}_{l+1}\\}$ be the $l$-th feature weaving layer. Here,
${\bf Z}^{A}_{l}=\\{{{\bf
z}^{A}_{l,n,m}}|\\{n,m\\}\in\bm{E}^{A\rightarrow\\!B}\\}$ is a set of
intermediate representations when $l>0$. We organize $\mathcal{M}$ by stacking
$L-1$ feature weaving layers ($l=0,\ldots,L-1$) and one output layer.

The feature weaving layer has two main components: a shared set-encoder and an
edge-selective cross-concatenation. The shared set-encoder serves to emphasize
the relative characteristics of each edge among the adjacent edges. The edge-
selective cross-concatenation is applied to the output of the shared set-
encoder to propagate the mutual stream information to the other stream. At the
output layer, we merge the outputs of the shared set-encoder into
$\hat{\bm{P}}$ instead of applying the cross-concatenation.

#### The shared set-encoder

seeks to extract representative edge features for matching. It is desirable
that edge features are highly discriminative for accurate, stable matching of
two elements, and then it is necessary to analyze each edge feature via the
deep network while considering the other edge features. PointNet [27] and
DeepSets [42] provide a suitable structure for this purpose, a combination of
max-pooling with two linear operations. However, these average aggregation
process is also known to smooth edge features as the network gets deeper [19,
25]. Accordingly, to alleviate this effect and extract discriminative features
with the deep network, we extend the conventional graph neural network
process.

Specifically, we first extract edge features via a trainable linear layer
$\phi^{1}_{l}(\cdot)$ and max-pooling $\max(\cdot)$ as follows:

$\displaystyle{{\bf h}^{A}_{l,n}}$
$\displaystyle=\max_{m_{k}\in\mathcal{N}(n)}(\phi^{1}_{l}({{\bf
z}^{A}_{l,n,m_{k}}})).$ (3)

Hence, ${{\bf h}^{A}_{l,n}}$ is expected to aggregate all the specific
characteristics of matching candidates. Next, we concatenate ${{\bf
h}^{A}_{l,n}}$ to the feature of each edge ${{\bf z}^{A}_{l,n,m}}$ and apply
another linear layer $\phi^{2}_{l}(\cdot)$ to compare them as follows:

$\displaystyle{{\bf g}^{A}_{l,n,m}}$
$\displaystyle=\text{act}(\text{BN}(\phi^{2}_{l}(cat({{\bf
z}^{A}_{l,n,m}},{{\bf h}^{A}_{l,n}})))),$ (4)

where $\text{act}(\cdot)$ and $\text{BN}(\cdot)$ represent pReLU activation
and batch normalization, respectively. Here, ${{\bf g}^{A}_{l,n,m}}$ is a
${D_{g}}$-dimensional feature. By feed-forwarding ${{\bf g}^{A}_{l,n,m}}$
without averaging, $\mathcal{M}$ can compare the candidates precisely.

#### The edge-selective cross-concatenation

seeks to propagate the mutual stream information to the other stream. To
accurately match the elements of two samples, it is important to
collaboratively consider valuable edges in the other stream. We achieve this
by concatenating the features in a stream-crossing manner (A to B and B to A).

In principle, the combined features should be relevant to each edge, and then
the concatenation of ${{\bf g}^{A}_{l,n,m}}$ and ${{\bf g}^{B}_{l,m,n}}$ is a
straightforward choice. However, because the similarity matrix
${\bm{S}}=[s_{n,m}]_{{N\times M}}$ is generally not symmetric, the $K$-nearest
selection also yields an asymmetric choice of $\bm{E}^{A\rightarrow\\!B}$ and
$\bm{E}^{B\rightarrow\\!A}$; that is, ${{\bf g}^{B}_{l,m,n}}$ is not always
calculated in stream B. In that case, we have the two options, eliminating
such edges from candidates or preparing a substitution of ${{\bf
g}^{B}_{l,m,n}}$ anyway. Since the former option could eliminate all
candidates in the worst case, we use the latter option, and we specify the
substitute for an absent ${{\bf g}^{B}_{l,m,n}}$ as the mean of the candidates
$\mathcal{N}(n)$. Overall, the edge-selective cross-concatenation is defined
as

$\displaystyle{{\bf z}^{A}_{l+1,n,m}}$ $\displaystyle={\rm cat}({{\bf
g}^{A}_{l,n,m}},{{{\bf g}^{\prime}}^{B}_{l,n,m}}),$ (5) $\displaystyle{{{\bf
g}^{\prime}}^{B}_{l,n,m}}$
$\displaystyle=\left\\{\begin{array}[]{c@{~}l@{}}{{\bf
g}^{B}_{l,m,n}}\hfil~{}&{\rm
if~{}}\\{m,n\\}\in\bm{E}^{B\rightarrow\\!A}\vspace{0.5em}\\\
\frac{1}{K}\mspace{-9.0mu}\sum\limits_{m_{k}\in\mathcal{N}(n)}\mspace{-9.0mu}{{\bf
g}^{B}_{l,m_{k},n}}\hfil~{}&{\rm otherwise}\end{array}\right.,$ (8)

where ${{\bf z}^{A}_{l+1,n,m}}$ is the output of $\text{FW}_{l}$, which is
feed-forwarded to the next feature weaving layer $\text{FW}_{l+1}$.

Note that, to stack the feature weaving layers deeply, we include shortcut
paths every two feature weaving layers to create a residual structure. For
this purpose, we use an identical channel width $D_{z}$ as the output of every
feature weaving layer (i.e., ${{\bf z}^{A}_{l,n,m}}\in\mathcal{R}^{{D_{z}}}$
for $l=1,\ldots,L$), where $D_{z}=2D_{g}$.

#### The output layer

finalizes ${\bf Z}^{A}_{L}$, the output of the ($L-1$) stacked feature weaving
layers, to be an estimate of the correspondence matrix $\hat{\bm{P}}$. This
layer first applies the same calculation as in Eq. (3), and then calculates

$\displaystyle{{\bf g}^{A}_{L,n,m}}$
$\displaystyle=\text{BN}(\phi^{2}_{l}(cat({{\bf z}^{A}_{L,n,m}},{{\bf
h}^{A}_{L,n}}))).$ (9)

This is almost the same as Eq. (4) but without activation. Batch normalization
is necessary to prevent the output from falling into a trivial local minimum,
which would yield quite small values at all points and make $\hat{\bm{P}}$ a
uniform matrix.

Figure 3: Experimental results on the rigid and non-rigid datasets. The
vertical axis indicates the matching performance (Eq. 14), and the horizontal
axis indicates the tolerant errors described in Subsec. 4.1.

Figure 4: The samples of the experimental results with the non-rigid dataset
in unsupervised setting. The green line shows the correctly matching pair, and
the red line shows the incorrectly matching pair. Each sample is randomly
selected for fair comparison, and we show correspondence relationships from
random 600 points for visibility. Note that we select the same points in each
method. The othere results can be seen in the supplemental materials.

Next, we calculate each element $\hat{p}^{A}_{n,m}$ in a matrix
$\hat{\bm{P}}^{A}$ as follows:

$\displaystyle\hat{p}^{A}_{n,m}$
$\displaystyle=\left\\{\begin{array}[]{ll}{{\bf g}^{A}_{L,n,m}}&{\rm
if~{}}\\{n,m\\}\in\bm{E}^{A\rightarrow\\!B}\\\ -\infty&{\rm
otherwise}\end{array}\right.$ (12)

By merging the two streams’ outputs, we obtain the final output as

$\displaystyle\hat{\bm{P}}=\frac{\hat{\bm{P}}^{A}+{{\hat{\bm{P}}}^{B^{\top}}}}{2}.$
(13)

##  4 Experiments

To evaluate the effectiveness of our edge-selective feature weaving (ESFW)
module, we conducted experiments with both supervised and unsupervised point
cloud registration tasks. Because the CorrNet3D model obviously outperforms
the other recent 3D point cloud registration methods on both rigid and non-
rigid datasets in both supervised and unsupervised settings, we adopted it as
a baseline. The performance of other methods can be seen in graphs of its
paper[44].

###  4.1 Experimental settings

Dataset. Following the experiments in [44], we conducted experiments in both
rigid and non-rigid settings. Note that, in each experiment, we used the data
split provided by [44]. For the rigid setting, we used the Surreal dataset
[33], which contains 230K point cloud samples for training and 100 samples for
testing. The 230K point clouds were randomly paired into 115K training
samples, and 100 test pairs were created by randomly rotating and translating
the test samples. For the non-rigid setting, we used the Surreal dataset and
the SHREC dataset [18] for training and test, which respectively contain 230K
and 860 samples. The training samples were the same as in the rigid setting,
while the 860 test samples were randomly combined into 430 test pairs. For all
the above datasets, each point cloud contained 1024 points.

Evaluation metrics. For the evaluation metrics, following [44], we used the
corresponding percentage (Corr) under various tolerant errors. The Corr and
tolerant errors are defined as follows:

Corr $\displaystyle=\frac{1}{N}||\bm{P}\odot\bm{P}^{\text{gt}}||_{1},$ (14)
tolerant errors $\displaystyle=r/{\text{dist}_{\text{max}}},$ (15)
$\displaystyle\text{dist}_{\text{max}}$
$\displaystyle:=\text{max}\\{||\bm{x}^{\text{A}}_{i}-\bm{x}^{\text{A}}_{j}||,\forall{i,j}\\},$
(16)

where $\odot$, $\bm{P}^{\text{gt}}$, $||\cdot||_{1}$ and $r$ respectively
represent the Hadamard product, ground-truth correspondence matrix, L1-norm
and tolerant radius.

Implementation details. In this experiment, we replaced the DeSmooth module in
the CorrNet3D model with our ESFW module, which was optimized according to the
loss function, as in CorrNet3D. For the parameters in the equations, we set
$K=128$, $D_{f}=128$, $L=10$ and $D_{g}=16$. Also, following [44], we trained
the model with the Adam optimizer, a learning rate of 1e-4, a batch size of 10
and 100 epochs. We implemented the model with the PyTorch framework on a Tesla
V100 GPU.

###  4.2 Experimental results

Experimental results on rigid dataset. We show the experimental results on the
rigid dataset in the “Rigid dataset” of Fig. 3. Here, the result labeled
“CorrNet3D+ESFW in supervised (resp. unsupervised)” shows the matching
performance of CorrNet3D with our ESFW module in the supervised (resp.
unsupervised) setting. As shown in the figure, the CorrNet3D+ESFW model
outperforms the baseline CorrNet3D model in both supervised and unsupervised
settings. From these results, we confirmed that our ESFW module could
significantly improve the performance of the CorrNet3D model simply by
including it in the matching procedure. In other words, for the rigid point
cloud matching, the trainable GNN-based function is more effective than the
non-trainable function. Furthermore, the performance improvement especially
for smaller tolerant errors reveals that our ESFW module discriminatively
finds the corresponding point from similar points. From these results, the
effectiveness of our ESFW module is confirmed in the rigid point cloud
matching setting.

Experimental results on non-rigid dataset. Figure 3 also shows the
experimental results on the non-rigid dataset. Again, the CorrNet3D+ESFW model
outperforms the CorrNet3D model in both the supervised and unsupervised
settings, which confirmed the effectiveness of our ESFW module in non-rigid
point cloud matching setting. On the other hand, the performance improvement
for smaller tolerant errors seems not as high as in the rigid case. We assume
that these results are due to the difficulties of the non-rigid point cloud
matching setting. For further performance improvement, investigation about
sophisticated collaboration between the baseline model and ESFW module is
required.

Qualitative results. We show the samples of the experimental results in Fig.
4. In each sample, the CorrNet3D+ESFW model can accurately match points
compared with the CorrNet3D model. Furthermore, the results of the
CorrNet3D+ESFW model in tolerant errors = $0.00$ is similar to the results of
the CorrNet3D model in tolerant errors = $0.06$. From these results, it is
considered that our ESFW module helps finer matching of CorrNet3D model.

Figure 5: Experimental results from varying the hyperparameter $N$ with the
rigid and non-rigid datasets. The vertical axis indicates the matching
performance (Eq. 14), and the horizontal axis indicates the hyperparameter
$N$.

###  4.3 Ablation study

For further understanding of our ESFW module, we extensively conducted various
extensive ablation studies. Specifically, we examined the effects of the
hyperparameters $N$, $K$, $L$, and $D_{g}$. Each experiment follows the same
settings described in Subsec. 4.1, and we report the experimental results of
changing each hyperparameter. Furthermore, we compared the matching
performance with that of the CorrNet3D model in both the supervised and
unsupervised settings. Because the CorrNet3D model does not include $K$, $L$,
and $D_{g}$, these performances are plotted as horizontal lines. For the
experiments with $K$, $L$, and $D_{g}$, we show the experimental results with
the tolerant errors = $0.00$.

Figure 6: Experimental results from varying the hyperparameter $K$ with the
rigid and non-rigid dataset. The vertical axis indicates the matching
performance (Eq. 14), and the horizontal axis indicates the hyperparameter
$K$.

Effect of $N$. The hyperparameter $N$ represents the number of input points,
and its value depends on the point cloud. To confirm the robustness of our
ESFW module with respect to $N$, we conducted experiments with
$N=64,~{}128,~{}256,~{}512$, and $1024$. In these experiments, we subsampled
the point clouds by furthest point sampling.

Figure 5 shows the experimental results. For each $N$, the CorrNet3D+ESFW
model outperforms the CorrNet3D model in both the supervised and unsupervised
settings, which confirmed the robustness of our ESFW module with respect to
$N$. On the other hand, we observed that the tendencies of each experimental
result are differed according to the test dataset. For example, with $N>256$,
although the matching performances are not decreased in the rigid dataset,
they are decreased on the non-rigid dataset. For practical application,
further analysis is required to understand the differences between the
datasets as the next step of our research.

Effects of $K$. The hyperparameter $K$ represents the number of edges coming
from each node. Accordingly, $K$ affects memory usage, and it is desirable to
adopt a smaller $K$. To examine the relationship between the matching
performance and $K$, we conducted experiments with $K=64$, $K=128$, and
$K=256$.

Experimental results are shown in Fig. 6. As shown in Fig. 6, the performance
of the CorrNet3D+ESFW model decreased with $K=64$. These results reveal that
our network is sensitive to the value of $K$, and we regard this as a
limitation of our approach. On the other hand, such sensitivity to a
hyperparameter can be observed universally in the field of point cloud
processing [29, 20, 11], and we believe that this limitation should be tackled
as a general challenge for the field.

Figure 7: Experimental results from varying the hyperparameter $L$ with the
rigid and non-rigid dataset. The vertical axis indicates the matching
performance (Eq. 14), and the horizontal axis indicates the hyperparameter
$L$.

Effects of $L$. The hyperparameter $L$ represents the number of feature
weaving layers in our ESFW module. To confirm that the ESFW module can extract
valuable edge features while alleviating the smoothing problem (as described
in “the shared set-encoder” of SubSec. 3.2), we conducted experiments with
$L=4$, $L=6$, $L=8$, and $L=10$.

Experimental results are shown in Fig. 7. As shown in Fig. 7, the performance
of the CorrNet3D+ESFW model improved in proportion to the number of feature
weaving layers $L$. From these results, we confirmed that the ESFW module can
extract valuable edge features while alleviating the smoothing problem.
Furthermore, even with smaller $L$, the CorrNet3D+ESFW model outperforms the
CorrNet3D model in both supervised and unsupervised settings. From these
results, we further confirmed the effectiveness of our ESFW module towards the
3D point cloud registration.

Figure 8: Experimental results from varying the hyperparameter $D_{g}$ with
the rigid and non-rigid dataset. The vertical axis indicates the matching
performance (Eq. 14), and the horizontal axis indicates the hyperparameter
$D_{g}$.

Effects of $D_{g}$. The hyperparameter $D_{g}$ represents the dimension of the
extracted features in our network, and the number of trainable parameters also
proportionately increases with $D_{g}$. Accordingly, $D_{g}$ directly effects
the memory usage, and it is desirable to adopt a smaller $D_{g}$. To examine
the relationship between the matching performance and $D_{g}$, we conducted
experiments with $D_{g}=8$, $D_{g}=16$, $D_{g}=32$, and $D_{g}=64$.

Experimental results are shown in Fig. 8. As shown in Fig. 8, the
CorrNet3D+ESFW model achieves the best performance with $D_{g}=16$ in both the
supervised and unsupervised settings. From these results, we confirmed that
adjustment of $D_{g}$ leads to further performance improvement. It is
considered that excessive and insufficient parameters respectively cause
overfitting and underfitting in general neural networks. On the other hand,
for each $D_{g}$, the CorrNet3D+ESFW model outperforms the CorrNet3D model in
both supervised and unsupervised settings. From these results, we confirmed
the robustness with respect to $D_{g}$.

##  5 Limitation

The first remarkable limitation of our approach is the sensitivity to the
hyperparameter $K$. As with the other point cloud processing methods [29, 20,
11], our ESFW module is also subject to negative effects in a particular
hyperparameter setting. Specifically, as shown in Fig. 6, it is crucial to
select a suitable value of $K$, which is the number of selected edges. Too
small a $K$ value can eliminate the edge toward the corresponding node,
whereas too large a $K$ value can lead to leave valueless edges. To overcome
this limitation, further architectural improvements are needed.

The second limitation is memory usage as shown in point cloud processing [34,
21] and graph neural networks [37, 16]. Our ESFW module process graphs with
edge information, and its memory usage thus becomes larger in proportion with
the numbers of input points $N$ and selected edges $K$. Because this
limitation is closely related to the first limitation, it is important to
consider an ideal architecture in terms of both limitations. Note that, since
the selection of $K$ edges is conducted based on the extracted features, we
believe that the architecture should be considered not only for our ESFW
module but also for the point cloud feature extractor.

##  6 Conclusion

In this paper, we have proposed a differentiable matching network for point
cloud matching that can be jointly learned with the feature extraction
procedure. Our experimental results demonstrate that our idea of joint
optimization between feature extraction and point-to-point matching expands
the limits of 3D point cloud registration and will cultivate further
performance improvement. In addition to the 3D point cloud registration, the
proposed network may also improve the performance of other matching tasks. In
our future works, we will tackle both architecture improvement and
verification on other tasks.

## References

  * [1] Yasuhiro Aoki, Hunter Goforth, Rangaprasad Arun Srivatsan, and Simon Lucey.  PointNetLK: Robust & efficient point cloud registration using pointnet.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7163–7172, 2019.
  * [2] Yaodong Cui, Ren Chen, Wenbo Chu, Long Chen, Daxin Tian, Ying Li, and Dongpu Cao.  Deep learning for image and point cloud fusion in autonomous driving: A review.  IEEE Transactions on Intelligent Transportation Systems, 2021.
  * [3] Nicolas Donati, Abhishek Sharma, and Maks Ovsjanikov.  Deep geometric functional maps: Robust feature learning for shape correspondence.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8592–8601, 2020.
  * [4] Guoguang Du, Kai Wang, Shiguo Lian, and Kaiyong Zhao.  Vision-based robotic grasping from object localization, object pose estimation to grasp estimation for parallel grippers: a review.  Artificial Intelligence Review, 54(3):1677–1734, 2021.
  * [5] Gil Elbaz, Tamar Avraham, and Anath Fischer.  3D point cloud registration for localization using a deep neural network auto-encoder.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4631–4640, 2017.
  * [6] Guillaume Gamelin, Amine Chellali, Samia Cheikh, Aylen Ricca, Cedric Dumas, and Samir Otmane.  Point-cloud avatars to improve spatial communication in immersive collaborative virtual environments.  Personal and Ubiquitous Computing, pages 1–18, 2020.
  * [7] Yulan Guo, Mohammed Bennamoun, Ferdous Sohel, Min Lu, Jianwei Wan, and Ngai Ming Kwok.  A comprehensive performance evaluation of 3D local feature descriptors.  International Journal of Computer Vision, 116(1):66–89, 2016.
  * [8] Yulan Guo, Hanyun Wang, Qingyong Hu, Hao Liu, Li Liu, and Mohammed Bennamoun.  Deep learning for 3D point clouds: A survey.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2020\.
  * [9] Xian-Feng Hana, Jesse S Jin, Juan Xie, Ming-Jie Wang, and Wei Jiang.  A comprehensive review of 3D point cloud descriptors.  arXiv preprint arXiv:1802.02297, 2, 2018.
  * [10] Xiaoshui Huang, Guofeng Mei, Jian Zhang, and Rana Abbas.  A comprehensive survey on point cloud registration.  arXiv preprint arXiv:2103.02690, 2021.
  * [11] Haiyong Jiang, Feilong Yan, Jianfei Cai, Jianmin Zheng, and Jun Xiao.  End-to-end 3D point cloud instance segmentation without detection.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 12796–12805, 2020.
  * [12] Andrew E Johnson and Martial Hebert.  Using spin images for efficient object recognition in cluttered 3D scenes.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 21(5):433–449, 1999.
  * [13] Pedro R Mendes Júnior, Roberto M De Souza, Rafael de O Werneck, Bernardo V Stein, Daniel V Pazinato, Waldir R de Almeida, Otávio AB Penatti, Ricardo da S Torres, and Anderson Rocha.  Nearest neighbors distance ratio open-set classifier.  Machine Learning, 106(3):359–386, 2017.
  * [14] Zhao Kang, Yiwei Lu, Yuanzhang Su, Changsheng Li, and Zenglin Xu.  Similarity learning via kernel preserving embedding.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 4057–4064, 2019.
  * [15] Zhao Kang, Chong Peng, and Qiang Cheng.  Kernel-driven similarity learning.  Neurocomputing, 267:210–219, 2017.
  * [16] Johannes Klicpera, Aleksandar Bojchevski, and Stephan Günnemann.  Predict then propagate: Graph neural networks meet personalized pagerank.  arXiv preprint arXiv:1810.05997, 2018.
  * [17] Chengcai Leng, Hai Zhang, Bo Li, Guorong Cai, Zhao Pei, and Li He.  Local feature descriptor for image matching: A survey.  IEEE Access, 7:6424–6434, 2018.
  * [18] Bo Li, Yijuan Lu, Chunyuan Li, Afzal Godil, Tobias Schreck, Masaki Aono, Martin Burtscher, Qiang Chen, Nihad Karim Chowdhury, Bin Fang, et al.  A comparison of 3D shape retrieval methods based on a large-scale benchmark supporting multimodal queries.  Computer Vision and Image Understanding, 131:1–27, 2015.
  * [19] Qimai Li, Zhichao Han, and Xiao-Ming Wu.  Deeper insights into graph convolutional networks for semi-supervised learning.  In Proceedings of the AAAI Conference on Artificial Intelligence, pages 3538–3545, 2018.
  * [20] Hao Liu, Yulan Guo, Yanni Ma, Yinjie Lei, and Gongjian Wen.  Semantic context encoding for accurate 3D point cloud segmentation.  IEEE Transactions on Multimedia, 2020.
  * [21] Jinxian Liu, Bingbing Ni, Caiyuan Li, Jiancheng Yang, and Qi Tian.  Dynamic points agglomeration for hierarchical point sets learning.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 7546–7555, 2019.
  * [22] Weiping Liu, Jia Sun, Wanyi Li, Ting Hu, and Peng Wang.  Deep learning on point clouds and its application: A survey.  Sensors, 19(19):4188, 2019.
  * [23] Daniel Maturana and Sebastian Scherer.  Voxnet: A 3D convolutional neural network for real-time object recognition.  In Proceedings of IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 922–928, 2015.
  * [24] Gonzalo Mena, David Belanger, Scott Linderman, and Jasper Snoek.  Learning latent permutations with Gumbel-Sinkhorn networks.  arXiv preprint arXiv:1802.08665, 2018.
  * [25] Kenta Oono and Taiji Suzuki.  Graph neural networks exponentially lose expressive power for node classification.  In Proceedings of the International Conference on Learning Representations, 2020.
  * [26] Pedro O Pinheiro.  Unsupervised domain adaptation with similarity learning.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 8004–8013, 2018.
  * [27] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas.  PointNet: Deep learning on point sets for 3D classification and segmentation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 652–660, 2017.
  * [28] Charles R Qi, Li Yi, Hao Su, and Leonidas J Guibas.  PointNet++: Deep hierarchical feature learning on point sets in a metric space.  arXiv preprint arXiv:1706.02413, 2017.
  * [29] Siamak Ravanbakhsh, Jeff Schneider, and Barnabas Poczos.  Deep learning with sets and point clouds.  arXiv preprint arXiv:1611.04500, 2016.
  * [30] Vinit Sarode, Xueqian Li, Hunter Goforth, Yasuhiro Aoki, Rangaprasad Arun Srivatsan, Simon Lucey, and Howie Choset.  PCRNet: Point cloud registration network using pointnet encoding.  arXiv preprint arXiv:1908.07906, 2019.
  * [31] Yang-Tian Sun, Qian-Cheng Fu, Yue-Ren Jiang, Zitao Liu, Yu-Kun Lai, Hongbo Fu, and Lin Gao.  Human motion transfer with 3D constraints and detail enhancement.  arXiv preprint arXiv:2003.13510, 2020.
  * [32] Federico Tombari, Samuele Salti, and Luigi Di Stefano.  Unique signatures of histograms for local surface description.  In Proceedings of European Conference on Computer Vision, pages 356–369. Springer, 2010.
  * [33] Gül Varol, Javier Romero, Xavier Martin, Naureen Mahmood, Michael J. Black, Ivan Laptev, and Cordelia Schmid.  Learning from synthetic humans.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017.
  * [34] Weiyue Wang, Ronald Yu, Qiangui Huang, and Ulrich Neumann.  SGPN: Similarity group proposal network for 3D point cloud instance segmentation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2569–2578, 2018.
  * [35] Yue Wang and Justin M Solomon.  Deep closest point: Learning representations for point cloud registration.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 3523–3532, 2019.
  * [36] Yue Wang, Yongbin Sun, Ziwei Liu, Sanjay E Sarma, Michael M Bronstein, and Justin M Solomon.  Dynamic graph cnn for learning on point clouds.  ACM Transactions on Graphics, 38(5):1–12, 2019.
  * [37] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  IEEE Transactions on Neural Networks and Learning Systems, 32(1):4–24, 2020.
  * [38] Yifan Xu, Tianqi Fan, Mingye Xu, Long Zeng, and Yu Qiao.  Spidercnn: Deep learning on point sets with parameterized convolutional filters.  In Proceedings of the European Conference on Computer Vision, pages 87–102, 2018.
  * [39] Jiaqi Yang, Qian Zhang, Ke Xian, Yang Xiao, and Zhiguo Cao.  Rotational contour signatures for both real-valued and binary feature representations of 3D local shape.  Computer Vision and Image Understanding, 160:133–147, 2017.
  * [40] Jiaqi Yang, Qian Zhang, Yang Xiao, and Zhiguo Cao.  Toldi: An effective and robust approach for 3D local shape description.  Pattern Recognition, 65:175–187, 2017.
  * [41] Zi Jian Yew and Gim Hee Lee.  RPM-Net: Robust point matching using learned features.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11824–11833, 2020.
  * [42] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexander J Smola.  Deep sets.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [43] Andy Zeng, Shuran Song, Matthias Nießner, Matthew Fisher, Jianxiong Xiao, and Thomas Funkhouser.  3DMatch: Learning local geometric descriptors from rgb-d reconstructions.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1802–1811, 2017.
  * [44] Yiming Zeng, Yue Qian, Zhiyu Zhu, Junhui Hou, Hui Yuan, and Ying He.  CorrNet3D: Unsupervised end-to-end learning of dense correspondence for 3D point clouds.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6052–6061, 2021.
  * [45] Zhiyuan Zhang, Yuchao Dai, and Jiadai Sun.  Deep learning based point cloud registration: an overview.  Virtual Reality & Intelligent Hardware, 2(3):222–246, 2020.
  * [46] Yu Zhong.  Intrinsic shape signatures: A shape descriptor for 3D object recognition.  In Proceedings of IEEE International Conference on Computer Vision Workshops, pages 689–696. IEEE, 2009.


# On growth, buckling, and rupture of aneurysms

Masoud Hejazi  A. Srikantha Phani  Department of Mechanical Engineering, 6250
Applied Science Lane, University of British Columbia, Vancouver, B.C, Canada
V6T 1Z4

###### Abstract

Aneurysms are localized bulges of arteries, and they can rupture with fatal
consequences. Complex mechanobiological factors preclude _in vivo_ testing to
assess the rupture risk of an aneurysm, and size based criteria are often used
in clinical practice to guide surgical interventions. It is often found that
tortuous and buckled aneurysms can exceed in size recommended for surgical
intervention, and yet do not rupture. This study addresses why buckled
aneurysms exhibit this intriguing behaviour by combining _in vitro_ inflation
experiments on hyperelastic tubes with finite element calculations. Using a
biologically relevant material model for an arterial wall, we show that
buckled aneurysms can grow in size without rupturing under favourable arterial
pre-tensions. Stretch reversal phenomenon exhibited by arteries governs
whether buckling or bulging occurs first. Exponential stiffening favours the
axial propagation of an aneurysm instead of radial growth in size. The choice
of failure criteria based on Ogden’s strain energy function, Gent’s first
stretch invariant, and Cauchy stress are discussed. Failure maps incorporating
post bifurcation (bulging and buckling) response are constructed to delineate
the regimes of growth, buckling and rupture of an aneurysm.

###### keywords:

Aneurysm, Buckling, Rupture, Finite Element Analysis

##  1 Introduction

An aneurysm is a balloon-like bulge within the wall of a blood vessel,
commonly observed in the abdominal and the thoracic regions of aorta, due to
local degradation in mechanical properties of an artery [Sakalihasan et al.,
2018]. Biomechanics of aneurysms has been reviewed in [Vorp, 2007, Sakalihasan
et al., 2018], arterial wall modelling in [Holzapfel and Ogden, 2010,
Humphrey, 2013, Fung, 2013]. An increasing incidence rate of aneurysms,
ranging from 5.6 to 10.6 per 100,000 persons/year, has been observed in the
past decade in developed countries [Acosta et al., 2006]. Mortality rate for a
ruptured aortic aneurysm can be up to 90%, especially when it occurs outside
of a hospital [Sampson et al., 2014]. Hence, estimating the rupture risk for
AAAs is critical after an aneurysm is detected. Clinical guidelines rely on
the maximum diameter of an aneurysm to estimate the rupture risk [Muluk et
al., 2017, Lazaris, 2019]. However, some AAAs are found unruptured despite
their maximum size exceeding the minimum recommended size (5.5 cm) for
surgical interventions [Lazaris, 2019, Fillinger et al., 2004]. Tortuous and
buckled AAAs whose size far exceed the guidelines have been observed without
rupture [Fillinger et al., 2004, Hejazi et al., 2021, Boyd et al., 2016, Ducas
et al., 2020, Choi et al., 2021]. Arterial stability with aneurysm [Lee et
al., 2014] and buckling without aneurysms [Liu et al., 2019] has also been
studied, among others. However, the rupture of buckled aneurysms remains to be
addressed.

Recently, the hypothesis that buckling offers a fail safe mode for AAAs has
been tested through controlled _in vitro_ experiments on hyperelastic tubes
and numerical modelling [Hejazi et al., 2021] by treating the walls as
isotropic and hyperelastic membranes. Pre-tension has been identified as the
key factor that governs whether buckling occurs or not. It has been reported
that, under low pre-tensions (arterial pre-stress) an aneurysm buckles and
grows in size without rupture. Buckling of an aneurysm is shown to arise from
the compressive stresses developed outside the bulged region which overcome
the arterial pre-tension. Hyperelastic models may be justified on the basis of
the degradation of elastic fibres due to aneurysm’s pathogenesis at incipient
rupture [Niestrawska et al., 2019]. Bulge formation and other bifurcation
phenomena in hyperelastic cylindrical tubes have been studied extensively
using approximate membrane theories for thin walled tubes [Haughton and Ogden,
1979b, a, Fu et al., 2012], diffuse interface models [Lestringant and Audoly,
2018, 2020], and for thick walled cylindrical tubes [Fu et al., 2016, Wang and
Fu, 2018, Wang et al., 2019, Ye et al., 2020]. Hyperelastic models have also
been used in finite element calculations [Länne et al., 1992, Gonçalves et
al., 2008] to study bulge formation in cylindrical tubes. While a
phenomenological and qualitative understanding of aneurysm formation and
growth can be obtained using hyperelastic models, they may not be suitable to
describe the complex material response of a fiber-reinforced, multi-layered,
and anisotropic arterial wall [Holzapfel et al., 2019]. Attempts have been
made to use more complex material models in finite element calculations
[Alhayani et al., 2013, Demirkoparan and Merodio, 2017, Font et al., 2021,
Topol et al., 2021, Dehghani et al., 2019]. These studies indicate the
asymmetry of aneurysm shape [Topol et al., 2021, Font et al., 2021] due to
residual stresses, buckling without bulging [Dehghani et al., 2019],
competition between radial expansion and axial propagation of a bulge
[Alhayani et al., 2014], and the effect of fiber reinforcement on the
initiation of the aneurysm [Demirkoparan and Merodio, 2017]. However, these
computational studies do not address the rupture in the post bulging regime;
the role of buckling in preventing it; the influence of arterial pre-tension
on the occurrence of buckling, bulging and rupture. Significantly, features
such as stretch reversal and exponential stiffening are absent in hyperelastic
material models and their influence on aneurysm evolution needs to be
clarified.

This study aims to understand the influence of biologically relevant
constitutive material laws on the biomechanics of growth and rupture of the
abdominal aortic aneurysm (AAA). To do so, we use the Gasser-Ogden-Holzapfel
(GOH) constitutive model of arteries [Gasser et al., 2006] to account for
exponential stiffening and stretch reversal [Sommer et al., 2010].
Experimental protocols and finite element analyses are described in Section 2.
Section 3 discusses results and factors governing the growth and rupture of an
aneurysm, culminating in failure maps. Clinical implications and limitations
are discussed. Cncluding remarks and avenues for further investigation are
provided in Section 4.

##  2 Methodology

We first present the experimental methodology, followed by finite element
calculations. Experiments on hyperelastic tube models will be used to validate
the finite element methodology, which will then be extended to include
biologically relevant material models.

###  2.1 Experiments on hyperelastic cylindrical tubes

Inflation tests on finite latex rubber (hyperelastic) tubes are conducted as
illustrated in Fig. 1. The purpose of these tests is to understand the
deformation phenomena of buckling and kinking, and their influence on the
axial and radial growth of a bulge. Here, the emphasis is on post bulge
regime, and the influence of material models on this regime that has not been
investigated in the literature. Material and geometric properties of the tubes
used are given in Table. 1. A nominally straight cylindrical rubber tube is
fixed at its bottom, and the top end is attached to an Instron (model 5969)
tensile tester’s jaw. A known pre-tension (pre-stretch) is applied by moving
the jaw upward and then holding it fixed at a prescribed displacement
throughout the inflation test. The load cell reads the axial force ($F$) in
this displacement control test. The air inlet is attached to the top end of
the tube and a pressure gauge reads the injection pressure. A DSLR camera
(Nikon D7000) records the deformation of the tube during air injection and a
pressure transducer records the internal pressure. An edge detection algorithm
in MATLAB© image processing toolbox is used to obtain the bulged profile for
further analysis. A number of experiments (39 in total) are performed for
different initial pre-tensions, $F_{i}$, and aspect ratios by changing the
initial lengths of the tube. The initial diameter and the wall thickness are
held fixed in all the experiments. Depending on the initial pre-tension
($F_{i}$) two possible outcomes of the inflation tests are shown in Fig. 1(c)
and(d), to be discussed later in Section 3.1.

Figure 1: Experimental setup for inflation test of pre-streched rubber tube
with fixed-fixed boundary conditions, $L$ (initial half-length), $D$ (initial
diameter), $H$ (initial thickness), $F$ (force applied/measured by Intron),
$l$ (stretched half-length), $l_{b}$ (bulge length), $d$ (deformed diameter),
$c$ (buckling amplitude), $F_{i}$ (longitudinal tension), and $h$ (deformed
thickness). (a) Longitudinal tension is applied at the upper end to reach the
required pre-stretch ($\lambda_{i}$). (b) Deformed (dashed) and undeformed
configurations during inflation subjected to fixed-fixed ends. (c) Bulge
initiation and buckling at low longitudinal tension (pre-stretch). (d) Bulge
initiation and rupture at high longitudinal tension. The scale bar in white is
2 cm.

###  2.2 Finite element analyses

Finite element analysis (FEA) offers the ability to systematically vary the
material and geometric parameters of the inflation test. We conduct FEA in two
parts. In the first part, virtual inflation tests on cylindrical hyperelastic
tubes are performed. Initial tension is varied in these simulations. The
purpose of these simulations is to verify the experimental observations of
buckling in the post bulge regime as well as validate the modelling procedure.
In the second part, we use this FEA procedure with a more biologically
relevant GOH material model for the wall.

The FEA has been conducted in three stages. In experiments, all the
displacement degrees of freedom have been constrained at the lower end, and at
the top end of the tube, as shown in Fig. 1(b). We can exploit symmetry of the
boundary conditions by analyzing only one half of the tube. Note that this
symmetry still allows the aneurysm to to grow radially asymmetric. The top end
of the tube with half length, $l$, as shown in Fig. 1(b), coincides with the
crown of the bulge and it is only constrained for the axial displacement to
reflect symmetry. The plane of symmetry is also indicated. We use ABAQUS©
Linear Perturbation Analysis (LPA) and Riks arc length method solvers to
obtain the critical bifurcation modes and post bifurcation solution,
respectively [Hibbit et al., 2012, Crisfield, 1983]. In stage 1 of the FEA,
LPA has been performed to obtain the unstable mode shapes of a cylindrical
tube subjected to a normalized unit value of prescribed axial stretch and
internal pressure. The global nodal imperfection has been introduced based on
the first two circumferential modes (global bulge shape) and the first axial
mode. It should be noted that due to manufacturing defects, these
imperfections naturally exist in the tubes subjected to the inflation test in
Section 2.1. In stage 2, the desired pre-stretch is applied and held constant
for various prescribed values, identical to that in experiments on rubber
tubes. Finally, the inflation test is simulated by increasing the internal
pressure. Both stages 2 and 3 use Riks method.

Linear shell elements with reduced integration (S4R) were used [Hibbit et al.,
2012]. We used 31 elements for the circumference and 199 elements along the
length of the tube (6169 total). Element size and amplitude of imperfection
are two important simulation parameters for the bifurcation analysis. For the
Riks method, global nodal imperfection of order of $10^{-3}$, normalized with
respect to element size, has been introduced for the first two circumferential
modes (global bulge shape) and the first axial mode (buckling mode). The
number of elements and the normalized (with respect element size) modal
amplitude of imperfection ($10^{-3}$) have been chosen based on the
convergence of the critical pressure associated with bulging instability.

###  2.3 Material models

Two material models are considered in the finite element analyses: isotropic
hyperelastic wall material and anisotropic GOH model. For the hyperelastic
model, we use two failure criteria, respectively based on Ogden’s strain
energy function (SEF) and Gent’s criterion based on the first stretch
invariant. For the GOH model we Cauchy stress based failure criterion
following [Sherifova and Holzapfel, 2019].

The Ogden’ SEF ($W_{\text{O}}$) for a cylindrical coordinate with plane stress
can be expressed as

$\begin{array}[]{c}W_{\text{O}}=\displaystyle\sum_{i=1}^{3}\frac{s_{i}}{\alpha_{i}}(\lambda_{1}^{\alpha_{i}}+\lambda_{2}^{\alpha_{i}}+\lambda_{3}^{\alpha_{i}}-3),\\\
\lambda_{1}=\lambda_{z},~{}\lambda_{2}=\lambda_{\theta},~{}\lambda_{3}=(\lambda_{z}\lambda_{\theta})^{-1},\end{array}$
(1)

where ($\lambda_{z}$ and $\lambda_{\theta}$) are principal stretches in
cylindrical coordinate system for an axisymmetric membrane. The material
parameters are given in Table. 1 based on our previous experimental
calibration [Hejazi et al., 2021]. We used two different rupture criteria: (1)
The maximum strain energy density based on Ogden’s model in Eq. (1); (2) the
maximum in-plane first stretch invariant based on Gent’s material model [Gent,
1996], which reads

$\begin{array}[]{c}I=\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}-3,\\\
\lambda_{1}=\lambda_{z},~{}\lambda_{2}=\lambda_{\theta},~{}\lambda_{3}=(\lambda_{z}\lambda_{\theta})^{-1}.\end{array}$
(2)

First, we use the calibrated material parameters for Ogden’s strain energy
function (SEF) in Eq. (1) to validate our FEA with experiments. Next, with
identical assumptions as above, namely, axisymmetric membrane and plane
stress, the Gasser-Ogden-Holzapfel (GOH) model’s SEF ($W_{\text{H}}$) can be
expressed as

$\begin{array}[]{c}W_{\text{H}}=\displaystyle\frac{c_{10}}{2}(I_{1}-3)+\displaystyle\frac{k_{1}}{k_{2}}\left(\text{exp}[k_{2}({\kappa}I_{1}+(1-3\kappa)I_{4}-1)^{2}]-1\right),\\\
I_{1}=\lambda_{z}^{2}+\lambda_{\theta}^{2}+(\lambda_{z}\lambda_{\theta})^{-2},~{}I_{4}=\lambda_{z}^{2}\sin^{2}{\gamma}+\lambda_{\theta}\cos^{2}{\gamma},\end{array}$
(3)

where, $c_{10}$, $k_{1}$, and $k_{2}$ are the material parameters regarding
the contribution of strain energy in the matrix and fibers. The structural
parameters are $\kappa$ and $\gamma$ represent the in-plane fiber gaussian
distribution and fiber orientation with respect to circumferential direction
($\theta$), respectively [Gasser et al., 2006]. The first term in the above
equation containing $c_{10}$ can represent the strain energy stored in the
extra cellular matrix and the second term is the strain energy stored in
reinforcing fibers (collagen). $k_{1}$ and $k_{2}$ respectively correspond to
the volume fraction and stiffness of collagen fibers. The parameters used in
the computations are given in Table. 1 baed on _ex-vivo_ calibration test on
AAA tissue samples [Niestrawska et al., 2019].

Table 1: Geometric and material properties of the tubes in experiment and FEA based on Eq. (1) and Eq. (3) Internal diameter, $D$ |  6.5 $mm$
---|---
Length, $2L$ |  10 - 20 $cm$
Tube wall thickness, $H$ |  1.5 $mm$
Ultimate uniaxial stretch, $\lambda_{ut}$ |  $7.85\pm 0.8$
Ogden’s SEF parameters |  $\alpha_{1}=1.04$; $\alpha_{2}=1.33$; $\alpha_{3}=4.46$
|  $s_{1}=704.8$; $s_{2}=373.2$; $s_{3}=2.8$ ($s_{i}$ in kPa)
GOH’s SEF parameters |  $c_{10}=5.9$ kPa; $k_{1}=18.2$ kPa; $k_{2}=17.7$
|  $\kappa=0.242$; $\gamma=60^{\circ}$

##  3 Results and discussion

In this section we first validate the FEA modelling procedure by comparing it
with experiments and then use the FEA subsequently for different material
models, geometric parameters, and pre-tensions.

###  3.1 Rubber tube: modeling and experiments

With fixed-fixed boundary conditions, a pre-stretched tube subjected to
internal pressure undergoes three stages of deformation as shown in Fig. 2(a):
(1) uniform expansion, (2) radial expansion and axial growth of the aneurysm,
(3) buckling and/or rupture. During the uniform expansion, the internal
pressure increases monotonically until a global long-wavelength bulge forms
with a relatively small amplitude. The formation of the long-wavelength bulge
is associated with a reduction of pressure and it rapidly localizes forming
the so-called aneurysm, or a bulge. Formation of the localized bulge leads to
a sudden pressure relaxation. As the bulge region invades the straight
portions of the tube, the axial force drops due to load relaxation. A
compressive region forms outside the bulged region, since both ends are
constrained, lateral buckling occurs. If the pre-stretch is high enough the
tube may rupture without buckling. This is because the decrease in axial force
due to an axially propagating bulge is insufficient to overcome the pre-
tension. For the case where buckling follows bulge formation, as shown in Fig.
2(a), the pressure-force curves show a good agreement between the experiments
and the FEA. The FEA shows a slightly stiffer response due to slight
differences in the material properties and discretization. In the post
buckling regime, shown in Fig. 2(b), the normalized buckling amplitude ratio
($c/d$) versus bulge length fraction ($l_{b}/l$) plot shows a similar level of
agreement. In the early stages of buckling, the lateral displacement increases
rapidly, but with the increasing buckling amplitude ($c$), we observe a
hardening behaviour that leads to a plateau (indicated by 3 in Fig. 2(b)) in
FEA. This deviation arises from the formation of a localized kink at the
bottom end and the neck of the bulge (See Fig. 2(c)). Compared to FEA, the
kink observed in the experiment, seeFig. 2(d)), is not as strong. Having
verified the validity of FEA, we now can vary material model and other
parameters.

Figure 2: Finite element analysis versus experiments on hyperelastic tubes.
(a) Longitudinal force versus internal pressure during inflation. The onset of
bulging and buckling are associated with a drop in the internal pressure and a
drop in the longitudinal force, respectively. (b) Buckling amplitude growth
versus bulge axial growth showing three regimes: (1-2) buckling initiation,
(2-3) flexural stiffening due to bulge axial growth and reduction in length of
uniform section, (3-end) buckling amplitude limited by local kinking. (c) FEA
stress distribution associated with points 1 to 3 in (b) and shear stress at
the neck of the the bulge and bottom fixed end subjected to kink. (d) Small
kink formation at the bottom fixed end of the tube during experiment.

The effect of pre-stretch on the progression of the buckling and bulge growth
is shown in Fig. 3. Here, the normalized buckling amplitude is plotted against
the fraction of the length of the tube, $l_{b}/l$, occupied by the bulge. With
larger axial propagation, the uniform section of the tube becomes shorter, and
the end kink has a more predominant effect on hardening transition in buckling
amplitude development (see the transition region in Fig. 3(a)). The axial and
radial propagation of a bulge is shown in Fig. 3(b). At the onset of bulge
propagation, the uniform expansion is followed by the formation of a long-
wavelength bulge, which then becomes localized by growing in radial size. The
radial propagation of the localized bulge gradually turns into axial
propagation until the buckling sets in. With buckling, the radial propagation
declines and the bulge length development accelerates (Fig. 3(b)). This
behaviour is a result of relaxation of axial force due to buckling, which
allows the bulge to easily propagate along the length of the tube. However,
with a higher pre-stretch, the bulge needs to invade longer until the buckling
(and hence relaxation) develops. This means that the onset of buckling is
delayed in Fig. 3(b).

Figure 3: Effect of longitudinal pre-stretch on the growth of $l_{b}$ (axial
length of the bulge) and $c$ (buckling amplitude);
$\lambda_{\theta}^{\text{max}}=d_{\text{max}}/d$. (a) By increasing tension,
buckling initiation requires larger bulge length and flexural stiffening is
more significant due to a shorter length of the uniform section. (b) The long
wave bulge profile ($l_{b}/l<0.1$) is followed by localized bulge formation
and rapid axial growth due to buckling. With higher longitudinal pre-stretch
the length of the long wave is larger and the buckling initiates at larger
diameter/circumferential stretch.

A rupture portrait based on buckling amplitude ratio ($c/d$) and maximum
circumferential stretch ($\lambda_{\theta}^{\text{max}}=d_{\text{max}}/d$) is
shown in Fig. 4(a). Two rupture criteria, of Ogden and Gent, are compared.
Using FEA, we first obtained the Gent’s first invariant stretch, Eq. (2), and
Ogden’s SEF, Eq. (1), for all the data points throughout the inflation of the
tube with different pre-stretch values ($1.75\leq\lambda_{i}\leq 1.4$). Then,
the contour lines, shown in Fig. 4(a), have been obtained for all the data
points to delineate the safe-unsafe boundary with respect to rupture. Along
these curves the pre-stretch varies as indicated in Fig. 4(a). By increasing
the pre-stretch, the buckling amplitude falls as well as the maximum
circumferential stretch. With a large enough pre-stretch (initial tension,
$F_{i}$), rupture without buckling (see the diamond data points in Fig. 4(b))
occurs. Analytical buckling and rupture envelopes are obtained by following
[Hejazi et al., 2021]. A modified post-buckling rupture line based on FEA is
shown for both Gent’s stretch invariant and Ogden’s strain energy density
based on the contour lines in Fig. 4(a). At higher values of pre-stretch,
Ogden’s criterion is conservative, and at lower pre-stretch values Gent’s
criterion is conservative. This is due to the additional tensile strain energy
induced by buckling.

Figure 4: Rupture behaviour of the rubber tube subjected to an internal
pressure and fixed longitudinal pre-stretch, $F_{\text{i}}$ (longitudinal
tension at the beginning of the inflation), $A_{0}$ (tube cross section area
in deformed configuration). (a) Rupture contour line in the plane of buckling
amplitude and maximum circumferential stretch based on Ogden’s maximum strain
energy density and Gent’s maximum stretch invariant. (b) Rupture and buckling
phase portrait showing the effect of pre-stretch and buckling on the rupture
behaviour.

###  3.2 GOH material model

We implemented the GOH material model from Eq. (3) with the parameters for the
control case are given in Table. 1.

Comparing the inflation tests on a GOH tube shown in Fig. 5(a) with Fig. 2(a)
for a hyperelastic tube, it can be observed that the decrease in pressure
after the formation of the aneurysm is less pronounced for GOH. Hence, the
propagation pressure for the aneurysm to spread along the length of the artery
is lower. The two instability modes, bulging and buckling can be identified in
the pressure-force curve in Fig. 5(b). Sudden drops in pressures and axial
force follow the bulging and buckling, respectively.

Figure 5: Inflation of a tube with arterial material properties based on fiber
reinforced hyperelastic GOH model, showing the onset of bulging and buckling.
(a) In contrast to a rubber tube, the internal pressure drop after bulge
formation is not significant during the propagation phase. (b) The pressure-
force curve shows the onset of bulging (softening behaviour of the pressure
($P$)) and onset of buckling (significant drop in the axial force ($F$)).
Here, buckling follows the formation of the aneurysm.

A 3D bifurcation map for the tube with the aortic wall is shown in Fig. 6(a).
In GOH model, compressive forces may arise in the straight regions of the tube
without a bulge formation due to stretch inversion. Buckling before bulging
can occur, which is not possible in a hyperelastic tube. In Fig. 6(b), we can
see a small drop in the pressure as the indicator of the bulge formation and
an immediate increase in buckling amplitude indicating the lateral buckling.
By comparing two different cases with small ($\lambda_{i}=1.15$) and large
($\lambda_{i}=1.4$) pre-stretches, we see how the order of bulging and
buckling may reverse (see also FEA results in Fig. 6(c) and (d)).

Figure 6: Arterial tube buckling and bulging bifurcation points (limit point
instabilities) in the plane of inflation volume versus internal pressure: $P$
(internal pressure), $V$ (internal volume), $V_{0}$ (internal volume at the
beginning of inflation),
$S_{\text{ini}}=\displaystyle\frac{1}{2}\left[\displaystyle\frac{\partial}{\partial\lambda_{\theta}}\left(\displaystyle\frac{\partial
W}{\partial\lambda_{z}}\right)\right]$ is the initial shear modulus. (a) The
bulge and buckling initiation are evident by increase in ($c$/$d$ and limit
point instability of $Pd$/$hS_{\text{ini}}$. (b) The effect of tension on the
buckling and bulge initiation indicates that with lower pre-stretch tube may
buckle before bulge initiation, but as the pre-stretch increases, bulge forms
before buckling. Similar to rubber tubes, the propagation pressure does not
change due to a change in the pre-stretch. (c) Post bulge buckling of the
arterial tube with high pre-stretch. (d) Post buckling bulge of the arterial
tube with low pre-stretch.

The growth rate of an AAA is one of the important factors in clinical decision
making. Here, we investigated the axial and radial growth and the effect of
buckling on these two parameters. Fig. 7(a) shows the progression of the
buckling amplitude ($c$) versus the bulge radial growth
($\lambda_{\theta}^{\text{max}}$). We see both cases of post bulge buckling
($\lambda_{i}=1.15$) and post-buckling bulge ($\lambda_{i}=1.4$). In both
cases, the initial stage of the buckling is associated with a surge in
amplitude ($c$), followed by a small decrease, and another increase of the
amplitude. The buckling and long-wavelength bulge initiate concurrently at the
early stages. By the formation of the localized bulge, the buckling amplitude
declines slightly, followed by significant radial growth of the bulge. With
the transition of radial growth into axial propagation, the buckling amplitude
increases further until it reaches the hardening as a result of the reduction
in uniform length and pinching effects at the end. Also, with a longer bulge
length, the maximum diameter tends to shrink as a result of relaxation of the
axial force due to buckling. The axial propagation of the bulge ($l_{b}$) is
shown in Fig. 7(b). In the first stage, we see the formation of the long-
wavelength bulge. The formation of the localized bulge results in a minor
reduction in the bulge length ($l_{b}$) and a significant increase in radial
size. In the next stage of the propagation, the maximum radial size
approximately remains constant (marginal reduction is observed) while the
bulge propagates axially.

Figure 7: Growth of a bulge in an arterial tube.
$\lambda_{\theta}^{\text{max}}=d_{\text{max}}$/$d$ (radial size), $l_{b}$/$l$
(axial size), $c$/$d$ (lateral buckling growth). (a) Bucking before
($\lambda_{\text{i}}$ =1.15) and after ($\lambda_{\text{i}}$ =1.4) bulge
formation. The buckling amplitude ($c$) slightly reduces during propagation
phase due to stress relaxation in the uniform section. (b) Initial bulge with
a long wave length transforms into a localized bulge and results in a slight
reduction in the bulge length ($l_{b}$). At the late stages of propagation,
the axial growth takes over radial growth.

Finally, the rupture map of a tube with the aortic wall with GOH material
model, subjected to axial pre-tension and internal pressure is shown in Fig.
8(a). By reducing the $c_{10}$ and increasing $k_{1}$, we can model the
elastin degradation, which is commonly observed for AAAs [Niestrawska et al.,
2019]. We showed two rupture thresholds by changing these parameters
($c_{10}$, $k_{1}$) parameter in the SEF function. The rupture criterion is
based on the maximum Cauchy stress, which is located at the crown of the bulge
on the outer wall with respect to buckling curvature. Similar to Fig. 4(a), we
obtained the failure points using FEA for various pre-stretch values and
sketched the failure boundary on buckling amplitude-circumferential stretch
plot. Both of these parameters are easy to measure using common imaging
techniques in clinical practice (e.g. CT scans and ultrasound). The control
case has a qualitatively similar failure envelope compared to the hyperelastic
tube in Fig. 4(a). With a larger buckling amplitude for low pre-tension tubes,
a smaller circumferential stretch is needed for rupture to occur. However, in
general, buckling acts as a fail safe mode, and postpones rupture. When
buckling occurs, rupture takes place at a much larger size of the aneurysm, as
shown in Fig. 8(a). The exponential stiffening is more pronounced at lower
values of $c_{10}$ and higher values of $k_{1}$, which prevents the radial
growth of the aneurysm and favours buckling. This is also evident in Fig. 8(b)
which shows higher amplitudes of buckling (lower size of aneurysm) for
decreasing $c_{10}$ and increasing $k_{1}$.

Figure 8: Rupture map for an arterial tube. (a) The effect of material
parameters on the rupture of bulged-buckled arterial tube in the plane of
maximum radial size and buckling amplitude. (b) Corresponding stress
distribution figures for the rupture contours in (a).

###  3.3 Clinical implications

In this study, we showed two mechanisms that may govern the arterial buckling.
The first one, buckling without bulge, is a result of stretch inversion, which
is sensitive to microscopic fiber orientation and composition. The second
mechanism, probably more relevant to AAAs, is the formation of the compressive
force due to axial propagation of the aneurysm. It is not clear whether the
buckling always follows the bulge or vice versa. However, many AAAs are
commonly observed in a buckled state [Fillinger et al., 2004, Ducas et al.,
2020, Hejazi et al., 2021] do not rupture at a much larger size [Choi et al.,
2021]. This study suggests that a straight aneurysm has a higher risk of
rupture compared to a buckled aneurysm.

###  3.4 Limitations

There are limitations of this study. The most significant is the lack of
haemodynamics and fluid-structure interactions. Such formidable calculations
need to be attempted in the future, specially when the artery undergoes a
structural instability like buckling. Arterial branches and the curvature
imposed by spine posterior support are not considered. Nonetheless, the
rupture maps for both the rubber tube and the tube with an aortic wall can
provide useful information regarding the rupture of tortuous AAAs by
correlating the 2D buckling amplitude and maximum diameter to the rupture
risk. We anticipate such studies in the future.

##  4 Conclusion

We performed a series of inflation tests on rubber tubes with axial pre-
tension under fixed-fixed boundary conditions. We performed finite element
calculations using ABAQUS© commercial solvers including linear perturbation
analysis and Riks method solver. The FEA code has shown reasonable agreement
with _in vitro_ experiments on hyperelastic rubber tubes. With a validated FEA
code, we implemented the properties of the aortic wall using fiber-reinforced
hyperelastic material with collagen and elastin fibers. We investigated the
bulging and buckling instabilities, and explored the bulge growth and rupture
in post bulge/buckling regime. We implemented different AAA pathogenesis
scenarios by changing the material parameters associated with fiber deposition
and degradation. We presented a rupture map in the plane of buckling amplitude
plotted against maximum circumferential stretch, both of which can be measured
using 2D clinical images (CT and ultrasound) with minimal technical
requirements (e.g. 3D image reconstruction). We concluded that the post bulge
buckling is a protective mechanism against rupture.

## Acknowledgment

We gratefully acknowledge the funding provided for this research by the
National Science and Engineering Research Council of Canada (NSERC) through
Discovery Grant. Discussions with Prof. Y. Hsiang from the Vancouver General
Hospital are thankfully acknowledged.

## References

  * Acosta et al. [2006] Acosta, S., Ögren, M., Bengtsson, H., Bergqvist, D., Lindblad, B., Zdanowski, Z., 2006\.  Increasing incidence of ruptured abdominal aortic aneurysm: a population-based study.  Journal of vascular surgery 44, 237–243.
  * Alhayani et al. [2013] Alhayani, A., Giraldo, J., Rodríguez, J., Merodio, J., 2013\.  Computational modelling of bulging of inflated cylindrical shells applicable to aneurysm formation and propagation in arterial wall tissue.  Finite elements in analysis and design 73, 20–29.
  * Alhayani et al. [2014] Alhayani, A., Rodríguez, J., Merodio, J., 2014.  Competition between radial expansion and axial propagation in bulging of inflated cylinders with application to aneurysms propagation in arterial wall tissue.  International journal of engineering science 85, 74–89.
  * Boyd et al. [2016] Boyd, A.J., Kuhn, D.C., Lozowy, R.J., Kulbisky, G.P., 2016\.  Low wall shear stress predominates at sites of abdominal aortic aneurysm rupture.  Journal of vascular surgery 63, 1613–1619.
  * Choi et al. [2021] Choi, S.H., Hejazi, M., Phani, S.A., Hsiang, Y., 2021\.  Evaluation of aortic tortuosity as a negative predictor of abdominal aortic aneurysm rupture, in: 2021 Vascular annual meeting, (SVS).
  * Crisfield [1983] Crisfield, M., 1983.  An arc-length method including line searches and accelerations.  International journal for numerical methods in engineering 19, 1269–1289.
  * Dehghani et al. [2019] Dehghani, H., Desena-Galarza, D., Jha, N., Reinoso, J., Merodio, J., 2019.  Bifurcation and post-bifurcation of an inflated and extended residually-stressed circular cylindrical tube with application to aneurysms initiation and propagation in arterial wall tissue.  Finite elements in analysis and design 161, 51–60.
  * Demirkoparan and Merodio [2017] Demirkoparan, H., Merodio, J., 2017\.  Bulging bifurcation of inflated circular cylinders of doubly fiber-reinforced hyperelastic material under axial loading and swelling.  Mathematics and Mechanics of Solids 22, 666–682.
  * Ducas et al. [2020] Ducas, A.A., Kuhn, D.C., Bath, L.C., Lozowy, R.J., Boyd, A.J., 2020.  Increased matrix metalloproteinase 9 activity correlates with flow-mediated intraluminal thrombus deposition and wall degeneration in human abdominal aortic aneurysm.  Journal of vascular surgery: vascular science 1, 190–199.
  * Fillinger et al. [2004] Fillinger, M.F., Racusin, J., Baker, R.K., Cronenwett, J.L., Teutelink, A., Schermerhorn, M.L., Zwolak, R.M., Powell, R.J., Walsh, D.B., Rzucidlo, E.M., 2004\.  Anatomic characteristics of ruptured abdominal aortic aneurysm on conventional ct scans: implications for rupture risk.  Journal of vascular surgery 39, 1243–1252.
  * Font et al. [2021] Font, A., Jha, N.K., Dehghani, H., Reinoso, J., Merodio, J., 2021.  Modelling of residually stressed, extended and inflated cylinders with application to aneurysms.  Mechanics research communications 111, 103643.
  * Fu et al. [2016] Fu, Y., Liu, J., Francisco, G., 2016.  Localized bulging in an inflated cylindrical tube of arbitrary thickness–the effect of bending stiffness.  Journal of the mechanics and physics of solids 90, 45–60.
  * Fu et al. [2012] Fu, Y.B., Rogerson, G.A., Zhang, Y.T., 2012.  Initiation of aneurysms as a mechanical bifurcation phenomenon.  International journal of non-linear mechanics 47, 179–184.
  * Fung [2013] Fung, Y.c., 2013.  Biomechanics: mechanical properties of living tissues.  Springer Science & Business Media.
  * Gasser et al. [2006] Gasser, T.C., Ogden, R.W., Holzapfel, G.A., 2006.  Hyperelastic modelling of arterial layers with distributed collagen fibre orientations.  Journal of the royal society interface 3, 15–35.
  * Gent [1996] Gent, A.N., 1996.  A new constitutive relation for rubber.  Rubber chemistry and technology 69, 59–61.
  * Gonçalves et al. [2008] Gonçalves, P.B., Pamplona, D., Lopes, S.R.X., 2008.  Finite deformations of an initially stressed cylindrical shell under internal pressure.  International journal of mechanical sciences 50, 92–103.
  * Haughton and Ogden [1979a] Haughton, D., Ogden, R., 1979a.  Bifurcation of inflated circular cylinders of elastic material under axial loading—ii. exact theory for thick-walled tubes.  Journal of the mechanics and physics of solids 27, 489–512.
  * Haughton and Ogden [1979b] Haughton, D.M., Ogden, R.W., 1979b.  Bifurcation of inflated circular cylinders of elastic material under axial loading—i. membrane theory for thin-walled tubes.  Journal of the mechanics and physics of solids 27, 179–212.
  * Hejazi et al. [2021] Hejazi, M., Hsiang, Y., Srikantha Phani, A., 2021.  Fate of a bulge in an inflated hyperelastic tube: theory and experiment.  Proceedings of the royal society A: mathematical, physical and engineering sciences 477, 20200837\.
  * Hibbit et al. [2012] Hibbit, H., Karlsson, B., Sorensen, E., 2012.  Abaqus user manual, version 6.12.  Simulia, providence, RI 545.
  * Holzapfel and Ogden [2010] Holzapfel, G.A., Ogden, R.W., 2010\.  Constitutive modelling of arteries.  Proceedings of the royal society A: mathematical, physical and engineering sciences 466, 1551–1597.
  * Holzapfel et al. [2019] Holzapfel, G.A., Ogden, R.W., Sherifova, S., 2019.  On fibre dispersion modelling of soft biological tissues: a review.  Proceedings of the royal society A: mathematical, physical and engineering sciences 475, 20180736\.
  * Humphrey [2013] Humphrey, J.D., 2013.  Cardiovascular solid mechanics: cells, tissues, and organs.  Springer Science & Business Media.
  * Länne et al. [1992] Länne, T., Sonesson, B., Bergqvist, D., Bengtsson, H., Gustafsson, D., 1992.  Diameter and compliance in the male human abdominal aorta: influence of age and aortic aneurysm.  European journal of vascular surgery 6, 178–184.
  * Lazaris [2019] Lazaris, A.M., 2019.  Regarding “the society for vascular surgery practice guidelines on the care of patients with an abdominal aortic aneurysm”.  Journal of vascular surgery 69, 975\.
  * Lee et al. [2014] Lee, A.Y., Sanyal, A., Xiao, Y., Shadfan, R., Han, H.C., 2014.  Mechanical instability of normal and aneurysmal arteries.  Journal of biomechanics 47, 3868–3875.
  * Lestringant and Audoly [2018] Lestringant, C., Audoly, B., 2018\.  A diffuse interface model for the analysis of propagating bulges in cylindrical balloons.  Proceedings of the royal society A: mathematical, physical and engineering sciences 474, 0313–0333.
  * Lestringant and Audoly [2020] Lestringant, C., Audoly, B., 2020\.  Asymptotically exact strain-gradient models for nonlinear slender elastic structures: a systematic derivation method.  Journal of the mechanics and physics of solids 136, 103730.
  * Liu et al. [2019] Liu, Y., Ye, Y., Althobaiti, A., Xie, Y.X., 2019\.  Prevention of localized bulging in an inflated bilayer tube.  International journal of mechanical sciences 153, 359–368.
  * Muluk et al. [2017] Muluk, S.L., Muluk, P.D., Shum, J., Finol, E.A., 2017\.  On the use of geometric modeling to predict aortic aneurysm rupture.  Annals of vascular surgery 44, 190–196.
  * Niestrawska et al. [2019] Niestrawska, J.A., Regitnig, P., Viertler, C., Cohnert, T.U., Babu, A.R., Holzapfel, G.A., 2019\.  The role of tissue remodeling in mechanics and pathogenesis of abdominal aortic aneurysms.  Acta biomaterialia 88, 149–161.
  * Sakalihasan et al. [2018] Sakalihasan, N., Michel, J.B., Katsargyris, A., Kuivaniemi, H., Defraigne, J.O., Nchimi, A., Powell, J.T., Yoshimura, K., Hultgren, R., 2018.  Abdominal aortic aneurysms.  Nature reviews disease primers 4, 1–22.
  * Sampson et al. [2014] Sampson, U.K., Norman, P.E., Fowkes, F.G.R., Aboyans, V., Song, Y., Harrell Jr, F.E., Forouzanfar, M.H., Naghavi, M., Denenberg, J.O., McDermott, M.M., et al., 2014\.  Estimation of global and regional incidence and prevalence of abdominal aortic aneurysms 1990 to 2010.  Global heart 9, 159–170.
  * Sherifova and Holzapfel [2019] Sherifova, S., Holzapfel, G.A., 2019\.  Biomechanics of aortic wall failure with a focus on dissection and aneurysm: a review.  Acta biomaterialia 99, 1–17.
  * Sommer et al. [2010] Sommer, G., Regitnig, P., Költringer, L., Holzapfel, G.A., 2010\.  Biaxial mechanical properties of intact and layer-dissected human carotid arteries at physiological and supraphysiological loadings.  American journal of physiology-heart and circulatory physiology 298, H898–H912.
  * Topol et al. [2021] Topol, H., Al-Chlaihawi, M.J., Demirkoparan, H., Merodio, J., 2021\.  Bulging initiation and propagation in fiber-reinforced swellable mooney–rivlin membranes.  Journal of engineering mathematics 128, 1–15.
  * Vorp [2007] Vorp, D.A., 2007.  Biomechanics of abdominal aortic aneurysm.  Journal of biomechanics 40, 1887–1902.
  * Wang and Fu [2018] Wang, J., Fu, Y., 2018.  Effect of double-fibre reinforcement on localized bulging of an inflated cylindrical tube of arbitrary thickness.  Journal of engineering mathematics 109, 21–30.
  * Wang et al. [2019] Wang, S., Guo, Z., Zhou, L., Li, L., Fu, Y., 2019\.  An experimental study of localized bulging in inflated cylindrical tubes guided by newly emerged analytical results.  Journal of the mechanics and physics of solids 124, 536–554.
  * Ye et al. [2020] Ye, Y., Liu, Y., Fu, Y., 2020\.  Weakly nonlinear analysis of localized bulging of an inflated hyperelastic tube of arbitrary wall thickness.  Journal of the mechanics and physics of solids 135, 103804.