# Second Order Kinematic Surface Fitting in Anatomical Structures

Wilhelm Wimmer  Hervé Delingette  Technical University of Munich, Germany; TUM
School of Medicine, Klinikum rechts der Isar, Department of
Otorhinolaryngology  Université Côte d’Azur, France; Inria Sophia Antipolis,
Epione

###### 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.

###### keywords:

stationary velocity field, symmetry, cochlea, shape classification

##  1 Introduction

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.

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, and aiding in surgical planning or implant design, e.g.,
for the cardiovascular system, bones, or the inner ear [34, 22, 6].

Kinematic surface fitting enables to extract stationary velocity fields from
oriented point clouds [11, 1]. 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 [1, 11] or used to approximate more complex structures [34]. 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 [25].

Motivated by the goal of extracting and classifying intrinsic shape properties
of the human cochlea [33], 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.
Forth, 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. [34], 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 $\mathbf{v}(\mathbf{p},\mathbf{m})\in\mathbb{R}^{3}$
with $\mathbf{m}$ being a compact set of parameters and
$\mathbf{p}\in\mathbb{R}^{3}$ being a point in 3-dimensional space. A smooth
surface $\mathcal{S}\in\mathbb{R}^{3}$ is a kinematic surface if it is
tangential to a stationary velocity field $\mathbf{v}(\mathbf{p},\mathbf{m})$,
i.e., if for each point lying on the surface $\mathbf{p}\in\mathcal{S}$ having
a normal vector $\mathbf{n}(\mathbf{p})$, the velocity field is in the tangent
plane of the surface at $\mathbf{p}$:
$\mathbf{v}(\mathbf{p},\mathbf{m})\cdot\mathbf{n}(\mathbf{p})=0$.

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

###  2.2 Link with previous work

#### Symmetric Surfaces

The notion of kinematic surfaces is closely linked with the notion of group
symmetry in 3-dimensional surfaces. Indeed, a parametric surface
$\mathcal{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}\longrightarrow\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}^{\star}(\mathbf{u},\boldsymbol{\theta}))$
where $\mathbf{u}^{\star}(\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}^{\star}(\mathbf{u},\boldsymbol{\theta}(t)))$,
$\forall t\in[0,b]$ and
$\mathbf{u}^{\star}(\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}^{\star}(\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 [15]. The proposed methods can be
posed as robust optimization problems [29], as the local estimation of texture
fractal dimensions [12] or by using convolutional neural networks [24, 35].
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 [31]. Centerline extraction of blood vessels is often
based on the local multiscale filtering [27] 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 [32]. 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. [2] 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 [16, 7] 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 $\mathcal{S}$ and its corresponding velocity field
$\mathbf{v}(\mathbf{p},\mathbf{m})$, the parameter set $\mathbf{m}$
characterizes the shape of the surface. However, this parameter set is not
unique since any co-linear velocity field
$\mu(\mathbf{p})\mathbf{v}(\mathbf{p},\mathbf{m})$,
$\mu(\mathbf{p})\in\mathbb{R}$ is also a valid velocity field for the same
kinematic surface. Furthermore, the geometric interpretation of the parameter
set $\mathbf{m}$ may not be obvious, as seen for the higher order velocity
field.

The analysis of the parametric velocity field
$\mathbf{v}(\mathbf{p},\mathbf{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:
$\mathbf{v}(\mathbf{p}_{0})=0$. These are points of accumulation or repulsion,
depending on the sign of the velocity divergence
$\nabla_{p}\cdot\mathbf{v}(\mathbf{p}_{0},\mathbf{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
$\mathbf{J}(\mathbf{p})=\nabla^{T}_{p}\mathbf{v}(\mathbf{p},\mathbf{m})$ [10].

Furthermore, an _axis of symmetry_ defined by a point $\mathbf{d}_{0}$ and a
direction vector $\mathbf{d}$ corresponds to a straight line, where the
velocity is directed along the direction $\mathbf{d}$ of this line:
$\mathbf{v}(\mathbf{d}_{0}+\lambda\mathbf{d})\times\mathbf{d}=\mathbf{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 [19], but we rely on
the approach proposed by Roth and Peikert [21] 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 [30] 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 [1], i.e., fields
for which the velocity is a linear function of the position:
$\mathbf{v}(\mathbf{p},\mathbf{m})=\mathbf{A}\mathbf{p}+\mathbf{c}$, where
$\mathbf{A}$ is a $3\times 3$ matrix and $\mathbf{c}$ a 3-dimensional vector.
The considered matrix $\mathbf{A}$ is restricted to be the sum of an
antisymmetric matrix and a matrix proportional to the identity matrix, such
that:

$\mathbf{v}(\mathbf{p},\mathbf{m})=\mathbf{A}\mathbf{p}+\mathbf{c}=\mathbf{r}\times\mathbf{p}+\mathbf{c}+\gamma\mathbf{p},$
(1)

where $\mathbf{m}=[\mathbf{r},\mathbf{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$ [11]. 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 [20].

#### Convergence Point

We can find the convergence center $\mathbf{p}_{0}$ by solving
$\mathbf{v}(\mathbf{p}_{0},\mathbf{m})\\\ =0$ [11]:

$\mathbf{p}_{0}=\frac{1}{\gamma(\mathbf{r}^{2}+\gamma^{2})}\Big{(}\gamma\mathbf{r}\times\mathbf{c}-\gamma^{2}\mathbf{c}-(\mathbf{r}\cdot\mathbf{c})\mathbf{r}\Big{)}.$
(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}).$
(3)

#### Streamlines

The continuous symmetry is defined by the streamlines $\mathcal{C}(u)$ defined
by the integration of the velocity $\frac{\partial\mathcal{C}}{\partial
u}=\mathbf{v}(\mathcal{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)\left(\cos(r(u-u_{0})\mathbf{e}_{x}+\sin(r(u-u_{0})\mathbf{e}_{y}+b\mathbf{e}_{z}\right)$

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
$\mathcal{S}=\\{\mathbf{p}_{i},\mathbf{n}_{i}\\}$:

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

Andrews and Séquin [1] provide a generalized and basis-independent formulation
of the kinematic surface fitting problem using an approximate maximum
likelihood approach [5]. 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_{\mathbf{p}}\|\nabla_{\mathbf{p}}(\mathbf{v}(\mathbf{p}_{i},\mathbf{m})\cdot\mathbf{n}_{i})\|^{2}}},$
(5)

where $w_{\mathbf{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_{\mathbf{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})={\cal N}(d_{i}|0,\Sigma)$. Then, we want to determine
$\mathbf{m}$ to maximize the log-likelihood:

$\displaystyle\log p(\mathcal{S}|\mathbf{m})$
$\displaystyle=\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}$
$\displaystyle=-\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 [1] 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}\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_{\mathbf{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}}$
(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 [5]:

$\mathbf{B}_{m}\mathbf{m}=\mathbf{C}_{m}\mathbf{m},$ (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}_{\mathbf{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}$ [1]:

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

and

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

##  4 Second Order Velocity Field

###  4.1 Definition

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 $\mathbf{p}$:

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

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

#### Convergence Point

There is no closed-form solution for isolating convergence points
$\mathbf{p}_{0}$ 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 ($\|\mathbf{v}(\mathbf{p}_{i})\|$). In analogy to Eq. 3,
we can find an alternative expression of the velocity field centered at
$\mathbf{p}_{0}$:

$\mathbf{v}(\mathbf{p})=\Big{(}\mathbf{t}\times(\mathbf{p}-\mathbf{p}_{0})+\mathbf{r}\Big{)}\times(\mathbf{p}-\mathbf{p}_{0})+\gamma(\mathbf{p}-\mathbf{p}_{0}).$
(11)

#### Streamlines

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

#### Core Line

When $\mathbf{t}\neq 0$, no straight axis of rotational symmetry can be found,
but we can extract core lines as defined in section2.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 [21, 30].

###  4.2 Kinematic Surface Fitting

To estimate the 10-dimensional vector
$\mathbf{m}=[\mathbf{t},\mathbf{r},\mathbf{c},\gamma]$ 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 $\mathbf{m}$. For
instance, the 10-dimensional vector $\mathbf{f}(\mathbf{p},\mathbf{n})$ now
writes as
$[(\mathbf{n}\times\mathbf{p})\times\mathbf{p},\mathbf{p}\times\mathbf{n},\mathbf{n},\mathbf{p}\cdot\mathbf{n}]$
and the matrix $\mathbf{M}_{i}$ can be obtained using Eq. 8.

To compute the normalization matrix
$\mathbf{N}_{i}=\mathbf{N}^{(1)}_{i}+w_{\mathbf{p}}\mathbf{N}^{(2)}_{i}$, we
first express the term $\|\mathbf{v}(\mathbf{p}_{i})\|^{2}$ as a
multiplication of the matrix $\mathbf{N}_{i}^{(1)}$ and $\mathbf{m}$ where :

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

with
$\mathbf{F}_{\mathbf{p}_{i}}=\mathbf{p}_{i}\otimes\mathbf{p}_{i}-(\mathbf{p}_{i}\cdot\mathbf{p}_{i})\mathbf{I}$,
where $\otimes$ denotes the tensor product such that
$\mathbf{a}\otimes\mathbf{b}=\mathbf{a}\mathbf{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}+(\mathbf{t}\cdot\mathbf{p})\mathbf{n}-2(\mathbf{p}\otimes\mathbf{t})\mathbf{n}+\mathbf{n}\times\mathbf{r}+\gamma\mathbf{n}.$
(12)

Then, we can express
$\|\nabla_{\mathbf{p}}(\mathbf{v}(\mathbf{p}_{i})\cdot\mathbf{n}_{i})\|^{2}$
in terms of the matrix $\mathbf{N}_{i}^{(2)}$ and $\mathbf{m}$:

$\mathbf{N}^{(2)}_{i}=\begin{bmatrix}\mathbf{G}_{\mathbf{p}_{i},\mathbf{n}_{i}}^{T}\mathbf{G}_{\mathbf{p}_{i},\mathbf{n}_{i}}&\mathbf{G}_{\mathbf{p}_{i},\mathbf{n}_{i}}\mathbf{A}_{\mathbf{n}_{i}}&\mathbf{0}&\mathbf{G}_{\mathbf{p}_{i},\mathbf{n}_{i}}\mathbf{n}_{i}\\\
\mathbf{A}_{\mathbf{n}_{i}}^{T}\mathbf{G}_{\mathbf{p}_{i},\mathbf{n}_{i}}^{T}&\mathbf{A}_{\mathbf{n}_{i}}^{T}\mathbf{A}_{\mathbf{n}_{i}}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{n}_{i}^{T}\mathbf{G}_{\mathbf{p}_{i},\mathbf{n}_{i}}^{T}&\mathbf{0}&\mathbf{0}&\mathbf{n}_{i}\cdot\mathbf{n}_{i}=1\end{bmatrix},$

with
$\mathbf{G}_{\mathbf{p}_{i},\mathbf{n}_{i}}=(\mathbf{p}_{i}\cdot\mathbf{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.

##  5 Robust Kinematic Surface Fitting

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[13] need to estimate precisely extra variables such
as the percentage of outliers.

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.

More precisely, we now assume that
$p(d_{i}|\mathbf{m})={\mathrm{St}}(d_{i}|0,\Sigma,\nu)=\\\
\int_{z_{i}}\mathcal{N}(d_{i}|0,\Sigma/z_{i})\,\mathrm{Ga}(z_{i}|\nu/2,\nu/2)\mathrm{d}z_{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\rightarrow+\infty$, the Student-t distribution is equivalent to the
Gaussian distribution and the variable $\nu$ is inversely proportional to the
number of outliers.

The robust estimation of a kinematic surface is now achieved with an
Expectation-Maximization scheme, where $z_{i}$ is the latent variable [23],
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)=\mathrm{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.

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:

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

where $\psi(x)$ denotes the digamma function. The value of $\nu$ after
performing a few iterations using an implementation of the Levenberg-Marquardt
algorithm [14].

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}^{T}\mathbf{N}_{i}\mathbf{m}}$
and
$\mathbf{C}_{m}=\sum_{i=1}^{n}z_{i}\frac{\mathbf{m}^{T}\mathbf{M}_{i}\mathbf{m}}{(\mathbf{m}^{T}\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.

##  6 Results

###  6.1 Validation

Figures 1 and 2 illustrate examples of non-robust and robust kinematic surface
fitting in synthetic test shapes. The fitting results are listed in Table 1.
In the first example, a logarithmic helical spiral, the fit second order
velocity field becomes degenerate, and is identical to the first order field
solution (Fig. 1a) as $\mathbf{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 a curved helix. It is obvious that the first order
velocity field is incapable of capturing the bent core line of the structure,
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).

Table 1: Validation of fitting on test surfaces (15 iterations,
$w_{\mathbf{p}}=0.001$).

First order rotational scale velocity field Second order rotational scale
velocity field Shape $\nu$ RMSE $\nu$ RMSE Logarithmic spiral, no outliers
(Fig. 1a) - 2.2$\cdot 10^{-3}$ - 2.2$\cdot 10^{-3}$ \- with outliers, non-
robust fitting (Fig. 1b) - 9.5$\cdot 10^{-2}$ - 1.9$\cdot 10^{-1}$ \- with
outliers, robust fitting (Fig. 1c) 0.96 9.5$\cdot 10^{-2}$ 0.99 9.5$\cdot
10^{-2}$ Bent helix, no outliers (Fig. 2a) - 2.8$\cdot 10^{-1}$ - 3.9$\cdot
10^{-2}$ \- with outliers, non-robust fitting (Fig. 2b) - 3.2$\cdot 10^{-1}$ -
2.0$\cdot 10^{-1}$ \- with outliers, robust fitting (Fig. 2c) 1.75 3.2$\cdot
10^{-1}$ 1.30 1.6$\cdot 10^{-1}$

abcd

Figure 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 $\mathbf{p}_{0}$ are shown as
circles. Continuous curves represent example center lines $\mathcal{C}$ traced
from a central seed point $\mathbf{p}_{\mathrm{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.

abcd

Figure 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. Continuous curves represent example center lines $\mathcal{C}$ traced
from a central seed point $\mathbf{p}_{\mathrm{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.

###  6.2 Center Line and Core Line Extraction

Figures 3 to 7 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 $\mathcal{C}$ of the
best fitting fields traced from selected central seed points
$\mathbf{p}_{\mathrm{s}}$ (black diamonds). The seed points were found by
computing the centroid of vertices of a cross-section of the surface aligned
with the rotational axis [11]. The parameters found for both velocity fields
are summarized in Table 2.

Table 2: Kinematic surface fitting parameters for anatomical structures (15
iterations, $w_{\mathbf{p}}=0.001$).

First order rotational scale velocity field Second order rotational scale
velocity field Anatomical structure $\mathbf{r}$ $\mathbf{c}$ $\gamma$ $\nu$
RMSE $\mathbf{t}$ $\mathbf{r}$ $\mathbf{c}$ $\gamma$ $\nu$ RMSE Aortic arch
(Fig. 3) $[0.83,-0.43,-0.27]^{T}$ $[0.03,0.04,0.11]^{T}$ -0.20 4.68 2.2$\cdot
10^{-1}$ $[0.05,0.76,0.46]^{T}$ $[0.36,-0.17,-0.16]^{T}$
$[-0.04,0.07,0.06]^{T}$ -0.09 5.52 1.9$\cdot 10^{-1}$ Left ventricle (Fig. 4)
$[0.02,-0.16,0.99]^{T}$ $[-0.02,0.008,0.003]^{T}$ -0.002 6.89 1.2$\cdot
10^{-1}$ $[0.24,-0.18,-0.03]^{T}$ $[0.03,-0.14,0.98]^{T}$
$[0.001,-0.005,0.0004]^{T}$ -0.002 8.77 1.1$\cdot 10^{-1}$ Rib (Fig. 5)
$[0.12,-0.11,-0.72]^{T}$ $[0.06,-0.33,0.24]^{T}$ -0.53 3.17 1.1$\cdot 10^{-1}$
$[-0.37,-0.70,0.20]^{T}$ $[0.06,-0.16,-0.43]^{T}$ $[0.02,-0.17,0.09]^{T}$
-0.27 2.74 8.3$\cdot 10^{-2}$ Cochlea (Fig. 6) $[-0.14,-0.09,0.97]^{T}$
$[-0.09,0.13,0.03]^{T}$ -0.08 2.51 1.8$\cdot 10^{-1}$ $[-0.46,0.23,0.17]^{T}$
$[-0.14,-0.03,0.81]^{T}$ $[-0.10,0.12,0.02]^{T}$ -0.13 2.19 1.6$\cdot 10^{-1}$
Semicircular canal (Fig. 7) $[-0.53,0.76,-0.36]^{T}$ $[-0.03,-0.13,0.01]^{T}$
-0.005 1.84 1.5$\cdot 10^{-1}$ $[-0.40,0.43,0.39]^{T}$
$[-0.47,0.42,-0.31]^{T}$ $[-0.08,-0.04,-0.04]^{T}$ -0.03 1.80 1.7$\cdot
10^{-1}$

Figure 3: Kinematic surface fitting of a human aortic arch [17]. 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 $\mathbf{p}_{0}$ are shown as
circles. Continuous curves represent example center lines $\mathcal{C}$ traced
from a central seed point $\mathbf{p}_{\mathrm{s}}$ (black diamonds).

Figure 4: Kinematic surface fitting of a human left ventricle [4]. 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 $\mathbf{p}_{0}$ are shown as
circles. Continuous curves represent example center lines $\mathcal{C}$ traced
from a central seed point $\mathbf{p}_{\mathrm{s}}$ (black diamonds).

Figure 5: Kinematic surface fitting of a human rib [18]. 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 $\mathbf{p}_{0}$ are shown as
circles. Continuous curves represent example center lines $\mathcal{C}$ traced
from a central seed point $\mathbf{p}_{\mathrm{s}}$ (black diamonds)

Figure 6: Example of a right human cochlea [26]. 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
$\mathbf{p}_{0}$ are shown as circles. Continuous curves represent example
center lines $\mathcal{C}$ traced from a central seed point
$\mathbf{p}_{\mathrm{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.

Figure 7: Kinematic surface fitting of a human semicircular canal [33].
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 $\mathbf{p}_{0}$ are
shown as circles. Continuous curves represent example center lines
$\mathcal{C}$ traced from a central seed point $\mathbf{p}_{\mathrm{s}}$
(black diamonds).

###  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. 6. 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. [25] (see Fig. 8). As an intrinsic shape parameter, we calculate the
magnitude of the projection of $\mathbf{t}$ on $\mathbf{r}$, that is,
$\|\mathrm{proj}_{\mathbf{r}}\mathbf{t}\|=\frac{\mathbf{t}\cdot\mathbf{r}}{\|\mathbf{r}\|}$.
The correlation between the parameter
$\|\mathrm{proj}_{\mathbf{r}}\mathbf{t}\|$ and the measured tilt angle is
demonstrated in a set of 31 human cochleae [33, 26]. This association could be
used to classify a cochlea in terms of coiling and tilting characteristics,
i.e., ’sloping’ vs. ’rollercoaster’ shapes [3].

Figure 8: Relationship between the intrinsic shape parameter
$\|\mathrm{proj}_{\mathbf{r}}\mathbf{t}\|$ and the measured tilt angle
$\theta$ [25] for 31 samples [33, 26]. The linear regression line and the
coefficient of determination are shown.

##  7 Discussion

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 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. 6), 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. 7). Similarly, a
first order rotational field can be used as an approximation for symmetry
detection, e.g., in a heart ventricle (Fig. 4) or for modiolar axis detection
in cochleae as shown by Wimmer et al. [34]).

It is noteworthy that the recognition of revolutional 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 [7]. 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 [8, 9,
28].

Additional applications lie in segmentation and surface generation, where
velocity fields can be used to trace out kinematic surfaces using generator
curves [11]. It is further possible to explore how useful our approach is for
modeling the biological growth of tissue, e.g., bones [36].

The global kinematic characterization of shapes may be desirable for the
classification of suitable structures, e.g., the cochlea (Fig. 6). 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.

##  8 Conclusions

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.

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.

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. 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. 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.

## CRediT Authorship Contribution Statement

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.

## Acknowledgments

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.

## References

  * Andrews and Séquin [2013] Andrews, J., Séquin, C.H., 2013.  Generalized, basis-independent kinematic surface fitting.  Computer-Aided Design 45, 615–620.
  * Arsigny et al. [2009] 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 et al. [2014] 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 et al. [2015] 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 [2007] Chernov, N., 2007.  On the convergence of fitting algorithms in computer vision.  Journal of Mathematical Imaging and Vision 27, 231–239.
  * Eckhoff et al. [2016] 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.
  * Günther and Theisel [2018] 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 et al. [2023a] Harte, N., Obrist, D., Caversaccio, M., Lajoinie, G.P., Wimmer, W., 2023a.  Transverse flow under oscillating stimulation in helical square ducts with cochlea-like geometrical curvature and torsion.  arXiv preprint arXiv:2303.15603 .
  * Harte et al. [2023b] Harte, N.C., Obrist, D., Caversaccio, M.D., Lajoinie, G.P., Wimmer, W., 2023b.  Wall Shear Stress and Pressure Fluctuations under Oscillating Stimulation in Helical Square Ducts with Cochlea-like Geometrical Curvature and Torsion, in: 2023 45th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), IEEE, Sydney, Australia. pp. 1–7.  URL: https://ieeexplore.ieee.org/document/10340844/, doi:10.1109/EMBC40787.2023.10340844.
  * Helman and Hesselink [1989] 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.
  * Hofer et al. [2005] 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.
  * Jayasuriya et al. [2013] 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.
  * Jung [2007] 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.
  * Liu and Rubin [1995] Liu, C., Rubin, D.B., 1995.  Ml estimation of the t distribution using em and its extensions, ecm and ecme.  Statistica Sinica , 19–39.
  * Liu [2009] 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.
  * Merzkirch [2012] Merzkirch, W., 2012.  Flow visualization.  Elsevier.
  * Online [2023a] Online, 2023a.  Aorta model - ultimaker thingiverse.  https://www.thingiverse.com/thing:2277464.  Accessed on 2023-07-26.
  * Online [2023b] Online, 2023b.  Rib cage model - nih 3d database.  https://3d.nih.gov/entries/3dpx-016836.  Accessed on 2023-07-26.
  * Post et al. [2003] 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.
  * Pottmann and Wallner [2006] Pottmann, B.O.H., Wallner, J., 2006.  Equiform kinematics and the geometry of line elements.  Contributions to Algebra and Geometry 47, 567–582.
  * Roth and Peikert [1998] 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.
  * Rueckert et al. [1997] 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.
  * Scheffler [2008] 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 .
  * Seo et al. [2022] 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.  doi:10.1109/CVPR52688.2022.00932.
  * Shin et al. [2013] 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.
  * Sieber et al. [2019] 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.
  * Sironi et al. [2015] Sironi, A., Türetken, E., Lepetit, V., Fua, P., 2015.  Multiscale centerline detection.  IEEE Transactions on Pattern Analysis and Machine Intelligence 38, 1327–1341.
  * Stokes et al. [2023] 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.
  * Thirion et al. [2000] 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.
  * Van Gelder and Pang [2009] 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.
  * van Stralen et al. [2008] 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.  doi:https://doi.org/10.1016/j.ultrasmedbio.2007.07.016.
  * Vercauteren et al. [2008] 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.
  * Wimmer et al. [2019a] 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.
  * Wimmer et al. [2019b] 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.
  * Wu et al. [2021] 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.
  * Xu et al. [2014] 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.


Thereby, the phase convention dependence is cancelled in $\epsilon_{K}$.

The simplest way $\epsilon_{K}$ can be related to the Hamiltonian off diagonal
elements is as follows BurasGuadagnoliIsidori2010

$\epsilon_{K}=\frac{\kappa_{\epsilon}e^{i\phi_{\epsilon}}}{\sqrt{2}\Delta
M_{K}}\text{Im}[M_{12}^{\rm SD}]\ ,$ (90)

where $\phi_{\epsilon}=(43.51\pm 0.05)^{o}$ and the phenomenological factor,
$\kappa_{\epsilon}=0.94\pm 0.02$ BurasGirrbach014 , summarizes the corrections
to $\epsilon_{K}$ due to long distance (LD) effects BurasGuadagnoliIsidori2010
. $M_{12}^{SD}$ is the short distance (SD) contribution to $M_{12}$ defined in
80. This SD contribution is calculated from the box diagrams of the
$K^{0}-\bar{K}^{0}$ mixing.202020The box diagrams in DQCD are those of Figure
4. They are similar in the SM, but replacing $X$ by the $W$ boson and $Q$ by
the up-type quarks. Equation (90) clearly shows how $\epsilon_{K}$ represents
$CP$ violation due to mixing ($\epsilon_{K}$ vanishes if $M_{12}^{\rm SD}=0$).
The SD matrix element is given by

$M_{12}^{\rm SD}=\frac{1}{2m_{K}}\bra{\bar{K}^{0}}{H_{\Delta S=2}^{\rm
Box}}\ket{K^{0}}^{\ast}.$ (91)

Here, $H_{\Delta S=2}^{\rm Box}$ is the effective Hamiltonian calculated from
the box diagrams of the $K^{0}-\bar{K}^{0}$ mixing. Finally, for the case of
the $K$ mesons only, since $\epsilon_{K}$ is very small,
$\mathcal{O}(10^{-3})$, the imaginary parts of $M_{12}$ and $\Gamma_{12}$ are
much smaller than their real parts, thus, from equation (84), we have the
approximation

$\Delta M_{K}\simeq 2\text{Re}\left[M_{12}\right]\ ,$ (92)

where $M_{12}$ contains both SD and LD contributions.

###  A.2 The $B-$Systems

We now explore the basic meson mixing formalism for $B_{d}$ and $B_{s}$. The
equations are the same as equations (80)$-$(88) developed for $K$ mesons,
except that the mass eigenstates are called $B_{\text{L}}$ and $B_{\text{H}}$
for the light and heavy state, respectively. Since $\Gamma_{12}\ll M_{12}$ for
the $B$-system BurasEtAlBBbar , we have, from equation (84) Buras2005 ,

$\Delta M_{B}=M_{H}-M_{L}\approx 2\mid M_{12}^{B}\mid\ .$ (93)

Like equation (91), we have for the B$-$system

$M_{12}^{B}=\frac{1}{2m_{B}}\bra{\bar{B}^{0}}{H_{\Delta
B=2}}\ket{B^{0}}^{\ast}$ (94)

Because the mass of the $B$ mesons is significantly above the QCD scale, the
LD effects are negligible in the calculation of $\Delta M_{B}$ Buras1998 .

Now expanding equation (82), into powers of $|\Gamma_{12}/M_{12}|$, we can
keep the leading term only, since $|\Gamma_{12}/M_{12}|$ is estimated to be of
order $10^{-4}$ Buras2005 ; NirCPV , giving

$\left(\frac{q}{p}\right)_{B}\approx\frac{M_{12}^{\ast}}{\mid M_{12}\mid}\ \ \
\stackrel{{\scriptstyle\mathclap{\mbox{\text{(\ref{DeltaMB})}}}}}{{\simeq}}\ \
\ \frac{2M_{12}^{\ast}}{\Delta M_{B}}\ .$ (95)

Thus $(q/p)_{B}$ is a pure phase to an excellent approximation, which means

$\left(\frac{q}{p}\right)_{B}=e^{2i\phi_{M}}\simeq\ \ \
\frac{2M_{12}^{\ast}}{\Delta M_{B}}\ ,$ (96)

where $\phi_{M}$ is a CPV phase due to mixing. By convention,
$\phi_{M}=-\beta$ for $B^{0}_{d}$, and $\phi_{M}=-\beta^{0}_{s}$ for $B_{s}$.

$\left(\frac{q}{p}\right)_{B_{d}}=e^{-2i\beta}\approx\frac{2M_{12}^{\ast}(B_{d})}{\Delta
M_{d}}\ ,$ (97)
$\left(\frac{q}{p}\right)_{B_{s}}=e^{-2i\beta_{s}}\approx\frac{2M_{12}^{\ast}(B_{s})}{\Delta
M_{s}}\ .$ (98)

Within the Standard Model, the contribution to the amplitudes of the
$B-\bar{B}$ mixing are dominated by the box diagrams with two $W$ propagators
and 2 internal $t$ quarks, so the SM weak phase of mixing comes dominantly
from the arguments of the CKM matrix elements. For $B_{d}^{0}$ this comes from
$V_{td}$:

$\text{Arg}\left(\frac{q}{p}(\text{SM})\right)_{B_{d}}=2\text{Arg}(V_{td})=-2\beta_{\text{SM}}\
,$ (99)

and for $B_{s}^{0}$ this comes from $V_{ts}$:

$\text{Arg}\left(\frac{q}{p}(\text{SM})\right)_{B_{s}}=2\text{Arg}(V_{ts})=-2(\beta_{s})_{\text{SM}}\
.$ (100)

Both $\beta_{\text{SM}}$ and $(\beta_{s})_{\text{SM}}$ are related to the
unique phase of the CKM matrix. Their predicted SM values are reported in
Table 4 from UTfit collaboration Bona2016Utfit2 . This same reference provides
the experimental measurements of these angles, reported in Table 4 as well.

##  Appendix B Mixing Induced $CP-$Asymmetry

In what follows, we will briefly introduce a very interesting class of
$CP-$asymmetries which result from $CP-$violation in the interference between
mixing and decay Buras2005 ; NirFPandCPV ; NirCPV ; CPasymExp  which can be
defined for neutral meson decays. Let $B^{0}$ be any neutral meson; knowing
that it oscillates in time between the flavour eigenstates $B^{0}$ and
$\bar{B}^{0}$, we define $B^{0}(t)$ as the meson at a time $t$, starting as a
$B^{0}$ at $t=0$. Similarly, $\bar{B}^{0}(t)$ is defined such as
$\bar{B}^{0}(0)=\bar{B}^{0}$, the $CP-$conjugate of $B^{0}$. In a decay
$B^{0}\rightarrow f$, where $f$ is a $CP-$eigenstate, the time-dependent
mixing induced $CP-$asymmetry is defined as:

$\mathcal{A}_{CP}(t)=\frac{\Gamma(B^{0}(t)\rightarrow
f)-\Gamma(\bar{B}^{0}(t)\rightarrow f)}{\Gamma(B^{0}(t)\rightarrow
f)+\Gamma(\bar{B}^{0}(t)\rightarrow f)}\ .$ (101)

This quantity can be separated into two terms, a decay only contribution with
coefficient $C_{f}$, and a term due to the interference with coefficient
$S_{f}$:

$\mathcal{A}_{CP}(t)=C_{f}\cos(\Delta M_{B}t)-S_{f}\sin(\Delta M_{B}t)\ ,$
(102)

where $\Delta M_{B}$ is the mass difference between the two $B^{0}$ mass
eigenstates. The expressions of the decay and mixing contributions can be
calculated and found to be Buras2005 ; NirFPandCPV ; NirCPV :

$\displaystyle C_{f}=$ $\displaystyle\frac{1-|\xi_{f}|^{2}}{1+|\xi_{f}|^{2}}\
,$ (103) $\displaystyle S_{f}=$
$\displaystyle\frac{-2\text{Im}(\xi_{f})}{1+|\xi_{f}|^{2}}\ ,$
$\displaystyle\xi_{f}=$
$\displaystyle\frac{q}{p}\frac{A(\bar{B}^{0}\rightarrow f)}{A(B^{0}\rightarrow
f)}\ ,$

where $A$ is the decay amplitude of the flavour eigenstate at $t=0$, and $p$
and $q$ are the mixing parameters defined in Appendix A.

In general, $A(B^{0}\rightarrow f)$ has different terms coming from different
Feynman diagrams. For example, if only two diagrams contribute to the decay
then

$\frac{A(\bar{B}^{0}\rightarrow f)}{A(B^{0}\rightarrow f)}=-\ \zeta_{f}\
\frac{A_{1}e^{i(\delta_{1}-\phi_{1})}+A_{2}e^{i(\delta_{2}-\phi_{2})}}{A_{1}e^{i(\delta_{1}+\phi_{1})}+A_{2}e^{i(\delta_{2}+\phi_{2})}}\
,$ (104)

where $\phi_{i}$ are the CPV (or the weak) phases, $\delta_{i}$ are the strong
phases,212121The strong phases are complex phases that are not present in the
Lagrangian. They appear in decay amplitudes when intermediate states go on-
shell via $CP-$conserving processes, so they do not change sign under $CP$.
They are called “strong phases” because they are usually due to strong
interactions. $\zeta_{f}$ is the $CP-$parity of the decay product $f$, and the
$A_{i}$ contain the hadronic matrix elements. If $\phi_{1}\neq\phi_{2}$, and
the two diagrams contribute with similar strengths, then the asymmetry will
suffer from hadronic uncertainties due to the presence of the strong phases
$\delta_{i}$, and the hadronic matrix elements in $A_{i}$. However, a very
interesting scenario happens when the decay is either dominated by one
diagram, or all diagrams have the same CPV phase
($\phi_{1}=\phi_{2}\equiv\phi_{D}$ in the example). In these cases, the strong
phases as well as the hadronic matrix elements cancel completely from
$A(\bar{B}^{0}\rightarrow f)/A(B^{0}\rightarrow f)$, and we are left with:

$\left[\frac{A(\bar{B}^{0}\rightarrow f)}{A(B^{0}\rightarrow
f)}\right]_{\text{sp}}=-\zeta_{f}e^{-i2\phi_{D}}\ ,$ (105)

where we added the index “$D$” to the weak phase $\phi$ to express its decay
origin, and we added the index “sp” to show that we are in the special case
where either one diagram dominates the decay or all diagrams have the same CPV
phase. In this special case, where equation (105) applies, since $(q/p)$ is
almost a pure phase for both $B^{0}_{d}$ and $B^{0}_{s}$ (equations (97, 98)),
then we have $|\xi_{f}|^{2}=1$ and thus, $C_{f}=0$ and
$S_{f}=-\zeta_{f}\sin(2\phi_{D}-2\phi_{M})$. Thus we get from equations (102)
and (103):

$[\mathcal{A}_{CP}(t)]_{\text{sp}}=\zeta_{f}\sin(2\phi_{D}-2\phi_{M})\sin(\Delta
M_{B}t)$ (106)

We note that the difference between two weak phases, $(\phi_{D}-\phi_{M})$, is
convention-independent.

The decays $(B_{d}\rightarrow\psi K_{\text{S}})$ and
$(B_{s}\rightarrow\psi\phi)$, have the same inclusive mode $b\rightarrow
sc\bar{c}$. In the SM, it is dominated by the tree diagram with an internal
$W$ boson. In NP, our DQCD model’s contribution starts only at the 2-loop
level, so it is negligible. Therefore, we are in the special case of one
dominant diagram, where equation (106) applies. For the considered decays,
$\phi_{D}$ is the weak phase of the SM tree diagram, that is,
$\phi_{D}=\text{Arg}(V_{cs}V_{cb}^{\ast})\simeq 0$. So the asymmetries in
these decays depend only on the weak phase of the mixing $\phi_{M}$, which is
$(-\beta)$ for $B_{d}$ and $(-\beta_{s})$ for $B_{s}$. Therefore

$\displaystyle\mathcal{A}_{CP}(B_{d}\rightarrow\psi K_{\text{S}})$
$\displaystyle=\zeta_{\psi K_{\text{S}}}\sin(2\beta)\sin[(\Delta M_{d})t]$
(107) $\displaystyle\mathcal{A}_{CP}(B_{s}\rightarrow\psi\phi)$
$\displaystyle=\zeta_{\psi\phi}\sin(2\beta_{s})\sin[(\Delta M_{s})t]\ ,$

where $\Delta M_{d(s)}$ is the mass difference between the two $B_{d(s)}$
states, and the $CP-$parity $\zeta_{\psi K_{\text{S}}}=-1$. The final state
$(\psi\phi)$ is a mixture of $CP-$even and $CP-$odd states, nevertheless, its
$CP-$parity can be resolved experimentally LinnPhD .

##  Appendix C Calculation of The $\Delta F=2$ Constraints

In this section we show the explicit calculation of the $\Delta F=2$
constraints on the DQCD scenario that are shown in Table 1. In the
calculations of the new physics (NP) contributions to the flavor observables,
the off-diagonal elements of the mixing mass matrix (see equations (91) and
(94)) are involved. For all three mesons ($M=K,\ B_{d},\ B_{s}$), they can be
written as

$M_{12}^{\rm NP}=\frac{1}{2m_{M}}\bra{\bar{M}^{0}}H_{\rm eff}^{\rm NP}(\Delta
F=2)\ket{M^{0}}^{\ast}=\frac{N_{c}}{384{\pi}^{2}M_{X}^{2}}m_{M}F_{M}^{2}\hat{B}_{M}\eta_{M}(\xi_{M}^{\ast})^{2}\
,$ (108)

where $m_{M}$ is the experimental value of the meson’s mass, $F_{M}$ is its
decay constant, $\hat{B}_{M}$ is its bag parameter, and $\eta_{M}$ includes
any RG and QCD corrections, and DQCD non-perturbative uncertainties, which, as
mentioned in Section 4, we take to be $1\pm 0.5$.

Table 4 contains all the experimental input data we use as well as the derived
$\Delta F=2$ constraints, along with the corresponding references. These
constraints will impose limits on the parameters of our model, whether the
real ones, or the $CP-$violating phases, that we will abbreviate to “CPV
phases.” The $C_{M}$ parameters in Table 4 are defined by

$C_{M}=\Delta M_{M}^{\rm exp}/\Delta M_{M}^{\rm SM}\ ,\ M=K,B_{d},B_{s}\ ,$
(109)

that is, the ratio of the experimental mass difference to the prediction from
the SM.

Experimental values PDG  | Lattice QCD values
---|---
$m_{K}=497.611\pm 0.013\text{ MeV}$ |  $F_{K}=155.8\pm 1.7\text{ MeV}$ QCDlatticeUT
$m_{B_{d}}=5279.62\pm 0.15\text{ MeV}$ |  $\hat{B}_{K}=0.7625\pm 0.0097$ FLAG016
$m_{B_{s}}=5366.82\pm 0.22\text{ MeV}$ |  $F_{B_{d}}\sqrt{\hat{B}_{B_{d}}}=219\pm 14\text{ MeV}$ FLAG016
$\Delta M_{K}=(3.484\pm 0.006)\times 10^{-12}\text{ MeV}$ |  $F_{B_{s}}\sqrt{\hat{B}_{B_{s}}}=270\pm 16\text{ MeV}$ FLAG016
$\Delta M_{d}=(3.354\pm 0.022)\times 10^{-10}\text{ MeV}$ |  $F_{B_{d}}=192.0\pm 4.3\text{ MeV}$ FLAG016
$\Delta M_{s}=(1.1688\pm 0.0014)\times 10^{-8}\text{ MeV}$ |  $F_{B_{s}}=228.4\pm 3.7\text{ MeV}$ FLAG016
$\Gamma_{K_{\text{L}}}=(1.286\pm 0.005)\times 10^{-17}$ GeV |  $\Delta F=2$ constraints
$\Gamma_{K^{+}}=(5.315\pm 0.009)\times 10^{-17}$ GeV |  $C_{K}\in[0.51,2.07]$ UTFit2008
$\Gamma_{B_{s}}=(4.358\pm 0.014)\times 10^{-13}$ GeV |  $C_{B_{d}}\in[0.77,1.35]$ UTfit2016
$\Gamma_{B_{d}}=(4.329\pm 0.011)\times 10^{-13}$ GeV |  $C_{B_{s}}\in[0.87,1.30]$ UTfit2016
$|\epsilon_{K}|=(2.228\pm 0.011)\times 10^{-3}$ |  $10^{3}|\epsilon_{K}|^{\rm SM}\in[1.76,2.33]$ BurasGirrbach014
$\sin 2\beta_{d}\in[0.634,0.726]$ Bona2016Utfit2  |  $\sin 2\beta_{d}^{\rm SM}\in[0.665,0.785]$ Bona2016Utfit2
$\beta_{s}[^{o}]\in[-0.91,2.85]$ Bona2016Utfit2  |  $\beta_{s}^{\rm SM}[^{o}]\in[0.97,1.13]$ Bona2016Utfit2
Table 4: Summary of quantities and constraints used. All values are taken at
95% confidence level (CL).

###  C.1 Meson Mass Differences

We show in this section how we get the constraints on the $\xi_{M}$ terms, due
to the mass difference data, for the $K$, $B$, and $B_{s}$ mesons.

The $CP$ conserving constraint, at the 2$\sigma$ level, from the $K-\bar{K}$
mixing is UTFit2008

$0.51\leq\ C_{K}\equiv\Delta M_{K}^{\rm exp}/\Delta M_{K}^{\rm SM}\ \leq 2.07\
,$ (110)

from which we can derive the allowed interval for the new physics
contribution:

$\Delta M_{K}^{\rm NP}=\Delta M_{K}^{\rm exp}-\Delta M_{K}^{\rm
SM}\pm\sqrt{\sigma_{\rm exp}^{2}+\sigma_{SM}^{2}}\ .$ (111)

From equation (92), we have the relation

$\Delta M_{K}^{\rm NP}=2{\rm Re}(M_{12}^{\rm NP})\ ,$ (112)

then we use equation (108) for $M_{12}^{\rm NP}$. Using the values given in
Table 4, we can constrain $\Delta M_{K}^{NP}$, and obtain:

$-6.93\times 10^{-4}\leq
N_{c}\text{Re}\left[({\xi}_{K}^{*})^{2}\right]\left(\frac{1\text{
TeV}}{M_{X}}\right)^{2}\leq 3.77\times 10^{-4},$ (113)

where ${\xi}_{K}$ is defined in equation (7), and $N_{c}$ is the number of
dark colors.

Moving to the $B_{d}$ system, the $B_{d}-\bar{B}_{d}$ mixing mass difference
constraint, at the 2$\sigma$ level, is UTfit2016

$0.77\leq C_{B_{d}}\leq 1.35\ .$ (114)

From equation (93), the relation between the new physics contribution to
$\Delta M_{B_{d(s)}}$ and $M_{12}^{NP}$, the NP’s contribution to $M_{12}$, is
given by:

$\Delta M_{B_{d(s)}}^{\rm NP}=2\mid M_{12}^{\rm
NP}\mid\,=\,\mid\bra{\bar{B}_{d(s)}^{0}}H_{\rm eff}^{\rm NP}(\Delta
B=2)\ket{B_{d(s)}^{0}}\mid/m_{B_{d(s)}}.$ (115)

Using the same strategy as in the $K-$system, we find

$N_{c}\mid{\xi}_{B_{d}}^{2}\mid\left(\frac{1\text{ TeV}}{M_{X}}\right)^{2}\leq
6.55\times 10^{-4}.$ (116)

Similarly, for the $B_{s}$ system the bound is given by

$0.87\leq C_{B_{s}}\leq 1.3\ .$ (117)

which then leads to

$N_{c}\mid{\xi}_{B_{s}}^{2}\mid\left(\frac{1\text{ TeV}}{M_{X}}\right)^{2}\leq
13.15\times 10^{-3}.$ (118)

###  C.2 The $\Delta F=2\ CP-$Violating Constraints

Now we briefly show how we get the constraints from the $CP$ violating
observables in the $K-\bar{K}$, $B_{d}-\bar{B}_{d}$, and $B_{s}-\bar{B}_{s}$
mixings.

We start with the $K-\bar{K}$ mixing, where we consider the $CP$ violating
observable $\mid\epsilon_{K}\mid$ which is the modulus of equation (90). From
the SM prediction at 95% CL BurasGirrbach014

$1.76\times 10^{-3}\leq|\epsilon_{K}^{\rm SM}|\leq 2.33\times 10^{-3}\ ,$
(119)

and the experimental value PDG  is $|\epsilon_{K}^{\rm exp}|=2.228(11)\times
10^{-3}$. We evaluate the constraint on the new physics contribution to
$|\epsilon_{K}|$ at 95% CL is

$\mid\epsilon_{K}\mid_{\rm NP}\ \leq\ 0.75\times 10^{-3}.$ (120)

Taking the modulus of $\epsilon_{K}$ in equation (90) and using equation
(108), we find

$\displaystyle\mid\epsilon_{K}\mid_{\rm NP}$
$\displaystyle\approx\frac{\kappa_{\epsilon}}{\sqrt{2}(\Delta M_{K})_{\rm
exp}}\mid\text{Im}[M_{12}^{\rm NP}]\mid$ (121)
$\displaystyle=\frac{\kappa_{\epsilon}m_{K}F_{K}^{2}\hat{B}_{K}\eta_{K}}{384\sqrt{2}{\pi}^{2}(\Delta
M_{K})_{\rm
exp}M_{X}^{2}}N_{c}\mid\text{Im}\left[(\xi_{K}^{\ast})^{2}\right]\mid\ ,$
$\displaystyle\ \ \kappa_{\epsilon}=0.94\pm 0.02\ .$

With the data of Table 4 we then get

$N_{c}\mid\text{Im}\left[({\xi}_{K}^{*})^{2}\right]\mid\left(\frac{1\text{
TeV}}{M_{X}}\right)^{2}\leq 1.64\times 10^{-6}.$ (122)

Let us now move to the constraints from the $CP-$violating observables in the
$B-$ system, that is, $\beta$ and $\beta_{s}$. These parameters can be
extracted from global fits HFAG  and are also measured directly from the
mixing-induced $CP$ asymmetries (defined in the previous section) in the
decays of $B$ mesons. This allows one to place a constraint on new physics
contributing to these processes. For example, the asymmetry in
$(B_{d}\rightarrow\psi K_{\text{S}})$ measures $\beta$ and the asymmetry in
$(B_{s}\rightarrow\psi\phi)$ measures $\beta_{s}$, as shown in equation (107)
. These decay modes provide “clean” measurements of those phases, because the
hadronic matrix elements and the strong phases drop out from their
$CP-$asymmetries, so the measurements of the CPV phases are free from hadronic
uncertainties. From equations (97)-(100), we can write the relations

$S_{\psi K_{\text{S}}}=\sin 2\beta_{\text{exp}}\approx-\frac{2}{\Delta
M_{d}}\text{Im}\left[M_{12}^{\ast}({\rm NP})+M_{12}^{\ast}({\rm
SM})\right]_{B_{d}}\ ,$ (123) $S_{\psi\phi}=-\zeta_{\psi\phi}\sin
2(\beta_{s})_{\text{exp}}\approx-\frac{2}{\Delta
M_{s}}\text{Im}\left[M_{12}^{\ast}({\rm NP})+M_{12}^{\ast}({\rm
SM})\right]_{B_{s}}\ ,$ (124)

and the SM terms in equations (123) and (124) represent respectively $\sin
2\beta_{\text{SM}}$ and $\sin 2(\beta_{s})_{\text{SM}}$. Thus from equations
(123) and (108), we can easily find

$\displaystyle\sin[2\beta_{\text{exp}}]-\sin[2\beta_{\text{SM}}]=$
$\displaystyle-\frac{2}{\Delta M_{d}}\text{Im}\left[M_{12}^{\ast}({\rm
NP})\right]$ (125) $\displaystyle=$ $\displaystyle-\frac{1}{m_{B_{d}}\Delta
M_{d}}\text{Im}\left[\bra{B_{d}^{0}}H_{\rm eff}^{\rm NP}(\Delta
B=2)\ket{B_{d}^{0}}^{\ast}\right]$ $\displaystyle=$
$\displaystyle-\frac{m_{B_{d}}F_{B_{d}}^{2}\hat{B}_{B_{d}}\eta_{B_{d}}}{192{\pi}^{2}\Delta
M_{d}M_{X}^{2}}N_{c}\text{Im}\left[(\xi_{B_{d}}^{\ast})^{2}\right]\ ,$

where we take the experimental values for $m_{B_{d}}$ and $\Delta M_{d}$.
Using the numerical values of Table 4, we find the constraint

$-0.86\times 10^{-4}\leq
N_{c}\text{Im}\left[({\xi}_{B_{d}}^{*})^{2}\right]\left(\frac{1\text{
TeV}}{M_{X}}\right)^{2}\leq 3.12\times 10^{-4}.$ (126)

We apply the same procedure for the $B_{s}-\bar{B}_{s}$ mixing, and we find
the following constraint from the difference between the experimental and the
SM values of the $CP-$violating observable $\beta_{s}$:

$-3.54\times 10^{-3}\leq
N_{c}\text{Im}\left[({\xi}_{B_{s}}^{*})^{2}\right]\left(\frac{1\text{
TeV}}{M_{X}}\right)^{2}\leq 3.88\times 10^{-3}.$ (127)

All the constraints are summarized in Table 1 in Section 4. From these results
we can infer model independent bounds on any new-physics model, contributing
to flavor observables, and which $\Delta F=2$ effective Hamiltonian, at some
scale $\Lambda$, contains (V+A)(V+A) (or (V-A)(V-A), since strong interactions
conserve $CP$) operators in the form:

$H_{\rm eff}^{\rm NP}(\Delta
F=2)=\frac{c_{qq^{\prime}}^{NP}}{\Lambda^{2}}\left[(\bar{q}\gamma_{\mu}P_{R(L)}q^{\prime})(\bar{q}\gamma^{\mu}P_{R(L)}q^{\prime})\right]+{\rm
h.c.}.$ (128)

The scale $M_{X}$ can be mapped to the scale suppressing the dimension 6
operator, for example $c_{ds}^{\text{NP}}=N_{c}{\xi}_{K}^{2}/(128{\pi}^{2})$.

Before closing this section, we compare our results to the work in
IsidoriEtAL010  which performed a similar analysis. At a scale $\Lambda\approx
M_{X}\approx 1\text{ TeV}$, equation (113) yields:

$\text{Re}[c_{ds}^{\rm NP}]\left(\Lambda=1\text{TeV}\right)\leq 2.98\times
10^{-7},$ (129)

which is at the same order of magnitude, but three times stronger than that of
reference IsidoriEtAL010 . In addition, the $CP-$violating constraint we find
in equation (122) is also about three times stronger than the constraint of
reference IsidoriEtAL010 . We believe these disagreements originate from the
fact that this reference took the SM leading contribution to $K-\bar{K}$
mixing amplitude with only the top quark running inside the loop Isidori015 ;
NirFPandCPV . This is only true for the $B-$ system, while for the $K-$system,
the dominant contribution to the real part of the amplitude is from an
internal $c$ quark, and the dominant contribution to the imaginary part is
from both $c$ and $t$ internal quarks (See section 10.2 and 10.4 of reference
Buras1998 ).

The constraints we find in equations (116), (118), and (126), are all about
one order of magnitude stronger than the ones of reference IsidoriEtAL010 . We
believe this is due to the improvement in the experimental data, which have
significantly changed for the $B-$system. We find a result that is the same
order of magnitude as the one in IsidoriEtAL010  when using the same input
data.

## References

  * (1) S. Nussinov, _TECHNOCOSMOLOGY: COULD A TECHNIBARYON EXCESS PROVIDE A ’NATURAL’ MISSING MASS CANDIDATE?_ , _Phys. Lett. B_ 165 (1985) 55.
  * (2) S. M. Barr, R. S. Chivukula and E. Farhi, _Electroweak Fermion Number Violation and the Production of Stable Particles in the Early Universe_ , _Phys. Lett. B_ 241 (1990) 387.
  * (3) S. M. Barr, _Baryogenesis, sphalerons and the cogeneration of dark matter_ , _Phys. Rev. D_ 44 (1991) 3062.
  * (4) S. Dodelson, B. R. Greene and L. M. Widrow, _Baryogenesis, dark matter and the width of the Z_ , _Nucl. Phys. B_ 372 (1992) 467.
  * (5) D. B. Kaplan, _A Single explanation for both the baryon and dark matter densities_ , _Phys. Rev. Lett._ 68 (1992) 741.
  * (6) V. A. Kuzmin, _A Simultaneous solution to baryogenesis and dark matter problems_ , _Phys. Part. Nucl._ 29 (1998) 257 [hep-ph/9701269].
  * (7) R. Foot and R. R. Volkas, _Was ordinary matter synthesized from mirror matter? An Attempt to explain why Omega(Baryon) approximately equal to 0.2 Omega(Dark)_ , _Phys. Rev. D_ 68 (2003) 021304 [hep-ph/0304261].
  * (8) R. Foot and R. R. Volkas, _Explaining Omega(Baryon) approximately 0.2 Omega(Dark) through the synthesis of ordinary matter from mirror matter: A More general analysis_ , _Phys. Rev. D_ 69 (2004) 123510 [hep-ph/0402267].
  * (9) D. Hooper, J. March-Russell and S. M. West, _Asymmetric sneutrino dark matter and the Omega(b) / Omega(DM) puzzle_ , _Phys. Lett. B_ 605 (2005) 228 [hep-ph/0410114].
  * (10) R. Kitano and I. Low, _Dark matter from baryon asymmetry_ , _Phys. Rev. D_ 71 (2005) 023510 [hep-ph/0411133].
  * (11) S. B. Gudnason, C. Kouvaris and F. Sannino, _Towards working technicolor: Effective theories and dark matter_ , _Phys. Rev. D_ 73 (2006) 115003 [hep-ph/0603014].
  * (12) D. E. Kaplan, M. A. Luty and K. M. Zurek, _Asymmetric Dark Matter_ , _Phys. Rev. D_ 79 (2009) 115016 [0901.4117].
  * (13) H. Davoudiasl and R. N. Mohapatra, _On Relating the Genesis of Cosmic Baryons and Dark Matter_ , _New J. Phys._ 14 (2012) 095011 [1203.1247].
  * (14) K. Petraki and R. R. Volkas, _Review of asymmetric dark matter_ , _Int. J. Mod. Phys. A_ 28 (2013) 1330028 [1305.4939].
  * (15) K. M. Zurek, _Asymmetric Dark Matter: Theories, Signatures, and Constraints_ , _Phys. Rept._ 537 (2014) 91 [1308.0338].
  * (16) Y. Bai and P. Schwaller, _Scale of dark QCD_ , _Phys. Rev._ D89 (2014) 063522 [1306.4676].
  * (17) N. Craig, A. Katz, M. Strassler and R. Sundrum, _Naturalness in the Dark at the LHC_ , _JHEP_ 07 (2015) 105 [1501.05310].
  * (18) D. Curtin and C. B. Verhaaren, _Discovering Uncolored Naturalness in Exotic Higgs Decays_ , _JHEP_ 12 (2015) 072 [1506.06141].
  * (19) H.-C. Cheng, E. Salvioni and Y. Tsai, _Exotic electroweak signals in the twin Higgs model_ , _Phys. Rev. D_ 95 (2017) 115035 [1612.03176].
  * (20) C. Kilic, S. Najjari and C. B. Verhaaren, _Discovering the Twin Higgs Boson with Displaced Decays_ , _Phys. Rev. D_ 99 (2019) 075029 [1812.08173].
  * (21) M. J. Strassler and K. M. Zurek, _Echoes of a hidden valley at hadron colliders_ , _Phys. Lett. B_ 651 (2007) 374 [hep-ph/0604261].
  * (22) M. J. Strassler and K. M. Zurek, _Discovering the Higgs through highly-displaced vertices_ , _Phys. Lett. B_ 661 (2008) 263 [hep-ph/0605193].
  * (23) T. Han, Z. Si, K. M. Zurek and M. J. Strassler, _Phenomenology of hidden valleys at hadron colliders_ , _JHEP_ 07 (2008) 008 [0712.2041].
  * (24) P. Schwaller, D. Stolarski and A. Weiler, _Emerging Jets_ , _JHEP_ 05 (2015) 059 [1502.05409].
  * (25) T. Cohen, M. Lisanti and H. K. Lou, _Semivisible Jets: Dark Matter Undercover at the LHC_ , _Phys. Rev. Lett._ 115 (2015) 171804 [1503.00009].
  * (26) C. Csaki, E. Kuflik, S. Lombardo and O. Slone, _Searching for displaced Higgs boson decays_ , _Phys. Rev. D_ 92 (2015) 073008 [1508.01522].
  * (27) S. Knapen, S. Pagan Griso, M. Papucci and D. J. Robinson, _Triggering Soft Bombs at the LHC_ , _JHEP_ 08 (2017) 076 [1612.00850].
  * (28) A. Pierce, B. Shakya, Y. Tsai and Y. Zhao, _Searching for confining hidden valleys at LHCb, ATLAS, and CMS_ , _Phys. Rev. D_ 97 (2018) 095033 [1708.05389].
  * (29) G. Burdman and G. Lichtenstein, _Displaced Vertices from Hidden Glue_ , _JHEP_ 08 (2018) 146 [1807.03801].
  * (30) H.-C. Cheng, L. Li, E. Salvioni and C. B. Verhaaren, _Light Hidden Mesons through the Z Portal_ , _JHEP_ 11 (2019) 031 [1906.02198].
  * (31) D. Linthorne and D. Stolarski, _Triggering on Emerging Jets_ , 2103.08620.
  * (32) S. Knapen, J. Shelton and D. Xu, _Perturbative benchmark models for a dark shower search program_ , _Phys. Rev. D_ 103 (2021) 115013 [2103.01238].
  * (33) J. Alimena et al., _Searching for long-lived particles beyond the Standard Model at the Large Hadron Collider_ , _J. Phys. G_ 47 (2020) 090501 [1903.04497].
  * (34) CMS collaboration, _Search for new particles decaying to a jet and an emerging jet_ , _JHEP_ 02 (2019) 179 [1810.10069].
  * (35) H. Mies, C. Scherb and P. Schwaller, _Collider constraints on dark mediators_ , 2011.13990.
  * (36) J. Kile and A. Soni, _Flavored Dark Matter in Direct Detection Experiments and at LHC_ , _Phys. Rev. D_ 84 (2011) 035016 [1104.5239].
  * (37) J. F. Kamenik and J. Zupan, _Discovering Dark Matter Through Flavor Violation at the LHC_ , _Phys. Rev. D_ 84 (2011) 111502 [1107.0623].
  * (38) B. Batell, J. Pradler and M. Spannowsky, _Dark Matter from Minimal Flavor Violation_ , _JHEP_ 08 (2011) 038 [1105.1781].
  * (39) P. Agrawal, S. Blanchet, Z. Chacko and C. Kilic, _Flavored Dark Matter, and Its Implications for Direct Detection and Colliders_ , _Phys. Rev. D_ 86 (2012) 055002 [1109.3516].
  * (40) B. Batell, T. Lin and L.-T. Wang, _Flavored Dark Matter and R-Parity Violation_ , _JHEP_ 01 (2014) 075 [1309.4462].
  * (41) L. Lopez-Honorez and L. Merlo, _Dark matter within the minimal flavour violation ansatz_ , _Phys. Lett. B_ 722 (2013) 135 [1303.1087].
  * (42) P. Agrawal, M. Blanke and K. Gemmler, _Flavored dark matter beyond Minimal Flavor Violation_ , _JHEP_ 10 (2014) 072 [1405.6709].
  * (43) F. Bishara and J. Zupan, _Continuous Flavor Symmetries and the Stability of Asymmetric Dark Matter_ , _JHEP_ 01 (2015) 089 [1408.3852].
  * (44) A. Hamze, C. Kilic, J. Koeller, C. Trendafilova and J.-H. Yu, _Lepton-Flavored Asymmetric Dark Matter and Interference in Direct Detection_ , _Phys. Rev. D_ 91 (2015) 035009 [1410.3030].
  * (45) L. Calibbi, A. Crivellin and B. Zaldívar, _Flavor portal to dark matter_ , _Phys. Rev. D_ 92 (2015) 016004 [1501.07268].
  * (46) P. Agrawal, Z. Chacko, C. Kilic and C. B. Verhaaren, _A Couplet from Flavored Dark Matter_ , _JHEP_ 08 (2015) 072 [1503.03057].
  * (47) F. Bishara, A. Greljo, J. F. Kamenik, E. Stamou and J. Zupan, _Dark Matter and Gauged Flavor Symmetries_ , _JHEP_ 12 (2015) 130 [1505.03862].
  * (48) B. Bhattacharya, D. London, J. M. Cline, A. Datta and G. Dupuis, _Quark-flavored scalar dark matter_ , _Phys. Rev. D_ 92 (2015) 115012 [1509.04271].
  * (49) M.-C. Chen, J. Huang and V. Takhistov, _Beyond Minimal Lepton Flavored Dark Matter_ , _JHEP_ 02 (2016) 060 [1510.04694].
  * (50) P. Agrawal, Z. Chacko, E. C. F. S. Fortes and C. Kilic, _Skew-Flavored Dark Matter_ , _Phys. Rev. D_ 93 (2016) 103510 [1511.06293].
  * (51) M. Blanke, S. Das and S. Kast, _Flavoured Dark Matter Moving Left_ , _JHEP_ 02 (2018) 105 [1711.10493].
  * (52) S. Renner and P. Schwaller, _A flavoured dark sector_ , _JHEP_ 08 (2018) 052 [1803.08080].
  * (53) M. Blanke and S. Kast, _Top-Flavoured Dark Matter in Dark Minimal Flavour Violation_ , _JHEP_ 05 (2017) 162 [1702.08457].
  * (54) T. Jubb, M. Kirk and A. Lenz, _Charming Dark Matter_ , _JHEP_ 12 (2017) 010 [1709.01930].
  * (55) M. Blanke, S. Das and S. Kast, _Flavoured Dark Matter Moving Left_ , _JHEP_ 02 (2018) 105 [1711.10493].
  * (56) D. E. Kaplan and M. D. Schwartz, _Constraining Light Colored Particles with Event Shapes_ , _Phys. Rev. Lett._ 101 (2008) 022002 [0804.2477].
  * (57) D. S. M. Alves, J. Galloway, J. T. Ruderman and J. R. Walsh, _Running Electroweak Couplings as a Probe of New Physics_ , _JHEP_ 02 (2015) 007 [1410.6810].
  * (58) OPAL collaboration, _Search for scalar top and scalar bottom quarks at LEP_ , _Phys. Lett. B_ 545 (2002) 272 [hep-ex/0209026].
  * (59) J. M. Pendlebury et al., _Revised experimental upper limit on the electric dipole moment of the neutron_ , _Phys. Rev. D_ 92 (2015) 092003 [1509.04411].
  * (60) T. Chupp, P. Fierlinger, M. Ramsey-Musolf and J. Singh, _Electric dipole moments of atoms, molecules, nuclei, and particles_ , _Rev. Mod. Phys._ 91 (2019) 015001 [1710.02504].
  * (61) T. Liu, Z. Zhao and H. Gao, _Experimental constraint on quark electric dipole moments_ , _Phys. Rev._ D97 (2018) 074018 [1704.00113].
  * (62) ACME collaboration, _Improved limit on the electric dipole moment of the electron_ , _Nature_ 562 (2018) 355.
  * (63) M. Le Dall and A. Ritz, _CP-violation and electric dipole moments_ , _Hyperfine Interact._ 214 (2013) 87.
  * (64) M. Pospelov and A. Ritz, _Electric dipole moments as probes of new physics_ , _Annals Phys._ 318 (2005) 119 [hep-ph/0504231].
  * (65) A. Hook, _TASI Lectures on the Strong CP Problem and Axions_ , _PoS_ TASI2018 (2019) 004 [1812.02669].
  * (66) S. Weinberg, _Larger Higgs Exchange Terms in the Neutron Electric Dipole Moment_ , _Phys. Rev. Lett._ 63 (1989) 2333.
  * (67) C. Hamzaoui and M. E. Pospelov, _Electric dipole moment of neutron in the Kobayashi-Maskawa model with four generations of quarks_ , _Phys. Lett._ B357 (1995) 616 [hep-ph/9503468].
  * (68) M. E. Peskin and D. V. Schroeder, _An Introduction to quantum field theory_. Addison-Wesley, Reading, USA, 1995.
  * (69) S. M. Barr and A. Zee, _Electric Dipole Moment of the Electron and of the Neutron_ , _Phys. Rev. Lett._ 65 (1990) 21.
  * (70) M. Fujiwara, J. Hisano and T. Toma, _Vanishing or non-vanishing rainbow? Reduction formulas of electric dipole moment_ , 2106.03384.
  * (71) M. Fujiwara, J. Hisano, C. Kanai and T. Toma, _Electric dipole moments in the extended scotogenic models_ , _JHEP_ 04 (2021) 114 [2012.14585].
  * (72) F. Farchioni, I. Montvay, G. Munster, E. E. Scholz, T. Sudmann and J. Wuilloud, _Hadron masses in QCD with one quark flavour_ , _Eur. Phys. J. C_ 52 (2007) 305 [0706.1131].
  * (73) Belle collaboration, _Search for $B^{0}$ decays to invisible final states_, _Phys. Rev._ D86 (2012) 032002 [1206.5948].
  * (74) BaBar collaboration, _Improved Limits on $B^{0}$ Decays to Invisible Final States and to $\nu\bar{\nu}\gamma$_, _Phys. Rev._ D86 (2012) 051105 [1206.2543].
  * (75) Belle-II collaboration, _The Belle II Physics Book_ , _PTEP_ 2019 (2019) 123C01 [1808.10567].
  * (76) A. Badin and A. A. Petrov, _Searching for light Dark Matter in heavy meson decays_ , _Phys. Rev._ D82 (2010) 034005 [1005.1277].
  * (77) B. Bhattacharya, C. M. Grant and A. A. Petrov, _Invisible widths of heavy mesons_ , _Phys. Rev._ D99 (2019) 093010 [1809.04606].
  * (78) A. Ali and C. Greub, _An Analysis of two-body nonleptonic B decays involving light mesons in the standard model_ , _Phys. Rev._ D57 (1998) 2996 [hep-ph/9707251].
  * (79) Particle Data Group collaboration, _Review of Particle Physics_ , _Chin. Phys._ C40 (2016) 100001.
  * (80) A. Ali, P. Ball, L. T. Handoko and G. Hiller, _A Comparative study of the decays $B\to$ ($K$, $K^{*)}\ell^{+}\ell^{-}$ in standard model and supersymmetric theories_, _Phys. Rev._ D61 (2000) 074024 [hep-ph/9910221].
  * (81) D. Melikhov and B. Stech, _Weak form-factors for heavy meson decays: An Update_ , _Phys. Rev._ D62 (2000) 014006 [hep-ph/0001113].
  * (82) J. M. Flynn, T. Izubuchi, T. Kawanai, C. Lehner, A. Soni, R. S. Van de Water et al., _$B\to\pi\ell\nu$ and $B_{s}\to K\ell\nu$ form factors and $|V_{ub}|$ from 2+1-flavor lattice QCD with domain-wall light quarks and relativistic heavy quarks_, _Phys. Rev._ D91 (2015) 074510 [1501.05373].
  * (83) Fermilab Lattice, MILC collaboration, _$B\to\pi\ell\nu$ semileptonic form factors from unquenched lattice QCD and determination of $|V_{ub}|$_, _PoS_ LATTICE2014 (2014) 385 [1411.6038].
  * (84) J. A. Bailey et al., _$B\to Kl^{+}l^{-}$ Decay Form Factors from Three-Flavor Lattice QCD_, _Phys. Rev. D_ 93 (2016) 025026 [1509.06235].
  * (85) Belle collaboration, _Search for $\boldsymbol{B\to h\nu\bar{\nu}}$ decays with semileptonic tagging at Belle_, _Phys. Rev._ D96 (2017) 091101 [1702.03224].
  * (86) C. Hambrock, A. Khodjamirian and A. Rusov, _Hadronic effects and observables in $B\to\pi\ell^{+}\ell^{-}$ decay at large recoil_, _Phys. Rev._ D92 (2015) 074020 [1506.07760].
  * (87) A. J. Buras, J. Girrbach-Noe, C. Niehoff and D. M. Straub, _$B\to{K}^{\left(\ast\right)}\nu\overline{\nu}$ decays in the Standard Model and beyond_, _JHEP_ 02 (2015) 184 [1409.4557].
  * (88) A. J. Buras, _Flavor physics and CP violation_ , in _2004 European School of High-Energy Physics, Sant Feliu de Guixols, Spain, 30 May - 12 June 2004_ , pp. 95–168, 2005, hep-ph/0505175, http://doc.cern.ch/yellowrep/2006/2006-003/p95.pdf.
  * (89) A. J. Buras, _Flavor dynamics: CP violation and rare decays_ , _Subnucl. Ser._ 38 (2002) 200 [hep-ph/0101336].
  * (90) Belle II collaboration, _Search for $B^{+}\to K^{+}\nu\bar{\nu}$ decays using an inclusive tagging method at Belle II_, 2104.12624.
  * (91) A. J. Buras, D. Buttazzo, J. Girrbach-Noe and R. Knegjens, _${K}^{+}\to{\pi}^{+}\nu\overline{\nu}$ and ${K}_{L}\to{\pi}^{0}\nu\overline{\nu}$ in the Standard Model: status and perspectives_, _JHEP_ 11 (2015) 033 [1503.02693].
  * (92) KOTO collaboration, _Study of the $K_{L}\\!\to\\!\pi^{0}\nu\overline{\nu}$ Decay at the J-PARC KOTO Experiment_, _Phys. Rev. Lett._ 126 (2021) 121801 [2012.07571].
  * (93) NA62 collaboration, _An investigation of the very rare ${K}^{+}\to{\pi}^{+}\nu\overline{\nu}$ decay_, _JHEP_ 11 (2020) 042 [2007.08218].
  * (94) N. Carrasco, P. Lami, V. Lubicz, L. Riggio, S. Simula and C. Tarantino, _$K\to\pi$ semileptonic form factors with $N_{f}=2+1+1$ twisted mass fermions_, _Phys. Rev. D_ 93 (2016) 114512 [1602.04113].
  * (95) Y. Grossman and Y. Nir, _K(L) — > pi0 neutrino anti-neutrino beyond the standard model_, _Phys. Lett._ B398 (1997) 163 [hep-ph/9701313].
  * (96) H. Nanjo, “Rare kaon decay experiments.” 2020.
  * (97) L.-L. Chau, W.-Y. Keung and M. D. Tran, _An Updated Analysis of the Quark Mixing Matrix and Its Phenomenological Implications_ , _Phys. Rev._ D27 (1983) 2145.
  * (98) A. J. Buras, D. Guadagnoli and G. Isidori, _On $\epsilon_{K}$ Beyond Lowest Order in the Operator Product Expansion_, _Phys. Lett._ B688 (2010) 309 [1002.3612].
  * (99) Y. Nir, _CP violation in and beyond the standard model_ , in _Proceedings, 27th SLAC Summer Institute on Particle Physics: CP Violation in and Beyond the Standard Model (SSI 99): Stanford, USA, July 7-16, 1999_ , pp. 165–243, 1999, hep-ph/9911321, http://www-public.slac.stanford.edu/sciDoc/docMeta.aspx?slacPubNumber=slac-r-719.
  * (100) A. J. Buras and J. Girrbach, _Towards the Identification of New Physics through Quark Flavour Violating Processes_ , _Rept. Prog. Phys._ 77 (2014) 086201 [1306.3775].
  * (101) A. J. Buras, W. Slominski and H. Steger, _B0 anti-B0 Mixing, CP Violation and the B Meson Decay_ , _Nucl. Phys._ B245 (1984) 369.
  * (102) A. J. Buras, _Weak Hamiltonian, CP violation and rare decays_ , in _Probing the standard model of particle interactions. Proceedings, Summer School in Theoretical Physics, NATO Advanced Study Institute, 68th session, Les Houches, France, July 28-September 5, 1997. Pt. 1, 2_ , pp. 281–539, 1998, hep-ph/9806471.
  * (103) UTfit collaboration, _Unitarity Triangle analysis in the Standard Model from the UTfit collaboration_ , _PoS_ ICHEP2016 (2016) 554.
  * (104) Y. Nir, _Flavour Physics and CP Violation_ , in _Proceedings, 7th CERN–Latin-American School of High-Energy Physics (CLASHEP2013): Arequipa, Peru, March 6-19, 2013_ , pp. 123–156, 2015, 1605.00433, DOI.
  * (105) A. J. Bevan, G. Inguglia and B. Meadows, _Time-dependent $CP$ asymmetries in $D$ and $B$ decays_, _Phys. Rev._ D84 (2011) 114009 [1106.5075].
  * (106) C. Linn, _Measurement of the CP-violating phase $\phi_{s}$ using $B_{s}\rightarrow J/\psi\phi$ and $B_{s}^{0}\rightarrow J/\psi\pi^{+}\pi^{-}$ decays with the LHCb experiment_, Ph.D. thesis, Heidelberg U., 4, 2013.
  * (107) J. Laiho, E. Lunghi and R. S. Van de Water, _Lattice QCD inputs to the CKM unitarity triangle analysis_ , _Phys. Rev._ D81 (2010) 034503 [0910.2928].
  * (108) S. Aoki et al., _Review of lattice results concerning low-energy particle physics_ , _Eur. Phys. J._ C77 (2017) 112 [1607.00299].
  * (109) UTfit collaboration, _Model-independent constraints on $\Delta F=2$ operators and the scale of new physics_, _JHEP_ 03 (2008) 049 [0707.0636].
  * (110) UTfit collaboration, _Unitarity Triangle analysis beyond the Standard Model from UTfit_ , _PoS_ ICHEP2016 (2016) 149.
  * (111) Heavy Flavor Averaging Group (HFAG) collaboration, _Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of summer 2014_, 1412.7515.
  * (112) G. Isidori, Y. Nir and G. Perez, _Flavor Physics Constraints for Physics Beyond the Standard Model_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 355 [1002.0900].
  * (113) G. Isidori, _Flavour Physics and Implication for New Phenomena_ , _Adv. Ser. Direct. High Energy Phys._ 26 (2016) 339 [1507.00867].