# Second order kinematic surface fitting in anatomical structures

Wilhelm Wimmer, Hervé Delingette

<figure><img src="image_1.png" /></figure>

## To cite this version:

Wilhelm Wimmer, Hervé Delingette. Second order kinematic surface fitting in anatomical structures. Medical Image Analysis, 2025, pp.103488. 10.1016/j.media.2025.103488. hal-04937619

HAL Id: hal-04937619

https://inria.hal.science/hal-04937619v1

Submitted on 10 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ELSEVIER

<figure><img src="image_4.png" /></figure>

2

Second Order Kinematic Surface Fitting in Anatomical Structures

3

Wilhelm Wimmer<sup>a,*</sup>, Hervé Delingette<sup>b</sup>

4

5

6

ARTICLE INFO

7

Article history:

Received -

Received in final form -

Accepted -

Available online -

8

Keywords:

stationary velocity field

symmetry

core line detection

cochlea

shape classification

9

10

# 1. Introduction

11

Symmetry is a fundamental concept in mathematics, shape analysis, and various scientific disciplines. It refers to the property of a feature (points, lines, surfaces, or other mathematical objects) that remains unchanged under specific transformations, such as reflections, rotations, translations, or combinations of these operations. Symmetry plays a crucial role in understanding the organization and order found in the world, and is widely observed in both biological and non-biological systems.

12

13

14

15

For biomedical purposes, geometrical symmetry becomes particularly significant. Anatomical structures in the human body often exhibit symmetrical properties. Detecting and analyzing symmetrical features can assist in the classification of the structure into a known taxonomy, identifying abnormalities, assessing asymmetries caused by injuries, diseases, or developmental issues,

16

17

ABSTRACT

Symmetry detection and morphological classification of anatomical structures play pivotal roles in medical image analysis. The application of kinematic surface fitting, a method for characterizing shapes through parametric stationary velocity fields, has shown promising results in computer vision and computer-aided design. However, existing research has predominantly focused on first order rotational velocity fields, which may not adequately capture the intricate curved and twisted nature of anatomical structures. To address this limitation, we propose an innovative approach utilizing a second order velocity field for kinematic surface fitting. This advancement accommodates higher rotational shape complexity and improves the accuracy of symmetry detection in anatomical structures. We introduce a robust fitting technique and validate its performance through testing on synthetic shapes and real anatomical structures. Our method not only enables the detection of curved rotational symmetries (*core lines*) but also facilitates morphological classification by deriving intrinsic shape parameters related to curvature and torsion. We illustrate the usefulness of our technique by categorizing the shape of human cochleae in terms of the intrinsic velocity field parameters. The results showcase the potential of our method as a valuable tool for medical image analysis, contributing to the assessment of complex anatomical shapes.

© 2025 Elsevier B. V. All rights reserved.

<footnote><sup>a</sup>Technical University of Munich, Germany; TUM School of Medicine, Klinikum rechts der Isar, Department of Otorhinolaryngology</footnote>

<footnote><sup>b</sup>Inria, Université Côte d'Azur, Epione team, Sophia Antipolis, France</footnote>

<footnote>*Corresponding author at: Department of Otorhinolaryngology, Technical University of Munich, Trogerstr. 32, 81675 Munich, Germany
E-mail: wilhelm.wimmer@tum.de</footnote>
and aiding in surgical planning or implant design, e.g., for the cardiovascular system, bones, or the inner ear (Wimmer et al., 2019b; Rueckert et al., 1997; Eckhoff et al., 2016).

Kinematic surface fitting enables to extract stationary velocity fields from oriented point clouds (Hofer et al., 2005; Andrews and Séquin, 2013). For shape analysis, the method offers two interesting features. First, stationary velocity fields can be expressed as a compact set of parameters, allowing the global characterization of geometries analogous to primitive surfaces such as spheres, cylinders, cones, helices, conical spirals, and other surfaces of revolution. Second, stationary velocity fields can be used to compute critical points (e.g., zero-velocity convergence centers) and extract invariants to detect translational or rotational symmetry in meshes.

In the field of medical image analysis, kinematic surface fitting has been of limited use. Since the velocity fields studied contain constant rotational components, they can only be applied to shapes with straight rotationally symmetric axes (Andrews and Séquin, 2013; Hofer et al., 2005) or used to approximate more complex structures (Wimmer et al., 2019b). In contrast, naturally evolved structures often have shapes with curved rotational symmetry axes or *core lines*. An illustrative example is the human cochlea, the organ of hearing. The spiral shape consists of successive turns that are not parallel, but twisted against each other (Shin et al., 2013).

Motivated by the goal of extracting and classifying intrinsic shape properties of the human cochlea (Wimmer et al., 2019a), here we investigate the applicability of kinematic surface fitting using a nonlinear velocity field. First, we provide a brief overview of conventional (i.e., first order) kinematic surface fitting. Then, we derive a robust fitting scheme for a second order velocity field and address the identification of core lines and critical points. In the final section of the paper, we discuss the capabilities and limitations of the method using selected examples of test geometries and anatomical structures.

The paper presents several novel contributions. First, we introduce second order velocity fields to the field of kinematic surface fitting. These velocity fields allow for the representation of more complex shapes, including those with curved rotational axes or core lines. Second, we provide a robust fitting method for the proposed second order velocity field that also accounts for outliers using the heavy tailed Student-t distribution. Third, we demonstrate the applicability of our method to extracting convergence points and core lines in anatomical structures, including the aorta and the cochlea. Fourth, we illustrate how the method can be used for morphological classification of anatomical structures based on global intrinsic shape properties.

The presented work extends our methods introduced in a conference paper Wimmer et al. (2019b), which applied kinematic surface fitting using a linear velocity field to detect a straight rotational symmetry axis in human cochleae. This paper presents in more details the concept of kinematic surfaces and its links with other concepts in computational anatomy. Besides, we introduce a non-linear velocity field (i.e., quadratic with respect to position) for kinematic surface fitting which remains a linear function of its parameters. Finally, we demonstrate the first application of morphometric classification of human cochlea based on intrinsic global shape properties derived from the parameters of that velocity field.

## 2. Kinematic Surfaces

Our work aims at extracting global geometric characteristics of anatomical structures by estimating the parameters of velocity fields that can be associated with those structures when considered as kinematic surfaces. In this section, we introduce the concept of kinematic surfaces, and describe their geometric properties based on the analysis of the velocity fields.

### 2.1. *Definition*

A kinematic surface is a surface which is at every point tangent to a parameterizable velocity field. We define a 3-dimensional parametric stationary velocity field $v(p, m) \in \mathbb{R}^3$ with **m** being a compact set of parameters and **p** $\in \mathbb{R}^3$ being a point in 3-
dimensional space. A smooth surface $S \in \mathbb{R}^3$ is a kinematic surface if it is tangential to a stationary velocity field $v(\mathbf{p}, \mathbf{m})$, i.e., if for each point lying on the surface $\mathbf{p} \in S$ having a normal vector $\mathbf{n}(\mathbf{p})$, the velocity field is in the tangent plane of the surface at $\mathbf{p}$: $v(\mathbf{p}, \mathbf{m}) \cdot \mathbf{n}(\mathbf{p}) = 0$.

The vector field $v(\mathbf{p}, \mathbf{m})$ is a tangent field of a surface, but also a velocity field, since a curve $C(u)$ can be grown from any seed point $\mathbf{p}_s$ lying on a surface by solving the differential equation $\frac{\partial C}{\partial u} = v(C(u))$, with $C(0) = \mathbf{p}_s$. This curve generated from $\mathbf{p}_s$ is likely to lie on the surface $S$ unless $\mathbf{p}_s$ is on the border of $S$.

## 2.2. Link with previous work

**Symmetric Surfaces.** The notion of kinematic surfaces is closely linked with the notion of group symmetry in 2-dimensional surfaces. Indeed, a parametric surface $S(\mathbf{u}) \in \mathbb{R}^3$, $\mathbf{u} \in \Omega_S$ is considered to be symmetric if it is invariant by the application of a transformation (considered to be a diffeomorphism) of the 3D Euclidean space. More formally, if we write $T(\mathbf{p}; \boldsymbol{\theta}) : \mathbb{R}^3 \to \mathbb{R}^3$, a symmetry transformation parameterized by a parameter $\boldsymbol{\theta}$, then $\forall \mathbf{u} \in \Omega_S$, $T(S(\mathbf{u}); \boldsymbol{\theta}) \in S$. The invariance can be written more explicitly as $T(S(\mathbf{u}); \boldsymbol{\theta}) = S(\mathbf{u}^*(\mathbf{u}, \boldsymbol{\theta}))$ where $\mathbf{u}^*(\mathbf{u}, \boldsymbol{\theta})$ is the mapping associating the parameter $\mathbf{u}$ with the parameter of its symmetric point.

The notion of kinematic surfaces is related to that of continuous symmetry where there is a one dimensional continuous function $\boldsymbol{\theta}(t)$, $t \in [0, b] \in \mathbb{R}$ for which the symmetry is defined. This implies that $T(S(\mathbf{u}); \boldsymbol{\theta}(t)) = S(\mathbf{u}^*(\mathbf{u}, \boldsymbol{\theta}(t)))$, $\forall t \in [0, b]$ and $\mathbf{u}^*(\mathbf{u}, \boldsymbol{\theta}(0)) = \mathbf{u}$. In other words, there is a continuous set of symmetric transformations that maps each point to a curve on the surface which starts from that point. Then we can consider the derivative of the symmetric transformation at $t = 0$, $\frac{\partial T}{\partial t}(S(\mathbf{u}); \boldsymbol{\theta}(0)) = \frac{\partial S}{\partial t}(\mathbf{u}^*(\mathbf{u}, \boldsymbol{\theta}(0)))$ and show that this is a vector field lying on the tangent plane of the surface at $S(\mathbf{u})$. Conversely, given a stationary velocity field and a kinematic surface, it is easy to define a set of symmetric transformations by integrating the velocity field over a fixed amount of time.

The concept of continuous symmetry can be opposed to that of discrete symmetry such as plane or rotational symmetries, where there exists only a discrete set of symmetric transformations for which the surface is invariant. Besides, continuous symmetry is related to key concepts in fundamental physics (like the definition of the action of a dynamic system) and to the notion of one-parameter subgroups of Lie groups in mathematics, both of which fall outside the scope of this paper.

**Symmetry Detection in Anatomical Structures.** Detection of global plane symmetries on anatomical structures such as the brain has been explored in previous work (Liu, 2009). The proposed methods can be posed as robust optimization problems (Thirion et al., 2000), as the local estimation of texture fractal dimensions (Jayasuriya et al., 2013) or by using convolutional neural networks (Seo et al., 2022; Wu et al., 2021). The estimation of the long axis of the left ventricle has also been achieved in various ways, including fitting a line through the circle centers detected with Hough transform (van Stralen et al., 2008). Centerline extraction of blood vessels is often based on the local multiscale filtering (Sironi et al., 2015) of the image, revealing the axial symmetry of the structure.

Thus, prior work attempts to extract either global simple primitives (planes or lines) of symmetry or non-parametric local centers of symmetry. The method proposed in this paper, based on kinematic surfaces, has a complementary objective by extracting global non-trivial shapes primitives in a robust way.

**Parametric Stationary Velocity Field for Image Registration.** Stationary velocity fields have been introduced to parameterize diffeomorphisms for solving image registration problems (Vercauteren et al., 2008). However, those velocity fields often consist of multidimensional images where a speed vector is stored at each voxel. In the proposed method, we consider very compact
parametric fields consisting of a dozen parameters instead of millions. Arsigny et al. (2009) have, however, introduced polyaffine registration approaches leading to a combination of affine registrations.

This paper considers the motion invariance of shapes with a second order velocity field that generalizes a single affine transformation.

*Computation Fluid Dynamics (CFD) Flow Analysis.* The CFD community has developed a large body of techniques for visualizing flow fields (Merzkirch, 2012; Günther and Theisel, 2018) and extracting quantitative features (e.g., vortices and vortex core lines) from those flow fields. Equivalent problems arise when analyzing velocity fields associated with kinematic surfaces, especially when dealing with higher order velocity fields. Therefore, the same tools can be applied for both domains.

## 2.3. Properties of Kinematic Surfaces

Given a kinematic surface $S$ and its corresponding velocity field $v(p, m)$, the parameter set $m$ characterizes the shape of the surface. However, this parameter set is not unique since any co-linear velocity field $\mu(p)v(p, m)$, $\mu(p) \in \mathbb{R}$ is also a valid velocity field for the same kinematic surface. Furthermore, the geometric interpretation of the parameter set $m$ may not be obvious, as seen for the higher order velocity field.

The analysis of the parametric velocity field $v(p, m)$ provides a framework for extracting meaningful geometric features. For instance, points of convergence (i.e., critical points) are points where the velocity is null: $v(p_0) = 0$. These are points of accumulation or repulsion, depending on the sign of the velocity divergence $\nabla_p \cdot v(p_0, m)$ at this location. The nature of the convergence points can be analyzed by looking at the eigenvalues (that can be complex or real) of the Jacobian matrix $J(p) = \nabla_p^T v(p, m)$ (Helman and Hesselink, 1989).

Furthermore, an axis of symmetry defined by a point $d_0$ and a direction vector $d$ corresponds to a straight line, where the velocity is directed along the direction $d$ of this line: $v(d_0 + \lambda d) \times d = 0$. Any seed point located in this axis will remain on this axis when moved according to the velocity field.

A generalized notion of curve of symmetry is given by core lines (corresponding to vortex core lines in the CFD community). There exists many competing definitions of vortices and vortex core lines (Post et al., 2003), but we rely on the approach proposed by Roth and Peikert (1998) that defines those lines as the points where the torsion of streamlines is null (i.e., where the derivative of the acceleration is parallel to the velocity field). Besides, it is possible to eliminate weak or short vortices by thresholding on two local attributes. We use the parallel vector operator defined in Van Gelder and Pang (2009) for a practical implementation. The detected core lines correspond geometrically to a symmetry curve around which the surface is wrapped.

## 2.4. First Order Kinematic Surfaces

To the best of our knowledge, the concept of kinematic surfaces has only been developed so far with linear or first order velocity fields (Andrews and Séquin, 2013), i.e., fields for which the velocity is a linear function of the position: $v(p, m) = Ap + c$, where $A$ is a $3 \times 3$ matrix and $c$ a 3-dimensional vector. The considered matrix $A$ is restricted to be the sum of an antisymmetric matrix and a matrix proportional to the identity matrix, such that:

$$v(p, m) = Ap + c = r \times p + c + \gamma p, \qquad (1)$$

where $m = [r, c, \gamma]$. Those linear velocity fields lead to a large family of kinematic surfaces that have simple geometric characteristics. We detail below some cases associated with this field for increasing levels of complexity.
**Constant Velocity Field.** The constant field $\mathbf{v}(\mathbf{p}, \mathbf{m}) = \mathbf{c}$ only allows translational motions along $\mathbf{c}$ and describes a linear extrusion of a generator curve.

**Scale Velocity Field.** The scale field $\mathbf{v}(\mathbf{p}, \mathbf{m}) = \mathbf{c} + \gamma\mathbf{p}$ also includes scaling by factor $\gamma$.

**Rotational Velocity Field.** The first order rotational field $\mathbf{v}(\mathbf{p}, \mathbf{m}) = \mathbf{r} \times \mathbf{p} + \mathbf{c}$ contains a rotational axis along the vector $\mathbf{r}$ and can be used to express symmetric rotational surfaces or a straight helix.

**Rotational Scale Velocity Field.** The first order rotational scale field can be expressed as $\mathbf{v}(\mathbf{p}, \mathbf{m}) = \mathbf{r} \times \mathbf{p} + \gamma\mathbf{p} + \mathbf{c}$, encoding cylindrical, conical, rotational, and helical motions, i.e., motions composed of rotation $\mathbf{r}$, translation $\mathbf{c}$ and scaling $\gamma$ (Hofer et al., 2005). Since scaling is contained within the magnitude of the rotation vector, effectively, this velocity field has only one true degree of freedom to describe spiral shape invariance (Pottmann and Wallner, 2006).

**Convergence Point.** We can find the convergence center $\mathbf{p}_0$ by solving $\mathbf{v}(\mathbf{p}_0, \mathbf{m}) = 0$ (Hofer et al., 2005):

$$ \mathbf{p}_0 = \frac{1}{\gamma(\mathbf{r}^2 + \gamma^2)} (\gamma\mathbf{r} \times \mathbf{c} - \gamma^2\mathbf{c} - (\mathbf{r} \cdot \mathbf{c})\mathbf{r}). \quad (2) $$

The velocity field can then be alternatively expressed with respect to $\mathbf{p}_0$:

$$ \mathbf{v}(\mathbf{p}, \mathbf{m}) = \mathbf{r} \times (\mathbf{p} - \mathbf{p}_0) + \gamma(\mathbf{p} - \mathbf{p}_0). \quad (3) $$

**Streamlines.** The continuous symmetry is defined by the streamlines $C(u)$ defined by the integration of the velocity $\frac{\partial C}{\partial u} = \mathbf{v}(C(u))$.

With first order velocity field, such streamlines are conical-spiral (aka concho-spiral) curves and can be computed in closed form as:

$$ \mathbf{p}(u) = \mathbf{p}_0 + a \exp(\gamma u) (\cos(r(u - u_0))\mathbf{e}_x + \sin(r(u - u_0))\mathbf{e}_y + b\mathbf{e}_z) $$

where $u_0, a, b$ are any scalar and $\mathbf{r} = r\mathbf{e}_z$ and $(\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$ make an orthonormal frame.

**Symmetry Axis.** With the velocity of Eq.1, the straight line passing through $\mathbf{p}_0$ and directed by $\mathbf{r}$ is an axis of symmetry around which each streamline turns around.

# 3. Gaussian Kinematic Surface Fitting

## 3.1. General Approach

We consider the problem of fitting a kinematic surface given a set of points with position $\mathbf{p}_i$ and their unit surface normal $\mathbf{n}_i$. This implies finding a parametric velocity field $\mathbf{v}(\mathbf{p}, \mathbf{m})$ for which the velocity field lies in the tangent plane at each sample point of the discrete surface $S = \{\mathbf{p}_i, \mathbf{n}_i\}$:

$$ \mathbf{v}(\mathbf{p}_i, \mathbf{m}) \cdot \mathbf{n}_i = 0, \quad \forall i \quad (4) $$

Andrews and Séquin (2013) provide a generalized and basis-independent formulation of the kinematic surface fitting problem using an approximate maximum likelihood approach (Chernov, 2007). With the proposed framework, the optimal parameter vector $\mathbf{m}$ is the one minimizing the following distance $d_i(\mathbf{m})$:

$$ d_i(\mathbf{m}) = \frac{\mathbf{v}(\mathbf{p}_i, \mathbf{m}) \cdot \mathbf{n}_i}{\sqrt{\|\mathbf{v}(\mathbf{p}_i, \mathbf{m})\|^2 + w_p \|\nabla_p(\mathbf{v}(\mathbf{p}_i, \mathbf{m}) \cdot \mathbf{n}_i)\|^2}}, \quad (5) $$
where $w_p$ is a scalar. This distance measures how well the normalized velocity field lies in the tangent plane. The velocity normalization avoids selecting the trivial solution $\mathbf{v}(\mathbf{p}, \mathbf{m}) = \mathbf{0}$ and is composed of 2 terms: the velocity norm and the norm of the gradient of the velocity-normal dot product. This second term scaled by $w_p$ aims to control the gradient norm of the distance (somewhat enforcing its Lipschitz continuity), thus to make the square distance as convex as possible.

We then assume that for a set of oriented points, the distance follows a Gaussian distribution with variance $\Sigma$, that is $p(d_i|\mathbf{m}) = N(d_i|0, \Sigma)$. Then, we want to determine $\mathbf{m}$ to maximize the log-likelihood:

$$ \log p(\mathcal{S}|\mathbf{m}) = \log \prod_{i=1}^{n} p(d_i|\mathbf{m}) = -\frac{n}{2} \log 2\pi\Sigma - \frac{1}{2} \sum_{i=1}^{n} \frac{d_i^2(\mathbf{m})}{\Sigma} = -\frac{n}{2} \log 2\pi\Sigma + \mathcal{L}(\mathbf{m}). $$

Note that it is not necessary to provide a consistent orientation of the points $\mathbf{p}_i$ since the criterion is invariant to a change of orientation when replacing $\mathbf{n}_i$ with $-\mathbf{n}_i$.

Andrews and Séquin (2013) make the additional assumption that the velocity field is a linear function of its parameters $\mathbf{m}$: $\frac{\partial^2 \mathbf{v}(\mathbf{p}, \mathbf{m})}{\partial \mathbf{m}_j \partial \mathbf{m}_k} = \mathbf{0}$. This implies that there exists a matrix $\mathbf{H}(\mathbf{p})$ such that $\mathbf{v}(\mathbf{p}, \mathbf{m}) = \mathbf{H}(\mathbf{p})\mathbf{m}$ and a vector $\mathbf{f}(\mathbf{p}, \mathbf{n})$ such that $\mathbf{v}(\mathbf{p}, \mathbf{m}) \cdot \mathbf{n} = \mathbf{m} \cdot \mathbf{f}(\mathbf{p}, \mathbf{n})$. As a consequence, the square distance $d_i^2$ is a ratio between two quadratic forms of $\mathbf{m}$:

$$ d_i^2(\mathbf{m}) = \frac{(\mathbf{v}(\mathbf{p}_i, \mathbf{m}) \cdot \mathbf{n}_i)^2}{\|\mathbf{v}(\mathbf{p}_i, \mathbf{m})\|^2 + w_p \|\nabla_{\mathbf{p}}(\mathbf{v}(\mathbf{p}_i, \mathbf{m}) \cdot \mathbf{n}_i)\|^2} = \frac{\mathbf{m}^T \mathbf{M}_i \mathbf{m}}{\mathbf{m}^T \mathbf{N}_i \mathbf{m}} \quad (6) $$

where $\mathbf{M}_i$ and $\mathbf{N}_i$ are two symmetric matrices.

Then, maximizing the log-likelihood $\mathcal{L}(\mathbf{m})$ is equivalent to minimizing $\sum_{i=1}^{n} \frac{\mathbf{m}^T \mathbf{M}_i \mathbf{m}}{\mathbf{m}^T \mathbf{N}_i \mathbf{m}}$, which leads to the generalized eigenvalue problem (Chernov, 2007):

$$ \mathbf{B}_m \mathbf{m} = \mathbf{C}_m \mathbf{m}, \quad (7) $$

with $\mathbf{B}_m = \sum_{i=1}^{n} \frac{\mathbf{M}_i}{\mathbf{m}^T \mathbf{N}_i \mathbf{m}}$ and $\mathbf{C}_m = \sum_{i=1}^{n} \frac{\mathbf{m}^T \mathbf{M}_i \mathbf{m}}{(\mathbf{m}^T \mathbf{N}_i \mathbf{m})^2} \mathbf{N}_i$. The non-linear problem can be solved by iteratively computing the matrices $\mathbf{B}_m$ and $\mathbf{C}_m$ for a given estimation of $\mathbf{m}$ and then estimating $\mathbf{m}$ as the eigenvector associated with the eigenvalue closest to zero.

## 3.2. Fitting First Order Velocity Fields

We consider the first order velocity field defined in Eq.1 where $\mathbf{m} = [\mathbf{r}, \mathbf{c}, \gamma]$. It is easy to verify that the velocity field is a linear function of parameters $\mathbf{m}$.

To compute the matrices $\mathbf{M}_i$ and $\mathbf{N}_i$, we find that $\nabla_{\mathbf{p}}(\mathbf{v}(\mathbf{p}_i) \cdot \mathbf{n}_i) = \mathbf{n}_i \times \mathbf{r} + \gamma \mathbf{n}_i = [\mathbf{A}_r + \gamma \mathbf{I}] \mathbf{n}_i$ where $\mathbf{A}_r$ is the skew-symmetric matrix associated with $\mathbf{r}$, and $\mathbf{I}$ is the identity matrix. Furthermore, we have $\|\nabla_{\mathbf{p}}(\mathbf{v}(\mathbf{p}_i) \cdot \mathbf{n}_i)\|^2 = \|\mathbf{n}_i \times \mathbf{r}\|^2 + \gamma^2$ and $\|\mathbf{v}(\mathbf{p}_i)\|^2 = \|\mathbf{r} \times \mathbf{p}_i\|^2 + \gamma^2 \|\mathbf{p}_i\|^2 + 2\gamma(\mathbf{p}_i \cdot \mathbf{c}) + 2[\mathbf{r}, \mathbf{p}_i, \mathbf{c}] + \|\mathbf{c}\|^2$. With a first order velocity field, the seven-dimensional vector $\mathbf{f}(\mathbf{p}, \mathbf{n})$ is $\mathbf{f}(\mathbf{p}, \mathbf{n}) := [\mathbf{p} \times \mathbf{n}, \mathbf{n}, \mathbf{p} \cdot \mathbf{n}]$.

To find the parameter vector $\mathbf{m}$, we solve equation (7) after computing the $7 \times 7$ matrices $\mathbf{M}_i$ and $\mathbf{N}_i$ (Andrews and Séquin, 2013):

$$ \mathbf{M}_i = \mathbf{f}(\mathbf{p}_i, \mathbf{n}_i) \mathbf{f}(\mathbf{p}_i, \mathbf{n}_i)^T \quad (8) $$

and

$$ \mathbf{N}_i = \begin{bmatrix} \mathbf{A}_{\mathbf{p}_i}^T \mathbf{A}_{\mathbf{p}_i} + w_p \mathbf{A}_{\mathbf{n}_i}^T \mathbf{A}_{\mathbf{n}_i} & -\mathbf{A}_{\mathbf{p}_i} & \mathbf{0} \\ -\mathbf{A}_{\mathbf{p}_i}^T & \mathbf{I} & \mathbf{p}_i \\ \mathbf{0} & \mathbf{p}_i^T & \mathbf{p}_i \cdot \mathbf{p}_i + w_p \end{bmatrix}. \quad (9) $$
176

# 4. Second Order Velocity Field

177

## 4.1. Definition

178

179

180

It is possible to identify a rotational symmetry axis of kinematic surfaces using a linear velocity field. However, for many anatomical structures, we hypothesize that such symmetry line may be curved and not straight. This is why we propose to consider a more sophisticated velocity field that is quadratic with respect to the position **p**:

$$v(\mathbf{p}) = (\mathbf{t} \times \mathbf{p}) \times \mathbf{p} + \mathbf{r} \times \mathbf{p} + \gamma \mathbf{p} + \mathbf{c}, \qquad (10)$$

181

182

183

184

with the second order kinematic parameter **t**. With this parametric velocity field, we replace the fixed axis **r** with the moving direction **t** × **p** + **r**, which is rotating around the direction **t**. Besides, the parameter vector **m** = [**t**, **r**, **c**, **γ**] is now a 10-dimensional vector and the velocity field remains a linear function of the parameters in **m**. Obviously, this second order velocity field generalizes the first order field, which corresponds to the case **t** = **0**.

185

186

187

188

**Convergence Point.** There is no closed-form solution for isolating convergence points **p**₀ in the case of the second order velocity field. Therefore, to find the convergence point, we employ a nonlinear programming solver ("fminsearch") minimizing the objective function defined as the norm of the velocity vector field (||**v**(**p**ᵢ)||). In analogy to Eq. 3, we can find an alternative expression of the velocity field centered at **p**₀:

$$v(\mathbf{p}) = (\mathbf{t} \times (\mathbf{p} - \mathbf{p}_0) + \mathbf{r}) \times (\mathbf{p} - \mathbf{p}_0) + \gamma(\mathbf{p} - \mathbf{p}_0). \qquad (11)$$

189

190

**Streamlines.** The streamlines cannot be computed in closed form, but must be integrated through Euler or more sophisticated integration methods.

191

192

193

**Core Line.** When **t** ≠ 0, no straight axis of rotational symmetry can be found, but we can extract core lines as defined in section 2.3 that can be seen as curves of symmetry around which particles swirl. In practice, we follow vorticity extrema in the vector field, implemented using the parallel vector operator (Roth and Peikert, 1998; Van Gelder and Pang, 2009).

194

## 4.2. Kinematic Surface Fitting

195

196

197

To estimate the 10-dimensional vector **m** = [**t**, **r**, **c**, **γ**] from a set of oriented points, we proceed as described in Section 3 since we have a linear relation between the velocity field and the parameter vector **m**. For instance, the 10-dimensional vector **f**(**p**, **n**) now writes as [(**n** × **p**) × **p**, **p** × **n**, **n**, **p** ⋅ **n**] and the matrix **M**ᵢ can be obtained using Eq. 8.

198

199

To compute the normalization matrix **N**ᵢ = **N**ᵢ⁽¹⁾ + wᵢ**N**ᵢ⁽²⁾, we first express the term ||**v**(**p**ᵢ)||² as a multiplication of the matrix **N**ᵢ⁽¹⁾ and **m** where:

$$\mathbf{N}_{i}^{(1)} = \begin{bmatrix} \mathbf{F}_{p_i}^T \mathbf{F}_{p_i} & -\mathbf{F}_{p_i} \mathbf{A}_{p_i} & \mathbf{F}_{p_i} & \mathbf{0} \\ -\mathbf{A}_{p_i}^T \mathbf{F}_{p_i}^T & \mathbf{A}_{p_i}^T \mathbf{A}_{p_i} & -\mathbf{A}_{p_i} & \mathbf{0} \\ \mathbf{F}_{p_i}^T & -\mathbf{A}_{p_i}^T & \mathbf{I} & \mathbf{p}_i \\ \mathbf{0} & \mathbf{0} & \mathbf{p}_i^T & \mathbf{p}_i \cdot \mathbf{p}_i \end{bmatrix},$$

200

with **F**ᵢₚᵢ = **p**ᵢ ⊗ **p**ᵢ - (**p**ᵢ ⋅ **p**ᵢ)I, where ⊗ denotes the tensor product such that **a** ⊗ **b** = **a****b**^T. Moreover, we find that:

$$\nabla_{\mathbf{p}}(\mathbf{v}(\mathbf{p}_i) \cdot \mathbf{n}_i) = (\mathbf{t} \otimes \mathbf{p})\mathbf{n}_i + (\mathbf{t} \cdot \mathbf{p})\mathbf{n}_i - 2(\mathbf{p} \otimes \mathbf{t})\mathbf{n}_i + \mathbf{n}_i \times \mathbf{r} + \gamma\mathbf{n}_i. \qquad (12)$$

201

Then, we can express ||∇**p**(**v**(**p**ᵢ) ⋅ **n**ᵢ) ||² in terms of the matrix **N**ᵢ⁽²⁾ and **m**:
$$ N_i^{(2)} = \begin{bmatrix} G_{p_i, n_i}^T & G_{p_i, n_i} & G_{p_i, n_i} & A_{n_i} & 0 & G_{p_i, n_i} & n_i \\ A_{n_i}^T & G_{p_i, n_i}^T & A_{n_i}^T & A_{n_i} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ n_i^T & G_{p_i, n_i}^T & 0 & 0 & 0 & n_i \cdot n_i = 1 \end{bmatrix} $$

with $G_{p_i, n_i} = (p_i \cdot n_i)\mathbf{I} - 2\mathbf{p}_i \otimes \mathbf{n}_i + \mathbf{n}_i \otimes \mathbf{p}_i$. As before, we obtain our parameters $\mathbf{m}$ by assembling matrices $\mathbf{B}_m$ and $\mathbf{C}_m$ and finally solving the generalized eigenvalue problem of Eq. 7.

202

203

204

## 5. Robust Kinematic Surface Fitting

205

206

207

208

The method for estimating kinematic surface parameters in Section 3 is a least square fitting method, which is sensitive to outliers. Such outliers may come either from the erroneous estimation of surface normals $\mathbf{n}_i$ at points $\mathbf{p}_i$ or by including points that are not lying on the kinematic surfaces. Classical robust least square methods such as m-estimators or least trimmed squares(Jung, 2007) need to estimate precisely extra variables such as the percentage of outliers.

209

210

211

We adopt a parametric approach by replacing the Gaussian likelihood with a heavy tailed Student-t distribution. The benefit of this distribution is that Student-t is a Gaussian Scale Mixture which makes it amenable to a data-driven iterative estimation with closed form solutions.

212

213

214

More precisely, we now assume that $p(d_i|\mathbf{m}) = \text{St}(d_i|0, \Sigma, \nu) = \int_{z_i} \mathcal{N}(d_i|0, \Sigma/z_i) \text{Ga}(z_i|\nu/2, \nu/2)dz_i$, where $z_i$ is the variance scale variable which has a prior given by the Gamma distribution parameterized by the degrees of freedom $\nu$. When $\nu \to +\infty$, the Student-t distribution is equivalent to the Gaussian distribution and the variable $\nu$ is inversely proportional to the number of outliers.

215

216

217

The robust estimation of a kinematic surface is now achieved with an Expectation-Maximization scheme, where $z_i$ is the latent variable (Scheffler, 2008), and where $\nu$, $z_i$ and $\Sigma$ are iteratively estimated. In the E-step, we estimate the posterior distribution of $z_i$ as a Gamma distribution: $p(z_i|d_i, \nu, \Sigma) = \text{Ga}(z_i|\frac{\nu+1}{2}, \frac{\nu}{2} + \frac{d_i^2(\mathbf{m})}{2\Sigma})$. This implies that the mean value of the latent variable $z_i$ is $\frac{(\nu+1)}{(\nu+d_i^2(\mathbf{m})/\Sigma)}$. Thus, for large of $\nu$ (i.e. when the distribution is close to be Gaussian), $z_i$ approaches 1, irrespective of the distance $d_i^2(\mathbf{m})$. Conversely, if $\nu$ is small, then $z_i$ is close to zero when the distance $d_i^2(\mathbf{m})$ is much larger than the variance $\Sigma$, i.e., when dealing with outliers. The quantity $\frac{d_i^2(\mathbf{m})}{\Sigma}$ acts here as a Mahalanobis distance.

218

219

220

221

222

In the M-step, the variance $\Sigma$ and degree of freedom $\nu$ are optimized with a fixed value of $z_i$ set to its mean value. The variance is then estimated as $\Sigma = \frac{1}{n} \sum_{i=1}^n z_i d_i^2(\mathbf{m})$. The value of $\nu$ is not obtained in closed form, but is the solution of the following equation:

$$ -\psi(\frac{\nu}{2}) + \log(\frac{\nu}{2}) + 1 + \psi(\frac{\nu+1}{2}) - \log(\frac{\nu+1}{2}) + \frac{1}{n} \sum_{i=1}^{n} (\log z_i - z_i) = 0, $$

223

224

where $\psi(x)$ denotes the digamma function. The value of $\nu$ after performing a few iterations using an implementation of the Levenberg-Marquardt algorithm (Liu and Rubin, 1995).

225

226

227

228

229

230

231

Finally, after convergence of the EM algorithm, we get the estimation of the latent $z_i$ of each data point, which also indicates if that point is an inlier ($z_i \approx 1$) or an outlier ($z_i \approx 0$). The maximization of the log-likelihood with the Student-t distribution is then equivalent to solving equation (7), but with $\mathbf{B}_m = \sum_{i=1}^n z_i \frac{\mathbf{M}_i}{\mathbf{m}'\mathbf{N}_i\mathbf{m}}$ and $\mathbf{C}_m = \sum_{i=1}^n z_i \frac{\mathbf{n}_i^T \mathbf{M}_i \mathbf{m}}{(\mathbf{m}'\mathbf{N}_i\mathbf{m})^2} \mathbf{N}_i$. Thus, the robust estimation of kinematic surfaces involves two nested iterative processes, the inner one being the estimation of the variable $z_i$ softly characterizing each oriented point and the outer one estimating the velocity field parameter $\mathbf{m}$ with weighted contributions from each data point. The pseudocode for the implementation of robust kinematic surface fitting using the second order velocity field is provided in Algorithm 1.
**Algorithm 1** Robust kinematic surface fitting using second order velocity field.

**input:** $S = \{\mathbf{p}_i, \mathbf{n}_i\} \in \mathbb{R}^{3 \times n}$ (surface points and normals); $k_{max}$ (# iterations); $w_p$ (weighting factor)

**output:** $\mathbf{m} \in \mathbb{R}^{10 \times 1}$ (kinematic parameters: $[\mathbf{t}, \mathbf{r}, \mathbf{c}, \gamma]$)

1: **procedure** INITIAL FITTING

2: Centralize and normalize $\{\mathbf{p}_i\}$

3: **for** $i \le n$ **do**

4: Compute matrix $\mathbf{M}_i \in \mathbb{R}^{10 \times 10}$

5: Compute matrix $\mathbf{N}_i \in \mathbb{R}^{10 \times 10}$ using weight $w_p$ (Andrews and Séquin, 2013)

6: **end for**

7: Compute $\mathbf{N} = \sum_{i=1}^{n} \mathbf{N}_i$ and $\mathbf{M} = \sum_{i=1}^{n} \mathbf{M}_i$

8: Solve generalized eigenvalue problem: $\mathbf{M}\mathbf{m} = \lambda\mathbf{N}\mathbf{m}$

9: Select eigenvector $\mathbf{m}$ associated with the smallest eigenvalue $\lambda$

10: **end procedure**

11: **procedure** ITERATIVE ROBUST FITTING

12: Initialize confidence weights $\{z_i\}$ with 1

13: **for** $k \le k_{max}$ **do**

14: Compute $\mathbf{B}_m = \sum_{i=1}^{n} z_i \frac{\mathbf{M}_i}{m_i^T \mathbf{N}_i \mathbf{m}}$ and $\mathbf{C}_m = \sum_{i=1}^{n} z_i \frac{\mathbf{m}_i^T \mathbf{M}_i \mathbf{m}_i}{(m_i^T \mathbf{N}_i \mathbf{m}_i)^2} \mathbf{N}_i$ using $\mathbf{m}$ found from the previous iteration

15: Solve generalized eigenvalue problem: $\mathbf{B}_m\mathbf{m} = \lambda\mathbf{C}_m\mathbf{m}$

16: Select eigenvector $\mathbf{m}$ associated with the smallest eigenvalue $\lambda$

17: Compute distance metric $\{d_i(\mathbf{m})\}$

18: Fit t-distribution to distances $\{d_i(\mathbf{m})\}$ and find $\Sigma$ and $\nu$

19: Compute new weights using $z_i = (\nu + 1)/(\nu + d_i^2(\mathbf{m})/\Sigma)$

20: **end for**

21: **end procedure**

22: **procedure** POST-PROCESSING

23: Compute velocity field $\mathbf{v}(\mathbf{p}) = (\mathbf{t} \times \mathbf{p}) \times \mathbf{p} + \mathbf{r} \times \mathbf{p} + \gamma\mathbf{p} + \mathbf{c}$

24: Find convergence point $\mathbf{p}_0 = \underset{\mathbf{p}}{\operatorname{argmin}}\|\mathbf{v}(\mathbf{p})\|$

$\mathbf{p}$

25: Compute core line by starting from $\mathbf{p}_0$ by solving $\mathbf{v} \parallel \nabla \times \mathbf{v}$ (Peikert and Roth, 1999)

26: Compute streamlines by numerical integration for any seed point $\mathbf{p}_s$

27: **end procedure**
# 6. Results

## 6.1. Validation

Figures 1 to 3 illustrate examples of non-robust and robust kinematic surface fitting in synthetic test shapes. The shapes were generated by sweeping a circular cross-section along parametric centerlines. To generate the logarithmic helical spiral in Figure 1, the centerline was defined by the expression $\mathbf{p}(s) = [Re^{\gamma s} \cos s, Re^{\gamma s} \sin s, He^{\gamma s}]^T$, with $R = 1.0$, $H = 1.5$, $\gamma = -0.08$, and the centerline starting point $\mathbf{p}_s = [1, 0, -1.5]^T$ (10'000 vertices). The circular cross-section had a radius of 0.5 and is scaled by $e^{\gamma s}$. In this case, the symmetry core line corresponds to the vector $\mathbf{r} = [0, 0, 1]^T$, converging at $\mathbf{p}_0 = [0, 0, 0]^T$.

For the bent helix, the circular cross-section was swept along the centerline

$$ \mathbf{p}(s) = [(R + r \cos 18s) \cos (18s), (R + r \cos (18s)) \sin s, r \sin (18s)]^T, $$

with $R = 7.0$ and $r = 0.5$, for a sector representing 1/6 of a full circle (22'500 vertices). The centerline starting point is $\mathbf{p}_s = [7.5, 0, 0]^T$. The bent helix does not have any convergence point of minimum velocity, but rather is winding around a circular core line with radius $R$.

As third shape, a full toroidal helix was generated (90'000 vertices) using the centerline equation:

$$ \mathbf{p}(s) = [1.5r \sin (15s), (R + r \cos 15s) \cos (15s), (R + r \cos (15s)) \sin s]^T, $$

with the major bending radius $R = 5.0$ and minor radius $r = 0.5$ (with elliptical eccentricity). The centerline starting point was $\mathbf{p}_s = [0, 0, 5.5]^T$. The solution for the core line is a circle with radius $R$ (no convergence point $p_0$). For all scenarios, an outlier structure comprised of a cylinder with a radius of 0.3 and a height of 3, containing 10% of the vertices from the test shapes, was included to illustrate the impact of robust fitting.

The fitting results are listed in Tables 1 and 2. In the first example, the logarithmic helical spiral, the fit second order velocity field becomes degenerate, and is identical to the first order field solution (Fig. 1a) as $t \approx 0$. Introducing outliers (cylindrical structure) distorts the solution of the second order velocity field (Fig. 1b). However, it can be recovered with the robust fitting procedure (Fig. 1c and Fig. 1d).

The second example shows the bent helix. It is obvious that the first order velocity field is incapable of capturing the curved core line of the shape, while the second order velocity field is correctly detecting the curvature (Fig. 2a). Adding outliers (cylindrical structure), results in a degraded solution of the second order velocity field fit (Fig. 2b). The robust estimation scheme can identify the outliers in terms of a different kinematic and assigns low confidence to the data points of the outliers (Fig. 2c and Fig. 2d). The toroidal helix shows how the second order velocity field is able to capture the higher order bending of the shape and detect the symmetries of the core and center lines 3a). Again, adding the outliers worsens results, as indicated in comparison with the ideal parametric reference curves 3c). Analogously, the robust scheme assigns lower confidence to the outlier vertices (Fig. 3c,d). The outcomes derived from the constant rotational velocity field are not shown, as they do not yield a meaningful result. It has to be noted, that the used parametric curve for generating the test shape represents only an approximate solution to the field.

## 6.2. Center Line and Core Line Extraction

Figures 4 to 8 illustrate the extraction of center lines (streamlines extracted from a central seed point) and core lines in different human anatomical structures. In all figures, we show the kinematic surface fitting results using the first order (yellow curves) and second order (red curves) velocity fields. The dashed curves indicate the detected core lines for each velocity field. Continuous curves represent streamlines $C$ of the best fitting fields traced from selected central seed points $\mathbf{p}_s$ (black diamonds). The seed
<table><caption>Table 1. First order velocity field fitting on synthetic test surfaces (15 iterations, $w_p = 0.001$). Root mean square errors are indicated for the residuals $d_i$, as well for the distances between reference lines and generated lines (core and center lines). If applicable, the Euclidean distance for the theoretical and generated center points $p_0$ are also indicated.</caption><thead><tr><th rowspan="2">Shape</th><th colspan="4">Parameters</th><th colspan="3">Fitting errors</th></tr><tr><th>r</th><th>c</th><th>γ</th><th>ν</th><th>√Σ<sub>i</sub>d<sup>2</sup><sub>i</sub>(m)/n</th><th>p<sub>0</sub></th><th>Core line</th><th>Center line</th></tr></thead><tbody><tr><td>Logarithmic spiral, no outliers (Fig. 1a)</td><td>[0.00, 0.00, 0.99]<sup>T</sup></td><td>[-0.06, 0.06, 0.04]<sup>T</sup></td><td>-0.08</td><td>-</td><td>2.2·10<sup>-3</sup></td><td>5.8·10<sup>-3</sup></td><td>1.3·10<sup>-4</sup></td><td>8.0·10<sup>-2</sup></td></tr><tr><td>- with outliers, non-robust fitting (Fig. 1b)</td><td>[-0.01, 0.00, 0.99]<sup>T</sup></td><td>[-0.06, 0.06, 0.02]<sup>T</sup></td><td>-0.07</td><td>-</td><td>9.5·10<sup>-2</sup></td><td>1.8·10<sup>-1</sup></td><td>1.8·10<sup>-2</sup></td><td>9.5·10<sup>-2</sup></td></tr><tr><td>- with outliers, robust fitting (Fig. 1c)</td><td>[0.00, 0.00, 0.99]<sup>T</sup></td><td>[-0.06, 0.06, 0.02]<sup>T</sup></td><td>-0.08</td><td>0.94</td><td>9.6·10<sup>-2</sup></td><td>5.8·10<sup>-3</sup></td><td>1.2·10<sup>-4</sup></td><td>8.1·10<sup>-2</sup></td></tr><tr><td>Bent helix, no outliers (Fig. 2a)</td><td>[-0.37, 0.63, -0.06]<sup>T</sup></td><td>[-0.01, 0.01, 0.68]<sup>T</sup></td><td>0.00</td><td>-</td><td>2.8·10<sup>-1</sup></td><td>-</td><td>not detected</td><td></td></tr><tr><td>Toroidal helix, no outliers (Fig. 3a)</td><td>[0.46, 0.00, 0.85]<sup>T</sup></td><td>[0.13, 0.00, 0.20]<sup>T</sup></td><td>0.00</td><td>-</td><td>5.7·10<sup>-1</sup></td><td>-</td><td>not detected</td><td></td></tr></tbody></table>

<table><caption>Table 2. Second order velocity field fitting on synthetic test surfaces (15 iterations, $w_p = 0.001$). Root mean square errors are indicated for the residuals $d_i$ (overall velocity field fit), as well for the distances between reference lines and generated lines (core and center lines). If applicable, the Euclidean distance for the theoretical and generated center points $p_0$ are also indicated.</caption><thead><tr><th rowspan="2">Shape</th><th colspan="4">Parameters</th><th colspan="3">Fitting errors</th></tr><tr><th>t</th><th>r</th><th>c</th><th>γ</th><th>ν</th><th>√Σ<sub>i</sub>d<sup>2</sup><sub>i</sub>(m)/n</th><th>p<sub>0</sub></th><th>Core line</th><th>Center line</th></tr></thead><tbody><tr><td>Logarithmic spiral, no outliers (Fig. 1a)</td><td>[0.00, 0.00, 0.00]<sup>T</sup></td><td>[0.00, 0.00, 0.99]<sup>T</sup></td><td>[-0.06, 0.06, 0.04]<sup>T</sup></td><td>-0.08</td><td>-</td><td>2.2·10<sup>-3</sup></td><td>5.8·10<sup>-3</sup></td><td>4.5·10<sup>-5</sup></td><td>9.1·10<sup>-2</sup></td></tr><tr><td>- with outliers, non-robust fitting (Fig. 1b)</td><td>[-0.28, -0.75, 0.14]<sup>T</sup></td><td>[-0.15, 0.06, 0.55]<sup>T</sup></td><td>[-0.03, 0.03, 0.01]<sup>T</sup></td><td>-0.04</td><td>-</td><td>1.9·10<sup>-1</sup></td><td>1.8·10<sup>0</sup></td><td>1.1·10<sup>0</sup></td><td>2.0·10<sup>0</sup></td></tr><tr><td>- with outliers, robust fitting (Fig. 1c)</td><td>[0.00, 0.00, 0.00]<sup>T</sup></td><td>[0.00, 0.00, 0.99]<sup>T</sup></td><td>[-0.06, 0.06, 0.02]<sup>T</sup></td><td>-0.08</td><td>0.96</td><td>9.6·10<sup>-2</sup></td><td>4.0·10<sup>-3</sup></td><td>2.3·10<sup>-4</sup></td><td>9.2·10<sup>-2</sup></td></tr><tr><td>Bent helix, no outliers (Fig. 2a)</td><td>[0.02, -0.04, -0.59]<sup>T</sup></td><td>[0.13, -0.22, 0.02]<sup>T</sup></td><td>[0.00, 0.00, -0.77]<sup>T</sup></td><td>0.00</td><td>-</td><td>3.9·10<sup>-2</sup></td><td>-</td><td>5.0·10<sup>-2</sup></td><td>2.9·10<sup>-1</sup></td></tr><tr><td>- with outliers, non-robust fitting (Fig. 2b)</td><td>[-0.01, -0.01, -0.80]<sup>T</sup></td><td>[-0.01, 0.16, 0.00]<sup>T</sup></td><td>[-0.03, 0.05, -0.58]<sup>T</sup></td><td>-0.01</td><td>-</td><td>3.2·10<sup>-1</sup></td><td>-</td><td>1.7·10<sup>0</sup></td><td>4.2·10<sup>0</sup></td></tr><tr><td>- with outliers, robust fitting (Fig. 2c)</td><td>[0.04, -0.05, -0.63]<sup>T</sup></td><td>[0.11, -0.18, 0.01]<sup>T</sup></td><td>[0.01, 0.01, -0.75]<sup>T</sup></td><td>-0.01</td><td>1.30</td><td>2.9·10<sup>-1</sup></td><td>-</td><td>7.4·10<sup>-2</sup></td><td>2.9·10<sup>-1</sup></td></tr><tr><td>Toroidal helix, no outliers (Fig. 3a)</td><td>[0.98, 0.00, 0.00]<sup>T</sup></td><td>[-0.04, 0.00, 0.01]<sup>T</sup></td><td>[0.19, 0.00, 0.00]<sup>T</sup></td><td>0.00</td><td>-</td><td>8.2·10<sup>-2</sup></td><td>-</td><td>6.5·10<sup>-2</sup></td><td>7.3·10<sup>-3</sup></td></tr><tr><td>- with outliers, non-robust fitting (Fig. 3b)</td><td>[0.98, -0.02, -0.02]<sup>T</sup></td><td>[-0.04, -0.06, 0.01]<sup>T</sup></td><td>[0.16, 0.00, 0.00]<sup>T</sup></td><td>-0.01</td><td>-</td><td>2.4·10<sup>-1</sup></td><td>-</td><td>8.8·10<sup>-2</sup></td><td>4.1·10<sup>-1</sup></td></tr><tr><td>- with outliers, robust fitting (Fig. 3c)</td><td>[0.98, -0.01, -0.01]<sup>T</sup></td><td>[-0.04, -0.06, 0.01]<sup>T</sup></td><td>[0.16, 0.00, 0.00]<sup>T</sup></td><td>0.00</td><td>1.23</td><td>2.4·10<sup>-1</sup></td><td>-</td><td>4.3·10<sup>-2</sup></td><td>1.5·10<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Fig. 1. Kinematic surface fitting of a logarithmic spiral. Displayed are the results for the first order rotational (yellow curves) and second order rotational (red curves) velocity fields. The dashed curves indicate the detected core lines, and convergence points $p_0$ are shown as circles. Continuous curves represent example center lines C traced from a central seed point $p_s$ (black diamonds). a) Ideal case: both velocity fields converge to the same solution. b) Non-robust fitting with added cylindrical structure as outlier (arrow). c) Robust fitting of shape with outlier (arrow). d) Visualization of confidence weights $z_i$ estimated for robust fitting of the second order velocity field. Note that the reference axis and center lines are not separately depicted because they coincide with the outcome of the first-order fit.</figcaption></figure>

points were found by computing the centroid of vertices of a cross-section of the surface aligned with the rotational axis (Hofer et al., 2005). The parameters found for both velocity fields are summarized in Table 3.

Figure 9 shows an application of the robust fitting scheme for symmetry detection using a surface model of a human cochlea obtained through micro-slicing. Because the procedure involves sequential grinding of the specimens, slight displacements can occur and distort the result (Sieber et al., 2019).

For kinematic shapes, including swept surfaces along equiform motions (Hofer et al., 2005), the approach enables to infer symmetries even in cases with limited information, e.g., incomplete surfaces or loss of detail because of low resolution (see Fig. 10a,b).
<figure><img src="image_6.png" /><figcaption>Fig. 2. Kinematic surface fitting of a bent helix. Displayed are the results for the first order rotational (yellow curves) and second order rotational (red curves) velocity fields. The dashed curves indicate the detected core lines. In addition, black and blue curves are shown for the parametric reference lines. Continuous curves represent example center lines C traced from a central seed point $p_s$ (black diamonds). a) Ideal case: Second order fitting can capture the bending, in contrast to the first order approach. b) Non-robust fitting with added cylindrical structure as outlier (arrow): the second order fit does not converge properly. c) Robust fitting with outliers (arrow): the second order fitting method converges to the correct solution. d) Visualization of confidence weights $z_i$ estimated for robust fitting of the second order velocity field.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 3. Kinematic surface fitting of a toroidal helix. Displayed are the second order rotational velocity fields (red curves). Magnified views are presented in the square insets. In addition, black and blue curves are shown for the parametric reference lines. Continuous curves represent example center lines C traced from a central seed point $p_s$ (black diamonds). a) Ideal case, b) non-robust fitting with added cylindrical structure as outlier, c) robust fitting with outliers, and d) visualization of confidence weights $z_i$ estimated for robust fitting.</figcaption></figure>

## 6.3. Morphological Classification

Kinematic surface fitting can be used to extract intrinsic global shape parameters for the morphological classification of anatomical structures. Here, we demonstrate an example of human cochlea classification. The cochlea has a spiral shape with successive turns inclined towards each other, as depicted in Fig. 7. To quantify the tilt, we measure the angle $\theta$ between the basal turn and the middle turn of the cochlea, as defined by Shin et al. (2013) (see Fig. 11a). As an intrinsic shape parameter, we calculate the
<table><caption>Table 3. Kinematic surface fitting parameters for anatomical structures (15 iterations, $w_p = 0.001$). Root mean square errors are indicated for the residuals $d_i$.</caption><thead><tr><th rowspan="2">Anatomical structure</th><th colspan="5">First order rotational scale velocity field</th><th colspan="5">Second order rotational scale velocity field</th></tr><tr><th>r</th><th>c</th><th>γ</th><th>ν</th><th>√Σd<sub>i</sub><sup>2</sup>(m)/n</th><th>t</th><th>r</th><th>c</th><th>γ</th><th>ν</th><th>√Σd<sub>i</sub><sup>2</sup>(m)/n</th></tr></thead><tbody><tr><td>Aortic arch (Fig. 4)</td><td>[0.83, -0.43, -0.27]<sup>T</sup></td><td>[0.03, 0.04, 0.11]<sup>T</sup></td><td>-0.20</td><td>4.68</td><td>2.2·10<sup>-1</sup></td><td>[0.05, 0.76, 0.46]<sup>T</sup></td><td>[0.36, -0.17, -0.16]<sup>T</sup></td><td>[0.04, 0.07, 0.06]<sup>T</sup></td><td>-0.09</td><td>5.52</td><td>1.9·10<sup>-1</sup></td></tr><tr><td>Left ventricle (Fig. 5)</td><td>[0.02, -0.16, 0.99]<sup>T</sup></td><td>[0.02, 0.008, 0.003]<sup>T</sup></td><td>-0.002</td><td>6.89</td><td>1.2·10<sup>-1</sup></td><td>[0.24, -0.18, -0.03]<sup>T</sup></td><td>[0.03, -0.14, 0.98]<sup>T</sup></td><td>[0.001, -0.005, 0.0004]<sup>T</sup></td><td>-0.002</td><td>8.77</td><td>1.1·10<sup>-1</sup></td></tr><tr><td>Rib (Fig. 6)</td><td>[0.12, -0.11, -0.72]<sup>T</sup></td><td>[0.06, -0.33, 0.24]<sup>T</sup></td><td>-0.53</td><td>3.17</td><td>1.1·10<sup>-1</sup></td><td>[0.37, -0.70, 0.20]<sup>T</sup></td><td>[0.06, -0.16, -0.43]<sup>T</sup></td><td>[0.02, -0.17, 0.09]<sup>T</sup></td><td>-0.27</td><td>2.74</td><td>8.3·10<sup>-2</sup></td></tr><tr><td>Cochlea (Fig. 7)</td><td>[0.14, -0.09, 0.97]<sup>T</sup></td><td>[0.09, 0.13, 0.03]<sup>T</sup></td><td>-0.08</td><td>2.51</td><td>1.8·10<sup>-1</sup></td><td>[0.46, 0.23, 0.17]<sup>T</sup></td><td>[0.14, -0.03, 0.81]<sup>T</sup></td><td>[0.10, 0.12, 0.02]<sup>T</sup></td><td>-0.13</td><td>2.19</td><td>1.6·10<sup>-1</sup></td></tr><tr><td>Semicircular canal (Fig. 8)</td><td>[0.53, 0.76, -0.36]<sup>T</sup></td><td>[0.03, -0.13, 0.01]<sup>T</sup></td><td>-0.005</td><td>1.84</td><td>1.5·10<sup>-1</sup></td><td>[0.40, 0.43, 0.39]<sup>T</sup></td><td>[0.47, 0.42, -0.31]<sup>T</sup></td><td>[0.08, -0.04, -0.04]<sup>T</sup></td><td>-0.03</td><td>1.80</td><td>1.7·10<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 4. Kinematic surface fitting of a human aortic arch (Online, 2023a). Displayed are the results for the first order rotational (yellow curves) and second order rotational (red curves) velocity fields. The dashed curves indicate the detected core lines, and convergence points $p_0$ are shown as circles. Continuous curves represent example center lines C traced from a central seed point $p_s$ (black diamonds).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. Kinematic surface fitting of a human left ventricle (Bai et al., 2015). Displayed are the results for the first order rotational (yellow curves) and second order rotational (red curves) velocity fields. The dashed curves indicate the detected core lines, and convergence points $p_0$ are shown as circles. Continuous curves represent example center lines C traced from a central seed point $p_s$ (black diamonds).</figcaption></figure>

magnitude of the projection of $t$ on $r$, that is, $||proj_r t|| = \frac{t \cdot r}{||r||}$. The correlation between the parameter $||proj_r t||$ and the measured tilt angle is demonstrated in a set of 31 human cochleae (Wimmer et al., 2019a; Sieber et al., 2019) (see Fig. 11b). This association could be used to classify a cochlea in terms of coiling and tilting characteristics, i.e., 'sloping' vs. 'rollercoaster' shapes (Avci et al., 2014). For comparison, we compute an approximation of the total curvature of the surface mesh (sum of Gauss curvatures in all vertices) to relate with the surface topology (Rusinkiewicz, 2004). Similar to the kinematic parameters, total curvature is an intrinsic shape property that can be used to classify the shapes in our example, however, no correlation can be observed (see Fig. 11c).
<figure><img src="image_10.png" /><figcaption>Fig. 6. Kinematic surface fitting of a human rib (Online, 2023b). Displayed are the results for the first order rotational (yellow curves) and second order rotational (red curves) velocity fields. The dashed curves indicate the detected core lines, and convergence points $p_0$ are shown as circles. Continuous curves represent example center lines $C$ traced from a central seed point $p_s$ (black diamonds)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7. Example of a right human cochlea (Sieber et al., 2019). Displayed are the results for the first order (yellow curves) and second order (red curves) velocity fields. The dashed curves indicate the detected core lines, and convergence points $p_0$ are shown as circles. Continuous curves represent example center lines $C$ traced from a central seed point $p_s$ (black diamonds). Note how the first order rotational scale velocity field (corresponding to the generation of a logarithmic spiral) fails to follow the twist between turns. The second order velocity field, however, can capture this degree of freedom.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 8. Kinematic surface fitting of a human semicircular canal (Wimmer et al., 2019a). Displayed are the results for the first order rotational (yellow curves) and second order rotational (red curves) velocity fields. The dashed curves indicate the detected core lines, and convergence points $p_0$ are shown as circles. Continuous curves represent example center lines $C$ traced from a central seed point $p_s$ (black diamonds).</figcaption></figure>

## 6.4. Limitations

It is important to note that the second-order velocity field can only describe shapes that belong to its solution set. Consequently, further work is necessary to better comprehend the solutions available to shape fitting. If the surfaces intended to fit deviate from the set of solutions, the method will fail.

Another limitation of the second-order velocity fitting technique is its need for an adequate number of vertices to reach conver-
a

<figure><img src="image_13.png" /></figure>

b

<figure><img src="image_14.png" /></figure>

c

<figure><img src="image_15.png" /></figure>

d

<figure><img src="image_16.png" /><figcaption>Fig. 9. Comparison between robust and non-robust second order kinematic surface fitting of a human cochlea surface found from micro-slicing (Sieber et al., 2019). a) Side view and b) top view of the cochlea with non-robust (blue curves) and robust (red curves) fitted centerlines. The non-robust fitting performs worse because it includes artifacts introduced by the microslicing procedure during preparation of the data set. The corresponding maps c) and d) visualize the assigned confidence weights $z_i$. The red spots on the surface indicate zones that indicate displacement caused by microslicing. In addition, the structures at the base and the apex of the cochlea are identified as different kinematic motions by the robust algorithm.</figcaption></figure>

a

<figure><img src="image_17.png" /></figure>

b

<figure><img src="image_18.png" /></figure>

c

<figure><img src="image_19.png" /><figcaption>Fig. 10. Examples of second-order kinematic fitting in a surface model of a human cochlea found from magnetic-resonance-imaging. a) Kinematic surface fitting and centerline detection (red curve) compared to curved skeleton algorithm (blue curve) (Gao et al., 2024). The skeletonization algorithm is lacking detail for the middle and apical turns of the cochlea (1'500 vertices). In contrast, the course of the centerline is inferred from the basal turn using the kinematic approach. b) Same mesh as in a) but cropped in half (600 vertices). The overall kinematic motion is still detectable using the algorithm. c) Same mesh as in a) but with reduced detail (250 vertices). The algorithm fails to converge to a reasonable solution. The dashed curves indicate the detected core lines, convergence points $p_0$ are shown as circles, and seed points $p_s$ are indicated as black diamonds.</figcaption></figure>

gence, differing from the first-order approach. Because of the additional torsion parameter, the solution is more prone to over-fitting.

An example is demonstrated in Figure 10c). However, more in-depth research is necessary to fully comprehend the constraints.

291

292

## 7. Discussion

293

We introduce a new approach for fitting kinematic surfaces to anatomical structures using a second order velocity field. This advancement permits for the incorporation of more intricate shapes and enhances the precision of symmetry detection. Our method

294

295
<figure><img src="image_20.png" /><figcaption>Fig. 11. Relationship between a) the measured tilt angle $\theta$ (Shin et al., 2013) and b) the intrinsic shape parameter $||proj_r t||$ for 31 samples (Wimmer et al., 2019a; Sieber et al., 2019). For comparison, c) the relationship between the total curvature of the surfaces and the tilt angle is displayed. The linear regression lines and the coefficients of determination are shown.</figcaption></figure>

not only permits the recognition of curved rotational symmetries (core lines), but also facilitates morphological categorization by computing intrinsic shape parameters related to curvature and torsion.

The nature of the velocity field enables to identify rotational symmetries of first and second orders. For instance, in the human cochlea (Fig. 7), the highest order core line is the modiolar axis, while a generator curve traced out from the center of the duct represents the symmetry of the duct. Depending on the purpose, detection of a first order rotational symmetry may be sufficient, e.g., in the semicircular canal the detection of the first order rotational symmetry axis represents its orientation (Fig. 8). Similarly, a first order rotational field can be used as an approximation for symmetry detection, e.g., in a heart ventricle (Fig. 5) or for modiolar axis detection in cochleae as shown by Wimmer et al. (2019b)).

It is noteworthy that the recognition of revolutionary or spiraling core lines is a frequent problem that arises when attempting to identify vortices (i.e., rotary motions) in the examination of fluid flows in various areas such as medicine, technology, meteorology, and more (Günther and Theisel, 2018). We could therefore envision applications of the technique in biomedical fluid mechanics to extract parametric core lines to define axial flow and link duct geometries to flow phenomena, which is relevant in the cardiovascular system and inner ear (?Harte et al., 2023; Stokes et al., 2023).

Additional applications lie in segmentation and surface generation, where velocity fields can be used to trace out kinematic surfaces using generator curves (Hofer et al., 2005). It is further possible to explore how useful our approach is for modeling the biological growth of tissue, e.g., bones (Xu et al., 2014).

The global kinematic characterization of shapes may be desirable for the classification of suitable structures, e.g., the cochlea (Fig. 7). However, our method is limited in that it does not enable for a local or piecewise extraction of surface features, including multiscale representations. The algorithm should be adapted to be responsive to local variations in symmetry.

The introduced second order velocity field produces a class of kinematic surfaces that extends the possibilities of the first order field. However, the number of described shapes is by no means exhaustive and limited to a certain set of solutions. For example, a toroidal helix can only be approximated by the velocity field. In a broader sense, we are not aware of a closed-form solution to detection core lines and convergence points.
319

# 8. Conclusions

320

321

322

323

This paper presents a novel approach for the fitting of kinematic surfaces to anatomical structures using a second order rotational scale velocity field. This advancement allows for the detection of more complex shapes, including those with curved rotational axes or core lines, enhancing the precision of symmetry detection and morphological classification. The method is capable of identifying rotational symmetries of both first and second orders, making it versatile for various applications in medical image analysis.

324

325

326

The results of our validation experiments demonstrate the effectiveness of the proposed approach, particularly in cases where first order velocity fields fail to capture the underlying shape characteristics. The introduction of a robust fitting scheme further enhances the method's resilience to outliers, making it suitable for real-world biomedical data.

327

328

329

330

331

332

333

By extending the capabilities of kinematic surface fitting, this approach opens up new possibilities for the characterization and classification of anatomical structures based on their intrinsic shape properties. It has the potential to be applied in various biomedical contexts, including the analysis of cardiovascular structures, the inner ear, and other complex biological systems (Harte et al., 2024; Suito et al., 2013). Our approach could be of use in radiology, orthopedics, and neurology, where anatomical shape examination and symmetry core line detection is necessary for diagnosis, treatment preparation, and surgical procedure planning (Weber et al., 2017). There are potential uses of biomedical fluid mechanics that are worth exploring, as well as the application to model tissue growth, e.g., bone formation.

334

# CRediT Authorship Contribution Statement

335

336

337

Wilhelm Wimmer: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing – original draft, Writing – review & editing, Visualization, Funding acquisition. Hervé Delingette: Conceptualization, Methodology, Formal analysis, Writing – original draft, Writing – review & editing.

338

# Acknowledgments

339

340

341

This work was supported by the Swiss National Science Foundation under project grant number 205321_200850 and by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.

342

# References

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

Andrews, J., Séquin, C.H., 2013. Generalized, basis-independent kinematic surface fitting. Computer-Aided Design 45, 615–620.

Arsigny, V., Commowick, O., Ayache, N., Pennec, X., 2009. A fast and log-euclidean polyaffine framework for locally linear registration. Journal of Mathematical Imaging and Vision 33, 222–238.

Avci, E., Nauwelaers, T., Lenarz, T., Hamacher, V., Kral, A., 2014. Variations in microanatomy of the human cochlea. Journal of Comparative Neurology 522, 3245–3261.

Bai, W., Shi, W., de Marvao, A., Dawes, T.J., O'Regan, D.P., Cook, S.A., Rueckert, D., 2015. A bi-ventricular cardiac atlas built from 1000+ high resolution mr images of healthy subjects and an analysis of shape and motion. Medical image analysis 26, 133–145.

Chernov, N., 2007. On the convergence of fitting algorithms in computer vision. Journal of Mathematical Imaging and Vision 27, 231–239.

Eckhoff, D.G., Jacofsky, D.J., Springer, B.D., Dunbar, M., Cherian, J.J., Elmallah, R.K., Mont, M.A., Greene, K.A., 2016. Bilateral symmetrical comparison of femoral and tibial anatomic features. The Journal of Arthroplasty 31, 1083–1090.

Gao, X., Loriot, S., Tagliasacchi, A., 2024. Triangulated surface mesh skeletonization, in: CGAL User and Reference Manual. 5.6.1 ed.. CGAL Editorial Board. URL: https://doc.cgal.org/5.6.1/Manual/packages.html#PkgSurfaceMeshSkeletonization.

Günther, T., Theisel, H., 2018. The state of the art in vortex extraction, in: Computer Graphics Forum, Wiley Online Library. pp. 149–173.

Harte, N., Obrist, D., Caversaccio, M., Lajoinie, G., Wimmer, W., 2024. Transverse flow under oscillating stimulation in helical square ducts with cochlea-like geometrical curvature and torsion. European Journal of Mechanics - B/Fluids 107, 165–174. URL: https://www.sciencedirect.com/science/article/pii/S099775462400089X. doi: https://doi.org/10.1016/j.euromechflu.2024.07.001.

Harte, N., Obrist, D., Caversaccio, M., Lajoinie, G.P., Wimmer, W., 2023. Wall shear stress and pressure fluctuations under oscillating stimulation in helical square ducts with cochlea-like geometrical curvature and torsion. EMBC, Sydney.

Harte, N., Obrist, D., Versluis, M., Jebbink, E.G., Caversaccio, M., Wimmer, W., Lajoinie, G., Second order and transverse flow visualization through three-dimensional particle image velocimetry in millimetric ducts 159, 111296. URL: https://linkinghub.elsevier.com/retrieve/pii/S0894177724001651, doi: 10.1016/j.expthermflusci.2024.111296.
364 Helman, J., Hesselink, L., 1989. Representation and display of vector field topology in fluid flow data sets. Computer 22, 27–36. doi:10.1109/2.35197.

365 Hofer, M., Odehnal, B., Pottmann, H., Steiner, T., Wallner, J., 2005. 3d shape recognition and reconstruction based on line element geometry, in: Tenth IEEE International Conference on Computer Vision (ICCV’05) Volume 1, IEEE, pp. 1532–1538.

367 Jayasuriya, S.A., Liew, A.W.C., Law, N.F., 2013. Brain symmetry plane detection based on fractal analysis. Computerized Medical Imaging and Graphics 37, 568–580.

368 Jung, K.M., 2007. Least trimmed squares estimator in the errors-in-variables model. Journal of Applied Statistics 34, 331–338. doi:10.1080/02664760601004973.

369 Liu, C., Rubin, D.B., 1995. ML estimation of the t distribution using em and its extensions, ecm and ecme. Statistica Sinica, 19–39.

371 Liu, S.X., 2009. Symmetry and asymmetry analysis and its implications to computer-aided diagnosis: A review of the literature. Journal of biomedical informatics 42, 1056–1064.

373 Merzkirch, W., 2012. Flow visualization. Elsevier.

374 Online, 2023a. Aorta model - ultimaker thingiverse. https://www.thingiverse.com/thing:2277464. Accessed on 2023-07-26.

375 Online, 2023b. Rib cage model - nih 3d database. https://3d.nih.gov/entries/3dpx-016836. Accessed on 2023-07-26.

376 Peikert, R., Roth, M., 1999. The "parallel vectors" operator-a vector field visualization primitive, in: Proceedings Visualization'99 (Cat. No. 99CB37067), Ieee. pp. 263–532.

378 Post, F.H., Vrolijk, B., Hauser, H., Laramee, R.S., Doleisch, H., 2003. The state of the art in flow visualisation: Feature extraction and tracking. Computer Graphics Forum 22, 775–792. doi:https://doi.org/10.1111/j.1467-8659.2003.00723.x.

380 Pottmann, B.O.H., Wallner, J., 2006. Equiform kinematics and the geometry of line elements. Contributions to Algebra and Geometry 47, 567–582.

381 Roth, M., Peikert, R., 1998. A higher-order method for finding vortex core lines, in: Proceedings Visualization'98 (Cat. No. 98CB36276), IEEE. pp. 143–150.

382 Rueckert, D., Burger, P., Forbat, S., Mohiaddin, R., Yang, G.Z., 1997. Automatic tracking of the aorta in cardiovascular mr images using deformable models. IEEE Transactions on medical imaging 16, 581–590.

384 Rusinkiewicz, S., 2004. Estimating curvatures and their derivatives on triangle meshes, in: Proceedings. 2nd International Symposium on 3D Data Processing, Visualization and Transmission, 2004. 3DPVT 2004., IEEE. pp. 486–493.

386 Scheffler, C., 2008. A derivation of the em updates for finding the maximum likelihood parameter estimates of the student'st distribution. Technical note. URL www.inference.phy.cam.ac.uk/cs482/publications/scheffler2008derivation.pdf.

388 Seo, A., Kim, B., Kwak, S., Cho, M., 2022. Reflection and rotation symmetry detection via equivariant learning, in: 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), IEEE Computer Society, Los Alamitos, CA, USA. pp. 9529–9538. URL: https://doi.ieeecomputersociety.org/10.1109/CVPR52688.2022.00932, doi:10.1109/CVPR52688.2022.00932.

391 Shin, K.J., Lee, J.Y., Kim, J.N., Yoo, J.Y., Shin, C., Song, W.C., Koh, K.S., 2013. Quantitative analysis of the cochlea using three-dimensional reconstruction based on microcomputed tomographic images. The Anatomical Record 296, 1083–1088.

393 Sieber, D., Erfurt, P., John, S., Santos, G.R.D., Schurzig, D., Sørensen, M.S., Lenarz, T., 2019. The openear library of 3d models of the human temporal bone based on computed tomography and micro-slicing. Scientific data 6, 1–9.

395 Sironi, A., Türetken, E., Lepetit, V., Fua, P., 2015. Multiscale centerline detection. IEEE Transactions on Pattern Analysis and Machine Intelligence 38, 1327–1341.

396 Stokes, C., Ahmed, D., Lind, N., Haupt, F., Becker, D., Hamilton, J., Muthurangu, V., von Tengg-Kobligk, H., Papadakis, G., Balabani, S., et al., 2023. Aneurysmal growth in type-b aortic dissection: assessing the impact of patient-specific inlet conditions on key haemodynamic indices. Journal of the Royal Society Interface 20, 20230281.

399 Suito, H., Ueda, T., Sze, D., 2013. Numerical simulation of blood flow in the thoracic aorta using a centerline-fitted finite difference approach. Japan Journal of Industrial and Applied Mathematics 30, 701–710.

401 Thirion, J.P., Prima, S., Subsol, G., Roberts, N., 2000. Statistical analysis of normal and abnormal dissymmetry in volumetric medical images. Medical Image Analysis 4, 111–121.

403 Van Gelder, A., Pang, A., 2009. Using pvsolve to analyze and locate positions of parallel vectors. IEEE Transactions on Visualization and Computer Graphics 15, 682–695.

405 van Stralen, M., Leung, K., Voormolen, M., de Jong, N., van der Steen, A., Reiber, J., Bosch, J., 2008. Time continuous detection of the left ventricular long axis and the mitral valve plane in 3-d echocardiography. Ultrasound in Medicine and Biology 34, 196–207. URL: https://www.sciencedirect.com/science/article/pii/S0301562907003869, doi:https://doi.org/10.1016/j.ultrasmedbio.2007.07.016.

408 Vercauteren, T., Pennec, X., Perchant, A., Ayache, N., 2008. Symmetric log-domain diffeomorphic registration: A demons-based approach, in: International conference on medical image computing and computer-assisted intervention, Springer. pp. 754–761.

410 Weber, S., Gavaghan, K., Wimmer, W., Williamson, T., Gerber, N., Anso, J., Bell, B., Feldmann, A., Rathgeb, C., Matulic, M., et al., 2017. Instrument flight to the inner ear. Science robotics 2, eaal4916.

412 Wimmer, W., Anschuetz, L., Weder, S., Wagner, F., Delingette, H., Caversaccio, M., 2019a. Human bony labyrinth dataset: Co-registered ct and micro-ct images, surface models and anatomical landmarks. Data in brief 27, 104782.

414 Wimmer, W., Vandersteen, C., Guevara, N., Caversaccio, M., Delingette, H., 2019b. Robust cochlear modiolar axis detection in ct, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 3–10.

416 Wu, H., Chen, X., Li, P., Wen, Z., 2021. Automatic symmetry detection from brain mri based on a 2-channel convolutional neural network. IEEE Transactions on Cybernetics 51, 4464–4475. doi:10.1109/TCYB.2019.2952937.

418 Xu, P.F., Houssin, N., Ferri-Lagneau, K.F., Thisse, B., Thisse, C., 2014. Construction of a vertebrate embryo from two opposing morphogen gradients. Science 344, 87–89.