Bugeau -Patrick Aurélie

Pérez

Aurélie Bugeau

Patrick Pérez

COM et COG INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET EN AUTOMATIQUE Detection and segmentation of moving objects in complex scenes

Keywords: motion detection, segmentation, mean shift clustering, graph cuts Détection détection de mouvement, segmentation, partitionnement mean shift, coupure de graphe

Detecting and segmenting moving objects in dynamic scenes is a hard but essential task in a large number of applications such as surveillance. Most existing methods only give good results in the case of persistent or slowly changing background, or if both the objects and the background can be characterized by simple parametric motions. This paper aims at detecting and segmenting foreground moving objects in the absence of such constraints. The sequences we consider have highly dynamic backgrounds, illumination changes and low contrasts, and can have been shot by a moving camera.

Three main steps compose the proposed method. First, moving points are selected within a sub-grid of image pixels. A descriptor is associated to each of these points. Clusters of points are then formed using a variable bandwidth mean shift with automatic bandwidth selection. Finally, segmentation of the object associated to a given cluster is performed using Graph cuts. Experiments and comparison to other motion detection methods on challenging sequences show the performance of the proposed method and its utility for video analysis in complex scenes.

Introduction

Detection of moving objects in sequences is an essential step for video analysis. It becomes a difficult task in the presence of a dynamic background. We are interested in very challenging sequences containing a complex motion in the background. This motion can also have a high amplitude. Furthermore, contrast between background and interesting objects can be small. Like in many applications in video analysis, the algorithms have to be robust to illumination and point of view changes. The last particularity about our sequences is that the camera taking the scene can move.

Existing methods

Different kinds of methods exist to solve the problem of motion detection and motion segmentation.

Good but incomplete reviews on motion detection methods can be found in [START_REF] Mitiche | Computation and analysis of image motion: a synopsis of current problems and methods[END_REF][START_REF] Radke | Image change detection algorithms: a systematic survey[END_REF]. Here, we divide these methods into four categories: thresholding, background modeling, layers extraction and finally saliency based methods.

Tresholding methods First works on motion detection were based on adjacent frames difference [START_REF] Jain | On the analysis of accumulative difference pictures from image sequence of real world scenes[END_REF]. The most obvious algorithm is to simply threshold the difference image. The choice of this threshold highly depends on the sequence, its noise, its motion. Furthermore there is no reason for this value to be constant on the whole image. Indeed, different objects and motions give different variations of luminosity. Many methods have been developed to decide whether or not a pixel has moved. The decision can be directly made independently on each pixel [START_REF] Konrad | Handbook of image and Video processing[END_REF] or on small blocks of pixels [START_REF] Fablet | Moving object detection in color image sequences using region-level graph labeling[END_REF]. With independent pixel-wise detections, detection maps are usually corrupted by holes in the mask of moving objects and false detections due to noise. These errors can be attenuated using regularization constraints and contextual information via Markov Random Fields [START_REF] Aach | Bayesian algorithms for change detection in image sequences using markov random fields[END_REF]. In [START_REF] Odobez | Separation of moving regions from background in an image sequence acquired with a mobile camera[END_REF],

Markov Random fields are applied after a step of motion compensation. This step is described in section 2. This method is well adapted to videos taken by a moving camera. The decision rule in many change detection algorithms is cast as a statistical test hypotheses [START_REF] Hsu | New likelihood test methods for change detection in image sequences[END_REF]. More complex methods are proposed in [START_REF] Rosin | Thresholding for change detection[END_REF] for modeling the spatial distribution of either noise or signal and selecting the appropriate threshold. In [START_REF] Veit | A maximality principle applied to a contrario motion detection[END_REF] an a contrario method is proposed. It is based on a perceptual grouping principle named the Helmholtz principle. It consists in defining an image model in the absence of moving objects instead of modeling the moving objects.

Background modeling and subtraction methods Methods based on adjacent frame difference are mostly sensitive to noise and to illuminations changes. When the number of frames in the sequence is high and there is not much change between consecutive frames, another solution to motion detection is background modeling. This technique is routinely used in the context of surveillance applications, when the camera is fixed. Background modeling methods can be classified as predictive or non predictive methods. Non-predictive methods build a probability density function of the intensity at an individual pixel. In static environment, the statistical distribution of a pixel can be represented by a single Gaussian [START_REF] Cavallaro | Video object extraction based on adaptive background and statistical change detection[END_REF] [26] [START_REF] Kanade | Advances in cooperative multi-sensor video surveillance[END_REF]. The foreground pixels are determined as those for which the intensity value is far from the mean background model and are clustered into objects. A variable number of Gaussian distributions corresponding to each different foreground object can be added. It was used by [START_REF] Hötter | Detection of moving objects using a robust displacement estimation including a statistical error analysis[END_REF] for generic objects and by [START_REF] Wren | Pfinder: Real-time tracking of the human body[END_REF] for people tracking. In the presence of dynamic background the use of a single Gaussian becomes inappropriate and a mixture of several Gaussians is preferred to model the background [START_REF] Friedman | Image segmentation in video sequences: A probabilistic approach[END_REF][START_REF] Grimson | Using adaptive tracking to classify and monitor activities in a site[END_REF]. When changes in the background are too fast, the variance of the Gaussians becomes too large and non parametric approaches are more suited. In [START_REF] Elgammal | Non-parametric model for background subtraction[END_REF], Gaussian kernels calculated on the past frames are used to model the density at a particular pixel. Contrary to previous approaches, this method addresses the uncertainty of spatial location. Until recently methods were almost all based on photometric properties. A lot of outdoor scenes exhibit a persistent INRIA motion which is well modeled by optical flow and a non parametric algorithm that combines color and flow features can be used [START_REF] Mittal | Motion-based background subtraction using adaptative kernel density estimation[END_REF]. As optical flow can not be computed when there is no intensity gradient, the authors have chosen to use kernels with variable bandwidth. In [START_REF] Ren | Statistical background modeling for non-stationary camera[END_REF] the authors extended a statistical background modeling technique to cope with non stationary camera. The current image is registered to the estimated background image using an affine or projective transformation. The foreground information can also be used as in [START_REF] Paragios | Adaptive detection and localization of moving objects in image sequences[END_REF] in which the background and foreground maps are forced to be a Markov random field. All these pixel-wise approaches allow an accurate detection of moving objects but are memory and possibly computationally expensive. Also, they can be sensitive to noise and they do not take into account spatial correlation. Spatial consistency can be added [START_REF] Sheikh | Bayesian modeling of dynamic scenes for object detection[END_REF] with a MAP-MRF modeling of both foreground and background. This method has been extended to novelty detection in [START_REF] Mahamud | Comparing belief propagation and graph cuts for novelty detection[END_REF].

Predictive methods use a dynamical model to predict the value of a pixel from previous observations. A Kalman filter based approach that models the dynamics of the intensity at a particular pixel can be used [START_REF] Karmann | Time-varying image processing and moving object recognition[END_REF][START_REF] Koller | Robust multiple car tracking with occlusion reasoning[END_REF]. In [START_REF] Toyama | Wallflower: Principles and practice of background maintenance[END_REF] an algorithm called wallflower is described. It uses a simpler version of the Kalman filter called Weiner filter to predict a pixel's current value from its k previous values.

Pixels whose prediction error is high are classified as changed pixels. Recent methods are based on more complicated models. For example, in [START_REF] Doretto | Dynamic textures[END_REF] and [START_REF] Zhong | Segmenting foreground objects from a dynamic textured background via a robust kalman filter[END_REF], an autoregressive model was proposed to capture background properties.

Background subtraction and tresholding methods are a preliminary step to moving object detection and subsequent processing is necessary to get the masks of moving objects.

Layer approaches Motion segmentation can be seen as the problem of fitting a collection of motion models to the spatio-temporal image data. This leads to the layer approach [START_REF] Darrel | Robust estimation of a multi-layered motion representation[END_REF] that tries to fit a mixture of motion models to the entire image. Layers are then found by associating each pixel to the model it belongs to. In many papers [START_REF] Ayer | Layered representation of motion video using robust maximumlikelihood estimation of mixture models and mdl encoding[END_REF][START_REF] Jepson | Mixture models for optical flow computation[END_REF][START_REF] Wang | Representing moving images with layers[END_REF][START_REF] Weiss | Smoothness in layers: Motion segmentation using nonparametric mixture estimation[END_REF] a mixture of probabilistic models is iteratively built with an Expectation-Maximization algorithm (EM). A major drawback of such approaches is that they are very sensitive to the initialization and are computationally expensive. In [START_REF] Xiao | Accurate motion layer segmentation and matting[END_REF] graph cuts have been used to extract these models or layers. After a number of seed regions are determined using two frames correspondences, these seed regions are first extended thanks to a graph cuts segmentation method. The resulting initial regions are then merged into layers according to motion similarities. This method requires the scene to be induced by multiple planar regions having an apparent affine motion. Other approaches aim at fitting a polynomial model to all the image measurements and then factorizing this polynomial to obtain the parameters of each 2D motion models (multi-body factorization) [START_REF] Vidal | A unified algebraic approach to 2-d and 3-d motion segmentation[END_REF]. It was adapted to both static and dynamic scenes [START_REF] Vidal | A closed form solution to direct motion segmentation[END_REF]. Recently, in [START_REF] Pundlik | Motion segmentation at any speed[END_REF], an incremental approach to layer extraction has been introduced. Feature points are detected, tracked and then merged into groups based on their motion. Objects are detected incrementally when enough evidence allows them to be distinguished from their background.

In [START_REF] Kim | Background updating for visual surveillance[END_REF] a combination of background modeling with a layer technique is proposed. The layers are called short-term backgrounds. The idea is to assign a layer to each moving object and to keep this object as a single layer even when it stops moving.

Saliency based methods

A last approach is to define moving objects as areas having salient motion.

In [START_REF] Wildes | A measure of motion salience for surveillance applications[END_REF], salient motion was defined as motion that is likely to result from a typical surveillance target (e.g. person, vehicle). This definition was used in [START_REF] Wixson | Detecting salient motion by accumulating directionally-consistent flow[END_REF] and [START_REF] Tian | Robust salient motion detection with complex background for real-time video surveillance[END_REF] to detect moving objects. The assumption made is that an object with salient motion moves in an approximate consistent direction during a time period. Therefore moving objects are searched as localized image regions that have moved in the same direction during a time period. In [START_REF] Tian | Robust salient motion detection with complex background for real-time video surveillance[END_REF], the accumulation of flow motions is done during 10 frames. A fusion of background modeling and saliency was proposed in [START_REF] Zhu | Learning a sparse, corner-based representation for timevarying background modeling[END_REF]. A specificity in this paper is that the background is only sparsely modeled on corners, and moving objects are then found by the clustering of foreground corner trajectories.

Overview of the paper

In this paper, we are interested in challenging sequences containing complex motions, with possible high amplitude and sudden changes in the background. For example, in the context of driver surveillance, the motions visible through the windows are often hard to characterize. The "background" is composed of both the passenger compartment and what is behind the windows. Furthermore, contrast between background and interesting objects (face, hands) can be low. Also, the motion of the "interesting" moving objects can be close to the one of the moving background. For example, on driver sequences, depending on the speed of the car, an arm going from the steering wheel to the face can have the same speed as some trees behind the window. Finally, the sequences we consider can be shot by a moving camera. Few frames of such a sequence are shown in figure 1. Most existing methods would fail to detect only the moving arm because of all the motions present behind the window and the low contrast between the harm and the guardrail. Our work does not aim at modeling the background or at finding every layer but only at detecting moving foreground objects. We define these objects as groups of pixels salient for both motion and color. Our algorithm can be divided in three main steps. First, in section 2, the camera motion is computed and a subgrid of "moving" pixels, i.e not belonging to camera motion, is selected. A descriptor is defined to characterize each points of this grid. They are then merged into clusters consistent for both color and motion (section 3), using mean shift filtering algorithm [START_REF] Comaniciu | Mean shift: A robust approach toward feature space analysis[END_REF]. From the clusters, the complete pixel-wise segmentation of moving objects is found using a MAP-MRF framework (section 4). Finally, section 5 presents some experimental results.

t = 15 t = 16 t = 17

Point selection and description

In this paper, similarly to [START_REF] Wixson | Detecting salient motion by accumulating directionally-consistent flow[END_REF] and [START_REF] Tian | Robust salient motion detection with complex background for real-time video surveillance[END_REF], the moving objects detection will only be performed on several points and their neighborhood. We have preferred the computational lightness and the noise robustness of these methods to the accuracy of the pixel-wise approaches. The first step of our algorithm is then to build a sub-grid of moving pixels and to select a descriptor for these pixels. This section is organized as follows. First we explain what are the moving pixels. Next we discuss the construction of the grid of points and finally we present the selected descriptors.

First, let us introduce the notations. In all the paper, P will denote the set of N pixels of the frame I t at time t from an input sequence of images. To each pixel s = (x, y) of P is associated a feature vector z t (s). In case of color sequences, z t (s) = (z

(G) t (s), z (C) t (s), z (M) t (s)),
(G) t (s), z (M) t (s)).

Sensor motion

Moving pixels are the pixels not belonging to the camera motion. We only work on moving pixels because we aim at detecting moving objects. We assume that the apparent motion induced by the physical motion of the camera is dominant in the image (among various movements in the image field, it is the one that concerns the larger number of pixels) and is well approximated by a 2D affine motion field. The parametric flow vector w θ (s) at location s = (x, y) reads:

w θ (s) = a 1 a 4 + a 2 a 3 a 5 a 6 . x y . (1)
Different methods are available to estimate the parameters of such a model [START_REF] Black | The robust estimation of multiple motions: parametric and piecewise-smooth flow fields[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF][START_REF] Sawhney | Compact representations of videos through dominant and multiple motion estimation[END_REF]. We use the robust real-time multiresolution algorithm described in [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]. The parameter vector θ = (a 1 , a 2 , a 3 , a 4 , a 5 , a 6 ) is estimated between two consecutive frames I t+1 and I t as follows:

θ = arg min θ s ρ z (G) t+1 (s + w θ (s)) -z (G) t (s) + ζ t , (2)
where ρ(x) is an M-estimator and ζ t is a global intensity shift that accounts for global illumination changes. The minimization is done through a multiscale Gauss-Newton method that yields a succession of reweighted least-squares problems. The auxiliary weight maps of the M-estimator will be denoted as W t (W t (s) ∈ [0, 1]). The final map indicates if a pixel participates to the final robust motion estimate (W t (s) close to 1) or is more considered as an outlier (W t (s) close to 0). A simple pixel-wise motion detector can be built using this map. A pixel is considered as "moving" at time t if it is an outlier to the dominant motion at times t and t -1:

M t (s) = 1 if W t (s + w θ,t-1 (s))) + W t-1 (s) = 0, 0 else . (3)
If M t (s) = 0, pixel s is considered as a motionless pixel, and it will not be used for the clustering step of the algorithm. The choice of pure outliers to dominant motion for moving pixels can seem drastic. However experiments have shown that no moving object is lost using such method. This choice has been made to avoid false labeling of pixels. Furthermore, it permits to deal with occlusion and disocclusion of the scene background. Figure 2 shows a result of this pixel-wise motion detector.

It is easy to see that this preliminary step is not sufficient to detect interesting moving objects. In this sequence, we are willing to detect the water skier. The motion in the water is complex, not repetitive, and so can not be estimated directly. A drawback of dominant motion estimation can appear in the presence of large uniform objects moving slowly. Inner portions of such objects often appear to belong to the dominant motion (see the person on figure 3). Our algorithm will handle this loss of moving pixels thanks to the segmentation step where the whole pixel grid is considered (see section 4).

Grid construction

The goal of the algorithm is to build groups of pixels consistent both for motion and for some photometric or colorimetric features. These groups must correspond to interesting moving objects. In [START_REF] Zhu | Learning a sparse, corner-based representation for timevarying background modeling[END_REF], moving objects are found using corners, detected with the Harris corner detector. The authors justify the use of corners by claiming that a moving object contains a large number of corners. In our experiments, we have observed that the number of corners belonging to a moving object can be much lower than the number of corners belonging to the background (figure 4). Besides, if variations in the background are fast and if parallax changes, the number of corners and their neighborhoods can be signicatively different from one frame to the other. Finally, corner detection adds one stage of calculation and requires two thresholds. As no a priori is assumed on the shape and texture of objects, we have chosen to use points of arbitrary type. Hence, we only use a grid of points regularly spread on the image. As the purpose is to detect moving objects, the simple pixel-wise motion detector from previous subsection is used to restrict this step to the grid subset: where w and h are the dimensions of the image and N 2 G the size of the grid before pruning. We have arbitrarily chosen to use the same number of points in the x and y axes. The value of the parameter N G is important. It controls the balance between computational cost (regional methods) and accuracy (local methods). Next step of the algorithm can become computationally expensive if the number of points of the grid is too large. An important thing to note is that N G may depend on the number m of "moving" pixels in the image, m = s∈P Mt(s). To limit the computational cost for clusters creation, we fix the number of moving points to n G (500 in our experiments) that will be approximatively kept in further steps of the algorithm. The size N G of the grid is then set as

INRIA (a) (b) (c) (d) (e) (f)
G = s = k.w N G , l.h N G , k = 0 • • • N G , l = 0 • • • N G | M t (s) = 1 , (4)
N G = ⌊ w * h * n G /m⌋,
where ⌊•⌋ is the integer part of •.

Description of the selected points

Now that the points are chosen, the features that will be used to create clusters corresponding to objects need to be defined. It is necessary to chose only few discriminant features. An object is defined as a moving and compact area over which the values of displacement and photometry are nearly constant. Color is not sufficient because the contrast between an object and the background can be small, nor is flow in case of similar motion between an object and the background. Hence the descriptor is formed by three different groups of features. The first group is composed of the coordinates of the point. The second group contains its motion, and the last one contains discriminant photometric features.

Motion features

As we try to detect moving objects, an essential feature is the displacement of the selected points.

Several types of methods exist for computing optical flow. Good reviews are given in [START_REF] Barron | Performance of optical flow techniques[END_REF][START_REF] Beauchemin | The computation of optical flow[END_REF]. We concentrate on the Lucas and Kanade algorithm [START_REF] Lucas | An iterative technique of image registration and its application to stereo[END_REF], with an incremental multiscale implementa-INRIA tion. The reasons for such a choice are that we want an efficient approach and that we do not aim at computing the flow on the whole image but only on several points (which are not particular feature points). In order to use classical notations, in this subsection I denotes the intensity function, I(x, y, t) is the intensity value at point s = (x, y) and at time t, and d x and d y are the horizontal and vertical components of the apparent displacement between instants t and t + 1. The initial assumptions are that image intensities are approximatively constant under motion for at least a short period of time and that motions are small:

I(x + d x , y + d y , t + 1) ≃ I(x, y, t) . (5)
Using a Taylor approximation, previous equation leads to the brightness constancy equation:

∂I ∂x d x + ∂I ∂y d y + ∂I ∂t = 0 . (6)
This equation has two unknowns and is thus insufficiently constrained. This is the aperture problem:

edges parallel to the motion do not convey motion information. Only the motion component in the direction of the local gradient of the intensity function can be computed. To find the optical flow some additional constraint must be added. The solution of Lucas and Kanade is to use a local neighborhood V(x, y) and to assume that the flow is constant in this small region of size n. This leads to an overdetermined system of equations solved using a least squares method. The flow is then obtained by solving:

argmin (dx,dy) (xi,yi)∈V(x,y) ∂I ∂x (x i , y i , t)d x + ∂I ∂y (x i , y i , t)d y + ∂I ∂t (x i , y i , t) 2 . ( 7
)
Unfortunately the associated 2x2 linear system is regular and well conditioned only for "textured" neighborhoods. Therefore we will not use the points for which the intensity gradient is very low in all the region. These points will not be considered in next steps of the algorithm and are therefore left apart from the grid:

G = s = k.w N G , l.h N G | M t (s) = 1 & ∃(x i , y i ) ∈ V k.w N G , l.h N G , |∇I(x i , y i )| = 0 . (8)
The brightness constancy equation is only valid for small motion. In order to access large motions, an incremental multiscale approach can be used. A Gaussian pyramid is derived from each image.

The optical flow constraint holds for coarsest scales of the pyramids even when the motion is large.

Estimated motion for large scales are used as initial estimates for lower scales.

The whole procedure can be applied, independently, to the different locations of any pixel set. No spatial consistency is enforced over estimates at neighboring locations. We could instead have used Horn and Schunk algorithm [START_REF] Horn | Determining optical flow[END_REF] that adds a smoothness term to regularize over the whole image or the robust estimation of Black and Anandan [START_REF] Black | A framework for the robust estimation of optical flow[END_REF] to get a better estimation. However these algorithms are more computationally expensive and we do not aim at having a dense and very precise estimation over the whole image.

The computation of flow using gradient based methods fails when the brightness constancy is not satisfied and when a point does not move as its close neighbors. Indeed, the two basic assumptions are not valid in such cases. Instead of using a more robust method to compute optical flow, we have [START_REF] Neyman | On the use and interpretation of certain test criteria for purposes of statistical inference[END_REF] argued that rather than having a single hypothesis that should be rejected or not, it is better to choose between two hypotheses. One of them is the null hypothesis and the other one is the alternative hypothesis. The goal of statistical tests is then to decide if the null hypothesis should be rejected. An alpha risk can be attached to the null hypothesis. This risk is defined as the risk of rejecting the null hypothesis when in fact it is true. It is stated in terms of probability and corresponds to the confidence level of a statistical test. In our case, the null hypothesis H 0 asserts that the two data samples X and Y are uncorrelated.

chosen
A statistical test related to the linear correlation is the so-called "p-value" (or observed significance level). To get the p-value, the T-statistics (or Student's statistics) has to be computed:

T = γ n -2 1 -γ 2 . ( 9
)
The parameter n is the size of the sample X and Y and γ the normalized cross correlation. A Tstatistics near zero is evidence under the null hypothesis that there is no correlation between the two samples. Assuming that H 0 is true, T is well represented by a Student's distribution defined by:

A(x) = 1 √ n -2 B( 1 2 , n-2 2 ) x -x 1 + y 2 n -2 3-n 2 dy , ( 10
)
where B is the beta function given by:

B(a, b) = B(b, a) = 1 0 x a-1 (1 -x) b-1 dx . (11)
The p-value is finally the probability, when the null hypothesis is true, that the absolute value of the T-statistics would exceed the alpha risk. It is equal to A(|T |). A small p-value means that the null hypothesis is false and that the two data samples are in fact correlated. If one wants to limit to 5% the alpha risk, then data are assumed correlated if the p-value is lower than 0.05. For more information on the theory and the implementation of the p-value, we refer to [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF].

If the p-value obtained for a point s is larger than 0.05, the flow for this point is considered as non valid. The point s will then not be processed in next steps of the algorithm. Finally, recalling that s ′ = (x + dx, y + dy), a new grid

G = s = k.w NG , l.h NG | Mt(s) = 1 & ∃(xi, yi) ∈ V k.w NG , l.h NG , |∇I(xi, yi)| = 0 & p-value(s, s ′ ) < 0.05 (12)
is obtained with a flow vector z The influence of the p-value is shown on figure 5. The first row shows the original image and the set of moving pixels. The second row shows the final grid obtained without using a validation step (Figure 5c), when keeping only vector flows for which the correlation coefficient is higher than 0.5 (Figure 5d) and the final grid obtained with the pvalue test (Figure 5e). The third row presents the associated motion fields, as estimated by multiscale

(a) (b) (c) (d) (e) (f) (g) (h)
Lucas-Kanade technique at each point of the various grids. The parameters that yield these grids are the following: the number of moving pixels m is 10845, the size of the image is 240x320 and the parameter N G is equal to 5. The grid finally contains 420 pixels when no validation is used, 184 pixels for a correlation test and 277 with the use of p-value. The p-value enables to keep more points than the correlation, especially important points on the arm of the driver.

Photometric features

The last features concern the photometry at selected locations. They are different for grayscale and color sequences. To be robust to noise, the features are computed over the neighborhood of each point of the grid defined in previous subsection.

Grayscale sequences:

In case of grayscale sequences, the first photometric feature is the intensity itself, and more precisely the mean z (G) t (s) of the luminance on a 3x3 window around the point s = (x, y). As the contrast between the object and the background can be small, this feature is not sufficient. Another interesting information is the texture. A lot of definitions and characterizations of texture can be found in literature. Here it will be simply associated to the quantity or the strength of edges present in the studied area. The feature that we introduce to capture the texture is the standard deviation σ ∆z (G)

t (s)
to the mean of the Laplacian of intensity on a 3x3 window around the point. To include some simple temporal consistency, we add image intensity and texture values at time t + 1 for the displaced point s ′ = s + (d x , d y ). Finally, for grayscale sequences, the descriptor at each individual valid point s = (x, y) of the grid, indexed by i (i = 1 . . . |G|), is:

x (i) = (x (i) 1 , x (i) 2 , x (i) 3 ), (13)
where

x (i) 1 = (x, y) , x (i) 2 = (d x , d y ) , x (i) 3 = (z (G) t (s), σ ∆z (G) t (s) , z (G) t+1 (s ′ ), σ ∆z (G) t+1 (s ′ ) )
, and s ′ = (x + dx, y + dy).

Color sequences:

For color sequences, the information given by the three channels have appeared to be sufficient during our experiments. Hence, no texture information will be added. However, we observed that the three color channels RGB do not give the best representation of our images. Indeed, they are highly correlated and this representation is sensitive to illumination changes. Therefore it is better to separate the luminance and the chrominance informations. Furthermore, most of our test sequences contain human skin, which has a specific signature in the space of chrominance [START_REF] Kjeldsen | Finding skin in color images[END_REF][START_REF] Singh | A robust skin color based face detection algorithm[END_REF]. Therefore, it is interesting to use a color system representing the chrominance instead of the classical RGB system. Most chrominance spaces give the same results on skin detection [START_REF] Terrillon | Comparative performance of different chrominance spaces for color segmentation and detection of human faces in complex scene images[END_REF]. In this paper we use the YUV color space. This choice proves appropriate for various types of sequences. Here again, to include some simple temporal consistency, we add image t + 1 chrominance values of the corresponding point.

Finally, for color sequences, the descriptor at each individual valid point indexed by i (i = 1 . . . |G|) of the grid is:

x (i) = (x (i) 1 , x (i) 2 , x (i) 3 ), (14)
where

x (i) 1 = (x, y) , x (i) 2 = (d x , d y ) , x (i) 3 = (z (C) t (s), z (C) t+1 (s ′ )) ,
where the 3-dimensional color feature vector is: z (C) t (s) = (Yt(s), Ut(s), Vt(s)).

Clustering points: mean shift for mixed feature spaces

Now that a grid of valid points has been chosen and described, we address the problem of grouping the points into clusters. Many data clustering methods have been described in the literature. A good review on classic clustering techniques can be found in [START_REF] Jain | Data clustering: a review[END_REF]. An appealing technique to extract the clusters is the mean shift algorithm, which does not require to fix the (maximum) number of clusters. On the other hand the kernel bandwidth and its shape have to be chosen or estimated for each dimension.

Fixed bandwidth mean shift partitioning

Mean shift is an iterative gradient ascent method used to locate the density modes of a cloud of points, i.e. the local maximum of its density. This technique, which we now summarize, is well described in [START_REF] Comaniciu | Mean shift: A robust approach toward feature space analysis[END_REF].

Given a set of n points {x (i) }i=1..n in the d-dimensional space R d , the non-parametric density estimation at each point x is given by:

f (x) = 1 n|H| 1/2 n i=1 K(H -1/2 (x -x (i) )) = c k n|H| 1/2 n i=1 k( H -1/2 (x -x (i) ) 2 ) (15)
where K is a kernel with associated profile k, H is the bandwidth matrix and c k is a strictly positive normalization constant which makes K(x) integrate to one. Introducing the notation

g(x) = -k ′ (x)
the density gradient reads:

∇ f (x) = H -1 f (x) m(x) ( 16
)
where m is the "mean shift" vector,

m(x) = n i=1 x (i) g H -1/2 (x -x (i) ) 2 n i=1 g H -1/2 (x -x (i) ) 2 -x . ( 17
)
Using exactly this displacement vector at each step of an iterative search guaranties convergence to the local maximum of the density [START_REF] Comaniciu | Mean shift: A robust approach toward feature space analysis[END_REF]. With a d-variate Gaussian kernel, equation [START_REF] Elgammal | Non-parametric model for background subtraction[END_REF] becomes

m(x) = n i=1 x (i) exp(-1 2 D 2 (x, x (i) , H)) n i=1 exp(-1 2 D 2 (x, x (i) , H)) -x (18)
where

D 2 (x, x (i) , H) ≡ (x -x (i) ) T H -1 (x -x (i) ) (19)
is the Mahalanobis distance from x to x (i) .

A mode seeking algorithm, or mean shift filtering can be derived by iteratively computing the mean shift vector. The final partition of the feature space is obtained by grouping together all the data points that converged to the same mode.

INRIA

Variable bandwidth mean shift

Usual mean shift procedures use a fixed bandwidth for all the data set. However, variable bandwidths are essential when local characteristics of the data vary significantly across the feature space. Two density kernel estimators have been introduced to take into account a bandwidth that is not fixed for each point [START_REF] Wand | Kernel Smoothing[END_REF]. The first one is called "balloon estimator" and makes the bandwidth vary at each estimation point. It was first introduced by Loftsgaarden and Quensberry [START_REF] Loftsgaarden | A nonparametric estimate of a multivariate density function[END_REF]. It is defined as:

f (x) = 1 n n i=1 K H(x) (x -x (i) ) = 1 n n i=1 1 |H(x)| 1/2 K(H(x) -1/2 (x -x (i) )) (20)
where H(x) is the bandiwth matrix varying at each estimation point. When applied globally, this estimator typically does not integrate to 1 and thus is usually not itself a density, even when K is.

In [START_REF] Terrell | Variable kernel density estimation[END_REF] the authors have investigated the degree of improvement that this estimator allows over fixed kernel estimates. For data up to 3 dimensions, the improvement seems to be very modest. However the balloon estimator becomes very efficient as soon as the number of dimensions becomes larger than 3.

The second density kernel estimator is the "sample point estimator" that makes the bandwidth vary at each data point. In [START_REF] Terrell | Variable kernel density estimation[END_REF] the advantages and drawbacks of this estimator have been studied. The major advantages are that it is a density and that a particular choice of the bandwidth can considerably reduce the bias [START_REF] Hall | Improved variable window kernel estimates of probability densities[END_REF]. However, finding this value for multivariate data is a hard problem not yet solved.

A disadvantage is that the estimate at a point may be influenced by observations very far away and not just by points nearby. In [START_REF] Terrell | Variable kernel density estimation[END_REF] simulations have shown that this estimator has a very good behavior for small-to-moderate sample sizes, but deteriorates in performance compared to fixed estimates as the sample size grows.

In this paper we are working on data having more than 3 dimensions. Therefore, we will concentrate only on the balloon estimator. The mean shift technique using the balloon estimator has been introduced in [START_REF] Bugeau | Bandwidth selection for kernel estimation in mixed multi-dimensional spaces[END_REF] for univariate and multivariate data.

Bandwidth selection for mean shift

A mode seeking algorithm can be derived by iteratively computing and applying the mean shift vector.

The partition of the feature space is obtained by grouping together all the data points whose associated mean shift procedures converged to the same mode. The quality of the results highly depends on the choice of the bandwidth matrix H.

The bandwidth selection can be based on statistical analysis or task-oriented. Statistical methods compute the best bandwidth by balancing the bias and the variance of the density estimate. Taskoriented methods rely on the stability of the feature space partitioning. In [START_REF] Comaniciu | An algorithm for data-driven bandwidth selection[END_REF] a method that combines the two different approaches has been introduced. This method is based on the maximization of the normalized mean shift vector. However, there is no theoretical result permitting to relate formally this maximization to the quality of the clustering. In [START_REF] Bugeau | Bandwidth selection for kernel estimation in mixed multi-dimensional spaces[END_REF], a method based on the task oriented step of [START_REF] Comaniciu | An algorithm for data-driven bandwidth selection[END_REF] has been developed. It is dedicated to the variable bandwidth selection in mixed multi-dimensional spaces. Here we briefly remind its principle.

The main idea underlying this bandwidth selection method is that if one cluster can be represented by a normal distribution, then the best cluster is the one for which the normal distribution is the most stable. If few points are added to the cluster or if some are left apart, the associated distribution should not change. Therefore by partitioning the data for several predefined bandwidths, one can identify the best bandwidth as the one that gave the most stable partitions.

Assume that the d-dimensional data can be decomposed as the Cartesian product of P independent spaces associated to different types of information (e.g. position, color), also called feature spaces or domains, with dimension dρ, ρ = 1 . . . P (where P ρ=1 dρ = d). The method in [START_REF] Bugeau | Bandwidth selection for kernel estimation in mixed multi-dimensional spaces[END_REF] proposes to find iteratively for each feature space ρ the best bandwidth Υ (i) ρ for each estimation point i. The iterative estimation of bandwidths is shown in algorithm 1. The final partition of the data is obtained by applying a last time the mean shift partitioning using the balloon estimator with the selected variable bandwidths

Υ (i) = diag[Υ (i) ρ , ρ = 1 . . . B].
The last point to discuss is how to choose the range of predefined bandwidths for our application of moving object detection. It is obvious that this choice influences the results. One could use a large sample of bandwidths but this would be very expensive as the number of mean shift partitioning procedures would be very large. In our application of moving object detection the data space is divided in P = 3 feature spaces: the position, the motion and the color. In all our experiments we have selected the same predefined bandwidth matrices. For all the feature spaces we have taken Bρ = 9 bandwidths. For the position they depend on the size of the grid:

H (b) 1 = w N G (1 + 3b B 1 -1 ) H (b) 2 = h N G (1 + 3b B 1 -1 ) , b = 1 . . . B 1 , (26)
where B1 is the number of predefined bandwidths for the first feature space (here B1 = 9). The range of predefined matrices for color and motion is directly computed from images noises. We define color, respectively motion, noises as the standard deviation over the whole grid of the color values difference, respectively motion, between two neighboring points of the grid. Introducing VG the set of pairs of neighboring points of the grid G, |VG | its cardinal, the mean and standard deviation vectors are:

α ρ = 1 |V G | (i,j)∈VG |x (i) ρ -x (j) ρ | (27)
and

β ρ = 1 |V G | (i,j)∈VG (|x (i) ρ -x (j) ρ | -α ρ ) 2 . ( 28
)
Some experimentations have shown the range of predefined matrices for color (ρ = 2) and motion (ρ = 3) can be defined as

H (b) ρ = β ρ (0.5 + 1.5b B ρ -1 )I dρ , b = 1 . . . B ρ , (29)
where I dρ is the identity matrix of size dρ. The clustering step decomposes the data into several clusters, each corresponding to a moving object or a moving part. We retain only large enough clusters (e.g., with more than 15 grid points). Finally we obtain k t clusters Cu,t, u = 1 . . . kt.

INRIA 4 Segmenting the objects

In order to get the complete masks of objects, a final step is necessary. Segmenting the object associated to a given cluster amounts to assigning a label l s , either "background" or "object", to each pixel s of the image. This problem can be reformulated into the graph cut framework as a bi-partitioning problem. In recent years graph cuts have been increasingly used in image segmentation. The reason for such a popularity is that the exact maximum a posteriori (MAP) of a two label pairwise Markov Random Field (MRF) can be computed in polynomial time using min-cut/max-flow algorithms [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF].

In seminal paper [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], Boykov et al. introduce an iterative foreground/background segmentation system based on this principle, using hard constraints provided by the user. Here we can directly learn some properties of the object from the points belonging to its cluster. These points are called inliers.

The energy function is defined so that its minimum should correspond to a good segmentation, in the sense that it is guided both by the motion and color of the inliers (observed foreground) and by distributions of color and motion built on the whole image. Also, boundaries of the segmentations should preferably coincide with large photometric gradients. These various specifications are captured by the following objective function:

E t (L) = -γ c s∈P ln(Pr(z (C) t (s)|l s )) -γ m s∈G ln(Pr(z (M) t (s)|l s ))+ λ (s,r)∈V exp - z (G) t (s) -z (G) t (r) 2 σ 2 1 dist(s, r) (1 -δ(l s , l r )) (30)
where L is the set of all the labels ls, s ∈ P, dist is a distance measure, and V is the set of unordered pairs (s, r) of neighboring elements of P. The parameters γm, γc, λ are some weight constants discussed below. Since we are segmenting each object independently, there is one such energy function to be minimized per cluster.

The two first terms of the cost function are based on pixel-wise modeling of color and motion features distributions. Motion term only concerns the points of the grid. For both color and motion, the object distribution is a mixture of Gaussians on the inliers. For the background, the mixture is built as follows. For color it is computed on the whole image whereas for motion it is only computed on the grid. In [START_REF] Boykov | Interactive graph cuts for optimal boundary and region segmentation of objects in n-d images[END_REF], authors have shown that it is possible to force some pixels to belong to the object or to the background. Here we force inliers to belong to the object. This is done by rewriting the energy function as:

E t (L) = -γ c s∈P\Cu,t ln(Pr(z (C) t (s)|l s )) -γ m s∈G\Cu,t ln(Pr(z (M) t (s)|l s ))- s∈Cu,t 1 + max s ′ ∈P λ r|(s ′ ,r)∈V V {s ′ ,r} δ(l s , "object") + λ (s,r)∈V V {s,r} (1 -δ(l s , l r )) (31)
where u is the object or cluster we are trying to segment and V {s,r} = exp -

z (G) t (s)-z (G) t (r) 2 σ 2 1 dist(s,r) .

Bugeau & Pérez

Because for motion we only consider points of the grid, we chose to set the parameters γc and γm such that γc = 1 and

γ m = N 2 G . (32)
This choice permits to give more influence to the points having a valid motion. If γm and γc were equal, the motion would have only a very small influence to the segmentation result.

The parameter σ in the third energy term can be related to noise [START_REF] Rother | grabcut": interactive foreground extraction using iterated graph cuts[END_REF]:

σ = 2 * (z (G) t (s) -z (G) t (r)) 2 (33)
where . denotes expectation over the whole image. The value of parameter λ has not been really studied in literature. To avoid a possible saturation of all binary edges in the max-flow procedure, we fix here its value as:

λ = 1 N -γ c
s∈P ls="bkg","obj" ln(Pr(z

(C) t (s)|l s )) -γ m
s∈G ls="bkg","obj" ln(Pr(z

(M) t (s)|l s )) . (34)
After the minimization, all pixels labeled as "object" form the final mask of the moving object.

INRIA

Experimental results

This section presents several experimental results on different kinds of sequences. These videos have been shot by a moving or a non moving camera, and they exhibit rather different types of motions.

Also we tested our method on both grayscale and color sequences. For all the tests, exactly the same parameters were used. The size of the grid as well as the predefined matrices for bandwidths selection are set as detailed in the paper. This section is decomposed in two parts. First we present the results obtained for both the clustering and the segmentation step. In a second part we see the final output (set of all segmented objects) of our algorithm as a motion detection mask in order to compare our approach with other methods mainly based on background modeling.

Sequences taken by a moving camera

We start by showing the results of the clustering and the segmentation on three different video sequences. For visualization purpose, an arbitrary individual color is assigned in each frame to each segmented object. This color only depends on the arbitrary order in which the objects were handled.

Hence there is no temporal consistency and the same object can be represented by different colors along the sequence.

The two following results are on color sequences shot by a moving camera. The first one is the water skier sequence already discussed in section 2. This sequence is hard because a lot of moving pixels are present in the water. On figure 6 we show the original images (first column), the moving pixels maps (second column), the results of the clustering steps (third column) and finally the segmented moving objects (fourth column). The water skier is most of the time well detected and segmented while no objects are found on the water. On the frame at time t = 124, the skier is not detected since, at this instant, his apparent motion is similar to the dominant motion estimated on the scene (see the mask of moving pixels). Note also that on the same frame, part of the water was detected as a moving object at the end of the clustering step. This happened several times during the sequence. However most of the time there is no corresponding moving object. The reason is that the water occupies a large part of the image, making the probability to belong to the background almost the same as the one to belong to the object. At the end of the energy minimization, all the pixels in the water are labeled as background, except for the inliers that were forced to belong to the object. These inliers are only visible by zooming on the segmentation image. In frame number 214, the cluster covers both the skin (arm, face, legs) and the body of the skier. This comes from the automatic bandwidth selection for the mean shift partitioning. The chrominance values corresponding to these different parts are very close. However it is not the case for the Y channel, corresponding to the intensity. The range of predefined bandwidths we are using is large, and the preferred bandwidth for the intensity turned out to be large.

The next sequence is the difficult driver sequence presented in the introduction. In this sequence we have to deal with the high dynamic background (behind the window), the low contrast between the shirt and the guardrail, some drastic illumination changes. Furthermore the motion of the hand and the trees behind the window are sometimes rather similar. The clusters corresponding to the moving objects and the segmentations can be seen on figure 7. Despite the number of moving pixels found behind the window, no objects are detected in this part of the image. The hand is very well detected in the three frames shown. However the arm is not detected in the last frame because not enough moving pixels with similar motion are found. Unfortunately a part of the passenger compartment and the guardrail are detected as moving objects. The clusters found there were small but as there is no strong enough contours, the segmentation spills over. We believe that detecting this area is not such a big problem. Indeed, it is not moving during several successive frames. Hence, adding some temporal consistency and/or tracking that would reject static objects would probably allow the removal of such spurious detection. On figure 8, we show another part of this sequence with a different background.

The results obtained, although not perfect, are rather encouraging for such a difficult sequence.

The last results presented in this subsection were obtained on a grayscale sequence with a non moving camera showing some cars and pedestrians moving in a street. The difficulty here comes from the high level of noise and the low contrast present along the sequence. Furthermore, there is a large number of small moving objects. If, for these reasons, the segmentation (not shown) are disappointing, the detected clusters are nonetheless interesting. In figure 9 we show the bounding boxes corresponding to these clusters. Another difficulty faced by the segmentation part lies in the small number of pixels available in each cluster for learning the foreground distributions used in the energy functions. Most of the largest pedestrians are detected. The one on the right and the middle is not most of the time, because there are not enough moving pixels with valid optical flow vectors on him. Also, some clusters corresponding to noise are sometimes detected. However, as there is not any constancy in their detection, once again, a temporal consistency or a tracking scheme would probably flag these objects as false detections.

In this subsection we have shown promising results on three different kinds of sequences. To extend the validation of our method, comparisons with other motion detection methods are presented in the sequel.

Comparison with other motion detection methods

We now compare our algorithm to the background modeling method of Grimson and Stauffer [START_REF] Grimson | Using adaptive tracking to classify and monitor activities in a site[END_REF] and the non parametric background modeling of Elgammal et al. [START_REF] Elgammal | Non-parametric model for background subtraction[END_REF]. We also show the masks of moving pixels resulting from the robust estimation and the compensation of the dominant motion with the technique of Odobez and Bouthemy [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]. In order to get motion detection results as binary maps, all the objects segmented by our approach are set as the moving areas of the image (white INRIA pixels). Also note that the sequences used here were taken by a fixed camera so that the background modeling methods can be applied.

The first sequence is a grayscale sequence of a pedestrian walking in front of water (figure 10). This is not a very difficult sequence for our algorithm as there are not many moving pixels detected in the water (second row of figure 10). While our method detects and segments well the person, the bike and the small cars moving on the bridge, the other methods detect several moving pixels in the water while there are many holes on the person masks. The problem of holes in the detection mask (moving pixels) was mentioned in section 2. Here they result from the fact that the coat of the person is not highly textured and the person is walking slowly. Despite this lack of motion information at places, the segmentation step enables to recover the whole person as one unique moving object.

The second sequence on which we compare our algorithm to others is a color sequence of two pedestrians walking behind some waving trees. This sequence is hard because of the complex movements of the trees and the frequent occlusions of the two persons. Grimson and Stauffer's method does not permit to distinguish the trees from the building behind. The results obtained by the non parametric method and the detection of moving pixels method are similar. These two methods do not detect the buildings but find many spurious motions within the trees. Our method also detects few objects within the trees but the segmentation permits to drop most of them. The masks of objects obtained by our algorithm are not perfect though. The partial occlusions lead to many high contours on the pedestrians which stop the flow. Also, due to these occlusions, the clusters are usually very small, which does not permit to construct enough good distributions of the interesting objects. Still the results obtained by our method, as compared to other motion detection techniques, are very promising.

Conclusion

A new technique to detect and segment moving objects in complex dynamic scenes shot by possibly moving cameras has been presented in this paper. The algorithm can be divided into three main steps.

First, a set of points are selected and described in terms of color and motion, then these points are clustered according to their descriptors and finally a segmentation is obtained from these clusters.

Each cluster corresponds to a moving area of the image. Several salient aspects of the contribution can be emphasized. First we only work on "moving" pixels belonging to a grid regularly spread on the whole image and with a motion estimate that passed a statistical test. Second, we use position, color and motion to describe each point. Third, we use a variable bandwidth mean shift using the balloon estimator and an automatic bandwidth selection to create the clusters. And finally, we use sparse motion data in an optimization framework to get the final segmentations of moving objects.

Until now, we have segmented each object independently. We are now trying to segment all the objects jointly by using a multilabel energy function that could be minimized using the α-expansion algorithm [START_REF] Boykov | Markov random fields with efficient approximations[END_REF][START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF].

It is worth emphasizing that the parameters involved in the preliminary motion computations (optic flow and parametric dominant motion) are fixed to the same values in all experiments, while the other parameters (for clustering and segmentation) are automatically selected.

Finally, the proposed method does not make use of any temporal consistency. We are now studying the introduction of such a consistency either on a frame-to-frame basis or within a tracker whose (re)initialization would rely on detection maps. To that end, we proposed in [START_REF] Bugeau | Joint tracking and segmentation of objects using graph cuts[END_REF] a first method to combine the tracking and the segmentation phase.

Figure 1 :

 1 Figure 1: Exemple of a driver sequence

Figure 2 :Figure 3 :

 23 Figure 2: Motion maps for the water skier sequence. (a)-(c) Initial grayscale frames 107 to 109. (d) Displaced frame difference at time t. (e) Displaced frame difference at time t + 1. (f) Result M t of the pixel-wise motion detector.

Figure 4 :

 4 Figure 4: Result of the Harris corner detector on frame 16 of the driver sequence and frame 108 of the water skier sequence.

  = (d x , d y ) associated to each of its point s. The size of the grid G will be denoted as |G|.

Figure 5 :

 5 Figure 5: Result of grid construction (second row) and associated motion fields (third row) on frame 16 of the driver sequence. (a) Original image. (b) Moving pixels. (c) Final grid without flow vectors validation (associated motion field is shown on image (f)). (d) Final grid obtained by a correlation validation with correlation threshold at 0.5. (associated motion field is shown on image (g)). (e) Final grid obtained by using the p-value (associated motion field is shown on image (h)).

Figure 6 :Figure 7 :Figure 8 : 50 INRIA

 67850 Figure 6: Results on the water skier sequence for frames 74, 124, 144, 214, 232, 236 and 242

Figure 9 :Figure 11 :

 911 Figure 9: Results on the traffic sequence. Bounding boxes of detected clusters are shown.

where z

		(G) t (s) is a one
	dimensional vector of grayscale value, z (C) t (s) a 3-dimensional vector of color values and z (M) t	(s) a
	2-dimensional vector characterizing the apparent motion. The computation of z (M) t	(s) is explained
	in section 2.3.1 and the one of z (C) t (s) in section 2.3.2. In case of grayscale sequences, the feature
	vector only contains two elements: z t (s) = (z

  to keep Lucas and Kanade algorithm, which is easy to implement, not expensive, robust to noise, and we verify afterwards if the flow vectors are good or not. To validate values of displacement, a comparison is done between the neighborhood of pixel s = (x, y) in image at time t (data sample X), and the neighborhood of the corresponding point s ′ = (x + dx, y + dy) at time t + 1 (data sample Y ). The linear relationship between intensity values of X and Y is often estimated by computing the normalized cross correlation γ. This coefficient is also known as Pearson product-moment correlation coefficient by statisticians. A value near zero indicates that the two samples are uncorrelated. Unfortunately, the correlation does not take into account the individual distributions of X and Y . Hence it is a poor statistics for deciding whether or not two distributions are really correlated.

In statistics, a result is called significant if it is unlikely to have occurred by chance. Statistical tests exist to assess this correlation. They are based on two principal notions: the null hypotheses and the alpha risk (or fixed level testing). In 1928, Neyman and Pearson

RR n°6282

Bugeau & Pérez

Algorithm 1 Iterative estimation of mean shift bandwidths

Given n data points x (i) = (x (i) ρ )ρ=1...P , i = 1 . . . n. Given a set of Bρ predefined bandwidths

ρ , b = 1 . . . B} for each feature space ρ, the bandwidth selection is as follows. For ρ ′ = 1, . . . , P

• Evaluate the bandwidth at the partition level: For all b = 1, . . . , B

1. For all ρ = 1, . . . , P, ρ = ρ ′ , i = 1, . . . , n compute H(i) ρ :

3. Partition the data using the balloon mean shift partitioning. The result is b) . We introduce the function c that associates the i -th data point to its cluster:

u ) of each cluster using:

and

5. Associate to each point x (i) the mean µ (b) c(i,b) and covariance Σ (b) c(i,b) of the cluster it belongs to. Denote p (i,b) the corresponding normal distribution.

• Evaluate the bandwidth at the data level: For each point x (i) 1. Select the scale b ′ giving the most stable normal distribution by solving:

where JS is the Jensen Shanon divergence defined by:

2. The best bandwidth

Original images

Moving pixels (section 2)

Non parametrical method [START_REF] Elgammal | Non-parametric model for background subtraction[END_REF] Grimson and Stauffer method [START_REF] Grimson | Using adaptive tracking to classify and monitor activities in a site[END_REF] Our method t = 34 t = 59 t = 84


R Boutteau
email: boutteau@esigelec.fr.

X Savatier
email: savatier@esigelec.fr.

J Y Ertaud
email: ertaud@esigelec.fr.

B Mazari
email: mazari@esigelec.fr.

An omnidirectional stereoscopic system for mobile robot navigation

Keywords: omnidirectional sensor, calibration, stereovision

This paper proposes a scheme for a 3D metric reconstruction of the environment of a mobile robot. We first introduce the advantages of a catadioptric stereovision sensor for autonomous navigation and how we have designed it with respect to the Single Viewpoint constraint. For applications such as path generation, the robot needs a metric reconstruction of its environment, so calibration of the sensor is required. After justification of the chosen model, a calibration method to obtain the model parameters and the relative pose of the two catadioptric sensors is presented. Knowledge of all the sensor parameters yields the 3D metric reconstruction of the environment by triangulation. Tools for calibration and relative pose estimation are presented and are available on the author's web page. The entire process has been evaluated using real data.

I. INTRODUCTION

Interest in autonomous robots has been growing over the past few years in many applications: intervention in hostile environments, preparation of military intervention, mapping, etc. In many case, the navigation should be done in an unknown environment and can be helped by 3D reconstruction.

Autonomous navigation requires a large field of view to provide a complete map of the environment. Therefore, interest for omnidirectional vision, which provides a 360-degree field of view, has consequently grown up significantly. To obtain a panoramic image, several methods are being explored: rotating cameras [START_REF] Benosman | Panoramic stereovision sensor[END_REF], multicamera systems and catadioptric sensors [START_REF] Baker | A theory of single-viewpoint catadioptric image formation[END_REF]. We chose to work on catadioptric sensors (camera/mirror combination) because they are able to provide a panoramic image instantaneously and without moving parts.

To achieve a 3D reconstruction, it is necessary to have two images of the environment from two different viewpoints either by moving a single camera or by using a multicamera system. The second way is more suitable for dynamic environments because it provides stereo-pair images instantaneously and the 3D reconstruction can be initialized without any motion.

Our work focuses on the development of a catadioptric stereovision sensor and of a complete toolchain, from calibration to 3D reconstruction.

The following section describes our catadioptric stereovision system and how we designed it to respect the Single View-Point constraint. Section III is devoted to the choice of the model. In section IV, we describe the calibration tool developed to determine the model parameters, while in section V we present the relative pose estimation of the two sensors. In section VI we explain how to obtain the 3D coordinates of points and our experimental results are shown in section VII. Finally, in sections VIII and IX, we draw some conclusions and establish future directions for research.

II. SENSOR DESCRIPTION A. Sensor configuration

It is a well-known fact that a 360-degree field of view presents advantages for navigation since there are interesting optical flow properties [START_REF] Mouaddib | Introduction à la vision panoramique catadioptrique[END_REF] and more visual features to track. Although it is possible to reconstruct the environment with only one camera, a stereoscopic sensor can produce a 3D reconstruction instantaneously (without displacement) and will give better results in dynamic scenes. Among all possible configurations of central catadioptric sensors describes by [START_REF] Baker | A theory of single-viewpoint catadioptric image formation[END_REF], we have chosen to combine two hy-perbolic mirrors with two cameras (see Fig. 1) for the sake of compactness (a parabolic mirror needs a bulky telecentric lens).

B. How to respect the Single-Viewpoint (SVP) constraint

The formation of images with catadioptric sensors is based on the Single-Viewpoint (SVP) theory [START_REF] Baker | A theory of single-viewpoint catadioptric image formation[END_REF]. When the Single-Viewpoint constraint is respected, sensed images are geometrically correct (pure perspective) and the epipolar geometry is applicable. In the case of a hyperbolic mirror (see Fig. 2), the optical center of the camera has to coincide with the second focus F of the hyperbola located at a distance of 2e from the mirror focus. The eccentricity e is a parameter of the mirror given by the manufacturer. Fig. 2. Image formation with a hyperbolic mirror. The camera center has to be located at 2e from the mirror focus to respect the SVP constraint.

A key step in designing a catadioptric sensor is to respect this constraint as much as possible. To achieve this, we first calibrate our camera with a standard calibration tool such as that of J.Y. Bouguet [START_REF] Bouguet | Jean-yves bouguet's www homepage[END_REF] to determine the central point and the focal length. Knowing the parameters of both the mirror and the camera, the image of the mirror on the image plane can be easily predicted if the SVP constraint is respected (see Fig. 2). The expected mirror boundaries are drawn and the mirror has then to be moved to fit this estimation (see Fig. 3). III. PROJECTION MODEL The modelling of the sensor is a necessary step before 3D reconstruction can take place because it establishes the relation between the 3D points of the scene and their projections in the image (pixel coordinates). Although there are many calibration methods, they can be classified into two main categories: parametric and non-parametric. The first family consists in finding an appropriate model for the projection of a 3D point onto the image plane. Non-parametric methods associate one projection ray to each pixel [START_REF] Ramalingram | Towards complete generic camera calibration[END_REF], [START_REF] Ragot | Calibration of a panoramic stereovision sensor : Analytical vs interpolation-based methods[END_REF], and provide a "black box model" of the sensor. They are well adapted for general purposes but they restrict the number of suitable 3D reconstruction algorithms. We consequently chose a parametric calibration method. Using a parametric method requires the choice of the model, which is very important because it has an effect on the complexity and the precision of the calibration process. Several models are available for catadioptric sensors: complete model, polynomial approximation of the projection function and generic model.

The complete model relies on the mirror equation, the camera parameters and the rigid transformation between them to calculate the projection function [START_REF] Gonzalez-Barbosa | Fast dense panoramic stereovision[END_REF]. The large number of parameters to be estimated leads to an error function which is difficult to minimize because of numerous local minima [START_REF] Mei | Single view point omnidirectional camera calibration from planar grids[END_REF]. The polynomial approximation of the projection function was introduced by Scaramuzza [START_REF] Scaramuzza | A flexible technique for accurate omnidirectional camera calibration and structure from motion[END_REF], who proposed a calibration toolbox for his model. The generic model, also known as the unified model, was introduced by Geyer [START_REF] Geyer | A unifying theory for central panoramic systems and practical implications[END_REF] and Barreto [START_REF] Barreto | A unifying geometric representation for central projection systems[END_REF], who proved its validity for all central catadioptric systems. This model was then modified by Mei, who generalized the projection matrix and also took into account the distortions (see [START_REF] Mei | Single view point omnidirectional camera calibration from planar grids[END_REF] for more details). We chose to work with the unified model described by Mei (see Fig. 4) because any catadioptric system can be used and the number of parameters to be estimated is quite reasonable.

A. Projection of a 3D point

A detailed description of the model can be found in [START_REF] Mei | Single view point omnidirectional camera calibration from planar grids[END_REF]. The model contains 11 intrinsic parameters: one mirror parameter (ξ), 5 distortion parameters (k 1 , k 2 , k 3 , k 4 , k 5 ) and 5 parameters for the generalised camera projection matrix (α, γ u , γ v , u 0 , v 0 ). To summarize, the projection p of a 3D point X can be computed using the following steps (see Fig. 4):

• The world point X in the mirror frame is projected onto the unit sphere:

X → X S • This point is then changed to a new reference frame cen- tered in C p • It is then projected onto the normalized plane: X S → m • Distortions are added: m → m d • A perspective projection is then applied: p = Km d =   γ u γ u α u 0 0 γ v v 0 0 0 1   m d (1)

B. Lifting

The lifting step is the calculation of the point X S on the unit sphere corresponding to a pixel, which is very useful for triangulation. Given m = [x y 1] T the coordinates of the pixel on the normalized plane, we have:

X S =      ξ+ √ 1+(1-ξ 2 )(x 2 +y 2 )
x 2 +y 2 +1

x

ξ+ √ 1+(1-ξ 2 )(x 2 +y 2 ) x 2 +y 2 +1 y ξ+ √ 1+(1-ξ 2 )(x 2 +y 2 ) x 2 +y 2 +1 -ξ      (2)
IV. CALIBRATION The calibration step is very easy to achieve because it only requires the catadioptric sensor to observe a planar pattern at different positions. The pattern can be freely moved (the motion does not need to be known) and the user only needs to select the four corners of the pattern.

The calibration process is similar to that of Christopher Mei [START_REF] Mei | Single view point omnidirectional camera calibration from planar grids[END_REF]. It consists of a minimization over all the model parameters of an error function between the estimated projection of the pattern corners and the measured projection thanks to the Levenberg-Marquardt algorithm.

V. RELATIVE POSE ESTIMATION

Once the sensor calibration is completed, we have to know the relative pose of the two sensors to achieve a 3D reconstruction by triangulation. There are two ways to determine the relative pose of the two sensors : either by finding the essential matrix using the epipolar geometry then decomposing it to obtain the translation and the rotation between the two sensors, or by using a pattern to find the relative pose directly. Both methods have been implemented.

A. Essential matrix estimation using epipolar geometry

Epipolar geometry describes the relationship between two cameras. Since it depends only on the parameters of the cameras and their relative pose (it is independent of the scene structure), it can be computed from the correspondence of a few points. The main aim of determining the epipolar geometry is to compute the relative pose of the two sensors (or the displacement of one sensor) but it is also useful for simplifying the search for corresponding points in the two images.

Epipolar geometry is well-known for classical cameras [12] [13]. Given a point on the first image, a line on which the corresponding point lies on the second image can be determined. This line is known as an epipolar line and corresponds to the intersection of an epipolar plane with the image plane. Epipolar geometry of catadioptric sensors is more complicated because of the complex shape of the mirror. It is nevertheless possible to argue from analogy with a classical camera by working with points on the unit sphere (lifted points) rather than image points.

Let the projection of a 3D point X onto the unit spheres be denoted X S1 and X S2 (see Fig. 5), and R and t be the rotation and the translation between the two sensors. The coplanarity constraint of the points X, X S1 , X S2 , C 1 and C 2 can be expressed as follows: The coplanarity constraint (3) can be expressed in matrix form:

X S2 R(t ∧ X S1 ) = 0 (3)
X T S2 EX S1 = 0 (4)
where

E = RS (5)
is the essential matrix first introduced by Longuet-Higgins [START_REF] Longuet | A computer algorithm for reconstructing a scene from two projections[END_REF] and S is an antisymmetric matrix characterizing the transla-

tion: S =   0 -t z t y t z 0 -t x -t y t x 0   .
Given two lifted points X S1 = x 1 y 1 z 1 T and

X S2 = x 2 y 2 z 2 T corresponding to the same 3D point X, (4)
becomes for each pair of matched points: With n matched points (n ≥ 8), we can build a system of equations of the form [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]:

Ae = 0 (7)
where e = e 11 e 12 e 13 • • • e 31 e 32 e 33 T .

The essential matrix is computed through the following steps:

• detection of Harris corners • matching of the corners with backward correlation in order to reject the false matches as described in [START_REF] Werth | A novel bidirectional framework for control and refinement of areabased correlation techniques[END_REF] • resolution of [START_REF] Gonzalez-Barbosa | Fast dense panoramic stereovision[END_REF]

by Singular Value Decomposition of A B. Calibration with pattern

To plan its path, a robot must have a metric model of its environment, which implies knowledge of the baseline. Nevertheless, the essential matrix is estimated up to a scale factor with the eight-point algorithm [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], which leads only to a projective reconstruction. To solve this problem, we developed another method for the relative pose estimation problem.

This method uses patterns shown at different positions as with calibration. Let X be a point with coordinates

X 1 = x 1 y 1 z 1
T in the first frame associated with the first sensor (see Fig. 6). Its coordinates in the second sensor reference frame can be expressed by: With n control points, we can build the following system using (8):

  x 2 y 2 z 2   =   r 11 r 12 r 13 r 21 r 22 r 23 r 31 r 32 r 33     x 1 y 1 z 1   +   t x t y t z   (8)
           x 1 1 y 1 1 z 1 1 0 0 0 0 0 0 1 0 0 0 0 0 x 1 1 y 1 1 z 1 1 0 0 0 0 1 0 0 0 0 0 0 0 x 1 1 y 1 1 z 1 1 0 0 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x n 1 y n 1 z n 1 0 0 0 0 0 0 1 0 0 0 0 0 x n 1 y n 1 z n 1 0 0 0 0 1 0 0 0 0 0 0 0 x n 1 y n 1 z n 1 0 0 1                               r 11 r 12 r 13 r 21 r 22 r 23 r 31 r 32 r 33 t x t y t z                    =            x 1 2 y 1 2 z 1 2 . . . x n 2 y n 2 z n 2            (9)
The 3D coordinates of the pattern corners are estimated using the Levenberg-Marquardt algorithm. The resolution of [START_REF] Scaramuzza | A flexible technique for accurate omnidirectional camera calibration and structure from motion[END_REF] gives the relative pose (R, t) of the two sensors including the scale factor.

Once the relative pose is computed, the essential matrix can be easily computed thanks to [START_REF] Ramalingram | Towards complete generic camera calibration[END_REF].

VI. 3D RECONSTRUCTION

Reconstruction consists in determining the 3D coordinates of a point X, given its projections x 1 and x 2 onto the two images. Consequently, the first step is to match pixels.

A. Epipolar geometry helps pixel matching

In addition to the essential matrix estimation, epipolar geometry possesses some helpful properties for pixel matching.

For each point X S1 on the first unit sphere, the epipolar curve is characterized by the normal of the plane n 2 = EX S1 . The epipolar curve on the second image is obtained by determining the intersection of the plane and the second equivalence sphere, which is a great circle. In the same way, n 1 = E T X S2 characterizes the epipolar curve corresponding to X S2 on the second sphere.

The epipoles are the points of intersection of the line joining the center with the unit sphere. All the epipolar curves intersect at the two epipoles. They can be directly computed from the essential matrix by singular value decomposition as in the classical case [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF].

B. Mid-point method

In theory, the rays corresponding to two matched pixels must intersect at X. In practice, however, various types of noise (distortions, small errors in the model parameters, etc) lead to lines generated by corresponding image points which do not always intersect. The problem is to find a 3D point which optimally fits the measured image points [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF].

This step is achieved by using the mid-point method (see Fig. 7). Lifted points X S1 and X S2 corresponding to the pixels are first computed using [START_REF] Baker | A theory of single-viewpoint catadioptric image formation[END_REF]. These points are then used to define the direction vectors v 1 and v 2 of the two rays. The 3D point X is chosen as the mid-point of the shortest transversal between the two rays. VII. EXPERIMENTAL RESULTS Each step of the 3D reconstruction was evaluated and a global result is shown at the end of this section. Experiments were implemented on real images and without prior knowledge to evaluate our system in real conditions.

A. Calibration

The calibration was evaluated by computing the Root of Mean Squares (RMS) distances (in pixels) between the extracted corners (77 per image) of the pattern and the reprojected ones. Figure 8 shows the RMS error with respect to the number of images of the pattern used for the calibration. We can notice that the error decreases when more images are used but tends to stabilize from 6 images. The minimum RMS value obtained is then about 0.02 pixels. The model parameters obtained by the calibration allows the lifting, i.e. the projection of the pixels onto the unit sphere as shown in Figure 9.

B. Relative pose estimation

The relative pose estimation was evaluated on real images. The sensor was mounted on a graduated rail and was moved 10cm by 10cm. At each position, an omnidirectional image was acquired with the aim of computing the displacement of the sensor according to the first position using five calibration patterns placed in the room (see Fig. 10).

Table I summarizes the results. The average error is less than 0.9%.

TABLE I ESTIMATION OF THE DISPLACEMENTS

Once the relative pose is estimated, we can compute the essential matrix and use it to check epipolar geometry properties. In figure 10, for each selected pixel on the left image (red crosses), the corresponding epipolar curve (green curves) is drawn on the right image and vice versa.

C. 3D reconstruction

The entire process was evaluated by a piecewise planar 3D reconstruction of our laboratory (see Fig. 11). The four corners of the ceiling were manually selected to check the size of the room. After triangulation, we estimated the ceiling dimension at 6.39 x 6.34 meters and the actual size is 6.4 x 6.4 m.

Fig. 1 .

 1 Fig. 1. View of our catadioptric stereovision sensor mounted on a Pioneer robot. Baseline is approximately 20cm for indoor environments and can be extended for outdoor environments. The overall height of the sensor is 40cm.

Fig. 3 .

 3 Fig. 3. Adjustment of the mirror position to respect the SVP constraint. The mirror border has to fit the estimation (green circle).

Fig. 4 .

 4 Fig. 4. Unified projection model.

Fig. 5 .

 5 Fig. 5. Epipolar geometry of spherical sensors.

x 2 x 1 e

 1 11 + x 2 y 1 e 12 + x 2 z 1 e 13 + • • • + z 2 z 1 e 33 = 0 (6) where E =   e 11 e 12 e 13 e 21 e 22 e 23 e 31 e 32 e 33   .

Fig. 6 .

 6 Fig. 6. Relative pose estimation principle.

Fig. 7 .

 7 Fig. 7. Mid-point method. The 3D point X is chosen as the mid-point of the shortest transversal between the two rays.

Fig. 8 .

 8 Fig. 8. RMS error versus the number of images of the pattern.

Fig. 9 .

 9 Fig. 9. An omnidirectional image and two views of its projection onto the unit sphere.

Fig. 10 .

 10 Fig. 10. Epipolar curves (green) corresponding to selected pixels (red crosses).

Fig. 11 .

 11 Fig. 11. Piecewise planar 3D reconstruction of our laboratory.

ACKNOWLEDGMENT

This research is supported by the DGA (Délégation Générale pour l'Armement), IPSIS (Ingénierie pour Signaux et Systèmes) and IRSEEM (Institut de Recherche en Systèmes Electroniques Embarqués).

VIII. DISCUSSION

In this paper, we have used a stereovision system based on two catadioptric sensors (see Fig. 1) mounted co-axially one above the other. Indeed, such a configuration greatly simplifies the epipolar geometry since epipolar curves become radial lines. By this way, matching horizontal lines is simplified while a motion allows vertical lines to be matched. A stereoscopic system also presents advantages compared to a structure from motion method because the baseline is well-known and not estimated by odometry which introduces many errors.

We have chosen to work with a unified model which is valid for all central sensors hence we have to respect the Single-Viewpoint constraint. Our method, presented in section 2, allows us to be as near as possible to the single viewpoint. The choice of the Mei's unified model leads to a very flexible calibration step because the number of parameters to be estimated is quite reasonable and easy to be initialized. A calibration tool was developed in C++ using the computer vision library OpenCV and can be freely downloaded from our website [START_REF] Boutteau | 3d reconstruction for autonomous navigation[END_REF]. The calibration tool does not require any commercial software and optimizes the computing time. A calibration with 10 images does not exceed 2 minutes.

Two methods have been developed for the relative pose estimation and have different fields of use. The first one relies on the eight-point algorithm and can be used in line because the process is entirely automatic but gives only the relative pose up to a scale factor. The other one uses calibration patterns and is more adapted for measurement thanks to the entire knowledge of the relative pose. Table II summarizes the main characteristics of the two methods used for the relative pose estimation and their field of use. The evaluation of our method using calibration patterns was done by estimating the motion of the sensor. The average error between estimation and real displacements is small (less than 2%) and probably induced by image noise which disturbs the corner extraction process. The precision is sufficient for applications such as navigation.

The piecewise planar 3D reconstruction is introduced for the validation of the calibration and triangulation steps. The good results obtained imply the reliable estimation of the parameters. The 3D reconstruction was done on manually se-lected points to avoid possible errors due to an automatic matching.

IX. CONCLUSION

In this paper, we have presented a scheme for 3D reconstruction of the environment of a mobile robot based on a catadioptric stereovision sensor. This scheme can be decomposed into three stages : calibration, estimation of the relative pose of the two sensors and triangulation.

Thanks to an accurate calibration and relative pose estimation, the metric of the scene is well respected and evaluated using real images.

Our future work will focus on automatic reconstruction, thanks to pixel matching or more elaborated primitives (lines, planes, etc). Once 3D reconstruction is achieved, the motion estimation of the robot will be addressed to be able to merge several local reconstructions.