# Fractional Order Periodic Maps: Stability Analysis and Application to the
Periodic-2 Limit Cycles in the Nonlinear Systems

Sachin Bhalekar<EMAIL_ADDRESS><EMAIL_ADDRESS>(Corresponding Author) Prashant M. Gade<EMAIL_ADDRESS>School of
Mathematics and Statistics, University of Hyderabad, Hyderabad, 500046 India
Department of Physics, Rashtrasant Tukadoji Maharaj Nagpur University, Nagpur

###### Abstract

We consider the stability of periodic map with period-$2$ in linear fractional
difference equations where the function is $f(x)=ax$ at even times and
$f(x)=bx$ at odd times. The stability of such a map for an integer order map
depends on product $ab$. The conditions are much complex for fractional maps
and depend on $ab$ as well as $a+b$. There are no superstable period-2 orbits.
These conditions are useful in obtaining stability conditions of
asymptotically periodic orbits with period-$2$ in the nonlinear case. The
stability conditions are demonstrated numerically. The formalism can be
generalized to higher periods.

††journal:

##  1 Introduction

The stability analysis of the fixed point is extremely helpful in the study of
nonlinear maps of integer order. The stability of periodic orbits is no less
important. In fact, one of the definitions of chaos is as follows: if $V$ is a
set and $F:V\rightarrow V$ is chaotic if it has a) sensitive dependence on
initial conditions, b) is topologically transitive and c) the periodic points
of $F$ are dense on $V$[1]. It has been shown that conditions b) and c) imply
a). Thus periodic points are crucial in the theory of chaos. Several invariant
properties of the chaotic attractor can be computed using unstable periodic
orbits[2]. The multifractal spectrum of the attractor can be computed using
unstable periodic orbits[3]. In a very striking result, it has been shown that
the statistical properties of turbulence can be computed using only one
unstable periodic orbit[4, 5]. We can compute Lyapunov exponent of the system
using unstable periodic orbit[6]. Several dynamic quantities can be computed
using eigenvalues of a few short fundamental cycles because they are
structured hierarchically. Longer cycles only offer higher order corrections
to these quantities[2]. Apart from invariant density, fractal dimension and
Lyapunov exponents, we can compute the topological and metric entropy of the
attractor[7]. Periodic orbits form a skeleton of chaotic attractors and
control schemes such as the Ott-Grebogi-Yorke scheme have been used in
controlling chaos and stabilizing a particular periodic orbit[8]. Several
bifurcations in the system, such as crisis, can be explained by understanding
the periodic orbits and their stable and unstable manifold[9]. (We note that
unstable periodic orbits are an important theoretical tool in studying quantum
chaos as well[10].) They have important applications, including control of
cardiac chaos[11]. Methods based on the detection of unstable periodic orbits
have been used to establish low dimensional chaos in crayfish caudal
photoreceptor [12]. They are used in the characterization, control and
prediction of experimental systems [13, 14]. In short, the importance of
periodic orbit in the theory of nonlinear dynamics and chaos cannot be
overemphasized.

In fractional order systems, the chaos theory is not as well developed as in
integer order systems. However, specific results about the stability of the
fixed point are obtained. We linearize around the fixed point and ensure
stability if the eigenvalues are inside the unit circle. In integer order
maps, the chain rule is applicable. Thus, we can study the stability of fixed
points of the function $f^{n}(x)$ for n-period orbits. We again linearize and
the stability of the n-period orbit $(x_{1},x_{2},\ldots,x_{n})$ of 1-d map is
dictated by the condition that $|f^{\prime}(x_{n})f^{\prime}(x_{n-1})\ldots
f^{\prime}(x_{1})|<1$.

Unfortunately, these conditions do not work for fractional order maps even for
fixed points and the stability is ensured if the eigenvalues are inside the
cardioid-shaped stability region in the complex plane[fixed]. This work shows
that the conditions are even more complicated for periodic points. (In
fractional maps, we have only asymptotically periodic points, not strictly
periodic ones). While it is indeed true that the fixed points are given by
$f(x)=x$ even for fractional order map, the period-n orbit is not given by
roots of equation $f^{n}(x)=x$. The fixed points of the twice iterated
nonlinear map $f^{2}(x)$ do not give the 2-period orbit for the fractional
order map, which is reached asymptotically. We will explicitly solve the
system for an asymptotically 2-period orbit. Besides, the stability of orbit
$(x_{1},x_{2})$ is not given by product $f^{\prime}(x_{1})f^{\prime}(x_{2})$.
It is also dependent in $f^{\prime}(x_{1})+f^{\prime}(x_{2})$. Thus the
stability analysis of higher period orbits for fractional maps is much more
complicated than integer order maps. The bifurcation diagrams can be
complicated as well. In integer order maps, we observe a period doubling
cascade. When certain period becomes unstable, we obtain the stable solution
with twice the period. In fractional order maps, the fixed point and period
two orbit can both be stable at same parameter value.

In this work, we first derive the analytic conditions for the stability of the
periodic map. The map is linear. However, it is different for odd and even
times. We find that the same conditions work for linearized asymptotically
period-two orbits of fractional nonlinear maps.

##  2 Preliminaries

In this section, we present some basic definitions and results. Let
$h>0,\;a\in\mathbb{R}$, $(h\mathbb{N})_{a}=\\{a,a+h,a+2h,\ldots\\}$ and
$\mathbb{N}_{a}=\\{a,a+1,a+2,\ldots\\}$.

###### Definition 2.1

(see [15, 16, 17]). For a function $x:(h\mathbb{N})_{a}\rightarrow\mathbb{R}$,
the forward h-difference operator if defined as

$(\Delta_{h}x)(t)=\frac{x(t+h)-x(t)}{h},$

where t $\in(h\mathbb{N})_{a}$.

Throughout this paper, we take $a=0$ and $h=1$.

###### Definition 2.2

[17] For a function $x:\mathbb{N_{\circ}}\rightarrow\mathbb{R}$ the fractional
sum of order $\alpha>0$ is given by

$(\Delta^{-\alpha}x)(t)=\frac{1}{\Gamma(\alpha)}\sum_{s=0}^{n}\frac{\Gamma(\alpha+n-s)}{\Gamma(n-s+1)}x(s),$
(1)

where, $t=\alpha+n,\;n\in\mathbb{N_{\circ}}$.

###### Definition 2.3

[17, 18] Let $\mu>0$ and $m-1<\mu<m$, where $m\in\mathbb{N}$,
$m=\lceil\mu\rceil$. The $\mu$th fractional Caputo-like difference is defined
as

$\Delta^{\mu}x(t)=\Delta^{-(m-\mu)}\left(\Delta^{m}x(t)\right),$ (2)

where $t\in\mathbb{N}_{m-\mu}$ and

$\Delta^{m}x(t)=\sum_{k=0}^{m}\binom{m}{k}(-1)^{m-k}x(t+k).$ (3)

Note: The discrete dynamical system $x(t+1)=f(x(t))$ can be written
equivalently as a difference equation $\Delta x(t)=f(x(t))-x(t)$ by
subtracting the term $x(t)$ from both sides. Further, we can generalize this
difference equation by replacing the operator $\Delta$ with the operator
$\Delta^{\alpha}$, where $0<\alpha<1$. If we consider the difference equations
in the higher dimensions and if the function $f$ is linear, then we get the
systems of the form Equation (LABEL:linCMLSys). This is the motivation behind
the term $A-I$. Of course, we can write the matrix $A^{\prime}=A-I$ and use
the results discussed in the literature (e.g., [19]) to analyze this system.

###### Definition 2.4

[17] The Z-transform of a sequence $\\{y(n)\\}_{n=0}^{\infty}$ is a complex
function given by $Y(z)=Z[y](z)=\sum_{k=0}^{\infty}y(k)z^{-k}$ where
$z\in\mathbb{C}$ is a complex number for which the series converges
absolutely.

###### Definition 2.5

[17] Let $\tilde{\phi}_{\alpha}(n)$ be a family of binomial functions defined
on $\mathbb{Z}$, parametrized by $\alpha$ defined by

$\displaystyle\tilde{\phi}_{\alpha}(n)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(n+\alpha-1)}{\Gamma(\alpha)\Gamma(n)}$ (8)
$\displaystyle=$ $\displaystyle\left(\begin{array}[]{c}n+\alpha-1\\\ n\\\
\end{array}\right)=(-1)^{n}\left(\begin{array}[]{c}-\alpha\\\
n\end{array}\right).$

Then

$Z(\tilde{\phi}_{\alpha}(t))=\frac{1}{(1-z^{-1})^{\alpha}},\quad|z|>1.$

###### Definition 2.6

[17] The convolution $\phi*x$ of the functions $\phi$ and $x$ defined on
$\mathbb{N}$ is defined as

$\left(\phi*x\right)(n)=\sum_{s=0}^{n}\phi(n-s)x(s)=\sum_{s=0}^{n}\phi(s)x(n-s).$

Then the Z-transform of this convolution is

$Z\left(\phi*x\right)(n)=\left(Z\left(\phi\right)(n)\right)\left(Z\left(x\right)(n)\right).$
(9)

##  3 The Model and Characteristic Equation

Let $x:\mathbb{N_{\circ}}\rightarrow\mathbb{R}$ and $f$ be a map defined by

$f\left(x(t)\right)=\begin{cases}ax(t),&\text{if $t$ is even}\\\
bx(t),&\text{if $t$ is odd},\end{cases}$ (10)

where $a$ and $b$ are real numbers.
We define the fractional order discrete dynamical system using this map as

$x(t+1)=x(0)+\sum_{j=0}^{t}\frac{\Gamma(t-j+\alpha)}{\Gamma(\alpha)\Gamma(t-j+1)}\left[f\left(x(j)\right)-x(j)\right].$
(11)

The traditional way to obtain the characteristic equation of the systems of
the form (11) is to take Z-transform and equate the coefficient of $Z[x(t)]$
to zero. The map $f$ defined in (10) can also be written as

$f\left(x(t)\right)=\frac{a+b+(-1)^{t}(a-b)}{2}x(t).$

Applying Z-transform to (11), we get

$X(z)-zx(0)=\frac{x(0)}{1-z^{-1}}+\frac{1}{(1-z^{-1})^{\alpha}}\left(\frac{a+b}{2}-1\right)X(z)+\frac{1}{(1-z^{-1})^{\alpha}}\left(\frac{a-b}{2}-1\right)X(-z).$
(12)

Note that the equation (12) cannot be used to find the characteristic equation
because of $X(-z)$ in the last term.
An elegant way to get the solution to this problem is to separate the terms
$x(t)$ with even $t$ from odd values of $t$. Let us define $p(t)=x(2t)$ and
$q(t)=x(2t+1)$. The system (11) can now be written in an equivalent form as

$\displaystyle p(t+1)$ $\displaystyle=$ $\displaystyle
x(0)+\sum_{k=0}^{t}\frac{\Gamma(2t+1-2k+\alpha)}{\Gamma(\alpha)\Gamma(2t-2k+2)}\left[(a-1)p(k)\right]$
$\displaystyle\,\,+\sum_{k=0}^{t}\frac{\Gamma(2t-2k+\alpha)}{\Gamma(\alpha)\Gamma(2t-2k+1)}\left[(b-1)q(k)\right],$
$\displaystyle q(t+1)$ $\displaystyle=$ $\displaystyle
x(0)+\sum_{k=0}^{t}\frac{\Gamma(2t+2-2k+\alpha)}{\Gamma(\alpha)\Gamma(2t-2k+3)}\left[(a-1)p(k)\right]$
(13)
$\displaystyle\,\,+(a-1)p(t+1)+\sum_{k=0}^{t}\frac{\Gamma(2t+1-2k+\alpha)}{\Gamma(\alpha)\Gamma(2t-2k+2)}\left[(b-1)q(k)\right].$

If we define $\phi_{1}(t)=\tilde{\phi}_{\alpha}(2t)$,
$\phi_{2}(t)=\tilde{\phi}_{\alpha}(2t+1)$ and
$\phi_{3}(t)=\tilde{\phi}_{\alpha}(2t+2)$ then the system (13) can be written
as

$\displaystyle p(t+1)$ $\displaystyle=$ $\displaystyle
x(0)+(a-1)\left(\phi_{2}*p\right)(t)+(b-1)\left(\phi_{1}*q\right)(t),$ (14)
$\displaystyle q(t+1)$ $\displaystyle=$ $\displaystyle
x(0)+(a-1)\left(\phi_{3}*p\right)(t)+(a-1)p(t+1)+(b-1)\left(\phi_{2}*q\right)(t).$

With a few computations, we get the Z-transforms as

$\displaystyle Z\left(\phi_{1}(t)\right)$ $\displaystyle=$
$\displaystyle\frac{\left(\sqrt{z}-1\right)^{-\alpha}+\left(\sqrt{z}+1\right)^{-\alpha}}{2z^{-\alpha/2}},$
$\displaystyle Z\left(\phi_{2}(t)\right)$ $\displaystyle=$
$\displaystyle\frac{\left(\sqrt{z}-1\right)^{-\alpha}-\left(\sqrt{z}+1\right)^{-\alpha}}{2z^{(-\alpha-1)/2}},$
$\displaystyle Z\left(\phi_{3}(t)\right)$ $\displaystyle=$
$\displaystyle-z+\frac{\left(\sqrt{z}-1\right)^{-\alpha}+\left(\sqrt{z}+1\right)^{-\alpha}}{2z^{-1-\alpha/2}},$
$\displaystyle Z\left(p(t+1)\right)$ $\displaystyle=$ $\displaystyle
zP(z)-zp(0),\,\,Z\left(q(t+1)\right)=zQ(z)-zq(0),$ (15)

where $Z(p(t))=P(z)$, $Z(q(t))=Q(z)$, $p(0)=x(0)$ and $q(0)=x(1)=ax(0)$.
Applying Z-transform to the system (14) and using (15), we get

$\displaystyle\left[z-(a-1)\frac{\left(\sqrt{z}-1\right)^{-\alpha}-\left(\sqrt{z}+1\right)^{-\alpha}}{2z^{(-\alpha-1)/2}}\right]P(z)$
$\displaystyle-(b-1)\frac{\left(\sqrt{z}-1\right)^{-\alpha}+\left(\sqrt{z}+1\right)^{-\alpha}}{2z^{-\alpha/2}}Q(z)$
$\displaystyle=$ $\displaystyle\frac{-z^{2}}{1-z}x(0),$
$\displaystyle(a-1)\frac{\left(\sqrt{z}-1\right)^{-\alpha}+\left(\sqrt{z}+1\right)^{-\alpha}}{2z^{-1-\alpha/2}}P(z)$
$\displaystyle+\left[(b-1)\frac{\left(\sqrt{z}-1\right)^{-\alpha}-\left(\sqrt{z}+1\right)^{-\alpha}}{2z^{(-\alpha-1)/2}}-z\right]Q(z)$
$\displaystyle=$ $\displaystyle\frac{z^{2}}{1-z}x(0).$ (16)

The characteristic equation of the system (13) (and hence of the system (11))
can now be obtained by equating the determinant of coefficients of the terms
$P(z)$ and $Q(z)$ in the system (16) to zero as below:

$-z(z-1)^{\alpha}-\frac{1}{2}(a+b-2)z^{\frac{1+\alpha}{2}}\left[\left(\sqrt{z}-1\right)^{-\alpha}-\left(\sqrt{z}+1\right)^{-\alpha}\right]+(a-1)(b-1)z^{\alpha}=0.$
(17)

##  4 Stable Region

The zero solution of system (11) is locally asymptotically stable if and only
if all the roots $z$ of the characteristic equation (17) satisfy $|z|<1$.
Therefore, the boundary of the stable region of the system (11) can be
obtained by substituting $z=e^{\iota t}$ in the characteristic equation (17).
We have,

$\displaystyle z-1$ $\displaystyle=$ $\displaystyle e^{\iota
t}-1=2\sin(t/2)e^{\iota(\pi+t)/2},$ $\displaystyle\sqrt{z}-1$ $\displaystyle=$
$\displaystyle e^{\iota t/2}-1=2\sin(t/4)e^{\iota(2\pi+t)/4},$
$\displaystyle\sqrt{z}+1$ $\displaystyle=$ $\displaystyle e^{\iota
t/2}+1=2\cos(t/4)e^{\iota t/4}.$ (18)

Using (18), we can rewrite the characteristic equation (17) as

$\displaystyle-2^{\alpha}\left(\sin(t/2)\right)^{\alpha}e^{\iota\left[t+\alpha(\pi+t)/2\right]}+(a-1)(b-1)e^{\iota\alpha
t}$ (19) $\displaystyle-\frac{1}{2}(a+b-2)e^{\iota
t(1+\alpha)/2}2^{\alpha}\left[\left(\sin(t/4)\right)^{\alpha}e^{\iota\alpha(2\pi+t)/4}-\left(\cos(t/4)\right)^{\alpha}e^{\iota\alpha
t/4}\right]$ $\displaystyle=$ $\displaystyle 0.$

Separating real and imaginary parts in (19), we get

$\displaystyle-2^{\alpha}\left(\sin(t/2)\right)^{\alpha}\cos\left(\frac{\alpha(\pi+t)}{2}+t\right)+(a-1)(b-1)\cos(\alpha
t)$
$\displaystyle-2^{\alpha-1}(a+b-2)\left(\sin(t/4)\right)^{\alpha}\cos\left(\frac{\alpha\pi}{2}+t\left(\frac{1}{2}+\frac{3\alpha}{4}\right)\right)$
$\displaystyle+2^{\alpha-1}(a+b-2)\left(\cos(t/4)\right)^{\alpha}\cos\left(t\left(\frac{1}{2}+\frac{3\alpha}{4}\right)\right)$
$\displaystyle=$ $\displaystyle 0,$ (20)
$\displaystyle-2^{\alpha}\left(\sin(t/2)\right)^{\alpha}\sin\left(\frac{\alpha(\pi+t)}{2}+t\right)+(a-1)(b-1)\sin(\alpha
t)$
$\displaystyle-2^{\alpha-1}(a+b-2)\left(\sin(t/4)\right)^{\alpha}\sin\left(\frac{\alpha\pi}{2}+t\left(\frac{1}{2}+\frac{3\alpha}{4}\right)\right)$
$\displaystyle+2^{\alpha-1}(a+b-2)\left(\cos(t/4)\right)^{\alpha}\sin\left(t\left(\frac{1}{2}+\frac{3\alpha}{4}\right)\right)$
$\displaystyle=$ $\displaystyle 0.$ (21)

Equation (21) is identically satisfied for $t=0$. For this value of $t$, the
equation (20) gives

$b=\frac{2\left(2^{\alpha}-1\right)+2a\left(1-2^{\alpha-1}\right)}{2\left(a-1+2^{\alpha-1}\right)}.$
(22)

This boundary curve (22) can also be written as

$\left(a-\left[1-2^{\alpha-1}\right]\right)\left(b-\left[1-2^{\alpha-1}\right]\right)=4^{\alpha-1}.$
(23)

Let us call this boundary curve as $\Gamma_{1}$. The lines $a=1-2^{\alpha-1}$
and $b=1-2^{\alpha-1}$ are asymptotes for $\Gamma_{1}$.
If we substitute $t=\pi$, then the equations (20) and (21) generate the
following common boundary curve $\Gamma_{2}$

$b=-\frac{1+2^{\alpha}-a+2^{\alpha/2}(a-2)\sin\left(\alpha\pi/4\right)}{a-1+2^{\alpha/2}\sin\left(\alpha\pi/4\right)}.$
(24)

Equivalently,

$\left(a-\left[1-2^{\alpha/2}\sin\left(\alpha\pi/4\right)\right]\right)\left(b-\left[1-2^{\alpha/2}\sin\left(\alpha\pi/4\right)\right]\right)=-2^{\alpha}\left(\cos(\alpha\pi/4)\right)^{2}.$
(25)

The lines $a=1-2^{\alpha/2}\sin\left(\alpha\pi/4\right)$ and
$b=1-2^{\alpha/2}\sin\left(\alpha\pi/4\right)$ are asymptotes for
$\Gamma_{2}$.
Furthermore, the system (20)-(21) can be solved for $a$ and $b$ as parametric
functions of $t\in[0,2\pi]$. We proceed as below:
Let us define
$s_{1}=\sin^{\alpha}\left(t/2\right),s_{2}=\sin\left(t+\alpha(\pi-t)/2\right)$,
$s_{3}=\sin^{\alpha}\left(t/4\right)$, $s_{4}=\sin\left(t(\alpha-2)/4\right)$,
$s_{5}=\cos^{\alpha}\left(t/4\right)$,
$s_{6}=\sin\left(((\alpha-2)t-2\alpha\pi)/4\right)$. Then

$\displaystyle a(t)$ $\displaystyle=$ $\displaystyle
1+\frac{-s_{1}s_{2}+\sqrt{s_{1}\left(2^{\alpha}(-s_{3}s_{4}+s_{5}s_{6})(s_{5}s_{4}-s_{3}s_{6})+s_{1}s_{2}^{2}\right)}}{s_{5}s_{4}-s_{3}s_{6}},$
$\displaystyle b(t)$ $\displaystyle=$ $\displaystyle
1-\frac{s_{1}s_{2}+\sqrt{s_{1}\left(2^{\alpha}(-s_{3}s_{4}+s_{5}s_{6})(s_{5}s_{4}-s_{3}s_{6})+s_{1}s_{2}^{2}\right)}}{s_{5}s_{4}-s_{3}s_{6}}$
(26)

is required parametric representation of the boundary curve, which we call
$\Gamma_{3}$.

###### Theorem 4.1

For $0<\alpha<1$, the region inside the boundary curves $\Gamma_{1}$,
$\Gamma_{2}$ and $\Gamma_{3}$ is bounded in the $ab$-plane. For any pair
$(a,b)$ in this bounded region, the zero solution of system (11) is locally
asymptotically stable.

Proof: We assume that $0<\alpha<1$. The intersection points between the curves
$\Gamma_{1}$ and $\Gamma_{2}$ are obtained by equating the right sides of the
equations (22) and (24). This gives the following two points:

$\displaystyle(a_{1},b_{1})$ $\displaystyle=$
$\displaystyle\left(\frac{2^{1+\alpha/2}-2\sin(\alpha\pi/4)-\sqrt{2^{\alpha}+4^{\alpha}-2^{1+3\alpha/2}\sin(\alpha\pi/4)}}{2^{\alpha/2}-2\sin(\alpha\pi/4)},\right.$
$\displaystyle\left.\frac{2^{1+\alpha/2}-2^{1+\alpha}\sin(\alpha\pi/4)-(2-2^{\alpha})\sqrt{2^{\alpha}+4^{\alpha}-2^{1+3\alpha/2}\sin(\alpha\pi/4)}}{2^{1+\alpha/2}+2^{3\alpha/2}-2^{1+\alpha}\sin(\alpha\pi/4)-2\sqrt{2^{\alpha}+4^{\alpha}-2^{1+3\alpha/2}\sin(\alpha\pi/4)}}\right),$
$\displaystyle(a_{2},b_{2})$ $\displaystyle=$
$\displaystyle\left(\frac{2^{1+\alpha/2}-2\sin(\alpha\pi/4)+\sqrt{2^{\alpha}+4^{\alpha}-2^{1+3\alpha/2}\sin(\alpha\pi/4)}}{2^{\alpha/2}-2\sin(\alpha\pi/4)},\right.$
$\displaystyle\left.\frac{2^{1+\alpha/2}-2^{1+\alpha}\sin(\alpha\pi/4)+(2-2^{\alpha})\sqrt{2^{\alpha}+4^{\alpha}-2^{1+3\alpha/2}\sin(\alpha\pi/4)}}{2^{1+\alpha/2}+2^{3\alpha/2}-2^{1+\alpha}\sin(\alpha\pi/4)+2\sqrt{2^{\alpha}+4^{\alpha}-2^{1+3\alpha/2}\sin(\alpha\pi/4)}}\right).$

As the curves $\Gamma_{j}$ are symmetric about the line $a=b$, so are these
intersection points. The point $(a_{1},b_{1})$ is above whereas
$(a_{2},b_{2})$ is below the line $a=b$ and both are in the first quadrant.
This also shows that there is no intersection between the curves $\Gamma_{1}$
and $\Gamma_{2}$ for the negative values of $a$ or $b$.
The curves $\Gamma_{2}$ and $\Gamma_{3}$ intersects each other at the points
$(a_{3},b_{3})$ and $(b_{3},a_{3})$, where

$\displaystyle a_{3}$ $\displaystyle=$ $\displaystyle
1+2^{\alpha/2}\frac{2-\alpha+\sqrt{2(2-2\alpha+\alpha^{2})\cos^{2}(\alpha\pi/4)}}{\alpha\cos(\alpha\pi/4)+(\alpha-2)\sin(\alpha\pi/4)},$
$\displaystyle b_{3}$ $\displaystyle=$ $\displaystyle
1+2^{\alpha/2}\frac{2-\alpha-\sqrt{2(2-2\alpha+\alpha^{2})\cos^{2}(\alpha\pi/4)}}{\alpha\cos(\alpha\pi/4)+(\alpha-2)\sin(\alpha\pi/4)}.$
(29)

Note that $a_{3}<0$ and $b_{3}>0$.
Furthermore, for any negative values of $a$ or $b$, the curve $\Gamma_{3}$
lies between the corresponding branches of the curves $\Gamma_{1}$ and
$\Gamma_{2}$. This shows that the region inside the boundary curves
$\Gamma_{1}$, $\Gamma_{2}$ and $\Gamma_{3}$ is bounded in the $ab$-plane.
Since the change in stability can occur only at these boundary curves and the
system (11) is stable at the origin, the bounded region mentioned above is the
stable region for the system (11). This proves the result.

The curves $\Gamma_{1}$ (blue color), $\Gamma_{2}$ (red color), $\Gamma_{3}$
(black color) and the stable region for $\alpha=0.5$ is shown in Figure 1.

Figure 1: The stable region of system (11) with $\alpha=0.5$

The stable orbit of the system (11) with $\alpha=0.5$ and $(a,b)=(0.6,0.7)$
inside the stable region in Figure 1 is shown in Figure 2. On the other hand,
the unstable orbit of this system with $(a,b)=(-2.5,3.6)$ outside the stable
region in Figure 1 is shown in Figure 3.

Figure 2: The stable orbit of system (11) with $\alpha=0.5$ and
$(a,b)=(0.6,0.7)$ Figure 3: The unstable orbit of system (11) with
$\alpha=0.5$ and $(a,b)=(-2.5,3.6)$

Note: If $\alpha=1$, then the curve $\Gamma_{3}$ merges with the branch of
$\Gamma_{1}$ in the third quadrant of $ab$-plane. Furthermore, the curves
$\Gamma_{1}$ and $\Gamma_{2}$ do not intersect each other in this case and
they get reduced to the curve $|ab|=1$. Note that $|ab|=1$ is the boundary of
the stable region for the classical map i.e., the system
$x(t+1)=f\left(x(t)\right)$, where $f$ is defined by (10). This shows that our
system and stability analysis are the continuous generalization to the
classical map and the corresponding stability.

##  5 Application to the nonlinear systems with period-2 limit cycles

It is proved that [20] the continuous-time fractional order autonomous systems
of differential equations cannot have periodic solutions. However, such
systems can have “asymptotic” periodic solutions or a limit cycle [20, 21]. In
this section, we show that the discrete-time fractional order systems

$x(t+1)=x(0)+\sum_{j=0}^{t}\frac{\Gamma(t-j+\alpha)}{\Gamma(\alpha)\Gamma(t-j+1)}\left[f\left(x(j)\right)-x(j)\right]$
(30)

also have the same property. Further, we propose a necessary and sufficient
condition for the existence of a period-2 limit cycle in the system (30).

###### Theorem 5.1

The system (30) cannot have a period-2 orbit.

Proof: If there exists the points $u$ and $v$ such that ${u,v}$ is a period-2
orbit of system (30) then $x(2k)=u,x(2k+1)=v$ for $k=0,1,2,\cdots$. Therefore,
for $t=0,1$ and $2$, the system (30) gives

$\displaystyle f(u)$ $\displaystyle=$ $\displaystyle v,$ (31) $\displaystyle
f(v)$ $\displaystyle=$ $\displaystyle(1-\alpha)v+\alpha u,$ (32)
$\displaystyle v$ $\displaystyle=$ $\displaystyle
u+\left(\frac{\alpha(\alpha+1)}{2}+1\right)\left(f(u)-u\right)+\alpha\left(f(v)-v\right),$
(33)

respectively. Using (31) and (32) in (33), we get

$\left(\alpha-\frac{\alpha+1}{2}\right)\left(u-v\right)=0.$ (34)

This implies either $u=v$ or $\alpha=1$. This contradiction shows that there
cannot be a period-2 orbit of the system (30).

###  5.1 The necessary and sufficient condition for the period-2 limit cycle
${u,v}$ in the system (30)

As in (13), we can split the system (30) as

$\displaystyle p(t+1)$ $\displaystyle=$ $\displaystyle
x(0)+\sum_{k=0}^{t}\frac{\Gamma(2t+1-2k+\alpha)}{\Gamma(\alpha)\Gamma(2t-2k+2)}\left[f(p(k))-p(k)\right]$
(35)
$\displaystyle\,\,+\sum_{k=0}^{t}\frac{\Gamma(2t-2k+\alpha)}{\Gamma(\alpha)\Gamma(2t-2k+1)}\left[f(q(k))-q(k)\right],$
$\displaystyle q(t+1)$ $\displaystyle=$ $\displaystyle
x(0)+\sum_{k=0}^{t}\frac{\Gamma(2t+2-2k+\alpha)}{\Gamma(\alpha)\Gamma(2t-2k+3)}\left[f(p(k))-p(k)\right]$
(36)
$\displaystyle\,\,+(a-1)p(t+1)+\sum_{k=0}^{t}\frac{\Gamma(2t+1-2k+\alpha)}{\Gamma(\alpha)\Gamma(2t-2k+2)}\left[f(q(k))-q(k)\right],$

where $p(t)=x(2t)$ and $q(t)=x(2t+1)$. If there exists a period-2 limit cycle
${u,v}$ in the system (30) then

$\lim_{t\to\infty}p(t)=u$

and

$\lim_{t\to\infty}q(t)=v$

.
If $t$ is very large and $k$ is very small, then the ratios of Gamma functions
in (35) and (36) become zero. On the other hand, if $k$ is very large in such
cases, then $p(k)\approx u$ and $q(k)\approx v$. Therefore, taking limit as
$t\to\infty$ and subtracting (35) from (36), we get

$\displaystyle u-v$ $\displaystyle=$
$\displaystyle(f(u)-u)\lim_{t\to\infty}\left(\sum_{k=0}^{t-1}\frac{\Gamma(2t-2k-1+\alpha)}{\Gamma(2t-2k)\Gamma(\alpha)}-\sum_{k=0}^{t}\frac{\Gamma(2t-2k+\alpha)}{\Gamma(2t-2k+1)\Gamma(\alpha)}\right)$
(37)
$\displaystyle+(f(v)-v)\lim_{t\to\infty}\left(\sum_{k=0}^{t-1}\frac{\Gamma(2t-2k-2+\alpha)}{\Gamma(2t-2k-1)\Gamma(\alpha)}-\sum_{k=0}^{t-1}\frac{\Gamma(2t-2k-1+\alpha)}{\Gamma(2t-2k)\Gamma(\alpha)}\right)$
$\displaystyle=$ $\displaystyle-[(f(u)-u)-(f(v)-v)]2^{-\alpha}.$

Similarly, taking limit as $t\to\infty$ and adding (35) and (36), we get

$u+v=2x(0)+2[(f(u)-u)+(f(v)-v)]\times[\lim_{t\to\infty}\sum_{k=0}^{2t}\frac{\Gamma(2t-k+\alpha)}{\Gamma(2t-k+1)\Gamma(\alpha)}].$
(38)

Since the limit in the equation (38) tends to infinity, and all other terms
are finite, we must have

$u+v=f(u)+f(v).$ (39)

Solving equations (37) and (39), we get

$\displaystyle f(u)$ $\displaystyle=$ $\displaystyle u+2^{\alpha-1}(v-u)$ (40)
$\displaystyle f(v)$ $\displaystyle=$ $\displaystyle v+2^{\alpha-1}(u-v).$
(41)

Note that for $\alpha=1$, the conditions (40)–(41) get reduced
$f(u)=v,\,f(v)=u$, the conditions for classical map $x(t+1)=f(x(t))$ to have a
period-2 orbit.
Recall that the period-2 orbit ${u,v}$ of the classical map is stable if
$|ab|<1$, where $a=f^{\prime}(u),\,b=f^{\prime}(v)$. Furthermore, there are
only an asymptotic period-2 orbits in the fractional order system (30) and
$f(p(k))\approx f(u)+ap(k)$, $f(q(k))\approx f(v)+bq(k)$. Therefore,
heuristically we can use the linearized stability analysis and expect that the
point $(a,b)=(f^{\prime}(u),f^{\prime}(v))$ should lie inside the stable
region of the system (13).

Thus the necessary and sufficient condition for the period-2 limit cycle
${u,v}$ in the system (30) is
(i) $f(u)=u+2^{\alpha-1}(v-u)$, $f(v)=v+2^{\alpha-1}(u-v)$, and
(ii) The point $(a,b)=(f^{\prime}(u),f^{\prime}(v))$ lies inside the stable
region bounded by the curves $\Gamma_{1}$, $\Gamma_{2}$ and $\Gamma_{3}$
defined in the Section 4. We verify this result with various well-known
systems in the subsection below.

###  5.2 Examples

###### Example 5.1

Consider the fractional order logistic map. In this case, we take the equation
(30) with $f(x)=\lambda x(1-x),$ $\lambda$ is a real parameter.

The conditions (40) and (41) give

$\displaystyle(u,v)$ $\displaystyle=$
$\displaystyle\left(\frac{(2^{\alpha}-1)+\lambda+\sqrt{(\lambda-(1-2^{\alpha}))(\lambda-(1+2^{\alpha}))}}{2\lambda},\right.$
$\displaystyle\left.\frac{(2^{\alpha}-1)+\lambda-\sqrt{(\lambda-(1-2^{\alpha}))(\lambda-(1+2^{\alpha}))}}{2\lambda}\right).$

These points are real if $\lambda>1+2^{\alpha}$.

Now, the point

$\displaystyle(a,b)$ $\displaystyle=$
$\displaystyle(f^{\prime}(u),f^{\prime}(v))$ $\displaystyle=$
$\displaystyle\left(1-2^{\alpha}+\sqrt{(\lambda-(1-2^{\alpha}))(\lambda-(1+2^{\alpha}))},1-2^{\alpha}-\sqrt{(\lambda-(1-2^{\alpha}))(\lambda-(1+2^{\alpha}))}\right)$

forms a straight line $L_{1}$, a parametric curve in $\lambda$ that intersects
the curve $\Gamma_{1}$ at $\lambda=1+2^{\alpha}$ and $\Gamma_{2}$ at
$\lambda=1+\sqrt{2^{\alpha}+2^{1+2\alpha}-2^{1+3\alpha/2}\sin(\alpha\pi/4)}$
(cf. Figure 4 for $\alpha=0.4$).

Figure 4: The points $(a,b)=(f^{\prime}(u),f^{\prime}(v))$ on the line $L_{1}$
indicates the period-2 limit cycle in the fractional order logistic map with
$\alpha=0.4$.

This shows that the fractional order logistic map has period-2 limit cycle if
and only if
$1+2^{\alpha}<\lambda<1+\sqrt{2^{\alpha}+2^{1+2\alpha}-2^{1+3\alpha/2}\sin(\alpha\pi/4)}$.
If $\alpha=0.4$ then we need $\lambda\in(2.31951,2.96595)$ for period-2 limit
cycle. For $\lambda<2.31951$, the trajectory settles down to an equilibrium
point. We can observe period-2 limit cycles when $\lambda\in(2.31951,2.96595)$
as expected (cf. Figure 5). Period-doubling is observed for $\lambda>2.96595$.

Figure 5: Period-2 limit cycle in the fractional order logistic map with
$\alpha=0.4$ and $\lambda=2.8$.

###### Example 5.2

Now, we consider the fractional order cubic map. We take the equation (30)
with $f(x)=\beta x(6-x^{2}),$ where $\beta<0$ is a real parameter.

We get the three expressions $(u_{0},v_{0})$, $(u_{1},v_{1})$ and
$(-u_{1},-v_{1})$ for the points $(u,v)$ by using the conditions (40) and
(41), where

$\displaystyle u_{0}$ $\displaystyle=$
$\displaystyle\sqrt{(6\beta+2^{\alpha}-1)/\beta},\,v_{0}=-u_{0},$
$\displaystyle u_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{12+\frac{2^{\alpha}-2}{\beta}-z},$
$\displaystyle v_{1}$ $\displaystyle=$
$\displaystyle\frac{2^{-1-\alpha}u_{1}\left((2^{\alpha}-2)\beta+(12+z)\beta^{2}\right)}{\beta},$
$\displaystyle z$ $\displaystyle=$
$\displaystyle\frac{-\sqrt{4-2^{2+\alpha}-3\times
4^{\alpha}+24(2^{\alpha}-2)\beta+144\beta^{2}}}{\beta}.$

Note that, the points $(a_{0},b_{0})=(f^{\prime}(u_{0}),f^{\prime}(v_{0}))$
form a straight-line $L_{2}$ in the $ab-$plane, where
$a_{0}=b_{0}=-12\beta-3(2^{\alpha}-1)$. This line $L_{2}$ intersects both the
branches of boundary curve $\Gamma_{1}$ at the parameter values
$\beta_{0}=(2-3\times 2^{\alpha})/12$ and $\beta_{1}=(1-2^{\alpha})/6$. For
$\alpha=0.7$, we have $\beta_{0}=-0.23946$ and $\beta_{1}=-0.104084$
respectively. If $\beta>\beta_{1}$ then the numbers $a_{0}$ and $b_{0}$ are
not real; whereas if $\beta<\beta_{0}$ then the points $(a_{0},b_{0})$ are
outside the stable region and hence the corresponding points $(u_{0},v_{0})$
don’t form period-2 limit cycle for this system.

The points $(a_{1},b_{1})=(f^{\prime}(u_{1}),f^{\prime}(v_{1}))$ form a curve
$L_{3}$ defined by $4\left(a^{2}+b^{2}\right)+\left(9\times
2^{\alpha}-18\right)(a+b)+10ab+18-9\times 2^{1+\alpha}=0$ in the $ab-$plane.
This curve $L_{3}$ intersects the boundary curves $\Gamma_{1}$ and
$\Gamma_{2}$ at the parameter values $\beta_{0}$ and $\beta_{2}$, where

$\displaystyle\beta_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{384}\left(64-33\times
2^{\alpha}+2^{1+\alpha/2}\sin(\alpha\pi/4)+2^{\alpha/2}\nu_{1}-3\sqrt{2}\nu_{2}\right.$
$\displaystyle-3\times 2^{(2+\alpha)/4}\left(179\times 2^{\alpha/2}+189\times
2^{1+3\alpha/2}-19\times 2^{\alpha/2}\cos(\alpha\pi/2)\right.$
$\displaystyle-25\times
3^{\alpha+3}\sin(\alpha\pi/4)+2^{1+\alpha/2}\sin^{2}(\alpha\pi/4)+5(2^{2+\alpha/2}-2^{1+\alpha})\nu_{1}$
$\displaystyle\left.\left.+3\times
2^{(1+\alpha)/2}\nu_{2}-3\sqrt{2}\nu_{1}\nu_{2}-6\sqrt{2}\sin(\alpha\pi/4)\nu_{2}\right)^{1/2}\right),$
$\displaystyle\nu_{1}$ $\displaystyle=$ $\displaystyle\left(34+81\times
2^{\alpha}-2\cos(\alpha\pi/2)-9\times
2^{2+\alpha/2}\sin(\alpha\pi/4)\right)^{1/4},$ $\displaystyle\nu_{2}$
$\displaystyle=$ $\displaystyle 2^{\alpha/2}\left(18+9\times
2^{\alpha}-2\cos(\alpha\pi/2)-2^{\alpha/2}\nu_{1}+2(-5\times
2^{1+\alpha/2}+\nu_{1})\sin(\alpha\pi/4)\right)^{1/2}.$

For $\alpha=0.7$, $\beta_{2}=-0.277584$.

Thus, the condition for the existence of period-2 limit cycle in the
fractional order cubic map is $\beta\in(\beta_{2},\beta_{1})$. Figure 6 shows
the curves $L_{2}$ and $L_{3}$ in the stable region for $\alpha=0.7$ and
$-0.277584<\beta<-0.104084$.

For $\alpha=0.7$, the trajectory of fractional order cubic map converges to
the period-2 point $(u_{0},v_{0})=(1.6963,-1.6963)$ when
$\beta=-0.20\in(\beta_{0},\beta_{1})$ (cf. Figure 7) and to the period-2 point
$(u_{1},v_{1})=(1.45647,-2.14495)$ when $\beta=-0.26\in(\beta_{2},\beta_{0})$
(cf. Figure 8). Note that the point $(-u_{1},-v_{1})$ indicates the existence
of “coexisting” asymptotic period-2 orbits as shown in Figure 9.

Figure 6: The curves $L_{2}$ and $L_{3}$ in the stable region for $\alpha=0.7$
and $-0.277584<\beta<-0.104084$. Figure 7: Asymptotic period-2 orbit in the
fractional order cubic map for $\alpha=0.7$ and $\beta=-0.20$. Figure 8:
Asymptotic period-2 orbit in the fractional order cubic map for $\alpha=0.7$
and $\beta=-0.26$. Figure 9: Coexisting asymptotic period-2 orbits in the
fractional order cubic map for $\alpha=0.7$ and $\beta=-0.24$.

###### Example 5.3

In this example, we discuss the conditions for asymptotic period-2 orbits in
the fractional order Gauss map (30), where $f(x)=e^{-7.5x^{2}}+\beta$ and
$\beta$ is a real parameter.

Due to the transcendental nature of the function $f$, we cannot have the exact
expressions for the points $(u,v)$ and $(a,b)$, unlike the previous examples.
Therefore, we verify the results using numerical approximations. It is
observed that, for $-0.05\leq\beta\leq 0.56$ and $\alpha=0.6$, the points
$(a,b)=(f^{\prime}(u),f^{\prime}(v))$ form a curve $L_{4}$ that remains inside
the stable region as shown in Figure 10. The system shows asymptotic period-2
orbits for all these parameter values, as expected. We did not observe
period-2 limit cycles outside this range.

Figure 10: The points on the curve $L_{4}$ for $-0.05\leq\beta\leq 0.56$ show
asymptotic period-2 orbits in the fractional order Gauss map with
$\alpha=0.6$.

##  6 Discussion

We have obtained the analytic conditions for the stability of periodic linear
map in fractional difference equations. We show that the same conditions help
us infer the stability of asymptotically periodic orbits of period-$2$ in
nonlinear fractional difference equations. This formalism can be potentially
generalized to higher periods.

Unstable periodic orbits form the skeleton of chaotic attractors in integer
order systems. They are useful in characterization, prediction and control.
Analysis of stable and unstable manifolds of periodic orbits is an
indispensable tool in the theory of dynamical systems. The presence of chaos
or the presence of stable or unstable manifolds of periodic orbits are open
questions in fractional order systems. However, finding basic stability
conditions for periodic orbit can be a useful step in formulating an analogous
theory for fractional systems.

## Acknowledgement

P. M. Gade thanks DST-SERB for financial assistance (Ref. CRG/2020/003993).

## References

  * [1] R. L. Devaney, Chaotic dynamical systems: An introduction, Addison-Wesley, 1989\.
  * [2] P. Cvitanović, Invariant measurement of strange sets in terms of cycles, Phys. Rev. Lett. 61 (1988) 2729–2732.  doi:10.1103/PhysRevLett.61.2729.
URL https://link.aps.org/doi/10.1103/PhysRevLett.61.2729

  * [3] C. Grebogi, E. Ott, J. A. Yorke, Unstable periodic orbits and the dimension of chaotic attractors, Phys. Rev. A 36 (1987) 3522–3524.  doi:10.1103/PhysRevA.36.3522.
URL https://link.aps.org/doi/10.1103/PhysRevA.36.3522

  * [4] G. KAWAHARA, S. KIDA, Periodic motion embedded in plane couette turbulence: regeneration cycle and burst, Journal of Fluid Mechanics 449 (2001) 291–300.  doi:10.1017/S0022112001006243.
  * [5] M. Kawasaki, S.-i. Sasa, Statistics of unstable periodic orbits of a chaotic dynamical system with a large number of degrees of freedom, Phys. Rev. E 72 (2005) 037202.  doi:10.1103/PhysRevE.72.037202.
URL https://link.aps.org/doi/10.1103/PhysRevE.72.037202

  * [6] R. Franzosi, P. Poggi, M. Cerruti-Sola, Lyapunov exponents from unstable periodic orbits, Phys. Rev. E 71 (2005) 036218.  doi:10.1103/PhysRevE.71.036218.
URL https://link.aps.org/doi/10.1103/PhysRevE.71.036218

  * [7] H. Fujisaka, Statistical dynamics generated by fluctuations of local lyapunov exponents, Progress of theoretical physics 70 (5) (1983) 1264–1275.
  * [8] E. Ott, C. Grebogi, J. A. Yorke, Controlling chaos, Phys. Rev. Lett. 64 (1990) 1196–1199.  doi:10.1103/PhysRevLett.64.1196.
URL https://link.aps.org/doi/10.1103/PhysRevLett.64.1196

  * [9] C. Grebogi, E. Ott, J. A. Yorke, Chaotic attractors in crisis, Phys. Rev. Lett. 48 (1982) 1507–1510.  doi:10.1103/PhysRevLett.48.1507.
URL https://link.aps.org/doi/10.1103/PhysRevLett.48.1507

  * [10] H.-J. Stöckmann, Quantum Chaos: An Introduction, Cambridge University Press, 1999\.  doi:10.1017/CBO9780511524622.
  * [11] A. Garfinkel, M. L. Spano, W. L. Ditto, J. N. Weiss, Controlling cardiac chaos, Science 257 (5074) (1992) 1230–1235.
  * [12] X. Pei, F. Moss, Characterization of low-dimensional dynamics in the crayfish caudal photoreceptor, Nature 379 (1996) 618–621.
  * [13] R. Roy, T. W. Murphy, T. D. Maier, Z. Gills, E. R. Hunt, Dynamical control of a chaotic laser: Experimental stabilization of a globally coupled system, Phys. Rev. Lett. 68 (1992) 1259–1262.  doi:10.1103/PhysRevLett.68.1259.
URL https://link.aps.org/doi/10.1103/PhysRevLett.68.1259

  * [14] L. Flepp, R. Holzner, E. Brun, M. Finardi, R. Badii, Model identification by periodic-orbit analysis for nmr-laser chaos, Phys. Rev. Lett. 67 (1991) 2244–2247.  doi:10.1103/PhysRevLett.67.2244.
URL https://link.aps.org/doi/10.1103/PhysRevLett.67.2244

  * [15] R. A. Ferreira, D. F. Torres, Fractional h-difference equations arising from the calculus of variations, Applicable Analysis and Discrete Mathematics (2011) 110–121.
  * [16] N. R. Bastos, R. A. Ferreira, D. F. Torres, Discrete-time fractional variational problems, Signal Processing 91 (3) (2011) 513–524.
  * [17] D. Mozyrska, M. Wyrwas, The z-transform method and delta type fractional difference operators, Discrete Dynamics in Nature and Society 2015 (2015).
  * [18] C. Fulai, L. Xiannan, Z. Yong, Existence results for nonlinear fractional difference equation, Advances in Difference Equations (2011).
  * [19] R. Abu-Saris, Q. Al-Mdallal, On the asymptotic stability of linear system of fractional-order difference equations, Fractional Calculus and Applied Analysis 16 (3) (2013) 613–629.
  * [20] E. Kaslik, S. Sivasundaram, Non-existence of periodic solutions in fractional-order dynamical systems and a remarkable difference between integer and fractional-order derivatives of periodic functions, Nonlinear Analysis: Real World Applications 13 (3) (2012) 1489–1497.
  * [21] S. Bhalekar, M. Patil, Singular points in the solution trajectories of fractional order dynamical systems, Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (11) (2018) 113123.


# Segmentation Guided Deep HDR Deghosting

K Ram Prabhakar, Susmit Agrawal, and R Venkatesh Babu K. R. Prabhakar, S.
Agrawal and R. V. Babu are with the Department of Computational and Data
Sciences, Indian Institute of Science, Bengaluru, KA, 560012, INDIA. e-mail:
<EMAIL_ADDRESS>
###### Abstract

We present a motion segmentation guided convolutional neural network (CNN)
approach for high dynamic range (HDR) image deghosting. First, we segment the
moving regions in the input sequence using a CNN. Then, we merge static and
moving regions separately with different fusion networks and combine fused
features to generate the final ghost-free HDR image. Our motion segmentation
guided HDR fusion approach offers significant advantages over existing HDR
deghosting methods. First, by segmenting the input sequence into static and
moving regions, our proposed approach learns effective fusion rules for
various challenging saturation and motion types. Second, we introduce a novel
memory network that accumulates the necessary features required to generate
plausible details in the saturated regions. The proposed method outperforms
nine existing state-of-the-art methods on two publicly available datasets and
generates visually pleasing ghost-free HDR results. We also present a large-
scale motion segmentation dataset of 3683 varying exposure images to benefit
the research community. The code and dataset are available at
http://val.serc.iisc.ernet.in/HDR/shdr/.

###### Index Terms:

High Dynamic Range image fusion, Exposure Fusion, Deghosting, Computational
Photography, Convolutional Neural Networks.

##  I Introduction

Human eyes can recognize a wide range of illumination present in nature, a
range significantly higher than any standard digital camera’s capability. Our
eyes can discern details in both highlight and shadow areas of the scene.
However, cameras have difficulty in capturing both ends of the brightness
spectrum. The captured image will likely have highlight or shadow regions,
where some parts are correctly exposed while a few details will be missing in
the other. High Dynamic Range (HDR) imaging is a photography technique used to
generate images with a wide range of illumination than possible with digital
cameras. HDR imaging faithfully reconstructs details in both highlight and
shadow regions as perceived by our human eyes.

Using a specialized hardware devices or sensors ([1, 2, 3]) is one way to
generate HDR image. However, such a custom device is costly for widespread
consumer use. Instead, a widely adopted way is to shoot a stack of bracketed
varying exposure low dynamic range (LDR) images to capture both ends of the
brightness spectra. The exposure stack is then merged to generate a single HDR
image with well-exposed highlights and shadows. Many methods ([4, 5, 6]) have
been proposed in the past to generate HDR from static images. Nevertheless,
the images captured in a practical application contain both camera and object
motion. Applying such static fusion methods on dynamic scenes result in
ghosting artifacts.

HDR deghosting methods in the literature can be classified into four major
categories.  1) Pixel rejection methods: These methods align the input stack
using global alignment techniques. The static pixels across all images are
merged using standard static HDR fusion techniques. Whereas, the dynamic
pixels from the affected images are rejected from the fusion process ([7, 8,
9, 10, 11, 12, 13]). As these methods reject pixels in the dynamic region,
they suffer from having only LDR content in the moving regions.  2) Alignment
methods: These methods align images using non-rigid approaches and merge them
using simple fusion techniques ([14, 15, 16, 17, 18, 19]). However, as these
methods use simple fusion techniques, they cannot handle alignment errors.  3)
Synthesis methods: These methods synthesize missing details due to saturation
in the reference image from other input images using patch-based optimization
techniques ([20, 21, 22]). They generate static sequences from the given
dynamic sequence and merge them using standard static fusion methods. These
methods introduce artifacts in challenging sequences with heavily saturated
regions (see Fig. 1), and they are computationally expensive.  4) Deep
learning methods: With the rise of deep learning, many approaches have been
proposed recently to generate HDR images with Convolutional Neural Networks
(CNN) ([23, 24, 25, 26, 27]). These methods generate significantly better HDR
images than traditional non-deep methods. However, these methods still produce
artifacts for challenging scenes with heavily saturated regions. As
highlighted in examples Fig. 1, 5 and 7, the current state-of-the-art methods
have artifacts in their results and there is room for improvement in HDR
deghosting.

Figure 1: The input varying exposure images (a) contain gradual exposure
change and moving objects. The objective of our proposed method is to combine
these input images and generate a visually pleasing HDR image without any
ghosting. As shown in zoomed regions on the right, the existing methods
produce results with artifacts or unable to recover saturated regions for
challenging sequences. In our proposed method, we simplify the problem by
dividing the inputs into static and dynamic parts. Then, we fuse them
independently to generate the final result. As highlighted in the zoomed
regions, our method generates ghosting and artifact-free results even in
regions affected by motion and saturation. Image courtesy of Kalantari and
Ramamoorthi [23]. Image best viewed zoomed-in electronic display.

We propose a novel CNN-based HDR deghosting method to address the
aforementioned issues. In our approach, we simplify the problem by segmenting
the image into dynamic or moving and static regions. Then, we merge static and
dynamic features separately using two different CNN fusion sub-networks. The
fused dynamic and static features are combined using a decoder network. Unlike
Yan et al. [25], we train a CNN segmentation network with ground truth
segmentation images for accurate static and dynamic segmentation. Also,
different from existing state-of-the-art deep learning methods like [23, 28,
25, 24, 27], our method can be extended to fuse an arbitrary number of images
and achieve better performance compared to Prabhakar et al. [26] in the same
category. Furthermore, we present a novel global memory network to improve the
quality of images predicted by model.

The existing deep learning-based methods use a single CNN model to fuse both
static and dynamic regions. Conversely, we take inspiration from the
literature and follow the divide-and-conquer approach. We divide the given
sequence into static and dynamic regions and fuse them separately. However,
unlike rejection based methods, we use CNN to learn the efficient fusion rule,
thus avoiding the shortcomings of rejection based methods. Our approach
combines the divide-and-conquer advantage of rejection methods with better
generalization capability of CNNs. In Yan et al. [25], the authors use
attention mechanisms to select the useful features; however, they do not
enforce any loss on the predicted attention maps. Hence, the network is free
to choose features that minimize the final HDR loss; may not necessarily
identify static and dynamic features (see Fig. 1). In our approach, we train a
segmentation network dedicated to identifying the dynamic regions and simplify
the fusion process. In summary, the main contributions of our work are as
follows,

  * •

We propose a motion segmentation assisted CNN for artifact-free HDR
deghosting. Our proposed method identifies the dynamic regions and learns
different fusion rules for static and dynamic, thus effectively dealing with
challenging motion and saturation conditions.

  * •

We propose a new memory network module that enhances the network’s ability to
generate rich details in heavily saturated regions.

  * •

We present a new large-scale motion segmentation dataset with 3683 varying
exposure images and their corresponding human annotated motion segmentation
maps to benefit the research community.

  * •

We perform extensive experiments to demonstrate better performance of the
proposed approach. Additionally, through rigorous ablation experiments, we
justify the importance of each module in our proposed approach.

The rest of the paper is organized as follows. We present the literature
methods and their limitations in Section II. We discuss our proposed approach
elaborately in Section III. Then, we discuss the proposed dataset and present
the experimental results in Section IV. In Section V, we discuss the extension
of the proposed method to an arbitrary number of images and computational
complexity. Finally, we conclude the paper in Section VI.

##  II Related Works

HDR deghosting has a vibrant literature history spanning almost two decades. A
wide range of algorithms was proposed in the literature to address HDR
deghosting. They can be broadly clubbed into four categories.

The first category of algorithms assumes that input images are mostly static
and contains few regions with motion. Hence, they aim to identify dynamic
pixels (pixels affected by motion) and process them separately. The static
pixels are merged using the conventional HDR merging process. Many different
approaches use different techniques to identify dynamic regions. Grosch [11]
proposes a method to identify dynamic pixels by thresholding the difference
between images. Wu et al. [29] identify dynamic pixels that violate brightness
consistency among the inputs. Gallo et al. [12] threshold the difference
between patches in logarithmic domain instead of pixels to identify dynamic
pixels. Heo et al. [30] estimate the joint probability density function
between input images and threshold it to determine the dynamic pixels. Later,
they refine it further using an energy minimization formulation. Min et al.
[31] perform multi-level thresholding of intensity histograms between input
images to recognize dynamic pixels. Raman et al. [32] develop on top of [12]
by comparing super-pixels.

Pece et al. [9] generates a median threshold bitmap by thresholding the pixels
with median value and use it to find motion regions. Eden et al. [33] produce
HDR result in two steps. The first step uses a graph cut technique to generate
output with a similar structure as that of the reference image but does not
have a complete dynamic range in all regions. Then the dynamic range of the
output is increased in the second step by borrowing properly exposed details
from other input images. Khan et al. [7] estimate weight value for each pixel
in all input images using two criteria: the probability that it is well
exposed and it belongs to the background. Then, images are combined using
estimated weight maps. Zhang et al. [10] detect dynamic pixels by comparing
image gradients. Granados et al. [5] propose an optimization method to
identify a coherent subset of images that be combined to generate output
without ghosting artifacts.

After identifying dynamic regions, these methods ignore images with dynamic
content and use information from the rest of the static images. Whereas, for
the static regions, they make use of contents from all images. There are two
sub-classes of algorithms within this category. In the first subclass, the
dynamic regions are replaced by content from a selected reference image ([11,
31, 30]). In the second subclass, the dynamic regions are filled by
information from static images; thus, these methods will predict static
regions in all pixels ([34, 33, 7]). These methods are fast and work better
for mostly static scenes. Their main drawback, however, is producing LDR
content in large moving regions.

The second class of algorithms aligns the input images with a selected
reference image to generate a static sequence. Then the conventional HDR merge
algorithm is used to fuse them. Ward’s [35] method aligns images through
homography estimated by comparing median threshold bitmaps of the individual
images. Tomaszewska et al. [36] use SIFT features followed by RANSAC to
compute homography for aligning images. Bogoni’s [14] method handles camera
motion by registering images using a global affine transform and object motion
with the optical flow. Kang et al. [15] generate HDR video from alternating
exposure sequences. They estimate bidirectional optical flow between the
current frame and consecutive frames to align and merge images. Zimmer et al.
[17] propose a optimization-based optical flow method to align images with two
energy terms. The data term ensures that images are aligned properly, while
the smoothness term ensures smooth flow among neighboring pixels. Gallo et al.
[19] introduce a fast algorithm suited for real-time mobile processing. Their
method speeds up the operation by computing optical flow only at selected
sparse locations and interpolating it for the other pixels. Compared to the
first class of algorithms, these methods can generate HDR content in moving
regions. In general, registration-based methods [35, 36] may fail in case of
complex deformable motions. Additionally, optical flow can produce erroneous
flow in heavily saturated regions.

The third class of algorithms synthesizes static sequences from dynamic
sequences using patch-match ([37]) or similar techniques. They begin with
selecting one of the input images as a reference. They then synthesize a new
sequence with the same structure as the reference but with same exposure as
the corresponding image in the input sequence. In the synthesized result, the
structure will resemble a reference where the reference is well exposed.
Whereas, in poorly exposed regions of reference, the structure is borrowed
from other images. Sen et al. [21] use multi-scale bidirectional similarity
metric for locating similar patches. While Hu et al. [20] make use of a patch-
match algorithm ([37]) for the same. This class of algorithms can handle
complex deformable motion and hallucinate details in the reference saturated
regions. While [21, 20] methods produce high-quality results, they are
computationally demanding.

The last class algorithm uses deep learning-based techniques to perform HDR
deghosting. Kalantari and Ramamoorthi [23] begin with aligning images using
optical flow. Later, they merge the aligned images using a CNN instead of a
conventional HDR merge process. The CNN is trained to merge them by ignoring
the artifacts introduced by optical flow. Wu et al. [24] method proposes a
simple CNN-based method to generate the final HDR image directly from
misaligned input images. Yan et al. [28] method uses a multi-scale CNN model
to reconstruct the HDR result in a coarse-to-fine strategy. In another work,
Yan et al. [25] propose a method to remove ghosting artifacts using an
attention mechanism. Yan et al. [38] propose a CNN model with non-local
module. A non-local module identifies matching neighbor features to fill in
the saturated regions in their approach. Prabhakar et al. [26] introduced a
scalable CNN method that can fuse an arbitrary number of images. They achieve
scalability by aggregating the mean and max of all input feature maps.
Recently, Prabhakar et al. [27] proposed a method to process very high-
resolution images using Bilateral Guided Upsampler (BGU). While the recent
deep learning-based methods have produced better ghost-free results than the
traditional methods, their results are still inaccurate in challenging
conditions (see Fig. 1, 5, 6, 7).

##  III Proposed method

With $N$ varying exposure LDR images, $\\{\mathbf{X}_{k}\\}$, $\forall
k=(1,\dots,N)$ as input, the goal of our approach is to merge them into a
single HDR image ($\mathbf{Y}$) without any ghosting artifacts.

Motivation: Assuming that the input sequence is camera motion aligned, the
challenge is to merge them while considering the object motion. We begin by
choosing an image with the least number of saturated pixels as the reference
image. The final result will resemble the chosen reference image in the
dynamic regions while combining HDR contents from all other images in static
regions. Hence, different fusion rules are required depending on the presence
or absence of object motion. Existing deep learning-based methods employ a
monolithic CNN architecture to learn fusion rules applicable to both
scenarios. Despite their success in typical scenes, they introduce artifacts
in some challenging conditions with regions affected by both motion and
saturation.

To address this problem, we segment the input images into static and dynamic
regions (see Fig. 2) and fuse them separately. Hence, the model has the
freedom to learn different fusion rules for static and dynamic regions.
Finally, we merge the fused static and dynamic features to generate fused HDR
output.

Formulation: For ease of understanding, we explain our approach with three
input varying exposure images, $N=3$; however, we can easily extend our method
to fuse arbitrary-length sequence as explained in Section V. Following [23,
24, 25, 26, 27], we choose the middle image ($\mathbf{X}_{2}$) as reference.
The generated $\mathbf{Y}$ will have the same structural details as of
reference in well-exposed regions of $\mathbf{X}_{2}$, and borrow details from
other images in ill-exposed regions of $\mathbf{X}_{2}$.

Figure 2: Overview of our proposed approach. We use a simple four-layer
encoder to generate $\mathbf{E}_{i}$ features for inputs $\mathbf{X}_{i}$. We
then segment $\mathbf{E}_{i}$ into static ($\mathbf{S}_{i}$) and dynamic
regions ($\mathbf{D}_{i}$) with the help of a segmentation map generated by
$\mathcal{S}$. $\mathbf{S}_{i}$ and $\mathbf{D}_{i}$ of different input images
are concatenated and fused by $\mathcal{F}_{S}$ and $\mathcal{F}_{D}$.
Finally, the fused static and dynamic features are merged by a decoder to
generate the final HDR image.

####  III-1 Feature encoder ($\mathcal{E}$)

We begin by extracting convolutional features for input varying exposure
images using a encoder $\mathcal{E}$. As the HDR domain data of input LDR
images help in locating misalignments, we concatenate both HDR and LDR domain
of $\mathbf{X}_{k}$ to form a six-channel input to $\mathcal{E}$. The LDR
images can be converted to HDR domain by,
$\mathbf{H}_{k}=\mathbf{X}^{2.2}_{k}/t_{k}$, $\forall k=(1,2,3)$, where
$t_{k}$ denote the exposure time of $\mathbf{X}_{k}$. Each input image is fed
to a shared encoder to extract individual feature maps $\\{\mathbf{E}_{k}\\}$,

$\mathbf{E}_{k}=\mathcal{E}\big{(}\\{\mathbf{X}_{k},\mathbf{H}_{k}\\}\big{)},\forall
k=(1,2,3),$ (1)

####  III-2 Dynamic region segmentation ($\mathcal{S}$)

In order to identify dynamic regions with respect to the reference image, we
use a separate segmentation network $\mathcal{S}$. We transfer the exposure of
the reference image $\mathbf{X}_{2}$ to non-reference images $\mathbf{X}_{1}$
and $\mathbf{X}_{3}$:

$\mathbf{H^{ref}_{i}}=X^{2.2}_{i}/t_{ref},\forall i=(1,3)$ (2)

and concatenate non-reference image ($k=1,3$) with the reference image and
pass as input to the segmentation network $\mathcal{S}$,

$\mathbf{W}_{i}=\mathcal{S}\big{(}\\{\mathbf{H^{ref}_{i}},\mathbf{H_{ref}}\\}\big{)},\forall
i=(1,3),$ (3)

The output of $\mathcal{S}$ is a single channel segmentation map
($\mathbf{W}_{i}$) at the same resolution as input, with values in
$\\{0,1\\}$. While the value ‘0’ indicates the absence of motion between
$\mathbf{X}_{i}$ and $\mathbf{X}_{2}$, value ‘1’ indicates the presence of
motion. During this process, we generate one segmentation map for each one of
the non-reference images. Then, we obtain the segmentation map for the
reference image by taking a union of the non-reference image segmentation
maps,

$\mathbf{W}_{2}=\mathbf{W}_{1}\displaystyle\bigcup\mathbf{W}_{3},$ (4)

$\mathbf{W}_{2}$ contains the motion in all sequences with respect to the
reference image $\mathbf{X}_{2}$. Using the predicted segmentation maps
$\mathbf{W}_{k}$, we divide the input image feature maps $\mathbf{E}_{k}$ into
static ($\mathbf{S}_{k}$) and dynamic features ($\mathbf{D}_{k}$),

$\displaystyle\mathbf{D}_{k}$
$\displaystyle=\mathbf{W}_{k}\times\mathbf{E}_{k},$ $\displaystyle\forall
k=1,2,3$ (5) $\displaystyle\mathbf{S}_{k}$
$\displaystyle=\big{[}1-\mathbf{W}_{k}\big{]}\times\mathbf{E}_{k},$
$\displaystyle\forall k=1,2,3$ (6)

For $\mathcal{S}$, we use two separate U-net [39] styled architectures for
predicting $\mathbf{W}_{1}$ and $\mathbf{W}_{3}$. It consists of an encoder
with four convolution blocks, with max-pooling after each block. We pass the
encoded features to a decoder with convolutional and upsampling layers.
Additionally, we also concatenate the encoder features with the decoder using
skip connections.

####  III-3 Fusion networks

As the static features do not have any moving regions, they can be fused by
combining best-exposed image features across all images. For dynamic features,
the structure should resemble the reference image in non-saturated regions and
should have borrowed HDR details from other images in saturated regions. We
achieve this objective by merging static and dynamic features using two
separate convolutional blocks: $\mathcal{F}_{S}$ and $\mathcal{F}_{D}$,
respectively.

The concatenated static features is fed as input to $\mathcal{F}_{S}$ and
concatenated dynamic features to $\mathcal{F}_{D}$. As both $\mathcal{F}_{S}$
and $\mathcal{F}_{D}$ are separate convolutional blocks, they learn two
different fusion rules for static and dynamic feature fusion.
$\mathcal{F}_{S}$ and $\mathcal{F}_{D}$ outputs the fused static feature map
($\mathbf{S}_{f}$) and a fused dynamic feature map ($\mathbf{D}_{f}$),

$\displaystyle\mathbf{S}_{f}$
$\displaystyle=\mathcal{F}_{S}\bm{(}\mathrm{CONCAT}(\\{\mathbf{S}_{k}\\})\bm{)}$
(7) $\displaystyle\mathbf{D}_{f}$
$\displaystyle=\mathcal{F}_{D}\bm{(}\mathrm{CONCAT}(\\{\mathbf{D}_{k}\\})\bm{)}$
(8)

####  III-4 Global Memory

The proposed approach’s performance can be bolstered and generalized to $N$
images using a concept we call Global Memory. The architecture features an
external memory, which is essentially a collection of concatenated encoder
feature maps. This memory is global in that it can be accessed using special
convolutional blocks anywhere in the network. We use these blocks for writing
to and reading from memory. The memory used here is unbounded; thus the size
grows with the number of images.

Write operation: The individual image features are added to the global memory
using a write operation. Each write operation appends 32 dimension feature map
to the memory (See Fig. 2). The operation consists of two parts: convolution
of the input features with trainable kernels and concatenation of the result
to the global memory. We apply the write operation to the features extracted
during the encoding process.

Input: Varying exposure images $\\{\mathbf{X}_{k}\\}$ and their exposure times
$t_{k}$, $\forall k=(1,2,3)$

Output: Fused HDR image, $\mathbf{Y}$

/* Extract features with Encoder($\mathcal{E}$) */

1 $\\{\mathbf{E}_{k}\\}$ =
$\mathcal{E}\bm{(}\\{\mathbf{X}_{k},\mathbf{H}_{k}\\}\bm{)},where\>\mathbf{H}_{k}=\mathbf{X}^{2.2}_{k}/t_{k}$

/* Segment dynamic regions with $\mathcal{S}$ */

2 $\\{\mathbf{W}_{i}\\}$ =
$\mathcal{S}\bm{(}\\{\mathbf{X}_{i},\mathbf{X}_{2}\\}\bm{)},\forall i=(1,3)$

$\mathbf{W}_{2}$ = $\mathbf{W}_{1}\displaystyle\bigcup\mathbf{W}_{3}$

/* Extract dynamic and static parts */

3 $\\{\mathbf{D}_{k}\\}$ = $\\{\mathbf{W}_{k}\\}\times\\{\mathbf{E}_{k}\\}$

$\\{\mathbf{S}_{k}\\}$ =
$\bm{[}1-\\{\mathbf{W}_{k}\\}\bm{]}\times\\{\mathbf{E}_{k}\\}$

/* Fuse dynamic and static parts */

4 $\mathbf{D}_{f}$ =
$\mathcal{F}_{D}\bm{(}\mathrm{CONCAT}(\\{\mathbf{D}_{k}\\})\bm{)}$

$\mathbf{S}_{f}$ =
$\mathcal{F}_{S}\bm{(}\mathrm{CONCAT}(\\{\mathbf{S}_{k}\\})\bm{)}$

/* Merge fused dynamic and static features */

5 $\mathbf{Z}$ =
$\mathbf{W}_{2}\times\mathbf{D}_{f}+\bm{(}1-\mathbf{W}_{2}\bm{)}\times\mathbf{S}_{f}$

/* Generate fused HDRI with decoder($\mathcal{D}$) */

6 $\mathbf{Y}$ = $\mathcal{D}\bm{(}\mathbf{Z}\bm{)}$

Algorithm 1 Overview of proposed method

Read operation: The aggregated global memory features are tapped during the
fusion process with a read operation. Each read from memory derives a fixed
number of channels from memory. This number is also set to 32 channels in our
implementation. The operation consists of a convolution of memory features
with trainable kernels, followed by a concatenation of the results with the
network’s main branch. Since the number of channels in the memory is variable,
we use grouped convolutions to reduce the channels to a fixed size. Each
group’s size depends on the number of channels in memory and the number of
channels to be appended to the main network. Instead of concatenating the
memory features during a read operation, we chose to distill the memory and
extract only the necessary features required for the current stage. We achieve
this by applying a 1$\times$1 convolution filter on the global memory. As
shown in Fig. 2, multiple individual read operations are applied during
$\mathcal{F}_{S}$ and $\mathcal{F}_{D}$.

Figure 3: Architecture of SDC block.

####  III-5 Feature decoder ($\mathcal{D}$)

The fused static ($\mathbf{S}_{f}$) and dynamic ($\mathbf{D}_{f}$) features
are combined using the reference image segmentation map, $\mathbf{W}_{2}$. The
non-zero pixels in $\mathbf{W}_{2}$ indicate the dynamic regions in all input
images, and zero pixels indicate otherwise. Hence, $\mathbf{D}_{f}$ is chosen
in the dynamic regions and $\mathbf{S}_{f}$ in the static regions by,

$\mathbf{Z}=\mathbf{W}_{2}\times\mathbf{D}_{f}+\bm{(}1-\mathbf{W}_{2}\bm{)}\times\mathbf{S}_{f}$
(9)

We pass the combined feature $\mathbf{Z}$ as input to a series of three
densely connected Stacked and Dilated Convolutional (SDC) blocks. SDC blocks
introduced by Schuster et al. [40] consists of four parallel convolutional
layers with dilation rate of $(1,2,3,4)$. We modify the original SDC block to
have three stacked convolutional layers with dilation rate of $(1,2,4)$ (see
Fig. 3). Additionally, to reduce the number of parameters, the stacked
convolutional layers share the weights between them. An SDC block consists of
three concurrent stacked convolutional layers. The stacked convolutional
layers contain 32 filters with 3$\times$3 kernel. The last three stacked
layers’ features are concatenated and further processed by a 1$\times$1
convolution layer.

The individual outputs of each SDC block are concatenated and passed to a
1$\times$1 convolutional layer to generate the final output $\mathbf{Y}$.

####  III-6 Loss functions

The loss to train our model consists of two sub losses. The first sub-loss is
computed between the predicted ($\mathbf{W}_{i}$) and ground truth
($\hat{\mathbf{W}}_{i}$) segmentation maps as,

$\mathcal{L}_{seg}=E_{D}(\mathbf{W}_{1},\hat{\mathbf{W}}_{1})+E_{D}(\mathbf{W}_{3},\hat{\mathbf{W}}_{3})$
(10)

where $E_{D}$ denotes the Dice loss error function ([41]). The second sub-loss
consists of $\ell_{2}$ loss between predicted $\mathbf{Y}$ and ground truth
HDR ($\hat{\mathbf{Y}}$) images. As it is conventional to display the HDR
images after tone mapping, we compute the loss after tone mapping the HDR
images with $\mu$-law operator. The $\mu$-law tone mapping function is defined
as,

$T(\mathbf{Y})=\dfrac{log(1+\mu\mathbf{Y})}{log(1+\mu)}$ (11)

where $\mu$ = 5000. Then, the second sub-loss is defined as,

$\mathcal{L}_{HDR}=\ell_{2}(T(\mathbf{Y}),T(\hat{\mathbf{Y}}))$ (12)

where $\ell_{2}$ denotes the mean squared error function. The final loss is
defined as weighted sum of two sub-losses,

$\mathcal{L}_{final}=\alpha\times\mathcal{L}_{seg}+\beta\times\mathcal{L}_{HDR}$
(13)

To ensure convergence, we initially use high $\alpha$ of 1 with low $\beta$
value of $10^{-4}$. Then after 50 epochs, we change $\alpha$ to $10^{-2}$ and
$\beta$ to 1.

##  IV Evaluation and Results

###  IV-A Multi-Exposure Dynamic Motion Segmentation dataset

In the field of HDR deghosting, since most approaches involve non-deep methods
for solving the problem, there is a lack of a large data repository of multi-
exposure image sequences needed for training deep neural networks. To tackle
this issue, we create MEDS - a Multi-Exposure Dynamic motion Segmentation
dataset. MEDS consists of 3683 multi-exposure images with segmentation
annotations of moving regions.

Figure 4: Example images and their annotated motion maps from MEDS dataset.
(a) Images captured with GoPro Hero 5 camera, (b) Images captured with Canon
EOS 500D and 600D cameras.

Among the 3683 images, 2660 images were captured with a GoPro Hero 5 camera.
Each exposure stack in this set has five images with $\\{-2,-1,0,+1,+2\\}$EV
bias. The remaining 1023 images were captured using Canon EOS 500D and 600D
cameras with each exposure stack having sequences of three images with
$\\{-2,0,+2\\}$EV bias. The images were captured with a tripod-mounted camera
in Auto-exposure bracketing mode by varying only the exposure time while the
ISO and aperture values are fixed. The dataset captures a wide variety of
scenes (indoor, outdoor, landscapes), different lighting conditions with
varying natural motion (slow-moving persons to fast-moving vehicles). Sample
images from the MEDS dataset are shown in Fig. 4.

Annotation: For addressing the task of dynamic moving object segmentation from
the registered input images, we provide human-annotated motion masks for 3683
scenes. The annotators visually inspect a GIF image created from the sequence
for a few seconds to identify the dynamic regions. Then, the annotators draw a
boundary around the moving regions. For each image, we obtain two different
annotations by two different individuals. The final annotation for the image
is obtained by taking the union of the two. In mismatch regions between the
two, a final annotator inspects the difference and chooses the correct one. In
total, six annotators were involved in creating the segmentation dataset. We
plan to release this dataset of images and their corresponding dynamic moving
object annotations publicly to benefit the HDR community. In Fig. 6, we
present example images, their ground truth, and predicted segmentation maps.

###  IV-B Implementation

We have implemented our model in TensorFlow ([42]). The weights are
initialized with Glorot uniform initialization ([43]) and was trained on a
machine with Intel core i7-9700F CPU and a NVIDIA Quadro RTX 6000 GPU. We use
Adam optimizer ([44]) with $10^{-4}$ learning rate and a batch size of four to
train the model for 200 epochs. The learning rate is reduced by a factor of
0.96 after every epoch. We train the model using Kalantari et al. [23] and
Prabhakar et al. [26] datasets. [23] dataset consists of 74 training and 15
validation sequences, while [26] dataset consists of 466 training and 116
validation sequences. In addition to [23] and [26], we have also validated our
model on publicly available datasets like Sen et al. [21] and Tursun et al.
[45].

As the ground truth deghosted HDR image is unavailable for MEDS dataset
images, they were used to pre-train the segmentation network only. Later, the
pre-trained segmentation network weights are used as initialization while
training with [23] dataset. We train the model with random patches of size
128$\times$128 cropped from random locations in the full image.

###  IV-C Quantitative results

TABLE I: Quantitative comparison between proposed method against nine state-
of-the-art methods. Best score highlighted in Blue with bold and second-best
score is highlighted in gray. Refer Section IV-C for more details.

| PSNR | SSIM |
---|---|---|---
| T1 | T2 | T3 | T4 | T5 | T1 | T2 | T3 | T4 | T5 |  HDR- VDP-2
M1 | 31.25 | 35.75 | 30.18 | 28.87 | 31.69 | 0.941 | 0.963 | 0.940 | 0.942 | 0.944 | 62.07
M2 | 38.57 | 40.94 | 27.81 | 30.33 | 31.35 | 0.971 | 0.978 | 0.955 | 0.966 | 0.969 | 64.74
M3 | 8.846 | 21.33 | 12.26 | 12.10 | 20.65 | 0.107 | 0.622 | 0.715 | 0.787 | 0.738 | 54.00
M4 | 14.21 | 14.13 | 25.23 | 20.26 | 26.82 | 0.350 | 0.882 | 0.925 | 0.923 | 0.935 | 57.95
M5 | 41.27 | 42.74 | 34.12 | 33.70 | 32.99 | 0.981 | 0.987 | 0.980 | 0.979 | 0.974 | 66.10
M6 | 40.91 | 41.65 | 34.98 | 34.54 | 33.69 | 0.986 | 0.986 | 0.982 | 0.982 | 0.977 | 67.44
M7 | 39.68 | 40.47 | 33.56 | 34.08 | 32.60 | 0.980 | 0.975 | 0.966 | 0.975 | 0.961 | 66.50
M8 | 41.33 | 42.82 | 35.18 | 36.94 | 36.29 | 0.986 | 0.989 | 0.984 | 0.985 | 0.982 | 67.15
M9 | 41.08 | 41.21 | 29.83 | 33.28 | 29.51 | 0.989 | 0.989 | 0.962 | 0.978 | 0.974 | 67.53
M10 | 41.47 | 42.03 | 35.70 | 38.48 | 36.67 | 0.992 | 0.992 | 0.988 | 0.989 | 0.986 | 67.26
(a)

| PSNR | SSIM |
---|---|---|---
| T1 | T2 | T3 | T4 | T5 | T1 | T2 | T3 | T4 | T5 |  HDR- VDP-2
M1 | 29.47 | 32.58 | 27.61 | 26.84 | 26.20 | 0.954 | 0.949 | 0.912 | 0.919 | 0.917 | 63.50
M2 | 32.93 | 33.43 | 30.22 | 29.15 | 30.91 | 0.972 | 0.964 | 0.950 | 0.948 | 0.950 | 65.47
M3 | 9.760 | 8.980 | 11.18 | 13.07 | 20.60 | 0.132 | 0.641 | 0.675 | 0.763 | 0.711 | 55.76
M4 | 14.19 | 15.66 | 22.47 | 22.04 | 24.97 | 0.442 | 0.869 | 0.879 | 0.897 | 0.904 | 58.74
M5 | 32.50 | 35.63 | 30.08 | 28.81 | 30.45 | 0.969 | 0.961 | 0.938 | 0.943 | 0.940 | 65.40
M6 | 34.40 | 38.03 | 33.35 | 31.82 | 32.66 | 0.977 | 0.971 | 0.962 | 0.957 | 0.956 | 66.59
M7 | 32.74 | 36.08 | 30.66 | 29.83 | 30.54 | 0.967 | 0.959 | 0.942 | 0.935 | 0.939 | 66.10
M8 | 34.98 | 38.30 | 32.99 | 31.76 | 32.53 | 0.978 | 0.970 | 0.960 | 0.952 | 0.953 | 66.25
M9 | 35.28 | 38.65 | 33.82 | 32.08 | 33.31 | 0.980 | 0.973 | 0.963 | 0.961 | 0.957 | 66.88
M10 | 37.21 | 39.71 | 34.99 | 32.83 | 33.25 | 0.984 | 0.977 | 0.969 | 0.965 | 0.962 | 68.03
(b)

In Table Ib, we report the quantitative results comparison between proposed
method and nine state-of-the-art approaches for [23] and [26] datasets. The
compared methods are:  (i) Hu13 - Hu et al. [20],  (ii) Sen12 - Sen et al.
[21],  (iii) Endo17 - Endo et al. [48],  (iv) Eilertsen17 - Eilertsen et al.
[49],  (v) Kalantari17 - Kalantari et al. [23],  (vi) Wu18 - Wu et al. [24],
(vii) Prabhakar19 - Prabhakar et al. [26],  (viii) Prabhakar20 - Prabhakar et
al. [27] (ix) Yan19 - Yan et al. [25].  Among the nine comparison methods,
Hu13 and Sen12 are classical non-deep approaches, whereas Endo17 and
Eilertsen17 are deep learning-based single image HDR generation methods. The
rest of them are deep learning-based state-of-the-art approaches for HDR
deghosting.

For evaluation, we use three full-reference metrics: PSNR, SSIM [50], and HDR-
VDP-2 [51]. HDR-VDP-2 is a full reference metric designed specifically to
evaluate HDR images in the linear domain. We compute PSNR and SSIM in the
linear domain and after applying four different tone mapping operators. One of
them is the $\mu$-law tonemapper, which was used to train the network in Eqn.
11. Hence, we evaluate with three different tone mappers other than the
$\mu$-law tone mapper for more representative or realistic evaluation. Those
three tone mappers are: Krawczyk et al. [46], Reinhard et al. [52] and Durand
et al. [47]. It should be noted that the network was not trained to optimize
these three tone mappers. In total, we compare our proposed approach against
nine state-of-the-art methods using eleven metrics.

TABLE II: Ablation experiments. Refer Section IV-D for more details. Method | Description | PSNR-L | PSNR-T
---|---|---|---
A1 | End-to-end training | 41.07 | 41.48
A2 | End-to-End training + segmentation loss | 40.87 | 41.38
A3 | No segmentation map | 40.58 | 40.69
A4 | Simple difference as segmentation map | 40.37 | 41.13
A5 | Simple difference + Triangle fusion | 34.40 | 33.50
A6 | CNN segmentation + Triangle fusion | 36.75 | 34.30
A7 | Arbitrary size memory network | 41.05 | 41.91
A8 | Single R/W submodule | 40.54 | 40.63
Decoder architecture ablations
A9 | Vanilla | 41.14 | 41.31
A10 | Res-Net [53] | 41.11 | 41.95
A11 | SDC without dense connection | 41.36 | 42.02
A12 | SDC with dense connection | 41.47 | 42.03

In Table Iba, we present the quantitative comparison on Kalantari et al. [23]
dataset. In terms of PSNR, the proposed method outperforms others in all tone
mappers except $\mu$-law tone mapper. In terms of SSIM, the proposed method
outperforms others in all categories. In Table Ibb, we present the
quantitative comparison on Prabhakar et al. [26] dataset. The proposed method
outperforms all other approaches in ten out of eleven metrics.

Figure 5: Qualitative comparison between different methods for examples from
[26] dataset. The results are shown after tonemapping with Durand et al. [47]
method. Image best viewed zoomed in electronic monitors.

###  IV-D Ablation analysis

In Table II, we present results for various baseline ablation experiments. All
these ablation experiments were trained and tested on Kalantari et al. [23]
dataset.

TABLE III: Comparison of two models trained with and without memory network for shared and unshared weights between $\mathcal{F}_{S}$ and $\mathcal{F}_{D}$. Refer Section IV-D for more details. | PSNR-L |  | HDR-VDP-2
---|---|---|---
| Shared | Unshared |  | Shared | Unshared
With memory | 41.10 | 41.47 |  | 67.11 | 67.26
Without memory | 40.54 | 40.68 |  | 67.01 | 66.90

  * (A1)

We train our proposed method in an end-to-end fashion. Similar to Yan et al.
[38], we don’t enforce loss for segmentation network. The final HDR
reconstruction loss is used to update both segmentation as well as fusion
network.

  * (A2)

With the same setting as A1, we train the model to minimize both HDR
reconstruction and the segmentation loss. This baseline experiment performs
slightly poorer than A1. This is because the fusion network is dependent on
the segmentation map for its operation. However, by updating both segmentation
and fusion network, the training becomes unstable and stuck at local minima.

  * (A3)

The segmentation map is not used to divide the features into dynamic and
static features. $\\{\mathbf{E}_{k}\\}$ features are fused by a single fusion
network and further processed by SDC blocks to generate the final result. As
shown in Fig. 10, without the segmentation maps, the output still has some
mild ghosting artifact around the edges.

  * (A4)

Instead of using a CNN to predict the segmentation map, we threshold the
difference between brightness normalized source and reference image with 0.1
value to generate the motion map. Such a motion segmentation map has
artifacts; even the illumination changes are detected as dynamic regions. Due
to this, the overall PSNR is lower than predicting the motion map with a CNN.

  * (A5)

The motion map generated by the simple difference approach from A4 is used to
compensate for motion among the input images. For the dynamic regions, we
replace source image pixels with exposure compensated reference image pixels.
This results in three static varying exposure images. Then, the static images
are fused using standard triangle function [4]. As expected, the segmentation
error propagates to the final result without any correction. Hence,
performance is inferior.

  * (A6)

We replace the simple difference map in A5 with the motion map predicted by
the trained CNN network. This method performs better than A5; however, the CNN
model still produces output with few artifacts.

  * (A7)

We extend the proposed approach to fuse an arbitrary number of input images.
See Section V for more details.

  * (A8)

Instead of using different convolution weights for each instance of reading
and writing operation, we share the weights across all read operations.
Similarly, we share the weights across all write operations.

  * (A9)

We use nine vanilla convolution layers instead of dense SDC blocks.

  * (A10)

We use Res-Net [53] for decoder model.

  * (A11)

SDC without dense connection is used.

  * (A12)

SDC with dense connection is used for decoder architecture.

In Table III, we present ablations results for using memory and sharing fusion
networks. By sharing the fusion networks $\mathcal{F}_{S}$ and
$\mathcal{F}_{D}$, a single fusion rule is used to merge both static and
dynamic regions. As evidenced by both PSNR-L and HDR-VDP-2 metrics, having a
memory network with unshared fusion networks outperforms other ablation
methods. In Fig. 11, we present qualitative results generated by models
trained with and without memory features. As seen in the figure, the network
trained without memory fails to hallucinate viable details in saturated and
occluded regions. With the support of memory features, the model can choose
details for filling in saturated regions and avoid unnatural artifacts.

Figure 6: (a)-(c) Input sequence, (d) - combined ground-truth motion
segmentation map of (a) and (b), (e) - predicted motion segmentation map for
(a) and (b), (f) - combined ground-truth motion segmentation map of (b) and
(c), and (g) - predicted motion segmentation map of (b) and (c). The
segmentation model has learned to differentiate between motion and saturation,
and accurately segment only the moving regions with precise boundaries.

In Table IV, we show ablation results with using combinations of different
loss functions on the tonemapped HDR images.

TABLE IV: Loss function ablations. Refer Section IV-D for more details. Methods | PSNR-L | PSNR-T
---|---|---
$\ell_{2}$ | 41.47 | 42.03
$\ell_{1}$ | 41.01 | 41.25
$\ell_{2}$ \+ $\ell_{1}$ | 40.85 | 41.32
$\ell_{1}$ \+ MS-SSIM | 41.07 | 41.40
$\ell_{2}$ \+ MS-SSIM | 40.13 | 39.83
$\ell_{1}$ \+ $\ell_{2}$ \+ MS-SSIM | 40.61 | 40.63
Figure 7: Qualitative comparison between proposed and state-of-the art methods
for a hard validation example from [23] dataset. Figure 8: Qualitative results
by different methods for an example from METU dataset [45]. (a) Hu13 [20], (b)
Sen12 [21], (c) Prabhakar19 [26], (d) Yan19 [25], (e) Prabhakar20 [27], and
(f) Proposed method. The highlighted zoomed regions show that existing methods
have artifacts in both moving and saturated regions. Comparatively, our
proposed method generates a ghosting-free result.

###  IV-E Qualitative results

In Fig. 1, we compare proposed method against state-of-the-art HDR deghosting
methods for a difficult validation image from [23] dataset111Please refer to
supplementary material for additional results.. It is a challenging sequence,
as the moving arm occludes the reference image saturated regions in the low
exposure image. The traditional methods such as Sen12 and Hu13 were unable to
hallucinate correct textures in those regions. In those regions, the reference
image does not possess any concrete structural information, and additionally,
those regions are occluded by the moving arm in the low exposure image. Thus,
the patch-match based synthesis methods fail to find correct details and fill
those regions with false-color details (black zoom box in Fig. 1).
Comparatively, our proposed method can produce textures with plausible
details, even in such challenging regions. In Fig. 5, Sen12 method exhibits
similar false color artifacts in saturated regions as highlighted by the red
box. In addition to false color generation, Sen12 method exhibits a texture
smoothening artifact in the tile region (red dotted zoom box). Similarly, Hu13
method also produces false colors in the saturated regions (blue dotted zoom
box). In the second result of Fig. 5, as highlighted in blue and red zoom box,
Hu13 smoothens details in the cloud region. Similar artifacts by Sen12 and
Hu13 are observed in another example from METU dataset [45] in Fig. 8. In
contrast, the results generated by the proposed method are void of such
artifacts. Furthermore, with the help of memory network, it fills in plausible
details and generate natural looking results.

Figure 9: Our proposed method can easily be extended to fuse arbitrary length
sequence by replacing the feature concatenation with Mean+Max fusion strategy
followed by Prabhakar et al. [26]. Top row - the seven varying exposure images
captured with Canon EOS-5D Mark III camera in AEB setting with EV-3 to EV+3
bias. Results by (a) Hu13 [20], (b) Sen12 [21], (c) Prabhakar19 [26], (d)
Proposed method, and (e) Ground truth. The zoomed regions shown in the last
row are in linear domain to distinguish the artifacts easily. Figure 10:
Qualitative results to highlight the importance of segmentation maps. (a)
Cropped and zoomed region from the reference image, (b) Result by the A3 model
trained without segmentation maps, (c) Result by the proposed model using
segmentation maps, and (d) Ground truth. A3 model suffers from mild ghosting
artifacts prominent around the edges of moving regions.

Kalantari17 method trains a CNN to remove optical flow artifacts in the
aligned images. Despite this, some of their results still have flow artifacts.
As shown in the black zoomed box of Fig. 1, the optical flow artifacts in the
region affected by motion and saturation are left uncorrected by their model.
Similar artifact is observed in Fig. 5 and 7 as well. Prabhakar19 [26] method
also follows a similar optical flow error correction setup as that of
Kalantari17; however, it refines the aligned images using separate networks
before fusing them. Despite this, Prabhakar19 method still produces results
with visible artifacts in moving saturated regions. This phenomenon is visible
in their results in Fig. 1, 5 and 7. Both of these methods use optical flow
for alignment, followed by fusion with CNN. As evidenced by their results,
these two methods find it hard to handle difficult situations because the
network treats saturated regions as moving irrespective of motion.

Figure 11: Qualitative results to highlight the importance of memory network.
(a) Cropped and zoomed region from the reference image, (b) Result by the
model trained without memory network, (c) Result by the proposed model using
memory network, and (d) Ground truth. As highlighted, using memory network
aides the network to generate plausible textures in regions affected by
saturation and motion.

By comparison, our proposed method segments images into static and dynamic
before fusion, thus avoiding any such confusions.

Unlike Kalantari17 and Prabhakar19, Yan19 and Wu18 train a CNN model to fuse
unaligned images. While this approach aids to avoid false color generation or
optical flow artifacts, they present challenges of their own. In their method,
a single CNN is trained to handle motion, saturation, and fill in saturated
details. Hence, in heavily saturated static regions, mistaking them for moving
regions, they choose the reference image. It can be observed in the green zoom
box of Fig. 1. The window is unaffected by motion but saturated in the
reference image. Instead of obtaining sharper details from the low exposure
image, Yan19 and Wu18 generate blurred details from the reference image.
Similar artifacts are observed in the moving saturated regions (blue, magenta
and red boxes) of Fig. 1 and in Fig. 5, 7 and 8. Our proposed method is void
of such artifacts by employing different fusion sub-networks for static and
dynamic parts.

##  V Discussion

####  V-1 Arbitrary size fusion

The individual image features are concatenated for fusion in the proposed
approach. Hence, similar to [23, 24, 25], it requires retraining for a
different number of input images. We address this challenge by using the mean
+ max operation used by Prabhakar19. For $\mathbf{N}$ input images, we obtain
$\\{\mathbf{E}_{N}\\}$ feature maps using a shared encoder. As the encoder is
shared, it can be reused for an arbitrary number of images without retraining.
The features are then divided into dynamic and static using a single
segmentation network. Instead of concatenating the segmented features, we
concatenate the mean and max of all $\\{\mathbf{E}_{N}\\}$ features. The
fusion networks further process the concatenated dynamic and static features.

In Table II A7, we present the quantitative result for this approach on
Kalantari17 dataset. Compared with the Prabhakar19 approach, our proposed
method offers 1.4 dB improvement in PSNR-T, while offering the same
scalability advantage. In Fig. 9, we show example results by our approach with
seven input images and compare with state-of-the-art scalable methods like
Sen12, Hu13, and Prabhakar19. While Hu13 introduces false color in saturated
regions, Sen12 over smoothens the necessary details. Whereas, the proposed
method has successfully retained the structure with correct color details.

####  V-2 Moving region segmentation

The segmentation network can accurately detect moving regions across image
pairs. A few qualitative results are shown in Fig. 6. Quantitatively, the mean
intersection-over-union (IoU) achieved by our segmentation network with
respect to annotated masks provided by [27] is 0.686 on 15 validation examples
of the [23] dataset.

####  V-3 Running time

Our proposed method can process three input images of 1000$\times$1500
resolution in about 3.31 seconds on an NVIDIA Quadro RTX 6000 GPU.
Comparatively, the non-deep learning methods such as Hu13 and Sen12 methods
take around 450 and 320 seconds. Among the deep learning methods, Kalantari17
takes up to 57 seconds, and Wu18 takes 6.8 seconds. Unlike Kalantari17 method,
our proposed approach avoids the optical flow correction overhead by dealing
with dynamic and static regions separately.

##  VI Conclusion

In this paper, we proposed a motion segmentation assisted CNN for HDR
deghosting. With the help of predicted motion segmentation maps, the proposed
method fuses static and dynamic regions independently. Our method learns
different and dedicated fusion rules suited for static and dynamic regions by
identifying moving regions. Hence, it is better equipped to handle complex
motion and saturation seen in practice. The fused static and dynamic features
are combined to generate final high-quality ghost-free HDR results.
Furthermore, the proposed memory network accumulates the essential features
required to generate plausible textures in the saturated regions. Extensive
experimental evaluation exemplifies that the proposed method produces visually
appealing results with clear texture and faithful color reconstruction without
ghosting artifacts. Additionally, we hope that the presented MEDS dataset will
help the research community to advance further in a similar direction.

## References

  * [1] M. D. Tocci, C. Kiser, N. Tocci, and P. Sen, “A versatile HDR video production system,” in _ACM Transactions on Graphics (TOG)_ , vol. 30, no. 4. ACM, 2011, p. 41.
  * [2] M. McGuire, W. Matusik, H. Pfister, B. Chen, J. F. Hughes, and S. K. Nayar, “Optical splitting trees for high-precision monocular imaging,” _IEEE Computer Graphics and Applications_ , vol. 27, no. 2, pp. 32–42, 2007.
  * [3] H. Zhao, B. Shi, C. Fernandez-Cull, S.-K. Yeung, and R. Raskar, “Unbounded high dynamic range photography using a modulo camera,” in _IEEE International Conference on Computational Photography (ICCP)_. IEEE, 2015, pp. 1–10.
  * [4] P. E. Debevec and J. Malik, “Recovering high dynamic range radiance maps from photographs,” in _ACM SIGGRAPH 2008 classes_. ACM, 2008, p. 31.
  * [5] M. Granados, B. Ajdin, M. Wand, C. Theobalt, H.-P. Seidel, and H. P. Lensch, “Optimal HDR reconstruction with linear digital cameras,” in _Computer Vision and Pattern Recognition (CVPR)_. IEEE, 2010, pp. 215–222.
  * [6] K. R. Prabhakar, V. S. Srikar, and R. V. Babu, “Deepfuse: A deep unsupervised approach for exposure fusion with extreme exposure image pairs,” in _IEEE International Conference on Computer Vision (ICCV)_. IEEE, 2017, pp. 4724–4732.
  * [7] E. A. Khan, A. Akyiiz, and E. Reinhard, “Ghost removal in high dynamic range images,” in _IEEE International Conference on Image Processing_ , 2006.
  * [8] K. Jacobs, C. Loscos, and G. Ward, “Automatic high-dynamic range image generation for dynamic scenes,” _IEEE Computer Graphics and Applications_ , no. 2, pp. 84–93, 2008.
  * [9] F. Pece and J. Kautz, “Bitmap movement detection: HDR for dynamic scenes,” in _2010 Conference on Visual Media Production_. IEEE, 2010, pp. 1–8.
  * [10] W. Zhang and W.-K. Cham, “Reference-guided exposure fusion in dynamic scenes,” _Journal of Visual Communication and Image Representation_ , vol. 23, no. 3, pp. 467–475, 2012.
  * [11] T. Grosch, “Fast and robust high dynamic range image generation with camera and object movement,” _Vision, Modeling and Visualization, RWTH Aachen_ , pp. 277–284, 2006.
  * [12] O. Gallo, N. Gelfandz, W.-C. Chen, M. Tico, and K. Pulli, “Artifact-free high dynamic range imaging,” in _ICCP_. IEEE, 2009, pp. 1–7.
  * [13] S. Raman and S. Chaudhuri, “Reconstruction of high contrast images for dynamic scenes,” _The Visual Computer_ , vol. 27, pp. 1099–1114, 2011, 10.1007/s00371-011-0653-0. [Online]. Available: http://dx.doi.org/10.1007/s00371-011-0653-0
  * [14] L. Bogoni, “Extending dynamic range of monochrome and color images through fusion,” in _Proceedings. 15th International Conference on Pattern Recognition._ , 2000.
  * [15] S. B. Kang, M. Uyttendaele, S. Winder, and R. Szeliski, “High dynamic range video,” in _ACM Transactions on Graphics (TOG)_ , vol. 22, no. 3. ACM, 2003, pp. 319–325.
  * [16] T. Jinno and M. Okuda, “Motion blur free hdr image acquisition using multiple exposures,” in _2008 15th IEEE International Conference on Image Processing_. IEEE, 2008, pp. 1304–1307.
  * [17] H. Zimmer, A. Bruhn, and J. Weickert, “Freehand HDR imaging of moving scenes with simultaneous resolution enhancement,” in _Computer Graphics Forum_ , vol. 30, no. 2. Wiley Online Library, 2011, pp. 405–414.
  * [18] J. Hu, O. Gallo, and K. Pulli, “Exposure stacks of live scenes with hand-held cameras,” in _ECCV_. Springer, 2012, pp. 499–512.
  * [19] O. Gallo, A. Troccoli, J. Hu, K. Pulli, and J. Kautz, “Locally non-rigid registration for mobile HDR photography,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2015, pp. 49–56.
  * [20] J. Hu, O. Gallo, K. Pulli, and X. Sun, “HDR deghosting: How to deal with saturation?” in _IEEE Conference on Computer Vision and Pattern Recognition_ , 2013.
  * [21] P. Sen, N. K. Kalantari, M. Yaesoubi, S. Darabi, D. B. Goldman, and E. Shechtman, “Robust patch-based HDR reconstruction of dynamic scenes.” _ACM Trans. Graph._ , vol. 31, no. 6, p. 203, 2012.
  * [22] K. Ma, H. Li, H. Yong, Z. Wang, D. Meng, and L. Zhang, “Robust multi-exposure image fusion: A structural patch decomposition approach,” _IEEE Transactions on Image Processing_ , vol. 26, no. 5, pp. 2519–2532, 2017.
  * [23] N. K. Kalantari and R. Ramamoorthi, “Deep high dynamic range imaging of dynamic scenes,” _ACM Transactions on Graphics (Proceedings of SIGGRAPH 2017)_ , vol. 36, no. 4, 2017.
  * [24] S. Wu, J. Xu, Y.-W. Tai, and C.-K. Tang, “Deep high dynamic range imaging with large foreground motions,” in _European Conference on Computer Vision_ , 2018, pp. 120–135.
  * [25] Q. Yan, D. Gong, Q. Shi, A. van den Hengel, C. Shen, I. Reid, and Y. Zhang, “Attention-guided network for ghost-free high dynamic range imaging,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 1751–1760.
  * [26] K. R. Prabhakar, R. Arora, A. Swaminathan, K. P. Singh, and R. V. Babu, “A fast, scalable, and reliable deghosting method for extreme exposure fusion,” in _2019 IEEE International Conference on Computational Photography (ICCP)_. IEEE, 2019, pp. 1–8.
  * [27] K. R. Prabhakar, S. Agrawal, D. Singh, B. Ashwath, and R. V. Babu, “Towards practical and efficient high-resolution HDR deghosting with CNN,” in _European Conference on Computer Vision (ECCV)_ , 2020.
  * [28] Q. Yan, D. Gong, P. Zhang, Q. Shi, J. Sun, I. Reid, and Y. Zhang, “Multi-scale dense networks for deep high dynamic range imaging,” in _IEEE Winter Conference on Applications of Computer Vision (WACV)_. IEEE, 2019, pp. 41–50.
  * [29] S. Wu, S. Xie, S. Rahardja, and Z. Li, “A robust and fast anti-ghosting algorithm for high dynamic range imaging,” in _2010 IEEE International Conference on Image Processing_. IEEE, 2010, pp. 397–400.
  * [30] Y. S. Heo, K. M. Lee, S. U. Lee, Y. Moon, and J. Cha, “Ghost-free high dynamic range imaging,” in _Asian Conference on Computer Vision_. Springer, 2010, pp. 486–500.
  * [31] T.-H. Min, R.-H. Park, and S. Chang, “Histogram based ghost removal in high dynamic range images,” in _2009 IEEE International Conference on Multimedia and Expo_. IEEE, 2009, pp. 530–533.
  * [32] S. Raman and S. Chaudhuri, “Reconstruction of high contrast images for dynamic scenes,” _The Visual Computer_ , vol. 27, no. 12, pp. 1099–1114, 2011.
  * [33] A. Eden, M. Uyttendaele, and R. Szeliski, “Seamless image stitching of scenes with large motions and exposure differences,” in _Conference on Computer Vision and Pattern Recognition_. IEEE, 2006.
  * [34] E. Reinhard and K. Devlin, “Dynamic range reduction inspired by photoreceptor physiology,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 11, no. 1, pp. 13–24, 2005.
  * [35] G. Ward, “Fast, robust image registration for compositing high dynamic range photographs from hand-held exposures,” _Journal of Graphics Tools_ , vol. 8, no. 2, pp. 17–30, 2003.
  * [36] A. Tomaszewska and R. Mantiuk, “Image registration for multiexposure high dynamic range image acquisition,” in _Proceedings of the International Conference on Computer Graphics, Visualization and Computer Vision_ , 2007.
  * [37] C. Barnes, E. Shechtman, A. Finkelstein, and D. B. Goldman, “Patchmatch: A randomized correspondence algorithm for structural image editing,” _ACM Transactions on Graphics-TOG_ , vol. 28, no. 3, p. 24, 2009.
  * [38] Q. Yan, L. Zhang, Y. Liu, Y. Zhu, J. Sun, Q. Shi, and Y. Zhang, “Deep hdr imaging via a non-local network,” _IEEE Transactions on Image Processing_ , vol. 29, pp. 4308–4322, 2020.
  * [39] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [40] R. Schuster, O. Wasenmuller, C. Unger, and D. Stricker, “Sdc-stacked dilated convolution: A unified descriptor network for dense matching tasks,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 2556–2565.
  * [41] R. Deng, C. Shen, S. Liu, H. Wang, and X. Liu, “Learning to predict crisp boundaries,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018, pp. 562–578.
  * [42] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard _et al._ , “Tensorflow: A system for large-scale machine learning,” in _12th $\\{$USENIX$\\}$ symposium on operating systems design and implementation ($\\{$OSDI$\\}$ 16)_, 2016, pp. 265–283.
  * [43] X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” in _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , 2010, pp. 249–256.
  * [44] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [45] O. T. Tursun, A. O. Akyüz, A. Erdem, and E. Erdem, “An objective deghosting quality metric for HDR images,” in _Computer Graphics Forum_ , vol. 35, no. 2. Wiley Online Library, 2016, pp. 139–152.
  * [46] G. Krawczyk, K. Myszkowski, and H.-P. Seidel, “Lightness perception in tone reproduction for high dynamic range images,” in _Computer Graphics Forum_ , vol. 24, no. 3. Citeseer, 2005, pp. 635–646.
  * [47] F. Durand and J. Dorsey, “Fast bilateral filtering for the display of high-dynamic-range images,” in _ACM transactions on graphics (TOG)_ , vol. 21, no. 3. ACM, 2002, pp. 257–266.
  * [48] Y. Endo, Y. Kanamori, and J. Mitani, “Deep reverse tone mapping.” _ACM Trans. Graph._ , vol. 36, no. 6, pp. 177–1, 2017.
  * [49] G. Eilertsen, J. Kronander, G. Denes, R. K. Mantiuk, and J. Unger, “HDR image reconstruction from a single exposure using deep CNNs,” _ACM Transactions on Graphics (TOG)_ , vol. 36, no. 6, p. 178, 2017.
  * [50] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” _IEEE Transactions on Image Processing_ , vol. 13, no. 4, pp. 600–612, 2004.
  * [51] R. Mantiuk, K. J. Kim, A. G. Rempel, and W. Heidrich, “Hdr-vdp-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions,” _ACM Transactions on graphics (TOG)_ , vol. 30, no. 4, p. 40, 2011.
  * [52] E. Reinhard, M. Stark, P. Shirley, and J. Ferwerda, “Photographic tone reproduction for digital images,” in _ACM Transactions on Graphics (TOG)_ , vol. 21, no. 3. ACM, 2002, pp. 267–276.
  * [53] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016.