# Vocal Fold Reconstruction from Optical Velocity and Displacement
Measurements

Daniel Zieger
Department of Computer Science
Friedrich-Alexander-Universität Erlangen-Nürnberg
Erlangen, Germany
<EMAIL_ADDRESS>
&Christoph Näger
Institute of Fluid Mechanics
Friedrich-Alexander-Universität Erlangen-Nürnberg
Erlangen, Germany
<EMAIL_ADDRESS>
&Stefan Becker
Institute of Fluid Mechanics
Friedrich-Alexander-Universität Erlangen-Nürnberg
Erlangen, Germany
<EMAIL_ADDRESS>
&Tobias Günther
Department of Computer Science
Friedrich-Alexander-Universität Erlangen-Nürnberg
Erlangen, Germany
<EMAIL_ADDRESS>

###### Abstract

The three-dimensional reconstruction of vocal folds in medicine usually
involves endoscopy and an approach to extract depth information like
structured light or stereo matching of images. The resulting mesh can
accurately represent the superior area of the vocal folds, while new
approaches also try to reconstruct the inferior area. We propose a novel
approach to extract the time-dependent 3D geometry of the vocal fold from
optical measurements on both the superior and inferior side, requiring optical
measurements only from the superior side. First, a time-dependent, tri-variate
surface velocity vector field is reconstructed using a high-speed camera and a
laser Doppler vibrometer in an experimental environment. This vector field
serves as target in an inverse finite-element simulation that optimizes the
forces applied to a deformable vocal fold model such that the resulting
movement after FEM simulation matches the velocity observations on the
superior side. The required forces for the finite element method simulation
are treated as unknowns and are assembled using multiple scalar fields. We use
tensor products in Bézier Bernstein basis for our scalar fields to reduce the
degrees of freedom for our optimization. We use gradient descent to optimize
the control points of the force field polynomials. Our utilized error metric
for gradient descent consists of two terms. The first term is used to match
the simulated velocities to the observed measurements, while the second term
measures the silhouette difference between observation and simulation.

##  1 Introduction

Human communication is primarily based on the voice. Consequently, the voice
has a great impact on the daily life of a person, as the quality of life of an
individual is significantly decreased if they suffer from voice impairments
[1]. Physically speaking, the human voice is generated in a complex interplay
of fluid flow, structural vibration, and acoustics. This process is called
phonation. During phonation, the vocal folds are excited by a flow of air
$\dot{V}$ from the lungs, causing them to vibrate. The vocal fold oscillation
in return entails a modulation of this airflow, resulting in a pulsating jet
flow in the vocal tract. This modulated jet flow is the main acoustic source
in phonation, as discussed by Kaltenbacher et al. [2]. The generated sound is
filtered through the vocal tract and radiated from the mouth, resulting in the
voice. Even though the vocal fold vibration is not the main direct sound
source in phonation, it still plays a very important part by its interaction
with the air flow. As it is difficult to perform extensive measurements of the
glottal air flow in vivo due to a restricted accessibility, clinicians often
rely on measurements of the vocal fold vibration for the diagnosis of speech
impairments. Endoscopic two-dimensional imaging is thereby the most
established method. This has been extended to high-speed videoendoscopy with
frame rates up to 20 kHz in the past to be able to resolve the full
oscillation cycle of the vocal folds in time [3, 4]. However, the work by
Döllinger and Berry [5] showed that there is a non-negligible out-of-plane
movement present in the vocal fold vibration, which is not captured with
conventional two-dimensional endoscopy. To be able to obtain a three-
dimensional vocal fold movement, a stereo triangulation technique has been
developed. For this, often a glass prism is used to generate two different
viewpoints with one endoscope, cf. [5]. Also two independent optical systems
that are displaced by a known distance have been used [6]. To facilitate
feature matching, marker points can be generated by projecting a structured
light pattern onto the vocal fold surface [7, 8, 9, 10]. To help practitioners
with diagnosing laryngeal disorders, Semmler et al. [10] developed a method of
reconstructing the superior vocal fold surface as a 3D-model. This has been
further developed by Henningson et al. [11], who reconstructed the inferior
vocal fold surface via B-spline interpolation, resulting in a full 3D vocal
fold model based on the M5-model by Scherer et al. [12]. The inferior vocal
fold surface movement here is however not necessarily physically valid, as no
physical constraints are used when applying the B-spline model. Nevertheless,
this is still a useful representation for helping practitioners making a
diagnosis. For a more detailed physical analysis, however, it would certainly
be more advantageous to apply a more physics-based approach of vocal fold
reconstruction, which is the contribution of this paper.

One application, in which 3D vocal fold models are currently used is the
simulation of the phonation process via computational fluid dynamics (CFD). In
fully coupled fluid-structure interaction simulations, the vocal fold movement
is directly simulated by coupling CFD with a lumped element model or finite
element simulation. As an example for the use of lumped elements, Erath et al.
[13] used a two mass model of the vocal folds in combination with boundary
layer methods for the flow. Using lumped elements is a strong simplification
of the vocal fold movement. One example for coupling a finite element solver
to a Navier-Stokes based flow solver is the work by Zheng et al. [14]. They
applied the immersed boundary method to analyse the aerodynamics of the
phonation process. This approach can generate very accurate results, while at
the same time a significant amount of computational resources is necessary to
perform such kind of simulations. For this reason, one-way coupled fluid-
structure simulations have emerged as a useful method for simulating human
phonation. Here, the vocal fold movement is being prescribed and therefore not
influenced by the surface pressure arising from the intraglottal flow field.
With this, the computationally costly fluid-structure interaction simulation
is reduced to a pure CFD-simulation. This method was applied first by Zörner
et al. [15] and was developed further in several publications with the goal of
making forward-coupled CFD-simulations clinically applicable for voice
treatment [16, 17, 18, 19]. Furthermore, an extensive analysis of the
aeroacoustic sources in phonation could be performed by Schoder et al. [20],
applying the perturbed convective wave equation. These simulations used a two-
parameter model for the vocal fold movement, such that the vocal fold motion
has a translational part in medial-lateral direction and a rotational part
generating the typical convergent-divergent glottal gap. The velocities of
opening and closing the gap were chosen to reproduce the glottal area waveform
of an experimental model. However, as stated before, a physiologically
realistic vocal fold movement features a significant out-of-plane movement
[5], which is not included in these models. Therefore, one-way coupled fluid-
structure simulations could benefit from an improved prescribed vocal fold
movement model.

Considering these two main applications, the present work introduces a
procedure to reconstruct a physically valid three-dimensional vocal fold model
from three-dimensional velocity measurements. The velocity data is obtained
experimentally from a synthetic vocal fold model via high-speed-camera imaging
and laser-Doppler-vibrometry. Nonlinear optimization is applied to fit the
movement of a finite-element-model of the vocal folds to the measured velocity
data in an optimal sense. Additionally, the glottal area shape is extracted
from the high-speed-camera images and included in the optimization procedure.
The result is a three-dimensional vocal fold model that can be used as an
input for future CFD-simulations of the phonation process with prescribed
vocal fold movement. Since a finite-element-model serves as shape prior, the
optimization of a full three-dimensional model (inferior and superior side)
can be performed from optical measurements of the superior side alone. Thus,
the same procedure can also be used in the medical diagnosis of voice
impairments via videoendoscopy by applying it to the patient specific vocal
fold movement. In general, the proposed procedure can be applied to reproduce
any three-dimensional displacement field of a solid body from three-
dimensional velocity data while at the same time satisfying the underlying
physical equations.

(a)  Schematic illustration of the experimental setup.

(b)  Shape parameters of the vocal fold model.

Figure 1: Illustration of the experimental setup. In (a), the images show 1)
the labeled parts of the entire setup, 2) the components of the high-speed
imaging, and 3) the components for the Laser-Doppler-Vibrometry. The vocal
folds are highlighted in red. In (b), shape parameters of the used vocal folds
model are shown with the units given in millimeters.

##  2 Experimental Setup

The experimental setup is a slightly modified version of the setup presented
by Lodermeyer et al. [21] and is shown in Fig. 1. A mass flow generator
containing a supercritical valve (see [22]) is responsible for the flow
generation through the test rig. The total mass flow enters a silencer that
dampens acoustic fluctuations coming from the supply hose. After the silencer,
a special curved subglottal channel is added to allow for optical access to
the inferior side of the vocal folds. After the vocal folds, a uniform vocal
tract with a length of 200 mm is attached. The length is chosen such that the
acoustic resonance frequencies of the vocal tract do not evoke acoustically
driven vocal folds as proposed by Zhang et al. [23]. The material of the vocal
folds is an elastomer with a Young’s modulus of $E=4.4$ kPa ([21]) and a
Poisson’s ratio of $\nu=0.499$ ([24]). Two different measurement methods were
chosen to obtain the three-dimensional velocities at the superior as well as
at the inferior vocal fold surfaces. They were conducted sequentially and
synchronized in a post processing step by making use of the simultaneously
acquired subglottal pressure signal. The subglottal pressure signal was
recorded using a Kulite XCQ-093 pressure sensor (Kulite Semiconductor Products
Inc., Leonia, NJ, USA) that was located 50mm upstream of the glottal exit at a
measurement frequency of 50 kHz. Planar high-speed imaging of the superior
vocal fold surface was performed to obtain information of the vocal fold
movement in the transversal plane. The images were acquired by a Photron SA-X2
high speed camera (Photron, Tokyo, Japan) in combination with a AF-S DX NIKKOR
18-300 mm 1:3.5-6.3G ED VR (Nikon, Tokyo, Japan) at a rate of 10,000 frames
per second. For the out-of-plane movement, Laser-Doppler-Vibrometry was
applied, using a Polytec PSV-500 Scanning Vibrometer (Polytec, Waldbronn,
Germany). A grid of $17\times 22$ measurement points on the superior vocal
fold surface was chosen, resulting in a total of 374 measurement points. The
measurement frequency of the vibrometer was set to 50 kHz.

Figure 2: Our experimental setup delivers high-speed camera images
$I(\mathbf{y}_{i},t_{i})$ and discrete vertical velocity measurements
$w_{i}(\mathbf{y}_{i},t_{i})$. In a precomputation, we reconstruct the
horizontal components $(u_{I}(\mathbf{y},t),v_{I}(\mathbf{y},t))$ and the
vertical component $w(\mathbf{y},t)$ of the surface velocity field
$\mathbf{v}$, as well as an unsigned distance field $\sigma(\mathbf{y},t)$
that measures the distance to the vocal fold opening. We use a finite element
(FEM) simulation that takes an initial vocal fold model $\mathbf{x}_{0}(u,v)$
and a surface force field $\mathbf{f}(u,v,t)$ as input and calculates a
deformed model with positions $\overline{\mathbf{x}}(u,v,t)$ and velocities
$\overline{\mathbf{v}}(u,v,t)$. Using gradient descent, we optimize for the
force field $\mathbf{f}(u,v,t)$ such that the simulated vocal fold motion and
the opening match the observations, i.e., the measurements that were processed
in the precomputation.

##  3 Vocal Fold Reconstruction

###  3.1 Overview

Our computational vocal fold reconstruction pipeline is divided into two
steps: (1) precomputation and (2) optimization. A schematic illustration is
provided in Fig. 2. In the precomputation, the discrete vibrometer
measurements are interpolated in space and time, resulting in the vertical
component of the vocal fold target surface velocity field. Further, the
horizontal component of the target surface velocity field is estimated using
an optical flow algorithm from the high-speed camera footage. In addition, a
time-dependent unsigned distance field is computed that measures the distance
to the silhouette of the opening, i.e., the gap. The optimization takes a
deformable vocal fold model, the reconstructed fields, and a force vector
field as input. The latter is applied in a finite-element method simulation to
the initial vocal fold model, which results in a deformed vocal fold model.
From the deformed model, the resulting surface vector field and the resulting
opening are compared to the optical measurements that were precomputed. Via
gradient descent, the force field is gradually improved to reduce the distance
between the deformed vocal fold model and the measurement data. In the
following sections, the individual steps are described formally in more
detail.

###  3.2 Reconstruction of Surface Vector Field from Optical Measurements

In this section, we reconstruct a continuous, time-dependent velocity vector
field
$\mathbf{v}(\mathbf{y},t):\mathbb{R}^{2}\times\mathbb{R}\rightarrow\mathbb{R}^{3}$
that describes the screen-space motion of the observed vocal fold surface from
the experimental setup described in Section 2. We combine this vector field
later with a strong shape prior of the vocal folds using the finite-element
method (FEM).

$6.5\;$\mathrm{m}\mathrm{s}$$

$6.6\;$\mathrm{m}\mathrm{s}$$

$6.7\;$\mathrm{m}\mathrm{s}$$

$6.8\;$\mathrm{m}\mathrm{s}$$

$6.9\;$\mathrm{m}\mathrm{s}$$

$7.0\;$\mathrm{m}\mathrm{s}$$

$7.1\;$\mathrm{m}\mathrm{s}$$

$7.2\;$\mathrm{m}\mathrm{s}$$

$7.3\;$\mathrm{m}\mathrm{s}$$

$7.4\;$\mathrm{m}\mathrm{s}$$

$7.5\;$\mathrm{m}\mathrm{s}$$

$7.6\;$\mathrm{m}\mathrm{s}$$

$7.7\;$\mathrm{m}\mathrm{s}$$

$7.8\;$\mathrm{m}\mathrm{s}$$

Figure 3:  Time series of images of the superior surface of the vocal folds
captured by a high-speed camera.

#### High-speed Camera

The high-speed camera provides a sequence of mono-chromatic images, as shown
in Fig. 3. These images can be described as a function of intensities
$I(\mathbf{y},t_{i}):\mathbb{R}^{2}\times\mathbb{N}_{0}\rightarrow\mathbb{R}$
with an image resolution $w\times h$, the screen-space coordinates
$\mathbf{y}\in[0,w-1]\times[0,h-1]\subset\mathbb{R}^{2}$ and with
$t_{i}\in[0,F)\subset\mathbb{N}_{0}$ being the currently sampled frame of the
video with $F$ many frames in total. Then, $u_{I}(\mathbf{y},t)$ and
$v_{I}(\mathbf{y},t)$ denote the estimated motion of the screen-space
coordinate $\mathbf{y}$ at time $t$ along the longitudinal and lateral axis,
respectively, resulting in the motion vector field:

$\displaystyle\mathbf{v}_{I}(\mathbf{y},t)=\begin{pmatrix}u_{I}(\mathbf{y},t)\\\
v_{I}(\mathbf{y},t)\end{pmatrix}.$ (1)

We applied the sparse to dense optical flow algorithm based on the Lucas-
Kanade feature tracker with pyramids by [25], which is available in OpenCV
[26]. Such dense optical flow algorithm estimates a motion for every pixel of
the input image, but tends to oversmooth flows (see Fig. 4(a)), which produces
artifacts at edges. To mitigate this problem, we apply a mask (see Fig. 4(b)),
which helps to distinguish the vocal folds from the gap between them. We use
the intensity values for the mask segmentation with a subsequent median filter
for more robustness to noise. Using the mask (see Fig. 4(c)), we set the
motion vectors inside the gap explicitly to zero vectors, since there is no
vocal fold movement. For consistency with the remaining data, the motion in
the image plane is converted to physical domain units. The high-speed images
and the resulting optical flows are discrete in space (individual pixels) and
in time (individual frames). We interpolate the resulting discrete vector
piecewise bilinearly in space with $C^{0}$ continuity and we interpolate
quadratically with $C^{1}$ continuity in time to extract a continuous time-
dependent two-dimensional vector field.

(a) Oversmoothed optical flow

(b) Binary silhouette mask

(c)  Applied silhouette mask to optical flow

Figure 4:  (a) OpenCV’s [26] sparse to dense optical flow algorithm
interpolates a flow in the gap, (b) shows the binary silhouette mask that
identifies the gap, (c) shows the cleaned optical flow with velocities set to
zero in the gap.

#### Vibrometer Measurements

Vibrometer measurements are taken at discrete grid points (see Fig. 5) with
screen-space coordinates $\mathbf{y}_{i}$, which remain the same across all
time steps. For each measurement point, we receive an associated velocity
component $w_{i}(\mathbf{y}_{i},t_{i})$ along the vertical axis. Unlike the
optical flow vector fields, the measurement points cannot be interpolated
bilinearly in space, because the measurement grid is not exactly regular.
Instead, we interpolate spatially using radial basis functions:

$\displaystyle
w(\mathbf{y},t_{i})=\sum_{i}\omega_{i}\cdot\phi(\|\mathbf{y}-\mathbf{y}_{i}\|)$
(2)

where $\phi(r)=r^{2}\ln(r)$ is the thin plate spline kernel and the weights
$\omega_{i}$ are linearly fitted per time step $t_{i}$ to meet all constraints
in a least-squares sense, i.e.,
$\sum_{i}\|w(\mathbf{y}_{i},t_{i})-w_{i}(\mathbf{y}_{i},t_{i})\|^{2}\rightarrow\min$.
Formally, we interpolate the vibrometer measurements linearly in time to
obtain a continuous field $w(\mathbf{y},t)$. With this, $w(\mathbf{y},t)$
gives us the velocity for an arbitrary screen-space coordinate $\mathbf{y}$ at
time $t$. Using the optical flow $\mathbf{v}_{I}(\mathbf{y},t)$ from Eq. (1),
we combine the longitudinal and lateral velocity components with the vertical
component from the vibrometer in Eq. (2) to obtain the tri-variate surface
vector field $\mathbf{v}(\mathbf{y},t)$ that describes the observed three-
dimensional motion of the vocal fold surface:

$\displaystyle\mathbf{v}(\mathbf{y},t)=\begin{pmatrix}u_{I}(\mathbf{y},t)\\\
v_{I}(\mathbf{y},t)\\\ \,w(\mathbf{y},t)\end{pmatrix}$ (3)

In the following, we optimize for a vocal fold model using FEM simulations
that match the observed motion.

Figure 5:  Schematic view of the vibrometer measurement points marked in red
on the vocal fold model.

###  3.3 Vocal Fold Model

We parameterized the inferior and superior surface of the M5 vocal fold model
with $(u,v)$ coordinates. Thus, given is the M5 vocal fold model with a
surface $\mathbf{x}_{0}(u,v):\mathbb{R}^{2}\rightarrow\mathbb{R}^{3}$,
parameterized by $(u,v)$ coordinates. For each time $t$, the model is subject
to a stationary FEM simulation, resulting in a deformed shape
$\overline{\mathbf{x}}(u,v,t)$. Let $\mathbf{f}(u,v,t)$ be the unknown time-
dependent force field acting on the surface. We compose the force field from a
tensor-product with unknown control forces $\mathbf{f}_{i,j,k}$. Thereby, the
polynomial degrees are $l$, $m$, and $n$, and the Bézier Bernstein basis
functions are $B_{i}^{n}(t)=\binom{n}{i}\;t^{i}\;(1-t)^{n-i}$. A force field
$\mathbf{f}(u,v,t)$ is given as tensor product via:

$\displaystyle\mathbf{f}(u,v,t)=\sum_{i=0}^{l}\sum_{j=0}^{m}\sum_{k=0}^{n}\mathbf{f}_{i,j,k}\;B_{i}^{l}(u)\;B_{j}^{m}(v)\;B_{k}^{n}(t)$
(4)

Discontinuities in the force field are required due to discontinuities in the
observable motion. For example, a discontinuity in the velocity along the
lateral axis can be observed when the vocal folds touch and are pushed apart
by the air flow. Thus, we separate our force field into 10 individual scalar
fields to have greater control over each individual component and to model
discontinuities in our force field. We employ five of them to act on the
inferior surface and five of them on the superior surface of the vocal folds.
Two of the five scalar fields represent the longitudinal component of the
force, two scalar fields represent the lateral component of the force and one
scalar field represents the vertical component of the force. The scalar fields
for the longitudinal and lateral axis cover only half of the vocal folds as
shown in Fig. 6.

(a) Scalar fields for lateral component

(b) Scalar fields for longitudinal component

(c) Scalar field for vertical component

Figure 6:  To model discontinuities, the force field is split in the lateral
and longitudinal direction into two pieces with $C^{0}$ continuity. a) shows
the two lateral pieces, b) shows the two longitudinal pieces, and c) shows
that there is only one piece used for the vertical component. The same layout
is used on the superior and inferior side of the vocal folds.

#### Finite Element Method

We utilize FEM simulations for the deformation of the initial vocal fold
model. We discretize the vocal fold model and apply the sampled forces from
the force field function $\mathbf{f}$ to the inferior and superior areas of
the model. The FEM simulation requires Lamé constants $\lambda$ and $\mu$,
which determine the stretching and stiffness of the material. $\lambda$ and
$\mu$ can be calculated using the Poisson ratio $\nu$ and the Young’s modulus
$E$ of the used elastomer:

$\displaystyle\lambda$
$\displaystyle=\frac{E\nu}{\left(1+\nu\right)\left(1-2\nu\right)}$
$\displaystyle\mu$ $\displaystyle=\frac{E}{2\left(1+\nu\right)}$ (5)

We define the stress
$\sigma(\mathbf{x}):\mathbb{R}^{3}\rightarrow\mathbb{R}^{3\times 3}$ as tensor
field for position $\mathbf{x}$ and denote the resulting displacement vector
field as $\mathbf{u}(\mathbf{x}):\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$. The
FEM simulation solves the linear elasticity partial differential equation
after the forces are applied with $\mathbf{I}$ being the identity matrix and
$\mathbf{f}(\mathbf{x}):\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ being the
externally applied forces:

$\displaystyle\mathbf{f}(\mathbf{x})+\left(\nabla\cdot\sigma\left(\mathbf{x}\right)^{\textrm{T}}\right)^{\textrm{T}}$
$\displaystyle=0$ (6) $\displaystyle\sigma\left(\mathbf{x}\right)$
$\displaystyle=\lambda\;\textrm{div}\left(\mathbf{u}(\mathbf{x})\right)\mathbf{I}+\mu\left(\nabla\mathbf{u}(\mathbf{x})+\nabla\mathbf{u}(\mathbf{x})^{\textrm{T}}\right)$
(7)

The resulting displacement field creates a deformation of the vocal folds (see
Fig. 7(b)), where the externally applied forces are in an equilibrium with the
internal forces of the deformed material. We additionally set boundary
conditions on the sides that are in contact with the fixture. These sides are
immovable, because the sides of the model are glued to the fixture in the
experimental setup and therefore cannot move.

(a) Rest state (before applying forces)

(b) Equilibrium state (after applying forces)

Figure 7: Comparison of one vocal fold before and after the FEM simulation,
where the displacement magnitude is mapped to color. (a) shows the initial
rest position without forces applied to it, and (b) shows the vocal fold in
its equilibrium state after the forces were applied.

#### Vocal Fold Motion

Our goal is to match the motion of the simulated vocal fold with the motion
observed in the measurement data. Next, we describe how the motion of the
simulated vocal fold is calculated in screen-space, which is the space in
which simulation and measurement are compared. For the screen-space
coordinates $\mathbf{y}$, we determine the $(u,v)$ coordinates of the visible
surface of the deformed model $\overline{\mathbf{x}}(u,v,t)$ by ray casting.
Thereby, the binary indicator $m(\mathbf{y})\in\\{0,1\\}$ is created, which
denotes whether the ray cast was successful, i.e., whether a vocal fold
surface is visible in the pixel at screen-space coordinate $\mathbf{y}$. The
screen-space motion is then determined by finite differences:

$\displaystyle\overline{\mathbf{v}}(\mathbf{y},t):=\frac{\mathrm{d}\overline{\mathbf{x}}}{\mathrm{d}t}(u,v,t)\approx\frac{\overline{\mathbf{x}}(u,v,t+\Delta
t)-\overline{\mathbf{x}}(u,v,t)}{\Delta t}$ (8)

where the mesh points with the same $(u,v)$ coordinates are taken, namely the
coordinates observed at time $t$ for pixel $\mathbf{y}$. With this, we observe
where a mesh point with a certain $(u,v)$ coordinate moves to over time in
screen-space.

#### Gap

The gap is the opening between the two individual vocal fold halves. The gap
is a defining feature of the vocal fold motion and we thus aim to preserve the
gap in our FEM-based reconstruction. To measure the extent of the gap, we
introduce the gap opening function
$o(\mathbf{y},t):\mathbb{R}^{2}\times\mathbb{R}\rightarrow\mathbb{R}^{+}$,
from which the silhouette of the gap can be determined, as illustrated in Fig.
8. For a given screen-space coordinate $\mathbf{y}$ at time $t$ the gap
opening function $o(\mathbf{y},t)$ denotes the unsigned lateral distance to
the closest position on the gap silhouette, which is determined by casting a
ray. If the ray does not intersect the gap, then the distance to the last
known opening location is taken.

Figure 8:  The gap opening function $o(\mathbf{y},t)$ returns for a given
screen-space coordinate $\mathbf{y}$ (red) the unsigned distance $d$ to the
closest point on the silhouette (blue) along the lateral axis.

###  3.4 Optimization

In the following, we describe the energy minimization approach that determines
the forces to apply in the FEM simulation. We initialize the forces with zero
and optimize them iteratively according to the error metric.

#### Error Metric

We minimize an error $\mathrm{E}$, which penalizes a mismatch in the
reproduced motion and a difference in the silhouette of the gap. For the
latter, we apply an arctangent to enhance the penalty on smaller differences.
We integrate the error over a selected time domain, from $T_{1}$ to $T_{2}$.

$\displaystyle\mathrm{E}=\frac{1}{T_{2}-T_{1}}\int_{T_{1}}^{T_{2}}\frac{1}{N}\sum_{i=1}^{N}\underbrace{m(\mathbf{y}_{i})\;\|\mathbf{v}\left(\mathbf{y}_{i},t\right)-\overline{\mathbf{v}}(\mathbf{y}_{i},t)\|^{2}}_{\textrm{motion}}+\underbrace{\arctan\left(o(\mathbf{y}_{i},t)\right)}_{\textrm{gap}}\;\mathrm{d}t$
(9)

The motion term of the error sums up the $L^{2}$ norm of the difference in the
measured velocity $\mathbf{v}$ from Eq. (3) and the simulated velocity
$\overline{\mathbf{v}}$ from Eq. (8) over all measurement points, if a
measurement point hits both the deformed model as well as the model in the
experiment ($m(\mathbf{y}_{i})=1$). The measurement point is ignored if the
deformed model and the ground truth measurement point do not have a valid
velocity, i.e., when $m(\mathbf{y}_{i})=0$. If only the vocal folds in the
experiment have a valid velocity for the measurement point, then the $L^{2}$
norm of the ground truth velocity is added. The second term of the error
metric aims to match the silhouette of the gap. We sum up the unsigned
distance returned from the gap opening function $o(\mathbf{y}_{i},t)$ after we
adjusted their influence using the $\arctan$ function. Finally, the
accumulated error is normalized by the number of measurement points $N$.

#### Gradient Descent

The error metric $\mathrm{E}$ is evaluated in each optimization step for a
given set of control forces $\mathbf{f}_{i,j,k}$. Finite differences are
applied to each individual control point of the force field to construct the
gradient of the error metric. Afterwards, gradient descent is used to find the
control forces that minimize the error:

$\displaystyle\mathbf{f}_{i,j,k}^{n+1}=\mathbf{f}_{i,j,k}^{n}-\gamma\frac{\partial}{\partial\mathbf{f}_{i,j,k}}\mathrm{E}(\mathbf{f}_{i,j,k}^{n})$
(10)

A gradient descent with step size $\gamma$ is able to find a local minimum of
a function by descending along the steepest descent direction. To accelerate
convergence, we perform a line search with increasing step size.

$3\times 3\times 10$ (one-sided forces) 🌑

$3\times 3\times 10$ (two-sided forces) 🌑

$5\times 5\times 15$ (two-sided forces) 🌑

$t_{0}$

$t_{150}$

$t_{310}$

$t_{450}$

$t_{499}$

Figure 9: Visual comparison of the three setups for different time steps of
the animation. Left: optimization for a force field with $3\times 3\times 10$
control points (longitudinal, lateral, time) on the inferior side only.
Middle: Adding a second force field on the superior side allows for a better
capture of the bending of vocal folds. Right: Increasing the resolution to
$5\times 5\times 15$ allows for bending while maintaining symmetry. The
displacement magnitude is color-coded.

##  4 Results

In the following section, we show the results of our vocal fold surface
reconstruction, and we compare the motion of the FEM model with the velocity
measurements. The unknown in our optimization is the force field $\mathbf{f}$
that is applied to the vocal fold surface in the FEM simulation. To validate
our choice of optimizing for a force field on both the inferior and the
superior side of the vocal model, we performed two experiments: first, we
applied the force field only on the inferior side, and second we applied the
force field on both sides. Another degree of freedom is the number of control
points that are used to model the force field as a tensor product. Here, we
compared a $3\times 3\times 10$ resolution with a $5\times 5\times 15$
resolution (longitudinal, lateral, time).

#### Time-Series Comparison with Optical Measurements

In Figure 9, the optimized FEM models are overlayed on top of the optical
high-speed camera measurements. The color encodes the displacement magnitude
$\|\mathbf{u}\|$. The vocal fold shape obtained from the two-sided force field
optimization matches the optical measurements better. In particular, the
opening of the gap and the folding of the vocal folds match the optical
measurements more closely. In the simulation results, however, an asymmetry in
the vocal folds is apparent, which could indicate a pathology. The asymmetry
reduces with increasing tensor product resolution. Closer inspection of the
high-speed camera footage revealed that an asymmetry is also present in the
experimental data, since the vocal folds were not glued in perfectly
symmetric.

#### Convergence

The error convergence plots of the gradient descent schemes are shown for all
experiments in Figure 10, where the optimization was terminated once the error
minimization converged. While the one-sided force field optimization converged
after about 25 iterations, the two-sided force field optimization converged
after 85 iterations and 115 iterations, respectively. Both two-sided force
field optimizations had a similar error decay. The model with higher
resolution, however, was able to reduce the error a bit further before
reaching convergence.

#### Agreement with Observed Velocities for Selected Surface Points

Aside from matching the silhouette of the gap, the objective function aims to
reconstruct the observed velocities on the superior side, considering both the
vibrometer measurement (vertical component) and the optical flow (horizontal
components). In Figure 11, we show the velocity signal over time for two
selected points on the vocal fold surface. The velocity signals are shown for
the vertical, longitudinal, and lateral component, separately. For the first
observation point (left column), the two-sided force fields were closer to the
measured data than the one-sided model. An increase in resolution of the force
field tensor product resulted in marginally lower errors, since the fold had
more degrees of freedom to match high-frequency fluctuations. Note, however,
that an exact reproduction of all high-frequent details is not desirable,
since some can be considered measurement noise. The second observation point
(right column) exhibited for the longitudinal component a temporal variation
that could not be matched by any of the three optimizations, due to an
insufficient number of degrees of freedom. The two-sided high-resolution
optimization reproduced the lateral component noticeably better, whereas the
differences in the vertical components were marginally improved.

Figure 10: Convergence plots of error E over the course of the gradient
descent minimization. With one-sided forces, the error reaches a plateau after
25 gradient descent iterations. With two-sided forces, the error can be
reduced further.

$\mathbf{y}_{i}=(7.45,4.54)$

$\mathbf{y}_{i}=(12.73,15.25)$

Figure 11: Plots of vertical, longitudinal, and lateral velocity components
(in rows) for vibrometer measurements  _🌑_ , our $3\times 3\times 10$ (one-
sided forces) _🌑_ , $3\times 3\times 10$ (two-sided forces) _🌑_ , and $5\times
5\times 15$ (two-sided forces) _🌑_ for two different selected image
coordinates (in columns). An increase in force field resolution and a two-
sided force field application result in the lowest error. Residuals are
visible in particular for high-frequent temporal changes due to insufficient
degrees of freedom in the optimization.

##  5 Conclusion

In this paper, we proposed a novel FEM-based approach to reconstruct the
three-dimensional motion of a complete vocal fold geometry based on optical
measurements that were taken from the superior side only. First, we
reconstructed a surface velocity field from the optical flow of a high-speed
camera recording (horizontal components) and a laser Doppler vibrometer
(vertical component). We then applied a time-dependent force field in tensor
product basis to the inferior and/or superior side of the vocal model and let
an FEM simulation deform the surface into an equilibrium state. The force
field control points are then optimized via gradient descent to minimize an
error metric that consisted of two terms. The first term measures the
difference between simulated velocities and observed measurements, while the
second term measures the silhouette difference between observation and
simulation. The results showed that our approach is able to qualitatively
reproduce the underlying ground truth data of surface velocity and glottal
gap. Increasing the resolution of the model resulted in a further reduction of
the error. We saw that for the surface velocities, our models were able to
reproduce the measurement data roughly. They were, however, not able to
reproduce high-frequency components of the data. This has on the one hand the
advantage, that it filters out measurement noise; on the other hand, this
might also cause an unwanted low-pass filtering of the physical vocal fold
surface velocities.

In future works, the behavior of our model in the frequency domain could be
analyzed by optimization over several oscillation periods, followed by a
transformation of the velocity data into the frequency domain. Furthermore, an
increased accuracy might be achieved by a further increase of the model
resolution. In both cases, an improved computation of the gradient needed for
optimization like automatic differentiation would benefit the procedure by
enhancing the computational efficiency of the numerical optimization,
rendering more complex models possible. As another step, also dynamic effects
of the vocal folds caused by e.g. viscoelasticity could be incorporated.

### Acknowledgements

This work is funded by the German Research Foundation (DFG) through the
projects ”Tracing the mechanisms that generate tonal content in voiced speech”
(project number: 446965891) and "End-to-End Optimization for Energy-Driven
Scientific Visualization" (project number: 517157369).

## References

  * [1] Ray M. Merrill, Nelson Roy, and Jessica Lowe.  Voice-related symptoms and their effects on quality of life.  Annals of Otology, Rhinology & Laryngology, 122(6):404–411, jun 2013.
  * [2] M. Kaltenbacher, S. Zörner, and A. Hüppe.  On the importance of strong fluid-solid coupling with application to human phonation.  Progress in Computational Fluid Dynamics, An International Journal, 14(1):2, 2014.
  * [3] Rita R. Patel, Denis Dubrovskiy, and Michael Döllinger.  Measurement of glottal cycle characteristics between children and adults: Physiological variations.  Journal of Voice, 28(4):476–486, jul 2014.
  * [4] Matthias Echternach, Michael Döllinger, Johan Sundberg, Louisa Traser, and Bernhard Richter.  Vocal fold vibrations at high soprano fundamental frequencies.  The Journal of the Acoustical Society of America, 133(2):EL82–EL87, feb 2013.
  * [5] Michael Döllinger and David A. Berry.  Computation of the three-dimensional medial surface dynamics of the vocal folds.  Journal of Biomechanics, 39(2):369–374, 2006.
  * [6] David E. Sommer, Isao T. Tokuda, Sean D. Peterson, Ken-Ichi Sakakibara, Hiroshi Imagawa, Akihito Yamauchi, Takaharu Nito, Tatsuya Yamasoba, and Niro Tayama.  Estimation of inferior-superior vocal fold kinematics from high-speed stereo endoscopic data in vivo.  The Journal of the Acoustical Society of America, 136(6):3290–3300, dec 2014.
  * [7] Georg Luegmair, Stefan Kniesburges, Maik Zimmermann, Alexander Sutor, Ulrich Eysholdt, and Michael Döllinger.  Optical reconstruction of high-speed surface dynamics in an uncontrollable environment.  IEEE Transactions on Medical Imaging, 29(12):1979–1991, dec 2010.
  * [8] Rita R. Patel, Kevin D. Donohue, Daniel Lau, and Harikrishnan Unnikrishnan.  In vivo measurement of pediatric vocal fold motion using structured light laser projection.  Journal of Voice, 27(4):463–472, jul 2013.
  * [9] Georg Luegmair, Daryush D. Mehta, James B. Kobler, and Michael Dollinger.  Three-dimensional optical reconstruction of vocal fold kinematics using high-speed video with a laser projection system.  IEEE Transactions on Medical Imaging, 34(12):2572–2582, dec 2015.
  * [10] Marion Semmler, Stefan Kniesburges, Veronika Birk, Anke Ziethe, Rita Patel, and Michael Dollinger.  3d reconstruction of human laryngeal dynamics based on endoscopic high-speed recordings.  IEEE Transactions on Medical Imaging, 35(7):1615–1624, jul 2016.
  * [11] Jann-Ole Henningson, Marc Stamminger, Michael Döllinger, and Marion Semmler.  Real-time 3D reconstruction of human vocal folds via high-speed laser-endoscopy.  In Lecture Notes in Computer Science, pages 3–12. Springer Nature Switzerland, 2022.
  * [12] Ronald C. Scherer, Daoud Shinwari, Kenneth J. De Witt, Chao Zhang, Bogdan R. Kucinschi, and Abdollah A. Afjeh.  Intraglottal pressure profiles for a symmetric and oblique glottis with a divergence angle of 10 degrees.  The Journal of the Acoustical Society of America, 109(4):1616–1630, apr 2001.
  * [13] Byron D. Erath, Sean D. Peterson, Matías Zañartu, George R. Wodicka, and Michael W. Plesniak.  A theoretical model of the pressure field arising from asymmetric intraglottal flows applied to a two-mass model of the vocal folds.  The Journal of the Acoustical Society of America, 130(1):389–403, jul 2011.
  * [14] X. Zheng, R. Mittal, Q. Xue, and S. Bielamowicz.  Direct-numerical simulation of the glottal jet and vocal-fold dynamics in a three-dimensional laryngeal model.  The Journal of the Acoustical Society of America, 130(1):404–415, jul 2011.
  * [15] S. Zörner, M. Kaltenbacher, and M. Döllinger.  Investigation of prescribed movement in fluid–structure interaction simulation for the human phonation process.  Computers & Fluids, 86:133–140, nov 2013.
  * [16] Hossein Sadeghi, Stefan Kniesburges, Manfred Kaltenbacher, Anne Schützenberger, and Michael Döllinger.  Computational models of laryngeal aerodynamics: Potentials and numerical costs.  Journal of Voice, 33(4):385–400, jul 2019.
  * [17] Stefan Schoder, Michael Weitz, Paul Maurerlehner, Alexander Hauser, Sebastian Falk, Stefan Kniesburges, Michael Döllinger, and Manfred Kaltenbacher.  Hybrid aeroacoustic approach for the efficient numerical simulation of human phonation.  The Journal of the Acoustical Society of America, 147(2):1179–1194, feb 2020.
  * [18] Sebastian Falk, Stefan Kniesburges, Stefan Schoder, Bernhard Jakubaß, Paul Maurerlehner, Matthias Echternach, Manfred Kaltenbacher, and Michael Döllinger.  3d-FV-FE aeroacoustic larynx model for investigation of functional based voice disorders.  Frontiers in Physiology, 12, mar 2021.
  * [19] Paul Maurerlehner, Stefan Schoder, Clemens Freidhager, Andreas Wurzinger, Alexander Hauser, Florian Kraxberger, Sebastian Falk, Stefan Kniesburges, Matthias Echternach, Michael Döllinger, and Manfred Kaltenbacher.  Efficient numerical simulation of the human voice.  e & i Elektrotechnik und Informationstechnik, 138(3):219–228, apr 2021.
  * [20] Stefan Schoder, Paul Maurerlehner, Andreas Wurzinger, Alexander Hauser, Sebastian Falk, Stefan Kniesburges, Michael Döllinger, and Manfred Kaltenbacher.  Aeroacoustic sound source characterization of the human voice production-perturbed convective wave equation.  Applied Sciences, 11(6):2614, mar 2021.
  * [21] Alexander Lodermeyer, Stefan Becker, Michael Döllinger, and Stefan Kniesburges.  Phase-locked flow field analysis in a synthetic human larynx model.  Experiments in Fluids, 56(4), mar 2015.
  * [22] F Durst, U Heim, B nsal, and G Kullik.  Mass flow rate control system for time-dependent laminar and turbulent flow investigations.  Measurement Science and Technology, 14(7):893–902, may 2003.
  * [23] Zhaoyan Zhang, Juergen Neubauer, and David A. Berry.  The influence of subglottal acoustics on laboratory models of phonation.  The Journal of the Acoustical Society of America, 120(3):1558–1569, sep 2006.
  * [24] Stefan J. Rupitsch, Jürgen Ilg, Alexander Sutor, Reinhard Lerch, and Michael Döllinger.  Simulation based estimation of dynamic mechanical properties for viscoelastic materials used for vocal fold models.  Journal of Sound and Vibration, 330(18-19):4447–4459, aug 2011.
  * [25] Jean-Yves Bouguet et al.  Pyramidal implementation of the affine lucas kanade feature tracker description of the algorithm.  Intel corporation, 5(1-10):4, 2001.
  * [26] G. Bradski.  The OpenCV Library.  Dr. Dobb’s Journal of Software Tools, 2000.


# BrainMetDetect: Predicting Primary Tumor from Brain Metastasis MRI Data
Using Radiomic Features and Machine Learning Algorithms

Hamidreza Sadeghsalehi
Iran University of Medical Sciences
<EMAIL_ADDRESS>

###### Abstract

Objective: Brain metastases (BMs) are common in cancer patients and
determining the primary tumor site is crucial for effective treatment. This
study aims to predict the primary tumor site from BM MRI data using radiomic
features and advanced machine learning algorithms.
Methods: We utilized a comprehensive dataset from Ocaña-Tienda et al. (2023)
comprising MRI and clinical data from 75 patients with BMs. Radiomic features
were extracted from post-contrast T1-weighted MRI sequences. Feature selection
was performed using the GINI index, and data normalization was applied to
ensure consistent scaling. We developed and evaluated Random Forest and
XGBoost classifiers, both with and without hyperparameter optimization using
the FOX (Fox optimizer) algorithm. Model interpretability was enhanced using
SHAP (SHapley Additive exPlanations) values.
Results: The baseline Random Forest model achieved an accuracy of 0.85, which
improved to 0.93 with FOX optimization. The XGBoost model showed an initial
accuracy of 0.96, increasing to 0.99 after optimization. SHAP analysis
revealed the most influential radiomic features contributing to the models’
predictions. The FOX-optimized XGBoost model exhibited the best performance
with a precision, recall, and F1-score of 0.99.
Conclusion: This study demonstrates the effectiveness of using radiomic
features and machine learning to predict primary tumor sites from BM MRI data.
The FOX optimization algorithm significantly enhanced model performance, and
SHAP provided valuable insights into feature importance. These findings
highlight the potential of integrating radiomics and machine learning into
clinical practice for improved diagnostic accuracy and personalized treatment
planning.

_K_ eywords Brain metastasis $\cdot$ MRI $\cdot$ Radiomics $\cdot$ Machine
learning $\cdot$ FOX optimization $\cdot$ Random Forest $\cdot$ XGBoost

##  1 Introduction

Brain metastases (BMs) are a significant medical challenge, representing the
most common type of intracranial tumor in adults, with the incidence rising
due to improved systemic control of primary cancers and increased longevity of
cancer patients [1]. Early and accurate identification of the primary tumor
site is crucial for determining the appropriate treatment strategy and
improving patient outcomes. Despite advancements in neuroimaging techniques,
differentiating the primary tumor origin based solely on brain metastasis
imaging remains complex and elusive [2].

Radiomics, a developing discipline that produces high-dimensional, mineable
data from medical images, has shown promise in enhancing diagnostic accuracy
in oncology [3]. Radiomic features capture the underlying heterogeneity of
tumors that may not be visually discernible, thus providing a robust
quantitative basis for disease characterization [4, 5]. Integrating these
features with advanced machine learning algorithms offers a compelling
approach to predict the primary tumor site from brain metastasis MRI data.

This study leverages a comprehensive dataset of clinical information and
annotated brain metastases MR images, as described by Ocaña-Tienda et al.
(2023) [6]. The dataset includes imaging studies from 75 patients diagnosed
with BMs, collected from five distinct medical institutions, ensuring a
diverse and representative sample. The clinical data encompass various patient
demographics, treatment regimens, and survival statistics, while the imaging
data consist of high-resolution post-contrast T1-weighted MRI sequences,
crucial for accurate feature extraction and analysis.

Our research aims to develop and evaluate the performance of several machine
learning models for predicting the primary tumor site using radiomic features
extracted from BM MRI data. We employ Random Forest and XGBoost classifiers,
both with and without hyperparameter optimization using the FOX (Fox
optimizer) algorithm. The FOX algorithm, inspired by the foraging behavior of
foxes, has demonstrated superior performance in optimizing complex functions
and is particularly suited for enhancing the predictive accuracy of our
models. Additionally, we utilize SHAP (SHapley Additive exPlanations) to
interpret and elucidate the contributions of individual radiomic features to
the model predictions, providing valuable insights into the clinical relevance
of these features.

Previous studies have explored various machine learning techniques for cancer
prediction and classification using radiomic features. For instance, Parmar et
al. (2015) demonstrated the utility of radiomics in predicting overall
survival in lung cancer patients using machine learning models [7]. Similarly,
Macyszyn et al. (2015) employed radiomic features to predict the molecular
subtypes of glioblastomas, showcasing the potential of this approach in neuro-
oncology [8]. However, the application of radiomics combined with advanced
optimization algorithms like FOX for predicting primary tumors from brain
metastasis MRI data remains underexplored.

In this study, we aim to fill this gap by systematically evaluating the
effectiveness of radiomic features and optimized machine learning models in
classifying primary tumor sites. Our work contributes to the growing body of
literature on radiomics and machine learning in medical imaging, offering a
novel approach that enhances predictive accuracy and model interpretability.
The findings of this research have significant implications for clinical
practice, potentially aiding in the development of personalized treatment
plans and improving patient outcomes in cases of brain metastasis.

##  2 Materials and Methods

###  2.1 Dataset

Source: The dataset used in this study is from Ocaña-Tienda et al. (2023),
titled "A comprehensive dataset of annotated brain metastasis MR images with
clinical and radiomic data." This dataset is publicly available and can be
accessed at https://doi.org/10.1038/s41597-023-02123-0 [6].

Description:

  * •

Patient Selection: The dataset comprises 75 patients diagnosed with brain
metastases (BMs), with data collected from five different medical
institutions. Patients included in the dataset were adults with a
pathologically confirmed diagnosis of BM, and data spanned from January 1,
2005, to December 31, 2021.

  * •

Imaging Data: The imaging data includes post-contrast T1-weighted high-
resolution MRI sequences with pixel spacing $\leq$ 2 mm and slice thickness
$\leq$ 2 mm. Only images devoid of noise or artifacts were included, ensuring
high-quality imaging data.

  * •

Clinical Data: Clinical data encompass age at diagnosis, gender, treatment
regimens, and survival statistics.

  * •

Primary Tumor Distribution: The dataset includes various primary tumor types:
non-small cell lung cancer (NSCLC) (38 cases), small cell lung cancer (SCLC)
(5 cases), breast cancer (22 cases), melanoma (6 cases), ovarian cancer (2
cases), kidney cancer (1 case), and uterine cancer (1 case).

  * •

Brain Metastases Data: A total of 75 patients with 260 distinct BMs were
included, with 637 imaging studies, 593 of which underwent semi-automatic
segmentation.

###  2.2 Data Preprocessing

Feature Selection: Feature selection was performed using the GINI index to
enhance model performance and reduce dimensionality. The GINI index is a
measure of impurity used in decision tree-based algorithms. In this study, the
GINI index was employed to evaluate the importance of each radiomic feature in
predicting the primary tumor from BM MRI data. The top 50 radiomic features
with the highest GINI index scores were selected for further analysis [9, 10].

Data Normalization: After feature selection, data normalization was performed
to ensure all features were on a consistent scale. Normalization techniques
such as z-score normalization (standardization) and min-max scaling were
applied, transforming the data such that each feature had a mean of 0 and a
standard deviation of 1 (z-score normalization) or a range between 0 and 1
(min-max scaling) [11]. This step was crucial to prevent features with larger
numerical ranges from disproportionately influencing the machine learning
models.

Data Splitting: The normalized dataset was split into training and test sets
with an 80-20 split ratio. This allowed for model training on a subset of the
data while preserving a separate subset for evaluating model performance [12].

###  2.3 Machine Learning Models

Random Forest:

  * •

Baseline Model: A Random Forest classifier was initially used without any
hyperparameter optimization to establish a baseline performance for primary
tumor classification [9].

  * •

FOX Optimized Model: The Random Forest classifier was then optimized using the
FOX (Fox optimizer) algorithm to improve performance. The FOX algorithm is a
nature-inspired optimization technique that mimics the foraging behavior of
foxes, using distance measurements between the fox and its prey to execute
efficient jumps and enhance the model’s performance [13].

XGBoost:

  * •

Baseline Model: An XGBoost (Extreme Gradient Boosting) classifier was also
applied without hyperparameter tuning to evaluate its initial performance
[14].

  * •

FOX Optimized Model: The XGBoost classifier underwent hyperparameter tuning
using the FOX algorithm to further enhance its predictive accuracy.

###  2.4 FOX Optimization Algorithm

Overview: The FOX (Fox optimizer) algorithm is a novel nature-inspired
optimization technique inspired by the foraging behavior of foxes. It measures
the distance between a fox and its prey to execute efficient jumps, optimizing
the process [15].

Application: In this study, the FOX algorithm was applied to fine-tune the
hyperparameters of both the Random Forest and XGBoost models. Specifically,
the algorithm searched for the optimal hyperparameters that maximized the
models’ predictive performance. The FOX algorithm’s effectiveness was
initially validated using benchmark functions and the CEC2019 benchmark test
functions, demonstrating superior performance compared to other optimization
algorithms such as Particle Swarm Optimization (PSO), Genetic Algorithm (GA),
and Grey Wolf Optimization (GWO) [16].

###  2.5 Model Explainability with SHAP

SHAP (SHapley Additive exPlanations): SHAP values provide a unified measure of
feature importance by assigning each feature a value representing its
contribution to the model’s prediction [17].

Application: SHAP was used to interpret the machine learning models’
predictions by explaining the impact of each radiomic feature on the
classification outcome. This approach helped to identify the most influential
features, providing insights into the clinical relevance of the radiomic
features used in the study. SHAP summary plots and bar charts were generated
to visualize the contributions of individual features to the model’s
predictions.

##  3 Results

###  3.1 Performance Metrics

The performance of each machine learning model in classifying primary tumors
from brain metastasis MRI data was evaluated using several metrics: precision,
recall, F1-score, and accuracy. These metrics provide a comprehensive
understanding of the models’ effectiveness.

###  3.2 Random Forest Results

Without Optimization: The baseline Random Forest classifier was implemented
without hyperparameter optimization. The performance metrics are as follows:

  * •

Precision: 0.86 (weighted average)

  * •

Recall: 0.85 (weighted average)

  * •

F1-score: 0.84 (weighted average)

  * •

Accuracy: 0.85

Figure 1: Confusion matrix for baseline Random Forest model.

With FOX Optimization: The Random Forest classifier was then optimized using
the FOX algorithm. The optimized model showed significant improvements in
performance:

  * •

Precision: 0.94 (weighted average)

  * •

Recall: 0.93 (weighted average)

  * •

F1-score: 0.93 (weighted average)

  * •

Accuracy: 0.93

Figure 2: Confusion matrix for FOX-optimized Random Forest model.

###  3.3 XGBoost Results

Without Tuning: The XGBoost classifier was applied without hyperparameter
tuning, yielding the following performance metrics:

  * •

Precision: 0.96 (weighted average)

  * •

Recall: 0.96 (weighted average)

  * •

F1-score: 0.96 (weighted average)

  * •

Accuracy: 0.96

Figure 3: Confusion matrix for baseline XGBoost model.

With FOX Optimization: Hyperparameter tuning of the XGBoost classifier using
the FOX algorithm resulted in further performance enhancement:

  * •

Precision: 0.99 (weighted average)

  * •

Recall: 0.99 (weighted average)

  * •

F1-score: 0.99 (weighted average)

  * •

Accuracy: 0.99

Figure 4: Confusion matrix for FOX-optimized XGBoost model.

###  3.4 SHAP Analysis

To understand the contributions of individual radiomic features to the model
predictions, SHAP (SHapley Additive exPlanations) values were calculated for
the optimized models. SHAP values provide insights into which features are
most influential in the classification task.

Random Forest Model: The SHAP summary plot for the FOX-optimized Random Forest
model shows the top features contributing to the model’s predictions. The SHAP
values indicate that certain radiomic features have a substantial impact on
differentiating primary tumor sites from brain metastasis MRI data.

Figure 5: SHAP summary plot for FOX-optimized Random Forest model.

XGBoost Model: Similarly, the SHAP summary plot for the FOX-optimized XGBoost
model highlights the most influential features. The SHAP analysis confirms the
importance of selected radiomic features and provides a transparent
explanation of the model’s decision-making process.

Figure 6: SHAP summary plot for FOX-optimized XGBoost model.

###  3.5 Comparative Analysis

A comparative analysis of the four models (baseline and optimized versions of
Random Forest and XGBoost) is presented in Table 1. The results demonstrate
that the FOX-optimized models significantly outperform their baseline
counterparts. The XGBoost model, particularly when optimized with the FOX
algorithm, achieves the highest accuracy and F1-score among all models
evaluated.

Table 1: Performance metrics of baseline and FOX-optimized machine learning models Model | Precision | Recall | F1-score | Accuracy
---|---|---|---|---
Random Forest (Baseline) | 0.86 | 0.85 | 0.84 | 0.85
Random Forest (FOX Optimized) | 0.94 | 0.93 | 0.93 | 0.93
XGBoost (Baseline) | 0.96 | 0.96 | 0.96 | 0.96
XGBoost (FOX Optimized) | 0.99 | 0.99 | 0.99 | 0.99

##  4 Discussion

The results of this study demonstrate the significant potential of using
radiomic features extracted from brain metastasis (BM) MRI data in conjunction
with advanced machine learning algorithms to predict the primary tumor site.
The application of the FOX (Fox optimizer) algorithm for hyperparameter
optimization proved to be highly effective in enhancing the performance of
both the Random Forest and XGBoost models.

###  4.1 Model Performance

The baseline Random Forest model provided a reasonable starting point,
achieving an accuracy of 0.85. However, when optimized with the FOX algorithm,
the accuracy increased to 0.93, demonstrating the substantial impact of
hyperparameter tuning. Similarly, the XGBoost model showed impressive baseline
performance with an accuracy of 0.96, which further improved to 0.99 after FOX
optimization. These findings highlight the importance of using robust
optimization techniques to achieve the best possible model performance [13].

###  4.2 Feature Importance and Model Interpretability

The use of SHAP (SHapley Additive exPlanations) for model interpretability
provided valuable insights into the contributions of individual radiomic
features to the classification outcomes. SHAP summary plots for both the
Random Forest and XGBoost models revealed that certain features significantly
influenced the predictions. This interpretability is crucial for clinical
applications, as it helps in understanding which aspects of the MRI data are
most informative and potentially linked to specific primary tumors [17].

###  4.3 Clinical Relevance

Identifying the primary tumor site in patients with brain metastases is
critical for tailoring treatment strategies and improving patient outcomes.
The high accuracy achieved by the optimized XGBoost model suggests that
machine learning models, combined with radiomic analysis, can serve as
powerful tools in clinical decision-making. These models could potentially be
integrated into clinical workflows to assist radiologists and oncologists in
diagnosing and planning treatments more effectively.

###  4.4 Limitations

Despite the promising results, several limitations should be noted. The
dataset, although comprehensive, is limited to 75 patients, which may not
capture the full heterogeneity of brain metastases and their primary tumors.
Future studies should aim to include larger and more diverse datasets to
validate and generalize the findings. Additionally, while the FOX optimization
algorithm proved effective, exploring other optimization techniques could
further enhance model performance.

###  4.5 Future Work

Future research should focus on expanding the dataset and exploring additional
radiomic features and machine learning algorithms. Integrating genomic and
molecular data with radiomic features could provide a more comprehensive
understanding of tumor biology and improve prediction accuracy. Moreover,
developing real-time predictive tools that can be seamlessly integrated into
clinical practice will be essential for translating these findings into
practical applications.

##  5 Conclusion

This study demonstrates the feasibility and effectiveness of using radiomic
features from brain metastasis MRI data to predict the primary tumor site with
high accuracy. The combination of advanced machine learning models,
particularly the XGBoost classifier, and the FOX optimization algorithm
significantly improved predictive performance. The use of SHAP for model
interpretability provided insights into the most influential features,
enhancing the clinical relevance of the findings.

The results underscore the potential of radiomics and machine learning in
advancing cancer diagnosis and treatment. By continuing to refine these models
and validate them with larger datasets, we can move closer to integrating
these technologies into routine clinical practice, ultimately improving
patient outcomes in cases of brain metastasis.

##  6 Data and Code Availability

Data Availability: The dataset used in this study is publicly available and
can be accessed through the following link:

  * •

Ocaña-Tienda, B., Pérez-Beteta, J., Villanueva-García, J.D. et al. "A
comprehensive dataset of annotated brain metastasis MR images with clinical
and radiomic data." Scientific Data, 10, 208 (2023).
https://doi.org/10.1038/s41597-023-02123-0

This dataset includes detailed imaging and clinical data for 75 patients with
brain metastases, collected from multiple medical institutions, ensuring a
robust and comprehensive sample for radiomic and machine learning analyses.

Code Availability: The code used to preprocess the data, extract radiomic
features, train machine learning models, optimize hyperparameters, and
interpret model predictions using SHAP is available in a public GitHub
repository. Researchers and practitioners can access the complete codebase
through the following link:

  * •

GitHub Repository: BrainMetDetect: https://github.com/hamidreza-s-
salehi/BrainMetDetect

This repository includes all necessary scripts, documentation, and
dependencies to reproduce the study’s results and facilitate further research
and development in this domain. By providing access to both the data and code,
we aim to promote transparency, reproducibility, and collaboration within the
scientific community.

## References

  * [1] Xiang Zhang, Wei Zhang, Wei-Dong Cao, Gang Cheng, Bolin Liu, and Jinxiang Cheng.  A review of current management of brain metastases.  Annals of surgical oncology, 19:1043–1050, 2012.
  * [2] Krutika Deshpande, Ian Buchanan, Vahan Martirosian, and Josh Neman.  Clinical perspectives in brain metastasis.  Cold Spring Harbor Perspectives in Medicine, 10(6):a037051, 2020.
  * [3] Robert J Gillies, Paul E Kinahan, and Hedvig Hricak.  Radiomics: images are more than pictures, they are data.  Radiology, 278(2):563–577, 2016.
  * [4] Vishwa Parekh and Michael A Jacobs.  Radiomics: a new application from established techniques.  Expert review of precision medicine and drug development, 1(2):207–226, 2016.
  * [5] Michele Avanzo, Lise Wei, Joseph Stancanello, Martin Vallieres, Arvind Rao, Olivier Morin, Sarah A Mattonen, and Issam El Naqa.  Machine and deep learning methods for radiomics.  Medical physics, 47(5):e185–e202, 2020.
  * [6] Beatriz Ocaña-Tienda, Julián Pérez-Beteta, José D Villanueva-García, José A Romero-Rosales, David Molina-García, Yannick Suter, Beatriz Asenjo, David Albillo, Ana Ortiz de Mendivil, Luis A Pérez-Romasanta, et al.  A comprehensive dataset of annotated brain metastasis mr images with clinical and radiomic data.  Scientific data, 10(1):208, 2023.
  * [7] Chintan Parmar, Patrick Grossmann, Derek Rietveld, Michelle M Rietbergen, Philippe Lambin, and Hugo JWL Aerts.  Radiomic machine-learning classifiers for prognostic biomarkers of head and neck cancer.  Frontiers in oncology, 5:272, 2015.
  * [8] Luke Macyszyn, Hamed Akbari, Jared M Pisapia, Xiao Da, Mark Attiah, Vadim Pigrish, Yingtao Bi, Sharmistha Pal, Ramana V Davuluri, Laura Roccograndi, et al.  Imaging patterns predict patient survival and molecular subtype in glioblastoma via machine learning techniques.  Neuro-oncology, 18(3):417–425, 2015.
  * [9] Leo Breiman.  Random forests.  Machine learning, 45:5–32, 2001.
  * [10] Mahesh Pal.  Random forest classifier for remote sensing classification.  International journal of remote sensing, 26(1):217–222, 2005.
  * [11] Jiawei Han, Jian Pei, and Hanghang Tong.  Data mining: concepts and techniques.  Morgan kaufmann, 2022.
  * [12] Ron Kohavi et al.  A study of cross-validation and bootstrap for accuracy estimation and model selection.  In Ijcai, volume 14, pages 1137–1145. Montreal, Canada, 1995.
  * [13] Seyedali Mirjalili, Seyed Mohammad Mirjalili, and Andrew Lewis.  Grey wolf optimizer.  Advances in engineering software, 69:46–61, 2014.
  * [14] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pages 785–794, 2016.
  * [15] Hardi Mohammed and Tarik Rashid.  Fox: a fox-inspired optimization algorithm.  Applied Intelligence, 53(1):1030–1050, 2023.
  * [16] Gai-Ge Wang, Suash Deb, and Zhihua Cui.  Monarch butterfly optimization.  Neural computing and applications, 31:1995–2014, 2019.
  * [17] Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  Advances in neural information processing systems, 30, 2017.