Sajid Javed

Arif Mahmood

Thierry Bouwmans

Soon Ki

Background-Foreground Modeling Based on Spatiotemporal Sparse Subspace Clustering

Keywords: Background modeling, Foreground Detection, Graph Regularization, Subspace Clustering, Robust Principal Component Analysis

Background estimation and foreground segmentation are important steps in many high-level vision tasks. Many existing methods estimate background as a low-rank component and foreground as a sparse matrix without incorporating the structural information. Therefore, these algorithms exhibit degraded performance in the presence of dynamic backgrounds, photometric variations, jitter, shadows, and large occlusions. We observe that these backgrounds often span multiple manifolds. Therefore, constraints that ensure continuity on those manifolds will result in better background estimation. Hence, we propose to incorporate the spatial and temporal sparse subspace clustering into the RPCA framework. To that end, we compute a spatial and temporal graph for a given sequence using motion-aware correlation coefficient. The information captured by both graphs is utilized by estimating the proximity matrices using both the normalized Euclidean and geodesic distances. The low-rank component must be able to efficiently partition the spatiotemporal graphs using these Laplacian matrices. Embedded with the RPCA objective function, these Laplacian matrices constrain the background model to be spatially and temporally consistent, both on linear and nonlinear manifolds. The solution of the proposed objective function is computed by using the LADMAP optimization scheme. Experiments are performed on challenging sequences from five publicly available datasets and are compared with 23 existing state-of-the-art methods. The results demonstrate excellent performance of the proposed algorithm for both background estimation and foreground segmentation.

I. INTRODUCTION

B ACKGROUND estimation and foreground segmentation originate in numerous applications in computer vision including moving object detection [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF], video surveillance [START_REF] Bouwmans | Robust PCA via principal component pursuit: a review for a comparative evaluation in video surveillance[END_REF], visual object tracking [START_REF] Zhang | Robust visual tracking via consistent low-rank sparse learning[END_REF], and salient motion detection [START_REF] Gao | Block-sparse RPCA for salient motion detection[END_REF]. Background modeling is mainly intended to efficiently and accurately extract a model which describes the scene without foreground objects. On the other hand, foreground detection is intended for segmenting moving objects from the known background model [START_REF] Bouwmans | Scene background initialization: A taxonomy[END_REF], [START_REF] Maddalena | Towards benchmarking scene background initialization[END_REF], [START_REF] Oreifej | Simultaneous video stabilization and moving object detection in turbulence[END_REF]. Both of these methods become challenging in the presence of dynamic background, changing lighting conditions, and jitter induced by the sensor. Background modeling also suffers in the presence of occlusion because of foreground objects. A number of techniques have been proposed in literature that mostly address relatively S. Javed and S. K. Jung are with Virtual Reality Laboratory, the School of Computer Science and Engineering, Kyungpook National University, 80 Daehak-ro, Buk-gu, Daegu, 41566, Republic of Korea. (e-mail: sajid@vr.knu.ac.kr, skjung@knu.ac.kr). (corresponding author: Soon Ki Jung) A. Mahmood is with Department of Computer Science and Engineering, Qatar University, Qatar. (email: arif.mahmood@qu.edu.qa) T. Bouwmans is with Universite de La Rochelle, Laboratoire Mathematiques, Image et Applications (MIA), [START_REF] Johnson | Efficient algorithms for shortest paths in sparse networks[END_REF] Avenue Albert Einstein, La Rochelle, BP 33060-17301, France. (email: tbouwman@univ-lr.fr)

Sequence of input images (a) (b) (c) (d) (e)
Background model estimation using RPCA methods [START_REF] Oreifej | Simultaneous video stabilization and moving object detection in turbulence[END_REF], [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF] Background model estimation using Proposed MSCL algorithm simple scenarios for the estimation of background and/or foreground component [START_REF] Bouwmans | Scene background initialization: A taxonomy[END_REF].

As an example, most methods perform well when the pixels or regions of the background are visible for longer periods in the training data. Furthermore, the performance of many methods [START_REF] Diego | Rejection based multipath reconstruction for background estimation in video sequences with stationary objects[END_REF], [START_REF] Liu | Background subtraction based on low-rank and structured sparse decomposition[END_REF], [START_REF] Ye | Foreground-Background Separation From Video Clips via Motion-Assisted Matrix Restoration[END_REF] degrades significantly in real-time scenarios, especially if the background scene is visible for short periods. Moreover, the performance of foreground detection methods is also effected by the dynamic pixels, such as those caused by the rippling of water surface, swaying of trees, and sudden variations in lighting conditions. Recently, Robust Principal Component Analysis (RPCA) has been shown to be an efficient framework for segregating background and foreground components [START_REF] Candès | Robust principal component analysis?[END_REF], [START_REF] Wright | Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[END_REF]. In RPCA, Wright et al. [START_REF] Wright | Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[END_REF] considered background-foreground modeling as a matrix decomposition problem:

B,F min||B|| * +λ 1 ||F|| 1 such that X = B + F, (1)
where X = [x 1 , x 2 , ..., x n ] ∈ R p×n is the input video sequence of n frames, and each x i ∈ R p denotes i-th frame. The low-rank component B corresponds to the background model whereas the locally deforming regions called foreground constitute the sparse matrix F. This decomposition is achieved by solving a convex optimization method [START_REF] Candès | Robust principal component analysis?[END_REF]. In [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], ||•|| * denotes the nuclear norm, which is defined as the sum of all singular values, and ||•|| 1 is the l 1 -norm. Nuclear norm enforces the background images to be linearly correlated, and the l 1 -norm constrains the foreground matrix to be sparse.

In many real scenarios, traditional RPCA methods [START_REF] Diego | Rejection based multipath reconstruction for background estimation in video sequences with stationary objects[END_REF], [START_REF] Oreifej | Simultaneous video stabilization and moving object detection in turbulence[END_REF], [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF] suffer from some obvious limitations. For example, many real-world scenarios contain redundant information in the form of motionless frames, where foreground objects remain static for long durations and then start moving. If the redundant information is not properly handled, the outlier regions will also appear in the estimated background (Fig. 1 (a)). Furthermore, sudden photometric variations of the background increase its rank; therefore, the low-rank background model estimation fails (Fig. 1 (b)). If both background and foreground coexist in each frame, the low-rank component cannot efficiently capture the background model because of abrupt changes (Fig. 1 (c)-(d)). In the case of some dynamic background sequences, the low-rank component fails to estimate the background model (Fig. 1 (e)). Given these common real-world scenarios, the existing methods fail to achieve good quality background-foreground models because of unrealistic assumptions.

We observe that the real-world background sequences, as shown in Fig. 1, may span one or more linear or nonlinear manifolds. To investigate this fact, we applied sparse subspace clustering on multiple dynamic background scenes. Nonlinear shapes of the clusters in 3D principal-space of geodesic distance may be observed in Fig. 2. In high dimensional space, the frames in each cluster span different nonlinear manifolds. Therefore, constraining the estimated low-rank B to remain consistent with the subspace structure of the original background significantly improves the background modeling performance.

Compared to the traditional RPCA techniques, following modifications are proposed in this study:

• Two graphs are constructed, namely, a temporal graph over the background frames and a spatial graph over the spatial background locations using motion-aware Correlation Coefficient (CC) as a similarity measure to handle the background modeling efficiently. • For background models spanning linear subspaces, the normalized Euclidean distance based proximity matrices are directly computed from these graphs. Sparse coding is used to linearly decompose each column of each proximity matrix using the remaining columns as dictionary items. The linear coefficients are used to recalculate a new proximity matrix which is then used to compute the sparse subspace based normalized Laplacian matrix [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF]. • For backgrounds spanning nonlinear manifolds, Geodesic Distance (GD) based proximity matrices are computed for both graphs. For each matrix, we apply sparse coding and compute the geodesic subspace based Laplacian matrix [START_REF] Mahmood | Semi-supervised spectral clustering for image set classification[END_REF]. The normalized Laplacian matrices computed over both Euclidean and Geodesic distances are embedded into the RPCA framework. These graph-based regularization's help us to overcome the limitations of the existing methods [START_REF] Candès | Robust principal component analysis?[END_REF], [START_REF] Liu | Background subtraction based on low-rank and structured sparse decomposition[END_REF], [START_REF] Oreifej | Simultaneous video stabilization and moving object detection in turbulence[END_REF], [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF]. To the best of our knowledge, this is the first study that integrates the spatiotemporal clustering information into low-rank component for improved background modeling.

We name the resulting algorithm as Motion-assisted Spatiotemporal Clustering of Low-rank (MSCL). In the proposed algorithm, we first detect a set of dynamic frames D within the sequence X by eradicating motionless frames, which help us in dealing with redundant data problem, as shown in Fig. 1 (a). For this purpose, we estimate the optical flow between the consecutive frames and use it to generate a binary motion mask, which assists us to obtain approximate knowledge on the pixels of background and foreground components. The spatial and temporal graphs are computed using motion-aware CC [START_REF] Mahmood | Exploiting transitivity of correlation for fast template matching[END_REF], which includes pixels exhibiting motion less than a predefined threshold. The use of CC-based regularization enforces the background model to be invariant to lighting condition variations and to enhance continuity on linear manifolds, while GD-based regularization ensures the continuity of the background model on nonlinear manifolds.

We solve our proposed objective function using the Linearized Alternating Direction Method with Adaptive Penalty (LADMAP) optimization because of its efficiency in using less auxiliary variables with convergence guarantee [START_REF] Lin | Linearized alternating direction method with adaptive penalty for low-rank representation[END_REF], [START_REF] Liu | Linearized alternating direction method with parallel splitting and adaptive penalty for separable convex programs in machine learning[END_REF]. We evaluated our proposed algorithm on a new dataset Scene Background Modeling.net (SBM.net) [START_REF] Pierre-Marc | Scene Background Modeling.net[END_REF], which contains 8 challenging categories and 80 videos. We compared the proposed MSCL algorithm with several existing state-of-theart methods. Our proposed algorithm outperforms these earlier methods on almost all categories of the SBM.net dataset 1 .

Compared to existing methods, our algorithm generates better foreground model; however, we do not explicitly constrain the foreground, which may result in a noisy and dis-contiguous

(a) (b) (c) (d) (e)
Highly dynamic input images

Ground truth images

Foreground estimation by using traditional RPCA [START_REF] Oreifej | Simultaneous video stabilization and moving object detection in turbulence[END_REF] Foreground estimation by using our proposed MSCL algorithm Foreground estimation by using our proposed MSCL-FL algorithm foreground (Fig. 3). Figs. 3 (a)-(e) show that although our proposed algorithm estimated an outlier-free background model in the case of highly dynamic background scenes, the foreground pixels are not contiguous. Therefore, we propose to enforce a post-processing step using Fused Lasso (FL) [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF], [START_REF] Xin | Background subtraction via generalized fused lasso foreground modeling[END_REF], [START_REF] Xue | Foreground estimation based on linear regression model with fused sparsity on outliers[END_REF] to strengthen the similarity among the neighboring foreground pixels. FL penalizes the pixels not only by l 1norm, but it also preserves the structure of the neighboring pixels. We name this extended algorithm as MSCL-FL for the detection of foreground objects. The proposed MSCL-FL algorithm has shown better performance compared to stateof-the-art methods including TVRPCA [START_REF] Cao | Total Variation Regularized RPCA for Irregularly Moving Object Detection Under Dynamic Background[END_REF], 2P-RPCA [START_REF] Gao | Block-sparse RPCA for salient motion detection[END_REF], DP-GMM [START_REF] Haines | Background subtraction with DirichletProcess mixture models[END_REF], GRASTA [START_REF] He | Incremental gradient on the grassmannian for online foreground and background separation in subsampled video[END_REF], TLSFSD [START_REF] Hu | Moving object detection using tensor-based low-rank and saliently fused-sparse decomposition[END_REF], SRPCA [START_REF] Sajid | Spatiotemporal low-rank modeling for complex scene background initialization[END_REF], LSD [START_REF] Liu | Background subtraction based on low-rank and structured sparse decomposition[END_REF], 3TD [START_REF] Oreifej | Simultaneous video stabilization and moving object detection in turbulence[END_REF], MODSM [START_REF] Pang | Incremental learning with saliency map for moving object detection[END_REF], MLRBS [START_REF] Qin | A background basis selectionbased foreground detection method[END_REF], PAWCS [START_REF] St-Charles | Universal Background Subtraction Using Word Consensus Models[END_REF], SuBSENSE [START_REF] St-Charles | Subsense: A universal change detection method with local adaptive sensitivity[END_REF], BMTDL [START_REF] Staglianò | Online space-variant background modeling with sparse coding[END_REF], GFL [START_REF] Xin | Background subtraction via generalized fused lasso foreground modeling[END_REF], GOSUS [START_REF] Xu | Gosus: Grassmannian online subspace updates with structured-sparsity[END_REF], LR-FSO [START_REF] Xue | Foreground estimation based on linear regression model with fused sparsity on outliers[END_REF], RMAMR [START_REF] Ye | Foreground-Background Separation From Video Clips via Motion-Assisted Matrix Restoration[END_REF], BRTF [START_REF] Zhao | Bayesian Robust Tensor Factorization for Incomplete Multiway Data[END_REF], GoDec [START_REF] Zhou | Godec: Randomized low-rank and sparse matrix decomposition in noisy case[END_REF],

and DECOLOR [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF] on publicly available datasets, such as Change Detection (CDnet) 2014 [START_REF] Wang | Cdnet 2014: an expanded change detection benchmark dataset[END_REF], I2R [START_REF] Li | Statistical modeling of complex backgrounds for foreground object detection[END_REF], Background Models Challenge (BMC) [START_REF] Vacavant | A benchmark dataset for outdoor foreground/background extraction[END_REF], and Wallflower [START_REF] Toyama | Wallflower: Principles and practice of background maintenance[END_REF].

The remaining content of this paper is organized as follows. In Section II, related work is reviewed. In Section III, we describe our method in detail. The experimental results are discussed in Section IV. Finally, our conclusions and future research directions are presented in Section V.

II. RELATED WORK

During the past few years, many research studies have been carried out on background subtraction or foreground detection [START_REF] Cao | Total Variation Regularized RPCA for Irregularly Moving Object Detection Under Dynamic Background[END_REF], [START_REF] Diego | Rejection based multipath reconstruction for background estimation in video sequences with stationary objects[END_REF], [START_REF] Haines | Background subtraction with DirichletProcess mixture models[END_REF], [START_REF] Javed | Motion-Aware Graph Regularized RPCA for Background Modeling of Complex Scenes[END_REF], [START_REF] Lee | Background Subtraction Using Background Sets with Image-and Color-Space Reduction[END_REF], [START_REF] St-Charles | Universal Background Subtraction Using Word Consensus Models[END_REF], [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF] as well as background initialization [START_REF] Bouwmans | Scene background initialization: A taxonomy[END_REF], [START_REF] Bouwmans | Robust PCA via principal component pursuit: a review for a comparative evaluation in video surveillance[END_REF], [START_REF] Erichson | Randomized low-rank Dynamic Mode Decomposition for motion detection[END_REF], [START_REF] Maddalena | Towards benchmarking scene background initialization[END_REF], [START_REF] Ye | Foreground-Background Separation From Video Clips via Motion-Assisted Matrix Restoration[END_REF]. In background subtraction, the emphasis is to improve the accuracy of foreground detection. On the other hand, the task of estimating a foreground-free image is called background modeling. Many surveys have also contributed to these topics [START_REF] Bouwmans | Scene background initialization: A taxonomy[END_REF], [START_REF] Bouwmans | Decomposition into Low-rank plus additive matrices for background/foreground separation: A review for a comparative evaluation with a large-scale dataset[END_REF], [START_REF] Bouwmans | Robust PCA via principal component pursuit: a review for a comparative evaluation in video surveillance[END_REF], [START_REF] Maddalena | Towards benchmarking scene background initialization[END_REF]. Gaussian Mixture Model (GMM) [START_REF] Stauffer | Adaptive background mixture models for real-time tracking[END_REF], [START_REF] Zhang | Mining semantic context information for intelligent video surveillance of traffic scenes[END_REF] uses a mixture of Gaussian probability density functions to model color intensity variations at each pixel. Latest GMM enhancements include bidirectional analysis [START_REF] Shimada | Background modeling based on bidirectional analysis[END_REF] and minimum spanning tree GMM [START_REF] Chen | Spatiotemporal gmm for background subtraction with superpixel hierarchy[END_REF]. Most GMM based methods suffer performance degradation in complex dynamic scenes. Recently, deep Convolutional Neural Network (CNN) based methods have also been proposed for foreground segmentation [START_REF] Braham | Deep background subtraction with scene-specific convolutional neural networks[END_REF], [START_REF] Wang | Interactive deep learning method for segmenting moving objects[END_REF], [START_REF] Zhang | Deep learning driven blockwise moving object detection with binary scene modeling[END_REF]. For instance, Wang et al. [START_REF] Wang | Interactive deep learning method for segmenting moving objects[END_REF] proposed a simple and effective CNN based method for estimating foreground regions. CNN based methods work well in many complex situations however, these methods require significant amount of labelled training data which may not always be available. In contrast, our proposed algorithm is unsupervised therefore do not require labelled training data.

In the current study, we propose a novel algorithm for background modeling, which draws inspiration from subspace learning methods. Wright et al. [START_REF] Wright | Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[END_REF] presented the first proposal of RPCA to handle the outliers in data. Candeś et al. [START_REF] Candès | Robust principal component analysis?[END_REF] used RPCA for background-foreground separation. RPCA-based approaches for background-foreground separation are not ideal for surveillance applications because these approaches suffer from high computational complexity. Moreover, traditional RPCA implementations processed data in batches. Many studies attempting to make the batch methods faster have been reported in literature [START_REF] Cao | Total Variation Regularized RPCA for Irregularly Moving Object Detection Under Dynamic Background[END_REF], [START_REF] Zhou | Godec: Randomized low-rank and sparse matrix decomposition in noisy case[END_REF], [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF]. However, batch methods are not real-time and mostly work offline. Some online methods have also been reported to handle this problem, while global optimality is still the challenging issue in these approaches [START_REF] He | Incremental gradient on the grassmannian for online foreground and background separation in subsampled video[END_REF], [START_REF] Javed | Robust background subtraction to global illumination changes via multiple features-based online robust principal components analysis with markov random field[END_REF], [START_REF] Xu | Gosus: Grassmannian online subspace updates with structured-sparsity[END_REF].

Many authors have contributed interesting studies in the direction of enhancing only foreground detection2 . For this purpose, a number of constraints have been suggested [START_REF] Cao | Total Variation Regularized RPCA for Irregularly Moving Object Detection Under Dynamic Background[END_REF], [START_REF] Javed | Robust background subtraction to global illumination changes via multiple features-based online robust principal components analysis with markov random field[END_REF], [START_REF] Liu | Background subtraction based on low-rank and structured sparse decomposition[END_REF], [START_REF] Xin | Background subtraction via generalized fused lasso foreground modeling[END_REF], [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF]. For example, Cao et al. [START_REF] Cao | Total Variation Regularized RPCA for Irregularly Moving Object Detection Under Dynamic Background[END_REF] improved the performance of foreground detection by proposing the total variation regularized RPCA method. Zhao et al. [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF] proposed a markov random field constraint on the foreground matrix to eliminate noise and small background movements. Although the segmentation performance improved, the foreground regions tend to be over-smoothed [START_REF] Liu | Background subtraction based on low-rank and structured sparse decomposition[END_REF] because of neighboring pixels smoothing constraints.

Unfortunately, research has not been focused in improving low-rank background modeling3 [START_REF] Bouwmans | Decomposition into Low-rank plus additive matrices for background/foreground separation: A review for a comparative evaluation with a large-scale dataset[END_REF]. Therefore, there is a need to design a robust algorithm to recover background components in real-life challenging scenarios [START_REF] Maddalena | Towards benchmarking scene background initialization[END_REF], [START_REF] Pierre-Marc | Scene Background Modeling.net[END_REF]. We attempt to fill this research gap by proposing a novel notion of encoding spatiotemporal similarity information in the component of the background model. Our idea is motivated by the recently proposed subspace clustering methods in [START_REF] Shahid | Robust principal component analysis on graphs[END_REF], [START_REF] Yin | Dual graph regularized latent low-rank representation for subspace clustering[END_REF], [START_REF] Zhang | Low-rank matrix approximation with manifold regularization[END_REF]. However, unlike these methods, we incorporate a motion-aware CC and GD based spatiotemporal regularization for clustering low-rank subspace. Our proposed method has some similarity with RMAMR [START_REF] Ye | Foreground-Background Separation From Video Clips via Motion-Assisted Matrix Restoration[END_REF] because both methods use motion mask. However, RMAMR has not suggested the spatiotemporal consistency constraints as we proposed using Laplacian matrices of graph structures. Our proposed objective function is novel and to the best of our knowledge, no similar work has been proposed previously.

III. PROPOSED ALGORITHM In this section, we describe each step of the proposed MSCL algorithm in detail. The system diagram of the proposed MSCL algorithm is presented in Fig. 4. We first compute the optical flow and create a binary motion mask. This motion mask enables the removal of motionless video frames from the input sequence. Then, spatial and temporal graphs are constructed using the correlation coefficient. The sparse linear decomposition is applied on each graph to estimate the normalized graph Laplacian matrices both in geodesic and correlation space. The four Laplacian matrices computed from these graphs are then incorporated in the RPCA framework. The proposed objective function is solved using the LADMAP optimization scheme. Finally, the Fused Lasso (FL) regularization is applied on the sparse component to enhance the foreground detection.

A. Motion Estimation

We incorporate motion information in our proposed algorithm by computing the dense optical flow [START_REF] Liu | Beyond pixels: exploring new representations and applications for motion analysis[END_REF] between each pair of consecutive frames in the given sequence X. Using motion information, we compute a motion mask M and a sequence of dynamic frames D. Let x i and x i-1 be the two consecutive frames in X at time t and t -1, respectively. Let v x i,k be the horizontal component of the motion vector and v y i,k be the vertical component at position k computed between the frames x i and x i-1 . Let m i ∈ {0, 1} be the corresponding motion mask, which is computed as

m i,k = 1, if (v x i,k ) 2 + (v y i,k ) 2 < τ, 0, otherwise, (2)
where τ is the threshold of motion magnitude, which is computed adaptively as the average of all pixels in the motion field. We select the threshold τ such that all pixels in X exhibiting motion larger than τ definitely belong to the foreground. The threshold τ is selected to be large enough so that the motion should not result because of noise in the background.

To prepare matrix D by eradicating the motionless frames in X using (2), the i-th frame in X is considered to be redundant or motionless if all entries are 1 in the corresponding i-th column of M; otherwise, if some entries are 0, then the frame is considered as dynamic and is appended in matrix D. Using this technique, the dimension of matrix X may be reduced significantly depending on the number of motionless frames.

B. Mathematical Formulation

Given the sequence D ∈ R p×c , where c is the number of dynamic frames, we require that the corresponding matrix B with singular vectors lie in a low-dimensional subspace by minimizing the loss function defined in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. Here, we define a matrix M ∈ R p×c4 as a concatenation of [m 1 , m 2 , .., m c ] as computed by [START_REF] Bouwmans | Scene background initialization: A taxonomy[END_REF]. We consider the pixel positions in D corresponding to 0 values in M as missing data and the remaining pixels as the observed data. The main objective is to estimate these missing values using optimization as:

B,F min||B|| * +λ 1 ||F|| 1 such that M • D = M • (B + F), (3)
where '•' denotes the element-wise product. The observed M • D data still contains a significant number of outliers, which actually does not belong to the background. To handle the missing values and outliers effectively in MSCL, we incorporate spatiotemporal smoothness constraints into (3) by encoding the pairwise similarities among the video frames and pixels. We compute a temporal graph across columns of D and a spatial graph across spatial locations of D using motion-aware Correlation-Coefficient (CC). Based on CCbased adjacency matrix, we compute another representation called Geodesic Distance (GD)-based adjacency matrix.

For each adjacency matrix, we compute a normalized Laplacian matrix, L c t ∈ R c×c corresponding to the CC-based adjacency matrix and L g t ∈ R c×c corresponding to the GD-based adjacency matrix. Similarly, for the spatial graph, we compute the normalized Laplacian matrices L c s ∈ R p×p corresponding to the CC-based adjacency matrix and L g s ∈ R p×p corresponding to the GD-based adjacency matrix. The CC-based regularization enforces the background model to be invariant to lighting condition variations and to enhance continuity on linear manifolds. The GD-based regularization ensures the continuity of the background model on the nonlinear manifolds. We observed that the integration of both types of regularization significantly improves the performance of the background model. The proposed MSCL model is then reformulated as

min B,F ||B|| * +Θ c (B, F) + Θ g (B, F) + λ 1 ||F|| 1 such that M • D = M • (B + F), (4)
where the functions Θ c and Θ g denote the spatiotemporal CC and GD-based regularization. It is important to note that the functions in (4) are explicitly related with D and B and can be shown to be related with B and F as:

Θ c (B, D) = Θ c (B, F + B) = Θ c (B, F). Θ c (B, F) = γ 1 2 Tr(B L c s B) + γ 2 2 Tr(BL c t B ), Θ g (B, F) = γ 3 2 Tr(B L g s B) + γ 4 2 Tr(B L g t B). (5)
We constrain the background model B to minimize tr(B L c s B), tr(BL c t B ), tr(B L g s B), and tr(BL g t B ), which is a spatiotemporal graph regularization. The parameters λ 1 , γ 1 , γ 2 , γ 3 , γ 4 > 0 assign relative importance to each of the terms while optimizing (4). If parameters γ 1 , γ 2 , γ 3 , and γ 4 , are equal to zero, then the model will degenerate to motionaware RPCA.

C. Temporal Graph Construction

For the sequence of dynamic frames D, we construct an undirected weighted graph G t = (V t , A t ), such that the vertices V t correspond to the columns of D and A t is the edge weight matrix. First, we compute the pairwise motion-aware correlation coefficient between frames (d i , d j ) [START_REF] Da | Pearson's Correlation Coefficient: A More Realistic Threshold for Applications on Autonomous Robotics[END_REF] as

ρ t (i, j) = p k=1 m i,j,k (d i,k -di )(d j,k -dj ) p k=1 m i,j,k [(d i,k -di ) 2 + (d j,k -dj ) 2 ] , ( 6
)
where k is the pixel position and m i,j,k = m i,k ∧ m j,k is a motion indicator computed from the motion masks of the two frames, m i,j,k = 1 if both m i,k = 1 and m j,k = 1; otherwise, it will be zero. di is the mean of the frame d i . The value of ρ t (i, j) ranges from -1 to +1, where +1 indicates a perfect linear relationship and -1 is a perfect inverse relationship, which may happen because of abrupt light intensity variations. We compute the weight matrix A t as normalized Euclidean distance matrix:

A t (i, j) = 2(1 -ρ t (i, j)). (7)
The proximity matrix A c t of the graph G t is computed as

A c t = exp -A t /2σ 2 t , (8)
where σ 2 t is the smoothing factor that strengthens the fusion between the adjacent pixels. We use σ t as an average distance among the nodes in G t [START_REF] Shahid | Robust principal component analysis on graphs[END_REF]. If A c t (i, j) > 0, then there is an edge between nodes d i and d j in G t .

D. Spatial Graph Construction

Using motion-aware correlation coefficient, we compute a spatial graph G s = (V s , A s ) such that V s corresponds to the rows of matrix D and A s is the edge weight matrix computed similar to A t . The spatial graph G s complements the information captured by the temporal graph G t . Given that G s encodes the notion of similarity among the spatial locations; therefore, it will enforce smoothness on the patch level in spatial dimensions in the low-rank model. Moreover, computing the correlation coefficient among the rows of patches allows the use of local information in the image sequence.

For each pixel in the first frame in D, we consider a patch of size u×u pixels in the corresponding 2-D image. Therefore, for each row of pixels in D, we obtain a row of patches or patch-row having the size c × u × u. The motion-aware correlation coefficient ρ s (i, j) is computed between all pairs of patch-rows resulting in a p×p matrix of correlations, which is transformed to normalized Euclidean distance A s by using

A s (i, j) = 2(1 -ρ s (i, j)). Spatial proximity matrix A c s is obtained by A c s (i, j) = exp -A s (i, j)/2σ 2 s
, where σ 2 s acts as a spatial smoothing factor computed as an average distance among nodes in G s .

E. Geodesic Distance Based Proximity Matrices

For the spatial and temporal graphs G t = (V t , A t ) and G s = (V s , A s ), we compute all pairs of shortest distance or Geodesic Distance (GD) [START_REF] Jiang | A new approach to shortest paths on networks based on the quantum bosonic mechanism[END_REF], [START_REF] Johnson | Efficient algorithms for shortest paths in sparse networks[END_REF]. We obtain two more proximity matrices A g t ∈ R c×c and A g s ∈ R p×p such that the i-th column in these matrices corresponds to the geodesic distances of the i-th node from the rest of the nodes in that graph. These matrices will be used to construct regularization, which makes the low-rank robust to the nonlinear variations in spatial and temporal dimensions.

For the temporal graph G t , the geodesic distance matrix P g t is computed using A t . For the spatial graph G s , the geodesic distance matrix P g s is computed using A s . The geodesic distance between two nodes is more meaningful than the Euclidean distance computed in the image space because two frames may have smaller Euclidean distance, while geodesic distance is larger along the nonlinear manifold. Using the geodesic distance matrix P g ∈ {P g t , P g s }, we compute the geodesic proximity matrices A g ∈ {A g t , A g s },

A g = exp(-P g • P g /2σ 2 g ), (9)
where σ g ∈ {σ s g , σ t g } provides smoothing in the geodesic space and computed as the average distance between the nodes in the corresponding graphs. In general, if a node d i is not reachable from node d j , then A g (i, j) = ∞ and P g (i, j) = 0.

F. Sparse Linear Coefficients Proximity Matrices

We observe that the columns of a proximity matrix A ∈ {A c t , A c s , A g t , A g s } span low dimensional subspaces compared to the overall matrix dimensions because the size of the temporal proximity matrices A c t and A g t is the number of frames × number of frames. In case there is no change in the background, ideally, the rank of both of these matrices will be 1. However, given the dynamic nature of the background, we observe that the rank of these matrices will be larger than 1, but much smaller than the number of frames. In Fig. 3, we plot the temporal and spatial graphs of the dynamic background sequences. Each graph contains multiple nonlinear clusters showing the presence of multiple manifolds. Therefore, a better background model can be computed by using subspace-based spectral clustering [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF]. We compute the proximity matrices by applying sparse linear decomposition on the matrix A. The computation of normalized Laplacian using sparse linear coefficients has been proposed in subspace sparse spectral clustering [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF], which may be reviewed for more theoretical details.

Let a i be the i-th column of A, which is decomposed as a linear combination of the rest of the matrix  = A \ a i .  is the same as A except that the a i column is replaced by a column of zeros. We find an optimal solution for sparse linear coefficients α α α i :

α α α * i := arg min α α αi (||a i -Âα α α i || 2 2 +λ 2 ||α α α i || 1 ), (10)
where λ 2 gives a relative importance to the 1 norm factor, which approximates the o norm actually required for sparsity. λ 2 is automatically computed from the data [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF]. The sparse linear decomposition with l 1 -norm constraints given by ( 10) is an unconstrained convex optimization problem also known as the Least Absolute Selection and Shrinkage Operator (LASSO) [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. We solve this problem using a fast solution available in the Least Angle Regression (LARS) framework [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF], which has the same asymptotic complexity as the simple least squares regression. Sparse linear decomposition is applied on each column of A and the resulting α α α i are arranged as columns resulting in a matrix of sparse linear coefficients E where E ∈ {E c

t , E c s , E g t , E g s }.

G. Computing Normalized Laplacian Matrices

The matrices of linear coefficients may not be symmetric E = E , which is partially because the set of vectors used to represent each vector a i ∈ A is slightly different from the set used to represent a j ∈ A. We make this relationship normalized and symmetric by taking the average of both coefficients normalized by the maximum value of its own set [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF]:

E(i, j) = E(j, i) = 1 2 α i,j max(α α α i ) + α j,i max(α α α j ) . (11)
Using each normalized symmetric matrix E of sparse linear coefficients, we compute a normalized graph Laplacian matrix

L where L ∈ {L c t , L c s , L g t , L g s } L = I -W -1/2 EW -1/2 , ( 12
)
where I is an identity matrix of the same size as E, W is the degree matrix whose diagonal entries are given by (j, j) = i E(i, j), and non-diagonal values are zero. The normalized Laplacian matrices are used to constrain the background model B, minimizing either the trace of B L B or BLB depending on the dimensions of L, as shown in the objective function (4). In the following section, we present the solution of the proposed objective function.

H. Proposed LADMAP Optimization

Model ( 4) is essentially a convex optimization problem, which is solved using the Linearized Alternating Direction Method with Adaptive Penalty (LADMAP) [START_REF] Lin | The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices[END_REF], [START_REF] Lin | Linearized alternating direction method with adaptive penalty for low-rank representation[END_REF], [START_REF] Wen | Alternating direction augmented Lagrangian methods for semidefinite programming[END_REF]. For this purpose, the linear equality constraints in (4) are removed by employing following augmented Lagrangian formulation:

L(B, F, Y, µ) = min B,F ||B|| * +Θ c (B, F) + Θ g (B, F) + λ 1 ||F|| 1 +tr(Y(M • (D -B -F))) + µ 2 ||M • (D -B -F)|| 2 F , (13)
where Y ∈ R p×c is a Lagrangian multiplier matrix and µ > 0 controls the penalty for violating the linear constraints.

Optimizing directly the primary variables B, F, and Y, the LADMAP method solves each variable iteratively one after another. LADMAP uses less auxiliary variables without matrix inversions and converges faster than the original ADM [START_REF] Wen | Alternating direction augmented Lagrangian methods for semidefinite programming[END_REF].

We derive the solutions to use the proximity operator of the nuclear norm effectively and l 1 -norm in solving subproblems including Laplacian terms.

Updating B: According to LADMAP, the linearization only over the augmented quadratic penalty term in [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF] solves subproblems B k+1 and F k+1 (in this case, k is the iteration index). By fixing F, the update for background B k+1 at the (k + 1)-th iteration is

B k+1 = arg min B L(B, F k , Y k , µ k ) = arg min B ||B|| * +Θ c (B, F) + Θ g (B, F) + Tr(Y k (M • (D -B -F k ))) + µ k 2 ||M • (D -B -F k )|| 2 F (14) If γ 1 , γ 2 , γ 3 , γ 4 > 0 in functions Θ c
and Θ g , then optimizing subproblem B k+1 does not lead to the exact solution. Therefore, the convergence analysis provided in [START_REF] Lin | Linearized alternating direction method with adaptive penalty for low-rank representation[END_REF] is not applicable. To use the closed-form solution to the proximity operator of the nuclear norm, which is given by Singular Value Thresholding (SVT) operator [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF], we further linearize the graph regularization terms to simplify the subproblem B k+1 . According to Lemma 2.1 in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], the trace norm of any graph regularization terms, e.g., γ1 2 tr(B L c s B) in ( 12) can be upper bounded by its proximal approximation which describes the local linearization of γ1 2 tr(B L c s B) at B k . Then, the solution of ( 14) leads to the following problem after applying simple linear algebra as

B k+1 = argmin B ||B|| * + η 1 2 × M • (B -B k ) + H η 1 2 F ( 15
)
For our convenience, we denote

                   η 1 = µ k + ρ 1 , ρ 1 = γ 1 ||L c s ||+γ 2 ||L c t ||+γ 3 ||L g s ||+γ 4 ||L g t || H = -(M • D) Ỹk + ρ ρ ρ 2 , ρ ρ ρ 2 = γ 2 B k L c t + γ 4 B k L g t + [γ 1 B k L c s + γ 3 B k L g s ] , Bk = (M • k ) - 1 η 1 [-(M • D) Ỹk + ρ ρ ρ 2 ], (16)
To solve B k+1 , first we need to compute the parameter µ k and matrix Ỹk as

             µ k+1 = µ k + ρ 0 µ max µ max = max{γ 2 ||L c t ||+γ 4 ||L g t ||, γ 1 ||L c s ||+γ 3 ||L g s ||} ||M • D|| 2 -λ 1 Ỹk = Y k + µ k M • (D -B k -F k ), Y k+1 = Y k + µ k M • (D -B k+1 -F k+1 ) (17)
where ρ 0 is constant and the regularization e.g., L c s is known as the spectral norm of a matrix L c s , which is the largest singular value of matrix L c s . The B k+1 has the following closed-form solution as

B k+1 = UT τ (Σ)V , where(U, Σ, V ) = SVD( Bk ) (18)
where Σ is the singular value matrix of Bk . The operator T τ (•) is the SVT [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF], which is defined by element-wise τ thresholding of Σ. Especially, let σ i be the i-th diagonal element of Σ, then T τ (Σ) is a diagonal matrix defined by T τ (Σ) = diag({(σ i ) -τ } + ), where a + is the positive part of a + = max(0, a). Updating F: When B is fixed, to update F k+1 , we derive from (13) the following sub-problem:

F k+1 = arg min F L(B k+1 , F, Y k , µ k ) = arg min F λ||F|| 1 +Tr(Y(M • (D -B k+1 -F))) + µ k 2 ||M • (D -B k+1 -F)|| 2 F ( 19
)
The closed-form solution of the above problem is

F k+1 = S λ/µ k M • (D -B k+1 ) + Ỹk µ k (20)
where S τ (•) is the shrinkage operation [START_REF] Lin | The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices[END_REF] defined by

S τ (F) = sgn(F) • max{|F|-τ 1 , 0}. (21)
Algorithm 1: Proposed MSCL for B modeling Input: X ∈ R p×n Initialization: D, M ∈ R p×c using (2), L c t , L g t ∈ R c×c , and L c s , L g s ∈ R p×p using ( 12), B 0 = F 0 = 0, γ = 10, 1 = 10 -4 , 2 = 10 -5 , µ 0 = 0.1, ρ 0 = 1.1, µ max using [START_REF] Haines | Background subtraction with DirichletProcess mixture models[END_REF] while not converged (k = 0, 1, ..) do 1. Compute Ỹk using (17) 2. Update B k+1 , F k+1 in parallel using ( 18), (20) 3. Update Y k+1 and µ k+1 according to [START_REF] Haines | Background subtraction with DirichletProcess mixture models[END_REF]. [START_REF] Bouwmans | Robust PCA via principal component pursuit: a review for a comparative evaluation in video surveillance[END_REF]. Check convergence: according to [START_REF] Jiang | A new approach to shortest paths on networks based on the quantum bosonic mechanism[END_REF]. end Output: B * , F * Convergence Condition: According to the KKT condition, the following criterion is defined for the sub-optimality of the solution of problem ( 4)

M • (D -B k+1 -F k+1 ) / M • D < 1 and max{ρ 3 , ρ 4 } ≤ 2 , where ρ 3 = µ k M • (F k+1 -F k ) , ρ 4 = µ k + 2ρ 1 √ λ 1 M • (B k+1 -B k ) , (22)
where 1 and 2 are the appropriate tolerances. Based on the stopping criteria defined in [START_REF] Jiang | A new approach to shortest paths on networks based on the quantum bosonic mechanism[END_REF], the sequences (B, F, Y) generated by the revised LADMAP converges to an optimal solution of problem (4). Algo. 1 describes the summary of

I. Foreground Detection

The additional constraints on the low-rank matrix enhanced the quality of the background model and improved the foreground detection. Often, the foreground is relatively smaller and sparser compared to the background. Consequently, minimizing the foreground term by imposing the l 1 norm on F makes the estimation even more sparser. We observe that without considering the spatial connectivity among the adjacent foreground pixels, especially in dynamic backgrounds, it is difficult to obtain precise foreground segmentation (Fig. 3).

To improve the foreground segmentation, we constrain the foreground pixels by employing Fused Lasso (FL) [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] as a post-processing step. We name the resulting algorithm as MSCL-FL for foreground detection, which is robust against dynamic background scenes. In MSCL-FL, it is assumed that if a pixel belongs to the foreground, then its neighbors would also belong to the foreground. The overall goal of MSCL-FL is to minimize the following energy function as: [START_REF] Johnson | Efficient algorithms for shortest paths in sparse networks[END_REF] where f (k) is the k-th vector of F and N is the spatial neighborhood system, i.e., (i, j) ∈ N when both pixels i and j are spatially connected. The first term ||f|| 1 is the observed data term that connects each pixel to background and foreground nodes. The second term is known as the smoothness term, which represents the relationship between Fig. 5. AGE measure (minimum is best) using different values of parameters γ 1 ,γ 2 ,γ 3 , and γ 4 . An experiment performed on "CameraParameter" sequence of the SBM.net dataset [START_REF] Pierre-Marc | Scene Background Modeling.net[END_REF].

||F|| FL = c k=1 {||f (k) || 1 +β (i,j)∈N w (k) ij |f (k) i -f (k) j |},
adjacent pixels i and j. β is a scalar parameter and w ij is the adaptive weighting factor that stabilizes the fusion between pixels as w

(k) ij = exp -||d (k) i -d (k) j || 2 2 /2σ 2
, where d is the pixel intensity of k-th frame and σ is the spread of the Gaussian function. The smoothness term encourages the same labeling to similar neighboring pixels. The energy function [START_REF] Johnson | Efficient algorithms for shortest paths in sparse networks[END_REF], which is nonconvex, is difficult to solve especially if the dimensionality is large. Therefore, we use the fast parametric flow method [START_REF] Gallo | A fast parametric maximum flow algorithm and applications[END_REF], [START_REF] Xin | Background subtraction via generalized fused lasso foreground modeling[END_REF].

IV. EXPERIMENTS

We performed extensive experiments on five publicly available video datasets including Scene background modeling.net (SBM.net) [START_REF] Pierre-Marc | Scene Background Modeling.net[END_REF], Change Detection.net (CDnet) 2014 [START_REF] Wang | Cdnet 2014: an expanded change detection benchmark dataset[END_REF], Background Models Challenge (BMC) [START_REF] Vacavant | A benchmark dataset for outdoor foreground/background extraction[END_REF], I2R [START_REF] Li | Statistical modeling of complex backgrounds for foreground object detection[END_REF], and Wallflower [START_REF] Toyama | Wallflower: Principles and practice of background maintenance[END_REF]. The results are compared with 23 state-ofthe-art methods, such as TVRPCA [START_REF] Cao | Total Variation Regularized RPCA for Irregularly Moving Object Detection Under Dynamic Background[END_REF], RMR [START_REF] Diego | Rejection based multipath reconstruction for background estimation in video sequences with stationary objects[END_REF], 2P-RPCA [START_REF] Gao | Block-sparse RPCA for salient motion detection[END_REF], RFSA [START_REF] Guo | Robust foreground detection using smoothness and arbitrariness constraints[END_REF], DP-GMM [START_REF] Haines | Background subtraction with DirichletProcess mixture models[END_REF], GRASTA [START_REF] He | Incremental gradient on the grassmannian for online foreground and background separation in subsampled video[END_REF], TLSFSD [START_REF] Hu | Moving object detection using tensor-based low-rank and saliently fused-sparse decomposition[END_REF], our previous algorithm MAGRPCA [START_REF] Javed | Motion-Aware Graph Regularized RPCA for Background Modeling of Complex Scenes[END_REF], LSD [START_REF] Liu | Background subtraction based on low-rank and structured sparse decomposition[END_REF], 3TD [START_REF] Oreifej | Simultaneous video stabilization and moving object detection in turbulence[END_REF], MODSM [START_REF] Pang | Incremental learning with saliency map for moving object detection[END_REF], MLRBS [START_REF] Qin | A background basis selectionbased foreground detection method[END_REF], SRPCA [START_REF] Sajid | Spatiotemporal low-rank modeling for complex scene background initialization[END_REF], PAWCS [START_REF] St-Charles | Universal Background Subtraction Using Word Consensus Models[END_REF], SuBSENSE [START_REF] St-Charles | Subsense: A universal change detection method with local adaptive sensitivity[END_REF], BMTDL [START_REF] Staglianò | Online space-variant background modeling with sparse coding[END_REF], GFL [START_REF] Xin | Background subtraction via generalized fused lasso foreground modeling[END_REF], GOSUS [START_REF] Xu | Gosus: Grassmannian online subspace updates with structured-sparsity[END_REF], LR-FSO [START_REF] Xue | Foreground estimation based on linear regression model with fused sparsity on outliers[END_REF], RMAMR [START_REF] Ye | Foreground-Background Separation From Video Clips via Motion-Assisted Matrix Restoration[END_REF], BRTF [START_REF] Zhao | Bayesian Robust Tensor Factorization for Incomplete Multiway Data[END_REF], GoDec [START_REF] Zhou | Godec: Randomized low-rank and sparse matrix decomposition in noisy case[END_REF], and DECOLOR [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF], using original author implementations. Background models are compared using Average Gray-level Error (AGE), percentage of Error Pixels (pEPs), Percentage of Clustered Error Pixels (pCEPs), Multi Scale Structural Similarity Index (MSSSIM), Color image Quality Measure (CQM), and PSNR [START_REF] Maddalena | Towards benchmarking scene background initialization[END_REF], [START_REF] Pierre-Marc | Scene Background Modeling.net[END_REF]. We aim to minimize AGE, pEPs, and pCEPs for more accurate background model recovery while maximizing MSSSIM, PSNR, and CQM.

The solution of the proposed model (4) requires parameters λ 1 , γ 1 , γ 2 , γ 3 , and γ 4 . We use λ 1 = 1/ max(p, c), where p is the number of pixels and c is the number of dynamic frames in matrix D [START_REF] Candès | Robust principal component analysis?[END_REF]. In our experiments, we used γ 1 = γ 2 = γ 3 = γ 4 = 10. We empirically verified that further tuning of these parameters may yield some improvement in accuracy (Fig. 5). However, to make the experiments repeatable, all results are reported with the same value given previously. The parameter µ k in ( 12) is initialized with 0.1 and then adjusted adaptively in the later iterations. For the construction of G s on image patches, we used the patch size of 5×5 pixels. The parameters required by the LADMAP are used as recommended by the original authors [START_REF] Lin | Linearized alternating direction method with adaptive penalty for low-rank representation[END_REF], [START_REF] Liu | Linearized alternating direction method with parallel splitting and adaptive penalty for separable convex programs in machine learning[END_REF], [START_REF] Yin | Dual graph regularized latent low-rank representation for subspace clustering[END_REF]. In addition, we used the PROPACK library [START_REF] Larsen | Lanczos bidiagonalization with partial reorthogonalization[END_REF] for LADMAP implementation.

A. Evaluation of MSCL on SBMI Dataset

We evaluated our proposed MSCL algorithm on the SBM.net dataset, which contains 8 challenging categories and 80 videos (Fig. 6). In contrast to the existing datasets, the SBM.net comprises complex scenes in which no foregroundfree images are available and the background is largely occluded by foreground objects. More details can be found at http://scenebackgroundmodeling.net. On the average, our algorithm outperformed all existing methods in terms of 6 accuracy measures included in the study [START_REF] Pierre-Marc | Scene Background Modeling.net[END_REF] (Figs. 7 and8). The MSCL outperformed all existing methods in three categories, namely, 'Basic', 'Clutter', and 'Illumination changes' (Fig. 8). A visual comparison of the estimated background model with best performing methods over 8 selected sequences (one sequence per category) is illustrated in Fig. 6. In Fig. 8, we present only 5 noteworthy methods and analyze their average performance on each category using the AGE metric, which is the l 1 -norm of the difference of ground truth and the estimated background image. Lower AGE value denotes better background estimation. For the 'Basic' category, only MSCL and RMR produced good results, while the remaining methods exhibit some discrepancies in the form of outliers because background pixels are mostly visible all the time in these sequences, although there are no frames containing complete background. In the challenging sequences of 'Clutter' category, the background remains largely occluded by the foreground objects. Fig. 8 shows that the proposed algorithm attains, on the average, the best accuracy in terms of AGE value. In the compared methods, overwhelming outliers of foreground objects were incorporated into the estimated background model (Fig. 6, 3 rd row). For 'Jitter' sequences, most of the methods perform good estimation of the background component. In this case, the proposed MSCL algorithm has shown the best performance. RMAMR and MAGRPCA also use motion information; there- fore, these methods also show better performance. For complex dynamics as presented by the sequences in the 'Background Motion' category, all methods have shown relatively poor performance because of the sudden change in the background in the presence of foreground objects, which always remain visible. In this category, MSCL is in the middle of the spectrum obtaining an AGE score of 11.21, while the best performing competitor DECOLOR obtained 10.5. Note that the previous version of our proposed algorithm MAGCPCA has obtained the best accuracy at 10.07, which shows that the inclusion of the correlation coefficient-based regularization has caused some degradation in this case, while the overall accuracy has significantly increased. In the 'Intermittent Motion' category, some foreground objects remain motionless and then start moving, creating outlier in the background. Our algorithm considered the motionless frames as redundant and removed them, which helped in reducing the outliers. This step is only useful for this category, while in other cases, this step has no effect. Fig. 8 demonstrates that majority of the methods failed to cope with incorporated outliers because of motionless foreground objects. Only MSCL and RMR methods effectively handled these outliers in the final estimation of background. The proposed MSCL provides the best performance in terms of AGE measure. For sequences that belong to the 'Illumination Changes' category, Fig. 8 shows that the proposed algorithm has exhibited significant improvement over current methods, which have shown degraded performance. The performance of the proposed algorithm is improved because of the use of the correlation coefficient for the computation of adjacency matrices, which are robust to illumination variations. All of the compared methods have shown relatively good performance on 'Very Long' and 'Very Short' categories. Considering all categories of the SBM.net dataset, on the average, the MSCL obtained an AGE score of 5.95 and has been ranked as the best performer (http://scenebackgroundmodeling.net). The improved performance is attributed to the spatiotemporal regularization, which was able to handle the outliers effectively and ensured smoothness of the background model on temporal and spatial manifolds.

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m)

B. Importance of Different Regularizations

To find the relative contribution of different regularizations in (4), we conduct a series of experiments. Setting all 0 results in LADMAP based RPCA, γ 2 = γ 4 = 0 results in spatially regularized RPCA, γ 1 = γ 3 = 0 results in temporally regularized RPCA, γ 1 = γ 2 = 0 results in nonlinear spatiotemporal regularized RPCA, and γ 3 = γ 4 = 0 results in spatiotemporal regularized RPCA. In addition, we also implemented the proposed objective function without using the Sparse Subspace Clustering (SSC) [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF]. Fig. 9 shows the comparison of these variants with the proposed algorithm with the full objective function on six different sequences from the SBM.net dataset. In all sequences, the performance of RPCA has remained the least (shown in black) and the performance of the proposed MSCL has remained the best (shown in red). The performance of MSCL without SSC (shown in green) has remained poor compared to the proposed MSCL with only spatial, only temporal, and both spatiotemporal constraints.

In the case of the cluttered scene (Fig. 9(a)), the temporal regularization (shown in blue) has performed better compared to spatial regularization (shown in magenta) because of the presence of more outliers along the rows of matrix D. In the simple case (Fig. 9(b)), the spatial constraints have mostly performed better. In varying lighting conditions (Fig. 9(c)), the spatial constraints performed better because they were not affected by the abrupt intensity variations. In the case of jitter (Fig. 9(d)), the spatial constraints performed better because of the abrupt motion of the camera, which degraded the interframe similarity, but retained the high intra-frame similarity.

In the case of linear light intensity variations (Fig. 9(e)),the correlation performed better, while in the case of nonlinear variations with abrupt changes (Fig. 9(f)), the geodesic distance performed better because in this case, the correlation coefficient was not able to handle the nonlinear changes. Therefore, in some cases, spatial regularization performed better, while in other cases, the performance of temporal regularization was better. However, the proposed MSCL algorithm integrated the strengths of different regularization schemes and achieved better performance compared to the individual regularization.

C. Evaluation of MSCL-FL for Foreground Detection

The main objective of the MSCL algorithm is to extract the background model for complex scenes. However, we have also observed improvements in the foreground detection. In this section, we compare the performance of MSCL and MSCL-FL algorithms with 20 existing methods on three publicly available datasets using the F 1 score as the performance measure [START_REF] Haines | Background subtraction with DirichletProcess mixture models[END_REF], [START_REF] Gao | Block-sparse RPCA for salient motion detection[END_REF], [START_REF] Zhou | Moving object detection by detecting contiguous outliers in the low-rank representation[END_REF]. 1) CDnet 2014 dataset [START_REF] Wang | Cdnet 2014: an expanded change detection benchmark dataset[END_REF] consists of eleven categories of videos and provides ground truth for each sequence. We only test seven challenging categories, namely, 'Baseline', 'Dynamic Background', 'Camera Jitter', 'Intermittent Object Motion', 'Thermal', 'Shadows', and 'Bad Weather'. These videos were captured using low-resolution IP cameras as well as thermal cameras. The spatial resolution varies from 320 × 240 to 720 × 576. The level of noise and compression artifacts also vary across videos. Table I shows the quantitative results of all the compared methods. Fig. 10 presents the visual results of the proposed MSCL-FL algorithm only.

The category Baseline contains four simple videos (Fig. 10a). As shown in Table I, all of the compared methods (excluding GRASTA) produce an average F 1 score of around 90%. Therefore, the Baseline category does not pose a challenge for most of the compared methods.

The Dynamic Background category contains six challenging videos (Fig. 10b) depicting outdoor scenes. This is the most difficult among all categories for mounted camera object detection, which contains sequences exhibiting dynamic background motions. Most state-of-the-art methods including DE-COLOR, 3TD, 2P-RPCA, LSD, and BMDTL, generate noisy foreground segments because of highly dynamic background regions. Therefore, the performance of these methods degraded compared to the proposed MSCL-FL algorithm (Table I). Given that no additional constraints are considered on the low-rank matrix; therefore, the FL-based methods, such as LR-FSO and GFL, also show a degraded performance (Table I), while very few methods, such as DP-GMM and RMAMR, attain a comparable performance of 80%. The PAWCS and the proposed algorithms MSCL and MSCL-FL are the best performers for these dynamic sequences. The proposed spatiotemporal continuity in the MSCL algorithm improves the foreground detection among all the compared algorithms. However, the fusion of neighboring information in the proposed MSCL-FL algorithm improves by 4% of the F 1 score. The proposed MSCL-FL algorithm achieves 90% average F 1 which is larger than the current state-ofthe-art methods.

The category Camera Jitter contains one indoor and three outdoor videos (Fig. 10a). The compared methods generate noisy foreground mask resulting in a low F 1 score because of jitter. One of the strengths of the MSCL-FL is to handle effectively the fast motion of the background caused by camera jitter and fast moving foreground objects. The comparison between MSCL-FL and GOSUS (structured sparsity constraints with RPCA) demonstrates that the performance gained by the proposed method is not a general feature of RPCA, but it is a combined result of the spatiotemporal regularization, the optimization method, and the foreground pixel labeling process using FL. The Shadows category comprises six videos (Fig. 10c) exhibiting both strong and faint shadows. For 3TD, LR-FSO, GoDec, and GRASTA methods, this category poses a big challenge (TableI). SuBSENSE, PAWCS, and MSCL-FL have achieved promising performance as compared to other methods. We observe that some hard shadows on the ground (Fig. 10c, Cubicle sequence) are still a major limitation of the top performing algorithms. The Thermal category consists of five sequences captured by the far-IR camera (Fig. 10c).

Color saturation is the main challenge in this category, which degrades the performance of DECOLOR, 2P-RPCA, TVR-PCA, RMAMR, GoDec, LSD, RMAMR, GFL, and GRASTA methods (TableI). The spatiotemporal regularization combined with the FL constraint in MSCL-FL are able to discriminate the background-foreground pixels effectively in the presence of color saturation.

The Intermittent Object Motion category includes six videos (Fig. 10b), which contain ghosting artifacts in the detected motion. All compared methods except SRPCA were not able to handle this challenge and obtained a low F 1 score. The proposed MSCL-FL was able to handle the challenge of ghosting 10a). The moving rain streaks or snowflakes must remain part of the background model. All compared methods except SuBSENSE and PAWCS were not able to handle this challenge efficiently. In the proposed MSCL-FL algorithm, the information of neighboring pixels was useful in handling the isolated moving parts of the background, such as snowflakes and rain streaks. On the average, on CDnet dataset, the MSCL-FL obtained 88% F 1 score, which is significantly larger than the existing best-performing methods including GOSUS, SuBSENSE, and PAWCS. The performance of other methods has remained low (Table I).

2) Background Models Challenge (BMC) 2012 Dataset:

In BMC dataset [START_REF] Vacavant | A benchmark dataset for outdoor foreground/background extraction[END_REF], we compare our results on 9 real videos including challenges, such as the presence of dynamic backgrounds, cast shadows, lighting conditions, intermittent foreground object motion, and bootstrapping. For fair comparison, the automatic evaluation tool by the original authors is used to compute the F 1 score. Fig. 11 presents the visual results of the proposed MSC-FL algorithm. Table II shows the average F 1 score of 11 compared methods. The MSCL-FL algorithm has achieved the best performance over videos 001, 004, 005, 006, 007, and 009, while the MSCL algorithm obtained the best results for 003 and 008. These results are attributed to the small size of the foreground objects, which hindered FL from distinguishing these pixels efficiently. The average F 1 score obtained by the proposed MSCL-FL is 86%, LSD is 83%, and MSCL is 82%, while all other compared methods are significantly lower.

3) I2R Dataset: We have also reported the performance of the proposed algorithms on the I2R dataset [START_REF] Li | Statistical modeling of complex backgrounds for foreground object detection[END_REF]. This dataset comprises 9 complex background videos including crowded foreground (Bootstrap, Shopping Mall, Escalator, and Airport Hall), dynamic background (Campus, Curtain, Fountain, and Water Surface), and drastic illumination variations (lobby).

The proposed algorithms are compared with 14 state-of-theart methods in Table III. On the average, GFL obtained the best performance of 0.85 while our proposed MSCL-FL obtained 0.84. All other compared methods exhibited degraded performance for these videos. In some cases, the MSCL performs better than MSCL-FL because of the over-smoothing induced by the FL. 4) Wallflower Dataset [START_REF] Toyama | Wallflower: Principles and practice of background maintenance[END_REF] consists of 6 challenging videos including Waving Trees, Moved Object, Light Switch, Time of Day, Camouflage, and Foreground Aperture. We compare our proposed algorithms with 12 state-of-the-art methods. On the average, the MSCL-FL outperformed the existing algorithms by a significant margin (TableIII). The GFL obtained significantly low performance of 0.84 compared to 0.92 obtained by MSCL-FL. It is because in the presence of camouflage sequence GFL failed to accurately detect moving foreground objects (Fig. 3 in [START_REF] Xin | Background subtraction via generalized fused lasso foreground modeling[END_REF]). Also in the bootstrap sequence because of cluttered foreground objects GFL suffered from oversmoothing. The closely located distinct objects were merged into single object. While the proposed MSCL-FL algorithm was able to obtain good performance in all cases.

D. Execution Time Comparison

Execution times are compared on a machine with Intel core i7 processor and 8GB RAM. The computational cost of our proposed algorithm is mainly determined by the LADMAP [START_REF] Lin | Linearized alternating direction method with adaptive penalty for low-rank representation[END_REF]. Let k denote the number of iterations and r be the lowest rank for matrix B. The construction of spatiotemporal CC-based graph Laplacians needs O(p 2 c + pc 2 ). For spatiotemporal GD-based graph Laplacians, we use Johnson's algorithm [START_REF] Johnson | Efficient algorithms for shortest paths in sparse networks[END_REF], which requires O(c log(c)+p log(p)+(c+p)e) time complexity, where e denotes edges. In each iteration for solving Alg. (1), SVT is applied to update the low-rank matrix B whose total complexity is O(rc 2 ) when we use partial SVD. Similarly, the soft thresholding operation to update the matrix F has a complexity of O(pc). Therefore, the cost of LADMAP iterations is O(krc 2 + kpc). The dominant computational complexity is O(p 2 c) if p > c and O(pc 2 ) if c > p. We have also compared the computational time of the proposed MSCL algorithm. For this purpose, we selected a highway sequence with 50 frames of size 320 × 240. The overall time taken by MSCL is 39.8 sec, which includes 20 sec taken by the optical flow computation and 19.8 sec by the rest of the algorithm. On the same sequence, GRASTA, DECOLOR, RFSA, RMAMR, 3TD, BSLSD, and GFL took 86, 291, 72, 30, 106, 71, and 66 sec, respectively. The execution time of the proposed algorithm can be further enhanced by using GPU or CPU based parallel implementation. Moreover, a faster optical flow computation method will also make the overall process faster.

V. CONCLUSION

In this study, an algorithm based on RPCA with spatiotemporal sparse spectral clustering based regularization is presented for efficient background modeling. To reduce outliers, motionless frames are removed. The regularization enforced the background model to be continuous on the low dimensional multiple manifolds both in the spatial and temporal dimensions, which is achieved by constructing two graphs to encode the temporal and spatial similarities. The objective function is efficiently solved using LADMAP. The major advantage of the proposed algorithm is its capacity to generate an accurate background model even in the presence of occlusions, clutter, jitter, and abrupt intensity variations. Large scale experimental evaluations on five datasets demonstrated that the proposed algorithm achieved the best performance compared to existing methods. However, videos captured from moving and PTZ cameras and online processing remain open challenges. We plan to investigate the possibility of extending the proposed algorithm to scenes that are more crowded and to those recorded using a moving camera by further extending the notion of data similarity using coarse-to-fine strategy.

Fig. 1 .

 1 Fig. 1. Estimated background examples in the SBM.net dataset [37]: (a) Sequence 'IMB01 from the category 'intermittent object motion contains redundant frames. (b) Sequence 'camera parameter from the category 'illumination changes. (c) Sequence 'People & Foliage from the category 'clutter has large occlusions. (d) Sequence 'boulvardJam from the category 'clutter has dominant foreground. (e) Sequence 'overpass from the category 'dynamic motion has always moving objects. In all of these cases, the proposed MSCL algorithm has produced significant improvement.

Fig. 2 .

 2 (a) and (c) Temporal graphs of the dynamic background sequences 'Water Surface' and 'Campus' from the I2R dataset [26]. (b) and (d) Corresponding spatial graphs of the same sequences. The edges are removed to avoid clutter. These graphs are plotted using three major principal components of the geodesic space. The different colors indicate different clusters computed using subspace spectral clustering in the geodesic space. The shapes of different clusters as well as the individual graphs are nonlinear, indicating that the background sequences span multiple nonlinear manifolds.

Fig. 3 .

 3 Fig. 3. Results of foreground detection. (a)-(b) 'Fall' sequence from the CDnet dataset, (c) 'Overpass' from the CDnet dataset, (d) 'Waving Trees' from the Wallflower dataset, and (e) 'Water Surface' the from I2R dataset. From top to bottom: (1) input images, (2) ground-truth images of the foreground, (3) foreground estimation using RPCA methods[START_REF] Oreifej | Simultaneous video stabilization and moving object detection in turbulence[END_REF], (4) estimated foreground mask by using the proposed MSCL algorithm, and (5) results of the foreground estimation by using our proposed MSCL-FL algorithm.

Fig. 4 .

 4 Fig. 4. Schematic illustration of the proposed MSCL algorithm for background estimation and foreground segmentation. Step (1) describes the motion estimation module of sequence X, the creation of motion mask M, and detection of dynamic sequence D. Step (2) describes the construction of spatial and temporal graphs and computation of correlation based and geodesic distance based proximity matrices. In step (3), sparse coding is applied on the proximity matrices computed in step (2) and computation of normalized graph Laplacian matrices. Step (4) shows the decomposition of sequence D into B and F components based on LADMAP optimization. In step (5), the FL constraints are employed on the foreground to ensure continuity.

Fig. 7 .

 7 Fig. 7. Average performance on each metric of MSCL with existing methods on the overall SBM.net dataset.

Fig. 6 .Fig. 8 .

 68 Fig. 6. Qualitative results of the proposed method. (a) 8 images from the input sequences, (b) ground truth images, (c) estimated background model by the proposed MSCL, (d) RMAMR, (e) RMR, (f) DECOLOR, (g) 3TD, (h) MAGRPCA, (i) RFSA, (j) GRASTA, (k) GOSUS, (l) BRTF, and (m) GoDec. From top to bottom: each input sequence is selected from each category. (1)sequence '511' from 'Basic', (2) 'advertisementBoard' from 'Background Motion', (3) 'boulevardJam' from 'Clutter', (4) 'badminton' from 'Jitter', (5) 'AVSS2007' from 'Intermittent Motion', (6) 'CameraParameter' from 'Illumination Changes', (7) 'BusStopMorning' from 'Very Long', and (8) 'CUHK Square' from 'Very Short'.

  Cluttered Scene Average Value of Gray-level Error (AGE) (a) Simple Scene (background pixels are more dominant) Average Value of Gray-level Error (AGE) Linear Manifold (linear change in lighting condition) Average Value of Gray-level Error (AGE) (e) Non-linear Manifold (cluttered scene: abrupt changes) Average Value of Gray-level Error (AGE) Proposed without SSC (34.84) Spatiotemporal GD Regularizations (11.11) Spatiotemporal CC Regularizations (23.72) Proposed MSCL (10.19)

9 .

 9 Performance of estimated background model with only RPCA, spatial regularization, only temporal regularization, only correlation coefficient based spatiotemporal regularizations, only geodesic distance based spatiotemporal regularizations, and without spatiotemporal sparse subspace clustering.

Fig. 10 .

 10 Qualitative results of the proposed MSCL-FL algorithm on each category of CDnet 2014 dataset [50]. From left to right: the visual results of (a) Baseline, Camera Jitter, and Bad Weather categories, (b) Dynamic Background and Intermittent Object Motion categories, and (c) Thermal and Shadow category.

Fig. 11 .

 11 Fig. 11. Visual results of the proposed MSCL-FL algorithm over 9 sequences of BMC dataset [49]. From left to right: 9 videos are presented. From top to bottom: input frames, ground truth images, and visual results of the proposed algorithm are presented.

TABLE I COMPARISON

 I OF AVERAGE F 1 SCORE ON CDNET DATASET.

														GFL GoDec GRASTA GOSUS MSCL MSCL-FL
	Baseline	0.92	0.88	0.92	0.92	0.92	0.95 *	0.94	0.88	0.84	0.82	0.89	0.80	0.83	0.77	0.66	0.90	0.87	0.94
	Dynamic Background	0.70	0.75	0.81	0.78	0.71	0.81	0.89	0.75	0.55	0.84	0.82	0.74	0.74	0.58	0.35	0.79	0.85	0.90 *
	Camera Jitter	0.77	0.72	0.74	0.81	0.78	0.81	0.81	0.72	0.63	0.78	0.75	0.76	0.78	0.48	0.43	0.82	0.83	0.86 *
	Shadow	0.83	0.68	0.81	0.80	0.81	0.89 *	0.89 *	0.81	0.71	0.77	0.73	0.69	0.82	0.51	0.52	0.84	0.82	0.86
		0.70	0.78	0.81	0.76	0.75	0.81	0.83	0.79	0.69	0.79	0.75	0.80	0.76	0.62	0.42	0.80	0.82	0.86 *
	Intermittent Object Motion	0.59	0.55	0.54	0.65	0.67	0.65	0.77	0.69	0.57	0.80	0.66	0.63	0.59	0.38	0.35	0.74	0.80	0.84 *
	Bad Weather	0.76	0.79	0.82	0.75	0.79	0.86	0.81	0.77	0.78	0.75	0.70	0.79	0.76	0.66	0.68	0.77	0.83	0.88 *
	Average	0.75	0.73	0.77	0.78	0.77	0.82	0.84	0.77	0.68	0.79	0.75	0.74	0.75	0.57	0.48	0.80	0.83	0.88 *

Categories

DECOLOR 3TD DP-GMM 2P-RPCA LSD SuBSENSE PAWCS BMTDL TVRPCA SRPCA RMAMR LR-FSO

TABLE II COMPARISON

 II OF AVERAGE F 1 SCORE ON EACH VIDEO OF BMC 2012 DATASET[START_REF] Vacavant | A benchmark dataset for outdoor foreground/background extraction[END_REF] (SEE FIG.11) Videos DECOLOR 3TD DP-GMM LSD TVRPCA SRPCA RMAMR LR-FSO GFL GoDec GRASTA MSCL MSCL-FL TABLE III COMPARISON OF AVERAGE F 1 SCORE ON I2R [26] AND WALLFLOWER [48] DATASETS. BOLD FACE NUMBER REPRESENTS BEST PERFORMING METHOD.

	001	0.77	0.79		0.72	0.79	0.76	0.79	0.78		0.71	0.78	0.75	0.81		0.80	0.84 *
	002	0.71	0.76		0.69	0.80	0.67	0.74	0.71		0.66	0.74	0.73	0.85 *		0.78	0.84
	003	0.65	0.70		0.75	0.94	0.68	0.83	0.78		0.70	0.61	0.93	0.78		0.96 *	0.88
	004	0.85	0.83		0.80	0.88	0.82	0.81	0.79		0.72	0.88	0.83	0.76		0.86	0.90 *
	005	0.74	0.79		0.71	0.73	0.77	0.80	0.76		0.66	0.80	0.70	0.68		0.79	0.83 *
	006	0.80	0.82		0.68	0.80	0.69	0.69	0.65		0.78	0.74	0.78	0.59		0.74	0.80 *
	007	0.70	0.73		0.65	0.81	0.71	0.70	0.64		0.54	0.69	0.67	0.75		0.76	0.78 *
	008	0.78	0.81		0.78	0.84	0.79	0.84	0.80		0.80	0.81	0.76	0.84		0.89 *	0.85
	009	0.90	0.85		0.79	0.92	0.88	0.86	0.82		0.82	0.83	0.89	0.87		0.86	0.94 *
	Average	0.76	0.78		0.73	0.83	0.75	0.78	0.74		0.71	0.76	0.78	0.77		0.82	0.86 *
	Datasets	DECOLOR 3TD DP-GMM LSD TLSFSD TVRPCA SRPCA RMAMR LR-FSO MLRSBS MODSM GFL RFSA GRASTA MSCL MSCL-FL
	I2R	0.74	0.72	0.70	0.75	0.76	0.69	0.80	0.75	0.69	0.76	0.76	0.85 *	0.71	0.54	0.82	0.84
	Wallflower	0.59	0.75	0.78	0.75	N/A	0.61	0.85	0.80	0.74	N/A	0.73	0.84	0.54	0.33	0.86	0.92 *
	artifact effectively. The large margin obtained by the proposed
	algorithms is attributed to the removal of the motionless frames
	and the encoding of spatiotemporal regularization in the low-
	rank background model. The Bad Weather category contains
	six videos (Fig.

http://pione.dinf.usherbrooke.ca/results/

http://wordpress-jodoin.dmi.usherb.ca/results2014/

http://pione.dinf.usherbrooke.ca/results/

For simplicity, we use the same notation M ∈ R p×c for representing the motion mask of only dynamic frames.

ACKNOWLEDGMENTS

This work was supported by the National Research Foundation of Korea funded by the Korean Government (NRF-20170915).