1
Simultaneous Video Stabilization and Moving
Object Detection in Turbulence
Omar Oreifej, Member, IEEE, Xin Li, Member, IEEE, and Mubarak Shah, Fellow, IEEE
Abstract—Turbulence mitigation refers to the stabilization of videos with non-uniform deformations due to the influence of optical
turbulence. Typical approaches for turbulence mitigation follow averaging or de-warping techniques. Although these methods
can reduce the turbulence, they distort the independently moving objects which can often be of great interest. In this paper,
we address the novel problem of simultaneous turbulence mitigation and moving object detection. We propose a novel three-
term low-rank matrix decomposition approach in which we decompose the turbulence sequence into three components: the
background, the turbulence, and the object. We simplify this extremely difficult problem into a minimization of nuclear norm,
Frobenius norm, and `
1
norm. Our method is based on two observations: First, the turbulence causes dense and Gaussian
noise, and therefore can be captured by Frobenius norm, while the moving objects are sparse and thus can be captured by `
1
norm. Second, since the object’s motion is linear and intrinsically different than the Gaussian-like turbulence, a Gaussian-based
turbulence model can be employed to enforce an additional constraint on the search space of the minimization. We demonstrate
the robustness of our approach on challenging sequences which are significantly distorted with atmospheric turbulence and
include extremely tiny moving objects.
Index Terms—Three-Term Decomposition, Turbulence Mitigation, Rank Optimization, Moving Object Detection, Particle Advec-
tion, Restoring Force.
F
1 INTRODUCTION
T
HE refraction index of the air varies based on
several atmospheric characteristics including the
air’s temperature, humidity, pressure, carbon dioxide
level, and dust density. Such conditions are typically
not homogeneous; for instance, a non-uniform tem-
perature distribution might be observed above a sur-
face receiving sunlight. Therefore, light rays travelling
through the air with such non-uniform changes in its
relative refraction index, will go through a complex
series of refraction and reflection causing extreme
spatially and temporally varying deformations to the
captured images [1], [2], [3], [4], [5].
On the other hand, if the objects of interest are
additionally moving in the scene, their motion will
be mixed up with the turbulence deformation in the
captured images, rendering the problem of detecting
the moving objects extremely difficult. In this paper,
we are interested in the dual problem of turbulence
mitigation (stabilizing the sequence) and moving ob-
ject detection under the turbulent medium. To the
best of our knowledge, such a problem has never
been explored before. Relevant previous approaches
• O. Oreifej, and Mubarak Shah are with the Department of Electrical
Engineering and Computer Science at University of Central
Florida, 4000 Central Florida Blvd., Orlando, FL 32816. E-mail:
oreifej@eecs.ucf.edu, shah@eecs.ucf.edu.
• Xin Li is with the Math Department at University of Central
Florida, 4000 Central Florida Blvd., Orlando, FL 32816. E-mail:
xli@math.ucf.edu.
have either focused on detecting moving objects or
de-warping a deformed sequence, but not on both
tasks concurrently. Note that other than image defor-
mation, atmospheric turbulence may cause blur if the
camera exposure time is not sufficiently short. In this
paper, however, we focus only on image deformation
because of the inherent confusion between the motion
of the object and the motion caused by the turbulence.
Given a sequence of frames {I
1
, ..., I
T
} acquired
from a stationary camera residing in a turbulent
medium while observing relatively tiny moving ob-
jects, we decompose the sequence into background,
turbulence, and object components. More precisely,
consider the frames matrix F = [vec{I
1
} · · · vec{I
T
}]
for I
k
∈ R
W ×H
(k = 1, 2, ..., T ), where W × H
denotes the frame resolution (width by height), and
vec : R
W ×H
→ R
M
is the operator which stacks the
image pixels as a column vector. We formulate our
decomposition of F as:
min
A,O,E
Rank(A) s.t. F = A + O + E, (1)
||O||
0
≤ s, ||E||
F
≤ σ,
where F, A, O, and E are the matrices of frames, back-
ground, object, and error (turbulence), respectively.
Here, the k · k
0
−norm counts the number of nonzero
entries, k · k
F
−norm is the Frobenius norm which
is equal to the square root of the sum of squared
elements in the matrix, s represents an upper bound
of the total number of moving objects’ pixels across
all images, and σ is a constant which reflects our
knowledge of the maximum total variance due to
2
Frames
Intensity
Model
Pre-
Processing
Processed
Frames
I
1
..I
n
Turbulence Noise Model
Three-Term
Low-Rank
Opmizaon
Object
Confidence
Map
Background
Object
Turbulence
Moon
Model
Fig. 1: The various steps of the proposed algorithm.
corrupted pixels across all images.
Our decomposition is based on the intrinsic prop-
erties of each of the components:
1) The Background: The scene in the background is
presumably static; thus, the corresponding com-
ponent in the frames of the sequence has linearly
correlated elements. Therefore, the background
component is expected to be the part of the
matrix which is of low rank. Minimizing the
rank of the low-rank component of the frames
matrix F emphasizes the structure of the linear
subspace containing the column space of the
background, which reveals the background.
2) The Turbulence: Previous work dealing with
turbulence, such as [6], [7], [8], [9], [10], demon-
strated that the fluctuations of fluids (for in-
stance air and water) attain Gaussian-like char-
acteristics such as being unimodal, symmetric,
and locally repetitive; therefore, the projected
deformations in the captured sequence often
approach a Gaussian distribution (we discuss
this in more detail in Section 3.2). For this reason,
the turbulence component can be captured by
minimizing its Frobenius norm. The Frobenius
norm of a matrix is the same as the Euclidean
norm of the vector obtained from the matrix by
stacking its columns. Therefore, as in the well-
known vector case, constraining the error in the
Euclidean norm is equivalent to controlling the
sample variance of the error. Furthermore, theo-
retically, the estimate obtained by the Frobenius
norm has several desirable statistical properties
[11].
3) The Moving Objects: We assume that the moving
objects are sparse in the sequence. This means
that the number of pixels occupied by the mov-
ing objects is small (or can be considered as
outliers) compared to the total number of pixels
in the frames. This is a reasonable assumption
for most realistic surveillance videos. For this
reason, the moving objects are best captured
by restricting the number of nonzero entries
(denoted by the `
0
norm of the matrix), which
is desirable for finding outliers.
In practice, parts of the turbulence could also ap-
pear as sparse errors in the object matrix O. There-
fore, an additional constraint needs to be enforced on
the moving objects. We employ a simple turbulence
model to compute an object confidence map which is
used to encourage the sparse solutions to be located
on regions exhibiting linear motion that is dissimilar
from the fluctuations of the turbulence. Under the
new constraint, the optimization problem (1) must be
reformulated as:
min
A,O,E
Rank(A) s.t. F = A + O + E, (2)
||Π(O)||
0
≤ s, ||E||
F
≤ σ,
where Π : R
M×T
→ R
M×T
is the object confidence
map, which is a linear operator that weights the
entries of O according to their confidence of cor-
responding to a moving object such that the most
probable elements are unchanged and the least are
set to zero.
Figure 1 shows a diagram of the proposed ap-
proach. We first apply a pre-precessing step to im-
prove the contrast of the sequence, and reduce the
spurious and random noise. Consequently, we obtain
an object confidence map using a turbulence model
which utilizes both the intensity and the motion cues.
Finally, we decompose the sequence into its compo-
nents using three-term rank minimization.
This paper makes three main contributions: First,
we propose a new variant of matrix decomposition
based on low-rank optimization and employ it to
solve the novel problem of simultaneous moving
object detection and turbulence mitigation in videos
3
distorted by atmospheric turbulence. Second, we pro-
pose a turbulence model based on both intensity and
motion cues, where the motion distribution is derived
from the Lagrangian particle advection framework
[12]. The turbulence model is used to enforce an
additional constraint on the decomposition which en-
courages the sparse solutions to be located in areas
with non-Gaussian motion. Finally, we propose an
additional force component in the particle advection
framework in order to stabilize the particles in the
turbulent medium and handle long sequences without
discontinuities.
The rest of the paper is organized as follows: In the
next section, we discuss the related works. In Section
3, we present our three-term rank minimization, fol-
lowed by our probabilistic formulation of the object
confidence map using particle advection. Section 4
discusses the details of implementation and the tech-
nical challenges of the three-term optimization. The
experiments and the results are described in Section
5. Finally, Section 6 concludes the paper.
2 RELATED WORK
Rank optimization-based video de-noising has re-
cently flourished with several successful works re-
ported, from which we will only discuss the most
related articles. Robust PCA was proposed in [13],
where a low rank matrix was recovered from a small
set of corrupted observations through convex pro-
gramming. Similar concepts were later employed in
[14] for video de-noising, where serious mixed noise
was extracted by grouping similar patches in both
spatial and temporal domains, and solving a low-
rank matrix completion problem. Additionally, in [15],
linear rank optimization was employed to align faces
with rigid transformations, and concurrently detect
noise and occlusions. In [16], Yu et. al. proposed
an efficient solution to subspace clustering problems
which involved the optimization of unitarily invariant
norms. Another variant of such space-time optimiza-
tion techniques is the total variation minimization,
where for instance in [17], Chan et. al. posed the
problem of video restoration as a minimization of
anisotropic total variation given in terms of `
1
-norm,
or isotropic variation given in terms of `
2
-norm. Con-
sequently, the Lagrange multiplier method was used
to solve the optimization function.
On the other hand, moving object detection is a
widely investigated problem. When the scene is static,
moving objects can be easily detected using frame
differencing. A better approach would be to use the
mean, the median, or the running average as the
background [18]. The so-called eigenbackground [19]
can also be obtained using PCA. However, when the
scene is constantly changing because of noise, light
changes, or camera shake, the intensities of image
pixels can be considered as independent random vari-
ables, which can be represented using a statistical
model such as a Gaussian, a mixture of Gaussians, or a
kernel density estimator. The model can then be used
to compute the probability for each pixel to belong to
either the background or the foreground. Examples
of such approaches include [20], [21], [22]. Addition-
ally, the correlation between spatially proximal pixels
could also be employed to improve the background
modelling using a joint domain (location) and range
(intensity) representation of image pixels such as in
[23].
Approaches for turbulence mitigation focused
mainly on registration-based techniques. In [6], [5],
[24], both the turbulence deformation parameters and
a super-resolution image were recovered using area-
based B-Spline registration. Moreover, in [25], Tian
and Narasimhan proposed recovering the large non-
rigid turbulence distortions through a “pull-back”
operation that utilizes several images with known
deformations. In model-based tracking [26], the char-
acteristics of the turbulence caused by water waves
were employed to estimate the water basis using PCA.
More recently, in [7], turbulence caused by water was
overcome by iteratively registering the sequence to
its mean followed by RPCA to extract the sparse er-
rors. Averaging-based techniques are also popular for
video de-noising and turbulence mitigation, includ-
ing pixel-wise mean/median, non-local means (NLM)
[27], [28], fourier-based averaging [29], and speckle
imaging [1], [3], [30]. Another category of methods
for turbulence mitigation is the lucky region approach
[4], [31], [32], where the least distorted patches of the
video are selected based on several quality statistics,
then those selected patches are fused together to
compose the recovered video.
Clearly, previous work in moving object detection
in dynamic scenes mostly focused on detecting the ob-
jects and did not consider recovering the background.
Inversely, previous work in turbulence mitigation did
not consider the possible interest in detecting moving
objects in the scene. In this paper, we pose the two
problems of moving object detection and turbulence
mitigation as one application for our proposed three-
term low-rank decomposition. We demonstrate how
to decompose a turbulent video into separate back-
ground, foreground, and turbulence components. It
is important to note that our method is not directly
comparable to background subtraction or turbulence
mitigation approaches; though, we do provide com-
petitive results on each task separately.
3 PROPOSED APPROACH
We decompose the matrix which contains the frames
of the turbulence video, into its components: the back-
ground, the turbulence, and the objects. The decompo-
sition is performed by solving the rank optimization
in equation (2), which enforces relevant constraints on
each component. In the next subsection we describe
the details of the decomposition approach.
4
3.1 Three-Term Decomposition
When solving equation (2), it is more convenient to
consider the Lagrange form of the problem:
min
A,O,E
Rank(A) + τ ||Π(O)||
0
+ λ||E||
2
F
(3)
s.t. F = A + O + E,
where τ and λ are weighting parameters. The opti-
mization of (3) is not directly tractable since the matrix
rank and the `
0
-norm are nonconvex and extremely
difficult to optimize. However, it was recently shown
in [13] that when recovering low-rank matrices from
sparse errors, if the rank of the matrix A to be
recovered is not too high and the number of non-
zero entries in O is not too large, then minimizing the
nuclear norm of A (sum of singular values
P
i
σ
i
(A))
and the `
1
-norm of O can recover the exact matrices.
Therefore, the nuclear norm and the `
1
-norm are the
natural convex surrogates for the rank function and
the `
0
-norm, respectively. Applying this relaxation,
our new optimization becomes:
min
A,O,E
||A||
∗
+ τ||Π(O)||
1
+ λ||E||
2
F
s.t. F = A + O + E,
(4)
where ||A||
∗
denotes the nuclear norm of matrix A. We
adopt the Augmented Lagrange Multiplier method
(ALM) [33] to solve the optimization problem (4). De-
fine the augmented Lagrange function for the problem
as:
L(A, O, E, Y ) = ||A||
∗
+ τ ||Π(O)||
1
+ λ||E||
2
F
+ (5)
hY, F − A − O − Ei +
β
2
||F − A − O − E||
2
F
,
where Y ∈ R
M×T
is a Lagrange multiplier matrix,
β is a positive scalar, and h, i denotes the matrix in-
ner product (trace(A
T
B)). Minimizing the function in
equation (5) can be used to solve the constrained op-
timization problem in equation (4). We use the ALM
algorithm to iteratively estimate both the Lagrange
multiplier and the optimal solution by iteratively
minimizing the augmented Lagrangian function:
(A
k+1
, O
k+1
, E
k+1
) = arg min
A,O,E
L(A, O, E, Y
k
), (6)
Y
k+1
= Y
k
+ β
k
(F
k+1
− A
k+1
− O
k+1
− E
k+1
).
When β
k
is a monotonically increasing positive se-
quence, the iterations converge to the optimal solution
of problem (4) [34]. However, solving equation (6)
directly is difficult; therefore, the solution is approxi-
mated using an alternating strategy minimizing the
augmented Lagrange function with respect to each
component separately:
A
k+1
= arg min
A
L(A, O
k
, E
k
, Y
k
), (7)
O
k+1
= arg min
O
L(A
k+1
, O, E
k
, Y
k
),
E
k+1
= arg min
E
L(A
k+1
, O
k+1
, E, Y
k
).
Following the idea of the singular value threshold-
ing algorithm [35], we derive the solutions for the
update steps in equation (7) for each of the nuclear,
Frobenius, and `
1
norms. Please refer to the sup-
plementary appendix for the complete derivations.
Consequently, a closed form solution for each of the
minimization problems is found:
UW V
T
= svd(F − O
k
− E
k
+ β
−1
k
Y
k
), (8)
A
k+1
= US
1/β
k
(W )V
T
,
O
k+1
= S
τ/β
k
Π
(F − A
k+1
− E
k
+ β
−1
k
Y
k
),
E
k+1
= (1 +
2λ
β
k
)
−1
(β
−1
k
Y
k
+ F − A
k+1
− O
k+1
),
where svd(M) denotes a full singular value decompo-
sition of matrix M , and S
α
(·) is the soft-thresholding
operator defined for a scalar x as:
S
α
(x) = sign(x) · max{|x| − α, 0}, (9)
and for two matrices A = (a
ij
) and B = (b
ij
) of the
same size, S
A
(B) applies the soft-thresholding entry-
wise outputting a matrix with entries S
a
ij
(b
ij
) .
The steps of our decomposition are summarized
in Algorithm 1. In the next subsection, we describe
our method to obtain the moving object confidence
map Π, which is employed as a prior in the rank
minimization problem.
3.2 Turbulence Model
We employ a turbulence model to enforce an ad-
ditional constraint on the rank minimization such
that moving objects are encouraged to be detected
in locations with non-Gaussian deformations. Exact
modelling of the turbulence is in fact ill-posed as
it follows a non-uniform distribution which varies
significantly in time, besides having an additional
complexity introduced during the imaging process;
thus, rendering the problem of modelling turbulence
extremely difficult. Although the refraction index of
the turbulent medium is often randomly changing,
it is also statistically stationary [3], [30], [36]; thus,
the deformations caused by turbulence are generally
repetitive and locally centered [6], [7], [8], [9], [37];
this encourages the use of Gaussian-based models
as approximate distributions that are general enough
to avoid overfitting, but rather capture significant
portion of the turbulent characteristics.
We use a Gaussian function to model the inten-
sity distribution of a pixel going through turbulence.
This is similar to [20] which employs a mixture of
5
Gaussians; however, we found that a single Gaussian
worked better since more complicated models often
require a period of training which is not available in
our sequences. Therefore, the intensity of a pixel at
location x is modelled using a Gaussian distribution:
I(x) ∼ N (µ
I
, σ
I
), (10)
where µ
I
and σ
I
are the mean and the standard
deviation at x, respectively. On the other hand, the
deformation caused by turbulence can be captured in
the motion domain besides the intensity. Therefore,
we combine the intensity and the motion features
to obtain a better model of turbulence. In order to
capture the ensemble motion in the scene, we use
the concept of a “particle” in a Lagrangian particle
trajectory acquisition approach. We assume that a
grid of particles is overlaid onto a scene where each
particle corresponds to a single pixel (the granularity
is controllable). The basic idea is to quantify the
scene’s motion in terms of the motion of the particles
which are driven by dense optical flow. A so-called
particle advection [12], [38], [39] procedure is applied
to produce the particle trajectories. Given a video
clip ∈ R
W ×H×T
, we denote the corresponding optical
flow by (U
t
w
, V
t
h
), where w ∈ [1, W ], h ∈ [1, H], and
t ∈ [1, T − 1]. The position vector (x
t
w
, y
t
h
) of the
particle at grid point (w, h) at time t is estimated by
solving the following differential equations:
dx
t
w
dt
= U
t
w
, (11)
dy
t
h
dt
= V
t
h
.
We use Euler’s method to solve them, similar to
[38]. By performing advection for the particles at all
grid points with respect to each frame of the clip,
we obtain the clip’s particle trajectory set, denoted by
{(x
t
w
, y
t
h
)|w ∈ [1, W ], h ∈ [1, H], t ∈ [1, T ]}.
We employ the spatial locations of the particle tra-
jectories (i.e. (x
t
w
, y
t
h
)) to model the turbulence motion
in the scene. The locations visited by a particle moving
due to the fluctuations of the turbulence have a uni-
modal and symmetric distribution which approaches
a Gaussian [6], [30], [10], [9]. This is dissimilar from
the linear motion of the particles driven by moving
objects. Therefore, we associate each particle with a
Gaussian with mean µ
M
and covariance matrix Σ
M
:
x ∼ N (µ
M
, Σ
M
). (12)
By augmenting the intensity model in equation (10)
with the motion model in equation (12), the total
confidence of corresponding to the turbulence versus
the moving objects for a particle at location x is
expressed as a linear opinion pooling of the motion
and the intensity cues
C(x) = wP(I(x)|µ
I
, σ
I
) + (1 − w)P(x|µ
M
, Σ
M
). (13)
The parameters of our model {w, µ
I
, σ
I
, µ
M
, Σ
M
}
can be learned by optimization using training se-
quences or set to constant values selected empirically.
In the context of our three-term decomposition, the
obtained confidence provides a rough prior knowl-
edge of the moving objects’ locations, which can be
incorporated into the matrix optimization problem in
equation (4). Interestingly, this prior employs motion
information; therefore, it is complementary to the
intensity-based rank optimization, and can signifi-
cantly improve the result.
At frame t, we evaluate all the particles’ locations
against their corresponding turbulence models and
obtain the turbulence confidence map C
t
∈ R
W ×H
.
While C
t
corresponds to the confidence of a particle
to belong to turbulence, the desired Π in equation
(4) corresponds to the confidence of belonging to
the moving objects; therefore, we define the object
confidence map Π as the complement of the stacked
turbulence confidence maps:
Π = 1 − [vec{C
1
} · · · vec{C
T
}]. (14)
3.3 Restoring Force
The particles carrying the object’s motion typically
drift far from their original locations leaving several
gaps in the sequence. In the presence of turbulence,
the drifting also occurs as a result of the turbulent
motion. Therefore, the particles need to be reinitial-
ized every certain number of frames which, however,
creates discontinuities. This is a typical hurdle in the
Lagrangian framework of fluid dynamics [12], [38],
[39], which constitutes a major impediment for the
application of particle flow to turbulence videos. In
order to handle the drifting and the discontinuity
problems associated with the particle flow, we use a
new force component in the advection equation:
dx
t
w
dt
= U
t
w
+ G(x, x
o
), (15)
dy
t
h
dt
= V
t
h
+ G(y, y
o
).
We refer to the new force as “Restoring Force” -
a reference to a local restoration force acting in the
direction of the original location of each particle. We
use a simple linear function to represent the restoring
force:
G(x, x
o
) =
x − x
o
s
, (16)
where s is a scaling factor which trades off the de-
tection sensitivity and the speed of recovery for the
particles. In other words, if s is set to a high value,
the effect of the restoring force will be negligible,
and therefore the particles will require a relatively
longer time to return to their original positions. In this
case, the sensitivity of moving object detection will be
![Fig. 6: Example frames illustrating the contribution of the turbulence model and the low-rank optimization in the total moving object detection performance, separately. The first row shows the frames, the second row shows the object confidence map obtained from the turbulence model (confidence values are mapped to [0−255], the highest confidence is black), and the third row shows the final object blob after the three-term low-rank decomposition. Clearly, the turbulence model provides a rough estimation of the object location, while the low-rank optimization refines the result to obtain an accurate detection.](/figures/fig-6-example-frames-illustrating-the-contribution-of-the-2cnqpo9p.webp)
![Fig. 5: Our moving object detection results on sample frames from sequences 3 and 4 compared to [22], [19], and [23]. For every sequence, the first row shows the frames, the second row shows the result from our method (taken from matrix O after the decomposition), the third, forth, and fifth rows show the background subtraction result obtained using [22], [19], and [23], respectively. Please zoom in to see the details.](/figures/fig-5-our-moving-object-detection-results-on-sample-frames-3mvdmpk7.webp)
