smmfstream#3dict/Type/Filespec/F (#3)/UF (#3)/EF ¡¡/F
¿¿dict/Type/RichMediaInstance/Subtype/Video/Asset /Params
¡¡/FlashVars(source=#3&skin=SkinOverAllNoFullNoCaption.swfskinAutoHide=trueskinBackgroundColor=0x5F5F5FskinBackgroundAlpha=0.75loop=True)¿¿dict/Type/RichMediaConfiguration/Subtype/Video/Instances
[]dict/Type/RichMediaContent/Assets ¡¡/Names [(#3) ]¿¿/Configurations
[]dict/Activation ¡¡/Condition/#1PVXA/Presentation
¡¡/Style/Embedded¿¿¿¿/Deactivation
¡¡/Condition/PI¿¿1113.33337pt6.94444pt1.94443pt/Subtype/RichMedia/BS ¡¡/W
0/S/S¿¿/Contents (embedded video file:#3)/NM (rma:#3)/AP ¡¡/N
¿¿/RichMediaSettings /RichMediaContent #2

# The Dirac-Delta Rogue Wave

Miguel C. N. Fiolhais1,2,3 and Andrea Ferroglia3,4
1 Science Department, Borough of Manhattan Community College,
The City University of New York, 199 Chambers St, New York, NY 10007, USA
2 LIP, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra,
Portugal
3 The Graduate School and University Center, The City University of New York,
365 Fifth Avenue, New York, NY 10016 USA
4 Physics Department, New York City College of Technology, The City University
of New York, 300 Jay Street, Brooklyn, NY 11201, USA

###### Abstract

This paper presents a double spatio-temporal localized Dirac-delta solution
for the linear wave equation. The solution arises from the interference of
sinusoidal waves with frequencies that vary as a function of the time of
emission. It is shown that the time-dependent frequency function required to
produce a localized Dirac-delta wave is exclusively determined by the
dispersion relation of the medium in which the wave is traveling. Both
numerical and exact analytical solutions for a typical physical scenario are
obtained and displayed. The application of this result to possible practical
uses is discussed.

##  I Introduction

A rogue wave, also known as a freak wave or monster wave, is an exceptionally
large ocean wave that occurs unexpectedly amid relatively smaller waves. Its
immense size, steepness, and unpredictability make it extremely dangerous for
ships, offshore structures and coastal areas. Rogue waves appear suddenly and
can reach heights of more than twice the significant wave height (the average
height of the highest third of waves in a given area.) Consequently, these
waves can have devastating impacts on vessels, capable of capsizing or
severely damaging even large ships PhysRevA.80.043818 ; PhysRevLett.114.213901
.

Once considered mythical or exaggerated accounts by sailors, the existence of
these waves has been confirmed through scientific observations and satellite
data christou ; benetazzo ; onoratoObservation . One possible explanation
stems from a solution of the non-linear Schrödinger equation - the Peregrine
soliton Peregrine1983 . The non-linear Schrödinger equation, is a fundamental
equation in quantum mechanics and nonlinear optics kaertner . It extends the
Schrödinger equation to describe the behavior of complex wave phenomena in
nonlinear systems and plays a crucial role in understanding the dynamics of
Bose-Einstein condensates, optical solitons, and other nonlinear wave
phenomena bronski ; ahmad . In the context of water waves, it accounts for
phenomena such as wave steepening, wave modulation, and wave interactions,
allowing for a more accurate representation of the complex dynamics of waves
in various conditions, such as those encountered in the ocean onorato . In
contrast to a typical soliton, which can retain its profile without alteration
while propagating, the Peregrine soliton exhibits a unique double spatio-
temporal localization. Consequently, when initiated with a modest oscillation
against a continuous background, the Peregrine soliton experiences a gradual
amplification of its intensity. Following this stage of maximum compression,
the soliton’s amplitude diminishes and its width expands. In other words,
large amplitudes that seemingly materialize out of thin air, vanish without a
trace just like a rogue wave. This phenomenon has been experimentally
reproduced in water tanks experiments PhysRevLett.106.204502 ;
PhysRevX.2.011015  and has also been observed in a variety of physical systems
ONORATO201347 . In the linear regime however, for example in the approximate
Airy wave theory for sea waves or in any other physical system that follows a
linear wave equation, rogue waves may result purely from the constructive
interference of smaller waves in dispersive media 2023PNAS..12006275H . In
certain conditions, waves traveling at different speeds may align and combine,
leading to the formation of a larger, more powerful wave for a short period of
time before disappearing onto the background.

The objective of this paper is to explore the possibility of generating such
effect with the emission of a sinusoidal wave with a varying frequency, such
as a chirp, in a dispersive medium. In order to do so, a general solution of
the classical linear wave equation is obtained for the extreme case of a
wavepacket that converges into a Dirac-delta function at a fixed location for
an arbitrary dispersive medium. The analytical solution is obtained by
considering the superposition of an infinite succession of infinitesimal
wavefronts generated at the source with a changing frequency. The same exact
solution is also derived by using the short-time Fourier transform of the wave
at the source, and by propagating each sinusoidal component of the transform
as a plane wave according to the dispersion relation of the medium. In
addition, it is shown that the time-dependent frequency function is determined
by the dispersion relation, making the production of a Dirac-delta rogue wave
possible in any dispersive medium, as long as the varying frequency of the
emitted sinusoidal wave is adjusted accordingly.

The results for the analytical solution are confirmed numerically for the
particular case of a dispersion relation where the angular frequency is
proportional to the square of the wavenumber. The numerical approach is
employed by considering an emission signal at the source comprising a
continuous succession of very short but finite sinusoidal signals, in order to
simulate a chirp. These plane waves are propagated into the medium according
to the dispersion relation, and their superposition is directly compared with
the analytical solution. The relevance and possible application of the result
presented in this paper is discussed in different contexts.

##  II Localized Dirac-Delta Wave

Consider a plane sinusoidal wave that is emitted along an $x$-axis with a
time-varying frequency. For simplicity, the source of the wave is made to
coincide with the origin of the $x$-axis. In addition, the amplitude at the
source at the instant $t_{s}$ is given the following value:

$F(x=0,t_{s})=\sin\varphi(t_{s})\,,$ (1)

where the instantaneous angular frequency at the source is given by the time
derivative of the angular phase,

$\omega(t_{s})=\frac{\partial\varphi(t_{s})}{\partial t_{s}}\,.$ (2)

In a dispersive medium, if the frequency at the source is varied over time,
successive fronts of the wave will be transmitted with a changing
instantaneous frequency at different phase velocities. If the variation of the
frequency is continuous, each wavefront emitted through an infinitesimal time
interval can be made to arrive simultaneously at a point $x=L$ at a time
instant $t_{0}$. By considering the wavefront emitted from the source at the
time instant $t_{s}$, the required velocity for it to arrive at $x=L$ can be
determined as,

$v=\frac{\omega}{k}=\frac{L}{t_{0}-t_{s}}\,.$ (3)

The convergence of phases is perhaps better observed in the diagrams shown in
Fig. 1, where the position and time coordinates are represented in the
vertical and horizontal axes, respectively. For convenience, dimensionless
variables are used for position and time, defined as:

$\widetilde{x}=\frac{x}{L}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\textrm{and}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\widetilde{t}=\frac{t}{t_{0}}.$
(4)

Figure 1: Dimensionless position-time diagrams for the emission of a wave
during the interval $\Delta\widetilde{t}$ and converging at
$\widetilde{x}=\widetilde{t}=1$. The left panel shows the path of a phase
emitted from the source at the instant $\widetilde{t_{s}}$ as a dashed line,
while the right panel shows the same phase emitted with a fixed frequency in
the finite interval between $\widetilde{t}_{s}-\Delta\widetilde{t}_{s}/2$ and
$\widetilde{t}_{s}+\Delta\widetilde{t}_{s}/2$.

If the dispersion relation of the medium, $\omega(k)$, is known, then Eq. (3)
fully determines the time dependency of frequency at the source
$\omega(\widetilde{t}_{s})$. Note that, for convenience, Eq. (3) can also be
expressed in terms of a dimensionless angular frequency,
$\widetilde{\omega}=\omega/\omega_{0}$, and wavenumber,
$\widetilde{k}=k/k_{0}$,

$\frac{\widetilde{\omega}}{\widetilde{k}}=\frac{k_{0}L}{\omega_{0}t_{0}}\frac{1}{1-\tilde{t}_{s}}\,,$
(5)

where $k_{0}$ is the wavenumber corresponding to the angular frequency
$\omega_{0}$ in a given dispersive medium. For simplicity, the special case
where the values of $\omega_{0}t_{0}$ and $k_{0}L$ are equal to one will be
considered for the remainder of this section.

As the time-dependent frequency at the origin is found directly from the
dispersion relation, the next step is to calculate the overall wave as a
result of the superposition of the infinite number of wavefronts emitted
during infinitesimal time intervals $\textrm{d}t_{s}$ between the initial time
instant and an arbitrary time interval $\Delta t$, with $\Delta t<t_{0}$.

The equation of a wavefront emitted from the origin at a time
$\widetilde{t}_{s}$ can be written as,

$f(\,\widetilde{x},\widetilde{t}\,)=\mathcal{U}\left(\,\widetilde{t}_{s}-\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,\right)\sin\left[\widetilde{\omega}(\widetilde{t}_{s})\widetilde{t}-\widetilde{k}(\widetilde{t}_{s})\widetilde{x}+\phi(\widetilde{t}_{s})\right]\,,$
(6)

for an arbitrary point $(\widetilde{t},\widetilde{x})$ in the left diagram of
Fig. 1, where

$\mathcal{U}(z)=\begin{cases}1&\text{if }z=0\,,\\\ 0&\text{if }z\neq
0\,.\end{cases}$ (7)

Note that the relation,

$\widetilde{t}_{s}=\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,,$ (8)

corresponds to the equation of the dashed line in the left graph of Fig. 1,
where $f$ takes a constant value along the line and is zero elsewhere. In
addition, as there is no superposition at the source, the phase
$\phi(\widetilde{t}_{s})$ is determined by the boundary condition,

$\displaystyle f(\,0,\widetilde{t}_{s}\,)$ $\displaystyle=$ $\displaystyle
F(\,0,\widetilde{t}_{s}\,)$ (9)
$\displaystyle\sin\left[\widetilde{\omega}(\widetilde{t}_{s})\widetilde{t}_{s}+\phi(\widetilde{t}_{s})\right]$
$\displaystyle=$ $\displaystyle\sin\varphi(\widetilde{t}_{s})\,,$ (10)

which leads to,

$\phi(\widetilde{t}_{s})=\varphi(\widetilde{t}_{s})-\widetilde{\omega}(\widetilde{t}_{s})\widetilde{t}_{s}\,.$
(11)

According to the superposition principle, the overall wave can then be
determined as the infinite sum of individual wavefronts emitted from the
source over the interval $\Delta\widetilde{t}$. To this effect, it is
convenient to perform this sum over small but finite intervals
$\Delta\widetilde{t}_{s}$, as represented on the right diagram of Fig. 1, and
then take the limit when these intervals become infinitesimal
($\Delta\widetilde{t}_{s}\to\textrm{d}\widetilde{t}_{s}$), so that the
interference is described by an integral. For instance, if a wave is being
emitted with a fixed angular frequency ${\omega}(\widetilde{t}_{s})$ in the
interval between $\widetilde{t}_{s}-\Delta\widetilde{t}_{s}/2$ and
$\widetilde{t}_{s}+\Delta\widetilde{t}_{s}/2$, the function $\mathcal{U}$ can
be expressed as,

$\mathcal{U}\left(\,\widetilde{t}_{s}-\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,\right)=\begin{cases}1&\text{if
}-\frac{\Delta\widetilde{t}_{s}}{2|1-\widetilde{x}|}<\widetilde{t}_{s}-\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}<+\frac{\Delta\widetilde{t}_{s}}{2|1-\widetilde{x}|}\,,\\\
0&\text{elsewhere,}\end{cases}$ (12)

where the condition expressed by the inequality in the equation above
corresponds to the region between the two dashed lines on the right panel of
Fig. 1. In the limit $\Delta\widetilde{t}_{s}\rightarrow 0$, these two
parallel lines merge and the function $\mathcal{U}$ can be written in terms of
the Dirac delta function as:

$\lim_{\Delta\widetilde{t}_{s}\to
0}\frac{1}{\Delta\widetilde{t}_{s}}\,\mathcal{U}\left(\,\widetilde{t}_{s}-\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,\right)=\frac{1}{|1-\widetilde{x}|}\delta\left(\,\widetilde{t}_{s}-\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,\right)\,.$
(13)

As a result, the final superposed wave in the continuous limit is,

$F(\,\widetilde{x},\widetilde{t}\,)=\frac{1}{|1-\widetilde{x}|}\int_{0}^{\Delta\widetilde{t}}\mathcal{\delta}\left(\,\widetilde{t}_{s}-\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,\right)\sin\left[\widetilde{\omega}(\widetilde{t}_{s})\widetilde{t}-\widetilde{k}(\widetilde{t}_{s})\widetilde{x}+\phi(\widetilde{t}_{s})\right]\textrm{d}\widetilde{t}_{s}\,.$
(14)

which in turn leads to the exact solution in a closed form:

$F(\,\widetilde{x},\widetilde{t}\,)=\frac{1}{|1-\widetilde{x}|}\sin\left[\widetilde{\omega}\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)\widetilde{t}-\widetilde{k}\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)\widetilde{x}+\phi\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)\right]\,.$
(15)

By using Eqs. (5) and (11), the superposition of waves in the continuous
limit, expressed in Eq. (14), can also be rewritten as,

$F(\,\widetilde{x},\widetilde{t}\,)=\frac{1}{|1-\widetilde{x}|}\int_{0}^{\Delta\widetilde{t}}\mathcal{\delta}\left(\,\widetilde{t}_{s}-\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,\right)\sin\left[\widetilde{\omega}(\widetilde{t}_{s})\left(\widetilde{t}-\widetilde{x}-(1-\widetilde{x})\widetilde{t}_{s}\right)+\varphi(\widetilde{t}_{s})\right]\textrm{d}\widetilde{t}_{s}\,,$
(16)

which upon integration leads to a more compact form,

$F(\,\widetilde{x},\widetilde{t}\,)=\frac{1}{|1-\widetilde{x}|}\sin{\varphi}\left(\,{\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}}\,\right)\,.$
(17)

It is clear from this expression that the amplitude diverges at
$\widetilde{x}=1$. It is important to stress that due to the superposition
principle, the integral of the function over time at any fixed location must
remain constant, _i.e._ independent of the position - every wavefront emitted
from the source over an infinitesimal time, must eventually cross any position
in front of the source. Therefore, as the wave is radiated from the source, it
converges into a Dirac delta function at a fixed location before dispersing
into infinity - a Dirac-delta rogue wave. Naturally, if the integral is
negative, the delta function is also negative at the divergence point, leading
to a narrow deep depression instead of a towering wave.

Indeed, in the Appendix, it is proven that $F$ behaves as a Dirac-Delta
function in both $\widetilde{x}$ and $\widetilde{t}$, so that for a smooth
test function $h$ one finds that

$\displaystyle\lim_{\widetilde{t}\to
1}\frac{1}{G_{x}(\,\Delta\widetilde{t}\,)}$
$\displaystyle\int_{\frac{\widetilde{t}-\Delta\widetilde{t}}{1-\Delta\widetilde{t}}}^{\widetilde{t}}\textrm{d}\widetilde{x}\,\,\,h(\,\widetilde{x}\,)F(\,\widetilde{x},\widetilde{t}\,)=\int_{-\infty}^{\infty}\textrm{d}\widetilde{x}\,\,\,h(\,\widetilde{x}\,)\delta(\,1-\widetilde{x}\,)=h(1)\,,$
$\displaystyle\lim_{\widetilde{x}\to
1}\frac{1}{G_{t}(\,\Delta\widetilde{t}\,)}$
$\displaystyle\int_{\widetilde{x}}^{\widetilde{x}+\Delta\widetilde{t}(1-\widetilde{x})}\\!\\!\\!\\!\textrm{d}\widetilde{t}\,\,\,h(\,\widetilde{t}\,)F(\,\widetilde{x},\widetilde{t}\,)=\int_{-\infty}^{\infty}\textrm{d}\widetilde{t}\,\,\,h(\,\widetilde{t}\,)\delta(\,1-\widetilde{t}\,)=h(1)\,,$
(18)

where $G_{x}$ and $G_{t}$ are suitable normalization factors whose analytic
expression is given in the Appendix.

##  III Short-Time Fourier Transform

The same analytical solution presented in Section II can also be derived using
a short-time Fourier transform (STFT) sejdic . In order to do so, the signal
at the source can be divided into shorter segments, and the Fourier transform
can be calculated separately on each segment. As such, in the limit in which
the short time intervals become infinitesimal, the signal at the source at a
given time $\tau$ can be expressed as,

$F(x=0,\tau)=\sin\varphi(\tau)=\int_{0}^{\Delta\widetilde{t}}\delta(\tau-\widetilde{t}_{s})\sin\varphi(\widetilde{t}_{s})\,\textrm{d}\widetilde{t}_{s}\,,$
(19)

where the Dirac delta $\delta(\tau-\widetilde{t}_{s})$ can be interpreted as
the window function centered around $\widetilde{t}_{s}$. The Fourier transform
of the integrand function,

$f(\,\tau,\widetilde{t}_{s}\,)=\delta(\tau-\widetilde{t}_{s})\sin\varphi(\widetilde{t}_{s})\,,$
(20)

with respect to $\tau$ is therefore,

$\mathcal{F}(\widetilde{\omega},\widetilde{t}_{s})=\frac{1}{\sqrt{2\pi}}\int_{0}^{\Delta\widetilde{t}}\delta(\tau^{\prime}-\widetilde{t}_{s})\sin\varphi(\widetilde{t}_{s})e^{-i\widetilde{\omega}\tau^{\prime}}\,\textrm{d}\tau^{\prime}=\frac{1}{\sqrt{2\pi}}\sin\varphi(\widetilde{t}_{s})e^{-i\widetilde{\omega}\widetilde{t}_{s}}\,,$
(21)

and the function can be recovered as,

$f(\,\tau,\widetilde{t}_{s}\,)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}\mathcal{F}(\widetilde{\omega},\widetilde{t}_{s})e^{i\widetilde{\omega}\tau}\textrm{d}\widetilde{\omega}\,.$
(22)

As each individual sinusoidal signal at a time $\tau$ at the source, that is
multiplied to the Fourier transform, must propagate as a plane wave according
to the dispersion relation the medium, the integrand function can be written
in terms of $\widetilde{x}$ and $\widetilde{t}$ as,

$f(\widetilde{x},\widetilde{t},\widetilde{t}_{s})=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}\mathcal{F}(\widetilde{\omega},\widetilde{t}_{s})e^{i\widetilde{\omega}\widetilde{t}-i\widetilde{k}\widetilde{x}}\textrm{d}\widetilde{\omega}\,,$
(23)

where $\widetilde{k}$ is the wave number associated to the frequency
$\widetilde{\omega}$. Consequently, the overall wave can also be calculated
for any pair of $\widetilde{x}$ and $\widetilde{t}$ as,

$F(\,\widetilde{x},\widetilde{t}\,)=\int_{0}^{\Delta\widetilde{t}}f(\widetilde{x},\widetilde{t},\widetilde{t}_{s})\,\textrm{d}\widetilde{t}_{s}=\frac{1}{{2\pi}}\int_{0}^{\Delta\widetilde{t}}\int_{-\infty}^{+\infty}\sin\varphi(\widetilde{t}_{s})e^{-i\widetilde{\omega}\widetilde{t}_{s}+i\widetilde{\omega}\widetilde{t}-i\widetilde{k}\widetilde{x}}\,\textrm{d}\widetilde{\omega}\,\textrm{d}\widetilde{t}_{s}$
(24)

This double integral can be rearranged as,

$F(\,\widetilde{x},\widetilde{t}\,)=\int_{0}^{\Delta\widetilde{t}}\left(\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{-i\widetilde{\omega}\widetilde{t}_{s}+i\widetilde{\omega}\widetilde{t}-i\widetilde{k}\widetilde{x}}\,\textrm{d}\widetilde{\omega}\right)\sin\varphi(\widetilde{t}_{s})\,\textrm{d}\widetilde{t}_{s}\,.$
(25)

By imposing the condition of Eq. (5), the final wave becomes,

$F(\,\widetilde{x},\widetilde{t}\,)=\frac{1}{|1-\widetilde{x}|}\int_{0}^{\Delta\widetilde{t}}\mathcal{\delta}\left(\,\widetilde{t}_{s}-\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,\right)\sin\varphi(\widetilde{t}_{s})\textrm{d}\widetilde{t}_{s}=\,\frac{1}{|1-\widetilde{x}|}\sin{\varphi}\left(\,{\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}}\,\right)\,,$
(26)

which is the exact same result obtained in Eq. (17).

##  IV Analytical vs Numerical Solution

For the physically common dispersion relation
$\widetilde{\omega}={\widetilde{k}}^{2}$, the dimensionless angular
frequencies at the source can be obtained using Eq. (5) as,

$\widetilde{\omega}(\widetilde{t}_{s})=\left(\frac{k_{0}L}{\omega_{0}t_{0}}\frac{1}{1-\widetilde{t}_{s}}\right)^{2}\,,$
(27)

which corresponds to the angular phase

${\varphi}(\widetilde{t}_{s})=\frac{k_{0}^{2}L^{2}}{\omega_{0}t_{0}}\frac{\widetilde{t}_{s}}{1-\widetilde{t}_{s}}\,.$
(28)

where the integration constant is fixed in such a way that $\varphi(0)=0$.
This result can be extended to a more general dispersion relation,
$\widetilde{\omega}=\widetilde{k}^{\alpha}$ ($\alpha\neq 1$), in which the
angular phase becomes:

$\varphi(\widetilde{t}_{s})=\omega_{0}t_{0}(\alpha-1)\left(\frac{k_{0}L}{\omega_{0}t_{0}}\right)^{\frac{\alpha}{\alpha-1}}\left[(1-\widetilde{t}_{s})^{\frac{1}{1-\alpha}}-1\right]\,.$
(29)

It must be stressed that the method fails for $\alpha=1$, which corresponds to
the case of a non-dispersive medium.

Assuming once more the special case where the values of $\omega_{0}t_{0}$ and
$k_{0}L$ are equal to one, the analytical solution can be presented in a
compact form as,

$\widetilde{\omega}=\widetilde{k}^{2}\Rightarrow
F(\,\widetilde{x},\widetilde{t}\,)=\frac{1}{|1-\widetilde{x}|}\sin\left(\,\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{t}}\,\right)\,.$
(30)

This result can be visualized in Figure 2, which shows the time profile of the
wave as it evolves through the spatial coordinate, compressing and increasing
its amplitude. The value of $\Delta\widetilde{t}$ is set at approximately
0.9263 for visual effect purposes, _i.e._ so the signal includes two full
cycles, starting and ending at zero.

Figure 2: Time profile of $F(\,\widetilde{x},\widetilde{t}\,)$ for the
dispersion relation $\widetilde{\omega}=\widetilde{k}^{2}$, as it evolves
through the spatial coordinate, compressing and increasing its amplitude, from
the source at $\widetilde{x}=0$ (blue) through $\widetilde{x}=1/3$ (orange)
and $\widetilde{x}=2/3$ (green), until the diverging point at
$\widetilde{x}=1$ (red), _i.e._ the Dirac-delta function. The value of
$\Delta\widetilde{t}$ is set at approximately 0.9263.

For any point $\widetilde{x}$ between $0$ and $1$, it can be seen that the
time dependence of the amplitude has the exact same form as the time profile
at $\widetilde{x}=0$, within the allowed time interval. After all, this is the
only possible way that such point could operate as a source and lead to the
same localized wave at $\widetilde{x}=1$. In other words, the amplitude at the
source is being mapped to any position $\widetilde{x}$ with the following
transform:

$\widetilde{t}_{s}\rightarrow\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,,$
(31)

shifting the function to the right by an amount $\widetilde{x}$ and
compressing it by a factor $1-\widetilde{x}$. As the integral remains constant
due to the superposition principle, the amplitude also scales with the same
factor.

At this point, it is convenient to revisit the progression of the wave in
space as shown in Fig. 2. For the sake of clarity, one can assume the wave in
the example is mechanical and the amplitude corresponds to a physical
displacement. At the source, _i.e._ at the fixed location $\widetilde{x}=0$,
an object riding the wave will complete two cycles in the dimensionless time
interval $\Delta\widetilde{t}$. However, at any other fixed location, for
example, $\widetilde{x}=1/3$ or $\widetilde{x}=2/3$, an object riding the wave
will complete the very same cycles in a shorter and shorter time,
respectively. In the extreme case of $\widetilde{x}=1$, the cycles would be
completed instantaneously. This leads to an apparent contradiction. When a
sinusoidal wave with a fixed frequency travels into a dispersive medium, the
frequency must be preserved while the wavenumber changes according to the
dispersion relation. Nonetheless, in this example, the frequency appears to
increase as the wave progresses into the medium. This is not the case,
however, as each component of the STFT maintains its own frequency while
traveling through the dispersive medium. The apparent increase in the
frequency of the object riding the wave is no more than a collective effect of
the superposition of these components.

###  IV.1 Numerical solution

Despite the fact that this analytical solution is supported by two different
formal derivations, it is critical to test this result with an independent
brute force numerical algorithm for confirmation. The procedure for the
numerical approach is presented below and compared with the analytical
solution of the physical case considered before.

The mathematical function described by Eq. (1) represents the source of a
sinusoidal plane wave with a varying frequency. In the discrete limit, this
can be approximated as a series of sinusoidal segments, where each segment has
a constant frequency,

$F_{n}(\,\widetilde{x}=0,\widetilde{t}_{s}\,)=\sin\left[\widetilde{\omega}(\widetilde{t}_{n})\widetilde{t}_{s}+\phi_{n}\right]\,,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\widetilde{t}_{n}-\frac{\Delta\widetilde{t}_{s}}{2}<\widetilde{t}_{s}<\widetilde{t}_{n}+\frac{\Delta\widetilde{t}_{s}}{2}\,,$
(32)

where $\Delta\widetilde{t}_{s}$ is the dimensionless time interval of each
segment and,

$\widetilde{t}_{n}=\left(n-\frac{1}{2}\right)\Delta\widetilde{t}_{s}\,,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,n=1\,,\,2\,,\,...\,,N$
(33)

corresponds to the time instant at the center of the $n$th segment, with
$N=\Delta\widetilde{t}/\Delta\widetilde{t}_{s}$, where
$\Delta\widetilde{t}_{s}$ is chosen in such a way that $N$ is an integer. The
dimensionless angular frequency is determined by Eq. (5) and by the dispersion
relation specific to medium. These relations ensure that each segment arrives
at $\widetilde{x}=1$ simultaneously at $\widetilde{t}=1$. In order to maintain
continuity, the phase of each sinusoidal segment can be determined recursively
as,

$\widetilde{\omega}(\widetilde{t}_{n})n\Delta\widetilde{t}_{s}+\phi_{n}=\widetilde{\omega}(\widetilde{t}_{n+1})n\Delta\widetilde{t}_{s}+\phi_{n+1}\,.$
(34)

As an example, the agreement between the analytical continuum limit of the
varying frequency wave at the source with the discrete sum of sinusoidal
segments can be observed in Fig. 3 for the dispersion relation
$\widetilde{\omega}=\widetilde{k}^{2}$ and $\Delta\widetilde{t}_{s}=0.025$.
The blue dots represent the beginning and end of each segment.

Figure 3: Comparison between the analytical function at the source,
$F(\,\widetilde{x}=0,\widetilde{t}_{s}\,)$, represented by the orange line,
with the discrete sum of sinusoidal segments shown in blue. The blue dots
represent the beginning and end of each segment for
$\Delta\widetilde{t}_{s}=0.025$.

As each time interval corresponds to a sinusoidal plane wave with a fixed
frequency being radiated from the source, the overall wave in the numerical
approximation can be obtained as the superposition of each segment moving
through space with their corresponding constant velocity,

$F(\,\widetilde{x},\widetilde{t}\,)=\sum_{n=1}^{N}F_{n}(\,\widetilde{x},\widetilde{t}\,)=\sum_{n=1}^{N}\mathcal{U}_{n}\sin\left[\widetilde{\omega}(\widetilde{t}_{n})\widetilde{t}-\widetilde{k}(\widetilde{t}_{n})\widetilde{x}+\phi_{n}\right]\,,$
(35)

where,

$\displaystyle\mathcal{U}_{n}=\begin{cases}1&\text{if
}\widetilde{x}\left[1-\left(n-\frac{1}{2}\right)\Delta\widetilde{t}_{s}\right]+\left(n-1\right)\Delta\widetilde{t}_{s}<\widetilde{t}<\widetilde{x}\left[1-\left(n-\frac{1}{2}\right)\Delta\widetilde{t}_{s}\right]+n\Delta\widetilde{t}_{s}\,,\\\
0&\text{otherwise.}\end{cases}$ (36)

It should be noted that in this numerical approach, as the signal is divided
into shorter segments of equal length, and approximated as a sinusoidal
function in each segment, $\mathcal{U}_{n}$ can be interpreted as a
rectangular window function. As the segments progress in space and time as
plane waves, their overlap will naturally lead to discontinuities at the edges
of the segments as the overall wave starts compressing. In order to address
this issue, a Constant Overlap-Add (COLA) window function, such as the Hann
function smith3 , is employed in the numerical calculation leading to Fig. 4,
in order to ensure the overall wave’s smoothness. It is worth stressing that
as the time interval of the segments approaches zero, the choice of the window
function becomes inconsequential.

The comparison between the numerical and the analytical solution is presented
in Fig. 4 for $\Delta\widetilde{t}_{s}=10^{-4}$. The comparison is shown for
different time instants right before and after the formation of the Dirac-
delta peak, presenting a nearly perfect agreement. The dimensionless time
instants in the panels range from 0.96 to 1.04, corresponding to the interval
where the wave compresses above tenfold the initial amplitude and diverges at
1.00. The complete time progression of the wave can also be seen in the
following video video , again showing the excellent agreement between the
analytical and numerical results. In both representations, it is visible that
the wave presents symmetric shapes before and after $\widetilde{t}=1$. For
example, the shape at $\widetilde{t}=0.96$ is mirrored with respect to
$\widetilde{x}=1$ when compared with the shape at $\widetilde{t}=1.04$.

Figure 4: Comparison between the numerical and the analytical solutions for
dimensionless time instants ranging from 0.96 to 1.04. The figure refers to
the case of the dispersion relation $\widetilde{\omega}=\widetilde{k}^{2}$.

##  V Discussion and Conclusions

Even though the previous numerical analysis focused on the case of a
particular dispersion relation, the same approach can be applied to any
physical scenario where dispersion occurs. In fact, the angular phase required
to produce a localized Dirac-delta wave was obtained for a general power
function dispersion relation in Eq. (29). This is the case, for example, of
water waves in the deep water regime where the dispersion relation is
approximately $\widetilde{\omega}={\widetilde{k}}^{1/2}$, according to the
Airy wave theory. As the frequency of water waves depends primarily on the
wind speed, it is plausible that a varying gust in the right conditions may
lead to the formation of a rogue wave in a similar fashion as the localized
Dirac-delta wave.

Creating a focused wave targeting a specific region of space holds immense
potential in a number of applications. In the realm of cancer treatment for
example, the ability to direct waves precisely to the affected area allows for
more effective and less invasive therapies. As several types of waves
experience dispersion in biological tissues schwan1 ; schwan2 ; guedert ;
sugimoto , a targeted approach like the one proposed in this paper could
potentially minimize collateral damage to surrounding healthy tissues in the
medical treatment of tumors, enhancing the precision of medical interventions.
Another appealing physical scenario is the case of ion-acoustic waves in
fusion plasmas. Ion acoustic waves are a distinctive form of plasma
oscillations that propagate through ionized gases, or plasmas, found in a
variety of astrophysical and laboratory settings gekelman . These waves arise
from the collective motion of ions in response to perturbations, generating
compression disturbances in the plasma. Dominated by the inertia of ions and
the restoring force due to thermal pressure, ion acoustic waves play a crucial
role in fusion experiments. While the dispersion of ion-acoustic waves is more
complex than the case analyzed in this paper fitzpatrick ; liu , the linear
nature of such waves allows for the possibility of producing a localized
Dirac-delta wave with an appropriate choice for the time-dependent angular
phase at the source. The production of such wave could possibly lead to a
localized increase of density above the threshold for nuclear fusion reactions
to occur, enabling a controlled release of energy, leaving the rest of the
plasma mostly unaffected.

A further compelling consequence of localized Dirac-delta waves is the ability
to produce an impulse signal. It is widely known that a chirp shares the same
power spectrum of an impulse signal smith . Therefore, when an impulse signal
is transmitted into a dispersive medium, it may undergo a conversion into a
chirp. This, in fact, is the principle behind Chirped Pulse Amplification
(CPA) strickland , where an ultra-short laser pulse is temporally extended
(chirped) to reduce its peak power using diffraction gratings. Subsequently,
it is amplified without causing damage to the laser medium. Following
amplification, the pulse is compressed back to its original duration,
resulting in an intense ultrashort laser pulse. This technique ultimately led
to the Nobel Prize in Physics in 2018, awarded to Strickland and Mourou.
Analogously, a chirp can also be propagated into a dispersive medium in a
manner that momentarily converges the signal into a Dirac-delta function
before dispersing into the background, much like a rogue wave. As such, if a
chirp signal is sent from a non-dispersive medium into a dispersive slab with
the exact thickness to converge into a Dirac-delta at the end of the slab, the
signal would be transmitted as a Dirac-delta pulse back into the non-
dispersive medium, as exemplified in Fig. 5. Therefore, the production of the
localized Dirac-delta wave as detailed in this paper can also be interpreted
as the time reversal ($t\rightarrow-t$) of a Dirac-delta pulse emitted into a
dispersive medium. As a result, it shall be possible to use this time symmetry
to find out experimentally which shape a signal should have at the source in
order to produce a Delta-function pulse at the desirable distance from the
source. This may be particularly useful in media with complicated dispersion
relations where the analytical calculation is challenging.

Figure 5: Progression of a chirped sinusoidal wave from a non-dispersive
medium into a dispersive slab. The thickness is such that the wave compresses
into a Dirac-delta at the end of the slab, and exits back into the non-
dispersive medium as a Dirac-delta pulse wave.

In summary, this study presents a novel solution to the linear wave equation
in dispersive media: a localized Dirac-delta wave. The proposed solution
arises from the interference of a sinusoidal wave characterized by a time-
dependent frequency, uniquely dictated by the dispersion relation of the
medium. Both numerical and exact analytical solutions are presented for a
particular physical scenario, showing a nearly perfect agreement between them.
While a thorough exploration of the practical applications stemming from this
result is not within the scope of this paper, the authors are optimistic that
the proposal may inspire and contribute to potential real-world
implementations.

## Acknowledgments

The authors would like to thank Dr. Giovanni Ossola for carefully reading the
manuscript and giving valuable suggestions.

## Appendix: Dirac Delta function properties of the function $F$

In this subsection it is shown that the integral of $F$ at fixed
$\widetilde{t}$ does not depend on $\widetilde{t}$, as long as
$\Delta\widetilde{t}<\widetilde{t}<1$. The integral of interest is

$G_{x}=\int_{\frac{\widetilde{t}-\Delta\widetilde{t}}{1-\Delta\widetilde{t}}}^{\widetilde{t}}\,\frac{1}{1-\widetilde{x}}\sin{\varphi\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)}\,\textrm{d}\widetilde{x}\,.$
(37)

It is now sufficient to change integration variable according to

$\widetilde{t}_{s}=\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,,\qquad\widetilde{x}=\frac{\widetilde{t}-\widetilde{t}_{s}}{1-\widetilde{t}_{s}}\,,\qquad\textrm{d}\widetilde{x}=-\frac{1-\widetilde{t}}{(1-\widetilde{t}_{s})^{2}}\textrm{d}\widetilde{t}_{s}\,,$
(38)

to obtain

$G_{x}(\,\Delta\widetilde{t}\,)=\int_{0}^{\Delta\widetilde{t}}\frac{1}{1-\widetilde{t}_{s}}\sin{\varphi(\widetilde{t}_{s})}\,\textrm{d}\widetilde{t}_{s}\,,$
(39)

where it was made explicit that $G$ is a function of the parameter
$\Delta\widetilde{t}$. It is therefore obvious that

$\frac{1}{G_{x}(\,\Delta\widetilde{t}\,)}\int_{\frac{\widetilde{t}-\Delta\widetilde{t}}{1-\Delta\widetilde{t}}}^{\widetilde{t}}\frac{1}{1-\widetilde{x}}\sin{\varphi\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)}\,\textrm{d}\widetilde{x}=1\,.$
(40)

With the same change of variables discussed above, it is possible to prove
that for a sufficiently smooth function $h(\widetilde{x})$

$I=\lim_{\widetilde{t}\to
1}\frac{1}{G_{x}(\,\Delta\widetilde{t}\,)}\int_{\frac{\widetilde{t}-\Delta\widetilde{t}}{1-\Delta\widetilde{t}}}^{\widetilde{t}}\,h(\widetilde{x})\frac{1}{1-\widetilde{x}}\sin{\varphi\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)}\,\textrm{d}\widetilde{x}=h(1)\,.$
(41)

Indeed, after the change in integration variables
$\widetilde{x}\to\widetilde{t}_{s}$ introduced above, one finds

$I=\lim_{\widetilde{t}\to
1}\frac{1}{G_{x}(\,\Delta\widetilde{t}\,)}\int_{0}^{\Delta\widetilde{t}}h\left(\frac{\widetilde{t}-\widetilde{t}_{s}}{1-\widetilde{t}_{s}}\right)\frac{1}{1-\widetilde{t}_{s}}\sin{\varphi(\widetilde{t}_{s})}\,\textrm{d}\widetilde{t}_{s}\,.$
(42)

At this stage, it is possible to take the limit of the integrand before
integration, to find

$I=\frac{1}{G_{x}(\,\Delta\widetilde{t}\,)}\int_{0}^{\Delta\widetilde{t}}\,h\left(1\right)\frac{1}{1-\widetilde{t}_{s}}\sin{\left(\varphi(\widetilde{t}_{s})\right)}\,\textrm{d}\widetilde{t}_{s}=h(1).$
(43)

It is therefore possible to conclude that

$\lim_{\widetilde{t}\to
1}\frac{1}{G_{x}(\,\Delta\widetilde{t}\,)}\int_{\frac{\widetilde{t}-\Delta\widetilde{t}}{1-\Delta\widetilde{t}}}^{\widetilde{t}}h(\widetilde{x})\frac{1}{1-\widetilde{x}}\sin{\varphi\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)}\,\textrm{d}\widetilde{x}=\int_{-\infty}^{\infty}h(\widetilde{x})\delta(1-\widetilde{x})\,\textrm{d}\widetilde{x}=h(1)\,.$
(44)

Similar considerations apply to the integrals of the function $F$ with respect
to $\widetilde{t}$. For instance, one can also show that the integral of $F$
at fixed $\widetilde{x}$ does not depend on $\widetilde{t}$. In this case, the
integral of interest is

$G_{t}=\int_{\widetilde{x}}^{\widetilde{x}+\Delta\widetilde{t}(1-\widetilde{x})}\frac{1}{1-\widetilde{x}}\sin{\varphi\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)}\,\textrm{d}\widetilde{t}\,.$
(45)

It is now sufficient to change integration variable according to

$\widetilde{t}_{s}=\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\,,\qquad\widetilde{t}=\widetilde{x}+\widetilde{t}_{s}(1-\widetilde{x})\,,\qquad\textrm{d}\widetilde{t}=(1-\widetilde{x})\textrm{d}\widetilde{t}_{s}\,,$
(46)

to obtain

$G_{t}(\,\Delta\widetilde{t}\,)=\int_{0}^{\Delta\widetilde{t}}\sin{\varphi(\widetilde{t}_{s})}\,\textrm{d}\widetilde{t}_{s}\,,$
(47)

where it was made explicit that $G$ is a function of the parameter
$\Delta\widetilde{t}$. It is therefore obvious that

$\frac{1}{G_{t}(\,\Delta\widetilde{t}\,)}\int_{\widetilde{x}}^{\widetilde{x}+\Delta\widetilde{t}(1-\widetilde{x})}\frac{1}{1-\widetilde{x}}\sin{\varphi\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)}\,\textrm{d}\widetilde{t}=1\,.$
(48)

With the same change of variables discussed above, it is possible to prove
that for a sufficiently smooth function $h(\widetilde{x})$

$I=\lim_{\widetilde{x}\to
1}\frac{1}{G_{t}(\,\Delta\widetilde{t}\,)}\int_{\widetilde{x}}^{\widetilde{x}+\Delta\widetilde{t}(1-\widetilde{x})}h(\,\widetilde{t}\,)\frac{1}{1-\widetilde{x}}\sin{\varphi\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)}\,\textrm{d}\widetilde{t}=h(1)\,.$
(49)

Indeed, after the change in integration variables
$\widetilde{x}\to\widetilde{t}_{s}$ introduced above, one finds

$I=\lim_{\widetilde{x}\to
1}\frac{1}{G_{t}(\,\Delta\widetilde{t}\,)}\int_{0}^{\Delta\widetilde{t}}h\left(\widetilde{x}+\widetilde{t}_{s}(1-\widetilde{x})\right)\sin{\varphi(\widetilde{t}_{s})}\,\textrm{d}\widetilde{t}_{s}\,.$
(50)

It is again possible to take the limit of the integrand before integration, to
find

$I=\frac{1}{G_{t}(\,\Delta\widetilde{t}\,)}\int_{0}^{\Delta\widetilde{t}}h\left(1\right)\sin{\varphi(\widetilde{t}_{s})}\,\textrm{d}\widetilde{t}_{s}=h(1).$
(51)

Consequently, one finds that

$\lim_{\widetilde{x}\to
1}\frac{1}{G_{t}(\,\Delta\widetilde{t}\,)}\int_{\widetilde{x}}^{\widetilde{x}+\Delta\widetilde{t}(1-\widetilde{x})}h(\,\widetilde{t}\,)\frac{1}{1-\widetilde{x}}\sin{\varphi\left(\frac{\widetilde{t}-\widetilde{x}}{1-\widetilde{x}}\right)}\,\textrm{d}\widetilde{t}=\int_{-\infty}^{\infty}h(\,\widetilde{t}\,)\delta(\,1-\widetilde{t}\,)\,\textrm{d}\widetilde{t}=h(1)\,.$
(52)

## References

  * (1) N. Akhmediev, J. M. Soto-Crespo and A. Ankiewicz, “How to excite a rogue wave,” Phys. Rev. A 80, 043818 (2009).
  * (2) S. Birkholz, C. Brée, A. Demircan, and G. Steinmeyer, “Predictability of Rogue Events,” Phys. Rev. Lett. 114, 213901 (2015).
  * (3) M. Christou and K. Ewans, “Field Measurements of Rogue Water Waves,” Journal of Physical Oceanography 44, 2317–2335 (2014).
  * (4) A. Benetazzo  _et al._ , “On the shape and likelihood of oceanic rogue waves,” Sci. Rep. 7, 8276 (2017).
  * (5) M. Onorato  _et al._ , “Observation of a giant nonlinear wave-packet on the surface of the ocean,” Sci. Rep. 11, 23606 (2021).
  * (6) D. H. Peregrine, “Water waves, nonlinear Schrödinger equations and their solutions,” J. Austral. Math. Soc. B 25, 16-43 (1983).
  * (7) F. X. Kaertner, _Ultrafast Optics_ (LibreTexts, Massachusetts Institute of Technology, 2023). See Chapter 3.3: The Nonlinear Schrödinger Equation.
  * (8) J. C. Bronski, L. D. Carr, B. Deconinck, and J. N. Kutz, “Bose-Einstein Condensates in Standing Waves: The Cubic Nonlinear Schrödinger Equation with a Periodic Potential,” Phys. Rev. Lett. 86, 1402 (2001).
  * (9) J. Ahmad  _et al._ , “Soliton solutions of fractional extended nonlinear Schrödinger equation arising in plasma physics and nonlinear optical fiber,” Sci Rep 13, 10877 (2023).
  * (10) M. Onorato, D. Proment, G. Clauss, M. Klein, “Rogue waves: from nonlinear Schrödinger breather solutions to sea-keeping test,” PLoS One 8, e54629 (2013).
  * (11) A. Chabchoub, N. P. Hoffmann. N. Akhmediev, “Rogue Wave Observation in a Water Wave Tank,” Phys. Rev. Lett. 106, 204502 (2011).
  * (12) A. Chabchoub, N. P. Hoffmann, M. Onorato, N. Akhmediev, “Super Rogue Waves: Observation of a Higher-Order Breather in Water Waves,” Phys. Rev. X 2, 011015 (2012).
  * (13) M. Onorato, S. Residori, U. Bortolozzo, A. Montina and F. T. Arecchi, “Rogue waves and their generating mechanisms in different physical contexts,” Physics Reports 528, 47-89 (2013).
  * (14) D. Häfner, J. Gemmrich and M. Jochum, “Machine-guided discovery of a real-world rogue wave model,” Proc. Natl. Acad. Sci. 120, e2306275120 (2023).
  * (15) E. Sejdić, I. Djurović and J. Jiang, “Time–frequency feature representation using energy concentration: An overview of recent advances,” Digital Signal Processing 19, 153-183 (2009).
  * (16) J. O. Smith III, _Spectral Audio Signal Processing_ (W3K Publishing, Stanford, 2011). See Chapter 8: Overlap-Add STFT Processing.
  * (17) https://youtu.be/R8hH_csY-zw
  * (18) H. P. Schwan, “Electrical properties of tissue and cell suspensions,” Advances in Biological and Medical Physics 5, 147–209 (1957).
  * (19) H. P. Schwan, “Electrical properties of tissues and cell suspensions: mechanisms and models,” Proceedings of 16th Annual International Conference of the IEEE Engineering in Medicine and Biology Society 1, A70-A71 (1994).
  * (20) R. Guedert  _et al._ , “Biological dispersion in the time domain using finite element method software,” Sci. Rep. 13, 22868 (2023).
  * (21) K. Sugimoto  _et al._ , “Clinical utilization of shear wave dispersion imaging in diffuse liver disease,” Ultrasonography 39, 3–10 (2020).
  * (22) W. Gekelman, “Ion acoustic wave experiments in a high school plasma physics laboratory,” American Journal of Physics 75, 103–110 (2007).
  * (23) R. Fitzpatrick, _Plasma Physics: An Introduction_ (CRC Press, Taylor & Francis Group, 2022). See subsection Ion Acoustic Waves in Chapter 7: Waves in Warm Plasmas.
  * (24) L. Liu and J. Du, “Ion acoustic waves in the plasma with the power-law q-distribution in nonextensive statistics,” Physica A 387, 4821-4827 (2008).
  * (25) S. W. Smith, _The Scientist and Engineer’s Guide to Digital Signal Processing_ (California Technical Publishing, San Diego, 1997). See subsection Chirp Signals in Chapter 11: Fourier Transform Pairs.
  * (26) D. Strickland and G. Mourou, “Compression of amplified chirped optical pulses,” Optics Communications 56, 219-221 (1985).


# Color image segmentation based on a convex K-means approach

Tingting Wu1, Xiaoyu Gu1, Jinbo Shao1, Ruoxuan Zhou2, Zhi Li3,∗ ∗The
corresponding author of this paper.1School of Science, Nanjing University of
Posts and Telecommunications, Nanjing, 210023, China.2College of Oceanography
and Space Informatics, China University of Petroleum, Qingdao, 266580,
China.3The Department of Computer Science and Technology, Shanghai Key
Laboratory of Multidimensional Information Processing, East China Normal
University, Shanghai, 200241, China. E-mail<EMAIL_ADDRESS>
###### Abstract

Image segmentation is a fundamental and challenging task in image processing
and computer vision. The color image segmentation is attracting more attention
due to the color image provides more information than the gray image. In this
paper, we propose a variational model based on a convex K-means approach to
segment color images. The proposed variational method uses a combination of
$l_{1}$ and $l_{2}$ regularizers to maintain edge information of objects in
images while overcoming the staircase effect. Meanwhile, our one-stage
strategy is an improved version based on the smoothing and thresholding
strategy, which contributes to improving the accuracy of segmentation. The
proposed method performs the following steps. First, we specify the color set
which can be determined by human or the K-means method. Second, we use a
variational model to obtain the most appropriate color for each pixel from the
color set via convex relaxation and lifting. The Chambolle-Pock algorithm and
simplex projection are applied to solve the variational model effectively.
Experimental results and comparison analysis demonstrate the effectiveness and
robustness of our method.

##  I Introduction

(a) Given image (b) SLaT [1] (c) Proposed

Figure 1: Segmentation results of Flowers with Gaussian noise (mean
$\tilde{\mu}=0$ and variance $\sigma^{2}=0.1$). Method [1] used classical
K-means framework (Figure (b)), where the pixels around the stalks are wrongly
classified; our proposed method in Figure (c) achieves a relatively better
classification, see the rectangle part.

(a) Given image (b) LC [2] (c) Proposed

Figure 2: Segmentation results of Parrot. Figure (b) produces some unwanted
tiny structures [2]; our result in Figure (c) can eliminate the unwanted tiny
structures near the boundary by introducing the square of the first-order
derivative term.

Image segmentation is a fundamental and significant problem in image
processing and computer vision, which has been widely applied in many fields
such as object detection and medical diagnosis [3, 4, 5, 6]. Image
segmentation aims at dividing an image of $N$ pixels into $K$ regions with
similar characteristics together (edges, intensities, colors or textures).
Various models and algorithms have been extensively applied for image
segmentation, including level-set methods [7, 8, 9, 10, 11, 12, 13], active
contours [14, 15, 16, 17], variational models [18, 19, 20, 21] and clustering
methods [22, 23, 24], etc. Particularly, many successful methods for image
segmentation are based on variational models and clustering methods.
Variational models are flexible in structure, diverse in form, and superior in
performance. The clustering methods are suitable and useful for both low- and
high-dimensional data.

The Mumford-Shah (MS) model [18] is one of the most popular and influential
variational segmentation models. The MS model proposed an energy minimization
problem which approximates the true solution by finding optimal piecewise
smooth approximations. However, it is very challenging to find or approximate
its minimizer due to the MS model is nonconvex. It is worth noting that Cai et
al. proposed a two-stage strategy, i.e., smoothing and thresholding (SaT), for
image segmentation [25]. In the first stage, a convex optimization problem
based on the MS model is solved to get a smooth approximation $u$ of the given
image $f$. This stage could be regarded as the process of smoothing. Then in
the second stage, they used the K-means algorithm or other clustering
algorithms to obtain thresholds for segmentation, which can be seen as the
process of thresholding. The segmentation stage is independent of the
optimization problem in the first stage, one can change the number of phases
$K$ or thresholds $\rho$ without resolving the optimization problem.

Since the SaT strategy is a flexible methodology for image segmentation, it
has been adopted and followed by many papers [26, 19, 20, 27, 28, 29, 30, 31,
21]. And the strategy was extended to images corrupted by Poisson or Gamma
noises [32], color images [1] and classification for high-dimensional data and
point clouds. For example, literature [1] proposed a three-stage method, i.e.,
Smoothing, Lifting and Thresholding, achieving a good result for segmenting
images corrupted by noise. In [21], Cai et al. revealed the connection between
image restoration and segmentation by using the SaT strategy.

In the above segmentation methods based on the SaT strategy, thresholding is
recognized as an effective method of image segmentation. A great number of
thresholding techniques have been widely used for grayscale and color image
segmentation. Furthermore, some clustering algorithms have been used by many
image segmentation models to determine thresholds for segmentation, such as
the fuzzy clustering [23, 24, 33, 5, 34] and the K-means clustering [1, 25,
26, 32, 35, 36], etc. In [23], a thresholding method based on fuzzy set was
proposed to deal with medical images, and the method got good performance
visually and numerically. Jia et al. [34] presented a robust self-sparse fuzzy
clustering algorithm for solving over-segmentation.

The K-means method is popular due to its simple principle and high calculation
rate. In 2017, Condat proposed two methods for image segmentation based on a
convex approach to the K-means method. The paper [2] proposed a definition of
discrete total variation (TV) and applied it to the segmentation of color
images. Meanwhile, the method [37] showed a convex formulation of data
clustering and image segmentation based on the “Potts” model (see Section II-A
for details), convex relaxation and lifting. However, the K centroids are pre-
assigned by the author before the segmentation takes place. The above methods
[37, 2] introduced a large number of candidates which did not improve the
image quality enough to be worth the extra computing cost, although avoided
dealing with a non-convex problem. Furthermore, models in [37, 2] are both
based on TV regularization, which could cause the staircase effect, especially
for noisy images.

For the SaT strategy, the segmentation is carried in the thresholding stage
and it uses the smoothed image obtained in the smoothing stage. The SaT
strategy is quite important, and it can keep main features of the image while
clear out the noise and meaningless features so that in all cases it can
obtain good segmentation results. Besides, the K-means method has been widely
used in segmentation method to determine the thresholds for the thresholding
stage in the SaT strategy. However, Figure 1(b) is the result of SLaT [1]
which used classical K-means framework to perform the segmentation. We can see
that the pixels around the stalks are wrongly classified, there are many “red
dots” around it (see the marked area).

Thus the goals of this paper are summarized as follows:

  * a.

to avoid solving non-convex and combinatorial problem for the sake of
stability;

  * b.

to reduce the staircase effect in the model based on the TV regularization;

  * c.

to inherit the advantages of the SaT strategy especially the smoothing stage;

  * d.

to adopt a new strategy for segmentation to determine the most appropriate
color for each pixel.

In summary, we propose a new method to achieve the above goals. The proposed
method combines convex relaxation based on the K-means method and lifting. We
perform an improved convex relaxation for the classical K-means method and
transform the problem into an assignment problem by the technique of lifting
[1, 38, 37, 2, 39]. On the one hand, we add a square of the first-order
derivative term based on the variational method [37] to reduce the staircase
effect. The segmentation performance can be well guaranteed by the combination
of $l_{1}$ and $l_{2}$ regularizers. On the other hand, we use the strategy of
lifting to determine the most appropriate color for each pixel from the color
set, which could overcome the drawback of the classical K-means method in
color image segmentation.

The organization of this paper is as follows. In Section 2, we give the basic
knowledge used throughout the paper. In Section 3, we propose our model for
image segmentation and present the Chambolle-Pock algorithm to solve the
suggested model. In Section 4, we present numerical experiments and some
comparisons to show the superiority of our one-stage method. Conclusions are
drawn in Section 5.

(a) Three-phases (b) Palette $M=4$ (c) Six-phases

(d) Palette $M=6$

Figure 3: The original synthetic images and corresponding Palettes with Three-
phases and Six-phases, respectively.

(a) Airplane (b) Hill (c) Flowers (d) Parrot

(e) Ladybug  (f) Sunflowers

Figure 4: The original real-world images with different color sets and
structures.

(a) Palette $M=2$

(b) Palette $M=3$

(c) Palette $M=4$

(d) Palette $M=6$

(e) Palette $M=6$

(f) Palette $M=6$

Figure 5: The used palettes in Figure 4(a)-(f).

##  II Model

###  II-A Potts model

Firstly, we give some definitions for the classical K-means method in image
processing. The data $f=(f_{n})_{n\in\Omega}$ is a 2D image with values
$f_{n}$, where domain
$\Omega=\\{1,\cdot\cdot\cdot,N_{1}\\}\times\\{1,\cdot\cdot\cdot,N_{2}\\}$, and
$n=N_{1}\times N_{2}$. The classical K-means problem minimizes the sum of
Euclidean distance from every points $f_{n}$ to the nearest centroid $c_{k}$,
aiming at partitioning $N$ points of $\mathcal{R}^{d}$ into $K$ groups, i.e.,
it solves the following problem to partition points into different clusters

$\underset{(\Omega_{k})_{k=1}^{K},(c_{k})_{k=1}^{K}}{\operatorname{minimize}}\frac{1}{2}\sum_{k=1}^{K}\sum_{n\in\Omega_{k}}\left\|f_{n}-c_{k}\right\|^{2}.$
(1)

Secondly, we consider a general problem based on the classical K-means
problem. Given an integer $K\geq 2$, one wants to partition $\Omega$ into $K$
regions $\Omega_{k}$ (so $\bigcup^{K}_{k=1}\Omega_{k}=\Omega$ and
$\Omega_{k^{{}^{\prime}}}\bigcap\Omega_{k}=\emptyset$, for all $k\neq
k^{{}^{\prime}}$), and to find the corresponding centroids
$c_{k}\in\mathcal{R}^{d}$, so as to

$\underset{(\Omega_{k})_{k=1}^{K},(c_{k})_{k=1}^{K}}{\operatorname{minimize}}\frac{1}{2}\sum_{k=1}^{K}\sum_{n\in\Omega_{k}}\left\|f_{n}-c_{k}\right\|^{2}+\frac{\lambda}{2}\sum_{k=1}^{K}\operatorname{per}\left(\Omega_{k}\right).$
(2)

In the above model (2), the first term is the data fidelity term, the second
term is the penalization of the sum of the region perimeter (the region
perimeter is $\operatorname{per}\left(\Omega_{k}\right)$, and
$\operatorname{per}$ denotes the perimeter). The parameter $\lambda$ is to
control the level of spatial regularization, which has been adopted in the NP-
hard piecewise constant Mumford-Shah problem [40, 18] to obtain spatial
homogeneity by penalizing the sum of the region perimeter. The model is named
as “Potts” (Penalization of the region perimeter) model. The classical K-means
problem is the special case of the equation (2), then we can partition the
point cloud $f=(f_{n})_{n\in\Omega}$ whose $c_{k}$ are the means. Please see
[37, 2] for more details.

###  II-B Hill-climbing method

Traditional hill-climbing segmentation [41, 42] is a simple, fast algorithm
that clusters colors of an image without any hand-tuning of parameters. Li et
al. made a few modifications to get a more stable hill-climbing procedure
[20]. Moreover, the most significant advantage of hill-climbing is that it
does not need prior knowledge of the number of clusters or the content of the
given image and can detect the number of segments. As papers [43, 20]
mentioned, the number of peaks is equal to the number of clusters, i.e., the
number of segments ($K$). This process is quite simple to implement: we just
group the matrix $f$ into a 3D histogram, then the hill-climbing algorithm
attempts to locate the peaks on the color histogram [20]. For more details,
see [43, 41, 20, 42] and references therein.

###  II-C The smoothed Potts model

As discussed in the above subsection, we first use the hill-climbing procedure
to detect the segments ($K$) in advance and then we use the K-means algorithm
to get the set $\Sigma={c_{k}\in[0,1]^{3}}$, where $k=1,\ldots,K$ of $K\geq 2$
colors $c_{k}$, denoting the triplets of R, G, B values. We can use the Matlab
K-means function “kmeans” to get the $c_{k}$. Here, it is worth mentioning
that there is no doubt that the color set could be specified by human or the
K-means method. We can choose the most suitable color for each pixel by the
segmentation accuracy to get a finer result.

We can define the segmented RGB color image $u$, whose R, G, B values at every
pixel is one of the $c_{k}$, close to $f$ (also known as piecewise constant MS
problem).

In the segmented image $u$, every element $u_{n}$ is one of the centroids
$c_{k}$, thus we can use the lifting strategy [1, 38, 37, 2, 39] to
reformulate the problem (2) as an equivalent assignment problem, then the
unknown variable is the assignment array $z$ indexed by $k=1,\ldots,K$ with
one more dimension than $f$. This reformulation meets the goal (d): to adopt a
new strategy for segmentation to determine the most appropriate color for each
pixel.

The assignment array $z$ has the following properties:

  * a.

The assignment array $z$ belongs to the set $\mathcal{A}$ of binary assignment
vectors.

  * b.

For every $n\in\Omega$ and $k=1,\ldots,K$, if $u_{n}=c_{k}$ then
$z_{k,n}\in\mathcal{A}^{\Omega}$ is equal to 1 and to 0 else, i.e., we want to
choose a $c_{k}$ from the color set to represent the pixel of the
$n\in\Omega$.

  * c.

Each vector $z_{:,n}=\left(z_{k,n}\right)_{k=1}^{K}$ with elements in ${0,1}$
and its sum is 1.

Based on the above properties, we can get $u$ by a simple summation

$u_{n}=\sum_{k=1}^{K}z_{k,n}c_{k},\quad\forall n\in\Omega.$ (3)

Then, in short $\langle z,w\rangle=\frac{1}{2}\left\|f-u\right\|^{2}$, where
$w_{k,n}=\frac{1}{2}\left\|f_{n}-c_{k}\right\|^{2}$.

The TV of the indicator function of set (1 inside, 0 outside) is equal to the
perimeter of that set [38]. We can rewrite the second term in (2) as

$\displaystyle\frac{\lambda}{2}\sum_{k=1}^{K}\operatorname{per}\left(\Omega_{k}\right)=\frac{\lambda}{2}\|\nabla
z\|_{1}.$ (4)

Here we have some discrete form of the TV, and $\|\nabla z\|_{1}$ denotes the
discrete total variation of $z$, i.e., $\|\nabla
z\|_{1}=\sum_{k}\sqrt{\left(\nabla_{x}z\right)_{k}^{2}+\left(\nabla_{y}z\right)_{k}^{2}}$
with $\nabla_{x}$ and $\nabla_{y}$ corresponding to the discrete derivative
operators in the $x$-direction and $y$-direction, respectively.

As we know, image processing models based on TV regularization may preserve
important edge information, but the staircase effect may be introduced. To
reduce the staircase effect [25], we add a square of the first-order
derivative in (4). In short, we denote

$\lambda\sum_{k=1}^{K}\operatorname{per}\left(\Omega_{k}\right):=\lambda\|\nabla
z\|_{1}+\frac{\mu}{2}\|\nabla z\|_{2}^{2}.$ (5)

The effect of the square of the first-order derivative term can be seen in
Figure 2.

As we can see, when we set $\mu=0$, model 5 becomes the model of [2].
Comparing Figure 2(b) with (c), we can find the (b) has some unwanted tiny
structures (see the marked area). We add the square of the first-order
derivative term, which is a smooth term. This term works and smoothes out the
individual dots, see Figure 2(c). More precisely, Figure 2(b) produces
unwanted tiny structures as the method of [2] based on TV regularization;
Figure 2(c) can eliminate the unwanted tiny structures near the boundary by
introducing the square of the first-order derivative term. So, it is of great
significance to add this smooth term. This modification meets the goal (b) and
(c): to reduce the staircase effect and inherit the advantages of the
smoothing stage.

(a) $\tilde{\mu}=0$, $\sigma^{2}=0.1$ (b) $\tilde{\mu}=0$, $\sigma^{2}=0.3$
(c) $\tilde{\mu}=0$, $\sigma^{2}=0.5$ (d) $\tilde{\mu}=0$, $\sigma^{2}=0.1$

Figure 6: The degraded images for Figure 3(a) and (c). For Figure (a)-(c), we
add the Gaussian noise with mean $\tilde{\mu}=0$ and variance
$\sigma^{2}=0.1$, $\sigma^{2}=0.3$, and $\sigma^{2}=0.5$, respectively. For
Figure (d), we add the Gaussian noise with mean $\tilde{\mu}=0$ and variance
$\sigma^{2}=0.1$.

(a) FRC [44] (b) DPP [45] (c) SLaT [1] (d) LC [2]

(e) RSSFC [34] (f) Proposed

Figure 7: The segmentation results of Three-shapes in Figure 6(a). We add the
Gaussian noise with mean $\tilde{\mu}=0$, variance $\sigma^{2}=0.1$.

(a) FRC [44] (b) DPP [45] (c) SLaT [1] (d) LC [2]

(e) RSSFC [34] (f) Proposed

Figure 8: The segmentation results of Three-shapes in Figure 6(b). We add the
Gaussian noise with mean $\tilde{\mu}=0$, variance $\sigma^{2}=0.3$.

(a) FRC [44] (b) DPP [45] (c) SLaT [1] (d) LC [2]

(e) RSSFC [34] (f) Proposed

Figure 9: The segmentation results of Three-shapes in Figure 6(c). We add the Gaussian noise with mean $\tilde{\mu}=0$, variance $\sigma^{2}=0.5$. TABLE I: The runtime in seconds of the color set by the K-means method and human specified (including the time taken in deciding the number of colors). Images | K-means | Human
---|---|---
Figure 6(a) | 8.36 | 157.67
Figure 6(b) | 8.36 | 159.54
Figure 6(c) | 8.39 | 159.93
Figure 6(d) | 6.18 | 275.27

(a) FRC [44] (b) DPP [45] (c) SLaT [1] (d) LC [2]

(e) RSSFC [34] (f) Proposed

Figure 10: The segmentation results of Six-phases in Figure 6(d). In Figure
(a) and (b), the FRC and the DPP fail for the case of the image with Gaussian
noise; the SLaT (Figure (c)) and RSSFC (Figure (e)) misclassify the color of
the area between the two neighbor circles; our proposed (Figure (f)) is better
than the LC (Figure (d)) as it not only tells the five different colors but
also preserves the shape of the circle. TABLE II: The comparisons of SA and
CPU-time in seconds by different methods in Figure 4.
(The best results in the simulation are given in bold.) Images | Methods | SA | Time(s) | Images | Methods | SA | Time(s) | Images | Methods | SA | Time(s)
---|---|---|---|---|---|---|---|---|---|---|---
Airplane | FRC | 0.8674 | 3.33 | Hill | FRC | 0.7690 | 2.95 | Flowers | FRC | 0.7965 | 4.33
DPP | 0.6807 | 4.89 | DPP | 0.8646 | 4.80 | DPP | 0.5986 | 5.91
SLaT | 0.9067 | 7.82 | SLaT | 0.8891 | 9.77 | SLaT | 0.8997 | 5.58
LC | 0.9095 | 72.11 | LC | 0.9065 | 73.04 | LC | 0.9038 | 90.06
CKC | 0.9101 | 3546.64 | CKC | 0.9099 | 3094.92 | CKC | 0.9048 | 3915.16
CQaS | 0.9059 | 2109.49 | CQaS | 0.9056 | 2098.49 | CQaS | 0.9104 | 2589.49
RSSFC | 0.8754 | 44.61 | RSSFC | 0.7805 | 54.99 | RSSFC | 0.8505 | 71.93
Proposed | 0.9170 | 14.34 | Proposed | 0.9109 | 16.05 | Proposed | 0.9151 | 16.55
Parrot | FRC | 0.8032 | 4.47 | Ladybug | FRC | 0.8245 | 2.00 | Sunflowers | FRC | 0.8216 | 1.51
DPP | 0.5846 | 6.10 | DPP | 0.6940 | 1.55 | DPP | 0.6712 | 1.48
SLaT | 0.7539 | 6.08 | SLaT | 0.7984 | 1.51 | SLaT | 0.8568 | 1.30
LC | 0.8376 | 79.22 | LC | 0.8582 | 25.36 | LC | 0.8769 | 22.23
CKC | 0.7827 | 2040.88 | CKC | 0.8389 | 2576.99 | CKC | 0.8798 | 2191.62
CQaS | 0.8269 | 1135.65 | CQaS | 0.8513 | 1424.72 | CQaS | 0.8809 | 1264.82
RSSFC | 0.6904 | 86.41 | RSSFC | 0.7225 | 20.27 | RSSFC | 0.6387 | 16.86
Proposed | 0.8389 | 15.52 | Proposed | 0.8624 | 7.55 | Proposed | 0.8824 | 7.05

##  III The proposed model

Now, we propose our model as follows

$\min E(z)=\underset{z\in\mathcal{A}^{\Omega}}{\min}\left\\{\lambda\|\nabla
z\|_{1}+\frac{\mu}{2}\|\nabla z\|_{2}^{2}+\langle z,w\rangle\right\\},$ (6)

where $\langle z,w\rangle$ is defined in Section II-C. The set $\mathcal{A}$
is non-convex, so the problem (6) is non-convex. Here we perform a convex
relaxation, replacing $\mathcal{A}$ with its convex hull (the simplex
$\Delta$), which is the set of vectors with nonnegative elements whose sum is
1. This strategy meets the goal (a): to avoid solving non-convex and
combinatorial problem for the sake of stability, then we reformulate the
problem (6) as follows

$\min_{z}\left\\{\lambda\|\nabla z\|_{1}+\frac{\mu}{2}\|\nabla
z\|_{2}^{2}+\langle z,w\rangle+\imath_{\Delta}(z)\right\\}.$ (7)

Since model (7) is convex, there are many methods to solve it such as the
split-Bregman algorithm [46], which is used widely to solve many problems
based on TV regularization and the alternating direction method of multipliers
(ADMM) [47, 48, 49]. Specially, the Chambolle-Pock algorithm [50] proves a
convergence rate.

In the following, we use the Chambolle-Pock algorithm to solve (7), first we
introduce new variable $v=\nabla z$, then we consider the following primal-
dual optimization problem

$\begin{split}\min_{z,v}\max_{q}\bigg{\\{}&\lambda\|v\|_{1}+\frac{\mu}{2}\|v\|_{2}^{2}+\langle
z,w\rangle+\\\ &\langle v-\nabla
z,q\rangle+\imath_{\Delta}(z)\bigg{\\}}.\end{split}$ (8)

Then the Chambolle-Pock algorithm is defined through the following iterations

$\displaystyle
q^{(k+1)}=\arg\max_{q}\bigg{\\{}\left\langle\overline{v}^{(k)}-\nabla\overline{z}^{(k)},q\right\rangle$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{1}{2\sigma}\left\|q-q^{(k)}\right\|_{2}^{2}\bigg{\\}},$
(9) $\displaystyle z^{(k+1)}=\arg\min_{z}\bigg{\\{}-\left\langle\nabla
z,q^{(k+1)}\right\rangle+\langle z,w\rangle+\imath_{\Delta}(z)$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{1}{2\tau}\left\|z-z^{(k)}\right\|_{2}^{2}\bigg{\\}},$
(10) $\displaystyle
v^{(k+1)}=\arg\min_{v}\bigg{\\{}\lambda\|v\|_{1}+\frac{\mu}{2}\|v\|_{2}^{2}+\left\langle
v,q^{(k+1)}\right\rangle$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{1}{2\tau}\left\|v-v^{(k)}\right\|_{2}^{2}\bigg{\\}},$
(11) $\displaystyle\overline{z}^{(k+1)}=2z^{(k+1)}-z^{(k)},$ (12)
$\displaystyle\overline{v}^{(k+1)}=2v^{(k+1)}-v^{(k)}.$ (13)

Since the objective functions (9) is quadratic, the update of $q$ can be
computed efficiently

$q^{(k+1)}=\sigma\left(\overline{v}^{(k)}-\nabla\overline{z}^{(k)}\right)+q^{(k)}.$
(14)

The solution of (11) can be easily obtained by applying the soft thresholding
operator. Denoting
$t^{(k)}=\frac{\tau}{\mu\tau+1}\left(\frac{1}{\tau}v^{(k)}-q^{(k+1)}\right)$,
we have

$v_{x}^{(k+1)}=\max\left\\{\left|t^{(k)}\right|-\frac{\lambda\tau}{\mu\tau+1},0\right\\}\cdot\frac{t_{x}^{(k)}}{\left|t^{(k)}\right|},$
(15)
$v_{y}^{(k+1)}=\max\left\\{\left|t^{(k)}\right|-\frac{\lambda\tau}{\mu\tau+1},0\right\\}\cdot\frac{t_{y}^{(k)}}{\left|t^{(k)}\right|}.$
(16)

Denoting $z_{0}^{k}=z^{k}-\tau(w+\operatorname{div}q^{(k+1)})$, we have

$z^{(k+1)}=\arg\min_{z}\left\\{\frac{1}{2\tau}\left\|z-z_{0}^{k}\right\|_{2}^{2}+\imath_{\Delta}(z)\right\\},$
(17)

then

$z^{(k+1)}=P_{\Delta}(z_{0}^{k}),$ (18)

where $P_{\Delta}$ is the simplex projection, which could be computed
efficiently by the code from https://lcondat.github.io/.

Our one-stage method has the following three steps:

  * 1.

First, we use the hill-climbing procedure to detect the number of segments
($K$) and use the K-means algorithm to get $K$ colors $c_{k}$.

  * 2.

Second, we solve our one-stage model (merging the process of smoothing and
segmentation) via the Chambolle-Pock algorithm to get the assignment array
$z$.

  * 3.

Third, we get the segmented image $u$ via simple summation.

Then the whole process of our one-stage strategy image segmentation method for
clean or noisy images can be summarized by the following Algorithm 1.

Algorithm 1: Chambolle-Pock algorithm to solve (7)
---
1\. Initialize: given the $z^{0},\overline{z}^{0}$, $c_{k}$, $q^{0}=0$,
$v^{0}=\overline{v}^{0}=0$.
2\. Do k=0, 1,…, maxiteation:
a. compute $q^{k+1}$ by(14);
b. compute $v^{k+1}$ by (15) and (16);
c. compute $z^{k+1}$ by (18);
d. update $\overline{z}^{k+1}$ and $\overline{v}^{k+1}$ by (12) and (13).
3\. Output: the segmented $u$ by (3).

(a) FRC [44] (b) DPP [45] (c) SLaT [1] (d) LC [2]

(e) RSSFC [34] (f) Proposed

Figure 11: The segmentation results of Airplane in Figure 4(a). The FRC
(Figure (a)) and the RSSFC (Figure (e)) do not clearly cut the airplane out as
there is a small area that is not separated into the background; the DPP
(Figure (b)) has a poor performance on this task; The SLaT (Figure (c)), the
LC (Figure (d)) and our proposed (Figure (f), $\lambda=0.2$) cut the plane
from the sky (background), and segment the image into two reasonable parts.

(a) FRC [44] (b) DPP [45] (c) SLaT [1] (d) LC [2]

(e) RSSFC [34] (f) Proposed

Figure 12: The segmentation results of Hill in Figure 4(b). The FRC (Figure
(a)) and RSSFC (Figure(e)) fail to separate the hills as an entirety; the DPP
(Figure(b)) is somewhat under-segmented that some small features have
vanished; the LC (Figure (d)) does not cut the hill clearly out; the SLaT
(Figure (c)) and our proposed (Figure (f), $\lambda=0.2$) get relatively
better results.

(a) FRC [44] (b) DPP [45] (c) SLaT [1] (d) LC [2]

(e) RSSFC [34] (f) Proposed

Figure 13: The segmentation results of Flowers in Figure 4(c) with Gaussian
noise ($\tilde{\mu}=0$, $\sigma^{2}=0.1$). The FRC (Figure (a)) fails to
distinguish the colors; the DPP (Figure (b)) does not make out the objects in
the flowers; the SLaT (Figure (c)) misclassifies the area around the
receptacle of the first and third flowers; the LC (Figure (d)) produces some
yellow “dots”; the RSSFC (Figure (e)) misclassifies the second flowers; our
proposed (Figure (f), $\lambda=0.4$) gets a relatively better result. Figure
14: The 279 colors candidates in [37] obtained by sampling the CIELAB space.

(a) SLaT [1] (b) LC [2] (c) CKC [37] (d) CQaS [37]

(e) RSSFC [34] (f) Proposed

Figure 15: The segmentation results of Parrot in Figure 4(d). We test the
K-means palette to specify the color set shown in Figure 5(d). The SLaT in
Figure (a) and CKC in Figure (c) undercut the image; the LC in Figure (b) gets
a more smooth result than the former methods; the CQaS in Figure (d) and RSSFC
in Figure (e) have a poor performance obviously; our proposed gets a
relatively good result at the object (parrot) in Figure (f) ($\lambda=0.09$).

(a) SLaT [1] (b) LC [2] (c) CKC [37] (d) CQaS [37]

(e) RSSFC [34] (f) Proposed

Figure 16: Segmentation results of Ladybug in Figure 4(e). The SLaT (Figure
(a)), LC (Figure (b)), CKC (Figure (c)) and our method (Figure (f),
$\lambda=0.05$) all have a good result; the RSSFC (Figure (e)) has some under-
segmentation to the objects in the image; our method is less complicated than
the CKC and CQaS, and achieves a comparable good result ($\lambda=0.02$).

(a) SLaT [1] (b) LC [2] (c) CKC [37] (d) CQaS [37]

(e) RSSFC [34] (f) Proposed

Figure 17: Segmentation results of Sunflowers in Figure 4(f). The candidates
of palette in this example is shown in Figure 5(f). The CQaS (Figure (d))
overcuts the sky; the RSSFC (Figure(e)) undercuts the sunflowers; the SLaT
(Figure (a)), LC (Figure (b)), CKC (Figure (c)) and our method (Figure (f))
all have a good result; our method is less complicated than the CKC and
achieves a comparable good result ($\lambda=0.09$).

##  IV Experiments

In this section, to show the superior performance of our proposed method, we
compare our method with some methods. The names of these models are the Fuzzy
Region Competition (FRC) model [44], the Dynamic Programming for the Potts
(DPP) [45], the Smoothing, Lifting and Thresholding (SLaT) model [1], the
Potts model for Segmentation (LC) [2], the Convex K-means Clustering (CKC)
model [37] and the Robust Self-Sparse Fuzzy Clustering (RSSFC) model [34]. The
FRC used fuzzy membership functions to approximate the piecewise constant
Mumford-Shah model. The DPP used the dynamic programming and the ADMM
algorithm to solve the Potts model which is non-convex. The SLaT uses a three-
stage approach for Segmenting color images (smoothing, lifting and
thresholding). The LC proposed a convex approach to K-means without limitation
for the segmentation of Parrot. The CKC presented a lifting strategy to solve
a piecewise constant Mumford-Shah problem. The RSSFC presented a method of
robust self-sparse fuzzy clustering for image segmentation. For all methods in
experiments, the scope of the main parameter is provided as follows:

  * a.

FRC [44]: the regularization parameter $\lambda\in$ [0.0001, 1];

  * b.

DPP [45]: the regularization parameter $\lambda\in$ [0.1, 5];

  * c.

SLaT [1]: the regularization parameter $\lambda\in$ [1, 10];

  * d.

LC [2]: the regularization parameter $\lambda\in$ [0.001, 10];

  * e.

CKC [37]: the regularization parameter $\lambda\in$ [1, 500];

  * f.

RSSFC [34]: the regularization parameter $\lambda\in$ [0.01, 1];

  * g.

Proposed: the regularization parameter $\lambda\in$ [0.0001, 1].

All the experiments were performed under Matlab R2020a and Windows 10(x64) on
a PC with an Intel Core (TM) i7 10750H CPU at 2.60GHz and 16.0GB of memory.

To quantitative comparison of the segmentation performance, we use one
objective indicator, i.e., the segmentation accuracy (SA) [34]. The SA is
computed as follows:

$SA=\sum_{i=1}^{k}\dfrac{S_{i}\cap G_{i}}{n},$ (19)

where the $S_{i}$ represents a segmentation result, the $G_{i}$ denotes the
corresponding Ground Truth, the $k$ is the number of clusters and $n$ is the
total number of pixels of images.

###  IV-A Segmentation of synthetic images

In this subsection, we show the effectiveness of our method on some simple
segmentation tasks. These images have obvious different colors and have a
reasonable criterion of the segmentation, e.g., we cut the objects from the
background. Meanwhile, we compare the robustness of our method with the SLaT
on different noise levels. We use the Matlab “imnoise” command to add Gaussian
noise.

As aforementioned, the color set could be determined by the K-means method or
human. In examples of Three-phases and Six-phases, we use the above two
methods to specify the color set. The runtime comparisons show in Table I
between the color set specified by the K-means method and human (including the
time taken in deciding the number of colors). We can learn that the runtime of
the color set specified by K-means method faster than human, and two methods
can arrow the same segmentation accuracy. Thus, we use the K-means method to
choose the color set in the following examples.

Example 1: Three-phases segmentation. Figure 3(a) is the three-phase synthetic
image consisting of three types of shapes. In Figure 7, 8 and 9, we add the
Gaussian noise with $\tilde{\mu}=0$, $\sigma^{2}=0.1$, 0.3 and 0.5
($\tilde{\mu}$: mean, $\sigma^{2}$: variance) respectively. With the variance
of noise increases, the performance of the FRC, SLaT and RSSFC becomes poorer
and poorer. We can see that the SLaT has misclassified colors. What’s more,
our method gets more smooth and clear results than the LC and DPP. We can
observe that our method is robust for different levels of noise.

Example 2: Six-phases segmentation. Figure 10 shows the segmentation results
of the six-phase synthetic image consisted of five overlapping circles with
different colors. It is natural for us to cut the five overlapping circles
with different colors from the background. From the results, we can observe
that the FRC, DPP and RSSFC fail for the case of the image with Gaussian
noise. The FRC is unable to tell the five colors, and the DPP can not keep the
shape of the black and green circle and the RSSFC fails to maintain the edges
of all circles. Although the SLaT gets a relatively good result, it
misclassifies the color of the area between the two neighbor circles. Our
method is better than the LC as it not only tells the five different colors
(see Figure 3(d)) but also preserves the shape of the circle.

###  IV-B Segmentation of real-world images

In this subsection, we test our method on six real-world images to show that
our method can get a comparable segmentation result to the CKC without complex
calculation. The choice of six real-world images is according to the
comparison methods [44, 45, 1, 2, 37, 34].

Example 3: Airplane segmentation. We consider a common and simple task of
segmenting an image into two parts: foreground (object) and background
regions. We use the airplane image to show the superior performance of our
method on this task. The SLaT, LC and our method cut the plane from the sky
(background) and segment the image into two reasonable parts. Although this is
a simple example, the DPP Figure 11(b) has a poor performance on this task.
The FRC and RSSFC do not clearly cut the airplane out as there is a small area
that is not separated into the background.

Example 4: Hill segmentation. Figure 12 shows the segmentation results of the
hill. The FRC and RSSFC fail to separate the hills as an entirety, the DPP is
somewhat under-segmented that some small features have vanished, and the LC
does not cut the hill clearly out. The SLaT and our method get relatively
better results.

Example 5: Flowers segmentation. Figure 4(c) consists of four flowers with
red, yellow and green color, see Figure 5(c). There is no doubt that we want
to cut the flowers with different colors from the sky (background). As we can
see from Figure 13, for the image with Gaussian noise, the FRC fails to
distinguish the colors. The DPP does not make out the objects in the flowers.
The SLaT misclassifies the area around the receptacle of the first and third
flowers, see the red “dots”. The LC produces some yellow “dots”. The RSSFC
misclassifies the second flower. Our method gets a relatively good result.

In the above examples, we find that the FRC, FRC and RSSFC perform badly, the
SLaT and LC are much superior to the former methods. In the following
examples, we mainly compare our method with the SLaT, LC, RSSFC and CKC for
three color images shown in [37] including the Parrot, the Ladybug and the
Sunflowers.

Here, we emphasize that the method [37] specified the $M$ (279 in [37])
candidates which takes its values in only $K$ among the $M$ candidates ($K\ll
M$). The $M$ colors are shown in Figure 14. Meanwhile, the method [37] has
shown that it is far superior to a two-step strategy: first estimates the $K$
colors using quantization and then solves the segmentation problem restricted
to these $M=K$ colors. For convenience, we call this as Coupling Quantization
and Segmentation (CQaS) method [37]. Although [37] proposed an accelerated
algorithm and the segmentation results are good, the computational time is
still long, as shown in Table II.

In the next subsections, we use the examples used in the [37] to show that our
simple modification could get a comparable result to the CKC without complex
calculation. For the sake of simplicity, we set $K=6$ in the following
examples.

Example 6: Parrot segmentation. In Figure 15, the CQaS and RSSFC have a poor
performance obviously. The LC gets a more smooth result than the former
method, but it yields a different color. The SLaT undercuts the image, and the
area of the eye has been removed. Although our method yields a wrong area in
the background, it gets a relatively good result at the object (parrot).
What’s more, our method with a suitable color set gets a comparable good
result to the CKC.

Example 7: Ladybug segmentation. In this example, Figure 5(e) shows the
palette used. Figure 16 shows the segmentation results of the ladybug. The
SLaT, LC, CKC, CQaS and our method all have a good result. Here, we emphasize
that our method is less complicated than the CKC and achieves a comparable
good result to the CKC.

Example 8: Sunflowers segmentation. The candidate of the palette in this
example is shown in Figure 5(f). Sunflowers have some flowers and the green
area at the down of the image, we just want to detect the flowers and take the
green area as an entirety. The CQaS overcuts the sky for the given image, and
the RSSFC undercuts the sky for the given image. The SLaT, LC, CKC and our
method all have a good result.

Here, we highlight that our method is less complicated than the CKC and
achieves a comparable good result to the CKC, as shown in Table II.

##  V Conclusion

In this paper, we propose a new variational segmentation model based on an
improved K-means method. We first give a color set determined by human or the
K-means method. Then we use a variational model to select the most appropriate
colors by SA for each pixel from the color set to get a finer result. The
experiments show the effectiveness and robustness of our proposed method.
Compared with some state-of-the-art methods, our method has a well running
speed while ensuring higher segmentation accuracy. Recently, deep learning has
attracted wide attention. In future work, we intend to use data-driven deep
learning to train the regularization and combine the variational model for
image segmentation.

## Acknowledgment

The research was supported by East China Normal University through startup
funding grant number 15901-120215-10671, the Natural Science Foundation of
China (Grant No. 61971234, 11501301, 62001167), the “1311 Talent Plan” of
NUPT, Hunan Provincial Key Laboratory of Mathematical Modeling and Analysis in
Engineering (Changsha University of Science & Technology) and Postgraduate
Research & Practice Innovation Program of Jiangsu Province (Grant No.
$\mathrm{SJCX19}\\_0267,\mathrm{SJCX20}\\_0229$).

## References

  * [1] X. Cai, R. Chan, M. Nikolova, and T. Zeng, “A three-stage approach for segmenting degraded color images: Smoothing, lifting and thresholding (slat),” _Journal of Scientific Computing_ , vol. 72, no. 3, pp. 1313–1332, 2017.
  * [2] L. Condat, “Discrete total variation: New definition and minimization,” _SIAM Journal on Imaging Sciences_ , vol. 10, no. 3, pp. 1258–1290, 2017\.
  * [3] W. R. Crum, C. Oscar, and D. L. G. Hill, “Generalized overlap measures for evaluation and validation in medical image analysis,” _IEEE Transactions on Medical Imaging_ , vol. 25, no. 11, pp. 1451–1461, 2006.
  * [4] H. Liu, M. Yan, E. Song, Y. Qian, X. Xu, R. Jin, L. Jin, and C.-C. Hung, “Label fusion method based on sparse patch representation for the brain mri image segmentation,” _IET Image Processing_ , vol. 11, no. 7, pp. 502–511, 2017.
  * [5] A. Namburu, S. K. Samayamantula, and S. R. Edara, “Generalised rough intuitionistic fuzzy c-means for magnetic resonance brain image segmentation,” _IET Image Processing_ , vol. 11, no. 9, pp. 777–785, 2017\.
  * [6] G. McConnon, F. Deravi, S. Hoque, K. Sirlantzis, and G. Howells, “Impact of common ophthalmic disorders on iris segmentation,” _2012 5th IAPR International Conference on Biometrics (ICB)_ , pp. 277–282, 2012.
  * [7] Y. Liu, C. He, P. Gao, Y. Wu, and Z. Ren, “A binary level set variational model with $l^{1}$ data term for image segmentation,” _Signal Processing_ , vol. 155, pp. 193–201, 2019.
  * [8] S. Luo, X.-C. Tai, L. Huo, Y. Wang, and R. Glowinski, “Convex shape prior for multi-object segmentation using a single level set function,” _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 613–621, 2019.
  * [9] L. Wang, L. Zhang, X. Yang, P. Yi, and H. Chen, “Level set based segmentation using local fitted images and inhomogeneity entropy,” _Signal Processing_ , vol. 167, p. 107297, 2020.
  * [10] S. Yan, X.-c. Tai, J. Liu, and H.-y. Huang, “Convexity shape prior for level set based image segmentation method,” _IEEE Transactions on Image Processing_ , vol. 29, pp. 7141–7152, 2020.
  * [11] R. Rahmat and D. Harris-Birtill, “Comparison of level set models in image segmentation,” _IET Image Processing_ , vol. 12, no. 12, pp. 2212–2221, 2018\.
  * [12] D. Wang and X.-P. Wang, “The iterative convolution-thresholding method (ictm) for image segmentation,” _arXiv: Computer Vision and Pattern Recognition (2019)_.
  * [13] D. Wang, H. Li, X. Wei, and X.-P. Wang, “An efficient iterative thresholding method for image segmentation,” _Journal of Computational Physics_ , vol. 350, pp. 657–667, 2017.
  * [14] T. F. Chan and L. A. Vese, “Active contours without edges,” _IEEE Transactions on image processing_ , vol. 10, no. 2, pp. 266–277, 2001.
  * [15] K. Ding, L. Xiao, and G. Weng, “Active contours driven by region-scalable fitting and optimized laplacian of gaussian energy for image segmentation,” _Signal Processing_ , vol. 134, pp. 224–233, 2017.
  * [16] R. Jin and G. Weng, “Active contours driven by adaptive functions and fuzzy c-means energy for fast image segmentation,” _Signal Processing_ , vol. 163, pp. 1–10, 2019.
  * [17] L. Wang, G. Chen, D. Shi, Y. Chang, S. Chan, J. Pu, and X. Yang, “Active contours driven by edge entropy fitting energy for image segmentation,” _Signal Processing_ , vol. 149, pp. 27–35, 2018.
  * [18] D. Mumford and J. Shah, “Optimal approximations by piecewise smooth functions and associated variational problems,” _Communications on Pure and Applied Mathematics_ , vol. 42, no. 5, pp. 577–685, 1989.
  * [19] Y. Duan, W. Huang, J. Zhou, H. Chang, and T. Zeng, “A two-stage image segmentation method using euler’s elastica regularized mumford-shah model,” _2014 22nd International Conference on Pattern Recognition_ , pp. 118–123, 2014.
  * [20] Z. Li and T. Zeng, “A two-stage image segmentation model for multi-channel images,” _Communications in Computational Physics_ , vol. 19, no. 4, pp. 904–926, 2016.
  * [21] X. Cai, R. Chan, C.-B. Schönlieb, G. Steidl, and T. Zeng, “Linkage between piecewise constant mumford–shah model and rudin–osher–fatemi model and its virtue in image segmentation,” _SIAM Journal on Scientific Computing_ , vol. 41, no. 6, pp. B1310–B1340, 2019.
  * [22] S. Zeng, X. Wang, H. Cui, C. Zheng, and D. Feng, “A unified collaborative multikernel fuzzy clustering for multiview data,” _IEEE Transactions on Fuzzy Systems_ , vol. 26, no. 3, pp. 1671–1687, 2017.
  * [23] V. Ananthi and P. Balasubramaniam, “A new thresholding technique based on fuzzy set as an application to leukocyte nucleus segmentation,” _Computer Methods and Programs in Biomedicine_ , vol. 134, pp. 165–177, 2016\.
  * [24] V. Ananthi, P. Balasubramaniam, and T. Kalaiselvi, “A new fuzzy clustering algorithm for the segmentation of brain tumor,” _Soft Computing_ , vol. 20, no. 12, pp. 4859–4879, 2016.
  * [25] X. Cai, R. Chan, and T. Zeng, “A two-stage image segmentation method using a convex variant of the mumford–shah model and thresholding,” _SIAM Journal on Imaging Sciences_ , vol. 6, no. 1, pp. 368–390, 2013.
  * [26] R. Chan, A. Lanza, S. Morigi, and F. Sgallari, “Convex non-convex image segmentation,” _Numerische Mathematik_ , vol. 138, no. 3, pp. 635–680, 2018\.
  * [27] C. Liu, M. K.-P. Ng, and T. Zeng, “Weighted variational model for selective image segmentation with application to medical images,” _Pattern Recognition_ , vol. 76, pp. 367–379, 2018.
  * [28] Q. Ma, J. Peng, and D. Kong, “Image segmentation via mean curvature regularized mumford-shah model and thresholding,” _Neural Processing Letters_ , vol. 48, no. 2, pp. 1227–1241, 2018.
  * [29] B. Shi, Z.-F. Pang, and J. Xu, “Image segmentation based on the hybrid total variation model and the k-means clustering strategy,” _Inverse Problems and Imaging_ , vol. 10, pp. 807–828, 2016.
  * [30] Z. Zhi, Y. Sun, and Z.-F. Pang, “Two-stage image segmentation scheme based on inexact alternating direction method,” _Numerical Mathematics: Theory, Methods and Applications_ , vol. 9, no. 3, pp. 451–469, 2016.
  * [31] Y. Yang, Q. Zhong, Y. Duan, and T. Zeng, “A weighted bounded hessian variational model for image labeling and segmentation,” _Signal Processing_ , p. 107564, 2020.
  * [32] R. Chan, H. Yang, and T. Zeng, “A two-stage image segmentation method for blurry images with poisson or multiplicative gamma noise,” _SIAM Journal on Imaging Sciences_ , vol. 7, no. 1, pp. 98–127, 2014.
  * [33] T. Lei, X. Jia, Y. Zhang, S. Liu, H. Meng, and A. K. Nandi, “Superpixel-based fast fuzzy c-means clustering for color image segmentation,” _IEEE Transactions on Fuzzy Systems_ , vol. 27, no. 9, pp. 1753–1766, 2018.
  * [34] X. Jia, T. Lei, X. Du, S. Liu, H. Meng, and A. K. Nandi, “Robust self-sparse fuzzy clustering for image segmentation,” _IEEE Access_ , vol. 8, pp. 146 182–146 195, 2020.
  * [35] R. Jin, C. Kou, R. Liu, and Y. Li, “A color image segmentation method based on improved k-means clustering algorithm,” _Proceedings of the International Conference on Information Engineering and Applications (IEA) 2012_ , pp. 499–505, 2013.
  * [36] M. Singh, P. Patel, D. Khosla, and T. Kim, “Segmentation of functional mri by k-means clustering,” _IEEE Transactions on Nuclear Science_ , vol. 43, no. 3, pp. 2030–2036, 1996.
  * [37] L. Condat, “A convex approach to k-means clustering and image segmentation,” _International Workshop on Energy Minimization Methods in Computer Vision and Pattern Recognition_ , pp. 220–234, 2017.
  * [38] A. Chambolle, D. Cremers, and T. Pock, “A convex approach to minimal partitions,” _SIAM Journal on Imaging Sciences_ , vol. 5, no. 4, pp. 1113–1158, 2012.
  * [39] T. Pock, D. Cremers, H. Bischof, and A. Chambolle, “Global solutions of variational models with convex regularization,” _SIAM Journal on Imaging Sciences_ , vol. 3, no. 4, pp. 1122–1145, 2010.
  * [40] L. Bar, T. F. Chan, G. Chung, M. Jung, N. Kiryati, N. Sochen, and L. A. Vese, “Mumford and shah model and its applications to image segmentation and image restoration,” _Handbook of Mathematical Methods in Imaging_ , pp. 1–52, 2014\.
  * [41] Z. Ding, J. Jia, and D. LI, “Fast clustering segmentation method combining hill-climbing for color image,” _Journal of Information and Computational Sciences_ , vol. 8, no. 14, pp. 2949–2957, 2011.
  * [42] T. Ohashi, Z. Aghbari, and A. Makinouchi, “Hill-climbing algorithm for efficient color-based image segmentation,” _IASTED International Conference on Signal Processing, Pattern Recognition, and Applications_ , pp. 17–22, 2003.
  * [43] B. Chandana, K. Srinivas, and R. K. Kumar, “Clustering algorithm combined with hill climbing for classification of remote sensing image,” _International Journal of Electrical and Computer Engineering_ , vol. 4, no. 6, p. 923, 2014.
  * [44] F. Li, M. K. Ng, T. Y. Zeng, and C. Shen, “A multiphase image segmentation method based on fuzzy region competition,” _SIAM Journal on Imaging Sciences_ , vol. 3, no. 3, pp. 277–299, 2010.
  * [45] M. Storath and A. Weinmann, “Fast partitioning of vector-valued images,” _SIAM Journal on Imaging Sciences_ , vol. 7, no. 3, pp. 1826–1852, 2014.
  * [46] T. Goldstein and S. Osher, “The split bregman method for l1-regularized problems,” _SIAM Journal on Imaging Sciences_ , vol. 2, no. 2, pp. 323–343, 2009.
  * [47] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein _et al._ , “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Foundations and Trends® in Machine learning_ , vol. 3, no. 1, pp. 1–122, 2011\.
  * [48] M. A. Figueiredo and J. M. Bioucas-Dias, “Restoration of poissonian images using alternating direction optimization,” _IEEE transactions on Image Processing_ , vol. 19, no. 12, pp. 3133–3145, 2010.
  * [49] K. Guo, D. Han, and T. Wu, “Convergence of alternating direction method for minimizing sum of two nonconvex functions with linear constraints,” _International Journal of Computer Mathematics_ , pp. 1–17, 2016.
  * [50] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” _Journal of Mathematical Imaging and Vision_ , vol. 40, no. 1, pp. 120–145, 2011.