Isaac Cohen

Isabelle Herlin

Sonia Bouzidi

Detection and Tracking of Vortices on Oceanographic Images

Keywords: Oceanographic images, Optical flow, Polynomial phase portrait, Geometric modelling, Nonlinear regularization

This paper deals with the problem of automatic interpretation of oceanographic images for vortices detection, modelization and tracking. We present a framework allowing the vortices detection. The processing is split into three parts: apparent motion computation from the image sequence, local interpretation of apparent motion and geometric modelization of the vortices. This scheme allows an efficient approach for vortices segmentation on very large image sequences. A set of experimental data shows the use of such framework for processing Advanced Very High Resolution Radiometer (AVHRR) and Coastal Zone Color Scanner (CZCS) image sequences.

Introduction

Oceanographics images obtained from environmental satellite platforms present a new challenge in computer science. The huge amount of data collected each day and the needs of characterizing some specific structures on these images for oceanographic monitoring justify our approach for the detection and tracking of vortices on oceanographic images.

This paper deals with the problem of automatic interpretation of oceanographic images for vortices detection, modelization and tracking.

We present a framework allowing the vortices detection. The processing is split into three parts:

• Apparent motion computation from the image sequence.

• Local interpretation of apparent motion.

• Geometric modelization of the interesting regions. This scheme allows an efficient approach for vortices segmentation. The first step concerns the localization of interesting structures on sequences of satellite images. These images are too large (typically 2048 × 1024 images) for modeling the vortices structures over the whole image. Instead we use an apparent motion computation on oceanographic temporal sequences and characterize the vortices center through a local polynomial approximation of optical flow.

In the regions that have been identified to show a turbulent motion we apply geometric modelization to obtain simultaneously a segmentation of the structures (i.e. the external boundary of the vortex), and a quantitative values (according to the model) allowing a study of the temporal evolution.

The aim of this paper, as it will be explained all over the sections, is to show that dynamic satellite images need a hierarchical approach with two steps, localization and segmentation, due to the large amount data. The paper will also explain our choices for this particular application.

We emphasize that this paper presents fundamental problems arising from applicative and context dependent research in an oceanographic application. Although basic tools like optical flow, vector field interpretation and geometric modelization are presented, we show how these methods can be reformulated and improved to incorporate knowledge of the particular application, like the type of movement occurring on a temperature front or on a vortex.

The paper is organized as follows: section 2 deals with the computation of apparent movement on oceanographic temporal sequences, according to the optical flow equation. In this section, we will discuss the choice of a variational formulation allowing flow discontinuities. These discontinuities characterize oceanographic patterns: temperature fronts and vortices on sea surface temperature images obtained by the Advanced Very High Resolution Radiometer (AVHRR) sensor. In section 3 we present a local polynomial approximation of the flow field, allowing an interpretation of the apparent movement and its fixed points. These points allow to initialize a geometric model of external boundary of the vortex as described in section 4. Finally section 5 presents some results obtained on AVHRR and Coastal Zone Color Scanner (CZCS) images.

Optical Flow Computation

Usually [START_REF] Berthold | Determining optical flow[END_REF][START_REF] Schunck | The image flow constraint equation[END_REF], the techniques used for computation of optical flow are based on the constant grey level value hypothesis:

I = constant ( 1
)
where I is the grey level value of a pixel. By differentiating this equation over time when a pixel is apparently moving on an image over a sequence, we obtain:

dI dt = 0 (2)
Equation 2 may be rewritten with partial differentiation:

dI dt = ∂I ∂x dx dt + ∂I ∂y dy dt + ∂I ∂t = 0 ⇐⇒ I x u + I y v + I t = 0 (3)
where I x , I y and I t are the partial derivatives of the image irradiance function I(x, y, t) at point (x, y) and (u, v) is the velocity of this pixel inside the image. Equation 3 is evidently not sufficient for computing the image velocity (u, v) at pixel (x, y) since the two components u and v are constrained by only one equation, called "optical flow equation". Therefore, most of the techniques make use of a regularity constraint which restrains the space of admissible solutions of optical flow equation, and allows to compute the velocity field on each image of the studied sequence. Usually, the regularity constraint is quadratics and therefore enforces the optical flow field to be continuous and smooth. For instance, to guarantee a unique solution to equation 3, one may use a regularizing term by minimizing the energy function:

E(u, v) = Ω (I x u + U y v + I t ) 2 dx + α Ω |∇u| 2 + |∇v| 2 dx dy (4)
where ∇ denotes the gradient operator. Here, the regularity constraint is applied by a quadratic regularizer Ω |∇u| 2 +|∇v| 2 dx dy that constrains the optical flow to be usually continuous.

Unfortunately the interesting structures on oceanographic images like temperature fronts and vortices; correspond to discontinuities in the optical flow. This can also be true for other applications where discontinuities occur on the boundary between two surfaces representing two different objects.

A number of authors [START_REF] Aisbett | Optical flow with an intensity-weighted smoothing[END_REF][START_REF] Nagel | An investigation of smoothness constraints for the estimation of displacement vector fields from image sequence[END_REF] have studied different methods to allow discontinuities of the optical flow. But most of these approaches require the precomputation of the discontinuity points and a specific processing at these points. For example, equation 4 is used with α = 0., at edge point and with a smooth increase of α in stationary regions to insure the smoothness of the resulting optical flow field. In our application, it is impossible to compute these discontinuity points: on the contrary these points are to be found by use of optical flow.

Naturally, this quadratic regularization is not the only one choice possible. Its use in most cases is due to its simplicity in computation since searching for a minimum of E given by Eq. 4 leads to a set of linear equations.

We have defined a new method for constraining equation 3 (see [START_REF] Cohen | Nonlinear variational method for optical flow computation[END_REF] for a complete description), which is more appropriate since it does not enforce the optical flow to be smooth but guarantees a unique solution. Figure 1 describes on the left, the type of real flow fields occurring on an oceanographic front and on the right, the type of optical flow that would be computed with a continuity constraint. We are interested in modelling temperature fronts by use of the discontinuity property. That is not possible if the optical flow is computed with a continuity constraint.

Instead of a quadratic regularizer with L 2 norm (|u| 2 = u 2 ) we use a

L 1 norm (|u| 1 = |u|) regularity constraint. L 1 regularizers preserve sharp
signals as well as discontinuities. Another advantage is the ability of this model to incorporate additional knowledge such as the characteristics of the noise or the confidence on the stationarity constraint.

We solve the problem of optical flow in the space of functions of bounded variations in IR 2 ; i.e. the space of functions

f = (f 1 , f 2 ) such that Ω |∇f 1 | + |∇f 2 |dx dy < +∞.
Computation of optical flow (u, v) can now be stated as a minimization of the following energy functional:

Ω ∂u ∂x 2 + ∂u ∂y 2 + ∂v ∂x 2 + ∂v ∂y 2 dx dy + Ω (I x u + I y v + I t ) 2 dx dy
(5) This energy minimization problem can be translated into Euler Lagrange form:

                   -∂ ∂x   u x u 2 x + u 2 y   -∂ ∂y   u y u 2 x + u 2 y   + (u I 2 x + v I x I y + I x I t ) = 0 -∂ ∂x   v x v 2 x + v 2 y   -∂ ∂y   v y v 2 x + v 2 y   + (u I x I y + v I 2 y + I y I t = 0
+ Boundary conditions (6) Equations in system 6 are non linear and therefore must be processed in a particular way by considering the associated evolution equation, or equivalently, the gradient descend method:

                           ∂u ∂t -∂ ∂x   u x u 2 x + u 2 y   -∂ ∂y   u y u 2 x + u 2 y   + (u I 2 x + v I x I y + I x I t ) = 0 ∂v ∂t -∂ ∂x   v x v 2 x + v 2 y   -∂ ∂y   v y v 2 x + v 2 y   + (u I x I y + v I 2 y + I y I t = 0 + Boundary conditions + Initial estimation : u(0, x, y) = u 0 (x, y) v(0, x, y) = v 0 (x, y) (7)
In [START_REF] Cohen | Nonlinear variational method for optical flow computation[END_REF], a solution of system 7 is described by a finite element method. This solution is done in two steps, first the temporal discretization of Eq. 7 and then the spatial discretization. This can be done since the temporal variable t and the space variables (x, y) are independent.

Motion Interpretation

As explained in the introduction, our aim is to make use of the optical flow, we described in the previous section, to localize interesting regions that have to be further modelized. The extraction of higher level description is necessary to characterize different types of movement.

Different studies have been done on linear phase portrait [START_REF] Ford | Image models for 2-D flow visualization and compression[END_REF][START_REF] Shu | Vector field analysis for oriented patterns[END_REF][START_REF] Rao | Computerized flow field analysis: Oriented textures fields[END_REF] and their use for characterizing oriented texture field. Zhong et al [START_REF] Zhong | Salient structure analysis of fluid flow[END_REF] propose an application of linear phase portrait interpretation of structures in a fluid flow. The main drawback of this method is that it can handle only one critical point. This is inconsistent for our application as we can observe different vortices on the data. Recently Ford and Strickland [START_REF] Ford | Nonlinear phase portrait models for oriented textures[END_REF] proposed a nonlinear phase portrait model allowing multiple critical points, but the model seems computationally expensive, and not very stable.

In this paper, we propose a new method for approximating an orientation field and for characterizing stationary points of trajectories obtained from an arbitrary phase portrait.

Approximation by a polynomial phase portrait model

Suppose that we wish to approximate our optical flow field f = (u, v) by a polynomial phase portrait model

g(x, y) =      dx dt = P (x, y) dy dt = Q(x, y) (8)
where P and

Q ∈ Q n (IR 2 ) = {p, st p(x, y) = i,j≤n a ij x i y j }.
We fit this model into the original field by minimizing locally an energy function

S(g) = 1 2 k,l∈W |f × g| 2 (9)
where f × g denotes the cross product of the two vectors f and g and (k, l) ∈ W is the considered neighborhood.

In the following we show that recovering the coefficients of the two polynomials P and Q by minimizing S amounts to an eigenvalue problem. This approach represents a major advantage since it leads to a linear problem independing of the chosen polynomials P and Q.

Let X = (x i y j ) i,j=0...n be a basis of Q n (IR 2 ), we can write:

g(x, y) =        P (x, y) = i,j≤n a i,j x i y j = Ω T a X Q(x, y) = i,j≤n b i,j x i y j = Ω T b X (10)
S can be rewritten as:

S(g) = 1 2 W (Ω T a Xf 2 -Ω T b Xf 1 ) 2 (11)
Let B = (Xf 2 , -Xf 1 ) T and L = (Ω a , Ω b ), we have:

S(g) = S(L) = 1 2 W (BL) 2 (12)
To obtain the unicity of the solution we add a normalization constraint L T L = 1 and finally the similarity measure between the optical flow field f = (u, v) and the model g = (Ω a X, Ω b X) is obtained by the minimization under constraint of S(L):

S(L) = 1 2 W L T B T BL + 1 2 λ(L T L -1) (13)
A minimum is characterized by the derivatives equations:

     ∂S ∂L = B T BL + λ L = 0 ∂S ∂λ = 1 2 (L T L -1) = 0 (14)
which lead to the following eigenvector problem:

-B T BL = λ L (15)
along with the normalization constraint:

L T L = 1 (16)
The vector L, representing the coefficients of the model in a polynomial basis, is characterized as the eigenvector associated to the largest eigenvalue of the matrix -B T B.

Interpretation

After approximation of the optical flow field by a polynomial flow field, it becomes possible to characterize some high-level information in the vector field.

In the following we report the classification of fixed points for a linear phase model and in the general case of polynomial models. Linear phase models can be used in the neighborhood of fixed points of arbitrary polynomial model as a first order approximation. Indeed, the behavior is the same as for the first order approximation (i.e. linear case) except that when the characteristic roots are pure complex, we have a center or a focus [START_REF] Lefschetz | Differential Equations: Geometric Theory[END_REF]. Nevertheless, the experiments shown in this paper were obtained with a third order polynomial model (i.e Q 3 ).

Linear case

Consider a linear phase model where the polynomials P and Q are in P 1 (IR 2 ) = {p st p(x, y) = a 00 + a 10 x + a 01 y}.

The vectors base may be written X = (1, x, y) and the model may be rewritten: , and we obtain the following classification:

• if λ 1 and λ 2 are real:

-if λ 1 and λ 2 are distinct and of the same sign, then x f is a node,

-if λ 1 = λ 2 < 0 or λ 1 = λ 2 > 0, then x f is a degenerate node, -if λ 1 λ 2 < 0, then x f is a saddle point, -if λ 1 = λ 2 = 0, then x f is a singular node,
• if λ 1 and λ 2 are complex:

-if Re(λ 1 ) = 0, then x f is a center, -if Re(λ 1 ) < 0, then x f is a stable focus,
-if Re(λ 1 ) > 0, then x f is an unstable focus.

In the case of a polynomial model, the previous linear approximation holds in the neighborhood of a fixed point and hence can be used for its classification but it requires the computation of the fixed points i.e. the solution of a set of polynomial equations of arbitrary degree. Instead we make use of the index of a flow field.

Index of a flow field

Let F = (P, Q) be a vector field defined over a Jordan curve J in the Euclidean plane, with no critical point on J. The index of F over J is proportional to the angular variation of the vector F (M ) (applied at M ∈ J) as M describes J. For the system (??), the index over an oriented Jordan curve J is given by:

Index(J) = 1 2π J d arctan Q P = 1 2π J P dQ -QdP P 2 + Q 2 . ( 18
)
A classification of the flow field F can be obtained by computing the index over a small circle surrounding an isolated critical point. Since it is hard to compute the stationary points of the system (??), it leads to solve a system of polynomials with arbitrary degree. We have chosen to compute locally the index of F = (P, Q) over the whole flow field. At each point we consider a circle contained in the centered window W and we use the following classification:

• The index of a focus, a center or a node is equal to +1,

• The index of a saddle point is -1.

Although this characterization is compendious, it characterizes the most important structures in a fluid flow field. The index measure computed over all the flow allows us also to obtain critical points locations without computing the roots of the system (??). Once we located these points we may use the linearization of polynomial model at the neighborhood of fixed points to obtain a complete description of the flow field in these points.

In this approach, we have proposed a two step method based on first computing an optical flow field preserving motion discontinuities and then fit locally a polynomial phase model to the obtained field. We could have used a polynomial model for flow computation [?] and its characterization. Using such a model has some limitations. Indeed, flow discontinuities and non uniform tessellation which are two important properties in oceanographic image sequences, cannot be used with a polynomial model.

Vortex segmentation and quantification

Optical flow and polynomial phase portrait methods allows to characterize some regions of interest. These regions correspond to node and focus points of the phase portrait model and we will modelize them as vortices in the following.

This modelization is based on the following properties of vortices:

• The region including the vortex is the result of a cold and a warm front meeting: so it is characterized by an important variation of gray level values,

• Edge detection within the image shows that a lot of points on the structure's boundary have a high gradient norm value,

• Variance calculus with application of hysteresis thresholding, allows to localize vortex rolls,

• Geometrically, the external shape of the vortex may be approximated by an arc of a circle.

Segmentation of the vortex rolls

The rolls of the vortex are small circle-like regions characterized by an important variation of gray level values. We have chosen to extract them by using the variance image, which is defined for each pixel of the data image, by the variance of a window centered around that point. A window of radius 10 pixels leads to fair results regarding to the size of the studied phenomena. Regions of interest are then extracted by hysteresis thresholding of that image. Thresholding values are automatically computed as corresponding to the classical 95 percent threshold on the variance histogram.

After this process, we still have to eliminate remaining false regions. We then use the fact that the regions we are searching for, have a small extension and have approximately the shape of a circle -or an ellipse with eccentricity close to -1-. This is done by computing principal directions of the region, and its extensions along those directions. A simple thresholding on the ratio of these extensions is enough to retain only good regions.

Segmentation Model

We modelize vortices by computing their center regions and by an arc of a circle delimiting their extension. We there use a half of a circle. We then have to determinate, firstly, the center and the orientation of this arc. This depends on the type of vortex and we distinguish two main cases :

• Symmetric case : the vortex is made of two rotating regions evolving in the same way. The extreme points of the arc will be on the axis between centers of gravity of the two regions, and the center will be the center of symmetry of the regions.

• Other cases : only one vortex has evolved enough. We then center the arc over the center of these regions.

The arc must be such that gradient norm is maximal along it. We express that with the help of the following energy functional :

E(r) = 1 2r ∇I(x(r), y(r)) 2 (19)
where r represents the radius of the circle's arc of our model.

Application on oceanographic images

Image processing tools have been applied quite successfully these past few years in the domain of land cartography. The shapes to be extracted are numerous and complex (roads, rivers, airports,..). On the other hand, oceanographic "shapes" have received very few attention up to now. Structures such as fronts and vortices can often be found on sea surface temperatures images (AVHRR). These structures may vary a lot with time, and therefore a specific process is required to detect them. Several problems arise in vortices detection and tracking :

• AVHRR images represent a huge amount of data to process daily. We then need to localize area of interest within a wide image, in order to work on a smaller one. This suppose that a vortex is never missed, and also that we minimize the number of false hits;

• a mathematical framework for segmenting the structure within the smaller image is required ;

• at last, we need a modelization of structure's shape that allows wide deformations.

This paper has focused on the first and the third points above. We will now present some results obtained on temperature and CZCS images .

We have started our process by a optical flow computation over an AVHRR image sequence. Figure 2 shows an overlay of the obtained displacement field on a frame. Some particular flow patterns can be seen over the whole image. We can distinguish node points, vortex points and saddle points. The purpose of the motion interpretation explained in section 3 is to localize automatically those regions. We have used a cubic polynomial phase portrait to classify the different flow patterns. Figure 3 shows the obtained regions of interest, where grey regions represent nodes and white regions, vortices.

The geometrical modelization is then applied on a small window centered on the region classified as vortices points. The maximization of the energy function described in equation 19 yields on a CZCS image (figure 4), an efficient representation of the vortex. An overlay of the model on the image data is shown in figure 5.

Finally, we have applied our geometric model on CZCS image sequence to characterize the evolution of the spatio-temporal deformable structure. This variation involves position, extension and geometry (and even topology, if two vortices merge for instance) as can be seen in figure 6 showing three successive frames of the image sequence. Those structures have a temporal evolution which depends on their location, size and the period of the year. We expect to use the obtained geometric model to characterize their evolution.

Conclusion

This paper has presented a global methodology to localize, segment and track vortices on oceanographic images. We make use of some information on the physical processes occurring on the images. In the future we plan to:

• reconsider the validity of optical flow equation for different applications :

would it be possible to make use of a different assumption on the variation of the gray level value?

• compare image processing and geophysical modelling for vortices detection and tracking in order to study if image processing of dynamic satellite image may be an alternative way of processing these data for a daily analysis. • consider the problem of data fusion. Information coming from Infra-red data (AVHRR) and color data (CZCS) are highly correlated but contrast and statistical properties of grey level values vary from one sensor to another. Moreover, the stability of the process is different for the different types of data.

Figure 1 :

 1 Figure 1: Left: Flow field on a temperature front. Right: Optical flow computed with a continuity constraint.

  + a 10 x + a 01 y = (a 00 , a 10 , a 01 )X dy dt = b 00 + b 10 x + b 01 y = (b 00 , b 10 , b 01 )X (17) Once we have recovered at each point the coefficients of the linear model a fixed point x f can be classified according to the eigenvalues, λ 1 and λ 2 , of the matrix A = a 10 a 01 b 10 b 01

Figure 2 :

 2 Figure 2: An overlay of the optical flow field and the corresponding AVHRR image.

Figure 3 :

 3 Figure 3: Interpretation of optical flow by a cubic polynomial phase portrait. Gray regions represent nodes while white regions represent vortices.

Figure 4 :

 4 Figure 4: A CZCS image of the western french coast (given by the OCEAN project).

Figure 5 :

 5 Figure 5: Vortex model applied to a CZCS image after localizing vortex rolls and fitting the arc representing vortex boundaries by maximizing Eq. 19.

Figure 6 :

 6 Figure 6: Results obtained on three consecutive CZCS images. An overlay of the obtained model is displayed on each frame.

ACKNOWLEDGEMENTS

We thank ACRI company for providing us AVHRR images under a research contract, and also for giving us information about the data and oceanic structures. The OCEAN (Ocean Color European Archive Network) project is the data source of the CZCS images.


Charles Huyghues-Despointes

Binh-Minh Bui-Xuan

Clémence Magnien

Clémence Magnien Forte

Forte ∆-connexité dans les flots de liens

Keywords: flots de liens, connexité, algorithmique

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

Introduction

The question of connectivity in graphs is well-defined and understood: it is possible to partition a graph's nodes into connected components, such that there is a path from any node to any other node of the same component. Moreover, the computation of this partition can be done in linear time in the size of the graph, with simple algorithms such as a breadth-first search.

In the case where we do not deal with a classical, static graph, but where the links evolve with time, we are faced with two problems. First, though it is straightforward to extend the classical notion of path to that of a temporal path, the corresponding reachability relation is not transitive. This means that the natural extension of the notion of connected components will not lead to a node partition, and instead the components will overlap. Second, computing these connected components is an NP-hard problem [START_REF] Bhadra | Complexity of connected components in evolving graphs and the computation of multicast trees in dynamic networks[END_REF].

The question of connectivity in the dynamic case has been addressed [BF03, NTM + 12, TMML10]. In these cases, a strongly connected component is defined as a set of nodes such that there exist a temporal path from any node to any other node in the set.

However, temporal paths may span a small or large fraction of the total duration of the considered link streams. Hence, the fact that a link stream is strongly connected can describe vastly different situations. For instance, there could exist a small time interval such that there are temporal paths from any node to any other node in this interval; in this case no description of the dataset outside of this time interval is given; or the temporal paths may have a time length comparable to that of the whole dataset; or other possibilities.

In order to provide a better description of a link stream with respect to connectivity, we are interested in finding stretches of time during which paths exist recurrently. In an email network, this would for instance amount to: every day, an information can be iteratively forwarded from any user X to any user Y of usergroup group. This corresponds to the notion of strongly ∆-connected components, introduced in [GCCLL15] as ∆-components.

In this paper, when addressing evolving network data, we will capture timestamped links between nodes under the formalism of a link stream. A link stream is a triple L = (T,V, E) where T = [0,t max ] represents a time interval, V is a finite set of nodes, and

E ⊆ {(t, uv) : t ∈ T ∧ u ∈ V ∧ v ∈ V ∧ u = v} is a set of
timestamped unordered pairs of nodes called timelinks. A timelink is a pair (t, uv) where t is a time instant and uv is a notation which represents the unordered link uv = vu between two distinct nodes u and v † .

We show that computing the maximal strongly ∆-connected components is a NP-hard problem, and we propose algorithms for computing upper bounds and lower bounds on their size.

2 Bounding strongly ∆-connected components in link streams Let L = (T,V, E) be a link stream where T = [0,t max ]. For any subset S ⊆ V of nodes we call L[S] = (T, S, E[S]) the substream of L induced by S, where

E[S] = E ∩ {(t, uv) : t ∈ T ∧ u ∈ S ∧ v ∈ S}.
Definition 1 (Strongly ∆-connected component) A temporal path from u to v, starting at t 0 and at t k , for k ≥ 0, is a sequence (t 0 , x 0 x 1 ), (t 1 , x 1 x 2 ), . . . , (t k , x k x k+1 ) of timelinks of E where t 0 < t 1 < . . . < t k , x 0 = u and x k+1 = v. Link stream L is strongly ∆-connected if for any pair of distinct nodes u and v, for any time instant t ∈ [0,t max -∆], there exists a temporal path from u to v starting after t and ending before t + ∆. A subset S ⊆ V of nodes is strongly ∆-connected if the induced substream L[S] is strongly ∆-connected, and it is called a component when maximal by inclusion.

For instance, consider the link stream L s = (T s ,V s , E s ) where T s = [0, 3], V s = {a, b, c, d, e}, and E s = {(0, ab), (1, bc), (1, be), (2, cd), (2, ed), (3, ad)}. For ∆ = 3, there are only three missing temporal paths: from c to e; from e to c; and from d to b. Hence, there are exactly four strongly ∆-connected components of L s , which are {a, b, c}, {a, b, e}, {a, c, d}, and {a, d, e}. They all overlap each other in a non-trivial way.

We claim that finding the largest strongly ∆-connected component of link stream L is NP-hard. The best way to see this is to consider the ∆-reachability digraph of L.

Definition 2 (∆-reachability digraph) The ∆-reachability digraph of link stream L = (T,V, E) is the di- graph D L = (V, A)
where an arc (u, v) belongs to the arc set A if and only if for any time instant t ∈ [0,t max -∆] there exists a temporal path from u to v starting after t and ending before t + ∆.

By extension, if S ⊆ V , we denote by D S the ∆-reachability subgraph of L[S].

Strongly ∆-connected components of link stream L are exactly maximal cliques of digraph D L . We can then prove the NP-hardness by making a reduction from the problem of cliques in classical graphs to the problem of cliques in ∆-reachability digraphs, in a similar way as what has been done in [START_REF] Bhadra | Complexity of connected components in evolving graphs and the computation of multicast trees in dynamic networks[END_REF].

Lower bounding strongly ∆-connected components

We are interested in the size of the largest strongly ∆-connected component. We present here a partition refinement approach for finding large strongly ∆-connected sets, which provides a lower bound for the size of the largest component.

Algorithm 1 Input: link stream L = (T,V, E). First partition V into strongly connected components of its ∆-reachability digraph D L . For each partition part S that is not strongly ∆-connected in L, successively perform both following splitting operations. First split S into two parts S \ {u} and {u} following some RULE. Second, split S \ {u} into strongly connected components of the ∆-reachability digraph D S\{u} . Keep on refining as long as there are partition parts that are not strongly ∆-connected. After refinement stability, scan through the partition parts and output their maximum size as a lower bound of the maximum strongly ∆-connected component of link stream L.

Claim 1 Algorithm 1 first computes a partition of V into strongly ∆-connected sets of L before outputting the maximum size of these sets, hence, correctly outputs a lower bound of the maximum strongly ∆connected component of link stream L.

Proof: The main idea here is that cliques of digraphs are subsets of their strongly connected component. 2 This process is an inversely incremental approach whose efficiency depends heavily on the step when we split S into S \ {u} and {u} according to RULE. The exact choice of RULE will be discussed below. Concerning complexity issues, it would be time consuming to compute the many ∆-reachability digraphs † As in most timestamped log data, we suppose that the input of timelink set E is already sorted by time instants. involved in the process. However, we are not really interested in these digraphs, except for their strongly connected components. We define the (undirected) underlying graph of link stream L = (T,V, E) as G = (V, E) where E = {uv : ∃t ∈ T, (t, uv) ∈ E}.

Claim 2 For any link stream, the strongly connected components of its ∆-reachability digraphs are exactly the connected components of its underlying graph.

Proving this claim is a straightforward exercise. Note that computing connected components of an undirected graph can be done in linear time, and the graph itself can be obtained implicitly from the corresponding link stream.

We now address the test condition whether a partition part S is strongly ∆-connected in link stream L. By lack of space, we skip the details and claim that it is possible to modify the algorithm computing temporal distance vector in [START_REF] Kossinets | The Structure of Information Pathways in a Social Communication Network[END_REF] and devise an algorithm deciding strong ∆-connectedness.

Eventually, we are left with defining the RULE for choosing the node u ∈ S in the splitting of S into S \ {u} and {u}. We pick u as the node of the ∆-reachability digraph of L[S] having the smallest sum of in and out degrees. We experimented with other criteria and this one provided the best results. We leave a precise evaluation of the impact of this RULE for future work. It is possible to find this node by modifying the above-mentioned variant of Algorithm [START_REF] Kossinets | The Structure of Information Pathways in a Social Communication Network[END_REF]. The overall complexity is therefore at most that of the routine for testing whether S is strongly ∆-connected in L.

Refining the timeline and upper bounding strongly ∆-connected components

Link streams may or may not be strongly ∆-connected, and this may evolve with time (and the sizes of components may also evolve). For instance, professional email exchanges are not to be expected during weekends and holidays. It is therefore important to be able to compute time intervals where some special phenomenon occurs, e.g. when the network records heavy activities or, on the contrary, a very quiet period. Formally, for any time interval I ⊆ T we call L[I] = (I,V, E[I]) the substream of L induced by I, where

E[I] = E ∩ {(t, uv) : t ∈ I}.
Here we cover the timeline with overlapping time intervals using the formalism of ∆-existence.

Definition 3 (∆-existence) A node u ∈ V of link stream L = (T,V, E) is said to ∆-exist at time instant t ∈ T if there exists a node v ∈ V and a time instant t ∈ [t,t + ∆] such that (t , uv) ∈ E. We say that an interval The computation of ∆-existence of a node can be done by scanning the link stream. From that computation, outputting consistent intervals (and hence candidate intervals) is a straightforward exercise.

I = [a, b] of time I ⊆ T is a consistent interval if
Finally, for every candidate interval I, the number of nodes which ∆-exist throughout I is an upper bound of the size of the maximum strongly ∆-connected component of L[I].

Application to the Rollernet dataset

We considered the Rollernet dataset [TLD + 09], collected during a rollerblade tour in Paris. 62 participants were equipped with wireless sensors recording when they are at a communication distance from one another. The total dataset duration is approximately 2 hours and 45 minutes (with a break of approximately 30 minutes). The candidate intervals, as well as the corresponding upper and lower bounds for the size of the largest strongly ∆-connected component are presented in Figure 1. The computation for each value of ∆ took just under 3s on a laptop.

We observe a very long candidate interval (minutes 20 to 65) where the bounds show that the maximum strongly ∆-connected component includes over 90% of the nodes, for both values of ∆. The following candidate interval also lasts for a long period of time (> 25 minutes). For ∆ = 6 minutes the bounds again indicate that the largest component spans over 90%of the nodes. For ∆ = 4 minutes, the lower bound shows that the largest component includes over 75% of the total number of nodes for the same interval. We believe that this bound can be improved to be closer to the upper bound, and leave this for future work.

Our results strongly suggest that, after a short initial phase (20 first minutes), the Rollernet network records dense and regular (every 5 minutes) exchanges for a long lasting period, until a drop of activity (between 90th and 120th minute, probably corresponding to the break) after which point a regular communication density can resume again.

Conclusion

We presented an algorithm for computing upper and lower bounds for the size of the largest strongly ∆connected components. We also introduced a method for splitting a linkstream's timeline into relevant subintervals. Illustrated on a real-world dataset, our method highlights some interesting features and allows a better understanding of the dataset.

Future work includes the computation of the exact time complexity of the algorithm, as well as formal results on the bounds' tightness. We think it is possible to improve the tightness of the lower bound. Finally we aim at making more detailed analysis of this and other datasets.

  every node of the link stream satisfies one of the following: either the node ∆-exists at all time instant t ∈ [a, b]; or, the node does not ∆-exist at any time instant t ∈ [a, b]. Definition 4 (Candidate interval) For any consistent interval I = [a, b] as described above, the candidate interval associated to I is the time interval [a, b + ∆]. Notice that, as consistent intervals form a partition of the timeline, candidate intervals may overlap.

Fig. 1 :

 1 Fig. 1: Candidate time interval covering of Rollernet, with upper and lower bounds of the size of the maximum strongly ∆-connected components for ∆ = 6 minutes (left) and ∆ = 4 minutes (right). Red (resp. blue) segments correspond to lower (resp. upper) bounds during the corresponding candidate interval. The middle points of the segments are joined for a better interpretation.