Rainer Kümmerle

Rudolph Triebel
email: rudolph.triebel@mavt.ethz.ch

Patrick Pfaff

Wolfram Burgard

Monte Carlo Localization in Outdoor Terrains using Multi-Level Surface Maps

In this paper we consider the problem of mobile robot localization with range sensors in outdoor environments. Our approach applies a particle filter to estimate the full six-dimensional state of the robot. To represent the environment we utilize multi-level surface maps which allow the robot to represent vertical structures and multiple levels in the environment. We describe probabilistic motion and sensor models to calculate the proposal distribution and to evaluate the likelihood of observations. Experimental results obtained with a mobile robot in an outdoor environment indicate that our approach can be used to robustly and accurately localize an outdoor vehicle. The experiments also demonstrate that multi-level surface maps lead to a significantly better localization performance than standard elevation maps.

Introduction

The problem of mobile robot localization with range sensors in outdoor environments arises whenever GPS signals are missing because of occlusions caused by buildings, bridges, or trees. In such situations, a mobile robot typically has to estimate its position in the environment using its other sensors and a map of the environment. In this paper, we consider the problem of localizing a mobile robot in outdoor environments by matching laser range measurements to a given map of the environment. One of the most popular representations for outdoor environments are elevation maps [START_REF] Bares | Ambler: An autonomous rover for planetary exploration[END_REF][START_REF] Hebert | Terrain mapping for a roving planetary explorer[END_REF][START_REF] Lacroix | Autonomous rover navigation on unknown terrains: Functions and integration[END_REF][START_REF] Parra | 3-d modelling and robot localization from visual and range data in natural scenes[END_REF]. The key idea underlying elevation maps is to store the 2 1 2dimensional height information of the terrain in a two-dimensional grid, which corresponds to a representation of the horizontal surfaces of the environment. Whereas the knowledge about the horizontal surfaces is well suited to support traversability analysis and path planning, it provides only weak support for the localization of the vehicle. Modeling only the horizontal surfaces, namely, means that vertical structures, which are frequently perceived by ground based vehicles, cannot be used to support localization. To avoid this problem, multi-level surface (MLS) maps [START_REF] Triebel | Multi level surface maps for outdoor terrain mapping and loop closing[END_REF] have been introduced. They can be regarded as an extension of the classical elevation maps as they additionally represent intervals corresponding to vertical objects in the environment. A further disadvantage of elevation maps is that they cannot represent multiple levels. This, for example, is important when mobile robots are deployed in environments with bridges or underpasses.

Fig. 1 depicts examples of an elevation map (left) and the corresponding MLS map (right) of the campus at the University of Freiburg. As can be seen from the images, the MLS map is able to represent the environment more accurately than the elevation map. In the MLS map, objects like trees and walls are represented properly. Another example is shown in Fig. 2. Here the robot is located in front of a wall. In the elevation map, the wall is not represented correctly, because the height values obtained from beams reflected by the wall are averaged which is the typical approach in elevation maps. This can lead to a poor estimate of the measurement likelihood at the particular robot position. In contrast, when the MLS map is used, one obtains a better value of the likelihood, because the wall is modeled correctly.

In this paper we present an approach to use the multi-level surface maps for localization. We present probabilistic motion and observation models and describe how these models can be utilized in a probabilistic localization scheme. We furthermore evaluate how the localization performance changes when standard elevation maps are used instead of MLS maps.

This paper is organized as follows. After discussing related work in the next section, we briefly describe Monte Carlo localization in Section 3. Our approach to apply the MLS maps is presented in Sections 4 and 5. Finally, in Section 6, we present experimental results illustrating the advantages of applying MLS maps for localization in outdoor environments.

Related Work

The problem of localization of mobile robots in outdoor environments with range sensors or cameras has been studied intensively in the past. For example, Adams et al. [START_REF] Adams | Particle filter based outdoor robot localization using natural features extracted from laser scanners[END_REF] extract predefined features from range scanners and apply a particle filter for localization. Davison and Kita [START_REF] Davison | 3d simultaneous localisation and map-building using active vision for a robot moving on undulating terrain[END_REF] utilize a Kalman filter for vision-based localization with point features on non-flat surfaces. Recently, Agrawal and Konolige [START_REF] Agrawal | Real-time localization in outdoor environments using stereo vision and inexpensive gps[END_REF] presented an approach to robot localization in outdoor terrains based on feature points that are tracked across frames in stereo images. Lingemann et al. [START_REF] Lingemann | High-speed laser localization for mobile robots[END_REF] recently described a method for fast localization in in-and outdoor environments. Their system operates on raw data sets, which results in huge memory requirements. Additionally, they apply a scan-matching routine for localization, which does not facilitate global localization. To reduce the memory requirements of outdoor terrain representations, several researchers applied elevation maps [START_REF] Bares | Ambler: An autonomous rover for planetary exploration[END_REF][START_REF] Hebert | Terrain mapping for a roving planetary explorer[END_REF][START_REF] Lacroix | Autonomous rover navigation on unknown terrains: Functions and integration[END_REF][START_REF] Parra | 3-d modelling and robot localization from visual and range data in natural scenes[END_REF]. A probabilistic approach to localize a planetary rover in such elevation maps has been described by Olson [START_REF] Olson | Probabilistic self-localization for mobile robots[END_REF]. In this system, elevation maps were sufficient to robustly localize the vehicle, mainly because the number of vertical and overhanging objects is negligible in environments like on Mars. However, environments on earth contain many objects like buildings or trees which have vertical or even overhanging surfaces. To address this issue, Pfaff et al. [START_REF] Pfaff | An efficient extension to elevation maps for outdoor terrain mapping and loop closing[END_REF] extended the elevation map approach by a cell classification. In this approach the individual cells are divided into three classes: cells which have been observed from above, cells with vertical objects, and cells with overhanging objects. This extension has been utilized to improve the data association during the scan matching process and enabled the robot to traverse cells with overhanging objects. In contrast to MLS maps, however, this approach still lacks the ability to store multiple surfaces at one position, which prevents the robot from dealing with situations, in which it has to traverse a bridge and move through the corresponding underpass. The goal of this paper is to develop a probabilistic localization method based on MLS maps and to demonstrate that the more accurate representation of the environment results in improved localization capabilities.

Monte Carlo Localization

To estimate the pose x = (x, y, z, ϕ, ϑ, ψ) of the robot in its environment, we consider probabilistic localization, which follows the recursive Bayesian filtering scheme. The key idea of this approach is to maintain a probability density p(x t | z 1:t , u 0:t-1 ) of the robot's location x t at time t given all observations z 1:t up to time t and all control inputs u 0:t-1 up to time t -1. This posterior is updated as follows:

p(x t | z 1:t , u 0:t-1 ) = α • p(z t | x t ) • p(x t | u t-1 , x t-1 ) • p(x t-1 ) dx t-1 . (1)
Here, α is a normalization constant ensuring that p(x t | z 1:t , u 0:t-1 ) sums up to one over all x t . The terms to be described in Eqn. ( 1) are the prediction model p(x t | u t-1 , x t-1 ) and the sensor model p(z t | x t ). One major contribution of this paper is an appropriate computation of these models in the case that an MLS map is given.

For the implementation of the described filtering scheme, we use a sample-based approach which is commonly known as Monte Carlo localization [START_REF] Dellaert | Monte Carlo localization for mobile robots[END_REF]. Monte-Carlo localization is a variant of particle filtering [START_REF] Doucet | Sequential Monte-Carlo Methods in Practice[END_REF] where each particle corresponds to a possible robot pose and has an assigned weight w i . The belief update from Eqn. ( 1) is performed by the following two alternating steps:

1. In the prediction step, we draw for each particle with weight w i a new particle according to w i and to the prediction model p(x t | u t-1 , x t-1 ). 2. In the correction step, a new observation z t is integrated. This is done by assigning a new weight w i to each particle according to the sensor model p(z t | x t ).

Furthermore, the particle set needs to be re-sampled according to the assigned weights to obtain a good approximation of the pose distribution with a finite number of particles. However, the resampling step can remove good samples from the filter which can lead to particle impoverishment. Thus, we need to find a criterion for deciding when to perform the resampling step. In this paper, we calculate the number N eff of effective particles according to the formula proposed by Doucet et al. [START_REF] Doucet | Sequential Monte-Carlo Methods in Practice[END_REF]

N eff = 1 N i=1 wi 2 , (2)
where wi refers to the normalized weight of sample i. In particular, we only resample if N eff drops below the threshold of N 2 where N is the number of samples. In the past, this approach has already successfully been applied in the context of the simultaneous mapping and localization (SLAM) problem [START_REF] Grisetti | Improving grid-based slam with raoblackwellized particle filters by adaptive proposals and selective resampling[END_REF].

Prediction Model for MLS Maps

The prediction model p(x t | u t-1 , x t-1 ) we use is based on an approach introduced by Eliazar et al. [START_REF] Eliazar | Learning probabilistic motion models for mobile robots[END_REF]. It reflects systematic errors such as drift, as well as the uncertainty in the execution of an action u = (x u , y u , θ u ), where (x u , y u ) is the translation and θ u the rotation angle. To incorporate this 2D motion into our 3D map we proceed as follows. First, we obtain a possible outcome (x v , y v , θ v ) of the action by applying the probabilistic model. Then, we adapt the motion vector v = (x v , y v ) to the shape of the 3D surface traversed by the robot. This surface is obtained from the given MLS map and consists of planar square patches. To adapt the motion vector, we discretize it into segments of length c, which is the cell size of the MLS map, in our case 0.1 m. For each segment, we determine the corresponding surface patch S and rotate the segment according to the orientation (ϕ S , ϑ S ) of the patch, where ϕ S is the rotation about the x-axis and ϑ S the rotation about the y-axis. The patch orientation is computed from the normal vector n S of the patch S, which in turn is obtained by fitting a plane into the local vicinity of S. The normal vector computation is done beforehand and constitutes an extension to the framework of MLS maps. In general, it is not robust against noise and small errors in the MLS map, which results in an uncertainty of the patch orientation. In our approach, we model this uncertainty by adding Gaussian noise to the orientation parameters ϕ S and ϑ S . Thus, our prediction model expresses the uncertainty in 5 out of 6 position parameters -x, y and ψ by the 2D motion model and ϕ and ϑ by our 3D extension. For the last one -the height value z -we have the constraint that the robot must stay on the ground. Therefore, we adjust the z-value manually whenever it is too high or too low. This is illustrated in Fig. 3. Finally, after concatenating all transformed motion vector segments, we obtain a new 3D motion vector v which is added to the current estimate of the robot position x t-1 to obtain a new position estimate x t .

Endpoint Sensor Model for MLS Maps

In our sensor model, we treat each beam independently and determine the likelihood of a whole laser scan by factorizing over all beams. Thus, we have

p(z | x) = K k=1 p(z k | x) (3)
where K is the number of beams in each laser measurement z. In Eqn.

(3) and in the following, we drop the index t for convenience. Our sensor model p(z k | x) is based on an approach that has been introduced by Thrun [START_REF] Thrun | A probabilistic online mapping algorithm for teams of mobile robots[END_REF] as likelihood fields (LF) or end point model. In particular, we formulate the sensor model p(z k | x) for each particular beam as a mixture of three different distributions:

p(z k | x) = α hit p hit (z k | x) + α rand p rand (z k | x) + α max p max (z k | x), (4)
where p hit is a normal distribution N(0, σ 2 ) that models situations in which the sensor detects an obstacle. Random measurements are modeled using a uniform distribution p rand (z k | x). Maximum range measurements are covered by a point mass distribution p max (z k | x). These three distributions are weighted by the non-negative parameters α hit , α rand , and α max , which sum up to one. The values for α hit , α rand , α max , and σ 2 used in our current implementation have been determined empirically.

In the end point model, the probability p hit (z k | x) only depends on the distance d k between the end point of the k-th laser beam and the closest obstacle in the map. Thus, the physical property of the laser beam is ignored, because the model just uses the end point and does not consider the beam characteristic of the laser. Therefore, we need to calculate the global coordinates for a beam end point. If we denote the angle of the k-th beam relative to the zero angle with ζ k , then the end point pk = ( xk , ỹk , zk ) T of that beam in the robot's own coordinate frame is calculated as

          xk ỹk zk           =           x ŷ ẑ          + Rz k           cos(ζ k ) sin(ζ k ) 0           , (5)
where ( x, ŷ, ẑ) T denotes the position of the sensor at time t and R is a rotation matrix that expresses the 3D sensor orientation in the robot's coordinate frame. For a given robot pose x = (x, y, z, ϕ, ϑ, ψ) at time t we can compute the global coordinates

p k = (x k , y k , z k ) T of the k-th beam end point p k as follows           x k y k z k           = R(ϕ, ϑ, ψ)           xk ỹk zk           +           x y z           , (6)
where R(ϕ, ϑ, ψ) denotes the rotation matrix for the given Euler angles ϕ, ϑ, and ψ.

In MLS maps, obstacles are represented as vertical surface patches, which can be seen as vertical segments of occupied space. Unfortunately, there is no efficient way to find the closest of all vertical segments to a given beam end point. Therefore, we use an approximation by uniformly sampling a set P of 3D points from all vertical patches. The distance d k of the k-th beam end point p k to the closest obstacle is then approximated as the Euclidean distance d(p k , P) between p k and P. This distance can be efficiently calculated by storing all points from P in a kD-tree. Equations. ( 5) and ( 6) describe a 3D transform T (z k ; x) of the measurement z k at position x. Using this and the fact that p hit is Gaussian, we can compute p hit as

p hit (z k | x) ≈ 1 √ 2πσ 2 exp        - 1 2 d(p k , P) σ 2        , (7)
where p k = T (z k ; x). Plugging this into Eqn. ( 4) and the result into Eqn. (3), we obtain the entire sensor model.

Experimental Results

The sensor and prediction models have been implemented in a particle filter algorithm and evaluated on real data acquired with a mobile robot. The robot is a Pioneer II AT system equipped with a SICK LMS laser range scanner and an AMTEC wrist unit, which is used as a pan/tilt device for the laser (see Figure 4). During the experiments described in this section, the laser was directed horizontally. The experiments are designed to investigate if the MLS map approach facilitates mobile robot localization and whether it yields better localization performance than the elevation maps.

Global Localization

The first set of experiments is designed to evaluate the performance of the MLS map approach in the context of a global localization task. Figure 5 depicts the convergence of the particles to the true position of the robot with 500,000 and 1,000,000 particles. Whereas the x-axis corresponds to the resampling step, the y-axis shows the number of particles in percent that are closer than 1 m to the true position, which has been

MLS Map Elevation Map

Fig. 6. The left image depicts the number of successful localizations after 15 resampling steps for the two different map representations for particle numbers from 250,000 up to 1,000,000. The right image shows the average localization error over all particles for a tracking experiment with 1,000 particles. In average the use of the MLS maps leads to smaller errors. computed by a tracking experiment with 100,000 particles. Shown are the evolutions of these numbers when the MCL is applied on standard elevation maps and on MLS maps. Note that the elevation map does not reach 100%. This is due to the fact that the sensor model for the standard elevation map relies on a highly smoothed likelihood function, which is good for global localization but does not achieve maximal accuracy during tracking. The application of a more peaked sensor model in the case of the standard elevation map would lead to much higher divergence rates. In both cases, a t-test showed that it is significantly better to apply the MLS maps than the standard elevation maps for the global localization task. Experiments with 250,000 and 750,000 particles showed the same behavior. Figure 6 shows the number of successful localizations for the two different map representations and for different numbers of particles. Here, we assumed that the localization was achieved when every particle differed by at most 1 m from the true location of the robot. We can see that the global localization performs more robust on the MLS map than on the standard elevation map.

Tracking

We also carried out experiments, in which we analyzed the accuracy of the MLS map approach in the context of a position tracking task. To obtain the corresponding data set, we steered along a loop in our campus environment. The traversed trajectory has a length of 284 meters. Figure 7 depicts a top view of the MLS map of our test environment. The blue / dark grey line shows the localized robot poses. The yellow / light grey line shows the pure odometry. Figure 6 depicts the average localization error for a tracking experiment with 1,000 particles. As can be seen from the figure, the MLS map approach outperforms the standard elevation map approach. The tracking experiments have been computed online on a standard PC with an AMD Athlon 64 3200+ processor. In the practical experiments we found that the use of the MLS maps results in a computational overhead of no more than 10% compared to elevation maps.

Conclusions

In this paper, we presented an approach to Monte-Carlo localization (MCL) using multi-level surface (MLS) maps. We applied the main concepts of MCL to this new 3D map representation by adapting the prediction model and the sensor model to MLS maps. Furthermore, we showed in experiments that MLS maps are better suited for the task of local and global localization than standard elevation maps. The slightly increased runtime due to the higher accuracy of MLS maps is compensated with a significantly faster convergence of the particle filter. This makes MLS maps useful for outdoor localization in cases in which no GPS signal is available.

Fig. 1 .

 1 Fig. 1. Elevation Map (left) and multi-level surface (MLS) map (right) of the Freiburg campus. The MLS map represents vertical structures more accurately and can deal with multiple surfaces that can be traversed by the robot.

Fig. 2 .

 2 Fig. 2. Advantage of the MLS map approach in comparison to the standard elevation maps. In contrast to the MLS map (right) the elevation map (left) lacks the ability to model vertical structures, because it averages over all measured height values. Since the distance of the endpoint of a laser beam to the closest point in the elevation map can have substantial deviations from the true distance, localization becomes harder.

Fig. 3 .

 3 Fig. 3. Application of our prediction model to a of 2D motion vectors (black). They are rotated to estimate the 3D motion vectors (red). The dashed line indicates the tolerance interval for the z-coordinate.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Robot Herbert used for the experiments.

Fig. 7 .

 7 Fig. 7. MLS map used for the localization experiments. The area represented by this map spans approximately 195 by 146 meters. The blue / dark grey line shows the localized robot poses. The yellow / light grey line shows the pure odometry. The traversed trajectory has a length of 284 meters.

Acknowledgment

This work has partly been supported by the German Research Foundation (DFG) within the Research Training Group 1103 and under contract number SFB/TR-8 and by the EC under contract number FP6-IST-027140-BACS and FP6-IST-034120, Action line: Micro/Nano Based Sub-Systems.


domain: math info

text:
Adaptive Space-time Examplar-Based Method for Image Sequence Restoration

We present a novel space-time exemplar-based method for image sequence restoration. We have designed a locally adaptive estimation framework which relies on the analysis of the bias-variance trade-off. It allows us to select on-line the best adapted filtering space-time neighborhood at each pixel. Furthermore, it involves an original space-time exemplar-based restoration which improves the overall performance. This method is unsupervised and requires no motion estimation. It nevertheless can be combined with motion estimation to cope with very large displacements due to camera motion. Experiments show that this method is able to drastically improve the quality of highly corrupted image sequences. Quantitative results on standard image sequences demonstrate that our approach outperforms other recent methods.

Introduction

Image sequence restoration takes a crucial place in several important domains. In biology, confocal microscopes provide noisy image sequences due to the limited number of photons. In astronomy, the same phenomenon occurs in telescopes and image sequences are also corrupted with an important noise level. In medicine, image sequence denoising is critical especially in ultra-sound and X-ray imaging. Video surveillance also provides some noisy image sequences. Finally, we can notice that the noise increases the required bandwidth for video sequence broadcasting and that a denoising algorithm can then produce a more compact image sequence digital representation [START_REF] Dugad | Noise reduction in video by joint spatial and temporal processing[END_REF].

The main difficulties in image sequence restoration arises from the nonstationarity of the spatio-temporal signal. Denoising methods must preserve these space-time discontinuities while minimizing the error between the unknown noise-free image sequence and the denoised sequence. The later is usually stated as the minimization of a local quadratic objective function.

In order to recover the correct intensity at every pixel of the image sequence, we can exploit the local correlation of the signal and use the neighboring pixels. A method based on a pre-determined shape and size (scale) will not be able to take into account the spatio-temporal discontinuities of the intensity function. This would result in a blurred image sequence. As for video sequences, temporal discontinuities are mainly due to motion, image sequence denoising methods have to be able to filter along motion trajectories. Thus, the local neighborhood can be oriented according to the motion direction; or conversely, the motion can be compensated for. This implies a motion estimation stage. However, motion estimation is usually not that easy specially in noisy sequences.

Another way to remove noise is to find similar samples in the image sequence. Exemplar-based approach has proved its efficiency in texture synthesis [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF], image in-painting [START_REF] Criminisi | Region filling and object removal by exemplar-based inpainting[END_REF] and also video completion [START_REF] Wexler | Space-time video completion[END_REF]. In [START_REF] Buades | Denoising image sequences does not require motion estimation[END_REF], an exemplar-based image sequence filtering method is also proposed: the non-local mean algorithm. In these methods the similarity between two points is measured by using patches that describe their local context. Examples are usually searched in the whole image or in an area around the considered pixel.

The original method we propose is an exemplar-based adaptive statistical method for image sequence restoration. It is related to the framework already described for still images [START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF][START_REF] Katkovnik | Adaptive window size image denoising based on intersection of confidence intervals (ICI) rule[END_REF][START_REF] Kervrann | An adaptive window approach for image smoothing and structures preserving[END_REF] and image sequences [START_REF] Ercole | Spatio-temporal pointwise adaptive denoising in video: 3d non parametric approach[END_REF]. Methods designed for still image like wavelet shrinkage [START_REF] Rajpoot | Adaptive wavelet restoration of noisy video sequences[END_REF], Wiener filtering [START_REF] Dekeyser | Spatio-temporal Wiener filtering of image sequences using a parametric motion model[END_REF] or PDE based filters [START_REF] Lee | Spatio-temporal video filtering algorithm based on 3-D anisotropic diffusion equation[END_REF] have been extended to image sequences. Unlike robust anisotropic diffusion [START_REF] Black | Robust anisotropic diffusion[END_REF] and non-linear Gaussian filtering [START_REF] Barash | A fundamental relationship between bilateral filtering, adaptive smoothing, and the nonlinear diffusion equation[END_REF], our local adaptive estimation approach supplies scale selection for each pixel by locally estimating the appropriate spacetime filtering window. Moreover, if required, the proposed method is able to process in a decoupled way the space and the time dimensions which is not the case for methods that simply consider the image sequence as a 3D volume [START_REF] Lee | Spatio-temporal video filtering algorithm based on 3-D anisotropic diffusion equation[END_REF][START_REF] Rajpoot | Adaptive wavelet restoration of noisy video sequences[END_REF]. Unlike [START_REF] Katkovnik | Adaptive window size image denoising based on intersection of confidence intervals (ICI) rule[END_REF][START_REF] Ercole | Spatio-temporal pointwise adaptive denoising in video: 3d non parametric approach[END_REF] our approach is not based on a geometrical partition of the neighborhood in sectors but use an appropriate and more flexible weighted sum of data points in a neighborhood. distance between a patch centered on a given pixel and examples patches in the estimated space-time neighborhood. Additionally, a confidence level (i.e. variance) attached to each restored pixel is provided. Lets us finally point out that the Total Variation minimization method [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] cannot be easily extended to space-time domain.

The remainder of the paper is organized as follows. Section 2 describes the general framework of our method and the intensity estimation is defined. In Section 3, the adaptive choice of the local space-time neighborhood is introduced. Section 4 deals with the similarity measure involved in the selection of the patches in the space-time neighborhood of the pixel to be restored and Section 5 gives details of the algorithm implementation. In Section 6, we report an important set of experimental results. Intensive comparisons with other recent methods has been carried out, demonstrating that our method outperforms them. We also present how our denoising method can be combined with a motion estimation method and we comment the usefulness of a motion compensation stage. Finally, Section 7 contains concluding remarks.

Proposed exemplar-based space-time approach

We consider the following statistical image model:

Y i = u(x i ) + ξ i , (1)
where x i ∈ Ω denotes the pixel location in the space-time volume Ω ⊂ IR 3 . The image function u i = u(x i ) has to be recovered from observations Y i . The errors ξ i are assumed to be independent zero-mean Gaussian variables with unknown variance τ 2 . We need minimal assumptions on the structure of the image for recovering u. We assume that u(x i ) is a locally piecewise constant function in a spacetime neighborhood of pixel x i . However, the size and shape of the neighborhood should varie in the image sequence because of non-stationarities and presence of discontinuities. If we can determine the adequate neighborhood for each pixel, then the regression function u can be estimated by optimizing a local maximum likelihood (ML) criterion. The proposed method addresses these two issues as described below.

We construct a sequence of increasing nested space-time neighborhoods (W i,n ) n∈[0:N ] centered at each point x i , i.e. W i,n ⊂ W i,n+1 with N indicating the largest window. Further details about the neighborhood design are described in Section 3. At the initialization, we choose a small neighborhood (the 8 nearest neighbors in the 2D space domain) as the pilot (starting) window W i,0 at point x i , ∀x i ∈ Ω. Then, we compute an initial estimate ûi,0 of u(x i ) and its associated variance σ2

i,0 as:

ûi,0 = xj ∈Wi,0 ω ij Y j and σ2 i,0 = τ 2 xj ∈Wi,0 ω ij (2)
Comment citer ce document : Boulanger, J., Kervrann Patch ûi,n-1 of size 3 × 3 associated to x i (in green).

x j location and corresponding 3 × 3 patch.

Space-time neighborhood W i,n (in blue).

Current pixel located at x i . where τ 2 is an empirical estimate of the noise variance τ 2 as described in Section 5. At the initialization step, the weights ω ij are defined as a function of the distance between two 2D image patches of p 2 pixels centered in x i and x j . At the next iteration, we consider a larger space-time neighborhood W i,1 such that W i,0 ⊂ W i,1 and calculate new estimates ûi,1 and σ2 i,1 over W i,1 . We continue this way, and at iteration n, we define the estimator as:

ûi,n = xj ∈Wi,n ω ij Y j and σ2 i,n = τ 2 xj ∈Wi,n ω 2 ij , (3)
The estimate ûi,n corresponds to a weighted average of the data points located in the space-time neighborhood. We propose to define weights ω ij as a function of the distance between two 2D regularized image patches ûi,n-1 and ûj,n-1 defined at iteration n -1 and centered in x i and x j as shown in Figure 1. In the two next sections, we address the problems of the definition of the neighborhood W i,n and propose the design of a distance to compare image patches.

3 Space-time neighborhood adaptation space-time volume as neighborhood shape. Nevertheless, we separate the space dimensions from the time dimension. Then, the space-time neighborhood is characterized by its extent in the space domain and its extent in the time domain. The use of two distinct extents in space and in time allows us to respect spacetime discontinuities and the space-time domain is not considered as an isotropic 3D volume. In Figure 2(a) the increasing neighborhood sequence is illustrated. At each iteration, the spatial extent and the temporal extent are alternatively increased until a stopping rule is satisfied.

Space-time neighborhood selection

A point-wise rules is used to guide the space-time neighborhood estimation process. This rule aims at estimating the optimal neighborhood at x i and is based on the measure of the closeness of the estimator ûi to the unknown function u i given by the local L 2 risk. This measure of performance can be decomposed in two terms, that is the squared bias and the variance as

E [û i,n -u i ] 2 = [bias (û i,n )] 2 + σ2 i,n , (4)
where E(.) denotes the mathematical expectation. In the sequel, we can reasonably assume that the squared bias is an increasing function of the neighborhood size and the variance is a decreasing function of the neighborhood size [START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF][START_REF] Katkovnik | Adaptive window size image denoising based on intersection of confidence intervals (ICI) rule[END_REF][START_REF] Kervrann | An adaptive window approach for image smoothing and structures preserving[END_REF][START_REF] Lepski | Asymptotically minimax adaptive estimation 1: upper bounds[END_REF]. Then, the optimal neighborhood will be the one that achieves an optimal compromise between these two terms. A closed-form optimal solution for the ideal window is not available for such a non-linear estimator but we can assume that the optimal neighborhood is the one for which the squared bias and the variance are nearly the same: [START_REF] Lepski | Asymptotically minimax adaptive estimation 1: upper bounds[END_REF].

E u * i,n -u i 2 ≈ 2σ 2 * i,n
Comment citer ce document : Boulanger, J., Kervrann A practical rule based on pairwise comparison of successive estimates can be derived to select the optimal neighborhood. It amounts to define the largest neighborhood satisfying the following point-wise statistical rule [START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF][START_REF] Kervrann | An adaptive window approach for image smoothing and structures preserving[END_REF][START_REF] Lepski | Asymptotically minimax adaptive estimation 1: upper bounds[END_REF]

|û i,n -ûi,n | < η σi,n , ∀n < n, (5)
as the optimal neighborhood. This rule can be interpreted as follows. While the succesive estimates ûi,n are sufficiently close to each other, we continue the estimation process. As shown in Figure 2, the estimation process is continued, while new estimates belong to the intersection of estimated confidence intervals [û i,n -ησ i,n , ûi,n + ησ i,n ]. Besides, it shows that we do not need to store all the previous estimates (û i,n ) n ≤n but only the intersection of confidence intervals, the last previous estimate and its variance for each pixel. Finally, the factor η can be easily chosen in the range [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF][START_REF] Wexler | Space-time video completion[END_REF] as justified in [START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF][START_REF] Katkovnik | Adaptive window size image denoising based on intersection of confidence intervals (ICI) rule[END_REF][START_REF] Kervrann | An adaptive window approach for image smoothing and structures preserving[END_REF].

Local similarity measure for patch selection

As an alternative to a geometry-based approach [START_REF] Katkovnik | Adaptive window size image denoising based on intersection of confidence intervals (ICI) rule[END_REF][START_REF] Ercole | Spatio-temporal pointwise adaptive denoising in video: 3d non parametric approach[END_REF], we prefer to use weights that allow us to select the data points in the neighborhood for filtering. This selection is based on the similarity between a given small image patch at point x i and patches at points x j belonging to the space-time neighborhood W i,n .

Considering small patches enables to capture the local intensity structure. The sum of square image differences is widely used for similarity measure between image patches. However, in order to take into account the local variance of the estimator we define the following symmetric distance between image patches ûi,n-1 and ûj,n-1 :

δ 2 ij = 1 2 (û i,n-1 -ûj,n-1 )V -1 i,n-1 (û i,n-1 -ûj,n-1 ) + (û i,n-1 -ûj,n-1 )V -1 j,n-1 (û i,n-1 -ûj,n-1 ) (6)
where V i,n-1 denotes the estimator covariance matrix. The two vectors ûi,n-1 and ûj,n-1 denote p × p patches respectively centered on point x i and point x j . The two matrices V i,n-1 and V j,n-1 are diagonal with the diagonal elements equal to σ2 i,n and σ2 j,n respectively. We decide that the two vectors ûi,n-1 and ûj,n-1 are similar with a probability of false alarm 1 -α, under the hypothesis that they are Gaussian distributed, using a classical χ 2 test with p 2 degrees of freedom. In other words, when δ 2 ij /λ < 1, with λ choosen as a quantile of a χ 2 p 2 ,1-α distribution we can decide that the two patches are similar. In our experiments, we use a confidence level of 99% for setting α to 0.01 .

The distance δ 2 ij is transformed into a similarity measure using the exponential kernel. We compute the similarity measure for all the points of the neighborhood and normalize it to obtain the following expression for the weights:

ω ij = e -δ 2 ij /2λ xj ∈Wi,n e -δ 2 ij /2λ . ( 7
)
Comment citer ce document : Boulanger, J., Kervrann Then, if the distance δ ij between two patches is important then the weight ω ij associated to pixel x j is small and the pixel will not participate in the intensity estimation at point x i . Consequently, the weights provide an efficient and flexible way to implicitly choose the appropriate pixels contributing to the intensity estimation in the adaptive space-time neighborhood while effectively preserving space-time discontinuities. Note that the process is entirely data-driven and does not require particular geometry adapted to image content.

Algorithm implementation

At the initialization of the algorithm, we have to estimate the noise variance τ 2 . It can be robustly estimated by calculating pseudo-residuals ε i as described in [START_REF] Gasser | Residual variance and residual pattern in nonlinear regression[END_REF]. If we choose a space-time 6-neighborhood (4 spatial neighbors and 2 temporal ones), pseudo-residuals are compactly represented by ε i = (8Y i -∆Y i ) / √ 42 where ∆Y i is the discrete Laplacian at x i and the constant √ 42 is introduced to ensure that E[ε 2 i ] = τ 2 . Given the residuals ε i , we can then robustly estimate the noise variance τ 2 by: τ = 1.4826 med i (| ε i -med j |ε j | |). In our experiments, λ is set to the 0.99 quantile of the χ 2 k,0.99 distribution with k ∈ {(2p + 1) 2 , p ∈ {1, 2, 3, 4}}. The last parameter η is set to 2 √ 2 to ensure a good accuracy of the estimation. During the estimation, we alternate the increasing of the spatial and temporal extents of the space-time neighborhoods. Furthermore, the algorithm can be easily parallelized. Actually, estimation steps use only local information and have been distributed over 8 CPUs dividing the computation time by 8. Finally, another possibility to speed up the algorithm is to use a dyadic scheme when increasing the extent of the neighborhood. To conclude, the proposed method is very simple to implement and do not require the fine adjustment of parameters λ and η which control the estimatation process.

Experimental results

In this section, we report a great number of experiments to show how our exemplar-based space-time adaptive filtering method performs on real image sequences with artificially added noise. The described method is compared to others commonly-used methods for denoising image sequences. The performance of the algorithm combined with a global affine motion compensation stage is also analyzed. For an objective performance evaluation, we use the Peak-Signal-to-Noise-Ratio global measure defined as P SN R = 20 log 10 (255/mse) where mse denotes the mean squared error between the original noise-free image sequence and the denoised image sequence.

Performance assessment and comparison with other recent methods

We first report experiments to evaluate the influence of the noise level and the patch size on the overall method performance. Figure 3 Fig. 3. Denoising performance for several noise levels and for several patch sizes. The P SN R is used to measure the overall performance of the filtering and the test sequence is Akyio ( 176 × 144 × 300). This sequence mainly exhibits low motion. Experiments show that the introduction of patches improve the P SN R of at least 2dB. We can also remark that the results for patch 7 × 7 and 9 × 9 are similar.

tained for 8 noise levels and 5 patch sizes. First, we can note that the exemplarbased method is stable while the noise level increases. The improvement gained by introducing patches (to be compared to pixel-wise version) is clearly demonstrated. As it could be expected, it is useless to consider too larges patches. When the size of the patch increases, the P SN R increases too, however, the result for size 7 × 7 and 9 × 9 are quite similar while the computation time is proportional to the number of points in the patch.

In Figure 4, we can observe the improvement brought by the progressive increasing neighborhood scheme in term of P SN R as well as in term of visual quality. Figure 4(c) contains denoised image sequence with the same weighted intensity estimation as in our method but with a fixed 11×11×11 neighborhood size which corresponds to the biggest one used to filter the image sequence with our method whose result is shown in Figure 4(d). The image sequence denoised using our method is smoother and its edges are better preserved compared to the sequence denoised with a fixed neighborhood which appears as blurred. This is due to the fact that the weights are calculated from a regularized version of the image sequence and not the original noisy sequence.

We have compared our method with 4 other recent methods: a combination of a spatial Wiener filter with a motion-compensated temporal Kalman filter [START_REF] Dugad | Noise reduction in video by joint spatial and temporal processing[END_REF], a space-time non linear adaptive K-NN filter [START_REF] Zlokolica | Motion and detail adaptive denoising in video[END_REF], a 3D wavelet-based method [START_REF] Rajpoot | Adaptive wavelet restoration of noisy video sequences[END_REF] and a 3D point-wise adaptive estimate using a geometrical design for the neighborhood [START_REF] Ercole | Spatio-temporal pointwise adaptive denoising in video: 3d non parametric approach[END_REF]. 8 test sequences corrupted with noise of different levels are used. For a fair evaluation, we have considered the results supplied by the authors themselves in the referenced papers. Therefore, we cannot provide the P SN R measure for all the test sequences. Table 1 contains all the available results. Our method clearly outperform all the other methods since it supplies the best P SN R results for all the test sequences sometimes with a quite significative improvement (up to 4dB). We have also to stress that the implementation of our method is easy and it involves no parameters tuning.

Comment
Figure 5 reports another experiments on the "Flower Garden" sequence. It involves an important apparent (varying) motion due to the camera motion. In order to give insights into the spatio-temporal behavior of the denoising method, we have displayed XT slices of the image sequence. Reported results demonstrate that our method can cope with the presence of motion while preserving as well as temporal discontinuities.

Comment citer ce document : Boulanger, J., Kervrann and Wiener denoising with motion compensation using dense motion field [START_REF] Dugad | Noise reduction in video by joint spatial and temporal processing[END_REF], (b) Adaptive K-NN space-time filter [START_REF] Zlokolica | Motion and detail adaptive denoising in video[END_REF],(c)Wavelet based method for image sequence denoising: TIWP3D [START_REF] Rajpoot | Adaptive wavelet restoration of noisy video sequences[END_REF], (d) 3D non-parametric regression approach [START_REF] Ercole | Spatio-temporal pointwise adaptive denoising in video: 3d non parametric approach[END_REF], (e) the proposed adaptive method with 7×7 patches and 6 iterations. Numerical results for the other method are taken from the related publications.

Filtering with motion compensation

In this section, we analyze the case of image sequences with a global camera motion and we investigate the addition of a motion compensation stage.

Global parametric motion can be accurately computed using a robust estimation method even if independent moving objects are present in the scene as described in [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF].

In Figure 6, we give the overall block-diagram of the resulting method forming by adding a motion compensated stage to our denoising method. Motion estimation and denoising are not performed alternatively. Actually, we consider that the global affine motion is correctly estimated with the method described in [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] (We used the software made available on the web). A similar exploitation of a parametric motion compensation was proposed in [START_REF] Dekeyser | Spatio-temporal Wiener filtering of image sequences using a parametric motion model[END_REF] and associated with a 3D Wiener filtering technique. After the motion estimation step, the whole sequence is back-warped in the frame (referential) of the first image. Then, we apply our denoising method to this motion-compensated sequence. The result of the filtering has to be warped again to recover the images into their original frames. This is accomplished by using a two-step procedure. First, the denoised images are warped by using the initially estimated motion. Then, resulting images undergo a second transformation inferred from motion estimated between each resulting denoised image and its corresponding image of the original sequence. The Comment citer ce document : Boulanger, J., Kervrann Fig. 6. Block diagram of the denoising method with motion compensation. A motion compensation is applied to the image sequence before the denoising stage. It corresponds to the dotted box on the left. This module is composed of the estimation of affine motion models and the projection of the image sequence in the coordinates system of the first image. After the denoising stage, a two-step warping procedure is applied to get the images in their initial frames. The second step estimates the motion between each image of the input sequence and its corresponding image in the denoised warped image sequence.

later step is made necessary by the error accumulation in combining the motion models estimated from frame to frame in the warping procedure. Figure 7(c) shows one image the sequence denoised without motion compensation and Figure 7(d) shows the denoising result with the proposed motion compensation stage. The processed image sequence contains two cars tracked by the camera from an helicopter. The warping stage involved in the motioncompensated version of the algorithm require bilinear interpolations that somewhat blur the image sequence. In this expermient, the P SN R value of the proposed motion-compensated method is thus slightly lower than the P SN R value of the output of the denoising method without motion compensation.

We have shown that interpolations involved in the warping stage can induce blur effects that reduce the P SN R value of the filtered image sequence. Then, if interpolations could be avoided, the results of the motion-compensated denoising method would be improved. This could be achieved by transforming the neighborhood instead of warping the data. Such an investigation is in progress. In Figure 8, we have simulated an image sequence by sliding a window over an aerial photography with translation of 10 pixels between successive frame. Since all the motions are translations of entire pixels, no interpolation is required in the warping steps. This experiment shows that the motion-compensated denoising method performs better in the case of large displacements. Consequently, the denoising method without motion-compensation is more appropriate for denoising image sequences with not too large motion magnitude.

Conclusion

We have described a novel and efficient unsupervised method for denoising image sequences framework. It can specify, in a simple data-driven way, the most appropriate space-time neighborhood and associate weights to select the data points involved in the intensity estimation at each pixel. Moreover, it involves an exemplar-based approach extended to the space-time domain. All the parameters of the algorithm are well calibrated and our method does not require any fine tuning. Quite satisfactory results have been obtained on several image sequences. Furthermore, it was experimentally demonstrated that our method outperforms other recent methods. The visual quality of the denoised image sequences is noticeable since noise is well smoothed out while spatial and temporal discontinuities are well preserved. Finally, some improvements are proposed to incorporate a motioncompensated stage It was shown that motion-compensation was really needed for large displacement only provided interpolations involved in the warping steps could be avoided which is under progress.

Comment citer ce document : Boulanger, J., Kervrann Gaussian white noise of variance 100, the P SN R is 28.11dB, (c) the same image of the sequence denoised with the proposed method without motion compensation, the P SN R is 31.59dB , (d) denoised sequence using the described motion-compensated denoising method with a P SN R = 32.97dB.

  The weights are defined by computing the Comment citer ce document : Boulanger, J., Kervrann, C., Bouthemy, P. (2006). Adaptive space-time patch-based method for image sequence restoration. In: Workshop on statistical methods in multi-image and video processing (SMVP'06) (p. 14). Presented at SMVP'06 Statistical Methods in Multi-Image and Video Processing, Graz, AUT (2006-05-07 -2006-05-13).

Fig. 1 .

 1 Fig. 1. Exemplar-based space-time approach. To each point of the image sequence is associated an estimated space-time neighborhood Wi,n. To each point xj of this neighborhood is associated a 3 × 3 patch. The weights ωij are defined as a function of the distance between the patch centered in point xi and patches centered in point xj ∈ Wi,n.

3. 1

 1 Space-time neighborhood geometry One important contribution of this work is the on-line adaptation of the spacetime neighborhood sequence (W i,n ) n∈[0:N ] . We consider a simple hyper-cubic Comment citer ce document : Boulanger, J., Kervrann, C., Bouthemy, P. (2006). Adaptive space-time patch-based method for image sequence restoration. In: Workshop on statistical methods in multi-image and video processing (SMVP'06) (p. 14). Presented at SMVP'06 Statistical Methods in Multi-Image and Video Processing, Graz, AUT (2006-05-07 -2006-05-13).

Fig. 2 .

 2 Fig. 2. (a) Spatio-temporal neighborhood: colors correspond to iterations plotted in (b); (b) confidence intervals: circles represent estimates ûi,n obtained at each iteration n. The gray rectangles represent the intersection between the current confidence interval and the previous one. As long as the estimate belongs to this intersection, the estimation process is updated.

  , C., Bouthemy, P. (2006). Adaptive space-time patch-based method for image sequence restoration. In: Workshop on statistical methods in multi-image and video processing (SMVP'06) (p. 14). Presented at SMVP'06 Statistical Methods in Multi-Image and Video Processing, Graz, AUT (2006-05-07 -2006-05-13).

  citer ce document : Boulanger, J., Kervrann, C., Bouthemy, P. (2006). Adaptive space-time patch-based method for image sequence restoration. In: Workshop on statistical methods in multi-image and video processing (SMVP'06) (p. 14). Presented at SMVP'06 Statistical Methods in Multi-Image and Video Processing, Graz, AUT (2006-05-07 -2006-05-13).

Fig. 4 .

 4 Fig. 4. (a) Original sequence (b) noisy sequence with a additive Gaussian white noise of standard deviation τ = 30 (P SN R = 18.58) (c) denoised sequence using a fixed 11 × 11 × 11 neighborhood and a 5 × 5 patch. The P SN R is 28.72dB and we can see that the sequence is quite blurred (d) denoised sequence with the proposed method. The P SN R is 32.14dB.

Fig. 5 .

 5 Fig. 5. Sequence "Garden". (a) original sequence, (b) noisy sequence with a additive Gaussian white noise of standard deviation τ = 30 ( P SN R = 18.58dB ) and (c) denoised sequence, P SN R = 23.59dB. (d), (e) and (f) represent corresponding XT slices in space-time domain. The camera motion appear as lines in the XT slices.

  . The proposed method is based upon an adaptive estimation statistical Comment citer ce document : Boulanger, J., Kervrann, C., Bouthemy, P. (2006). Adaptive space-time patch-based method for image sequence restoration. In: Workshop on statistical methods in multi-image and video processing (SMVP'06) (p. 14). Presented at SMVP'06 Statistical Methods in Multi-Image and Video Processing, Graz, AUT (2006-05-07 -2006-05-13).

Fig. 7 .

 7 Fig. 7. Sequence "Avenger". (a) original sequence (b) noisy sequence with a additive Gaussian white noise of variance 100 (c) denoised sequence with the proposed method P SN R = 31.36 (d) denoised with the proposed motion-compensated method P SN R = 31.04. (e) and (f) are details of respectively (c) and (d). Some spatial blurring effects due to interpolation are visible in (f)

  , C., Bouthemy, P. (2006). Adaptive space-time patch-based method for image sequence restoration. In: Workshop on statistical methods in multi-image and video processing (SMVP'06) (p. 14). Presented at SMVP'06 Statistical Methods in Multi-Image and Video Processing, Graz, AUT (2006-05-07 -2006-05-13).

Fig. 8 .

 8 Fig. 8. Image sequence simulated by sliding a window in an aerial photography. (a) one image of original sequence, (b) the same image of the noisy sequence with a additiveGaussian white noise of variance 100, the P SN R is 28.11dB, (c) the same image of the sequence denoised with the proposed method without motion compensation, the P SN R is 31.59dB , (d) denoised sequence using the described motion-compensated denoising method with a P SN R = 32.97dB.

  , C., Bouthemy, P. (2006). Adaptive space-time patch-based method for image sequence restoration. In: Workshop on statistical methods in multi-image and video processing (SMVP'06) (p. 14). Presented at SMVP'06 Statistical Methods in Multi-Image and Video Processing, Graz, AUT (2006-05-07 -2006-05-13).

Table 1 .

 1 PSNR results for 8 test sequences and 5 denoising methods. (a) Join Kalman

	Sequence name P SN R (a) (b) (c) (d)	(e)
	Akiyo	22	-	-	-33.86 34.31
	Salesman	28 34.4 32.5 -	-35.13
		24 31.1		-	-32.60
	Garden	28	-28.2 -	-31.33
	Miss America 1 28	-35.3 -	-39.39
	Miss America 2	7	-	-26.36 -26.69
	(128 3 )	12	-	-28.16 -29.63
		17	-	-30.46 -32.05
		22	-	-32.66 -34.20
	Suzie	28 34.8 -	-	-37.07
		24 32.0 -	-	-35.11
	Trevor	28 33.9 34.1 -	-36.68
		24 31.3 -	-	-34.79
	Foreman	28 33.9 -	-	-34.94
		24 31.1 -	-	-32.90

, C., Bouthemy, P. (2006). Adaptive space-time patch-based method for image sequence restoration. In: Workshop on statistical methods in multi-image and video processing (SMVP'06) (p. 14). Presented at SMVP'06 Statistical Methods in Multi-Image and Video Processing, Graz, AUT (2006-05-07 -2006-05-13).

  , C., Bouthemy, P. (2006). Adaptive space-time patch-based method for image sequence restoration. In: Workshop on statistical methods in multi-image and video processing (SMVP'06) (p. 14). Presented at SMVP'06 Statistical Methods in Multi-Image and Video Processing, Graz, AUT (2006-05-07 -2006-05-13).

					image sequence
	Input	Motion estimation			Affine motion parameters
		Back-warping	Denoising	Warping	Motion estimation
		Motion compensation			Back-warping	Output
			1	2	Warping