Silvia Zuffi

Javier Romero

Cordelia Schmid

Michael J Black

Estimating Human Pose with Flowing Puppets

Figure 1. Flowing puppets. (a) Frame with a hypothesized human "puppet" model; (b) Dense flow between frame (a) and its neighboring frames; (c) The flow of the puppet is approximated by a part-based affine motion model; (d) The prediction of the puppet from (a) into the adjacent frames using the estimated flow.

Introduction

We address the problem of estimating the 2D pose of a person in a monocular video sequence captured under uncontrolled conditions, without manual initialization. In a single frame, pose estimation is challenging and current methods tend to do poorly at estimating the pose of the limbs. Arms and hands are relatively small and can be difficult to localize due to occlusion, motion blur, accidental alignment, ambiguities, poor contrast, etc. Incorporating information from multiple frames may alleviate some of these problems but the question is how best to achieve this? Previous approaches use image evidence in individual frames and then try to infer a coherent sequence of poses by imposing priors that encode smooth motion over time. Such approaches can work well for tracking where an initial pose is given but, as we describe below, are difficult to use for the general pose inference problem.

Like previous work we want to exploit the inherent consistency of appearance and motion over time. However, we take a novel approach that does not rely on human motion priors. Instead, we exploit optical flow in three ways: 1) to exploit image evidence from adjacent frames, 2) to propagate information over time, and 3) to provide richer cues for pose estimation.

Our approach is enabled by recent advances in methods for dense optical flow computation and by a recently introduced 2D model of articulated human body shape, the Deformable Structures model (DS) [24]. The availability of accurate estimates of dense optical flow allows us to consider the flow as an observation, while the articulated model of 2D body shape provides a tool for modeling the regions of motion of a moving person. The question is: How can optical flow be incorporated to make the pose inference problem simpler and more accurate?

Consider the problem of estimating body pose in Fig. 1. Assume we have a hypothesis for the body at frame t (Fig. 1(a), red). In any given frame, the image evidence may be ambiguous and we would like to combine evidence from multiple frames to more robustly infer pose. Due to the complexity of human pose, we perform inference using a distribution of "particles" at each frame, where each particle represents the pose of the body. If we are lucky and have particles at frame t and t + 1 that are both correct, then the poses in each frame explain the image evidence and the change in pose between frames is consistent with the flow. This should increase our confidence in the solution.

Unfortunately such an approach is not practical. Estimating the pose of the body in two frames simultaneously effectively doubles the size of the state space which, for articulated body models, is already high. Alternatively, if we independently estimate the pose in both frames then, given the high-dimensional space and a small set of particles, we will have to be extremely lucky to have two poses that are consistent with the image evidence in both frames and the optical flow. We need a different approach.

Our first solution is to estimate the pose of the body only at one frame (keeping the dimensionality under control) and to use the optical flow to check how good this solution is in neighboring frames. We refer to the body model as a "puppet" because it can be "puppeteered" by the optical flow. Given a pose at frame t we use the computed dense optical flow (Fig. 1(b)) to predict how the puppet should move into the next frame, forwards and backwards in time. The puppet flow (Fig. 1(c)), estimated from the dense optical flow, provides the prediction of the puppet in the next fame (Fig. 1(d)). We then extend our image likelihood to take into account evidence from the neighboring frames. The advantage is that inference takes place for a single puppet at a time but we are able to incorporate information from multiple frames.

Image evidence is computed in each frame using a DS model [24]. We describe upper body pose estimation but the method should be applicable to full body pose as well. This model captures the rough shape of a person and how the shape of the body parts deform with pose. We go beyond previous work to train a multi-scale image likelihood that captures the statistics of image gradients along the contour. By learning the model at multiple scales, we capture information about how real people deviate from the model.

Our second use of optical flow is in search. Our optimization uses a particle-based stochastic search method [START_REF] Ivekovič | Human body pose estimation with particle swarm optimisation[END_REF]. We initialize particles on each frame of the video sequence using a state-of-the-art single-frame pose estimation method [23]. We take the most likely particles in a given frame and use the puppet flow to predict their poses in adjacent frames. This enriches the particle set at neighboring frames. Inference always happens in a single frame but the two methods above serve to incorporate information from adjacent frames. We generate additional pose proposals that incorporate information about the possible location of hands based on image and flow evidence; this is our third use of flow.

We compare our method with [23] and [20] on the Video-Pose2.0 dataset [20]. VideoPose2.0 is a complex and challenging benchmark for pose estimation methods in video sequences that includes very difficult sequences where the appearance of the people can be easily confounded with the background. Until now Sapp et al. [20] had the best results on this dataset but they rely on knowing the correct scale of the person and a bounding box used to rescale and crop the video frames. Here we remove these restrictions while obtaining more accurate estimates of the wrists.

In summary our work proposes a new way of integrating information over time to do human pose estimation in video. The key idea is to use the optical flow field to define "puppets" that "flow" from one time to the next, allowing us to integrate image evidence from multiple frames in a principled and effective way and to propagate good solutions in time. A good pose is one that is good in multiple frames and agrees with the optical flow.

Background and Related Work

Human pose tracking in video. Much of the early work in the field addresses the tracking problem. It assumes either that there is a known template [4] or the first pose is known [15]. The tracking literature is extensive and beyond our scope.

Human pose estimation in still images. There is a similarly large literature on 2D human pose estimation in static images. Again a full review here is not feasible but the most relevant recent work is based on pictorial structures (PS) models [START_REF] Andriluka | Pictorial structures revisited: People detection and articulated pose estimation[END_REF]7,[START_REF] Felzenszwalb | Pictorial structures for object recognition[END_REF]18,[START_REF] Sapp | Adaptive pose priors for pictorial structures[END_REF]. Such models are widely used but still have trouble accurately localizing human arms. Here we use the Flexible Mixtures of Parts (FMP) model [23] on each still video frame to provide an initialization. FMP is one of the most adopted methods for human pose estimation, due to its computational efficiency and ability to detect people at different scales.

Human pose from video. Surprisingly little work has addressed the combination of monocular pose estimation with tracking in uncontrolled environments. The problem is sometimes referred to as articulated motion parsing [20]. In early work, Ramanan et al. [17] assume there is at least one frame in the video sequence where the person is in an easy to detect pose. Based on this detection, they build a person-specific appearance model and perform independent pose estimation on each frame using this appearance model. A similar approach is used by Buehler et al. [3] who introduce temporal information by identifying key frames with reliable poses. These effectively become anchor frames and they exploit temporal coherence to link poses in the intermediate frames. Ferrari et al. [START_REF] Ferrari | Progressive search space reduction for human pose estimation[END_REF] formulate a spatiotemporal parsing problem and perform simultaneous pose estimation over a set of frames. In addition to single-frame potentials, the model includes temporal dependencies that represent continuity of appearance and pose. These methods rely on static image likelihoods and do not use optical flow.

Sapp et al.

[20] exploit optical flow information to locate foreground contours. This integrates well with their pose estimation method, which exploits image contours and region-based likelihoods. The idea of using flow discontinuities as a cue for pose estimation dates at least to [START_REF] Sminchisescu | Estimating articulated human motion with covariance scaled sampling[END_REF] on 3D body pose estimation in monocular video. In a recent work, Fragkiadaki et al. [START_REF] Fragkiadaki | Pose from flow and flow from pose estimation[END_REF] exploit optical flow for segmenting body parts and propagating segmentations over time.

Body representation. The representation of the body for monocular articulated motion parsing has not received much attention but, we argue, is critically important. While there are a wide range of 3D body representations, here we focus on 2D, where the range of models is remarkably limited. The most common approach represents the body as rectangular (or polygonal) regions [START_REF] Andriluka | Pictorial structures revisited: People detection and articulated pose estimation[END_REF]4,7,[START_REF] Felzenszwalb | Pictorial structures for object recognition[END_REF]15,18] (Fig. 2). The pose of the body can then be represented either explicitly by a kinematic tree [4,15] or by a probabilistic collection of parts [START_REF] Andriluka | Pictorial structures revisited: People detection and articulated pose estimation[END_REF][START_REF] Felzenszwalb | Pictorial structures for object recognition[END_REF].

In contrast to polygonal parts, the Contour People [START_REF] Freifeld | Contour people: A parametrized model of 2D articulated human shape[END_REF] and Deformable Structures [24] models are derived from a realistic 3D model of body shape and better capture the 2D shape of the person including perspective effects, foreshortening, and non-rigid deformations with pose. The DS model is like a pictorial structures model in that its parts are connected only probabilistically and the configuration of parts is inferred using belief propagation.

Note that in [6] body pose is displayed in a way that looks like a DS model but is not. They take the rectangular body parts of a standard PS model and smooth them with the probability distribution for the part. This produces attractive contours but these are unrelated to the model or the inference.

Alternatively Andriluka et al. [START_REF] Andriluka | Pictorial structures revisited: People detection and articulated pose estimation[END_REF] work with a traditional rectangular part model but use shape contexts to learn a model that captures information about the shape of the parts within these regions. This is not an explicit model of shape and is consequently not appropriate for our task. Guan et al. [12] propose a 2D model of clothed body shape but the model is not articulated. For our application, there are two properties that a representation must satisfy. First, it must be able to represent occlusion. Pure part-based models like pictorial structures, have no notion of what is in front of what. To make sense of the optical flow in the scene this is necessary and is supported by our model (though crudely here). Second, it should approximate the shape of the body. As we will see, to "flow" the puppet in time requires that we associate observed optical flow with body parts. If the parts match the size and shape of parts in the image, this is easy. Rectangular parts would make this much harder. Consequently, here we take the idea of the DS model but, instead of a distributed collection of parts, our state space represents the full pose of the model; this allows us to model occlusions. We augment the DS model with a scale parameter to capture the overall size of the person in the image. The DS model is learned for random poses and cameras. We introduce prior knowledge on camera location and pose in TV shows by redefining the mean DS shape as the average shape in the Buffy training set [START_REF] Ferrari | Progressive search space reduction for human pose estimation[END_REF] annotated with the DS model. We also go beyond previous work to learn a new, multi-scale, image likelihood that captures image statistics along the contour of the puppet.

Model

We briefly summarize the DS model and refer the reader to [24] for details. DS is a gender-specific part-based probabilistic model, where contour points of body parts are represented in local coordinate systems by linear models of the form

p i y i = B i z i + m i (1)
where p i are contour points, y i are joint points, z i are PCA (Principal Component Analysis) coefficients, B i is a matrix of principal components, and m i is the mean part contour. Let l i = (c i , θ i , z i ), where c i is the center of the part i and θ i is the part orientation. The correlation between the shape coefficients, z i , and body pose parameters captures how shape varies with pose and is modeled with pairwise Multivariate Gaussian distributions over the relative pose and shape coefficients of connected body parts. The probability of a model instance is factored as:

p(l|π DS ) ∝ (i,j)∈E p ij (l i , l j |π ij ) (2)
where E is the set of pairs of connected parts and π DS represents the model parameters. The DS model does not in- clude a scale variable in the potentials, but a scale factor can be specified and it is used to convert the model from the DS model space to image pixel coordinates. Let x t be a vector of DS model variables and the scale at time t (i.e. x t = [l t , s t ]), let I t be the image frame at time t, and U t,t+1 the dense optical flow between images I t and I t+1 . We define the posterior distribution over the DS model variables and scale for each frame in the sequence of N frames as:

p(X|I, U, π DS ) ∝ N -1 t=1 p(I t+1 |x t+1 )p(x t+1 |x t , U t,t+1 ) N t=1 p(I t |x t ) N t=1 p(l t |π DS ) N t=1 p(s t |π s ) (3)
where ) is the static image likelihood for the frame at t+1, evaluated for xt+1 , which is the "flowing puppet" of x t given the flow U t,t+1 (see below). Here our likelihood uses flowing puppets in the forward direction, but our formulation is general and can be extended to consider flowing puppets generated with backward flow and for more than one time step.

X = [x 1 , ..., x N ], I = [I 1 , ..., I N ] , U = [U 1,

Flowing puppets

Given a DS puppet defined by the variables x t , and given the dense flow U t,t+1 , the corresponding flowing puppet for frame t + 1 is generated by propagating x t to x t+1 through the flow. The conditional probability distribution p(x t+1 |x t , U t,t+1 ) expresses the noisy generative process for the flowing puppet xt+1 .

Exploiting the part-based representation of the DS model, we define a layered model of the body parts with a fixed depth ordering. We assume the torso is the most distant part, then comes the head, the right and the left upper arms, then the right and left lower arms. Figure 3 shows an example of the layer map, where warm colors indicate parts that are closer to the camera. Given the visibility mask for each body part, we consider the corresponding pixels in the optical flow map U t,t+1 . This is where the DS body shape representation becomes important. Figure 1(c) shows a puppet, x t , overlaid on the forward and backward optical flow fields (i.e., computed from t to t + 1 and from t to t-1). We fit an affine motion model to the optical flow vectors within each body part. The resulting puppet flow field is illustrated in Fig. 1(c); this is our estimate for how the puppet should move from frame to frame. We then apply the estimated affine motion to the joints of each part, resulting in predicted puppets, xt-1 and xt+1 , at the adjacent frames (Fig. 1(d), white). Our current process of generating the flowing puppet does not include a noise model, thus the probability distribution p(x t+1 |x t , U t,t+1 ) is simply a delta function centered on the predicted puppet.

Image likelihood

The static image likelihood is defined by three terms: a contour-based term p c (I t |x t ) that encourages alignment of the DS puppet contour with the edges in the image, a color term p s (I t |x t ) that encodes the knowledge that lower arms and hands are likely to be skin, and a hand likelihood p h (I t |x t ) computed from a hand probability map generated from a hand detector:

p(I t |x t ) = p c (I t |x t )p s (I t |x t )p h (I t |x t ). (4)
The DS model we use is learned from a 3D model that does not include hand pose variations, consequently our 2D model does not have separate hand parts with their own articulation parameters. Instead, hands are included in the model of the shape of the lower arm. To define a regionbased likelihood for hand color, we simply consider the image area surrounded by the points in the lower arm that correspond to hand contours.

Similar to [24], we define the contour-based likelihood p c (I t |x t ) as the output of an SVM classifier with a feature vector of HOG descriptors [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] that are steered to the contour orientation and computed along the model contour. In order to obtain a likelihood model that is more robust to scale variations, we compute the features at different scales. Figure 4 shows an example for the upper arm: We use a 3-level pyramid; HOG cells are placed at contour points (blue), inside the contour (red), and outside (green). We learn the contour-based image likelihood from the Buffy dataset [START_REF] Ferrari | Progressive search space reduction for human pose estimation[END_REF].

The skin-color likelihood p s (I t |x t ) is defined by a histogram of skin chroma and hue color coordinates in the CIE L*C*h* color space. We define a log probability map for the image, and compute log skin color likelihood for the hands, and the lower arms without hands, as the average log probability value in the part region. The skin color likelihood is learned from a dataset of skin colors downloaded from the web as in [24].

The hand likelihood p h (I t |x t ) is based on a hand probability map generated by a hand detector using optical flow. The hand detector is defined as in [20] by computing the gradient magnitude of the flow, then learning an SVM classifier. The hand probability map is built as the max response from the detector at each image location over a discrete set of hand orientations. The detector is learned from training images from the VideoPose2.0 dataset [20] where we have manually annotated hands with oriented bounding boxes. Figure 5 shows an example of the hand probability map.

Inference

The DS model state space consists of pose and shape parameters. To reduce the number of variables during inference, we define the configuration of a body as only the location of the joint points y t , marginalizing out the shape parameters, and thus reducing the number of variables. We use a set of 11 joints points, namely nose, neck, right and left shoulders, belly button, right and left hips, elbows and wrists. From these joints, we can easily compute the DS puppet parameters, l i = (c i , θ i , z i ), where the shape coefficients z i represent the expected shape for each part. The state space for a puppet in a frame is then x t = [y t , s t ], where s t is the puppet scale.

We convert the probabilistic formulation, Eq. (3), into an energy

E(X) = λ DS N t=1 E DS (y t ) + λ c N t=1 E c (x t )+ λ s N t=1 E s (x t ) + λ h N t=1 E h (x t ) + λ c N -1 t=1 E c (x t+1 ) + λ s N -1 t=1 E s (x t+1 ) + λ h N -1 t=1 E h (x t+1 )
where E DS (y t ) = -log p(l t |π DS ) (see Eq. ( 2)), E c (x t ), E s (x t ) and E h (x t ) are the energy terms associated with the contour-based, the skin-color based, and the hand-detectorbased image likelihood on the current frame, respectively. E c (x t+1 ), E s (x t+1 ) and E h (x t+1 ) are the negative log likelihoods of the puppet in frame t propagated to the frame t + 1 through the dense optical flow U t,t+1 . We use a uniform prior for scale, as bounds for the scale parameters are set in the optimizer. The weights for the energy terms are estimated with grid search on a validation set.

We seek a maximum a posteriori estimate and minimize the energy using a novel iterative approach based on frame optimization and propagation. For the frame optimization, we adopt Particle Swarm Optimization (PSO) [START_REF] Ivekovič | Human body pose estimation with particle swarm optimisation[END_REF]. PSO searches the parameter space by perturbing the particles. Perturbing the vertices can produce implausible puppets so we first convert the pose into a joint angle representation, do this perturbation in joint angle space, convert back to joint positions, and then to the expected DS model to obtain contours and regions. Inspired by [START_REF] John | Markerless human articulated tracking using hierarchical particle swarm optimisation[END_REF] we employ a hierarchical strategy optimizing first the torso, head and right arm, then the whole puppet. In addition, we resample the particles according to their energy. PSO requires setting bounds for the variables to optimize: We estimated bounds for scale, limb angles, and limb lengths from the VideoPose2.0 training set.

The process of optimization and propagation is visually described in Fig. 6. We start by initializing a set of P particles on each frame (Fig. 6, first row). Then the video sequence is scanned forward and backward to propagate the best M particles from a frame to the next using the flow (Fig. 6, second row). Each frame in the sequence is then optimized in turn, using PSO, starting from the first frame, and proceeding forward for all the frames then backwards. Figure 6, third row, shows examples of particles after PSO. After optimizing pose in each frame, the best M particles are propagated to the neighbors, forward and backward, using the flow (Fig. 6, fourth row). After propagation, each frame has P+2M particles, but only the best P particles are retained for the frame. This process of optimization and propagation iterates for a defined number of runs, R. In our experiments we used P=40, M=5, R=8. We run the optimization with 3 different seeds for the random number generator and select the solution with the minimum energy.

We use the Flexible Mixtures of Parts (FMP) model [23] to provide the optimization with a good starting guess for the solution. We use the code provided by the authors with their model trained on the Buffy dataset. The FMP model generates a "stickman" as output. In order to map the FMP stickman to the DS model at a proper scale, we learn a regression function between the scale of the stickman and DS from manually annotated frames of the VideoPose2.0 training set. To further help the optimizer, we generate additional initial puppets relocating the wrists of the FMP solution to likely hands locations. We exploit both image cues and motion cues for hand detection. We train and use a hand detector based on the method described in [START_REF] Mittal | Hand detection using multiple proposals[END_REF], which uses image features like statistics of image gradients and colors, to provide initial guesses for the hand location. We also exploited the hand detector trained on optical flow described in Sec. 3.2 to generate proposals for wrist positions.

Experiments

The VideoPose2.0 dataset [20] Results are reported as in [20] as the percentage of joints that have a distance lower than a threshold in pixels from the ground truth (Fig. 7). As a baseline we report results for Yang et al. [23], which we use during initialization. Note that it performs rather poorly but was not trained for this dataset. We compare different variants of our method (FP for Flowing Puppets). First, we report results without and with median filtering as applied in [20]. Second, to show the benefit of our optimization strategy, we show results obtained without exploiting the dense flow for propagation and likelihood (FP, -flow). We report performance better than the state-of-the-art for wrists, significantly improving over [23], and at the performance level of [20] for elbows. We can observe that our approach performs significantly better than the static image detector of [23]. A recent paper [START_REF] Fragkiadaki | Pose from flow and flow from pose estimation[END_REF] also performs pose estimation on the VideoPose2.0 dataset, but for testing they select a set of the clips that is different from the one specified in the dataset; a direct comparison is not possible. Figure 8 shows several examples of correctly predicted body pose, with the DS puppet overlaid on the image and on the optical flow. Figure 9 shows some representative failure cases.

Conclusions

Given recent improvements in the accuracy of optical flow estimation, we argue that it is a useful source of information for human pose estimation in video. Here we use flow in a novel way to make predictions about the pose of the body in neighboring frames. Given a representation of body shape, we use the optical flow forwards and backwards in time from a given frame to predict how the body should move, creating what we call a flowing puppet. If the body pose is correctly estimated in the current frame, and the flow is accurate, then our method should accurately predict the pose in neighboring frames. We use this to construct an extended energy function that incorporates image evidence from multiple frames. We also use our flowing puppets to propagate good candidate poses during optimization and to hypothesize putative hand locations.

The approach improves accuracy and robustness relative to a baseline method that does not use puppet flow. If the pose in one frame is ambiguous, it may not be in neigh- boring frames. Thus by integrating evidence over time we may be able to resolve such ambiguities. This represents a novel approach to temporal estimation of body pose in video. In particular there is no temporal prior to restrict the motion. Instead we rely on flow to "link" information across time. We tested the method on the challenging VideoPose2.0 dataset and found that we improved over previous results in the difficult problem of localizing the wrists.

This work opens up many possibilities for future research. Other properties of the flow, such as motion boundaries could be used. Reasoning about depth ordering of the parts should be added. Some sort of temporal reasoning could be included to propagate information beyond the neighboring frames, for example using particle smoothing or a Viterbi-like algorithm. More accurate models of body shape, hair, and clothing might also improve the results. Note that, while VideoPose2.0 contains background motions, more work should be done to evaluate robustness to multiple moving people and other scene motion. Finally, flowing puppets could be used to build a temporally consistent appearance model across several frames, which could provide stronger image evidence.

Figure 2 .

 2 Figure 2. 2D models of human shape. From the left: Cardboard People (in two viewpoints), Pictorial Structures (PS), Contour People (CP) and Deformable Structures (DS).

Figure 3 .

 3 Figure 3. DS puppet layer. (1) Frame; (2) Corresponding puppet layer with parts ordered by fixed order. The warmer the color, the closer to the camera.

Figure 4 .

 4 Figure 4. Contour-based likelihood features. HOG descriptors are steered to the contour orientation and computed at contour points (blue), inside (red) and outside the contour (green) in a 3level pyramid.

Figure 5 .

 5 Figure 5. Hand detection. Example of output from the hand detector trained on optical flow. Image (left), optical flow (center), and hand probability map defined from running a flow-based hand detector on the flow (right).

Figure 6 .

 6 Figure 6. Particle-based optimization. Particles are initialized on each frame (first row), then the M best are propagated through the flow forward and backward (second row).For a defined number of iterations particles are then locally optimized, then the M best are propagated to the neighbors (third and fourth row). Then the best particle on each frame is returned as the solution (last row).

  contains 1225 frames from two popular TV shows (Friends and Lost) corresponding to 44 clips. The dataset is divided into 706 training frames and 519 test frames in 18 clips. For consistency with [20], we use the dataset with only every other frame of the original video sequences. It contains frames at the original size, and frames that have been cropped and rescaled to have the person in the middle of the frame to meet the needs of the pose estimation method of [20]. In contrast to [20] we use the original frames, since we model scale explicitly and estimate during optimization. We annotate the clips for gender and use a DS model of the appropriate gender. The dense flow is computed with the method of [22] in both the forward and backward time direction.

Figure 7 .

 7 Figure 7. Results. Accuracy of elbow (left) and wrist detection (right) for different threshold distances from ground truth. We significantly improve over our baseline (Yang and Ramanan [23]) and over the state-of-the-art (Sapp et al. [20]) in wrist detection. FP stands for Flowing Puppet.

Figure 8 .

 8 Figure 8. Estimated body pose. Successful detection results from 9 test clips are shown (2 frames per clip). Images are shown with the estimated puppet overlaid in white. Below each image is the estimated forward flow field color coded as in [2] with the puppet overlaid in black.

Figure 9 .

 9 Figure 9. Estimated body pose. Examples of failure cases. In all cases the image evidence supports incorrect poses.

Acknowledgements. CS was supported in part by the ERC advanced grant Allegro. We thank Ben Sapp for helpful feedback on his code.


Mohammed El

Delphine Rivasseau

Stéphanie Jehan-Besson

Marinette Revenu

David Tschumperlé

Luc Brun

Marc Duranton

A Time-Consistent Video Segmentation Algorithm designed for Real-Time Implementation

In this paper, we propose a time consistent video segmentation algorithm designed for real-time implementation. Our segmentation algorithm is based on a region merging process that combines both spatial and motion information. The spatial segmentation takes benefit of an adaptive decision rule and a specific order of merging. Our method has proven to be efficient for the segmentation of natural images (flat or textured regions) with few parameters to be set. Temporal consistency of the segmentation is ensured by incorporating motion information through the use of an improved change detection mask. This mask is designed using both illumination differences between frames, and region segmentation of the previous frame. By considering both pixel and region levels, we obtain a particularly efficient algorithm at a low computational cost, allowing its implementation in real-time on the TriMedia processor for CIF image sequences.

I. INTRODUCTION

Segmentation of videos into homogeneous regions is an important issue for many video applications such as regionbased motion estimation, image enhancement (since different processing may be applied for different regions), 2D to 3D conversion (segmentation can be used for depth estimation). These applications require two main features from segmentation: accuracy of regions boundaries in the spatial segmentation and temporal stability of the segmentation from frame to frame. Spatial segmentation can be classified into two main categories, namely contour-based and region-based methods. In the first category, edges are computed and connected components are extracted [START_REF] Iannizzotto | Fast and accurate edge-based segmentation with no contour smoothing in 2-d real images[END_REF]. The main drawback of such an approach is that the computation of the gradient is prone to large errors especially in noisy images. Moreover we cannot take benefit of statistical properties of the considered image regions. The second category of methods, i.e. region-based, is less sensitive to noise. We are interested here in a bottom-up segmentation approach. In these methods, two important points must be considered: the first one is the order of merging and the second one is the similarity criterion, see for example [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF], [START_REF] Shi | Normalized cuts and image segmentation[END_REF], [START_REF] Sharon | Fast multiscale image segmentation[END_REF]. When dealing with video segmentation, the temporal dimension must be added. Various approaches have been tested. Some authors consider video data as a volume [START_REF] Wang | Automatic video object segmentation via 3d structure tensor[END_REF], while others take benefit of motion information, such as scene change detection or motion field [START_REF] Deng | Unsupervised segmentation of colourtexture regions in images and video[END_REF], [START_REF] Wang | Unsupervised video segmentation based on watersheds and temporal tracking[END_REF]. Other works propose the object tracking [START_REF] Moscheni | Spatio-temporal segmentation based on region merging[END_REF] which is beyond the topic of this paper.

In this paper, we propose a spatial segmentation that takes benefit of both an adaptive decision rule and an original order of merging. This method gives good results for spatial segmentation with few parameters to be set. We use motion information to improve the temporal consistency. Motion estimation is a real bottleneck for real-time implementation, so we rather propose to combine an improved Change Detection Mask (CDM ) with spatial segmentation in order to improve the temporal consistency of our segmentation. By making comparisons both at a pixel level and at a region level, we obtain an efficient algorithm for video segmentation at a low computational cost. Our method runs in real-time on the TriMedia processor for CIF image sequences. Moreover, experimental results conducted on real video sequences demonstrate a good temporal consistency.

The paper is organized as follows. The spatial segmentation method is detailed in section II. The temporal consistency improvement is detailed in section III. In section IV, we discuss the implementation of the algorithm. Experimental results and measures are given in section V.

II. SPATIAL SEGMENTATION

Let us consider an image I, the notation |.| represents the cardinal and I(p, n) the pixel intensity at position p = (x, y) T in the frame n.

A region-based segmentation problem aims at finding a relevant partition of the image domain in m regions {S 1 , S 2 , .., S m }. The algorithm of segmentation detailed here uses an implicit encoding of the initial 4 connected planar sampling grid of pixels. It combines a specific order of fusion with an adaptive threshold. Both steps are detailed thereafter.

A. ORDER OF MERGING

The order of merging is built on the edges weights as in [START_REF] Nock | Statistical region merging[END_REF], [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF]. The idea behind this order of merging is to merge first similar regions rather than different ones. An edge e denotes a couple of pixels (p, p ) in a 4-connectivity scheme. The similarity between pixels is measured by computing the distance between two pixel colors as follows:

w(p, p , n) = I∈{Y,U,V } (I(p, n) -I(p , n)) 2 . (1)
For our algorithm, we consider the Y U V color space, since this is the native color space of CIF sequences. The color space (L * a * b * ) provides partitions with a little greater subjective quality but with a higher computational cost.

The edges are sorted in increasing order of their weights and corresponding couples of pixels are processed in this order for merging. As far as the implementation is concerned, the image is only scanned twice for this sorting: One time, in order to compute the number of edges with same weights and store this number in a table (edge weight histogram). Second, in order to store each edge in the memory part corresponding to its weight.

B. THE CRITERION OF MERGING

Given two regions S 1 and S 2 , we want to know if these two regions have to be merged. A similarity criterion between the two regions must then be chosen and evaluated. Couple of adjacent regions will be considered as similar and merged if their merging criterion is below a given threshold. The choice of the corresponding threshold is often difficult. In this paper, we propose an adaptive threshold for a similarity criterion based on the means of the intensities of these two regions. This adaptive threshold is justified using statistical inequalities as in [START_REF] Nock | Statistical region merging[END_REF] but we here propose a simpler statistical interpretation of the image which leads to a more adapted criterion for real-time implementation. Such a predicate gives very good segmentation results as shown in Fig. 2.

1) Merging predicate: The mean of the intensities of the region S i is computed as follows:

S i = 1 |S i | k=|Si| k=0 I i (p k , n)
where I i (p k , n) is the intensity of the k th pixel of region S i . The merging predicate is then :

P (S 1 , S 2 ) = true if (S 1 -S 2 ) 2 ≤ Qg 2 ( 1 |S1| + 1 |S2| ) f alse otherwise
(2) where g is the maximum level of I (g = 255 for grayscale images). The parameter Q allows to tune the coarseness of the segmentation. In the experiments, we choose Q = 2 which gives good results for CIF (Common Input Format = 352 × 288) images.

2) Statistical justification of the predicate: Classically, the image I is considered to be an observation of a perfect image I * and then pixel intensities are considered as observations of a vector of random variable (r.v.) noted X = (X 1 , .., X n ) T . The merging predicate is based on the use of statistical inequalities and especially the McDiarmid inequality as in [START_REF] Nock | Statistical region merging[END_REF].

Theorem 2.1: (The independent bounded difference inequality [START_REF] Habib | Probabilistic methods for algorithmic discrete math[END_REF]) Let X = (X 1 , X 2 , ..., X n ) be a family of n independent r.v. with X k taking values in a set A k for each k. Suppose that the real valued function f defined on k A k satisfies |f (x) -f (x )| ≤ c k whenever vectors x and x differ only in the k th coordinate. Let E(f (X)) be the expected value of the r.v. f (X), then for any τ ≥ 0,

P r(|f (X) -E(f (X))| > τ ) ≤ 2 exp -2τ 2 / k c 2 k (3)
When testing the deviation between two adjacent regions S 1 , S 2 , we consider the following vector of random variables:

(X 1 , ..., X n ) = I * 1 (p 1 ), ..., I * 1 (p |S1| ), I * 2 (p 1 ), ..., I * 2 (p |S2|
) where I i (p j ) is the intensity of the j th pixel of S i corresponding to the observation of the random variable I * i (p j ). In this case, the size of the random vector is n = |S 1 | + |S 2 |. In order to apply Theorem 2.1, we then choose f (x) = S 1 -S 2 . By inversion of the theorem, we have with a probability of at least 1 -δ (0 < δ ≤ 1):

|(S 1 -S 2 ) -E S 1 -S 2 | ≤ g Q 1 |S 1 | + 1 |S 2 | where Q = 1 2 * ln( 2 δ )
, and E(S) is the expected value of S. If S 1 and S 2 belong to the same region in I * , the expectation E S 1 -S 2 will be null and the predicate follows.

C. SEGMENTATION ALGORITHM

Our spatial segmentation could be divided in three steps. In the first one, we compute the weights of edges and their histogram. In the second step we sort edges in an increasing way of their weights. In the last step we merge pixels or regions connected by edges following their order. The implementation of this last step is given in section IV.

III. TIME CONSISTENCY IMPROVEMENT

In video segmentation, the quality of the spatial segmentation is not the only requirement, time consistency is also a very important one. If in two successive frames, one region is segmented very differently because of noise, occlusion or deocclusion, results of segmentation would be very difficult to exploit for any application like image enhancement, depth estimation and motion estimation. Many works, see for example [START_REF] Wang | Unsupervised video segmentation based on watersheds and temporal tracking[END_REF], use motion estimation to improve time consistency in video segmentation. However, motion estimation is a real bottleneck for real-time implementation and is even sometimes unreliable. In this paper, we combine an improved Change Detection Mask (CDM ) with spatial segmentation in order to improve the temporal consistency of our segmentation.

The CDM is designed using both illumination differences between frames and region segmentation of the previous frame. We first detect changing pixels using the frame difference. Then, we take benefit of the region segmentation of the previous frame in order to classify the pixels not only at a pixel level but also at a region level. Given the current frame I(:, n) and the previous one I(:, n -1), the frame difference F D is given by F D(p) = |I(p, n) -I(p, n -1)|. Classically, F D is thresholded in order to distinguish changing pixels from noise. A pixel p, with L(p) = 1 is denoted a changing pixel. We then use the previous segmentation in order to convert the CDM from the pixel level to a region level which is more reliable. For each region in the previous segmentation S i , we compute τ (S i ) which represents the ratio of changing pixels in the region S i . We have τ

(S i ) = N i,changing |Si|
where N i,changing is the number of changing pixels in the region S i . Pixels are then classified using three categories:

where tr 2 is a positive constant. In the experiments, we take tr 2 = 0.01 (i.e. a region is a moving region when it contains at least 1% of moving pixels). The value of the threshold is chosen so that we don't miss any moving region. Every pixel of regions qualified as static are labelled using CDM (p) = 0. The two other labels concern pixels within moving regions. Depending on the value of the frame difference, the pixel is qualified as a changing one (CDM (p) = 1) or a no changing one (CDM (p) = 2). Such a classification is then used to segment the current frame. Firstly, static regions are kept as they were segmented in the previous frame. Secondly we apply a connected component labelling (CCL) algorithm [START_REF] Wu | Optimizing connected component labeling algorithms[END_REF] to extract connected component of pixels with CDM (p) = 2. This second step builds seeds from the segmentation of the previous frame. These seeds link the current segmentation to the previous one in a time consistent way. Thirdly, we apply the spatial segmentation only on edges (p, p ) connecting changing pixels, and those connecting changing pixels with seeds built in the second step. We show here an example of segmentation results with and without time consistency in Fig. 1(d) and 1(c). In Fig. 1(b), the different labels (i.e. values of CDM ) are given using three values of intensity (black for CDM (p) = 0, gray for CDM (p) = 2 and white for CDM (p) = 1).

In section V, we propose the computation of an objective measure for temporal consistency. The results show a real improvement for time consistency for different sequences. Moreover, the way we exploit the CDM decreases also the amount of computation of the algorithm since the edges in static area are not reconsidered, and those linking no changing pixels in moving area are simply processed by a CCL algorithm.

IV. IMPLEMENTATION CONSIDERATIONS

In this section, we propose to describe optimizations that have been made to allow a real time treatment. The whole algorithm of video segmentation is summarized in figure 3. The algorithm 1 describes more particularly the merging loop which is the critical part for a real time implementation. The term N e represents the number of edges within the image I in the 4-connectivity. In the merging process, we use the UNION-FIND data structure [START_REF] Fiorio | Two linear time union-find strategies for image processing[END_REF]. The UNION function fuses Fig. 2.

Results of segmentation of frames from the sequences P aris, Bridge, Akiyo and M obile. Apart from this merging loop, all other functions access pixels data in a predictable way (for example from top to bottom left to right). The cache memory benefits from this regularity, since it exploits spatial and temporal locality of data, and consequently causes less cache misses. In the merging loop, the UNION-FIND data structure is unpredictable, and consequently causes an important data cache stalls. To reduce the data cache stalls cycles, we investigate two optimizations, which are independent of the cache memory type (size, mapping, block size, replacement algorithm...).

The first optimization concerns the organization of the data. Each pixel in the image could be a representative of one region so the data are stored with an array A of size |I|, where each element A(p) stores the mean colors (Y , U , V ) and the cardinal |S| of the region whose representative is p. We verify that this approach reduces the number of data cache stalls. The second optimization concerns the UNION-FIND data structure. It is encoded separately thanks to an array F of size |I|, where F (p) stores the address of the father of p. Each region is thus encoded by a tree using such fatherchild relationships. The alternative solution, which consists in storing both the labels and the data in the same structure, would induce the storage of useless data within the cache during the FIND operations.

We experimented these optimizations on the TriMedia processor [START_REF]pnx1500 databook[END_REF]. The cache memory of this particular TriMedia is 128 KByte, 4 way associative, with block of 128 Byte. The replacement algorithm used is LRU . We get a reduction of 3 M Cycles per frame. Notice that these optimizations are useful for many segmentation algorithms that use the UNION-FIND data structure.

There is an important amount of DLP (data level parallelism) in our algorithm (computation of edge's weight, frame difference, classification of pixels in CDM ). This allows to increase the throughput (i.e. amount of pixels processed per unit time), by processing data in parallel when it is possible. The core of TriMedia is a VLIW architecture with 5 issues slots. Each slot has some functional unit, and each functional unit could process 4 bytes in parallel (SIMD mode). The ILP (Instruction Level Parallelism) is extracted by the compiler, while the DLP could be exploited through the use of custom operations, loop unrolling, and grafting. So we use these optimizations to exploit the DLP available in our algorithm.

V. EXPERIMENTAL RESULTS

In this section we present experimental results of our algorithm run on TriMedia with many very known CIF sequences. The table I shows the following measures:

1) We use a classical measure to evaluate time consistency.

Given the segmentation of previous frame SEG(n -1) and the segmentation of the current one SEG(n), we find the correspondence between regions in SEG(n -1) and SEG(n). Two regions (S i,n-1 ∈ SEG(n -1) and S j,n ∈ SEG(n) correspond if they have the most overlapping area Overlap(i, j) = |S i,n-1 ∩ S j,n |. We sum the number of pixels in the overlapping areas of corresponding regions over the image. The consistency measure is the percentage of this number to the size of the image. 2) We give the number of Mcycles the algorithm takes on TriMedia. The exploitation of the CDM reduces the computational cost. This reduction depends on the correlation between two successive frames. Since the execution of spatial segmentation without CDM is the critical case of our algorithm, we show the number of Mcycles corresponding to this critical case.

When enforcing consistency through the CDM , time consistency is higher, and visually, segmentation is more stable from frame to frame and still fit very well regions boundaries as shown in Fig. 2. With a 450 M Hz TriMedia, we are able to

Fig. 1 .0

 1 Fig. 1. (a) Segmentation of a frame n-1. (b) CDM computed from the segmentation of the frame n-1 and frame difference between frame n-1 and frame n. (c) Segmentation of the frame n without CDM . (d) Segmentation of the frame n with CDM .

Fig. 3 .

 3 Fig. 3. The general diagram of video segmentation

Algorithm 1

 1 The merging loop for i := 1 to N e do Read the i th edge: (p 1 , p 2 ); S 1 = F IN D(p 1 ); S 2 = F IN D(p 2 ); if P (S 1 , S 2 ) = T rue then U N ION (S 1 , S 2 ) end if end for

  for video processing requires low computational methods. Directed by this constraint, we propose here an efficient time consistent algorithm for video segmentation. This method gives accurate results with high temporal consistency. It is based on a statistical spatial segmentation and on the computation of an improved change detection mask. It runs in real-time on only one processor. Experimental results demonstrate the applicability of this method.

			TABLE I
	EXPERIMENTAL MEASURES OF OUR ALGORITHM
	Sequence	Akiyo	Tennis Table	Paris	Mobile
	Time Consistency
	(without CDM )	0.88	0.73	0.89	0.84
	Time Consistency
	(with CDM )	0.98	0.92	0.97	0.92
	Mcycles per frame
	(without CDM )	15.68	15.84	16.59	16.20
	process more than 25 frames per second, which is sufficient
	for real-time execution.
		VI. CONCLUSION
	Designing usable algorithms

ACKNOWLEDGEMENT

The authors would like to thank Patrick Meuwissen, O.P. Gangwal and Zbigniew Chamski for their constructive suggestions.


J Zhu

R Le Bouquin Jeannès

C Yang

J J Bellanger

H Shu

PARTIAL MUTUAL INFORMATION FOR SIMPLE MODEL ORDER DETERMINATION IN MULTIVARIATE EEG SIGNALS AND ITS APPLICATION TO TRANSFER ENTROPY

Keywords: Effective connectivity, partial mutual information, transfer entropy

The context of this work is the analysis of depth electroencephalographic signals recorded with depth electrodes during seizures in patients with drug-resistant epilepsy. Usually, different phases are observed during the seizure process and we aim to determine how cerebral structures get involved during these phases, in particular whether some structures can "drive" other ones. To this end, we consider a pair of signals and use transfer entropy which needs beforehand to choose efficiently the size of two conditioning vectors built on the past values of these signals.

In this contribution, we extend a partial mutual information based technique, first developed for monochannel prediction models, to the case of two channels. Experimental results on signals generated either by a linear autoregressive model or by a physiology-based model of coupled neuronal populations support the relevance of the proposed approach.

INTRODUCTION

Over the last decade, some measures have been developed to deal with functional connectivity and effective connectivity between different cortical sites, such as in chaotic systems and multivariate neurobiological signals [START_REF] Engel | Outcome with respect to epileptic seizures[END_REF][START_REF] Schreiber | Measuring information transfer[END_REF][START_REF] Sabesan | Information flow in couple nonlinear system: Application to the epileptic human brain[END_REF][START_REF] Sabesan | Information flow and application to epileptogenic focus localization from intracranial EEG[END_REF], which is an important topic in neuroscience. In [START_REF] Schreiber | Measuring information transfer[END_REF], a theoretical information measure named Transfer Entropy (TE) was proposed to identify the direction of the information flow and to quantify the strength of coupling between complex systems. This model-free technique can be considered as an extension of the linear Granger causality index, which is based on a parametric linear modeling of observations, to a non parametric index well suited to non linear signals. When defined from signal Y to X , TE depends on two transition (conditional) probabilities measures characterizing stochastic dynamical links. One of them considers transitions from past values of X to the current value, and the other one considers transitions from past values of X and to the same current value. Its estimation can be strongly dependent on the sizes (named k and in the sequel) of the two vectors containing the past values of

Y l

X and Y respectively. The problem addressed in this contribution is to choose k and for pairs of signals that could be not well characterized by linear AR models. If we assume AR modeling, the standard Akaike methodology can be used. If we do not assume it, a possible approach is to select nonlinear predictors in a given class. This implies order determination and parameters estimation. For one channel (signal l X conditioned only on its own past), a simpler approach based on partial mutual information (PMI) was first proposed by Sharma in [5] to build a robust and effective probabilistic forecast model. The corresponding algorithm requires only computation of mutual information for 2D distributions. Compared to the Mutual Information criterion (MI) [6], PMI takes interdependences between candidate variables into account and is suitable for both linear and nonlinear models. In this paper, transfer entropy is first presented in section 2.1, before giving some details on PMI (section 2.2.1). In section 2.2.2, we extend PMI to select predictors in two-channel signals using two approaches for characterizing models orders and we introduce them in TE estimation algorithm. Experimental results are presented in section 3 before drawing some conclusions.

METHODS AND MATERIALS

Transfer Entropy

In the sequel,

p n u denotes 1 1 , ..., n n n p u u u .
Considering a k -th order Markov process X , we have

' ' 1 1 1 1 / / ' : k k k k n n n n n n x X X x k k dP x dP x n n X X (1)
By considering the auxiliary random process Y , relation [START_REF] Engel | Outcome with respect to epileptic seizures[END_REF] can be extended to formalize the absence of information flow from Y to X (for given and values):

k l 1 1 1 1 / / , k k k k l l n n n n n n n n n n X X x X X x Y y dP x dP x (2)
The deviation from this assumption can be quantified using the Kullback pseudo-metric, which leads to define the transfer entropy from Y to X

1 1 1 , , 1 / , 2 1 / 1 log , , k k l l n n n n n k l k k n n n k l X X Y n n n X X x Y y Y X n X X x k l n n n dP TE x dP dP x x y (3)
where the ratio in (3) corresponds to the Radon-Nikodym derivative of the numerator conditional probability measure with respect to the denominator's one [7]. TE is not symmetric. It can be estimated from observations , , 1,...,

n n x y n N , using a kernel discrete estimation of 1 , , k l n n n X X Y distribution [2]: 1 1 1 , ˆ, , , 1,..., 1 n m k l k k n n n n m m l l n m n r x x P x x y x x r y y C n k N (4)
which depends on a neighborhood size (radius r ). Then, it can be used to compute the estimation

1 1 2 1 ˆ| , ˆ, , log ˆ| k l n n n k l Y X n n n k n n n P x x y TE P x x y P x x (5)
where conditional probabilities in (5) are obtained from estimated joint probabilities in [START_REF] Sabesan | Information flow and application to epileptogenic focus localization from intracranial EEG[END_REF]. is defined by . X Z is estimated as a kernel estimation from realizations ( , )

i i x z ,
. With a Gaussian kernel whose covariance is chosen to be equal to the sample joint covariance of ( ,

1,.., i n ) X Z , then | E X Z can be written as 1 [ | ] 1 i i i ZZ XZ n T E X Z w x Z z S S i (7)
where is the sample cross-covariance between XZ S X and Z , ZZ S is the sample covariance of Z , and is given by: i w 1,.., ( )/ ( )

i i l n w GK Z z GK Z z l (8) with 1 1/2 2 2 ( ) 2 exp 2 T K ZZ ZZ u S u GK u S .
is a smoothing parameter chosen as in [5]. Let us consider a univariate stationary process n X and suppose we want to choose an efficient subset Z of regressors in the set max , 1,..,

c n i Z X i n to predict n X (i.e.
n Y X in Figure 1), where is the maximal predictor lag value. The algorithm proposed by Sharma to build max n Z is described in Figure 1.

Extension of PMI to 2 signals

We proposed two approaches (Method 1 and Method 2) to extend the original PMI to a bivariate time series to determine the orders and l needed to compute as in section 2.

EXPERIMENTAL RESULTS

In this section, we used Method 1 and Method 2 to compute X Y TE and on two types of signals. The first one was a toy linear AR model and the other one was a realistic EEG model. Each experiment was repeated 100 times and the following results correspond to averaged values.

Y X TE

Linear Model

For the linear stochastic system we considered, the following two signals were generated:

1 4 ' 1 0.6 0.4 0.6 n n n n n n n X X X Y X W W (9)
where and were independent white Gaussian noises with zero means and unit variances.

Y X

Physiology-based Model

We used a physiology-based time continuous SDE (Stochastic Differential Equation) model to represent the electrical activity of two distant -and possibly coupledneuronal populations and . This model produces two outputs X Pop Y Pop X and which can be interpreted as two iEEG (intracranial EEG) signals similar to those recorded with proximal electrodes [8]. In the model, each population contains three subpopulations of neurons that mutually interact via excitatory or inhibitory feedback linking main pyramidal cells with two other types of local interneurons. Since pyramidal cells are excitatory neurons that project their axons to other areas of the brain, the model accounts for this organization by using the average pulse rate of action potentials from the main cells of one population,

, as an excitatory input to main cells inputs of another population, . The connection from the first population to the second one is represented by a parameter which is proportional to the number of corresponding axonal links. An appropriate setting of this

CONCLUSIONS

In this paper, we focused on information propagation between two observations using TE and employed the extended PMI method in the determination of models orders.

The proposed approach appeared efficient to reveal the information flow direction between two signals in unidirectional case. In a future work, we will test it on multivariate systems including bidirectional flows. We also plan (i) to conduct statistical analysis on the estimators we proposed to compare them with other order determination methods, and (ii) to take into account supplementary EEG observations for conditioning surrounding information. [6] A.M. Fraser and H.L. Swinney, "Independent coordinates for strange attractors from mutual information", Physical Review A 33(2): 1134-1140, 1986.

[7] R. M. Gray, "Entropy and information theory", New Work, USA, Springer Science+business Media, 2011.

[8] F. Wendling, A. Hernandez, J.J. Bellanger, P. Chauvel, and F. Bartolomei, "Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG", Journal of Clinical Neurophysiology, vol. 22, pp. 343-356, 2005.

X 2 . 2 .

 22 of the values. In[START_REF] Sabesan | Information flow and application to epileptogenic focus localization from intracranial EEG[END_REF], the value of k is determined heuristically from the autoand the parameter is set to 1.l Determination of k and l using PMI

ZStartFig. 1 .

 1 Fig. 1. Algorithm flow chart of PMI

  the corresponding TE values vs. ln r (Figure2.B). According to Figure2, transfer entropy is well estimated using both methods. The flow propagation from signal Y X to signal Y is correctly established, the highest dynamics being obtained using Method 1. It is also clear that there was no influence of signal to Y X .

2 Fig. 2 .

 22 Results of Transfer Entropy between X and Y in AR4 model (Eq. (9)). Top: Results for method 1. Bottom: Results for method 2. A. Plot of the mean of C r (in ln scale) over the trials. B. Plot of TE (bits) vs. radius (in ln scale), for direction of flow r X Y (blue solid line) and (red dotted line).

  effective connectivity. Other introduced parameters are intra-population parameters whose tuning allows modifying the type of activity (normal/epileptic). This model was used to simulate 100 blocks of 4-second length signals. Sampling rate was equal to 256 Hz. Model parameters were such that a fast quasisinusoidal (25 Hz) activity (similar to that observed at seizure onset) was generated by the two populations which were unidirectionally coupled as already mentioned.Referencing to [8], we used parameters values for ( A , ,) equal respectively to (5, 3, 20) and (3.5, 3.5, 84) in and with (unidirectional flow). Experimental results are shown in Figure3. Even if the dynamics of transfer entropy is reduced compared to the linear model, the second method appears relevant and outperforms the first one in detecting the flow propagation from

2 Fig. 3 .

 23 Simulation results on physiological model. Top: Results for method 1. Bottom: Results for method 2. A. Plot of the mean of C r (in ln scale) over the trials. B. Plot of TE (bits) vs. radius r (in ln scale), for direction of flow X Y (blue solid line) and Y (red dotted line). X [5] A. Sharma, "Seasonal to interannual rainfall probabilistic forecasts for improved water supply management: Part 1-A strategy for system predictor identification", Journal of Hydrology 239(1-4): 232-239, 2000.

  1.

						k	TE	Y X
	Method 1: At the first step, the candidate set
	Z	c	Z	cx	cy Z initially contains lag versions of both X
	( cx Z ) and Y ( cy Z ), and the forecasting target is n X .
	Running PMI algorithm (see Figure 1), we get
	Z		Z	xx	x Z where xx y Z and xy Z are subsets of selected
	predictors extracted from past values of X and	Y
	respectively. The orders k and are finally equal to the l
	maximal lags in	xx Z and xy Z .
	Method 2: We run PMI algorithm twice. Firstly, we do not
	use Y . From the initial lag versions of X , cx Z , we obtain
	a new set xx Z . In a second step, we keep this set xx Z and
	the algorithm searches in the set

ACKNOWLEDGMENTS

The authors thank F. Wendling for providing the database on the physiological model and A. Sharma for fruitful discussions on PMI.


Nazim Fatès
email: nazim.fates@loria.fr

Is there something like "modellability

Invitation by Pr. Julio Aracena.

In his seminal paper On computable numbers with an application to the Entscheidungsproblem of 1936, Turing opened the way to consider computability by showing the limits of computation. In this talk we will ask if there is anything that could be called "modellability" and would be analogous to computability. To this end, we propose to examine the dynamical behaviour of various cellular automata and multi-agent systems and endeavour to evaluate the robustness of these systems to various perturbations.

1 About the robustness of the most simple rules α-synchronous CA One easy way to obtain a probabilistic cellular automata (CA) is to alter the global transition function and to consider α-synchronous CA in which cells are updated according to a random local choice with probability α, the synchrony rate. This means that for each cell, independently, the transition results in applying the rule with probability α and staying in the same state with probability 1α.

A first experimental study of such CA was carried out on the set of the 88 minimal ECA. This study revealed that there is no straightforward links between the classical synchronous behaviour of the rules and their robustness. In particular, there seems to be no obvious correlation between the well-known empirical classification of Wolfram and the classes of robustness that can be defined.

Fully asynchronous CA

In order to start studying CA with an analytical approach, one has to start with the most simple systems. Fully asynchronous CA are defined as applying one cell chosen randomly and uniformly at each time step. In other word, the choice of the cell is memoryless and the local transition function is applied on only one cell, while the state of the other cells is left unchanged. Note that the definition only applies for finite set of cells (|L| < ∞) and can not be extended directly for infinite CA. A first classification of the rules with two quiescent states (that is, such that f (0, 0, 0) = 0 and f (1, 1, 1) = 1) was proposed. The classes are defined with respect to the worst average expected convergence time to a fixed point. The interesting point is that there does exist strong correlation between the behaviour of a CA and its time of convergence (more exactly to the scaling time of convergence). It is still an open problem to extend this classification to the 88 ECAs or to the two-dimensional rules. See Ref. [START_REF] Fatès | A note on the classification of the most simple asynchronous cellular automata[END_REF] for a recent synthesis on this question.

So at this point the question we face is : what is the origin of this difficulty to understand the robustness of such simple models of complex systems?

Is it because of the non-existence of objects, entities that carry out the information from one place to the other, like in most "physical models" ? In cellular automata, transmission of information is "abstract" and even when we do α-asynchronous updating, we still assume that each cells reads simultaneously all the states of the neighbours. For an systematic numerical study of the relaxation of this hypothesis, see the work on so-called β-and γ-asynchronism by Bouré et al. [START_REF] Bouré | Probing robustness of cellular automata through variations of asynchronous updating[END_REF].

2 Are multi-agent systems "different" ?

We will now consider different models of multi-agent systems. The first one, the multi-turmite model, is the most closely related to cellular automata. The second one is a model of the swarming phenomenon ("collective motion" is currently a fashionable topic). The third one shows how a swarm of virtual amoebae may exploit an active environment in order to realise a decentralised gathering.

Multi-turmite model (Langton's ants)

Langton ants, or turmites (name given by Dewdney), are one of the most simple multi-agent system where the local rule is simple and the global dynamics complex ; they were initially introduced by Langton to study "artificial life". We refer to the work of Gajardo et al. for the study of the dynamics of a single ant. When we take the simultaneous presence of multiple ants on the grid, very few works are found in literature. We propose a way to model such a simple multi-agent system as a discrete dynamical system [START_REF] Chevrier | Multi-agent systems as discrete dynamical systems: Influences and reactions as a modelling principle[END_REF]. We show that even for the synchronous update, various updating schemes can be defined, depending on the way the rules are interpreted. We present examples where starting from the same initial conditions with two turmites, various updating schemes lead to qualitatively different evolutions. Novel observations are made such as the formation of cyclic behaviours, infinitely extending patterns, deadlocks, etc.

The case of synchronous vs. asynchronous updating was also tackled recently. This requires to extend the formalism and also leads to interesting new observations such as the existence of gliders that are robust to a change from synchronous to sequential (deterministic) updating [START_REF] Belgacem | Robustness of multi-agent models: the example of collaboration between turmites with synchronous and asynchronous updating[END_REF].

In this case, we see that it is difficult to find robust phenomena. Hence the question, is it because the system is somewhat "too discrete" (but the majority rule is robust) ? because it is deterministic ? etc.

LGCA model of swarming Swarming refers to the capacity of for huge ensembles of animals to form a coherent group and move together. This phenomenon is often seen with birds (starlings), fishes, amoebae, etc. We are interested in re-visiting a model by A. Deutsch in order to assess its robustness. The model is simple enough as it has only two parameters : ρ, the initial density and σ, the sensitivity to control to which extent a particle may align with its neighbours. The main studied phenomenon is the existence of a phase transition from a disorder state where no pattern appears to an order state where particles tend to align and form diagonal stripes. However, a more extensive studies reveals that there exist (at least) new phenomena : (a) if the sensitivity σ is made sufficiently high, one no longer observes stripes but rather clusters, (b) if the density of particles ρ is high enough, we observe an astounding new behaviour as particles tend to anti-align instead of aligning! [START_REF] Bouré | A robustness approach to study metastable behaviours in a lattice-gas model of swarming[END_REF]. This patterns disappears with an asynchronous updating of the LGCA, while other patterns seem robust, at least up to a certain degree [START_REF] Bouré | First steps on asynchronous lattice-gas models with an application to a swarming rule[END_REF]. It is an open problem to characterise the nature of the phase transitions between patterns, in particular to determine if the transition from order to disorder is continuous or not.

Discussion

It is now time to come back to our three initial question: "Is there something like modellability ?" Clearly, the time has not yet come where an analogue to computability can be defined in the context of modelling ; however, there seems to be the need for it... In the context of cellular automata and simple reactive multi-agent systems, we asked the following three questions:

How can we think about the limits of modelling ?

We examined some examples of simple models where a simple variation, namely the updating scheme, often resulted in qualitative change of behaviour. The case of binary cellular automata shows that this change is not systematic and that there exist rules which are robust. The study focuses on cases of sensibility as this challenges the validity of our models of complex systems? This leads us to question: what is a good model for synchrony in such models? Do the limits shown by Turing also apply to the models of physical situations?

In many models of "natural" systems, which are stochastic, the undecidability results do not apply. Nevertheless, we see very puzzling phenomena by numerical simulations and there are precise cases as with non-equilibrium phase transitions (e.g. directed percolation) where we know that analytical results are hard to obtain.

What advantages can we expect from the examination of the robustness of discrete models ?

Showing the limits allows us to replace models in a restricted context of validity. One phenomenon becomes valid only in the region where it is robust. In the other regions, new phenomena can be observed. These new phenomena are often interesting by themselves. An image to conclude about broadening the perspective: the difference between the Mona Lisa by da Vinci and the Marriage at Cana by Veronese...