# Singlet Dirac dark matter streamlined

Carlos E. Yaguna  Escuela de Física, Universidad Pedagógica y Tecnológica de
Colombia,
Avenida Central del Norte # 39-115, Tunja, Colombia  Óscar Zapata  Instituto
de Física, Universidad de Antioquia,
Calle 70 # 52-21, Apartado Aéreo 1226, Medellín, Colombia

###### Abstract

We propose a new and compact realization of singlet Dirac dark matter within
the WIMP framework. Our model replaces the standard $Z_{2}$ stabilizing
symmetry with a $Z_{6}$, and uses spontaneous symmetry breaking to generate
the dark matter mass, resulting in a much simplified scenario for Dirac dark
matter. Concretely, we extend the Standard Model (SM) with just two new
particles, a Dirac fermion (the dark matter) and a real scalar, both charged
under the $Z_{6}$ symmetry. After acquiring a vacuum expectation value, the
scalar gives mass to the dark matter and mixes with the Higgs boson, providing
the link between the dark sector and the SM particles. With only four free
parameters, this new model is extremely simple and predictive. We study the
dark matter density as a function of the model’s free parameters and use a
likelihood approach to determine its viable parameter space. Our results
demonstrate that the dark matter mass can be as large as $6$ TeV while
remaining consistent with all known theoretical and experimental bounds. In
addition, a large fraction of viable models turns out to lie within the
sensitivity of future direct detection experiments, furnishing a promising way
to test this appealing scenario.

##  I Introduction

Weakly interacting massive particles (WIMPs) still are one of the best
motivated candidates to explain the dark matter observed throughout the
Universe Roszkowski et al. (2018); Arcadi et al. (2018). They can be easily
incorporated into many Standard Model (SM) extensions, naturally give rise to
a relic abundance of the right order of magnitude Aghanim et al. (2020), and
can usually be probed in current and future dark matter Schumann (2019) or
collider experiments Buchmueller et al. (2017). Within this WIMP framework,
many different dark matter models can be envisaged, depending on the number
and type of new fields that are included and on the additional symmetries that
may be imposed.

In this paper we propose a new and simplified scenario for Dirac dark matter.
In it, the SM is extended with just two additional fields –a Dirac fermion
(the dark matter) and a real scalar, both singlets under the SM gauge group.
The main novelty of our model is the use of a $Z_{6}$ discrete symmetry rather
than the standard $Z_{2}$ as a stabilizing symmetry for the dark matter
particle. Under it, only the new fields transform non-trivially. The resulting
scenario is very predictive as it admits only four free parameters,
significantly fewer than similar models for Dirac dark matter. After
introducing the model and describing its dark matter phenomenology, we use a
likelihood approach to identify its viable parameter space –the regions
consistent with all known theoretical and experimental bounds. Our results
demonstrate that this scenario is viable over a wide range of dark matter
masses extending from about $60$ GeV to $6$ TeV. Moreover, the new scalar
turns out to play an essential role: it generates the dark matter mass via its
vacuum expectation value; it allows the dark matter to annihilate efficiently
in the early Universe so as to satisfy the relic density constraint; and it
links the dark matter with the SM fields, providing a way to probe this
scenario in current and planned experiments. This simple setup, therefore,
provides a minimal, consistent and testable realization of singlet Dirac dark
matter.

Numerous studies have investigated singlet fermion dark matter models,
primarily employing a $Z_{2}$ symmetry and introducing additional free
parameters (e.g., Kim et al. (2008); Baek et al. (2012a, b); Esch et al.
(2013); Freitas et al. (2015)). These scenarios have also been analyzed from
an effective field theory standpoint (e.g., Kim and Lee (2007); Fedderke et
al. (2014); Matsumoto et al. (2014); Athron et al. (2019)), and the prospect
of electroweak baryogenesis has gained recent attention (e.g., Fairbairn and
Hogan (2013); Li and Zhou (2014); Ghorbani (2017); Beniwal et al. (2019);
Gould and Xie (2023)). Over the past few years, $Z_{N}$ symmetries have been
integrated into diverse dark matter models, often involving multi-component
dark matter (e.g., Yaguna and Zapata (2020); Bélanger et al. (2020); Belanger
et al. (2022); Yaguna and Zapata (2021, 2022); Choi et al. (2021); Díaz Sáez
and Ghorbani (2023); Jurčiukonis and Lavoura (2023); Bélanger et al. (2023);
Dupuis (2016); Bhattacharya et al. (2018); Ho et al. (2022); Yaguna and Zapata
(2024)). WIMP models using a $Z_{4}$ with a singlet scalar mediator featuring
a Dirac and a Majorana dark matter candidate were analyzed in Ref. Dupuis
(2016) and Ref. Bhattacharya et al. (2018), respectively (in Refs. Ho et al.
(2022) and Yaguna and Zapata (2024) alternative dark matter production
mechanisms were implemented). To our knowledge, the $Z_{6}$ model for Dirac
dark matter we focus on is novel.

The rest of the paper is organized as follows. In the next section the model
is introduced and its free parameters are identified. The dark matter
phenomenology is discussed in section III, where the variation of the relic
density with the model’s parameters is analyzed. Section IV presents our main
results. There, the viable regions of this model are determined and analyzed,
and the detections prospects are investigated. In section V we contrast this
model against related scenarios considered previously in the literature, and
some of its possible extensions are mentioned. Finally, our conclusions are
presented in section VI .

##  II The Model

To account for the dark matter, we propose to extend the SM particle content
with just two chiral fermions ($\psi_{L}$ and $\psi_{R}$) and a real scalar
field ($\phi$), all singlets under the SM gauge symmetry. In addition, a
$Z_{6}$ discrete symmetry is assumed under which only the new fields transform
non-trivially and as follows111Notice that both fermions must be charged under
the $Z_{6}$ to avoid a Yukawa interaction term with the lepton doublet, e.g.,
$\bar{L}\tilde{H}\psi_{R}$.

$\displaystyle\phi\rightarrow\omega_{6}^{3}\,\phi,\,\,\,\psi_{L}\rightarrow\omega_{6}\,\psi_{L},\,\,\,\psi_{R}\rightarrow\omega_{6}^{4}\,\psi_{R};\hskip
28.45274pt\omega_{6}=e^{\frac{2\pi i}{6}}.$ (1)

This particle content and symmetry transformations specify our setup.

A consequence of the above transformations is that bare mass terms for the
fermions are forbidden by the $Z_{6}$ symmetry, either of Dirac or of Majorana
type. It is possible, however, to write the following interaction term between
$\psi_{L}$, $\psi_{R}$ and $\phi$:

$\displaystyle\mathcal{L}_{\psi}$
$\displaystyle=\,\,-y_{s}\overline{\psi}_{L}\psi_{R}\,\phi+\rm{h.c.},$ (2)

which will induce a Dirac mass term for the fermions once $\phi$ acquires a
vacuum expectation value.

The tree-level scalar potential for this model is

$\displaystyle\mathcal{V}_{0}$
$\displaystyle=\,\,-\mu_{H}^{2}|H|^{2}+\lambda_{H}|H|^{4}-\frac{1}{2}\mu_{\phi}^{2}\phi^{2}+\frac{1}{4}\lambda_{\phi}\phi^{4}+\frac{1}{2}\lambda_{SH}|H|^{2}\phi^{2},$
(3)

where $H$ is the SM Higgs doublet. Notice that terms such as $\phi^{3}$ or
$\phi|H|^{2}$ are forbidden by the $Z_{6}$ symmetry because $\phi$ transforms
to $-\phi$ under it. As explained above, $\phi$ must acquire a vacuum
expectacion value so we write

$\displaystyle H$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}0\\\
H^{0}_{R}+v_{H}\end{array}\right),\hskip
28.45274pt\phi=\left(\phi_{R}+v_{\phi}\right).$ (6)

After the $Z_{6}$ and electroweak symmetry breaking the real components
$\phi_{R}$ and $H^{0}_{R}$ get mixed through the $\lambda_{SH}$ term.
Specifically, the $2\times 2$ mass matrix in the $(H^{0}_{R},\phi_{R})$ basis
has the entries

$\displaystyle M^{2}_{11}$ $\displaystyle=2\lambda_{H}v_{H}^{2},$ (7)
$\displaystyle M^{2}_{22}$ $\displaystyle=2\lambda_{\phi}v_{\phi}^{2},$ (8)
$\displaystyle M^{2}_{12}$
$\displaystyle=M^{2}_{21}=\lambda_{SH}v_{\phi}v_{H}.$ (9)

The mass eigenstates $h,S$ are defined through the rotation matrix

$\displaystyle\begin{pmatrix}H^{0}_{R}\\\ \phi_{R}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}h\\\ S\end{pmatrix},$ (10)

with a mixing angle

$\displaystyle\tan 2\theta$
$\displaystyle=\frac{\lambda_{SH}v_{\phi}v_{H}}{\lambda_{\phi}v_{\phi}^{2}-\lambda_{H}v_{H}^{2}},$
(11) $\displaystyle\sin 2\theta$
$\displaystyle=\frac{2\lambda_{SH}v_{\phi}v_{H}}{M_{S}^{2}-M_{h}^{2}},$ (12)
$\displaystyle\cos 2\theta$
$\displaystyle=\frac{2\lambda_{\phi}v_{\phi}^{2}-2\lambda_{H}v_{H}^{2}}{M_{h}^{2}-M_{S}^{2}},$
(13)

such that ${\rm diag}(M_{h}^{2},M_{S}^{2})=R(\theta)^{T}M^{2}R(\theta)$. The
scalar masses are

$\displaystyle M_{h,S}^{2}$
$\displaystyle=\frac{1}{2}\left[2\lambda_{H}v_{H}^{2}+2\lambda_{\phi}v_{\phi}^{2}\mp\sqrt{\left[2\lambda_{H}v_{H}^{2}-2\lambda_{\phi}v_{\phi}^{2}\right]^{2}+4\lambda_{SH}^{2}v_{\phi}^{2}v_{H}^{2}}\right],$
(14) $\displaystyle M_{h}^{2}$
$\displaystyle=2(v_{H}^{2}\lambda_{H}\cos^{2}\theta-\lambda_{SH}v_{\phi}v_{H}\cos\theta\sin\theta+\lambda_{\phi}v_{\phi}^{2}\sin^{2}\theta),$
(15) $\displaystyle M_{S}^{2}$
$\displaystyle=2(v_{H}^{2}\lambda_{H}\sin^{2}\theta+\lambda_{SH}v_{\phi}v_{H}\cos\theta\sin\theta+\lambda_{\phi}v_{\phi}^{2}\cos^{2}\theta)$
(16)

where we assumed, without loss of generality, that $M_{h}<M_{S}$. $h$ can be
identified with the Higgs boson observed at the LHC with a mass of about $125$
GeV Aad et al. (2015).

From the minima conditions we can obtain $\mu_{H}$ and $\mu_{\phi}$, which are
replaced back in $\mathcal{V}$:

$\displaystyle\frac{\partial\mathcal{V}}{\partial\phi_{R}}=0$
$\displaystyle\to
v_{\phi}\left(-\mu_{\phi}^{2}+\lambda_{\phi}v_{\phi}^{2}+\frac{1}{2}\lambda_{SH}v_{H}^{2}\right)=0,$
(17) $\displaystyle\frac{\partial\mathcal{V}}{\partial H^{0}_{R}}=0$
$\displaystyle\to
v_{H}\left(-\mu_{H}+\lambda_{H}^{2}v_{H}^{2}+\frac{1}{2}\lambda_{SH}v_{\phi}^{2}\right)=0,$
(18)

leading to

$\displaystyle\lambda_{H}$
$\displaystyle=(M_{h}^{2}\cos^{2}\theta+M_{S}^{2}\sin^{2}\theta)/(2v_{H}^{2}),$
(19) $\displaystyle v_{\phi}$
$\displaystyle=(M_{S}^{2}-M_{h}^{2})\sin\theta\cos\theta/(\lambda_{SH}v_{H}),$
(20) $\displaystyle\lambda_{\phi}$
$\displaystyle=\left(M_{h}^{2}\sin^{2}\theta+M_{S}^{2}\cos^{2}\theta\right)/(2v_{\phi}^{2}).$
(21)

In this way, the parameters of the scalar potential can be written in terms of
the physical parameters.

We can now write the fermion interaction term in the mass eigenstate basis,

$\displaystyle\mathcal{L}_{\psi}$
$\displaystyle=\,\,-y_{s}\overline{\psi}_{L}\psi_{R}\,\phi+\rm{h.c.}$ (22)
$\displaystyle=\,\,-m_{\psi}\overline{\psi}\psi\left[1+(-h\,\sin\theta+S\,\cos\theta)/v_{\phi}\right],$
(23)

where

$\displaystyle y_{s}=m_{\psi}/v_{\phi}.$ (24)

Thus, the chiral fields $\psi_{L}$ and $\psi_{R}$ form a massive Dirac
fermion, $\psi=\psi_{L}+\psi_{R}$. Notice that $\psi$ is automatically stable
and thus the dark matter particle in this model. It couples only to $h$ and
$S$ and with an strength proportional to its mass. In addition, since $\theta$
must be small (see section IV), $\psi$ couples more strongly to $S$ than to
$h$.

Altogether, this model includes only 4 free parameters, which can be taken to
be

$\displaystyle m_{\psi},M_{S},\sin\theta,\lambda_{SH}.$ (25)

All other parameters are either known or can be determined from these using
the equations above. Let us emphasize that no extra assumptions have been made
to simplify the Lagrangian or the parameter space –this is the most general
setup consistent with the particle content and the symmetries imposed. Given
this restricted set of parameters, one may wonder whether it is actually
possible to explain the observed density of dark matter while remaining
consistent with known experimental and theoretical bounds. That is the
question we will address in the next two sections.

##  III Dark matter phenomenology

In this model, the dark matter particle couples directly only to $h$ and $S$,
and through them, via mixing, to the other SM particles. For illustration,
figure 1 displays the Feynman diagrams that contribute to dark matter
annihilation into two scalars. Throughout this work, the dark matter
observables are calculated by micrOMEGAs version 5.0 Bélanger et al. (2018)
via LanHEP Semenov (2016). Let us stress that micrOMEGAs numerically solves
the Boltzmann equation for $\psi$, automatically taking into account all the
different processes that contribute to the dark matter relic density,
including decays and $2\to 2$ scatterings.

Figure 1: Some of the Feynman diagrams that contribute to dark matter
annihilation in this model. In the left diagram, the final state scalars can
be replaced by $W,Z$ or SM fermions.

Figure 2 displays, as a function of the dark matter mass, the relative
contributions of different final states to the relic density. As anticipated,
the dominant annihilation final state strongly depends on $M_{\psi}$, being
$b\bar{b}$ for dark matter masses smaller than $M_{h}/2$ –below the Higgs
resonance. Right above the Higgs resonance, other final states become relevant
including $W^{+}W^{-}$, $HH$ and $HS$. Finally, for $M_{\psi}\gtrsim M_{S}$ it
is the $SS$ final state that entirely dominates the relic density. Even if
this figure was obtained for a specific set of parameters –$M_{S}=400$ GeV,
$\lambda_{SH}=\sin\theta=0.1$– the pattern it reveals is quite generic.
Notice, from the right diagram in figure 1, that the cross section for the
annihilation into $SS$ goes like $y_{s}^{4}$ and since $M_{\psi}\propto
y_{s}$, that process is expected to dominate at large dark matter masses. In
fact, as we will see in the next section, for _all_ the viable models that lie
outside the resonances, the process that sets the dark matter relic density is
indeed $\psi\bar{\psi}\to SS$.

Figure 2: The relative contributions of different final states to the total
relic density.

Let us now study how the relic density varies with the free parameters of the
model. Figure 3 shows $\Omega h^{2}$ versus the dark matter mass for different
values of the Higgs-portal coupling, $\lambda_{SH}$. The other two parameters
were taken as $\sin\theta=0.1$ and $M_{S}=300$ GeV. Notice that $\lambda_{SH}$
has a direct impact on the dark matter interactions because, according to
equations (20) and (24), $y_{s}\propto\lambda_{SH}$. Thus, we expect
$\Omega\propto\lambda_{SH}^{-2}$ at low dark matter masses whereas
$\Omega\propto\lambda_{SH}^{-4}$ at high dark matter masses, in agreement with
what is observed in the figure. Two additional features are clearly observed
in this plot. First, the suppression of the relic density at the Higgs and $S$
resonances –$M_{\psi}=62.5,150$ GeV respectively. Second, the drop in the
relic density once the $SS$ annihilation channel becomes kinematically
accessible ($M_{\psi}\sim 300$ GeV). Finally, notice that agreement with the
observed dark matter density Aghanim et al. (2020) (horizontal cyan band) can
be obtained, for this set of parameters, only if $\lambda_{SH}$ is not too
small and if the dark matter mass lies either at the resonances or above
$M_{S}$.

Figure 3: The relic density as a function of the dark matter mass for
different values of $\lambda_{SH}$. From top to bottom the lines correspond to
$\lambda_{SH}=0.003$ (solid, red), $0.01$ (dashed, green), $0.03$ (dotted,
blue), and $0.1$ (dash-dotted, orange). For this figure $M_{S}=300$ GeV and
$\sin\theta=0.1$. The (cyan) horizontal band corresponds to the observed dark
matter density.

Figure 4 shows instead the dark matter density for different values of
$M_{S}$. For this plot the dimensionless parameters $\lambda_{SH}$ and
$\sin\theta$ were both set to $0.1$. In this case, the position of the $S$
resonance moves from line to line, and the relic density tends to increases
with $M_{S}$. Notice that for $M_{S}=1.6$ TeV the predicted relic density is
always larger than the observed value (cyan horizontal band), whereas for the
other cases it is not so and it becomes possible to satisfy the relic density
constraint for some values of $M_{DM}$. Finally, the variation of the relic
density with $\sin\theta$ is illustrated in figure 5 for $M_{S}=300$ GeV and
$\lambda_{SH}=0.01$. $\sin\theta$ affects the dark matter interactions mainly
through $v_{\phi}$ –see equation (20)–, which in turn determines $y_{s}$. At
low dark matter masses, the dependence on $\sin\theta$ cancels out among the
two vertices but at high masses we expect $\Omega h^{2}\propto\sin\theta^{4}$,
in agreement with the plot. Note that the observed value of the dark matter
density tends to prefer smaller values of $\sin\theta$.

Figure 4: The relic density as a function of the dark matter mass for
different values of $M_{S}$. From bottom to top the lines correspond to
$M_{S}=200$ GeV (solid, red), $400$ GeV (dashed, green), $800$ GeV (dotted,
blue), and $1.6$ TeV (dash-dotted, orange). For this figure
$\lambda_{SH}=\sin\theta=0.1$. The (cyan) horizontal band corresponds to the
observed dark matter density. Figure 5: The relic density as a function of the
dark matter mass for different values of $\sin\theta$. From top to bottom the
lines correspond to $\sin\theta=0.3$ (solid, red), $0.1$ (dashed, green),
$0.03$ (dotted, blue), and $0.01$ (dash-dotted, orange). For this figure
$M_{S}=300$ GeV and $\lambda_{SH}=0.01$. The (cyan) horizontal band
corresponds to the observed dark matter density.

Besides the relic density, another crucial observable in this model is the
dark matter scattering cross section, which is subject to strong limits from
direct detection experiments such as XENON-nT Aprile et al. (2022) and LZ
Aalbers et al. (2023). Dark matter-quark scattering proceeds via a t-channel
diagram mediated by $h$ and $S$ –see Figure 6. The resulting dark matter-
nucleon cross section is spin-independent and given by

$\sigma_{SI}=\frac{\mu^{2}}{\pi}\left(\frac{y\sin\theta\cos\theta}{v_{H}}\right)^{2}\left(\frac{1}{M_{h}^{2}}-\frac{1}{m_{S}^{2}}\right)^{2}m_{\mathcal{N}}^{2}f_{\mathcal{N}}^{2}$
(26)

with $m_{\mathcal{N}}=939$ MeV, $f_{\mathcal{N}}\approx 0.3$ and
$\mu=m_{\psi}m_{\mathcal{N}}/(m_{\psi}+m_{\mathcal{N}})$. From this
information alone, however, it is not possible to say much about the
importance of current direct detection bounds. To that end, we must first find
the models that are consistent with the observed value of the dark matter
density. That is precisely what we are going to do in the next section.

Figure 6: The Feynman diagrams that contribute to direct detection in this
model.

##  IV The viable parameter space

In this section we present the main results of the exploration of the
parameter space: the regions compatible with all the constraints are defined
and characterized, and the direct and indirect detection prospects are also
analyzed.

###  IV.1 Sampling

Free parameter | Range | Prior
---|---|---
$M_{\psi}$ |  $[10,10^{4}]$ GeV | Log
$M_{S}$ |  $[150,10^{4}]$ GeV | Log
$\lambda_{SH}$ | $[10^{-4},10]$ | Log
$\sin\theta$ | $[10^{-5},1]$ | Log
Table 1: Intervals for the model’s free parameters used with the differential
evolution scan, alongside their corresponding prior likelihoods.

Our numerical analysis aims to make a statistically sound inference on the
viable parameter space AbdusSalam et al. (2022). For this purpose, we adopt a
frequentist approach and perform 4-dimensional scans using the differential
evolution algorithm Storn and Price (1997); Virtanen et al. (2019). The
relevant outputs for a frequentist analysis, such as best-fit points and 68.3%
(1$\sigma$) and 95.4% (2$\sigma$) confidence level (C.L.) regions, are created
by maximizing the likelihood via Wilks’ theorem Wilks (1938). We assume the
following likelihood function

$\displaystyle\mathcal{L}(\vec{\alpha})=\mathcal{L}_{\Omega}(\vec{\alpha})\times\mathcal{L}_{\text{DD}}(\vec{\alpha}),$
(27)

that accounts for the relic density measurement and the direct detection upper
bound as a function of the free parameters
$\vec{\alpha}=(M_{\psi},M_{S},\lambda_{SH},\sin\theta)$ (see Table 1).
Moreover, we rely on the profile likelihood ratio (PLR), which is defined as

$\displaystyle\Lambda=\frac{\mathcal{L}(\alpha_{i})}{\mathcal{L}_{\rm max}},$
(28)

where

$\displaystyle\mathcal{L}(\alpha_{i})=\underset{\alpha_{i}\notin\vec{\alpha}}{{\rm
max}}\,[\mathcal{L(\vec{\alpha})}],\hskip 28.45274pt\mathcal{L}_{\rm max}={\rm
max}\,[\mathcal{L(\vec{\alpha})}],$ (29)

We assume Gaussian likelihoods in such a way that their corresponding
expressions read

$\displaystyle\chi_{\mathcal{O}}^{2}(\vec{\alpha})=-2\ln\mathcal{L}_{\mathcal{O}}(\vec{\alpha})=\left(\frac{\mathcal{O}^{\text{th}}-\mathcal{O}^{\text{obs}}}{\Sigma_{\mathcal{O}}}\right)^{2},$
(30)

where $\mathcal{O}^{\text{th}}$ and $\mathcal{O}^{\text{obs}}$ are the
theoretical prediction and experimental measurement of the observable
$\mathcal{O}$, respectively, and $\Sigma_{\mathcal{O}}$ is the associated
total standard deviation. From the dark matter abundance measurement reported
by PLANCK Aghanim et al. (2020),

$\displaystyle\Omega_{\text{DM}}h^{2}=0.1198\pm 0.0012.$ (31)

we associate the the central value to $\Omega^{\text{obs}}=0.1198$. Moreover,
we consider a model to be compatible with this measurement if the relic
density obtained from micrOMEGAs, $\Omega^{\text{th}}$, lies roughly between
0.11 and 0.13, that is, within an estimated theoretical and experimental
uncertainty of 10%, $\Sigma_{\Omega}=0.1\Omega^{\text{obs}}$. Instead, to
account for the current upper bounds reported by the LZ collaboration Aalbers
et al. (2023) a null mean value is adopted,
$\sigma_{\text{SI}}^{\text{obs}}=0$, whereas $\sigma_{\text{SI}}^{\text{th}}$
is calculated through micrOMEGAS. Finally,
$\Sigma^{2}_{\sigma_{\text{SI}}}=[0.2\sigma_{\text{SI}}^{\text{th}}]^{2}+[\Sigma_{\sigma_{\text{SI}}}^{\text{LZ}}/1.64]^{2}$,
where $0.2\sigma_{\text{SI}}^{\text{th}}$ corresponds to an estimated
theoretical uncertainty of 20% and $\Sigma_{\sigma_{\text{SI}}}^{\text{LZ}}$
is the 90% confidence level upper bound (within $2\sigma$ sensitivity band) on
the spin-independent cross section reported by LZ collaboration.

We impose the following theoretical constraints on the scalar potential
parameters:

  * •

Perturbatibity:

$\displaystyle|\lambda_{H}|<\frac{2}{3}\pi,\,\,\,|\lambda_{\phi}|<\pi,\,\,\,|\lambda_{SH}|<4\pi,$
(32)

  * •

Perturbative unitarity:

$\displaystyle|\lambda_{H}|<8\pi,\,\,\,|\lambda_{\phi}|<8\pi,\,\,\,\left|6\lambda_{H}+\lambda_{\phi}\pm\sqrt{(6\lambda_{H}-\lambda_{\phi})^{2}+8\lambda^{2}_{SH}}\right|<16\pi.$
(33)

  * •

Vacuum stability: For this case the necessary and sufficient vacuum stability
conditions are

$\displaystyle\lambda_{H}>0,\,\,\,\lambda_{\phi}>0,\,\,\,\lambda_{SH}+2\sqrt{\lambda_{H}\lambda_{\phi}}>0.$
(34)

  * •

Global vacuum $(v_{H},v_{S})$:

$\displaystyle\lambda_{SH}v_{\phi}^{2}v_{H}^{2}+\lambda_{\phi}v_{\phi}^{4}>0,\,\,\,\,\,\lambda_{H}v_{H}^{4}+\lambda_{SH}v_{\phi}^{2}v_{H}^{2}>0.$
(35)

Besides, the Yukawa coupling is subject to the perturbativity constraint

$\displaystyle y_{s}<\sqrt{4\pi}.$ (36)

Concerning collider constraints, since the Higgs signal strength gets modified
by a factor of $\cos^{2}\theta$, it follows that $|\sin\theta|\lesssim 0.3$
for $M_{S}>M_{h}$ Falkowski et al. (2015); Arcadi et al. (2020); Ferber et al.
(2023) from the most recent LHC results ATL (2022); Tumasyan et al. (2022).
The viable parameter space of the model is determined after imposing all the
constraints and selecting those samples fulfilling
$-2\ln\mathcal{L}(\vec{\alpha})\leq 5.99$ ($2\sigma$ of C.L. for a single
degree of freedom).

###  IV.2 Results

Figure 7: Viable parameter space in terms of the free parameters of the model.

The viable regions of this model are shown in figure 7, projected onto the
plane $(M_{\psi},M_{S})$ in the left panel, and $(M_{\psi},\lambda_{SH})$ in
the right panel, with the value of $\sin\theta$ color coded. Remarkably, these
viable points span a range of two orders of magnitude in the dark matter mass,
from $60$ GeV to $6$ TeV approximately. In the left panel, three regions can
be clearly differentiated, depending on how the relic density constraint is
satisfied: the vertical line at $M_{\psi}\sim M_{h}/2$ is due to the Higgs
resonance; the diagonal line extending up to dark matter masses of $5$ TeV
corresponds to the $S$ resonance; and the remaining triangular region
comprises the models where the relic density is obtained via the annihilations
$\psi\bar{\psi}\to SS$. Notice that this region prefers small values of
$\sin\theta$, as explained in the previous section, whereas the two resonances
favor larger values. From the right panel, we also see that, as expected,
$\lambda_{SH}$ tends to be larger at the resonances. In fact, outside the
resonances we find $\lambda_{HS}\lesssim 0.5$, a fact that, as shown in the
next subsection, has profound implications for the detection prospects in this
model.

It is also useful to know the allowed values for the parameters
$(\lambda_{H},\lambda_{\phi},v_{\phi},y_{s})$, even if they are not free –Eqs.
(19)-(24). In figure 8 we display these parameters and at the same time
classify the viable points according to the mechanism responsible for their
relic density: green for the Higgs resonance, blue for the $S$ resonance, and
red for outside the resonances (equivalently $\psi\bar{\psi}\to SS$). Notice
from the top-left panel that for the resonant regions $\lambda_{H}$ may exceed
significantly its SM value ($0.127$), providing an indirect way to probe this
scenario. From the bottom-right panel, instead, we can see that $y_{s}$ must
be large when the relic density is determined via $\psi\bar{\psi}\to SS$, as
it is the coupling that fixes that annihilation cross section.

Figure 8: Non-free parameters $(\lambda_{H},\lambda_{\phi},v_{\phi},y_{s})$ as
a function of the dark matter mass according to their belongingness to the
region I (red plus symbols), II (blue dots) or III (green crosses).

Figure 9 shows the relative contributions of different final states to the
total relic density, for points in the $S$ resonance (colored plus symbols)
and outside the resonances (black dots). Within the resonance we see that the
dark matter annihilates mostly into gauge ($\sim 75\%$) and Higgs ($\sim
25\%$) bosons. In contrast, the process that sets the dark matter relic
density for all viable models lying outside the resonances is
$\psi+\bar{\psi}\to S+S$ (black points), as anticipated in the previous
section –see Fig. 2.

Figure 9: Relative contributions of the dark matter annihilation channels to
the relic abundance. The colored plus symbols correspond to viable models
inside $S$-resonance (region II) while the black dots to those outside this
(region I).

The results for the one-dimensional PLR from the 4-dimensional scans are
displayed in Fig. 10. The other two 1D PLR for $M_{S}$ and $\sin\theta$ are
not shown since they are not constrained by the fit (the PLR is basically flat
and equal to 1 within the whole range for these parameters). From these
results it follows that the global fit slightly prefers dark matter masses
above the Higgs mass and Higgs portal couplings below $0.5$. The $1\sigma$ and
$2\sigma$ C.L. regions in the 2-dimensional planes involving free parameters
are mostly unconstrained by our global fit except in the white zones in Fig.
7.

Figure 10: One-dimensional PLR for $M_{\psi}$ (left panel) and $\lambda_{SH}$
(right panel).

###  IV.3 Dark matter detection prospects

Now that the viable parameter space of this model has been determined, it is
time to look at its implications, in particular regarding dark matter
detection searches. The most important dark matter annihilation channels are
displayed in figure 11, following the same convention as in figure 9. Being
$p$-wave suppressed, the annihilation into $SS$ features a cross section about
five orders of magnitude below the thermal value, at around $3\times
10^{-31}$cm3/s. For models at the $S$ resonance, the cross section can be
larger but they all lie well below the thermal value and below the sensitivity
of current experiments. To illustrate this point, the limits reported by the
Fermi collaboration from $\gamma$-ray observations of dShps Ackermann et al.
(2015) and the limits Reinert and Winkler (2018) from antiproton searches
Aguilar et al. (2016) for $WW$ final states are shown as dashed and dotted
lines respectively. Indirect detection searches, therefore, do not currently
constraint the viable parameter space of this model and are not expected to do
so in the near future.

Figure 11: More important dark matter annihilation rates as a function of
$M_{\psi}$. The colored plus symbols correspond to viable models inside
$S$-resonance while the black dots to those outside this. The dashed and
dotted lines are the limits reported by the Fermi collaboration from
$\gamma$-ray observations of dShps Ackermann et al. (2015) and the limits
Reinert and Winkler (2018) from antiproton searches Aguilar et al. (2016) for
$WW$ final states, respectively.

More interesting is the case for direct detection. Figure 12 shows the viable
models projected onto the plane ($M_{\psi}$, $\sigma^{SI}$), with the color
indicating the mechanism responsible for the relic density, as before. By
construction, all the viable models are consistent with the most recent
results from the LZ experiment, which sets the upper bound on the spin-
independent cross section. The dotted line corresponds to the expected
sensitivity of the DARWIN experiment Aalbers et al. (2016). Remarkably, this
model can give rise to direct detection signals over the whole range of viable
dark matter masses. In fact, the entire Higgs-resonance region could be
excluded if no events are observed at DARWIN. And most of the $S$-resonance
region will be probed by that experiment. The fact that the resonances are
easier to test in direct detection experiments is related to the value of
$\lambda_{HS}$, which is higher in such regions, as shown previously.
Regarding the out of resonances region, we see from the figure that, even if
not guaranteed, a significant fraction of models may also induce an observable
signal in DARWIN. Hence, dark matter direct detection in this model looks very
encouraging and it provides the most promising way to test it in the near
future.

Figure 12: Spin-independent cross-sections for elastic scattering of $\psi$
with nuclei as a fucntion of $M_{\psi}$ in regions I, II and III. The dashed
line shows the projected sensitivity of DARWIN experiment Aalbers et al.
(2016).

##  V Discussion

As we have seen, the $Z_{6}$ model we proposed is a simple, viable and
testable scenario for Dirac dark matter. Let us briefly contrast this scenario
against the $Z_{2}$-invariant singlet fermion model often considered in the
literature –e.g., Kim et al. (2008); Baek et al. (2012a, b); Freitas et al.
(2015); Lopez-Honorez et al. (2012); Esch et al. (2013). In it, the fermion is
odd while the scalar is even under the $Z_{2}$ so that the Lagrangian includes
the following extra terms

$\displaystyle\mathcal{L}_{Z_{2}}\supset-
m_{D}\bar{\psi}\psi-\left[\mu^{3}_{1}\phi+\frac{1}{2}\mu_{2}\phi|H|^{2}+\frac{1}{3}\mu_{3}\phi^{3}+\rm{h.c.}\right].$
(37)

Even if the $\mu_{1}$ term can be set to zero under a shift of the singlet
field, the extra scalar parameters $\mu_{2}$ and $\mu_{3}$ induce additional
scalar interactions that give rise to fermion self-annihilations and elastic
scatterings with nuclei. Moreover, a bare mass term is now allowed for the
dark matter ($m_{D}$ above), which means that $m_{\psi}$ and $y_{s}$ become
independent of each other (see Eq. (36)). Thus, the $Z_{6}$ model studied here
is simpler and more predictive.

One may wonder why do we impose a $Z_{6}$ and not some other $Z_{N}$ symmetry.
It turns out that, for our particle content, the $Z_{6}$ is the lowest $Z_{N}$
symmetry that can generate Dirac dark matter via spontaneous symmetry breaking
–a $Z_{2}$ or a $Z_{4}$ would instead induce Majorana masses. Having three
non-equivalent and non-trivial charges, each for the chiral fields
$\psi_{L,R}$ and the scalar field $\phi$, the $Z_{6}$ forbids bare mass terms
for the fermions, either of Dirac or of Majorana type, while allowing a Yukawa
interaction among all three fields. This interaction not only induces a Dirac
mass after symmetry breaking but it also links the dark matter to the known
particles, as shown before. The $Z_{6}$ is thus quite special in this context.

An important feature of the $Z_{N}$ symmetries is that they may naturally
appear as remnants of spontaneously broken $U(1)$ gauge symmetries. It is
natural, therefore, to ask whether the $Z_{6}$ model investigated here can be
incorporated into a $U(1)$ gauge extension of the SM. Besides a complex scalar
(rather than real), such a scenario would demand an extended fermion sector to
ensure the cancellation of $U(1)$ anomalies. Due to the new gauge interactions
and the extra fermion fields, which may constitute a second dark matter sector
Bernal et al. (2019), the phenomenology of this gauge model is expected to
differ from the $Z_{6}$ scenario we discussed. Another interesting aspect we
have not addressed is the fact that the reduction in the number of free
parameters (with respect to the $Z_{2}$-variant) may substantially modify
those regions of the parameter space suitable for a strong first-order
electroweak phase transition Fairbairn and Hogan (2013); Li and Zhou (2014);
Ghorbani (2017); Beniwal et al. (2019); Gould and Xie (2023). Hence, the
possibility of realizing electroweak baryogenesis Ghorbani (2017); Beniwal et
al. (2019) as well as the signatures of this model at current and future
gravitational wave experiments Beniwal et al. (2019) need to be studied in a
future work.

##  VI Conclusions

One of the most pressing problems in fundamental physics today is determining
the correct Standard Model extension that accounts for the dark matter. In
this work we analyzed a new concise model in which the dark matter particle is
a Dirac fermion, singlet under the SM gauge group. This model includes only
two new particles –the dark matter and a real scalar field– and four free
parameters, which renders it very predictive. Two features set this scenario
apart from related WIMP models. First, the dark matter stabilizing symmetry is
a $Z_{6}$ rather than the standard $Z_{2}$. Second, the mass term for the dark
matter particle is generated through the spontaneous breaking of the $Z_{6}$
symmetry. After studying its dark matter phenomenology, the viable parameter
space of this model was explored using a profile likelihood approach. Our
results indicate that this scenario remains viable over a wide range of dark
matter masses, extending well into the multi-TeV region. Three mechanisms can
successfully account for the observed dark matter abundance: annihilation at
the Higgs resonance ($M_{\psi}\approx 62.5$ GeV), annihilation at the new
scalar resonance ($M_{\psi}\approx M_{S}/2$), and annihilations into the new
scalars, $\psi\bar{\psi}\to SS$, ($M_{\psi}>M_{S}$). All of them will be
probed, either completely or partially, by future direct detection
experiments. Thus, this $Z_{6}$ scenario not only offers one of the minimal
dark matter extensions of the Standard Model, but also presents the most
compelling realization of singlet Dirac dark matter.

## Acknowledgments

OZ and CY received funding from the Patrimonio Autónomo - Fondo Nacional de
Financiamiento para la Ciencia, la Tecnología y la Innovación Francisco José
de Caldas (MinCiencias - Colombia) grant 82315-2021-1080. OZ has been
partially supported by Sostenibilidad-UdeA and the UdeA/CODI Grants
2020-33177, 2022-52380 and 2023-59130.

## References

  * Roszkowski et al. (2018) L. Roszkowski, E. M. Sessolo, and S. Trojanowski, Rept. Prog. Phys. 81, 066201 (2018), eprint 1707.06277.
  * Arcadi et al. (2018) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, Eur. Phys. J. C 78, 203 (2018), eprint 1703.07364.
  * Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], eprint 1807.06209.
  * Schumann (2019) M. Schumann, J. Phys. G 46, 103003 (2019), eprint 1903.03026.
  * Buchmueller et al. (2017) O. Buchmueller, C. Doglioni, and L. T. Wang, Nature Phys. 13, 217 (2017), eprint 1912.12739.
  * Kim et al. (2008) Y. G. Kim, K. Y. Lee, and S. Shin, JHEP 05, 100 (2008), eprint 0803.2932.
  * Baek et al. (2012a) S. Baek, P. Ko, and W.-I. Park, JHEP 02, 047 (2012a), eprint 1112.1847.
  * Baek et al. (2012b) S. Baek, P. Ko, W.-I. Park, and E. Senaha, JHEP 11, 116 (2012b), eprint 1209.4163.
  * Esch et al. (2013) S. Esch, M. Klasen, and C. E. Yaguna, Phys. Rev. D 88, 075017 (2013), eprint 1308.0951.
  * Freitas et al. (2015) A. Freitas, S. Westhoff, and J. Zupan, JHEP 09, 015 (2015), eprint 1506.04149.
  * Kim and Lee (2007) Y. G. Kim and K. Y. Lee, Phys. Rev. D 75, 115012 (2007), eprint hep-ph/0611069.
  * Fedderke et al. (2014) M. A. Fedderke, J.-Y. Chen, E. W. Kolb, and L.-T. Wang, JHEP 08, 122 (2014), eprint 1404.2283.
  * Matsumoto et al. (2014) S. Matsumoto, S. Mukhopadhyay, and Y.-L. S. Tsai, JHEP 10, 155 (2014), eprint 1407.1859.
  * Athron et al. (2019) P. Athron et al. (GAMBIT), Eur. Phys. J. C 79, 38 (2019), eprint 1808.10465.
  * Fairbairn and Hogan (2013) M. Fairbairn and R. Hogan, JHEP 09, 022 (2013), eprint 1305.3452.
  * Li and Zhou (2014) T. Li and Y.-F. Zhou, JHEP 07, 006 (2014), eprint 1402.3087.
  * Ghorbani (2017) P. H. Ghorbani, JHEP 08, 058 (2017), eprint 1703.06506.
  * Beniwal et al. (2019) A. Beniwal, M. Lewicki, M. White, and A. G. Williams, JHEP 02, 183 (2019), eprint 1810.02380.
  * Gould and Xie (2023) O. Gould and C. Xie, JHEP 12, 049 (2023), eprint 2310.02308.
  * Yaguna and Zapata (2020) C. E. Yaguna and O. Zapata, JHEP 03, 109 (2020), eprint 1911.05515.
  * Bélanger et al. (2020) G. Bélanger, A. Pukhov, C. E. Yaguna, and O. Zapata, JHEP 09, 030 (2020), eprint 2006.14922.
  * Belanger et al. (2022) G. Belanger, A. Mjallal, and A. Pukhov, Phys. Rev. D 105, 035018 (2022), eprint 2108.08061.
  * Yaguna and Zapata (2021) C. E. Yaguna and O. Zapata, JHEP 10, 185 (2021), eprint 2106.11889.
  * Yaguna and Zapata (2022) C. E. Yaguna and O. Zapata, Phys. Rev. D 105, 095026 (2022), eprint 2112.07020.
  * Choi et al. (2021) S.-M. Choi, J. Kim, P. Ko, and J. Li, JHEP 09, 028 (2021), eprint 2103.05956.
  * Díaz Sáez and Ghorbani (2023) B. Díaz Sáez and K. Ghorbani, JCAP 02, 002 (2023), eprint 2203.09282.
  * Jurčiukonis and Lavoura (2023) D. Jurčiukonis and L. Lavoura, PTEP 2023, 023B02 (2023), eprint 2210.12133.
  * Bélanger et al. (2023) G. Bélanger, A. Pukhov, C. E. Yaguna, and O. Zapata, JHEP 03, 100 (2023), eprint 2212.07488.
  * Dupuis (2016) G. Dupuis, JHEP 07, 008 (2016), eprint 1604.04552.
  * Bhattacharya et al. (2018) S. Bhattacharya, I. de Medeiros Varzielas, B. Karmakar, S. F. King, and A. Sil, JHEP 12, 007 (2018), eprint 1806.00490.
  * Ho et al. (2022) S.-Y. Ho, P. Ko, and C.-T. Lu, JHEP 03, 005 (2022), eprint 2201.06856.
  * Yaguna and Zapata (2024) C. E. Yaguna and O. Zapata, Phys. Rev. D 109, 015002 (2024), eprint 2308.05249.
  * Aad et al. (2015) G. Aad et al. (ATLAS, CMS), Phys. Rev. Lett. 114, 191803 (2015), eprint 1503.07589.
  * Bélanger et al. (2018) G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov, and B. Zaldivar, Comput. Phys. Commun. 231, 173 (2018), eprint 1801.03509.
  * Semenov (2016) A. Semenov, Comput. Phys. Commun. 201, 167 (2016), eprint 1412.5016.
  * Aprile et al. (2022) E. Aprile et al. (XENON), Phys. Rev. Lett. 129, 161805 (2022), eprint 2207.11330.
  * Aalbers et al. (2023) J. Aalbers et al. (LZ), Phys. Rev. Lett. 131, 041002 (2023), eprint 2207.03764.
  * AbdusSalam et al. (2022) S. S. AbdusSalam et al., Rept. Prog. Phys. 85, 052201 (2022), eprint 2012.09874.
  * Storn and Price (1997) R. Storn and K. Price, Journal of Global Optimization 11, 341 (1997).
  * Virtanen et al. (2019) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al., CoRR abs/1907.10121 (2019), eprint 1907.10121, URL http://arxiv.org/abs/1907.10121.
  * Wilks (1938) S. S. Wilks, Annals Math. Statist. 9, 60 (1938).
  * Falkowski et al. (2015) A. Falkowski, C. Gross, and O. Lebedev, JHEP 05, 057 (2015), eprint 1502.01361.
  * Arcadi et al. (2020) G. Arcadi, A. Djouadi, and M. Raidal, Phys. Rept. 842, 1 (2020), eprint 1903.03616.
  * Ferber et al. (2023) T. Ferber, A. Grohsjean, and F. Kahlhoefer (2023), eprint 2305.16169.
  * ATL (2022) Nature 607, 52 (2022), [Erratum: Nature 612, E24 (2022)], eprint 2207.00092.
  * Tumasyan et al. (2022) A. Tumasyan et al. (CMS), Nature 607, 60 (2022), eprint 2207.00043.
  * Ackermann et al. (2015) M. Ackermann et al. (Fermi-LAT), Phys. Rev. Lett. 115, 231301 (2015), eprint 1503.02641.
  * Reinert and Winkler (2018) A. Reinert and M. W. Winkler, JCAP 01, 055 (2018), eprint 1712.00002.
  * Aguilar et al. (2016) M. Aguilar et al. (AMS), Phys. Rev. Lett. 117, 091103 (2016).
  * Aalbers et al. (2016) J. Aalbers et al. (DARWIN), JCAP 11, 017 (2016), eprint 1606.07001.
  * Lopez-Honorez et al. (2012) L. Lopez-Honorez, T. Schwetz, and J. Zupan, Phys. Lett. B 716, 179 (2012), eprint 1203.2064.
  * Bernal et al. (2019) N. Bernal, D. Restrepo, C. Yaguna, and O. Zapata, Phys. Rev. D 99, 015038 (2019), eprint 1808.03352.


IEEEexample:BSTcontrol

# Echocardiography Segmentation Using Neural ODE-based Diffeomorphic
Registration Field

Phi V. Nguyen  Hieu H. Pham  Long Q. Tran Manuscript received November xx,
2021; revised April xx, 2022. This work was supported in part by the xx
*Correspondence authorsThis work was supported by the Vingroup Innovation
Foundation (VINIF) under project code VINIF.2019.DA02 and is also supported by
the VinUni-Illinois Smart Health Center at VinUniversity. Phi V. Nguyen, Long
Q. Tran are with the Institute for Artificial Intelligent, University of
Engineering and Technology, Vietnam National University, Vietnam.Hieu H. Pham
are with the College of Engineering & Computer Science, VinUni-Illinois Smart
Health Center, VinUniversity, Vietnam and the Coordinated Science Laboratory,
University of Illinois Urbana-Champaign, USA.Copyright ©2022 IEEE. Personal
use of this material is permitted. However, permission to use this material
for any other purposes must be obtained from the IEEE by sending a request to
<EMAIL_ADDRESS>
###### Abstract

Convolutional neural networks (CNNs) have recently proven their excellent
ability to segment 2D cardiac ultrasound images. However, the majority of
attempts to perform full-sequence segmentation of cardiac ultrasound videos
either rely on models trained only on keyframe images or fail to maintain the
topology over time. To address these issues, in this work, we consider
segmentation of ultrasound video as a registration estimation problem and
present a novel method for diffeomorphic image registration using neural
ordinary differential equations (Neural ODE). In particular, we consider the
registration field vector field between frames as a continuous trajectory ODE.
The estimated registration field is then applied to the segmentation mask of
the first frame to obtain a segment for the whole cardiac cycle. The proposed
method, Echo-ODE, introduces several key improvements compared to the previous
state-of-the-art. Firstly, by solving a continuous ODE, the proposed method
achieves smoother segmentation, preserving the topology of segmentation maps
over the whole sequence (Hausdorff distance: 3.7–4.4). Secondly, it maintains
temporal consistency between frames without explicitly optimizing for temporal
consistency attributes, achieving temporal consistency in 91% of the videos in
the dataset. Lastly, the proposed method is able to maintain the clinical
accuracy of the segmentation maps (MAE of the LVEF: 2.7–3.1). The results show
that our method surpasses the previous state-of-the-art in multiple aspects,
demonstrating the importance of spatial-temporal data processing for the
implementation of Neural ODEs in medical imaging applications. These findings
open up new research directions for solving echocardiography segmentation
tasks.

CNN, neural ordinary differential equation, cardiac segmentation, ultrasound,
left ventricle.

##  1 Introduction

Echocardiography is one of the most popular non-invasive procedures for
identifying structural abnormalities of the heart. It is inexpensive, requires
no incisions, and can be performed quickly, enabling cardiologists to diagnose
heart disease or assess a patient’s risk of future cardiac events [1]. Despite
its ease of use and cost-effectiveness, ultrasound is well known for its noisy
scans and variable image quality, which can depend on factors such as the
skill of the person operating the device [2], making it challenging for
automated systems to utilize the imaging data for patient diagnosis.

However, with the development of deep neural networks, particularly
convolutional neural networks (CNNs), outstanding results have been achieved
on multiple cardiovascular assessments, such as left ventricle segmentation or
ejection fraction [3, 4, 5]. Despite obtaining acceptable results on these
evaluations, these models are still hindered in accurately defining two
crucial phases of the cardiac cycle, namely end-diastole (ED) and end-systole
(ES), due to the inherent noise present in ultrasound images and the high cost
associated with annotating datasets. This limitation adversely affects the
algorithm’s ability to accurately segment the whole cardiac cycle, which is
crucial for evaluating other parameters, such as global longitudinal strain
[6, 7].

While current methods for segmenting the cardiac cycle in echocardiography
sequences have been introduced in the literature [4, 8, 9], they have several
major limitations. For instance, approaches presented in [4, 9] either process
individual frames or regularize the segmentation maps generated by models,
failing to consider the crucial temporal relationship between consecutive
frames in a video sequence. In contrast, [8] addresses the temporal
relationship between consecutive frames and utilizes a displacement field to
regularize the segmentation maps, but it does not ensure the consistency and
topology of the segmentation maps across the entire cardiac cycle.

To address these limitations, we present in this paper a novel approach for
segmenting the cardiac cycle in echocardiography sequences. This approach
allows generating a displacement field and propagating the segmentation map
using a Neural ODE [10]. The motivation behind using the Neural ODE for
segmentation in echocardiography video is to generate a smooth and topology-
preserving registration from an ODE and to take advantage of the ability to
learn dynamic and complex data from deep neural networks. Furthermore,
incorporating segmentation from a 2D CNN model enables the production of a
registration that includes spatial features from a conventional 2D
segmentation model while also benefiting from the improved temporal
representation provided by the Neural ODE. Additionally, by incorporating
clinical loss functions, such as the mean absolute error (MAE) of the left-
ventricular ejection fraction (LVEF), into the Neural ODE model, this approach
can produce estimations of clinical parameters that are as accurate as those
produced by the conventional 2D segmentation model. Finally, learning to
estimate registration for propagating segmentation maps represents the
underlying motion of the left ventricle chamber, which is beneficial for other
clinical diagnoses in echocardiography imaging.

Our contributions: We summarize our main contributions as follows:

  * •

We present a novel echocardiography video segmentation method that propagates
the segmentation map through the registration field obtained from a Neural
ODE, Echo-ODE. The proposed approach optimizes both spatial and temporal
information to effectively perform echocardiography segmentation;

  * •

We introduce collection of loss function and regularization techniques that
not only enhance segmentation performance but also improve the approach’s
effectiveness in terms of anatomical, cardiac motion, topology preservation,
and clinical aspects;

  * •

The proposed approach serves as post-processing for a trained segmentation
model and optimizes a single video sequence, thereby eliminating the need for
any additional training data;

  * •

Extensive experiments on the CAMUS dataset demonstrate that our approach
achieves state-of-the-art performance while showing interesting properties,
e.g., preserving desirable topology and registration.

The rest of this paper is organized as follows: Previous works related to
echocardiography segmentation are reviewed and discussed in Section 2. Section
4 explains the proposed network architecture that is able to perform the
diffeomorphic registration for echocardiography segmentation. Section 5
describes the dataset used in this study and the experimental settings.
Experimental results and discussions are presented in Section 6. We also
conduct in this section various experiments to investigate the model’s
properties. Finally, Section 7 concludes the paper with future research
directions.

##  2 Previous Works

This section reviews recent literature on full cardiac cycle segmentation,
especially on left ventricle segmentation and the Neural ODE. Details are
below.

###  2.1 Echocardiography Segmentation

####  2.1.1 Keyframe Segmentation

Because of the expense of collecting and accurately annotating
echocardiography data, previous works only focused on performing
echocardiography segmentation using ES and ED images [3, 5]. The first notable
data collected were CAMUS by Leclerc et al. [3] in 2019, which has 500
patients and was annotated by expert cardiologists. Within the same year,
Ouyang et al. [5] published EchoNet which was designed for segmenting the left
ventricle in A4C view sequences and estimating the ejection fraction (EF).

With the availability of large datasets, many conventional segmentation
architectures were being assessed. In the same CAMUS work, Leclerc et al. [3]
evaluated various U-Net architect [11] in terms of segmentation of the left
ventricle (LV), myocardium (MYO) and left atrium (LA) for the ED and ES frames
of the cardiac cycle for apical 4 chamber (A4C) and apical 2 chamber (A2C)
views. Ouyang et al. [5] also provided the segmentation result of DeepLabv3
[12] architecture to the A4C view sequence. Beside experimenting with
universal segmentation architectures, many methods were also proposed to
specifically tackle problems on echocardiography datasets. Veni et al. [13]
used a shape-driven deformable model with features provided by deep learning
networks to regulate the segmentation mask. Thomas et al. [14] proposed a
novel model based on Graph Convolutional Networks, which was used to detect
keypoints from ultrasound images. The result was then used to segment the LV
and estimate the LVEF.

Despite their success, the limitations of the above segmentation methods for
the LV were inherent in their sole consideration of the ED and ES phases,
ignoring the valuable temporal information that is crucial for accurate
estimation of the EF, as demonstrated in work by Wei et al. [8].

####  2.1.2 Spatial Temporal Segmentation

Recent papers [4, 15, 8] proposed specialized frameworks to focus on the
temporal aspect of echocardiography segmentation. PV-LVNet [4] used temporal
information extracted from recurrent networks to locate and crop the LV across
whole US sequences. MV-RAN [15] used both temporal and multi-view information
to perform segmentation of the LV with a convolutional LSTM to capture
temporal information. Wei et al. [8] introduced CLAS, a co-learning shape and
appearance model for temporal-consistent heart segmentation of
echocardiography sequences with sparsely labeled data. The framework learned
the 2D+time apical long-axis cardiac shape such that the segmented sequences
can benefit from temporal and spatial consistency. While relying solely on ED
and ES annotations, CLAS [8] also predicted deformation fields and used this
information to regularize segmentation during training instead of propagating
the segmentation to the whole sequence.

Another approach is the refinement of segmentation masks. Painchaud et al. [9]
utilized AR-VAE [16] to post-process the latent representation of segmentation
for every frame in the cycle to achieve the desired temporal consistency. The
methods helped achieve SOTA results in LV segmentation and interpretable
embedding space. However, the method did not account for temporal features
from ultrasound but only used segmentation results from a pretrained model,
missed crucial temporal information presented in ultrasound sequence.

We address in this work the limitations of previous approaches by introducing
a novel method, Echo-ODE, for echocardiography segmentation that leverages
Neural ODE to model the temporal feature between frames, resulting in smooth
and accurate segmentation. While prior research has focused primarily on
keyframe segmentation and neglected valuable temporal information from
ultrasound sequences, the method exploits the unique capabilities of Neural
ODE to effectively capture the dynamic motion of the LV throughout the entire
cardiac cycle.

###  2.2 Deformable Registration Learning

In the field of medical imaging, deformable image registration is used to find
a dense, non-linear correspondence transform that can smoothly change one
image or video into another [17]. Deformable registration for echocardiography
can be used to register 3D models to 3D points, quantify myocardial motion and
strain from multiple 3D ultrasound sequences, and segment the LV in 3D
echocardiography [18].

With the development of deep learning algorithms, researchers in computer
vision and medical imaging are turning to deep neural networks for solving
problems in the field of deformable image registration [17]. One of the
primary deep learning architectures used for deformable registration tasks is
VoxelMorph [19], which used a convolutional neural network to automatically
create an accurate deformation field between two medical scans. Adversarial
similarity networks were also used to improve the quality of the registration
field by calculating the appearance loss between a warped image and a fixed
image. In echocardiography processing, deformable image registration is also
used to estimate the flow of motion between ultrasound images. A notable work
by Zhu et al. [20] that used U-Net [11] to measure blood flow between
ultrasound images at test-time training.

Solving for diffeomorphisms by putting the registration problem into the form
of an ordinary differential equation (ODE) whose solution is a smooth
deformation field is a promising direction. Dalca et al. [21] extended
VoxelMorph by forming the registration field as a diffeomorphic field of an
ODE. By constraining the transformations to diffeomorphisms, the model from
Dalca et al. [21] could ensure certain desirable properties when transforming
the image, such as the preservation of topology. The ODE itself could later be
approximated by a neural network, which is more flexible and tractable for
large-scale problems [10]. Xu et al. [22] introduced a multi-scale Neural ODE
for estimating registration of multi-contrast 3D MR images. Joshi et al. [23]
utilized Lipschitz continuous ResBlocks to solve the ODE that governs
diffeomorphic deformations. Despite the success of these methods on different
modalities, as far as we know, there is no work that uses diffeomorphic
deformations from Neural ODE to solve the problem of echocardiography
segmentation.

Our proposed method is an attempt to overcome the limitations of existing SOTA
techniques. Ultrasound images are inherently noisy, which makes solving the
Neural ODE without regularization and guidance inadequate for ensuring smooth
and accurate registration while preserving the topology of segmentation maps
throughout the entire sequence. To address this issue, we employ a balance
between temporal consistency and segmentation map accuracy by incorporating
both temporal features from the ultrasound sequence and segmentation maps from
a pretrained segmentation model. Moreover, the proposed method integrates a
clinical loss function, which facilitates the preservation of clinical
accuracy for relevant metrics.

Figure 1:  An overview of Echo-ODE for segmenting echocardiography. Taking the
echocardiography sequence $I$ and two segmentation map $m_{0}$ and $m_{1}$ as
a starting point. At time step $t$, the the model encodes the ultrasound
image, $I_{t}$, and decode it with an encoding of $t$ by a positional encoding
in latent space. The model then generates an initial vector field by spatially
adaptive normalizing the decoder features with positional embedding of $t$ and
segmentation map $m_{0}$. The ODE solves for the deformation vector field
$\phi_{t}$, used to generate a pseudo-image $\hat{I}_{t}$ and a pseudo-
segmentation map $\hat{m}_{t}$.

##  3 Preliminaries

In this section, we provide an overview of the video segmentation and
diffeomorphic registration methods that are the foundation of the proposed
approach.

###  3.1 Echocardiography Segmentation

In echocardiography segmentation, given a sequence of ultrasound images from
ED to ES, we have $I=\\{I_{0},\dots,I_{t},\dots,I_{1}|t\in[0,1]\\}$, where
$I_{0}$ denotes the B-mode ultrasound image at the beginning of the cardiac
cycle (ED), $I_{1}$ denotes the B-mode ultrasound image at the end of the
cardiac cycle (ES), and $I_{t}$ denotes the intermediate image between ED and
ES. Specifically, when $t\rightarrow 1$ the heart goes through systolic phase,
undergoes contraction and ejects blood into the arteries. Each B-mode
ultrasound image $I_{t}:\Omega\rightarrow\mathbb{R}$, where $\Omega$ is the
image domain, has the segmentation map
$\hat{m}_{t}:\Omega\rightarrow\mathbb{Z}$, in which each pixel
$\hat{m}_{t}(p)=k$ corresponds to an object of interest, e.g., $k=1$
represents the LV region.

###  3.2 Deformation Registration

Deformable registration aims to find an optimal transformation function that
maps every point in the source image to a corresponding point in the target
image, under some constraints of smoothness and preservation of topological
properties. The transfer of segmentation maps is then accomplished by applying
this deformation field to the source segmentation, effectively repositioning
and reshaping the segmented structures to align with the corresponding
structures in the target image.

The deformable registration in our study could be formulated as follows. Given
two B-mode ultrasound images at two time frames, beginning of cardiac cycle
and arbitrary time, 0 and $t$, represented as $I_{0},I_{t}$ . Let each pixel
location be represented by $p\in\Omega\subseteq\mathbb{R}^{2}$. The goal of
deformation registration is to find a transform
$\phi_{t}:\Omega\rightarrow\Omega$ that minimizes a loss function
$\mathcal{L}(I_{t},\Phi_{t}(I_{0}))$, where
$\Phi_{t}(I_{0})(p)=I_{0}(\phi_{t}(p))$. Here $\mathcal{L}$ measures the
dissimilarity between the image $I_{0}$ transformed by $\phi_{t}$ and the
image $I_{t}$ under some regularization constrains $\mathcal{R}$. Given the
segmentation map $m_{0}$ of $I_{0}$, our goal is to apply the previously
computed transformation $\phi_{t}$ to transfer the segmentation map from the
beginning of the cardiac cycle to an arbitrary time $t$. This process results
in a new segmentation map $\hat{m_{t}}$ that aligns with the image $I_{t}$.

###  3.3  Diffeomorphic Registration Using Neural ODE

To ensure the smoothness and topology preservation properties of the
deformation, we define the deformation field as an ODE. Following the works of
Rousseau et al.[24] and Ashburner [25], we parameterize the deformation
registration $\phi_{t}$ by the ODE:

$\frac{\mathrm{d}\phi_{t}(p)}{\mathrm{d}t}=V_{\theta}(p,t),$ (1)

where $\phi_{t}$ is the diffeomorphism at time $t$, and $V_{\theta}$ is the
stationary velocity vector field parameterized by $\theta$, which is computed
by a 2D convolutional neural network. By integrating the velocity field from
time 0 to $t$, we obtain the deformation at any time $t$ ($\phi_{t}$),
starting with the identity $\phi_{0}(p)=p$:

$\phi_{t}(p)=\phi_{0}(p)+\int_{0}^{t}{V_{\theta}(\phi_{s}(p),s)}ds.$ (2)

In group theory, integrating this ODE to compute $\phi_{t}$ is analogous to
computing the exponential map of the velocity field $V_{\theta}$ (which can be
considered an element of a Lie algebra associated with a Lie group [26]). The
inverse transform $\phi_{t}^{-1}$ is obtained by integrating the negative of
the velocity field.

##  4 Proposed Method

Our method focuses on echocardiography segmentation, where we employ
diffeomorphic registration as a tool to improve the segmentation process. The
method measures a diffeomorphic registration between first frame and arbitrary
frame and uses the registration to propagate the initial segmentation map
throughout the entire sequence. An overview of our method is visually
represented in Fig. 1.

###  4.1 Learning registration field

The deformation vector field $\phi_{t}$ must satisfy three criteria: (1) the
deformation should be smooth; (2) the deformation vector field should
represent the motion present in the ultrasound sequence; and (3) the topology
of segmentation must be preserved after individual transformations, while the
semantic information of the target segmentation, as learned by pre-trained
segmentation models, must be accurately reflected in the applied
transformations. Our method is illustrated in Fig. 1 and detailed in Algorithm
1. We proposed a auto encoder model based on UNet [11] to effectively estimate
the velocity field $V_{\theta}(p,t)$. The input for this model is the image
$I_{t}$, segmentation map of the first frame $m_{0}$, and timestep $t$. The
objective of the model is to estimate the registration field $V$ for ODE
solver.

Encoder-Decoder: At the encoder module, the model encodes the image into a
latent vector $v$ as a conventional auto-encoder model. At the decoder module,
we augment the latent vector $v$ with cosine positional encoding for timestep
$t$, as described in [27], in order to maintain the order of the frame
sequence, ensuring that the temporal relationships between frames are
preserved. We also normalize decoder features to guide the flow field to
follow LV, as described in the following section.

Algorithm 1 Learning diffeomorphism $\phi:I_{0}\rightarrow I_{1}$

1:Sequence images $I=\\{I_{t}\\}_{t\in[0,1]}$, 2D model’s segmentation maps
$m=\\{m_{0},m_{1}\\}$

2:Registration fields $\phi=\\{\phi_{t}\\}_{t\in[0,1]}$ $\implies$ full
cardiac segmentation $\hat{m}=\\{\hat{m}_{t}\\}_{t\in[0,1]}$

3:parameterized by neural networks: $V_{\theta}$ in (1).

4:while V has not converged do

5: Set $\phi_{0}$ as identical transformation (vector field).

6: for $i=0,\ldots,1$ do

7: Using Neural ODE($V_{\theta},\phi_{0})$ find a solution $\phi_{t}$.

8: Given $\phi_{t}$, warp $I_{0}$ to $\hat{I}_{t}$

9: end for

10: Given $\phi_{1}$, warp $m_{0}$ to $\hat{m}_{1}$

11: update $V_{\theta}$ by minimizing ODE loss:
$\mathcal{L}_{sim}(I,\hat{I})+\mathcal{L}_{seg}(m_{1},\hat{m_{1}})+\mathcal{L}_{EF}+\mathcal{R}_{Jdet}(\phi)+\mathcal{R}_{smt}(\phi)$

12:end while

Temporal Spatial Adaptive Normalization: Due to the fact that motion change
between frames happens around the edge of anatomical structure and the
correctly segments only happen on keyframe images, we extend SPADE [28] toward
spatially-adaptive and time-dependent normalization by adding positional
encoding of $t$, so the model can dynamically be guided by segmentation maps
in the early part of the cardiac cycle, and use features from ultrasound later
on. Specifically, in order to estimate the registration $V$, for each decoder
feature map $f_{n}\subseteq\mathbb{R}^{H_{n}\times W_{n}\times C_{n}}$, we use
the segmentation maps of the first frame $m_{0}$ and position encoding of
timestep $t$ as inputs to normalize feature map $f_{n}$. We estimate two
normalization parameters $\gamma_{n},\beta_{n}\subseteq\mathbb{R}^{H_{n}\times
W_{n}\times C_{n}}$. The normalized feature is given by:

$f_{n}=\gamma_{n}\otimes f_{n}\oplus\beta_{n},$

where $\otimes$ is element-wise multiplication and $\oplus$ is element-wise
addition. This procedure is applied to each frame of the cardiac cycle, and
the results are integrated by ODE solver, resulting in the following sequence
of diffeomorphic registrations: $\phi=\\{\phi_{t}|i=0,...,1\\}$.

###  4.2 Loss Function and Regularization

Following the work of Wu et al. [29], we design Algorithm 1 in order to find
the optimal parameter $\theta$ that is satisfied:

$\begin{array}[]{r}\theta=\operatorname*{argmin}_{\theta}(\mathcal{L}_{sim}+\mathcal{L}_{seg}+\mathcal{L}_{EF}+\mathcal{R}_{Jdet}+\mathcal{R}_{smt}),\end{array}$
(3)

Similarity loss is used to minimize the discrepancy between the pseudo-images
$\hat{I}$ and the ground truth $I$. Given the prevalence of speckle noise in
ultrasound images, we apply Gaussian kernels $\mathcal{K}$ (for individual
images in the sequence) to both pseudo and ground truth images to smooth out
images. The mean square error for each pixel is employed to minimize the
difference between pseudo images $\hat{I}$ and the ground truth $I$, which is
calculated as follows:

$\mathcal{L}_{sim}=\mathbb{E}_{t\sim\mathcal{U}(0,1)}(I_{t}-\hat{I}_{t})^{2}.$
(4)

Segmentation loss To enable the registration process to concentrate on the
segmentation problem, we propose a secondary segmentation loss between the
propagated segmentation maps and the segmentation maps generated by pre-
trained segmentation models on keyframe only. The rationale for incorporating
this additional loss is to guide the registration estimation to focus on the
keyframe, where pre-trained model performed well, instead of whole sequence.
The loss is calculated using the Dice coefficient as follows:

$\mathcal{L}_{seg}=2\frac{|\hat{m}_{1}\cap m_{1}|}{|\hat{m}_{1}|+|m_{1}|},$
(5)

where $\hat{m}_{1}$ is the segmentation map from the registration model and
$m_{1}$ is the segmentation map from the segmentation model at the ES frame.

Ejection Fraction loss The main purpose of doing segmentation is to measure
clinical parameters such as EF, LV mass, and LV volume. Painchaud et al. [9]
showed that forcing the temporal consistency between segmentation maps caused
the degradation of EF estimation in comparison to the segmentation maps from
the pretrained segmentation model. To mitigate this problem, we propose a
clinical attribute loss to minimize the difference between the EF obtained
from segmentation maps from the registration model and EF obtained from
segmentation maps from the pretrained segmentation model. The EF loss is
calculated by:

$\displaystyle\mathcal{L}_{EF}$ $\displaystyle=\left|\frac{EDV-
ESV}{EDV}-\frac{EDV_{\text{pre}}-ESV_{\text{pre}}}{EDV_{\text{pre}}}\right|$
$\displaystyle=1-\frac{\sum_{i}\sum_{j}\mathbbm{1}[\hat{m}_{1}(i,j)=\text{LV}]}{\sum_{i}\sum_{j}\mathbbm{1}[m_{1}(i,j)=\text{LV}]+\epsilon},$

where $EDV$ and $ESV$ are the end-diastolic volumes and end-systolic volumes
from the registration model, and $EDV_{pre}$ and $ESV_{\text{pre}}$ are the
end-diastolic volumes and end-systolic volumes from the pretrained
segmentation model. Since only a 4 chamber view is available for EF estimation
and $EDV=ESV_{\text{pred}}$, the volume is then approximated by the area (in
pixels) of the LV segmentation map. Beside loss functions, we further add
regularization terms to make sure the diffeomorphism vector field achieves
smoothness and topology-preserving properties.

Topology preserving regularizer The first regularization term is to ensure the
Jacobian determinant of each deformable vector field $\phi_{t}$ is non-
negative. The Jacobian determinant is positively making vector field
$\phi_{t}$ create no fold during transformation, thus preserving the topology
of the segmentation map after transformation by enforcing local orientation
consistency of each pixel concerning neighbor pixels. Similar to NODEO [29],
we minimize the negative elements of $\phi_{t}$ by apply ReLU activation
function and use L2 norm to avoid sparsity of L1 norm, the regularization
$\mathcal{R}_{Jdet}$ is given by:

$\mathcal{R}_{Jdet}=\mathbb{E}_{t\sim\mathcal{U}(0,1)}\left\|\textit{ReLU}\left(-\left(\left|\mathcal{D}_{\phi_{t}}(\mathrm{p})\right|\right)\right)\right\|_{2}^{2},$
(6)

Smooth regularizer The second term of regularization, $\mathcal{L}_{smt}$,
serves to promote the smoothness of the transformation. This is achieved by
minimizing the spatial gradients of the transformed image, thereby ensuring a
seamless transition between adjacent pixels:

$\mathcal{R}_{smt}=\mathbb{E}_{t,p}\left(\left\|\nabla_{\phi}(\mathrm{p})\right\|_{2}^{2}\right).$
(7)

##  5 EXPERIMENTS

###  5.1 Datasets and experimental settings

Due to the lack of echocardiography datasets contains segmentation maps
annotated for a full cardiac cycle. We evaluate the effectiveness of the
proposed method for echocardiography segmentation tasks using only the CAMUS
dataset [3, 9]. The first version of the CAMUS dataset provides annotations
for the ED and ES frames for each cardiac cycle. These frames were chosen to
represent the peak of systole and diastole, respectively. The annotations in
this version of the dataset encompass the locations of the left ventricle
endocardium ($LV_{endo}$), myocardium ($LV_{epi}$), and left atrium ($LA$)
boundaries. Designed for training segmentation models, the two-frame dataset
comprises a total of 450 videos, each capturing a complete cardiac cycle (from
ED to ES and from ES to ED). The dataset was subsequently partitioned into
training and validation sets at a 90% to 10% ratio. All single-frame
segmentation models employed to guide the ODE model were trained using this
dataset.

The second version of the CAMUS dataset [9] features full-cycle annotations
for echocardiography videos, with an emphasis on temporal segmentation of the
left ventricle. In particular, this dataset includes annotations for the
endocardium ($LV_{endo}$) and myocardium ($LV_{epi}$) of the left ventricle
for every frame from ED to ES. This version of the CAMUS dataset is utilized
to assess the performance of segmentation methods developed using the two-
frame dataset.

###  5.2 Pretrained Segmentation Model

We employed three distinct architectures to guide the learning of the
registration field, namely DeepLabv3 [12], U-Net[11], and ENet [30]. DeepLabv3
[12] is a fully convolutional neural network (FCN) that utilizes atrous
convolutions and a multi-scale feature fusion strategy to enhance semantic
segmentation performance. U-Net [11] is an FCN with a symmetric architecture
that employs skip connections to combine information from various levels of
resolution. ENet [30] is a lightweight network that leverages an encoder-
decoder architecture to strike a balance between accuracy and computational
cost.

###  5.3 Implementation Details and Training Procedure

The training process for the registration field was conducted in two stages.
The first stage involved training the segmentation models using the two-frame
dataset, aiming to generate accurate segmentation maps of the LV. The
segmentation maps created in the first stage served as a guide for estimating
the registration field in the second stage on the full-cycle dataset.

In the initial stage, we trained the segmentation models using the two-frame
dataset. All segmentation models were trained with the Dice loss function, a
commonly used loss function for segmentation tasks. The parameters of the
segmentation models were optimized using the Adam optimizer with a learning
rate of 0.01 and a batch size of 8. Training was terminated when the
validation loss ceased to improve. The segmentation models were trained for
100 epochs, utilizing a total of 720 training samples and 450 validation
samples. The segmentation models were also evaluated using the Dice similarity
coefficient (DSC) metric.

For the second stage, we estimated the registration field for each sample in
the full-cycle dataset. This means that each model implicitly represents only
a single echocardiography sequence. We used U-Net model architecture with an
encoder consisting of four convolutional layers with channel sizes [16, 32,
32, 32] and a decoder comprising seven convolutional layers with channel sizes
[32, 32, 32, 32, 32, 16, 16]. Each sample was optimized for 1000 steps. We
used Euler’s method with a step size of 0.1 to solve the ODE. The parameters
of the ODE model were optimized using the Adam optimizer with a learning rate
of 0.01.

Table 1: Performance of the methods on temporal consistency, anatomical accuracy, and clinical metric. The Hausdorff distance (HD), mean absolute error of ejection fraction (EF) are reported in terms of the mean and the standard deviation across 98 test sequences, we also report mean Dice score and omit the standard deviation since we notice no significant difference between methods. The temporal error (Temp. Error) is the number of sequences that have at least one temporal inconsistency. The mDice is the mean Dice score of the LV myocardium and LV epicardium. The Dice LVepi and Dice LVmyo are the Dice scores of the LV epicardium and LV myocardium, respectively. The best results are highlighted in bold. Methods |  | Temporal | Anatomical | Clinical
---|---|---|---|---
| HD | Temp. Error | mDice | Dice LVepi | Dice LVmyo | EF
CLAS [8] |  | 7.2 $\pm$ 5.2 | 98 | 0.911 | 0.898 | 0.925 | 6.0 $\pm$ 5.1
Baseline | U-Net [11] | 4.6 $\pm$ 4.7 | 98 | 0.953 | 0.943 | 0.964 | 3.2 $\pm$ 2.6
ENet [30] | 4.6 $\pm$ 3.1 | 98 | 0.951 | 0.943 | 0.960 | 2.9 $\pm$ 2.5
DeepLabv3 [12] | 5.6 $\pm$ 4.7 | 98 | 0.945 | 0.935 | 0.955 | 2.6 $\pm$ 2.3
Cardiac - VAE [9] | U-Net [11] | 4.2 $\pm$ 1.6 | 0 | 0.938 | 0.931 | 0.945 | 3.6 $\pm$ 3.3
ENet [30] | 4.3 $\pm$ 1.6 | 1 | 0.943 | 0.937 | 0.949 | 3.0 $\pm$ 2.4
DeepLabv3 [12] | 4.8 $\pm$ 1.5 | 0 | 0.943 | 0.937 | 0.949 | 3.3 $\pm$ 2.5
Echo-ODE (ours) | U-Net [11] | 3.7 $\pm$ 1.6 | 12 | 0.957 | 0.948 | 0.967 | 3.1 $\pm$ 2.6
ENet [30] | 3.9 $\pm$ 1.6 | 14 | 0.956 | 0.948 | 0.963 | 2.8 $\pm$ 2.5
DeepLabv3 [12] | 4.4 $\pm$ 1.8 | 9 | 0.950 | 0.940 | 0.959 | 2.7 $\pm$ 2.3

###  5.4 Evaluation Metrics

In this study, we employed a combination of three distinct metrics to evaluate
the effectiveness of our proposed method: (1) temporal smoothness and
consistency, (2) anatomical accuracy, and (3) clinical relevance. To assess
segmentation quality, we report DSC for left ventricle endocardium
($LV_{endo}$), left ventricle myocardium ($LV_{epi}$), and mean DSC of those
classes. In order to assess the smoothness of the segmentation maps, we
calculate the average Hausdorff Distance (HD) between consecutive frames
across the entire sequence. Although our model was not explicitly optimized to
minimize temporal inconsistency, we also reported the seven inconsistency
indicators proposed by [9] to enable fair comparisons with other methods. The
attributes used to characterize a 2D cardiac shape include: i) LV area, ii) LV
width at the valves, iii) LV length from valves to apex, iv) orientation of
the LV principal axis, v) myocardium (MYO) area, vi) horizontal center of mass
of the epicardium (EPI), and vii) vertical center of mass of the EPI. To
evaluate the clinical relevance of our approach, we used MAE of EF between
segmentation maps generated by our method and ground truth, to assess the
performance of our method in providing clinically relevant information. The
temporal error indicates the number of sequences with at least one temporal
attribute inconsistency.

###  5.5 Baselines

We conducted a comparative analysis of the proposed method against two
categories of models: 1) popular 2D segmentation models, and 2) specialized
temporal segmentation models for echocardiography. The popular 2D segmentation
models that we compared against are U-Net [11], ENet [30], and DeepLabv3 [12].
In contrast, the specialized segmentation models we compared against are CLAS
[8] and Cardiac AR VAE. CLAS [8] is a 2D segmentation model that is
specifically designed for echocardiography segmentation. It employs an
auxiliary registration loss to guide the segmentation model. On the other
hand, Cardiac AR VAE is a temporal segmentation model that uses an Attribute-
Regularized VAE (AR-VAE) to encode segmentation maps, utilizes a latent space
that contains seven interpretable attributes to model the temporal dynamics of
the left ventricle and nine residual attributes to encode additional
information. The smooth segmentation maps is obtained by smooth out those
attributes.

To evaluate the performance of the proposed and baseline methods, we trained
the models on a two-frame dataset and evaluated them on a full-cycle dataset.
The evaluation is conducted on 98 test sequences, and the results are reported
in terms of mean and standard deviation over all test samples. We used default
parameters from the original paper to train CLAS [8] and Cardiac AR VAE. All
the baseline and proposed models were implemented using PyTorch and trained on
a single NVIDIA Ampere A100 GPU.

##  6 Results

Table 1 presents the performance of five different methods on three evaluation
metrics: temporal smoothness and consistency, anatomical accuracy, and
clinical metric. In general, the Echo-ODE method outperforms the other methods
in terms of temporal consistency, with the lowest Hausdorff distance of 3.7.
Moreover, it achieves the highest mean Dice score of 0.967 for the LV
myocardium and 0.957 for the LV epicardium, respectively. In terms of clinical
metrics, our method achieved the lowest mean absolute error of ejection
fraction of 2.7, and while not surpassing baseline models, our method still
performs better than the current temporal SOTA. Overall, the Echo-ODE method
shows promising results in accurately segmenting cardiac structures and
ensuring temporal consistency in echocardiography sequences. The following
sections discuss more details on different aspects of the results.

###  6.1 Improvement over Anatomical Segmentation

The results of methods on segmentation metrics are presented in Table 1. While
accurate segmentation is important, we believe that temporal smoothness and
accuracy are equally critical factors that contribute to the effectiveness of
our approach, and maintaining a consistent level of segmentation quality
across different time points and images is essential for ensuring reliable
diagnoses and treatments [6]. The proposed method achieved an average Dice
coefficient at best of a mean Dice score of 0.957 for using U-Net [11] as
segmentation guidance, which is slightly better than the 2D segmentation
models, which achieved a mean Dice score ranging from 0.945 to 0.953. Compared
to models that utilized temporal information, our method outperforms Cardiac
AR-VAE and CLAS [8] by a substantial margin. Furthermore, we also evaluated
the results on different sub-regions of the $LV_{epi}$ and $LV_{myo}$, our
model obtained an average Dice coefficient at best of 0.948 for $LV_{epi}$ and
0.967 for $LV_{myo}$, respectively, which further demonstrate the improvement
in anatomical segmentation over baselines and previous SOTA using temporal
information with only 0.943 for $LV_{epi}$ and 0.964 for $LV_{myo}$ at best.
We also notice that the increase in segmentation accuracy is proportional to
the accuracy of pre-trained segmentation models. For example, 2D U-Net [11]
achieves a better mean Dice score of 0.953 when compared to ENet [30] and
Deeplabv3, and our method that utilizes U-Net [11] also achieves better
performance than ENet [30] and Deeplabv3. In a later section, we further
investigate the impact of initial segmentation.

Figure 2: Qualitative comparisons of different methods on the typical frames
from CAMUS datasets. Different color represents different semantic classes:
$LV_{myo}$ (teal), $LV_{epi}$ (yellow) and background (purple). Both Cardiac
AR-VAE and Echo-ODE use the same U-Net [11] as input. More results can be
found in the supplementary materials.

###  6.2 Segmentation Smoothness and Temporal Consistency

Assessing the temporal consistency of the segmentation maps is important for
the clinical application of the segmentation models. The results are presented
in the first column of Table 1. The Echo-ODE method outperforms the other
methods in terms of temporal smoothness with the lowest Hausdorff distance.
Our method and Cardiac AR-VAE also produce segmentation maps with a lower
standard deviation of Hausdorff distance of around 1.6, compared to 3.1-4.7 in
baseline models, indicating that the temporal models produce consistent
results. The results show that the proposed method achieved a mean Hausdorff
distance of 3.7 to 4.4, which is lower than 4.2-4.8 for Cardiac AR-VAE and 7.2
for CLAS [8], indicating the effectiveness of modeling the motion using neural
ODE. Furthermore, we also detected 9 to 12 sequences that have at least one
temporally inconsistent frame out of 98 sequences. While the temporal
consistency of our method is not as good as Cardiac AR-VAE in the definition
of Painchaud et al. [9], it is still better than 2D segmentation methods and
CLAS [8]. It is important to note that Cardiac VAE achieves these results by
explicitly optimizing for the temporal consistency attribute, while optimizing
the whole latent space of VAE does not yield such good results. The proposed
method, on the other hand, does not require any assumptions from the data
modality, which is better suited for a broader range of use cases.

(a) Highlight shape difference of segmentation maps, where the left half of
endocardium of Cardiac AR-VAE appears to have different shape.

(b) LV length from apex point to mitral valve of different models.

Figure 3: Visually and qualitative demonstration of different segmentation
maps from different segmentation models.

###  6.3 Improvement over Clinical Metrics

The data in Table 1 indicates that our proposed methods surpass CLAS [8] and
Cardiac AR-VAE, which are specialized designed models, in terms of clinical
metric. Specifically, our approach achieved a mean absolute error of 2.7 to
3.1 for EF, which is notably lower than the 0.6-0.7 range of Cardiac AR-VAE
and 2.9-3.1 for CLAS [8], demonstrating the accuracy of our method for
estimating cardiac parameters. Additionally, our method does not experience a
decline in clinical metrics when smoothing segmentation maps, which is a
limitation of the Cardiac AR-VAE model. Nonetheless, we observed that our
method does not yield significantly higher scores than 2D segmentation models
in terms of clinical metrics, except for the U-Net [11] and ENet [30] cases,
which had a 0.1 higher score than the original 2D segmentation models. The
other two models, ENet [30] and Deeplabv3 [12] have slightly lower scores than
2D segmentation models. This is because our approach relies on segmentation
maps and EF of baseline models to provide guidance, which are not as precise
as ground truth segmentation maps from the outset. Future work will involve
developing a better registration method that does not require segmentation
maps for regularization.

###  6.4 Shape correctness and localization

When qualitatively evaluating the results of the previous SOTA, Cardiac AR-
VAE, we observed a difference in the shape of the segmentation maps before and
after being passed through the auto-encoder. This issue only occurred in the
VAE model and was not present in the proposed method or other segmentation
models. For instance, in Fig. 3(a), we can see that while the segmentation
results from the original 2D segmentation models and the Echo-ODE produced
shapes similar to the ground truth, the segmentation results from Cardiac AR-
VAE tended to smooth out the contour of the segmentation maps. This caused the
segmentation maps to diverge from the characteristic anatomical structure. To
further illustrate the differences, we also plotted the length of the LV’s
apex and center of the mitral valve from the segmentation maps over time for
different models in Fig. 3(b). The figure indicates that Cardiac AR-VAE had a
differing length of the LV, whereas our method had similar lengths and also
had a smooth attribute over time. The cause of this issue may be due to the
fact that the latent space of Cardiac AR-VAE encodes side information related
to the shape of segmentation maps, and the smoothing action on the latent
space over time may have inadvertently smoothed out the shape of the
anatomical structure. On the other hand, other methods, including ours, do not
require manipulation over the latent space and therefore do not suffer from
this phenomenon.

###  6.5 Impact of Initial Segmentation on Anatomical and Clinical Accuracy

Our method’s segmentation is obtained by propagating through the whole
sequence from the first ED frame. Concerns may arise about the accuracy of
segmentation maps from pre-trained models for the entire sequence. To assess
the impact of the initial segmentation on the final results, we compared our
method’s segmentation using the initial frame from pre-trained segmentation
models versus ground truth segmentation. Despite using the same registration
field, we reported both anatomical and temporal metrics in Table 2. The
findings show that using ground truth as the initial segmentation map leads to
better segmentation smoothness with a higher Hausdorff score of 0.1-0.6.
However, the impact on anatomical accuracy and clinical metrics is minimal.
The largest difference was observed with Deeplabv3 [12], where the Dice score
was 0.08 lower when using a pre-trained segmentation map compared to the
ground truth segmentation, while the MAE of EF remained the same. This
suggests that the registration field is more accurate and smooth when close to
the ground truth region, but using pre-trained segmentation maps as the
initial segmentation is still effective for our method.

Table 2: Performance of the proposed method on different segmentation model
after replace initial segmentation map by the ground truth. all the results
are reported in terms of mean across test set. The difference with using
segmentation map of correspond segmentation model is reported in parenthesis.

Segmentator | HD | mDice | Dice LVepi | Dice LVmyo | EF
---|---|---|---|---|---
U-Net [11] | 3.6 (-0.1) $\pm$ 1.7 | 0.959(+0.02) | 0.950(+0.02) | 0.968(+0.01) | 3.1(+0.00) $\pm$ 2.6
ENet [30] | 3.6 (-0.3) $\pm$ 1.7 | 0.959(+0.03) | 0.951(+0.03) | 0.967(+0.04) | 2.8(+0.00) $\pm$ 2.6
DeepLabv3 [12] | 3.8 (-0.6) $\pm$ 1.9 | 0.956(+0.06) | 0.948(+0.08) | 0.965(+0.04) | 2.7(+0.00) $\pm$ 2.4

###  6.6 Topology Preserving and Registration

This section presents the results of our method’s ability to perform
registration and preserve topology.

Fig. 2 shows qualitative outcomes of different approaches on different
sequence frames. U-Net [11], ENet [30], and Deeplabv3 2D segmentation
algorithms preserve the LV’s shape but not its topology or smoothness. The
CLAS [8] model maintains the shape correctness and smoothness of segmentation
maps in early cardiac cycle but fails to preserve its topology later on. The
Cardiac AR-VAE is more stable when the sequence is long but loses its LV shape
after encoding and decoding the segmentation map. Our method, in the other
hand, both preserves LV topology and segmentation map smoothness. Fig. 4
demonstrates the registration field and how it transforms the image. For the
first frame, the registration field is identity and the grid is not deformed,
making the ultrasound identical to the original image. In the last frame, the
ultrasound is deformed towards the real ultrasound image. The grid shows the
topology preservation of the registration, where there is minimal folds or
holes in the grid. It is important to note that the topology preservation of
the registration field are obtained by minimizing the Jacobian determinant of
the registration field. Therefore, there is no guarantee that the registration
field will preserve the segmentation map’s topology at all times, but it
trying to preserve the topology of the registration field as much as possible.
The method only ensure the topology of propagation, the correctness of
segmentation maps would depend on the correctness of initial segmentation map.
By using guidance from segmentation, we can see from the grid that the
registration field follows the motion pattern of the LV, and the model
attempts to register the ultrasound by following the speckle of the LV. This
property allows us to use the registration field for other tasks, such as
tracking the motion of the LV.

Figure 4: Topology preserving. The registration results of our method (left)
when applied to the first ultrasound frame of the sequence. The deformed grid
(right) when applying the registration field (middle). The ultrasound frame is
jet-colored for easier visualization of the speckle noise.

###  6.7 Potential Limitations and Weaknesses

Our method demonstrates promising results but has limitations, as visualized
in Fig. 5. Firstly, there’s a trade-off between learning complex registration
fields and contour continuity, with noise in fast-motion cases. This issue
stems from pixel-level loss functions; future work will explore alternative
loss functions for better contour capture. Secondly, the proposed method
inconsistently localizes motion, as seen in Fig. 5(b). The global loss
function minimizes mean absolute error, potentially leading to incorrect
motion locations. Incorporating prior knowledge like optical flow could help,
but the feasibility of deep learning remains uncertain. Lastly, the evaluation
of registration and topology preservation warrants closer examination using
quantitative and qualitative measures. We focused on video segmentation, which
may not fully capture our method’s performance in these aspects. Future work
should involve comprehensive evaluation, the development of specific metrics,
and the incorporation of diverse datasets for robust effectiveness assessment.

(a) An example of an incorrect LV segmentation map with a non-smooth contour
due to noise.

(b) In cases of rapid cardiac motion, the neural ODE solution may not
accurately follow the heart’s movements.

Figure 5: Illustration of the two most common errors observed in the results
of our method.

##  7 Conclusion

In this study, we introduced the Echo-ODE method, leveraging Neural ODE to
enhance echocardiography segmentation. Our approach produces accurate and
temporally consistent segmentation maps without explicit temporal
optimization, outperforming state-of-the-art methods on the CAMUS dataset in
terms of smoothness, anatomical accuracy, and clinical metrics. We also
explored the impact of initial segmentation, finding minimal effects on
accuracy and clinical metrics. Despite limitations in handling complex
deformations and reliance on pre-trained models, the Echo-ODE approach shows
promise for broader medical imaging analysis. Future directions include
investigating advanced Neural ODE architectures, addressing limitations, and
exploring uncertainty estimation techniques for more reliable segmentation
maps. Additionally, exploring the method’s applicability to other
echocardiography processing and synthetic tasks could advance cardiac
research.

## References

  * [1] M. D. Cheitlin, J. S. Alpert, W. F. Armstrong _et al._ , “ACC/AHA guidelines for the clinical application of echocardiography: a report of the American College of Cardiology/American Heart Association Task Force on Practice Guidelines,” _Circulation_ , vol. 95, no. 6, pp. 1686–1744, 1997\.
  * [2] G. Varoquaux and V. Cheplygina, “Machine learning for medical imaging: methodological failures and recommendations for the future,” _NPJ Digital Medicine_ , vol. 5, no. 1, p. 48, 2022.
  * [3] S. Leclerc, E. Smistad, J. Pedrosa _et al._ , “Deep Learning for Segmentation Using an Open Large-Scale Dataset in 2D Echocardiography”,” _IEEE Trans Med Imaging_ , vol. 38, no. 9, pp. 2198–2210, 2019.
  * [4] R. Ge, G. Yang, Y. Chen _et al._ , “V-LVNet: Direct left ventricle multitype indices estimation from 2D echocardiograms of paired apical views with deep neural networks,” _Med Image Anal_ , vol. 58, p. 101554, 2019\.
  * [5] D. Ouyang, B. He, A. Ghorbani, M. P. Lungren, E. A. Ashley, D. H. Liang, and J. Y. Zou, “EchoNet-Dynamic: a Large New Cardiac Motion Video Data Resource for Medical Machine Learning,” in _NeurIPS ML4H Workshop: Vancouver, BC, Canada_ , 2019.
  * [6] O. A. Smiseth, H. Torp, A. Opdahl _et al._ , “Myocardial strain imaging: how useful is it in clinical decision making?” _European Heart Journal_ , vol. 37, no. 15, pp. 1196–1207, 2016.
  * [7] M. Cikes and S. D. Solomon, “Beyond ejection fraction: an integrative approach for assessment of cardiac structure and function in heart failure,” _European Heart Journal_ , vol. 37, no. 21, pp. 1642–1650, 2016.
  * [8] H. Wei, H. Cao, Y. Cao _et al._ , “Temporal-Consistent Segmentation of Echocardiography with Co-learning from Appearance and Shape,” in _Proc. MICCAI, LNCS_ , 2020, pp. 623–632.
  * [9] N. Painchaud, N. Duchateau, O. Bernard, and P.-M. Jodoin, “Echocardiography Segmentation With Enforced Temporal Consistency,” _IEEE Trans Med Imaging_ , vol. 41, no. 10, pp. 2867–2878, 2022.
  * [10] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, “Neural Ordinary Differential Equations,” in _Proc. NeurIPS_ , vol. 31, 2018.
  * [11] O. Ronneberger, P. Fischer, and T. Brox, “U-Net: Convolutional Networks for Biomedical Image Segmentation,” in _Proc. MICCAI, LNCS_ , 2015, pp. 234–241.
  * [12] L.-C. Chen, Y. Zhu, G. Papandreou _et al._ , “Encoder-Decoder with Atrous Separable Convolution for Semantic Image Segmentation,” in _Proc. ECCV_ , 2018, pp. 801–818.
  * [13] G. Veni, M. Moradi, H. Bulu _et al._ , “Echocardiography segmentation based on a shape-guided deformable model driven by a fully convolutional network prior,” in _Proc. ISBI_ , 2018, pp. 898–902.
  * [14] S. Thomas, A. Gilbert, and G. Ben-Yosef, “Light-weight Spatio-Temporal Graphs for Segmentation and Ejection Fraction Prediction in Cardiac Ultrasound,” in _Proc. MICCAI, LNCS_ , 2022, pp. 380–390.
  * [15] M. Li, C. Wang, H. Zhang, and G. Yang, “MV-RAN: Multiview recurrent aggregation network for echocardiographic sequences segmentation and full cardiac cycle analysis,” _Comput Biol and Med_ , vol. 120, p. 103728, 2020\.
  * [16] A. Pati and A. Lerch, “Attribute-based regularization of latent spaces for variational auto-encoders,” _Neural Computing and Applications_ , vol. 33, pp. 4429–4444, 2021.
  * [17] A. Sotiras, C. Davatzikos, and N. Paragios, “Deformable Medical Image Registration: A Survey,” _IEEE Trans Med Imaging_ , vol. 32, no. 7, pp. 1153–1190, 2013.
  * [18] S. M. Eick, E. A. Enright, S. D. Geiger, and _et al._ , “Associations of maternal stress, prenatal exposure to per-and polyfluoroalkyl substances (pfas), and demographic risk factors with birth outcomes and offspring neurodevelopment: an overview of the echo. ca. il prospective birth cohorts,” _International journal of environmental research and public health_ , vol. 18, no. 2, p. 742, 2021.
  * [19] G. Balakrishnan, A. Zhao, M. R. Sabuncu, J. Guttag, and A. V. Dalca, “VoxelMorph: A Learning Framework for Deformable Medical Image Registration,” _IEEE Trans Med Imaging_ , vol. 38, no. 8, pp. 1788–1800, 2019.
  * [20] W. Zhu, Y. Huang, D. Xu _et al._ , “Test-Time Training for Deformable Multi-Scale Image Registration,” in _Proc. ICRA_. IEEE, 2021, pp. 13 618–13 625.
  * [21] A. V. Dalca, G. Balakrishnan, J. Guttag, and M. R. Sabuncu, “Unsupervised Learning for Fast Probabilistic Diffeomorphic Registration,” in _Proc. MICCAI, LNCS_ , 2018, pp. 729–738.
  * [22] J. Xu, E. Z. Chen, X. Chen, T. Chen, and S. Sun, “Multi-scale Neural ODEs for 3D Medical Image Registration,” in _Proc. MICCAI, LNCS_ , 2021, pp. 213–223.
  * [23] A. Joshi and Y. Hong, “Diffeomorphic Image Registration using Lipschitz Continuous Residual Networks,” in _Proc. MIDL_ , 2021.
  * [24] F. Rousseau, L. Drumetz, and R. Fablet, “Residual Networks as Flows of Diffeomorphisms,” _Journal of Mathematical Imaging and Vision_ , vol. 62, no. 3, pp. 365–375, 2020.
  * [25] J. Ashburner, “A fast diffeomorphic image registration algorithm,” _Neuroimage_ , vol. 38, no. 1, pp. 95–113, 2007.
  * [26] A. Iserles, H. Z. Munthe-Kaas, S. P. Nørsett, and A. Zanna, “Lie-group methods,” _Acta Numerica_ , vol. 9, pp. 215–365, 2000.
  * [27] A. Vaswani, N. Shazeer, N. Parmar _et al._ , “Attention is All you Need,” in _Proc. NeurIPS_ , vol. 30. Curran Associates, Inc., 2017.
  * [28] T. Park, M.-Y. Liu, T.-C. Wang, and J.-Y. Zhu, “Semantic image synthesis with spatially-adaptive normalization,” in _Proc. CVPR_ , 2019.
  * [29] Y. Wu, T. Z. Jiahao, J. Wang _et al._ , “NODEO: A Neural Ordinary Differential Equation Based Optimization Framework for Deformable Image Registration,” in _Proc. CVPR_ , 2022, pp. 20 772–20 781.
  * [30] A. Paszke, A. Chaurasia, S. Kim _et al._ , “ENet: A Deep Neural Network Architecture for Real-Time Semantic Segmentation,” _arXiv preprint arXiv:1606.02147_ , 2016.