Antoine Manzanera
email: antoine.manzanera@ensta.fr

Σ-∆ background subtraction and the Zipf law

Keywords: Image processing, Motion detection, Background subtraction, Σ-∆ modulation, Vector data parallelism for(t=1, t<=tstop, t+=1) {// Time loop

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Background subtraction is a very popular class of methods for detecting moving objects in a scene observed by a stationary camera [START_REF] Karmann | Moving Object Recognition Using an Adaptive Background Memory[END_REF] [2] [3] [4] [START_REF] Cheung | Robust techniques for background subtraction in urban traffic video[END_REF]. In every pixel p of the image, the observed data is a time series I t (p), corresponding to the values taken by p in the video I, as a function of time t. As only temporal processing will be considered in this paper, the argument p will be discarded. The principle of background subtraction methods is to discriminate the pixels of moving objects (the foreground) from those of the static scene (the background), by detecting samples which are significantly deviating from the statistical behaviour of the time series. To do this, one needs to estimate the graylevel distribution with respect to time, i.e. f t (x) = P (I t = x). As the conditions of the static scene are subject to changes (lighting conditions typically), f t is not stationary, and must be constantly re-estimated. For the sake of computational efficiency, which is particularly critical for video processing, it is desirable that f t should be represented by a small number of estimates which can be computed recursively. In this paper, we focus on recursive estimation of mono-modal distributions, which means that we assume that the time series corresponding to the different values possibly taken by the background along time, presents one single mode. This may not be a valid assumption for every pixel, but it does not affect the interest of the principle since the technique presented can be extended to multi-modal background estimation.

The Σ-∆ background estimation [START_REF] Mcfarlane | Segmentation and tracking of piglets in images[END_REF] [7] [START_REF] Manzanera | A new motion detection algorithm based on Σ-∆ background estimation[END_REF] is a simple and powerful non linear background subtraction technique, which consists in incrementing (resp. decrementing) the current estimate by an constant value if it is smaller (resp greater) than the current sample. Our objective is to discuss the foundations of this method, with regards to statistical estimation. We show the relation between the Σ-∆ estimation and the probabilistic model of Zipf-Mandelbrot, and compare it with two other recursive methods: exponential smoothing and Gaussian estimation. Section 2 presents the general framework of recursive estimation. Section 3 presents the particular case of Σ-∆ estimation, and provides the full numerical algorithm to compute it in the mono-modal case. Section 4 shows some results and discuss the computational advantages of Σ-∆ background subtraction, proposing in particular a complete vector data parallel implementation adapted to the SIMD-within-register framework. Section 5 concludes and presents the possible extension of the primitive algorithm.

Recursive estimation

If one should represent f t by one single scalar estimate, one of the most natural would be the average M t of the time series I t . The naive recursive computation: M t = 1 t I t + t-1 t M t-1 can be used as initialisation for the small values of t, but is not numerically realisable for long series. So one common solution is to use a constant weight (or learning rate) α ∈]0, 1[ for the current sample:

M t = αI t + (1 -α)M t-1
. This is sometimes referred to as running average, and corresponds to the recursive implementation of exponential smoothing.

One way of generalising this is to write the updating equation in an incremental form: M t = M t-1 + δ t (I t ), where δ t is the increment function, depending on the current sample I t . In the case of exponential smoothing, δ t is the affine function α(I t -M t-1 ) (Fig. 1(1)). This linear dependence is not satisfying, since, in most cases, a sample which is far from the average is out of the background model and should have little effect on the estimate updating. This problem can still be addressed in the exponential smoothing framework, by using two distinct constants α 1 and α 2 such that α 1 > α 2 , and by defining δ t (I t ) = α 1 (I t -M t-1 ) if I t is in the background model, and δ t (I t ) = α 2 (I t -M t-1 ) if I t is foreground. This results in a discontinuous increment function δ t , as shown in Fig. 1(2), where the decision background/foreground is done by simply thresholding the absolute difference: The pixel is foreground if |I t -M t-1 | > T h. It appears however, that the discontinuity of δ t makes the choice of T h critical.

To get a more continuous behaviour, we shall follow [START_REF] Stauffer | Learning patterns of activity using real-time tracking[END_REF], who suggests that the weight α attached to the current sample I t should depend on its probability f t (I t ). But, as noted by [START_REF] Power | Understanding background mixture models for foreground segmentation[END_REF], the mere product

δ t (I t ) = f t (I t ) × α(I t -M t-1 )
suggested by [START_REF] Stauffer | Learning patterns of activity using real-time tracking[END_REF] is not usable in practise because of increments too small in general to be numerically operative. A proper way to achieve this, if α max is the maximal desired weight, and as M t-1 is the mode of the current distribution f t , is to use:

δ t = α max f t (I t ) f t (M t-1 ) × (I t -M t-1 ). (1)
If we use a Gaussian distribution as density model like in [START_REF] Stauffer | Learning patterns of activity using real-time tracking[END_REF], we get the following increment function:

δ t = α max × exp( -(I t -M t-1 ) 2 2V t-1 ) × (I t -M t-1 ). ( 2
)
The model needs the temporal variance V t . In [START_REF] Stauffer | Learning patterns of activity using real-time tracking[END_REF], it is computed as the Gaussian estimation of the series (I t -M t ) 2 . But this leads to a double auto-reference in the definitions of M t and V t , which is prejudicial to the adaptation capability of the algorithm. We recommend rather to compute V t as the exponential smoothing of (I t -M t ) 2 , using a fixed learning rate α V .

One of the interest of computing V t is that it provides a natural criterion of decision background/foreground through the condition |I t -M t | > N × √ V t , with N typically between 2 and 4. Note that the increment function (Fig. 1(3)) is very different from the previous ones, and has a derivative-of-Gaussian shape.

This Gaussian estimation provides some attractive features compared to the exponential smoothing: the update of the estimates depends on the probability of the current sample, and the increment values are globally higher when the background density is more dispersed. Nevertheless, it is less used than exponential smoothing because of the computational load much higher. Now, what does the increment function look like if we take the Zipf law as the probabilistic model ?

Zipfian background estimation

Originally the Zipf law is an empirical principle [START_REF] Zipf | Human behavior and the principle of least-effort[END_REF] at the crossroads of linguistic and information theory, stating that, in any sense-making text, the probability of occurrence of the n th most frequent word is 1/n the probability of occurrence of the (first) most frequent word. So the Zipf distribution is a hyperbolic decreasing function. Recently, it has been used in several applications of image processing [START_REF] Caron | A method for detecting artificial objects in natural environments[END_REF], in particular as a model for the distribution of local spatial features. We use it here as a model for (pixel-wise) temporal distribution.

Because of the divergence of the sum 1/n, the Zipf density function includes a power factor: 1/n s , with s > 1. The general expression of the continuous symmetric Zipf-Mandelbrot distribution can be written:

Z (µ,k,s) (x) = (s -1)k s-1 2(|x -µ| + k) s . (3)
In this expression, the parameter µ represents the mode of the distribution, and k determines its dispersion. The remarkable property of Z, taken as the density model f t of the time series I t (and then, replacing eq. 3 in eq. 1), is the shape of the increment function δ t (Fig. 1(3)), which is close to the Heaviside shaped function:

H (µ,κ) (x) = -κ if x < µ, +κ if x > µ, with κ = α max k s .
Thus it is naturally related to the Σ-∆ modulation, classically used in Analog to Digital Conversion:

Results

Figure 2 shows the background estimation for all the time indexes, and one particular pixel. This is a difficult case for pixel-level motion detection: an outdoor scene where the background signal (high grass meadow with wind) is corrupted by the passage of two foreground objects. The Boolean condition "D t > V t " is used as foreground classification. Figure 3 shows the result for all the pixels, at 4 different time indexes of the classical Hall sequence, in which two people are moving in radial motion, i.e. in the direction of the optical axis. This is a difficult case too, since the values in the centre of the moving objects do not change much (aperture problem). For reference, the last row of Figure 3 displays the hand drawn ground truth for the 4 frames.

This ground truth is used for the quantitative evaluation (detection and false alarm rates are averaged on these 4 key frames) shown on Table 1, for different values of the amplification constant N, and of the updating period T V . Those results are resumed on Figure 4, where the 9 Σ-∆ algorithms are compared with 6 different Gaussian algorithms. Note that these figures relate to pixel-level methods, and should not be interpreted in absolute, since a simple spatial regularisation appreciably improves the two measures, in all cases.

N=1

N=2 N=4 TV = 1 (0.74,0.25) (0.62,0.10) (0.53,0.02) TV = 8 (0.91,0.38) (0.87,0.23) (0.85,0.12) TV = 32 (0.95,0.45) (0.94,0.38) (0.94,0.33) Table 1. (Detection, False alarm) rates for 9 Σ-∆ background subtraction algorithms. Measures are averaged on the 4 key frames of the Hall sequence.

The relevance and power of the Σ-∆ estimation, as a pixel-level temporal filter, is comparable to that of the Gaussian estimation, whereas its computational cost is even inferior to that of exponential smoothing. Indeed, the algorithm proposed in Section 3 is straightforward to compute in any fixed-point arithmetic, using an instruction set limited to: absolute difference, comparison, increment/decrement. Thus, it is well adapted to real-time implementation using dedicated or programmable hardware.

Another important computational property of Σ-∆ background subtraction, is that, once chosen the number of bits used to represent the estimates M t and V t , every computation can be made at full precision without increasing the data width. This allows in particular to make the most of the data parallelism provided by the SIMD-WAR (Single Instruction Multiple Data Within A Register) paradigm, which consists in concatenating many short operands in one single very long variable, and then applying scalar operations on the long variables. This implementation is available on most personal computers, using for example the SSE-2 (Streaming SIMD Extensions 2) instructions of the Intel R C++ compiler [START_REF] Intel | Intel R C++ Compiler for Linux Systems -User's Guide[END_REF]. We provide hereunder the vectorised pseudo-code of the Σ-∆ background subtraction. Here, a 16-times acceleration is achieved by performing the operations on a 128-bit register made of 16 8-bit data. the background updating to the temporal dispersion. We have shown the computational advantages of the Σ-∆ estimation, illustrated by the vector SIMD implementation.

The limitations of this algorithm -used "as is" in a motion detection systemare inherent to its mono-modal nature: first, one single mode in the density can be inefficient to discriminate moving objects over a complicated background, and second, one single dispersion estimate, related to one time constant, may not be sufficient for certain kind of motion such as remote objects with radial velocity w.r.t. the optical centre. Nevertheless the basic model can be enriched, either by using a multi-modal Zipfian distribution like it is done in [START_REF] Stauffer | Learning patterns of activity using real-time tracking[END_REF] for Gaussian estimation, or by using several time magnitudes, as shown in [START_REF] Manzanera | A new motion detection algorithm based on Σ-∆ background estimation[END_REF].

  , Th=30, α 1 =1/16, α 2 =1/32

( 4 )Fig. 1 .

 41 Fig. 1. The different increment functions δt associated to the different distribution models: (X axis: graylevel It, Y axis: increment value δt(It)). (1) Exponential smoothing (plain: α = 1/32; dashed: α = 1/16) (1) Bi-level exponential smoothing (m = 100, T h = 30, α1 = 1/16, α2 = 1/32) (3) Gaussian laws, αmax = 1/4 (plain: µ = 100, σ = 30; dashed: µ = 150, σ = 50) (4) Zipf laws, αmax = 1/4 (plain: µ = 100, k = 1, s = 1.1; dashed: µ = 150, k = 5, s = 1.1).

Fig. 2 .

 2 Fig. 2. Σ-∆ background estimation running on a given pixel. (X axis: time index, Y axis: graylevel). All values are 8-bit coded. Amplification factor N is 2. Variance updating period TV is 1. Plain line: It, Dotted line: Mt, Impulses: Dt, Dashed line: Vt.

Fig. 3 .

 3 Fig. 3. Background subtraction shown at different frames of the Hall image sequence. Row 1: Original sequence, frames 38, 90, 170 and 250. Rows 2 and 3: Σ-∆ Background and foreground, with N=2, and TV = 8. Row 4: (Fore)ground truth.

Fig. 4 .

 4 Fig. 4. Detection / false alarm rates diagram, for 9 Σ-∆ and 6 Gaussian background subtraction algorithms. Measures are averaged on the 4 key frames of the Hall sequence.

For every time step ∆t: If M t-∆t > I t then M t = M t-∆t -ε ; If M t-∆t < I t then M t = M t-∆t + ε ; Here, the average increment per time unit is κ = ε ∆t . Digitally, the elementary increment ε is the least significant bit of the representation, i.e. 1 if the values are integer-coded. Adaptation to the dispersion of the model can then be done by tuning the updating period ∆t: the greater the variance, the smaller ∆t should be. The following algorithm reproduces such behaviour. The principle is to attach to every pixel, in addition to the mode estimator M t , a dispersion estimator V t . Suppose that V t ∈]0, 2 m -1[, which means that it is coded on m bits:

For every frame t:

Here x%y is x modulo y. The purpose of the two first lines of the algorithm (which are computed once for all the pixels at every frame) is to find the greatest power of two (pow2) that divides the time index modulo 2 m (rank). Once this has been determined, it is used to compute the minimal value of V t-1 for which the Σ-∆ estimate M t will be updated. Thus the (log-)period of update of M t is inversely proportional to the (log-)dispersion: if V t > 2 m-1 , M t will be updated every frame, if 2 m-2 ≤ V t < 2 m-1 , M t will be updated every 2 frames, and so on.

The dispersion factor V t is computed here as the Σ-∆ estimation of the absolute differences D t , amplified by a parameter N . Like in Gaussian estimation, we avoid double auto-reference by updating V t at a fixed period T V . V t can be used as a foreground criterion directly: the sample I t is considered foreground if D t > V t . V min and V max are simply used to control the overflows ; 2 and 2 m -1 are their typical values.

Note that the time constants, which represent the period response of the background estimation algorithm, are related here to the dynamics (the number of significant bits) of V t , and to its updating period T V . For the other methods, the time constants were associated to the inverse of the learning rates: 1/α i for the exponential smoothing and 1/α max and 1/α V for Gaussian estimation.


Paul Nadrag
email: paul.nadrag@ensta.fr

Antoine Manzanera
email: antoine.manzanera@ensta.fr

Nicolas Burrus
email: nicolas.burrus@ensta.fr

Smart retina as a contour-based visual interface

The purpose of this work is to provide a robust visionbased input device. In our system, a programmable retina is looking at the user who sends commands by moving his hand. The fusion between the acquisition and the processing functions of the retina allows a close adaptation to the lighting conditions and to the dynamic range of the scene. Thanks to its optical input and massive parallelism, the retina computes efficiently the contours of the moving objects. This feature has nice properties in terms of motion detection capabilities and allows a dramatic reduction in the volume of data to be output of the retina. An external low-power processor then performs global computations on the output data, such as extreme points or geometric moments, which are temporally filtered to generate a command.

Introduction

A vision-based interface is a system allowing to send symbolic commands to a computer just as a mouse or a keyboard does, but with visual interaction only. Using a camera as an input device is an attractive issue for numerous applications, from disabled persons assistance to intelligent home environments. To be really useful, such systems must not be invasive at all: it must neither impose to the user the wearing of any equipment, nor induce any constraint in the decoration of the scene observed by the camera. But detecting or recognizing the visual features independently of the scene content and its variation is a difficult task.

We are using in this work a vision system composed of a smart camera (programmable retina) and an external lowpower processor (ARM). The interest of this approach is twofold:

1. The intimate combination of acquisition and processing functions of the retina allows to maintain a high fitness between the scene lighting properties and the further processing: in our case, the acquisition is adapted according to a local contrast criterion, in order to enhance the stability of the contours, which are the visual features chosen to detect the interaction.

2. The massive data reduction which is made between the optical input of the retina and the contour descriptor manipulated by the external processor is in perfect adequation with the I/O flow of the visual interface: image in, symbol out. By minimizing data transfers thanks to the massive parallelism and to the global measure output, the system is designed to optimize this flow from an energetic point of view. Section 2 presents the Programmable retina within the vision system, its instruction set, acquisition and programming paradigms. Section 3 deals with image processing algorithms: the moving contour detection is detailed, and its properties are discussed. In Section 4, the higher level processing is presented and our first experiments in a drawing interface application are exposed. Finally, Section 5 outlines the contributions and presents some perspectives of this work.

Programmable retina

The programmable retina concept originates from the NCP (Neighborhood Combinatorial Processing) retinas [START_REF] Zavidovique | Bilevel processing of multilevel images[END_REF], which were SIMD Boolean machines. The NSIP (Near Sensor Image Processing) concept [START_REF] Forchheimer | Near-sensor image processing: a new paradigm[END_REF] then allowed to process graylevel images. Now, the highly submicronic level of CMOS technology allows to put more and more powerful processing circuitry aside the photoreceptors while preserving good acquisition performance and resolution. The most recent circuit that has been fabricated and validated is a 200 × 200 retina called Pvlsar 34, with an elementary digital processor and 48 bits of memory within every pixel.

The algorithm presented in this paper was implemented on the architecture presented in Figure 1). The retina Pvlsar 34 is both a CMOS sensor and a parallel machine. It is a grid of 200 × 200 pixels/processors connected by a regular 4-connected rectangular mesh. The processors operate synchronously, on their local data, a sequence of instructions being sent by the controller, which is the NIOS processor IP core of an Excalibur FPGA chip. The host processor or cortex is the ARM processor hardware core of the Excalibur. It can exchange data with the retina, modify the program sent by the NIOS to the retina, and is in charge of the higher levels of computation (i.e. non-image processing).

Every pixel/processor p of the grid G is composed of:

• one photoreceptor,

• one Analog to Digital Converter,

• 48 bits of digital memory {p i }, 1 ≤ i ≤ 48.

• one Boolean Unit (BU), which can read some bits of the digital memory, compute a Boolean operation and write its output on one bit of the digital memory. Regarding data extraction, there are two ways to output information from the retina:

0 0 0 0 0 1 1
• by transfering sequentially (translating) the image, to get the exact content of one or more bit planes of the digital memory.

• by using the Analog Global Summer, which provides in constant time an approximate measure of the number of 1s in any bit plane of the digital memory. Although simple, this last feature is important as it provides efficiently global measures that are very useful to get spatial statistics or to detect the convergence of relaxation algorithms.

The Analog to Digital Conversion is performed during the acquisition, and is a plain part of the retinal program. Indeed, the acquisition of a n graylevels image is done by multiple readings of the photoreceptor which provide (n -1) binary images (level sets) along the time (Figure 2). The photoreceptor (photodiode) is loaded to tension V max and then unloads linearly with time, with a slope proportional to the amount of photons received. The graylevel is then computed in every pixel by counting (on log 2 (n) bits), the number of times the tension at the bounds of the photodiode is below a threshold tension V s . Thus, an acquisition is defined by the set of time indices {x i } i=0..n-1 , corresponding to the reading instants of the photodiode.

The acquisition and digital coding of a n graylevels image is performed in n + log(n) operations, corresponding to the sequence of Boolean operations performed to add the (n-1) level sets. Of course, the resulting digital image occu- pies log(n) bits of memory in every pixel. Then, the retinal program applies a sequence of arithmetic and logic operations that are completely written in software, at the bit level. However, if the number of operations (and thus the energy consumption), is not critical, the possibility of processing during acquisition allows a lot of elegant computations of image operators. In particular, many operators of mathematical morphology are naturally computed by combining their binary counterparts over the level sets of the image (this is a consequence of the extension scheme of the mathematical morphology from sets to functions) [START_REF] Serra | Image analysis and mathematical morphology[END_REF]. In our application for example, the morphological gradient (see Section 3) can be computed by summing the level lines (i.e. the contours of every level set) instead of the level sets.

But a more immediate consequence of this mode of acquisition is that it allows to control very simply the histogram of the image, by modifying the reading instants of the photoreceptor (e.g. logarithmic compression, gamma correction). Furthermore, it can be coupled with the analog global summer to adapt the reading instants to the graylevels repartition (e.g. histogram equalization). In the following of this section, we describe the adaptive acquisition we actually perform in our system.

As we wish to describe the static background of the scene (and then, by discrimination, the moving objects) by their contours, we need to compute contours that are both rich (to get all the structure from the objects, without regard to their contrast) and temporally stable. Meeting these two requirements is very hard, in particular when coding a high dynamic scene with a limited number of bits, because in this case, the contrast within the extreme (i.e. very light or very dark) regions is inexistent (saturation effect).

On the retina, a triple acquisition is performed: at each frame time t, 3 graylevels images are coded B t , M t and D t corresponding respectively to bright, medium and dark regions of the graylevel distribution. Then 3(n -1) level sets are computed instead of n -1, corresponding to the 3 sets of time indices

{b t i } i=1..n-1 , {m t i } i=1..n-1
, and {d t i } i=1..n-1 respectively. These time indices are related as follows. First, the first index of the medium acquisition is temporally updated by a small increment:

if |M n/2 t | < (w × h)/2 then m t+1 1 = m t 1 + ε else m t+1 1 = m t 1 -ε where M n/2 t
is the (binary) n/2-level set of M t (or equivalently, the most significant bit of the binary code of M t ), |X| is the cardinality (the number of 1s) of X, w and h are the dimensions of the image (here w = h = 200). This first condition corresponds to adapting incrementally the time reference such that the median value of the medium image M t is half of the n graylevels dynamics. |M n/2 t |, which is the only global measure needed by this algorithm, is computed thanks to the analog summer.

Second, the first indices of the two other acquisitions are placed in the time axis according to the relation:

b t 1 + δ bm = m t 1 = d t
1 -δ md where δ bm (resp. δ md ) represents the time offset -manually fixed -between the bright and medium (resp. between the medium and dark) regions.

Finally, every image is coded according to a logarithmic representation of the graylevels:

x t i+1 = αx t i with x = b, m, and d respectively.

Figure 3(1-3) shows the 3 resulting images B t , M t and D t . Recomposing a high dynamic range image from these 3 images, or reducing the representation of the image (e.g. from 3log(n) to log(n) bits) can be done using different techniques that have been addressed in the literature [START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF] [START_REF] Fattal | Gradient domain high dynamic range compression[END_REF]. But in our application, we do not wish to recompose an image, but only to estimate the local contrast to compute the contour. We then simply compute a reduced representation of the contrast by computing the maximum of the local contrast in the 3 images:

G t B = max(g B (B t ), g B (M t ), g B (D t ))
where g B is the morphological gradient (see Section 3). Some other results that have been achieved by using the artificial retina can be found in the [START_REF] Burrus | Adaptive vision leveraging digital retinas: Extracting meaningful segments[END_REF] and [START_REF] Richefeu | A morphological dominant points detection and its cellular implementation[END_REF] papers.

Moving Contour detection

The static part of the algorithm used in the detection of the moving contours is based on the morphological edge detector described by Lee et al in [START_REF] Lee | Morphologic edge detection[END_REF]. This non-linear framework, using essentially min and max operations, is well suited to the Boolean computational model of the programmable retina. The morphological edge detection provides a measure of local contrast within every pixel, thanks to the dilation residue operator: let f (x, y) be a grayscale image and b the structuring element, which is a set containing the origin (0, 0). The dilation of f by b if defined as:

( f ⊕ b)(x, y) = max (i, j)∈b f (x -i, y -j)
The edges obtained by using the dilation residue, or (dilation-) gradient, noted by g b (x, y), will be:

g b ( f )(x, y) = ( f ⊕ b)(x, y) -f (x, y) Like [9]
, we combine different gradients obtained using several structuring elements:

g B ( f ) = min(g b 1 ( f ), ..., g b n ( f ))
, in order to be less sensitive to noise and edge orientation. And according to Section 2, we maximize the gradient over the three acquired images:

G t B = max(g B (B t ), g B (M t ), g B (D t )): see Figure 3(4).
Once the (grayscale) edge map computed, we get the (binary) contours by employing a non-maximal suppresion technique, adapted from [START_REF] Chanda | A multiscale morphologic edge detector[END_REF] to suit the massive parallelism of the retina. The method of suppresion is the following: let p be the pixel that is currently being processed (do not forget that the other 39,999 pixels of the retina are simultaneously undergoing the same sequence of operations) and his twelve neighbours, denoted using cardinal directions indices as follows:

p After this local maxima selection two last operations are completed to compute the contour image C t : first, the contour with a gradient value of 1 are discarded to avoid the fake edges induced by the quantization effect on a grayscale ramp (those for example, due to the vignetting effect of the lens). Then, to eliminate salt-and-pepper-like detections, we also impose the condition that a pixel can be taken into account as a contour pixel only if it has at least two neighbours classified as contours too. This gives a lower rate of false detections, although it lowers the quantity of information theoretically obtainable from an image. Figure 3(5) shows an example of contours image.

To discriminate the moving contours from those of the static scene, we need to compare the current occurrences with the past edges. The basic approach is the following one: we keep track of the number of times when a pixel was a contour, by computing, for any instant t, two descriptors attached to every pixel, and related to two different timescales: the short term history (STH) and the long term history (LTH). At every frame, if the pixel belongs to an edge, then STH is increased. If not, it is decreased. The increment/decrement is adjusted to not exceed the interval allowed by the dynamics of STH. The update for the LTH is almost the same, except that instead of the contour C t , the most significant bit of STH is used, and the timescale (i.e. the update period) is larger.

The computation of the moving contours is obtained from a logical AND between the current contour C t and the pixels with a low LTH. As, obviously, a pixel can belong to both the contour of a moving object and the contour of the static scene, one iteration of reconstruction (or geodesic dilation) is performed: a dilation is computed on the intersection set, followed by a logical AND with C t . This step allows to limit the disconnection on the resulting moving contours image E t . The relevance of this basic approach is based on the fact that -by construction -the contours are one-pixel-thick, and then an object must be very still to let his contours be taken into account by the LTH. To enhance this behavior, however, and avoid the blurry, low-value disturbances which may be caused by the contours of a moving object in STH, we lower the increment frequency in STH for the pixels belonging to

E t .
Another problem in the basic approach is generated by the fact that, generally, the moving objects are not transparent, which makes possible for a part of the background contours to be erased by the overlapping objects. The solution has been provided by the introduction of discrimination between the increment and the decrement of the LTH, the latter being discouraged. Thus, even if an object is obstructing the background scenery for a large amount of time (adjustable through parameterisation), it is possible for the object's contours to integrate themselves in the background (if the object is standing still) and, at the same time, for the background contours to stay longer in the scene, as if the new object was transparent. This property, together with the intrinsic thinness of the contours, significantly reduces the "ghost effect" (part of the static scene appearing as moving contours after being covered for a long time by a moving object). See Figure 3(6) for an example of moving contours image.

The measured times of execution for the different parts of the algorithm are the following: 10ms are needed by the retina to do its job (determination of the gradients, update of STH and LTH, extraction of moving contours, etc.), 8.6ms are needed for the ARM to compute the new trace (see Section 4) and 50ms are necessary to acquire the three initial images. The time needed to capture the three images is largely scene-dependent, a darker scene will demand a longer exposure time. As one can observe, the bottleneck is not induced by the processors in the retina themselves, but by the images' capture time. From an energetic point of view, the power consumption of the system is only a few tens milliwatts, 3 to 4 orders of magnitude as small as that of a PC.

Because of the high sensitivity induced by this algorithm, we get some false detections in addition to the true contours of the moving objects. Since some local filtering has already been applied, those false detections are gathered into small connected components. As objects of interest generally have important areas, we reduce the number of false alarms by filtering out small regions. This requires non-local computations on the image, surpassing the present possibilities of the retina. Thus, this step is performed by the host (ARM) processor. The further analysis, consisting in extracting and exploiting regional and global quantities on the moving contours is discussed in the next section.

Vision-based interface

The purpose of the image processing techniques presented above is to allow a stable and complete extraction of the contours of the moving objects, without imposing constraint of uniformity in the background. With respect to the input data flow, this feature represents a small amount of data which can be quickly extracted from the retina and processed by the external processor. We are currently beginning the exploitation of this visual feature for higher level processing. We present in this section the first experiments of visual interaction we have performed.

In these experiments, we are testing a drawing interface: the user draws a bidimensional curve by moving his finger or a pen above the retina, which is oriented toward the ceiling. Such interface can be used to input hand drawn char-

(1) (2)
(3) (4)

(5) (6) acters or sketches [START_REF] Erdem | Vision-based continuous graffiti-like text entry system[END_REF], but also some simple hand gestures [START_REF] Huang | Hand gesture modeling, analysis, and synthesis[END_REF], for which the curve represents a 2d model to be further processed and recognized.

One difficulty in developping visual interface is the necessity of several modalities in the interaction, just like the action of a mouse stops if it is dropped, or changes if its button is depressed. We wish to use the visual interface alone, without interaction of the mouse or keyboard. This implies that several parameters of the visual feature must be exploited, in order to define not only the position of the "tool", but also the action modality. In our example, the 3 modalities are (1) Draw, (2) Do nothing, (3) Erase.

So the first parameter computed from the visual feature is a global one ; it is the orientation of the moving object. Like [START_REF] Ahmad | A usable real-time 3d hand tracker[END_REF], we compute the second order moments (v x , v y ), i.e. the variances of the two coordinates of the moving contour points. The rough orientation is then used to define the modality: if v x /v y is superior to a given threshold σ D , then the "Draw" action is selected. If v y /v x is superior to a given threshold σ E , then the "Erase" action is selected. Otherwise, the "Do nothing" action is selected.

The second parameter extracted from the moving contour is a local one: the extreme left position of the moving objects, if any (there must be a connected curve with at least N min points) is used as the candidate position of the tool. This candidate position is checked according to a predictive filtering (it must not be too far from the previous position of the same tool). If the position is confirmed, then the following action is applied on the drawing widget: if "draw" is selected then create a line between previous and current position ; if "erase" is selected, then erase all points situated at the right of the current position.

On Figure 4, we can see the tracking system in action. Special attention should be paid to the transitions over the very textured areas (the neon grid, for instance). Because we work with contours only, we could expect some difficulties when crossing these areas, but by using a wisely chosen set of filters (enclosure and selection based on the elements' size), we are able to obtain a very high detection rate of the fingertip.

(1)

(

(3) (4) Figure 4. Drawing example using the visual mouse

Conclusion

We have presented in this paper a project of visual interface based on the extraction and processing of the contours of the moving objects. At the moment, we have mainly developped the low-level part of the project, which corresponds to the extraction of the moving contours, using a programmable retina, both as camera and parallel computer. The combination of acquisition and image processing allows a smart adaptation to the lighting properties of the scene. The visual interaction can be seen as a huge data reduction, turning video flow into symbols. The computational model of the programmable retina limits the data transfers to the minimum, and then reduces the energy consumed by this data reduction.

The higher level (interpretation of the moving contours) is being investigated, and will be pursued in the following months. We have presented our first experiments, based on single point/multiple actions interface. The perspectives of this work are also: global hand posture estimation for controlling a 3d mouse, and traking a 3d model to recover the complete posture of the arm/hand/fingers.

Figure 2 .

 2 Figure 2. Graylevel acquisition by multiple reading of the photodiode.

Figure 3 .

 3 Acquisition and processing within the programmable retina: (1) short exposure acquisition B t . (2) medium exposure acquisition M t . (3) long exposure acquisition D t . (4) maximum of the morphological gradient over the three acquisitions G B t . (5) binary contours C t . (6) moving contours detected E t .

  NW + p W + p SW < p N + p + p S > p NE + p E + p SE (corresponding to a horizontal ridge) 2. p NW + p N + p NE < p W + p + p E > p SW + p N + p SE (corresponding to a vertical ridge) 3. p NN + p NW + p WW < p NE + p+ p SW > p EE + p SE + p SS (corresponding to a ridge with an orientation of 135 • ) 4. p WW + p SW + p SS < p NW + p+ p SE > p NN + p NE + p EE (corresponding to a ridge with an orientation of 45 • )

	p WW	p NW p W p SW	NN p N p p S p SS	p NE p E p SE	p EE
	p will be selected as a contour point if at least one of the
	following four conditions is true:
	1. p


M Julien

M Antoine Manzanera

Détection de mouvement par capteur intelligent

Keywords: détection du mouvement, système temps réel, rétine artificielle, moyenne récursive, morphologie mathématique, modulation Σ-∆ motion detection, real-time system, artificial retina, recursive mean, mathematical morphology, Σ-∆ modulation

Cet article présente le portage de trois méthodes de détection du mouvement, temps réel, dans une séquence vidéo sur un système de vision à base de rétine numérique. Basé sur une architecture massivement parallèle, ce capteur intelligent permet d'effectuer des opérations élémentaires y compris pendant la phase d'acquisition de l'image. Les trois méthodes de détection du mouvement, ainsi que leur implémentation sur la rétine sont ensuite présentées. Pour chaque méthode, un effort d'explication est fait, et l'algorithme complet est fourni explicitement.

Introduction

Alors que les systèmes classiques de vision artificielle séparent complètement capture et traitement de l'image, ces deux fonctions s'avèrent en réalité intimement associées dans les systèmes visuels biologiques. Le caractère fortement submicronique des techniques d'intégration CMOS actuelles permet justement d'introduire des circuits de traitement de l'information au sein de chaque pixel d'un imageur, qui devient alors une rétine artificielle. La rétine numérique [START_REF] Paillet | Second Generation Programmable Artificial Retina[END_REF] [START_REF] Komuro | A digital vision chip specialized for high-speed target tracking[END_REF] est un circuit combinant fonctions d'acquisition et de traitement des images. Elle est constituée d'une grille bidimensionnelle de "cellules", chaque cellule étant formée de l'association d'un élément photosensible (photorécepteur), et d'un élément de calcul (processeur). Du point de vue algorithmique, la rétine numérique est une machine massivement parallèle de type réseau d'automates cellulaires, caractérisée par des éléments de calcul de capacité limitée (en terme de puissance de calcul et de mémorisation). Dans un système de vision, la rétine est associée à un hôte appelé cortex, typiquement composé d'un processeur scalaire (microprocesseur basse puissance DSP) et/ou d'un circuit de logique programmable (FPGA). Le cortex a deux fonctions : (1) contrôler la rétine, en lui envoyant les séquences de commandes; (2) effectuer les tâches de haut niveau du problème de vision à partir des informations extraites de la rétine. Un système de vision à base de rétine programmable doit permettre d'atteindre de grandes performances en terme de vitesse, grâce au parallélisme massif de la rétine, et aussi en terme de consommation d'énergie, de par la réduction drastique du flux d'information au sein du système. L'utilisation d'un tel système nous amène à une réflexion de fond sur les méthodes à employer : Comment sortir de l'information du circuit rétinien ? Sous quelles formes ? Quels sont les descripteurs appropriés à de tels calculs ? Comment adapter les algorithmes d'analyse d'images ? Comment mettre en place les mécanismes de communication entre la rétine et le cortex ? De par ses avantages, le couple (rétine-cortex) est particulièrement adapté aux systèmes de vision portables, mobiles, autonomes. Cela implique que la perception du mouvement doit être une caractéristique forte du système à base de rétine. Les applications sont nombreuses et variées : systèmes de télésurveillance, capteurs de robots mobiles, systèmes de poursuite automatique, caméras compressives... Nous présentons ici une étude sur l'implantation sur rétine programmable de trois algorithmes différents de détection du mouvement dans un contexte de caméra fixe. Les contraintes liées à ce système ont orienté les choix algorithmiques vers des méthodes récursives afin de réduire le coût en terme de Z(p) la plus grande puissance de 2 qui divise un naturel p, le codage de Gray du niveau de gris est défini par ses chiffres binaires (G i ) 0≤i≤log(n-1) avec :

G j = Z(i)=2 j C i (1)
Où représente l'imparité. On peut en outre passer facilement d'un codage Gray au codage naturel comme le montre la figure 2.

//Acquisition de l'image et codage en code Gray Le fait que l'algorithmique intervient aussi précocement dans les rétines programmables permet naturellement d'effectuer des traitements pendant l'acquisition de l'image, en effectuant des calculs à chaque ensemble de niveaux. Le cadre formel le plus adapté est donc celui des opérations morphologiques, grâce à l'équivalence entre morphologie en niveau de gris et morphologie ensembliste sur chaque ensemble de niveau. En outre, on peut intervenir sur la dynamique de l'image, en modifiant les instants de lecture de la photodiode (compression logarithmique par exemple) ou encore en adaptant les instants de lecture à la répartition des niveaux de gris (i.e. égalisation d'histogramme). de détection du mouvement Le défi de la détection du mouvement dans le cas de capteurs fixes tient dans la capacité d'effectuer une bonne segmentation des objets en mouvement indépendamment de leurs tailles, de leurs vélocités, ou de leurs contrastes par rapport au fond. Dans un tel contexte, l'élaboration des algorithmes de détection du mouvement se fait en tenant compte des contraintes suivantes :

Pour (i = 1; i < n; i + +) { //Seuillage du niveau i X i = ACQ(i) //Calcul de l'index: plus grande puissance de 2 qui divise i index = 0 x = i Tantque ((x%2) == 0) faire {index + +; x = x/2} G index = XOR(G index , X i ) } //Passage en code naturel N log(n) = G log(n) ; Pour (i = log(n) -1; i ≥ 0; i --) X i = XOR(G i , X i+1 );
• Le système doit être en mesure de fonctionner sans intervention humaine pendant un long moment, et être capable de prendre en compte des changements graduels ou soudains tels que les variations d'illumination ou la présence de nouveaux objets statiques dans la scène. Le système doit donc être temporellement adaptatif;

M t = αI t + (1 -α)M t-1 (2) = α(I t -M t-1 ) + M t-1 (3)
En pratique, pour éviter de multiplier les erreurs d'arrondis dues aux multiplications par α, on utilise la seconde formule. La Figure 3 montre les résultats du calcul de fond récursif pour un pixel particulier correspondant à une zone de passage d'un objet en mouvement. On remarque que la moyenne suit le signal selon une croissance ou une décroissance exponentielle.

Implémentation

A des fins d'implantation rapide, on limite la division à un simple décalage en prenant α = 2 -k . Ainsi pour une image acquise en n bits, l'intervalle de valeur pour lequel le calcul ait un sens est α ∈ [2 -(n-1) , 1]. La multiplication par α revient alors à effectuer un décalage à droite de k bits. Notons que grâce au bit de signe, lorsqu'on ajoute un nombre positif inférieur à 2 k , l'effet est nul alors que lorsqu'on ajoute un nombre négatif quelconque, on décrémente toujours au moins de 1. Nous menons donc le calcul une trame sur deux en ajoutant α * (X -M ) et une trame sur deux en soustrayant α * (M -X) (voir figure 4). Ce calcul est mathématiquement identique mais permet de tenir compte du problème lié aux arrondis. Le fait d'intervertir le sens de la différence une trame sur deux permet statistiquement de pouvoir ajuster au mieux le calcul du fond. En effet, à chaque instant t, la probabilité que l'image acquise X soit supérieure à la moyenne récursive M est égale à celle que X soit inférieure à M .

Gradient temporel morphologique oublieux

Principe

Le second exemple présenté est un opérateur différentiel basé sur un filtre hybride combinant opérations linéaires et morphologiques [START_REF] Richefeu | A new hybrid differential filter for motion detection[END_REF]. Il estime récursivement une amplitude de variation en chaque pixel. Cette méthode permet de détecter efficacement les mouvements de faible amplitude (objets lents, petits ou faiblement contrastés) avec une relativement bonne résistance au bruit grâce à la partie linéaire du filtre. En utilisant le paramètre α, un nombre réel compris entre 0 et

Implémentation

Afin de mener le calcul du gradient morphologique oublieux sur la rétine numérique, nous reprenons le principe du filtre précedent. En effet, le calcul des minima et maxima revient à calculer deux moyennes récursives en utilisant en plus un outil de comparaison. La constante α a la même dimension que dans le calcul de la moyenne récursive. De plus, ici on connaît le signe des différences entre l'image acquise X et les minima m et maxima M . 6 Moyenne Σ -∆

ε τ (I t ) δ τ (I t ) γ τ (I t ) m t M t Γ t

Principe

Dans le dernier algorithme [START_REF] Manzanera | A robust and computationally efficient motion detection algorithm based on Σ-∆ background estimation[END_REF], l'estimation du fond est vue comme la simulation d'une conversion numérique d'un signal analogique en utilisant la modulation Σ-∆ qui utilise seulement des comparaisons et des incrémentations/décrémentations élémentaires (Table 1(1). Comme la précision de la modulation Σ-∆ est limitée aux signaux dont la dérivée temporelle absolue est inférieure à un, l'erreur de modulation est proportionnelle au taux de variation du signal, correspondant ici à la mesure de vraisemblance en chaque pixel. Nous utilisons donc les différences absolues entre I t et M t comme premier estimateur: la différence ∆ t (Table 1(2)). On utilise aussi ce filtre pour calculer la variance temporelle du pixel, représentant sa mesure d'activité, utilisée pour décider si le pixel est plutôt mobile ou fixe. Donc le second estimateur V t (Tableau 1(3)) utilisé dans la méthode a la dimension d'un écart-type temporel calculé comme le filtre Σ-∆ de la séquence des différences ∆ t . Comme l'on s'intéresse aux pixels dont le taux de variation est significativement supérieur à leur activité temporelle, on applique le filtre Σ-∆ à N fois les différences non nulles. La Figure 9 montre les résultats du filtrage Σ-∆, toujours sur le même pixel. On peut remarquer comment le comportement non linéaire du filtre lui confère une plus grande robustesse. Finalement, la détection au niveau pixel est effectué par comparaison entre ∆ t et V t (Tableau 1(4)). La figure 10 montre le résultat de cette méthode pour une trame prise dans une scène de traffic urbain.

Implémentation

La première étape de la moyenne Σ-∆ consiste à effectuer la différence signée entre l'image acquise X et la moyenne M . Ensuite, comme pour le calcul des minima et maxima du gradient morphologique oublieux, on regarde le signe Nous utilisons le dialogue rétine-cortex pour adapter dynamiquement ce seuil en fonction d'une estimation de la distribution spatiale des différences. Pour réduire au maximum le flux d'information, nous utilisons une somme sur une image binaire, en appliquant récursivement l'ajustement suivant : après avoir seuillé l'image à la valeur τ , on calcule sur la rétine l'ensemble des points isolés (voir code rétine ci-dessous), puis on somme sur le cortex l'ensemble ainsi obtenu. Si la proportion de points isolés est inférieure à un taux fixé ρ (typiquement de l'ordre de 1%), on décrémente τ , sinon on l'incrémente. La justification de cet ajustement repose sur l'hypothèse qu'un point isolé est forcément du bruit et donc que le comptage de tels points fournit une estimation de ce bruit.

Un tel mécanisme peut également être utilisé dans le cas du troisième algorithme (moyenne Σ-∆) pour positionner automatiquement le paramètre N du calcul de V t . Enfin, il convient de préciser que de telles mesures globales peuvent êtres calculées efficacement sur la rétine en utilisant des opérateurs intégraux analogiques tels que des mesures de courant consommé [START_REF] Bernard | Output methods for an associative operation of programmable artificial retinas[END_REF]. Un fois la détection purement temporelle D t calculée, on la régularise en exploitant les corrélations spatiales grâce au filtrage suivant : on élimine les petites composantes connexes de D t en utilisant l'ouverture par reconstruction :

L (0) t = Rec Dt (ε B λ (D t )) (4)
où ε est l'érosion morphologique, B λ (l'élément structurant), est une boule de rayon λ et Rec Y (X) est la reconstruction géodésique de X dans Y . Finalement, on effectue une confirmation temporelle en calculant une autre reconstruction : Et enfin, ce dernier code permet d'effectuer la reconstruction du résultat P dans la détection courante (F ) :

L t = Rec L (0) t (L (0) t-1 ) (5)
// Reconstruction P = AND(P,F); for (i=0;i<N;i++) { P = OR(P,P.Ouest); P = OR(P,P.Est); P = OR(P,P.Sud); P = OR(P,P.Nord); P = AND(P,F); } Le modèle de rétine existant actuellement étant une machine parallèle purement synchrone, on utilise dans la reconstruction un nombre fixe d'itérations (N assez grand) au lieu d'une véritable relaxation. Néanmoins, la prochaine génération de rétine intègrera des opérateurs asynchrones qui permettront le calcul de la reconstruction géodésique en temps quasi constant.

Conclusion

Nous avons présenté trois exemples de calcul de détection du mouvement sur système de vision à base de rétine numérique. Nous avons mis en avant cette méthodologie particulière de conception d'algorithmes sur un tel système qui prend en compte la quantité limitée de mémoire disponible ainsi que les contraintes du calcul cellulaire. Une évaluation comparative de la qualité de détection des différents algorithmes reste à faire. Néanmoins nous avons noté expérimentalement les faits suivants :

• La méthode basée sur le calcul d'une moyenne récursive est moins coûteuse mais moins adaptative et moins précise en terme de localisation;

• La méthode basée sur les opérateurs de morphologie oublieux offre une plus grande capacité de détection mais est moins précise dans la localisation des objets détectés;

• La méthode basée sur l'estimation Σ-∆ est plus adaptative localement et temporellement, la localisation est plus précise mais elle est plus sensible aux modifications brutales du fond.

Figure 2 :

 2 Figure 2: Algorithme d'acquisition d'une image en niveau de gris.

Figure 3 :

 3 Figure 3: Calcul du fond par moyenne récursive. La courbe rouge correspond à la séquence originale, la courbe pointillée bleue au fond et la courbe en tiret vert à la différence au fond (α = 16).

InitialisationFigure 5 :

 5 Figure 5: Algorithme de calcul des opérateurs morphologiques oublieux.

Figure 6 :

 6 Figure 6: Filtrage morphologique oublieux. La courbe rouge correspond à la séquence originale, la courbe violette au maximum oublieux, la courbe pointillée bleue au minimum oublieux et la courbe en tiret vert le gradient morphologique oublieux (α = 8).

  En effet, X -M et m -X sont négatifs et l'on peut donc effectuer les opérations sans problèmes d'arrondis comme le montre la figure 8. La comparaison entre l'image acquise X, le maxima M AX et le minima M IN consiste à utiliser deux variables E et F ; E (resp. F ) vaut 1 si et seulement si X est inférieure à M IN (resp. X est supérieure à M AX). On utilise deux bascules D et G afin de savoir si le bit précédent était déjà inférieur ou supérieur à celui de X. Dans le code rétine suivant, A, B et C sont des variables temporaires : // Comparaison entre X, MAX et MIN D=0; E=0; F=0; G=0; for (i=7;i>=0;i--) { A = AND(NOT(X_i),MIN_i); B = AND(X_i,NOT(MIN_i)); C = AND(A,NOT(D)); E = OR(E,C); A = OR(A,B); D = OR(A,D); A = AND(X_i,NOT(MAX_i)); B = AND(NOT(X_i),MAX_i); C = AND(A,NOT(G)); F = OR(F,C); A = OR(A,B); G = OR(A,B); }

Figure 7 :Figure 8 :

 78 Figure 7: Application des opérateurs morphologique oublieux (en bas), comparés aux opérateurs morphologiques classiques (en haut), calculés à la trame t = 19 de la séquence "Tennis" (Berkeley). Pour comparaison, les opérateurs morphologiques classiques utilisent pour élément structurant le segment dans la dimension temporelle τ = [-8, 0]; et les opérateurs morphologiques oublieux α = 1/9. Les gradients sont affichés en mode vidéo inversé.

Figure 9 :

 9 Figure 9: Filtrage Σ-∆. La courbe rouge correspond à la séquence originale, la courbe pointillée bleue à la moyenne Σ-∆ et la courbe en tiret vert à la différence à la moyenne.

/

  / Detecteur de Points isoles P = NOT(P) P1 = AND(P.Nord, P.Est) P2 = P1.Sud P1 = AND(P1, P2.Ouest) P2 = AND(P.Nord, P.Sud) P2 = AND(P2.Est, P2.Ouest) P1 = AND(P1, P2) D = AND(NOT(P), P1)

  mémoire utilisée et vers des méthodes disposant d'une grande adaptabilité (globale (temporelle) pour les deux premières et locale (spatiotemporelle) pour la dernière méthode). Après avoir présenté brièvement la rétine programmable et le cadre général du problème de détection du mouvement, nous verrons trois exemples de ces techniques, leurs principes et leurs intérêts puis une étude de leur implémentation rétinienne. Nous présenterons d'abord une technique basée sur une moyenne récursive classique, puis une nouvelle méthode de différentiation fondée sur le filtrage morphologique oublieux, une combinaison d'opérations linéaires et morphologiques. Nous verrons son principe et son implémentation sur système à base de rétine numérique. Nous présenterons ensuite une technique s'appuyant sur une analogie de la modulation Σ-∆ en électronique. Ce filtrage très simple dans son principe nous permet d'aborder des problèmes de traitement d'images de plus haut niveau. Enfin, nous donnerons quelques exemples de techniques d'adaptation automatique d'un niveau de seuillage et de traitements spatiaux et nous concluerons sur les perspectives qu'apportent un tel système et ses performances.

2 Algorithmique de la rétine numérique

	V
	V s
	t
	Figure 1: Acquisition d'un niveau de gris par la lecture différents instants au cours de la décharge de la photodiode. à
	Comme nous l'avons vu précédemment les rétines sont adaptées à des applications portables, mobiles ou au-tonomes mais la petite capacité des éléments de cal-cul, et en particulier leur faible mémoire (typiquement quelques dizaines de bits par processseur cellulaire) con-stitue un réel défi pour la perception du mouvement. Cependant, parce que dans un système à base de rétine l'algorithmique débute dès l'acquisition et se poursuit dans la manière de coder une donnée et de l'extraire, c'est toute la méthodologie de conception des algorithmes qui se trouve transformée, en particulier par l'adaptation très fine à la dynamique tonale et temporelle rendue possible par la rétine programmable. La rétine est une grille de processeurs à topologie régulière dotée d'une architecture massivement parallèle. Chaque cellule de la rétine exécute à chaque cycle d'horloge la même instruction élémentaire : soit une opération booléenne entre deux de ses registres mémoire; soit une translation élementaire vers un de ses quatre voisins. L'acquisition d'une image numérique en n niveaux de gris s'effectue sur la rétine grâce à la lecture multiple du pho-torécepteur qui fournit (n -1) coupes binaires au cours du temps (Figure 1). En effet, le photorécepteur (photo-diode) se décharge d'autant plus vite (croissance hyper-bolique) qu'il rec ¸oit de photons. On effectue alors pour chaque pixel, un comptage sur log(n) bits du nombre de coupes où le pixel vaut 1. La technique de comptage la plus rapide sur le circuit con-siste à utiliser un compteur de Gray. En effet, le principe d'un tel compteur est que les codages de deux nombres consécutifs ne diffèrent que d'un seul bit et donc chaque incrément ne coûte qu'un seul XOR. Si l'image est représentée par ses coupes (C i ) 1≤i≤n-1 et si l'on note

  1, la dilatation (resp. érosion) temporelle oublieuse M t (resp. m t ) est définie comme le montre la Figure5. Comme pour la moyenne récursive classique, l'inverse de α a la dimension du temps. La sémantique de M

t (x) (resp. m t (x)) est alors la valeur maximale (resp. minimale) (estimée) observée au pixel x pendant les 1/α dernières trames. Pour α tendant vers un, M t (resp. m t ) tend vers I t , et pour α tendant vers zéro, M t (resp. m t ) tend vers la valeur maximale (resp. minimale) observée durant toute la séquence.

  MORPHO pour gradient temporel morphologique oublieux et Σ-∆ pour la moyenne Σ-∆) présentés en terme de temps de calcul (nombre d'instructions (INS)) et d'utilisation de l'espace mémoire (nombre de registres utilisés (REG)). Ainsi pour une rétine numérique cadencée à 1 Mhz (fréquence d'instructions), le temps de calcul par trame nécessaire pour toute la détection est inférieur à 0,45 ms pour la moyenne récursive et inférieur à 0,65 ms pour le gradient temporel morphologique oublieux et pour la moyenne Σ-∆.

		ALGORITHME Acquisition
	INS		262
	REG		9
		MOYREC	MORPHO	Σ-∆
	INS	88	303	312
	REG 18	36	28
		SEUIL GLOBAL
	INS REG	8 9
		FILTRAGE SPATIAL
	INS		66
	REG		5
	TOTAL MOYREC MORPHO Σ-∆ 424 639 640 REG 23 INS 41 32
	Figure 11: Tableau comparatif des trois algorithmes (MOYREC pour moyenne récursive,
	L t représente l'étiquette finale, le résultat de la détection. Sa sémantique est : les objets "plus gros que" λ qui appa-raissent sur deux trames consécutives. Pour effectuer une érosion par la boule élémentaire 4-connexe sur la rétine, il suffit de prendre la conjonction des quatre voisins de chaque pixel P :
	// Erosion 8-connexe
	P = AND(P,P.Nord)
	P = AND(P,P.Est)
	P = AND(P,P.Sud
	P = AND(P,P.Ouest)

vironnements dynamiques est une tâche bien plus ardue, pénalisant la détection des mouvements de faible amplitude (des objets petits ou très lents ou très peu contrastés). La troisième approche est basée sur le calcul de la vitesse locale apparente (flot optique)[START_REF] Beauchemin | The computation of optical flow[END_REF] qui est utilisé en entrée de la segmentation spatiale[START_REF] Ranchin | Moving objects segmentation using optical flow estimation[END_REF]. Ces méthodes engendrent des informations riches mais sont généralement plus complexes à calculer et sont aussi plus sensibles à la fiabilité des résultats issus du flot optique. Ainsi un compromis doit-être trouvé entre le calcul du flot optique et la précision de la segmentation. Enfin, plus récemment, des filtres morphologiques ont été employés[START_REF] Decenciere Ferrandiere | Application of the morphological geodesic reconstruction to image sequence analysis[END_REF] [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF] [1] afin d'analyser des séquences vidéos. En utilisant un élément structurant spatiotemporel, une amplitude de variation locale peut-être calculée comme mesure de vraisemblance du mouvement. Une telle mesure peut être utile pour détecter de petites amplitudes de mouvement, mais comme elle est sensible au bruit spatial, elle est souvent intégrée sur des régions entière en utilisant des opérateurs connexes. Nous allons maintenant présenter les trois méthodes de détection du mouvement (moyenne récursive, gradient morphologique oublieux et moyenne Σ-∆) que nous avons implantées sur la rétine numérique.

Actuellement, nous étudions la possibilité de combiner ces différentes méthodes afin de tirer parti au mieux des possibilités de chaque filtre. Un autre exemple d'implémentation d'algorithmes se trouve en [15] où l'on montre comment le calcul de points d'intérêt peut être mené de fac ¸on optimale sur de telles architectures.

Initialisation pour chaque pixel x:

pour chaque trame t pour chaque pixel x:

(3) (4)

Adaptation et traitements spatiaux

Pour les deux premiers algorithmes (moyenne récursive et gradient morphologique oublieux), nous utilisons un seuil global (i.e. le même pour tous les pixels). La valeur de ce seuil doit nécessairement être ajustée dynamiquement si l'on souhaite que la détection soit adaptative.


Julien Richefeu

Antoine Manzanera

A NEW HYBRID DIFFERENTIAL FILTER FOR MOTION DETECTION

Keywords: hybrid filter, temporal morphology, motion detection

A new operator to compute time differentiation in an image sequence is presented. It is founded on hybrid filters combining morphological and linear recursive operations. It estimates recursively the amplitude of time-variation within a certain interval. It combines the change detection capability of the temporal morphological gradient, and the (exponential) smoothing effect of the linear recursive average. It is particularly suited to small and low amplitude motion. We show how to use this filter within an adaptive motion detection algorithm.

Introduction and Preliminaries

In some applications such as video databases, video compression, security monitoring or medical imaging, motion information is usually a more significant cue than color or texture. More specifically, the need to detect interesting moving objects in the scene is a fundamental low-level step in many vision systems. In this paper, we will concentrate on the case of video surveillance systems using a stationary camera.

The challenge of motion detection lies in the ability of performing an accurate segmentation of the moving objects independently on their size, velocity and contrast with respect to the background of the scene. A large set of motion detection algorithms has already been proposed in the literature. We can file them into four main categories according to the type of inter-frame computations, i.e. to the way the time differentiation is performed.

The first one is based on temporal gradient: a motion likelihood index is measured by the instantaneous change in the image intensity computed by differentiation of consecutive frames [START_REF] Bouthemy | Recovery of moving objects masks in an image sequence using local spatiotemporal contextual information[END_REF]. These methods are naturally adaptive to changing environments, but are also dependent on the velocity and size of moving objects. This drawback can be minimized using a multiple bank of spatiotemporal filters which is done at the price of an increased complexity.

The second category are the background subtraction techniques [START_REF] Toyoma | Wallflower: principles and practice of background maintenance[END_REF] [4] [START_REF] Piccardi | Background subtraction techniques: a review[END_REF], that use a reference image (background), representing the stationary elements in the scene. Here the motion likelihood measure is the difference between the current frame and the background. These methods are less dependent on the velocity and size of the objects. Nevertheless, the adaptation to dynamic environment is a much more difficult task, which can penalize the detection of small amplitude motion (very slow or low contrast objects).

The third type of approach is based on the computation of the local apparent velocity (optical flow) [START_REF] Beauchemin | The computation of optical flow[END_REF] which is used as input of a spatial segmentation [START_REF] Ranchin | Moving objects segmentation using optical flow estimation[END_REF]. This method provides valuable information but it is in general more computationally complex and it is also sensitive to the reliability of the optical flow. Thus a trade-off has to be found between the smoothness of the flow field and the accuracy of segmentation.

More recently, morphological filters have been employed [START_REF] Decenciere Ferrandiere | Application of the morphological geodesic reconstruction to image sequence analysis[END_REF] [10] [START_REF] Agnus | Spatio-temporal segmentation using morphological tools[END_REF] for video sequence analysis. By using spatiotemporal structuring elements, a local amplitude of variation can be computed as motion likelihood index. Such measure can be useful to detect small amplitude motion, but as it is sensitive to outliers, it is usually integrated over regions using connected operators.

We propose in this paper a new differential operator based on a hybrid filter, combining morphological and linear operations. It computes a pixel-wise amplitude of time-variation over a recursively defined "temporal window". This method is designed to address the problem of small objects and slow motion while providing a certain noise immunity thanks to its linear part.

We first present the forgetting morphological temporal gradient in Section 2. Then we show how to use the output of this filter within a motion detection algorithm in Section 3. Results are presented and discussed in the same section.

The forgetting morphological temporal gradient

In this section, we introduce the forgetting morphological temporal gradient and show its interests for motion detection systems.

Considering an image sequence I t (x) where t is a time index and x a (bidimensional) space index, morphological temporal filters are defined using temporal structuring element

τ = [t 1 , t 2 ]. The temporal erosion (resp. dilation) is defined by ε τ (I t )(x) = min z∈τ {I t+z (x)} (resp. δ τ (I t )(x) = max z∈τ {I t+z (x)}).
The temporal (morphological) gradient γ is then defined by γ τ (I t ) = δ τ (I t )ε τ (I t ).

τ represents the temporal interval of interaction, which can be causal (e.g. [-3, 0]), anti-causal (e.g. [0, +5]) or both (e.g. [-1, +1]). Thus the temporal gradient corresponds to the amplitude of variation within this interval.

The use of this operator suffers from two major drawbacks: (1) it implies the use of a buffer with size corresponding to the diameter of the structuring element, which can be very memory consuming; (2) it is very sensitive to sudden large variations (like impulse noise or slight oscillations of the sensor).

To cope with these two problems, we use hybrid filters, which can be viewed as a recursive estimation of the values of the temporal erosion and dilation. Using parameter α, which is a real number between 0 and 1, the forgetting temporal dilation M t (resp. erosion m t ) is defined as shown in Figure 1. As in the classical running average (or exponential smoothing) defined by A t (x) = αI t (x) + (1 -α)A t-1 (x), the inverse of α has the dimension of time. The semantics of M t (x) (resp. m t (x)) is then the (estimated) maximal (resp. minimal) value observed at pixel x within the 1/α last frames. So as α tends to unity, M t (resp. m t ) tends to I t , and as α tends to zero, M t (resp. m t ) tends to the maximal (resp. minimal) value observed during the whole sequence. The use of the term "forgetting" is justified by the fact that these operators attach more importance to the near past than to the far past. Γ t , the forgetting morphological temporal gradient, is used in the following as a differentiation filter because of its interesting properties: (1) it has the dimension of an amplitude of time-variation, so it is able to integrate motion over a long period depending on 1/α, and then to detect small or slow moving objects; (2) it is less sensitive to impulse noise because of its forgetting term, corresponding to the exponentially decreasing weights attached to the past values; (3) it only requires the use of two buffers to compute the forgetting erosion and dilation. Figure 2 displays the forgetting morphological operations compared with their morphological counterparts.

The forgetting morphological gradient thus represents a relevant motion likelihood measure. We show on the next section how it can be used within a complete moving objects detection algorithm.

Motion detection algorithm

The filter presented above makes possible a good level of detection for motion whose amplitude is below the spatiotemporal discretization. Nevertheless, the forgetting term α needs to be adapted to the velocity of the moving object. As there are several objects with different sizes and velocities within the scene, it is necessary to adjust locally the value of α to the observations.

In addition to this, we need a decision rule to discriminate the moving objects from the background. Because the scene is constantly evolving -typically under illumination or weather condition changes -the decision criterion has to be temporally adaptive. Furthermore the temporal variation, possibly due to moving objects, but also to noise or irrelevant motion, is not uniformly distributed through the scene. So the decision must be locally differentiated.

In our algorithm, we compute, from the forgetting morphological temporal filter, a local estimation of the spatiotemporal activity. This estimation is used both for deciding the pixel-level motion label and for adjusting the value of the forgetting term α.

Following [START_REF] Manzanera | A robust and computationally efficient motion detection algorithm based on Σ-∆ background estimation[END_REF], we use the Σ-∆ filter to compute a second order statistics on the sequence: The Σ-∆ filter S t of a time series X t is a recursive approximation of the median defined by: S

t = S t-1 + 1 if S t < X t and S t = S t-1 -1 if S t > X t .
Now, what we compute exactly is V t , the Σ-∆ filter of N times the nonzero values of the forgetting morphological temporal filter. Then, the local estimation of the spatiotemporal activity is defined by Θ t = G σ (V t ), where G σ is the bidimensionnal Gaussian filter with standard deviation σ.

Θ t is used as the first step of decision for the moving label of each pixel: D t = 1 if Γ t > Θ t and D t = 0 elsewhere.

Θ t is also used to update the forgetting term for every pixel. Inspired by [START_REF] Pic | Active background estimation: Computing a pixelwise learning rate from local confidence and global correlation values[END_REF] who used a locally adaptive learning rate for recursive background estimation, we compute the local α for pixel x by a similar formula:

α t (x) = 2 -Θ t (x) 2 k 2
where Θ t is the complementary version of Θ t (e.g. Θ t = 255 -Θ t for images coded on 255 gray level) and k is a constant used to set the range value of α. α is then locally adjusted in such a way that the forgetting filters use long term memory in poor spatiotemporal activity areas, and short term memory in areas with high spatiotemporal activity. Thus, the detection of slow, small, or low-contrast moving objects will be enhanced, while large moving objects with high contrast will be better segmented. Figure 3 shows the motion detection algorithm steps applied on the Hamburg Taxi sequence. The last image in the figure represents D t , the detection result after removal of the smallest regions using an alternated filter by reconstruction with a ball of radius 1 [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: applications and efficient algorithms[END_REF]. It can be seen that the small amplitude motion like the pedestrian on the top left, as well as low-contrast moving object like the dark car on the bottom left are well detected while the high contrasted taxi at the center is better segmented. This is an effect of the adaptable memory of the forgetting filters, in order to fit their detecting ability to the amount of motion.

Conclusion

We have presented a new hybrid differential filter and shown how it can be used in a motion detection algorithm with a local adaptation of the forgetting terms. Like the recursive filters, the computation time and memory consumption does not depend on the size of the temporal window. At the present time, we are investigating more sophisticated spatial interaction, in order to improve the spatiotemporal adaptivity, and quantify precisely the validity range of the algorithm.

InitializationFigure 1 .

 1 Figure 1. The forgetting morphological temporal operators: Mt is the forgetting dilation, mt the forgetting erosion and Γt the forgetting morphological gradient.

Figure 2 .

 2 Figure 2. Application of the forgetting morphological operators (bottom), compared with the classical morphological temporal operators (top), computed on the frame t = 19 of the classical "Tennis" sequence (Berkeley). For comparison purposes the structuring element is τ = [-8, 0], and the forgetting term is α = 1/9. The gradients are displayed in reverse video mode. Note that we use the symmetrical gradient in order to treat dark and light objects the same way.

Figure 3 .

 3 Figure 3. Results for the Hamburg Taxi sequence (frame n. 20). The parameters used are N = 2 (number of deviations) for St, σ = 2.5 for the standard deviation of the Gaussian filter used to compute Θt, k = 140 as the constant used in the computation of αt.


A Manzanera

J C Richefeu

A robust and computationally efficient motion detection algorithm based on -¡ background estimation

This paper presents a new algorithm to detect moving objects within a scene acquired by a stationary camera. A simple recursive non linear operator, the £ -¤ filter, is used to estimate two orders of temporal statistics for every pixel of the image. The output data provide a scene characterization allowing a simple and efficient pixel-level change detection framework. For a more suitable detection, exploiting spatial correlation in these data is necessary. We use them as a multiple observation field in a Markov model, leading to a spatiotemporal regularization of the pixel-level solution. This method yields a good trade-off in terms of robustness and accuracy, with a minimal cost in memory and a low computational complexity.

Introduction

To detect the moving objects in an image sequence is a very important low-level task for many computer vision applications, such as video surveillance, traffic monitoring, sign language recognition. When the camera is stationary, a class of methods usually employed is background subtraction. The principle of these methods is to build a model of the static scene (i.e. without moving objects) called background, and then to compare every frame of the sequence to this background in order to discriminate the regions of unusual motion, called foreground (the moving objects).

Many algorithms have been developed for background subtraction: recent reviews and evaluations can be found in [START_REF] Lee | Background estimation for video surveillance[END_REF] [START_REF] Chalidabhongse | A perturbation method for evaluating background subtraction algorithms[END_REF] [3] [START_REF] Piccardi | Background subtraction techniques: a review[END_REF]. In this paper, we are more specifically interested in outdoor video surveillance systems with long autonomy. The difficulty in devising background subtraction algorithms in such context lies in the respect of several constraints:

¥ The system must keep working without human interaction for a long time, and then take into account gradual or sudden changes such as illumination variation or new static objects settling in the scene. This means that the background must be temporally adaptive. ¥ The system must be able to discard irrelevant motion such as waving bushes or flowing water. It should also be robust to slight oscillations of the camera. This means that there must be a local estimation for the confidence in the background value.

¥ The system must be real-time, compact and lowpower, so the algorithms must not use much resource, in terms of computing power and memory.

The two first conditions imply that statistical measures on the temporal activity must be locally available in every pixel, and constantly updated. This excludes any basic approach like using a single model such as the previous frame or a temporal average for the background, and global thresholding for decision.

Some background estimation methods are based on the analysis of the histogram of the values taken by each pixel within a fixed number ¦ of past frames. The mean, the median or the mode of the histogram can be chosen to set the background value, and the foreground can be discriminated by comparing the difference between the current frame and the background with the histogram variance. More sophisticated techniques are also based on the ¦ past frames history: linear prediction [START_REF] Toyoma | Wallflower: Principles and Practice of Background Maintenance[END_REF], kernel density estimation [START_REF] Elgammal | Non-parametric Model for Background Subtraction[END_REF] [10], or principal component analysis [START_REF] Oliver | A bayesian computer vision system for modeling human interactions[END_REF]. These methods require a great amount of memory, since ¦ needs to be large (usually more than 50) for robustness purposes. So they are not compatible with our third condition.

Much more attractive for our requirements are the recursive methods, that do not keep in memory a histogram for each pixel, but rather a fixed number of estimates computed recursively. These estimates can be the mean and variance of a Gaussian distribution [START_REF] Wren | Pfinder: Real-time tracking of the human body[END_REF], or different states of the background (e.g. its values and temporal derivatives) estimated by predictive (e.g. Kalman) filter [START_REF] Karmann | Time-Varying Image Processing and Moving Object Recognition, chapter Moving Object Recognition Using an Adaptive Background Memory[END_REF]. But it is difficult to get robust estimates of the background with linear recursive framework, unless a multi-modal distribution (e.g. multiple Gaussian [START_REF] Stauffer | Learning patterns of activity using real-time tracking[END_REF]) is explicitly used, which is done at the price of an increasing complexity and memory requirement. Furthermore, these methods rely on parameters such as the learning rates used in the recursive linear filters, setting the relative weights of the background states and the new observations, whose tuning can be tricky, which makes difficult the fulfillment of the first condition stated above.

A recursive approximation of the temporal median was proposed in [START_REF] Mcfarlane | Segmentation and tracking of piglets in images[END_REF] to compute the background. The interest of this method lies in the robustness provided by the non linearity compared to the linear recursive average, and in the very low computational cost. In this paper, we investigate some nice properties of this method, introducing the notion of £ -¤ filtering, and using it to obtain a non-parametric motion detection. In Section 2, we use the £ -¤ filter to compute two orders of temporal statistics for each pixel of the sequence, providing a multiple observation field, and a pixel-level decision framework. Then, in Section 3, we exploit the spatial correlation in these data using a Markov based spatial regularization algorithm. High level processing and feedback are then discussed in Section 4, where some results are displayed.

Temporal processing

Our first background estimate, shown on Table 1(1), is the same as [START_REF] Mcfarlane | Segmentation and tracking of piglets in images[END_REF], where § ©¨is the input sequence, and ¨the estimated background value. As noticed in [START_REF] Mcfarlane | Segmentation and tracking of piglets in images[END_REF], if § ¨is a discrete random signal, the most probable values of ¨lie in an interval such that there are as many indices ! " $# such that § &% " ' , as indices ! (" )# such that § &% 0 1 . So ¨is an approximation of the median of § ¨. But this filter has other interesting properties, relative to the change detection in time-varying signals. Indeed, we interpret this background estimation as the simulation of a digital conversion of a time-varying analog signal using £ -¤ modulation (A/D conversion using only comparison and elementary increment/decrement, hence the name £ -¤ filter).

As the precision of the £ -¤ modulation is limited to signals with absolute time-derivative less than unity, the modulation error is proportional to the variation rate of the signal, corresponding here to a motion likelihood measure of the pixels. We then use the absolute difference between § &ä nd 2¨as the first observation field: the difference ¤ 3¨(Table 1(2)).

Unlike [START_REF] Mcfarlane | Segmentation and tracking of piglets in images[END_REF], we also use this filter to compute the timevariance of the pixels, representing their motion activity measure, used to decide whether the pixel is more likely "moving" or "stationary". Then the second observation field 4 ¨(Table 1(3)) used in our method has the dimension of a temporal standard deviation. It is computed as a £ -¤ filter of the difference sequence ¤ 3¨. This provides a mea-sure of temporal activity of the pixels. As we are interested in pixels whose variation rate is significantly over its temporal activity, we apply the £ -¤ filter to the sequence of 5 times the non-zero differences.

Finally, the pixel-level detection is simply performed by comparing ¤ 3¨and 4 6¨(Table 1(

For each frame t for each pixel 7 : g P hA 7 6C ED (i 8 @P hA 7 6C pe `F P aA 7 6C i

Initialization for each pixel 7 : q r9 BA 7 6C sD g 9 BA 7 6C

For each frame t for each pixel 7 such that g P aA 7 IC utD 2v : if q rP RQ TS UA 7 6C wV fx 'y g P hA 7 6C , q rP aA 7 6C ED q rP RQ S A 7 6C 6b 2c if q rP RQ TS UA 7 6C wd fx 'y g P hA 7 6C , q rP aA 7 6C ED q rP RQ S A 7 6C e fc

For each frame t for each pixel 7 : if g P hA 7 6C wV q rP aA 7 6C then P aA 7 6C sD 2v

else P A 7 6C ED c (4)
Table 1: The £ -¤ background estimation: (1) Computation of the £ -¤ mean. (2) Computation of the difference between the image and the £ -¤ mean (motion likelihood measure). ( 3) Computation of the £ -¤ variance defined as the £ -¤ mean of 5 times the non-zero differences. (4) Computation of the motion label by comparison between the difference and the variance.

Figures 2 to 4 display an example of the evolution over time of the different values computed as above, for three pixels extracted from a country scene similar to the image shown on Figure 1, for a 1000 frames sequence. The red solid line represents the input image § ¨. The blue dashed line corresponds to the £ -¤ mean ¨. The green dashed line represents the difference ¤ ¨. Finally, the purple dotted line is the £ -¤ variance 4 ¨(using 5 ). The detection field 3¨is not represented explicitly, but corresponds to the Boolean indicator of the condition "the green dashed line is over the purple dotted one". The pixel used in Figure 2 is a pixel in a still zone, with flat temporal activity, such as a remote area of the static background (in this example, a sky lightly covered with slowly moving clouds). For such pixels, the high frequency variation corresponds to temporal noise due to the acquisition and digitization processes. The low frequency variations are due to illumination changes or slow motion of low contrast objects. The pixel used in Figure 3 is a pixel in a motion area, such as tracks or corridors (in this example, a country road with 4 vehicles passing away). In that case, the moving objects give rise to sharp changes that are not taken into account by the £ -¤ mean, and then the difference field shows a peak. Such peaks are discriminated thanks to the comparison with the £ -¤ variance.

Clutter area Motion area

The pixel used in Figure 4 is a pixel in a clutter area, i.e. a zone of physical changes due to intrinsic nature of the scene rather than moving objects. Examples of such areas are: trees moving with the wind, river, or crowd in a urban scene (in our example, high grass in the foreground of the scene). In that case, the difference field shows a repetition of peaks, and if these peaks are close enough from each other with respect to the delay induced by £ -¤ modulation, then they will be taken into account in the £ -¤ variance, in such a way that the difference will remain less that the variance. Figure 6 (1)-( 4) displays the result of this method for one frame of an urban traffic sequence. The four images represent respectively § ¨, ¨, 4 ¨and ¨. It can be seen that the discrimination of "moving" pixels corresponds to the detection of temporally salient pixels with respect to the temporal activity. This allows to discard irrelevant (clutter) motion, but also to be less sensitive to sensor oscillations, as it is shown on Figure 6 (5) and ( 6): A uniform random oscilla-tion of pixels has been simulated on the same sequence. In this case, 2¨converges to an approximation of the spatiotemporal median, and then 4 ¨(Fig. 6(5)) emphasizes the regions of high contrast, thus increasing the local threshold in these regions, and avoiding the detection of the whole scene contour in ¨(Fig. 6(6)).

The £ -¤ background estimation provides a simple and efficient method to detect the significantly changing pixels in a static scene, with respect to a time constant depending on the number of gray levels of the images. Nevertheless it is a pure temporal processing, which can only provide pixel-level detection. In the following section, we use a Markovian framework to improve the motion detection by exploiting the spatial correlation in the data.

Spatial regularization

We follow the model that has been used for real-time implementation in different architectures in [START_REF] Caplier | Mrf based motion detection algorithm image processing board implementation[END_REF] and in [START_REF] Lacassagne | Motion detection, labeling, data association and tracking in real-time on risc computer[END_REF]. This Markov model is based on the estimation of a binary (background/foreground) motion field given an observation field , by maximizing a Bayesian maximum a posteriori criterion, i.e. given a realization of the observation field , finding the realization of the motion label field that maximizes the conditional probability d e sf g 3 6h . Under the hypothesis that is a Markov field, and a probabilistic model linking and , this corresponds to finding the motion field that minimizes the global energy function defined over the set of pixels i as follows: typically 4). This algorithms is known to be very sensitive to the quality of the initial value of the estimated motion field. [START_REF] Caplier | Mrf based motion detection algorithm image processing board implementation[END_REF] and [START_REF] Lacassagne | Motion detection, labeling, data association and tracking in real-time on risc computer[END_REF] use a threshold of the observation (i.e. the absolute difference between two consecutive frames) as the initial estimation of .

In our algorithm, we use the same model, with the following exceptions (referring to the variables used in Table 1):

¥ for the observation , we use ¤ the difference between the background and the current frame.

¥ we use the £ -¤ variance 4 as a second observation field, to estimate locally the dispersion factor:

¥¤ ¦ h
is used instead of ¢ for weighting the relative importance of j u with respect to j q

.

¥ the initialization of the Markovian relaxation corresponds to the pixel-level detection .

What are the advantages of this algorithm compared to the original Markovian model ?

First, the difference with the £ -¤ background is more robust than the frame to frame difference, because it combines information over a large period of time instead of two frames. It is much less sensitive to the aperture problem, which makes difficult the detection of large homogeneous zones in motion. It is also less dependent on the velocity of the objects.

Next, for the same reasons, is in general much better to initialize the relaxation than a binarized frame difference, because it is closer to the expected solution. It must be emphasized that, for the ICM algorithm, once the initialization computed, the other parameters of the model are not critical, and have shown good behaviors on lots of different sequences.

Finally, the dispersion parameter is computed locally, then no global computation is needed at each frame. This allows the computation of the whole algorithm using only local memory sharing, thus permitting a massive spatial parallelism.

To increase the confidence in the detected objects, a higher level of processing, is needed, involving regional and global computation. This is discussed in the following section.

Region level processing

A precise description of the region level processing is not in the scope of this paper. We will only present its principles and discuss the feedback that can be done on the local computation in order to increase the robustness of the background estimation.

The output of the Markovian regularization is a lowlevel estimation of the foreground, consisting of the temporal salient pixels presenting spatial correlation. In order to enhance the quality of detection and lower the false alarm rate, some higher level processing is needed, using regional and global computations.

The pixels are grouped into regions representing objects. It is usually done by connected components labeling followed by fusion. The resulting objects then undergo morphological filtering, which can reject some objects under size or shape criteria. Kinematic filtering can also be employed in order to discriminate the objects whose motion is consistent with regard to the application (e.g. car, pedestrian,...).

In addition to the diminution of the false alarm rate, the interest of the global level processing is to allow a frame rate feedback on the low-level detection. One of the most straightforward example is the adaptation to a sudden change of background: if a global confidence index in the background (e.g. the relative surface area occupied by the foreground) decreases beneath a certain level, the decision can be made to re-initialize the background, in order to lower the re-adaptation time.

Another useful feedback is to attach a confidence index § to each filtered object, from 1 (lowest confidence) to ¨(absolute confidence). Those indices are then used as a period of update for the £ -¤ estimation: if the pixel s belongs to the foreground with confidence are updated only every § xs h frames. This enhances the quality of detection by increasing the robustness of the nonlinear filter, and avoiding the objects stopping temporarily (such as cars at red light) to enter in the background too quickly (thus generating "ghost" when they move again).

Figures 6(9) and [START_REF] Mittal | Motion-based background subtraction using adaptive kernel density estimation[END_REF] show an application of the relevance feedback, using a uniform confidence index of 6 over a mask which is simply a morphological filtering of the spatial regularization output.

Conclusions

We have presented a new algorithm allowing a robust and accurate detection of moving objects for a small cost in memory consumption and computational complexity. We have emphasized the nice properties of the £ -¤ filter for the detection of salient features in time-varying signal, showing that the interest of such filter goes well beyond its temporal median convergence property.

We have adapted a classical Markovian model to perform the spatial regularization that is needed to eliminate irrelevant salient pixels and aggregate the relevant ones in significant region. The multiple observation field produced by the £ -¤ estimation has shown relevant as input of the Markovian relaxation algorithm.

Because it only relies on pixel-wise or spatially limited interactions, the whole low-level processing (temporal processing and spatial regularization) is suited to a massively parallel implementation. We are at the present time implementing the algorithm on a programmable artificial retina [START_REF] Komuro | A digital vision chip specialized for high-speed target tracking[END_REF], which is a fine-grained parallel machine with optical input. The algorithm is indeed well adapted to the architecture, which consists in a mesh of tiny processors with limited memory and computation power. We have already implemented an alternate version of the algorithm (same temporal processing followed by a spatiotemporal morphological filtering) with the following performances, for a 200x200 retina array running at 25 Mhz, using 8 bits per pixels: 2.25 ms per frame, of which only 0.75 ms for the sole computation, and the rest for the acquisition.

The present limitation of our approach lies in the adaptation capability to certain complex scenes. It is the case for very low motion, that is likely to be taken as a background characteristics. A straightforward way to address this problem is to downsample the updating rate of the estimates, i.e. to use only every ¡ th frame to compute ¨and every © th frame to compute 4 ¨. See Figure 5 for an example of result (in this sequence, the scene is always full with persons moving then stopping for a while, and the camera is slightly oscillating). But this adds two learning rate parameters that can be critical, because reducing the updating rate increases the sensitivity to oscillations and the adaptation to a new scene. Another limitation case is the wide amplitude periodical motion (e.g sea surge), that will be classified as foreground if the period is too long. We are investigating the possibilities to go beyond these limits by combining different variance models using different sampling periods in the £ -¤ estimation, in order to get a richer quantitative estimation of the motion activity.

In future works, we will also focus on more sophisticated high-level filtering, in order to increase the robustness of the background estimation, and then of the whole detection.

Figure 1 :

 1 Figure 1: Example of observed scene, with 3 particular pixels from 3 different areas.

  Signal" "Still/Difference" "Still/Sigma-Delta_mean" "Still/Sigma-Delta_variance"

Figure 2 :

 2 Figure 2: Temporal variation of a pixel value in a still area.

  Difference" "Motion/Sigma-Delta_mean" "Motion/Sigma-Delta_variance"

Figure 3 :

 3 Figure 3: Temporal variation of a pixel value in a motion area.

Figure 4 :

 4 Figure 4: Temporal variation of a pixel value in a clutter area.

Figure 6 :

 6 Figure 6: Result of the proposed algorithm on a traffic sequence. (1) § ©¨(2) (3) 4 ¨(displayed with normalized histogram). (4) 3¨(N=2). (5) 4 ¨with simulated oscillations of the camera (Normalized histogram). (6) ¨for the oscillating camera (N=2). (7) Detection after Markovian regularization (5 iterations). (8) idem for the oscillating camera. (9) 4 ¨using relevance feedback (applying the same transformation as image (3). (10) Detection using relevance feedback on the pixel level processing and Markovian regularization.

  4)).

	Initialization
	for each pixel 7 : 8 @9 BA 7 6C ED GF 9 HA 7 IC
	For each frame t
	for each pixel 7 : if 8 @P RQ TS UA 7 6C WV XF P YA 7 6C , if 8 @P RQ TS UA 7 6C Wd XF P YA 7 6C ,	8 `P aA 7 6C ED 8 `P aA 7 6C ED	8 @P RQ TS UA 7 6C 6b 2c 8 @P RQ TS UA 7 6C e fc

  on the nature of the neighborhood. We use a uniform 6-connected spatiotemporal topology with 3 different values p p for the 4 spatial neighbors, T for the past neighbor, and p for the future neighbor.

												Rs use the absolute difference between two consecutive frames h ah £ if xs h has the foreground value. [1] and [7] as the observation field. They use a constant value for £ (20), and estimate ¢ by computing the spatial variance of the observation. j
												The minimization of the global energy the deterministic relaxation called iterated conditional mode is realized by (ICM): all the pixels are sequentially updated, and each pixel s is given the label xs j q j u h corresponding to the smallest local energy xs h ah vt xs h © a xs h ah . Usually, instead of a true relaxation, a limited number of scans is performed (
		j		j rq		j vu
		lk m n po I j q	Rs	h hh Wt	Rs	h U	Rs	h hh w ,
	j q	with j u and Rs	h U	xs	h hh v Rs h hh	k n {z E| }m w~4 T y g Rs	Rs h 2	h © a	xs	h hh , h hh	.
	Rs vide spatiotemporal regularity in the motion field. It is h hh is called model energy and is designed to pro-based on the Markovian modeling of as a Gibbs field, where is the set of neighbors of the pixel s , and the po-tential functions 4 p Rs h © a h hh equals m y if Rs h h , and t m y if Rs h h . The m y are positive constants, whose values depend j vu
	Rs to ensure a certain level of attachment to the input data, h © a Rs h ah is called fitness energy and is designed i.e. the observation . This term comes from the condi-tional probability of the observation field , with respect to the motion field , assuming that xs h Rs h hh t ¡ 6 ¢ h , with ¡ 6 a¢ h a centered Gaussian noise of vari-ance ¢ , xs h hh if xs h has the background value, and